Add csebins script (#1976)

Jürgen Knödlseder
1 parent 06ae0bbb
Showing 11 changed files with 619 additions and 6 deletions
ChangeLog
NEWS
cscripts/Makefile.am
cscripts/__init__.py.in
cscripts/csebins.par
cscripts/csebins.py
doc/source/users/reference_manual/csebins.rst
doc/source/users/reference_manual/cslightcrv.rst
doc/source/users/reference_manual/index.rst
test/test_csebins.py
test/test_python_cscripts.py
-2017-03-21
+2017-03-22
  
         * Version 1.3.0 released
         ========================
  
+        Add csebins script (#1976)
         Take into account RoI in ctskymap and use numerical integration for IRF (#1975)
         Set local language to english (#1370)
         Enhance logfile output of csroot2caldb script
 New Features and Important Changes in ctools 1.3.0.dev1
  
-21 March 2017
+22 March 2017
  
  
 Introduction
@@ -112,7 +112,10 @@ None
  
 cscripts
 --------
-None
+The csebins script has been added that generates an energy boundaries
+extension that can be used for a stacked analysis. The energy boundaries
+are set as function of the variation in the effective area and background
+rates of the Instrument Response Function (#1976).
  
  
 Examples
@@ -22,6 +22,7 @@
 # Scripts to be distributed in $(pyexecdir)/cscripts
 pythonscripts = $(top_builddir)/cscripts/__init__.py \
                 $(srcdir)/cscaldb.py \
+                $(srcdir)/csebins.py \
                 $(srcdir)/csinfo.py \
                 $(srcdir)/cslightcrv.py \
                 $(srcdir)/csmodelinfo.py \
@@ -53,6 +54,7 @@ pythonscripts = $(top_builddir)/cscripts/__init__.py \
  
 # Parfiles to be distributed in $(prefix)/syspfiles
 parfiles = $(srcdir)/cscaldb.par \
+           $(srcdir)/csebins.par \
            $(srcdir)/cslightcrv.par \
            $(srcdir)/csmodelinfo.par \
            $(srcdir)/csmodelmerge.par \
@@ -79,6 +81,7 @@ parfiles = $(srcdir)/cscaldb.par \
  
 # Scripts for which symbolic links will be created in bindir
 linkscripts = cscaldb \
+              csebins \
               csinfo \
               cslightcrv \
               csmodelinfo \
@@ -25,6 +25,7 @@ cur_version = sys.version_info
 # List all scripts in this module
 __all__ = [
     "cscaldb",
+    "csebins",
     "csinfo",
     "cslightcrv",
     "csmodelinfo",
@@ -59,6 +60,7 @@ if cur_version &gt; req_version:
 # Import
 from cscripts.cscaldb       import cscaldb
 from cscripts.csinfo        import csinfo
+from cscripts.csebins       import csebins
 from cscripts.cslightcrv    import cslightcrv
 from cscripts.csmodelinfo   import csmodelinfo
 from cscripts.csmodelmerge  import csmodelmerge
+#############################################################################
+#                                                                           #
+#                       CTA Science Analysis Tools                          #
+#                                                                           #
+# ------------------------------------------------------------------------- #
+#                                                                           #
+#  File:      csebins.par                                                   #
+#                                                                           #
+#  Author:    Juergen Knoedlseder                                           #
+#             Jurgen.Knodlseder@irap.omp.eu                                 #
+#             IRAP                                                          #
+#                                                                           #
+#  Purpose:   Parameter file for the csebins script                         #
+#                                                                           #
+#############################################################################
+#
+# Input and output parameters
+#============================
+inobs,   f, a, NONE,,, "Input observation definition XML file"
+caldb,   s, a, prod2,,, "Calibration database"
+irf,     s, a, South_0.5h,,, "Instrument response function"
+outfile, f, a, ebounds.fits,,, "Output energy boundaries file"
+
+#
+# Script parameters
+#==================
+emin,      r, a, 0.02,,, "Minimum energy (TeV)"
+emax,      r, a, 300.0,,, "Maximum energy (TeV)"
+aeffthres, r, h, 0.2,,, "Fractional effective area change threshold"
+bkgthres,  r, h, 0.5,,, "Fractional background rate change threshold"
+
+#
+# Standard parameters
+#====================
+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, csebins.log,,, "Log filename"
+#! /usr/bin/env python
+# ==========================================================================
+# Generates energy boundaries for stacked analysis
+#
+# Copyright (C) 2017 Juergen Knoedlseder
+#
+# 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 <http://www.gnu.org/licenses/>.
+#
+# ==========================================================================
+import sys
+import math
+import gammalib
+import ctools
+#from cscripts import obsutils
+
+
+# ============= #
+# csebins class #
+# ============= #
+class csebins(ctools.cscript):
+    """
+    Generates energy boundaries for stacked analysis
+    """
+
+    # Constructor
+    def __init__(self, *argv):
+        """
+        Constructor
+        """
+        # Set name and version
+        self._name    = 'csebins'
+        self._version = '1.0.0'
+
+        # Set members
+        self._ebounds = gammalib.GEbounds()
+
+        # Initialise observation container from constructor arguments.
+        self._obs, argv = self._set_input_obs(argv)
+
+        # Initialise script by calling the appropriate class constructor.
+        self._init_cscript(argv)
+
+        # Return
+        return
+
+
+    # Private methods
+    def _get_parameters(self):
+        """
+        Get parameters from parfile
+        """
+        # Initialise observation container if it is currently empty
+        if self._obs.size() == 0:
+
+            # Query "inobs" parameter
+            filename = self['inobs'].filename()
+
+            # If an observation definition file was provided then load it,
+            # otherwise build a single observation with response information
+            if self._is_valid_filename(filename):
+                self._obs = gammalib.GObservations(filename)
+            else:
+                cta   = gammalib.GCTAObservation()
+                caldb = gammalib.GCaldb('cta', self['caldb'].string())
+                rsp   = gammalib.GCTAResponseIrf(self['irf'].string(), caldb)
+                cta.response(rsp)
+                self._obs.append(cta)
+
+        # Query input parameters
+        self['emin'].real()
+        self['emax'].real()
+        self['aeffthres'].real()
+        self['bkgthres'].real()
+
+        # Query ahead output model filename
+        if self._read_ahead():
+            self['outfile'].filename()
+
+        #  Write input parameters into logger
+        self._log_parameters(gammalib.TERSE)
+
+        # Return
+        return
+
+    def _insert_energy(self, responses, energy, comment):
+        """
+        Set energy boundaries for one CTA observation
+
+        Parameters
+        ----------
+        responses : list of dict
+            List of response dictionaries
+        energy : float
+            Energy to insert (TeV)
+        comment : str
+            Reason for energy insertion
+        """
+        # Set thresholds of all response components
+        logE = math.log10(energy)
+        for rsp in responses:
+            rsp['thres'] = rsp['irf'](logE, 0.0, 0.0)
+
+        # Convert energy into gammalib energy
+        eng = gammalib.GEnergy(energy, 'TeV')
+  
+        # Set insertion message
+        msg = str(energy)+' TeV'
+        if len(comment) > 0:
+            msg += ' ('+comment+')'
+  
+        # If energy boundaries are empty then append an energy interval
+        # with zero size ...
+        if self._ebounds.size() == 0:
+
+            # Append energy interval with zero size
+            self._ebounds.append(eng, eng)
+
+            # Log first energy boundary
+            self._log_value(gammalib.NORMAL, 'First boundary', msg)
+
+        # ... otherwise append an interval
+        elif self._ebounds.size() > 0:
+
+            # If current boundaries have zero size then recover boundary energy,
+            # remove boundary, and append a boundary of non-zero size
+            if self._ebounds.emin() == self._ebounds.emax():
+
+                # Recover old energy, remove boundary and append new boundary
+                emin = self._ebounds.emin().copy() # Because emin becomes zero after remove
+                emax = eng
+                self._ebounds.remove(0)
+                self._ebounds.append(emin, emax)
+
+                # Log other energy boundary
+                self._log_value(gammalib.NORMAL, 'Boundary', msg)
+
+            else:
+
+                # Use old emax as new emin. Only append if interval is non zero
+                emin = self._ebounds.emax()
+                if emin < eng:
+                    self._ebounds.append(emin, eng)
+
+                    # Log other energy boundary
+                    self._log_value(gammalib.NORMAL, 'Boundary', msg)
+
+        # Return
+        return
+
+    def _set_ebounds(self, obs):
+        """
+        Set energy boundaries for one CTA observation
+
+        Parameters
+        ----------
+        obs : `~gammalib.GObservations`
+            CTA observation
+        """
+        # Get parameters
+        emin      = self['emin'].real()
+        emax      = self['emax'].real()
+        aeffthres = self['aeffthres'].real()
+        bkgthres  = self['bkgthres'].real()
+
+        # Build list of response dictionaries
+        responses = []
+        
+        # Loop over all observations
+        for run in obs:
+
+            # Get response name
+            mission    = run.response().caldb().mission()
+            instrument = run.response().caldb().instrument()
+            rspname    = run.response().rspname()
+            name       = '%s::%s::%s' % (mission, instrument, rspname)
+
+            # If response name exists already then examine next observation
+            exists = False
+            for rsp in responses:
+                if rsp['name'] == name:
+                    exists = True
+                    break
+            if exists:
+                continue
+
+            # Append effective area and background template to list
+            aeff = {'name': name, 'type': 'Aeff',
+                    'irf': run.response().aeff(), 'thres': 0.0}
+            bkg  = {'name': name, 'type': 'Background',
+                    'irf': run.response().background(), 'thres': 0.0}
+            responses.append(aeff)
+            responses.append(bkg)
+
+            # Log response name
+            self._log_value(gammalib.NORMAL, 'Append response', '%s (%s) [%s]' %
+                            (rspname, instrument, mission))
+
+        # Setup energy vector in log10(TeV)
+        logEs = [math.log10(1.0e-3*float(i)) for i in range(int(emin*1000),int(emax*1000))]
+
+        # Loop over all energies
+        for logE in logEs:
+
+            # Loop over all response components
+            for rsp in responses:
+
+                # Get IRF value
+                irf = rsp['irf'](logE, 0.0, 0.0)
+
+                # Get threshold
+                if rsp['type'] == 'Aeff':
+                    threshold = aeffthres
+                else:
+                    threshold = bkgthres
+
+                # If threshold is zero and response is positive then insert
+                # an energy ...
+                if rsp['thres'] == 0.0 and irf > 0.0:
+                    self._insert_energy(responses, math.pow(10.0, logE), rsp['type'])
+                    break
+
+                # ... otherwise if threshold is not zero then check whether
+                # the fractional change exceeds the threshold
+                elif rsp['thres'] != 0.0:
+                    f = irf/rsp['thres']
+                    if 1.0-f > threshold or f-1.0 > threshold:
+                        self._insert_energy(responses, math.pow(10.0, logE), rsp['type'])
+                        break
+
+        # Insert energy at end of range
+        self._insert_energy(responses, emax, 'End of range')
+
+        # Log number of energy boundaries
+        self._log_value(gammalib.NORMAL, 'Number of boundaries', self._ebounds.size())
+
+        # Return
+        return
+
+
+    # Public methods
+    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 observation into logger
+        self._log_observations(gammalib.NORMAL, self._obs, 'Input observation')
+
+        # Write header
+        self._log_header1(gammalib.TERSE, 'Define energy boundaries')
+
+        # Define energy boundaries for observation container
+        self._set_ebounds(self._obs)
+
+        # Return
+        return
+
+    def execute(self):
+        """
+        Execute the script
+        """
+        # Open logfile
+        self.logFileOpen()
+
+        # Read ahead output parameters
+        self._read_ahead(True)
+
+        # Run the scrip
+        self.run()
+
+        # Save energy boundaries
+        self.save()
+
+        # Return
+        return
+
+    def save(self):
+        """
+        Save energy boundaries
+        """
+        # Write header
+        self._log_header1(gammalib.TERSE, 'Save energy boundaries')
+
+        # Get energy boundaries filename
+        outfile = self['outfile'].filename()
+
+        # Log file name
+        self._log_value(gammalib.NORMAL, 'Energy boundaries file', outfile.url())
+
+        # Save energy boundaries
+        self._ebounds.save(outfile, self._clobber())
+
+        # Return
+        return
+
+    def ebounds(self):
+        """
+        Return energy boundaries
+
+        Returns
+        -------
+        ebounds : `~gammalib.GEbounds()`
+            Energy boundaries
+        """
+        # Return
+        return self._ebounds
+
+
+# ======================== #
+# Main routine entry point #
+# ======================== #
+if __name__ == '__main__':
+
+    # Create instance of application
+    app = csebins(sys.argv)
+
+    # Execute application
+    app.execute()
+.. _csebins:
+
+csebins
+=======
+
+Generates energy boundaries for stacked analysis.
+
+
+Synopsis
+--------
+
+This script generates the energy bounaries for a stacked analysis by inspecting
+the variation in the effective area and the background template of the response
+functions. The script takes either an observation definition XML file or a
+specific instrument response function on input.
+
+csebins writes an energy boundary extension into a FITS file.
+
+
+General parameters
+------------------
+
+``inobs [file]``
+    Input observation definition XML file.
+
+``caldb [string]``
+    Calibration database.
+ 	 	 
+``irf [string]``
+    Instrumental response function.
+
+``outfile [file]``
+    Name of the energy boundary output file.
+
+``emin [real]``
+    Lower energy limit of energy boundaries (in TeV).
+ 	 	 
+``emax [real]``
+    Upper energy limit of energy boundaries (in TeV).
+
+``(aeffthres = 0.2) [real]``
+    Fractional change in effective area that leads to insertion of a new energy
+    boundary.
+
+``(bkgthres = 0.5) [real]``
+    Fractional change in background rate that leads to insertion of a new energy
+    boundary.
+
+
+Standard parameters
+-------------------
+
+``(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 energy boundaries output 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 = csebins.log) [filename]``
+    Log filename.
+
+
+Related tools or scripts
+------------------------
+
+:doc:`ctbin`
+:doc:`ctexpcube`
+:doc:`ctpsfcube`
+:doc:`ctedispcube`
+:doc:`ctbkgcube`
@@ -46,7 +46,7 @@ General parameters
     Applies energy dispersion to response computation.
  
 ``outfile [file]``
-    Name of the source spectrum output file.
+    Name of the light curve output file.
  
 ``(tbinalg = GTI) <FILE|LIN|GTI> [string]``
     Algorithm for defining time bins.
@@ -125,7 +125,7 @@ Standard parameters
      chatter = 4: detailed report about the task execution
  
 ``(clobber = yes) [boolean]``
-    Specifies whether an existing source spectrum output file should be overwritten.
+    Specifies whether an existing light curve output file should be overwritten.
  
 ``(debug = no) [boolean]``
     Enables debug mode. In debug mode the executable will dump any log file output to the console.
@@ -40,6 +40,7 @@ cscripts
    :maxdepth: 1
  
    cscaldb --- Lists available instrument response functions <cscaldb>
+   csebins --- Generates energy boundaries for stacked analysis <csebins>
    cslightcrv --- Computes lightcurve <cslightcrv>
    csmodelinfo --- Shows model container content <csmodelinfo>
    csmodelmerge --- Merges several model containers into one file <csmodelmerge>
+#! /usr/bin/env python
+# ==========================================================================
+# This scripts performs unit tests for the cslightcrv script.
+#
+# Copyright (C) 2016-2017 Juergen Knoedlseder
+#
+# 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 <http://www.gnu.org/licenses/>.
+#
+# ==========================================================================
+import gammalib
+import cscripts
+from testing import test
+
+
+# ============================= #
+# Test class for csebins script #
+# ============================= #
+class Test(test):
+    """
+    Test class for csebins script
+
+    This test class makes unit tests for the csebins script by using it from
+    the command line and from Python.
+    """
+    
+    # Constructor
+    def __init__(self):
+        """
+        Constructor
+        """
+        # Call base class constructor
+        test.__init__(self)
+
+        # Return
+        return
+
+    # Set test functions
+    def set(self):
+        """
+        Set all test functions.
+        """
+        # Set test name
+        self.name('csebins')
+
+        # Append tests
+        self.append(self._test_cmd, 'Test csebins on command line')
+        self.append(self._test_python, 'Test csebins from Python')
+
+        # Return
+        return
+
+    # Test csebins on command line
+    def _test_cmd(self):
+        """
+        Test csebins on the command line.
+        """
+        # Set script name
+        cslightcrv = self._script('csebins')
+
+        # Setup cslightcrv command
+        cmd = cslightcrv+' inobs="NONE"'+ \
+                         ' caldb="'+self._caldb+'" irf="'+self._irf+'"'+ \
+                         ' emin=0.1 emax=100.0'+ \
+                         ' outfile="csebins_cmd1.fits"'+ \
+                         ' logfile="csebins_cmd1.log" chatter=1'
+
+        # Check if execution of wrong command fails
+        self.test_assert(self._execute('command_that_does_not_exist') != 0,
+             'Self test of test script')
+
+        # Check if execution was successful
+        self.test_assert(self._execute(cmd) == 0,
+             'Check successful execution from command line')
+
+        # Check energy boundaries
+        self._check_ebounds('csebins_cmd1.fits', 33)
+
+        # Setup cslightcrv command
+        cmd = cslightcrv+' inobs="NONE"'+ \
+                         ' caldb="Database_that_does_not_exist" irf="'+self._irf+'"'+ \
+                         ' emin=0.1 emax=100.0'+ \
+                         ' outfile="csebins_cmd2.fits"'+ \
+                         ' logfile="csebins_cmd2.log" chatter=1'
+
+        # Check if execution failed
+        self.test_assert(self._execute(cmd) != 0,
+             'Check invalid input file when executed from command line')
+
+        # Return
+        return
+
+    # Test csebins from Python
+    def _test_python(self):
+        """
+        Test csebins from Python
+        """
+        # Set-up csebins
+        ebins = cscripts.csebins()
+        ebins['inobs']   = 'NONE'
+        ebins['caldb']   = self._caldb
+        ebins['irf']     = self._irf
+        ebins['emin']    = 0.1
+        ebins['emax']    = 100.0
+        ebins['outfile'] = 'csebins_py1.fits'
+        ebins['logfile'] = 'csebins_py1.log'
+        ebins['chatter'] = 2
+
+        # Run csebins script and energy boundaries
+        ebins.logFileOpen()   # Make sure we get a log file
+        ebins.run()
+        ebins.save()
+
+        # Check energy boundaries
+        self._check_ebounds('csebins_py1.fits', 33)
+
+        # Return
+        return
+
+    # Check energy boundaries result file
+    def _check_ebounds(self, filename, bins):
+        """
+        Check energy boundaries file
+        """
+        # Load energy boundaries file
+        ebounds = gammalib.GEbounds(filename)
+
+        # Check number of energy boundaries
+        self.test_value(ebounds.size(), bins)
+        
+        # Return
+        return
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This scripts performs unit tests for cscripts
 #
-# Copyright (C) 2016 Juergen Knoedlseder
+# Copyright (C) 2016-2017 Juergen Knoedlseder
 #
 # 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
@@ -25,6 +25,7 @@ import cscripts
 sys.path.append(os.path.dirname(os.path.abspath(__file__)))
 import test_cscript
 import test_cscaldb
+import test_csebins
 import test_csfindobs
 import test_csinfo
 import test_cslightcrv
@@ -113,6 +114,7 @@ def test(installed=False, debug=False):
     # Define list of test suites
     tests = [test_cscript.Test(),
              test_cscaldb.Test(),
+             test_csebins.Test(),
              test_csinfo.Test(),
              test_cslightcrv.Test(),
              test_csmodelinfo.Test(),