Implement a new method to compute the butterfly diagram (#1498).

According to a study by Stephan Fegan done in the context of Fermi/LAT analysis, the best way to determine a butterfly diagram is to draw the envelope of the error ellipse of the model parameters. I implemented this method now in ctbutterfly. I also changed the output file, having now two columns giving as function of energy th minimum and maximum intensities.

Implement a new method to compute the butterfly diagram (#1498).
According to a study by Stephan Fegan done in the context of Fermi/LAT analysis, the best way to determine a butterfly diagram is to draw the envelope of the error ellipse of the model parameters. I implemented this method now in ctbutterfly. I also changed the output file, having now two columns giving as function of energy th minimum and maximum intensities.
Jürgen Knödlseder
1 parent c5b2ebf7
Showing 8 changed files with 396 additions and 187 deletions
doc/source/reference_manual/ctbutterfly.rst
doc/source/reference_manual/cterror.rst
doc/source/reference_manual/ctlike.rst
examples/show_butterfly.py
src/ctbutterfly/ctbutterfly.cpp
src/ctbutterfly/ctbutterfly.hpp
src/ctbutterfly/ctbutterfly.par
src/cterror/cterror.par
@@ -9,8 +9,15 @@ Computes butterfly.
 Synopsis
 --------
-Calculates the confidence band of a spectral model, taking into account the
-covariance matrix of a likelihood fit.
+Calculates the confidence band of a power law model, taking into account the
+covariance matrix of a likelihood fit. The tool derives the envelope of 
+power law models with parameters that lie within the error ellipse. It 
+outputs an ASCII file in column separated value (CSV) format containing 
+the results. Each row in the result file corresponds to an energy. The 
+columns give the energy in MeV, the fitted intensity in ph/cm2/s/MeV,
+and the minimum and maximum intensity that envelope the power law model 
+with a given confidence (by default a confidence level of 68% is used, but 
+the level can be adjusted using the ``confidence`` parameter).
 General parameters
@@ -45,7 +52,13 @@ General parameters
 ``outfile [file]``
     Output butterfly ASCII file name.
+
+``(refit = no) [boolean]``
+    Performs refitting of solution ignoring any provided covariance matrix.
+``(confidence = 0.68) [real]``
+    Confidence level for error computation.
+    
 ``(matrix = "NONE") [file]``
     Input covariance matrix file (not used)
@@ -89,11 +102,12 @@ Standard parameters
 ``(mode = ql) [string]``
     Mode of automatic parameters (default is "ql", i.e. "query and learn").
-``(logfile = ctskymap.log) [string]``
+``(logfile = ctbutterfly.log) [string]``
     Name of log file.
 Related tools
 -------------
-None
+:ref:`ctulimit`
+:ref:`cterror`
@@ -89,4 +89,5 @@ Standard parameters
 Related tools
 -------------
-:ref:`ctulimit`
 \ No newline at end of file
+:ref:`ctulimit`
+:ref:`ctbutterfly`
@@ -12,11 +12,10 @@ Synopsis
 Determines source model parameters, such as flux, spectral index, source 
 position, and source extent from a maximum likelihood analysis of CTA data.
 The analysis can be done in a binned or an unbinned formulation of the
-log-likelihood function.
-For binned analysis, a counts cube produced by :doc:`ctbin` is required.
-For unbinned analysis, an event list processed by :doc:`ctselect` is required.
-Based on the input file format, ctlike automatically selects between binned
-and unbinned maximum likelihood analysis.
+log-likelihood function. For binned analysis, a counts cube produced by
+:doc:`ctbin` is required. For unbinned analysis, an event list processed by
+:doc:`ctselect` is required. Based on the input file format, ctlike
+automatically selects between binned and unbinned maximum likelihood analysis.
 General parameters
@@ -18,7 +18,6 @@
 # along with this program.  If not, see <http://www.gnu.org/licenses/>.
 #
 # ==========================================================================
-
 import gammalib
 import sys
 try:
@@ -26,32 +25,34 @@ try:
 except:
     sys.exit("This script needs matplotlib")
-# =============== #
-# Script entry    #
-# =============== #    
-if __name__ == "__main__":
-    # check for given butterfly file
+# ============ #
+# Script entry #
+# ============ #    
+if __name__ == "__main__":
+    """
+    """
+    # Check for given butterfly file
     if not len(sys.argv) == 2:
         sys.exit("Usage: show_butterfly.py butterfly.txt")
     # Read given butterfly file    
     filename = sys.argv[1]
-    csv = gammalib.GCsv(filename)
+    csv      = gammalib.GCsv(filename)
-    # initialise arrays to be filled
+    # Initialise arrays to be filled
     butterfly_x = []
     butterfly_y = []
+    line_x      = []
+    line_y      = []
-    line_x = []
-    line_y = []
-
-    # loop over rows of the file
+    # Loop over rows of the file
     nrows = csv.nrows()
     for row in range(nrows):
+
         # Compute upper edge of confidence band
         butterfly_x.append(csv.real(row,0))
-        butterfly_y.append(csv.real(row,1)+csv.real(row,2))
+        butterfly_y.append(csv.real(row,2))
         # Set line values
         line_x.append(csv.real(row,0))
@@ -61,21 +62,19 @@ if __name__ == &quot;__main__&quot;:
     # of the confidence band    
     for row in range(nrows):
         index = nrows - 1 - row
-
         butterfly_x.append(csv.real(index,0))
-        low_error = csv.real(index,1)-csv.real(index,2)
+        low_error = csv.real(index,3)
         if low_error < 1e-26:
             low_error = 1e-26
         butterfly_y.append(low_error)   
-    # plot the butterfly and spectral line       
+    # Plot the butterfly and spectral line       
     plt.figure()
     plt.ylim([1e-26,1e-14])
     plt.loglog()
     plt.grid()
     plt.plot(line_x,line_y,color='black',ls='-')
     plt.fill(butterfly_x,butterfly_y,color='green',alpha=0.5)
-
-    plt.xlabel("Energy [MeV]")
-    plt.ylabel(r"dN/dE [MeV$^{-1}$ s$^{-1}$ cm$^{-2}$]")    
+    plt.xlabel("Energy (MeV)")
+    plt.ylabel(r"dN/dE (MeV$^{-1}$ s$^{-1}$ cm$^{-2}$)")    
     plt.show()
@@ -36,7 +36,7 @@
 /* __ Method name definitions ____________________________________________ */
 #define G_GET_PARAMETERS                      "ctbutterfly::get_parameters()"
 #define G_SAVE                                          "ctbutterfly::save()"
-#define G_GET_MODEL_PARAMETER            "ctbutterfly::get_model_parameter()"
+#define G_CHECK_MODEL                            "ctbutterfly::check_model()"
 /* __ Debug definitions __________________________________________________ */
@@ -241,11 +241,40 @@ void ctbutterfly::run(void)
         log << m_obs << std::endl;
     }
+    // If refit is selected then do maximum likelihood fit
+    if (m_refit) {
+
+        // Write header
+        if (logTerse()) {
+            log << std::endl;
+            log.header1("Compute best-fit likelihood");
+        }
+
+        // Optimize and compute errors
+        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.
+        GObservations::likelihood likelihood = m_obs.function();
+        m_covariance = likelihood.curvature()->invert();
+
+        // Write optimised model and covariance matrix into logger
+        if (logTerse()) {
+            log << m_opt << std::endl;
+            log << m_obs.models() << std::endl;
+            log << m_covariance << std::endl;
+        }
+
+    } // endif: refit selected
+
     // Get model instance for further computations
     GModels models = m_obs.models();
-    // Compute covariance matrix if necessary
-    if (m_covariance.size() == 0) {
+    // If no refit is selected and no covariance matrix has been specified
+    // then compute now the covariance matrix
+    if (!m_refit && m_covariance.size() == 0) {
         // Write header
         if (logTerse()) {
@@ -258,101 +287,195 @@ void ctbutterfly::run(void)
         GObservations::likelihood likelihood = m_obs.function();
         likelihood.eval(pars);
-        // Get covariance matrix by inverting the curvature matrix
+        // 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();
-    }
+        // Log covariance matrix
+        if (logTerse()) {
+            log << m_covariance << std::endl;
+        }
+
+    } // endif: covariance computation needed
     // Write header
     if (logTerse()) {
         log << std::endl;
-        log.header1("Generate butterfly");
+        log.header1("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;
+        }
     }
-    // Initialise dummy time to evaluate spectral model
-    GTime time = GTime();
+    // Check indices
-    // Initialise vector of gradients
-    GVector grad = GVector(m_obs.models().npars());
+    // 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");
-    // Loop over energy bins positions
-    for (int i = 0; i < m_ebounds.size(); ++i) {
+    // Show parameter indices
+    if (logTerse()) {
+        log << gammalib::parformat("Prefactor");
+        log << prefactor_mean << " ph/cm2/s/MeV (";
+        log << inx_prefactor << ")" << std::endl;
+        log << gammalib::parformat("Index");
+        log << index_mean << " (" << inx_index << ")" << std::endl;
+        log << gammalib::parformat("Pivot energy");
+        log << pivot << std::endl;
+    }
-        // Initialise number of gradients
-        int num_gradient = 0;
+    // 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);
-        // Get the energy of current bin
-        GEnergy energy = m_ebounds.elogmean(i);
+    // Initialise min and max intensities
+    m_min_intensities.assign(m_ebounds.size(), 1e30);
+    m_max_intensities.assign(m_ebounds.size(), 0.0);
-        // Header for verbose logging
-        if (logExplicit()) {
-            std::string msg = "Computing butterfly for bin number "+
-                              gammalib::str(i)+" at "+energy.print();
-            log << std::endl;
-            log.header2(msg);
-        }
+    // 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);
-        // Initialise flux value
-        double model_flux = 0.0;
+    // Show eigenvectors and eigenvalues
+    if (logTerse()) {
+        log << gammalib::parformat("Eigenvalue 1")  << lambda1 << std::endl;
+        log << gammalib::parformat("Eigenvalue 2")  << lambda2 << std::endl;
+        log << gammalib::parformat("Eigenvector 1") << vector1 << std::endl;
+        log << gammalib::parformat("Eigenvector 2") << vector2 << std::endl;
+    }
-        // Loop over models
-        for (int j = 0; j < models.size(); ++j) {
+    // Confidence scaling
+    double sigma = gammalib::erfinv(m_confidence) * gammalib::sqrt_two;
+    double scale = (sigma*sigma);
-            // If source name is found, retrieve the nominal flux
-            if (models[j]->name() == m_srcname) {
+    // Show confidence stuff
+    if (logTerse()) {
+        log << gammalib::parformat("Confidence level");
+        log << m_confidence << std::endl;
+        log << gammalib::parformat("Corresponding scaling");
+        log << scale << std::endl;
+    }
-                // Note that we have to compute the gradients here
-                // repetively to update the gradients in the vector
-                // No check for valid pointers necessary here as this 
-                // was covered by get_model_parameter()
-                GModelSky*      skymodel = dynamic_cast<GModelSky*>(models[j]);
-                GModelSpectral* spectral = skymodel->spectral();
-                model_flux = spectral->eval_gradients(energy, time);
+    // Compute minor and major axes of error ellipse
+    double major = scale*std::sqrt(lambda1);
+    double minor = scale*std::sqrt(lambda2);
-            }
+    // Write header
+    if (logTerse()) {
+        log << std::endl;
+        log.header1("Generate butterfly");
+    }
-            // Loop over model parameters
-            for (int i = 0; i < models[j]->size(); ++i) {
+    // 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;
+
+        // Show prefactor and index
+        if (logTerse()) {
+            log << gammalib::parformat("Angle "+gammalib::str(t*gammalib::rad2deg));
+            log << prefactor << " ph/cm2/s/MeV; ";
+            log << index << std::endl;
+        }
-                // Get model parameter
-                GModelPar parameter = models[j]->at(i);
+        // Setup power law model
+        GModelSpectralPlaw plaw(prefactor, index, pivot);
-                // Set gradient value to the vector
-                grad[num_gradient] = parameter.gradient();
+        // Evaluate power law at each energy and extract minimum and
+        // maximum intensity
+        for (int k = 0; k < m_ebounds.size(); ++k) {
-                // increment vector index
-                num_gradient++;
+            // Get the energy of current bin
+            GTime   time;
+            GEnergy energy = m_ebounds.elogmean(k);
-            } // endfor: Looped over model parameters
+            // 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
+
+    // Write header
+    if (logTerse()) {
+        log << std::endl;
+        log.header1("Results");
+    }
-        } // endfor: Looped over models
+    // Setup power law model
+    GModelSpectralPlaw plaw(prefactor_mean, index_mean, pivot);
-        // Multiply covariance to the gradient vector
-        GVector vector = m_covariance * grad;
+    // Show results
+    for (int k = 0; k < m_ebounds.size(); ++k) {
-        // Get the error from the scalar product
-        double error = std::sqrt(grad * vector);
+        // Get the energy of current bin
+        GTime   time;
+        GEnergy energy = m_ebounds.elogmean(k);
-        // Multiply in scale of spectral parameter
-        error *= m_par_scale;
+        // Evaluate power law
+        double intensity = plaw.eval(energy, time);
-        // Store flux, value and energy for saving
-        m_fluxes.push_back(model_flux);
+        // Store results
         m_energies.push_back(energy.MeV());
-        m_errors.push_back(error);
+        m_intensities.push_back(intensity);
         // Log information
-        if (logExplicit()) {
-            log << " Flux at ";
-            log << energy.print();
-            log << ": ";
-            log << model_flux;
-            log << " +- ";
-            log << error;
-            log << " ph/cm2/s/MeV";
+        if (logTerse()) {
+            log << gammalib::parformat("Intensity at "+energy.print());
+            log << intensity;
+            log << " (";
+            log << m_min_intensities[k];
+            log << ", ";
+            log << m_max_intensities[k];
+            log << ") ph/cm2/s/MeV";
+            log << std::endl;
         }
-    } //endfor: loop over energy bin
+    } // endfor
     // Restore energy dispersion flag for all CTA observations
     for (int i = 0; i < m_obs.size(); ++i) {
@@ -384,14 +507,15 @@ void ctbutterfly::save(void)
     // Get output filename
     m_outfile = (*this)["outfile"].filename();
-    // Create CSV table with 3 columns
-    GCsv table(m_fluxes.size(), 3);
+    // Create CSV table with 4 columns
+    GCsv table(m_energies.size(), 4);
     // Fill CSV table
-    for (int i = 0; i < m_fluxes.size(); ++i) {
+    for (int i = 0; i < m_energies.size(); ++i) {
         table.real(i, 0, m_energies[i]);
-        table.real(i, 1, m_fluxes[i]);
-        table.real(i, 2, m_errors[i]);
+        table.real(i, 1, m_intensities[i]);
+        table.real(i, 2, m_min_intensities[i]);
+        table.real(i, 3, m_max_intensities[i]);
     }
     // Save CSV table
@@ -415,18 +539,25 @@ void ctbutterfly::init_members(void)
 {
     // Initialise members
     m_srcname.clear();
-    m_outfile.clear();
-    m_ebounds.clear();
+    m_confidence  = 0.68;
+    m_max_iter    = 50;
     m_apply_edisp = false;
+    m_refit       = false;
+    m_ebounds.clear();
+    m_outfile.clear();
     // Initialise protected members
-    m_skymodel  = NULL;
-    m_par_scale = 0.0;
+    m_opt.clear();
     m_obs.clear();
     m_covariance.clear();
     m_energies.clear();
-    m_fluxes.clear();
-    m_errors.clear();
+    m_intensities.clear();
+    m_min_intensities.clear();
+    m_max_intensities.clear();
+
+    // Set optimizer parameters
+    m_opt.max_iter(m_max_iter);
+    m_opt.max_stalls(10);
     // Return
     return;
@@ -441,20 +572,22 @@ void ctbutterfly::init_members(void)
 void ctbutterfly::copy_members(const ctbutterfly& app)
 {
     // Copy attributes
-    m_srcname = app.m_srcname;
-    m_outfile = app.m_outfile;
-    m_ebounds = app.m_ebounds;
+    m_srcname     = app.m_srcname;
+    m_confidence  = app.m_confidence;
+    m_max_iter    = app.m_max_iter;
     m_apply_edisp = app.m_apply_edisp;
-
+    m_refit       = app.m_refit;
+    m_ebounds     = app.m_ebounds;
+    m_outfile     = app.m_outfile;
+ 
     // Copy protected members
-    m_obs        = app.m_obs;
-    m_covariance = app.m_covariance;
-    m_energies   = app.m_energies;
-    m_fluxes     = app.m_fluxes;
-    m_errors     = app.m_errors;
-
-    // Extract model parameter
-    get_model_parameter();
+    m_obs             = app.m_obs;
+    m_opt             = app.m_opt;
+    m_covariance      = app.m_covariance;
+    m_energies        = app.m_energies;
+    m_intensities     = app.m_intensities;
+    m_min_intensities = app.m_min_intensities;
+    m_max_intensities = app.m_max_intensities;
     // Return
     return;
@@ -513,15 +646,9 @@ void ctbutterfly::get_parameters(void)
     // Get name of test source and check container for this name
     m_srcname = (*this)["srcname"].string();
-    if (!m_obs.models().contains(m_srcname)) {
-        std::string msg = "Source \""+m_srcname+"\" not found in model "
-                          "container. Please add a source with that name "
-                          "or check for possible typos.";
-        throw GException::invalid_value(G_GET_PARAMETERS, msg);
-    }
-    // Get relevant model and parameter for model computation
-    get_model_parameter();
+    // Check model name and type
+    check_model();
     // Create energy boundaries from user parameters
     m_ebounds = create_ebounds();
@@ -539,6 +666,10 @@ void ctbutterfly::get_parameters(void)
     // Read energy dispersion flag
     m_apply_edisp = (*this)["edisp"].boolean();
+    // Read other parameters
+    m_confidence = (*this)["confidence"].real();
+    m_refit      = (*this)["refit"].boolean();
+
     // Optionally read ahead parameters so that they get correctly
     // dumped into the log file
     if (read_ahead()) {
@@ -551,62 +682,114 @@ void ctbutterfly::get_parameters(void)
 /***********************************************************************//**
- * @brief Get application parameters
+ * @brief Check if sky model is valid
  *
  * @exception GException::invalid_value
- *            Did not find a valid model parameter
+ *            Did not find a valid model
  *
- * Extracts a pointer to the sky model (m_skymodel) and a pointer to the
- * relevant model parameter (m_model_par) from the model container.
+ * Checks whether the specified sky model exists and whether it is a power
+ * law.
  ***************************************************************************/
-void ctbutterfly::get_model_parameter(void)
+void ctbutterfly::check_model(void)
 {
     // Continue only if source model exists
     if (m_obs.models().contains(m_srcname)) {
         // Get relevant model and parameter for upper limit computation.
-        GModels& models = const_cast<GModels&>(m_obs.models());
-        m_skymodel      = dynamic_cast<GModelSky*>(models[m_srcname]);
-        if (m_skymodel == NULL) {
+        GModels&   models = const_cast<GModels&>(m_obs.models());
+        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.";
-            throw GException::invalid_value(G_GET_MODEL_PARAMETER, msg);
+            throw GException::invalid_value(G_CHECK_MODEL, msg);
         }
-        GModelPar* modelpar;
-
-        if (m_skymodel->spectral()->type() == "NodeFunction") {
-            std::string msg = "\"NodeFunction\" cannot be used as spectral "
-                              "model for an butterfly computation. "
-                              "Please specify another spectral model.";
-            throw GException::invalid_value(G_GET_MODEL_PARAMETER, msg);
-        }
-        if (m_skymodel->spectral()->has_par("Normalization")) {
-            modelpar = &(m_skymodel->spectral()->operator[]("Normalization"));
-        }
-        else if (m_skymodel->spectral()->has_par("Prefactor")) {
-            modelpar = &(m_skymodel->spectral()->operator[]("Prefactor"));
-        }
-        else if (m_skymodel->spectral()->has_par("Integral")) {
-            modelpar = &(m_skymodel->spectral()->operator[]("Integral"));
-        }
-        else if (m_skymodel->spectral()->has_par("Value")) {
-            modelpar = &(m_skymodel->spectral()->operator[]("Value"));
-        }
-        else {
-            std::string msg = "Require spectral parameter \"Normalization\", "
-                              "\"Prefactor\", \"Integral\" or \"Value\" for "
-                              "butterfly computation. The specified source "
-                              "\""+m_srcname+"\" does not have such a parameter.";
-            throw GException::invalid_value(G_GET_MODEL_PARAMETER, 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);
         }
-        // Set scale value from parameter
-        m_par_scale = modelpar->scale();
-
     } // endif: source model existed
+    // 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.";
+        throw GException::invalid_value(G_CHECK_MODEL, msg);
+    }
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Compute normalized eigenvectors and eigenvalues
+ *
+ * @param[in] a Covariance matrix element (0,0)
+ * @param[in] b Covariance matrix element (0,1)
+ * @param[in] c Covariance matrix element (1,0)
+ * @param[in] d Covariance matrix element (1,1)
+ * @param[out] lambda1 First eigenavlue
+ * @param[out] lambda2 Second eigenavlue
+ * @param[out] vector1 First eigenvector (normalized)
+ * @param[out] vector2 Second eigenvector (normalized)
+ *
+ * Computes the eigenvalues and eigenvectors for a 2x2 matrix. See
+ * 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)
+{
+    // Compute trace and determinant
+    double trace = a + d;
+    double det   = a*d - b*c;
+
+    // Compute eigenvalues
+    double arg1 = 0.5*trace;
+    double arg2 = std::sqrt(0.25*trace*trace - det);
+    *lambda1    = arg1 + arg2;
+    *lambda2    = arg1 - arg2;
+
+    // Compute eigenvectors
+    double vector1_1;
+    double vector1_2;
+    double vector2_1;
+    double vector2_2;
+    if (c != 0.0) {
+        (*vector1)[0] = *lambda1 - det;
+        (*vector1)[1] = c;
+        (*vector2)[0] = *lambda2 - det;
+        (*vector2)[1] = c;
+    }
+    else if (b != 0.0) {
+        (*vector1)[0] = b;
+        (*vector1)[1] = *lambda1 - a;
+        (*vector2)[0] = b;
+        (*vector2)[1] = *lambda2 - a;
+    }
+    else {
+        (*vector1)[0] = 1.0;
+        (*vector1)[1] = 0.0;
+        (*vector2)[0] = 0.0;
+        (*vector2)[1] = 1.0;
+    }
+
+    // Normalize eigenvectors
+    *vector1 /= norm(*vector1);
+    *vector2 /= norm(*vector2);
+    
     // Return
     return;
 }
@@ -78,22 +78,33 @@ protected:
     void copy_members(const ctbutterfly& app);
     void free_members(void);
     void get_parameters(void);
-    void get_model_parameter(void);
+    void check_model(void);
+    void eigenvectors(const double& a,
+                      const double& b,
+                      const double& c,
+                      const double& d,
+                      double*       lambda1,
+                      double*       lambda2,
+                      GVector*      vector1,
+                      GVector*      vector2);
     // User parameters
-    std::string m_srcname;     //!< Name of source to compute butterfly
-    std::string m_outfile;     //!< Output ascii file
-    GEbounds    m_ebounds;     //!< Energy binning definition
-    bool        m_apply_edisp; //!< Apply energy dispersion?
+    std::string m_srcname;      //!< Name of source to compute butterfly
+    double      m_confidence;   //!< Confidence level
+    int         m_max_iter;     //!< Maximum number of iterations
+    bool        m_apply_edisp;  //!< Apply energy dispersion?
+    bool        m_refit;        //!< Do refit?
+    GEbounds    m_ebounds;      //!< Energy binning definition
+    std::string m_outfile;      //!< Output ascii file
     // Protected members
-    GObservations       m_obs;        //!< Observation container
-    GMatrixSparse       m_covariance; //!< Covariance matrix
-    GModelSky*          m_skymodel;   //!< Pointer to sky model
-    double              m_par_scale;  //!< Scale of normalising parameter
-    std::vector<double> m_energies;   //!< Energy values for storage
-    std::vector<double> m_fluxes;     //!< Flux values per energy bin
-    std::vector<double> m_errors;     //!< Flux errors per energy bin
+    GOptimizerLM        m_opt;             //!< Optimizer
+    GObservations       m_obs;             //!< Observation container
+    GMatrixSparse       m_covariance;      //!< Covariance matrix
+    std::vector<double> m_energies;        //!< Energy values
+    std::vector<double> m_intensities;     //!< Power law intensity
+    std::vector<double> m_min_intensities; //!< Minimum intensities
+    std::vector<double> m_max_intensities; //!< Maximum intensities
 };
 /***********************************************************************//**
@@ -17,17 +17,19 @@
 #
 # General parameters
 #===================
-inobs,    f, a, "events.fits",,,"Input event list, counts cube or observation definition file"
-inmodel,  f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
-srcname,  s, a, "Crab",,,"Source of interest"
-expcube,  f, a, "NONE",,, "Exposure cube file (only needed for stacked analysis)"
-psfcube,  f, a, "NONE",,, "PSF cube file (only needed for stacked analysis)"
-bkgcube,  f, a, "NONE",,, "Background cube file (only needed for stacked analysis)"
-caldb,    s, a, "prod2",,, "Calibration database"
-irf,      s, a, "South_50h",,, "Instrument response function"
-edisp,    b, h, no,,, "Apply energy dispersion?"
-outfile,  f, a, "butterfly.txt",,, "Output ascii file"
-matrix,   f, h, "NONE",,, "Input covariance Matrix FITS file"
+inobs,      f, a, "events.fits",,,"Input event list, counts cube or observation definition file"
+inmodel,    f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
+srcname,    s, a, "Crab",,,"Source of interest"
+expcube,    f, a, "NONE",,, "Exposure cube file (only needed for stacked analysis)"
+psfcube,    f, a, "NONE",,, "PSF cube file (only needed for stacked analysis)"
+bkgcube,    f, a, "NONE",,, "Background cube file (only needed for stacked analysis)"
+caldb,      s, a, "prod2",,, "Calibration database"
+irf,        s, a, "South_50h",,, "Instrument response function"
+edisp,      b, h, no,,, "Apply energy dispersion?"
+outfile,    f, a, "butterfly.txt",,, "Output ascii file"
+refit,      b, h, no,,, "Do refitting?"
+confidence, r, h, 0.68,0.0,1.0, "Confidence level"
+matrix,     f, h, "NONE",,, "Input covariance Matrix FITS file"
 #
 # Energy binning parameters
@@ -6,7 +6,7 @@
 # ------------------------------------------------------------------------- #
 #                                                                           #
 #  File:      cterror.par                                                   #
-#  Version:   01-00-00                                                      #
+#  Version:   1.0.0                                                         #
 #                                                                           #
 #  Author:    Florent Forest                                                #
 #             florent.forest@isae.fr                                        #