Merge branch '1330-ctmodel-multiple-handling' into add-base-class

Jürgen Knödlseder
2 parents c600f360 28a6964e
Showing 8 changed files with 456 additions and 449 deletions
ChangeLog
NEWS
pyext/ctmodel.i
src/ctmodel/ctmodel.cpp
src/ctmodel/ctmodel.hpp
src/ctmodel/ctmodel.par
test/test_ctmodel.py
test/test_ctools.sh
@@ -6,6 +6,7 @@
         Clean-up public interface
 	Adapt to GammaLib 0.9.0 interface
         ctbin now always merges event lists into a single counts cube
+        ctmodel now always creates a single model cube
         Add ctexpcube tool (generates exposure cube)
         Add ctpsfcube tool (generates PSF cube)
         Add ctbkgcube tool (generates background cube)
@@ -98,8 +98,26 @@ extension in the FITS file.
 None
  
  
-10. ctmodel - CTA counts map modelling
---------------------------------------
+10. ctmodel - CTA counts cube modelling
+---------------------------------------
+WARNING: The philosophy behind the ctmodel tool changed fundamentally!
+
+ctmodel now always produces a single model cube on output. There are
+two important input parameters: "infile" and "obsfile". "infile"
+specifies a counts or model cube that is used to define the binning
+of the model cube. The functionnality for this parameter is as before.
+If "infile=NONE", the cube will be defined according to task parameters
+that specify the spatial and the energy binning. The "obsfile"
+parameter is new and allows specifying the observation definition
+for which the model cube should be computed. If an event list or a
+counts cube is presented, the pointing and Good Time Intervals will
+be extracted from the FITS file header. If an observation definition
+XML file is presented, ctmodel will now loop over all CTA observations
+in the file and add the model values to the model cube. If no
+response information is associated to the "obsfile" (which is typically
+the case for event lists or counts cubes), the "caldb" and "irf"
+parameters will be used to identify the response.
+
 ctmodel now has full access to the calibration database, allowing
 the use for 2D response functions in FITS files.
  
 /***************************************************************************
- *                     ctmodel - CTA counts model tool                     *
+ *                ctmodel - CTA model cube generation tool                 *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2012-2014 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
@@ -20,7 +20,7 @@
  ***************************************************************************/
 /**
  * @file ctmodel.i
- * @brief CTA counts model tool definition
+ * @brief CTA model cube tool definition
  * @author Juergen Knoedlseder
  */
 %{
@@ -51,6 +51,7 @@ public:
     void                 run(void);
     void                 save(void);
     const GObservations& obs(void) const;
+    const GCTAEventCube& cube(void) const;
 };
  
  
 /***************************************************************************
- *                     ctmodel - CTA counts model tool                     *
+ *                ctmodel - CTA model cube generation tool                 *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2012-2014 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
@@ -20,7 +20,7 @@
  ***************************************************************************/
 /**
  * @file ctmodel.cpp
- * @brief CTA counts model tool implementation
+ * @brief CTA model cube tool implementation
  * @author Juergen Knoedlseder
  */
  
@@ -33,8 +33,7 @@
 #include "GTools.hpp"
  
 /* __ Method name definitions ____________________________________________ */
-#define G_SETUP_OBS                                    "ctmodel::setup_obs()"
-#define G_MODEL_MAP                    "ctmodel::model_map(GCTAObservation*)"
+#define G_GET_EBOUNDS                                "ctmodel::get_ebounds()"
  
 /* __ Debug definitions __________________________________________________ */
  
@@ -212,10 +211,10 @@ void ctmodel::execute(void)
     // Signal that some parameters should be read ahead
     m_read_ahead = true;
  
-    // Create the model map(s)
+    // Create model cube
     run();
  
-    // Save the model map(s) into FITS file
+    // Save model cube into FITS file
     save();
  
     // Return
@@ -276,23 +275,15 @@ void ctmodel::run(void)
     // Write header
     if (logTerse()) {
         log << std::endl;
-        if (m_obs.size() > 1) {
-            log.header1("Generate model maps");
-        }
-        else {
-            log.header1("Generate model map");
-        }
+        log.header1("Generate model cube");
     }
  
-    // Initialise observation counter
-    int n_observations = 0;
+    // Initialise counts cube
+    init_cube();
  
     // Loop over all observations in the container
     for (int i = 0; i < m_obs.size(); ++i) {
  
-        // Initialise event input and output filenames
-        m_infiles.push_back("");
-
         // Get CTA observation
         GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
  
@@ -309,35 +300,18 @@ void ctmodel::run(void)
                 }
             }
  
-            // Increment number of observations
-            n_observations++;
-
-            // Save event file name (for possible saving)
-            m_infiles[i] = obs->eventfile();
-
-            // Generate model map
-            model_map(obs, m_obs.models());
+            // Fill the cube
+            fill_cube(obs);
  
         } // endif: CTA observation found
  
     } // endfor: looped over observations
  
-    // If more than a single observation has been handled then make sure
-    // that an XML file will be used for storage
-    if (n_observations > 1) {
-        m_use_xml = true;
-    }
-
-    // Write observation(s) into logger
+    // Log cube
     if (logTerse()) {
         log << std::endl;
-        if (m_obs.size() > 1) {
-            log.header1("Observations after model map generation");
-        }
-        else {
-            log.header1("Observation after model map generation");
-        }
-        log << m_obs << std::endl;
+        log.header1("Model cube");
+        log << m_cube << std::endl;
     }
  
     // Restore energy dispersion flag for all CTA observations
