Johan Bregeon / ctools

Commit 6762b992b9297f953c4537cb1053d8a08d6504a4

Authored by Michael Mayer
1 parent bab5390b

added structure for ctulimit

Inline Side-by-side
Showing 7 changed files with 675 additions and 2 deletions
pyext/ctools/tools.i
  View file @6762b99
... ... @@ -58,3 +58,4 @@
58 58 %include "ctcubemask.i"
59 59 %include "cttsmap.i"
60 60 %include "ctbutterfly.i"
  61 +%include "ctulimit.i"
... ...
src/Makefile.am
  View file @6762b99
... ... @@ -17,7 +17,8 @@ SUBDIRS = support \
17 17 ctbkgcube \
18 18 ctcubemask \
19 19 cttsmap \
20   - ctbutterfly
  20 + ctbutterfly \
  21 + ctulimit
21 22  
22 23 # Scripts to be distributed
23 24 scripts = ctools-setup \
... ... @@ -41,7 +42,8 @@ libctools_la_LIBADD = support/libsupport.la \
41 42 ctbkgcube/libctbkgcube.la \
42 43 ctcubemask/libctcubemask.la \
43 44 cttsmap/libcttsmap.la \
44   - ctbutterfly/libctbutterfly.la
  45 + ctbutterfly/libctbutterfly.la \
  46 + ctulimit/libctulimit.la
45 47  
46 48  
47 49 # Distribute scripts
... ...
src/ctulimit/Makefile.am 0 → 100644
  View file @6762b99
  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 = ctulimit
  8 +
  9 +# Sources
  10 +sources = ctulimit.cpp
  11 +
  12 +# Libtool library
  13 +noinst_LTLIBRARIES = libctulimit.la
  14 +libctulimit_la_SOURCES = $(sources)
  15 +
  16 +# Program definition
  17 +ctulimit_SOURCES = main.cpp
  18 +ctulimit_LDFLAGS = $(LIB_LINK_OPTION) @LDFLAGS@
  19 +ctulimit_CPPFLAGS = @CPPFLAGS@
  20 +ctulimit_LDADD = libctulimit.la $(top_srcdir)/src/support/libsupport.la
  21 +
  22 +# Parameter files
  23 +pfilesdir = $(prefix)/syspfiles
  24 +pfiles_DATA = ctulimit.par
