Add "addbounds" parameter to response cube tools

Jürgen Knödlseder
1 parent e8719183
Showing 19 changed files with 607 additions and 101 deletions
cscripts/cspull.py
doc/source/users/reference_manual/ctbkgcube.rst
doc/source/users/reference_manual/ctedispcube.rst
doc/source/users/reference_manual/ctexpcube.rst
doc/source/users/reference_manual/ctpsfcube.rst
src/ctbkgcube/ctbkgcube.cpp
src/ctbkgcube/ctbkgcube.hpp
src/ctbkgcube/ctbkgcube.par
src/ctedispcube/ctedispcube.cpp
src/ctedispcube/ctedispcube.hpp
src/ctedispcube/ctedispcube.par
src/ctexpcube/ctexpcube.cpp
src/ctexpcube/ctexpcube.hpp
src/ctexpcube/ctexpcube.par
src/ctpsfcube/ctpsfcube.cpp
src/ctpsfcube/ctpsfcube.hpp
src/ctpsfcube/ctpsfcube.par
src/support/ctool.cpp
src/support/ctool.hpp
@@ -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
@@ -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?
@@ -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?
@@ -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?
@@ -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?
@@ -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&amp; 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<GChatter>(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<GChatter>((*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<const GWcs*>(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<const GCTAObservation*>
+                                         (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;
+};
@@ -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
@@ -28,11 +28,12 @@ outmodel, f, a, &quot;NONE&quot;,,, &quot;Output model XML file&quot;
 #
 # 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
@@ -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<GChatter>(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<GChatter>((*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<const GWcs*>(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<const GCTAObservation*>
+                                         (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;
+};
@@ -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
@@ -26,11 +26,12 @@ outcube,  f, a, &quot;edispcube.fits&quot;,,, &quot;Output energy dispersion cube file&quot;
 #
 # 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
@@ -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&amp; 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<GChatter>(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<GChatter>((*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<const GWcs*>(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<const GCTAObservation*>
+                                         (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;
+};
@@ -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
@@ -26,11 +26,12 @@ outcube,  f, a, &quot;expcube.fits&quot;,,, &quot;Output exposure cube file&quot;
 #
 # 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
@@ -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<GChatter>(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<GChatter>((*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<const GWcs*>(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<const GCTAObservation*>
+                                         (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;
+};
@@ -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
@@ -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, &quot;psfcube.fits&quot;,,, &quot;Output PSF cube file&quot;
 #
 # 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
@@ -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
@@ -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