Merge branch '1051-create-ctbutterfly' of https://cta-git.irap.omp.eu/ctools int…

…o 1051-create-ctbutterfly

Merge branch '1051-create-ctbutterfly' of https://cta-git.irap.omp.eu/ctools int…
…o 1051-create-ctbutterfly
Michael Mayer
2 parents df8dd6d9 3d171380
Showing 12 changed files with 113 additions and 95 deletions
.gitignore
AUTHORS
ChangeLog
NEWS
pyext/ctbutterfly.i
src/ctbutterfly/Makefile.am
src/ctbutterfly/ctbutterfly.cpp
src/ctbutterfly/ctbutterfly.hpp
src/ctbutterfly/ctbutterfly.par
src/ctbutterfly/main.cpp
test/test_ctbutterfly.py
test/test_ctools.sh
@@ -57,6 +57,7 @@ src/ctexpcube/ctexpcube
 src/ctpsfcube/ctpsfcube
 src/ctbkgcube/ctbkgcube
 src/ctcubemask/ctcubemask
+src/ctbutterfly/ctbutterfly
 test/pfiles/
 test/*.fits
 test/*.dat
@@ -64,6 +65,7 @@ test/*.xml
 test/*.trs
 test/test_cscripts.sh
 test/reports/*.xml
+test/butterfly.txt
 examples/*.fits
 examples/*.xml
 examples/*.dat
 Juergen Knoedlseder (jurgen.knodlseder@irap.omp.eu)
 Jean-Baptiste Cayrou (cayrou@etud.insa-toulouse.fr)
 Christoph Deil (Christoph.Deil@mpi-hd.mpg.de)
-Michael Mayer (michael.mayer@desy.de)
+Michael Mayer (michael.mayer@physik.hu-berlin.de)
 Chia-Chun Lu (Chia-Chun.Lu@mpi-hd.mpg.de)
-2014-10-15  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
+2014-10-26  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
         * Version 0.8.0 released
         ========================
@@ -10,6 +10,7 @@
         Add ctbkgcube tool (generates background cube)
         Add ctcubemask tool (masks event bins in event cube)
         Add cttsmap tool (generates TS map)
+        Add ctbutterfly tool (generates butterfly diagram)
         Add csresmap script (generates residual map)
         Introduce ctool base class
         ctbin now always merges event lists into a single counts cube
@@ -2,7 +2,7 @@ New Features and Important Changes in ctools 00-08-00
 Author(s) : Juergen Knoedlseder
-15 October 2014
+26 October 2014
 1. Introduction
@@ -168,12 +168,20 @@ bin content to -1. Filtering can be done by energy, for a circular
 Region of Interest (ROI) or using ds9 exclusion regions.
-16. Examples
+16. ctbutterfly - Generate butterfly diagrams
+=============================================
+The first version of the ctbutterfly tool is included in the release.
+This tools derives a butterfly diagram from the fitted model parameters
+for one specified source. The result will be written into an ASCII
+file.
+
+
+17. Examples
 ------------
 None
-17. Scripts
+18. Scripts
 -----------
 Add logL and event statistics to cspull output file.
@@ -196,11 +204,11 @@ Add the csresmap script that uses ctbin and ctmodel to generate a
 residual counts cube for a given model.
-18. caldb
+19. caldb
 ---------
 None
-19. models
+20. models
 ----------
 None
 /***************************************************************************
- *                    ctbutterfly - butterfly calculation tool             *
+ *                  ctbutterfly - butterfly calculation tool               *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2014 by Michael Mayer                                    *
  * ----------------------------------------------------------------------- *
@@ -20,7 +20,7 @@
  ***************************************************************************/
 /**
  * @file ctbutterfly.i
- * @brief tool for butterfly calculations
+ * @brief Butterfly calculation tool
  * @author Michael Mayer
  */
 %{
@@ -32,7 +32,7 @@
 /***********************************************************************//**
  * @class ctbutterfly
  *
- * @brief butterfly calculation tool
+ * @brief Butterfly calculation tool
  ***************************************************************************/
 class ctbutterfly : public ctool {
 public:
@@ -48,7 +48,6 @@ public:
     void                 run(void);
     void                 save(void);
     const GObservations& obs(void) const;
-    
 };
+## Process this file with automake to produce Makefile.in
+
+# Set include directories
+INCLUDES = -I$(top_srcdir)/src/support
+
+# Program name
+bin_PROGRAMS = ctbutterfly
+
+# Sources
+sources = ctbutterfly.cpp
+
+# Libtool library
+noinst_LTLIBRARIES = libctbutterfly.la
+libctbutterfly_la_SOURCES = $(sources)
+
+# Program definition
+ctbutterfly_SOURCES = main.cpp
+ctbutterfly_LDFLAGS = $(LIB_LINK_OPTION) @LDFLAGS@
+ctbutterfly_CPPFLAGS = @CPPFLAGS@
+ctbutterfly_LDADD = libctbutterfly.la $(top_srcdir)/src/support/libsupport.la
+
+# Parameter files
+pfilesdir = $(prefix)/syspfiles
+pfiles_DATA = ctbutterfly.par
 /***************************************************************************
- *                    ctbutterfly - butterfly calculation tool                    *
+ *                 ctbutterfly - butterfly calculation tool                *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2014 by Michael Mayer                                    *
  * ----------------------------------------------------------------------- *
@@ -34,8 +34,8 @@
 #include "GTools.hpp"
 /* __ Method name definitions ____________________________________________ */
-#define G_GET_PARAMETERS                          "ctbutterfly::get_parameters()"
-#define G_SAVE                                                  "ctbutterfly::save()"
+#define G_GET_PARAMETERS                      "ctbutterfly::get_parameters()"
+#define G_SAVE                                          "ctbutterfly::save()"
 /* __ Debug definitions __________________________________________________ */
@@ -70,7 +70,7 @@ ctbutterfly::ctbutterfly(void) : ctool(CTBUTTERFLY_NAME, CTBUTTERFLY_VERSION)
  * observations container.
  ***************************************************************************/
 ctbutterfly::ctbutterfly(const GObservations& obs) :
-         ctool(CTBUTTERFLY_NAME, CTBUTTERFLY_VERSION)
+             ctool(CTBUTTERFLY_NAME, CTBUTTERFLY_VERSION)
 {
     // Initialise members
     init_members();
@@ -91,7 +91,7 @@ ctbutterfly::ctbutterfly(const GObservations&amp; obs) :
  * @param[in] argv Array of command line arguments.
  ***************************************************************************/
 ctbutterfly::ctbutterfly(int argc, char *argv[]) :
-         ctool(CTBUTTERFLY_NAME, CTBUTTERFLY_VERSION, argc, argv)
+             ctool(CTBUTTERFLY_NAME, CTBUTTERFLY_VERSION, argc, argv)
 {
     // Initialise members
     init_members();
@@ -240,11 +240,13 @@ void ctbutterfly::run(void)
             log.header1("Compute covariance matrix");
         }
-		GOptimizerPars pars = models.pars();
+        // Evaluate curvature matrix at the actual parameters
+		GOptimizerPars            pars       = models.pars();
 		GObservations::likelihood likelihood = m_obs.function();
 		likelihood.eval(pars);
-		GMatrixSparse* curvature = dynamic_cast<GMatrixSparse*>(likelihood.curvature());
-		m_covariance = curvature->invert();
+
+        // Get covariance matrix by inverting the curvature matrix
+        m_covariance = likelihood.curvature()->invert();
 	}
@@ -263,6 +265,7 @@ void ctbutterfly::run(void)
     // 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
@@ -270,8 +273,8 @@ void ctbutterfly::run(void)
     	// Header for verbose logging
     	if (logExplicit()) {
-            std::string msg = "Computing butterfly for bin number "+gammalib::str(i)+
-                              " at "+energy.print();
+            std::string msg = "Computing butterfly for bin number "+
+                              gammalib::str(i)+" at "+energy.print();
             log << std::endl;
             log.header2(msg);
     	}
@@ -280,7 +283,7 @@ void ctbutterfly::run(void)
     	double model_flux = 0.0;
     	// Loop over models
-    	for (int j = 0; j < models.size(); ++j ) {
+    	for (int j = 0; j < models.size(); ++j) {
     	    // Check wether model is a skymodel
     	    GModelSky* skymodel = dynamic_cast<GModelSky*>(models[j]);
@@ -296,17 +299,15 @@ void ctbutterfly::run(void)
     	        // Set flux value of source of interest
     	        if (skymodel->name() == m_srcname) {
-    	               model_flux = spectral->eval_gradients(energy, time);
+                    model_flux = spectral->eval_gradients(energy, time);
     	        } // endif: model was source of interest
     	        // Loop over model parameters, get gradients
     	        // and assign them to the vector
     	        for (int k = 0; k < spectral->size(); ++k) {
-
-    	            grad[num_gradient] =  (*spectral)[k].gradient();
-    	            num_gradient += 1;
-
-    	        }// endfor: loop over spectral parameters
+    	            grad[num_gradient] = (*spectral)[k].gradient();
+    	            num_gradient++;
+    	        } // endfor: looped over spectral parameters
     	        // Skip temporal models
     	        num_gradient += skymodel->temporal()->size();
@@ -352,9 +353,10 @@ void ctbutterfly::run(void)
 /***********************************************************************//**
- * @brief Save maps
+ * @brief Save butterfly diagram
  *
- * This method saves the butterfly into an ascii file
+ * Saves the butterfly diagram into an ascii file using a column separated
+ * value (CSV) format with blanks as separators.
  ***************************************************************************/
 void ctbutterfly::save(void)
 {
@@ -367,33 +369,18 @@ void ctbutterfly::save(void)
     // Get output filename
     m_outfile = (*this)["outfile"].filename();
-    if (gammalib::file_exists(m_outfile)) {
-
-        if ((*this)["clobber"].boolean()) {
-            std::remove(m_outfile.c_str());
-        }
-        else {
-            std::string msg = "File \""+m_outfile+"\" already exists and clobber is set to \"false\"";
-            throw GException::file_open_error(G_SAVE,msg);
-        }
-    }
-
-    // Create  file instance
-    std::fstream file;
+    // Create CSV table with 3 columns
+    GCsv table(m_fluxes.size(), 3);
-    file.open(m_outfile, std::fstream::out | std::fstream::app);
-
-    // Write short header
-    file << "# Energy [MeV] | flux [ph/cm2/s/MeV] | e_flux [ph/cm2/s/MeV]" << std::endl;
-
-    // Write values to file
+    // Fill CSV table
     for (int i = 0; i < m_fluxes.size(); ++i) {
-
-        file << m_energies[i] << " " << m_fluxes[i] << " " << m_errors[i] << std::endl;
+        table.real(i, 0, m_energies[i]);
+        table.real(i, 1, m_fluxes[i]);
+        table.real(i, 2, m_errors[i]);
     }
-    // close file
-    file.close();
+    // Save CSV table
+    table.save(m_outfile, " ", clobber());
     // Return
     return;
@@ -413,15 +400,15 @@ void ctbutterfly::init_members(void)
 {
     // Initialise members
     m_infile.clear();
-    m_outfile.clear();
     m_srcname.clear();
+    m_outfile.clear();
     m_ebounds.clear();
     // Initialise protected members
     m_obs.clear();
     m_covariance.clear();
-    m_fluxes.clear();
     m_energies.clear();
+    m_fluxes.clear();
     m_errors.clear();
     // Return
@@ -437,17 +424,17 @@ void ctbutterfly::init_members(void)
 void ctbutterfly::copy_members(const ctbutterfly& app)
 {
     // Copy attributes
-    m_infile   = app.m_infile;
-    m_outfile  = app.m_outfile;
-    m_srcname     = app.m_srcname;
-    m_ebounds     = app.m_ebounds;
+    m_infile  = app.m_infile;
+    m_srcname = app.m_srcname;
+    m_outfile = app.m_outfile;
+    m_ebounds = app.m_ebounds;
     // Copy protected members
     m_obs        = app.m_obs;
     m_covariance = app.m_covariance;
-    m_fluxes         = app.m_fluxes;
-    m_energies    = app.m_energies;
-    m_errors         = app.m_errors;
+    m_energies   = app.m_energies;
+    m_fluxes     = app.m_fluxes;
+    m_errors     = app.m_errors;
     // Return
     return;
@@ -535,11 +522,11 @@ void ctbutterfly::get_parameters(void)
     if (!m_obs.models().contains(m_srcname)) {
         std::string msg = "Source \""+m_srcname+"\" not found in model "
                           "container. Please add a source with that name "
-                          "or check for a possible typos.";
+                          "or check for possible typos.";
     	throw GException::invalid_value(G_GET_PARAMETERS, msg);
     }
-    // get binning information
+    // Get energy binning information
     m_ebounds  = get_ebounds();
     // Get matrix file name and load if possible
@@ -561,4 +548,3 @@ void ctbutterfly::get_parameters(void)
     // Return
     return;
 }
-
 /***************************************************************************
- *                    ctbutterfly - butterfly calculation tool                    *
+ *                ctbutterfly - butterfly calculation tool                 *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2014 by Michael Mayer                                    *
  * ----------------------------------------------------------------------- *
@@ -40,17 +40,18 @@
 /***********************************************************************//**
  * @class ctbutterfly
  *
- * @brief butterfly calculation tool
+ * @brief Butterfly calculation tool
  *
  * This class computes the confidence interval of a fitted spectrum
- * (butterfly) for a given set of observations
+ * (butterfly) for a given set of observations.
  *
  * The class operates on predefined observation containers, an individual
  * event list or an observation definition XML file.
  *
- * During the computation the covariance matrix of the fit is used to propagate
- * uncertainties and their correlations through the entire energy range.
- * The output is saved as an ascii files containing the confidence band boundaries
+ * During the computation the covariance matrix of the fit is used to
+ * propagate uncertainties and their correlations through the entire energy
+ * range. The output is saved as an ascii files containing the confidence
+ * band boundaries
  ***************************************************************************/
 class ctbutterfly : public ctool {
@@ -79,19 +80,17 @@ protected:
     void get_parameters(void);
     // User parameters
-    std::string              m_infile;     //!< Input file
-    std::string              m_srcname;    //!< Name of source to compute butterfly
-    std::string              m_outfile;    //!< Output ascii file
-
-    GEbounds            m_ebounds; //!< Energy binning definition
+    std::string m_infile;     //!< Input file
+    std::string m_srcname;    //!< Name of source to compute butterfly
+    std::string m_outfile;    //!< Output ascii file
+    GEbounds    m_ebounds;    //!< Energy binning definition
       // Protected members
-    GObservations            m_obs;        //!< Observation container
-    GMatrixSparse           m_covariance; //!< covariance matrix
-
-    std::vector<double>   m_energies; //!< energy values for storage
-    std::vector<double>   m_fluxes; //!< flux values per energy bin
-    std::vector<double>   m_errors; //!< flux errors per energy bin
+    GObservations       m_obs;        //!< Observation container
+    GMatrixSparse       m_covariance; //!< Covariance matrix
+    std::vector<double> m_energies;   //!< Energy values for storage
+    std::vector<double> m_fluxes;     //!< Flux values per energy bin
+    std::vector<double> m_errors;     //!< Flux errors per energy bin
 };
 /***********************************************************************//**
@@ -19,7 +19,7 @@
 #===================
 infile,   f, a, "events.fits",,,"Input event list, cube or observation definition file"
 srcmdl,   f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
-srcname,  s, a, "CrabNebula",,,"Source of interest"
+srcname,  s, a, "Crab",,,"Source of interest"
 caldb,    s, a, "dummy",,, "Calibration database"
 irf,      s, a, "cta_dummy_irf",,, "Instrument response function"
 outfile,  f, a, "butterfly.txt",,, "Output ascii file"
 /***************************************************************************
- *                      ctbutterfly - TS map calculation tool                  *
+ *                   ctbutterfly - TS map calculation tool                 *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2014 by Michael Mayer                                    *
  * ----------------------------------------------------------------------- *
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This scripts performs unit tests for the ctbutterfly tool.
 #
-# Copyright (C) 2014 Juergen Knoedlseder
+# Copyright (C) 2014 Michal Mayer
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU General Public License as published by
@@ -22,9 +22,9 @@ import gammalib
 import ctools
-# ========================== #
+# =============================== #
 # Test class for ctbutterfly tool #
-# ========================== #
+# =============================== #
 class Test(gammalib.GPythonTestSuite):
     """
     Test class for ctbutterfly tool.
@@ -545,19 +545,18 @@ fi
 $ECHO " ok"
-
 #
 # Test ctbutterfly
 # ===============
 $ECHO -n "Test ctbutterfly: "
 $ctbutterfly infile="data/crab_events.fits.gz" \
-srcmdl="data/crab.xml" \
-srcname="Crab" \
-outfile="butterfly.txt" \
-caldb="irf" \
-irf="cta_dummy_irf" \
-emin=0.1 \
-emax=100.0
+             srcmdl="data/crab.xml" \
+             srcname="Crab" \
+             outfile="butterfly.txt" \
+             caldb="irf" \
+             irf="cta_dummy_irf" \
+             emin=0.1 \
+             emax=100.0
 $ECHO -n "."
 if [ -s "butterfly.txt" ]
 then