Implement cslightcrv unit test.

Jürgen Knödlseder
1 parent d0045ac8
Showing 6 changed files with 349 additions and 52 deletions
cscripts/Makefile.am
cscripts/cslightcrv.py
test/test_cscaldb.py
test/test_cslightcrv.py
test/test_ctbin.py
test/test_python_cscripts.py
@@ -26,13 +26,17 @@ linkscripts = cspull cssens cstsdist cscaldb csresmap \
 # Test scripts to be distributed in $(pyexecdir)/cscripts/tests
 testscripts = tests/__init__.py \
               ../test/test_python_cscripts.py \
-              ../test/test_cscaldb.py
+              ../test/test_cscaldb.py \
+              ../test/test_cslightcrv.py
 # Test data to be distributed in $(pyexecdir)/cscripts/tests/data
-testdata    = 
+testdata    = ../test/data/crab.xml \
+              ../test/data/crab_events.fits.gz \
+              ../test/data/crab_cntmap.fits.gz \
+              ../test/data/exclusion.reg
 # Test IRFs to be distributed in $(pyexecdir)/cscripts/tests/irf
-testirf     = 
+#testirf     = 
 # Distribute scripts and parfiles
 EXTRA_DIST         = $(binscripts) $(pythonscripts)
@@ -46,8 +50,8 @@ testscriptsdir     = $(pyexecdir)/cscripts/tests
 testscripts_DATA   = $(testscripts)
 testdatadir        = $(pyexecdir)/cscripts/tests/data
 testdata_DATA      = $(testdata)
-testirfdir         = $(pyexecdir)/cscripts/tests/irf
-testirf_DATA       = $(testirf)
+#testirfdir         = $(pyexecdir)/cscripts/tests/irf
+#testirf_DATA       = $(testirf)
 # Add hook rule to create symbolic links for cscripts
 install-data-hook:
@@ -76,12 +76,15 @@ class cslightcrv(ctools.cscript):
         self._stacked = False
         # Initialise observation container from constructor arguments.
-        self._obs = self._set_input_obs(argv)
+        self._obs, argv = self._set_input_obs(argv)
-        # Initialise application by calling the appropriate class
-        # constructor.
+        # Initialise script by calling the appropriate class constructor.
         self._init_cscript(argv)
+        # Set logger properties
+        self._log_header()
+        self._log.date(True)
+
         # Return
         return
@@ -91,12 +94,10 @@ class cslightcrv(ctools.cscript):
         """
         Get parameters from parfile.
         """
-        # Set observation if not done before
-        if self._obs == None or self._obs.size() == 0:
-            self._require_inobs("cslightcrv::_get_parameters()")
-            self._obs = self._get_observations()
+        # Setup observations
+        self._setup_observations(self._obs)
-        # Set models if we have none
+        # Set models if there are none in the container
         if self._obs.models().size() == 0:
             self._obs.models(self["inmodel"].filename())
@@ -436,16 +437,26 @@ class cslightcrv(ctools.cscript):
         if self._logExplicit():
             self._log.header3("Creating point spread function cube")
-        # Create psf cube
+        # Compute spatial binning for point spread function and
+        # energy dispersion cubes
+        binsz = 0.1 * self["binsz"].real()
+        nxpix = self["nxpix"].integer() / 10
+        nypix = self["nypix"].integer() / 10
+        if nxpix < 2:
+            nxpix = 2
+        if nypix < 2:
+            nypix = 2
+
+        # Create point spread function cube
         psfcube = ctools.ctpsfcube(obs)
         psfcube["incube"]   = "NONE"
         psfcube["usepnt"]   = False
         psfcube["ebinalg"]  = "LOG"
         psfcube["xref"]     = self["xref"].real()
         psfcube["yref"]     = self["yref"].real()
-        psfcube["binsz"]    = self["binsz"].real()
-        psfcube["nxpix"]    = self["nxpix"].integer()
-        psfcube["nypix"]    = self["nypix"].integer()
+        psfcube["binsz"]    = binsz
+        psfcube["nxpix"]    = nxpix
+        psfcube["nypix"]    = nypix
         psfcube["enumbins"] = self["enumbins"].integer()
         psfcube["emin"]     = self["emin"].real()
         psfcube["emax"]     = self["emax"].real()    
@@ -460,16 +471,16 @@ class cslightcrv(ctools.cscript):
             if self._logExplicit():
                 self._log.header3("Creating energy dispersion cube")
