diff --git a/configure.ac b/configure.ac index 8c174ef..b4132dd 100644 --- a/configure.ac +++ b/configure.ac @@ -888,6 +888,7 @@ AC_CONFIG_FILES([Makefile src/ctcubemask/Makefile src/cttsmap/Makefile src/ctbutterfly/Makefile + src/ctulimit/Makefile src/support/Makefile src/ctools-setup test/Makefile diff --git a/doc/source/reference_manual/ctulimit.rst b/doc/source/reference_manual/ctulimit.rst new file mode 100644 index 0000000..170dc03 --- /dev/null +++ b/doc/source/reference_manual/ctulimit.rst @@ -0,0 +1,99 @@ +.. _ctulimit: + +ctulimit +=========== + +Computes upper limit for a specific sky model component. + + +Synopsis +-------- + +Computes the upper limit flux of a given source and stores the value +in an ASCII file. + + +General parameters +------------------ + +``inobs [file]`` + Input event list, counts cube or observation definition XML file. + +``inmodel [file]`` + Model XML file containing the source and background definitions. + +``srcname [string]`` + Name of model component for which upper limit should be computed. + +``expcube = NONE [file]`` + Exposure cube file (only needed for stacked analysis). + +``psfcube = NONE [file]`` + PSF cube file (only needed for stacked analysis). + +``caldb [string]`` + Calibration database. + +``irf [string]`` + Instrumental response function. + +``outfile [file]`` + Output ASCII file name. + +``(confidence = 0.95) [real]`` + Confidence Level of upper limit. + +``(sigma_min = 0.0) [real]`` + Minimum boundary to start searching for upper limit value. + Number of standard deviations above best fit value + +``(sigma_max = 10.0) [real]`` + Maximum boundary to start searching for upper limit value. + Number of standard deviations above best fit value + +``(emin = 1.0) [real])`` + Minimum energy of flux limits (in TeV). + +``(emax = 100) [real]`` + Maximum energy of flux limits (in TeV). + +``(tol = 1e-5) [real]`` + Computation tolerance. + +``(max_iter = 50) [integer]`` + Maximum number of iterations before stopping the upper + limit computations. + + +Standard parameters +------------------- + +``(chatter = 2) [integer]`` + Verbosity of the executable: + chatter = 0: no information will be logged + + chatter = 1: only errors will be logged + + chatter = 2: errors and actions will be logged + + chatter = 3: report about the task execution + + chatter = 4: detailed report about the task execution + +``(clobber = yes) [boolean]`` + Specifies whether an existing output counts cube should be overwritten. + +``(debug = no) [boolean]`` + Enables debug mode. In debug mode the executable will dump any log file output to the console. + +``(mode = ql) [string]`` + Mode of automatic parameters (default is "ql", i.e. "query and learn"). + +``(logfile = ctulimit.log) [string]`` + Name of log file. + + +Related tools +------------- + +None diff --git a/doc/source/reference_manual/index.rst b/doc/source/reference_manual/index.rst index 4495dee..eccacb8 100644 --- a/doc/source/reference_manual/index.rst +++ b/doc/source/reference_manual/index.rst @@ -27,6 +27,7 @@ ctools ctselect --- Selects event data ctskymap --- Generates sky map cttsmap --- Generates Test Statistics map + ctulimit --- Calculates upper limit cscripts diff --git a/pyext/ctools/tools.i b/pyext/ctools/tools.i index 7d75cf1..381daa8 100755 --- a/pyext/ctools/tools.i +++ b/pyext/ctools/tools.i @@ -58,3 +58,4 @@ %include "ctcubemask.i" %include "cttsmap.i" %include "ctbutterfly.i" +%include "ctulimit.i" diff --git a/pyext/ctulimit.i b/pyext/ctulimit.i new file mode 100644 index 0000000..e722ff8 --- /dev/null +++ b/pyext/ctulimit.i @@ -0,0 +1,64 @@ +/*************************************************************************** + * ctulimit - Upper limit calculation tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2015 by Michael Mayer * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file ctulimit + * @brief Upper limit calculation tool + * @author Michael Mayer + */ +%{ +/* Put headers and other declarations here that are needed for compilation */ +#include "ctulimit.hpp" +%} + + +/***********************************************************************//** + * @class ctulimit + * + * @brief Upper limit calculation tool + ***************************************************************************/ +class ctulimit : public ctool { +public: + // Constructors and destructors + ctulimit(void); + explicit ctulimit(const GObservations& obs); + ctulimit(int argc, char *argv[]); + ctulimit(const ctulimit& app); + virtual ~ctulimit(void); + + // Methods + void clear(void); + void run(void); + void save(void); + const GObservations& obs(void) const; + const double& diff_ulimit(void) const; + const double& flux_ulimit(void) const; + const double& eflux_ulimit(void) const; +}; + + +/***********************************************************************//** + * @brief Upper limit calculation tool Python extensions + ***************************************************************************/ +%extend ctulimit { + ctulimit copy() { + return (*self); + } +} diff --git a/src/Makefile.am b/src/Makefile.am index bb0cddf..76f52a1 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -17,7 +17,8 @@ SUBDIRS = support \ ctbkgcube \ ctcubemask \ cttsmap \ - ctbutterfly + ctbutterfly \ + ctulimit # Scripts to be distributed scripts = ctools-setup \ @@ -41,7 +42,8 @@ libctools_la_LIBADD = support/libsupport.la \ ctbkgcube/libctbkgcube.la \ ctcubemask/libctcubemask.la \ cttsmap/libcttsmap.la \ - ctbutterfly/libctbutterfly.la + ctbutterfly/libctbutterfly.la \ + ctulimit/libctulimit.la # Distribute scripts diff --git a/src/ctulimit/Makefile.am b/src/ctulimit/Makefile.am new file mode 100644 index 0000000..d5d844a --- /dev/null +++ b/src/ctulimit/Makefile.am @@ -0,0 +1,24 @@ +## Process this file with automake to produce Makefile.in + +# Set include directories +INCLUDES = -I$(top_srcdir)/src/support + +# Program name +bin_PROGRAMS = ctulimit + +# Sources +sources = ctulimit.cpp + +# Libtool library +noinst_LTLIBRARIES = libctulimit.la +libctulimit_la_SOURCES = $(sources) + +# Program definition +ctulimit_SOURCES = main.cpp +ctulimit_LDFLAGS = $(LIB_LINK_OPTION) @LDFLAGS@ +ctulimit_CPPFLAGS = @CPPFLAGS@ +ctulimit_LDADD = libctulimit.la $(top_srcdir)/src/support/libsupport.la + +# Parameter files +pfilesdir = $(prefix)/syspfiles +pfiles_DATA = ctulimit.par diff --git a/src/ctulimit/ctulimit.cpp b/src/ctulimit/ctulimit.cpp new file mode 100644 index 0000000..743344c --- /dev/null +++ b/src/ctulimit/ctulimit.cpp @@ -0,0 +1,685 @@ +/*************************************************************************** + * ctulimit - Upper limit calculation tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2015 by Michael Mayer * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file ctulimit.hpp + * @brief Upper limit calculation tool interface implementation + * @author Michael Mayer + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#include "ctulimit.hpp" +#include "GTools.hpp" +#include "GOptimizer.hpp" + +/* __ Method name definitions ____________________________________________ */ +#define G_GET_PARAMETERS "ctulimit::get_parameters()" +#define G_GET_MODEL_PARAMETER "ctulimit::get_model_parameter()" +#define G_UL_BISECTION "ctulimit::ul_bisection(double&, double&)" +#define G_EVALUATE "ctulimit::evaluate(double&)" + +/* __ Debug definitions __________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + ***************************************************************************/ +ctulimit::ctulimit(void) : ctool(CTULIMIT_NAME, CTULIMIT_VERSION) +{ + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Observations constructor + * + * param[in] obs Observation container. + * + * This method creates an instance of the class by copying an existing + * observations container. + ***************************************************************************/ +ctulimit::ctulimit(const GObservations& obs) : + ctool(CTULIMIT_NAME, CTULIMIT_VERSION) +{ + // Initialise members + init_members(); + + // Set observations + m_obs = obs; + + // Return + return; +} + + + +/***********************************************************************//** + * @brief Command line constructor + * + * @param[in] argc Number of arguments in command line. + * @param[in] argv Array of command line arguments. + ***************************************************************************/ +ctulimit::ctulimit(int argc, char *argv[]) : + ctool(CTULIMIT_NAME, CTULIMIT_VERSION, argc, argv) +{ + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] app Application. + ***************************************************************************/ +ctulimit::ctulimit(const ctulimit& app) : ctool(app) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(app); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +ctulimit::~ctulimit(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] app Application. + * @return Application. + ***************************************************************************/ +ctulimit& ctulimit::operator=(const ctulimit& app) +{ + // Execute only if object is not identical + if (this != &app) { + + // Copy base class members + this->ctool::operator=(app); + + // Free members + free_members(); + + // Initialise members + init_members(); + + // Copy members + copy_members(app); + + } // endif: object was not identical + + // Return this object + return *this; +} + + +/*========================================================================== + = = + = Public methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Clear instance + ***************************************************************************/ +void ctulimit::clear(void) +{ + // Free members + free_members(); + this->ctool::free_members(); + this->GApplication::free_members(); + + // Initialise members + this->GApplication::init_members(); + this->ctool::init_members(); + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Computes the upper limit + * + * This method calculates the upper limit depending on the given algorithm + ***************************************************************************/ +void ctulimit::run(void) +{ + // If we're in debug mode then all output is also dumped on the screen + if (logDebug()) { + log.cout(true); + } + + // Get task parameters + get_parameters(); + + // Write parameters into logger + if (logTerse()) { + log_parameters(); + log << std::endl; + } + + // Write observation(s) into logger + if (logTerse()) { + log << std::endl; + if (m_obs.size() > 1) { + log.header1("Observations"); + } + else { + log.header1("Observation"); + } + log << m_obs << std::endl; + } + + // Save original models + GModels models_orig = m_obs.models(); + + // Save original log-likelihood. If the value is zero it has never been + // computed hence we compute it now. + m_best_logL = m_obs.logL(); + if (m_best_logL == 0.0) { + + // Write header + if (logTerse()) { + log << std::endl; + log.header1("Compute best-fit likelihood"); + } + + // Reoptimize if likelihood was not given before + GOptimizerLM* opt = new GOptimizerLM(); + m_obs.optimize(*opt); + m_obs.errors(*opt); + m_best_logL = m_obs.logL(); + + // Write optimised model into logger + if (logTerse()) { + log << m_obs.models() << std::endl; + } + + } // endif: likelihood was zero + + // Compute parameter bracketing + double value = m_model_par->value(); + double error = m_model_par->error(); + double parmin = value + m_sigma_min * error; + double parmax = value + m_sigma_max * error; + + // Write header + if (logTerse()) { + log << std::endl; + log.header1("Compute upper limit"); + log << gammalib::parformat("Model name"); + log << m_skymodel->name() << std::endl; + log << gammalib::parformat("Parameter name"); + log << m_model_par->name() << std::endl; + log << gammalib::parformat("Initial parameter range"); + log << "["; + log << parmin; + log << ", "; + log << parmax; + log << "]" << std::endl; + } + + // Compute upper limit + ulimit_bisection(parmin, parmax); + + // Write final parameter range + if (logTerse()) { + log << gammalib::parformat("Final parameter"); + log << m_model_par->value() << std::endl; + } + + // Get reference energy for differential upper limit + GEnergy eref = GEnergy(m_eref, "TeV"); + + // Create energy range for flux limits + GEnergy emin = GEnergy(m_emin, "TeV"); + GEnergy emax = GEnergy(m_emax, "TeV"); + + // Compute upper limit intensity and fluxes + m_diff_ulimit = m_skymodel->spectral()->eval(eref, GTime()); + m_flux_ulimit = m_skymodel->spectral()->flux(emin, emax); + m_eflux_ulimit = m_skymodel->spectral()->eflux(emin, emax); + + // Write results to logfile + if (logTerse()) { + log << std::endl; + log.header1("Upper limit results"); + log << gammalib::parformat("Differential flux limit"); + log << m_diff_ulimit; + log << " ph/cm2/s/MeV at "; + log << m_eref << " TeV"; + log << std::endl; + log << gammalib::parformat("Integral flux limit"); + log << m_flux_ulimit; + log << " ph/cm2/s within ["; + log << m_emin << "-" << m_emax; + log << "] TeV"; + log << std::endl; + log << gammalib::parformat("Energy flux limit"); + log << m_eflux_ulimit; + log << " erg/cm2/s within ["; + log << m_emin << "-" << m_emax; + log << "] TeV"; + log << std::endl; + } + + // Recover original models + m_obs.models(models_orig); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Save maps + * + * This method saves the upper limit to an ascii file + ***************************************************************************/ +void ctulimit::save(void) +{ + // Write header + if (logTerse()) { + log << std::endl; + log.header1("Save upper limit"); + } + + // Get output filename + m_outfile = (*this)["outfile"].filename(); + + // Create CSV table with 3 columns + GCsv table(1, 3); + + table.real(0, 0, m_diff_ulimit); + table.real(0, 1, m_flux_ulimit); + table.real(0, 2, m_eflux_ulimit); + + // Save CSV table + table.save(m_outfile, " ", clobber()); + + // Return + return; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void ctulimit::init_members(void) +{ + // Initialise user parameters + m_outfile.clear(); + m_srcname.clear(); + m_sigma_min = 0.0; + m_sigma_max = 0.0; + m_eref = 0.0; + m_emin = 0.0; + m_emax = 0.0; + m_tol = 1.0e-6; + m_max_iter = 50; + + // Initialise protected members + m_obs.clear(); + m_models.clear(); + m_dlogL = 0.0; + m_skymodel = NULL; + m_model_par = NULL; + m_best_logL = 0.0; + m_flux_ulimit = 0.0; + m_diff_ulimit = 0.0; + m_eflux_ulimit = 0.0; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] app Application. + ***************************************************************************/ +void ctulimit::copy_members(const ctulimit& app) +{ + // Copy user parameters + m_outfile = app.m_outfile; + m_srcname = app.m_srcname; + m_sigma_min = app.m_sigma_min; + m_sigma_max = app.m_sigma_max; + m_eref = app.m_eref; + m_emin = app.m_emin; + m_emax = app.m_emax; + m_tol = app.m_tol; + m_max_iter = app.m_max_iter; + + + // Copy protected members + m_obs = app.m_obs; + m_models = app.m_models; + m_dlogL = app.m_dlogL; + m_best_logL = app.m_best_logL; + m_diff_ulimit = app.m_diff_ulimit; + m_flux_ulimit = app.m_flux_ulimit; + m_eflux_ulimit = app.m_eflux_ulimit; + + // Extract model parameter + get_model_parameter(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void ctulimit::free_members(void) +{ + // Return + return; +} + + +/***********************************************************************//** + * @brief Get application parameters + * + * @exception GException::invalid_value + * Test source not found. + * + * Get all task parameters from parameter file or (if required) by querying + * the user. + ***************************************************************************/ +void ctulimit::get_parameters(void) +{ + // If there are no observations in container then load them via user + // parameters + if (m_obs.size() == 0) { + + // Throw exception if no input observation file is given + require_inobs(G_GET_PARAMETERS); + + // Build observation container + m_obs = get_observations(); + + } // endif: there was no observation in the container + + + // If there is are no models associated with the observations then + // load now the model definition + if (m_obs.models().size() == 0) { + + // Get models XML filename + std::string filename = (*this)["inmodel"].filename(); + + // Setup models for optimizing. + m_obs.models(GModels(filename)); + + } // endif: no models were associated with observations + + // Get name of test source and check container for this name + m_srcname = (*this)["srcname"].string(); + if (!m_obs.models().contains(m_srcname)) { + std::string msg = "Source \""+m_srcname+"\" not found in model " + "container. Please add a source with that name " + "or check for possible typos."; + throw GException::invalid_value(G_GET_PARAMETERS, msg); + } + + // Get relevant model and parameter for upper limit computation + get_model_parameter(); + + // Get confidence level + double confidence = (*this)["confidence"].real(); + if (confidence != 0.95) { + std::string msg = "Confidence level different from 95% requested." + " Currently only 95% is supported."; + throw GException::invalid_value(G_GET_PARAMETERS, msg); + } + else { + // Set Likelihood difference for 95% CL. + // See Minuit Handbook + m_dlogL = 3.84 / 2.0; + } + + // Read starting boundaries for bisection + m_sigma_min = (*this)["sigma_min"].real(); + m_sigma_max = (*this)["sigma_max"].real(); + + // Read energy values + m_eref = (*this)["eref"].real(); + m_emin = (*this)["emin"].real(); + m_emax = (*this)["emax"].real(); + + // Read precision + m_tol = (*this)["tol"].real(); + m_max_iter = (*this)["max_iter"].integer(); + + // Optionally read ahead parameters so that they get correctly + // dumped into the log file + if (read_ahead()) { + m_outfile = (*this)["outfile"].filename(); + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Get application parameters + * + * @exception GException::invalid_value + * Did not find a valid model parameter + * + * Extracts a pointer to the sky model (m_skymodel) and a pointer to the + * relevant model parameter (m_model_par) from the model container. + ***************************************************************************/ +void ctulimit::get_model_parameter(void) +{ + // Get relevant model and parameter for upper limit computation. + GModels& models = const_cast(m_obs.models()); + m_skymodel = dynamic_cast(models[m_srcname]); + if (m_skymodel == NULL) { + std::string msg = "Source \""+m_srcname+"\" is not a sky model. " + "Please specify the name of a sky model for " + "upper limit computation."; + throw GException::invalid_value(G_GET_MODEL_PARAMETER, msg); + } + if (m_skymodel->spectral()->has_par("Normalization")) { + m_model_par = &(m_skymodel->spectral()->operator[]("Normalization")); + } + else if (m_skymodel->spectral()->has_par("Prefactor")) { + m_model_par = &(m_skymodel->spectral()->operator[]("Prefactor")); + } + else if (m_skymodel->spectral()->has_par("Integral")) { + m_model_par = &(m_skymodel->spectral()->operator[]("Integral")); + } + else { + std::string msg = "Require spectral parameter \"Normalization\", " + "\"Prefactor\" or \"Integral\" for upper limit " + "computation. The specified source \""+m_srcname+ + "\" does not have such a parameter."; + throw GException::invalid_value(G_GET_MODEL_PARAMETER, msg); + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Performs upper limit computation by using a bisection method + * + * @param[in] min Minimum parameter value + * @param[in] max Maximum parameter value + * + * This method calculates the upper limit using a bisection method. + ***************************************************************************/ +void ctulimit::ulimit_bisection(const double& min, const double& max) +{ + // Copy values to working values + double wrk_min = min; + double wrk_max = max; + + // Initialise iteration counter + int iter = 0; + + // Loop until breaking condition is reached + while (true) { + + // Log information + if (logExplicit()) { + log << gammalib::parformat("Iteration "+gammalib::str(iter)); + log << "["; + log << wrk_min; + log << ", "; + log << wrk_max; + log << "]" << std::endl; + } + + // Throw exception if maximum iterations are reached + if (iter > m_max_iter) { + std::string msg = "The maximum number of "+gammalib::str(m_max_iter)+ + " has been reached. You may consider to increase" + " the \"max_iter\" parameter and re-run ctulimit."; + throw GException::invalid_value(G_UL_BISECTION, msg); + } + + // Compute center of boundary + double mid = (wrk_min + wrk_max) / 2.0; + + // Calculate function value + double eval_mid = evaluate(mid); + + // Check for convergence inside tolerance + if (std::abs(eval_mid) < m_tol) { + break; + } + + // Change boundaries for further iteration + if (eval_mid > 0.0) { + wrk_max = mid; + } + else if (eval_mid < 0.0) { + wrk_min = mid; + } + + // Increment counter + iter++; + + } // endwhile + + // Return + return; + +} + + +/***********************************************************************//** + * @brief Evaluates the log-likelihood + * + * @param[in] value Parameter factor value + * @return Log-likelihood value + * + * This method evaluates the log-likelihood as a function of the parameter + * of interest. + ***************************************************************************/ +double ctulimit::evaluate(const double& value) +{ + // Initialise log-likelihood value + double logL = 0.0; + + // Check if given parameter is within boundaries + if (value > m_model_par->min() && value < m_model_par->max()) { + + // Change parameter factor + m_model_par->value(value); + + // Evaluate likelihood for new model container + m_obs.eval(); + + // Retrieve likelihood + logL = m_obs.logL(); + + } // endif: value was inside allowed range + + // ... otherwise signal that the parameter went outside the boundaries + else { + std::string msg = "Value of parameter \""+m_model_par->name()+"\" " + "outside of validity range requested. To omit " + "this error please enlarge the parameter range " + "in the model XML file."; + throw GException::invalid_value(G_EVALUATE, msg); + } + + // Compute function value + double logL_difference = logL - m_best_logL - m_dlogL; + + // Return + return logL_difference; +} diff --git a/src/ctulimit/ctulimit.hpp b/src/ctulimit/ctulimit.hpp new file mode 100644 index 0000000..c7bc326 --- /dev/null +++ b/src/ctulimit/ctulimit.hpp @@ -0,0 +1,153 @@ +/*************************************************************************** + * ctulimit - Upper limit calculation tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2015 by Michael Mayer * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file ctulimit.hpp + * @brief Upper limit calculation tool interface definition + * @author Michael Mayer + */ + +#ifndef CTULIMIT_HPP +#define CTULIMIT_HPP + +/* __ Includes ___________________________________________________________ */ +#include "GammaLib.hpp" +#include "GCTALib.hpp" +#include "ctool.hpp" + +/* __Definitions _________________________________________________________ */ +#define CTULIMIT_NAME "ctulimit" +#define CTULIMIT_VERSION "01-00-00" + + +/***********************************************************************//** + * @class ctulimit + * + * @brief Upper limit calculation tool + * + * This class computes and upper limit for a given source. + * + * The class operates on predefined observation containers, an individual + * event list or an observation definition XML file. + * + * During the computation the likelihood function is inspected to find the + * best value for the upper limit. + ***************************************************************************/ +class ctulimit : public ctool { + +public: + // Constructors and destructors + ctulimit(void); + explicit ctulimit(const GObservations& obs); + ctulimit(int argc, char *argv[]); + ctulimit(const ctulimit& app); + virtual ~ctulimit(void); + + // Operators + ctulimit& operator=(const ctulimit& app); + + // Methods + void clear(void); + void run(void); + void save(void); + const GObservations& obs(void) const; + const double& diff_ulimit(void) const; + const double& flux_ulimit(void) const; + const double& eflux_ulimit(void) const; + +protected: + // Protected methods + void init_members(void); + void copy_members(const ctulimit& app); + void free_members(void); + void get_parameters(void); + void get_model_parameter(void); + void ulimit_bisection(const double& parmin, const double& parmax); + double evaluate(const double& value); + + // User parameters + std::string m_srcname; //!< Name of source which is moved around + std::string m_outfile; //!< Output ascii file + double m_sigma_min; //!< Starting value minimum (multiple fit errors above fit values) + double m_sigma_max; //!< Starting value maximum (multiple fit errors above fit values) + double m_eref; //!< Reference energy for flux limits (TeV) + double m_emin; //!< Minimum energy for flux limits (TeV) + double m_emax; //!< Maximum energy for flux limits (TeV) + double m_tol; //!< Tolerance for limit determination + int m_max_iter; //!< Maximum number of iterations + + // Protected members + GObservations m_obs; //!< Observation container + GModels m_models; + double m_dlogL; //!< Likelihood difference for upper limit computation + GModelSky* m_skymodel; //!< Pointer to sky model + GModelPar* m_model_par; //!< Pointer to model parameter + double m_best_logL; //!< Best fit log likelihood of given model + double m_diff_ulimit; //!< Differential upper limit value + double m_flux_ulimit; //!< Flux upper limit value + double m_eflux_ulimit; //!< Energy flux upper limits +}; + + +/***********************************************************************//** + * @brief Return observation container + * + * @return Reference to observation container + ***************************************************************************/ +inline +const GObservations& ctulimit::obs(void) const +{ + return m_obs; +} + +/***********************************************************************//** + * @brief Return differential upper limit + * + * @return differential upper limit + ***************************************************************************/ +inline +const double& ctulimit::diff_ulimit(void) const +{ + return m_diff_ulimit; +} + +/***********************************************************************//** + * @brief Return flux upper limit + * + * @return flux upper limit + ***************************************************************************/ +inline +const double& ctulimit::flux_ulimit(void) const +{ + return m_flux_ulimit; +} + +/***********************************************************************//** + * @brief return energy flux upper limit + * + * @return energy flux upper limits + ***************************************************************************/ +inline +const double& ctulimit::eflux_ulimit(void) const +{ + return m_eflux_ulimit; +} + +#endif /* CTULIMIT_HPP */ diff --git a/src/ctulimit/ctulimit.par b/src/ctulimit/ctulimit.par new file mode 100644 index 0000000..0489eeb --- /dev/null +++ b/src/ctulimit/ctulimit.par @@ -0,0 +1,48 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: ctulimit.par # +# Version: 01-00-00 # +# # +# Author: Michael Mayer # +# michael.mayer@physik.hu-berlin.de # +# HU Berlin # +# # +# Purpose: Parameter file for the ctulimit application # +# # +############################################################################# +# +# General parameters +#=================== +inobs, f, a, "events.fits",,, "Input event list, counts cube or observation definition file" +inmodel, f, a, "$CTOOLS/share/models/crab.xml",,, "Source model" +srcname, s, a, "Crab",,, "Source of interest" +expcube, f, a, "NONE",,, "Exposure cube file (only needed for stacked analysis)" +psfcube, f, a, "NONE",,, "PSF cube file (only needed for stacked analysis)" +caldb, s, a, "dummy",,, "Calibration database" +irf, s, a, "cta_dummy_irf",,, "Instrument response function" +outfile, f, a, "upperlimit.txt",,, "Output ASCII file" + +# +# Upper limit calculation parameters +#=================================== +confidence, r, h, 0.95,0.0,1.0, "Confidence level" +sigma_min, r, h, 0.0,,, "Start minimum value (multiple sigmas above fit value)" +sigma_max, r, h, 10.0,,, "Start maximum value (multiple sigmas above fit value)" +eref, r, h, 1.0,,, "Reference energy for differential limits (TeV)" +emin, r, h, 1.0,,, "Minimum energy for flux limits (TeV)" +emax, r, h, 100.0,,, "Maximum energy for flux limits (TeV)" +tol, r, h, 1e-5,,, "Computation tolerance" +max_iter, i, h, 50,,, "Maximum number of iterations" + +# +# Standard parameters +#==================== +chatter, i, h, 2,0,4, "Chattiness of output" +clobber, b, h, yes,,, "Overwrite existing output files with new output files?" +debug, b, h, no,,, "Debugging mode activated" +mode, s, h, "ql",,, "Mode of automatic parameters" +logfile, f, h, "ctulimit.log",,, "Log filename" diff --git a/src/ctulimit/main.cpp b/src/ctulimit/main.cpp new file mode 100644 index 0000000..0afef65 --- /dev/null +++ b/src/ctulimit/main.cpp @@ -0,0 +1,79 @@ +/*************************************************************************** + * ctulimit - upper limit calculation tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2015 by Michael Mayer * + * ----------------------------------------------------------------------- * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + * * + ***************************************************************************/ +/** + * @file ctulimit/main.cpp + * @brief CTA upper limit generation + * @author Michael Mayer + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#include "ctulimit.hpp" + + +/***********************************************************************//** + * @brief Main entry point of application + * + * @param[in] argc Number of command line arguments. + * @param[in] argv Command line arguments. + * + * This is the main entry point of the ctulimit application. It allocates a + * ctupperlimit object and runs the application. + ***************************************************************************/ +int main (int argc, char *argv[]) +{ + // Initialise return code + int rc = 1; + + // Create instance of application + ctulimit application(argc, argv); + + // Run application + try { + // Execute application + application.execute(); + + // Signal success + rc = 0; + } + catch (std::exception &e) { + + // Extract error message + std::string message = e.what(); + std::string signal = "*** ERROR encounterted in the execution of" + " ctulimit. Run aborted ..."; + + // Write error in logger + application.log << signal << std::endl; + application.log << message << std::endl; + + // Write error on standard output + std::cout << signal << std::endl; + std::cout << message << std::endl; + + } // endcatch: catched any application error + + // Return + return rc; +} diff --git a/test/test_ctools.sh b/test/test_ctools.sh index 89597ae..0984d7c 100755 --- a/test/test_ctools.sh +++ b/test/test_ctools.sh @@ -46,12 +46,14 @@ ctpsfcube=../src/ctpsfcube/ctpsfcube ctbkgcube=../src/ctbkgcube/ctbkgcube ctcubemask=../src/ctcubemask/ctcubemask ctbutterfly=../src/ctbutterfly/ctbutterfly +ctulimit=../src/ctulimit/ctulimit # # Remove any existing result files # ================================ rm -rf *.fits *.xml +rm -rf ulimit.dat butterfly.txt # @@ -597,3 +599,25 @@ else exit 1 fi $ECHO " ok" + + +# +# Test ctulimit +# =============== +$ECHO -n "Test ctulimit: " +$ctulimit inobs="data/crab_events.fits.gz" \ + inmodel="data/crab.xml" \ + srcname="Crab" \ + outfile="ulimit.dat" \ + caldb="irf" \ + irf="cta_dummy_irf" +$ECHO -n "." +if [ -s "ulimit.dat" ] +then + $ECHO -n "." +else + $ECHO " ulimit.dat file is not found" + exit 1 +fi +$ECHO " ok" + diff --git a/test/test_ctulimit.py b/test/test_ctulimit.py new file mode 100755 index 0000000..49685f0 --- /dev/null +++ b/test/test_ctulimit.py @@ -0,0 +1,94 @@ +#! /usr/bin/env python +# ========================================================================== +# This scripts performs unit tests for the ctulimit tool. +# +# Copyright (C) 2015 Michael Mayer +# +# This program is free software: you can redistribute it and/or modify +# it under the terms of the GNU General Public License as published by +# the Free Software Foundation, either version 3 of the License, or +# (at your option) any later version. +# +# This program is distributed in the hope that it will be useful, +# but WITHOUT ANY WARRANTY; without even the implied warranty of +# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +# GNU General Public License for more details. +# +# You should have received a copy of the GNU General Public License +# along with this program. If not, see . +# +# ========================================================================== +import gammalib +import ctools + + +# ========================== # +# Test class for cttsmap tool # +# ========================== # +class Test(gammalib.GPythonTestSuite): + """ + Test class for cttsmap tool. + """ + # Constructor + def __init__(self): + """ + Constructor. + """ + # Call base class constructor + gammalib.GPythonTestSuite.__init__(self) + + # Set members + self.events_name = "data/crab_events.fits" + self.model_name = "data/crab.xml" + self.caldb = "irf" + self.irf = "cta_dummy_irf" + + # Return + return + + # Set test functions + def set(self): + """ + Set all test functions. + """ + # Set test name + self.name("ctulimit") + + # Append tests + self.append(self.test_functional, "Test ctulimit functionality") + + # Return + return + + # Test cttsmap functionnality + def test_functional(self): + """ + Test cttsmap functionnality. + """ + # Set-up cttsmap + ulimit = ctools.ctulimit() + ulimit["inobs"].filename(self.events_name) + ulimit["inmodel"].filename(self.model_name) + ulimit["srcname"].string("Crab") + ulimit["caldb"].string(self.caldb) + ulimit["irf"].string(self.irf) + ulimit["outfile"].filename("ulimit.dat") + + # Run tool + self.test_try("Run ctulimit") + try: + ulimit.run() + self.test_try_success() + except: + self.test_try_failure("Exception occured in ctulimit.") + + # Save results + self.test_try("Save results") + try: + ulimit.save() + self.test_try_success() + except: + self.test_try_failure("Exception occured in saving results.") + + # Return + return diff --git a/test/test_python.py b/test/test_python.py index 6dd5aaa..7da0400 100755 --- a/test/test_python.py +++ b/test/test_python.py @@ -34,6 +34,7 @@ import test_ctpsfcube import test_ctbkgcube import test_ctcubemask import test_ctbutterfly +import test_ctulimit import test_pipelines @@ -60,6 +61,7 @@ if __name__ == '__main__': suite_ctbkgcube = test_ctbkgcube.Test() suite_ctcubemask = test_ctcubemask.Test() suite_ctbutterfly = test_ctbutterfly.Test() + suite_ctulimit = test_ctulimit.Test() suite_pipelines = test_pipelines.Test() # Setup unit tests @@ -75,6 +77,7 @@ if __name__ == '__main__': suite_ctbkgcube.set() suite_ctcubemask.set() suite_ctbutterfly.set() + suite_ctulimit.set() suite_pipelines.set() # Append tests to container @@ -90,6 +93,7 @@ if __name__ == '__main__': suites.append(suite_ctbkgcube) suites.append(suite_ctcubemask) suites.append(suite_ctbutterfly) + suites.append(suite_ctulimit) suites.append(suite_pipelines) # Set PFILES environment variable