Add initial version of ctmapcube tool (so far the initial version generates an empty map cube)

Jürgen Knödlseder
1 parent d98a593d
Showing 13 changed files with 962 additions and 0 deletions
.gitignore
configure.ac
pyext/Makefile.am
pyext/ctmapcube.i
pyext/ctools/tools.i
src/Makefile.am
src/ctmapcube/Makefile.am
src/ctmapcube/ctmapcube.cpp
src/ctmapcube/ctmapcube.hpp
src/ctmapcube/ctmapcube.par
src/ctmapcube/main.cpp
test/test_ctmapcube.py
test/test_python_ctools.py
@@ -65,6 +65,7 @@ src/ctexpcube/ctexpcube
 src/ctpsfcube/ctpsfcube
 src/ctedispcube/ctedispcube
 src/ctbkgcube/ctbkgcube
+src/ctmapcube/ctmapcube
 src/ctcubemask/ctcubemask
 src/ctbutterfly/ctbutterfly
 src/ctulimit/ctulimit
@@ -1013,6 +1013,7 @@ AC_CONFIG_FILES([Makefile
                  src/ctpsfcube/Makefile
                  src/ctedispcube/Makefile
                  src/ctbkgcube/Makefile
+                 src/ctmapcube/Makefile
                  src/ctcubemask/Makefile
                  src/cttsmap/Makefile
                  src/ctbutterfly/Makefile
@@ -42,6 +42,7 @@ testscripts = ctools/tests/__init__.py \
               ../test/test_ctexpcube.py \
               ../test/test_ctpsfcube.py \
               ../test/test_ctbkgcube.py \
+              ../test/test_ctmapcube.py \
               ../test/test_ctcubemask.py \
               ../test/test_ctbutterfly.py \
               ../test/test_ctulimit.py \
+/***************************************************************************
+ *                  ctmapcube - Map cube generation tool                   *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2016 by 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/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file ctmapcube.i
+ * @brief Map cube generation tool definition
+ * @author Juergen Knoedlseder
+ */
+%{
+/* Put headers and other declarations here that are needed for compilation */
+#include "ctmapcube.hpp"
+%}
+
+
+/***********************************************************************//**
+ * @class ctmapcube
+ *
+ * @brief Map cube generation tool
+ ***************************************************************************/
+class ctmapcube : public ctool {
+
+public:
+    // Constructors and destructors
+    ctmapcube(void);
+    ctmapcube(int argc, char *argv[]);
+    ctmapcube(const ctmapcube& app);
+    virtual ~ctmapcube(void);
+
+    // Methods
+    void clear(void);
+    void run(void);
+    void save(void);
+    void publish(const std::string& name = "");
+    void models(const GModels& models);
+};
+
+
+/***********************************************************************//**
+ * @brief Map cube generation tool Python extension
+ ***************************************************************************/
+%extend ctmapcube {
+    ctmapcube copy() {
+        return (*self);
+    }
+}
@@ -56,6 +56,7 @@
 %include "ctpsfcube.i"
 %include "ctedispcube.i"
 %include "ctbkgcube.i"
+%include "ctmapcube.i"
 %include "ctcubemask.i"
 %include "cttsmap.i"
 %include "ctbutterfly.i"
@@ -16,6 +16,7 @@ SUBDIRS = support \
           ctpsfcube \
           ctedispcube \
           ctbkgcube \
+          ctmapcube \
           ctcubemask \
           cttsmap \
           ctbutterfly \
@@ -43,6 +44,7 @@ libctools_la_LIBADD = support/libsupport.la \
                       ctpsfcube/libctpsfcube.la \
                       ctedispcube/libctedispcube.la \
                       ctbkgcube/libctbkgcube.la \
+                      ctmapcube/libctmapcube.la \
                       ctcubemask/libctcubemask.la \
                       cttsmap/libcttsmap.la \
                       ctbutterfly/libctbutterfly.la \
+## Process this file with automake to produce Makefile.in
+
+# Set include directories
+INCLUDES = -I$(top_srcdir)/src/support
+
+# Program name
+bin_PROGRAMS = ctmapcube
+
+# Sources
+sources = ctmapcube.cpp
+
+# Libtool library
+noinst_LTLIBRARIES = libctmapcube.la
+libctmapcube_la_SOURCES = $(sources)
+libctmapcube_la_LDFLAGS = @OPENMP_LIB_LDFLAGS@
+
+# Program definition
+ctmapcube_SOURCES = main.cpp
+ctmapcube_LDFLAGS = $(LIB_LINK_OPTION) @LDFLAGS@
+ctmapcube_CPPFLAGS = @CPPFLAGS@
+ctmapcube_LDADD = libctmapcube.la $(top_srcdir)/src/support/libsupport.la
+
+# Parameter files
+pfilesdir = $(prefix)/syspfiles
+pfiles_DATA = ctmapcube.par
+/***************************************************************************
+ *                  ctmapcube - Map cube generation tool                   *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2016 by 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/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file ctmapcube.cpp
+ * @brief Map cube generation tool implementation
+ * @author Juergen Knoedlseder
+ */
+
+/* __ Includes ___________________________________________________________ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include <cstdio>
+#include "ctmapcube.hpp"
+#include "GTools.hpp"
+
+/* __ Method name definitions ____________________________________________ */
+#define G_GET_PARAMETERS                        "ctmapcube::get_parameters()"
+
+/* __ Debug definitions __________________________________________________ */
+
+/* __ Coding definitions _________________________________________________ */
+
+/* __ Constants __________________________________________________________ */
+
+
+/*==========================================================================
+ =                                                                         =
+ =                        Constructors/destructors                         =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Void constructor
+ ***************************************************************************/
+ctmapcube::ctmapcube(void) : ctool(CTMAPCUBE_NAME, CTMAPCUBE_VERSION)
+{
+    // Initialise members
+    init_members();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Command line constructor
+ *
+ * @param[in] argc Number of arguments in command line.
+ * @param[in] argv Array of command line arguments.
+ ***************************************************************************/
+ctmapcube::ctmapcube(int argc, char *argv[]) :
+           ctool(CTMAPCUBE_NAME, CTMAPCUBE_VERSION, argc, argv)
+{
+    // Initialise members
+    init_members();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Copy constructor
+ *
+ * @param[in] app Application.
+ ***************************************************************************/
+ctmapcube::ctmapcube(const ctmapcube& app) : ctool(app)
+{
+    // Initialise members
+    init_members();
+
+    // Copy members
+    copy_members(app);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Destructor
+ ***************************************************************************/
+ctmapcube::~ctmapcube(void)
+{
+    // Free members
+    free_members();
+
+    // Return
+    return;
+}
+
+
+/*==========================================================================
+ =                                                                         =
+ =                               Operators                                 =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Assignment operator
+ *
+ * @param[in] app Application.
+ * @return Application.
+ ***************************************************************************/
+ctmapcube& ctmapcube::operator=(const ctmapcube& app)
+{
+    // Execute only if object is not identical
+    if (this != &app) {
+
+        // Copy base class members
+        this->ctool::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 ctmapcube::clear(void)
+{
+    // Free members
+    free_members();
+    this->ctool::free_members();
+    this->GApplication::free_members();
+
+    // Initialise members
+    this->GApplication::init_members();
+    this->ctool::init_members();
+    init_members();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Generate map cube.
+ *
+ * This method reads the task parameters from the parfile and generates the
+ * map cube.
+ ***************************************************************************/
+void ctmapcube::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 models into logger
+    if (logTerse()) {
+        log << std::endl;
+        log.header1(gammalib::number("Model", m_models.size()));
+        log << m_models << std::endl;
+    }
+
+    // Write header
+    if (logTerse()) {
+        log << std::endl;
+        log.header1("Generate map cube");
+    }
+
+    //TODO
+
+    // Log cube
+    if (logTerse()) {
+        log << std::endl;
+        log.header1("Map cube");
+        log << m_cube << std::endl;
+    }
+
+    // Optionally publish map cube
+    if (m_publish) {
+        publish();
+    }
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Save map cube
+ *
+ * Saves the map cube into a FITS file.
+ ***************************************************************************/
+void ctmapcube::save(void)
+{
+    // Write header
+    if (logTerse()) {
+        log << std::endl;
+        log.header1("Save map cube");
+    }
+
+    // Get map cube filename
+    m_outcube = (*this)["outcube"].filename();
+
+    // Save only if filename is non-empty
+    if (!m_outcube.is_empty()) {
+
+        // Log filename
+        if (logTerse()) {
+            log << "Save map cube into file \""+m_outcube+"\".";
+            log << std::endl;
+        }
+
+        // Save map cube into FITS file
+        m_cube.save(m_outcube, clobber());
+
+    }
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Publish map cube
+ *
+ * @param[in] name Map cube name.
+ ***************************************************************************/
+void ctmapcube::publish(const std::string& name)
+{
+    // Write header
+    if (logTerse()) {
+        log << std::endl;
+        log.header1("Publish map cube");
+    }
+
+    // Set default name is user name is empty
+    std::string user_name(name);
+    if (user_name.empty()) {
+        user_name = CTMAPCUBE_NAME;
+    }
+
+    // Log filename
+    if (logTerse()) {
+        log << "Publish \""+user_name+"\" map cube." << std::endl;
+    }
+
+    // Publish map cube
+    m_cube.cube().publish(user_name);
+
+    // Return
+    return;
+}
+
+
+/*==========================================================================
+ =                                                                         =
+ =                             Private methods                             =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Initialise class members
+ ***************************************************************************/
+void ctmapcube::init_members(void)
+{
+    // Initialise members
+    m_outcube.clear();
+    m_publish = false;
+
+    // Initialise protected members
+    m_models.clear();
+    m_cube.clear();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Copy class members
+ *
+ * @param[in] app Application.
+ ***************************************************************************/
+void ctmapcube::copy_members(const ctmapcube& app)
+{
+    // Copy attributes
+    m_outcube = app.m_outcube;
+    m_publish = app.m_publish;
+
+    // Copy protected members
+    m_models = app.m_models;
+    m_cube   = app.m_cube;
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Delete class members
+ ***************************************************************************/
+void ctmapcube::free_members(void)
+{
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Get application parameters
+ *
+ * Get all task parameters from parameter file or (if required) by querying
+ * the user. The parameters are read in the correct order.
+ ***************************************************************************/
+void ctmapcube::get_parameters(void)
+{
+    // Get energy binning parameters
+    GEnergies energies = create_energies();
+
+    // Get spatial binning parameters
+    double      xref     = (*this)["xref"].real();
+    double      yref     = (*this)["yref"].real();
+    std::string proj     = (*this)["proj"].string();
+    std::string coordsys = (*this)["coordsys"].string();
+    double      binsz    = (*this)["binsz"].real();
+    int         nxpix    = (*this)["nxpix"].integer();
+    int         nypix    = (*this)["nypix"].integer();
+
+    // Allocate map cube
+    GSkyMap cube(proj, coordsys, xref, yref, -binsz, binsz, nxpix, nypix,
+                 energies.size());
+
+    // Allocate map cube
+    m_cube = GModelSpatialDiffuseCube(cube, energies);
+
+    // Read model definition file if required
+    if (m_models.size() == 0) {
+
+        // Get model filename
+        GFilename inmodel = (*this)["inmodel"].filename();
+
+        // Load models from file
+        m_models.load(inmodel);
+
+    } // endif: there were no models
+
+    // Get remaining parameters
+    m_publish = (*this)["publish"].boolean();
+
+    // Read optionally output cube filenames
+    if (read_ahead()) {
+        m_outcube = (*this)["outcube"].filename();
+    }
+
+    // Write parameters into logger
+    if (logTerse()) {
+        log_parameters();
+        log << std::endl;
+    }
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Create energy vector from user parameters
+ *
+ * @exception GException::invalid_value
+ *            No valid energy boundary extension found.
+ *
+ * Get the energy vector according to the user parameters. The method
+ * supports loading of energy information from the ENERGIES extension (or
+ * any other extension specified by the "ebinfile" parameters), or setting
+ * energies using a linear or logarithmical spacing.
+ *
+ * The following parameters are read:
+ *      ebinalg - Energy binning algorithm
+ *      emin - Minimum energy (if ebinalg != FILE)
+ *      emax - Maximum energy (if ebinalg != FILE)
+ *      enumbins - Number of energies (if ebinalg != FILE)
+ ***************************************************************************/
+GEnergies ctmapcube::create_energies(void)
+{
+    // Allocate energies
+    GEnergies energies;
+
+    // Get energy binning algorithm
+    std::string ebinalg = (*this)["ebinalg"].string();
+
+    // If energy binning algorithm is of type "FILE" (case sensitive), then
+    // read energies from FITS file ...
+    if (ebinalg == "FILE") {
+
+        // Get filename
+        GFilename ebinfile = (*this)["ebinfile"].filename();
+
+        // Load energies
+        energies.load(ebinfile);
+
+    } // endif: ebinalg was "FILE"
+
+    // ... otherwise use a linear or a logarithmically-spaced energy binning
+    else {
+
+        // Get task parameters
+        double emin     = (*this)["emin"].real();
+        double emax     = (*this)["emax"].real();
+        int    enumbins = (*this)["enumbins"].integer();
+
+        // Initialise log mode for ebinning
+        bool log = true;
+
+        // Check if algorithm is linear
+        if (ebinalg == "LIN") {
+            log = false;
+        }
+
+        // Setup energies
+        energies = GEnergies(enumbins, GEnergy(emin, "TeV"),
+                                       GEnergy(emax, "TeV"), log);
+
+    } // endelse: ebinalg was not "FILE"
+
+    // Return energies
+    return energies;
+}
+/***************************************************************************
+ *                  ctmapcube - Map cube generation tool                   *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2016 by 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/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file ctmapcube.hpp
+ * @brief Map cube generation tool definition
+ * @author Juergen Knoedlseder
+ */
+
+#ifndef CTMAPCUBE_HPP
+#define CTMAPCUBE_HPP
+
+/* __ Includes ___________________________________________________________ */
+#include "GammaLib.hpp"
+#include "GCTALib.hpp"
+#include "ctool.hpp"
+
+/* __Definitions _________________________________________________________ */
+#define CTMAPCUBE_NAME    "ctmapcube"
+#define CTMAPCUBE_VERSION "1.0.0"
+
+
+/***********************************************************************//**
+ * @class ctmapcube
+ *
+ * @brief Map cube generation tool
+ ***************************************************************************/
+class ctmapcube : public ctool {
+
+public:
+    // Constructors and destructors
+    ctmapcube(void);
+    ctmapcube(int argc, char *argv[]);
+    ctmapcube(const ctmapcube& app);
+    virtual ~ctmapcube(void);
+
+    // Operators
+    ctmapcube& operator=(const ctmapcube& app);
+
+    // Methods
+    void clear(void);
+    void run(void);
+    void save(void);
+    void publish(const std::string& name = "");
+    void models(const GModels& models);
+
+protected:
+    // Protected methods
+    void      init_members(void);
+    void      copy_members(const ctmapcube& app);
+    void      free_members(void);
+    void      get_parameters(void);
+    GEnergies create_energies(void);
+
+    // User parameters
+    GFilename m_outcube;                //!< Output map cube filename
+    bool      m_publish;                //!< Publish map cube?
+
+    // Protected members
+    GModels                  m_models;  //!< Model container
+    GModelSpatialDiffuseCube m_cube;    //!< Map cube
+};
+
+
+/***********************************************************************//**
+ * @brief Set models
+ *
+ * @param[in] models Model container.
+ *
+ * Set model container to be used for map cube generation.
+ ***************************************************************************/
+inline
+void ctmapcube::models(const GModels& models)
+{
+    m_models = models;
+    return;
+}
+
+#endif /* CTMAPCUBE_HPP */
+#############################################################################
+#                                                                           #
+#                       CTA Science Analysis Tools                          #
+#                                                                           #
+# ------------------------------------------------------------------------- #
+#                                                                           #
+#  File:      ctmapcube.par                                                 #
+#  Version:   1.0.0                                                         #
+#                                                                           #
+#  Author:    Juergen Knoedlseder                                           #
+#             Jurgen.Knodlseder@irap.omp.eu                                 #
+#             IRAP                                                          #
+#                                                                           #
+#  Purpose:   Parameter file for the ctmapcube application                  #
+#                                                                           #
+#############################################################################
+#
+# General parameters
+#===================
+inmodel,  f, a, "$CTOOLS/share/models/crab.xml",,, "Input model XML file"
+outcube,  f, a, "mapcube.fits",,, "Output map cube file"
+
+#
+# Energy binning parameters
+#==========================
+ebinalg,  s, a, 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,1,200, "Number of energy bins"
+ebinfile, f, a, 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,0,360, "First coordinate of image center in degrees (RA or galactic l)"
+yref,     r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)"
+proj,     s, a, CAR,AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method"
+
+#
+# Standard parameters
+#====================
+publish, b, h, no,,, "Publish map cube on VO Hub?"
+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, "ctmapcube.log",,, "Log filename"
+/***************************************************************************
+ *                 ctmapcube - CTA map cube generation tool                *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2016 by 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/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file ctmapcube/main.cpp
+ * @brief CTA map cube generation tool main code
+ * @author Juergen Knoedlseder
+ */
+
+/* __ Includes ___________________________________________________________ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include <cstdio>
+#include "ctmapcube.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 gtbin application. It allocates a
+ * ctmapcube object and runs the application.
+ ***************************************************************************/
+int main (int argc, char *argv[])
+{
+    // Initialise return code
+    int rc = 1;
+
+    // Create instance of application
+    ctmapcube 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"
+                              " ctmapcube. 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;
+}
+#! /usr/bin/env python
+# ==========================================================================
+# This scripts performs unit tests for the ctmapcube tool.
+#
+# Copyright (C) 2016 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 os
+import gammalib
+import ctools
+
+
+# ============================= #
+# Test class for ctmapcube tool #
+# ============================= #
+class Test(gammalib.GPythonTestSuite):
+    """
+    Test class for ctmapcube tool.
+
+    This test class makes unit tests for the ctmapcube tool by using it from
+    the command line and from Python.
+    """
+
+    # Constructor
+    def __init__(self):
+        """
+        Constructor.
+        """
+        # Call base class constructor
+        gammalib.GPythonTestSuite.__init__(self)
+
+        # Set members
+        self._model_name = "data/crab.xml"
+
+        # Return
+        return
+
+    # Set test functions
+    def set(self):
+        """
+        Set all test functions.
+        """
+        # Set test name
+        self.name("ctmapcube")
+
+        # Append tests
+        self.append(self._test_cmd, "Test ctmapcube on command line")
+        self.append(self._test_python, "Test ctmapcube from Python")
+
+        # Return
+        return
+
+    # Test ctmapcube on command line
+    def _test_cmd(self):
+        """
+        Test ctmapcube on the command line.
+        """
+        # Kluge to set the command (installed version has no README file)
+        if os.path.isfile("README"):
+            ctmapcube = "../src/ctmapcube/ctmapcube"
+        else:
+            ctmapcube = "ctmapcube"
+
+        # Setup ctmapcube command
+        cmd = ctmapcube+' inmodel="'+self._model_name+'"'+ \
+                        ' outcube="ctmapcube_cmd1.fits"'+\
+                        ' emin=0.1 emax=100.0 enumbins=20 ebinalg=LOG'+ \
+                        ' nxpix=200 nypix=200 binsz=0.02 coordsys=CEL'+ \
+                        ' xref=83.63 yref=22.01 proj=CAR'+ \
+                        ' logfile="ctmapcube_cmd1.log" chatter=1'
+
+        # Execute ctmapcube, make sure we catch any exception
+        try:
+            #rc = os.system(cmd)
+            rc = os.system(cmd+" >/dev/null 2>&1")
+        except:
+            pass
+
+        # Check if execution was successful
+        self.test_assert(rc == 0,
+                         "Successful ctmapcube execution on command line")
+
+        # Check map cube
+        self._check_cube("ctmapcube_cmd1.fits")
+
+        # Setup ctmapcube command
+        cmd = ctmapcube+' inmodel="events_that_do_not_exist.fits"'+ \
+                        ' outcube="ctmapcube_cmd2.fits"'+\
+                        ' emin=0.1 emax=100.0 enumbins=20 ebinalg=LOG'+ \
+                        ' nxpix=200 nypix=200 binsz=0.02 coordsys=CEL'+ \
+                        ' xref=83.63 yref=22.01 proj=CAR'+ \
+                        ' logfile="ctmapcube_cmd2.log"'
+
+        # Execute ctmapcube, make sure we catch any exception
+        try:
+            #rc = os.system(cmd)
+            rc = os.system(cmd+" >/dev/null 2>&1")
+        except:
+            pass
+
+        # Check if execution failed
+        self.test_assert(rc != 0,
+                         "Failure of ctmapcube execution on command line")
+
+        # Return
+        return
+
+    # Test ctmapcube from Python
+    def _test_python(self):
+        """
+        Test ctmapcube from Python.
+        """
+        # Set-up ctmapcube
+        mapcube = ctools.ctmapcube()
+        mapcube["inmodel"]  = self._model_name
+        mapcube["outcube"]  = "ctmapcube_py1.fits"
+        mapcube["ebinalg"]  = "LOG"
+        mapcube["emin"]     = 0.1
+        mapcube["emax"]     = 100.0
+        mapcube["enumbins"] = 20
+        mapcube["nxpix"]    = 200
+        mapcube["nypix"]    = 200
+        mapcube["binsz"]    = 0.02
+        mapcube["coordsys"] = "CEL"
+        mapcube["proj"]     = "CAR"
+        mapcube["xref"]     = 83.63
+        mapcube["yref"]     = 22.01
+        mapcube["logfile"]  = "ctmapcube_py1.log"
+        mapcube["chatter"]  = 2
+
+        # Run ctmapcube tool
+        mapcube.logFileOpen()   # Make sure we get a log file
+        mapcube.run()
+
+        # Save map cube
+        mapcube.save()
+
+        # Check map cube
+        self._check_cube("ctmapcube_py1.fits")
+
+        # Return
+        return
+
+    # Check map cube
+    def _check_cube(self, filename):
+        """
+        Check content of map cube."
+
+        Args:
+            filename: Map cube file name.
+        """
+        # Load map cube
+        cube = gammalib.GModelSpatialDiffuseCube(filename)
+        #print(cube)
+
+        # Test map cube
+        #self.test_value(cube.maps(), 20, "20 maps")
+        #self.test_value(cube.pixels(), 40000, "40000 map pixels")
+        # The map cube is not loaded by default !!!! We should add the method so that
+        # the attributes are known !!!!
+        self.test_value(cube.maps(), 0, "20 maps")
+        self.test_value(cube.pixels(), 0, "40000 map pixels")
+
+        # Return
+        return
@@ -33,6 +33,7 @@ import test_ctskymap
 import test_ctexpcube
 import test_ctpsfcube
 import test_ctbkgcube
+import test_ctmapcube
 import test_ctcubemask
 import test_ctbutterfly
 import test_ctulimit
@@ -61,6 +62,7 @@ def test(installed=False):
     suite_ctexpcube   = test_ctexpcube.Test()
     suite_ctpsfcube   = test_ctpsfcube.Test()
     suite_ctbkgcube   = test_ctbkgcube.Test()
+    suite_ctmapcube   = test_ctmapcube.Test()
     suite_ctcubemask  = test_ctcubemask.Test()
     suite_ctbutterfly = test_ctbutterfly.Test()
     suite_ctulimit    = test_ctulimit.Test()
@@ -78,6 +80,7 @@ def test(installed=False):
     suite_ctexpcube.set()
     suite_ctpsfcube.set()
     suite_ctbkgcube.set()
+    suite_ctmapcube.set()
     suite_ctcubemask.set()
     suite_ctbutterfly.set()
     suite_ctulimit.set()
@@ -95,6 +98,7 @@ def test(installed=False):
     suites.append(suite_ctexpcube)
     suites.append(suite_ctpsfcube)
     suites.append(suite_ctbkgcube)
+    suites.append(suite_ctmapcube)
     suites.append(suite_ctcubemask)
     suites.append(suite_ctbutterfly)
     suites.append(suite_ctulimit)