Add H.E.S.S. background lookup to models and enhance csbkgmodel unit test

Jürgen Knödlseder
1 parent 19567f9c
Showing 6 changed files with 177 additions and 128 deletions
cscripts/Makefile.am
cscripts/csbkgmodel.py
models/Makefile.am
models/hess_bkg_lookup.fits
test/data/hess_bkg_lookup.fits
test/test_csbkgmodel.py
@@ -235,7 +235,8 @@ testdata = $(top_srcdir)/test/data/crab.xml \
            $(top_srcdir)/test/data/onoff_rmf2.fits.gz \
            $(top_srcdir)/test/data/model_crab_radialacceptance.xml \
            $(top_srcdir)/test/data/sigdist_map.fits.gz \
-           $(top_srcdir)/test/data/variability.fits.gz
+           $(top_srcdir)/test/data/variability.fits.gz \
+           $(top_srcdir)/test/data/hess_bkg_lookup.fits
  
 # Test data to be distributed in $(pyexecdir)/cscripts/tests/data/iactdata
 testiactdata = $(top_srcdir)/test/data/iactdata/hdu-index.fits.gz \
@@ -30,7 +30,6 @@ class csbkgmodel(ctools.csobservation):
     """
     Generates background model
     """
-
     # Constructor
     def __init__(self, *argv):
         """
@@ -46,7 +45,7 @@ class csbkgmodel(ctools.csobservation):
         # Return
         return
  
-    # State methods por pickling
+    # __getstate__ method for pickling
     def __getstate__(self):
         """
         Extend ctools.csobservation getstate method to include some members
@@ -64,6 +63,7 @@ class csbkgmodel(ctools.csobservation):
         # Return pickled dictionary
         return state
  
