Integrate improvement of ctbutterfly (#1862) and remove tool specific version strings

Jürgen Knödlseder
1 parent f3c1b389
Showing 35 changed files with 365 additions and 354 deletions
ChangeLog
NEWS
doc/source/users/reference_manual/ctbutterfly.rst
src/ctbin/ctbin.cpp
src/ctbin/ctbin.hpp
src/ctbkgcube/ctbkgcube.cpp
src/ctbkgcube/ctbkgcube.hpp
src/ctbutterfly/ctbutterfly.cpp
src/ctbutterfly/ctbutterfly.hpp
src/ctcubemask/ctcubemask.cpp
src/ctcubemask/ctcubemask.hpp
src/ctedispcube/ctedispcube.cpp
src/ctedispcube/ctedispcube.hpp
src/cterror/cterror.cpp
src/cterror/cterror.hpp
src/ctexpcube/ctexpcube.cpp
src/ctexpcube/ctexpcube.hpp
src/ctlike/ctlike.cpp
src/ctlike/ctlike.hpp
src/ctmapcube/ctmapcube.cpp
-2017-04-05
+2017-04-10
         * Version 1.3.0 released
         ========================
+        Implement Gaussian error propagation in ctbutterfly (#1862)
         Copy cscripts into bin directory instead of creating a symbolic links (#1987)
         If code originates from Git, always compile the Python bindings (#1985)
         Raise exception instead of exiting with return code in Python tests (#1983)
 New Features and Important Changes in ctools 1.3.0.dev1
-5 April 2017
+10 April 2017
 Introduction
@@ -113,7 +113,11 @@ None
 ctbutterfly - Generate butterfly diagrams
 -----------------------------------------
-None
+The tool now uses Gaussian error propagation as the standard method for the
+computation of the butterfly diagram. The envelope method that existed before
+is still available for power law spectral models as an option. The "method"
+parameter has been added to allow the selection of the computation method
+(#1862).
 ctulimit - Compute upper limit
@@ -10,7 +10,7 @@ Synopsis
 --------
 This tool calculates a butterfly diagram for a specific source according to 
-its spectral model. The butterfly diagram is the envelope of all power law
+its spectral model. The butterfly diagram is the envelope of all spectral
 models that are within a given confidence limit compatible with the data.
 The default method used for the calculation is Gaussian error propagation
 using the covariance matrix from a maximum likelihood fit. By default a
@@ -77,7 +77,7 @@ General parameters
     Performs maximum likelihood fitting of input model ignoring any provided
     covariance matrix.
-``(method = GAUSSIAN|ENVELOPE) [string]``
+``(method = "GAUSSIAN") <GAUSSIAN|ENVELOPE> [string]``
     Computation method.
 ``(confidence = 0.68) [real]``
@@ -55,7 +55,7 @@
  *
  * Constructs an empty event binning tool.
  ***************************************************************************/
-ctbin::ctbin(void) : ctobservation(CTBIN_NAME, CTBIN_VERSION)
+ctbin::ctbin(void) : ctobservation(CTBIN_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -73,7 +73,7 @@ ctbin::ctbin(void) : ctobservation(CTBIN_NAME, CTBIN_VERSION)
  * Constructs event binning tool from an observation container.
  ***************************************************************************/
 ctbin::ctbin(const GObservations& obs) :
-       ctobservation(CTBIN_NAME, CTBIN_VERSION, obs)
+       ctobservation(CTBIN_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -94,7 +94,7 @@ ctbin::ctbin(const GObservations&amp; obs) :
  * parameter setting.
  ***************************************************************************/
 ctbin::ctbin(int argc, char *argv[]) : 
-       ctobservation(CTBIN_NAME, CTBIN_VERSION, argc, argv)
+       ctobservation(CTBIN_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -31,8 +31,7 @@
 #include "ctobservation.hpp"
 /* __Definitions _________________________________________________________ */
-#define CTBIN_NAME    "ctbin"
-#define CTBIN_VERSION "1.2.0"
+#define CTBIN_NAME "ctbin"
 /***********************************************************************//**
@@ -50,7 +50,7 @@
 /***********************************************************************//**
  * @brief Void constructor
  ***************************************************************************/
-ctbkgcube::ctbkgcube(void) : ctobservation(CTBKGCUBE_NAME, CTBKGCUBE_VERSION)
+ctbkgcube::ctbkgcube(void) : ctobservation(CTBKGCUBE_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -69,7 +69,7 @@ ctbkgcube::ctbkgcube(void) : ctobservation(CTBKGCUBE_NAME, CTBKGCUBE_VERSION)
  * observations container.
  ***************************************************************************/
 ctbkgcube::ctbkgcube(const GObservations& obs) :
-		ctobservation(CTBKGCUBE_NAME, CTBKGCUBE_VERSION, obs)
+		ctobservation(CTBKGCUBE_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -87,7 +87,7 @@ ctbkgcube::ctbkgcube(const GObservations&amp; obs) :
  * @param[in] argv Array of command line arguments.
  ***************************************************************************/
 ctbkgcube::ctbkgcube(int argc, char *argv[]) :
-		ctobservation(CTBKGCUBE_NAME, CTBKGCUBE_VERSION, argc, argv)
+		ctobservation(CTBKGCUBE_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -33,8 +33,7 @@
 #include "ctobservation.hpp"
 /* __Definitions _________________________________________________________ */
-#define CTBKGCUBE_NAME    "ctbkgcube"
-#define CTBKGCUBE_VERSION "1.1.0"
+#define CTBKGCUBE_NAME "ctbkgcube"
 /***********************************************************************//**
 /***************************************************************************
  *                 ctbutterfly - butterfly calculation tool                *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2014-2016 by Michael Mayer                               *
+ *  copyright (C) 2014-2017 by Michael Mayer                               *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -52,7 +52,7 @@
 /***********************************************************************//**
  * @brief Void constructor
  ***************************************************************************/
-ctbutterfly::ctbutterfly(void) : ctlikelihood(CTBUTTERFLY_NAME, CTBUTTERFLY_VERSION)
+ctbutterfly::ctbutterfly(void) : ctlikelihood(CTBUTTERFLY_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -71,7 +71,7 @@ ctbutterfly::ctbutterfly(void) : ctlikelihood(CTBUTTERFLY_NAME, CTBUTTERFLY_VERS
  * observations container.
  ***************************************************************************/
 ctbutterfly::ctbutterfly(const GObservations& obs) :
-                     ctlikelihood(CTBUTTERFLY_NAME, CTBUTTERFLY_VERSION, obs)
+             ctlikelihood(CTBUTTERFLY_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -89,7 +89,7 @@ ctbutterfly::ctbutterfly(const GObservations&amp; obs) :
  * @param[in] argv Array of command line arguments.
  ***************************************************************************/
 ctbutterfly::ctbutterfly(int argc, char *argv[]) :
-                     ctlikelihood(CTBUTTERFLY_NAME, CTBUTTERFLY_VERSION, argc, argv)
+             ctlikelihood(CTBUTTERFLY_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -245,7 +245,7 @@ void ctbutterfly::run(void)
         // Get covariance matrix.
         GObservations::likelihood likelihood = m_obs.function();
-        if(!m_gepmode) {
+        if (m_method == "GAUSSIAN") {
             m_covariance = likelihood.curvature()->invert();
         }
         else {   
@@ -275,7 +275,7 @@ void ctbutterfly::run(void)
         likelihood.eval(pars);
         // Get covariance matrix.
-        if(!m_gepmode) {
+        if (m_method == "GAUSSIAN") {
             m_covariance = likelihood.curvature()->invert();
         }
         else {   
@@ -287,262 +287,13 @@ void ctbutterfly::run(void)
     } // endif: covariance computation needed
-    // Write header into logger
-    log_header1(TERSE, "Prepare butterfly computation");
-
-    if (!m_gepmode) {
-
-        // Find parameter indices. We do this by the memory location of the
-        // model parameters.
-        GModelSky* skymodel = dynamic_cast<GModelSky*>(models[m_srcname]);
-        GModelPar* par_prefactor = &(*skymodel)["Prefactor"];
-        GModelPar* par_index     = &(*skymodel)["Index"];
-        int        inx_prefactor = -1;
-        int        inx_index     = -1;
-        for (int i = 0; i < models.npars(); ++i) {
-            if (models.pars()[i] == par_prefactor) {
-                inx_prefactor = i;
-            }
-            if (models.pars()[i] == par_index) {
-                inx_index = i;
-            }
-        }
-
-        // Store model values
-        double  prefactor_mean  = (*skymodel)["Prefactor"].value();
-        double  index_mean      = (*skymodel)["Index"].value();
-        double  pivot_mean      = (*skymodel)["PivotEnergy"].value();
-        double  prefactor_scale = (*skymodel)["Prefactor"].scale();
-        double  index_scale     = (*skymodel)["Index"].scale();
-        GEnergy pivot(pivot_mean, "MeV");
-
-        // Write parameter indices into logger
-        log_value(NORMAL, "Prefactor", gammalib::str(prefactor_mean)+
-                " ph/cm2/s/MeV ("+gammalib::str(inx_prefactor)+")");
-        log_value(NORMAL, "Index", gammalib::str(index_mean)+
-                " ("+gammalib::str(inx_index)+")");
-        log_value(NORMAL, "Pivot energy", pivot.print());
-
-        // Extract covariance elements
-        double cov_pp = m_covariance(inx_prefactor,inx_prefactor);
-        double cov_pg = m_covariance(inx_prefactor,inx_index);
-        double cov_gg = m_covariance(inx_index,inx_index);
-
-        // Compute eigenvectors and eigenvalues
-        double  lambda1;
-        double  lambda2;
-        GVector vector1(2);
-        GVector vector2(2);
-        eigenvectors(cov_pp, cov_pg, cov_pg, cov_gg,
-                &lambda1, &lambda2, &vector1, &vector2);
-
-        // Write eigenvectors and eigenvalues into logger
-        log_value(NORMAL, "Eigenvalue 1", lambda1);
-        log_value(NORMAL, "Eigenvalue 2", lambda2);
-        log_value(NORMAL, "Eigenvector 1", vector1.print());
-        log_value(NORMAL, "Eigenvector 2", vector2.print());
-
-        // Confidence scaling, calculated from a Chi-Sqr. 
-        // distribution with 2 dof
-        double scale = std::sqrt( -2. * std::log(1.-m_confidence) );
-
-        // Write confidence level information into logger
-        log_value(NORMAL, "Confidence level", m_confidence);
-        log_value(NORMAL, "Corresponding scaling", scale);
-
-        // Compute minor and major axes of error ellipse
-        double major = scale*std::sqrt(lambda1);
-        double minor = scale*std::sqrt(lambda2);
-
-        // Write header into logger
-        log_header1(TERSE, "Generate butterfly");
-
-        // Initialise result arrays
-        m_energies.assign(m_ebounds.size(), 0.0);
-        m_intensities.assign(m_ebounds.size(), 0.0);
-        m_min_intensities.assign(m_ebounds.size(), 1e30);
-        m_max_intensities.assign(m_ebounds.size(), 0.0);
-
-        // Walk around the error ellipse
-        int    nt = 360;
-        double dt = gammalib::twopi/double(nt);
-        double t  = 0.0;
-        for (int i = 0; i < nt; ++i, t += dt) {
-
-            // Pre-compute sin and cosine
-            double sin_t = std::sin(t);
-            double cos_t = std::cos(t);
-
-            // Compute prefactor and index
-            double prefactor = (major * cos_t * vector1[0] +
-                    minor * sin_t * vector2[0]) * prefactor_scale +
-                            prefactor_mean;
-            double index     = (major * cos_t * vector1[1] +
-                    minor * sin_t * vector2[1]) * index_scale +
-                            index_mean;
-
-            // Write prefactor and index into logger
-            log_value(EXPLICIT, "Angle "+gammalib::str(t*gammalib::rad2deg),
-                    gammalib::str(prefactor)+" ph/cm2/s/MeV; "+
-                    gammalib::str(index));
-
-            // Setup power law model
-            GModelSpectralPlaw plaw(prefactor, index, pivot);
-
-            // Evaluate power law at each energy and extract minimum and
-            // maximum intensity
-            for (int k = 0; k < m_ebounds.size(); ++k) {
-
-                // Get the energy of current bin
-                GTime   time;
-                GEnergy energy = m_ebounds.elogmean(k);
-
-                // Evaluate power law
-                double intensity = plaw.eval(energy, time);
-
-                // Get extremes
-                if (intensity < m_min_intensities[k]) {
-                    m_min_intensities[k] = intensity;
-                }
-                if (intensity > m_max_intensities[k]) {
-                    m_max_intensities[k] = intensity;
-                }
-
-            } // endfor: looped over all energies
-
-        } // endfor: looped over all angles
-
-        // Setup power law model for mean prefactor and index
-        GModelSpectralPlaw plaw(prefactor_mean, index_mean, pivot);
-
-        // Compute energy and intensity vectors for mean power law
-        for (int k = 0; k < m_ebounds.size(); ++k) {
-
-            // Get the energy of current bin
-            GTime   time;
-            GEnergy energy = m_ebounds.elogmean(k);
-
-            // Evaluate power law
-            double intensity = plaw.eval(energy, time);
-
-            // Store results
-            m_energies[k]    = energy.MeV();
-            m_intensities[k] = intensity;
-
-        } // endfor: looped over energies
-
-    } //endif: use butterfly envelope for computation 
-
+    // Compute butterfly depending on selected method
+    if (m_method == "GAUSSIAN") {
+        gaussian_error_propagation(models);
+    }
     else {
-
-        // Use gaussian error propagation scheme to calculate
-        // the butterfly
-
-        // Confidence scaling
-        double scale = gammalib::erfinv(m_confidence) * gammalib::sqrt_two;
-        
-        // Write confidence level information into logger
-        log_value(NORMAL, "Confidence level", m_confidence);
-        log_value(NORMAL, "Corresponding scaling", scale);
-
-        // Write header into logger
-        log_header1(TERSE, "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) {
-
-            // Initialise number of gradients
-            int num_gradient = 0;
-
-            // Get the energy of current bin
-            GEnergy energy = m_ebounds.elogmean(i);
-
-            // 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) {
-
-                    // Set flux value of source of interest, evaluate gradients
-                    if (skymodel->name() == m_srcname) {
-
-                        // Get spectral component of model
-                        GModelSpectral *spectral = skymodel->spectral();
-
-                        model_flux = spectral->eval(energy, time, true);
-
-                        // Loop over all model parameters
-                        for (int k = 0; k < skymodel->size(); ++k) {
-
-                            // set gradient if we deal with spectral component
-                            bool parIsSpectral = false;
-
-                            for (int l = 0; l < spectral->size(); ++l) {
-
-                                if ((&((*spectral)[l])) == (&((*skymodel)[k]))) {
-                                    grad[num_gradient] = (*spectral)[l].gradient();
-                                    num_gradient++;
-                                    parIsSpectral = true;
-                                    break;
-                                }
-                            } // endfor: looped over spectral parameters
-
-
-                            if (!parIsSpectral) {
-                                // If parameter is spatial or temporal, count up
-                                num_gradient++;
-                            }
-
-                        } // endfor: looped over all model parameters
-
-
-
-                    } // endif: model was source of interest
-                    else {
-                        num_gradient += skymodel->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);
-
-            // Multiply with confidence scaling
-            error *= scale;
-
-            // Store flux, value and energy for saving
-            m_intensities.push_back(model_flux);
-            m_energies.push_back(energy.MeV());
-            m_min_intensities.push_back(model_flux-error);
-            m_max_intensities.push_back(model_flux+error);
-
-        } //endfor: loop over energy bin        
-
-    }  //endelse: use gaussian error propagation
-
+        ellipsoid_boundary(models);
+    }
     // Write header into logger
     log_header1(TERSE, "Results");
@@ -558,8 +309,6 @@ void ctbutterfly::run(void)
     } //endfor: loop over energy bin  
-
-
     // Restore energy dispersion flag for all CTA observations
     for (int i = 0; i < m_obs.size(); ++i) {
         GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
@@ -627,10 +376,10 @@ void ctbutterfly::init_members(void)
 {
     // Initialise members
     m_srcname.clear();
+    m_method.clear();
     m_confidence  = 0.68;
     m_max_iter    = 50;
     m_apply_edisp = false;
-    m_gepmode     = false;
     m_fit         = false;
     m_chatter     = static_cast<GChatter>(2);
     m_ebounds.clear();
@@ -661,10 +410,10 @@ void ctbutterfly::copy_members(const ctbutterfly&amp; app)
 {
     // Copy attributes
     m_srcname     = app.m_srcname;
+    m_method      = app.m_method;
     m_confidence  = app.m_confidence;
     m_max_iter    = app.m_max_iter;
     m_apply_edisp = app.m_apply_edisp;
-    m_gepmode     = app.m_gepmode;
     m_fit         = app.m_fit;
     m_chatter     = app.m_chatter;
     m_ebounds     = app.m_ebounds;
@@ -730,7 +479,7 @@ void ctbutterfly::get_parameters(void)
     // Get remaining parameters
     m_apply_edisp = (*this)["edisp"].boolean();
-    m_gepmode     = (*this)["method"].string() == std::string("GAUSSIAN");
+    m_method      = (*this)["method"].string();
     m_confidence  = (*this)["confidence"].real();
     m_fit         = (*this)["fit"].boolean();
     m_chatter     = static_cast<GChatter>((*this)["chatter"].integer());
@@ -753,6 +502,273 @@ void ctbutterfly::get_parameters(void)
 /***********************************************************************//**
+ * @brief Compute butterfly using the ellipsoid boundary method
+ ***************************************************************************/
+void ctbutterfly::ellipsoid_boundary(GModels& models)
+{
+    // Write header into logger
+    log_header1(TERSE, "Prepare butterfly computation");
+
+    // Find parameter indices. We do this by the memory location of the
+    // model parameters.
+    GModelSky* skymodel = dynamic_cast<GModelSky*>(models[m_srcname]);
+    GModelPar* par_prefactor = &(*skymodel)["Prefactor"];
+    GModelPar* par_index     = &(*skymodel)["Index"];
+    int        inx_prefactor = -1;
+    int        inx_index     = -1;
+    for (int i = 0; i < models.npars(); ++i) {
+        if (models.pars()[i] == par_prefactor) {
+            inx_prefactor = i;
+        }
+        if (models.pars()[i] == par_index) {
+            inx_index = i;
+        }
+    }
+
+    // Store model values
+    double  prefactor_mean  = (*skymodel)["Prefactor"].value();
+    double  index_mean      = (*skymodel)["Index"].value();
+    double  pivot_mean      = (*skymodel)["PivotEnergy"].value();
+    double  prefactor_scale = (*skymodel)["Prefactor"].scale();
+    double  index_scale     = (*skymodel)["Index"].scale();
+    GEnergy pivot(pivot_mean, "MeV");
+
+    // Write parameter indices into logger
+    log_value(NORMAL, "Prefactor", gammalib::str(prefactor_mean)+
+                                   " ph/cm2/s/MeV ("+
+                                   gammalib::str(inx_prefactor)+")");
+    log_value(NORMAL, "Index", gammalib::str(index_mean)+
+                               " ("+gammalib::str(inx_index)+")");
+    log_value(NORMAL, "Pivot energy", pivot.print());
+
+    // Extract covariance elements
+    double cov_pp = m_covariance(inx_prefactor,inx_prefactor);
+    double cov_pg = m_covariance(inx_prefactor,inx_index);
+    double cov_gg = m_covariance(inx_index,inx_index);
+
+    // Compute eigenvectors and eigenvalues
+    double  lambda1;
+    double  lambda2;
+    GVector vector1(2);
+    GVector vector2(2);
+    eigenvectors(cov_pp, cov_pg, cov_pg, cov_gg,
+                 &lambda1, &lambda2, &vector1, &vector2);
+
+    // Write eigenvectors and eigenvalues into logger
+    log_value(NORMAL, "Eigenvalue 1", lambda1);
+    log_value(NORMAL, "Eigenvalue 2", lambda2);
+    log_value(NORMAL, "Eigenvector 1", vector1.print());
+    log_value(NORMAL, "Eigenvector 2", vector2.print());
+
+    // Confidence scaling, calculated from a Chi-Sqr. 
+    // distribution with 2 dof
+    double scale = std::sqrt( -2. * std::log(1.-m_confidence) );
+
+    // Write confidence level information into logger
+    log_value(NORMAL, "Confidence level", m_confidence);
+    log_value(NORMAL, "Corresponding scaling", scale);
+
+    // Compute minor and major axes of error ellipse
+    double major = scale*std::sqrt(lambda1);
+    double minor = scale*std::sqrt(lambda2);
+
+    // Write header into logger
+    log_header1(TERSE, "Generate butterfly");
+
+    // Initialise result arrays
+    m_energies.assign(m_ebounds.size(), 0.0);
+    m_intensities.assign(m_ebounds.size(), 0.0);
+    m_min_intensities.assign(m_ebounds.size(), 1e30);
+    m_max_intensities.assign(m_ebounds.size(), 0.0);
+
+    // Walk around the error ellipse
+    int    nt = 360;
+    double dt = gammalib::twopi/double(nt);
+    double t  = 0.0;
+    for (int i = 0; i < nt; ++i, t += dt) {
+
+        // Pre-compute sin and cosine
+        double sin_t = std::sin(t);
+        double cos_t = std::cos(t);
+
+        // Compute prefactor and index
+        double prefactor = (major * cos_t * vector1[0] +
+                           minor * sin_t * vector2[0]) * prefactor_scale +
+                           prefactor_mean;
+        double index     = (major * cos_t * vector1[1] +
+                           minor * sin_t * vector2[1]) * index_scale +
+                           index_mean;
+
+        // Write prefactor and index into logger
+        log_value(EXPLICIT, "Angle "+gammalib::str(t*gammalib::rad2deg),
+                gammalib::str(prefactor)+" ph/cm2/s/MeV; "+
+                gammalib::str(index));
+
+        // Setup power law model
+        GModelSpectralPlaw plaw(prefactor, index, pivot);
+
+        // Evaluate power law at each energy and extract minimum and
+        // maximum intensity
+        for (int k = 0; k < m_ebounds.size(); ++k) {
+
+            // Get the energy of current bin
+            GTime   time;
+            GEnergy energy = m_ebounds.elogmean(k);
+
+            // Evaluate power law
+            double intensity = plaw.eval(energy, time);
+
+            // Get extremes
+            if (intensity < m_min_intensities[k]) {
+                m_min_intensities[k] = intensity;
+            }
+            if (intensity > m_max_intensities[k]) {
+                m_max_intensities[k] = intensity;
+            }
+
+        } // endfor: looped over all energies
+
+    } // endfor: looped over all angles
+
+    // Setup power law model for mean prefactor and index
+    GModelSpectralPlaw plaw(prefactor_mean, index_mean, pivot);
+
+    // Compute energy and intensity vectors for mean power law
+    for (int k = 0; k < m_ebounds.size(); ++k) {
+
+        // Get the energy of current bin
+        GTime   time;
+        GEnergy energy = m_ebounds.elogmean(k);
+
+        // Evaluate power law
+        double intensity = plaw.eval(energy, time);
+
+        // Store results
+        m_energies[k]    = energy.MeV();
+        m_intensities[k] = intensity;
+
+    } // endfor: looped over energies
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Compute butterfly using Gaussian error propagation
+ *
+ * Computes the butterfly diagram using Gaussian error propagation. This
+ * works for any kind of model.
+ ***************************************************************************/
+void ctbutterfly::gaussian_error_propagation(GModels& models)
+{
+    // Write header into logger
+    log_header1(TERSE, "Prepare butterfly computation");
+
+    // Confidence scaling
+    double scale = gammalib::erfinv(m_confidence) * gammalib::sqrt_two;
+    
+    // Write confidence level information into logger
+    log_value(NORMAL, "Confidence level", m_confidence);
+    log_value(NORMAL, "Corresponding scaling", scale);
+
+    // Write header into logger
+    log_header1(TERSE, "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
+    for (int i = 0; i < m_ebounds.size(); ++i) {
+
+        // Initialise number of gradients
+        int num_gradient = 0;
+
+        // Get the mean logarithmic energy of current bin
+        GEnergy energy = m_ebounds.elogmean(i);
+
+        // 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) {
+
+                // Set flux value of source of interest, evaluate gradients
+                if (skymodel->name() == m_srcname) {
+
+                    // Get spectral component of model
+                    GModelSpectral *spectral = skymodel->spectral();
+
+                    // Get model flux
+                    model_flux = spectral->eval(energy, time, true);
+
+                    // Loop over all model parameters
+                    for (int k = 0; k < skymodel->size(); ++k) {
+
+                        // Set gradient if we deal with spectral component
+                        bool parIsSpectral = false;
+                        for (int l = 0; l < spectral->size(); ++l) {
+                            if ((&((*spectral)[l])) == (&((*skymodel)[k]))) {
+                                grad[num_gradient] = (*spectral)[l].gradient();
+                                num_gradient++;
+                                parIsSpectral = true;
+                                break;
+                            }
+                        } // endfor: looped over spectral parameters
+
+                        // If parameter is spatial or temporal, count up
+                        if (!parIsSpectral) {
+                            num_gradient++;
+                        }
+
+                    } // endfor: looped over all model parameters
+
+                } // endif: model was source of interest
+                else {
+                    num_gradient += skymodel->size();
+                }
+
+            } // endif: model was sky model
+
+            // Skip other models (e.g. GModelData instances)
+            else {
+                num_gradient += models[j]->size();
+            }
+
+        } // 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);
+
+        // Multiply with confidence scaling
+        error *= scale;
+
+        // Store flux, value and energy for saving
+        m_intensities.push_back(model_flux);
+        m_energies.push_back(energy.MeV());
+        m_min_intensities.push_back(model_flux-error);
+        m_max_intensities.push_back(model_flux+error);
+
+    } // endfor: Looped over energy bins
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Check if sky model is valid
  *
  * @exception GException::invalid_value
@@ -776,13 +792,15 @@ void ctbutterfly::check_model(void)
             throw GException::invalid_value(G_CHECK_MODEL, msg);
         }
-        // Check that the spectral model is a power law for default mode
-        if(!m_gepmode){
+        // Check that the spectral model is a power law for the "ENVELOPE"
+        // method
+        if (m_method == "ENVELOPE") {
             if (model->spectral()->type() != "PowerLaw") {
                 std::string msg = "\""+model->spectral()->type()+"\" cannot be "
-                        "used as spectral model for an butterfly "
-                        "computation with method=ENVELOPE. Please specify a "
-                        "power law model or switch to method GAUSSIAN.";
+                                  "used as spectral model for an butterfly "
+                                  "computation with method=ENVELOPE. Please "
+                                  "specify a power law model or switch to "
+                                  "method=GAUSSIAN.";
                 throw GException::invalid_value(G_CHECK_MODEL, msg);
             }
         }
 /***************************************************************************
  *                ctbutterfly - butterfly calculation tool                 *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2014-2016 by Michael Mayer                               *
+ *  copyright (C) 2014-2017 by Michael Mayer                               *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -33,8 +33,7 @@
 #include "ctlikelihood.hpp"
 /* __Definitions _________________________________________________________ */
-#define CTBUTTERFLY_NAME    "ctbutterfly"
-#define CTBUTTERFLY_VERSION "1.1.0"
+#define CTBUTTERFLY_NAME "ctbutterfly"
 /***********************************************************************//**
@@ -77,6 +76,10 @@ protected:
     void copy_members(const ctbutterfly& app);
     void free_members(void);
     void get_parameters(void);
+    void gaussian_error_propagation(GModels& models);
+    void ellipsoid_boundary(GModels& models);
+
+
     void check_model(void);
     void eigenvectors(const double& a,
                       const double& b,
@@ -89,10 +92,10 @@ protected:
     // User parameters
     std::string m_srcname;      //!< Name of source to compute butterfly
+    std::string m_method;       //!< Computation method
     double      m_confidence;   //!< Confidence level
     int         m_max_iter;     //!< Maximum number of iterations
     bool        m_apply_edisp;  //!< Apply energy dispersion?
-    bool        m_gepmode;      //!< Switch to gaussian error propagation scheme
     bool        m_fit;          //!< Do fit?
     GEbounds    m_ebounds;      //!< Energy binning definition
     GFilename   m_outfile;      //!< Output ASCII file name
@@ -51,7 +51,7 @@
 /***********************************************************************//**
  * @brief Void constructor
  ***************************************************************************/
-ctcubemask::ctcubemask(void) : ctobservation(CTCUBEMASK_NAME, CTCUBEMASK_VERSION)
+ctcubemask::ctcubemask(void) : ctobservation(CTCUBEMASK_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -70,7 +70,7 @@ ctcubemask::ctcubemask(void) : ctobservation(CTCUBEMASK_NAME, CTCUBEMASK_VERSION
  * an observation container.
  ***************************************************************************/
 ctcubemask::ctcubemask(const GObservations& obs) :
-            ctobservation(CTCUBEMASK_NAME, CTCUBEMASK_VERSION, obs)
+            ctobservation(CTCUBEMASK_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -89,7 +89,7 @@ ctcubemask::ctcubemask(const GObservations&amp; obs) :
  * @param[in] argv Array of command line arguments.
  ***************************************************************************/
 ctcubemask::ctcubemask(int argc, char *argv[]) : 
-            ctobservation(CTCUBEMASK_NAME, CTCUBEMASK_VERSION, argc, argv)
+            ctobservation(CTCUBEMASK_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -35,8 +35,7 @@
 #include "ctobservation.hpp"
 /* __Definitions _________________________________________________________ */
-#define CTCUBEMASK_NAME    "ctcubemask"
-#define CTCUBEMASK_VERSION "1.2.0"
+#define CTCUBEMASK_NAME "ctcubemask"
 /***********************************************************************//**
@@ -52,7 +52,7 @@
  *
  * Constructs an empty energy dispersion tool.
  ***************************************************************************/
-ctedispcube::ctedispcube(void) : ctobservation(CTEDISPCUBE_NAME, CTEDISPCUBE_VERSION)
+ctedispcube::ctedispcube(void) : ctobservation(CTEDISPCUBE_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -71,7 +71,7 @@ ctedispcube::ctedispcube(void) : ctobservation(CTEDISPCUBE_NAME, CTEDISPCUBE_VER
  * in an observation container @p obs.
  ***************************************************************************/
 ctedispcube::ctedispcube(const GObservations& obs) :
-             ctobservation(CTEDISPCUBE_NAME, CTEDISPCUBE_VERSION, obs)
+             ctobservation(CTEDISPCUBE_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -92,7 +92,7 @@ ctedispcube::ctedispcube(const GObservations&amp; obs) :
  * on the command line.
  ***************************************************************************/
 ctedispcube::ctedispcube(int argc, char *argv[]) :
-             ctobservation(CTEDISPCUBE_NAME, CTEDISPCUBE_VERSION, argc, argv)
+             ctobservation(CTEDISPCUBE_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -33,8 +33,7 @@
 #include "ctobservation.hpp"
 /* __Definitions _________________________________________________________ */
-#define CTEDISPCUBE_NAME    "ctedispcube"
-#define CTEDISPCUBE_VERSION "1.1.0"
+#define CTEDISPCUBE_NAME "ctedispcube"
 /***********************************************************************//**
@@ -52,7 +52,7 @@
  *
  * Constructs an empty cterror tool.
  ***************************************************************************/
-cterror::cterror(void) : ctlikelihood(CTERROR_NAME, CTERROR_VERSION)
+cterror::cterror(void) : ctlikelihood(CTERROR_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -70,7 +70,7 @@ cterror::cterror(void) : ctlikelihood(CTERROR_NAME, CTERROR_VERSION)
  * Constructs cterror tool from an observations container.
  ***************************************************************************/
 cterror::cterror(const GObservations& obs) :
-         ctlikelihood(CTERROR_NAME, CTERROR_VERSION, obs)
+         ctlikelihood(CTERROR_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -90,7 +90,7 @@ cterror::cterror(const GObservations&amp; obs) :
  * parameters that are provided as command line arguments.
  ***************************************************************************/
 cterror::cterror(int argc, char *argv[]) :
-         ctlikelihood(CTERROR_NAME, CTERROR_VERSION, argc, argv)
+         ctlikelihood(CTERROR_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -31,8 +31,7 @@
 #include "ctlikelihood.hpp"
 /* __Definitions _________________________________________________________ */
-#define CTERROR_NAME    "cterror"
-#define CTERROR_VERSION "1.2.0"
+#define CTERROR_NAME "cterror"
 /***********************************************************************//**
@@ -50,7 +50,7 @@
 /***********************************************************************//**
  * @brief Void constructor
  ***************************************************************************/
-ctexpcube::ctexpcube(void) : ctobservation(CTEXPCUBE_NAME, CTEXPCUBE_VERSION)
+ctexpcube::ctexpcube(void) : ctobservation(CTEXPCUBE_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -69,7 +69,7 @@ ctexpcube::ctexpcube(void) : ctobservation(CTEXPCUBE_NAME, CTEXPCUBE_VERSION)
  * observations container.
  ***************************************************************************/
 ctexpcube::ctexpcube(const GObservations& obs) :
-           ctobservation(CTEXPCUBE_NAME, CTEXPCUBE_VERSION, obs)
+           ctobservation(CTEXPCUBE_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -87,7 +87,7 @@ ctexpcube::ctexpcube(const GObservations&amp; obs) :
  * @param[in] argv Array of command line arguments.
  ***************************************************************************/
 ctexpcube::ctexpcube(int argc, char *argv[]) :
-           ctobservation(CTEXPCUBE_NAME, CTEXPCUBE_VERSION, argc, argv)
+           ctobservation(CTEXPCUBE_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -33,8 +33,7 @@
 #include "ctobservation.hpp"
 /* __Definitions _________________________________________________________ */
-#define CTEXPCUBE_NAME    "ctexpcube"
-#define CTEXPCUBE_VERSION "1.1.0"
+#define CTEXPCUBE_NAME "ctexpcube"
 /***********************************************************************//**
@@ -48,7 +48,7 @@
 /***********************************************************************//**
  * @brief Void constructor
  ***************************************************************************/
-ctlike::ctlike(void) : ctlikelihood(CTLIKE_NAME, CTLIKE_VERSION)
+ctlike::ctlike(void) : ctlikelihood(CTLIKE_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -66,7 +66,7 @@ ctlike::ctlike(void) : ctlikelihood(CTLIKE_NAME, CTLIKE_VERSION)
  * Constructs ctlike tool from an observations container.
  ***************************************************************************/
 ctlike::ctlike(const GObservations& obs) :
-        ctlikelihood(CTLIKE_NAME, CTLIKE_VERSION, obs)
+        ctlikelihood(CTLIKE_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -86,7 +86,7 @@ ctlike::ctlike(const GObservations&amp; obs) :
  * that are provided as command line arguments.
  ***************************************************************************/
 ctlike::ctlike(int argc, char *argv[]) : 
-        ctlikelihood(CTLIKE_NAME, CTLIKE_VERSION, argc, argv)
+        ctlikelihood(CTLIKE_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -31,8 +31,7 @@
 #include "ctlikelihood.hpp"
 /* __Definitions _________________________________________________________ */
-#define CTLIKE_NAME    "ctlike"
-#define CTLIKE_VERSION "1.2.0"
+#define CTLIKE_NAME "ctlike"
 /***********************************************************************//**
@@ -50,7 +50,7 @@
 /***********************************************************************//**
  * @brief Void constructor
  ***************************************************************************/
-ctmapcube::ctmapcube(void) : ctool(CTMAPCUBE_NAME, CTMAPCUBE_VERSION)
+ctmapcube::ctmapcube(void) : ctool(CTMAPCUBE_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -67,7 +67,7 @@ ctmapcube::ctmapcube(void) : ctool(CTMAPCUBE_NAME, CTMAPCUBE_VERSION)
  * @param[in] argv Array of command line arguments.
  ***************************************************************************/
 ctmapcube::ctmapcube(int argc, char *argv[]) :
-           ctool(CTMAPCUBE_NAME, CTMAPCUBE_VERSION, argc, argv)
+           ctool(CTMAPCUBE_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -33,8 +33,7 @@
 #include "ctool.hpp"
 /* __Definitions _________________________________________________________ */
-#define CTMAPCUBE_NAME    "ctmapcube"
-#define CTMAPCUBE_VERSION "1.3.0"
+#define CTMAPCUBE_NAME "ctmapcube"
 /***********************************************************************//**
@@ -53,7 +53,7 @@ const GEnergy g_energy_margin(1.0e-12, &quot;TeV&quot;);
 /***********************************************************************//**
  * @brief Void constructor
  ***************************************************************************/
-ctmodel::ctmodel(void) : ctobservation(CTMODEL_NAME, CTMODEL_VERSION)
+ctmodel::ctmodel(void) : ctobservation(CTMODEL_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -72,7 +72,7 @@ ctmodel::ctmodel(void) : ctobservation(CTMODEL_NAME, CTMODEL_VERSION)
  * observations container.
  ***************************************************************************/
 ctmodel::ctmodel(const GObservations& obs) :
-         ctobservation(CTMODEL_NAME, CTMODEL_VERSION, obs)
+         ctobservation(CTMODEL_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -90,7 +90,7 @@ ctmodel::ctmodel(const GObservations&amp; obs) :
  * @param[in] argv Array of command line arguments.
  ***************************************************************************/
 ctmodel::ctmodel(int argc, char *argv[]) :
-         ctobservation(CTMODEL_NAME, CTMODEL_VERSION, argc, argv)
+         ctobservation(CTMODEL_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -33,8 +33,7 @@
 #include "ctobservation.hpp"
 /* __Definitions _________________________________________________________ */
-#define CTMODEL_NAME    "ctmodel"
-#define CTMODEL_VERSION "1.1.0"
+#define CTMODEL_NAME "ctmodel"
 /***********************************************************************//**
@@ -67,7 +67,7 @@ const double g_roi_margin = 0.5;      //!&lt; Simulation radius margin (degrees)
  *
  * Constructs an empty ctobssim tool.
  ***************************************************************************/
-ctobssim::ctobssim(void) : ctobservation(CTOBSSIM_NAME, CTOBSSIM_VERSION)
+ctobssim::ctobssim(void) : ctobservation(CTOBSSIM_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -85,7 +85,7 @@ ctobssim::ctobssim(void) : ctobservation(CTOBSSIM_NAME, CTOBSSIM_VERSION)
  * Constructs ctobssim tool from an observation container.
  ***************************************************************************/
 ctobssim::ctobssim(const GObservations& obs) :
-          ctobservation(CTOBSSIM_NAME, CTOBSSIM_VERSION, obs)
+          ctobservation(CTOBSSIM_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -105,7 +105,7 @@ ctobssim::ctobssim(const GObservations&amp; obs) :
  * setting.
  ***************************************************************************/
 ctobssim::ctobssim(int argc, char *argv[]) :
-          ctobservation(CTOBSSIM_NAME, CTOBSSIM_VERSION, argc, argv)
+          ctobservation(CTOBSSIM_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -33,8 +33,7 @@
 #include "ctobservation.hpp"
 /* __Definitions _________________________________________________________ */
-#define CTOBSSIM_NAME    "ctobssim"
-#define CTOBSSIM_VERSION "1.2.0"
+#define CTOBSSIM_NAME "ctobssim"
 /***********************************************************************//**
@@ -49,7 +49,7 @@
 /***********************************************************************//**
  * @brief Void constructor
  ***************************************************************************/
-ctpsfcube::ctpsfcube(void) : ctobservation(CTPSFCUBE_NAME, CTPSFCUBE_VERSION)
+ctpsfcube::ctpsfcube(void) : ctobservation(CTPSFCUBE_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -68,7 +68,7 @@ ctpsfcube::ctpsfcube(void) : ctobservation(CTPSFCUBE_NAME, CTPSFCUBE_VERSION)
  * observations container.
  ***************************************************************************/
 ctpsfcube::ctpsfcube(const GObservations& obs) :
-           ctobservation(CTPSFCUBE_NAME, CTPSFCUBE_VERSION, obs)
+           ctobservation(CTPSFCUBE_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -86,7 +86,7 @@ ctpsfcube::ctpsfcube(const GObservations&amp; obs) :
  * @param[in] argv Array of command line arguments.
  ***************************************************************************/
 ctpsfcube::ctpsfcube(int argc, char *argv[]) :
-           ctobservation(CTPSFCUBE_NAME, CTPSFCUBE_VERSION, argc, argv)
+           ctobservation(CTPSFCUBE_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -33,8 +33,7 @@
 #include "ctobservation.hpp"
 /* __Definitions _________________________________________________________ */
-#define CTPSFCUBE_NAME    "ctpsfcube"
-#define CTPSFCUBE_VERSION "1.1.0"
+#define CTPSFCUBE_NAME "ctpsfcube"
 /***********************************************************************//**
@@ -53,7 +53,7 @@
 /***********************************************************************//**
  * @brief Void constructor
  ***************************************************************************/
-ctselect::ctselect(void) : ctobservation(CTSELECT_NAME, CTSELECT_VERSION)
+ctselect::ctselect(void) : ctobservation(CTSELECT_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -72,7 +72,7 @@ ctselect::ctselect(void) : ctobservation(CTSELECT_NAME, CTSELECT_VERSION)
  * provided in an observation container.
  ***************************************************************************/
 ctselect::ctselect(const GObservations& obs) :
-          ctobservation(CTSELECT_NAME, CTSELECT_VERSION, obs)
+          ctobservation(CTSELECT_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -91,7 +91,7 @@ ctselect::ctselect(const GObservations&amp; obs) :
  * @param[in] argv Array of command line arguments.
  ***************************************************************************/
 ctselect::ctselect(int argc, char *argv[]) : 
-          ctobservation(CTSELECT_NAME, CTSELECT_VERSION, argc, argv)
+          ctobservation(CTSELECT_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -35,8 +35,7 @@
 #include "ctobservation.hpp"
 /* __Definitions _________________________________________________________ */
-#define CTSELECT_NAME    "ctselect"
-#define CTSELECT_VERSION "1.2.0"
+#define CTSELECT_NAME "ctselect"
 /***********************************************************************//**
@@ -52,7 +52,7 @@
  *
  * Constructs an empty sky mapping tool.
  ***************************************************************************/
-ctskymap::ctskymap(void) : ctobservation(CTSKYMAP_NAME, CTSKYMAP_VERSION)
+ctskymap::ctskymap(void) : ctobservation(CTSKYMAP_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -70,7 +70,7 @@ ctskymap::ctskymap(void) : ctobservation(CTSKYMAP_NAME, CTSKYMAP_VERSION)
  * Constructs sky mapping tool from an observation container.
  ***************************************************************************/
 ctskymap::ctskymap(const GObservations& obs) :
-          ctobservation(CTSKYMAP_NAME, CTSKYMAP_VERSION, obs)
+          ctobservation(CTSKYMAP_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -91,7 +91,7 @@ ctskymap::ctskymap(const GObservations&amp; obs) :
  * parameter setting.
  ***************************************************************************/
 ctskymap::ctskymap(int argc, char *argv[]) : 
-          ctobservation(CTSKYMAP_NAME, CTSKYMAP_VERSION, argc, argv)
+          ctobservation(CTSKYMAP_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -32,8 +32,7 @@
 #include "ctobservation.hpp"
 /* __Definitions _________________________________________________________ */
-#define CTSKYMAP_NAME    "ctskymap"
-#define CTSKYMAP_VERSION "1.3.0"
+#define CTSKYMAP_NAME "ctskymap"
 /***********************************************************************//**
@@ -51,7 +51,7 @@
  *
  * Constructs an empty cttsmap tool.
  ***************************************************************************/
-cttsmap::cttsmap(void) : ctlikelihood(CTTSMAP_NAME, CTTSMAP_VERSION)
+cttsmap::cttsmap(void) : ctlikelihood(CTTSMAP_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -69,7 +69,7 @@ cttsmap::cttsmap(void) : ctlikelihood(CTTSMAP_NAME, CTTSMAP_VERSION)
  * Constructs cttsmap tool from an observations container.
  ***************************************************************************/
 cttsmap::cttsmap(const GObservations& obs) :
-         ctlikelihood(CTTSMAP_NAME, CTTSMAP_VERSION, obs)
+         ctlikelihood(CTTSMAP_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -89,7 +89,7 @@ cttsmap::cttsmap(const GObservations&amp; obs) :
  * that are provided as command line arguments.
  ***************************************************************************/
 cttsmap::cttsmap(int argc, char *argv[]) :
-         ctlikelihood(CTTSMAP_NAME, CTTSMAP_VERSION, argc, argv)
+         ctlikelihood(CTTSMAP_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -31,8 +31,7 @@
 #include "ctlikelihood.hpp"
 /* __Definitions _________________________________________________________ */
-#define CTTSMAP_NAME    "cttsmap"
-#define CTTSMAP_VERSION "1.2.0"
+#define CTTSMAP_NAME "cttsmap"
 /***********************************************************************//**
@@ -53,7 +53,7 @@
  *
  * Constructs an empty ctulimit tool.
  ***************************************************************************/
-ctulimit::ctulimit(void) : ctlikelihood(CTULIMIT_NAME, CTULIMIT_VERSION)
+ctulimit::ctulimit(void) : ctlikelihood(CTULIMIT_NAME, VERSION)
 {
     // Initialise members
     init_members();
@@ -71,7 +71,7 @@ ctulimit::ctulimit(void) : ctlikelihood(CTULIMIT_NAME, CTULIMIT_VERSION)
  * Constructs ctulimit tool from an observations container.
  ***************************************************************************/
 ctulimit::ctulimit(const GObservations& obs) :
-          ctlikelihood(CTULIMIT_NAME, CTULIMIT_VERSION, obs)
+          ctlikelihood(CTULIMIT_NAME, VERSION, obs)
 {
     // Initialise members
     init_members();
@@ -92,7 +92,7 @@ ctulimit::ctulimit(const GObservations&amp; obs) :
  * parameters that are provided as command line arguments.
  ***************************************************************************/
 ctulimit::ctulimit(int argc, char *argv[]) :
-          ctlikelihood(CTULIMIT_NAME, CTULIMIT_VERSION, argc, argv)
+          ctlikelihood(CTULIMIT_NAME, VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -31,8 +31,7 @@
 #include "ctlikelihood.hpp"
 /* __Definitions _________________________________________________________ */
-#define CTULIMIT_NAME    "ctulimit"
-#define CTULIMIT_VERSION "1.2.0"
+#define CTULIMIT_NAME "ctulimit"
 /***********************************************************************//**