Use ctool::create_energies() method in csbkgmodel

Jürgen Knödlseder
1 parent e35c6b25
Showing 3 changed files with 33 additions and 15 deletions
cscripts/csbkgmodel.par
cscripts/csbkgmodel.py
doc/source/users/reference_manual/csbkgmodel.rst
@@ -25,7 +25,12 @@ outmodel, f, a, model.xml,,, &quot;Output model definition XML file&quot;
 instrument, s, a, CTA,,, "Instrument name (for multi-instrument analysis)"
 spatial,    s, a, GAUSS,IRF|AEFF|GAUSS,, "Spatial model component"
 gradient,   b, a, yes,,, "Allow for spatial gradient?"
-spectral,   s, a, HESS,PLAW|HESS,, "Spectral model component"
+spectral,   s, a, NODES,PLAW|NODES,, "Spectral model component"
+ebinalg,    s, a, LOG,FILE|LIN|LOG,, "Algorithm for defining energy nodes"
+emin,       r, a, 0.1,,, "First node energy (TeV)"
+emax,       r, a, 100.0,,, "Last node energy (TeV)"
+enumbins,   i, a, 10,1,200, "Number of 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)"
@@ -65,6 +65,8 @@ class csbkgmodel(ctools.csobservation):
         if self['spatial'].string() == 'GAUSS':
             self['gradient'].boolean()
         self['spectral'].string()
+        if self['spatial'].string() == 'NODES':
+            self._create_energies()
         self['runwise'].boolean()
         self['rad'].real()
@@ -204,16 +206,8 @@ class csbkgmodel(ctools.csobservation):
         # Extract spectral model of background component
         spectrum = model.spectral()
-        # Define 8 node energies
-        energies = gammalib.GEnergies()
-        energies.append(gammalib.GEnergy(0.1, 'TeV'))
-        energies.append(gammalib.GEnergy(0.2, 'TeV'))
-        energies.append(gammalib.GEnergy(0.4, 'TeV'))
-        energies.append(gammalib.GEnergy(0.8, 'TeV'))
-        energies.append(gammalib.GEnergy(1.6, 'TeV'))
-        energies.append(gammalib.GEnergy(3.2, 'TeV'))
-        energies.append(gammalib.GEnergy(6.4, 'TeV'))
-        energies.append(gammalib.GEnergy(12.8, 'TeV'))
+        # Define node energies
+        energies  = self._create_energies()
         # Create node spectrum
         nodes = gammalib.GModelSpectralNodes(spectrum, energies)
@@ -221,7 +215,7 @@ class csbkgmodel(ctools.csobservation):
         # Set minimum node intensities (this avoids NaNs during model fitting)
         for i in range(nodes.size()):
-            nodes[i].min(1.0e-6*nodes[i].value())
+            nodes[i].min(1.0e-10*nodes[i].value())
         # Set node spectrum
         model.spectral(nodes)
@@ -270,9 +264,9 @@ class csbkgmodel(ctools.csobservation):
         # Extract fitted model
         model = like.obs().models()[0].copy()
-        # If a HESS model is requested then refit a node spectrum
-        if self['spectral'].string() == 'HESS':
-        
+        # If a NODES model is requested then refit a node spectrum
+        if self['spectral'].string() == 'NODES':
+
             # Create nodes spectrum from fitted initial model
             model = self._create_nodes(model)
@@ -42,6 +42,25 @@ General parameters
     ``HESS`` specifies a piecewise broken power-law using 8 energy nodes
     spread over the energy range from 100 GeV to 12.8 TeV.
+``ebinalg <FILE|LIN|LOG> [string]``
+    Algorithm for defining energy nodes. For ``FILE``, the energy nodes are
+    defined in a FITS file that is specified by the ``ebinfile`` parameter,
+    for ``LIN`` and ``LOG`` there will be ``enumbins`` energy nodes spaced
+    linearly or logarithmically between ``emin`` and ``emax``, respectively.
+
+``emin [real]``
+    First node energy (in TeV) if ``LIN`` or ``LOG`` energy algorithms are
+    used.
+
+``emax [real]``
+    Last node energy (in TeV) if ``LIN`` or ``LOG`` energy algorithms are used.
+
+``enumbins [integer]``
+    Number of energy nodes if ``LIN`` or ``LOG`` energy algorithms are used.
+
+``ebinfile [file]``
+    Name of the file containing the energy node definition if ``ebinalg=FILE``.
+
 ``runwise [boolean]``
     Generate runwise background model? If ``yes`` is specified a background
     model component for each observation in the observation definition XML