-            # Create edisp cube
+            # Create energy dispersion cube
             edispcube = ctools.ctedispcube(obs)
             edispcube["incube"]   = "NONE"
             edispcube["usepnt"]   = False
             edispcube["ebinalg"]  = "LOG"
             edispcube["xref"]     = self["xref"].real()
             edispcube["yref"]     = self["yref"].real()
-            edispcube["binsz"]    = self["binsz"].real()
-            edispcube["nxpix"]    = self["nxpix"].integer()
-            edispcube["nypix"]    = self["nypix"].integer()
+            edispcube["binsz"]    = binsz
+            edispcube["nxpix"]    = nxpix
+            edispcube["nypix"]    = nypix
             edispcube["enumbins"] = self["enumbins"].integer()
             edispcube["emin"]     = self["emin"].real()
             edispcube["emax"]     = self["emax"].real()    
@@ -18,6 +18,7 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 # ==========================================================================
+import os
 import gammalib
 import cscripts
@@ -28,7 +29,11 @@ import cscripts
 class Test(gammalib.GPythonTestSuite):
     """
     Test class for cscaldb script.
+
+    This test class makes unit tests for the cslightcrv script by using it
+    from the command line and from Python.
     """
+
     # Constructor
     def __init__(self):
         """
@@ -49,15 +54,15 @@ class Test(gammalib.GPythonTestSuite):
         self.name("cscaldb")
         # Append tests
-        self.append(self._test_functional, "Test cscaldb functionality")
+        self.append(self._test_python, "Test cscaldb from Python")
         # Return
         return
-    # Test cscaldb functionality
-    def _test_functional(self):
+    # Test cscaldb from Python
+    def _test_python(self):
         """
