diff --git a/ChangeLog b/ChangeLog index ee614f2..c89f6c3 100755 --- a/ChangeLog +++ b/ChangeLog @@ -3,6 +3,8 @@ * Bug fix version 1.7.4 released ================================ + Replace spectral model by file function in csspec (#3624) + 2020-12-02 diff --git a/NEWS b/NEWS index 945f13c..3cfba34 100755 --- a/NEWS +++ b/NEWS @@ -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 diff --git a/configure.ac b/configure.ac index 147026b..515280b 100755 --- a/configure.ac +++ b/configure.ac @@ -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) diff --git a/cscripts/csspec.py b/cscripts/csspec.py index cdd8b00..9606fda 100755 --- a/cscripts/csspec.py +++ b/cscripts/csspec.py @@ -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') diff --git a/doc/source/admin/release_history/1.7.4.rst b/doc/source/admin/release_history/1.7.4.rst index 0680714..131eae6 100644 --- a/doc/source/admin/release_history/1.7.4.rst +++ b/doc/source/admin/release_history/1.7.4.rst @@ -10,3 +10,5 @@ ctools 1.7.4 is a bug fix release for ctools release 1.7. The following bugs were fixed: +* [`3624 `_] - + Replace spectral model by file function in :ref:`csspec` diff --git a/doc/source/users/reference_manual/csspec.rst b/doc/source/users/reference_manual/csspec.rst index 2557e8d..2ee3ee5 100755 --- a/doc/source/users/reference_manual/csspec.rst +++ b/doc/source/users/reference_manual/csspec.rst @@ -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 ------------------