Add GAUSS(E) model to csbkgmodel (#2868)

Jürgen Knödlseder
1 parent c5d689f9
Showing 5 changed files with 79 additions and 15 deletions
ChangeLog
NEWS
cscripts/csbkgmodel.par
cscripts/csbkgmodel.py
doc/source/admin/release_history/1.6.rst
-2019-04-03
+2019-04-14
  
         * Version 1.6.0 released
         ========================
  
+        Add GAUSS(E) model to csbkgmodel (#2868)
         Move ctool::set_obs_bounds() to ctobservation::set_obs_bounds()
         Log warning if edisp=no (#2842)
         Fix ring background computation in ctskymap for small roiradius (#2839)
 New Features and Important Changes in ctools 1.6.0
  
-3 April 2019
+14 April 2019
  
  
 Introduction
@@ -270,6 +270,9 @@ None
  
 csbkgmodel - Generate background model for 3D analysis
 ------------------------------------------------------
+A model named GAUSS(E) was added to implement an energy-dependent Gaussian
+spatial function (#2868).
+
 Introduced the script that generates a background for a 3D analysis. This
 script is particularly useful to generate a background model for the analysis
 of the H.E.S.S. public data release (#2694).
@@ -25,13 +25,16 @@ 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,, "Spatial model component"
+spatial,    s, a, GAUSS,IRF|AEFF|GAUSS|GAUSS(E),, "Spatial model component"
+snumbins,   i, a, 3,1,200, "Number of energy nodes for sigma(E)"
+smin,       r, a, 0.1,,, "Lower energy limit for sigma(E) energy nodes (TeV)"
+smax,       r, a, 10.0,,, "Upper energy limit for sigma(E) energy nodes (TeV)"
 gradient,   b, a, yes,,, "Allow for spatial gradient?"
 spectral,   s, a, NODES,PLAW|NODES,, "Spectral model component"
-ebinalg,    s, a, LOG,FILE|LIN|LOG,, "Algorithm for defining energy nodes"
+ebinalg,    s, a, LOG,FILE|LIN|LOG,, "Algorithm for defining spectral energy nodes"
 emin,       r, a, 0.1,,, "Lower energy limit (TeV)"
 emax,       r, a, 100.0,,, "Upper energy limit (TeV)"
-enumbins,   i, a, 8,1,200, "Number of energy nodes"
+enumbins,   i, a, 8,1,200, "Number of spectral energy nodes"
 ebinfile,   f, a, NONE,,, "Name of the file containing the energy node definition"
 runwise,    b, a, yes,,, "Generate runwise background model?"
 rad,        r, h, 2.0,,, "Radius for event selection (degrees)"
@@ -95,10 +95,15 @@ class csbkgmodel(ctools.csobservation):
  
         # Query input parameters
         spatial = self['spatial'].string()
-        if spatial == 'GAUSS':
+        if spatial == 'GAUSS(E)':
+            snumbins = self['snumbins'].integer()
+            if snumbins > 1:
+                self['smin'].real()
+                self['smax'].real()
+        if spatial == 'GAUSS' or spatial == 'GAUSS(E)':
             self['gradient'].boolean()
-        self['spectral'].string()
-        if spatial == 'NODES':
+        spectral = self['spectral'].string()
+        if spectral == 'NODES':
             self._create_energies()
         self['runwise'].boolean()
         self['emin'].real()
@@ -142,10 +147,17 @@ class csbkgmodel(ctools.csobservation):
         # Return instrument
         return instrument
  
-    def _generate_initial_model(self):
+    def _generate_initial_model(self, sigma_min=2.0, sigma_max=1000.0):
         """
         Generate initial background model
  
+        Parameters
+        ----------
+        sigma_min : float, optional
+            Minimum sigma
+        sigma_max : float, optional
+            Maximum sigma
+
         Returns
         -------
         model : `~gammalib.GModelData()`
@@ -165,22 +177,65 @@ class csbkgmodel(ctools.csobservation):
  
         # 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
+
+            # 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(E) model
+        elif self['spatial'].string() == 'GAUSS(E)':
+
+            # Set spatial model
+            if self['snumbins'].integer() == 1:
                 factor1 = gammalib.GCTAModelRadialGauss(3.0)
-                factor1['Sigma'].min(1.5)
-                factor1['Sigma'].max(10.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 = gammalib.GCTAModelRadialGauss(3.0)
-                spatial['Sigma'].min(1.5)
-                spatial['Sigma'].max(10.0)
+                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)
-            model    = gammalib.GCTAModelBackground(spatial, spectral)
+
+            # Set background model
+            model = gammalib.GCTAModelBackground(spatial, spectral)
  
         # Any other strings (should never occur)
         else:
@@ -46,6 +46,8 @@ Bug fixes
 Improvements
 ------------
  
+* [`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>`_] -
   Log warning if ``edisp=no``
 * Add ``max_iter`` parameter to :ref:`ctbutterfly`, :ref:`cterror`,