Replace spectral model by file function in csspec (#3624)

Jürgen Knödlseder
1 parent d8f2d169
Showing 6 changed files with 49 additions and 15 deletions
ChangeLog
NEWS
configure.ac
cscripts/csspec.py
doc/source/admin/release_history/1.7.4.rst
doc/source/users/reference_manual/csspec.rst
@@ -3,6 +3,8 @@
         * Bug fix version 1.7.4 released
         ================================
+        Replace spectral model by file function in csspec (#3624)
+
 2020-12-02
@@ -256,7 +256,10 @@ None
 csspec - Generates Spectral Energy Distribution
 -----------------------------------------------
-None
+The script now transforms the spectral model of the source of interest into a
+file function, using the model parameters specified in the input model. This
+allows to determine spectral energy distributions for combined spectral models
+(#3624).
 cssrcdetect - Detects sources in sky map
@@ -54,7 +54,7 @@ AC_CONFIG_HEADER([config.h])
 # 5. If any interfaces have been removed or changed since the last public   #
 #    release, then set age to 0.                                            #
 #############################################################################
-CTOOLS_LT_VERSION="7:1:0"
+CTOOLS_LT_VERSION="7:2:0"
 AC_SUBST(CTOOLS_LT_VERSION)
@@ -2,7 +2,7 @@
 # ==========================================================================
 # Generates a spectrum.
 #
-# Copyright (C) 2014-2020 Michael Mayer
+# Copyright (C) 2014-2021 Michael Mayer
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -20,6 +20,7 @@
 # ==========================================================================
 import sys
 import math
+import tempfile    # Kludge for file function generation
 import gammalib
 import ctools
 from cscripts import mputils
@@ -366,6 +367,24 @@ class csspec(ctools.csobservation):
                                          ' Fixing "'+par.name()+'"')
                     par.fix()
+                # Kludge: convert spectral model into a file function
+                nbins       = self._ebounds.size()
+                num         = nbins * 20
+                energies    = gammalib.GEnergies(num, self._ebounds.emin(0),
+                                                      self._ebounds.emax(nbins-1))
+                _, filename = tempfile.mkstemp()
+                f = open(filename, 'wb')
+                spectral    = model.spectral()
+                for energy in energies:
+                    value = spectral.eval(energy)
+                    f.write('%20.6f  %20e\n' % (energy.MeV(), value))
+                f.close()
+                spectral = gammalib.GModelSpectralFunc(filename, 1.0)
+                model.spectral(spectral)
+                self._log_string(gammalib.EXPLICIT, ' Converting spectral model '
+                                 'into file function via temporary file "%s"' %
+                                 filename)
+
                 # Free the normalisation parameter which is assumed to be
                 # the first spectral parameter
                 normpar = model.spectral()[0]
@@ -895,6 +914,7 @@ class csspec(ctools.csobservation):
                 result['dloglike']  = dlogL
                 result['logL']      = loglike
+            # Compute upper flux limit
             if self['calc_ulim'].boolean():
                 # Logging information
@@ -1009,7 +1029,7 @@ class csspec(ctools.csobservation):
         if (self._method != 'NODES') or (not self['dll_freenodes'].boolean()):
             for par in spectral:
                 par.fix()
-        
+
         # Re-compute the log-likelihood
         like.run()
         loglike = like.obs().logL()
@@ -1130,7 +1150,7 @@ class csspec(ctools.csobservation):
         table.append(dnde)
         table.append(TSvalues)
         table.append(Npred_values)
-        
+
         # Define the SED type
         table.card('SED_TYPE', 'norm,e2dnde,dnde,npred', 'SED type')
@@ -10,3 +10,5 @@ ctools 1.7.4 is a bug fix release for ctools release 1.7.
 The following bugs were fixed:
+* [`3624 <https://cta-redmine.irap.omp.eu/issues/3624>`_] -
+  Replace spectral model by file function in :ref:`csspec`
@@ -20,14 +20,24 @@ 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
-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
+the ``method`` parameter, ``SLICE`` and ``NODES``.
+
+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. In order to cope with
+all kinds of spectral models, including composite spectra, input spectrum will
+be converted into a file function, with the normalisation of the file function
+being the fit parameter for each energy bin. This implies that spectral
+components are merged into a single function using the model parameters that
+are specified in the input model.
+
+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.
+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.
 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
@@ -54,9 +64,6 @@ the fitted flux and flux error, the Test Statistics value (option
 ``calc_ts=yes``), the upper flux limit (option ``calc_ulim=yes``) and the
 predicted number of events (only for unbinned data).
-.. warning::
-   The upper limit computation is not yet implemented for the ``NODES`` method.
-
 General parameters
 ------------------