+    # __setstate__ method for pickling
     def __setstate__(self, state):
         """
         Extend ctools.csobservation setstate method to include some members
@@ -102,8 +102,9 @@ class csbkgmodel(ctools.csobservation):
                 self['smax'].real()
         if spatial == 'LOOKUP':
             self['slufile'].filename()
-        if spatial == 'LOOKUP' or spatial == 'GAUSS' or \
-           spatial == 'GAUSS(E)' or spatial == 'PROFILE':
+        if spatial == 'LOOKUP'   or spatial == 'GAUSS'   or \
+           spatial == 'GAUSS(E)' or spatial == 'PROFILE' or \
+           spatial == 'POLYNOM':
             self['gradient'].boolean()
         spectral = self['spectral'].string()
         if spectral == 'NODES':
@@ -178,133 +179,74 @@ class csbkgmodel(ctools.csobservation):
             spectral = gammalib.GModelSpectralPlaw(1.0e-13, -2.5, epivot)
             model    = gammalib.GCTAModelAeffBackground(spectral)
  
-        # Handle LOOKUP model
-        elif self['spatial'].string() == 'LOOKUP':
-
-            # Set spatial model
-            factor1 = gammalib.GCTAModelSpatialLookup(self['slufile'].filename())
-
-            # Optionally add gradient
-            if self['gradient'].boolean():
-                spatial = gammalib.GCTAModelSpatialMultiplicative()
-                factor2 = gammalib.GCTAModelSpatialGradient()
-                spatial.append(factor1)
-                spatial.append(factor2)
-            else:
-                spatial = factor1
-
-            # Set spectral model
-            epivot   = gammalib.GEnergy(1.0, 'TeV')
-            spectral = gammalib.GModelSpectralPlaw(3.0e-4, -1.5, epivot)
-            spectral['Prefactor'].min(1.0e-8)
-
-            # Set background model
-            model = gammalib.GCTAModelBackground(spatial, spectral)
-
-        # Handle GAUSS model
-        elif self['spatial'].string() == 'GAUSS':
-
-            # Set spatial model
-            factor1 = gammalib.GCTAModelRadialGauss(3.0)
-            factor1['Sigma'].min(sigma_min)
-            factor1['Sigma'].max(sigma_max)
-
-            # Optionally add gradient
-            if self['gradient'].boolean():
-                spatial = gammalib.GCTAModelSpatialMultiplicative()
-                factor2 = gammalib.GCTAModelSpatialGradient()
-                spatial.append(factor1)
-                spatial.append(factor2)
-            else:
-                spatial = factor1
+        # Handle other models
+        else:
  
-            # Set spectral model
-            epivot   = gammalib.GEnergy(1.0, 'TeV')
-            spectral = gammalib.GModelSpectralPlaw(3.0e-4, -1.5, epivot)
-            spectral['Prefactor'].min(1.0e-8)
+            # Handle LOOKUP model
+            if self['spatial'].string() == 'LOOKUP':
  
-            # Set background model
-            model = gammalib.GCTAModelBackground(spatial, spectral)
+                # Set spatial model
+                factor1 = gammalib.GCTAModelSpatialLookup(self['slufile'].filename())
  
-        # Handle GAUSS(E) model
-        elif self['spatial'].string() == 'GAUSS(E)':
+            # Handle GAUSS model
+            elif self['spatial'].string() == 'GAUSS':
  
-            # Set spatial model
-            if self['snumbins'].integer() == 1:
+                # Set spatial model
                 factor1 = gammalib.GCTAModelRadialGauss(3.0)
                 factor1['Sigma'].min(sigma_min)
                 factor1['Sigma'].max(sigma_max)
-            else:
-                emin     = gammalib.GEnergy(self['smin'].real(), 'TeV')
-                emax     = gammalib.GEnergy(self['smax'].real(), 'TeV')
-                energies = gammalib.GEnergies(self['snumbins'].integer(), emin, emax)
-                spectrum = gammalib.GModelSpectralConst(3.0)
-                nodes    = gammalib.GModelSpectralNodes(spectrum, energies)
-                for i in range(nodes.nodes()):
-                    nodes[i*2+1].min(sigma_min)
-                    nodes[i*2+1].max(sigma_max)
-                nodes.autoscale()
-                factor1 = gammalib.GCTAModelSpatialGaussSpectrum(nodes)
  
-            # Optionally add gradient
-            if self['gradient'].boolean():
-                spatial = gammalib.GCTAModelSpatialMultiplicative()
-                factor2 = gammalib.GCTAModelSpatialGradient()
-                spatial.append(factor1)
-                spatial.append(factor2)
-            else:
-                spatial = factor1
-
-            # Set spectral model
-            epivot   = gammalib.GEnergy(1.0, 'TeV')
-            spectral = gammalib.GModelSpectralPlaw(3.0e-4, -1.5, epivot)
-            spectral['Prefactor'].min(1.0e-8)
-
-            # Set background model
-            model = gammalib.GCTAModelBackground(spatial, spectral)
-
-        # Handle PROFILE model
-        elif self['spatial'].string() == 'PROFILE':
-
-            # Set spatial model
-            factor1 = gammalib.GCTAModelRadialProfile(2.0, 4.0, 5.0)
-            factor1['Width'].min(1.0)
-            factor1['Width'].max(10.0)
-            factor1['Core'].min(1.0)
-            factor1['Core'].max(10.0)
-            factor1['Tail'].min(-10.0)
-            factor1['Tail'].max(10.0)
-            factor1['Tail'].free()
-
-            # Optionally add gradient
-            if self['gradient'].boolean():
-                spatial = gammalib.GCTAModelSpatialMultiplicative()
-                factor2 = gammalib.GCTAModelSpatialGradient()
-                spatial.append(factor1)
-                spatial.append(factor2)
+            # Handle GAUSS(E) model
+            elif self['spatial'].string() == 'GAUSS(E)':
+
+                # Set spatial model. If a single energy bin is specified then
+                # use the GAUSS model. Otherwise allocate a node spectrum for
+                # the energy nodes of GAUSS(E).
+                if self['snumbins'].integer() == 1:
+                    factor1 = gammalib.GCTAModelRadialGauss(3.0)
+                    factor1['Sigma'].min(sigma_min)
+                    factor1['Sigma'].max(sigma_max)
+                else:
+                    emin     = gammalib.GEnergy(self['smin'].real(), 'TeV')
+                    emax     = gammalib.GEnergy(self['smax'].real(), 'TeV')
+                    energies = gammalib.GEnergies(self['snumbins'].integer(), emin, emax)
+                    spectrum = gammalib.GModelSpectralConst(3.0)
+                    nodes    = gammalib.GModelSpectralNodes(spectrum, energies)
+                    for i in range(nodes.nodes()):
+                        nodes[i*2+1].min(sigma_min)
+                        nodes[i*2+1].max(sigma_max)
+                    nodes.autoscale()
+                    factor1 = gammalib.GCTAModelSpatialGaussSpectrum(nodes)
+
+            # Handle PROFILE model
+            elif self['spatial'].string() == 'PROFILE':
+
+                # Set spatial model
+                factor1 = gammalib.GCTAModelRadialProfile(2.0, 4.0, 5.0)
+                factor1['Width'].min(1.0)
+                factor1['Width'].max(10.0)
+                factor1['Core'].min(1.0)
+                factor1['Core'].max(10.0)
+                factor1['Tail'].min(1.0)
+                factor1['Tail'].max(10.0)
+                factor1['Tail'].free()
+
+            # Handle POLYNOM model
+            elif self['spatial'].string() == 'POLYNOM':
+
+                # Set spatial model
+                factor1 = gammalib.GCTAModelRadialPolynom([1.0, -0.1, +0.1])
+                factor1['Coeff0'].min(0.1)
+                factor1['Coeff0'].max(10.0)
+                factor1['Coeff0'].fix()
+                factor1['Coeff1'].min(-10.0)
+                factor1['Coeff1'].max(10.0)
+                factor1['Coeff2'].min(-10.0)
+                factor1['Coeff2'].max(10.0)
+
+            # Any other strings (should never occur)
             else:
-                spatial = factor1
-
-            # Set spectral model
-            epivot   = gammalib.GEnergy(1.0, 'TeV')
-            spectral = gammalib.GModelSpectralPlaw(3.0e-4, -1.5, epivot)
-            spectral['Prefactor'].min(1.0e-8)
-
-            # Set background model
-            model = gammalib.GCTAModelBackground(spatial, spectral)
-
-        # Handle POLYNOM model
-        elif self['spatial'].string() == 'POLYNOM':
-
-            # Set spatial model
-            factor1 = gammalib.GCTAModelRadialPolynom([1.0, -0.1, +0.1])
-            factor1['Coeff0'].min(0.1)
-            factor1['Coeff0'].max(10.0)
-            factor1['Coeff0'].fix()
-            factor1['Coeff1'].min(-10.0)
-            factor1['Coeff1'].max(10.0)
-            factor1['Coeff2'].min(-10.0)
-            factor1['Coeff2'].max(10.0)
+                model = None
  
             # Optionally add gradient
             if self['gradient'].boolean():
@@ -323,10 +265,6 @@ class csbkgmodel(ctools.csobservation):
             # Set background model
             model = gammalib.GCTAModelBackground(spatial, spectral)
  
-        # Any other strings (should never occur)
-        else:
-            model = None
-
         # Set background name and instrument
         if model is not None:
             model.name('Background')
@@ -37,7 +37,8 @@ models = 1dc_howto.xml \
          crab_profile.xml \
          disk.xml \
          gauss.xml \
-         shell.xml
+         shell.xml \
+         hess_bkg_lookup.fits
  
 # Install models
 modelsdir   = $(datadir)/models
@@ -45,6 +45,7 @@ class Test(test):
         # Set test datasets and parameters
         self._myevents1 = self._datadir + '/crab_offaxis1.fits'
         self._myevents2 = self._datadir + '/crab_offaxis2.fits'
+        self._lookup    = self._datadir + '/hess_bkg_lookup.fits'
  
         # Return
         return
@@ -250,6 +251,114 @@ class Test(test):
         # Check background model
         self._check_bkg_model('csbkgmodel_py5.xml', nmodels=2)
  
+        # Test GAUSS(E) spatial model
+        bkgmodel = cscripts.csbkgmodel()
+        bkgmodel['inobs']      = self._events
+        bkgmodel['caldb']      = self._caldb
+        bkgmodel['irf']        = self._irf
+        bkgmodel['instrument'] = 'CTA'
+        bkgmodel['spatial']    = 'GAUSS(E)'
+        bkgmodel['snumbins']   = 2
+        bkgmodel['smin']       = 1.0
+        bkgmodel['smax']       = 10.0
+        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_py6.xml'
+        bkgmodel['logfile']    = 'csbkgmodel_py6.log'
+
+        # Run csbkgmodel script and save background model
+        bkgmodel.logFileOpen()   # Make sure we get a log file
+        bkgmodel.execute()
+
+        # Check background model
+        self._check_bkg_model('csbkgmodel_py6.xml')
+
+        # Test LOOKUP spatial model
+        bkgmodel = cscripts.csbkgmodel()
+        bkgmodel['inobs']      = self._events
+        bkgmodel['caldb']      = self._caldb
+        bkgmodel['irf']        = self._irf
+        bkgmodel['instrument'] = 'CTA'
+        bkgmodel['spatial']    = 'LOOKUP'
+        bkgmodel['slufile']    = self._lookup
+        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_py7.xml'
+        bkgmodel['logfile']    = 'csbkgmodel_py7.log'
+
+        # Run csbkgmodel script and save background model
+        bkgmodel.logFileOpen()   # Make sure we get a log file
+        bkgmodel.execute()
+
+        # Check background model
+        self._check_bkg_model('csbkgmodel_py7.xml')
+
+        # Test PROFILE spatial model
+        bkgmodel = cscripts.csbkgmodel()
+        bkgmodel['inobs']      = self._events
+        bkgmodel['caldb']      = self._caldb
+        bkgmodel['irf']        = self._irf
+        bkgmodel['instrument'] = 'CTA'
+        bkgmodel['spatial']    = 'PROFILE'
+        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_py8.xml'
+        bkgmodel['logfile']    = 'csbkgmodel_py8.log'
+
+        # Run csbkgmodel script and save background model
+        bkgmodel.logFileOpen()   # Make sure we get a log file
+        bkgmodel.execute()
+
+        # Check background model
+        self._check_bkg_model('csbkgmodel_py8.xml')
+
+        # Test POLYNOM spatial model
+        bkgmodel = cscripts.csbkgmodel()
+        bkgmodel['inobs']      = self._events
+        bkgmodel['caldb']      = self._caldb
+        bkgmodel['irf']        = self._irf
+        bkgmodel['instrument'] = 'CTA'
+        bkgmodel['spatial']    = 'POLYNOM'
+        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_py9.xml'
+        bkgmodel['logfile']    = 'csbkgmodel_py9.log'
+
+        # Run csbkgmodel script and save background model
+        bkgmodel.logFileOpen()   # Make sure we get a log file
+        bkgmodel.execute()
+
+        # Check background model
+        self._check_bkg_model('csbkgmodel_py9.xml')
+
         # Return
         return