... ...
src/ctulimit/ctulimit.cpp 0 → 100644
  View file @6762b99
  1 +/***************************************************************************
  2 + * ctulimitt - upper limit 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 ctulimit.hpp
  23 + * @brief upper limit calculation tool interface implementation
  24 + * @author Michael Mayer
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include <cstdio>
  32 +#include "ctulimit.hpp"
  33 +#include "GTools.hpp"
  34 +
  35 +/* __ Method name definitions ____________________________________________ */
  36 +#define G_GET_PARAMETERS "ctulimit::get_parameters()"
  37 +
  38 +/* __ Debug definitions __________________________________________________ */
  39 +
  40 +/* __ Coding definitions _________________________________________________ */
  41 +
  42 +
  43 +/*==========================================================================
  44 + = =
  45 + = Constructors/destructors =
  46 + = =
  47 + ==========================================================================*/
  48 +
  49 +/***********************************************************************//**
  50 + * @brief Void constructor
  51 + ***************************************************************************/
  52 +ctulimit::ctulimit(void) : ctool(CTULIMIT_NAME, CTULIMIT_VERSION)
  53 +{
  54 + // Initialise members
  55 + init_members();
  56 +
  57 + // Return
  58 + return;
  59 +}
  60 +
  61 +
  62 +/***********************************************************************//**
  63 + * @brief Observations constructor
  64 + *
  65 + * param[in] obs Observation container.
  66 + *
  67 + * This method creates an instance of the class by copying an existing
  68 + * observations container.
  69 + ***************************************************************************/
  70 +ctulimit::ctulimit(const GObservations& obs) :
  71 + ctool(CTULIMIT_NAME, CTULIMIT_VERSION)
  72 +{
  73 + // Initialise members
  74 + init_members();
  75 +
  76 + // Set observations
  77 + m_obs = obs;
  78 +
  79 + // Return
  80 + return;
  81 +}
  82 +
  83 +
  84 +
  85 +/***********************************************************************//**
  86 + * @brief Command line constructor
  87 + *
  88 + * @param[in] argc Number of arguments in command line.
  89 + * @param[in] argv Array of command line arguments.
  90 + ***************************************************************************/
  91 +ctulimit::ctulimit(int argc, char *argv[]) :
  92 + ctool(CTULIMIT_NAME, CTULIMIT_VERSION, argc, argv)
  93 +{
  94 + // Initialise members
  95 + init_members();
  96 +
  97 + // Return
  98 + return;
  99 +}
  100 +
  101 +
  102 +/***********************************************************************//**
  103 + * @brief Copy constructor
  104 + *
  105 + * @param[in] app Application.
  106 + ***************************************************************************/
  107 +ctulimit::ctulimit(const ctulimit& app) : ctool(app)
  108 +{
  109 + // Initialise members
  110 + init_members();
  111 +
  112 + // Copy members
  113 + copy_members(app);
  114 +
  115 + // Return
  116 + return;
  117 +}
  118 +
  119 +
  120 +/***********************************************************************//**
  121 + * @brief Destructor
  122 + ***************************************************************************/
  123 +ctulimit::~ctulimit(void)
  124 +{
  125 + // Free members
  126 + free_members();
  127 +
  128 + // Return
  129 + return;
  130 +}
  131 +
  132 +
  133 +/*==========================================================================
  134 + = =
  135 + = Operators =
  136 + = =
  137 + ==========================================================================*/
  138 +
  139 +/***********************************************************************//**
  140 + * @brief Assignment operator
  141 + *
  142 + * @param[in] app Application.
  143 + * @return Application.
  144 + ***************************************************************************/
  145 +ctulimit& ctulimit::operator=(const ctulimit& app)
  146 +{
  147 + // Execute only if object is not identical
  148 + if (this != &app) {
  149 +
  150 + // Copy base class members
  151 + this->ctool::operator=(app);
  152 +
  153 + // Free members
  154 + free_members();
  155 +
  156 + // Initialise members
  157 + init_members();
  158 +
  159 + // Copy members
  160 + copy_members(app);
  161 +
  162 + } // endif: object was not identical
  163 +
  164 + // Return this object
  165 + return *this;
  166 +}
  167 +
  168 +
  169 +/*==========================================================================
  170 + = =
  171 + = Public methods =
  172 + = =
  173 + ==========================================================================*/
  174 +
  175 +/***********************************************************************//**
  176 + * @brief Clear instance
  177 + ***************************************************************************/
  178 +void ctulimit::clear(void)
  179 +{
  180 + // Free members
  181 + free_members();
  182 + this->ctool::free_members();
  183 + this->GApplication::free_members();
  184 +
  185 + // Initialise members
  186 + this->GApplication::init_members();
  187 + this->ctool::init_members();
  188 + init_members();
  189 +
  190 + // Return
  191 + return;
  192 +}
  193 +
  194 +
  195 +/***********************************************************************//**
  196 + * @brief Computes the upper limit
  197 + *
  198 + * This method calculates the upper limit depending on the given algorithm
  199 + ***************************************************************************/
  200 +void ctulimit::run(void)
  201 +{
  202 + // If we're in debug mode then all output is also dumped on the screen
  203 + if (logDebug()) {
  204 + log.cout(true);
  205 + }
  206 +
  207 + // Get task parameters
  208 + get_parameters();
  209 +
  210 + // Write parameters into logger
  211 + if (logTerse()) {
  212 + log_parameters();
  213 + log << std::endl;
  214 + }
  215 +
  216 + // Write observation(s) into logger
  217 + if (logTerse()) {
  218 + log << std::endl;
  219 + if (m_obs.size() > 1) {
  220 + log.header1("Observations");
  221 + }
  222 + else {
  223 + log.header1("Observation");
  224 + }
  225 + log << m_obs << std::endl;
  226 + }
  227 +
  228 + // Write header
  229 + if (logTerse()) {
  230 + log << std::endl;
  231 + log.header1("Compute upper limit");
  232 + }
  233 +
  234 +
  235 + // Compute upper limit here
  236 +
  237 + // Return
  238 + return;
  239 +}
  240 +
  241 +
  242 +/***********************************************************************//**
  243 + * @brief Save maps
  244 + *
  245 + * This method saves the upper limit to an ascii file
  246 + ***************************************************************************/
  247 +void ctulimit::save(void)
  248 +{
  249 + // Write header
  250 + if (logTerse()) {
  251 + log << std::endl;
  252 + log.header1("Save TS map");
  253 + }
  254 +
  255 + // Get output filename
  256 + m_outfile = (*this)["outfile"].filename();
  257 +
  258 +
  259 + // Return
  260 + return;
  261 +}
  262 +
  263 +
  264 +/*==========================================================================
  265 + = =
  266 + = Private methods =
  267 + = =
  268 + ==========================================================================*/
  269 +
  270 +/***********************************************************************//**
  271 + * @brief Initialise class members
  272 + ***************************************************************************/
  273 +void ctulimit::init_members(void)
  274 +{
  275 + // Initialise members
  276 + m_infile.clear();
  277 + m_outfile.clear();
  278 +
  279 +
  280 + // Initialise protected members
  281 + m_obs.clear();
  282 +
  283 +
  284 + // Return
  285 + return;
  286 +}
  287 +
  288 +
  289 +/***********************************************************************//**
  290 + * @brief Copy class members
  291 + *
  292 + * @param[in] app Application.
  293 + ***************************************************************************/
  294 +void ctulimit::copy_members(const ctulimit& app)
  295 +{
  296 + // Copy attributes
  297 + m_infile = app.m_infile;
  298 + m_outfile = app.m_outfile;
  299 + m_modelfile = app.m_modelfile;
  300 + m_srcname = app.m_srcname;
  301 + m_parname = app.m_parname;
  302 + m_algorithm = app.m_algorithm;
  303 +
  304 + // Copy protected members
  305 + m_obs = app.m_obs;
  306 +
  307 + // Return
  308 + return;
  309 +}
  310 +
  311 +
  312 +/***********************************************************************//**
  313 + * @brief Delete class members
  314 + ***************************************************************************/
  315 +void ctulimit::free_members(void)
  316 +{
  317 + // Return
  318 + return;
  319 +}
  320 +
  321 +
  322 +/***********************************************************************//**
  323 + * @brief Get application parameters
  324 + *
  325 + * @exception GException::invalid_value
  326 + * Test source not found or no RA/DEC parameters found for test
  327 + * source.
  328 + *
  329 + * Get all task parameters from parameter file or (if required) by querying
  330 + * the user. Most parameters are only required if no observation exists so
  331 + * far in the observation container. In this case, a single CTA observation
  332 + * will be added to the container, using the definition provided in the
  333 + * parameter file.
  334 + ***************************************************************************/
  335 +void ctulimit::get_parameters(void)
  336 +{
  337 + // If there are no observations in container then add a single CTA
  338 + // observation using the parameters from the parameter file
  339 + if (m_obs.size() == 0) {
  340 +
  341 + // Allocate CTA observation
  342 + GCTAObservation obs;
  343 +
  344 + // Get event file name
  345 + std::string filename = (*this)["infile"].filename();
  346 +
  347 + // Try first to open as FITS file
  348 + try {
  349 +
  350 + // Load data
  351 + obs.load(filename);
  352 +
  353 + // Set response
  354 + set_obs_response(&obs);
  355 +
  356 + // Append observation to container
  357 + m_obs.append(obs);
  358 +
  359 + }
  360 +
  361 + // ... otherwise try to open as XML file
  362 + catch (GException::fits_open_error &e) {
  363 +
  364 + // Load observations from XML file
  365 + m_obs.load(filename);
  366 +
  367 + // Check if all observations have response information. If
  368 + // not, get the calibration database parameters and set
  369 + // the response properly
  370 + set_response(m_obs);
  371 +
  372 + } // endcatch: file was an XML file
  373 +
  374 + } // endif: there was no observation in the container
  375 +
  376 + // If there is are no models associated with the observations then
  377 + // load now the model definition
  378 + if (m_obs.models().size() == 0) {
  379 +
  380 + // Get models XML filename
  381 + std::string filename = (*this)["srcmdl"].filename();
  382 +
  383 + // Setup models for optimizing.
  384 + m_obs.models(GModels(filename));
  385 +
  386 + } // endif: no models were associated with observations
  387 +
  388 + // Get name of test source and check container for this name
  389 + m_srcname = (*this)["srcname"].string();
  390 + if (!m_obs.models().contains(m_srcname)) {
  391 + std::string msg = "Source \""+m_srcname+"\" not found in model "
  392 + "container. Please add a source with that name "
  393 + "or check for a possible typos.";
  394 + throw GException::invalid_value(G_GET_PARAMETERS, msg);
  395 + }
  396 +
  397 +
  398 + // Optionally read ahead parameters so that they get correctly
  399 + // dumped into the log file
  400 + if (read_ahead()) {
  401 + m_outfile = (*this)["outfile"].filename();
  402 + }
  403 +
  404 + // TODO read other parameters here
  405 +
  406 + // Return
  407 + return;
  408 +}
  409 +
  410 +
