Michael Mayer / ctools

Commit 22700ce9e48afc2ff08259bb546466f682fc15d8

Authored by Jürgen Knödlseder
2 parents 9a22fd59 3d171380

Merge branch '1051-create-ctbutterfly' into integration 

Conflicts:
	ChangeLog
	NEWS
Inline Side-by-side
Showing 18 changed files with 1015 additions and 10 deletions
.gitignore
  View file @22700ce
... ... @@ -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 @22700ce
1   -2014-10-24 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   -24 October 2014
  5 +26 October 2014
6 6  
7 7  
8 8 1. Introduction
... ... @@ -183,12 +183,20 @@ bin content to -1. Filtering can be done by energy, for a circular
183 183 Region of Interest (ROI) or using ds9 exclusion regions.
184 184  
185 185  
186   -16. Examples
  186 +16. ctbutterfly - Generate butterfly diagrams
  187 +=============================================
  188 +The first version of the ctbutterfly tool is included in the release.
  189 +This tools derives a butterfly diagram from the fitted model parameters
  190 +for one specified source. The result will be written into an ASCII
  191 +file.
  192 +
  193 +
  194 +17. Examples
187 195 ------------
188 196 None
189 197  
190 198  
191   -17. Scripts
  199 +18. Scripts
192 200 -----------
193 201 Add logL and event statistics to cspull output file.
194 202  
... ... @@ -211,11 +219,11 @@ Add the csresmap script that uses ctbin and ctmodel to generate a
211 219 residual counts cube for a given model.
212 220  
213 221  
214   -18. caldb
  222 +19. caldb
215 223 ---------
216 224 None
217 225  
218 226  
219   -19. models
  227 +20. models
220 228 ----------
221 229 None
... ...
configure.ac
  View file @22700ce
