diff --git a/ChangeLog b/ChangeLog index 07f42ed..70709ba 100644 --- a/ChangeLog +++ b/ChangeLog @@ -9,6 +9,7 @@ Add keywords to sky maps produced by ctskymap (#2194) Enhance efficiency of computations in ctmodel (#2189) Add csmodelsois script (#2180) + Add "NODES" method to csspec to support non-CTA analysis (#1927) Use t-type parameters for all "tmin" and "tmax" parameters (#1864) diff --git a/NEWS b/NEWS index 3add308..96a53b5 100644 --- a/NEWS +++ b/NEWS @@ -140,6 +140,13 @@ for classical On/Off IACT analysis (#2208). The csspec script now can also handle On/Off spectra for SED generation (#2231). +The csspec script now can also handle non-CTA data. For that purpose a "method" +parameter has been introduced that allows to selected between slicing of the data +(method="SLICE"), which is only supported for CTA observations, and replacing of +the spectral model by a node function (method="NODES"). If method="AUTO" the tool +will use "SLICE" for observations containing only CTA data and "NODES" for observations +containing other data (#1927). + Examples -------- diff --git a/cscripts/csspec.par b/cscripts/csspec.par index ff3409d..0b78e9d 100644 --- a/cscripts/csspec.par +++ b/cscripts/csspec.par @@ -31,6 +31,7 @@ outfile, f, a, spectrum.fits,,, "Output spectrum file" # # Script parameters #================== +method, s, a, AUTO,SLICE|NODES|AUTO,, "Spectrum generation method" emin, r, a, 0.1,,, "Lower energy limit (TeV)" emax, r, a, 100.0,,, "Upper energy limit (TeV)" enumbins, i, a, 20,,, "Number of energy bins" diff --git a/cscripts/csspec.py b/cscripts/csspec.py index a4e4169..014ad6a 100755 --- a/cscripts/csspec.py +++ b/cscripts/csspec.py @@ -29,7 +29,7 @@ import ctools class csspec(ctools.cscript): """ Generates a spectrum - + This class implements the generation of a spectral energy distribution from data. """ @@ -48,6 +48,7 @@ class csspec(ctools.cscript): self._fits = None self._binned_mode = False self._onoff_mode = False + self._method = 'AUTO' # Initialise observation container from constructor arguments self._obs, argv = self._set_input_obs(argv) @@ -75,16 +76,6 @@ class csspec(ctools.cscript): if self._obs == None or self._obs.size() == 0: self._require_inobs('csspec::get_parameters()') self._obs = self._get_observations() - # if we have one binned CTA observation, i.e. if we are in - # binned mode - #self._binned_mode = False - #self._onoff_mode = False - #if self._obs.size() == 1: - # if self._obs[0].classname() == 'GCTAObservation': - # if self._obs[0].eventtype() == 'CountsCube': - # self._binned_mode = True - # elif self._obs[0].classname() == 'GCTAOnOffObservation': - # self._onoff_mode = True # Set models if we have none if self._obs.models().size() == 0: @@ -106,6 +97,8 @@ class csspec(ctools.cscript): n_unbinned += 1 elif obs.classname() == 'GCTAOnOffObservation': n_onoff += 1 + n_cta = n_unbinned + n_binned + n_onoff + n_other = self._obs.size() - n_cta # Set script mode according to CTA observations if n_unbinned == 0 and n_binned != 0 and n_onoff == 0: @@ -126,6 +119,9 @@ class csspec(ctools.cscript): # Set ebounds self._set_ebounds() + # Get method + self._method = self['method'].string() + # Query other parameeters self['edisp'].boolean() self['calc_ulim'].boolean() @@ -140,6 +136,35 @@ class csspec(ctools.cscript): # Write input parameters into logger self._log_parameters(gammalib.TERSE) + # Write spectrum method header and parameters + self._log_header1(gammalib.TERSE, 'Spectrum method') + self._log_value(gammalib.TERSE, 'Unbinned CTA observations', n_unbinned) + self._log_value(gammalib.TERSE, 'Binned CTA observations', n_binned) + self._log_value(gammalib.TERSE, 'On/off CTA observations', n_onoff) + self._log_value(gammalib.TERSE, 'Other observations', n_other) + + # If there are a mix of CTA and non-CTA observations and the method + # is 'SLICE' then log a warning that non-CTA observations will be + # ignored + if n_cta > 0 and n_other > 0 and self._method == 'SLICE': + self._log_string(TERSE, '\nWARNING: Only CTA observation can be ' + 'handled with the "SLICE" method, all ' + 'non-CTA observation will be ignored.') + + # If there are only non-CTA observations and the method is 'SLICE' + # then stop now + 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.' + raise RuntimeError(msg) + else: + self._method = 'NODES' + + # Log selected spectrum method + self._log_value(gammalib.TERSE, 'Selected spectrum method', self._method) + # Return return @@ -311,6 +336,56 @@ class csspec(ctools.cscript): # Return return + def _set_replace_src_spectrum_by_nodes(self): + """ + Replace source spectrum by node 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 node function + if model.name() == self['srcname'].string(): + + # Setup energies at log mean energies of bins + energies = gammalib.GEnergies() + for i in range(self._ebounds.size()): + energies.append(self._ebounds.elogmean(i)) + + # Setup spectral node function + spectrum = gammalib.GModelSpectralNodes(model.spectral(), energies) + spectrum.autoscale() + + # Make sure that all nodes are positive + for i in range(spectrum.nodes()): + par = spectrum[i*2+1] + value = par.value() + min = 1.0e-30 * value + if min <= 0.0: + min = 1.0e-50 + if min < value: + value = min + par.value(value) + par.min(min) + + # Set spectral component of source model + model.spectral(spectrum) + + # 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): """ @@ -444,10 +519,10 @@ class csspec(ctools.cscript): # Use ctselect for unbinned analysis else: - + # Write header self._log_header3(gammalib.EXPLICIT, 'Selecting events') - + # Select events select = ctools.ctselect(self._obs) select['ra'] = 'UNDEFINED' @@ -458,7 +533,7 @@ class csspec(ctools.cscript): select['tmin'] = 'UNDEFINED' select['tmax'] = 'UNDEFINED' select.run() - + # Retrieve observation obs = select.obs().copy() @@ -610,6 +685,95 @@ class csspec(ctools.cscript): # Return results return results + def _fit_model(self): + """ + Fit model to observations + + Returns + ------- + results : list of dict + List of dictionaries with fit results + """ + # Write header + self._log_header1(gammalib.TERSE, 'Generate spectrum') + + # Write header for fitting + self._log_header3(gammalib.EXPLICIT, 'Performing fit') + + # Perform maximum likelihood fit + like = ctools.ctlike(self._obs) + like['edisp'] = self['edisp'].boolean() + like.run() + + # Initialise fit results + results = [] + + # Extract fit results + model = like.obs().models()[self['srcname'].string()] + spectrum = model.spectral() + ts = model.ts() + + # Loop over all nodes + for i in range(spectrum.nodes()): + + # Get energy boundaries + emin = self._ebounds.emin(i) + emax = self._ebounds.emax(i) + elogmean = self._ebounds.elogmean(i) + + # Initialise dictionary + result = {'energy': elogmean.TeV(), + 'energy_low': (elogmean - emin).TeV(), + 'energy_high': (emax - elogmean).TeV(), + 'flux': 0.0, + 'flux_err': 0.0, + 'TS': 0.0, + 'ulimit': 0.0, + 'Npred': 0.0} + + # Convert differential flux and flux error to nuFnu + norm = elogmean.MeV() * elogmean.MeV() * gammalib.MeV2erg + result['flux'] = spectrum[i*2+1].value() * norm + result['flux_err'] = spectrum[i*2+1].error() * norm + + # Compute TS + if self['calc_ts'].boolean(): + + # Copy observation container + obs = like.obs().copy() + + # Set intensity of node to tiny value + par = obs.models()[self['srcname'].string()].spectral()[i*2+1] + value = 1.0e-30 * par.value() + if par.min() > value: + par.min(value) + par.value(value) + par.fix() + + # Perform maximum likelihood fit + tslike = ctools.ctlike(obs) + tslike['edisp'] = self['edisp'].boolean() + tslike.run() + + # Store Test Statistic + model = tslike.obs().models()[self['srcname'].string()] + result['TS'] = ts - model.ts() + + # Log information + value = '%e +/- %e' % (result['flux'], result['flux_err']) + if self['calc_ulim'].boolean() and result['ulimit'] > 0.0: + value += ' [< %e]' % (result['ulimit']) + value += ' erg/cm2/s' + if self['calc_ts'].boolean() and result['TS'] > 0.0: + value += ' (TS = %.3f)' % (result['TS']) + self._log_value(gammalib.TERSE, 'Bin '+str(i+1), value) + + # Append results + results.append(result) + + # Return results + return results + def _create_fits(self, results): """ Create FITS file @@ -688,14 +852,26 @@ class csspec(ctools.cscript): # Write input observation container into logger self._log_observations(gammalib.NORMAL, self._obs, 'Input observation') + # Write spectral binning into logger + self._log_spectral_binning() + # Adjust model parameters dependent on input user parameters self._adjust_model_pars() - # Write spectral binning into logger - self._log_spectral_binning() + # Case A: "SLICE" method + if self._method == 'SLICE': - # Fit energy bins - results = self._fit_energy_bins() + # Fit energy bins + results = self._fit_energy_bins() + + # Case B: "NODES" method + elif self._method == 'NODES': + + # Replace source spectrum by nodes function + self._set_replace_src_spectrum_by_nodes() + + # Fit model + results = self._fit_model() # Create FITS file self._create_fits(results) @@ -738,7 +914,7 @@ class csspec(ctools.cscript): # Get outmap parameter outfile = self['outfile'].filename() - + # Log file name self._log_value(gammalib.NORMAL, 'Spectrum file', outfile.url()) diff --git a/doc/source/admin/release_history/1.5.rst b/doc/source/admin/release_history/1.5.rst index a7bdd67..67b2ea1 100644 --- a/doc/source/admin/release_history/1.5.rst +++ b/doc/source/admin/release_history/1.5.rst @@ -31,6 +31,8 @@ Improvements Add keywords to sky maps produced by :ref:`ctskymap` * [`2189 `_] - Enhance efficiency of computations in :ref:`ctmodel` +* [`1927 `_] - + Add ``NODES`` method to :ref:`csspec` to support non-CTA analysis New features diff --git a/doc/source/users/reference_manual/csspec.rst b/doc/source/users/reference_manual/csspec.rst index 2f47c6b..f39a2da 100644 --- a/doc/source/users/reference_manual/csspec.rst +++ b/doc/source/users/reference_manual/csspec.rst @@ -79,8 +79,13 @@ General parameters Name of the source in the source model XML file which should be used for sensitivity computation. -``outfile [file]`` - Output source spectrum file. +``outfile [file]`` + Spectrum generation method. + ``SLICE``: Slice energy interval into ``enumbins`` independent bins. + ``NODES``: Replace spectral model by node function with ``enumbins`` nodes. + ``AUTO`` : Automatically select method based on input data. Will use ``SLICE`` + for CTA-only observations and ``NODES`` for observations including + non-CTA data. ``emin [real]`` Lower energy limit of events (in TeV).