From 0c70d44e9ef42b58a3ed881ec2b18a58da57607b Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Thu, 21 Jul 2016 18:31:44 +0200 Subject: [PATCH] Add "addbounds" parameter to response cube tools --- cscripts/cspull.py | 15 +++++++++++---- doc/source/users/reference_manual/ctbkgcube.rst | 13 ++++++++++--- doc/source/users/reference_manual/ctedispcube.rst | 16 +++++++++------- doc/source/users/reference_manual/ctexpcube.rst | 20 +++++++++++++------- doc/source/users/reference_manual/ctpsfcube.rst | 15 ++++++++------- src/ctbkgcube/ctbkgcube.cpp | 135 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------------- src/ctbkgcube/ctbkgcube.hpp | 5 ++++- src/ctbkgcube/ctbkgcube.par | 11 ++++++----- src/ctedispcube/ctedispcube.cpp | 113 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++----------- src/ctedispcube/ctedispcube.hpp | 5 ++++- src/ctedispcube/ctedispcube.par | 11 ++++++----- src/ctexpcube/ctexpcube.cpp | 122 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------- src/ctexpcube/ctexpcube.hpp | 7 +++++-- src/ctexpcube/ctexpcube.par | 11 ++++++----- src/ctpsfcube/ctpsfcube.cpp | 113 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------- src/ctpsfcube/ctpsfcube.hpp | 5 ++++- src/ctpsfcube/ctpsfcube.par | 13 +++++++------ src/support/ctool.cpp | 75 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/support/ctool.hpp | 3 +++ 19 files changed, 607 insertions(+), 101 deletions(-) diff --git a/cscripts/cspull.py b/cscripts/cspull.py index 655d33c..ade0ae2 100755 --- a/cscripts/cspull.py +++ b/cscripts/cspull.py @@ -20,6 +20,7 @@ # ========================================================================== import sys import csv +import math import gammalib import ctools from cscripts import obsutils @@ -163,6 +164,12 @@ class cspull(ctools.cscript): self._log('\n') self._log.header1('Compute stacked response') + # Set number of energy bins to at least 30 per energy decade + enumbins = int((math.log10(self['emax'].real()) - + math.log10(self['emin'].real())) * 30.0) + if self._enumbins > enumbins: + enumbins = self._enumbins + # Get stacked exposure expcube = ctools.ctexpcube(self._obs) expcube['incube'] = 'NONE' @@ -171,7 +178,7 @@ class cspull(ctools.cscript): expcube['binsz'] = self._binsz expcube['nxpix'] = self._npix expcube['nypix'] = self._npix - expcube['enumbins'] = self._enumbins + expcube['enumbins'] = enumbins expcube['emin'] = self['emin'].real() expcube['emax'] = self['emax'].real() expcube['coordsys'] = self._coordsys @@ -195,7 +202,7 @@ class cspull(ctools.cscript): psfcube['binsz'] = self._binsz*10.0 psfcube['nxpix'] = self._npix/10 psfcube['nypix'] = self._npix/10 - psfcube['enumbins'] = self._enumbins + psfcube['enumbins'] = enumbins psfcube['emin'] = self['emin'].real() psfcube['emax'] = self['emax'].real() psfcube['coordsys'] = self._coordsys @@ -220,7 +227,7 @@ class cspull(ctools.cscript): edispcube['binsz'] = self._binsz*10.0 edispcube['nxpix'] = self._npix/10 edispcube['nypix'] = self._npix/10 - edispcube['enumbins'] = self._enumbins + edispcube['enumbins'] = enumbins edispcube['emin'] = self['emin'].real() edispcube['emax'] = self['emax'].real() edispcube['coordsys'] = self._coordsys @@ -244,7 +251,7 @@ class cspull(ctools.cscript): bkgcube['binsz'] = self._binsz bkgcube['nxpix'] = self._npix bkgcube['nypix'] = self._npix - bkgcube['enumbins'] = self._enumbins + bkgcube['enumbins'] = enumbins bkgcube['emin'] = self['emin'].real() bkgcube['emax'] = self['emax'].real() bkgcube['coordsys'] = self._coordsys diff --git a/doc/source/users/reference_manual/ctbkgcube.rst b/doc/source/users/reference_manual/ctbkgcube.rst index babe5cf..94ab50d 100644 --- a/doc/source/users/reference_manual/ctbkgcube.rst +++ b/doc/source/users/reference_manual/ctbkgcube.rst @@ -25,10 +25,14 @@ cube as ``incube`` parameter. This will instruct ctbkgcube to extract the background cube definition (such as sky coordinates and projection, number of pixels, pixel scale, energy binning) from the counts cube. +Note that there should at least be 25 bins per energy decade in the background +cube to assure that the energy dependence of the background rate is sufficiently +well sampled (in particular at low energies). + ctbkgcube generates a background cube FITS file comprising two extensions. The primary extension contains a 3-dimensional image that contains the -background cube values. The next extension named ``EBOUNDS`` contains a -binary table that defines the energy boundaries of the background cube. +background cube values. The next extension named ``ENERGIES`` contains a +binary table that defines the energies of the background cube. ctbkgcube generates also an output model XML file that can serve as input for a maximum likelihood analysis. The output model XML file is a copy of @@ -82,7 +86,10 @@ General parameters ``ebinfile [file]`` Name of the file containing the energy bin definition. - + +``(addbounds = yes) [boolean]`` + Add energies to the background cube at the observation energy boundaries? + ``(usepnt = no) [boolean]`` Use CTA pointing direction for cube centre instead of xref/yref parameters? diff --git a/doc/source/users/reference_manual/ctedispcube.rst b/doc/source/users/reference_manual/ctedispcube.rst index bd9ff2b..1bf7127 100644 --- a/doc/source/users/reference_manual/ctedispcube.rst +++ b/doc/source/users/reference_manual/ctedispcube.rst @@ -29,13 +29,12 @@ rather coarse spatial binning to keep the file size manageable (with a typical value of ``binsz=1.0``). ctedispcube generates an energy dispersion cube FITS file comprising three -extensions. The primary extension contains a 3-dimensional image that contains -the energy disperison values. The energy and fraction dimensions of the enrgy -dispersion cube are folded into the 3rd dimension of the FITS image. The next -extension named ``ENERGIES`` contains a binary table that defines the energy -values of the energy dispersion cube. The last extension named ``MIGRAS`` -contains a binary table that defines the migration values of the cube which -is the ratio between reconstructed and true photon energy. +extensions. The primary extension contains a 4-dimensional image that contains +the energy disperison values. The next extension named ``ENERGIES`` contains +a binary table that defines the energies of the energy dispersion cube. The +last extension named ``MIGRAS`` contains a binary table that defines the +migration values of the cube which is the ratio between reconstructed and true +photon energy. General parameters @@ -71,6 +70,9 @@ General parameters ``ebinfile [file]`` Name of the file containing the energy bin definition. +``(addbounds = yes) [boolean]`` + Add energies to the energy dispersion cube at the observation energy boundaries? + ``(usepnt = no) [boolean]`` Use CTA pointing direction for cube centre instead of xref/yref parameters? diff --git a/doc/source/users/reference_manual/ctexpcube.rst b/doc/source/users/reference_manual/ctexpcube.rst index 2855d63..be60228 100644 --- a/doc/source/users/reference_manual/ctexpcube.rst +++ b/doc/source/users/reference_manual/ctexpcube.rst @@ -20,7 +20,7 @@ that has been used in the generation of the counts cube using :doc:`ctbin`. The binning of an exposure cube does not need to correspond to the binning of a counts cube, as exposure values will in any case be determined by linear -interpolation. For the user's convenience, howver, the counts cube can be +interpolation. For the user's convenience, however, the counts cube can be specified as ``incube`` parameter which will instruct ctexpcube to extract the exposure cube definition (such as sky coordinates and projection, number of pixels, pixel scale, energy binning) from the counts cube. Note, however, @@ -32,13 +32,16 @@ due to the point spread function and energy dispersion. By computing the exposure cube on the same grid as the counts cube, the spill over of events from sources at the edge of cube will not be handled correctly. +Note that there should at least be 25 bins per energy decade in the exposure +cube to assure that the energy dependence of the exposure is sufficiently +well sampled (in particular at low energies). + ctexpcube generates an exposure cube FITS file comprising three extensions. The primary extension contains a 3-dimensional image that contains the -exposure values. The next extension named ``EBOUNDS`` contains a binary table -that defines the energy boundaries of the exposure cube. The last extension -named ``GTI`` contains a binary table that provides the Good Time Intervals -of all observations that have been used for the computation of the exposure -cube. +exposure values. The next extension named ``ENERGIES`` contains a binary table +that defines the energies of the exposure cube. The last extension named ``GTI`` +contains a binary table that provides the Good Time Intervals of all +observations that have been used for the computation of the exposure cube. General parameters @@ -73,7 +76,10 @@ General parameters ``ebinfile [file]`` Name of the file containing the energy bin definition. - + +``(addbounds = yes) [boolean]`` + Add energies to the exposure cube at the observation energy boundaries? + ``(usepnt = no) [boolean]`` Use CTA pointing direction for cube centre instead of xref/yref parameters? diff --git a/doc/source/users/reference_manual/ctpsfcube.rst b/doc/source/users/reference_manual/ctpsfcube.rst index 14d30e6..929b371 100644 --- a/doc/source/users/reference_manual/ctpsfcube.rst +++ b/doc/source/users/reference_manual/ctpsfcube.rst @@ -29,13 +29,11 @@ recommended to use a rather coarse spatial binning to keep the file size manageable (with a typical value of ``binsz=1.0``). ctpsfcube generates a point spread function cube FITS file comprising three -extensions. The primary extension contains a 3-dimensional image that contains -the point spread function values. The energy and offset angle dimensions -of the point spread function cube are folded into the 3rd dimension of the -FITS image. The next extension named ``EBOUNDS`` contains a binary table -that defines the energy boundaries of the exposure cube. The last extension -named ``DELTAS`` contains a binary table that defines the offset angles -between true and measured arrival direction of the photon. +extensions. The primary extension contains a 4-dimensional image that contains +the point spread function values. The next extension named ``ENERGIES`` contains +a binary table that defines the energies of the point spread function cube. +The last extension named ``DELTAS`` contains a binary table that defines the +offset angles between true and measured arrival direction of the photon. General parameters @@ -71,6 +69,9 @@ General parameters ``ebinfile [file]`` Name of the file containing the energy bin definition. +``(addbounds = yes) [boolean]`` + Add energies to the point spread function cube at the observation energy boundaries? + ``(usepnt = no) [boolean]`` Use CTA pointing direction for cube centre instead of xref/yref parameters? diff --git a/src/ctbkgcube/ctbkgcube.cpp b/src/ctbkgcube/ctbkgcube.cpp index ad98167..6fd0cba 100644 --- a/src/ctbkgcube/ctbkgcube.cpp +++ b/src/ctbkgcube/ctbkgcube.cpp @@ -217,23 +217,46 @@ void ctbkgcube::run(void) log << std::endl; } + // Warn if there are not enough energy bins + int n = m_background.energies().size(); + double logEmin = std::log10(m_background.energies()[0].TeV()); + double logEmax = std::log10(m_background.energies()[n-1].TeV()); + int nrequired = int((logEmax - logEmin) * 25.0); + if ((n < nrequired) && logTerse()) { + log << std::endl; + log << "WARNING: "; + log << "Only " << n-1 << " energy bins have been requested. This may be " + "too few energy bins."; + log << std::endl << " "; + log << "At least 25 bins per decade in energy are recommended, which " + "in your case"; + log << std::endl << " "; + log << "would be " << nrequired-1 << " bins. Consider increasing the " + "number of energy bins."; + log << std::endl; + } + + // Write header + if (logTerse()) { + log << std::endl; + log.header1("Initialise background cube"); + } + + // Initialise exposure cube + init_cube(); + // Write observation(s) into logger if (logTerse()) { log << std::endl; - if (m_obs.size() > 1) { - log.header1("Observations"); - } - else { - log.header1("Observation"); - } - log << m_obs << std::endl; + log.header1(gammalib::number("Observation",m_obs.size())); + log << m_obs.print(m_chatter) << std::endl; } // Log input model if (logTerse()) { log << std::endl; log.header1("Input model"); - log << m_obs.models() << std::endl; + log << m_obs.models().print(m_chatter) << std::endl; } // Write header @@ -330,6 +353,11 @@ void ctbkgcube::run(void) // Append model to output container m_outmdl.append(model); + // Log background cube + if (logTerse()) { + log << m_background.print(m_chatter) << std::endl; + } + // Log output model if (logTerse()) { log << std::endl; @@ -445,10 +473,14 @@ void ctbkgcube::publish(const std::string& name) ***************************************************************************/ void ctbkgcube::init_members(void) { - // Initialise members + // Initialise user parameters m_outcube.clear(); m_outmodel.clear(); - m_publish = false; + m_addbounds = true; + m_publish = false; + m_chatter = static_cast(2); + + // Initialise protected members m_obs.clear(); m_background.clear(); m_bkgmdl.clear(); @@ -466,10 +498,14 @@ void ctbkgcube::init_members(void) ***************************************************************************/ void ctbkgcube::copy_members(const ctbkgcube& app) { - // Copy members - m_outmodel = app.m_outmodel; - m_outcube = app.m_outcube; - m_publish = app.m_publish; + // Copy user parameters + m_outmodel = app.m_outmodel; + m_outcube = app.m_outcube; + m_addbounds = app.m_addbounds; + m_publish = app.m_publish; + m_chatter = app.m_chatter; + + // Copy protected members m_obs = app.m_obs; m_background = app.m_background; m_bkgmdl = app.m_bkgmdl; @@ -558,7 +594,9 @@ void ctbkgcube::get_parameters(void) } // Get remaining parameters - m_publish = (*this)["publish"].boolean(); + m_addbounds = (*this)["addbounds"].boolean(); + m_publish = (*this)["publish"].boolean(); + m_chatter = static_cast((*this)["chatter"].integer()); // Read output filenames (if needed) if (read_ahead()) { @@ -569,3 +607,70 @@ void ctbkgcube::get_parameters(void) // Return return; } + + +/***********************************************************************//** + * @brief Initialise background cube + * + * Initialise the background cube. + ***************************************************************************/ +void ctbkgcube::init_cube(void) +{ + // Extract exposure cube definition + const GWcs* proj = static_cast(m_background.cube().projection()); + std::string wcs = m_background.cube().projection()->code(); + std::string coords = m_background.cube().projection()->coordsys(); + double x = proj->crval(0); + double y = proj->crval(1); + double dx = proj->cdelt(0); + double dy = proj->cdelt(1); + int nx = m_background.cube().nx(); + int ny = m_background.cube().ny(); + + // Extract energies + GEnergies energies = m_background.energies(); + + // If requested, insert energies at all event list energy boundaries + if (m_addbounds) { + + // Set logger + GLog* logger = NULL; + if (logTerse()) { + logger = &log; + } + + // Loop over all observations + for (int i = 0; i < m_obs.size(); ++i) { + + // Get observation and continue only if it is a CTA observation + const GCTAObservation* cta = dynamic_cast + (m_obs[i]); + + // Skip observation if it's not a CTA observation + if (cta == NULL) { + continue; + } + + // Skip observation if it does not contain an event list + if (cta->eventtype() != "EventList") { + continue; + } + + // Insert energy boundaries + energies = insert_energy_boundaries(energies, *cta, logger); + + } // endfor: looped over all observations + + } // endif: energy bin insertion requested + + // Setup background cube + m_background = GCTACubeBackground(wcs, coords, x, y, dx, dy, nx, ny, energies); + + // Log background cube + if (logTerse()) { + log << m_background.print(m_chatter) << std::endl; + } + + // Return + return; +}; diff --git a/src/ctbkgcube/ctbkgcube.hpp b/src/ctbkgcube/ctbkgcube.hpp index 8de0f5b..8e6c095 100644 --- a/src/ctbkgcube/ctbkgcube.hpp +++ b/src/ctbkgcube/ctbkgcube.hpp @@ -69,12 +69,15 @@ protected: void copy_members(const ctbkgcube& app); void free_members(void); void get_parameters(void); - void fill_cube(GCTAObservation* obs); + void init_cube(void); + //void fill_cube(GCTAObservation* obs); // User parameters GFilename m_outcube; //!< Filename of output cube GFilename m_outmodel; //!< Filename of output XML model + bool m_addbounds; //!< Add energies at boundaries? bool m_publish; //!< Publish background cube? + GChatter m_chatter; //!< Chattiness // Protected members GObservations m_obs; //!< Observation container diff --git a/src/ctbkgcube/ctbkgcube.par b/src/ctbkgcube/ctbkgcube.par index 0de9d88..82edd81 100644 --- a/src/ctbkgcube/ctbkgcube.par +++ b/src/ctbkgcube/ctbkgcube.par @@ -28,11 +28,12 @@ outmodel, f, a, "NONE",,, "Output model XML file" # # Energy binning parameters #========================== -ebinalg, s, h, "LOG",FILE|LIN|LOG,, "Algorithm for defining energy bins" -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" -ebinfile, f, h, "NONE",,, "Name of the file containing the energy bin definition" +ebinalg, s, h, "LOG",FILE|LIN|LOG,, "Algorithm for defining energy bins" +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" +ebinfile, f, h, "NONE",,, "Name of the file containing the energy bin definition" +addbounds, b, h, yes,,, "Add energies at observation energy boundaries?" # # Spatial binning parameters diff --git a/src/ctedispcube/ctedispcube.cpp b/src/ctedispcube/ctedispcube.cpp index 3bc2c8c..728ecb0 100644 --- a/src/ctedispcube/ctedispcube.cpp +++ b/src/ctedispcube/ctedispcube.cpp @@ -229,16 +229,20 @@ void ctedispcube::run(void) log << std::endl; } + // Write header + if (logTerse()) { + log << std::endl; + log.header1("Initialise energy dispersion cube"); + } + + // Initialise energy dispersion cube + init_cube(); + // Write observation(s) into logger if (logTerse()) { log << std::endl; - if (m_obs.size() > 1) { - log.header1("Observations"); - } - else { - log.header1("Observation"); - } - log << m_obs << std::endl; + log.header1(gammalib::number("Observation",m_obs.size())); + log << m_obs.print(m_chatter) << std::endl; } // Write header @@ -247,9 +251,14 @@ void ctedispcube::run(void) log.header1("Generate energy dispersion cube"); } - // Fill Edisp + // Fill energy dispersion cube m_edispcube.fill(m_obs, &log); + // Log energy dispersion cube + if (logTerse()) { + log << m_edispcube.print(m_chatter) << std::endl; + } + // Return return; } @@ -302,8 +311,12 @@ void ctedispcube::save(void) ***************************************************************************/ void ctedispcube::init_members(void) { - // Initialise members + // Initialise user parameters m_outcube.clear(); + m_addbounds = true; + m_chatter = static_cast(2); + + // Initialise protected members m_obs.clear(); m_edispcube.clear(); @@ -319,8 +332,12 @@ void ctedispcube::init_members(void) ***************************************************************************/ void ctedispcube::copy_members(const ctedispcube& app) { - // Copy members + // Copy user parameters m_outcube = app.m_outcube; + m_addbounds = app.m_addbounds; + m_chatter = app.m_chatter; + + // Copy protected members m_obs = app.m_obs; m_edispcube = app.m_edispcube; @@ -380,7 +397,7 @@ void ctedispcube::get_parameters(void) } - // ... otherwise setup the exposure cube from the counts map + // ... otherwise setup the energy disperison cube from the counts map else { // Load event cube from filename @@ -391,6 +408,10 @@ void ctedispcube::get_parameters(void) } // endelse: cube loaded from file + // Get remaining parameters + m_addbounds = (*this)["addbounds"].boolean(); + m_chatter = static_cast((*this)["chatter"].integer()); + // Read output filename (if needed) if (read_ahead()) { m_outcube = (*this)["outcube"].filename(); @@ -399,3 +420,73 @@ void ctedispcube::get_parameters(void) // Return return; } + + +/***********************************************************************//** + * @brief Initialise energy dispersion cube + * + * Initialise the energy dispersion cube. + ***************************************************************************/ +void ctedispcube::init_cube(void) +{ + // Extract energy dispersion definition + const GWcs* proj = static_cast(m_edispcube.cube().projection()); + std::string wcs = m_edispcube.cube().projection()->code(); + std::string coords = m_edispcube.cube().projection()->coordsys(); + double x = proj->crval(0); + double y = proj->crval(1); + double dx = proj->cdelt(0); + double dy = proj->cdelt(1); + int nx = m_edispcube.cube().nx(); + int ny = m_edispcube.cube().ny(); + double mmax = (*this)["migramax"].real(); + int nmbins = (*this)["migrabins"].integer(); + + // Extract energies + GEnergies energies = m_edispcube.energies(); + + // If requested, insert energies at all event list energy boundaries + if (m_addbounds) { + + // Set logger + GLog* logger = NULL; + if (logTerse()) { + logger = &log; + } + + // Loop over all observations + for (int i = 0; i < m_obs.size(); ++i) { + + // Get observation and continue only if it is a CTA observation + const GCTAObservation* cta = dynamic_cast + (m_obs[i]); + + // Skip observation if it's not a CTA observation + if (cta == NULL) { + continue; + } + + // Skip observation if it does not contain an event list + if (cta->eventtype() != "EventList") { + continue; + } + + // Insert energy boundaries + energies = insert_energy_boundaries(energies, *cta, logger); + + } // endfor: looped over all observations + + } // endif: energy bin insertion requested + + // Setup energy dispersion cube + m_edispcube = GCTACubeEdisp(wcs, coords, x, y, dx, dy, nx, ny, energies, + mmax, nmbins); + + // Log energy dispersion cube + if (logTerse()) { + log << m_edispcube.print(m_chatter) << std::endl; + } + + // Return + return; +}; diff --git a/src/ctedispcube/ctedispcube.hpp b/src/ctedispcube/ctedispcube.hpp index 62f206f..fce7c9f 100644 --- a/src/ctedispcube/ctedispcube.hpp +++ b/src/ctedispcube/ctedispcube.hpp @@ -67,9 +67,12 @@ protected: void copy_members(const ctedispcube& app); void free_members(void); void get_parameters(void); + void init_cube(void); // User parameters - GFilename m_outcube; //!< Output energy dispersion cube file name + GFilename m_outcube; //!< Output exposure cube file + bool m_addbounds; //!< Add energies at boundaries? + GChatter m_chatter; //!< Chattiness // Protected members GObservations m_obs; //!< Observation container diff --git a/src/ctedispcube/ctedispcube.par b/src/ctedispcube/ctedispcube.par index 21bb9b3..148c16c 100644 --- a/src/ctedispcube/ctedispcube.par +++ b/src/ctedispcube/ctedispcube.par @@ -26,11 +26,12 @@ outcube, f, a, "edispcube.fits",,, "Output energy dispersion cube file" # # Energy binning parameters #========================== -ebinalg, s, h, "LOG", FILE|LIN|LOG,, "Algorithm for defining energy bins" -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" -ebinfile, f, h, "NONE",,, "Name of the file containing the energy bin definition" +ebinalg, s, h, "LOG", FILE|LIN|LOG,, "Algorithm for defining energy bins" +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" +ebinfile, f, h, "NONE",,, "Name of the file containing the energy bin definition" +addbounds, b, h, yes,,, "Add energies at observation energy boundaries?" # # Spatial binning parameters diff --git a/src/ctexpcube/ctexpcube.cpp b/src/ctexpcube/ctexpcube.cpp index adda25d..1b88f0e 100644 --- a/src/ctexpcube/ctexpcube.cpp +++ b/src/ctexpcube/ctexpcube.cpp @@ -217,11 +217,39 @@ void ctexpcube::run(void) log << std::endl; } + // Warn if there are not enough energy bins + int n = m_expcube.energies().size(); + double logEmin = std::log10(m_expcube.energies()[0].TeV()); + double logEmax = std::log10(m_expcube.energies()[n-1].TeV()); + int nrequired = int((logEmax - logEmin) * 25.0); + if ((n < nrequired) && logTerse()) { + log << std::endl; + log << "WARNING: "; + log << "Only " << n-1 << " energy bins have been requested. This may be " + "too few energy bins."; + log << std::endl << " "; + log << "At least 25 bins per decade in energy are recommended, which " + "in your case"; + log << std::endl << " "; + log << "would be " << nrequired-1 << " bins. Consider increasing the " + "number of energy bins."; + log << std::endl; + } + + // Write header + if (logTerse()) { + log << std::endl; + log.header1("Initialise exposure cube"); + } + + // Initialise exposure cube + init_cube(); + // Write observation(s) into logger if (logTerse()) { log << std::endl; log.header1(gammalib::number("Observation",m_obs.size())); - log << m_obs << std::endl; + log << m_obs.print(m_chatter) << std::endl; } // Write header @@ -230,9 +258,14 @@ void ctexpcube::run(void) log.header1("Generate exposure cube"); } - // Fill exposure + // Fill exposure cube m_expcube.fill(m_obs, &log); + // Log exposure cube + if (logTerse()) { + log << m_expcube.print(m_chatter) << std::endl; + } + // Optionally publish exposure cube if (m_publish) { publish(); @@ -321,9 +354,11 @@ void ctexpcube::publish(const std::string& name) ***************************************************************************/ void ctexpcube::init_members(void) { - // Initialise members + // Initialise user parameters m_outcube.clear(); - m_publish = false; + m_addbounds = true; + m_publish = false; + m_chatter = static_cast(2); // Initialise protected members m_obs.clear(); @@ -341,9 +376,11 @@ void ctexpcube::init_members(void) ***************************************************************************/ void ctexpcube::copy_members(const ctexpcube& app) { - // Copy attributes - m_outcube = app.m_outcube; - m_publish = app.m_publish; + // Copy user parameters + m_outcube = app.m_outcube; + m_addbounds = app.m_addbounds; + m_publish = app.m_publish; + m_chatter = app.m_chatter; // Copy protected members m_obs = app.m_obs; @@ -413,7 +450,9 @@ void ctexpcube::get_parameters(void) } // endelse: cube loaded from file // Get remaining parameters - m_publish = (*this)["publish"].boolean(); + m_addbounds = (*this)["addbounds"].boolean(); + m_publish = (*this)["publish"].boolean(); + m_chatter = static_cast((*this)["chatter"].integer()); // Read output filename (if needed) if (read_ahead()) { @@ -423,3 +462,70 @@ void ctexpcube::get_parameters(void) // Return return; } + + +/***********************************************************************//** + * @brief Initialise exposure cube + * + * Initialise the exposure cube. + ***************************************************************************/ +void ctexpcube::init_cube(void) +{ + // Extract exposure cube definition + const GWcs* proj = static_cast(m_expcube.cube().projection()); + std::string wcs = m_expcube.cube().projection()->code(); + std::string coords = m_expcube.cube().projection()->coordsys(); + double x = proj->crval(0); + double y = proj->crval(1); + double dx = proj->cdelt(0); + double dy = proj->cdelt(1); + int nx = m_expcube.cube().nx(); + int ny = m_expcube.cube().ny(); + + // Extract energies + GEnergies energies = m_expcube.energies(); + + // If requested, insert energies at all event list energy boundaries + if (m_addbounds) { + + // Set logger + GLog* logger = NULL; + if (logTerse()) { + logger = &log; + } + + // Loop over all observations + for (int i = 0; i < m_obs.size(); ++i) { + + // Get observation and continue only if it is a CTA observation + const GCTAObservation* cta = dynamic_cast + (m_obs[i]); + + // Skip observation if it's not a CTA observation + if (cta == NULL) { + continue; + } + + // Skip observation if it does not contain an event list + if (cta->eventtype() != "EventList") { + continue; + } + + // Insert energy boundaries + energies = insert_energy_boundaries(energies, *cta, logger); + + } // endfor: looped over all observations + + } // endif: energy bin insertion requested + + // Setup exposure cube + m_expcube = GCTACubeExposure(wcs, coords, x, y, dx, dy, nx, ny, energies); + + // Log exposure cube + if (logTerse()) { + log << m_expcube.print(m_chatter) << std::endl; + } + + // Return + return; +}; diff --git a/src/ctexpcube/ctexpcube.hpp b/src/ctexpcube/ctexpcube.hpp index beac456..108c3c6 100644 --- a/src/ctexpcube/ctexpcube.hpp +++ b/src/ctexpcube/ctexpcube.hpp @@ -68,10 +68,13 @@ protected: void copy_members(const ctexpcube& app); void free_members(void); void get_parameters(void); + void init_cube(void); // User parameters - GFilename m_outcube; //!< Output exposure cube file - bool m_publish; //!< Publish exposure cube? + GFilename m_outcube; //!< Output exposure cube file + bool m_addbounds; //!< Add energies at boundaries? + bool m_publish; //!< Publish exposure cube? + GChatter m_chatter; //!< Chattiness // Protected members GObservations m_obs; //!< Observation container diff --git a/src/ctexpcube/ctexpcube.par b/src/ctexpcube/ctexpcube.par index 9cbb8f7..594d3b8 100644 --- a/src/ctexpcube/ctexpcube.par +++ b/src/ctexpcube/ctexpcube.par @@ -26,11 +26,12 @@ outcube, f, a, "expcube.fits",,, "Output exposure cube file" # # Energy binning parameters #========================== -ebinalg, s, h, "LOG", FILE|LIN|LOG,,"Algorithm for defining energy bins" -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" -ebinfile, f, h, "NONE",,, "Name of the file containing the energy bin definition" +ebinalg, s, h, "LOG", FILE|LIN|LOG,,"Algorithm for defining energy bins" +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" +ebinfile, f, h, "NONE",,, "Name of the file containing the energy bin definition" +addbounds, b, h, yes,,, "Add energies at observation energy boundaries?" # # Spatial binning parameters diff --git a/src/ctpsfcube/ctpsfcube.cpp b/src/ctpsfcube/ctpsfcube.cpp index 6d7b3a7..3e0c664 100644 --- a/src/ctpsfcube/ctpsfcube.cpp +++ b/src/ctpsfcube/ctpsfcube.cpp @@ -216,16 +216,20 @@ void ctpsfcube::run(void) log << std::endl; } + // Write header + if (logTerse()) { + log << std::endl; + log.header1("Initialise PSF cube"); + } + + // Initialise exposure cube + init_cube(); + // Write observation(s) into logger if (logTerse()) { log << std::endl; - if (m_obs.size() > 1) { - log.header1("Observations"); - } - else { - log.header1("Observation"); - } - log << m_obs << std::endl; + log.header1(gammalib::number("Observation",m_obs.size())); + log << m_obs.print(m_chatter) << std::endl; } // Write header @@ -234,9 +238,14 @@ void ctpsfcube::run(void) log.header1("Generate PSF cube"); } - // Fill PSF + // Fill PSF cube m_psfcube.fill(m_obs, &log); + // Log PSF cube + if (logTerse()) { + log << m_psfcube.print(m_chatter) << std::endl; + } + // Return return; } @@ -286,8 +295,10 @@ void ctpsfcube::save(void) ***************************************************************************/ void ctpsfcube::init_members(void) { - // Initialise members + // Initialise user parameters m_outcube.clear(); + m_addbounds = true; + m_chatter = static_cast(2); // Initialise protected members m_obs.clear(); @@ -305,12 +316,14 @@ void ctpsfcube::init_members(void) ***************************************************************************/ void ctpsfcube::copy_members(const ctpsfcube& app) { - // Copy attributes - m_outcube = app.m_outcube; + // Copy user parameters + m_outcube = app.m_outcube; + m_addbounds = app.m_addbounds; + m_chatter = app.m_chatter; // Copy protected members - m_obs = app.m_obs; - m_psfcube = app.m_psfcube; + m_obs = app.m_obs; + m_psfcube = app.m_psfcube; // Return return; @@ -379,6 +392,10 @@ void ctpsfcube::get_parameters(void) } // endelse: cube loaded from file + // Get remaining parameters + m_addbounds = (*this)["addbounds"].boolean(); + m_chatter = static_cast((*this)["chatter"].integer()); + // Read output filename (if needed) if (read_ahead()) { m_outcube = (*this)["outcube"].filename(); @@ -387,3 +404,73 @@ void ctpsfcube::get_parameters(void) // Return return; } + + +/***********************************************************************//** + * @brief Initialise PSF cube + * + * Initialise the PSF cube. + ***************************************************************************/ +void ctpsfcube::init_cube(void) +{ + // Extract exposure cube definition + const GWcs* proj = static_cast(m_psfcube.cube().projection()); + std::string wcs = m_psfcube.cube().projection()->code(); + std::string coords = m_psfcube.cube().projection()->coordsys(); + double x = proj->crval(0); + double y = proj->crval(1); + double dx = proj->cdelt(0); + double dy = proj->cdelt(1); + int nx = m_psfcube.cube().nx(); + int ny = m_psfcube.cube().ny(); + double dmax = (*this)["amax"].real(); + int ndbins = (*this)["anumbins"].integer(); + + // Extract energies + GEnergies energies = m_psfcube.energies(); + + // If requested, insert energies at all event list energy boundaries + if (m_addbounds) { + + // Set logger + GLog* logger = NULL; + if (logTerse()) { + logger = &log; + } + + // Loop over all observations + for (int i = 0; i < m_obs.size(); ++i) { + + // Get observation and continue only if it is a CTA observation + const GCTAObservation* cta = dynamic_cast + (m_obs[i]); + + // Skip observation if it's not a CTA observation + if (cta == NULL) { + continue; + } + + // Skip observation if it does not contain an event list + if (cta->eventtype() != "EventList") { + continue; + } + + // Insert energy boundaries + energies = insert_energy_boundaries(energies, *cta, logger); + + } // endfor: looped over all observations + + } // endif: energy bin insertion requested + + // Setup PSF cube + m_psfcube = GCTACubePsf(wcs, coords, x, y, dx, dy, nx, ny, energies, + dmax, ndbins); + + // Log PSF cube + if (logTerse()) { + log << m_psfcube.print(m_chatter) << std::endl; + } + + // Return + return; +}; diff --git a/src/ctpsfcube/ctpsfcube.hpp b/src/ctpsfcube/ctpsfcube.hpp index 7b469b3..3079e0f 100644 --- a/src/ctpsfcube/ctpsfcube.hpp +++ b/src/ctpsfcube/ctpsfcube.hpp @@ -67,9 +67,12 @@ protected: void copy_members(const ctpsfcube& app); void free_members(void); void get_parameters(void); + void init_cube(void); // User parameters - GFilename m_outcube; //!< Output PSF cube file name + GFilename m_outcube; //!< Output PSF cube file name + bool m_addbounds; //!< Add energies at boundaries? + GChatter m_chatter; //!< Chattiness // Protected members GObservations m_obs; //!< Observation container diff --git a/src/ctpsfcube/ctpsfcube.par b/src/ctpsfcube/ctpsfcube.par index b9f99eb..ee0ad01 100644 --- a/src/ctpsfcube/ctpsfcube.par +++ b/src/ctpsfcube/ctpsfcube.par @@ -5,7 +5,7 @@ # ------------------------------------------------------------------------- # # # # File: ctpsfcube.par # -# Version: 1.0.0 # +# Version: 1.1.0 # # # # Author: Chia-Chun Lu # # c.lu@mpi-hd.mpg.de # @@ -26,11 +26,12 @@ outcube, f, a, "psfcube.fits",,, "Output PSF cube file" # # Energy binning parameters #========================== -ebinalg, s, h, "LOG", FILE|LIN|LOG,, "Algorithm for defining energy bins" -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" -ebinfile, f, h, "NONE",,, "Name of the file containing the energy bin definition" +ebinalg, s, h, "LOG", FILE|LIN|LOG,, "Algorithm for defining energy bins" +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" +ebinfile, f, h, "NONE",,, "Name of the file containing the energy bin definition" +addbounds, b, h, yes,,, "Add energies at observation energy boundaries?" # # Spatial binning parameters diff --git a/src/support/ctool.cpp b/src/support/ctool.cpp index c8b8cde..d5362dd 100644 --- a/src/support/ctool.cpp +++ b/src/support/ctool.cpp @@ -1303,6 +1303,81 @@ std::string ctool::get_obs_header(const GObservation* obs) /***********************************************************************//** + * @brief Insert observation energy boundaries into list of energies + * + * @param[in] energies Energies. + * @param[in] obs Observation container. + * @return Energies. + ***************************************************************************/ +GEnergies ctool::insert_energy_boundaries(const GEnergies& energies, + const GCTAObservation& obs, + GLog* log) const +{ + // Create copy of input energies + GEnergies engs = energies; + + // Get the energy boundaries of the event list + GEbounds ebounds = obs.ebounds(); + + // Loop over all boundary energies + for (int iebin = 0; iebin <= ebounds.size(); ++iebin) { + + // Get boundary energy + GEnergy energy; + if (iebin < ebounds.size()) { + energy = ebounds.emin(iebin); + } + else { + energy = ebounds.emax(); + } + + // Get index before which the energy should be appended + bool insert = true; + int index = -1; + for (int k = 0; k < engs.size(); ++k) { + + // If energy exists already then skip the boundary and examine the + // next one. We consider here 1 MeV as being sufficiently close. + if (std::abs(engs[k].MeV()-energy.MeV()) < 1.0) { + insert = false; + break; + } + + // If energy is above the boundary energy we found the index + if (engs[k] > energy) { + index = k; + break; + } + + } // endfor: search index for insertion + + // Insert energy if requested + if (insert) { + + // Insert or append energy + if (index != -1) { + engs.insert(index, energy); + } + else { + engs.append(energy); + } + + // Log energy insertion + if (log != NULL) { + *log << gammalib::parformat("Insert energy"); + *log << energy.print() << std::endl; + } + + } // endif: energy insertion requested + + } // endfor: looped over all energy boundaries + + // Return energies + return engs; +} + + +/***********************************************************************//** * @brief Dumps help text in the console ***************************************************************************/ void ctool::provide_help(void) const diff --git a/src/support/ctool.hpp b/src/support/ctool.hpp index dfb92ee..0bb2e86 100644 --- a/src/support/ctool.hpp +++ b/src/support/ctool.hpp @@ -92,6 +92,9 @@ protected: GSkyDir get_mean_pointing(const GObservations& obs); size_t get_current_rss(void); std::string get_obs_header(const GObservation* obs); + GEnergies insert_energy_boundaries(const GEnergies& energies, + const GCTAObservation& obs, + GLog* log = NULL) const; // Protected members bool m_read_ahead; //!< Read ahead output parameters -- libgit2 0.21.4