Dmitry Malyshev / ctools

Commit e511b18718ffd6adc799c3c78a9458024b06cd5e

Authored by Jürgen Knödlseder
2 parents cd821f28 411b7316

Merge branch '1363-create-ctulimit' into integration

Inline Side-by-side
Showing 14 changed files with 1281 additions and 2 deletions
configure.ac
  View file @e511b18
... ... @@ -888,6 +888,7 @@ AC_CONFIG_FILES([Makefile
888 888 src/ctcubemask/Makefile
889 889 src/cttsmap/Makefile
890 890 src/ctbutterfly/Makefile
  891 + src/ctulimit/Makefile
891 892 src/support/Makefile
892 893 src/ctools-setup
893 894 test/Makefile
... ...
doc/source/reference_manual/ctulimit.rst 0 → 100644
  View file @e511b18
  1 +.. _ctulimit:
  2 +
  3 +ctulimit
  4 +===========
  5 +
  6 +Computes upper limit for a specific sky model component.
  7 +
  8 +
  9 +Synopsis
  10 +--------
  11 +
  12 +Computes the upper limit flux of a given source and stores the value
  13 +in an ASCII file.
  14 +
  15 +
  16 +General parameters
  17 +------------------
  18 +
  19 +``inobs [file]``
  20 + Input event list, counts cube or observation definition XML file.
  21 +
  22 +``inmodel [file]``
  23 + Model XML file containing the source and background definitions.
  24 +
  25 +``srcname [string]``
  26 + Name of model component for which upper limit should be computed.
  27 +
  28 +``expcube = NONE [file]``
  29 + Exposure cube file (only needed for stacked analysis).
  30 +
  31 +``psfcube = NONE [file]``
  32 + PSF cube file (only needed for stacked analysis).
  33 +
  34 +``caldb [string]``
  35 + Calibration database.
  36 +
  37 +``irf [string]``
  38 + Instrumental response function.
  39 +
  40 +``outfile [file]``
  41 + Output ASCII file name.
  42 +
  43 +``(confidence = 0.95) [real]``
  44 + Confidence Level of upper limit.
  45 +
  46 +``(sigma_min = 0.0) [real]``
  47 + Minimum boundary to start searching for upper limit value.
  48 + Number of standard deviations above best fit value
  49 +
  50 +``(sigma_max = 10.0) [real]``
  51 + Maximum boundary to start searching for upper limit value.
  52 + Number of standard deviations above best fit value
  53 +
  54 +``(emin = 1.0) [real])``
  55 + Minimum energy of flux limits (in TeV).
  56 +
  57 +``(emax = 100) [real]``
  58 + Maximum energy of flux limits (in TeV).
  59 +
  60 +``(tol = 1e-5) [real]``
  61 + Computation tolerance.
  62 +
  63 +``(max_iter = 50) [integer]``
  64 + Maximum number of iterations before stopping the upper
  65 + limit computations.
  66 +
  67 +
  68 +Standard parameters
  69 +-------------------
  70 +
  71 +``(chatter = 2) [integer]``
  72 + Verbosity of the executable:
  73 + chatter = 0: no information will be logged
  74 +
  75 + chatter = 1: only errors will be logged
  76 +
  77 + chatter = 2: errors and actions will be logged
  78 +
  79 + chatter = 3: report about the task execution
  80 +
  81 + chatter = 4: detailed report about the task execution
  82 +
  83 +``(clobber = yes) [boolean]``
  84 + Specifies whether an existing output counts cube should be overwritten.
  85 +
  86 +``(debug = no) [boolean]``
  87 + Enables debug mode. In debug mode the executable will dump any log file output to the console.
  88 +
  89 +``(mode = ql) [string]``
  90 + Mode of automatic parameters (default is "ql", i.e. "query and learn").
  91 +
  92 +``(logfile = ctulimit.log) [string]``
  93 + Name of log file.
  94 +
  95 +
  96 +Related tools
  97 +-------------
  98 +
  99 +None
... ...
doc/source/reference_manual/index.rst
  View file @e511b18
... ... @@ -27,6 +27,7 @@ ctools
27 27 ctselect --- Selects event data <ctselect>
28 28 ctskymap --- Generates sky map <ctskymap>
29 29 cttsmap --- Generates Test Statistics map <cttsmap>
  30 + ctulimit --- Calculates upper limit <ctulimit>
30 31  
31 32  
32 33 cscripts
... ...
pyext/ctools/tools.i
  View file @e511b18
... ... @@ -58,3 +58,4 @@
58 58 %include "ctcubemask.i"
59 59 %include "cttsmap.i"
60 60 %include "ctbutterfly.i"
  61 +%include "ctulimit.i"
... ...
pyext/ctulimit.i 0 → 100644
  View file @e511b18
  1 +/***************************************************************************
  2 + * ctulimit - Upper limit calculation tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2015 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
  23 + * @brief Upper limit calculation tool
  24 + * @author Michael Mayer
  25 + */
  26 +%{
  27 +/* Put headers and other declarations here that are needed for compilation */
  28 +#include "ctulimit.hpp"
  29 +%}
  30 +
  31 +
  32 +/***********************************************************************//**
  33 + * @class ctulimit
  34 + *
  35 + * @brief Upper limit calculation tool
  36 + ***************************************************************************/
  37 +class ctulimit : public ctool {
  38 +public:
  39 + // Constructors and destructors
  40 + ctulimit(void);
  41 + explicit ctulimit(const GObservations& obs);
  42 + ctulimit(int argc, char *argv[]);
  43 + ctulimit(const ctulimit& app);
  44 + virtual ~ctulimit(void);
  45 +
  46 + // Methods
  47 + void clear(void);
  48 + void run(void);
  49 + void save(void);
  50 + const GObservations& obs(void) const;
  51 + const double& diff_ulimit(void) const;
  52 + const double& flux_ulimit(void) const;
  53 + const double& eflux_ulimit(void) const;
  54 +};
  55 +
  56 +
  57 +/***********************************************************************//**
  58 + * @brief Upper limit calculation tool Python extensions
  59 + ***************************************************************************/
  60 +%extend ctulimit {
  61 + ctulimit copy() {
  62 + return (*self);
  63 + }
  64 +}
... ...
src/Makefile.am
  View file @e511b18
... ... @@ -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 @e511b18
  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 @e511b18
  1 +/***************************************************************************
  2 + * ctulimit - Upper limit calculation tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2015 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 +#include "GOptimizer.hpp"
  35 +
  36 +/* __ Method name definitions ____________________________________________ */
  37 +#define G_GET_PARAMETERS "ctulimit::get_parameters()"
  38 +#define G_GET_MODEL_PARAMETER "ctulimit::get_model_parameter()"
  39 +#define G_UL_BISECTION "ctulimit::ul_bisection(double&, double&)"
  40 +#define G_EVALUATE "ctulimit::evaluate(double&)"
  41 +
  42 +/* __ Debug definitions __________________________________________________ */
  43 +
  44 +/* __ Coding definitions _________________________________________________ */
  45 +
  46 +
  47 +/*==========================================================================
  48 + = =
  49 + = Constructors/destructors =
  50 + = =
  51 + ==========================================================================*/
  52 +
  53 +/***********************************************************************//**
  54 + * @brief Void constructor
  55 + ***************************************************************************/
  56 +ctulimit::ctulimit(void) : ctool(CTULIMIT_NAME, CTULIMIT_VERSION)
  57 +{
  58 + // Initialise members
  59 + init_members();
  60 +
  61 + // Return
  62 + return;
  63 +}
  64 +
  65 +
  66 +/***********************************************************************//**
  67 + * @brief Observations constructor
  68 + *
  69 + * param[in] obs Observation container.
  70 + *
  71 + * This method creates an instance of the class by copying an existing
  72 + * observations container.
  73 + ***************************************************************************/
  74 +ctulimit::ctulimit(const GObservations& obs) :
  75 + ctool(CTULIMIT_NAME, CTULIMIT_VERSION)
  76 +{
  77 + // Initialise members
  78 + init_members();
  79 +
  80 + // Set observations
  81 + m_obs = obs;
  82 +
  83 + // Return
  84 + return;
  85 +}
  86 +
  87 +
  88 +
  89 +/***********************************************************************//**
  90 + * @brief Command line constructor
  91 + *
  92 + * @param[in] argc Number of arguments in command line.
  93 + * @param[in] argv Array of command line arguments.
  94 + ***************************************************************************/
  95 +ctulimit::ctulimit(int argc, char *argv[]) :
  96 + ctool(CTULIMIT_NAME, CTULIMIT_VERSION, argc, argv)
  97 +{
  98 + // Initialise members
  99 + init_members();
  100 +
  101 + // Return
  102 + return;
  103 +}
  104 +
  105 +
  106 +/***********************************************************************//**
  107 + * @brief Copy constructor
  108 + *
  109 + * @param[in] app Application.
  110 + ***************************************************************************/
  111 +ctulimit::ctulimit(const ctulimit& app) : ctool(app)
  112 +{
  113 + // Initialise members
  114 + init_members();
  115 +
  116 + // Copy members
  117 + copy_members(app);
  118 +
  119 + // Return
  120 + return;
  121 +}
  122 +
  123 +
  124 +/***********************************************************************//**
  125 + * @brief Destructor
  126 + ***************************************************************************/
  127 +ctulimit::~ctulimit(void)
  128 +{
  129 + // Free members
  130 + free_members();
  131 +
  132 + // Return
  133 + return;
  134 +}
  135 +
  136 +
  137 +/*==========================================================================
  138 + = =
  139 + = Operators =
  140 + = =
  141 + ==========================================================================*/
  142 +
  143 +/***********************************************************************//**
  144 + * @brief Assignment operator
  145 + *
  146 + * @param[in] app Application.
  147 + * @return Application.
  148 + ***************************************************************************/
  149 +ctulimit& ctulimit::operator=(const ctulimit& app)
  150 +{
  151 + // Execute only if object is not identical
  152 + if (this != &app) {
  153 +
  154 + // Copy base class members
  155 + this->ctool::operator=(app);
  156 +
  157 + // Free members
  158 + free_members();
  159 +
  160 + // Initialise members
  161 + init_members();
  162 +
  163 + // Copy members
  164 + copy_members(app);
  165 +
  166 + } // endif: object was not identical
  167 +
  168 + // Return this object
  169 + return *this;
  170 +}
  171 +
  172 +
  173 +/*==========================================================================
  174 + = =
  175 + = Public methods =
  176 + = =
  177 + ==========================================================================*/
  178 +
  179 +/***********************************************************************//**
  180 + * @brief Clear instance
  181 + ***************************************************************************/
  182 +void ctulimit::clear(void)
  183 +{
  184 + // Free members
  185 + free_members();
  186 + this->ctool::free_members();
  187 + this->GApplication::free_members();
  188 +
  189 + // Initialise members
  190 + this->GApplication::init_members();
  191 + this->ctool::init_members();
  192 + init_members();
  193 +
  194 + // Return
  195 + return;
  196 +}
  197 +
  198 +
  199 +/***********************************************************************//**
  200 + * @brief Computes the upper limit
  201 + *
  202 + * This method calculates the upper limit depending on the given algorithm
  203 + ***************************************************************************/
  204 +void ctulimit::run(void)
  205 +{
  206 + // If we're in debug mode then all output is also dumped on the screen
  207 + if (logDebug()) {
  208 + log.cout(true);
  209 + }
  210 +
  211 + // Get task parameters
  212 + get_parameters();
  213 +
  214 + // Write parameters into logger
  215 + if (logTerse()) {
  216 + log_parameters();
  217 + log << std::endl;
  218 + }
  219 +
  220 + // Write observation(s) into logger
  221 + if (logTerse()) {
  222 + log << std::endl;
  223 + if (m_obs.size() > 1) {
  224 + log.header1("Observations");
  225 + }
  226 + else {
  227 + log.header1("Observation");
  228 + }
  229 + log << m_obs << std::endl;
  230 + }
  231 +
  232 + // Save original models
  233 + GModels models_orig = m_obs.models();
  234 +
  235 + // Save original log-likelihood. If the value is zero it has never been
  236 + // computed hence we compute it now.
  237 + m_best_logL = m_obs.logL();
  238 + if (m_best_logL == 0.0) {
  239 +
  240 + // Write header
  241 + if (logTerse()) {
  242 + log << std::endl;
  243 + log.header1("Compute best-fit likelihood");
  244 + }
  245 +
  246 + // Reoptimize if likelihood was not given before
  247 + GOptimizerLM* opt = new GOptimizerLM();
  248 + m_obs.optimize(*opt);
  249 + m_obs.errors(*opt);
  250 + m_best_logL = m_obs.logL();
  251 +
  252 + // Write optimised model into logger
  253 + if (logTerse()) {
  254 + log << m_obs.models() << std::endl;
  255 + }
  256 +
  257 + } // endif: likelihood was zero
  258 +
  259 + // Compute parameter bracketing
  260 + double value = m_model_par->value();
  261 + double error = m_model_par->error();
  262 + double parmin = value + m_sigma_min * error;
  263 + double parmax = value + m_sigma_max * error;
  264 +
  265 + // Write header
  266 + if (logTerse()) {
  267 + log << std::endl;
  268 + log.header1("Compute upper limit");
  269 + log << gammalib::parformat("Model name");
  270 + log << m_skymodel->name() << std::endl;
  271 + log << gammalib::parformat("Parameter name");
  272 + log << m_model_par->name() << std::endl;
  273 + log << gammalib::parformat("Initial parameter range");
  274 + log << "[";
  275 + log << parmin;
  276 + log << ", ";
  277 + log << parmax;
  278 + log << "]" << std::endl;
  279 + }
  280 +
  281 + // Compute upper limit
  282 + ulimit_bisection(parmin, parmax);
  283 +
  284 + // Write final parameter range
  285 + if (logTerse()) {
  286 + log << gammalib::parformat("Final parameter");
  287 + log << m_model_par->value() << std::endl;
  288 + }
  289 +
  290 + // Get reference energy for differential upper limit
  291 + GEnergy eref = GEnergy(m_eref, "TeV");
  292 +
  293 + // Create energy range for flux limits
  294 + GEnergy emin = GEnergy(m_emin, "TeV");
  295 + GEnergy emax = GEnergy(m_emax, "TeV");
  296 +
  297 + // Compute upper limit intensity and fluxes
  298 + m_diff_ulimit = m_skymodel->spectral()->eval(eref, GTime());
  299 + m_flux_ulimit = m_skymodel->spectral()->flux(emin, emax);
  300 + m_eflux_ulimit = m_skymodel->spectral()->eflux(emin, emax);
  301 +
  302 + // Write results to logfile
  303 + if (logTerse()) {
  304 + log << std::endl;
  305 + log.header1("Upper limit results");
  306 + log << gammalib::parformat("Differential flux limit");
  307 + log << m_diff_ulimit;
  308 + log << " ph/cm2/s/MeV at ";
  309 + log << m_eref << " TeV";
  310 + log << std::endl;
  311 + log << gammalib::parformat("Integral flux limit");
  312 + log << m_flux_ulimit;
  313 + log << " ph/cm2/s within [";
  314 + log << m_emin << "-" << m_emax;
  315 + log << "] TeV";
  316 + log << std::endl;
  317 + log << gammalib::parformat("Energy flux limit");
  318 + log << m_eflux_ulimit;
  319 + log << " erg/cm2/s within [";
  320 + log << m_emin << "-" << m_emax;
  321 + log << "] TeV";
  322 + log << std::endl;
  323 + }
  324 +
  325 + // Recover original models
  326 + m_obs.models(models_orig);
  327 +
  328 + // Return
  329 + return;
  330 +}
  331 +
  332 +
  333 +/***********************************************************************//**
  334 + * @brief Save maps
  335 + *
  336 + * This method saves the upper limit to an ascii file
  337 + ***************************************************************************/
  338 +void ctulimit::save(void)
  339 +{
  340 + // Write header
  341 + if (logTerse()) {
  342 + log << std::endl;
  343 + log.header1("Save upper limit");
  344 + }
  345 +
  346 + // Get output filename
  347 + m_outfile = (*this)["outfile"].filename();
  348 +
  349 + // Create CSV table with 3 columns
  350 + GCsv table(1, 3);
  351 +
  352 + table.real(0, 0, m_diff_ulimit);
  353 + table.real(0, 1, m_flux_ulimit);
  354 + table.real(0, 2, m_eflux_ulimit);
  355 +
  356 + // Save CSV table
  357 + table.save(m_outfile, " ", clobber());
  358 +
  359 + // Return
  360 + return;
  361 +}
  362 +
  363 +
  364 +/*==========================================================================
  365 + = =
  366 + = Private methods =
  367 + = =
  368 + ==========================================================================*/
  369 +
  370 +/***********************************************************************//**
  371 + * @brief Initialise class members
  372 + ***************************************************************************/
  373 +void ctulimit::init_members(void)
  374 +{
  375 + // Initialise user parameters
  376 + m_outfile.clear();
  377 + m_srcname.clear();
  378 + m_sigma_min = 0.0;
  379 + m_sigma_max = 0.0;
  380 + m_eref = 0.0;
  381 + m_emin = 0.0;
  382 + m_emax = 0.0;
  383 + m_tol = 1.0e-6;
  384 + m_max_iter = 50;
  385 +
  386 + // Initialise protected members
  387 + m_obs.clear();
  388 + m_models.clear();
  389 + m_dlogL = 0.0;
  390 + m_skymodel = NULL;
  391 + m_model_par = NULL;
  392 + m_best_logL = 0.0;
  393 + m_flux_ulimit = 0.0;
  394 + m_diff_ulimit = 0.0;
  395 + m_eflux_ulimit = 0.0;
  396 +
  397 + // Return
  398 + return;
  399 +}
  400 +
  401 +
  402 +/***********************************************************************//**
  403 + * @brief Copy class members
  404 + *
  405 + * @param[in] app Application.
  406 + ***************************************************************************/
  407 +void ctulimit::copy_members(const ctulimit& app)
  408 +{
  409 + // Copy user parameters
  410 + m_outfile = app.m_outfile;
  411 + m_srcname = app.m_srcname;
  412 + m_sigma_min = app.m_sigma_min;
  413 + m_sigma_max = app.m_sigma_max;
  414 + m_eref = app.m_eref;
  415 + m_emin = app.m_emin;
  416 + m_emax = app.m_emax;
  417 + m_tol = app.m_tol;
  418 + m_max_iter = app.m_max_iter;
  419 +
  420 +
  421 + // Copy protected members
  422 + m_obs = app.m_obs;
  423 + m_models = app.m_models;
  424 + m_dlogL = app.m_dlogL;
  425 + m_best_logL = app.m_best_logL;
  426 + m_diff_ulimit = app.m_diff_ulimit;
  427 + m_flux_ulimit = app.m_flux_ulimit;
  428 + m_eflux_ulimit = app.m_eflux_ulimit;
  429 +
  430 + // Extract model parameter
  431 + get_model_parameter();
  432 +
  433 + // Return
  434 + return;
  435 +}
  436 +
  437 +
  438 +/***********************************************************************//**
  439 + * @brief Delete class members
  440 + ***************************************************************************/
  441 +void ctulimit::free_members(void)
  442 +{
  443 + // Return
  444 + return;
  445 +}
  446 +
  447 +
  448 +/***********************************************************************//**
  449 + * @brief Get application parameters
  450 + *
  451 + * @exception GException::invalid_value
  452 + * Test source not found.
  453 + *
  454 + * Get all task parameters from parameter file or (if required) by querying
  455 + * the user.
  456 + ***************************************************************************/
  457 +void ctulimit::get_parameters(void)
  458 +{
  459 + // If there are no observations in container then load them via user
  460 + // parameters
  461 + if (m_obs.size() == 0) {
  462 +
  463 + // Throw exception if no input observation file is given
  464 + require_inobs(G_GET_PARAMETERS);
  465 +
  466 + // Build observation container
  467 + m_obs = get_observations();
  468 +
  469 + } // endif: there was no observation in the container
  470 +
  471 +
  472 + // If there is are no models associated with the observations then
  473 + // load now the model definition
  474 + if (m_obs.models().size() == 0) {
  475 +
  476 + // Get models XML filename
  477 + std::string filename = (*this)["inmodel"].filename();
  478 +
  479 + // Setup models for optimizing.
  480 + m_obs.models(GModels(filename));
  481 +
  482 + } // endif: no models were associated with observations
  483 +
  484 + // Get name of test source and check container for this name
  485 + m_srcname = (*this)["srcname"].string();
  486 + if (!m_obs.models().contains(m_srcname)) {
  487 + std::string msg = "Source \""+m_srcname+"\" not found in model "
  488 + "container. Please add a source with that name "
  489 + "or check for possible typos.";
  490 + throw GException::invalid_value(G_GET_PARAMETERS, msg);
  491 + }
  492 +
  493 + // Get relevant model and parameter for upper limit computation
  494 + get_model_parameter();
  495 +
  496 + // Get confidence level
  497 + double confidence = (*this)["confidence"].real();
  498 + if (confidence != 0.95) {
  499 + std::string msg = "Confidence level different from 95% requested."
  500 + " Currently only 95% is supported.";
  501 + throw GException::invalid_value(G_GET_PARAMETERS, msg);
  502 + }
  503 + else {
  504 + // Set Likelihood difference for 95% CL.
  505 + // See Minuit Handbook
  506 + m_dlogL = 3.84 / 2.0;
  507 + }
  508 +
  509 + // Read starting boundaries for bisection
  510 + m_sigma_min = (*this)["sigma_min"].real();
  511 + m_sigma_max = (*this)["sigma_max"].real();
  512 +
  513 + // Read energy values
  514 + m_eref = (*this)["eref"].real();
  515 + m_emin = (*this)["emin"].real();
  516 + m_emax = (*this)["emax"].real();
  517 +
  518 + // Read precision
  519 + m_tol = (*this)["tol"].real();
  520 + m_max_iter = (*this)["max_iter"].integer();
  521 +
  522 + // Optionally read ahead parameters so that they get correctly
  523 + // dumped into the log file
  524 + if (read_ahead()) {
  525 + m_outfile = (*this)["outfile"].filename();
  526 + }
  527 +
  528 + // Return
  529 + return;
  530 +}
  531 +
  532 +
  533 +/***********************************************************************//**
  534 + * @brief Get application parameters
  535 + *
  536 + * @exception GException::invalid_value
  537 + * Did not find a valid model parameter
  538 + *
  539 + * Extracts a pointer to the sky model (m_skymodel) and a pointer to the
  540 + * relevant model parameter (m_model_par) from the model container.
  541 + ***************************************************************************/
  542 +void ctulimit::get_model_parameter(void)
  543 +{
  544 + // Get relevant model and parameter for upper limit computation.
  545 + GModels& models = const_cast<GModels&>(m_obs.models());
  546 + m_skymodel = dynamic_cast<GModelSky*>(models[m_srcname]);
  547 + if (m_skymodel == NULL) {
  548 + std::string msg = "Source \""+m_srcname+"\" is not a sky model. "
  549 + "Please specify the name of a sky model for "
  550 + "upper limit computation.";
  551 + throw GException::invalid_value(G_GET_MODEL_PARAMETER, msg);
  552 + }
  553 + if (m_skymodel->spectral()->has_par("Normalization")) {
  554 + m_model_par = &(m_skymodel->spectral()->operator[]("Normalization"));
  555 + }
  556 + else if (m_skymodel->spectral()->has_par("Prefactor")) {
  557 + m_model_par = &(m_skymodel->spectral()->operator[]("Prefactor"));
  558 + }
  559 + else if (m_skymodel->spectral()->has_par("Integral")) {
  560 + m_model_par = &(m_skymodel->spectral()->operator[]("Integral"));
  561 + }
  562 + else {
  563 + std::string msg = "Require spectral parameter \"Normalization\", "
  564 + "\"Prefactor\" or \"Integral\" for upper limit "
  565 + "computation. The specified source \""+m_srcname+
  566 + "\" does not have such a parameter.";
  567 + throw GException::invalid_value(G_GET_MODEL_PARAMETER, msg);
  568 + }
  569 +
  570 + // Return
  571 + return;
  572 +}
  573 +
  574 +
  575 +/***********************************************************************//**
  576 + * @brief Performs upper limit computation by using a bisection method
  577 + *
  578 + * @param[in] min Minimum parameter value
  579 + * @param[in] max Maximum parameter value
  580 + *
  581 + * This method calculates the upper limit using a bisection method.
  582 + ***************************************************************************/
  583 +void ctulimit::ulimit_bisection(const double& min, const double& max)
  584 +{
  585 + // Copy values to working values
  586 + double wrk_min = min;
  587 + double wrk_max = max;
  588 +
  589 + // Initialise iteration counter
  590 + int iter = 0;
  591 +
  592 + // Loop until breaking condition is reached
  593 + while (true) {
  594 +
  595 + // Log information
  596 + if (logExplicit()) {
  597 + log << gammalib::parformat("Iteration "+gammalib::str(iter));
  598 + log << "[";
  599 + log << wrk_min;
  600 + log << ", ";
  601 + log << wrk_max;
  602 + log << "]" << std::endl;
  603 + }
  604 +
  605 + // Throw exception if maximum iterations are reached
  606 + if (iter > m_max_iter) {
  607 + std::string msg = "The maximum number of "+gammalib::str(m_max_iter)+
  608 + " has been reached. You may consider to increase"
  609 + " the \"max_iter\" parameter and re-run ctulimit.";
  610 + throw GException::invalid_value(G_UL_BISECTION, msg);
  611 + }
  612 +
  613 + // Compute center of boundary
  614 + double mid = (wrk_min + wrk_max) / 2.0;
  615 +
  616 + // Calculate function value
  617 + double eval_mid = evaluate(mid);
  618 +
  619 + // Check for convergence inside tolerance
  620 + if (std::abs(eval_mid) < m_tol) {
  621 + break;
  622 + }
  623 +
  624 + // Change boundaries for further iteration
  625 + if (eval_mid > 0.0) {
  626 + wrk_max = mid;
  627 + }
  628 + else if (eval_mid < 0.0) {
  629 + wrk_min = mid;
  630 + }
  631 +
  632 + // Increment counter
  633 + iter++;
  634 +
  635 + } // endwhile
  636 +
  637 + // Return
  638 + return;
  639 +
  640 +}
  641 +
  642 +
  643 +/***********************************************************************//**
  644 + * @brief Evaluates the log-likelihood
  645 + *
  646 + * @param[in] value Parameter factor value
  647 + * @return Log-likelihood value
  648 + *
  649 + * This method evaluates the log-likelihood as a function of the parameter
  650 + * of interest.
  651 + ***************************************************************************/
  652 +double ctulimit::evaluate(const double& value)
  653 +{
  654 + // Initialise log-likelihood value
  655 + double logL = 0.0;
  656 +
  657 + // Check if given parameter is within boundaries
  658 + if (value > m_model_par->min() && value < m_model_par->max()) {
  659 +
  660 + // Change parameter factor
  661 + m_model_par->value(value);
  662 +
  663 + // Evaluate likelihood for new model container
  664 + m_obs.eval();
  665 +
  666 + // Retrieve likelihood
  667 + logL = m_obs.logL();
  668 +
  669 + } // endif: value was inside allowed range
  670 +
  671 + // ... otherwise signal that the parameter went outside the boundaries
  672 + else {
  673 + std::string msg = "Value of parameter \""+m_model_par->name()+"\" "
  674 + "outside of validity range requested. To omit "
  675 + "this error please enlarge the parameter range "
  676 + "in the model XML file.";
  677 + throw GException::invalid_value(G_EVALUATE, msg);
  678 + }
  679 +
  680 + // Compute function value
  681 + double logL_difference = logL - m_best_logL - m_dlogL;
  682 +
  683 + // Return
  684 + return logL_difference;
  685 +}
... ...
src/ctulimit/ctulimit.hpp 0 → 100644
  View file @e511b18
  1 +/***************************************************************************
  2 + * ctulimit - Upper limit calculation tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2015 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 "01-00-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 source.
  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
  51 + * best value for the upper limit.
  52 + ***************************************************************************/
  53 +class ctulimit : public ctool {
  54 +
  55 +public:
  56 + // Constructors and destructors
  57 + ctulimit(void);
  58 + explicit ctulimit(const GObservations& obs);
  59 + ctulimit(int argc, char *argv[]);
  60 + ctulimit(const ctulimit& app);
  61 + virtual ~ctulimit(void);
  62 +
  63 + // Operators
  64 + ctulimit& operator=(const ctulimit& app);
  65 +
  66 + // Methods
  67 + void clear(void);
  68 + void run(void);
  69 + void save(void);
  70 + const GObservations& obs(void) const;
  71 + const double& diff_ulimit(void) const;
  72 + const double& flux_ulimit(void) const;
  73 + const double& eflux_ulimit(void) const;
  74 +
  75 +protected:
  76 + // Protected methods
  77 + void init_members(void);
  78 + void copy_members(const ctulimit& app);
  79 + void free_members(void);
  80 + void get_parameters(void);
  81 + void get_model_parameter(void);
  82 + void ulimit_bisection(const double& parmin, const double& parmax);
  83 + double evaluate(const double& value);
  84 +
  85 + // User parameters
  86 + std::string m_srcname; //!< Name of source which is moved around
  87 + std::string m_outfile; //!< Output ascii file
  88 + double m_sigma_min; //!< Starting value minimum (multiple fit errors above fit values)
  89 + double m_sigma_max; //!< Starting value maximum (multiple fit errors above fit values)
  90 + double m_eref; //!< Reference energy for flux limits (TeV)
  91 + double m_emin; //!< Minimum energy for flux limits (TeV)
  92 + double m_emax; //!< Maximum energy for flux limits (TeV)
  93 + double m_tol; //!< Tolerance for limit determination
  94 + int m_max_iter; //!< Maximum number of iterations
  95 +
  96 + // Protected members
  97 + GObservations m_obs; //!< Observation container
  98 + GModels m_models;
  99 + double m_dlogL; //!< Likelihood difference for upper limit computation
  100 + GModelSky* m_skymodel; //!< Pointer to sky model
  101 + GModelPar* m_model_par; //!< Pointer to model parameter
  102 + double m_best_logL; //!< Best fit log likelihood of given model
  103 + double m_diff_ulimit; //!< Differential upper limit value
  104 + double m_flux_ulimit; //!< Flux upper limit value
  105 + double m_eflux_ulimit; //!< Energy flux upper limits
  106 +};
  107 +
  108 +
  109 +/***********************************************************************//**
  110 + * @brief Return observation container
  111 + *
  112 + * @return Reference to observation container
  113 + ***************************************************************************/
  114 +inline
  115 +const GObservations& ctulimit::obs(void) const
  116 +{
  117 + return m_obs;
  118 +}
  119 +
  120 +/***********************************************************************//**
  121 + * @brief Return differential upper limit
  122 + *
  123 + * @return differential upper limit
  124 + ***************************************************************************/
  125 +inline
  126 +const double& ctulimit::diff_ulimit(void) const
  127 +{
  128 + return m_diff_ulimit;
  129 +}
  130 +
  131 +/***********************************************************************//**
  132 + * @brief Return flux upper limit
  133 + *
  134 + * @return flux upper limit
  135 + ***************************************************************************/
  136 +inline
  137 +const double& ctulimit::flux_ulimit(void) const
  138 +{
  139 + return m_flux_ulimit;
  140 +}
  141 +
  142 +/***********************************************************************//**
  143 + * @brief return energy flux upper limit
  144 + *
  145 + * @return energy flux upper limits
  146 + ***************************************************************************/
  147 +inline
  148 +const double& ctulimit::eflux_ulimit(void) const
  149 +{
  150 + return m_eflux_ulimit;
  151 +}
  152 +
  153 +#endif /* CTULIMIT_HPP */
... ...
src/ctulimit/ctulimit.par 0 → 100644
  View file @e511b18
  1 +#############################################################################
  2 +# #
  3 +# CTA Science Analysis Tools #
  4 +# #
  5 +# ------------------------------------------------------------------------- #
  6 +# #
  7 +# File: ctulimit.par #
  8 +# Version: 01-00-00 #
  9 +# #
  10 +# Author: Michael Mayer #
  11 +# michael.mayer@physik.hu-berlin.de #
  12 +# HU Berlin #
  13 +# #
  14 +# Purpose: Parameter file for the ctulimit application #
  15 +# #
  16 +#############################################################################
  17 +#
  18 +# General parameters
  19 +#===================
  20 +inobs, f, a, "events.fits",,, "Input event list, counts cube or observation definition file"
  21 +inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model"
  22 +srcname, s, a, "Crab",,, "Source of interest"
  23 +expcube, f, a, "NONE",,, "Exposure cube file (only needed for stacked analysis)"
  24 +psfcube, f, a, "NONE",,, "PSF cube file (only needed for stacked analysis)"
  25 +caldb, s, a, "dummy",,, "Calibration database"
  26 +irf, s, a, "cta_dummy_irf",,, "Instrument response function"
  27 +outfile, f, a, "upperlimit.txt",,, "Output ASCII file"
  28 +
  29 +#
  30 +# Upper limit calculation parameters
  31 +#===================================
  32 +confidence, r, h, 0.95,0.0,1.0, "Confidence level"
  33 +sigma_min, r, h, 0.0,,, "Start minimum value (multiple sigmas above fit value)"
  34 +sigma_max, r, h, 10.0,,, "Start maximum value (multiple sigmas above fit value)"
  35 +eref, r, h, 1.0,,, "Reference energy for differential limits (TeV)"
  36 +emin, r, h, 1.0,,, "Minimum energy for flux limits (TeV)"
  37 +emax, r, h, 100.0,,, "Maximum energy for flux limits (TeV)"
  38 +tol, r, h, 1e-5,,, "Computation tolerance"
  39 +max_iter, i, h, 50,,, "Maximum number of iterations"
  40 +
  41 +#
  42 +# Standard parameters
  43 +#====================
  44 +chatter, i, h, 2,0,4, "Chattiness of output"
  45 +clobber, b, h, yes,,, "Overwrite existing output files with new output files?"
  46 +debug, b, h, no,,, "Debugging mode activated"
  47 +mode, s, h, "ql",,, "Mode of automatic parameters"
  48 +logfile, f, h, "ctulimit.log",,, "Log filename"
... ...
src/ctulimit/main.cpp 0 → 100644
  View file @e511b18
  1 +/***************************************************************************
  2 + * ctulimit - upper limit calculation tool *
  3 + * ----------------------------------------------------------------------- *
  4 + * copyright (C) 2015 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 + " ctulimit. 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 +}
... ...
test/test_ctools.sh
  View file @e511b18
... ... @@ -46,12 +46,14 @@ ctpsfcube=../src/ctpsfcube/ctpsfcube
46 46 ctbkgcube=../src/ctbkgcube/ctbkgcube
47 47 ctcubemask=../src/ctcubemask/ctcubemask
48 48 ctbutterfly=../src/ctbutterfly/ctbutterfly
  49 +ctulimit=../src/ctulimit/ctulimit
49 50  
50 51  
51 52 #
52 53 # Remove any existing result files
53 54 # ================================
54 55 rm -rf *.fits *.xml
  56 +rm -rf ulimit.dat butterfly.txt
55 57  
56 58  
57 59 #
... ... @@ -597,3 +599,25 @@ else
597 599 exit 1
598 600 fi
599 601 $ECHO " ok"
  602 +
  603 +
  604 +#
  605 +# Test ctulimit
  606 +# ===============
  607 +$ECHO -n "Test ctulimit: "
  608 +$ctulimit inobs="data/crab_events.fits.gz" \
  609 + inmodel="data/crab.xml" \
  610 + srcname="Crab" \
  611 + outfile="ulimit.dat" \
  612 + caldb="irf" \
  613 + irf="cta_dummy_irf"
  614 +$ECHO -n "."
  615 +if [ -s "ulimit.dat" ]
  616 +then
  617 + $ECHO -n "."
  618 +else
  619 + $ECHO " ulimit.dat file is not found"
  620 + exit 1
  621 +fi
  622 +$ECHO " ok"
  623 +
... ...
test/test_ctulimit.py 0 → 100755
  View file @e511b18
  1 +#! /usr/bin/env python
  2 +# ==========================================================================
  3 +# This scripts performs unit tests for the ctulimit tool.
  4 +#
  5 +# Copyright (C) 2015 Michael Mayer
  6 +#
  7 +# This program is free software: you can redistribute it and/or modify
  8 +# it under the terms of the GNU General Public License as published by
  9 +# the Free Software Foundation, either version 3 of the License, or
  10 +# (at your option) any later version.
  11 +#
  12 +# This program is distributed in the hope that it will be useful,
  13 +# but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  15 +# GNU General Public License for more details.
  16 +#
  17 +# You should have received a copy of the GNU General Public License
  18 +# along with this program. If not, see <http://www.gnu.org/licenses/>.
  19 +#
  20 +# ==========================================================================
  21 +import gammalib
  22 +import ctools
  23 +
  24 +
  25 +# ========================== #
  26 +# Test class for cttsmap tool #
  27 +# ========================== #
  28 +class Test(gammalib.GPythonTestSuite):
  29 + """
  30 + Test class for cttsmap tool.
  31 + """
  32 + # Constructor
  33 + def __init__(self):
  34 + """
  35 + Constructor.
  36 + """
  37 + # Call base class constructor
  38 + gammalib.GPythonTestSuite.__init__(self)
  39 +
  40 + # Set members
  41 + self.events_name = "data/crab_events.fits"
  42 + self.model_name = "data/crab.xml"
  43 + self.caldb = "irf"
  44 + self.irf = "cta_dummy_irf"
  45 +
  46 + # Return
  47 + return
  48 +
  49 + # Set test functions
  50 + def set(self):
  51 + """
  52 + Set all test functions.
  53 + """
  54 + # Set test name
  55 + self.name("ctulimit")
  56 +
  57 + # Append tests
  58 + self.append(self.test_functional, "Test ctulimit functionality")
  59 +
  60 + # Return
  61 + return
  62 +
  63 + # Test cttsmap functionnality
  64 + def test_functional(self):
  65 + """
  66 + Test cttsmap functionnality.
  67 + """
  68 + # Set-up cttsmap
  69 + ulimit = ctools.ctulimit()
  70 + ulimit["inobs"].filename(self.events_name)
  71 + ulimit["inmodel"].filename(self.model_name)
  72 + ulimit["srcname"].string("Crab")
  73 + ulimit["caldb"].string(self.caldb)
  74 + ulimit["irf"].string(self.irf)
  75 + ulimit["outfile"].filename("ulimit.dat")
  76 +
  77 + # Run tool
  78 + self.test_try("Run ctulimit")
  79 + try:
  80 + ulimit.run()
  81 + self.test_try_success()
  82 + except:
  83 + self.test_try_failure("Exception occured in ctulimit.")
  84 +
  85 + # Save results
  86 + self.test_try("Save results")
  87 + try:
  88 + ulimit.save()
  89 + self.test_try_success()
  90 + except:
  91 + self.test_try_failure("Exception occured in saving results.")
  92 +
  93 + # Return
  94 + return
... ...
test/test_python.py
  View file @e511b18
... ... @@ -34,6 +34,7 @@ import test_ctpsfcube
34 34 import test_ctbkgcube
35 35 import test_ctcubemask
36 36 import test_ctbutterfly
  37 +import test_ctulimit
37 38 import test_pipelines
38 39  
39 40  
... ... @@ -60,6 +61,7 @@ if __name__ == &#39;__main__&#39;:
60 61 suite_ctbkgcube = test_ctbkgcube.Test()
61 62 suite_ctcubemask = test_ctcubemask.Test()
62 63 suite_ctbutterfly = test_ctbutterfly.Test()
  64 + suite_ctulimit = test_ctulimit.Test()
63 65 suite_pipelines = test_pipelines.Test()
64 66  
65 67 # Setup unit tests
... ... @@ -75,6 +77,7 @@ if __name__ == &#39;__main__&#39;:
75 77 suite_ctbkgcube.set()
76 78 suite_ctcubemask.set()
77 79 suite_ctbutterfly.set()
  80 + suite_ctulimit.set()
78 81 suite_pipelines.set()
79 82  
80 83 # Append tests to container
... ... @@ -90,6 +93,7 @@ if __name__ == &#39;__main__&#39;:
90 93 suites.append(suite_ctbkgcube)
91 94 suites.append(suite_ctcubemask)
92 95 suites.append(suite_ctbutterfly)
  96 + suites.append(suite_ctulimit)
93 97 suites.append(suite_pipelines)
94 98  
95 99 # Set PFILES environment variable
... ...