Merge branch 'ctbkgcube-2' into integration

Jürgen Knödlseder
2 parents 535bdb93 3992ea54
Showing 26 changed files with 1309 additions and 35 deletions
.gitignore
ChangeLog
NEWS
configure.ac
pyext/ctbkgcube.i
pyext/ctools/ctools.i
src/Makefile.am
src/ctbin/ctbin.cpp
src/ctbkgcube/Makefile.am
src/ctbkgcube/ctbkgcube.cpp
src/ctbkgcube/ctbkgcube.hpp
src/ctbkgcube/ctbkgcube.par
src/ctbkgcube/main.cpp
src/ctexpcube/ctexpcube.par
src/ctpsfcube/ctpsfcube.par
test/test_ctbin.py
test/test_ctbkgcube.py
test/test_ctcubemask.py
test/test_ctexpcube.py
test/test_ctlike.py
@@ -54,6 +54,7 @@ src/ctselect/ctselect
 src/ctskymap/ctskymap
 src/ctexpcube/ctexpcube
 src/ctpsfcube/ctpsfcube
+src/ctbkgcube/ctbkgcube
 src/ctcubemask/ctcubemask
 test/pfiles/
 test/*.fits
-2014-07-20  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
+2014-07-22  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
         * Version 0.8.0 released
         ========================
@@ -6,7 +6,8 @@
 	Adapt to GammaLib 0.9.0 interface
         ctbin now always merges event lists into a single counts cube
         Add ctexpcube tool (generates exposure cube)
-        Add ctpsfcube tool (generated PSF cube)
+        Add ctpsfcube tool (generates PSF cube)
+        Add ctbkgcube tool (generates background cube)
         Add ctcubemask tool (masks event bins in event cube)
         Support energy dispersion in ctobssim, ctlike and ctmodel	
         Correct ctobssim for thread safety
@@ -2,7 +2,7 @@ New Features and Important Changes in ctools 00-08-00
 Author(s) : Juergen Knoedlseder
-20 July 2014
+22 July 2014
 1. Introduction
@@ -97,7 +97,15 @@ definition file. The cube definition may be taken from a counts cube
 or may be specified by user parameters.
-12. ctcubemask - CTA cube filtering
+12. ctbkgcube - CTA background cube generation
+----------------------------------------------
+The first version of the ctbkgcube tool is included in the release.
+ctbkgcube computes the background cube for an observation definition 
+file. The cube definition may be taken from a counts cube or may be
+specified by user parameters.
+
+
+13. ctcubemask - CTA cube filtering
 -----------------------------------
 The first version of the ctcubemask tool is included in the release.
 ctcubemask allows filtering of bins in a counts cube by setting the
@@ -105,12 +113,12 @@ bin content to -1. Filtering can be done by energy, for a circular
 Region of Interest (ROI) or using ds9 exclusion regions.
-13. Examples
+14. Examples
 ------------
 None
-14. Scripts
+15. Scripts
 -----------
 Add logL and event statistics to cspull output file.
@@ -119,11 +127,11 @@ Convert all tabs into spaces.
 The obsutils script has now been moved into the ctools module namespace.
-15. caldb
+16. caldb
 ---------
 None
-16. models
+17. models
 ----------
 None
@@ -884,6 +884,7 @@ AC_CONFIG_FILES([Makefile
                  src/ctselect/Makefile
                  src/ctexpcube/Makefile
                  src/ctpsfcube/Makefile
+                 src/ctbkgcube/Makefile
                  src/ctcubemask/Makefile
                  src/ctools-setup
                  test/Makefile
+/***************************************************************************
+ *                  ctbkgcube - CTA background cube tool                   *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2014 by Chia-Chun Lu                                     *
+ * ----------------------------------------------------------------------- *
+ *                                                                         *
+ *  This program is free software: you can redistribute it and/or modify   *
+ *  it under the terms of the GNU General Public License as published by   *
+ *  the Free Software Foundation, either version 3 of the License, or      *
+ *  (at your option) any later version.                                    *
+ *                                                                         *
+ *  This program is distributed in the hope that it will be useful,        *
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
+ *  GNU General Public License for more details.                           *
+ *                                                                         *
+ *  You should have received a copy of the GNU General Public License      *
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file ctbkgcube.i
+ * @brief CTA background cube tool definition
+ * @author Chia-Chun Lu
+ */
+
+%{
+/* Put headers and other declarations here that are needed for compilation */
+#include "ctbkgcube.hpp"
+#include "GTools.hpp"
+%}
+
+/***********************************************************************//**
+ * @class ctbkgcube
+ *
+ * @brief CTA background cube tool
+ ***************************************************************************/
+class ctbkgcube : public GApplication  {
+
+public:
+    // Constructors and destructors
+    ctbkgcube(void);
+    explicit ctbkgcube(const GObservations& obs);
+    ctbkgcube(int argc, char *argv[]);
+    ctbkgcube(const ctbkgcube& app);
+    virtual ~ctbkgcube(void);
+
+    // Methods
+    void                 clear(void);
+    void                 execute(void);
+    void                 run(void);
+    void                 save(void);
+    const GObservations& obs(void) const;
+    void                 get_parameters(void);
+    void                 get_obs(void);
+    void                 get_ebounds(void);
+    void                 set_from_cntmap(const std::string& filename);
+    void                 fill_cube(GCTAObservation* obs);
+};
+
+/***********************************************************************//**
+ * @brief CTA background cube tool Python extension
+ ***************************************************************************/
+%extend ctbkgcube {
+    ctbkgcube copy() {
+        return (*self);
+    }
+}
@@ -53,4 +53,5 @@
 %include "ctskymap.i"
 %include "ctexpcube.i"
 %include "ctpsfcube.i"