-        Test cscaldb functionality.
+        Test cscaldb from Python.
         """
         # Set-up cscaldb
         caldb = cscripts.cscaldb()
+#! /usr/bin/env python
+# ==========================================================================
+# This scripts performs unit tests for the cslightcrv script.
+#
+# Copyright (C) 2016 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 os
+import gammalib
+import cscripts
+
+
+# ================================ #
+# Test class for cslightcrv script #
+# ================================ #
+class Test(gammalib.GPythonTestSuite):
+    """
+    Test class for cslightcrv script.
+
+    This test class makes unit tests for the cslightcrv script by using it
+    from the command line and from Python.
+    """
+    
+    # Constructor
+    def __init__(self):
+        """
+        Constructor.
+        """
+        # Call base class constructor
+        gammalib.GPythonTestSuite.__init__(self)
+
+        # Set members
+        self._events_name = "data/crab_events.fits"
+        self._model_name  = "data/crab.xml"
+
+        # Return
+        return
+
+    # Set test functions
+    def set(self):
+        """
+        Set all test functions.
+        """
+        # Set test name
+        self.name("cslightcrv")
+
+        # Append tests
+        self.append(self._test_cmd, "Test cslightcrv on command line")
+        self.append(self._test_python, "Test cslightcrv from Python")
+
+        # Return
+        return
+
+    # Test cslightcrv on command line
+    def _test_cmd(self):
+        """
+        Test cslightcrv on the command line.
+        """
+        # Kluge to set the command (installed version has no README file)
+        if os.path.isfile("README"):
+            cslightcrv = "../cscripts/cslightcrv.py"
+        else:
+            cslightcrv = "cslightcrv"
+
+        # Setup cslightcrv command
+        cmd = cslightcrv+' inobs="'+self._events_name+'"'+ \
+                         ' inmodel="'+self._model_name+'"'+ \
+                         ' srcname="Crab" caldb="prod2" irf="South_0.5h"'+ \
+                         ' tbinalg="LIN" tmin=51544.50 tmax=51544.53'+ \
+                         ' tbins=3 enumbins=0 emin=0.1 emax=100.0'+ \
+                         ' outfile="lightcurve_cmd1.fits"'
+
+        # Execute cslightcrv, make sure we catch any exception
+        try:
+            rc = os.system(cmd+" &>/dev/null")
+        except:
+            pass
+
+        # Check if execution was successful
+        self.test_assert(rc == 0,
+                         "Successful cslightcrv execution on command line")
+
+        # Check light curve
+        self._check_light_curve("lightcurve_cmd1.fits", 3)
+
+        # Setup cslightcrv command
+        cmd = cslightcrv+' inobs="events_that_do_not_exist.fits"'+ \
+                         ' inmodel="'+self._model_name+'"'+ \
+                         ' srcname="Crab" caldb="prod2" irf="South_0.5h"'+ \
+                         ' tbinalg="LIN" tmin=51544.50 tmax=51544.53'+ \
+                         ' tbins=3 enumbins=0 emin=0.1 emax=100.0'+ \
+                         ' outfile="lightcurve_cmd1.fits"'
+
+        # Execute cslightcrv, make sure we catch any exception
+        try:
+            rc = os.system(cmd+" &>/dev/null")
+        except:
+            pass
+
+        # Check if execution failed
+        self.test_assert(rc != 0,
+                         "Failure of cslightcrv execution on command line")
+
+        # Return
+        return
+
+    # Test cslightcrv from Python
+    def _test_python(self):
+        """
+        Test cslightcrv from Python.
+        """
+        # Set-up unbinned cslightcrv
+        lcrv = cscripts.cslightcrv()
+        lcrv["inobs"]    = self._events_name
+        lcrv["inmodel"]  = self._model_name
+        lcrv["srcname"]  = "Crab"
+        lcrv["caldb"]    = "prod2"
+        lcrv["irf"]      = "South_0.5h"
+        lcrv["tbinalg"]  = "LIN"
+        lcrv["tmin"]     = 51544.50
+        lcrv["tmax"]     = 51544.53
+        lcrv["tbins"]    = 3
+        lcrv["enumbins"] = 0
+        lcrv["emin"]     = 0.1
+        lcrv["emax"]     = 100.0
+        lcrv["outfile"]  = "lightcurve_py1.fits"
+
+        # Run cslightcrv script and save light curve
+        lcrv.run()
+        lcrv.save()
+
+        # Check light curve
+        self._check_light_curve("lightcurve_py1.fits", 3)
+
+        # Now use FILE as time bin algorithm. For this we need first to
+        # create an ASCII file. We use now 6 time bins. The ASCII file
+        # is saved into the file "lightcurve_py2.dat".
+        csv    = gammalib.GCsv(2,2)
+        tmin   = 51544.50
+        tdelta = 0.01
+        for i in range(csv.nrows()):
+            csv[i,0] = "%.5f" % (tmin +  i    * tdelta)
+            csv[i,1] = "%.5f" % (tmin + (i+1) * tdelta)
+        csv.save("lightcurve_py2.dat", " ", True)
+
+        # Set-up unbinned cslightcrv
+        lcrv = cscripts.cslightcrv()
+        lcrv["inobs"]    = self._events_name
+        lcrv["inmodel"]  = self._model_name
+        lcrv["srcname"]  = "Crab"
+        lcrv["caldb"]    = "prod2"
+        lcrv["irf"]      = "South_0.5h"
+        lcrv["tbinalg"]  = "FILE"
+        lcrv["tbinfile"] = "lightcurve_py2.dat"
+        lcrv["enumbins"] = 0
+        lcrv["emin"]     = 0.1
+        lcrv["emax"]     = 100.0
+        lcrv["outfile"]  = "lightcurve_py2.fits"
+
+        # Execute cslightcrv script
+        lcrv.execute()
+
+        # Check light curve
+        self._check_light_curve("lightcurve_py2.fits", 2)
+
+        # Now we setup an observation container on input. We attached the
+        # model to the observation container so that cslightcrv should
+        # no longer query for the parameter.
+        cta = gammalib.GCTAObservation(self._events_name)
+        obs = gammalib.GObservations()
+        obs.append(cta)
+        obs.models(self._model_name)
+
+        # Set-up unbinned cslightcrv from observation container. Now use
+        # the GTI algorithm so that we test all timing algorithms.
+        lcrv = cscripts.cslightcrv(obs)
+        lcrv["srcname"]  = "Crab"
+        lcrv["caldb"]    = "prod2"
+        lcrv["irf"]      = "South_0.5h"
+        lcrv["tbinalg"]  = "GTI"
+        lcrv["enumbins"] = 0
+        lcrv["emin"]     = 0.1
+        lcrv["emax"]     = 100.0
+        lcrv["outfile"]  = "lightcurve_py3.fits"
+
+        # Execute cslightcrv script
+        lcrv.execute()
+
+        # Check light curve
+        self._check_light_curve("lightcurve_py3.fits", 1)
+
+        # Finally we set-up a binned cslightcrv
+        lcrv = cscripts.cslightcrv()
+        lcrv["inobs"]    = self._events_name
+        lcrv["inmodel"]  = self._model_name
+        lcrv["srcname"]  = "Crab"
+        lcrv["caldb"]    = "prod2"
+        lcrv["irf"]      = "South_0.5h"
+        lcrv["tbinalg"]  = "LIN"
+        lcrv["tmin"]     = 51544.50
+        lcrv["tmax"]     = 51544.53
+        lcrv["tbins"]    = 2
+        lcrv["emin"]     = 0.1
+        lcrv["emax"]     = 100.0
+        lcrv["enumbins"] = 10
+        lcrv["coordsys"] = "CEL"
+        lcrv["proj"]     = "TAN"
+        lcrv["xref"]     = 83.63
+        lcrv["yref"]     = 22.01
+        lcrv["nxpix"]    = 20
+        lcrv["nypix"]    = 20
+        lcrv["binsz"]    = 0.02
+        lcrv["outfile"]  = "lightcurve_py4.fits"
+
+        # Execute cslightcrv script
+        lcrv.execute()
+
+        # Check light curve
+        self._check_light_curve("lightcurve_py4.fits", 2)
+
+        # Return
+        return
+
+    # Check light curve result file
+    def _check_light_curve(self, filename, bins):
+        """
+        Check light curve file.
+        """
+        # Expected column names
+        cols = ["MJD", "e_MJD", "Prefactor", "e_Prefactor",
+                "Index", "e_Index", "TS", "UpperLimit"]
+
+        # Open FITS file
+        fits = gammalib.GFits(filename)
+
+        # Check FITS file structure
+        self.test_value(fits.size(), 2,
+             'Check for 2 extensions in light curve FITS file')
+        self.test_assert(fits.contains("LIGHTCURVE"),
+             'FITS file contains "LIGHTCURVE" extension')
+
+        # Get LIGHTCURVE table
+        table = fits["LIGHTCURVE"]
+
+        # Check FITS table structure
+        self.test_value(table.ncols(), len(cols),
+             'Check for %d columns in light curve FITS table' % len(cols))
+        self.test_value(table.nrows(), bins,
+             'Check for %d rows in light curve FITS table' % bins)
+        for col in cols:
+            self.test_assert(table.contains(col),
+                 'FITS file contains "'+col+'" column')
+
+        # Close FITS file
+        fits.close()
+        
+        # Return
+        return
@@ -57,14 +57,14 @@ class Test(gammalib.GPythonTestSuite):
         self.name("ctbin")
         # Append tests
-        self.append(self._test_ctbin_cmd, "Test ctbin on command line")
-        self.append(self._test_ctbin_python, "Test ctbin from Python")
+        self.append(self._test_cmd, "Test ctbin on command line")
+        self.append(self._test_python, "Test ctbin from Python")
         # Return
         return
     # Test ctbin on command line
-    def _test_ctbin_cmd(self):
+    def _test_cmd(self):
         """
         Test ctbin on the command line.
         """
@@ -93,7 +93,7 @@ class Test(gammalib.GPythonTestSuite):
         # Load counts cube and check content.
         evt = gammalib.GCTAEventCube("cntmap_cmd1.fits")
-        self._test_cube(evt, 5542)
+        self._check_cube(evt, 5542)
         # Setup ctbin command
         cmd = ctbin+' inobs="events_that_do_not_exist.fits"'+ \
@@ -116,7 +116,7 @@ class Test(gammalib.GPythonTestSuite):
         return
     # Test ctbin from Python
-    def _test_ctbin_python(self):
+    def _test_python(self):
         """
         Test ctbin from Python.
         """
@@ -140,15 +140,15 @@ class Test(gammalib.GPythonTestSuite):
         bin.run()
         # Check content of observation and cube
-        self._test_observation(bin, 5542)
-        self._test_cube(bin.cube(), 5542)
+        self._check_observation(bin, 5542)
+        self._check_cube(bin.cube(), 5542)
         # Test copy constructor
         cpy_bin = bin.copy()
         # Check content of observation and cube
-        self._test_observation(cpy_bin, 5542)
-        self._test_cube(cpy_bin.cube(), 5542)
+        self._check_observation(cpy_bin, 5542)
+        self._check_cube(cpy_bin.cube(), 5542)
         # Run copy of ctbin tool again
         cpy_bin.run()
@@ -157,8 +157,8 @@ class Test(gammalib.GPythonTestSuite):
         # event cube as on input the observation is binned, and any binned
         # observation will be skipped, hence the counts cube should be
         # empty.
-        self._test_observation(cpy_bin, 0)
-        self._test_cube(cpy_bin.cube(), 0)
+        self._check_observation(cpy_bin, 0)
+        self._check_cube(cpy_bin.cube(), 0)
         # Clear and run copy of ctbin again
         #cpy_bin.clear()
@@ -172,7 +172,7 @@ class Test(gammalib.GPythonTestSuite):
         # Load counts cube and check content.
         evt = gammalib.GCTAEventCube("cntmap.fits")
-        self._test_cube(evt, 5542)
+        self._check_cube(evt, 5542)
         # Prepare observation container for stacked analysis
         cta = gammalib.GCTAObservation(self._events_name)
@@ -204,8 +204,8 @@ class Test(gammalib.GPythonTestSuite):
         # Check content of observation and cube (need multiplier=3 since
         # three identical observations have been appended)
-        self._test_observation(bin, 5542, multiplier=3)
-        self._test_cube(bin.cube(), 5542, multiplier=3)
+        self._check_observation(bin, 5542, multiplier=3)
+        self._check_cube(bin.cube(), 5542, multiplier=3)
         # Set-up ctbin using an observation container
         bin = ctools.ctbin(obs)
@@ -227,7 +227,7 @@ class Test(gammalib.GPythonTestSuite):
         # Load counts cube and check content.
         evt = gammalib.GCTAEventCube("cntmap2.fits")
-        self._test_cube(evt, 5542, multiplier=3)
+        self._check_cube(evt, 5542, multiplier=3)
         # Clear tool
         #bin.clear()
@@ -242,9 +242,9 @@ class Test(gammalib.GPythonTestSuite):
         return
     # Check observation
-    def _test_observation(self, ctbin, nevents, multiplier=1):
+    def _check_observation(self, ctbin, nevents, multiplier=1):
         """