... ...
src/ctulimit/ctulimit.hpp 0 → 100644
  View file @6762b99
  1 +/***************************************************************************
  2 + * ctulimit - Upper limit 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 ctulimit.hpp
  23 + * @brief upper limit calculation tool interface definition
  24 + * @author Michael Mayer
  25 + */
  26 +
  27 +#ifndef CTULIMIT_HPP
  28 +#define CTULIMIT_HPP
  29 +
  30 +/* __ Includes ___________________________________________________________ */
  31 +#include "GammaLib.hpp"
  32 +#include "GCTALib.hpp"
  33 +#include "ctool.hpp"
  34 +
  35 +/* __Definitions _________________________________________________________ */
  36 +#define CTULIMIT_NAME "ctulimit"
  37 +#define CTULIMIT_VERSION "00-01-00"
  38 +
  39 +
  40 +/***********************************************************************//**
  41 + * @class ctulimit
  42 + *
  43 + * @brief Upper limit calculation tool
  44 + *
  45 + * This class computes and upper limit for a given parameter.
  46 + *
  47 + * The class operates on predefined observation containers, an individual
  48 + * event list or an observation definition XML file.
  49 + *
  50 + * During the computation the likelihood function is inspected to find the best
  51 + * value for the upper limit.
  52 + *
  53 + ***************************************************************************/
  54 +class ctulimit : public ctool {
  55 +
  56 +public:
  57 + // Constructors and destructors
  58 + ctulimit(void);
  59 + explicit ctulimit(const GObservations& obs);
  60 + ctulimit(int argc, char *argv[]);
  61 + ctulimit(const ctulimit& app);
  62 + virtual ~ctulimit(void);
  63 +
  64 + // Operators
  65 + ctulimit& operator=(const ctulimit& app);
  66 +
  67 + // Methods
  68 + void clear(void);
  69 + void run(void);
  70 + void save(void);
  71 + const GObservations& obs(void) const;
  72 + const double& ulimit(void) const;
  73 +
  74 +protected:
  75 + // Protected methods
  76 + void init_members(void);
  77 + void copy_members(const ctulimit& app);
  78 + void free_members(void);
  79 + void get_parameters(void);
  80 +
  81 + // User parameters
  82 + std::string m_infile; //!< Input file
  83 + std::string m_modelfile;//!< Name of model optimised model file
  84 + std::string m_srcname; //!< Name of source which is moved around
  85 + std::string m_parname; //!< Name of parameter upper limit should be computed
  86 + std::string m_outfile; //!< Output counts map or XML file
  87 + std::string m_algorithm; //!< algorithm of upper limit computation
  88 +
  89 + // Protected members
  90 + GObservations m_obs; //!< Observation container
  91 +
  92 + double m_ulimit; //!<Upper limit value
  93 +
  94 +};
  95 +
  96 +
  97 +/***********************************************************************//**
  98 + * @brief Return observation container
  99 + *
  100 + * @return Reference to observation container
  101 + ***************************************************************************/
  102 +inline
  103 +const GObservations& ctulimit::obs(void) const
  104 +{
  105 + return m_obs;
  106 +}
  107 +
  108 +/***********************************************************************//**
  109 + * @brief Return observation container
  110 + *
  111 + * @return Reference to observation container
  112 + ***************************************************************************/
  113 +inline
  114 +const double& ctulimit::ulimit(void) const
  115 +{
  116 + return m_ulimit;
  117 +}
  118 +
  119 +
  120 +#endif /* CTULIMIT_HPP */
