Add BINS method to csspec (#3606)

Jürgen Knödlseder
1 parent 51307080
Showing 7 changed files with 125 additions and 34 deletions
ChangeLog
NEWS
README.md
cscripts/csspec.par
cscripts/csspec.py
doc/source/admin/release_history/2.0.rst
doc/source/users/reference_manual/csspec.rst
-2021-03-05
+2021-04-13
  
         * Version 2.0.0 released
         ========================
  
+        Add BINS method to csspec (#3606)
         Add hidden "fix_spat" parameter to cttsmap (#3562)
         Add optimizer setter to ctlikelihood base class
         Add Sphinx tutorials test script (#3510)
 New Features and Important Changes in ctools 2.0.0
  
-5 March 2021
+13 April 2021
  
  
 Introduction
@@ -285,7 +285,9 @@ None
  
 csspec - Generates Spectral Energy Distribution
 -----------------------------------------------
-None
+Added BINS method that replaces the spectral model by a bin function that is
+implemented by the GammaLib class GModelSpectralBins. Using a bin function
+assures that spectral points are statistically independent (#3606).
  
  
 cssrcdetect - Detects sources in sky map
 ctools information
 ==================
-* Version:             2.0.0.dev (5 March 2021)
+* Version:             2.0.0.dev (13 April 2021)
 * GammaLib dependency: 2.0.0.dev
  
 [![Build Status](https://cta-jenkins.irap.omp.eu/buildStatus/icon?job=ctools-integrate-os)](https://cta-jenkins.irap.omp.eu/job/ctools-integrate-os/)
@@ -31,7 +31,7 @@ outfile,   f, a, spectrum.fits,,, &quot;Output spectrum file&quot;
 #
 # Script parameters
 #==================
-method,        s, a, AUTO,SLICE|NODES|AUTO,, "Spectrum generation method"
+method,        s, a, AUTO,SLICE|NODES|BINS|AUTO,, "Spectrum generation method"
 ebinalg,       s, a, LOG,FILE|LIN|LOG|POW,, "Algorithm for defining energy bins"
 emin,          r, a, 0.1,,, "Start value for first energy bin in TeV"
 emax,          r, a, 100.0,,, "Stop value for last energy bin in TeV"
@@ -43,6 +43,7 @@ calc_ts,       b, h, yes,,, &quot;Compute TS value in each bin&quot;
 calc_ulim,     b, h, yes,,, "Compute upper limit in each bin"
 fix_srcs,      b, h, yes,,, "Fix other skymodel parameters"
 fix_bkg,       b, h, no,,, "Fix background parameters"
+bingamma,      r, h, -2.0,,, "Spectral index for BINS method"
 dll_sigstep,   r, h, 0.0,0.0,100.0, "Step size in standard deviations for log-like profiles"
 dll_sigmax,    r, h, 5.0,1.0,100.0, "Maximum number of standard deviations for log-like profiles"
 dll_freenodes, b, h, no,,, "Free nodes not being fit when computing log-like profiles"
@@ -144,9 +144,9 @@ class csspec(ctools.csobservation):
         n_cta   = n_unbinned + n_binned + n_onoff
         n_other = self.obs().size() - n_cta
  
-        # If spectrum method is not "NODES" then set spectrum method and
-        # script mode according to type of observations
-        if self._method != 'NODES':
+        # If spectrum method is not "NODES" or 'BINS" then set spectrum method
+        # and script mode according to type of observations
+        if self._method != 'NODES' and self._method != 'BINS':
             if n_other > 0:
                 self._method = 'NODES'
             else:
@@ -178,6 +178,7 @@ class csspec(ctools.csobservation):
         self['calc_ts'].boolean()
         self['fix_bkg'].boolean()
         self['fix_srcs'].boolean()
+        self['bingamma'].real()
  
         # Setup dlog-likelihood parameters
         self['dll_sigstep'].real()
@@ -213,10 +214,10 @@ class csspec(ctools.csobservation):
         elif n_other > 0:
             if self._method == 'SLICE':
                 msg = 'Selected "SLICE" method but none of the observations ' \
-                      'is a CTA observation. Please select "AUTO" or "NODES" ' \
-                      'if no CTA observation is provided.'
+                      'is a CTA observation. Please select "AUTO", "NODES" ' \
+                      'or "BINS" if no CTA observation is provided.'
                 raise RuntimeError(msg)
-            else:
+            elif self._method == 'AUTO':
                 self._method = 'NODES'
  
         # Log selected spectrum method
@@ -456,6 +457,64 @@ class csspec(ctools.csobservation):
         # Return
         return
  
+    def _set_replace_src_spectrum_by_bins(self):
+        """
+        Replace source spectrum by bin function
+        """
+        # Initialise model container
+        models = gammalib.GModels()
+
+        # Loop over model containers
+        for model in self.obs().models():
+
+            # If we deal with source model then replace the spectral model
+            # by a bin function
+            if model.name() == self['srcname'].string():
+
+                # Get spectral bins index
+                bingamma = self['bingamma'].real()
+
+                # Setup spectral bin function
+                spectrum = gammalib.GModelSpectralBins(model.spectral(),
+                                                       self._ebounds,
+                                                       bingamma)
+                spectrum.autoscale()
+
+                # Make sure that all bins are positive. Autoscale all
+                # parameters so that their nominal value is unity.
+                for i in range(spectrum.bins()):
+                    name    = 'Intensity%d' % i
+                    par     = spectrum[name]
+                    par.autoscale()
+                    value   = par.value()
+                    minimum = 1.0e-20 * value
+                    if minimum <= 0.0:
+                        minimum = 1.0e-20
+                        if minimum < value:
+                            value = minimum
+                    par.value(value)
+                    par.min(minimum)
+
+                # Set spectral component of source model
+                model.spectral(spectrum)
+
+                # Make sure that TS computation is disabled (makes computation
+                # faster)
+                model.tscalc(False)
+
+                # Append model
+                models.append(model)
+
+            # ... otherwise just append model
+            else:
+                models.append(model)
+
+        # Put new model in observation containers
+        self.obs().models(models)
+
+        # Return
+        return
+
     def _select_onoff_obs(self, obs, emin, emax):
         """
         Select an energy interval from one CTA On/Off observation
@@ -854,8 +913,8 @@ class csspec(ctools.csobservation):
         # 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()):
+        # Loop over all nodes or bins
+        for i in range(self._ebounds.size()):
  
             # Get energy boundaries
             emin      = self._ebounds.emin(i)
@@ -878,23 +937,26 @@ class csspec(ctools.csobservation):
  
             # Convert differential flux and flux error to nuFnu
             norm               = elogmean.MeV() * elogmean.MeV()  * gammalib.MeV2erg
-            result['flux']     = spectrum[i*2+1].value()
-            result['e2dnde']   = spectrum[i*2+1].value() * norm
-            result['flux_err'] = spectrum[i*2+1].error()
+            result['flux']     = spectrum.intensity(i)
+            result['e2dnde']   = spectrum.intensity(i) * norm
+            result['flux_err'] = spectrum.error(i)
  
             # Compute upper flux limit
             ulimit_value = -1.0
-            parname = 'Intensity%d' % i
+            parname      = 'Intensity%d' % i
  
             # Compute delta log-likelihood
             if self['dll_sigstep'].real() > 0.0:
  
                 # Extract the normalization scan values
-                (norm_scan, dlogL, loglike) = self._profile_logL(like.copy(), parname, elogmean)
+                (norm_scan, dlogL, loglike) = self._profile_logL(like.copy(),
+                                                                 parname,
+                                                                 elogmean)
                 result['norm_scan'] = norm_scan
                 result['dloglike']  = dlogL
                 result['logL']      = loglike
  
+            # Compute upper limit
             if self['calc_ulim'].boolean():
  
                 # Logging information
@@ -1006,7 +1068,8 @@ class csspec(ctools.csobservation):
         source.tscalc(False)
  
         # Fix all parameters in the spectral model
-        if (self._method != 'NODES') or (not self['dll_freenodes'].boolean()):
+        if (self._method != 'NODES' and self._method != 'BINS') or \
+           (not self['dll_freenodes'].boolean()):
             for par in spectral:
                 par.fix()
  
@@ -1197,6 +1260,15 @@ class csspec(ctools.csobservation):
             # Fit model
             results = self._fit_model()
  
+        # Case C: "BINS" method
+        elif self._method == 'BINS':
+
+            # Replace source spectrum by bins function
+            self._set_replace_src_spectrum_by_bins()
+
+            # Fit model
+            results = self._fit_model()
+
         # Create FITS file
         self._create_fits(results)
  
 .. _2.0:
  
-ctools 2.0 release (5 March 2021)
-=================================
+ctools 2.0 release (13 April 2021)
+==================================
  
 Overview
 --------
@@ -24,6 +24,8 @@ Bug fixes
 Improvements
 ------------
  
+* [`3606 <https://cta-redmine.irap.omp.eu/issues/3606>`_] -
+  Add ``BINS`` method to :ref:`csspec`
 * [`3562 <https://cta-redmine.irap.omp.eu/issues/3562>`_] -
   Add hidden ``fix_spat`` parameter to :ref:`cttsmap`
 * Add optimizer setter to ``ctlikelihood`` base class.
@@ -19,15 +19,20 @@ significance of the source detection. Optionally, it also computes an upper
 flux limit that is particularly useful in case that the source is not
 significantly detected within a spectral bin (hidden parameter ``calc_ulim``).
  
-There are two fundamental methods to run the script depending on the value of
+There are three fundamental methods to run the script depending on the value of
 the ``method`` parameter. If ``method=SLICE`` the energy interval defined by the
 ``emin`` and ``emax`` parameters is divided into a number of energy bins and an
-independent maximum likelihood fit is performed in each of the energy bins,
-while if ``method=NODES`` the spectral model will be replaced by a node function
-that will be fit to all data. The latter is particularly useful if non-CTA data
-should be fitted such as data from Fermi/LAT or COMPTEL. By default, ``method`` is
-set to ``AUTO`` which will automatically select the method based on the input data.
-For CTA-only observations ``SLICE`` will be used, otherwise ``NODES`` will be used.
+independent maximum likelihood fit is performed in each of the energy bins. If
+``method=NODES`` the spectral model will be replaced by a node function that will
+be fit to all data. Alternatively, if ``method=BINS`` the spectral model will be
+replaced by a bin function that will be fit to all data. While with the ``NODES``
+method the spectral points are in general not statistically independent, the
+``BINS`` method assures that the spectral points are statistically independent.
+The ``NODES`` and ``BINS`` methods are particularly useful if data from instruments
+other than CTA should be analyses, such as data from Fermi/LAT, COMPTEL or
+INTEGRAL/SPI. By default, ``method`` is set to ``AUTO`` which will automatically
+select the method based on the input data. For CTA-only observations ``SLICE`` will
+be used, otherwise ``NODES`` will be used.
  
 The spectral binning is either defined by a FITS file containing the energy
 boundaries of each bin (option ``ebinalg=FILE``) or as ``enumbins`` bins spread
@@ -95,13 +100,17 @@ General parameters
 ``outfile [file]``
     Output spectrum FITS file.
  
-``method <SLICE|NODES|AUTO> [string]``
+``method <SLICE|NODES|BINS|AUTO> [string]``
     Spectrum generation method.
-    ``SLICE`` will slice energy interval into ``enumbins`` independent bins,
-    ``NODES`` will replace spectral model by node function with ``enumbins``
-    nodes, and ``AUTO`` will automatically select the method based on the input
-    data. For CTA-only observations ``SLICE`` will be used, otherwise ``NODES``
-    will be used.
+    ``SLICE`` will slice the energy interval into ``enumbins`` independent
+    bins,
+    ``NODES`` will replace the spectral model by a node function with
+    ``enumbins`` nodes,
+    ``BINS`` will replace the spectral model by a bin function with
+    ``enumbins`` nodes,
+    and ``AUTO`` will automatically select the method based on the input
+    data. If ``AUTO`` is selected and there are only CTA observations
+    then ``SLICE`` will be used, otherwise ``NODES`` will be used.
  
 ``ebinalg <FILE|LIN|LOG|POW> [string]``
     Algorithm for defining energy bins. For ``FILE``, the energy bins are defined
@@ -151,6 +160,10 @@ General parameters
 ``(fix_bkg = no) [boolean]``
     Fix background model parameters?
  
+``(bingamma = -2.0) [real]``
+    Spectral index for ``BINS`` method. Within each spectral bin a power law
+    will be assumed with the specified spectral index.
+
 ``(dll_sigstep = 0) [real]``
     Step size in standard deviations for log-like profiles. Note this value 
     determines the number of steps to use.