Add more spatial models to csbkgmodel (#2878)

Jürgen Knödlseder
1 parent ea84b18f
Showing 7 changed files with 123 additions and 8 deletions
ChangeLog
NEWS
cscripts/csbkgmodel.par
cscripts/csbkgmodel.py
doc/source/admin/release_history/1.6.rst
doc/source/users/reference_manual/csbkgmodel.rst
test/test_csbkgmodel.py
-2019-04-14
+2019-04-26
         * Version 1.6.0 released
         ========================
+        Add more spatial models to csbkgmodel (#2878)
         Add GAUSS(E) model to csbkgmodel (#2868)
         Move ctool::set_obs_bounds() to ctobservation::set_obs_bounds()
         Log warning if edisp=no (#2842)
 New Features and Important Changes in ctools 1.6.0
-14 April 2019
+26 April 2019
 Introduction
@@ -270,6 +270,11 @@ None
 csbkgmodel - Generate background model for 3D analysis
 ------------------------------------------------------
+Three additional spatial models were added: LOOKUP which uses a 2D lookup table,
+PROFILE which uses an energy-independent profile, and POLYNOM which uses an
+energy-independent polynom of third order. The parameter "slufile" was added
+to allow specification of a lookup table file (#2878).
+
 A model named GAUSS(E) was added to implement an energy-dependent Gaussian
 spatial function (#2868).
@@ -25,7 +25,8 @@ outmodel, f, a, model.xml,,, &quot;Output model definition XML file&quot;
 # Script parameters
 #==================
 instrument, s, a, CTA,,, "Instrument name (for multi-instrument analysis)"
-spatial,    s, a, GAUSS,IRF|AEFF|GAUSS|GAUSS(E),, "Spatial model component"
+spatial,    s, a, GAUSS,IRF|AEFF|LOOKUP|GAUSS|GAUSS(E)|PROFILE|POLYNOM,, "Spatial model component"
+slufile,    f, a, lookup.fits,,, "Name of lookup table file for LOOKUP spatial model"
 snumbins,   i, a, 3,1,200, "Number of energy nodes for GAUSS(E) spatial model"
 smin,       r, a, 0.1,,, "Lower energy limit of energy nodes for GAUSS(E) spatial model (TeV)"
 smax,       r, a, 10.0,,, "Upper energy limit of energy nodes for GAUSS(E) spatial model (TeV)"
@@ -100,7 +100,10 @@ class csbkgmodel(ctools.csobservation):
             if snumbins > 1:
                 self['smin'].real()
                 self['smax'].real()
-        if spatial == 'GAUSS' or spatial == 'GAUSS(E)':
+        if spatial == 'LOOKUP':
+            self['slufile'].filename()
+        if spatial == 'LOOKUP' or spatial == 'GAUSS' or \
+           spatial == 'GAUSS(E)' or spatial == 'PROFILE':
             self['gradient'].boolean()
         spectral = self['spectral'].string()
         if spectral == 'NODES':
@@ -175,6 +178,29 @@ 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':
@@ -237,6 +263,66 @@ class csbkgmodel(ctools.csobservation):
             # 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)
+            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)
+
+            # 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)
+
         # Any other strings (should never occur)
         else:
             model = None
@@ -46,6 +46,8 @@ Bug fixes
 Improvements
 ------------
+* [`2878 <https://cta-redmine.irap.omp.eu/issues/2878>`_] -
+  Add more spatial models to :ref:`csbkgmodel`
 * [`2868 <https://cta-redmine.irap.omp.eu/issues/2868>`_] -
   Add ``GAUSS(E)`` model to :ref:`csbkgmodel`
 * [`2842 <https://cta-redmine.irap.omp.eu/issues/2842>`_] -
@@ -50,9 +50,29 @@ General parameters
     instrument, the instrument name will be extracted from the observation
     definition XML file.
-``spatial <IRF|AEFF|GAUSS|GAUSS(E)> [string]``
-    Spatial model component. So far IRF template, effective area, a radial
-    Gaussian and an energy-dependent radial Gaussian models are supported.
+``spatial <IRF|AEFF|LOOKUP|GAUSS|GAUSS(E)|PROFILE|POLYNOM> [string]``
+    Spatial model component. The following options exist:
+     ``IRF``: The background template included in the Instrument Response
+     Functions will be used (energy-dependent)
+
+     ``AEFF``: The effective area will be used (energy-dependent)
+
+     ``LOOKUP``: A lookup table, specified by the ``slufile`` parameter, will
+     be used (energy-dependent).
+
+     ``GAUSS``: A radial Gaussian in offset angle squared will be used.
+
+     ``GAUSS(E)``: An energy-dependent radial Gaussian in offset angle squared
+     will be used. The ``snumbins`` parameter specifies the number of energy
+     nodes, the ``smin`` and ``smax`` parameters specify the energy range for
+     the nodes.
+
+     ``PROFILE``: A radial profile will be used.
+
+     ``POLYNOM``: A third order polynomial will be used.
+
+``slufile [file]``
+    Name of the file containing the lookup table if ``spatial=LOOKUP``.
 ``snumbins [integer]``
     Number of energy nodes for GAUSS(E) spatial model.
@@ -198,7 +198,7 @@ class Test(test):
         bkgmodel['caldb']      = self._caldb
         bkgmodel['irf']        = self._irf
         bkgmodel['instrument'] = 'CTA'
-        bkgmodel['spatial']    = 'AEFF'
+        bkgmodel['spatial']    = 'IRF'
         bkgmodel['gradient']   = False
         bkgmodel['spectral']   = 'PLAW'
         bkgmodel['ebinalg']    = 'LOG'