... ...
src/ctulimit/ctulimit.par 0 → 100644
  View file @6762b99
  1 +#############################################################################
  2 +# #
  3 +# CTA Science Analysis Tools #
  4 +# #
  5 +# ------------------------------------------------------------------------- #
  6 +# #
  7 +# File: ctulimit.par #
  8 +# Version: 00-01-00 #
  9 +# #
  10 +# Author: Michael Mayer #
  11 +# michael.mayer@physik.hu-berlin.de #
  12 +# DESY #
  13 +# #
  14 +# Purpose: Parameter file for the ctulimit 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",,,"Test source"
  23 +parname, s, a, "Prefactor",,,"Parameter name for which upper limit should be computed"
  24 +algorithm, s, a, "Bayesian",,,"Limit calculation algorithm"
  25 +caldb, s, a, "dummy",,, "Calibration database"
  26 +irf, s, a, "cta_dummy_irf",,, "Instrument response function"
  27 +outfile, f, a, "crab_limit.txt",,, "Output file"
  28 +
  29 +
  30 +#
  31 +# Standard parameters
  32 +#====================
  33 +chatter, i, h, 2,0,4, "Chattiness of output"
  34 +clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
  35 +debug, b, h, no,,, "Debugging mode activated"
  36 +mode, s, h, "ql",,, "Mode of automatic parameters"
  37 +logfile, f, h, "cttsmap.log",,, "Log filename"
... ...
src/ctulimit/main.cpp 0 → 100644
  View file @6762b99
  1 +/***************************************************************************
  2 + * ctulimit - upper limit 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 ctulimit/main.cpp
  23 + * @brief CTA upper limit generation
  24 + * @author Michael Mayer
  25 + */
  26 +
  27 +/* __ Includes ___________________________________________________________ */
  28 +#ifdef HAVE_CONFIG_H
  29 +#include <config.h>
  30 +#endif
  31 +#include <cstdio>
  32 +#include "ctulimit.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 ctulimit application. It allocates a
  42 + * ctupperlimit 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 + ctulimit 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 + " cttsmap. 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 +}
... ...