From 31d34ecd3763ff910d5d3f6322291a99d13fced2 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Tue, 23 Apr 2019 15:18:51 +0200 Subject: [PATCH] Add "POW" energy binning to ctool::create_energies() method and csbkgmodel script --- cscripts/csbkgmodel.par | 3 ++- cscripts/csspec.py | 17 +++++++++++++---- doc/source/users/reference_manual/csbkgmodel.rst | 14 ++++++++++++-- src/support/ctool.cpp | 35 ++++++++++++++++++++++++++++++++++- 4 files changed, 61 insertions(+), 8 deletions(-) diff --git a/cscripts/csbkgmodel.par b/cscripts/csbkgmodel.par index f5f0593..552b67d 100755 --- a/cscripts/csbkgmodel.par +++ b/cscripts/csbkgmodel.par @@ -31,11 +31,12 @@ smin, r, a, 0.1,,, "Lower energy limit of energy nodes for GAUSS(E) spatia smax, r, a, 10.0,,, "Upper energy limit of energy nodes for GAUSS(E) spatial model (TeV)" gradient, b, a, yes,,, "Allow for spatial gradient?" spectral, s, a, NODES,PLAW|NODES,, "Spectral model component" -ebinalg, s, a, LOG,FILE|LIN|LOG,, "Algorithm for defining spectral energy nodes" +ebinalg, s, a, LOG,FILE|LIN|LOG|POW,, "Algorithm for defining spectral energy nodes" emin, r, a, 0.1,,, "Lower energy limit (TeV)" emax, r, a, 100.0,,, "Upper energy limit (TeV)" enumbins, i, a, 8,1,200, "Number of spectral energy nodes" ebinfile, f, a, NONE,,, "Name of the file containing the energy node definition" +ebingamma, r, a, 1.0,,, "Exponent for POW energy node definition" runwise, b, a, yes,,, "Generate runwise background model?" rad, r, h, 2.0,,, "Radius for event selection (degrees)" diff --git a/cscripts/csspec.py b/cscripts/csspec.py index c936ee3..9cdd52d 100755 --- a/cscripts/csspec.py +++ b/cscripts/csspec.py @@ -652,7 +652,7 @@ class csspec(ctools.csobservation): 'Npred': 0.0} # Write header for fitting - self._log_header3(gammalib.EXPLICIT, 'Performing fit') + self._log_header3(gammalib.EXPLICIT, 'Performing fit in energy bin') # Setup maximum likelihood fit like = ctools.ctlike(obs) @@ -667,6 +667,9 @@ class csspec(ctools.csobservation): # Perform maximum likelihood fit like.run() + # Write model results for explicit chatter level + self._log_string(gammalib.EXPLICIT, str(like.obs().models())) + # Continue only if log-likelihood is non-zero if like.obs().logL() != 0.0: @@ -688,7 +691,8 @@ class csspec(ctools.csobservation): if self['calc_ulim'].boolean(): # Logging information - self._log_header3(gammalib.EXPLICIT, 'Computing upper limit') + self._log_header3(gammalib.EXPLICIT, + 'Computing upper limit for energy bin') # Create upper limit object ulimit = ctools.ctulimit(like.obs()) @@ -801,7 +805,7 @@ class csspec(ctools.csobservation): self._log_header1(gammalib.TERSE, 'Generate spectrum') # Write header for fitting - self._log_header3(gammalib.EXPLICIT, 'Performing fit') + self._log_header3(gammalib.EXPLICIT, 'Performing model fit') # Perform maximum likelihood fit like = ctools.ctlike(self.obs()) @@ -816,6 +820,9 @@ class csspec(ctools.csobservation): spectrum = model.spectral() logL0 = like.obs().logL() + # 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()): @@ -844,7 +851,9 @@ class csspec(ctools.csobservation): if self['calc_ulim'].boolean(): # Logging information - self._log_header3(gammalib.EXPLICIT, 'Computing upper limit') + self._log_header3(gammalib.EXPLICIT, + 'Computing upper limit for node energy %f TeV' % + result['energy']) # Copy observation container obs = like.obs().copy() diff --git a/doc/source/users/reference_manual/csbkgmodel.rst b/doc/source/users/reference_manual/csbkgmodel.rst index 5f9a771..4a98de2 100755 --- a/doc/source/users/reference_manual/csbkgmodel.rst +++ b/doc/source/users/reference_manual/csbkgmodel.rst @@ -75,11 +75,17 @@ General parameters Spectral model component. ``PLAW`` specifies a simple power law model, ``NODES`` specifies a piecewise broken power-law. -``ebinalg [string]`` +``ebinalg [string]`` Algorithm for defining energy nodes. For ``FILE``, the energy nodes are defined in a FITS file that is specified by the ``ebinfile`` parameter, for ``LIN`` and ``LOG`` there will be ``enumbins`` energy nodes spaced linearly or logarithmically between ``emin`` and ``emax``, respectively. + For ``POW`` a spacing with ``enumbins`` nodes between ``emin`` and ``emax`` + will be defined so that the integral over a power law with a spectral index + defined by the parameter ``ebingamma`` will be constant. A value of + ``ebingamma=1`` correspond to the ``LOG`` algorithm, a smaller value provides + a finer binning for larger energies while a larger value provides a finer + binning for smaller energies. ``emin [real]`` Lower energy limit (in TeV). @@ -88,11 +94,15 @@ General parameters Upper energy limit (in TeV). ``enumbins [integer]`` - Number of energy nodes if ``LIN`` or ``LOG`` energy algorithms are used. + Number of energy nodes if ``LIN``, ``LOG`` or ``POW`` energy algorithms are + used. ``ebinfile [file]`` Name of the file containing the energy node definition if ``ebinalg=FILE``. +``ebingamma [real]`` + Power law index for energy node definition if ``ebinalg=POW``. + ``runwise [boolean]`` Generate runwise background model? If ``yes`` is specified a background model component for each observation in the observation definition XML diff --git a/src/support/ctool.cpp b/src/support/ctool.cpp index 9b206fa..a5bcb87 100755 --- a/src/support/ctool.cpp +++ b/src/support/ctool.cpp @@ -696,10 +696,11 @@ GEbounds ctool::create_ebounds(void) * The following parameters are read: * * ebinalg - Energy binning algorithm - * ebinfile - Name of file with energy binning + * ebinfile - Name of file with energy binning (if ebinalg == FILE) * emin - Minimum energy (if ebinalg != FILE) * emax - Maximum energy (if ebinalg != FILE) * enumbins - Number of energy bins (if ebinalg != FILE) + * ebingamma - Power law index (if ebinalg == POW) ***************************************************************************/ GEnergies ctool::create_energies(void) { @@ -783,6 +784,38 @@ GEnergies ctool::create_energies(void) } // endif: ebinalg was "FILE" + // ... otherwise, if energy binning algorithm is of type "POW" (case + // sensitive), then set the energy binning according to a power-law + else if (ebinalg == "POW") { + + // Get task parameters + double emin = (*this)["emin"].real(); + double emax = (*this)["emax"].real(); + int enumbins = (*this)["enumbins"].integer(); + double ebingamma = (*this)["ebingamma"].real(); + + // Precomputation + double a = 1.0 - ebingamma; + double c = (a == 0.0) + ? 1.0 / (std::log(emax) - std::log(emin)) + : a / (std::pow(emax,a) - std::pow(emin,a)); + double b = double(c*(enumbins-1.0)); + + // Initialise first energy + double e = emin; + energies.append(GEnergy(e, "TeV")); + + // Loop over energy bins + for (int i = 0; i < enumbins-1; ++i) { + double log_e_next = (a == 0.0) + ? 1.0/b + std::log(e) + : std::log(a/b + std::pow(e,a)) / a; + double e_next = std::exp(log_e_next); + energies.append(GEnergy(e_next, "TeV")); + e = e_next; + } + } + // ... otherwise use a linear or a logarithmically-spaced energy binning else { -- libgit2 0.21.4