-        Test content of an observation.
+        Check content of an observation.
         Args:
             ctbin:   ctbin instance
@@ -270,15 +270,15 @@ class Test(gammalib.GPythonTestSuite):
                         "Pointing Declination is 22.01 deg")
         # Test event cube
-        self._test_cube(obs.events(), nevents, multiplier=multiplier)
+        self._check_cube(obs.events(), nevents, multiplier=multiplier)
         # Return
         return
     # Check event cube
-    def _test_cube(self, cube, nevents, multiplier=1):
+    def _check_cube(self, cube, nevents, multiplier=1):
         """
-        Test content of event cube."
+        Check content of event cube."
         Args:
             cube:    Event cube
@@ -18,12 +18,13 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 # ==========================================================================
+import os
+import sys
 import gammalib
 import cscripts
-import sys
-import os
 sys.path.append(os.path.dirname(os.path.abspath(__file__)))
 import test_cscaldb
+import test_cslightcrv
 # ================== #
@@ -37,13 +38,16 @@ def test(installed=False):
     suites = gammalib.GTestSuites("cscripts Python module unit testing")
     # Allocate test suites and append them to the container
-    suite_cscaldb = test_cscaldb.Test()
+    suite_cscaldb    = test_cscaldb.Test()
+    suite_cslightcrv = test_cslightcrv.Test()
     # Setup unit tests
     suite_cscaldb.set()
+    suite_cslightcrv.set()
     # Append tests to container
     suites.append(suite_cscaldb)
+    suites.append(suite_cslightcrv)
     # If we have an installed version then create a temporary directory and
     # copy over all information that is needed
@@ -60,9 +64,13 @@ def test(installed=False):
         dirname = os.path.dirname(testdir)
         # Copy over test data and irf
-        #os.system("cp -r %s %s" % (dirname+"/data", "data"))
+        os.system("cp -r %s %s" % (dirname+"/data", "data"))
         #os.system("cp -r %s %s" % (dirname+"/irf",  "irf"))
+    # ... otherwise set the calibration database
+    else:
+        os.environ['CALDB'] = "../caldb"
+
     # Set PFILES environment variable
     try:
         os.mkdir("pfiles")
@@ -97,8 +105,6 @@ def test(installed=False):
 # Main routine entry point #
 # ======================== #
 if __name__ == '__main__':
-    """
-    Perform unit testing for cscripts used from Python.
-    """
+
     # Run tests
     test()