@@ -354,43 +328,135 @@ void ctmodel::run(void)
  
  
 /***********************************************************************//**
- * @brief Save model map(s)
- *
- * This method saves the model map(s) into FITS file(s). There are two
- * modes, depending on the m_use_xml flag.
+ * @brief Save model cube
  *
- * If m_use_xml is true, all model map(s) will be saved into FITS files,
- * where the output filenames are constructued from the input filenames by
- * prepending the m_prefix string to name. Any path information will be
- * stripped form the input name, hence event files will be written into the
- * local working directory (unless some path information is present
- * in the prefix). In addition, an XML file will be created that gathers
- * the filename information for the model map(s). If an XML file was present
- * on input, all metadata information will be copied from this input file.
- *
- * If m_use_xml is false, the model map will be saved into a FITS file.
+ * Saves the model cube into a FITS file specified using the "outfile"
+ * task parameter.
  ***************************************************************************/
 void ctmodel::save(void)
 {
     // Write header
     if (logTerse()) {
         log << std::endl;
-        if (m_obs.size() > 1) {
-            log.header1("Save observations");
-        }
-        else {
-            log.header1("Save observation");
-        }
+        log.header1("Save cube");
     }
  
-    // Case A: Save model map(s) and XML metadata information
-    if (m_use_xml) {
-        save_xml();
-    }
+    // Make sure we have the FITS filename
+    m_outfile = (*this)["outfile"].filename();
  
-    // Case B: Save model map as FITS file
-    else {
-        save_fits();
+    // Save model cube into FITS file
+    m_cube.save(m_outfile, clobber());
+
+    // Return
+    return;
+}
+
+
+/*==========================================================================
+ =                                                                         =
+ =                             Private methods                             =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Initialise class members
+ ***************************************************************************/
+void ctmodel::init_members(void)
+{
+    // Initialise members
+    m_infile.clear();
+    m_obsfile.clear();
+    m_outfile.clear();
+    m_caldb.clear();
+    m_irf.clear();
+    m_srcmdl.clear();
+    m_proj.clear();
+    m_coordsys.clear();
+    m_ebinalg.clear();
+    m_ebinfile.clear();
+    m_ra          = 0.0;
+    m_dec         = 0.0;
+    m_deadc       = 1.0;
+    m_tmin        = 0.0;
+    m_tmax        = 0.0;
+    m_emin        = 0.0;
+    m_emax        = 0.0;
+    m_enumbins    = 0;
+    m_xref        = 0.0;
+    m_yref        = 0.0;
+    m_binsz       = 0.0;
+    m_nxpix       = 0;
+    m_nypix       = 0;
+    m_apply_edisp = false;
+
+    // Initialise protected members
+    m_obs.clear();
+    m_cube.clear();
+    m_ebounds.clear();
+    m_gti.clear();
+    m_read_ahead = false;
+
+    // Set logger properties
+    log.date(true);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Copy class members
+ *
+ * @param[in] app Application.
+ ***************************************************************************/
+void ctmodel::copy_members(const ctmodel& app)
+{
+    // Copy attributes
+    m_infile      = app.m_infile;
+    m_obsfile     = app.m_obsfile;
+    m_outfile     = app.m_outfile;
+    m_caldb       = app.m_caldb;
+    m_irf         = app.m_irf;
+    m_srcmdl      = app.m_srcmdl;
+    m_proj        = app.m_proj;
+    m_coordsys    = app.m_coordsys;
+    m_ebinalg     = app.m_ebinalg;
+    m_ebinfile    = app.m_ebinfile;
+    m_ra          = app.m_ra;
+    m_dec         = app.m_dec;
+    m_deadc       = app.m_deadc;
+    m_tmin        = app.m_tmin;
+    m_tmax        = app.m_tmax;
+    m_emin        = app.m_emin;
+    m_emax        = app.m_emax;
+    m_enumbins    = app.m_enumbins;
+    m_xref        = app.m_xref;
+    m_yref        = app.m_yref;
+    m_binsz       = app.m_binsz;
+    m_nxpix       = app.m_nxpix;
+    m_nypix       = app.m_nypix;
+    m_apply_edisp = app.m_apply_edisp;
+
+    // Copy protected members
+    m_obs        = app.m_obs;
+    m_cube       = app.m_cube;
+    m_ebounds    = app.m_ebounds;
+    m_gti        = app.m_gti;
+    m_read_ahead = app.m_read_ahead;
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Delete class members
+ ***************************************************************************/
+void ctmodel::free_members(void)
+{
+    // Write separator into logger
+    if (logTerse()) {
+        log << std::endl;
     }
  
     // Return
@@ -406,62 +472,77 @@ void ctmodel::save(void)
  ***************************************************************************/
 void ctmodel::get_parameters(void)
 {
-    // If we do not have any observations in the container then get all
-    // information to build (at least) one
-    if (m_obs.size() == 0) {
+    // Read input and optionally output cube filenames
+    m_infile = (*this)["infile"].filename();
+    if (m_read_ahead) {
+        m_outfile = (*this)["outfile"].filename();
+    }
  
-        // Read general parameters. Only read the output parameters if
-        // required, and also only read the model XML filename if no
-        // model is yet in the container.
-        m_infile = (*this)["infile"].filename();
-        if (m_read_ahead) {
-            m_outfile = (*this)["outfile"].filename();
-            m_prefix  = (*this)["prefix"].string();
-        }
-        m_caldb  = (*this)["caldb"].string();
-        m_irf    = (*this)["irf"].string();
-        if (m_obs.models().size() == 0) {
-            m_srcmdl = (*this)["srcmdl"].filename();
+    // Read model definiton file if required
+    if (m_obs.models().size() == 0) {
+        m_srcmdl = (*this)["srcmdl"].filename();
+    }
+
+    // Get energy dispersion flag parameters
+    m_apply_edisp = (*this)["edisp"].boolean();
+
+    // If the observation container is empty then read the parameters that
+    // are necessary to fill it
+    if (m_obs.size() == 0) {
+    
+        // Read observation definition filename and response parameters
+        m_obsfile = (*this)["obsfile"].filename();
+        m_caldb   = (*this)["caldb"].string();
+        m_irf     = (*this)["irf"].string();
+
+        // If no observation definition file has been specified then read all
+        // parameters that are necessary to create an observation from scratch
+        // (see method setup_obs)
+        if ((m_obsfile == "NONE") || (gammalib::strip_whitespace(m_obsfile) == "")) {
+            m_ra    = (*this)["ra"].real();
+            m_dec   = (*this)["dec"].real();
+            m_deadc = (*this)["deadc"].real();
+            m_tmin  = (*this)["tmin"].real();
+            m_tmax  = (*this)["tmax"].real();
+            m_emin  = (*this)["emin"].real();
+            m_emax  = (*this)["emax"].real();
         }
  
-        // If there is no input filename then read all parameters that
-        // are required to build a model map from scratch
-        if ((m_infile == "NONE") || (gammalib::strip_whitespace(m_infile) == "")) {
-            m_ra       = (*this)["ra"].real();
-            m_dec      = (*this)["dec"].real();
-            m_deadc    = (*this)["deadc"].real();
-            m_tmin     = (*this)["tmin"].real();
-            m_tmax     = (*this)["tmax"].real();
-            m_emin     = (*this)["emin"].real();
-            m_emax     = (*this)["emax"].real();
-            m_enumbins = (*this)["enumbins"].integer();
-            m_proj     = (*this)["proj"].string();
-            m_coordsys = (*this)["coordsys"].string();
-            m_xref     = (*this)["xref"].real();
-            m_yref     = (*this)["yref"].real();
-            m_binsz    = (*this)["binsz"].real();
-            m_nxpix    = (*this)["nxpix"].integer();
-            m_nypix    = (*this)["nypix"].integer();
+    }
+
+    // Check if we need response information
+    for (int i = 0; i < m_obs.size(); ++i) {
+        GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
+        if (obs != NULL) {
+            if (!obs->has_response()) {
+                m_caldb = (*this)["caldb"].string();
+                m_irf   = (*this)["irf"].string();
+                break;
+            }
         }
     }
  
-    // ... otherwise, read only the parameters that are required for the
-    // model generation.
-    else {
-        if (m_read_ahead) {
-            m_outfile = (*this)["outfile"].filename();
-            m_prefix  = (*this)["prefix"].string();
+    // If no cube file has been specified then read all parameters that
+    // are necessary to create the cube from scratch (see method init_cube)
+    if ((m_infile == "NONE") || (gammalib::strip_whitespace(m_infile) == "")) {
+        m_ebinalg  = (*this)["ebinalg"].string();
+        if (m_ebinalg == "FILE") {
+            m_ebinfile = (*this)["ebinfile"].filename();
         }
-        m_caldb  = (*this)["caldb"].string();
-        m_irf    = (*this)["irf"].string();
-        if (m_obs.models().size() == 0) {
-            m_srcmdl = (*this)["srcmdl"].filename();
+        else {
+            m_emin     = (*this)["emin"].real();
+            m_emax     = (*this)["emax"].real();
+            m_enumbins = (*this)["enumbins"].integer();
         }
+        m_proj     = (*this)["proj"].string();
+        m_coordsys = (*this)["coordsys"].string();
+        m_xref     = (*this)["xref"].real();
+        m_yref     = (*this)["yref"].real();
+        m_binsz    = (*this)["binsz"].real();
+        m_nxpix    = (*this)["nxpix"].integer();
+        m_nypix    = (*this)["nypix"].integer();
     }
  
-    // Get other parameters
-    m_apply_edisp = (*this)["edisp"].boolean();
-
     // Return
     return;
 }
@@ -470,37 +551,33 @@ void ctmodel::get_parameters(void)
 /***********************************************************************//**
  * @brief Setup observation container
  *
- * @exception GException::no_cube
- *            No event cube found in CTA observation.
- *
- * This method sets up the observation container for processing. There are
- * two cases:
+ * This method makes sure that we have an observation container with an
+ * associated model and that all CTA observations have a valid response
+ * functions.
  *
- * If there are no observations in the actual observation container, the
- * method will check in "infile" parameter. If this parameter is "NONE" or
- * empty, the task parameters will be used to construct a model map.
- * Otherwise, the method first tries to interpret the "infile" parameter as
- * a counts map, and attemps loading of the file in an event cube. If this
- * fails, the method tries to interpret the "infile" parameter as an
- * observation definition XML file. If this also fails, an exception will
- * be thrown.
+ * If on input the observation container is empty, the method sets up an
+ * observation container by either loading the information from an
+ * observation definition XML file, by loading an event list or a counts
+ * cube as a single observation, or by constructing an observation from
+ * the scratch from the task parameters.
  *
- * If observations exist already in the observation container, the method
- * will simply keep them.
+ * If no model exists yet in the observation container, the method will load
+ * the model from the file given by the task parameter "srcmdl".
  *
- * Test if all CTA observations contain counts maps.
- *
- * Finally, if no models exist so far in the observation container, the
- * models will be loaded from the model XML file.
+ * The method then loops over all observations in the container and makes
+ * sure that all CTA observations have a valid response function. If no
+ * response function exists so far, the response function will be set
+ * using the "caldb" and "irf" task parameters.
  ***************************************************************************/
 void ctmodel::setup_obs(void)
 {
-    // If there are no observations in the container then try to build some
+    // If there are no observations in the container then get them either
+    // from the observation definition file or from the task parameters
     if (m_obs.size() == 0) {
  
-        // If no input filename has been specified, then create a model map
-        // from the task parameters
-        if ((m_infile == "NONE") || (gammalib::strip_whitespace(m_infile) == "")) {
+        // If no observation definition file has been specified, then create
+        // an observation container from the task parameters
+        if ((m_obsfile == "NONE") || (gammalib::strip_whitespace(m_obsfile) == "")) {
  
             // Set pointing direction
             GCTAPointing pnt;
@@ -508,27 +585,27 @@ void ctmodel::setup_obs(void)
             skydir.radec_deg(m_ra, m_dec);
             pnt.dir(skydir);
  
-            // Setup energy range covered by model
-            GEnergy  emin(m_emin, "TeV");
-            GEnergy  emax(m_emax, "TeV");
-            GEbounds ebds(m_enumbins, emin, emax);
-
             // Setup time interval covered by model
             GGti  gti;
             GTime tmin(m_tmin);
             GTime tmax(m_tmax);
             gti.append(tmin, tmax);
  
-            // Setup skymap
-            GSkymap map = GSkymap(m_proj, m_coordsys,
-                                  m_xref, m_yref, -m_binsz, m_binsz,
-                                  m_nxpix, m_nypix, m_enumbins);
-
-            // Create model cube from sky map
-            GCTAEventCube cube(map, ebds, gti);
+            // Setup energy range covered by model
+            GEnergy  emin(m_emin, "TeV");
+            GEnergy  emax(m_emax, "TeV");
+            GEbounds ebounds(1, emin, emax);
  
-            // Allocate CTA observation
+            // Allocate CTA observation and empty event list
             GCTAObservation obs;
+            GCTAEventList   list;
+
+            // Set event list GTI and energy boundaries and append event list
+            // to observation (we need in fact only the GTI for the
+            // computations)
+            list.gti(gti);
+            list.ebounds(ebounds);
+            obs.events(list);
  
             // Set CTA observation attributes
             obs.pointing(pnt);
@@ -536,34 +613,25 @@ void ctmodel::setup_obs(void)
             obs.livetime(gti.ontime()*m_deadc);
             obs.deadc(m_deadc);
  
-            // Set event cube in observation
-            obs.events(cube);
-
             // Append CTA observation to container
             m_obs.append(obs);
  
-            // Signal that no XML file should be used for storage
-            m_use_xml = false;
-
-        } // endif: created model map from task parameters
+        } // endif: created single CTA observation from task parameters
  
         // ... otherwise try to load information from the file
         else {
  
-            // First try to open the file as a counts map
+            // First try to open the file as an event list or counts map
             try {
  
                 // Allocate CTA observation
                 GCTAObservation obs;
  
-                // Load counts map in CTA observation
-                obs.load(m_infile);
+                // Load CTA observation from FITS file
+                obs.load(m_obsfile);
  
                 // Append CTA observation to container
                 m_obs.append(obs);
-
-                // Signal that no XML file should be used for storage
-                m_use_xml = false;
  
             }
  
@@ -572,10 +640,7 @@ void ctmodel::setup_obs(void)
  
                 // Load observations from XML file. This will throw
                 // an exception if it fails.
-                m_obs.load(m_infile);
-
-                // Signal that XML file should be used for storage
-                m_use_xml = true;
+                m_obs.load(m_obsfile);
  
             }
  
@@ -589,8 +654,7 @@ void ctmodel::setup_obs(void)
         m_obs.models(GModels(m_srcmdl));
     }
  
-    // Check if all CTA observations contain an event cube and setup response
-    // for all observations
+    // Setup response for all observations
     for (int i = 0; i < m_obs.size(); ++i) {
  
         // Is this observation a CTA observation?
@@ -599,12 +663,6 @@ void ctmodel::setup_obs(void)
         // Yes ...
         if (obs != NULL) {
  
-            // Throw an exception if this observation does not contain
-            // an event cube
-            if (dynamic_cast<const GCTAEventCube*>(obs->events()) == NULL) {
-                throw GException::no_cube(G_SETUP_OBS);
-            }
-        
             // Set response if we don't have one
             if (!obs->has_response()) {
  
@@ -635,211 +693,133 @@ void ctmodel::setup_obs(void)
  
  
 /***********************************************************************//**
- * @brief Generate model map
- *
- * @param[in] obs CTA observation pointer.
- * @param[in] models Model container.
+ * @brief Initialise model cube
  *
- * @exception GException::no_cube
- *            No event cube found in CTA observation.
+ * Initialises model cube by either reading the cube definition from an
+ * existing cube file specified by the "infile" parameter, or by constructing
+ * a cube from the energy and spatial binning information provided by the
+ * task parameters. The GTIs of the model cube will be undefined on
+ * exit.
  ***************************************************************************/
-void ctmodel::model_map(GCTAObservation* obs, const GModels& models)
+void ctmodel::init_cube(void)
 {
-    // Continue only if observation pointer is valid
-    if (obs != NULL) {
+    // Clear model cube and GTIs
+    m_cube.clear();
+    m_gti.clear();
  
-        // Get event cube pointer
-        GCTAEventCube* cube = 
-            const_cast<GCTAEventCube*>(dynamic_cast<const GCTAEventCube*>(obs->events()));
+    // If no cube file has been specified then create a cube from the task
+    // parameters
+    if ((m_infile == "NONE") || (gammalib::strip_whitespace(m_infile) == "")) {
  
-        // Throw an exception if the observation does not hold and event
-        // cube
-        if (cube == NULL) {
-            throw GException::no_cube(G_MODEL_MAP);
-        }
+        // Set dummy GTI that is needed for event cube creation
+        GGti gti;
+        gti.append(GTime(0.0), GTime(0.1234));
+    
+        // Set energy boundaries
+        get_ebounds();
  
-        // Initialise statistics
-        double sum = 0.0;
+        // Setup skymap
+        GSkymap map = GSkymap(m_proj, m_coordsys,
+                              m_xref, m_yref, -m_binsz, m_binsz,
+                              m_nxpix, m_nypix, m_enumbins);
  
-        // Loop over all events in counts map
-        for (int i = 0; i < cube->size(); ++i) {
+        // Create model cube
+        m_cube = GCTAEventCube(map, m_ebounds, gti);
  
-            // Get event bin
-            GCTAEventBin* bin = (*cube)[i];
-            
-            // Compute model value for event bin
-            double model = 
-                   models.eval(*(const_cast<const GCTAEventBin*>(bin)), *obs) *
-                   bin->size();
+    }
  
-            // Store value
-            bin->counts(model);
+    // ... otherwise load cube from file
+    else {
  
-            // Sum all events
-            sum += model;
-        }
+        // Load cube
+        m_cube.load(m_infile);
  
-        // Log results
-        if (logTerse()) {
-            log << gammalib::parformat("Model events in cube");
-            log << sum << std::endl;
+        // Set all cube bins to zero
+        for (int i = 0; i < m_cube.size(); ++i) {
+            m_cube[i]->counts(0.0);
         }
  
-    } // endif: observation pointer was not valid
+    }
  
     // Return
     return;
 }
-    
  
-/*==========================================================================
- =                                                                         =
- =                             Private methods                             =
- =                                                                         =
- ==========================================================================*/
  
 /***********************************************************************//**
- * @brief Initialise class members
+ * @brief Fill model into model cube
+ *
+ * @param[in] obs CTA observation.
+ *
+ * Adds the expected number of events for a given observation to the events
+ * that are already found in the model cube. The method also updates the
+ * GTI of the model cube so that cube GTI is a list of the GTIs of all
+ * observations that were used to generate the model cube.
  ***************************************************************************/
-void ctmodel::init_members(void)
+void ctmodel::fill_cube(const GCTAObservation* obs)
 {
-    // Initialise members
-    m_infile.clear();
-    m_outfile.clear();
-    m_prefix.clear();
-    m_caldb.clear();
-    m_irf.clear();
-    m_srcmdl.clear();
-    m_proj.clear();
-    m_coordsys.clear();
-    m_ra          = 0.0;
-    m_dec         = 0.0;
-    m_deadc       = 1.0;
-    m_tmin        = 0.0;
-    m_tmax        = 0.0;
-    m_emin        = 0.0;
-    m_emax        = 0.0;
-    m_enumbins    = 0;
-    m_xref        = 0.0;
-    m_yref        = 0.0;
-    m_binsz       = 0.0;
-    m_nxpix       = 0;
-    m_nypix       = 0;
-    m_apply_edisp = false;
-
-    // Initialise protected members
-    m_obs.clear();
-    m_infiles.clear();
-    m_use_xml    = false;
-    m_read_ahead = false;
-
-    // Set logger properties
-    log.date(true);
-
-    // Return
-    return;
-}
+    // Continue only if observation pointer is valid
+    if (obs != NULL) {
  
+        // Get energy boundaries and GTI references for observation
+        const GEbounds& ebounds = obs->events()->ebounds();
+        const GGti&     gti     = obs->events()->gti();
  
-/***********************************************************************//**
- * @brief Copy class members
- *
- * @param[in] app Application.
- ***************************************************************************/
-void ctmodel::copy_members(const ctmodel& app)
-{
-    // Copy attributes
-    m_infile      = app.m_infile;
-    m_outfile     = app.m_outfile;
-    m_prefix      = app.m_prefix;
-    m_caldb       = app.m_caldb;
-    m_irf         = app.m_irf;
-    m_srcmdl      = app.m_srcmdl;
-    m_proj        = app.m_proj;
-    m_coordsys    = app.m_coordsys;
-    m_ra          = app.m_ra;
-    m_dec         = app.m_dec;
-    m_deadc       = app.m_deadc;
-    m_tmin        = app.m_tmin;
-    m_tmax        = app.m_tmax;
-    m_emin        = app.m_emin;
-    m_emax        = app.m_emax;
-    m_enumbins    = app.m_enumbins;
-    m_xref        = app.m_xref;
-    m_yref        = app.m_yref;
-    m_binsz       = app.m_binsz;
-    m_nxpix       = app.m_nxpix;
-    m_nypix       = app.m_nypix;
-    m_apply_edisp = app.m_apply_edisp;
+        // Initialise statistics
+        double sum              = 0.0;
+        int    num_outside_ebds = 0;
  
-    // Copy protected members
-    m_obs        = app.m_obs;
-    m_infiles    = app.m_infiles;
-    m_use_xml    = app.m_use_xml;
-    m_read_ahead = app.m_read_ahead;
+        // Setup cube GTIs for this observation
+        m_cube.gti(obs->events()->gti());
  
-    // Return
-    return;
-}
+        // Loop over all cube bins
+        for (int i = 0; i < m_cube.size(); ++i) {
  
+            // Get cube bin
+            GCTAEventBin* bin = m_cube[i];
  
-/***********************************************************************//**
- * @brief Delete class members
- ***************************************************************************/
-void ctmodel::free_members(void)
-{
-    // Write separator into logger
-    if (logTerse()) {
-        log << std::endl;
-    }
+            // Skip bin if it is outside the energy range of the observation
+            if (!ebounds.contains(bin->energy())) {
+                num_outside_ebds++;
+                continue;
+            }
  
-    // Return
-    return;
-}
+            // Get actual bin value
+            double value = bin->counts();
+            
+            // Compute model value for cube bin
+            double model = m_obs.models().eval(*bin, *obs) * bin->size();
  
+            // Add model to actual value
+            value += model;
+            sum   += model;
  
-/***********************************************************************//**
- * @brief Set output file name.
- *
- * @param[in] filename Input file name.
- *
- * Converts an input filename into an output filename by prepending the
- * prefix stored in the member m_prefix to the input filename. Any path will
- * be stripped from the input filename. Also a trailing ".gz" will be
- * stripped.
- ***************************************************************************/
-std::string ctmodel::set_outfile_name(const std::string& filename) const
-{
-    // Split input filename into path elements
-    std::vector<std::string> elements = gammalib::split(filename, "/");
+            // Store value
+            bin->counts(value);
  
-    // The last path element is the filename
-    std::string outname = m_prefix + elements[elements.size()-1];
+        } // endfor: looped over all cube bins
  
-    // Strip any ".gz"
-    outname = gammalib::strip_chars(outname, ".gz");
-    
-    // Return output filename
-    return outname;
-}
+        // Append GTIs of observation to list of GTIs
+        m_gti.extend(gti);
  
+        // Update GTIs
+        m_cube.gti(m_gti);
  
-/***********************************************************************//**
- * @brief Save model map in FITS format.
- *
- * Save the model map as a FITS file. The filename of the FITS file is
- * specified by the m_outfile member.
- ***************************************************************************/
-void ctmodel::save_fits(void)
-{
-    // Get output filename
-    m_outfile = (*this)["outfile"].filename();
+        // Log results
+        if (logTerse()) {
+            log << gammalib::parformat("Model events in cube");
+            log << sum << std::endl;
+            log << gammalib::parformat("Cube bins outside energy range");
+            log << num_outside_ebds << std::endl;
+        }
  
-    // Get CTA observation from observation container
-    GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]);
+        // Log cube
+        if (logExplicit()) {
+            log.header2("Model cube");
+            log << m_cube << std::endl;
+        }
  
-    // Save model map
-    save_model_map(obs, m_outfile);
+    } // endif: observation was valid
  
     // Return
     return;
@@ -847,83 +827,76 @@ void ctmodel::save_fits(void)
  
  
 /***********************************************************************//**
- * @brief Save model map(s) in XML format.
+ * @brief Get the energy boundaries
+ *
+ * @exception GException::invalid_value
+ *            No valid energy boundary extension found.
  *
- * Save the model map(s) into FITS files and write the file path information
- * into a XML file. The filename of the XML file is specified by the
- * m_outfile member, the filename(s) of the model map(s) are built by
- * prepending the prefix given by the m_prefix member to the input model
- * map(s) filenames. Any path present in the input filename will be stripped,
- * i.e. the model map(s) will be written in the local working directory
- * (unless a path is specified in the m_prefix member).
+ * Get the energy boundaries according to the user parameters. The method
+ * supports loading of energy boundary information from the EBOUNDS or
+ * ENERGYBINS extension, or setting energy boundaries using a linear
+ * or logarithmical spacing.
+ *
+ * @todo Ultimately this method should go in a support library as it is
+ * used by several ctools.
  ***************************************************************************/
-void ctmodel::save_xml(void)
+void ctmodel::get_ebounds(void)
 {
-    // Get output filename and prefix
-    m_outfile = (*this)["outfile"].filename();
-    m_prefix  = (*this)["prefix"].string();
-
-    // Issue warning if output filename has no .xml suffix
-    std::string suffix = gammalib::tolower(m_outfile.substr(m_outfile.length()-4,4));
-    if (suffix != ".xml") {
-        log << "*** WARNING: Name of observation definition output file \""+
-               m_outfile+"\"" << std::endl;
-        log << "*** WARNING: does not terminate with \".xml\"." << std::endl;
-        log << "*** WARNING: This is not an error, but might be misleading."
-               " It is recommended" << std::endl;
-        log << "*** WARNING: to use the suffix \".xml\" for observation"
-               " definition files." << std::endl;
-    }
-
-    // Loop over all observation in the container
-    for (int i = 0; i < m_obs.size(); ++i) {
-
-        // Get CTA observation
-        GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
-
-        // Handle only CTA observations
-        if (obs != NULL) {
-
-            // Set event output file name
-            std::string outfile = set_outfile_name(m_infiles[i]);
-
-            // Store output file name in observation
-            obs->eventfile(outfile);
-
-            // Save event list
-            save_model_map(obs, outfile);
-
-        } // endif: observation was a CTA observations
+    // Initialse energy boundary information
+    m_ebounds.clear();
+
+    // Check whether energy binning information should be read from a FITS
+    // file ...
+    if (m_ebinalg == "FILE") {
+
+        // Open energy boundary file using the EBOUNDS or ENERGYBINS
+        // extension. Throw an exception if opening fails.
+        GFits file(m_ebinfile);
+        if (file.contains("EBOUNDS")) {
+            file.close();
+            m_ebounds.load(m_ebinfile,"EBOUNDS");
+        }
+        else if (file.contains("ENERGYBINS")) {
+            file.close();
+            m_ebounds.load(m_ebinfile,"ENERGYBINS");
+        }
+        else {
+            file.close();
+            std::string msg = "No extension with name \"EBOUNDS\" or"
+                              " \"ENERGYBINS\" found in FITS file"
+                              " \""+m_ebinfile+"\".\n"
+                              "An \"EBOUNDS\" or \"ENERGYBINS\" extension"
+                              " is required if the parameter \"ebinalg\""
+                              " is set to \"FILE\".";
+            throw GException::invalid_value(G_GET_EBOUNDS, msg);
+        }
+        
+        // Set enumbins parameter to number of ebounds
+        m_enumbins = m_ebounds.size();
  
-    } // endfor: looped over observations
+    } // endif: ebinalg was "FILE"
  
-    // Save observations in XML file
-    m_obs.save(m_outfile);
+    // ... otherwise use a linear or a logarithmically-spaced energy binning
+    else {
  
-    // Return
-    return;
-}
+        // Initialise log mode for ebinning
+        bool log = true;
  
+        // check if algorithm is linear
+        if (m_ebinalg == "LIN") {
+            log = false;
+        }
  
-/***********************************************************************//**
- * @brief Save a single model map into a FITS file
- *
- * @param[in] obs Pointer to CTA observation.
- * @param[in] outfile Output file name.
- *
- * This method saves a single model map into a FITS file. The method does
- * nothing if the observation pointer is not valid.
- ***************************************************************************/
-void ctmodel::save_model_map(const GCTAObservation* obs,
-                             const std::string&     outfile) const
-{
-    // Save only if observation is valid
-    if (obs != NULL) {
+        // todo: should we also check if m_ebinalg is "LOG"
+        // and throw an exception if neither LIN/LOG/FILE
+        // is given?
  
-        // Save observation into FITS file
-        obs->save(outfile, clobber());
+        // Setup energy range covered by data
+        GEnergy  emin(m_emin, "TeV");
+        GEnergy  emax(m_emax, "TeV");
+        m_ebounds = GEbounds(m_enumbins, emin, emax, log);
  
-    } // endif: observation was valid
+    } //endif: ebinalg was not "FILE"
  
     // Return
     return;
 /***************************************************************************
- *                     ctmodel - CTA counts model tool                     *
+ *                ctmodel - CTA model cube generation tool                 *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2012-2014 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
@@ -20,7 +20,7 @@
  ***************************************************************************/
 /**
  * @file ctmodel.hpp
- * @brief CTA counts model tool definition
+ * @brief CTA model cube tool definition
  * @author Juergen Knoedlseder
  */
  
@@ -69,25 +69,23 @@ public:
     void                 run(void);
     void                 save(void);
     const GObservations& obs(void) const;
+    const GCTAEventCube& cube(void) const;
  
 protected:
     // Protected methods
-    void           init_members(void);
-    void           copy_members(const ctmodel& app);
-    void           free_members(void);
-    void           get_parameters(void);
-    void           setup_obs(void);
-    void           model_map(GCTAObservation* obs, const GModels& models);
-    std::string    set_outfile_name(const std::string& filename) const;
-    void           save_fits(void);
-    void           save_xml(void);
-    void           save_model_map(const GCTAObservation* obs,
-                                  const std::string&     outfile) const;
-
+    void init_members(void);
+    void copy_members(const ctmodel& app);
+    void free_members(void);
+    void get_parameters(void);
+    void setup_obs(void);
+    void init_cube(void);
+    void fill_cube(const GCTAObservation* obs);
+    void get_ebounds(void);
+    
     // User parameters
-    std::string m_infile;      //!< Input counts map or XML file
-    std::string m_outfile;     //!< Output model map or XML file
-    std::string m_prefix;      //!< Prefix for multiple model maps
+    std::string m_infile;      //!< Input counts cube
+    std::string m_obsfile;     //!< Event list, counts map or observation definition XML file
+    std::string m_outfile;     //!< Output model cube
     std::string m_caldb;       //!< Calibration database
     std::string m_irf;         //!< Instrument response functions
     std::string m_srcmdl;      //!< Source model
@@ -96,6 +94,8 @@ protected:
     double      m_deadc;       //!< Deadtime correction factor
     double      m_tmin;        //!< Start time
     double      m_tmax;        //!< Stop time
+    std::string m_ebinalg;     //!< Algorithm for energy binning
+    std::string m_ebinfile;    //!< FITS-file containing energy binning
     double      m_emin;        //!< Lower energy
     double      m_emax;        //!< Upper energy
     int         m_enumbins;    //!< Number of energy bins
@@ -109,10 +109,11 @@ protected:
     bool        m_apply_edisp; //!< Apply energy dispersion?
  
     // Protected members
-    GObservations            m_obs;        //!< Observation container
-    std::vector<std::string> m_infiles;    //!< Input event filenames
-    bool                     m_use_xml;    //!< Use XML file instead of FITS file
-    bool                     m_read_ahead; //!< Read ahead parameters
+    GObservations m_obs;        //!< Observation container
+    GCTAEventCube m_cube;       //!< Model cube
+    GEbounds      m_ebounds;    //!< Energy boundaries
+    GGti          m_gti;        //!< Model cube GTIs
+    bool          m_read_ahead; //!< Read ahead parameters
 };
  
  
@@ -127,4 +128,16 @@ const GObservations&amp; ctmodel::obs(void) const
     return m_obs;
 }
  
+
+/***********************************************************************//**
+ * @brief Return model cube
+ *
+ * @return Reference to model cube
+ ***************************************************************************/
+inline
+const GCTAEventCube& ctmodel::cube(void) const
+{
+    return m_cube;
+}
+
 #endif /* CTMODEL_HPP */
@@ -5,7 +5,7 @@
 # ------------------------------------------------------------------------- #
 #                                                                           #
 #  File:      ctmodel.par                                                   #
-#  Version:   00-02-00                                                      #
+#  Version:   00-03-00                                                      #
 #                                                                           #
 #  Author:    Juergen Knoedlseder                                           #
 #             Jurgen.Knodlseder@irap.omp.eu                                 #
@@ -19,22 +19,19 @@
 #===================
 infile,   f, a, "NONE",,, "Input counts map or observation definition file"
 outfile,  f, a, "modmap.fits",,, "Output counts map or observation definition file"
-prefix,   s, h, "modmap_",,, "Prefix for counts maps in observation definition file"
-caldb,    s, a, "$CTOOLS/share/caldb/data/cta/dummy",,, "Calibration database"
-irf,      s, a, "cta_dummy_irf",,, "Instrument response function"
-edisp,    b, h, no,,, "Apply energy dispersion?"
 srcmdl,   f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
+edisp,    b, h, no,,, "Apply energy dispersion?"
 #
-# Instrument parameters
-#======================
-ra,       r, a,    83.63,   0,360, "RA of pointing (degrees)"
-dec,      r, a,    22.01, -90, 90, "Dec of pointing (degrees)"
-deadc,    r, h,    0.95,    0,  1, "Average deadtime correction factor"
-#
-# Time range parameters
-#======================
-tmin,     r, a,     0.0,    0,   , "Start time (MET in s)"
-tmax,     r, a,  1800.0,    0,   , "End time (MET in s)"
+# Observation parameters
+#=======================
+obsfile,  f, a, "NONE",,, "Observation definition file"
+caldb,    s, a, "$CTOOLS/share/caldb/data/cta/dummy",,, "Calibration database"
+irf,      s, a, "cta_dummy_irf",,, "Instrument response function"
+ra,       r, a, 83.63,0,360, "RA of pointing (degrees)"
+dec,      r, a, 22.01,-90,90, "Dec of pointing (degrees)"
+tmin,     r, a, 0.0,0,, "Start time (MET in s)"
+tmax,     r, a, 1800.0,0,, "End time (MET in s)"
+deadc,    r, h, 0.95,0,1, "Average deadtime correction factor"
 #
 # Energy binning parameters
 #==========================
@@ -68,9 +68,10 @@ class Test(gammalib.GPythonTestSuite):
         model = ctools.ctmodel()
         model["infile"].filename("NONE")
         model["outfile"].filename("modmap.fits")
+        model["srcmdl"].filename(self.model_name)
+        model["obsfile"].filename("NONE")
         model["caldb"].string(self.caldb)
         model["irf"].string(self.irf)
-        model["srcmdl"].filename(self.model_name)
         model["ra"].real(83.63)
         model["dec"].real(22.01)
         model["tmin"].real(0.0)
@@ -197,6 +197,7 @@ $ECHO -n &quot;Test ctmodel: &quot;
 #
 # Run 1
 $ctmodel infile="cntmap.fits" \
+         obsfile="cntmap.fits" \
          outfile="modmap1.fits" \
          caldb="irf" \
          irf="cta_dummy_irf" \
@@ -212,6 +213,7 @@ fi
 #
 # Run 3
 $ctmodel infile="NONE" \
+         obsfile="NONE" \
          outfile="modmap3.fits" \
          caldb="irf" \
          irf="cta_dummy_irf" \
@@ -241,6 +243,7 @@ fi
 #
 # Run 2
 $ctmodel infile="cntmap2.fits" \
+         obsfile="cntmap2.fits" \
          outfile="modmap2.fits" \
          caldb="irf" \
          irf="cta_dummy_irf" \