diff --git a/ChangeLog b/ChangeLog index ca67c15..1e03971 100644 --- a/ChangeLog +++ b/ChangeLog @@ -3,6 +3,7 @@ * Version 1.5.0 released ======================== + Add "statistic" parameter to all likelihood tools (#2237) Rename ctlike "stat" parameter to "statistic" (#2237) Enable handling of On/Off observations in csspec (#2231) Add csphagen script (#2208) diff --git a/NEWS b/NEWS index cfc9303..c15edcc 100644 --- a/NEWS +++ b/NEWS @@ -127,6 +127,11 @@ ctulimit - Compute upper limit Add "statistic" parameter to allow setting of optimizer statistic (#2237). +support - Support classes +------------------------- +Add ctobservation.obs() method to set an observation container. + + cscaldb - Inspect calibration database -------------------------------------- None @@ -238,6 +243,10 @@ None csspec - Generates Spectral Energy Distribution ----------------------------------------------- +The csspec script now derived from the cslikelihood base class. A "statistic" +parameter was added to allow setting the optimizer statistic for the spectrum +generation (#2237). + 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" diff --git a/cscripts/csspec.par b/cscripts/csspec.par index 0b78e9d..302478b 100644 --- a/cscripts/csspec.par +++ b/cscripts/csspec.par @@ -37,6 +37,7 @@ emax, r, a, 100.0,,, "Upper energy limit (TeV)" enumbins, i, a, 20,,, "Number of energy bins" ebinalg, s, a, LOG,LIN|LOG|FILE,, "Binning algorithm" ebinfile, f, a, NONE,,, "Name of the file containing the energy bin definition" +statistic, s, h, DEFAULT,DEFAULT|CSTAT|WSTAT|CHI2,, "Optimization statistic" 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" diff --git a/cscripts/csspec.py b/cscripts/csspec.py index d7e0df0..0b97b8c 100755 --- a/cscripts/csspec.py +++ b/cscripts/csspec.py @@ -26,12 +26,12 @@ import ctools # ============ # # csspec class # # ============ # -class csspec(ctools.cscript): +class csspec(ctools.cslikelihood): """ Generates a spectrum - This class implements the generation of a spectral energy distribution - from data. + This class implements the generation of a Spectral Energy Distribution + (SED) from gamma-ray observations. """ # Constructors and destructors @@ -39,9 +39,8 @@ class csspec(ctools.cscript): """ Constructor """ - # Set name - self._name = 'csspec' - self._version = ctools.__version__ + # Initialise application by calling the appropriate class constructor + self._init_cslikelihood('csspec', ctools.__version__, argv) # Initialise data members self._ebounds = gammalib.GEbounds() @@ -50,12 +49,6 @@ class csspec(ctools.cscript): self._onoff_mode = False self._method = 'AUTO' - # Initialise observation container from constructor arguments - self._obs, argv = self._set_input_obs(argv) - - # Initialise application by calling the appropriate class constructor - self._init_cscript(argv) - # Return return @@ -73,13 +66,16 @@ class csspec(ctools.cscript): Get parameters from parfile and setup the observation """ # Set observation if not done before - if self._obs == None or self._obs.size() == 0: + if self.obs().size() == 0: self._require_inobs('csspec::get_parameters()') - self._obs = self._get_observations() + self.obs(self._get_observations()) + + # Set observation statistic + self._set_obs_statistic(gammalib.toupper(self['statistic'].string())) # Set models if we have none - if self._obs.models().size() == 0: - self._obs.models(self['inmodel'].filename()) + if self.obs().models().size() == 0: + self.obs().models(self['inmodel'].filename()) # Query source name self['srcname'].string() @@ -92,7 +88,7 @@ class csspec(ctools.cscript): n_unbinned = 0 n_binned = 0 n_onoff = 0 - for obs in self._obs: + for obs in self.obs(): if obs.classname() == 'GCTAObservation': if obs.eventtype() == 'CountsCube': n_binned += 1 @@ -101,7 +97,7 @@ class csspec(ctools.cscript): elif obs.classname() == 'GCTAOnOffObservation': n_onoff += 1 n_cta = n_unbinned + n_binned + n_onoff - n_other = self._obs.size() - n_cta + n_other = self.obs().size() - n_cta # If spectrum method is not "NODES" then set spectrum method and # script mode according to type of CTA observations @@ -194,9 +190,9 @@ class csspec(ctools.cscript): # should work reasonably well since for each observation the # relevant energy bins will be selected. if self._binned_mode: - cube_ebounds = self._obs[0].events().ebounds() + cube_ebounds = self.obs()[0].events().ebounds() else: - cube_ebounds = self._obs[0].rmf().etrue() + cube_ebounds = self.obs()[0].rmf().etrue() # Loop over user energy boundaries and collect all energy bins # that overlap @@ -244,7 +240,7 @@ class csspec(ctools.cscript): ' Specify overlapping energy range.' raise RuntimeError(msg) - # Unbinned mode + # Unbinned mode else: self._ebounds = self._create_ebounds() @@ -260,15 +256,15 @@ class csspec(ctools.cscript): # Log counts cube energy range for binned mode if self._binned_mode: - cube_ebounds = self._obs[0].events().ebounds() + cube_ebounds = self.obs()[0].events().ebounds() value = '%s - %s' % (str(cube_ebounds.emin()), str(cube_ebounds.emax())) self._log_value(gammalib.TERSE, 'Counts cube energy range', value) # Log RMF energy range for On/Off mode elif self._onoff_mode: - etrue = self._obs[0].rmf().etrue() - ereco = self._obs[0].rmf().emeasured() + etrue = self.obs()[0].rmf().etrue() + ereco = self.obs()[0].rmf().emeasured() vtrue = '%s - %s' % (str(etrue.emin()), str(etrue.emax())) vreco = '%s - %s' % (str(ereco.emin()), str(ereco.emax())) self._log_value(gammalib.TERSE, 'RMF true energy range', vtrue) @@ -291,7 +287,7 @@ class csspec(ctools.cscript): self._log_header1(gammalib.TERSE, 'Adjust model parameters') # Adjust model parameters dependent on input user parameters - for model in self._obs.models(): + for model in self.obs().models(): # Initialise TS flag for all models to false model.tscalc(False) @@ -350,7 +346,7 @@ class csspec(ctools.cscript): models = gammalib.GModels() # Loop over model containers - for model in self._obs.models(): + for model in self.obs().models(): # If we deal with source model then replace the spectral model # by a node function @@ -360,7 +356,7 @@ class csspec(ctools.cscript): 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() @@ -388,7 +384,7 @@ class csspec(ctools.cscript): models.append(model) # Put new model in observation containers - self._obs.models(models) + self.obs().models(models) # Return return @@ -463,9 +459,11 @@ class csspec(ctools.cscript): rmf[idst_true, idst_reco] = obs.rmf()[isrc_true, isrc_reco] # Set On/Off observations - id = obs.id() + id = obs.id() + statistic = obs.statistic() obs = gammalib.GCTAOnOffObservation(pha_on, pha_off, arf, rmf) obs.id(id) + obs.statistic(statistic) # Return observation return obs @@ -493,7 +491,7 @@ class csspec(ctools.cscript): self._log_header3(gammalib.EXPLICIT, 'Filtering cube') # Select layers - cubemask = ctools.ctcubemask(self._obs) + cubemask = ctools.ctcubemask(self.obs()) cubemask['regfile'] = 'NONE' cubemask['ra'] = 'UNDEFINED' cubemask['dec'] = 'UNDEFINED' @@ -516,12 +514,12 @@ class csspec(ctools.cscript): # Loop over all input observations and select energy bins for # all On/Off observations - for run in self._obs: + for run in self.obs(): if run.classname() == 'GCTAOnOffObservation': obs.append(self._select_onoff_obs(run, emin, emax)) # Append models - obs.models(self._obs.models()) + obs.models(self.obs().models()) # Use ctselect for unbinned analysis else: @@ -530,7 +528,7 @@ class csspec(ctools.cscript): self._log_header3(gammalib.EXPLICIT, 'Selecting events') # Select events - select = ctools.ctselect(self._obs) + select = ctools.ctselect(self.obs()) select['ra'] = 'UNDEFINED' select['dec'] = 'UNDEFINED' select['rad'] = 'UNDEFINED' @@ -707,7 +705,7 @@ class csspec(ctools.cscript): self._log_header3(gammalib.EXPLICIT, 'Performing fit') # Perform maximum likelihood fit - like = ctools.ctlike(self._obs) + like = ctools.ctlike(self.obs()) like['edisp'] = self['edisp'].boolean() like.run() @@ -856,7 +854,7 @@ class csspec(ctools.cscript): self._get_parameters() # Write input observation container into logger - self._log_observations(gammalib.NORMAL, self._obs, 'Input observation') + self._log_observations(gammalib.NORMAL, self.obs(), 'Input observation') # Write spectral binning into logger self._log_spectral_binning() @@ -931,7 +929,7 @@ class csspec(ctools.cscript): # Set default name is user name is empty if not name: - user_name = self._name + user_name = self._name() else: user_name = name @@ -959,7 +957,7 @@ class csspec(ctools.cscript): Set model """ # Copy models - self._obs.models(models.clone()) + self.obs().models(models.clone()) # Return return diff --git a/doc/source/users/reference_manual/csspec.rst b/doc/source/users/reference_manual/csspec.rst index f39a2da..9c179c1 100644 --- a/doc/source/users/reference_manual/csspec.rst +++ b/doc/source/users/reference_manual/csspec.rst @@ -68,7 +68,7 @@ General parameters ``caldb [string]`` Calibration database. - + ``irf [string]`` Instrumental response function. @@ -89,19 +89,27 @@ General parameters ``emin [real]`` Lower energy limit of events (in TeV). - + ``emax [real]`` Upper energy limit of events (in TeV). - + ``enumbins [integer]`` Number of energy bins. - + ``ebinalg [string]`` Algorithm for defining energy bins. - + ``ebinfile [file]`` Name of the file containing the energy bin definition. +``(statistic = DEFAULT) [string]`` + Optimization statistic. ``DEFAULT`` uses the default statistic for all + observations, which is ``CSTAT`` or the statistic specified in the + observation definition XML file. ``CSTAT`` uses the C statistic for + all observations, ``WSTAT`` uses the W statistic for all On/Off + observations, and ``CHI2`` uses the Chi squared statistic for all + binned or stacked observations. + ``(calc_ts = yes) [boolean]`` Compute TS for each spectral point? @@ -124,21 +132,23 @@ Standard parameters ``(chatter = 2) [integer]`` Verbosity of the executable: ``chatter = 0``: no information will be logged - + ``chatter = 1``: only errors will be logged - + ``chatter = 2``: errors and actions will be logged - + ``chatter = 3``: report about the task execution - + ``chatter = 4``: detailed report about the task execution - + ``(clobber = yes) [boolean]`` - Specifies whether an existing source spectrum output file should be overwritten. - + Specifies whether an existing source spectrum output file should be + overwritten. + ``(debug = no) [boolean]`` - Enables debug mode. In debug mode the executable will dump any log file output to the console. - + Enables debug mode. In debug mode the executable will dump any log file + output to the console. + ``(mode = ql) [string]`` Mode of automatic parameters (default is "ql", i.e. "query and learn"). diff --git a/pyext/ctlikelihood.i b/pyext/ctlikelihood.i index c446cc5..fb0bd4c 100644 --- a/pyext/ctlikelihood.i +++ b/pyext/ctlikelihood.i @@ -141,15 +141,23 @@ public: # The function supports either an observation container, or an argument # list or no argument as "argv" parameter. The function also writes the # header in the log file and switches the date on for logging. -def _init_cslikelihood(self, argv): +def _init_cslikelihood(self, name, version, argv): if len(argv) > 0 and isinstance(argv[0],gammalib.GObservations): - cslikelihood.__init__(self, self._name, self._version, argv[0]) + cslikelihood.__init__(self, name, version, argv[0]) elif len(argv) > 0: - cslikelihood.__init__(self, self._name, self._version, *argv) + cslikelihood.__init__(self, name, version, *argv) else: - cslikelihood.__init__(self, self._name, self._version) + cslikelihood.__init__(self, name, version) # Set logger properties self._log_header() self._log.date(True) cslikelihood._init_cslikelihood = _init_cslikelihood + +# Execute the script +def _execute(self): + self.logFileOpen() + self._read_ahead(True) + self.run() + self.save() +cslikelihood.execute = _execute %} diff --git a/pyext/ctobservation.i b/pyext/ctobservation.i index 67e9966..23e3578 100644 --- a/pyext/ctobservation.i +++ b/pyext/ctobservation.i @@ -88,6 +88,7 @@ public: virtual void save(void) = 0; // Public methods + void obs(const GObservations& obs); const GObservations& obs(void) const; void write_ogip_keywords(GFitsHDU* hdu) const; @@ -142,19 +143,19 @@ public: # The function supports either an observation container, or an argument # list or no argument as "argv" parameter. The function also writes the # header in the log file and switches the date on for logging. -def _init_csobservation(self, argv): +def _init_csobservation(self, name, version, argv): if len(argv) > 0 and isinstance(argv[0],gammalib.GObservations): - csobservation.__init__(self, self._name, self._version, argv[0]) + csobservation.__init__(self, name, version, argv[0]) elif len(argv) > 0: - csobservation.__init__(self, self._name, self._version, *argv) + csobservation.__init__(self, name, version, *argv) else: - csobservation.__init__(self, self._name, self._version) + csobservation.__init__(self, name, version) # Set logger properties self._log_header() self._log.date(True) csobservation._init_csobservation = _init_csobservation -# Define an iterator over all observations +# Define an iterator over all unbinned observations def _unbinned_observations(self): obs = self._first_unbinned_observation() while obs != None: @@ -162,4 +163,12 @@ def _unbinned_observations(self): obs = self._next_unbinned_observation() ctool._unbinned_observations = _unbinned_observations cscript._unbinned_observations = _unbinned_observations + +# Execute the script +def _execute(self): + self.logFileOpen() + self._read_ahead(True) + self.run() + self.save() +csobservation.execute = _execute %} diff --git a/src/support/ctlikelihood.cpp b/src/support/ctlikelihood.cpp index 9604a1d..f7869a8 100644 --- a/src/support/ctlikelihood.cpp +++ b/src/support/ctlikelihood.cpp @@ -305,6 +305,9 @@ double ctlikelihood::evaluate(GModelPar& par, const double& value) * For On/Off CTA observations of type GCTAOnOffObservation the fit statistic * is always set. Possible values are @c POISSON, @c CSTAT or @c WSTAT (case * insensitive). + * + * If @p statistic is any other string the fit statistic of the observations + * will not be modified. ***************************************************************************/ void ctlikelihood::set_obs_statistic(const std::string& statistic) {