Olga Sergijenko / ctools

Commit 77328f2638f3afdcbe6df09a58e327b44f1f2696

Authored by Jürgen Knödlseder
2 parents 4b7d6042 fc5614e8

Merge branch 'pull-request-1991-ctprob-eval-event-probs' into integration

Inline Side-by-side
Showing 19 changed files with 1697 additions and 150 deletions
.gitignore
  View file @77328f2
... ... @@ -71,6 +71,7 @@ src/ctcubemask/ctcubemask
71 71 src/ctbutterfly/ctbutterfly
72 72 src/ctulimit/ctulimit
73 73 src/cterror/cterror
  74 +src/ctprob/ctprob
74 75 test/pfiles/
75 76 test/caldb/
76 77 test/prod3
... ...
configure.ac
  View file @77328f2
... ... @@ -1035,6 +1035,7 @@ AC_CONFIG_FILES([Makefile
1035 1035 src/ctbutterfly/Makefile
1036 1036 src/ctulimit/Makefile
1037 1037 src/cterror/Makefile
  1038 + src/ctprob/Makefile
1038 1039 src/support/Makefile
1039 1040 src/ctools-setup
1040 1041 test/Makefile
... ...
doc/source/users/reference_manual/ctprob.rst 0 → 100644
  View file @77328f2
  1 +.. _ctprob:
  2 +
  3 +ctprob
  4 +========
  5 +
  6 +Compute event probability for a given model.
  7 +
  8 +
  9 +Synopsis
  10 +--------
  11 +
  12 +This tool takes an event list file (or an XML observation definition file) and
  13 +a model as input and computes for each event the probability that it originates
  14 +from a specific model component. The probability is computed by by evaluating
  15 +the differential event probability for a given model and normalizing these
  16 +probabilities so that the sum of the probabilities for all model components is
  17 +unity for each event.
  18 +
  19 +The tool appends for each model component a single precision column to the event
  20 +list that contains the event probabilities. The name of the column is build from
  21 +the model name prefixed with ``PROB_``.
  22 +
  23 +
  24 +General parameters
  25 +------------------
  26 +
  27 +``inobs [file]``
  28 + Input event list or observation definition XML file
  29 +
  30 +``outobs [file]``
  31 + Output event list or observation definition XML file.
  32 +
  33 +``(prefix = "prob_") [string]``
  34 + Prefix for output event lists in observation definition file.
  35 +
  36 +``inmodel [string]``
  37 + Input model definition XML file.
  38 +
  39 +``caldb [string]``
  40 + Calibration database.
  41 +
  42 +``irf [string]``
  43 + Instrument response function.
  44 +
  45 +``(edisp = no) [boolean]``
  46 + Applies energy dispersion to response computation.
  47 +
  48 +
  49 +Standard parameters
  50 +-------------------
  51 +
  52 +``(publish = no) [boolean]``
  53 + Specifies whether the event list(s) should be published on VO Hub.
  54 +
  55 +``(chatter = 2) [integer]``
  56 + Verbosity of the executable:
  57 + chatter = 0: no information will be logged
  58 +
  59 + chatter = 1: only errors will be logged
  60 +
  61 + chatter = 2: errors and actions will be logged
  62 +
  63 + chatter = 3: report about the task execution
  64 +
  65 + chatter = 4: detailed report about the task execution
  66 +
  67 +``(clobber = yes) [boolean]``
  68 + Specifies whether existing output files should be overwritten.
  69 +
  70 +``(debug = no) [boolean]``
  71 + Enables debug mode. In debug mode the executable will dump any log file output to the console.
  72 +
  73 +``(mode = ql) [string]``
  74 + Mode of automatic parameters (default is "ql", i.e. "query and learn").
  75 +
  76 +``(logfile = ctprob.log) [string]``
  77 + Name of log file.
  78 +
  79 +
  80 +Related tools or scripts
  81 +------------------------
  82 +
  83 +:doc:`ctphase`
  84 +:doc:`csphasecrv`
  85 +
... ...
doc/source/users/reference_manual/index.rst
  View file @77328f2
... ... @@ -27,6 +27,7 @@ ctools
27 27 ctmodel --- Computes model counts cube <ctmodel>
28 28 ctobssim --- Simulate observations <ctobssim>
29 29 ctpsfcube --- Generates point spread function cube <ctpsfcube>
  30 + ctprob --- Computes event probability for a given model <ctprob>
30 31 ctselect --- Selects event data <ctselect>
31 32 ctskymap --- Generates sky map <ctskymap>
32 33 cttsmap --- Generates Test Statistic map <cttsmap>
... ...
pyext/Makefile.am
  View file @77328f2
... ... @@ -58,6 +58,7 @@ testscripts = ctools/tests/__init__.py \
58 58 $(top_srcdir)/test/test_ctbutterfly.py \
59 59 $(top_srcdir)/test/test_ctulimit.py \
60 60 $(top_srcdir)/test/test_cterror.py \
  61 + $(top_srcdir)/test/test_ctprob.py \
61 62 $(top_srcdir)/test/test_pipelines.py
62 63  
63 64 # Test data to be distributed in $(pyexecdir)/ctools/tests/data
... ...
pyext/ctools/tools.i
  View file @77328f2
... ... @@ -72,6 +72,7 @@ namespace std {
72 72 %include "ctbutterfly.i"
73 73 %include "ctulimit.i"
74 74 %include "cterror.i"
  75 +%include "ctprob.i"
75 76  
76 77 /* __ Test function ______________________________________________________ */
77 78 %pythoncode %{
... ...
pyext/ctprob.i 0 → 100644
  View file @77328f2
  1 +/***************************************************************************
  2 + * ctprob - Computes probability for a given model *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2012-2016 by Leonardo Di Venere *
  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 ctprob.i
  23 + * @brief Computes probability for a given model
  24 + * @author Leonardo Di Venere
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "ctprob.hpp"
  29 +%}
  30 +
  31 +
  32 +/***********************************************************************//**
  33 + * @class ctprob
  34 + *
  35 + * @brief Data selection tool
  36 + ***************************************************************************/
  37 +class ctprob : public ctobservation {
  38 +public:
  39 + // Constructors and destructors
  40 + ctprob(void);
  41 + explicit ctprob(const GObservations& obs);
  42 + ctprob(int argc, char *argv[]);
  43 + ctprob(const ctprob& app);
  44 + virtual ~ctprob(void);
  45 +
  46 + // Methods
  47 + void clear(void);
  48 + void run(void);
  49 + void save(void);
  50 + void publish(const std::string& name = "");
  51 +};
  52 +
  53 +
  54 +/***********************************************************************//**
  55 + * @brief Data selection tool Python extension
  56 + ***************************************************************************/
  57 +%extend ctprob {
  58 + ctprob copy() {
  59 + return (*self);
  60 + }
  61 +}
... ...
src/Makefile.am
  View file @77328f2
... ... @@ -40,7 +40,8 @@ SUBDIRS = support \
40 40 cttsmap \
41 41 ctbutterfly \
42 42 ctulimit \
43   - cterror
  43 + cterror \
  44 + ctprob
44 45  
45 46 # Scripts to be distributed
46 47 scripts = ctools-setup \
... ... @@ -68,7 +69,8 @@ libctools_la_LIBADD = support/libsupport.la \
68 69 cttsmap/libcttsmap.la \
69 70 ctbutterfly/libctbutterfly.la \
70 71 ctulimit/libctulimit.la \
71   - cterror/libcterror.la
  72 + cterror/libcterror.la \
  73 + ctprob/libctprob.la
72 74 libctools_la_LDFLAGS = -version-info $(CTOOLS_LT_VERSION) @OPENMP_LIB_LDFLAGS@
73 75  
74 76 # Distribute scripts
... ...
src/ctprob/Makefile.am 0 → 100644
  View file @77328f2
  1 +#############################################################################
  2 +# Makefile.am for ctprob tool #
  3 +# ------------------------------------------------------------------------- #
  4 +# Copyright (C) 2017 Leonardo Di Venere #
  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 +# Set include directories
  23 +INCLUDES = -I$(top_srcdir)/src/support
  24 +
  25 +# Program name
  26 +bin_PROGRAMS = ctprob
  27 +
  28 +# Sources
  29 +sources = ctprob.cpp ctprob.hpp
  30 +
  31 +# Libtool library
  32 +noinst_LTLIBRARIES = libctprob.la
  33 +libctprob_la_SOURCES = $(sources)
  34 +libctprob_la_LDFLAGS = @OPENMP_LIB_LDFLAGS@
  35 +
  36 +# Program definition
  37 +ctprob_SOURCES = main.cpp
  38 +ctprob_LDFLAGS = $(LIB_LINK_OPTION) @LDFLAGS@
  39 +ctprob_CPPFLAGS = @CPPFLAGS@
  40 +ctprob_LDADD = libctprob.la $(top_builddir)/src/support/libsupport.la
  41 +
  42 +# Parameter files
  43 +pfilesdir = $(prefix)/syspfiles
  44 +pfiles_DATA = $(srcdir)/ctprob.par
  45 +
  46 +# Extra files to be distributed in the tar file
  47 +EXTRA_DIST = $(pfiles_DATA)
... ...
src/ctprob/ctprob.cpp 0 → 100644
  View file @77328f2
  1 +/***************************************************************************
  2 + * ctprob - Computes event probability for a given model *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2017 by Leonardo Di Venere *
  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 ctprob.cpp
  23 + * @brief Computes event probability for a given model
  24 + * @author Leonardo Di Venere
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include <cstdlib>
  32 +#include "ctprob.hpp"
  33 +#include "GTools.hpp"
  34 +
  35 +/* __ Method name definitions ____________________________________________ */
  36 +#define G_GET_OBS "ctprob::get_obs()"
  37 +#define G_EVALUATE_PROB "ctprob::evaluate_probability(GCTAObservation*)"
  38 +
  39 +/* __ Debug definitions __________________________________________________ */
  40 +
  41 +/* __ Coding definitions _________________________________________________ */
  42 +
  43 +
  44 +/*==========================================================================
  45 + = =
  46 + = Constructors/destructors =
  47 + = =
  48 + ==========================================================================*/
  49 +
  50 +/***********************************************************************//**
  51 + * @brief Void constructor
  52 + ***************************************************************************/
  53 +ctprob::ctprob(void) : ctobservation(CTPROB_NAME, VERSION)
  54 +{
  55 + // Initialise members
  56 + init_members();
  57 +
  58 + // Return
  59 + return;
  60 +}
  61 +
  62 +
  63 +/***********************************************************************//**
  64 + * @brief Observations constructor
  65 + *
  66 + * param[in] obs Observation container.
  67 + *
  68 + * Creates an instance of the class that is initialised using the information
  69 + * provided in an observation container.
  70 + ***************************************************************************/
  71 +ctprob::ctprob(const GObservations& obs) :
  72 + ctobservation(CTPROB_NAME, VERSION, obs)
  73 +{
  74 + // Initialise members
  75 + init_members();
  76 +
  77 + // Return
  78 + return;
  79 +}
  80 +
  81 +
  82 +
  83 +
  84 +/***********************************************************************//**
  85 + * @brief Command line constructor
  86 + *
  87 + * @param[in] argc Number of arguments in command line.
  88 + * @param[in] argv Array of command line arguments.
  89 + ***************************************************************************/
  90 +ctprob::ctprob(int argc, char *argv[]) :
  91 + ctobservation(CTPROB_NAME, VERSION, argc, argv)
  92 +{
  93 + // Initialise members
  94 + init_members();
  95 +
  96 + // Return
  97 + return;
  98 +}
  99 +
  100 +
  101 +/***********************************************************************//**
  102 + * @brief Copy constructor
  103 + *
  104 + * @param[in] app Application.
  105 + ***************************************************************************/
  106 +ctprob::ctprob(const ctprob& app) : ctobservation(app)
  107 +{
  108 + // Initialise members
  109 + init_members();
  110 +
  111 + // Copy members
  112 + copy_members(app);
  113 +
  114 + // Return
  115 + return;
  116 +}
  117 +
  118 +
  119 +/***********************************************************************//**
  120 + * @brief Destructor
  121 + ***************************************************************************/
  122 +ctprob::~ctprob(void)
  123 +{
  124 + // Free members
  125 + free_members();
  126 +
  127 + // Return
  128 + return;
  129 +}
  130 +
  131 +
  132 +/*==========================================================================
  133 + = =
  134 + = Operators =
  135 + = =
  136 + ==========================================================================*/
  137 +
  138 +/***********************************************************************//**
  139 + * @brief Assignment operator
  140 + *
  141 + * @param[in] app Application.
  142 + * @return Application.
  143 + ***************************************************************************/
  144 +ctprob& ctprob::operator=(const ctprob& app)
  145 +{
  146 + // Execute only if object is not identical
  147 + if (this != &app) {
  148 +
  149 + // Copy base class members
  150 + this->ctobservation::operator=(app);
  151 +
  152 + // Free members
  153 + free_members();
  154 +
  155 + // Initialise members
  156 + init_members();
  157 +
  158 + // Copy members
  159 + copy_members(app);
  160 +
  161 + } // endif: object was not identical
  162 +
  163 + // Return this object
  164 + return *this;
  165 +}
  166 +
  167 +
  168 +/*==========================================================================
  169 + = =
  170 + = Public methods =
  171 + = =
  172 + ==========================================================================*/
  173 +
  174 +/***********************************************************************//**
  175 + * @brief Clear ctprob tool
  176 + *
  177 + * Clears ctprob tool.
  178 + ***************************************************************************/
  179 +void ctprob::clear(void)
  180 +{
  181 + // Free members
  182 + free_members();
  183 + this->ctobservation::free_members();
  184 + this->ctool::free_members();
  185 +
  186 + // Clear base class (needed to conserve tool name and version)
  187 + this->GApplication::clear();
  188 +
  189 + // Initialise members
  190 + this->ctool::init_members();
  191 + this->ctobservation::init_members();
  192 + init_members();
  193 +
  194 + // Write header into logger
  195 + log_header();
  196 +
  197 + // Return
  198 + return;
  199 +}
  200 +
  201 +
  202 +/***********************************************************************//**
  203 + * @brief Runs the ctprob tool
  204 + *
  205 + * This method reads in the application parameters and loops over all
  206 + * unbinned observations to compute the probability for each event that it
  207 + * comes from one of the model components.
  208 + ***************************************************************************/
  209 +void ctprob::run(void)
  210 +{
  211 + // Switch screen logging on in debug mode
  212 + if (logDebug()) {
  213 + log.cout(true);
  214 + }
  215 +
  216 + // Get parameters
  217 + get_parameters();
  218 +
  219 + // Set energy dispersion flags of all CTA observations and save old
  220 + // values in save_edisp vector
  221 + std::vector<bool> save_edisp = set_edisp(m_obs, m_apply_edisp);
  222 +
  223 + // Write input observation container into logger
  224 + log_observations(NORMAL, m_obs, "Input observation");
  225 +
  226 + // Write header into logger
  227 + log_header1(TERSE, "Compute event probabilities");
  228 +
  229 + // Initialise counters
  230 + int n_observations = 0;
  231 +
  232 + // Loop over all unbinned CTA observations in the container
  233 + for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL;
  234 + obs = next_unbinned_observation()) {
  235 +
  236 + // Increment counter
  237 + n_observations++;
  238 +
  239 + // Evaluate event probabilities
  240 + evaluate_probability(obs);
  241 +
  242 + } // endfor: looped over observations
  243 +
  244 + // If more than a single observation has been handled then make sure that
  245 + // an XML file will be used for storage
  246 + if (n_observations > 1) {
  247 + m_use_xml = true;
  248 + }
  249 +
  250 + // Write observation(s) into logger
  251 + log_observations(NORMAL, m_obs, "Output observation");
  252 +
  253 + // Restore energy dispersion flags of all CTA observations
  254 + restore_edisp(m_obs, save_edisp);
  255 +
  256 + // Optionally publish event list(s)
  257 + if ((*this)["publish"].boolean()) {
  258 + publish();
  259 + }
  260 +
  261 + // Return
  262 + return;
  263 +}
  264 +
  265 +
  266 +/***********************************************************************//**
  267 + * @brief Save the selected event list(s)
  268 + *
  269 + * This method saves the selected event list(s) into FITS file(s). There are
  270 + * two modes, depending on the m_use_xml flag.
  271 + *
  272 + * If m_use_xml is true, all selected event list(s) will be saved into FITS
  273 + * files, where the output filenames are constructued from the input
  274 + * filenames by prepending the m_prefix string to name. Any path information
  275 + * will be stripped form the input name, hence event files will be written
  276 + * into the local working directory (unless some path information is present
  277 + * in the prefix). In addition, an XML file will be created that gathers
  278 + * the filename information for the selected event list(s). If an XML file
  279 + * was present on input, all metadata information will be copied from this
  280 + * input file.
  281 + *
  282 + * If m_use_xml is false, the selected event list will be saved into a FITS
  283 + * file.
  284 + ***************************************************************************/
  285 +void ctprob::save(void)
  286 +{
  287 + // Write header into logger
  288 + log_header1(TERSE, gammalib::number("Save event list", m_obs.size()));
  289 +
  290 + // Case A: Save event file(s) and XML metadata information
  291 + if (m_use_xml) {
  292 + save_xml();
  293 + }
  294 +
  295 + // Case B: Save event file as FITS file
  296 + else {
  297 + save_fits();
  298 + }
  299 +
  300 + // Return
  301 + return;
  302 +}
  303 +
  304 +
  305 +/***********************************************************************//**
  306 + * @brief Publish event lists
  307 + *
  308 + * @param[in] name Event list name.
  309 + ***************************************************************************/
  310 +void ctprob::publish(const std::string& name)
  311 +{
  312 + // Write header into logger
  313 + log_header1(TERSE, gammalib::number("Publish event list", m_obs.size()));
  314 +
  315 + // Loop over all observation in the container
  316 + for (int i = 0; i < m_obs.size(); ++i) {
  317 +
  318 + // Get CTA observation
  319 + GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
  320 +
  321 + // Handle only CTA observations
  322 + if (obs != NULL) {
  323 +
  324 + // Continue only if there is an event list
  325 + if (obs->events()->size() != 0) {
  326 +
  327 + // Set default name if user name is empty
  328 + std::string user_name(name);
  329 + if (user_name.empty()) {
  330 + user_name = CTPROB_NAME;
  331 + }
  332 +
  333 + // If there are several event lists then add an index
  334 + if (m_use_xml) {
  335 + user_name += gammalib::str(i);
  336 + }
  337 +
  338 + // Write event list name into logger
  339 + log_value(NORMAL, "Event list name", user_name);
  340 +
  341 + // Write events into in-memory FITS file
  342 + GFits fits;
  343 + obs->write(fits);
  344 +
  345 + // Publish
  346 + fits.publish(gammalib::extname_cta_events, user_name);
  347 +
  348 + } // endif: there were events
  349 +
  350 + } // endif: observation was a CTA observation
  351 +
  352 + } // endfor: looped over observations
  353 +
  354 + // Return
  355 + return;
  356 +}
  357 +
  358 +
  359 +/*==========================================================================
  360 + = =
  361 + = Private methods =
  362 + = =
  363 + ==========================================================================*/
  364 +
  365 +/***********************************************************************//**
  366 + * @brief Initialise class members
  367 + ***************************************************************************/
  368 +void ctprob::init_members(void)
  369 +{
  370 + // Initialise parameters
  371 + m_outobs.clear();
  372 + m_prefix.clear();
  373 + /* m_usepnt = false;
  374 + m_roi.clear();
  375 + m_tmin = 0.0;
  376 + m_tmax = 0.0;
  377 + m_emin = 0.0;
  378 + m_emax = 0.0;
  379 + m_expr.clear();
  380 + m_usethres.clear();*/
  381 + m_chatter = static_cast<GChatter>(2);
  382 +
  383 + // Initialise protected members
  384 + m_infiles.clear();
  385 + m_evtname.clear();
  386 + m_gtiname.clear();
  387 + // m_timemin.clear();
  388 + //m_timemax.clear();
  389 + //m_select_energy = true;
  390 + //m_select_roi = true;
  391 + //m_select_time = true;
  392 +
  393 + // Return
  394 + return;
  395 +}
  396 +
  397 +
  398 +/***********************************************************************//**
  399 + * @brief Copy class members
  400 + *
  401 + * @param[in] app Application.
  402 + ***************************************************************************/
  403 +void ctprob::copy_members(const ctprob& app)
  404 +{
  405 + // Copy parameters
  406 + m_outobs = app.m_outobs;
  407 + m_prefix = app.m_prefix;
  408 + /* m_usepnt = app.m_usepnt;
  409 + m_roi = app.m_roi;
  410 + m_tmin = app.m_tmin;
  411 + m_tmax = app.m_tmax;
  412 + m_emin = app.m_emin;
  413 + m_emax = app.m_emax;
  414 + m_expr = app.m_expr;
  415 + m_usethres = app.m_usethres;*/
  416 + m_chatter = app.m_chatter;
  417 +
  418 + // Copy protected members
  419 + m_infiles = app.m_infiles;
  420 + m_evtname = app.m_evtname;
  421 + m_gtiname = app.m_gtiname;
  422 + /* m_timemin = app.m_timemin;
  423 + m_timemax = app.m_timemax;
  424 + m_select_energy = app.m_select_energy;
  425 + m_select_roi = app.m_select_roi;
  426 + m_select_time = app.m_select_time;*/
  427 +
  428 + // Return
  429 + return;
  430 +}
  431 +
  432 +
  433 +/***********************************************************************//**
  434 + * @brief Delete class members
  435 + ***************************************************************************/
  436 +void ctprob::free_members(void)
  437 +{
  438 + // Return
  439 + return;
  440 +}
  441 +
  442 +
  443 +/***********************************************************************//**
  444 + * @brief Get application parameters
  445 + *
  446 + * Get all user parameters from parameter file or (if required) by querying
  447 + * the user. Times are assumed to be in the native CTA MJD format.
  448 + *
  449 + * This method also loads observations if no observations are yet allocated.
  450 + * Observations are either loaded from a single CTA even list, or from a
  451 + * XML file using the metadata information that is stored in that file.
  452 + ***************************************************************************/
  453 +void ctprob::get_parameters(void)
  454 +{
  455 + // Setup observations from "inobs" parameter. Do not request response
  456 + // information and do not accept counts cubes.
  457 + if (m_obs.size() == 0) {
  458 + get_obs();
  459 + }
  460 + // ... otherwise add response information and energy boundaries in case
  461 + // that they are missing
  462 + else {
  463 + setup_observations(m_obs);
  464 + }
  465 +
  466 + // Read model definition file if required
  467 + if (m_obs.models().size() == 0) {
  468 +
  469 + // Get model filename
  470 + std::string inmodel = (*this)["inmodel"].filename();
  471 +
  472 + // Load models from file
  473 + m_obs.models(inmodel);
  474 +
  475 + } // endif: there were no models
  476 +
  477 + // Get energy dispersion flag parameters
  478 + m_apply_edisp = (*this)["edisp"].boolean();
  479 +
  480 +
  481 + // Get remaining parameters
  482 + m_publish = (*this)["publish"].boolean();
  483 + m_chatter = static_cast<GChatter>((*this)["chatter"].integer());
  484 +
  485 + // Optionally read ahead parameters so that they get correctly
  486 + // dumped into the log file
  487 + if (read_ahead()) {
  488 + m_outobs = (*this)["outobs"].filename();
  489 + m_prefix = (*this)["prefix"].string();
  490 + }
  491 +
  492 + // Write parameters into logger
  493 + log_parameters(TERSE);
  494 +
  495 + // Return
  496 + return;
  497 +}
  498 +
  499 +
  500 +/***********************************************************************//**
  501 + * @brief Get observation container
  502 + *
  503 + * @exception GException::invalid_value
  504 + * Invalid FITS file.
  505 + *
  506 + * Get an observation container according to the user parameters. The method
  507 + * supports loading of a individual FITS file or an observation definition
  508 + * file in XML format.
  509 + ***************************************************************************/
  510 +void ctprob::get_obs(void)
  511 +{
  512 + // Get the filename from the input parameters
  513 + std::string filename = (*this)["inobs"].filename();
  514 +
  515 + // If no observation definition file throw an exception
  516 + if (!is_valid_filename(filename)) {
  517 + throw GException::invalid_value(G_GET_OBS, "Input event list is not valid.");
  518 + } // endif: filename was not valid
  519 +
  520 + // ... otherwise we have a file name
  521 + else {
  522 +
  523 + // If file is a FITS file then create an empty CTA observation
  524 + // and load file into observation
  525 + if (GFilename(filename).is_fits()) {
  526 +
  527 + // Allocate empty CTA observation
  528 + GCTAObservation cta;
  529 +
  530 + // Load data
  531 + cta.load(filename);
  532 +
  533 + // Set response
  534 + set_obs_response(&cta);
  535 +
  536 + // Append observation to container
  537 + m_obs.append(cta);
  538 +
  539 + // Signal that no XML file should be used for storage
  540 + m_use_xml = false;
  541 +
  542 + }
  543 +
  544 + // ... otherwise load file into observation container
  545 + else {
  546 +
  547 + // Load observations from XML file
  548 + m_obs.load(filename);
  549 +
  550 + // For all observations that have no response, set the response
  551 + // from the task parameters
  552 + set_response(m_obs);
  553 +
  554 + // Set observation boundary parameters (emin, emax, rad)
  555 + set_obs_bounds(m_obs);
  556 +
  557 + // Signal that XML file should be used for storage
  558 + m_use_xml = true;
  559 +
  560 + } // endelse: file was an XML file
  561 +
  562 + }
  563 +
  564 + // Return
  565 + return;
  566 +
  567 +}
  568 +
  569 +
  570 +
  571 +/***********************************************************************//**
  572 + * @brief Evaluate probability for events
  573 + *
  574 + * @param[in,out] obs CTA observation.
  575 + *
  576 + * @exception GException::invalid_value
  577 + * Observation does not contain an event list.
  578 + *
  579 + * Computes for all events the probability that an event arises from a
  580 + * specific model component. The computation is done by by evaluating the
  581 + * differential event probability for a given model and normalizing these
  582 + * probabilities so that the sum of the probabilities for all model
  583 + * components is unity for each event.
  584 + *
  585 + * The method appends for each model component a single precision column to
  586 + * the event list that contains the event probabilities. The name of the
  587 + * column is build from the model name prefixed with `PROB_`.
  588 + ***************************************************************************/
  589 +void ctprob::evaluate_probability(GCTAObservation* obs)
  590 +{
  591 + // Write header into logger
  592 + log_header3(NORMAL, "Evaluate events probabilities");
  593 +
  594 + // Make sure that the observation holds a CTA event list. If this
  595 + // is not the case then throw an exception.
  596 + GCTAEventList* events = dynamic_cast<GCTAEventList*>(const_cast<GEvents*>
  597 + (obs->events()));
  598 + if (events == NULL) {
  599 + std::string msg = "CTA Observation does not contain an event list. An "
  600 + "event list is needed to compute the event "
  601 + "probabilities.";
  602 + throw GException::invalid_value(G_EVALUATE_PROB, msg);
  603 + }
  604 +
  605 + // Continue only if there are model components
  606 + int nmodels = m_obs.models().size();
  607 + if (nmodels > 0) {
  608 +
  609 + // Allocate single precision FITS table columns for all model components
  610 + std::vector<GFitsTableFloatCol*> columns;
  611 + for (int k = 0; k < nmodels; ++k) {
  612 + const std::string colname = "PROB_" + m_obs.models()[k] ->name();
  613 + GFitsTableFloatCol* col =
  614 + new GFitsTableFloatCol(colname, events->size());
  615 + columns.push_back(col);
  616 + }
  617 +
  618 + // Loop over all events in the event list
  619 + for (int i = 0; i < events->size(); ++i) {
  620 +
  621 + // Initialise total probability
  622 + double total = 0.0;
  623 +
  624 + // Get point to event
  625 + const GEvent* event = (*events)[i];
  626 +
  627 + // Initialise probability vector
  628 + std::vector<double> probabilities(nmodels, 0.0);
  629 +
  630 + // Loop over models
  631 + for (int k = 0; k < nmodels; ++k) {
  632 + double value = m_obs.models()[k]->eval(*event, *obs);
  633 + probabilities[k] = value;
  634 + total += value;
  635 + }
  636 +
  637 + // Normalize probabilities to unity
  638 + if (total > 0.0) {
  639 + for (int k = 0; k < nmodels; ++k) {
  640 + probabilities[k] /= total;
  641 + }
  642 + }
  643 +
  644 + // Store probabilities in FITS columns
  645 + for (int k = 0; k < nmodels; ++k) {
  646 + (*(columns[k]))(i) = probabilities[k];
  647 + }
  648 +
  649 + } // endfor: looped over all events
  650 +
  651 + // Append columns to event list
  652 + for (int k = 0; k < nmodels; ++k) {
  653 + events->append_column(*columns[k]);
  654 + }
  655 +
  656 + } // endif: there were model components
  657 +
  658 + // Return
  659 + return;
  660 +}
  661 +
  662 +
  663 +/***********************************************************************//**
  664 + * @brief Set output file name.
  665 + *
  666 + * @param[in] filename Input file name.
  667 + *
  668 + * Converts an input file name into an output filename by prepending the
  669 + * prefix stored in the member m_prefix to the input file name. Any path as
  670 + * well as extension will be stripped from the input file name. Also a
  671 + * trailing ".gz" will be stripped as one cannot write into gzipped files.
  672 + ***************************************************************************/
  673 +std::string ctprob::set_outfile_name(const std::string& filename) const
  674 +{
  675 + // Create filename
  676 + GFilename fname(filename);
  677 +
  678 + // Split input filename without any extensions into path elements
  679 + std::vector<std::string> elements = gammalib::split(fname.url(), "/");
  680 +
  681 + // The last path element is the filename
  682 + std::string outname = m_prefix + elements[elements.size()-1];
  683 +
  684 + // Strip any ".gz"
  685 + outname = gammalib::strip_chars(outname, ".gz");
  686 +
  687 + // Return output filename
  688 + return outname;
  689 +}
  690 +
  691 +
  692 +/***********************************************************************//**
  693 + * @brief Save event list in FITS format.
  694 + *
  695 + * Save the event list as a FITS file. The file name of the FITS file is
  696 + * specified by the `outobs` parameter.
  697 + ***************************************************************************/
  698 +void ctprob::save_fits(void)
  699 +{
  700 + // Save only if there are observations
  701 + if (m_obs.size() > 0) {
  702 +
  703 + // Get output filename
  704 + m_outobs = (*this)["outobs"].filename();
  705 +
  706 + // Loop over all unbinned CTA observations in the container
  707 + for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL;
  708 + obs = next_unbinned_observation()) {
  709 +
  710 + // Get input file name and default extension names
  711 + std::string infile = obs->eventfile();
  712 + std::string evtname = gammalib::extname_cta_events;
  713 + std::string gtiname = gammalib::extname_gti;
  714 +
  715 + // Create file name object
  716 + GFilename fname(m_outobs);
  717 +
  718 + // Extract filename and event extension name
  719 + std::string outfile = fname.url();
  720 +
  721 + // Append event extension name. We handle here the possibility
  722 + // to write the events into a different extension.
  723 + if (fname.has_extname()) {
  724 + outfile += "["+fname.extname()+"]";
  725 + }
  726 + else {
  727 + outfile += "["+gammalib::extname_cta_events+"]";
  728 + }
  729 +
  730 + // Log filename
  731 + log_value(NORMAL, "Event list file", outfile);
  732 +
  733 + // Save event list
  734 + save_event_list(obs, infile, evtname, gtiname, outfile);
  735 +
  736 + // Exit the loop (if this method is called we should only
  737 + // have a single unbinned observation in the observation
  738 + // container
  739 + break;
  740 +
  741 + } // endfor: looped over all unbinned observations
  742 +
  743 + } // endif: there were observations
  744 +
  745 + // Return
  746 + return;
  747 +}
  748 +
  749 +
  750 +/***********************************************************************//**
  751 + * @brief Save event list(s) in XML format.
  752 + *
  753 + * Save the event list(s) into FITS files and write the file path information
  754 + * into a XML file. The filename of the XML file is specified by the `outobs`
  755 + * parameter, the filename(s) of the event lists are built by prepending a
  756 + * prefix to the input event list filenames. Any path present in the input
  757 + * filename will be stripped, i.e. the event list(s) will be written in the
  758 + * local working directory (unless a path is specified in the prefix).
  759 + ***************************************************************************/
  760 +void ctprob::save_xml(void)
  761 +{
  762 + // Get output filename and prefix
  763 + m_outobs = (*this)["outobs"].filename();
  764 + m_prefix = (*this)["prefix"].string();
  765 +
  766 + // Issue warning if output filename has no .xml suffix
  767 + log_string(TERSE, warn_xml_suffix(m_outobs));
  768 +
  769 + // Loop over all unbinned CTA observations in the container
  770 + for (GCTAObservation* obs = first_unbinned_observation(); obs != NULL;
  771 + obs = next_unbinned_observation()) {
  772 +
  773 + // Get input file name and default extension names
  774 + std::string infile = obs->eventfile();
  775 + std::string evtname = gammalib::extname_cta_events;
  776 + std::string gtiname = gammalib::extname_gti;
  777 +
  778 + // Extract event and GTI extension names from input FITS file
  779 + GFilename fname(infile);
  780 + if (fname.has_extname()) {
  781 + evtname = fname.extname();
  782 + }
  783 + gtiname = get_gtiname(fname.url(), evtname);
  784 +
  785 + // Set event output file name
  786 + std::string outfile = set_outfile_name(infile);
  787 +
  788 + // Append event extension name
  789 + outfile += "["+evtname+"]";
  790 +
  791 + // Log filename
  792 + log_value(NORMAL, "Event list file", outfile);
  793 +
  794 + // Store output file name in observation
  795 + obs->eventfile(outfile);
  796 +
  797 + // Save event list
  798 + save_event_list(obs, infile, evtname, gtiname, outfile);
  799 +
  800 + } // endfor: looped over observations
  801 +
  802 + // Write observation definition XML file name into logger
  803 + log_value(NORMAL, "Obs. definition file", m_outobs);
  804 +
  805 + // Save observations in XML file
  806 + m_obs.save(m_outobs);
  807 +
  808 + // Return
  809 + return;
  810 +}
... ...
src/ctprob/ctprob.hpp 0 → 100644
  View file @77328f2
  1 +/***************************************************************************
  2 + * ctprob - Computes event probability for a given model *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2017 by Leonardo Di Venere *
  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 ctprob.hpp
  23 + * @brief Computes event probability for a given model
  24 + * @author Leonardo Di Venere
  25 + */
  26 +
  27 +#ifndef CTPROB_HPP
  28 +#define CTPROB_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include <vector>
  32 +#include <string>
  33 +#include "GammaLib.hpp"
  34 +#include "GCTALib.hpp"
  35 +#include "ctobservation.hpp"
  36 +
  37 +/* __Definitions _________________________________________________________ */
  38 +#define CTPROB_NAME "ctprob"
  39 +
  40 +
  41 +/***********************************************************************//**
  42 + * @class ctprob
  43 + *
  44 + * @brief Data selection tool
  45 + ***************************************************************************/
  46 +class ctprob : public ctobservation {
  47 +
  48 +public:
  49 + // Constructors and destructors
  50 + ctprob(void);
  51 + explicit ctprob(const GObservations& obs);
  52 + ctprob(int argc, char *argv[]);
  53 + ctprob(const ctprob& app);
  54 + virtual ~ctprob(void);
  55 +
  56 + // Operators
  57 + ctprob& operator=(const ctprob& app);
  58 +
  59 + // Methods
  60 + void clear(void);
  61 + void run(void);
  62 + void save(void);
  63 + void publish(const std::string& name = "");
  64 +
  65 +protected:
  66 + // Protected methods
  67 + void init_members(void);
  68 + void copy_members(const ctprob& app);
  69 + void free_members(void);
  70 + void get_parameters(void);
  71 + void get_obs(void);
  72 + void evaluate_probability(GCTAObservation* obs);
  73 + std::string set_outfile_name(const std::string& filename) const;
  74 + void save_fits(void);
  75 + void save_xml(void);
  76 +
  77 + // User parameters
  78 + std::string m_outobs; //!< Output event list or XML file
  79 + std::string m_prefix; //!< Prefix for multiple event lists
  80 + bool m_apply_edisp; //!< Apply energy dispersion?
  81 + bool m_publish; //!< Publish event list?
  82 + GChatter m_chatter; //!< Chattiness
  83 +
  84 + // Protected members
  85 + std::vector<std::string> m_infiles; //!< Input event filenames
  86 + std::vector<std::string> m_evtname; //!< Event extension names
  87 + std::vector<std::string> m_gtiname; //!< GTI extension names
  88 +};
  89 +
  90 +#endif /* CTPROB_HPP */
... ...
src/ctprob/ctprob.par 0 → 100644
  View file @77328f2
  1 +#############################################################################
  2 +# #
  3 +# CTA Science Analysis Tools #
  4 +# #
  5 +# ------------------------------------------------------------------------- #
  6 +# #
  7 +# File: ctprob.par #
  8 +# #
  9 +# Author: Leonardo Di Venere #
  10 +# leonardo.divenere@ba.infn.it #
  11 +# INFN #
  12 +# #
  13 +# Purpose: Parameter file for the ctprob application #
  14 +# #
  15 +#############################################################################
  16 +#
  17 +# General parameters
  18 +#===================
  19 +inobs, f, a, "events.fits",,, "Input event list or observation definition XML file"
  20 +outobs, f, a, "prob_events.fits",,, "Output event list or observation definition XML file"
  21 +prefix, s, h, "prob_",,, "Prefix for event lists in observation definition XML file"
  22 +inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Input model definition XML file"
  23 +caldb, s, a, "prod2",,, "Calibration database"
  24 +irf, s, a, "South_0.5h",,, "Instrument response function"
  25 +edisp, b, h, no,,, "Apply energy dispersion?"
  26 +
  27 +#
  28 +# Standard parameters
  29 +#====================
  30 +publish, b, h, no,,, "Publish event list on VO Hub?"
  31 +chatter, i, h, 2,0,4, "Chattiness of output"
  32 +clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
  33 +debug, b, h, no,,, "Debugging mode activated"
  34 +mode, s, h, "ql",,, "Mode of automatic parameters"
  35 +logfile, f, h, "ctprob.log",,, "Log filename"
... ...
src/ctprob/main.cpp 0 → 100644
  View file @77328f2
  1 +/***************************************************************************
  2 + * ctprob - Computes event probability for a given model *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2017 by Leonardo Di Venere *
  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 main.cpp
  23 + * @brief Computes event probability for a given model
  24 + * @author Leonardo Di Venere
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include "support.hpp"
  32 +#include "ctprob.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 ctprob application. It allocates a
  42 + * ctprob object and executes the application. Any exceptions that occur
  43 + * will be catched and corresponding error messages written in the
  44 + * application logger and into the standard output.
  45 + ***************************************************************************/
  46 +int main (int argc, char *argv[])
  47 +{
  48 + // Initialise return code
  49 + int rc = 1;
  50 +
  51 + // Initialise pointer on application
  52 + ctprob* application = NULL;
  53 +
  54 + // Execute application
  55 + try {
  56 +
  57 + // Create instance of application
  58 + application = new ctprob(argc, argv);
  59 +
  60 + // Execute ctool
  61 + rc = execute_ctool(application);
  62 +
  63 + }
  64 +
  65 + catch (std::exception &e) {
  66 +
  67 + // Report exception
  68 + report_ctool_failure("ctprob", e.what());
  69 +
  70 + }
  71 +
  72 + // Delete application
  73 + delete application;
  74 +
  75 + // Return
  76 + return rc;
  77 +}
... ...
src/ctselect/ctselect.cpp
  View file @77328f2
... ... @@ -1175,39 +1175,6 @@ std::string ctselect::set_outfile_name(const std::string&amp; filename) const
1175 1175  
1176 1176  
1177 1177 /***********************************************************************//**
1178   - * @brief Get Good Time Intervals extension name
1179   - *
1180   - * @param[in] filename Input file name.
1181   - * @param[in] evtname Events extension name.
1182   - *
1183   - * Extracts the Good Time Intervals extension name from the event file. We
1184   - * do this by loading the events and accessing the Good Time Intervals
1185   - * extension name using the GCTAEventList::gtiname() method. If the file name
1186   - * is empty, the method returns `GTI`.
1187   - ***************************************************************************/
1188   -std::string ctselect::get_gtiname(const std::string& filename,
1189   - const std::string& evtname) const
1190   -{
1191   - // Initialise GTI name
1192   - std::string gtiname = gammalib::extname_gti;
1193   -
1194   - // Continue only if the filename is not empty
1195   - if (!filename.empty()) {
1196   -
1197   - // Load events
1198   - GCTAEventList events(filename+"["+evtname+"]");
1199   -
1200   - // Get GTI name
1201   - gtiname = events.gtiname();
1202   -
1203   - }
1204   -
1205   - // Return GTI name
1206   - return (gtiname);
1207   -}
1208   -
1209   -
1210   -/***********************************************************************//**
1211 1178 * @brief Save event list in FITS format.
1212 1179 *
1213 1180 * Save the event list as a FITS file. The file name of the FITS file is
... ... @@ -1325,109 +1292,3 @@ void ctselect::save_xml(void)
1325 1292 // Return
1326 1293 return;
1327 1294 }
1328   -
1329   -
1330   -/***********************************************************************//**
1331   - * @brief Save event list into FITS file
1332   - *
1333   - * @param[in] obs Pointer to CTA observation.
1334   - * @param[in] infile Input file name.
1335   - * @param[in] evtname Event extension name.
1336   - * @param[in] gtiname GTI extension name.
1337   - * @param[in] outfile Output file name.
1338   - *
1339   - * Saves an event list and the corresponding Good Time Intervals into a FITS
1340   - * file and copy all others extensions from the input file to the output
1341   - * file.
1342   - *
1343   - * If an extension name is specified in the @p outfile argument, the events
1344   - * and eventually also the Good Time Intervals will be extracted from the
1345   - * argument and used for writing the events. The format is
1346   - *
1347   - * <filename>[<event extension name;GTI extension name>]
1348   - *
1349   - * where <filename> needs to be replaced by the name of the FITS file,
1350   - * and <event extension name;GTI extension name> by the name of the events
1351   - * and Good Time Intervals extensions. For example
1352   - *
1353   - * myfits.fits[EVENTS1;GTI1]
1354   - *
1355   - * will write the selected events into the "EVENTS1" extension and the
1356   - * Good Time Intervals into the "GTI1" extension of the "myfits.fits" FITS
1357   - * file. If the Good Time Intervals extension name is skipped, e.g.
1358   - *
1359   - * myfits.fits[EVENTS1]
1360   - *
1361   - * the original extension name for the Good Time Intervals will be kept.
1362   - * Analogously, only the Good Time Intervals extension name can be changed
1363   - * by specifying
1364   - *
1365   - * myfits.fits[;GTI1]
1366   - *
1367   - * In none of the cases will the original events and Good Time Intervals be
1368   - * copied over to the output file.
1369   - ***************************************************************************/
1370   -void ctselect::save_event_list(const GCTAObservation* obs,
1371   - const std::string& infile,
1372   - const std::string& evtname,
1373   - const std::string& gtiname,
1374   - const std::string& outfile) const
1375   -{
1376   - // Save only if we have an event list
1377   - if (obs->eventtype() == "EventList") {
1378   -
1379   - // Set output FITS file event extension names
1380   - GFilename outname(outfile);
1381   - std::string outevt = evtname;
1382   - std::string outgti = gtiname;
1383   - if (outname.has_extname()) {
1384   - std::vector<std::string> extnames =
1385   - gammalib::split(outname.extname(), ";");
1386   - if (extnames.size() > 0) {
1387   - std::string extname = gammalib::strip_whitespace(extnames[0]);
1388   - if (!extname.empty()) {
1389   - outevt = extname;
1390   - }
1391   - }
1392   - if (extnames.size() > 1) {
1393   - std::string extname = gammalib::strip_whitespace(extnames[1]);
1394   - if (!extname.empty()) {
1395   - outgti = extname;
1396   - }
1397   - }
1398   - }
1399   -
1400   - // Create output FITS file
1401   - GFits outfits;
1402   -
1403   - // Write observation into FITS file
1404   - obs->write(outfits, outevt, outgti);
1405   -
1406   - // Copy all extensions other than evtname and gtiname extensions
1407   - // from the input to the output event list. The evtname and
1408   - // gtiname extensions are written by the save method, all others
1409   - // that may eventually be present have to be copied over
1410   - // explicitly.
1411   - GFits infits(infile);
1412   - for (int extno = 1; extno < infits.size(); ++extno) {
1413   - GFitsHDU* hdu = infits.at(extno);
1414   - if (hdu->extname() != evtname &&
1415   - hdu->extname() != gtiname &&
1416   - hdu->extname() != outevt &&
1417   - hdu->extname() != outgti) {
1418   - outfits.append(*hdu);
1419   - }
1420   - }
1421   -
1422   - // Close input file
1423   - infits.close();
1424   -
1425   - // Save file to disk and close it (we need both operations)
1426   - outfits.saveto(outname.url(), clobber());
1427   - outfits.close();
1428   -
1429   - } // endif: observation was unbinned
1430   -
1431   - // Return
1432   - return;
1433   -}
... ...
src/ctselect/ctselect.hpp
  View file @77328f2
1 1 /***************************************************************************
2 2 * ctselect - Data selection tool *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2017 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -77,15 +77,8 @@ protected:
77 77 std::string check_infile(const std::string& filename,
78 78 const std::string& evtname) const;
79 79 std::string set_outfile_name(const std::string& filename) const;
80   - std::string get_gtiname(const std::string& filename,
81   - const std::string& evtname) const;
82 80 void save_fits(void);
83 81 void save_xml(void);
84   - void save_event_list(const GCTAObservation* obs,
85   - const std::string& infile,
86   - const std::string& evtname,
87   - const std::string& gtiname,
88   - const std::string& outfile) const;
89 82  
90 83 // User parameters
91 84 std::string m_outobs; //!< Output event list or XML file
... ...
src/support/ctool.cpp
  View file @77328f2
... ... @@ -1730,6 +1730,145 @@ std::vector&lt;bool&gt; ctool::cube_layer_usage(const GEbounds&amp; cube_ebounds,
1730 1730  
1731 1731  
1732 1732 /***********************************************************************//**
  1733 + * @brief Save event list into FITS file
  1734 + *
  1735 + * @param[in] obs Pointer to CTA observation.
  1736 + * @param[in] infile Input file name.
  1737 + * @param[in] evtname Event extension name.
  1738 + * @param[in] gtiname GTI extension name.
  1739 + * @param[in] outfile Output file name.
  1740 + *
  1741 + * Saves an event list and the corresponding Good Time Intervals into a FITS
  1742 + * file and copy all others extensions from the input file to the output
  1743 + * file.
  1744 + *
  1745 + * If an extension name is specified in the @p outfile argument, the events
  1746 + * and eventually also the Good Time Intervals will be extracted from the
  1747 + * argument and used for writing the events. The format is
  1748 + *
  1749 + * <filename>[<event extension name;GTI extension name>]
  1750 + *
  1751 + * where <filename> needs to be replaced by the name of the FITS file,
  1752 + * and <event extension name;GTI extension name> by the name of the events
  1753 + * and Good Time Intervals extensions. For example
  1754 + *
  1755 + * myfits.fits[EVENTS1;GTI1]
  1756 + *
  1757 + * will write the selected events into the `EVENTS1` extension and the
  1758 + * Good Time Intervals into the `GTI1` extension of the `myfits.fits` FITS
  1759 + * file. If the Good Time Intervals extension name is skipped, e.g.
  1760 + *
  1761 + * myfits.fits[EVENTS1]
  1762 + *
  1763 + * the original extension name for the Good Time Intervals will be kept.
  1764 + * Analogously, only the Good Time Intervals extension name can be changed
  1765 + * by specifying
  1766 + *
  1767 + * myfits.fits[;GTI1]
  1768 + *
  1769 + * In none of the cases will the original events and Good Time Intervals be
  1770 + * copied over to the output file.
  1771 + ***************************************************************************/
  1772 +void ctool::save_event_list(const GCTAObservation* obs,
  1773 + const std::string& infile,
  1774 + const std::string& evtname,
  1775 + const std::string& gtiname,
  1776 + const std::string& outfile) const
  1777 +{
  1778 + // Save only if we have an event list
  1779 + if (obs->eventtype() == "EventList") {
  1780 +
  1781 + // Set output FITS file event extension names
  1782 + GFilename outname(outfile);
  1783 + std::string outevt = evtname;
  1784 + std::string outgti = gtiname;
  1785 + if (outname.has_extname()) {
  1786 + std::vector<std::string> extnames =
  1787 + gammalib::split(outname.extname(), ";");
  1788 + if (extnames.size() > 0) {
  1789 + std::string extname = gammalib::strip_whitespace(extnames[0]);
  1790 + if (!extname.empty()) {
  1791 + outevt = extname;
  1792 + }
  1793 + }
  1794 + if (extnames.size() > 1) {
  1795 + std::string extname = gammalib::strip_whitespace(extnames[1]);
  1796 + if (!extname.empty()) {
  1797 + outgti = extname;
  1798 + }
  1799 + }
  1800 + }
  1801 +
  1802 + // Create output FITS file
  1803 + GFits outfits;
  1804 +
  1805 + // Write observation into FITS file
  1806 + obs->write(outfits, outevt, outgti);
  1807 +
  1808 + // Copy all extensions other than evtname and gtiname extensions
  1809 + // from the input to the output event list. The evtname and
  1810 + // gtiname extensions are written by the save method, all others
  1811 + // that may eventually be present have to be copied over
  1812 + // explicitly.
  1813 + GFits infits(infile);
  1814 + for (int extno = 1; extno < infits.size(); ++extno) {
  1815 + GFitsHDU* hdu = infits.at(extno);
  1816 + if (hdu->extname() != evtname &&
  1817 + hdu->extname() != gtiname &&
  1818 + hdu->extname() != outevt &&
  1819 + hdu->extname() != outgti) {
  1820 + outfits.append(*hdu);
  1821 + }
  1822 + }
  1823 +
  1824 + // Close input file
  1825 + infits.close();
  1826 +
  1827 + // Save file to disk and close it (we need both operations)
  1828 + outfits.saveto(outname.url(), clobber());
  1829 + outfits.close();
  1830 +
  1831 + } // endif: observation was unbinned
  1832 +
  1833 + // Return
  1834 + return;
  1835 +}
  1836 +
  1837 +
  1838 +/***********************************************************************//**
  1839 + * @brief Get Good Time Intervals extension name
  1840 + *
  1841 + * @param[in] filename Input file name.
  1842 + * @param[in] evtname Events extension name.
  1843 + *
  1844 + * Extracts the Good Time Intervals extension name from the event file. We
  1845 + * do this by loading the events and accessing the Good Time Intervals
  1846 + * extension name using the GCTAEventList::gtiname() method. If the file name
  1847 + * is empty, the method returns `GTI`.
  1848 + ***************************************************************************/
  1849 +std::string ctool::get_gtiname(const std::string& filename,
  1850 + const std::string& evtname) const
  1851 +{
  1852 + // Initialise GTI name
  1853 + std::string gtiname = gammalib::extname_gti;
  1854 +
  1855 + // Continue only if the filename is not empty
  1856 + if (!filename.empty()) {
  1857 +
  1858 + // Load events
  1859 + GCTAEventList events(filename+"["+evtname+"]");
  1860 +
  1861 + // Get GTI name
  1862 + gtiname = events.gtiname();
  1863 +
  1864 + }
  1865 +
  1866 + // Return GTI name
  1867 + return (gtiname);
  1868 +}
  1869 +
  1870 +
  1871 +/***********************************************************************//**
1733 1872 * @brief Set warning string if there are too few energies
1734 1873 *
1735 1874 * @param[in] energies Energies.
... ...
src/support/ctool.hpp
  View file @77328f2
1 1 /***************************************************************************
2 2 * ctool - ctool base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2014-2016 by Juergen Knoedlseder *
  4 + * copyright (C) 2014-2017 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -121,6 +121,13 @@ protected:
121 121 std::vector<bool> cube_layer_usage(const GEbounds& cube_ebounds,
122 122 const GEbounds& list_ebounds) const;
123 123 bool is_valid_filename(const GFilename& filename) const;
  124 + std::string get_gtiname(const std::string& filename,
  125 + const std::string& evtname) const;
  126 + void save_event_list(const GCTAObservation* obs,
  127 + const std::string& infile,
  128 + const std::string& evtname,
  129 + const std::string& gtiname,
  130 + const std::string& outfile) const;
124 131  
125 132 // Protected warning strings
126 133 std::string warn_too_few_energies(const GEnergies& energies) const;
... ...
test/test_ctprob.py 0 → 100755
  View file @77328f2
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This scripts performs unit tests for the ctprob tool.
  4 +#
  5 +# Copyright (C) 2017 by Leonardo Di Venere
  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 os
  22 +import gammalib
  23 +import ctools
  24 +from testing import test
  25 +
  26 +
  27 +# ============================ #
  28 +# Test class for ctprob tool #
  29 +# ============================ #
  30 +class Test(test):
  31 + """
  32 + Test class for ctprob tool
  33 + """
  34 + # Constructor
  35 + def __init__(self):
  36 + """
  37 + Constructor
  38 + """
  39 + # Call base class constructor
  40 + test.__init__(self)
  41 +
  42 + # Set test data
  43 + self._invalid_events = self._datadir + '/invalid_event_list.fits'
  44 +
  45 + # Return
  46 + return
  47 +
  48 + # Set test functions
  49 + def set(self):
  50 + """
  51 + Set all test functions
  52 + """
  53 + # Set test name
  54 + self.name('ctprob')
  55 +
  56 + # Append tests
  57 + self.append(self._test_cmd, 'Test ctprob on command line')
  58 + self.append(self._test_python, 'Test ctprob from Python')
  59 +
  60 + # Return
  61 + return
  62 +
  63 + # Test ctprob on command line
  64 + def _test_cmd(self):
  65 + """
  66 + Test ctprob on the command line
  67 + """
  68 + # Set tool name
  69 + ctprob = self._tool('ctprob')
  70 +
  71 + # Setup ctprob command
  72 + cmd = ctprob+' inobs="'+self._events+'"'+ \
  73 + ' outobs="ctprob_cmd1.fits"'+ \
  74 + ' inmodel='+self._model+ \
  75 + ' caldb='+self._caldb+' irf='+self._irf+ \
  76 + ' edisp=no'+ \
  77 + ' logfile="ctprob_cmd1.log" chatter=1'
  78 +
  79 + # Check if execution was successful
  80 + self.test_assert(self._execute(cmd) == 0,
  81 + 'Check successful execution from command line')
  82 +
  83 + # Check result file
  84 + self._check_result_file('ctprob_cmd1.fits')
  85 +
  86 + # Setup ctprob command
  87 + cmd = ctprob+' inobs="event_file_that_does_not_exist.fits"'+ \
  88 + ' outobs="ctprob_cmd2.fits"'+ \
  89 + ' inmodel='+self._model+ \
  90 + ' caldb='+self._caldb+' irf='+self._irf+ \
  91 + ' edisp=no'+ \
  92 + ' logfile="ctprob_cmd2.log" chatter=1'
  93 +
  94 + # Check if execution failed
  95 + self.test_assert(self._execute(cmd) != 0,
  96 + 'Check invalid input file when executed from command line')
  97 +
  98 + # Setup ctprob --help option
  99 + cmd = ctprob+' --help'
  100 +
  101 + # Check if execution was successful in case that the CTOOLS
  102 + # environment variable was set or failed otherwise
  103 + if 'CTOOLS' in os.environ:
  104 + self.test_value(self._execute(cmd), 0,
  105 + 'Check successful execution with --help option')
  106 + else:
  107 + self.test_assert(self._execute(cmd) != 0,
  108 + 'Check execution failure with --help option')
  109 +
  110 + # Return
  111 + return
  112 +
  113 + # Test ctprob from Python
  114 + def _test_python(self):
  115 + """
  116 + Test ctprob from Python
  117 + """
  118 + # Allocate empty ctprob tool
  119 + prob = ctools.ctprob()
  120 +
  121 + # Check that empty ctprob tool holds an empty observation
  122 + self._check_obs(prob.obs(), nobs=0)
  123 +
  124 + # Check that saving does nothing
  125 + prob['outobs'] = 'ctprob_py0.fits'
  126 + prob['logfile'] = 'ctprob_py0.log'
  127 + prob.logFileOpen()
  128 + prob.save()
  129 + self.test_assert(not os.path.isfile('ctprob_py0.fits'),
  130 + 'Check that no event list has been created')
  131 +
  132 + # Check that clearing does not lead to an exception or segfault
  133 + prob.clear()
  134 +
  135 + # Now set ctprob parameters
  136 + prob['inobs'] = self._events
  137 + prob['outobs'] = 'ctprob_py1.fits'
  138 + prob['inmodel'] = self._model
  139 + prob['caldb'] = self._caldb
  140 + prob['irf'] = self._irf
  141 + prob['edisp'] = False
  142 + prob['logfile'] = 'ctprob_py1.log'
  143 + prob['chatter'] = 2
  144 +
  145 + # Run ctprob tool
  146 + prob.logFileOpen() # Make sure we get a log file
  147 + prob.run()
  148 + prob.save()
  149 +
  150 + # Check result file
  151 + self._check_result_file('ctprob_py1.fits')
  152 +
  153 + # Setup ctprob tool for an invalid event file
  154 + prob = ctools.ctprob()
  155 + prob['inobs'] = self._invalid_events
  156 + prob['outobs'] = 'ctprob_py2.fits'
  157 + prob['inmodel'] = self._model
  158 + prob['caldb'] = self._caldb
  159 + prob['irf'] = self._irf
  160 + prob['edisp'] = False
  161 + prob['logfile'] = 'ctprob_py2.log'
  162 + prob['chatter'] = 3
  163 +
  164 + # Execute tool
  165 + prob.logFileOpen() # Needed to get a new log file
  166 + self.test_try('Test invalid event file')
  167 + try:
  168 + prob.execute()
  169 + self.test_try_failure('Exception not thrown')
  170 + except RuntimeError:
  171 + self.test_try_success()
  172 +
  173 + # Setup ctprob with a list of observations and test prefix
  174 + prob = ctools.ctprob()
  175 + prob['inobs'] = self._datadir + '/obs_unbinned.xml'
  176 + prob['prefix'] = 'ctprob_py3_'
  177 + prob['outobs'] = 'ctprob_py3_obs.xml'
  178 + prob['inmodel'] = self._model
  179 + prob['caldb'] = self._caldb
  180 + prob['irf'] = self._irf
  181 + prob['edisp'] = True
  182 + prob['logfile'] = 'ctprob_py3.log'
  183 + prob['chatter'] = 4
  184 +
  185 + # Run ctprob tool
  186 + prob.logFileOpen() # Make sure we get a log file
  187 + prob.execute()
  188 +
  189 + # Check result file
  190 + self._check_result_file('ctprob_py3_crab_events.fits')
  191 +
  192 + # Return
  193 + return
  194 +
  195 + # Check result file
  196 + def _check_result_file(self, filename, nevents=6141):
  197 + """
  198 + Check result file
  199 +
  200 + Parameters
  201 + ----------
  202 + filename : str
  203 + Event list file name
  204 + nevents : int, optional
  205 + Expected number of events
  206 + """
  207 + # Open result file
  208 + events = gammalib.GCTAEventList(filename)
  209 +
  210 + # Check event list
  211 + self._check_events(events, nevents=nevents)
  212 +
  213 + # Check that probability columns exist
  214 + colname1 = 'PROB_Background'
  215 + colname2 = 'PROB_Crab'
  216 + self._check_column(filename, colname1)
  217 + self._check_column(filename, colname2)
  218 +
  219 + # Check that the probability values are correctly normalized
  220 + self._check_normalization(filename, [colname1,colname2])
  221 +
  222 + # Return
  223 + return
  224 +
  225 + # Check observation and event list
  226 + def _check_obs(self, obs, nobs=1, nevents=6141):
  227 + """
  228 + Check observation and event list
  229 +
  230 + Parameters
  231 + ----------
  232 + obs : `~gammalib.GObservations`
  233 + Models
  234 + nobs : int, optional
  235 + Expected number of observations
  236 + nevents : int, optional
  237 + Expected number of events
  238 + """
  239 + # Check number of observations
  240 + self.test_value(obs.size(), nobs, 'Check number of observations')
  241 +
  242 + # If there is an observation then check the event list of the first
  243 + # one
  244 + if obs.size() > 0:
  245 + self._check_events(obs[0].events(), nevents=nevents)
  246 +
  247 + # Return
  248 + return
  249 +
  250 + # Check events
  251 + def _check_events(self, events, nevents=6141):
  252 + """
  253 + Check event list
  254 +
  255 + Parameters
  256 + ----------
  257 + events : `~gammalib.GCTAEventList`
  258 + Event list
  259 + nevents : int, optional
  260 + Expected number of events
  261 + """
  262 + # Check event list
  263 + self.test_value(events.size(), nevents, 'Check number of events')
  264 +
  265 + # Return
  266 + return
  267 +
  268 +
  269 + def _check_column(self, filename, colname):
  270 + """
  271 + Check that column colname exists.
  272 +
  273 + Parameters
  274 + ----------
  275 + filename : str
  276 + Event list file name
  277 + colname : str
  278 + Name of the column to check
  279 + """
  280 + # Check that column exists
  281 + gfits = gammalib.GFits(filename)
  282 + gtable = gfits[1]
  283 + check = int(gtable.contains(colname))
  284 + self.test_value(check, 1, 'Check column '+colname)
  285 +
  286 + # Return
  287 + return
  288 +
  289 +
  290 + def _check_normalization(self, filename, colnames):
  291 + """
  292 + Check that the probability values in columns identified by colnames are
  293 + correctly normalized.
  294 +
  295 + Parameters
  296 + ----------
  297 + filename : str
  298 + Event list file name
  299 + colnames : list
  300 + Names of the columns to check
  301 + """
  302 + # Check that the probability values are correctly normalized
  303 +
  304 + # Open fits table
  305 + gfits = gammalib.GFits(filename)
  306 + gtable = gfits[1]
  307 + ncols = gtable.ncols()
  308 + nevt = gtable.nrows()
  309 + srcs = []
  310 +
  311 + # Loop over colnames
  312 + for colname in colnames:
  313 +
  314 + # Get the correct column corresponding to colname from the table
  315 + for icol in range(ncols):
  316 + gcol = gtable[icol]
  317 + if gcol.name()==colname:
  318 + break
  319 +
  320 + # Read and save column content
  321 + src1 = []
  322 + for ievt in range(nevt):
  323 + src1.append(gcol[ievt])
  324 + srcs.extend(src1)
  325 +
  326 + # Compute sum
  327 + tot = sum(srcs)
  328 + self.test_value(int(tot+0.5), nevt, 'Check that probability columns '
  329 + 'are normalized')
  330 +
  331 + # Return
  332 + return
  333 +
... ...
test/test_python_ctools.py
  View file @77328f2
... ... @@ -39,6 +39,7 @@ import test_ctcubemask
39 39 import test_ctbutterfly
40 40 import test_ctulimit
41 41 import test_cterror
  42 +import test_ctprob
42 43 import test_pipelines
43 44  
44 45  
... ... @@ -116,6 +117,7 @@ def test(installed=False):
116 117 test_ctbutterfly.Test(),
117 118 test_ctulimit.Test(),
118 119 test_cterror.Test(),
  120 + test_ctprob.Test(),
119 121 test_pipelines.Test()]
120 122  
121 123 # Allocate test suite container
... ...