Add "POW" energy binning to ctool::create_energies() method and csbkgmodel script

Jürgen Knödlseder
1 parent f9342fcb
Showing 4 changed files with 61 additions and 8 deletions
cscripts/csbkgmodel.par
cscripts/csspec.py
doc/source/users/reference_manual/csbkgmodel.rst
src/support/ctool.cpp
@@ -31,11 +31,12 @@ smin,       r, a, 0.1,,, &quot;Lower energy limit of energy nodes for GAUSS(E) spatia
 smax,       r, a, 10.0,,, "Upper energy limit of energy nodes for GAUSS(E) spatial model (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 spectral energy nodes"
+ebinalg,    s, a, LOG,FILE|LIN|LOG|POW,, "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 spectral energy nodes"
 ebinfile,   f, a, NONE,,, "Name of the file containing the energy node definition"
+ebingamma,  r, a, 1.0,,, "Exponent for POW energy node definition"
 runwise,    b, a, yes,,, "Generate runwise background model?"
 rad,        r, h, 2.0,,, "Radius for event selection (degrees)"
@@ -652,7 +652,7 @@ class csspec(ctools.csobservation):
                   'Npred':       0.0}
         # Write header for fitting
-        self._log_header3(gammalib.EXPLICIT, 'Performing fit')
+        self._log_header3(gammalib.EXPLICIT, 'Performing fit in energy bin')
         # Setup maximum likelihood fit
         like = ctools.ctlike(obs)
@@ -667,6 +667,9 @@ class csspec(ctools.csobservation):
         # Perform maximum likelihood fit
         like.run()
+        # Write model results for explicit chatter level
+        self._log_string(gammalib.EXPLICIT, str(like.obs().models()))
+
         # Continue only if log-likelihood is non-zero
         if like.obs().logL() != 0.0:
@@ -688,7 +691,8 @@ class csspec(ctools.csobservation):
             if self['calc_ulim'].boolean():
                 # Logging information
-                self._log_header3(gammalib.EXPLICIT, 'Computing upper limit')
+                self._log_header3(gammalib.EXPLICIT,
+                                  'Computing upper limit for energy bin')
                 # Create upper limit object  
                 ulimit = ctools.ctulimit(like.obs())
@@ -801,7 +805,7 @@ class csspec(ctools.csobservation):
         self._log_header1(gammalib.TERSE, 'Generate spectrum')
         # Write header for fitting
-        self._log_header3(gammalib.EXPLICIT, 'Performing fit')
+        self._log_header3(gammalib.EXPLICIT, 'Performing model fit')
         # Perform maximum likelihood fit
         like          = ctools.ctlike(self.obs())
@@ -816,6 +820,9 @@ class csspec(ctools.csobservation):
         spectrum = model.spectral()
         logL0    = like.obs().logL()
+        # Write model results for explicit chatter level
+        self._log_string(gammalib.EXPLICIT, str(like.obs().models()))
+
         # Loop over all nodes
         for i in range(spectrum.nodes()):
@@ -844,7 +851,9 @@ class csspec(ctools.csobservation):
             if self['calc_ulim'].boolean():
                 # Logging information
-                self._log_header3(gammalib.EXPLICIT, 'Computing upper limit')
+                self._log_header3(gammalib.EXPLICIT,
+                                  'Computing upper limit for node energy %f TeV' %
+                                  result['energy'])
                 # Copy observation container
                 obs = like.obs().copy()
@@ -75,11 +75,17 @@ General parameters
     Spectral model component. ``PLAW`` specifies a simple power law model,
     ``NODES`` specifies a piecewise broken power-law.
-``ebinalg <FILE|LIN|LOG> [string]``
+``ebinalg <FILE|LIN|LOG|POW> [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.
+    For ``POW`` a spacing with ``enumbins`` nodes between ``emin`` and ``emax``
+    will be defined so that the integral over a power law with a spectral index
+    defined by the parameter ``ebingamma`` will be constant. A value of
+    ``ebingamma=1`` correspond to the ``LOG`` algorithm, a smaller value provides
+    a finer binning for larger energies while a larger value provides a finer
+    binning for smaller energies.
 ``emin [real]``
     Lower energy limit (in TeV).
@@ -88,11 +94,15 @@ General parameters
     Upper energy limit (in TeV).
 ``enumbins [integer]``
-    Number of energy nodes if ``LIN`` or ``LOG`` energy algorithms are used.
+    Number of energy nodes if ``LIN``, ``LOG`` or ``POW`` energy algorithms are
+    used.
 ``ebinfile [file]``
     Name of the file containing the energy node definition if ``ebinalg=FILE``.
+``ebingamma [real]``
+    Power law index for energy node definition if ``ebinalg=POW``.
+
 ``runwise [boolean]``
     Generate runwise background model? If ``yes`` is specified a background
     model component for each observation in the observation definition XML
@@ -696,10 +696,11 @@ GEbounds ctool::create_ebounds(void)
  * The following parameters are read:
  *
  *      ebinalg - Energy binning algorithm
- *      ebinfile - Name of file with energy binning
+ *      ebinfile - Name of file with energy binning (if ebinalg == FILE)
  *      emin - Minimum energy (if ebinalg != FILE)
  *      emax - Maximum energy (if ebinalg != FILE)
  *      enumbins - Number of energy bins (if ebinalg != FILE)
+ *      ebingamma - Power law index (if ebinalg == POW)
  ***************************************************************************/
 GEnergies ctool::create_energies(void)
 {
@@ -783,6 +784,38 @@ GEnergies ctool::create_energies(void)
     } // endif: ebinalg was "FILE"
+    // ... otherwise, if energy binning algorithm is of type "POW" (case
+    // sensitive), then set the energy binning according to a power-law
+    else if (ebinalg == "POW") {
+
+        // Get task parameters
+    	double emin      = (*this)["emin"].real();
+    	double emax      = (*this)["emax"].real();
+    	int    enumbins  = (*this)["enumbins"].integer();
+        double ebingamma = (*this)["ebingamma"].real();
+
+        // Precomputation
+        double a = 1.0 - ebingamma;
+        double c = (a == 0.0)
+                   ? 1.0 / (std::log(emax) - std::log(emin))
+                   : a   / (std::pow(emax,a) - std::pow(emin,a));
+        double b = double(c*(enumbins-1.0));
+
+        // Initialise first energy
+        double e = emin;
+        energies.append(GEnergy(e, "TeV"));
+
+        // Loop over energy bins
+        for (int i = 0; i < enumbins-1; ++i) {
+            double log_e_next = (a == 0.0)
+                                ? 1.0/b + std::log(e)
+                                : std::log(a/b + std::pow(e,a)) / a;
+            double e_next = std::exp(log_e_next);
+            energies.append(GEnergy(e_next, "TeV"));
+            e = e_next;
+        }
+    }
+
     // ... otherwise use a linear or a logarithmically-spaced energy binning
     else {