added first version of ctbutterfly including a unit test

Michael Mayer
1 parent a2954be5
Showing 13 changed files with 990 additions and 3 deletions
configure.ac
pyext/ctbutterfly.i
pyext/ctools/ctools.i
src/Makefile.am
src/ctbutterfly/ctbutterfly.cpp
src/ctbutterfly/ctbutterfly.hpp
src/ctbutterfly/ctbutterfly.par
src/ctbutterfly/main.cpp
src/ctobssim/ctobssim.cpp
test/test_ctbutterfly.py
test/test_ctools.sh
test/test_cttsmap.py
test/test_python.py
@@ -887,6 +887,7 @@ AC_CONFIG_FILES([Makefile
                  src/ctbkgcube/Makefile
                  src/ctcubemask/Makefile
                  src/cttsmap/Makefile
+                 src/ctbutterfly/Makefile
                  src/support/Makefile
                  src/ctools-setup
                  test/Makefile
+/***************************************************************************
+ *                    ctbutterfly - butterfly calculation tool             *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2014 by Michael Mayer                                    *
+ * ----------------------------------------------------------------------- *
+ *                                                                         *
+ *  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 ctbutterfly.i
+ * @brief tool for butterfly calculations
+ * @author Michael Mayer
+ */
+%{
+/* Put headers and other declarations here that are needed for compilation */
+#include "ctbutterfly.hpp"
+%}
+
+
+/***********************************************************************//**
+ * @class ctbutterfly
+ *
+ * @brief butterfly calculation tool
+ ***************************************************************************/
+class ctbutterfly : public ctool {
+public:
+    // Constructors and destructors
+    ctbutterfly(void);
+    explicit ctbutterfly(const GObservations& obs);
+    ctbutterfly(int argc, char *argv[]);
+    ctbutterfly(const ctbutterfly& app);
+    virtual ~ctbutterfly(void);
+    
+    // Methods
+    void                 clear(void);
+    void                 run(void);
+    void                 save(void);
+    const GObservations& obs(void) const;
+    
+};
+
+
+/***********************************************************************//**
+ * @brief butterfly calculation tool Python extensions
+ ***************************************************************************/
+%extend ctbutterfly {
+    ctbutterfly copy() {
+        return (*self);
+    }
+}
@@ -57,3 +57,4 @@
 %include "ctbkgcube.i"
 %include "ctcubemask.i"
 %include "cttsmap.i"
+%include "ctbutterfly.i"
@@ -16,7 +16,8 @@ SUBDIRS = support \
           ctpsfcube \
           ctbkgcube \
           ctcubemask \
-          cttsmap
+          cttsmap \
+          ctbutterfly
 # Scripts to be distributed
 scripts = ctools-setup \
@@ -39,7 +40,8 @@ libctools_la_LIBADD = support/libsupport.la \
                       ctpsfcube/libctpsfcube.la \
                       ctbkgcube/libctbkgcube.la \
                       ctcubemask/libctcubemask.la \
-                      cttsmap/libcttsmap.la
+                      cttsmap/libcttsmap.la \
+                      ctbutterfly/libctbutterfly.la              
 # Distribute scripts
+/***************************************************************************
+ *                    ctbutterfly - butterfly calculation tool                    *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2014 by Michael Mayer                                    *
+ * ----------------------------------------------------------------------- *
+ *                                                                         *
+ *  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 ctbutterfly.hpp
+ * @brief Butterfly computation tool
+ * @author Michael Mayer
+ */
+
+/* __ Includes ___________________________________________________________ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include <cstdio>
+#include <fstream>
+#include "ctbutterfly.hpp"
+#include "GTools.hpp"
+
+/* __ Method name definitions ____________________________________________ */
+#define G_GET_PARAMETERS                          "ctbutterfly::get_parameters()"
+#define G_SAVE                                                  "ctbutterfly::save()"
+
+/* __ Debug definitions __________________________________________________ */
+
+/* __ Coding definitions _________________________________________________ */
+
+
+/*==========================================================================
+ =                                                                         =
+ =                        Constructors/destructors                         =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Void constructor
+ ***************************************************************************/
+ctbutterfly::ctbutterfly(void) : ctool(CTBUTTERFLY_NAME, CTBUTTERFLY_VERSION)
+{
+    // Initialise members
+    init_members();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Observations constructor
+ *
+ * param[in] obs Observation container.
+ *
+ * This method creates an instance of the class by copying an existing
+ * observations container.
+ ***************************************************************************/
+ctbutterfly::ctbutterfly(const GObservations& obs) :
+         ctool(CTBUTTERFLY_NAME, CTBUTTERFLY_VERSION)
+{
+    // Initialise members
+    init_members();
+
+    // Set observations
+    m_obs = obs;
+
+    // Return
+    return;
+}
+
+
+
+/***********************************************************************//**
+ * @brief Command line constructor
+ *
+ * @param[in] argc Number of arguments in command line.
+ * @param[in] argv Array of command line arguments.
+ ***************************************************************************/
+ctbutterfly::ctbutterfly(int argc, char *argv[]) :
+         ctool(CTBUTTERFLY_NAME, CTBUTTERFLY_VERSION, argc, argv)
+{
+    // Initialise members
+    init_members();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Copy constructor
+ *
+ * @param[in] app Application.
+ ***************************************************************************/
+ctbutterfly::ctbutterfly(const ctbutterfly& app) : ctool(app)
+{
+    // Initialise members
+    init_members();
+
+    // Copy members
+    copy_members(app);
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Destructor
+ ***************************************************************************/
+ctbutterfly::~ctbutterfly(void)
+{
+    // Free members
+    free_members();
+
+    // Return
+    return;
+}
+
+
+/*==========================================================================
+ =                                                                         =
+ =                               Operators                                 =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Assignment operator
+ *
+ * @param[in] app Application.
+ * @return Application.
+ ***************************************************************************/
+ctbutterfly& ctbutterfly::operator=(const ctbutterfly& 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 ctbutterfly::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 Computes the butterfly
+ *
+ * This method takes the spectral fit and its covariance matrix to compute
+ * the confidence band via Gaussian error propagation
+ ***************************************************************************/
+void ctbutterfly::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;
+    }
+
+    // Get model instance for further computations
+    GModels models = m_obs.models();
+
+    // Compute covariance matrix if necessary
+	if (m_covariance.size() == 0) {
+
+        // Write header
+        if (logTerse()) {
+            log << std::endl;
+            log.header1("Compute covariance matrix");
+        }
+
+		GOptimizerPars pars = models.pars();
+		GObservations::likelihood likelihood = m_obs.function();
+		likelihood.eval(pars);
+		GMatrixSparse* curvature = dynamic_cast<GMatrixSparse*>(likelihood.curvature());
+		m_covariance = curvature->invert();
+
+	}
+
+    // Write header
+    if (logTerse()) {
+        log << std::endl;
+        log.header1("Generate butterfly");
+    }
+
+    // Initialise dummy time to evaluate spectral model
+    GTime time = GTime();
+
+    // Initialise vector of gradients
+    GVector grad = GVector(m_obs.models().npars());
+
+    // Loop over energy bins positions
+    for (int i = 0; i < m_ebounds.size(); ++i) {
+
+        int num_gradient = 0;
+
+    	// Get the energy of current bin
+        GEnergy energy = m_ebounds.elogmean(i);
+
+    	// Header for verbose logging
+    	if (logExplicit()) {
+            std::string msg = "Computing butterfly for bin number "+gammalib::str(i)+
+                              " at "+energy.print();
+            log << std::endl;
+            log.header2(msg);
+    	}
+
+    	// Initialise model flux value
+    	double model_flux = 0.0;
+
+    	// Loop over models
+    	for (int j = 0; j < models.size(); ++j ) {
+
+    	    // Check wether model is a skymodel
+    	    GModelSky* skymodel = dynamic_cast<GModelSky*>(models[j]);
+
+    	    // Yes ...
+    	    if (skymodel != NULL) {
+
+    	        // Skip spatial models
+    	        num_gradient += skymodel->spatial()->size();
+
+    	        // Get pointer to spectral model
+    	        GModelSpectral* spectral = skymodel->spectral();
+
+    	        // Set flux value of source of interest
+    	        if (skymodel->name() == m_srcname) {
+    	               model_flux = spectral->eval_gradients(energy, time);
+    	        } // endif: model was source of interest
+
+    	        // Loop over model parameters, get gradients
+    	        // and assign them to the vector
+    	        for (int k = 0; k < spectral->size(); ++k) {
+
+    	            grad[num_gradient] =  (*spectral)[k].gradient();
+    	            num_gradient += 1;
+
+    	        }// endfor: loop over spectral parameters
+
+    	        // Skip temporal models
+    	        num_gradient += skymodel->temporal()->size();
+
+    	    } // endif: model was sky model
+
+    	    else {
+
+    	        // Skip other models (e.g. GModelData instances)
+    	        num_gradient += models[j]->size();
+
+    	    } // endelse: Model was not GModelSky
+
+    	} // endfor: Looped over models
+
+    	// Multiply covariance to the gradient vector
+    	GVector vector = m_covariance * grad;
+
+    	// Get the error from the scalar product
+    	double error = std::sqrt(grad * vector);
+
+    	// Store flux, value and energy for saving
+    	m_fluxes.push_back(model_flux);
+    	m_energies.push_back(energy.MeV());
+    	m_errors.push_back(error);
+
+    	 // Log information
+        if (logExplicit()) {
+            log << " Flux at ";
+            log << energy.print();
+            log << ": ";
+            log << model_flux;
+            log << " +- ";
+            log << error;
+            log << " ph/cm2/s/MeV";
+        }
+
+    } //endfor: loop over energy bin
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Save maps
+ *
+ * This method saves the butterfly into an ascii file
+ ***************************************************************************/
+void ctbutterfly::save(void)
+{
+    // Write header
+    if (logTerse()) {
+        log << std::endl;
+        log.header1("Save Butterfly to file");
+    }
+
+    // Get output filename
+    m_outfile = (*this)["outfile"].filename();
+
+    if (gammalib::file_exists(m_outfile)) {
+
+        if ((*this)["clobber"].boolean()) {
+            std::remove(m_outfile.c_str());
+        }
+        else {
+            std::string msg = "File \""+m_outfile+"\" already exists and clobber is set to \"false\"";
+            throw GException::file_open_error(G_SAVE,msg);
+        }
+    }
+
+    // Create  file instance
+    std::fstream file;
+
+    file.open(m_outfile, std::fstream::out | std::fstream::app);
+
+    // Write short header
+    file << "# Energy [MeV] | flux [ph/cm2/s/MeV] | e_flux [ph/cm2/s/MeV]" << std::endl;
+
+    // Write values to file
+    for (int i = 0; i < m_fluxes.size(); ++i) {
+
+        file << m_energies[i] << " " << m_fluxes[i] << " " << m_errors[i] << std::endl;
+    }
+
+    // close file
+    file.close();
+
+    // Return
+    return;
+}
+
+
+/*==========================================================================
+ =                                                                         =
+ =                             Private methods                             =
+ =                                                                         =
+ ==========================================================================*/
+
+/***********************************************************************//**
+ * @brief Initialise class members
+ ***************************************************************************/
+void ctbutterfly::init_members(void)
+{
+    // Initialise members
+    m_infile.clear();
+    m_outfile.clear();
+    m_srcname.clear();
+    m_ebounds.clear();
+
+    // Initialise protected members
+    m_obs.clear();
+    m_covariance.clear();
+    m_fluxes.clear();
+    m_energies.clear();
+    m_errors.clear();
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Copy class members
+ *
+ * @param[in] app Application.
+ ***************************************************************************/
+void ctbutterfly::copy_members(const ctbutterfly& app)
+{
+    // Copy attributes
+    m_infile   = app.m_infile;
+    m_outfile  = app.m_outfile;
+    m_srcname     = app.m_srcname;
+    m_ebounds     = app.m_ebounds;
+
+    // Copy protected members
+    m_obs        = app.m_obs;
+    m_covariance = app.m_covariance;
+    m_fluxes         = app.m_fluxes;
+    m_energies    = app.m_energies;
+    m_errors         = app.m_errors;
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Delete class members
+ ***************************************************************************/
+void ctbutterfly::free_members(void)
+{
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Get application parameters
+ *
+ * @exception GException::invalid_value
+ *            Test source not found or no RA/DEC parameters found for test
+ *            source.
+ *
+ * Get all task parameters from parameter file or (if required) by querying
+ * the user. Most parameters are only required if no observation exists so
+ * far in the observation container. In this case, a single CTA observation
+ * will be added to the container, using the definition provided in the
+ * parameter file.
+ ***************************************************************************/
+void ctbutterfly::get_parameters(void)
+{
+    // If there are no observations in container then add a single CTA
+    // observation using the parameters from the parameter file
+    if (m_obs.size() == 0) {
+
+        // Allocate CTA observation
+        GCTAObservation obs;
+
+        // Get event file name
+        std::string filename = (*this)["infile"].filename();
+
+        // Try first to open as FITS file
+        try {
+
+            // Load data
+            obs.load(filename);
+
+            // Set response
+            set_obs_response(&obs);
+
+            // Append 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);
+
+            // Check if all observations have response information. If
+            // not, get the calibration database parameters and set
+            // the response properly
+            set_response(m_obs);
+
+        } // endcatch: file was an XML file
+
+    } // endif: there was no observation in the container
+
+    // If there is are no models associated with the observations then
+    // load now the model definition
+    if (m_obs.models().size() == 0) {
+
+        // Get models XML filename
+        std::string filename = (*this)["srcmdl"].filename();
+
+        // Setup models for optimizing.
+        m_obs.models(GModels(filename));
+
+    } // endif: no models were associated with observations
+
+    // Get name of test source and check container for this name
+    m_srcname = (*this)["srcname"].string();
+    if (!m_obs.models().contains(m_srcname)) {
+        std::string msg = "Source \""+m_srcname+"\" not found in model "
+                          "container. Please add a source with that name "
+                          "or check for a possible typos.";
+    	throw GException::invalid_value(G_GET_PARAMETERS, msg);
+    }
+
+    // get binning information
+    m_ebounds  = get_ebounds();
+
+    // Get matrix file name and load if possible
+    std::string matrixfilename = (*this)["matrix"].filename();
+    if (matrixfilename != "NONE") {
+        std::string msg = "Loading of matrix from file not implemented yet. "
+                           "Use filename = \"NONE\" to induce a recomputation of "
+                           "the matrix internally.";
+        throw GException::feature_not_implemented(G_GET_PARAMETERS, msg);
+        // m_covariance.load(matrixfilename);
+    }
+
+    // Optionally read ahead parameters so that they get correctly
+    // dumped into the log file
+    if (read_ahead()) {
+        m_outfile = (*this)["outfile"].filename();
+    }
+
+    // Return
+    return;
+}
+
+/***************************************************************************
+ *                    ctbutterfly - butterfly calculation tool                    *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2014 by Michael Mayer                                    *
+ * ----------------------------------------------------------------------- *
+ *                                                                         *
+ *  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 ctbutterfly.hpp
+ * @brief Butterfly calculation tool interface definition
+ * @author Michael Mayer
+ */
+
+#ifndef CTBUTTERFLY_HPP
+#define CTBUTTERFLY_HPP
+
+/* __ Includes ___________________________________________________________ */
+#include "GammaLib.hpp"
+#include "GCTALib.hpp"
+#include "ctool.hpp"
+
+/* __Definitions _________________________________________________________ */
+#define CTBUTTERFLY_NAME    "ctbutterfly"
+#define CTBUTTERFLY_VERSION "00-01-00"
+
+
+/***********************************************************************//**
+ * @class ctbutterfly
+ *
+ * @brief butterfly calculation tool
+ *
+ * This class computes the confidence interval of a fitted spectrum
+ * (butterfly) for a given set of observations
+ *
+ * The class operates on predefined observation containers, an individual
+ * event list or an observation definition XML file.
+ *
+ * During the computation the covariance matrix of the fit is used to propagate
+ * uncertainties and their correlations through the entire energy range.
+ * The output is saved as an ascii files containing the confidence band boundaries
+ ***************************************************************************/
+class ctbutterfly : public ctool {
+
+public:
+    // Constructors and destructors
+    ctbutterfly(void);
+    explicit ctbutterfly(const GObservations& obs);
+    ctbutterfly(int argc, char *argv[]);
+    ctbutterfly(const ctbutterfly& app);
+    virtual ~ctbutterfly(void);
+
+    // Operators
+    ctbutterfly& operator=(const ctbutterfly& app);
+
+    // Methods
+    void                 clear(void);
+    void                 run(void);
+    void                 save(void);
+    const GObservations& obs(void) const;
+
+protected:
+    // Protected methods
+    void init_members(void);
+    void copy_members(const ctbutterfly& app);
+    void free_members(void);
+    void get_parameters(void);
+
+    // User parameters
+    std::string              m_infile;     //!< Input file
+    std::string              m_srcname;    //!< Name of source to compute butterfly
+    std::string              m_outfile;    //!< Output ascii file
+
+    GEbounds            m_ebounds; //!< Energy binning definition
+
+      // Protected members
+    GObservations            m_obs;        //!< Observation container
+    GMatrixSparse           m_covariance; //!< covariance matrix
+
+    std::vector<double>   m_energies; //!< energy values for storage
+    std::vector<double>   m_fluxes; //!< flux values per energy bin
+    std::vector<double>   m_errors; //!< flux errors per energy bin
+};
+
+/***********************************************************************//**
+ * @brief Return observation container
+ *
+ * @return Reference to observation container
+ ***************************************************************************/
+inline
+const GObservations& ctbutterfly::obs(void) const
+{
+    return m_obs;
+}
+
+#endif /* CTBUTTERFLY_HPP */
+#############################################################################
+#                                                                           #
+#                       CTA Science Analysis Tools                          #
+#                                                                           #
+# ------------------------------------------------------------------------- #
+#                                                                           #
+#  File:      ctbutterfly                                                   #
+#  Version:   00-01-00                                                      #
+#                                                                           #
+#  Author:    Michael Mayer                                                 #
+#             michael.mayer@physik.hu-berlin.de                             #
+#             HU Berlin                                                     #
+#                                                                           #
+#  Purpose:   Parameter file for the ctbutterfly application                #
+#                                                                           #
+#############################################################################
+#
+# General parameters
+#===================
+infile,   f, a, "events.fits",,,"Input event list, cube or observation definition file"
+srcmdl,   f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
+srcname,  s, a, "CrabNebula",,,"Source of interest"
+caldb,    s, a, "dummy",,, "Calibration database"
+irf,      s, a, "cta_dummy_irf",,, "Instrument response function"
+outfile,  f, a, "butterfly.txt",,, "Output ascii file"
+matrix,   f, h, "NONE",,, "Input covariance Matrix FITS 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, h, 100,,,"Number of energy bins"
+ebinfile, f, h, "NONE",,,"Name of the file containing the energy bin definition"
+
+#
+# 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, "ctbutterfly.log",,, "Log filename"
+/***************************************************************************
+ *                      ctbutterfly - TS map calculation tool                  *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2014 by Michael Mayer                                    *
+ * ----------------------------------------------------------------------- *
+ *                                                                         *
+ *  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 ctbutterfly/main.cpp
+ * @brief CTA butterfly computation tool
+ * @author Michael Mayer
+ */
+
+/* __ Includes ___________________________________________________________ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include <cstdio>
+#include "ctbutterfly.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 ctbutterfly application. It allocates a
+ * ctbutterfly object and runs the application.
+ ***************************************************************************/
+int main (int argc, char *argv[])
+{
+    // Initialise return code
+    int rc = 1;
+
+    // Create instance of application
+    ctbutterfly 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"
+                              " ctbutterfly. 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;
+}
@@ -791,6 +791,7 @@ void ctobssim::simulate_source(GCTAObservation* obs, const GModels&amp; models,
                                                                photons[i],
                                                                *obs,
                                                                ran);
+
                                 if (event != NULL) {
                                     // Use event only if it falls within ROI
+#! /usr/bin/env python
+# ==========================================================================
+# This scripts performs unit tests for the ctbutterfly 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 ctbutterfly tool #
+# ========================== #
+class Test(gammalib.GPythonTestSuite):
+    """
+    Test class for ctbutterfly tool.
+    """
+    # Constructor
+    def __init__(self):
+        """
+        Constructor.
+        """
+        # Call base class constructor
+        gammalib.GPythonTestSuite.__init__(self)
+
+        # Set members
+        self.events_name = "data/crab_events.fits"
+        self.model_name  = "data/crab.xml"
+        self.caldb       = "irf"
+        self.irf         = "cta_dummy_irf"
+
+        # Return
+        return
+
+    # Set test functions
+    def set(self):
+        """
+        Set all test functions.
+        """
+        # Set test name
+        self.name("ctbutterfly")
+
+        # Append tests
+        self.append(self.test_functional, "Test ctbutterfly functionality")
+
+        # Return
+        return
+
+    # Test ctbutterfly functionnality
+    def test_functional(self):
+        """
+        Test ctbutterfly functionnality.
+        """
+        # Set-up ctbutterfly
+        butterfly = ctools.ctbutterfly()
+        butterfly["infile"].filename(self.events_name)
+        butterfly["srcmdl"].filename(self.model_name)
+        butterfly["srcname"].string("Crab")
+        butterfly["caldb"].string(self.caldb)
+        butterfly["irf"].string(self.irf)
+        butterfly["emin"].real(0.1)
+        butterfly["emax"].real(100)
+        butterfly["outfile"].filename("butterfly.txt")
+        
+        # Run tool
+        self.test_try("Run ctbutterfly")
+        try:
+            butterfly.run()
+            self.test_try_success()
+        except:
+            self.test_try_failure("Exception occured in butterfly.")
+
+        # Save results
+        self.test_try("Save results")
+        try:
+            butterfly.save()
+            self.test_try_success()
+        except:
+            self.test_try_failure("Exception occured in saving results.")
+        
+        # Return
+        return
@@ -45,6 +45,7 @@ ctexpcube=../src/ctexpcube/ctexpcube
 ctpsfcube=../src/ctpsfcube/ctpsfcube
 ctbkgcube=../src/ctbkgcube/ctbkgcube
 ctcubemask=../src/ctcubemask/ctcubemask
+ctbutterfly=../src/ctbutterfly/ctbutterfly
 #
@@ -542,3 +543,27 @@ else
   exit 1
 fi
 $ECHO " ok"
+
+
+
+#
+# Test ctbutterfly
+# ===============
+$ECHO -n "Test ctbutterfly: "
+$ctbutterfly infile="data/crab_events.fits.gz" \
+srcmdl="data/crab.xml" \
+srcname="Crab" \
+outfile="butterfly.txt" \
+caldb="irf" \
+irf="cta_dummy_irf" \
+emin=0.1 \
+emax=100.0
+$ECHO -n "."
+if [ -s "butterfly.txt" ]
+then
+    $ECHO -n "."
+else
+    $ECHO " butterfly.txt file is not valid"
+    exit 1
+fi
+$ECHO " ok"
@@ -40,7 +40,7 @@ class Test(gammalib.GPythonTestSuite):
         # Set members
         self.events_name = "data/crab_events.fits"
         self.model_name  = "data/crab.xml"
-        self.caldb       = "irf"
+        self.caldb       = "dummy"
         self.irf         = "cta_dummy_irf"
         # Return
@@ -34,6 +34,7 @@ import test_ctpsfcube
 import test_ctbkgcube
 import test_ctcubemask
 import test_pipelines
+import test_ctbutterfly
 # ======================== #
@@ -59,6 +60,7 @@ if __name__ == &#39;__main__&#39;:
     suite_ctbkgcube  = test_ctbkgcube.Test()
     suite_ctcubemask = test_ctcubemask.Test()
     suite_pipelines  = test_pipelines.Test()
+    suite_ctbutterfly  = test_ctbutterfly.Test()
     # Setup unit tests
     suite_ctobssim.set()
@@ -73,6 +75,7 @@ if __name__ == &#39;__main__&#39;:
     suite_ctbkgcube.set()
     suite_ctcubemask.set()
     suite_pipelines.set()
+    suite_ctbutterfly.set()
     # Append tests to container
     suites.append(suite_ctobssim)
@@ -87,6 +90,7 @@ if __name__ == &#39;__main__&#39;:
     suites.append(suite_ctbkgcube)
     suites.append(suite_ctcubemask)
     suites.append(suite_pipelines)
+    suites.append(suite_ctbutterfly)
     # Set PFILES environment variable
     try: