Merge branch 'pull-request-1862_butterfly' into integration

# Conflicts: # ChangeLog # NEWS

Merge branch 'pull-request-1862_butterfly' into integration
# Conflicts: # ChangeLog # NEWS
Jürgen Knödlseder
2 parents 13844fca 53181887
Showing 36 changed files with 422 additions and 273 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/ctbutterfly/ctbutterfly.par
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
-2017-04-06
+2017-04-10
         * Version 1.3.0 released
         ========================
         Use correct ROOT histograms in csroot2caldb
         Implement regression scheme in cssens and remove num_avg parameter (#1989)
+        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
-6 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
@@ -3,21 +3,25 @@
 ctbutterfly
 ===========
-Computes butterfly diagram for a power law model.
+Computes butterfly diagram for a given spectral model.
 Synopsis
 --------
-This tool calculates a butterfly diagram for a specific source with power law
-spectral model. The butterfly diagram is the envelope of all power law models
-that are within a given confidence limit compatible with the data. By default
-a confidence level of 68% is used, but the level can be adjusted using the 
-hidden ``confidence`` parameter. ctbutterfly computes this envelope by
-evaluating for each energy the minimum and maximum intensity of all power law
-models that fall within the error ellipse of the prefactor and index parameters.
-The error ellipse is derived from the covariance matrix of a maximum likelihood
-fit.
+This tool calculates a butterfly diagram for a specific source according to 
+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
+confidence level of 68% is used, but the level can be adjusted using the
+hidden ``confidence`` parameter. For power law models, an alternative
+calculation method can be specified by setting the hidden parameter
+``method=ENVELOPE``. By using this method, the envelope is computed by
+evaluating for each energy the minimum and maximum intensity of all power
+law models that fall within the error ellipse of the prefactor and index
+parameters. The error ellipse is derived from the covariance matrix of a
+maximum likelihood fit.
 ctbutterfly assumes that the input model (parameter ``inmodel``) has been 
 adjusted using :doc:`ctlike` to the data, but if this is not the case you 
@@ -68,11 +72,14 @@ General parameters
 ``outfile [file]``
     Output butterfly ASCII file name.
-
+    
 ``(fit = no) [boolean]``
     Performs maximum likelihood fitting of input model ignoring any provided
     covariance matrix.
- 	 	 
+
+``(method = "GAUSSIAN") <GAUSSIAN|ENVELOPE> [string]``
+    Computation method.
+    
 ``(confidence = 0.68) [real]``
     Confidence level for error computation.
@@ -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();
@@ -243,11 +243,14 @@ void ctbutterfly::run(void)
         m_obs.optimize(m_opt);
         m_obs.errors(m_opt);
-        // Get covariance matrix by inverting the curvature matrix. Note that
-        // the curvature matrix is for the unscaled parameters, hence we have
-        // to scale them before we can proceed.
+        // Get covariance matrix.
         GObservations::likelihood likelihood = m_obs.function();
-        m_covariance = likelihood.curvature()->invert();
+        if (m_method == "GAUSSIAN") {
+            m_covariance = likelihood.curvature()->invert();
+        }
+        else {   
+            m_covariance = likelihood.covariance();
+        }
         // Write optimizer, optimised model and covariance matrix into logger
         log_string(NORMAL, m_opt.print(m_chatter));
@@ -271,157 +274,26 @@ void ctbutterfly::run(void)
         GObservations::likelihood likelihood = m_obs.function();
         likelihood.eval(pars);
-        // Get covariance matrix by inverting the curvature matrix. Note that
-        // the curvature matrix is for the unscaled parameters, hence we have
-        // to scale them before we can proceed.
-        m_covariance = likelihood.curvature()->invert();
+        // Get covariance matrix.
+        if (m_method == "GAUSSIAN") {
+            m_covariance = likelihood.curvature()->invert();
+        }
+        else {   
+            m_covariance = likelihood.covariance();
+        }
         // Write covariance matrix into logger
         log_string(EXPLICIT, m_covariance.print(m_chatter));
     } // endif: covariance computation needed
-    // 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;
-        }
+    // Compute butterfly depending on selected method
+    if (m_method == "GAUSSIAN") {
+        gaussian_error_propagation(models);
+    }
+    else {
+        ellipsoid_boundary(models);
     }
-
-    // Store model values and scales
-    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
-    double sigma = gammalib::erfinv(m_confidence) * gammalib::sqrt_two;
-    double scale = (sigma*sigma);
-
-    // 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
     // Write header into logger
     log_header1(TERSE, "Results");
@@ -431,11 +303,11 @@ void ctbutterfly::run(void)
         // Write results into logger
         log_value(NORMAL, "Intensity at "+gammalib::str(m_energies[k])+" MeV",
-                  gammalib::str(m_intensities[k])+" ("+
-                  gammalib::str(m_min_intensities[k])+", "+
-                  gammalib::str(m_max_intensities[k])+") ph/cm2/s/MeV");
+                gammalib::str(m_intensities[k])+" ("+
+                gammalib::str(m_min_intensities[k])+", "+
+                gammalib::str(m_max_intensities[k])+") ph/cm2/s/MeV");
-    } // endfor: looped over energies
+    } //endfor: loop over energy bin  
     // Restore energy dispersion flag for all CTA observations
     for (int i = 0; i < m_obs.size(); ++i) {
@@ -504,6 +376,7 @@ 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;
@@ -537,6 +410,7 @@ 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;
@@ -544,7 +418,7 @@ void ctbutterfly::copy_members(const ctbutterfly&amp; app)
     m_chatter     = app.m_chatter;
     m_ebounds     = app.m_ebounds;
     m_outfile     = app.m_outfile;
- 
+
     // Copy protected members
     m_covariance      = app.m_covariance;
     m_energies        = app.m_energies;
@@ -583,16 +457,13 @@ void ctbutterfly::get_parameters(void)
 {
     // Setup observations from "inobs" parameter
     setup_observations(m_obs);
-    
+
     // Setup models from "inmodel" parameter
     setup_models(m_obs, (*this)["srcname"].string());
     // Get name of test source
     m_srcname = (*this)["srcname"].string();
-    // Check model name and type
-    check_model();
-
     // Create energy boundaries from user parameters
     m_ebounds = create_ebounds();
@@ -600,18 +471,22 @@ void ctbutterfly::get_parameters(void)
     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.";
+                "Use filename = \"NONE\" to induce a recomputation of "
+                "the matrix internally.";
         throw GException::feature_not_implemented(G_GET_PARAMETERS, msg);
         // m_covariance.load(matrixfilename);
     }
     // Get remaining parameters
     m_apply_edisp = (*this)["edisp"].boolean();
+    m_method      = (*this)["method"].string();
     m_confidence  = (*this)["confidence"].real();
     m_fit         = (*this)["fit"].boolean();
     m_chatter     = static_cast<GChatter>((*this)["chatter"].integer());
+    // Check model name and type
+    check_model();
+
     // Optionally read ahead parameters so that they get correctly
     // dumped into the log file
     if (read_ahead()) {
@@ -627,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
@@ -645,17 +787,22 @@ void ctbutterfly::check_model(void)
         GModelSky* model  = dynamic_cast<GModelSky*>(models[m_srcname]);
         if (model == NULL) {
             std::string msg = "Source \""+m_srcname+"\" is not a sky model. "
-                              "Please specify the name of a sky model for "
-                              "butterfly computation.";
+                    "Please specify the name of a sky model for "
+                    "butterfly computation.";
             throw GException::invalid_value(G_CHECK_MODEL, msg);
         }
-        // Check that the spectral model is a power law
-        if (model->spectral()->type() != "PowerLaw") {
-            std::string msg = "\""+model->spectral()->type()+"\" cannot be "
-                              "used as spectral model for an butterfly "
-                              "computation. Please specify a power law model.";
-            throw GException::invalid_value(G_CHECK_MODEL, msg);
+        // 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.";
+                throw GException::invalid_value(G_CHECK_MODEL, msg);
+            }
         }
     } // endif: source model existed
@@ -663,8 +810,8 @@ void ctbutterfly::check_model(void)
     // Otherwise throw an exception
     else {
         std::string msg = "Source \""+m_srcname+"\" not found in model "
-                          "container. Please add a source with that name "
-                          "or check for possible typos.";
+                "container. Please add a source with that name "
+                "or check for possible typos.";
         throw GException::invalid_value(G_CHECK_MODEL, msg);
     }
@@ -689,13 +836,13 @@ void ctbutterfly::check_model(void)
  * http://www.math.harvard.edu/archive/21b_fall_04/exhibits/2dmatrices/index.html
  ***************************************************************************/
 void ctbutterfly::eigenvectors(const double& a,
-                               const double& b,
-                               const double& c,
-                               const double& d,
-                               double*       lambda1,
-                               double*       lambda2,
-                               GVector*      vector1,
-                               GVector*      vector2)
+        const double& b,
+        const double& c,
+        const double& d,
+        double*       lambda1,
+        double*       lambda2,
+        GVector*      vector1,
+        GVector*      vector2)
 {
     // Compute trace and determinant
     double trace = a + d;
@@ -730,7 +877,7 @@ void ctbutterfly::eigenvectors(const double&amp; a,
     // Normalize eigenvectors
     *vector1 /= norm(*vector1);
     *vector2 /= norm(*vector2);
-    
+
     // Return
     return;
 }
 /***************************************************************************
  *                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,6 +92,7 @@ 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?
@@ -27,6 +27,7 @@ caldb,      s, a, &quot;prod2&quot;,,, &quot;Calibration database&quot;
 irf,        s, a, "South_0.5h",,, "Instrument response function"
 edisp,      b, h, no,,, "Apply energy dispersion?"
 outfile,    f, a, "butterfly.txt",,, "Output ASCII file"
+method,     s, h, GAUSSIAN, GAUSSIAN|ENVELOPE,,"Butterfly computation method"
 #
 # Butterfly parameters
@@ -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"
 /***********************************************************************//**