... ... @@ -887,6 +887,7 @@ AC_CONFIG_FILES([Makefile
887 887 src/ctbkgcube/Makefile
888 888 src/ctcubemask/Makefile
889 889 src/cttsmap/Makefile
  890 + src/ctbutterfly/Makefile
890 891 src/support/Makefile
891 892 src/ctools-setup
892 893 test/Makefile
... ...
pyext/ctbutterfly.i 0 → 100644
  View file @22700ce
  1 +/***************************************************************************
  2 + * ctbutterfly - butterfly calculation tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2014 by Michael Mayer *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file ctbutterfly.i
  23 + * @brief Butterfly calculation tool
  24 + * @author Michael Mayer
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "ctbutterfly.hpp"
  29 +%}
  30 +
  31 +
  32 +/***********************************************************************//**
  33 + * @class ctbutterfly
  34 + *
  35 + * @brief Butterfly calculation tool
  36 + ***************************************************************************/
  37 +class ctbutterfly : public ctool {
  38 +public:
  39 + // Constructors and destructors
  40 + ctbutterfly(void);
  41 + explicit ctbutterfly(const GObservations& obs);
  42 + ctbutterfly(int argc, char *argv[]);
  43 + ctbutterfly(const ctbutterfly& app);
  44 + virtual ~ctbutterfly(void);
  45 +
  46 + // Methods
  47 + void clear(void);
  48 + void run(void);
  49 + void save(void);
  50 + const GObservations& obs(void) const;
  51 +};
  52 +
  53 +
  54 +/***********************************************************************//**
  55 + * @brief butterfly calculation tool Python extensions
  56 + ***************************************************************************/
  57 +%extend ctbutterfly {
  58 + ctbutterfly copy() {
  59 + return (*self);
  60 + }
  61 +}
... ...
pyext/ctools/ctools.i
  View file @22700ce
... ... @@ -57,3 +57,4 @@
57 57 %include "ctbkgcube.i"
58 58 %include "ctcubemask.i"
59 59 %include "cttsmap.i"
  60 +%include "ctbutterfly.i"
... ...
src/Makefile.am
  View file @22700ce
... ... @@ -16,7 +16,8 @@ SUBDIRS = support \
16 16 ctpsfcube \
17 17 ctbkgcube \
18 18 ctcubemask \
19   - cttsmap
  19 + cttsmap \
  20 + ctbutterfly
20 21  
21 22 # Scripts to be distributed
22 23 scripts = ctools-setup \
... ... @@ -39,7 +40,8 @@ libctools_la_LIBADD = support/libsupport.la \
39 40 ctpsfcube/libctpsfcube.la \
40 41 ctbkgcube/libctbkgcube.la \
41 42 ctcubemask/libctcubemask.la \
42   - cttsmap/libcttsmap.la
  43 + cttsmap/libcttsmap.la \
  44 + ctbutterfly/libctbutterfly.la
43 45  
44 46  
45 47 # Distribute scripts
... ...
src/ctbutterfly/Makefile.am 0 → 100644
  View file @22700ce
  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 0 → 100644
  View file @22700ce
  1 +/***************************************************************************
  2 + * ctbutterfly - butterfly calculation tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2014 by Michael Mayer *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file ctbutterfly.hpp
  23 + * @brief Butterfly computation tool
  24 + * @author Michael Mayer
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include <cstdio>
  32 +#include <fstream>
  33 +#include "ctbutterfly.hpp"
  34 +#include "GTools.hpp"
  35 +
  36 +/* __ Method name definitions ____________________________________________ */
  37 +#define G_GET_PARAMETERS "ctbutterfly::get_parameters()"
  38 +#define G_SAVE "ctbutterfly::save()"
  39 +
  40 +/* __ Debug definitions __________________________________________________ */
  41 +
  42 +/* __ Coding definitions _________________________________________________ */
  43 +
  44 +
  45 +/*==========================================================================
  46 + = =
  47 + = Constructors/destructors =
  48 + = =
  49 + ==========================================================================*/
  50 +
  51 +/***********************************************************************//**
  52 + * @brief Void constructor
  53 + ***************************************************************************/
  54 +ctbutterfly::ctbutterfly(void) : ctool(CTBUTTERFLY_NAME, CTBUTTERFLY_VERSION)
  55 +{
  56 + // Initialise members
  57 + init_members();
  58 +
  59 + // Return
  60 + return;
  61 +}
  62 +
  63 +
  64 +/***********************************************************************//**
  65 + * @brief Observations constructor
  66 + *
  67 + * param[in] obs Observation container.
  68 + *
  69 + * This method creates an instance of the class by copying an existing
  70 + * observations container.
  71 + ***************************************************************************/
  72 +ctbutterfly::ctbutterfly(const GObservations& obs) :
  73 + ctool(CTBUTTERFLY_NAME, CTBUTTERFLY_VERSION)
  74 +{
  75 + // Initialise members
  76 + init_members();
  77 +
  78 + // Set observations
  79 + m_obs = obs;
  80 +
  81 + // Return
  82 + return;
  83 +}
  84 +
  85 +
  86 +
  87 +/***********************************************************************//**
  88 + * @brief Command line constructor
  89 + *
  90 + * @param[in] argc Number of arguments in command line.
  91 + * @param[in] argv Array of command line arguments.
  92 + ***************************************************************************/
  93 +ctbutterfly::ctbutterfly(int argc, char *argv[]) :
  94 + ctool(CTBUTTERFLY_NAME, CTBUTTERFLY_VERSION, argc, argv)
  95 +{
  96 + // Initialise members
  97 + init_members();
  98 +
  99 + // Return
  100 + return;
  101 +}
  102 +
  103 +
  104 +/***********************************************************************//**
  105 + * @brief Copy constructor
  106 + *
  107 + * @param[in] app Application.
  108 + ***************************************************************************/
  109 +ctbutterfly::ctbutterfly(const ctbutterfly& app) : ctool(app)
  110 +{
  111 + // Initialise members
  112 + init_members();
  113 +
  114 + // Copy members
  115 + copy_members(app);
  116 +
  117 + // Return
  118 + return;
  119 +}
  120 +
  121 +
  122 +/***********************************************************************//**
  123 + * @brief Destructor
  124 + ***************************************************************************/
  125 +ctbutterfly::~ctbutterfly(void)
  126 +{
  127 + // Free members
  128 + free_members();
  129 +
  130 + // Return
  131 + return;
  132 +}
  133 +
  134 +
  135 +/*==========================================================================
  136 + = =
  137 + = Operators =
  138 + = =
  139 + ==========================================================================*/
  140 +
  141 +/***********************************************************************//**
  142 + * @brief Assignment operator
  143 + *
  144 + * @param[in] app Application.
  145 + * @return Application.
  146 + ***************************************************************************/
  147 +ctbutterfly& ctbutterfly::operator=(const ctbutterfly& app)
  148 +{
  149 + // Execute only if object is not identical
  150 + if (this != &app) {
  151 +
  152 + // Copy base class members
  153 + this->ctool::operator=(app);
  154 +
  155 + // Free members
  156 + free_members();
  157 +
  158 + // Initialise members
  159 + init_members();
  160 +
  161 + // Copy members
  162 + copy_members(app);
  163 +
  164 + } // endif: object was not identical
  165 +
  166 + // Return this object
  167 + return *this;
  168 +}
  169 +
  170 +
  171 +/*==========================================================================
  172 + = =
  173 + = Public methods =
  174 + = =
  175 + ==========================================================================*/
  176 +
  177 +/***********************************************************************//**
  178 + * @brief Clear instance
  179 + ***************************************************************************/
  180 +void ctbutterfly::clear(void)
  181 +{
  182 + // Free members
  183 + free_members();
  184 + this->ctool::free_members();
  185 + this->GApplication::free_members();
  186 +
  187 + // Initialise members
  188 + this->GApplication::init_members();
  189 + this->ctool::init_members();
  190 + init_members();
  191 +
  192 + // Return
  193 + return;
  194 +}
  195 +
  196 +
  197 +/***********************************************************************//**
  198 + * @brief Computes the butterfly
  199 + *
  200 + * This method takes the spectral fit and its covariance matrix to compute
  201 + * the confidence band via Gaussian error propagation
  202 + ***************************************************************************/
  203 +void ctbutterfly::run(void)
  204 +{
  205 + // If we're in debug mode then all output is also dumped on the screen
  206 + if (logDebug()) {
  207 + log.cout(true);
  208 + }
  209 +
  210 + // Get task parameters
  211 + get_parameters();
  212 +
  213 + // Write parameters into logger
  214 + if (logTerse()) {
  215 + log_parameters();
  216 + log << std::endl;
  217 + }
  218 +
  219 + // Write observation(s) into logger
  220 + if (logTerse()) {
  221 + log << std::endl;
  222 + if (m_obs.size() > 1) {
  223 + log.header1("Observations");
  224 + }
  225 + else {
  226 + log.header1("Observation");
  227 + }
  228 + log << m_obs << std::endl;
  229 + }
  230 +
  231 + // Get model instance for further computations
  232 + GModels models = m_obs.models();
  233 +
  234 + // Compute covariance matrix if necessary
  235 + if (m_covariance.size() == 0) {
  236 +
  237 + // Write header
  238 + if (logTerse()) {
  239 + log << std::endl;
  240 + log.header1("Compute covariance matrix");
  241 + }
  242 +
  243 + // Evaluate curvature matrix at the actual parameters
  244 + GOptimizerPars pars = models.pars();
  245 + GObservations::likelihood likelihood = m_obs.function();
  246 + likelihood.eval(pars);
  247 +
  248 + // Get covariance matrix by inverting the curvature matrix
  249 + m_covariance = likelihood.curvature()->invert();
  250 +
  251 + }
  252 +
  253 + // Write header
  254 + if (logTerse()) {
  255 + log << std::endl;
  256 + log.header1("Generate butterfly");
  257 + }
  258 +
  259 + // Initialise dummy time to evaluate spectral model
  260 + GTime time = GTime();
  261 +
  262 + // Initialise vector of gradients
  263 + GVector grad = GVector(m_obs.models().npars());
  264 +
  265 + // Loop over energy bins positions
  266 + for (int i = 0; i < m_ebounds.size(); ++i) {
  267 +
  268 + // Initialise number of gradients
  269 + int num_gradient = 0;
  270 +
  271 + // Get the energy of current bin
  272 + GEnergy energy = m_ebounds.elogmean(i);
  273 +
  274 + // Header for verbose logging
  275 + if (logExplicit()) {
  276 + std::string msg = "Computing butterfly for bin number "+
  277 + gammalib::str(i)+" at "+energy.print();
  278 + log << std::endl;
  279 + log.header2(msg);
  280 + }
  281 +
  282 + // Initialise model flux value
  283 + double model_flux = 0.0;
  284 +
  285 + // Loop over models
  286 + for (int j = 0; j < models.size(); ++j) {
  287 +
  288 + // Check wether model is a skymodel
  289 + GModelSky* skymodel = dynamic_cast<GModelSky*>(models[j]);
  290 +
  291 + // Yes ...
  292 + if (skymodel != NULL) {
  293 +
  294 + // Skip spatial models
  295 + num_gradient += skymodel->spatial()->size();
  296 +
  297 + // Get pointer to spectral model
  298 + GModelSpectral* spectral = skymodel->spectral();
  299 +
  300 + // Set flux value of source of interest
  301 + if (skymodel->name() == m_srcname) {
  302 + model_flux = spectral->eval_gradients(energy, time);
  303 + } // endif: model was source of interest
  304 +
  305 + // Loop over model parameters, get gradients
  306 + // and assign them to the vector
  307 + for (int k = 0; k < spectral->size(); ++k) {
  308 + grad[num_gradient] = (*spectral)[k].gradient();
  309 + num_gradient++;
  310 + } // endfor: looped over spectral parameters
  311 +
  312 + // Skip temporal models
  313 + num_gradient += skymodel->temporal()->size();
  314 +
  315 + } // endif: model was sky model
  316 +
  317 + else {
  318 +
  319 + // Skip other models (e.g. GModelData instances)
  320 + num_gradient += models[j]->size();
  321 +
  322 + } // endelse: Model was not GModelSky
  323 +
  324 + } // endfor: Looped over models
  325 +
  326 + // Multiply covariance to the gradient vector
  327 + GVector vector = m_covariance * grad;
  328 +
  329 + // Get the error from the scalar product
  330 + double error = std::sqrt(grad * vector);
  331 +
  332 + // Store flux, value and energy for saving
  333 + m_fluxes.push_back(model_flux);
  334 + m_energies.push_back(energy.MeV());
  335 + m_errors.push_back(error);
  336 +
  337 + // Log information
  338 + if (logExplicit()) {
  339 + log << " Flux at ";
  340 + log << energy.print();
  341 + log << ": ";
  342 + log << model_flux;
  343 + log << " +- ";
  344 + log << error;
  345 + log << " ph/cm2/s/MeV";
  346 + }
  347 +
  348 + } //endfor: loop over energy bin
  349 +
  350 + // Return
  351 + return;
  352 +}
  353 +
  354 +
  355 +/***********************************************************************//**
  356 + * @brief Save butterfly diagram
  357 + *
  358 + * Saves the butterfly diagram into an ascii file using a column separated
  359 + * value (CSV) format with blanks as separators.
  360 + ***************************************************************************/
  361 +void ctbutterfly::save(void)
  362 +{
  363 + // Write header
  364 + if (logTerse()) {
  365 + log << std::endl;
  366 + log.header1("Save Butterfly to file");
  367 + }
  368 +
  369 + // Get output filename
  370 + m_outfile = (*this)["outfile"].filename();
  371 +
  372 + // Create CSV table with 3 columns
  373 + GCsv table(m_fluxes.size(), 3);
  374 +
  375 + // Fill CSV table
  376 + for (int i = 0; i < m_fluxes.size(); ++i) {
  377 + table.real(i, 0, m_energies[i]);
  378 + table.real(i, 1, m_fluxes[i]);
  379 + table.real(i, 2, m_errors[i]);
  380 + }
  381 +
  382 + // Save CSV table
  383 + table.save(m_outfile, " ", clobber());
  384 +
  385 + // Return
  386 + return;
  387 +}
  388 +
  389 +
  390 +/*==========================================================================
  391 + = =
  392 + = Private methods =
  393 + = =
  394 + ==========================================================================*/
  395 +
  396 +/***********************************************************************//**
  397 + * @brief Initialise class members
  398 + ***************************************************************************/
  399 +void ctbutterfly::init_members(void)
  400 +{
  401 + // Initialise members
  402 + m_infile.clear();
  403 + m_srcname.clear();
  404 + m_outfile.clear();
  405 + m_ebounds.clear();
  406 +
  407 + // Initialise protected members
  408 + m_obs.clear();
  409 + m_covariance.clear();
  410 + m_energies.clear();
  411 + m_fluxes.clear();
  412 + m_errors.clear();
  413 +
  414 + // Return
  415 + return;
  416 +}
  417 +
  418 +
  419 +/***********************************************************************//**
  420 + * @brief Copy class members
  421 + *
  422 + * @param[in] app Application.
  423 + ***************************************************************************/
  424 +void ctbutterfly::copy_members(const ctbutterfly& app)
  425 +{
  426 + // Copy attributes
  427 + m_infile = app.m_infile;
  428 + m_srcname = app.m_srcname;
  429 + m_outfile = app.m_outfile;
  430 + m_ebounds = app.m_ebounds;
  431 +
  432 + // Copy protected members
  433 + m_obs = app.m_obs;
  434 + m_covariance = app.m_covariance;
  435 + m_energies = app.m_energies;
  436 + m_fluxes = app.m_fluxes;
  437 + m_errors = app.m_errors;
  438 +
  439 + // Return
  440 + return;
  441 +}
  442 +
  443 +
  444 +/***********************************************************************//**
  445 + * @brief Delete class members
  446 + ***************************************************************************/
  447 +void ctbutterfly::free_members(void)
  448 +{
  449 + // Return
  450 + return;
  451 +}
  452 +
  453 +
  454 +/***********************************************************************//**
  455 + * @brief Get application parameters
  456 + *
  457 + * @exception GException::invalid_value
  458 + * Test source not found or no RA/DEC parameters found for test
  459 + * source.
  460 + *
  461 + * Get all task parameters from parameter file or (if required) by querying
  462 + * the user. Most parameters are only required if no observation exists so
  463 + * far in the observation container. In this case, a single CTA observation
  464 + * will be added to the container, using the definition provided in the
  465 + * parameter file.
  466 + ***************************************************************************/
  467 +void ctbutterfly::get_parameters(void)
  468 +{
  469 + // If there are no observations in container then add a single CTA
  470 + // observation using the parameters from the parameter file
  471 + if (m_obs.size() == 0) {
  472 +
  473 + // Allocate CTA observation
  474 + GCTAObservation obs;
  475 +
  476 + // Get event file name
  477 + std::string filename = (*this)["infile"].filename();
  478 +
  479 + // Try first to open as FITS file
  480 + try {
  481 +
  482 + // Load data
  483 + obs.load(filename);
  484 +
  485 + // Set response
  486 + set_obs_response(&obs);
  487 +
  488 + // Append observation to container
  489 + m_obs.append(obs);
  490 +
  491 + }
  492 +
  493 + // ... otherwise try to open as XML file
  494 + catch (GException::fits_open_error &e) {
  495 +
  496 + // Load observations from XML file
  497 + m_obs.load(filename);
  498 +
  499 + // Check if all observations have response information. If
  500 + // not, get the calibration database parameters and set
  501 + // the response properly
  502 + set_response(m_obs);
  503 +
  504 + } // endcatch: file was an XML file
  505 +
  506 + } // endif: there was no observation in the container
  507 +
  508 + // If there is are no models associated with the observations then
  509 + // load now the model definition
  510 + if (m_obs.models().size() == 0) {
  511 +
  512 + // Get models XML filename
  513 + std::string filename = (*this)["srcmdl"].filename();
  514 +
  515 + // Setup models for optimizing.
  516 + m_obs.models(GModels(filename));
  517 +
  518 + } // endif: no models were associated with observations
  519 +
  520 + // Get name of test source and check container for this name
  521 + m_srcname = (*this)["srcname"].string();
  522 + if (!m_obs.models().contains(m_srcname)) {
  523 + std::string msg = "Source \""+m_srcname+"\" not found in model "
  524 + "container. Please add a source with that name "
  525 + "or check for possible typos.";
  526 + throw GException::invalid_value(G_GET_PARAMETERS, msg);
  527 + }
  528 +
  529 + // Get energy binning information
  530 + m_ebounds = get_ebounds();
  531 +
  532 + // Get matrix file name and load if possible
  533 + std::string matrixfilename = (*this)["matrix"].filename();
  534 + if (matrixfilename != "NONE") {
  535 + std::string msg = "Loading of matrix from file not implemented yet. "
  536 + "Use filename = \"NONE\" to induce a recomputation of "
  537 + "the matrix internally.";
  538 + throw GException::feature_not_implemented(G_GET_PARAMETERS, msg);
  539 + // m_covariance.load(matrixfilename);
  540 + }
  541 +
  542 + // Optionally read ahead parameters so that they get correctly
  543 + // dumped into the log file
  544 + if (read_ahead()) {
  545 + m_outfile = (*this)["outfile"].filename();
  546 + }
  547 +
  548 + // Return
  549 + return;
  550 +}
... ...
src/ctbutterfly/ctbutterfly.hpp 0 → 100644
  View file @22700ce
  1 +/***************************************************************************
  2 + * ctbutterfly - butterfly calculation tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2014 by Michael Mayer *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file ctbutterfly.hpp
  23 + * @brief Butterfly calculation tool interface definition
  24 + * @author Michael Mayer
  25 + */
  26 +
  27 +#ifndef CTBUTTERFLY_HPP
  28 +#define CTBUTTERFLY_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include "GammaLib.hpp"
  32 +#include "GCTALib.hpp"
  33 +#include "ctool.hpp"
  34 +
  35 +/* __Definitions _________________________________________________________ */
  36 +#define CTBUTTERFLY_NAME "ctbutterfly"
  37 +#define CTBUTTERFLY_VERSION "00-01-00"
  38 +
  39 +
  40 +/***********************************************************************//**
  41 + * @class ctbutterfly
  42 + *
  43 + * @brief Butterfly calculation tool
  44 + *
  45 + * This class computes the confidence interval of a fitted spectrum
  46 + * (butterfly) for a given set of observations.
  47 + *
  48 + * The class operates on predefined observation containers, an individual
  49 + * event list or an observation definition XML file.
  50 + *
  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
  55 + ***************************************************************************/
  56 +class ctbutterfly : public ctool {
  57 +
  58 +public:
  59 + // Constructors and destructors
  60 + ctbutterfly(void);
  61 + explicit ctbutterfly(const GObservations& obs);
  62 + ctbutterfly(int argc, char *argv[]);
  63 + ctbutterfly(const ctbutterfly& app);
  64 + virtual ~ctbutterfly(void);
  65 +
  66 + // Operators
  67 + ctbutterfly& operator=(const ctbutterfly& app);
  68 +
  69 + // Methods
  70 + void clear(void);
  71 + void run(void);
  72 + void save(void);
  73 + const GObservations& obs(void) const;
  74 +
  75 +protected:
  76 + // Protected methods
  77 + void init_members(void);
  78 + void copy_members(const ctbutterfly& app);
  79 + void free_members(void);
  80 + void get_parameters(void);
  81 +
  82 + // User parameters
  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 +
  88 + // Protected members
  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
  94 +};
  95 +
  96 +/***********************************************************************//**
  97 + * @brief Return observation container
  98 + *
  99 + * @return Reference to observation container
  100 + ***************************************************************************/
  101 +inline
  102 +const GObservations& ctbutterfly::obs(void) const
  103 +{
  104 + return m_obs;
  105 +}
  106 +
  107 +#endif /* CTBUTTERFLY_HPP */
... ...
src/ctbutterfly/ctbutterfly.par 0 → 100644
  View file @22700ce
  1 +#############################################################################
  2 +# #
  3 +# CTA Science Analysis Tools #
  4 +# #
  5 +# ------------------------------------------------------------------------- #
  6 +# #
  7 +# File: ctbutterfly #
  8 +# Version: 00-01-00 #
  9 +# #
  10 +# Author: Michael Mayer #
  11 +# michael.mayer@physik.hu-berlin.de #
  12 +# HU Berlin #
  13 +# #
  14 +# Purpose: Parameter file for the ctbutterfly application #
  15 +# #
  16 +#############################################################################
  17 +#
  18 +# General parameters
  19 +#===================
  20 +infile, f, a, "events.fits",,,"Input event list, cube or observation definition file"
  21 +srcmdl, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
  22 +srcname, s, a, "Crab",,,"Source of interest"
  23 +caldb, s, a, "dummy",,, "Calibration database"
  24 +irf, s, a, "cta_dummy_irf",,, "Instrument response function"
  25 +outfile, f, a, "butterfly.txt",,, "Output ascii file"
  26 +matrix, f, h, "NONE",,, "Input covariance Matrix FITS file"
  27 +
  28 +#
  29 +# Energy binning parameters
  30 +#==========================
  31 +ebinalg, s, h, "LOG", FILE|LIN|LOG,,"Algorithm for defining energy bins"
  32 +emin, r, a, 0.1,,,"Start value for first energy bin in TeV"
  33 +emax, r, a, 100.0,,,"Stop value for last energy bin in TeV"
  34 +enumbins, i, h, 100,,,"Number of energy bins"
  35 +ebinfile, f, h, "NONE",,,"Name of the file containing the energy bin definition"
  36 +
  37 +#
  38 +# Standard parameters
  39 +#====================
  40 +chatter, i, h, 2,0,4, "Chattiness of output"
  41 +clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
  42 +debug, b, h, no,,, "Debugging mode activated"
  43 +mode, s, h, "ql",,, "Mode of automatic parameters"
  44 +logfile, f, h, "ctbutterfly.log",,, "Log filename"
... ...
src/ctbutterfly/main.cpp 0 → 100644
  View file @22700ce
  1 +/***************************************************************************
  2 + * ctbutterfly - TS map calculation tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2014 by Michael Mayer *
  5 + * ----------------------------------------------------------------------- *
  6 + * *
  7 + * This program is free software: you can redistribute it and/or modify *
  8 + * it under the terms of the GNU General Public License as published by *
  9 + * the Free Software Foundation, either version 3 of the License, or *
  10 + * (at your option) any later version. *
  11 + * *
  12 + * This program is distributed in the hope that it will be useful, *
  13 + * but WITHOUT ANY WARRANTY; without even the implied warranty of *
  14 + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
  15 + * GNU General Public License for more details. *
  16 + * *
  17 + * You should have received a copy of the GNU General Public License *
  18 + * along with this program. If not, see <http://www.gnu.org/licenses/>. *
  19 + * *
  20 + ***************************************************************************/
  21 +/**
  22 + * @file ctbutterfly/main.cpp
  23 + * @brief CTA butterfly computation tool
  24 + * @author Michael Mayer
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include <cstdio>
  32 +#include "ctbutterfly.hpp"
  33 +
  34 +
  35 +/***********************************************************************//**
  36 + * @brief Main entry point of application
  37 + *
  38 + * @param[in] argc Number of command line arguments.
  39 + * @param[in] argv Command line arguments.
  40 + *
  41 + * This is the main entry point of the ctbutterfly application. It allocates a
  42 + * ctbutterfly object and runs the application.
  43 + ***************************************************************************/
  44 +int main (int argc, char *argv[])
  45 +{
  46 + // Initialise return code
  47 + int rc = 1;
  48 +
  49 + // Create instance of application
  50 + ctbutterfly application(argc, argv);
  51 +
  52 + // Run application
  53 + try {
  54 + // Execute application
  55 + application.execute();
  56 +
  57 + // Signal success
  58 + rc = 0;
  59 + }
  60 + catch (std::exception &e) {
  61 +
  62 + // Extract error message
  63 + std::string message = e.what();
  64 + std::string signal = "*** ERROR encounterted in the execution of"
  65 + " ctbutterfly. Run aborted ...";
  66 +
  67 + // Write error in logger
  68 + application.log << signal << std::endl;
  69 + application.log << message << std::endl;
  70 +
  71 + // Write error on standard output
  72 + std::cout << signal << std::endl;
  73 + std::cout << message << std::endl;
  74 +
  75 + } // endcatch: catched any application error
  76 +
  77 + // Return
  78 + return rc;
  79 +}
... ...
src/ctobssim/ctobssim.cpp
  View file @22700ce
... ... @@ -963,6 +963,7 @@ void ctobssim::simulate_source(GCTAObservation* obs, const GModels&amp; models,
963 963 photons[i],
964 964 *obs,
965 965 ran);
  966 +
966 967 if (event != NULL) {
967 968  
968 969 // Use event only if it falls within ROI
... ...
test/test_ctbutterfly.py 0 → 100755
  View file @22700ce
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This scripts performs unit tests for the ctbutterfly tool.
  4 +#
  5 +# Copyright (C) 2014 Michal Mayer
  6 +#
  7 +# This program is free software: you can redistribute it and/or modify
  8 +# it under the terms of the GNU General Public License as published by
  9 +# the Free Software Foundation, either version 3 of the License, or
  10 +# (at your option) any later version.
  11 +#
  12 +# This program is distributed in the hope that it will be useful,
  13 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  15 +# GNU General Public License for more details.
  16 +#
  17 +# You should have received a copy of the GNU General Public License
  18 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  19 +#
  20 +# ==========================================================================
  21 +import gammalib
  22 +import ctools
  23 +
  24 +
  25 +# =============================== #
  26 +# Test class for ctbutterfly tool #
  27 +# =============================== #
  28 +class Test(gammalib.GPythonTestSuite):
  29 + """
  30 + Test class for ctbutterfly tool.
  31 + """
  32 + # Constructor
  33 + def __init__(self):
  34 + """
  35 + Constructor.
  36 + """
  37 + # Call base class constructor
  38 + gammalib.GPythonTestSuite.__init__(self)
  39 +
  40 + # Set members
  41 + self.events_name = "data/crab_events.fits"
  42 + self.model_name = "data/crab.xml"
  43 + self.caldb = "irf"
  44 + self.irf = "cta_dummy_irf"
  45 +
  46 + # Return
  47 + return
  48 +
  49 + # Set test functions
  50 + def set(self):
  51 + """
  52 + Set all test functions.
  53 + """
  54 + # Set test name
  55 + self.name("ctbutterfly")
  56 +
  57 + # Append tests
  58 + self.append(self.test_functional, "Test ctbutterfly functionality")
  59 +
  60 + # Return
  61 + return
  62 +
  63 + # Test ctbutterfly functionnality
  64 + def test_functional(self):
  65 + """
  66 + Test ctbutterfly functionnality.
  67 + """
  68 + # Set-up ctbutterfly
  69 + butterfly = ctools.ctbutterfly()
  70 + butterfly["infile"].filename(self.events_name)
  71 + butterfly["srcmdl"].filename(self.model_name)
  72 + butterfly["srcname"].string("Crab")
  73 + butterfly["caldb"].string(self.caldb)
  74 + butterfly["irf"].string(self.irf)
  75 + butterfly["emin"].real(0.1)
  76 + butterfly["emax"].real(100)
  77 + butterfly["outfile"].filename("butterfly.txt")
  78 +
  79 + # Run tool
  80 + self.test_try("Run ctbutterfly")
  81 + try:
  82 + butterfly.run()
  83 + self.test_try_success()
  84 + except:
  85 + self.test_try_failure("Exception occured in butterfly.")
  86 +
  87 + # Save results
  88 + self.test_try("Save results")
  89 + try:
  90 + butterfly.save()
  91 + self.test_try_success()
  92 + except:
  93 + self.test_try_failure("Exception occured in saving results.")
  94 +
  95 + # Return
  96 + return
... ...
test/test_ctools.sh
  View file @22700ce
... ... @@ -45,6 +45,7 @@ ctexpcube=../src/ctexpcube/ctexpcube
45 45 ctpsfcube=../src/ctpsfcube/ctpsfcube
46 46 ctbkgcube=../src/ctbkgcube/ctbkgcube
47 47 ctcubemask=../src/ctcubemask/ctcubemask
  48 +ctbutterfly=../src/ctbutterfly/ctbutterfly
48 49  
49 50  
50 51 #
... ... @@ -542,3 +543,26 @@ else
542 543 exit 1
543 544 fi
544 545 $ECHO " ok"
  546 +
  547 +
  548 +#
  549 +# Test ctbutterfly
  550 +# ===============
  551 +$ECHO -n "Test ctbutterfly: "
  552 +$ctbutterfly infile="data/crab_events.fits.gz" \
  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
  560 +$ECHO -n "."
  561 +if [ -s "butterfly.txt" ]
  562 +then
  563 + $ECHO -n "."
  564 +else
  565 + $ECHO " butterfly.txt file is not valid"
  566 + exit 1
  567 +fi
  568 +$ECHO " ok"
... ...
test/test_cttsmap.py
  View file @22700ce
... ... @@ -40,7 +40,7 @@ class Test(gammalib.GPythonTestSuite):
40 40 # Set members
41 41 self.events_name = "data/crab_events.fits"
42 42 self.model_name = "data/crab.xml"
43   - self.caldb = "irf"
  43 + self.caldb = "dummy"
44 44 self.irf = "cta_dummy_irf"
45 45  
46 46 # Return
... ...
test/test_python.py
  View file @22700ce
... ... @@ -34,6 +34,7 @@ import test_ctpsfcube
34 34 import test_ctbkgcube
35 35 import test_ctcubemask
36 36 import test_pipelines
  37 +import test_ctbutterfly
37 38  
38 39  
39 40 # ======================== #
... ... @@ -59,6 +60,7 @@ if __name__ == &#39;__main__&#39;:
59 60 suite_ctbkgcube = test_ctbkgcube.Test()
60 61 suite_ctcubemask = test_ctcubemask.Test()
61 62 suite_pipelines = test_pipelines.Test()
  63 + suite_ctbutterfly = test_ctbutterfly.Test()
62 64  
63 65 # Setup unit tests
64 66 suite_ctobssim.set()
... ... @@ -73,6 +75,7 @@ if __name__ == &#39;__main__&#39;:
73 75 suite_ctbkgcube.set()
74 76 suite_ctcubemask.set()
75 77 suite_pipelines.set()
  78 + suite_ctbutterfly.set()
76 79  
77 80 # Append tests to container
78 81 suites.append(suite_ctobssim)
... ... @@ -87,6 +90,7 @@ if __name__ == &#39;__main__&#39;:
87 90 suites.append(suite_ctbkgcube)
88 91 suites.append(suite_ctcubemask)
89 92 suites.append(suite_pipelines)
  93 + suites.append(suite_ctbutterfly)
90 94  
91 95 # Set PFILES environment variable
92 96 try:
... ...