+%include "ctbkgcube.i"
 %include "ctcubemask.i"
@@ -13,6 +13,7 @@ SUBDIRS = ctbin \
           ctskymap \
           ctexpcube \
           ctpsfcube \
+          ctbkgcube \
           ctcubemask 
 # Scripts to be distributed
@@ -33,6 +34,7 @@ libctools_la_LIBADD = ctbin/libctbin.la \
                       ctskymap/libctskymap.la \
                       ctexpcube/libctexpcube.la \
                       ctpsfcube/libctpsfcube.la \
+                      ctbkgcube/libctbkgcube.la \
                       ctcubemask/libctcubemask.la
 # Distribute scripts
@@ -777,7 +777,7 @@ void ctbin::fill_cube(GCTAObservation* obs)
         // Make sure that the observation holds a CTA event list. If this
         // is not the case then throw an exception.
         if (dynamic_cast<const GCTAEventList*>(obs->events()) == NULL) {
-            throw GException::no_list(G_BIN_EVENTS);
+            throw GException::no_list(G_FILL_CUBE);
         }
         // Initialise binning statistics
+## Process this file with automake to produce Makefile.in
+
+# Program name
+bin_PROGRAMS = ctbkgcube
+
+# Sources
+sources = ctbkgcube.cpp
+
+# Libtool library
+noinst_LTLIBRARIES = libctbkgcube.la
+libctbkgcube_la_SOURCES = $(sources)
+
+# Program definition
+ctbkgcube_SOURCES = main.cpp
+ctbkgcube_LDFLAGS = $(LIB_LINK_OPTION) @LDFLAGS@
+ctbkgcube_CPPFLAGS = @CPPFLAGS@
+ctbkgcube_LDADD = libctbkgcube.la
+
+# Parameter files
+pfilesdir = $(prefix)/syspfiles
+pfiles_DATA = ctbkgcube.par
+/***************************************************************************
+ *                   ctbkgcube - CTA background cube tool                  *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2014 by Chia-Chun Lu                                     *
+ * ----------------------------------------------------------------------- *
+ *                                                                         *
+ *  This program is free software: you can redistribute it and/or modify   *
+ *  it under the terms of the GNU General Public License as published by   *
+ *  the Free Software Foundation, either version 3 of the License, or      *
+ *  (at your option) any later version.                                    *
+ *                                                                         *
+ *  This program is distributed in the hope that it will be useful,        *
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
+ *  GNU General Public License for more details.                           *
+ *                                                                         *
+ *  You should have received a copy of the GNU General Public License      *
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file ctbkgcube.cpp
+ * @brief CTA background cube tool implementation
+ * @author Chia-Chun Lu
+ */
+
+/* __ Includes ___________________________________________________________ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include <cstdio>
+#include "ctbkgcube.hpp"
+#include "GTools.hpp"
+
+/* __ Method name definitions ____________________________________________ */
+#define G_GET_PARAMETERS                        "ctbkgcube::get_parameters()"
+#define G_GET_OBS                                      "ctbkgcube::get_obs()"
+#define G_GET_EBOUNDS                              "ctbkgcube::get_ebounds()"
+#define G_SET_FROM_CNTMAP          "ctbkgcube::set_from_cntmap(std::string&)"
+
+/* __ Debug definitions __________________________________________________ */
+
+/* __ Coding definitions _________________________________________________ */
+
+
+/*==========================================================================
+ =                                                                         =
+ =                        Constructors/destructors                         =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Void constructor
+ ***************************************************************************/
+ctbkgcube::ctbkgcube(void) : GApplication(CTBKGCUBE_NAME, CTBKGCUBE_VERSION)
+{
+    // Initialise members
+    init_members();
+
+    // Write header into logger
+    log_header();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Observations constructor
+ *
+ * param[in] obs Observation container.
+ *
+ * This method creates an instance of the class by copying an existing
+ * observations container.
+ ***************************************************************************/
+ctbkgcube::ctbkgcube(const GObservations& obs) :
+           GApplication(CTBKGCUBE_NAME, CTBKGCUBE_VERSION)
+{
+    // Initialise members
+    init_members();
+
+    // Set observations
+    m_obs = obs;
+
+    // Write header into logger
+    log_header();
+
+    // Return
+    return;
+}
+
+
+
+/***********************************************************************//**
+ * @brief Command line constructor
+ *
+ * @param[in] argc Number of arguments in command line.
+ * @param[in] argv Array of command line arguments.
+ ***************************************************************************/
+ctbkgcube::ctbkgcube(int argc, char *argv[]) :
+           GApplication(CTBKGCUBE_NAME, CTBKGCUBE_VERSION, argc, argv)
+{
+    // Initialise members
+    init_members();
+
+    // Write header into logger
+    log_header();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Copy constructor
+ *
+ * @param[in] app ctbkgcube application.
+ ***************************************************************************/
+ctbkgcube::ctbkgcube(const ctbkgcube& app) : GApplication(app)
+{
+    // Initialise members
+    init_members();
+
+    // Copy members
+    copy_members(app);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Destructor
+ ***************************************************************************/
+ctbkgcube::~ctbkgcube(void)
+{
+    // Free members
+    free_members();
+
+    // Return
+    return;
+}
+
+
+/*==========================================================================
+ =                                                                         =
+ =                               Operators                                 =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Assignment operator
+ *
+ * @param[in] app ctbkgcube application.
+ * @return Returns ctbkgcube application.
+ ***************************************************************************/
+ctbkgcube& ctbkgcube::operator=(const ctbkgcube& app)
+{
+    // Execute only if object is not identical
+    if (this != &app) {
+
+        // Copy base class members
+        this->GApplication::operator=(app);
+
+        // Free members
+        free_members();
+
+        // Initialise members
+        init_members();
+
+        // Copy members
+        copy_members(app);
+
+    } // endif: object was not identical
+
+    // Return this object
+    return *this;
+}
+
+
+/*==========================================================================
+ =                                                                         =
+ =                            Public methods                               =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Clear instance
+ ***************************************************************************/
+void ctbkgcube::clear(void)
+{
+    // Free members
+    free_members();
+    this->GApplication::free_members();
+
+    // Initialise members
+    this->GApplication::init_members();
+    init_members();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Execute application
+ *
+ * This is the main execution method of the ctbkgcube class. It is invoked
+ * when the executable is called from command line. The method generates
+ * the background cube and saves the result.
+ ***************************************************************************/
+void ctbkgcube::execute(void)
+{
+    // Signal that some parameters should be read ahead
+    m_read_ahead = true;
+
+    // Create the background cube
+    run();
+
+    // Save the background cube into the output FITS file
+    save();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Generate the background cube(s).
+ *
+ * This method reads the task parameters from the parfile, sets up the
+ * observation container, loops over all CTA observations in the container
+ * and generates a background cube from the CTA observations.
+ ***************************************************************************/
+void ctbkgcube::run(void)
+{
+    // If we're in debug mode then all output is also dumped on the screen
+    if (logDebug()) {
+        log.cout(true);
+    }
+
+    // Get task parameters
+    get_parameters();
+
+    // Write parameters into logger
+    if (logTerse()) {
+        log_parameters();
+        log << std::endl;
+    }
+
+    // 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;
+    }
+
+    // Write header
+    if (logTerse()) {
+        log << std::endl;
+        log.header1("Prepare model");
+    }
+
+    // Copy models from observation container
+    m_bkgmdl = m_obs.models();
+
+    // Remove all model that are not CTA background models from the
+    // container
+    int num = m_bkgmdl.size();
+    for (int i = num-1; i >= 0; --i) {
+
+        // Flag removal
+        bool remove = true;
+
+        // Do we have a CTA, HESS or VERITAS specific model?
+        if (m_bkgmdl[i]->is_valid("CTA", "") ||
+            m_bkgmdl[i]->is_valid("HESS", "") ||
+            m_bkgmdl[i]->is_valid("VERITAS", "")) {
+
+            // Do we have a background model?
+            if (dynamic_cast<GModelData*>(m_bkgmdl[i]) != NULL) {
+
+                // ... then keep model
+                remove = false;
+
+            } // endif: had a background model
+        } // endif: had a CTA, HESS or VERITAS model
+
+        // Log results
+        if (logTerse()) {
+            if (remove) {
+                log << gammalib::parformat("Remove model");
+            }
+            else {
+                log << gammalib::parformat("Keep model");
+            }
+            log << m_bkgmdl[i]->name();
+            log << " ";
+            log << m_bkgmdl[i]->type();
+            log << "(";
+            log << m_bkgmdl[i]->instruments();
+            log << ")";
+            log << std::endl;
+        }
+
+        // Remove model if requested
+        if (remove) {
+            m_bkgmdl.remove(i);
+        }
+
+    } // endfor: looped over all background models
+
+    // Write header
+    if (logTerse()) {
+        log << std::endl;
+        log.header1("Generate background cube");
+    }
+
+    // Initialise observation counter
+    int n_observations = 0;
+
+    // Loop over all observations in the container
+    for (int i = 0; i < m_obs.size(); ++i) {
+      
+        // Get CTA observation
+        GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
+
+        // Continue only if observation is a CTA observation
+        if (obs != NULL) {
+
+            // Write header for observation
+            if (logTerse()) {
+                if (obs->name().length() > 1) {
+                    log.header3("Observation "+obs->name());
+                }
+                else {
+                    log.header3("Observation "+gammalib::str(i));
+                }
+            }
+
+            // Increment number of observations
+            n_observations++;
+
+            // Fill the cube
+            fill_cube(obs);
+
+        } // endif: CTA observation found
+
+    } // endfor: looped over observations
+
+    // Log results
+    if (logTerse()) {
+        log << gammalib::parformat("Number of observations");
+        log << n_observations << std::endl;
+    }
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Generate background cube
+ *
+ * @param[in] obs Pointer to CTA observation.
+ *
+ * Fills
+ ***************************************************************************/
+void ctbkgcube::fill_cube(GCTAObservation* obs)
+{
+    // Continue only if observation pointer is valid
+    if (obs != NULL) {
+
+        // Initialise statistics
+        double sum = 0.0;
+
+        // Set GTI of actual observations as the GTI of the event cube
+        m_bkgcube.gti(obs->events()->gti());
+
+        // Loop over all bins in background cube
+        for (int i = 0; i < m_bkgcube.size(); ++i) {
+
+            // Get event bin
+            GCTAEventBin* bin = m_bkgcube[i];
+            
+            // Compute model value for event bin
+            double model = 
+                   m_bkgmdl.eval(*(const_cast<const GCTAEventBin*>(bin)), *obs) *
+                   bin->size();
+
+            // Add existing number of counts
+            model += bin->counts();
+
+            // Store cumulated value
+            bin->counts(model);
+
+            // Sum all events
+            sum += model;
+        }
+
+        // Log results
+        if (logTerse()) {
+            log << gammalib::parformat("Background events in cube");
+            log << sum << std::endl;
+        }
+
+    } // endif: observation pointer was not valid
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Save background cube
+ *
+ * Save the background cube into the file specified by the outfile parameter.
+ ***************************************************************************/
+void ctbkgcube::save(void)
+{
+    // Write header
+    if (logTerse()) {
+        log << std::endl;
+        log.header1("Save background cube");
+    }
+
+    // Create energies container from energy boundaries
+    GEnergies energies;
+    for (int i = 0; i < m_ebounds.size(); ++i) {
+        energies.append(m_ebounds.elogmean(i));
+    }
+
+    // Get output filename
+    std::string filename = (*this)["outfile"].filename();
+
+    // Create empty FITS file
+    GFits fits;
+
+    // Write background cube
+    m_bkgcube.map().write(fits);
+
+    // Write energies
+    energies.write(fits);
+    
+    // Save FITS file
+    fits.saveto(filename, clobber());
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Get application parameters
+ *
+ * @exception GException::invalid_value
+ *            No background model definition XML file specified.
+ *
+ * Get all task parameters from parameter file or (if required) by querying
+ * the user. The parameters are read in the correct order.
+ ***************************************************************************/
+void ctbkgcube::get_parameters(void)
+{
+    // If we do not have any observations in the container then load them
+    // using the "infile" parameter
+    if (m_obs.size() == 0) {
+        get_obs();
+    }
+
+    // If there are no models associated with the observations then load now
+    // the model definition from the XML file
+    if (m_obs.models().size() == 0) {
+        if ((*this)["bkgmdl"].is_undefined()) {
+            std::string msg = "No background model definition XML file "
+                              "specified. Please set the \"bkgmdl\" parameter "
+                              "to the background model definition XML "
+                              "filename.";
+            throw GException::invalid_value(G_GET_PARAMETERS, msg);
+        }
+        std::string bkgmdl = (*this)["bkgmdl"].filename();
+        GModels     models(bkgmdl);
+        m_obs.models(models);
+    }
+
+    // If no counts map is specified then setup the background cube from
+    // the user parameters
+    std::string cntmap = (*this)["cntmap"].filename();
+    if ((gammalib::toupper(cntmap) == "NONE") ||
+        (gammalib::strip_whitespace(cntmap) == "")) {
+    
+        // Get user parameters for counts map definition
+        std::string proj     = (*this)["proj"].string();
+        std::string coordsys = (*this)["coordsys"].string();
+        double      xref     = (*this)["xref"].real();
+        double      yref     = (*this)["yref"].real();
+        double      binsz    = (*this)["binsz"].real();
+        int         nxpix    = (*this)["nxpix"].integer();
+        int         nypix    = (*this)["nypix"].integer();
+
+        // Get energy definition
+        get_ebounds();
+
+        // Define skymap for background cube
+        GSkymap map(proj, coordsys, xref, yref,
+                    -binsz, binsz, nxpix, nypix,
+                    m_ebounds.size());
+
+        // Allocate background cube using a dummy GTI
+        GGti gti;
+        gti.append(GTime(0.0), GTime(1.0)); // Dummy GTI
+        m_bkgcube = GCTAEventCube(map, m_ebounds, gti);
+
+    }
+
+    // ... otherwise setup the background cube from the counts map
+    else {
+    
+        // Set background cube from counts map
+        set_from_cntmap(cntmap);
+    
+    }
+
+    // Read output filename (if needed)
+    if (m_read_ahead) {
+        std::string filename = (*this)["outfile"].filename();
+    }
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Get observation definition
+ *
+ * @exception GException::invalid_value
+ *            No event list, counts cube or observation definition file
+ *            specified.
+ *
+ * Get observation definition from the user parameters.
+ ***************************************************************************/
+void ctbkgcube::get_obs(void)
+{
+    // Get input filename
+    std::string filename = (*this)["infile"].filename();
+
+    // Check whether infile parameter is undefined
+    if (gammalib::toupper(filename) == "NONE") {
+        std::string msg = "No event list, counts cube or observation definition "
+                          "file specified. Please set the \"infile\" parameter "
+                          "to either an event list filename, a counts cube "
+                          "filename or an observation definition filename.";
+        throw GException::invalid_value(G_GET_OBS, msg);
+    }
+
+    // Try first to open as FITS file
+    try {
+
+        // Allocate CTA observation
+        GCTAObservation obs;
+        
+        // Load input file in CTA observation
+        obs.load(filename);
+
+        // Append CTA observation to container
+        m_obs.append(obs);
+            
+    }
+        
+    // ... otherwise try to open as XML file
+    catch (GException::fits_open_error &e) {
+
+        // Load observations from XML file
+        m_obs.load(filename);
+
+    }
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Get energy boundaries from parameters
+ *
+ * @exception GException::invalid_value
+ *            Energy bin definition file not defined or invalid extension
+ *            name encountered.
+ *
+ * Get the energy boundaries from the user parameters.
+ ***************************************************************************/
+void ctbkgcube::get_ebounds(void)
+{
+    // Determine the energy binning alogrithm
+    std::string ebinalg = (*this)["ebinalg"].string();
+
+    // If we have the binning given by a file then try to get the boundaries
+    // from that file
+    if (ebinalg == "FILE") {
+
+        // Get filename
+        std::string filename = (*this)["ebinfile"].filename();
+
+        // Check whether infile parameter is undefined
+        if (gammalib::toupper(filename) == "NONE") {
+            std::string msg = "No energy bin definition file specified. "
+                              "Please set the \"ebinfile\" parameter to a "
+                              "file containing the energy bin definition.";
+            throw GException::invalid_value(G_GET_EBOUNDS, msg);
+        }
+
+        // Open fits file to check which extension is given
+        GFits file(filename);
+
+        // Check first for EBOUNDS extension
+        if (file.contains("EBOUNDS")) {
+            file.close();
+            m_ebounds.load(filename, "EBOUNDS");
+        }
+
+        // ... then check for ENERGYBINS extension
+        else if (file.contains("ENERGYBINS")) {
+            file.close();
+            m_ebounds.load(filename, "ENERGYBINS");
+        }
+
+        // ... otherwise throw an exception
+        else {
+            file.close();
+            std::string msg = "No extension with name \"EBOUNDS\" or"
+                              " \"ENERGYBINS\" found in FITS file"
+                              " \""+filename+"\".\n"
+                              "An \"EBOUNDS\" or \"ENERGYBINS\" extension"
+                              " is required if the parameter \"ebinalg\""
+                              " is set to \"FILE\".";
+            throw GException::invalid_value(G_GET_EBOUNDS, msg);
+        }
+    }
+    
+    // ... otherwise read emin, emax and nebins
+    else {
+
+        // Get the relevant parameters
+    	double emin     = (*this)["emin"].real();
+    	double emax     = (*this)["emax"].real();
+    	int    enumbins = (*this)["enumbins"].integer();
+        bool   log      = ((*this)["ebinalg"].string() == "LIN") ? false : true;
+
+        // Create energy boundaries
+        m_ebounds = GEbounds(enumbins,
+                             GEnergy(emin, "TeV"),
+                             GEnergy(emax, "TeV"),
+                             log);
+
+    }
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Set background cube definition from counts map
+ *
+ * @exception GException::invalid_value
+ *            Invalid counts map projection or invalid events encountered.
+ *
+ * Set background cube definition from counts map.
+ ***************************************************************************/
+void ctbkgcube::set_from_cntmap(const std::string& filename)
+{
+    // Allocate CTA observation
+    GCTAObservation obs;
+        
+    // Load counts map in CTA observation
+    obs.load(filename);
+
+    // Set background cube from counts map
+    const GCTAEventCube* cube = dynamic_cast<const GCTAEventCube*>(obs.events());
+
+    // Continue only if cube is valid
+    if (cube != NULL) {
+
+        // Get sky map projection
+        const GWcs* wcs = dynamic_cast<const GWcs*>(cube->map().projection());
+        
+        // Continue only if projection is valid
+        if (wcs != NULL) {
+            
+            // Get user parameters for counts map definition
+            std::string proj     = wcs->code();
+            std::string coordsys = wcs->coordsys();
+            double      xref     = wcs->crval(0);
+            double      yref     = wcs->crval(1);
+            double      dx       = wcs->cdelt(0);
+            double      dy       = wcs->cdelt(1);
+            int         nx       = cube->map().nx();
+            int         ny       = cube->map().ny();
+
+            // Get energy definition
+            m_ebounds = cube->ebounds();
+
+            // Define skymap for background cube
+            GSkymap map(proj, coordsys, xref, yref,
+                        dx, dy, nx, ny,
+                        m_ebounds.size());
+
+            // Allocate background cube
+            m_bkgcube = GCTAEventCube(map, m_ebounds, obs.events()->gti());
+
+        } // endif: WCS projection was valid
+
+        // ... projection is not of WCS type
+        else {
+            std::string msg = "Counts map project is not of WCS type.";
+            throw GException::invalid_value(G_SET_FROM_CNTMAP, msg);
+        }
+
+    } // endif: observation contained an events cube
+
+    // ... there is not events cube
+    else {
+        std::string msg = "No events cube found in file \""
+                          ""+filename+"\".";
+        throw GException::invalid_value(G_SET_FROM_CNTMAP, msg);
+    }
+
+    // Return
+    return;
+}
+
+
+/*==========================================================================
+ =                                                                         =
+ =                             Private methods                             =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Initialise class members
+ ***************************************************************************/
+void ctbkgcube::init_members(void)
+{
+    // Initialise members
+    m_read_ahead = false;
+    m_obs.clear();
+    m_bkgcube.clear();
+    m_bkgmdl.clear();
+    m_ebounds.clear();
+
+    // Set logger properties
+    log.date(true);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Copy class members
+ *
+ * @param[in] app Application.
+ ***************************************************************************/
+void ctbkgcube::copy_members(const ctbkgcube& app)
+{
+    // Copy members
+    m_read_ahead = app.m_read_ahead;
+    m_obs        = app.m_obs;
+    m_bkgcube    = app.m_bkgcube;
+    m_bkgmdl     = app.m_bkgmdl;
+    m_ebounds    = app.m_ebounds;
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Delete class members
+ ***************************************************************************/
+void ctbkgcube::free_members(void)
+{
+    // Write separator into logger
+    if (logTerse()) {
+        log << std::endl;
+    }
+
+    // Return
+    return;
+}
+/***************************************************************************
+ *                   ctbkgcube - CTA background cube tool                  *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2014 by Chia-Chun Lu                                     *
+ * ----------------------------------------------------------------------- *
+ *                                                                         *
+ *  This program is free software: you can redistribute it and/or modify   *
+ *  it under the terms of the GNU General Public License as published by   *
+ *  the Free Software Foundation, either version 3 of the License, or      *
+ *  (at your option) any later version.                                    *
+ *                                                                         *
+ *  This program is distributed in the hope that it will be useful,        *
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
+ *  GNU General Public License for more details.                           *
+ *                                                                         *
+ *  You should have received a copy of the GNU General Public License      *
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file ctbkgcube.hpp
+ * @brief CTA background cube tool definition
+ * @author Chia-Chun Lu
+ */
+
+#ifndef CTBKGCUBE_HPP
+#define CTBKGCUBE_HPP
+
+/* __ Includes ___________________________________________________________ */
+#include "GammaLib.hpp"
+#include "GCTALib.hpp"
+
+/* __Definitions _________________________________________________________ */
+#define CTBKGCUBE_NAME    "ctbkgcube"
+#define CTBKGCUBE_VERSION "00-01-00"
+
+
+/***********************************************************************//**
+ * @class ctbkgcube
+ *
+ * @brief CTA background cube tool
+ ***************************************************************************/
+class ctbkgcube : public GApplication  {
+
+public:
+    // Constructors and destructors
+    ctbkgcube(void);
+    explicit ctbkgcube(const GObservations& obs);
+    ctbkgcube(int argc, char *argv[]);
+    ctbkgcube(const ctbkgcube& app);
+    virtual ~ctbkgcube(void);
+
+    // Operators
+    ctbkgcube& operator=(const ctbkgcube& app);
+
+    // Methods
+    void                 clear(void);
+    void                 execute(void);
+    void                 run(void);
+    void                 save(void);
+    const GObservations& obs(void) const;
+    void                 get_parameters(void);
+    void                 get_obs(void);
+    void                 get_ebounds(void);
+    void                 set_from_cntmap(const std::string& filename);
+    void                 fill_cube(GCTAObservation* obs);
+
+protected:
+    // Protected methods
+    void init_members(void);
+    void copy_members(const ctbkgcube& app);
+    void free_members(void);
+
+    // Protected members
+    bool          m_read_ahead;  //!< Read ahead parameters
+    GObservations m_obs;         //!< Observation container
+    GCTAEventCube m_bkgcube;     //!< Background cube
+    GModels       m_bkgmdl;      //!< CTA background models
+    GEbounds      m_ebounds;     //!< Energy boundaries
+};
+
+
+/***********************************************************************//**
+ * @brief Return observation container
+ *
+ * @return Observation container.
+ ***************************************************************************/
+inline
+const GObservations& ctbkgcube::obs(void) const
+{
+    // Return observation container
+    return m_obs;
+}
+
+#endif /* CTBKGCUBE_HPP */
+#############################################################################
+#                                                                           #
+#                       CTA Science Analysis Tools                          #
+#                                                                           #
+# ------------------------------------------------------------------------- #
+#                                                                           #
+#  File:      ctbkgcube.par                                                 #
+#  Version:   00-01-00                                                      #
+#                                                                           #
+#  Author:    Chia-Chun Lu                                                  #
+#             c.lu@mpi-hd.mpg.de                                            #
+#             MPIK                                                          #
+#                                                                           #
+#  Purpose:   Parameter file for the ctbkgcube application                  #
+#                                                                           #
+#############################################################################
+#
+# General parameters
+#===================
+infile,   f, a, "NONE",,, "Event list, counts cube or observation definition file"
+bkgmdl,   f, a, "NONE",,, "Background model XML file"
+cntmap,   f, a, "NONE",,, "Counts cube for background cube definition"
+outfile,  f, a, "bkgcube.fits",,, "Output background cube file"
+#
+# Energy binning parameters
+#==========================
+ebinalg,  s, h, "LOG", FILE|LIN|LOG,,"Algorithm for defining energy bins"
+emin,     r, a,  0.1,,,"Start value for first energy bin in TeV"
+emax,     r, a, 100.0,,,"Stop value for last energy bin in TeV"
+enumbins, i, a, 20,,,"Number of energy bins"
+ebinfile, f, h, "NONE",,,"Name of the file containing the energy bin definition"
+#
+# Spatial binning parameters
+#===========================
+nxpix,    i, a, 200,,, "Size of the X axis in pixels"
+nypix,    i, a, 200,,, "Size of the Y axis in pixels"
+binsz,    r, a, 0.02,,, "Image scale (in degrees/pixel)"
+coordsys, s, a, "CEL", CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)"
+xref,     r, a, 83.63,,, "First coordinate of image center in degrees (RA or galactic l)"
+yref,     r, a, 22.01,,, "Second coordinate of image center in degrees (DEC or galactic b)"
+axisrot,  r, h, 0.0,,, "Rotation angle of image axis, in degrees"
+proj,     s, a, "CAR", AIT|AZP|CAR|MER|STG|TAN,, "Projection method e.g. AIT|AZP|CAR|MER|STG|TAN"
+#
+# Standard parameters
+#====================
+chatter, i, h, 2,0,4, "Chattiness of output"
+clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
+debug,   b, h, no,,, "Debugging mode activated"
+mode,    s, h, "ql",,, "Mode of automatic parameters"
+logfile, f, h, "ctbkgcube.log",,, "Log filename"
+/***************************************************************************
+ *                  ctbkgcube - CTA background cube tool                   *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2014 by Chia-Chun Lu                                     *
+ * ----------------------------------------------------------------------- *
+ *                                                                         *
+ *  This program is free software: you can redistribute it and/or modify   *
+ *  it under the terms of the GNU General Public License as published by   *
+ *  the Free Software Foundation, either version 3 of the License, or      *
+ *  (at your option) any later version.                                    *
+ *                                                                         *
+ *  This program is distributed in the hope that it will be useful,        *
+ *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
+ *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
+ *  GNU General Public License for more details.                           *
+ *                                                                         *
+ *  You should have received a copy of the GNU General Public License      *
+ *  along with this program.  If not, see <http://www.gnu.org/licenses/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file ctbkgcube/main.cpp
+ * @brief CTA background cube tool main code
+ * @author Chia-Chun Lu
+ */
+
+/* __ Includes ___________________________________________________________ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include <cstdio>
+#include "ctbkgcube.hpp"
+
+
+/***********************************************************************//**
+ * @brief Main entry point of application
+ *
+ * @param[in] argc Number of command line arguments.
+ * @param[in] argv Command line arguments.
+ *
+ * This is the main entry point of the ctbkgcube application. It allocates a
+ * ctbkgcube object and runs the application.
+ ***************************************************************************/
+int main (int argc, char *argv[])
+{
+    // Initialise return code
+    int rc = 1;
+
+    // Create instance of application
+    ctbkgcube application(argc, argv);
+
+    // Run application
+    try {
+        // Execute application
+        application.execute();
+
+        // Signal success
+        rc = 0;
+    }
+    catch (std::exception &e) {
+
+        // Extract error message
+        std::string message = e.what();
+        std::string signal  = "*** ERROR encounterted in the execution of"
+                              " ctbkgcube. Run aborted ...";
+
+        // Write error in logger
+        application.log << signal  << std::endl;
+        application.log << message << std::endl;
+
+        // Write error on standard output
+        std::cout << signal  << std::endl;
+        std::cout << message << std::endl;
+
+    } // endcatch: catched any application error
+
+    // Return
+    return rc;
+}
@@ -17,7 +17,7 @@
 #
 # General parameters
 #===================
-infile,   f, a, "NONE",,, "Observation definition file"
+infile,   f, a, "NONE",,, "Event list, counts cube or observation definition file"
 cntmap,   f, a, "NONE",,, "Counts map for exposure cube definition"
 outfile,  f, a, "expcube.fits",,, "Output exposure cube file"
 caldb,    s, a, "dummy",,, "Calibration database"
@@ -17,7 +17,7 @@
 #
 # General parameters
 #===================
-infile,   f, a, "NONE",,, "Observation definition file"
+infile,   f, a, "NONE",,, "Event list, counts cube or observation definition file"
 cntmap,   f, a, "NONE",,, "Counts map for exposure cube definition"
 outfile,  f, a, "psfcube.fits",,, "Output exposure cube file"
 caldb,    s, a, "dummy",,, "Calibration database"
@@ -22,9 +22,9 @@ import gammalib
 import ctools
-# ============================ #
-# Test class for ctobssim tool #
-# ============================ #
+# ========================= #
+# Test class for ctbin tool #
+# ========================= #
 class Test(gammalib.GPythonTestSuite):
     """
     Test class for ctbin tool.
+#! /usr/bin/env python
+# ==========================================================================
+# This scripts performs unit tests for the ctbkgcube tool.
+#
+# Copyright (C) 2014 Juergen Knoedlseder
+#
+# This program is free software: you can redistribute it and/or modify
+# it under the terms of the GNU General Public License as published by
+# the Free Software Foundation, either version 3 of the License, or
+# (at your option) any later version.
+#
+# This program is distributed in the hope that it will be useful,
+# but WITHOUT ANY WARRANTY; without even the implied warranty of
+# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+# GNU General Public License for more details.
+#
+# You should have received a copy of the GNU General Public License
+# along with this program.  If not, see <http://www.gnu.org/licenses/>.
+#
+# ==========================================================================
+import gammalib
+import ctools
+
+
+# ============================= #
+# Test class for ctbkgcube tool #
+# ============================= #
+class Test(gammalib.GPythonTestSuite):
+    """
+    Test class for ctbkgcube tool.
+    """
+    # Constructor
+    def __init__(self):
+        """
+        Constructor.
+        """
+        # Call base class constructor
+        gammalib.GPythonTestSuite.__init__(self)
+
+        # Set members
+        self.events_name = "data/crab_events.fits"
+        self.bkg_model   = "data/crab.xml"
+
+        # Return
+        return
+
+    # Set test functions
+    def set(self):
+        """
+        Set all test functions.
+        """
+        # Set test name
+        self.name("ctbkgcube")
+
+        # Append tests
+        self.append(self.test_functional, "Test ctbkgcube functionality")
+
+        # Return
+        return
+
+    # Test ctbkgcube functionnality
+    def test_functional(self):
+        """
+        Test ctbkgcube functionnality.
+        """
+        # Set-up ctbkgcube
+        bkgcube = ctools.ctbkgcube()
+        bkgcube["infile"].filename(self.events_name)
+        bkgcube["bkgmdl"].filename(self.bkg_model)
+        bkgcube["cntmap"].filename("NONE")
+        bkgcube["outfile"].filename("psfcube.fits")
+        bkgcube["ebinalg"].string("LOG")
+        bkgcube["emin"].real(0.1)
+        bkgcube["emax"].real(100.0)
+        bkgcube["enumbins"].integer(20)
+        bkgcube["nxpix"].integer(10)
+        bkgcube["nypix"].integer(10)
+        bkgcube["binsz"].real(0.4)
+        bkgcube["coordsys"].string("CEL")
+        bkgcube["proj"].string("CAR")
+        bkgcube["xref"].real(83.63)
+        bkgcube["yref"].real(22.01)
+        
+        # Run tool
+        self.test_try("Run ctbkgcube")
+        try:
+            bkgcube.run()
+            self.test_try_success()
+        except:
+            self.test_try_failure("Exception occured in ctbkgcube.")
+
+        # Save background cube
+        self.test_try("Save background cube")
+        try:
+            bkgcube.save()
+            self.test_try_success()
+        except:
+            self.test_try_failure("Exception occured in saving background cube.")
+
+        # Return
+        return
@@ -22,9 +22,9 @@ import gammalib
 import ctools
-# ============================ #
-# Test class for ctobssim tool #
-# ============================ #
+# ============================== #
+# Test class for ctcubemask tool #
+# ============================== #
 class Test(gammalib.GPythonTestSuite):
     """
     Test class for ctcubemask tool.
@@ -22,9 +22,9 @@ import gammalib
 import ctools
-# ============================ #
-# Test class for ctobssim tool #
-# ============================ #
+# ============================= #
+# Test class for ctexpcube tool #
+# ============================= #
 class Test(gammalib.GPythonTestSuite):
     """
     Test class for ctexpcube tool.
@@ -22,9 +22,9 @@ import gammalib
 import ctools
-# ============================ #
-# Test class for ctobssim tool #
-# ============================ #
+# ========================== #
+# Test class for ctlike tool #
+# ========================== #
 class Test(gammalib.GPythonTestSuite):
     """
     Test class for ctlike tool.
@@ -22,9 +22,9 @@ import gammalib
 import ctools
-# ============================ #
-# Test class for ctobssim tool #
-# ============================ #
+# =========================== #
+# Test class for ctmodel tool #
+# =========================== #
 class Test(gammalib.GPythonTestSuite):
     """
     Test class for ctmodel tool.
@@ -42,6 +42,7 @@ ctselect=../src/ctselect/ctselect
 ctskymap=../src/ctskymap/ctskymap
 ctexpcube=../src/ctexpcube/ctexpcube
 ctpsfcube=../src/ctpsfcube/ctpsfcube
+ctbkgcube=../src/ctbkgcube/ctbkgcube
 ctcubemask=../src/ctcubemask/ctcubemask
@@ -444,6 +445,51 @@ $ECHO &quot; ok&quot;
 #
+# Test ctbkgcube
+# ==============
+$ECHO -n "Test ctbkgcube: "
+#
+# Run 1
+$ctbkgcube infile="data/crab_events.fits" \
+           bkgmdl="data/crab.xml" \
+           cntmap="NONE" \
+           outfile="bkgcube1.fits" \
+           emin=0.1 \
+           emax=100.0 \
+           enumbins=20 \
+           nxpix=10 \
+           nypix=10 \
+           binsz=0.4 \
+           coordsys="CEL" \
+           xref=83.63 \
+           yref=22.01 \
+           proj="CAR"
+$ECHO -n "."
+if [ -s "bkgcube1.fits" ]
+then
+  $ECHO -n "."
+else
+  $ECHO " bkgcube1.fits file is not valid"
+  exit 1
+fi
+#
+# Run 2
+$ctbkgcube infile="data/crab_events.fits" \
+           bkgmdl="data/crab.xml" \
+           cntmap="data/crab_cntmap.fits" \
+           outfile="bkgcube2.fits"
+$ECHO -n "."
+if [ -s "bkgcube2.fits" ]
+then
+  $ECHO -n "."
+else
+  $ECHO " bkgcube2.fits file is not valid"
+  exit 1
+fi
+$ECHO " ok"
+
+
+#
 # Test ctcubemask
 # ===============
 $ECHO -n "Test ctcubemask: "
@@ -22,9 +22,9 @@ import gammalib
 import ctools
-# ============================ #
-# Test class for ctobssim tool #
-# ============================ #
+# ============================= #
+# Test class for ctpsfcube tool #
+# ============================= #
 class Test(gammalib.GPythonTestSuite):
     """
     Test class for ctpsfcube tool.
@@ -23,7 +23,7 @@ import ctools
 # ============================ #
-# Test class for ctobssim tool #
+# Test class for ctselect tool #
 # ============================ #
 class Test(gammalib.GPythonTestSuite):
     """
@@ -23,7 +23,7 @@ import ctools
 # ============================ #
-# Test class for ctobssim tool #
+# Test class for ctskymap tool #
 # ============================ #
 class Test(gammalib.GPythonTestSuite):
     """
@@ -28,9 +28,10 @@ import test_ctbin
 import test_ctlike
 import test_ctmodel
 import test_ctskymap
-import test_ctcubemask
 import test_ctexpcube
 import test_ctpsfcube
+import test_ctbkgcube
+import test_ctcubemask
 import test_pipelines
@@ -51,9 +52,10 @@ if __name__ == &#39;__main__&#39;:
     suite_ctlike     = test_ctlike.Test()
     suite_ctmodel    = test_ctmodel.Test()
     suite_ctskymap   = test_ctskymap.Test()
-    suite_ctcubemask = test_ctcubemask.Test()
     suite_ctexpcube  = test_ctexpcube.Test()
     suite_ctpsfcube  = test_ctpsfcube.Test()
+    suite_ctbkgcube  = test_ctbkgcube.Test()
+    suite_ctcubemask = test_ctcubemask.Test()
     suite_pipelines  = test_pipelines.Test()
     # Setup unit tests
@@ -63,9 +65,10 @@ if __name__ == &#39;__main__&#39;:
     suite_ctlike.set()
     suite_ctmodel.set()
     suite_ctskymap.set()
-    suite_ctcubemask.set()
     suite_ctexpcube.set()
     suite_ctpsfcube.set()
+    suite_ctbkgcube.set()
+    suite_ctcubemask.set()
     suite_pipelines.set()
     # Append tests to container
@@ -75,9 +78,10 @@ if __name__ == &#39;__main__&#39;:
     suites.append(suite_ctlike)
     suites.append(suite_ctmodel)
     suites.append(suite_ctskymap)
-    suites.append(suite_ctcubemask)
     suites.append(suite_ctexpcube)
     suites.append(suite_ctpsfcube)
+    suites.append(suite_ctbkgcube)
+    suites.append(suite_ctcubemask)
     suites.append(suite_pipelines)
     # Set PFILES environment variable