From 9d50309d45e87b60995e06907ad4d9863a493ebc Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Sun, 26 Oct 2014 20:25:58 +0100 Subject: [PATCH] Integrate ctbutterfly tool. --- .gitignore | 2 ++ AUTHORS | 2 +- ChangeLog | 3 ++- NEWS | 18 +++++++++++++----- pyext/ctbutterfly.i | 7 +++---- src/ctbutterfly/Makefile.am | 24 ++++++++++++++++++++++++ src/ctbutterfly/ctbutterfly.cpp | 85 ++++++++++++++++++++++++++++++++++++++++++++++--------------------------------------- src/ctbutterfly/ctbutterfly.hpp | 33 ++++++++++++++++----------------- src/ctbutterfly/ctbutterfly.par | 2 +- src/ctbutterfly/main.cpp | 2 +- test/test_ctbutterfly.py | 6 +++--- test/test_ctools.sh | 15 +++++++-------- 12 files changed, 119 insertions(+), 80 deletions(-) create mode 100644 src/ctbutterfly/Makefile.am diff --git a/.gitignore b/.gitignore index 3edf77c..d1de7f9 100644 --- a/.gitignore +++ b/.gitignore @@ -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 diff --git a/AUTHORS b/AUTHORS index 7a3d3f6..fd1f82b 100644 --- a/AUTHORS +++ b/AUTHORS @@ -1,5 +1,5 @@ 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) diff --git a/ChangeLog b/ChangeLog index 60f5aee..a6bdba7 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,4 +1,4 @@ -2014-10-15 Juergen Knoedlseder +2014-10-26 Juergen Knoedlseder * 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 diff --git a/NEWS b/NEWS index e08cb1c..2185843 100644 --- a/NEWS +++ b/NEWS @@ -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 diff --git a/pyext/ctbutterfly.i b/pyext/ctbutterfly.i index 627ebce..cc9f7c9 100644 --- a/pyext/ctbutterfly.i +++ b/pyext/ctbutterfly.i @@ -1,5 +1,5 @@ /*************************************************************************** - * 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; - }; diff --git a/src/ctbutterfly/Makefile.am b/src/ctbutterfly/Makefile.am new file mode 100644 index 0000000..efd235c --- /dev/null +++ b/src/ctbutterfly/Makefile.am @@ -0,0 +1,24 @@ +## 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 diff --git a/src/ctbutterfly/ctbutterfly.cpp b/src/ctbutterfly/ctbutterfly.cpp index 1d66fbf..4844278 100644 --- a/src/ctbutterfly/ctbutterfly.cpp +++ b/src/ctbutterfly/ctbutterfly.cpp @@ -1,5 +1,5 @@ /*************************************************************************** - * 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& 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(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(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,7 +353,10 @@ void ctbutterfly::run(void) /***********************************************************************//** - * @brief Save maps + * @brief Save butterfly + * + * @exception GException::invalid_value + * File exists already and clobber is false * * This method saves the butterfly into an ascii file ***************************************************************************/ @@ -367,32 +371,36 @@ void ctbutterfly::save(void) // Get output filename m_outfile = (*this)["outfile"].filename(); + // If file exists already then delete it in case the clobber parameter + // is true 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); + std::string msg = "Output file \""+m_outfile+"\" already exists " + "and clobber parameters is set to \"false\"."; + throw GException::invalid_value(G_SAVE, msg); } } - // Create file instance + // Create file instance std::fstream file; - file.open(m_outfile, std::fstream::out | std::fstream::app); + // Open file + file.open(m_outfile.c_str(), 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; + file << "# Energy [MeV] | flux [ph/cm2/s/MeV] | e_flux [ph/cm2/s/MeV]"; + file << std::endl; // Write values to file for (int i = 0; i < m_fluxes.size(); ++i) { - - file << m_energies[i] << " " << m_fluxes[i] << " " << m_errors[i] << std::endl; + file << m_energies[i] << " " << m_fluxes[i] << " " << m_errors[i]; + file << std::endl; } - // close file + // Close file file.close(); // Return @@ -413,15 +421,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 +445,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 +543,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 +569,3 @@ void ctbutterfly::get_parameters(void) // Return return; } - diff --git a/src/ctbutterfly/ctbutterfly.hpp b/src/ctbutterfly/ctbutterfly.hpp index 81367c5..329175f 100644 --- a/src/ctbutterfly/ctbutterfly.hpp +++ b/src/ctbutterfly/ctbutterfly.hpp @@ -1,5 +1,5 @@ /*************************************************************************** - * 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 m_energies; //!< energy values for storage - std::vector m_fluxes; //!< flux values per energy bin - std::vector m_errors; //!< flux errors per energy bin + GObservations m_obs; //!< Observation container + GMatrixSparse m_covariance; //!< Covariance matrix + std::vector m_energies; //!< Energy values for storage + std::vector m_fluxes; //!< Flux values per energy bin + std::vector m_errors; //!< Flux errors per energy bin }; /***********************************************************************//** diff --git a/src/ctbutterfly/ctbutterfly.par b/src/ctbutterfly/ctbutterfly.par index dedb4af..7cc0c48 100644 --- a/src/ctbutterfly/ctbutterfly.par +++ b/src/ctbutterfly/ctbutterfly.par @@ -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" diff --git a/src/ctbutterfly/main.cpp b/src/ctbutterfly/main.cpp index 23f5f11..920caa1 100644 --- a/src/ctbutterfly/main.cpp +++ b/src/ctbutterfly/main.cpp @@ -1,5 +1,5 @@ /*************************************************************************** - * ctbutterfly - TS map calculation tool * + * ctbutterfly - TS map calculation tool * * ----------------------------------------------------------------------- * * copyright (C) 2014 by Michael Mayer * * ----------------------------------------------------------------------- * diff --git a/test/test_ctbutterfly.py b/test/test_ctbutterfly.py index 3c11324..a3e160b 100755 --- a/test/test_ctbutterfly.py +++ b/test/test_ctbutterfly.py @@ -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. diff --git a/test/test_ctools.sh b/test/test_ctools.sh index 3577bd3..146a395 100755 --- a/test/test_ctools.sh +++ b/test/test_ctools.sh @@ -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 -- libgit2 0.21.4