Add caldb and irf parameters to csbkgmodel, add min/max limits to "Sigma" model …

…parameter, and add csbkgmodel unit test

Add caldb and irf parameters to csbkgmodel, add min/max limits to "Sigma" model …
…parameter, and add csbkgmodel unit test
Jürgen Knödlseder
1 parent c0406cf2
Showing 6 changed files with 329 additions and 14 deletions
cscripts/Makefile.am
cscripts/csbkgmodel.par
cscripts/csbkgmodel.py
doc/source/users/reference_manual/csbkgmodel.rst
test/test_csbkgmodel.py
test/test_python_cscripts.py
@@ -129,6 +129,7 @@ testscripts = $(srcdir)/tests/__init__.py \
               $(top_srcdir)/test/test_python_cscripts.py \
               $(top_srcdir)/test/test_cscript.py \
               $(top_srcdir)/test/test_cscaldb.py \
+              $(top_srcdir)/test/test_csbkgmodel.py \
               $(top_srcdir)/test/test_csebins.py \
               $(top_srcdir)/test/test_csfindobs.py \
               $(top_srcdir)/test/test_csinfo.py \
@@ -17,6 +17,8 @@
 # Input and output parameters
 #============================
 inobs,    f, a, obs.xml,,, "Input observation definition XML file"
+caldb,    s, a, prod2,,, "Calibration database"
+irf,      s, a, South_0.5h,,, "Instrument response function"
 outmodel, f, a, model.xml,,, "Output model definition XML file"
  
 #
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Generates background model
 #
-# Copyright (C) 2018 Juergen Knoedlseder
+# Copyright (C) 2018-2019 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
@@ -53,9 +53,9 @@ class csbkgmodel(ctools.csobservation):
         """
         Get parameters from parfile
         """
-        # Setup observations (do not require response, allow event list and
-        # counts cube)
-        self._setup_observations(self.obs(), False, True, True)
+        # Setup observations (require response, allow event list and counts
+        # cube)
+        self._setup_observations(self.obs(), True, True, True)
  
         # Set instrument name
         self._instrument = self._get_instrument()
@@ -135,11 +135,15 @@ class csbkgmodel(ctools.csobservation):
             if self['gradient'].boolean():
                 spatial = gammalib.GCTAModelSpatialMultiplicative()
                 factor1 = gammalib.GCTAModelRadialGauss(3.0)
+                factor1['Sigma'].min(1.5)
+                factor1['Sigma'].max(10.0)
                 factor2 = gammalib.GCTAModelSpatialGradient()
                 spatial.append(factor1)
                 spatial.append(factor2)
             else:
                 spatial = gammalib.GCTAModelRadialGauss(3.0)
+                spatial['Sigma'].min(1.5)
+                spatial['Sigma'].max(10.0)
             epivot   = gammalib.GEnergy(1.0, 'TeV')
             spectral = gammalib.GModelSpectralPlaw(3.0e-4, -1.5, epivot)
             spectral['Prefactor'].min(1.0e-8)
@@ -303,15 +307,6 @@ class csbkgmodel(ctools.csobservation):
     def _generate_runwise_bkg(self):
         """
         Generate background models
-
-        Parameters
-        ----------
-        responses : list of dict
-            List of response dictionaries
-        energy : float
-            Energy to insert (TeV)
-        comment : str
-            Reason for energy insertion
         """
         # Loop over observations
         for run in self.obs():
@@ -20,6 +20,12 @@ General parameters
 ``inobs [file]``
     Event list, counts cube or input observation definition XML file.
  
+``caldb [string]``
+    Calibration database.
+
+``irf [string]``
+    Instrument response function.
+
 ``outmodel [file]``
     Output model XML file.
  
