diff --git a/ChangeLog b/ChangeLog index 57ad70b..dc48462 100755 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -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) diff --git a/NEWS b/NEWS index 1f43670..67c241a 100755 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ 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 diff --git a/README.md b/README.md index 5a87f6a..a12531b 100755 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ 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/) diff --git a/cscripts/csspec.par b/cscripts/csspec.par index e333e97..4edbbae 100755 --- a/cscripts/csspec.par +++ b/cscripts/csspec.par @@ -31,7 +31,7 @@ outfile, f, a, spectrum.fits,,, "Output spectrum file" # # 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,,, "Compute TS value in each bin" 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" diff --git a/cscripts/csspec.py b/cscripts/csspec.py index b8e854e..8a75044 100755 --- a/cscripts/csspec.py +++ b/cscripts/csspec.py @@ -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) diff --git a/doc/source/admin/release_history/2.0.rst b/doc/source/admin/release_history/2.0.rst index 5c1ab24..1d91353 100644 --- a/doc/source/admin/release_history/2.0.rst +++ b/doc/source/admin/release_history/2.0.rst @@ -1,7 +1,7 @@ .. _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 `_] - + Add ``BINS`` method to :ref:`csspec` * [`3562 `_] - Add hidden ``fix_spat`` parameter to :ref:`cttsmap` * Add optimizer setter to ``ctlikelihood`` base class. diff --git a/doc/source/users/reference_manual/csspec.rst b/doc/source/users/reference_manual/csspec.rst index 2557e8d..453ef7d 100755 --- a/doc/source/users/reference_manual/csspec.rst +++ b/doc/source/users/reference_manual/csspec.rst @@ -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 [string]`` +``method [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 [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.