Jürgen Knödlseder / ctools

Commit 9d50309d45e87b60995e06907ad4d9863a493ebc

Authored by Jürgen Knödlseder
1 parent 89895869

Integrate ctbutterfly tool.

Inline Side-by-side
Showing 12 changed files with 119 additions and 80 deletions
.gitignore
  View file @9d50309
... ... @@ -57,6 +57,7 @@ src/ctexpcube/ctexpcube
57 57 src/ctpsfcube/ctpsfcube
58 58 src/ctbkgcube/ctbkgcube
59 59 src/ctcubemask/ctcubemask
  60 +src/ctbutterfly/ctbutterfly
60 61 test/pfiles/
61 62 test/*.fits
62 63 test/*.dat
... ... @@ -64,6 +65,7 @@ test/*.xml
64 65 test/*.trs
65 66 test/test_cscripts.sh
66 67 test/reports/*.xml
  68 +test/butterfly.txt
67 69 examples/*.fits
68 70 examples/*.xml
69 71 examples/*.dat
... ...
1 1 Juergen Knoedlseder (jurgen.knodlseder@irap.omp.eu)
2 2 Jean-Baptiste Cayrou (cayrou@etud.insa-toulouse.fr)
3 3 Christoph Deil (Christoph.Deil@mpi-hd.mpg.de)
4   -Michael Mayer (michael.mayer@desy.de)
  4 +Michael Mayer (michael.mayer@physik.hu-berlin.de)
5 5 Chia-Chun Lu (Chia-Chun.Lu@mpi-hd.mpg.de)
... ...
ChangeLog
  View file @9d50309
1   -2014-10-15 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
  1 +2014-10-26 Juergen Knoedlseder <jurgen.knodlseder@irap.omp.eu>
2 2  
3 3 * Version 0.8.0 released
4 4 ========================
... ... @@ -10,6 +10,7 @@
10 10 Add ctbkgcube tool (generates background cube)
11 11 Add ctcubemask tool (masks event bins in event cube)
12 12 Add cttsmap tool (generates TS map)
  13 + Add ctbutterfly tool (generates butterfly diagram)
13 14 Add csresmap script (generates residual map)
14 15 Introduce ctool base class
15 16 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
2 2  
3 3 Author(s) : Juergen Knoedlseder
4 4  
5   -15 October 2014
  5 +26 October 2014
6 6  
7 7  
8 8 1. Introduction
... ... @@ -168,12 +168,20 @@ bin content to -1. Filtering can be done by energy, for a circular
168 168 Region of Interest (ROI) or using ds9 exclusion regions.
169 169  
170 170  
171   -16. Examples
  171 +16. ctbutterfly - Generate butterfly diagrams
  172 +=============================================
  173 +The first version of the ctbutterfly tool is included in the release.
  174 +This tools derives a butterfly diagram from the fitted model parameters
  175 +for one specified source. The result will be written into an ASCII
  176 +file.
  177 +
  178 +
  179 +17. Examples
172 180 ------------
173 181 None
174 182  
175 183  
176   -17. Scripts
  184 +18. Scripts
177 185 -----------
178 186 Add logL and event statistics to cspull output file.
179 187  
... ... @@ -196,11 +204,11 @@ Add the csresmap script that uses ctbin and ctmodel to generate a
196 204 residual counts cube for a given model.
197 205  
198 206  
199   -18. caldb
  207 +19. caldb
200 208 ---------
201 209 None
202 210  
203 211  
204   -19. models
  212 +20. models
205 213 ----------
206 214 None
... ...
pyext/ctbutterfly.i
  View file @9d50309
1 1 /***************************************************************************
2   - * ctbutterfly - butterfly calculation tool *
  2 + * ctbutterfly - butterfly calculation tool *
3 3 * ----------------------------------------------------------------------- *
4 4 * copyright (C) 2014 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
... ... @@ -20,7 +20,7 @@
20 20 ***************************************************************************/
21 21 /**
22 22 * @file ctbutterfly.i
23   - * @brief tool for butterfly calculations
  23 + * @brief Butterfly calculation tool
24 24 * @author Michael Mayer
25 25 */
26 26 %{
... ... @@ -32,7 +32,7 @@
32 32 /***********************************************************************//**
33 33 * @class ctbutterfly
34 34 *
35   - * @brief butterfly calculation tool
  35 + * @brief Butterfly calculation tool
36 36 ***************************************************************************/
37 37 class ctbutterfly : public ctool {
38 38 public:
... ... @@ -48,7 +48,6 @@ public:
48 48 void run(void);
49 49 void save(void);
50 50 const GObservations& obs(void) const;
51   -
52 51 };
53 52  
54 53  
... ...
src/ctbutterfly/Makefile.am 0 → 100644
  View file @9d50309
  1 +## Process this file with automake to produce Makefile.in
  2 +
  3 +# Set include directories
  4 +INCLUDES = -I$(top_srcdir)/src/support
  5 +
  6 +# Program name
  7 +bin_PROGRAMS = ctbutterfly
  8 +
  9 +# Sources
  10 +sources = ctbutterfly.cpp
  11 +
  12 +# Libtool library
  13 +noinst_LTLIBRARIES = libctbutterfly.la
  14 +libctbutterfly_la_SOURCES = $(sources)
  15 +
  16 +# Program definition
  17 +ctbutterfly_SOURCES = main.cpp
  18 +ctbutterfly_LDFLAGS = $(LIB_LINK_OPTION) @LDFLAGS@
  19 +ctbutterfly_CPPFLAGS = @CPPFLAGS@
  20 +ctbutterfly_LDADD = libctbutterfly.la $(top_srcdir)/src/support/libsupport.la
  21 +
  22 +# Parameter files
  23 +pfilesdir = $(prefix)/syspfiles
  24 +pfiles_DATA = ctbutterfly.par
... ...
src/ctbutterfly/ctbutterfly.cpp
  View file @9d50309
1 1 /***************************************************************************
2   - * ctbutterfly - butterfly calculation tool *
  2 + * ctbutterfly - butterfly calculation tool *
3 3 * ----------------------------------------------------------------------- *
4 4 * copyright (C) 2014 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
... ... @@ -34,8 +34,8 @@
34 34 #include "GTools.hpp"
35 35  
36 36 /* __ Method name definitions ____________________________________________ */
37   -#define G_GET_PARAMETERS "ctbutterfly::get_parameters()"
38   -#define G_SAVE "ctbutterfly::save()"
  37 +#define G_GET_PARAMETERS "ctbutterfly::get_parameters()"
  38 +#define G_SAVE "ctbutterfly::save()"
39 39  
40 40 /* __ Debug definitions __________________________________________________ */
41 41  
... ... @@ -70,7 +70,7 @@ ctbutterfly::ctbutterfly(void) : ctool(CTBUTTERFLY_NAME, CTBUTTERFLY_VERSION)
70 70 * observations container.
71 71 ***************************************************************************/
72 72 ctbutterfly::ctbutterfly(const GObservations& obs) :
73   - ctool(CTBUTTERFLY_NAME, CTBUTTERFLY_VERSION)
  73 + ctool(CTBUTTERFLY_NAME, CTBUTTERFLY_VERSION)
74 74 {
75 75 // Initialise members
76 76 init_members();
... ... @@ -91,7 +91,7 @@ ctbutterfly::ctbutterfly(const GObservations&amp; obs) :
91 91 * @param[in] argv Array of command line arguments.
92 92 ***************************************************************************/
93 93 ctbutterfly::ctbutterfly(int argc, char *argv[]) :
94   - ctool(CTBUTTERFLY_NAME, CTBUTTERFLY_VERSION, argc, argv)
  94 + ctool(CTBUTTERFLY_NAME, CTBUTTERFLY_VERSION, argc, argv)
95 95 {
96 96 // Initialise members
97 97 init_members();
... ... @@ -240,11 +240,13 @@ void ctbutterfly::run(void)
240 240 log.header1("Compute covariance matrix");
241 241 }
242 242  
243   - GOptimizerPars pars = models.pars();
  243 + // Evaluate curvature matrix at the actual parameters
  244 + GOptimizerPars pars = models.pars();
244 245 GObservations::likelihood likelihood = m_obs.function();
245 246 likelihood.eval(pars);
246   - GMatrixSparse* curvature = dynamic_cast<GMatrixSparse*>(likelihood.curvature());
247   - m_covariance = curvature->invert();
  247 +
  248 + // Get covariance matrix by inverting the curvature matrix
  249 + m_covariance = likelihood.curvature()->invert();
248 250  
249 251 }
250 252  
... ... @@ -263,6 +265,7 @@ void ctbutterfly::run(void)
263 265 // Loop over energy bins positions
264 266 for (int i = 0; i < m_ebounds.size(); ++i) {
265 267  
  268 + // Initialise number of gradients
266 269 int num_gradient = 0;
267 270  
268 271 // Get the energy of current bin
... ... @@ -270,8 +273,8 @@ void ctbutterfly::run(void)
270 273  
271 274 // Header for verbose logging
272 275 if (logExplicit()) {
273   - std::string msg = "Computing butterfly for bin number "+gammalib::str(i)+
274   - " at "+energy.print();
  276 + std::string msg = "Computing butterfly for bin number "+
  277 + gammalib::str(i)+" at "+energy.print();
275 278 log << std::endl;
276 279 log.header2(msg);
277 280 }
... ... @@ -280,7 +283,7 @@ void ctbutterfly::run(void)
280 283 double model_flux = 0.0;
281 284  
282 285 // Loop over models
283   - for (int j = 0; j < models.size(); ++j ) {
  286 + for (int j = 0; j < models.size(); ++j) {
284 287  
285 288 // Check wether model is a skymodel
286 289 GModelSky* skymodel = dynamic_cast<GModelSky*>(models[j]);
... ... @@ -296,17 +299,15 @@ void ctbutterfly::run(void)
296 299  
297 300 // Set flux value of source of interest
298 301 if (skymodel->name() == m_srcname) {
299   - model_flux = spectral->eval_gradients(energy, time);
  302 + model_flux = spectral->eval_gradients(energy, time);
300 303 } // endif: model was source of interest
301 304  
302 305 // Loop over model parameters, get gradients
303 306 // and assign them to the vector
304 307 for (int k = 0; k < spectral->size(); ++k) {
305   -
306   - grad[num_gradient] = (*spectral)[k].gradient();
307   - num_gradient += 1;
308   -
309   - }// endfor: loop over spectral parameters
  308 + grad[num_gradient] = (*spectral)[k].gradient();
  309 + num_gradient++;
  310 + } // endfor: looped over spectral parameters
310 311  
311 312 // Skip temporal models
312 313 num_gradient += skymodel->temporal()->size();
... ... @@ -352,7 +353,10 @@ void ctbutterfly::run(void)
352 353  
353 354  
354 355 /***********************************************************************//**
355   - * @brief Save maps
  356 + * @brief Save butterfly
  357 + *
  358 + * @exception GException::invalid_value
  359 + * File exists already and clobber is false
356 360 *
357 361 * This method saves the butterfly into an ascii file
358 362 ***************************************************************************/
... ... @@ -367,32 +371,36 @@ void ctbutterfly::save(void)
367 371 // Get output filename
368 372 m_outfile = (*this)["outfile"].filename();
369 373  
  374 + // If file exists already then delete it in case the clobber parameter
  375 + // is true
370 376 if (gammalib::file_exists(m_outfile)) {
371   -
372 377 if ((*this)["clobber"].boolean()) {
373 378 std::remove(m_outfile.c_str());
374 379 }
375 380 else {
376   - std::string msg = "File \""+m_outfile+"\" already exists and clobber is set to \"false\"";
377   - throw GException::file_open_error(G_SAVE,msg);
  381 + std::string msg = "Output file \""+m_outfile+"\" already exists "
  382 + "and clobber parameters is set to \"false\".";
  383 + throw GException::invalid_value(G_SAVE, msg);
378 384 }
379 385 }
380 386  
381   - // Create file instance
  387 + // Create file instance
382 388 std::fstream file;
383 389  
384   - file.open(m_outfile, std::fstream::out | std::fstream::app);
  390 + // Open file
  391 + file.open(m_outfile.c_str(), std::fstream::out | std::fstream::app);
385 392  
386 393 // Write short header
387   - file << "# Energy [MeV] | flux [ph/cm2/s/MeV] | e_flux [ph/cm2/s/MeV]" << std::endl;
  394 + file << "# Energy [MeV] | flux [ph/cm2/s/MeV] | e_flux [ph/cm2/s/MeV]";
  395 + file << std::endl;
388 396  
389 397 // Write values to file
390 398 for (int i = 0; i < m_fluxes.size(); ++i) {
391   -
392   - file << m_energies[i] << " " << m_fluxes[i] << " " << m_errors[i] << std::endl;
  399 + file << m_energies[i] << " " << m_fluxes[i] << " " << m_errors[i];
  400 + file << std::endl;
393 401 }
394 402  
395   - // close file
  403 + // Close file
396 404 file.close();
397 405  
398 406 // Return
... ... @@ -413,15 +421,15 @@ void ctbutterfly::init_members(void)
413 421 {
414 422 // Initialise members
415 423 m_infile.clear();
416   - m_outfile.clear();
417 424 m_srcname.clear();
  425 + m_outfile.clear();
418 426 m_ebounds.clear();
419 427  
420 428 // Initialise protected members
421 429 m_obs.clear();
422 430 m_covariance.clear();
423   - m_fluxes.clear();
424 431 m_energies.clear();
  432 + m_fluxes.clear();
425 433 m_errors.clear();
426 434  
427 435 // Return
... ... @@ -437,17 +445,17 @@ void ctbutterfly::init_members(void)
437 445 void ctbutterfly::copy_members(const ctbutterfly& app)
438 446 {
439 447 // Copy attributes
440   - m_infile = app.m_infile;
441   - m_outfile = app.m_outfile;
442   - m_srcname = app.m_srcname;
443   - m_ebounds = app.m_ebounds;
  448 + m_infile = app.m_infile;
  449 + m_srcname = app.m_srcname;
  450 + m_outfile = app.m_outfile;
  451 + m_ebounds = app.m_ebounds;
444 452  
445 453 // Copy protected members
446 454 m_obs = app.m_obs;
447 455 m_covariance = app.m_covariance;
448   - m_fluxes = app.m_fluxes;
449   - m_energies = app.m_energies;
450   - m_errors = app.m_errors;
  456 + m_energies = app.m_energies;
  457 + m_fluxes = app.m_fluxes;
  458 + m_errors = app.m_errors;
451 459  
452 460 // Return
453 461 return;
... ... @@ -535,11 +543,11 @@ void ctbutterfly::get_parameters(void)
535 543 if (!m_obs.models().contains(m_srcname)) {
536 544 std::string msg = "Source \""+m_srcname+"\" not found in model "
537 545 "container. Please add a source with that name "
538   - "or check for a possible typos.";
  546 + "or check for possible typos.";
539 547 throw GException::invalid_value(G_GET_PARAMETERS, msg);
540 548 }
541 549  
542   - // get binning information
  550 + // Get energy binning information
543 551 m_ebounds = get_ebounds();
544 552  
545 553 // Get matrix file name and load if possible
... ... @@ -561,4 +569,3 @@ void ctbutterfly::get_parameters(void)
561 569 // Return
562 570 return;
563 571 }
564   -
... ...
src/ctbutterfly/ctbutterfly.hpp
  View file @9d50309
1 1 /***************************************************************************
2   - * ctbutterfly - butterfly calculation tool *
  2 + * ctbutterfly - butterfly calculation tool *
3 3 * ----------------------------------------------------------------------- *
4 4 * copyright (C) 2014 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
... ... @@ -40,17 +40,18 @@
40 40 /***********************************************************************//**
41 41 * @class ctbutterfly
42 42 *
43   - * @brief butterfly calculation tool
  43 + * @brief Butterfly calculation tool
44 44 *
45 45 * This class computes the confidence interval of a fitted spectrum
46   - * (butterfly) for a given set of observations
  46 + * (butterfly) for a given set of observations.
47 47 *
48 48 * The class operates on predefined observation containers, an individual
49 49 * event list or an observation definition XML file.
50 50 *
51   - * During the computation the covariance matrix of the fit is used to propagate
52   - * uncertainties and their correlations through the entire energy range.
53   - * The output is saved as an ascii files containing the confidence band boundaries
  51 + * During the computation the covariance matrix of the fit is used to
  52 + * propagate uncertainties and their correlations through the entire energy
  53 + * range. The output is saved as an ascii files containing the confidence
  54 + * band boundaries
54 55 ***************************************************************************/
55 56 class ctbutterfly : public ctool {
56 57  
... ... @@ -79,19 +80,17 @@ protected:
79 80 void get_parameters(void);
80 81  
81 82 // User parameters
82   - std::string m_infile; //!< Input file
83   - std::string m_srcname; //!< Name of source to compute butterfly
84   - std::string m_outfile; //!< Output ascii file
85   -
86   - GEbounds m_ebounds; //!< Energy binning definition
  83 + std::string m_infile; //!< Input file
  84 + std::string m_srcname; //!< Name of source to compute butterfly
  85 + std::string m_outfile; //!< Output ascii file
  86 + GEbounds m_ebounds; //!< Energy binning definition
87 87  
88 88 // Protected members
89   - GObservations m_obs; //!< Observation container
90   - GMatrixSparse m_covariance; //!< covariance matrix
91   -
92   - std::vector<double> m_energies; //!< energy values for storage
93   - std::vector<double> m_fluxes; //!< flux values per energy bin
94   - std::vector<double> m_errors; //!< flux errors per energy bin
  89 + GObservations m_obs; //!< Observation container
  90 + GMatrixSparse m_covariance; //!< Covariance matrix
  91 + std::vector<double> m_energies; //!< Energy values for storage
  92 + std::vector<double> m_fluxes; //!< Flux values per energy bin
  93 + std::vector<double> m_errors; //!< Flux errors per energy bin
95 94 };
96 95  
97 96 /***********************************************************************//**
... ...
src/ctbutterfly/ctbutterfly.par
  View file @9d50309
... ... @@ -19,7 +19,7 @@
19 19 #===================
20 20 infile, f, a, "events.fits",,,"Input event list, cube or observation definition file"
21 21 srcmdl, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
22   -srcname, s, a, "CrabNebula",,,"Source of interest"
  22 +srcname, s, a, "Crab",,,"Source of interest"
23 23 caldb, s, a, "dummy",,, "Calibration database"
24 24 irf, s, a, "cta_dummy_irf",,, "Instrument response function"
25 25 outfile, f, a, "butterfly.txt",,, "Output ascii file"
... ...
src/ctbutterfly/main.cpp
  View file @9d50309
1 1 /***************************************************************************
2   - * ctbutterfly - TS map calculation tool *
  2 + * ctbutterfly - TS map calculation tool *
3 3 * ----------------------------------------------------------------------- *
4 4 * copyright (C) 2014 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
... ...
test/test_ctbutterfly.py
  View file @9d50309
... ... @@ -2,7 +2,7 @@
2 2 # ==========================================================================
3 3 # This scripts performs unit tests for the ctbutterfly tool.
4 4 #
5   -# Copyright (C) 2014 Juergen Knoedlseder
  5 +# Copyright (C) 2014 Michal Mayer
6 6 #
7 7 # This program is free software: you can redistribute it and/or modify
8 8 # it under the terms of the GNU General Public License as published by
... ... @@ -22,9 +22,9 @@ import gammalib
22 22 import ctools
23 23  
24 24  
25   -# ========================== #
  25 +# =============================== #
26 26 # Test class for ctbutterfly tool #
27   -# ========================== #
  27 +# =============================== #
28 28 class Test(gammalib.GPythonTestSuite):
29 29 """
30 30 Test class for ctbutterfly tool.
... ...
test/test_ctools.sh
  View file @9d50309
... ... @@ -545,19 +545,18 @@ fi
545 545 $ECHO " ok"
546 546  
547 547  
548   -
549 548 #
550 549 # Test ctbutterfly
551 550 # ===============
552 551 $ECHO -n "Test ctbutterfly: "
553 552 $ctbutterfly infile="data/crab_events.fits.gz" \
554   -srcmdl="data/crab.xml" \
555   -srcname="Crab" \
556   -outfile="butterfly.txt" \
557   -caldb="irf" \
558   -irf="cta_dummy_irf" \
559   -emin=0.1 \
560   -emax=100.0
  553 + srcmdl="data/crab.xml" \
  554 + srcname="Crab" \
  555 + outfile="butterfly.txt" \
  556 + caldb="irf" \
  557 + irf="cta_dummy_irf" \
  558 + emin=0.1 \
  559 + emax=100.0
561 560 $ECHO -n "."
562 561 if [ -s "butterfly.txt" ]
563 562 then
... ...