+#! /usr/bin/env python
+# ==========================================================================
+# This scripts performs unit tests for the csbkgmodel script.
+#
+# Copyright (C) 2019 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 csbkgmodel script #
+# ================================ #
+class Test(test):
+    """
+    Test class for csbkgmodel script
+
+    This test class makes unit tests for the csbkgmodel script by using it
+    from the command line and from Python.
+    """
+
+    # Constructor
+    def __init__(self):
+        """
+        Constructor
+        """
+        # Call base class constructor
+        test.__init__(self)
+
+        # Set test datasets and parameters
+        self._myevents1 = self._datadir + '/crab_offaxis1.fits'
+        self._myevents2 = self._datadir + '/crab_offaxis2.fits'
+
+        # Return
+        return
+
+    # Set test functions
+    def set(self):
+        """
+        Set all test functions.
+        """
+        # Set test name
+        self.name('csbkgmodel')
+
+        # Append tests
+        self.append(self._test_cmd, 'Test csbkgmodel on command line')
+        self.append(self._test_python, 'Test csbkgmodel from Python')
+        self.append(self._test_pickeling, 'Test csbkgmodel pickeling')
+
+        # Return
+        return
+
+    # Test csbkgmodel on command line
+    def _test_cmd(self):
+        """
+        Test csbkgmodel on the command line.
+        """
+        # Set script name
+        csbkgmodel = self._script('csbkgmodel')
+
+        # Setup csbkgmodel command
+        cmd = csbkgmodel+' inobs="'+self._events+'"'+ \
+                         ' caldb="'+self._caldb+'" irf="'+self._irf+'"' \
+                         ' instrument=CTA spatial=GAUSS gradient=yes'+ \
+                         ' spectral=NODES ebinalg=LOG emin=1.0 emax=100.0'+ \
+                         ' enumbins=8 runwise=yes rad=2.0'+ \
+                         ' outmodel="csbkgmodel_cmd1.xml"'+ \
+                         ' logfile="csbkgmodel_cmd1.log" chatter=4'
+
+        # Check if execution was successful
+        self.test_assert(self._execute(cmd) == 0,
+             'Check successful execution from command line')
+
+        # Check background model
+        self._check_bkg_model('csbkgmodel_cmd1.xml')
+
+        # Setup csbkgmodel command
+        cmd = csbkgmodel+' inobs="events_that_do_not_exist.fits"'+ \
+                         ' caldb="'+self._caldb+'" irf="'+self._irf+'"' \
+                         ' instrument=CTA spatial=GAUSS gradient=yes'+ \
+                         ' spectral=NODES ebinalg=LOG emin=1.0 emax=100.0'+ \
+                         ' enumbins=8 runwise=yes rad=2.0'+ \
+                         ' outmodel="csbkgmodel_cmd2.xml"' + \
+                         ' logfile="csbkgmodel_cmd2.log" chatter=1'
+
+        # Check if execution failed
+        self.test_assert(self._execute(cmd, success=False) != 0,
+             'Check invalid input file when executed from command line')
+
+        # Check csbkgmodel --help
+        self._check_help(csbkgmodel)
+
+        # Return
+        return
+
+    # Test csbkgmodel from Python
+    def _test_python(self):
+        """
+        Test csbkgmodel from Python
+        """
+        # Set-up csbkgmodel
+        bkgmodel = cscripts.csbkgmodel()
+        bkgmodel['inobs']      = self._events
+        bkgmodel['caldb']      = self._caldb
+        bkgmodel['irf']        = self._irf
+        bkgmodel['instrument'] = 'CTA'
+        bkgmodel['spatial']    = 'GAUSS'
+        bkgmodel['gradient']   = True
+        bkgmodel['spectral']   = 'NODES'
+        bkgmodel['ebinalg']    = 'LOG'
+        bkgmodel['emin']       = 1.0
+        bkgmodel['emax']       = 100.0
+        bkgmodel['enumbins']   = 8
+        bkgmodel['runwise']    = True
+        bkgmodel['rad']        = 2.0
+        bkgmodel['chatter']    = 2
+        bkgmodel['outmodel']   = 'csbkgmodel_py1.xml'
+        bkgmodel['logfile']    = 'csbkgmodel_py1.log'
+
+        # Run csbkgmodel script and save background model
+        bkgmodel.logFileOpen()   # Make sure we get a log file
+        bkgmodel.run()
+        bkgmodel.save()
+
+        # Check background model
+        self._check_bkg_model('csbkgmodel_py1.xml')
+
+        # Now test without gradient, power law and not runwise
+        bkgmodel = cscripts.csbkgmodel()
+        bkgmodel['inobs']      = self._events
+        bkgmodel['caldb']      = self._caldb
+        bkgmodel['irf']        = self._irf
+        bkgmodel['instrument'] = 'CTA'
+        bkgmodel['spatial']    = 'GAUSS'
+        bkgmodel['gradient']   = False
+        bkgmodel['spectral']   = 'PLAW'
+        bkgmodel['ebinalg']    = 'LOG'
+        bkgmodel['emin']       = 1.0
+        bkgmodel['emax']       = 100.0
+        bkgmodel['enumbins']   = 8
+        bkgmodel['runwise']    = False
+        bkgmodel['rad']        = 2.0
+        bkgmodel['chatter']    = 3
+        bkgmodel['outmodel']   = 'csbkgmodel_py2.xml'
+        bkgmodel['logfile']    = 'csbkgmodel_py2.log'
+
+        # Execute csbkgmodel script
+        bkgmodel.logFileOpen()   # Make sure we get a log file
+        bkgmodel.execute()
+
+        # Check background model
+        self._check_bkg_model('csbkgmodel_py2.xml')
+
+        # Now test AEFF model
+        bkgmodel = cscripts.csbkgmodel()
+        bkgmodel['inobs']      = self._events
+        bkgmodel['caldb']      = self._caldb
+        bkgmodel['irf']        = self._irf
+        bkgmodel['instrument'] = 'CTA'
+        bkgmodel['spatial']    = 'AEFF'
+        bkgmodel['gradient']   = False
+        bkgmodel['spectral']   = 'PLAW'
+        bkgmodel['ebinalg']    = 'LOG'
+        bkgmodel['emin']       = 1.0
+        bkgmodel['emax']       = 100.0
+        bkgmodel['enumbins']   = 8
+        bkgmodel['runwise']    = False
+        bkgmodel['rad']        = 2.0
+        bkgmodel['chatter']    = 4
+        bkgmodel['outmodel']   = 'csbkgmodel_py3.xml'
+        bkgmodel['logfile']    = 'csbkgmodel_py3.log'
+
+        # Execute csbkgmodel script
+        bkgmodel.logFileOpen()   # Make sure we get a log file
+        bkgmodel.execute()
+
+        # Check background model
+        self._check_bkg_model('csbkgmodel_py3.xml')
+
+        # Now test IRF model
+        bkgmodel = cscripts.csbkgmodel()
+        bkgmodel['inobs']      = self._events
+        bkgmodel['caldb']      = self._caldb
+        bkgmodel['irf']        = self._irf
+        bkgmodel['instrument'] = 'CTA'
+        bkgmodel['spatial']    = 'AEFF'
+        bkgmodel['gradient']   = False
+        bkgmodel['spectral']   = 'PLAW'
+        bkgmodel['ebinalg']    = 'LOG'
+        bkgmodel['emin']       = 1.0
+        bkgmodel['emax']       = 100.0
+        bkgmodel['enumbins']   = 8
+        bkgmodel['runwise']    = False
+        bkgmodel['rad']        = 2.0
+        bkgmodel['chatter']    = 4
+        bkgmodel['outmodel']   = 'csbkgmodel_py4.xml'
+        bkgmodel['logfile']    = 'csbkgmodel_py4.log'
+
+        # Execute csbkgmodel script
+        bkgmodel.logFileOpen()   # Make sure we get a log file
+        bkgmodel.execute()
+
+        # Check background model
+        self._check_bkg_model('csbkgmodel_py4.xml')
+
+        # Test with multiple input observations
+        obs = gammalib.GObservations()
+        for s, events in enumerate([self._myevents1, self._myevents2]):
+            run = gammalib.GCTAObservation(events)
+            run.id(str(s + 1))
+            run.response(self._irf, gammalib.GCaldb('cta', self._caldb))
+            obs.append(run)
+
+        # Set-up csbkgmodel
+        bkgmodel = cscripts.csbkgmodel(obs)
+        bkgmodel['instrument'] = 'CTA'
+        bkgmodel['spatial']    = 'GAUSS'
+        bkgmodel['gradient']   = True
+        bkgmodel['spectral']   = 'NODES'
+        bkgmodel['ebinalg']    = 'LOG'
+        bkgmodel['emin']       = 1.0
+        bkgmodel['emax']       = 100.0
+        bkgmodel['enumbins']   = 8
+        bkgmodel['runwise']    = True
+        bkgmodel['rad']        = 2.0
+        bkgmodel['chatter']    = 2
+        bkgmodel['outmodel']   = 'csbkgmodel_py5.xml'
+        bkgmodel['logfile']    = 'csbkgmodel_py5.log'
+
+        # Execute csbkgmodel script
+        bkgmodel.logFileOpen()   # Make sure we get a log file
+        bkgmodel.execute()
+
+        # Check background model
+        self._check_bkg_model('csbkgmodel_py5.xml', nmodels=2)
+
+        # Return
+        return
+
+    # Test csbkgmodel pickeling
+    def _test_pickeling(self):
+        """
+        Test csbkgmodel pickeling
+        """
+        # Perform pickeling tests of empty class
+        self._pickeling(cscripts.csbkgmodel())
+
+        # Set-up csbkgmodel
+        bkgmodel = cscripts.csbkgmodel()
+        bkgmodel['inobs']      = self._events
+        bkgmodel['caldb']      = self._caldb
+        bkgmodel['irf']        = self._irf
+        bkgmodel['instrument'] = 'CTA'
+        bkgmodel['spatial']    = 'GAUSS'
+        bkgmodel['gradient']   = True
+        bkgmodel['spectral']   = 'NODES'
+        bkgmodel['ebinalg']    = 'LOG'
+        bkgmodel['emin']       = 1.0
+        bkgmodel['emax']       = 100.0
+        bkgmodel['enumbins']   = 8
+        bkgmodel['runwise']    = True
+        bkgmodel['rad']        = 2.0
+        bkgmodel['chatter']    = 2
+        bkgmodel['outmodel']   = 'csbkgmodel_py1_pickle.xml'
+        bkgmodel['logfile']    = 'csbkgmodel_py1_pickle.log'
+
+        # Perform pickeling tests of filled class
+        obj = self._pickeling(bkgmodel)
+
+        # Run csbkgmodel script and save background model
+        obj.logFileOpen()   # Make sure we get a log file
+        obj.run()
+        obj.save()
+
+        # Check background model
+        self._check_bkg_model('csbkgmodel_py1_pickle.xml')
+
+        # Return
+        return
+
+    # Check background model
+    def _check_bkg_model(self, filename, nmodels=1):
+        """
+        Check background model
+        """
+        # Load model
+        models = gammalib.GModels(filename)
+
+        # Check model size
+        self.test_value(models.size(), nmodels,
+                        'Check for %d model(s) in XML file' % nmodels)
+
+        # Return
+        return
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This scripts performs unit tests for cscripts
 #
-# Copyright (C) 2016-2018 Juergen Knoedlseder
+# Copyright (C) 2016-2019 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_csbkgmodel
 import test_csebins
 import test_csfindobs
 import test_csinfo
@@ -119,6 +120,7 @@ def test(installed=False, debug=False):
  
     # Define list of test suites
     tests = [test_cscript.Test(),
+             test_csbkgmodel.Test(),
              test_cscaldb.Test(),
              test_csebins.Test(),
              test_csinfo.Test(),