diff --git a/cscripts/Makefile.am b/cscripts/Makefile.am
index b623a23..6712daf 100644
--- a/cscripts/Makefile.am
+++ b/cscripts/Makefile.am
@@ -16,6 +16,7 @@ pythonscripts = __init__.py \
cssens.py \
csspec.py \
cstsdist.py \
+ cstsmapsplit.py \
cstsmapmerge.py \
csworkflow.py \
csfindobs.py \
@@ -38,6 +39,7 @@ parfiles = cscaldb.par \
cssens.par \
csspec.par \
cstsdist.par \
+ cstsmapsplit.par \
cstsmapmerge.par \
csworkflow.par \
csfindobs.par \
@@ -60,6 +62,7 @@ linkscripts = cscaldb \
cssens \
csspec \
cstsdist \
+ cstsmapsplit \
cstsmapmerge \
csworkflow \
csfindobs \
@@ -84,6 +87,7 @@ testscripts = tests/__init__.py \
../test/test_cssens.py \
../test/test_csspec.py \
../test/test_cstsdist.py \
+ ../test/test_cstsmapsplit.py \
../test/test_cstsmapmerge.py \
../test/test_csworkflow.py \
../test/test_csfindobs.py \
diff --git a/cscripts/__init__.py.in b/cscripts/__init__.py.in
index 7d299d2..36ced7a 100644
--- a/cscripts/__init__.py.in
+++ b/cscripts/__init__.py.in
@@ -38,6 +38,7 @@ __all__ = [
"cssens",
"csspec",
"cstsdist",
+ "cstsmapsplit",
"cstsmapmerge",
"csworkflow",
"obsutils"
@@ -63,6 +64,7 @@ from cscripts.csroot2caldb import csroot2caldb
from cscripts.cssens import cssens
from cscripts.csspec import csspec
from cscripts.cstsdist import cstsdist
+from cscripts.cstsmapsplit import cstsmapsplit
from cscripts.cstsmapmerge import cstsmapmerge
from cscripts.csworkflow import csworkflow
from cscripts import obsutils
diff --git a/cscripts/cstsmapsplit.par b/cscripts/cstsmapsplit.par
new file mode 100644
index 0000000..db7065c
--- /dev/null
+++ b/cscripts/cstsmapsplit.par
@@ -0,0 +1,55 @@
+#############################################################################
+# #
+# CTA Science Analysis Tools #
+# #
+# ------------------------------------------------------------------------- #
+# #
+# File: cstsmapsplit.par #
+# Version: 1.1.0 #
+# #
+# Author: Michael Mayer #
+# michael.mayer@physik.hu-berlin.de #
+# HU Berlin #
+# #
+# Purpose: Parameter file for the cstsmapsplit application #
+# #
+#############################################################################
+#
+# General parameters
+#===================
+inobs, f, a, events.fits,,, "Input event list, counts cube or observation definition XML file"
+inmodel, f, a, $CTOOLS/share/models/crab.xml,,, "Input model XML file"
+srcname, s, a, Crab,,, "Test source name"
+expcube, f, a, NONE,,, "Input exposure cube file (only needed for stacked analysis)"
+psfcube, f, a, NONE,,, "Input PSF cube file (only needed for stacked analysis)"
+edispcube, f, a, NONE,,, "Input energy dispersion cube file (only needed for stacked analysis)"
+bkgcube, f, a, NONE,,, "Input background cube file (only needed for stacked analysis)"
+caldb, s, a, prod2,,, "Calibration database"
+irf, s, a, South_0.5h,,, "Instrument response function"
+edisp, b, h, no,,, "Apply energy dispersion?"
+outmap, f, a, tsmap.fits,,, "Output Test Statistic map file"
+bins_per_job, i, a, 5,,, "Number of TS map bins per running job"
+compute_null, b, a, yes,,,"Compute null hypothesis first?"
+outfile, f, a, commands.txt,,, "ASCII file containing all commands"
+
+#
+# Spatial binning parameters
+#===========================
+usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?"
+nxpix, i, a, 200,,, "Size of the X axis in pixels"
+nypix, i, a, 200,,, "Size of the Y axis in pixels"
+binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)"
+coordsys, s, a, CEL, CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
+xref, r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)"
+yref, r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)"
+proj, s, a, CAR, AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method"
+
+#
+# Standard parameters
+#====================
+publish, b, h, no,,, "Publish Test Statistic map on VO Hub?"
+chatter, i, h, 2,0,4, "Chattiness of output"
+clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
+debug, b, h, no,,, "Debugging mode activated"
+mode, s, h, "ql",,, "Mode of automatic parameters"
+logfile, f, h, cstsmapsplit.log,,, "Log filename"
diff --git a/cscripts/cstsmapsplit.py b/cscripts/cstsmapsplit.py
new file mode 100755
index 0000000..135ca08
--- /dev/null
+++ b/cscripts/cstsmapsplit.py
@@ -0,0 +1,309 @@
+#!/usr/bin/env python
+# ==========================================================================
+# Create commands to split TS map computation
+#
+# Copyright (C) 2016 Michael Mayer
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+#
+# ==========================================================================
+import sys
+import os
+import math
+import gammalib
+import ctools
+
+
+# ================== #
+# cstsmapsplit class #
+# ================== #
+class cstsmapsplit(ctools.cscript):
+ """
+ Generates commands to split TS map computation
+
+ This class implements the creation of file containing commands to split
+ the computation of a TS map into several machines. It derives from
+ the ctools.cscript class which provides support for parameter files,
+ command line arguments, and logging. In that way the Python script
+ behaves just as a regular ctool.
+ """
+
+ # Consructor
+ def __init__(self, *argv):
+ """
+ Constructor
+ """
+
+ # Set name and version
+ self._name = 'cstsmapsplit'
+ self._version = '1.1.0'
+
+ # Set data members
+ self._outmap = gammalib.GFilename()
+ self._bins_per_job = 0
+ self._compute_null = False
+ self._outfile = gammalib.GFilename()
+ self._map = gammalib.GSkyMap()
+ self._cmd = []
+
+ # Initialise observation container from constructor arguments.
+ self._obs, argv = self._set_input_obs(argv)
+
+ # Initialise application by calling the appropriate class
+ # constructor.
+ self._init_cscript(argv)
+
+ # Return
+ return
+
+
+ # Private methods
+ def _get_parameters(self):
+ """
+ Get parameters from parfile and setup the observation
+ """
+
+ # If there are no observations in container then get some ...
+ if self._obs.size() == 0:
+ self._obs = self._get_observations()
+
+ # Set models if we have none
+ if self._obs.models().size() == 0:
+ self._obs.models(self['inmodel'].filename())
+
+ # Get TS map paramters
+ self._map = self._create_map(self._obs)
+ self._srcname = self['srcname'].string()
+ self._outmap = self['outmap'].filename()
+
+ # Get additional parameters
+ self._bins_per_job = self['bins_per_job'].integer()
+ self._compute_null = self['compute_null'].boolean()
+ self._outfile = self['outfile'].filename()
+
+ # Write input parameters into logger
+ if self._logTerse():
+ self._log_parameters()
+ self._log('\n')
+
+ # Return
+ return
+
+
+ def _compute_null_hypothesis(self):
+ """
+ Compute null hypothesis
+
+ Returns
+ -------
+ logl : float
+ Log-likelihood of null hypothesis
+ """
+
+ # Store original models
+ models_orig = self._obs.models()
+
+ # Get models instance
+ models0 = self._obs.models()
+
+ # Remove test source
+ models0.remove(self._srcname)
+
+ # Initialise optimizer
+ opt = gammalib.GOptimizerLM()
+
+ # Set models for null hypothesis
+ self._obs.models(models0)
+
+ # Optimise observation container
+ self._obs.optimize(opt)
+
+ # Set original models to container
+ self._obs.models(models_orig)
+
+ # Return optimised log-likelihood of null hypothesis
+ return -(opt.value())
+
+
+ def run(self):
+ """
+ Run the script
+ """
+
+ # Switch screen logging on in debug mode
+ if self._logDebug():
+ self._log.cout(True)
+
+ # Get parameters
+ self._get_parameters()
+
+ # Write information into logger
+ if self._logTerse():
+ self._log('\n')
+ self._log.header1('Test source')
+ model = self._obs.models()[self._srcname]
+ self._log(str(model))
+ self._log('\n')
+
+ # Set log-likelihood to zero
+ logL0 = 0.0
+
+ # Pre-compute null hypothesis if requested
+ if self._compute_null:
+
+ # Write information into logger
+ if self._logTerse():
+ self._log('\n')
+ self._log.header1('Compute null hypothesis')
+
+ # Compute null hypothesis
+ logL0 = self._compute_null_hypothesis()
+
+ # Write likelihood into logger
+ if self._logTerse():
+ self._log(gammalib.parformat('Source removed'))
+ self._log(self._srcname)
+ self._log('\n')
+ self._log(gammalib.parformat('Log-likelihood'))
+ self._log(repr(logL0))
+ self._log('\n')
+
+ # Get parameters of TS map ctool
+ pars = gammalib.GApplicationPars('cttsmap.par')
+
+ # Compute total number of jobs
+ nbins = self._map.npix()
+ njobs = int(math.ceil(float(nbins) / float(self._bins_per_job)) + 0.1)
+
+ # Set parameters to be skipped now, we will deal with them later
+ skip_pars = ['binmin', 'binmax', 'logL0', 'outmap', 'logfile']
+
+ # Set tool name for computation
+ base_command = 'cttsmap'
+
+ # Write information into logger
+ if self._logTerse():
+ self._log('\n')
+ self._log.header1('Create commands')
+ self._log(gammalib.parformat('Number of cttsmap calls'))
+ self._log(str(njobs))
+ self._log('\n')
+
+ # Loop over TS map parameters
+ for par in pars:
+
+ # Skip if we deal with them later
+ if par.name() in skip_pars:
+ continue
+
+ # Skip if they need to be queried
+ # This way we ensure we pass only parameters
+ # that were queried before
+ if self[par.name()].is_query():
+ continue
+
+ # Set TS map parameter according to input from
+ # this script
+ par.value(self[par.name()].value())
+
+ # Append command to set parameter
+ base_command += ' ' + par.name() + '=' + par.value()
+
+ # Append null hypothesis parameter
+ base_command += ' logL0=' + repr(logL0)
+
+ # Set binning to start from zero
+ binmin = 0
+ binmax = 0
+
+ # Clear command sequence
+ self._cmd = []
+
+ # Loop over jobs and create commands
+ for job in range(njobs):
+
+ # Set specific outmap file name
+ outmap = self._outmap.url().replace('.fits', '_'+str(job)+'.fits')
+
+ # Set bin numbers to be computed in this job
+ binmin = binmax
+ binmax = binmin + self._bins_per_job
+ if binmax > nbins:
+ binmax = nbins
+
+ # Setup sliced command
+ sliced_command = '%s binmin=%s binmax=%s outmap=%s logfile=%s' % \
+ (base_command, str(binmin), str(binmax),
+ outmap, outmap.replace('.fits','.log'))
+
+ # Append command to list of commands
+ self._cmd.append(sliced_command)
+
+ # Write information into logger
+ if self._logExplicit():
+ self._log('\n')
+ self._log.header2('Commands')
+ for cmd in self._cmd:
+ self._log(cmd)
+ self._log('\n')
+ self._log('\n')
+
+ # Return
+ return
+
+
+ def execute(self):
+ """
+ Execute the script
+ """
+
+ # Run the script
+ self.run()
+
+ # Save results
+ self.save()
+
+ # Return
+ return
+
+
+ def save(self):
+ """
+ Save commands to ASCII file
+ """
+
+ # Open file
+ f = open(self._outfile.url(), 'w')
+
+ # Write commands to file
+ for cmd in self._cmd:
+ f.write(cmd + '\n')
+
+ # Close file
+ f.close()
+
+ # Return
+ return
+
+
+# ======================== #
+# Main routine entry point #
+# ======================== #
+if __name__ == '__main__':
+
+ # Create instance of application
+ app = cstsmapsplit(sys.argv)
+
+ # Execute application
+ app.execute()
diff --git a/doc/source/users/reference_manual/cstsmapsplit.rst b/doc/source/users/reference_manual/cstsmapsplit.rst
new file mode 100644
index 0000000..40eb71e
--- /dev/null
+++ b/doc/source/users/reference_manual/cstsmapsplit.rst
@@ -0,0 +1,127 @@
+.. _cstsmapsplit:
+
+cstsmapsplit
+============
+
+Generate commands to split the Test Statistic map computation
+
+
+Synopsis
+--------
+
+This script generates a sequence of commands to split the Test Statistic (TS)
+map computation into a distinct number of :ref:`cttsmap` computations. The
+sequence of commands is written into an ASCII file, where each line in the
+file corresponds to the command for an individual :ref:`cttsmap` run. Each TS
+map computed in that way will be suffixed with the command number, and the
+map can be combined into a single TS map using the :ref:`cstsmapmerge` script.
+
+
+General parameters
+------------------
+
+``inobs [file]``
+ Input event list, counts cube or observation definition XML file.
+
+``inmodel [file]``
+ Input model XML file.
+
+``srcname [string]``
+ Name of source model for which the TS map should be computed.
+
+``expcube [file]``
+ Input exposure cube file (only needed for stacked analysis).
+
+``psfcube [file]``
+ Input PSF cube file (only needed for stacked analysis).
+
+``bkgcube [file]``
+ Input background cube file (only needed for stacked analysis).
+
+``caldb [string]``
+ Calibration database.
+
+``irf [string]``
+ Instrumental response function.
+
+``(edisp = no) [boolean]``
+ Applies energy dispersion to response computation.
+
+``outmap [file]``
+ Output Test Statistic map file.
+
+``bins_per_job [integer]``
+ Number of bins of the TS map that should be computed within one command
+
+``compute_null [boolean]``
+ Pre-computes the null hypothesis globally and passes the value to the
+ individual commands. If the parameter is set to ``no``, the null
+ hypothesis will be computed in each run of :ref:`cttsmap` individually.
+
+``outfile [file]``
+ Output ASCII file name where the commands are written to.
+
+``(usepnt = no) [boolean]``
+ Use CTA pointing direction for map centre instead of xref/yref parameters?
+
+``nxpix [integer]``
+ Number of map pixels in Right Ascension or Galactic longitude.
+
+``nypix [integer]``
+ Number of map pixels in Declination or Galactic latitude.
+
+``binsz [real]``
+ Map pixel size (in degrees/pixel).
+
+``coordsys [string]``
+ Coordinate system (CEL - celestial, GAL - galactic).
+
+``xref [real]``
+ Right Ascension / Galactic longitude of map centre (J2000, in degrees).
+
+``yref [real]``
+ Declination / Galactic latitude of map centre (J2000, in degrees).
+
+``proj [string]``
+ Projection method.
+
+
+
+Standard parameters
+-------------------
+
+``(publish = no) [boolean]``
+ Specifies whether the TS map should be published on VO Hub.
+
+``(chatter = 2) [integer]``
+ Verbosity of the executable:
+ chatter = 0: no information will be logged
+
+ chatter = 1: only errors will be logged
+
+ chatter = 2: errors and actions will be logged
+
+ chatter = 3: report about the task execution
+
+ chatter = 4: detailed report about the task execution
+
+``(clobber = yes) [boolean]``
+ Specifies whether an existing output TS map file should be overwritten.
+
+``(debug = no) [boolean]``
+ Enables debug mode. In debug mode the executable will dump any log file
+ output to the console.
+
+``(mode = ql) [string]``
+ Mode of automatic parameters (default is "ql", i.e. "query and learn").
+
+``(logfile = cstsmapsplit.log) [string]``
+ Name of log file.
+
+
+Related tools or scripts
+------------------------
+
+:ref:`cttsmap`
+:ref:`cstsmapmerge`
+
diff --git a/doc/source/users/reference_manual/index.rst b/doc/source/users/reference_manual/index.rst
index 536ada0..b67ad95 100644
--- a/doc/source/users/reference_manual/index.rst
+++ b/doc/source/users/reference_manual/index.rst
@@ -49,8 +49,10 @@ cscripts
csresmap --- Generates residual map
cssens --- Computes CTA sensitivity
csspec --- Computes spectral points
- cstsdist --- Generates TS distribution
- cstsmapmerge --- Merges slices from ts map computations
+ cstsdist --- Generates Test Statistic distribution
+ cstsmapsplit --- Creates commands to split the Test Statistic map computations
+ cstsmapmerge --- Merges slices from Test Statistic map computations
+
Scripts to manage an IACT database
diff --git a/test/test_cstsmapsplit.py b/test/test_cstsmapsplit.py
new file mode 100755
index 0000000..d41fcb6
--- /dev/null
+++ b/test/test_cstsmapsplit.py
@@ -0,0 +1,182 @@
+#! /usr/bin/env python
+# ==========================================================================
+# This scripts performs unit tests for the cstsmapsplit script.
+#
+# Copyright (C) 2016 Michael Mayer
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program. If not, see .
+#
+# ==========================================================================
+import os
+import gammalib
+import cscripts
+
+
+# =============================== #
+# Test class for cstssplit script #
+# =============================== #
+class Test(gammalib.GPythonTestSuite):
+ """
+ Test class for cstssplit script
+
+ This test class makes unit tests for the csstsplit script by using it
+ from the command line and from Python.
+ """
+
+ # Constructor
+ def __init__(self):
+ """
+ Constructor
+ """
+ # Call base class constructor
+ gammalib.GPythonTestSuite.__init__(self)
+
+ # Return
+ return
+
+ # Set test functions
+ def set(self):
+ """
+ Set all test functions
+ """
+ # Set test name
+ self.name('cstsmapsplit')
+
+ # Append tests
+ self.append(self._test_cmd, 'Test cstsmapsplit on command line')
+ #self.append(self._test_python, 'Test cstsmapsplit from Python')
+
+ # Return
+ return
+
+ # Test cstsmapsplit on command line
+ def _test_cmd(self):
+ """
+ Test cstsmapsplit on the command line.
+ """
+ # Kluge to set the command (installed version has no README file)
+ if os.path.isfile("README.md"):
+ cstsmapsplit = "../cscripts/cstsmapsplit.py"
+ else:
+ cstsmapsplit = "cstsmapsplit"
+
+ # Setup cstsmapsplit command
+ cmd = cstsmapsplit+' inobs="selected_events.fits"'+ \
+ ' inmodel="data/crab.xml"'+ \
+ ' srcname="Crab"'+ \
+ ' caldb="prod2"'+ \
+ ' irf="South_0.5h"'+ \
+ ' outmap="tsmap.fits"'+ \
+ ' nxpix=5'+ \
+ ' nypix=5'+ \
+ ' binsz=0.05'+ \
+ ' coordsys="CEL"'+ \
+ ' xref=83.63'+ \
+ ' yref=22.01'+ \
+ ' proj="CAR"'+ \
+ ' compute_null=yes'+ \
+ ' bins_per_job=5'+ \
+ ' outfile="tsmap_cmds1.txt"'
+
+
+ # Execute cstssplit, make sure we catch any exception
+ try:
+ rc = os.system(cmd+" >/dev/null 2>&1")
+ except:
+ pass
+
+ # Check if execution was successful
+ self.test_assert(rc == 0,
+ "Successful cstsssplit execution on command line")
+
+ # Check observation definition XML file
+ self._check_cmdfile('tsmap_cmds1.txt')
+
+ # Setup cstsmapsplit command
+ cmd = cstsmapsplit+' inobs="selected_events.fits"'+ \
+ ' inmodel="data/crab.xml"'+ \
+ ' srcname="Does_not_exist"'+ \
+ ' caldb="prod2"'+ \
+ ' irf="South_0.5h"'+ \
+ ' outmap="tsmap.fits"'+ \
+ ' nxpix=5'+ \
+ ' nypix=5'+ \
+ ' binsz=0.05'+ \
+ ' coordsys="CEL"'+ \
+ ' xref=83.63'+ \
+ ' yref=22.01'+ \
+ ' proj="CAR"'+ \
+ ' compute_null=yes'+ \
+ ' bins_per_job=5'+ \
+ ' outfile="tsmap_cmds2.txt"'
+
+
+ # Execute cstssplit, make sure we catch any exception
+ try:
+ rc = os.system(cmd+" >/dev/null 2>&1")
+ except:
+ pass
+
+ # Check if execution failed
+ self.test_assert(rc != 0,
+ "Failure of cstssplit execution on command line")
+
+ # Return
+ return
+
+ # Test cstsmapsplit from Python
+ def _test_python(self):
+ """
+ Test cstsmapsplit from Python.
+ """
+ # Set-up cstsmapsplit
+ tsmapsplit = cscripts.cstsmapsplit()
+ tsmapsplit['inobs'] = 'obs_selected.xml'
+ tsmapsplit['inmodel'] = 'data/crab.xml'
+ tsmapsplit['xref'] = 83.6331
+ tsmapsplit['yref'] = 22.01
+ tsmapsplit['coordsys'] = 'CEL'
+ tsmapsplit['proj'] = 'CAR'
+ tsmapsplit['binsz'] = 0.05
+ tsmapsplit['nxpix'] = 5
+ tsmapsplit['nypix'] = 5
+ tsmapsplit['compute_null'] = False
+ tsmapsplit['bins_per_job'] = 5
+ tsmapsplit['outfile'] = 'tsmap_cmd_py1.txt'
+
+ # Run cstssplit script and save run list
+ tsmapsplit.logFileOpen() # Make sure we get a log file
+ tsmapsplit.run()
+ tsmapsplit.save()
+
+ # Check observation definition XML file
+ self._check_cmdfile('tsmap_cmd_py1.txt')
+
+ # Return
+ return
+
+ # Check observation definition XML file
+ def _check_cmdfile(self, filename):
+ """
+ Check ascii files containing commands
+ """
+
+ # Create file name instance
+ fname = gammalib.GFilename(filename)
+
+ # Check if execution was successful
+ self.test_assert(fname.exists(), 'Check of output filename exists')
+
+ # Return
+ return
diff --git a/test/test_python_cscripts.py b/test/test_python_cscripts.py
index bd8d68c..63c6150 100755
--- a/test/test_python_cscripts.py
+++ b/test/test_python_cscripts.py
@@ -37,6 +37,7 @@ import test_csroot2caldb
import test_cssens
import test_csspec
import test_cstsdist
+import test_cstsmapsplit
import test_cstsmapmerge
import test_csworkflow
import test_csiactcopy
@@ -111,6 +112,7 @@ def test(installed=False, debug=False):
suite_cssens = test_cssens.Test()
suite_csspec = test_csspec.Test()
suite_cstsdist = test_cstsdist.Test()
+ suite_cstsmapsplit = test_cstsmapsplit.Test()
suite_cstsmapmerge = test_cstsmapmerge.Test()
suite_csworkflow = test_csworkflow.Test()
@@ -128,6 +130,7 @@ def test(installed=False, debug=False):
suite_cssens.set()
suite_csspec.set()
suite_cstsdist.set()
+ suite_cstsmapsplit.set()
suite_cstsmapmerge.set()
suite_csworkflow.set()
@@ -145,6 +148,7 @@ def test(installed=False, debug=False):
suites.append(suite_cssens)
suites.append(suite_csspec)
suites.append(suite_cstsdist)
+ suites.append(suite_cstsmapsplit)
suites.append(suite_cstsmapmerge)
suites.append(suite_csworkflow)