From 6762b992b9297f953c4537cb1053d8a08d6504a4 Mon Sep 17 00:00:00 2001 From: Michael Mayer Date: Wed, 12 Nov 2014 11:55:51 +0100 Subject: [PATCH] added structure for ctulimit --- pyext/ctools/tools.i | 1 + src/Makefile.am | 6 ++++-- src/ctulimit/Makefile.am | 24 ++++++++++++++++++++++++ src/ctulimit/ctulimit.cpp | 410 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/ctulimit/ctulimit.hpp | 120 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/ctulimit/ctulimit.par | 37 +++++++++++++++++++++++++++++++++++++ src/ctulimit/main.cpp | 79 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 7 files changed, 675 insertions(+), 2 deletions(-) create mode 100644 src/ctulimit/Makefile.am create mode 100644 src/ctulimit/ctulimit.cpp create mode 100644 src/ctulimit/ctulimit.hpp create mode 100644 src/ctulimit/ctulimit.par create mode 100644 src/ctulimit/main.cpp 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/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..6925e1b --- /dev/null +++ b/src/ctulimit/ctulimit.cpp @@ -0,0 +1,410 @@ +/*************************************************************************** + * ctulimitt - upper limit calculation tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2014 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" + +/* __ Method name definitions ____________________________________________ */ +#define G_GET_PARAMETERS "ctulimit::get_parameters()" + +/* __ 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; + } + + // Write header + if (logTerse()) { + log << std::endl; + log.header1("Compute upper limit"); + } + + + // Compute upper limit here + + // 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 TS map"); + } + + // Get output filename + m_outfile = (*this)["outfile"].filename(); + + + // Return + return; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void ctulimit::init_members(void) +{ + // Initialise members + m_infile.clear(); + m_outfile.clear(); + + + // Initialise protected members + m_obs.clear(); + + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] app Application. + ***************************************************************************/ +void ctulimit::copy_members(const ctulimit& app) +{ + // Copy attributes + m_infile = app.m_infile; + m_outfile = app.m_outfile; + m_modelfile = app.m_modelfile; + m_srcname = app.m_srcname; + m_parname = app.m_parname; + m_algorithm = app.m_algorithm; + + // Copy protected members + m_obs = app.m_obs; + + // 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 or no RA/DEC parameters found for test + * source. + * + * Get all task parameters from parameter file or (if required) by querying + * the user. Most parameters are only required if no observation exists so + * far in the observation container. In this case, a single CTA observation + * will be added to the container, using the definition provided in the + * parameter file. + ***************************************************************************/ +void ctulimit::get_parameters(void) +{ + // If there are no observations in container then add a single CTA + // observation using the parameters from the parameter file + if (m_obs.size() == 0) { + + // Allocate CTA observation + GCTAObservation obs; + + // Get event file name + std::string filename = (*this)["infile"].filename(); + + // Try first to open as FITS file + try { + + // Load data + obs.load(filename); + + // Set response + set_obs_response(&obs); + + // Append observation to container + m_obs.append(obs); + + } + + // ... otherwise try to open as XML file + catch (GException::fits_open_error &e) { + + // Load observations from XML file + m_obs.load(filename); + + // Check if all observations have response information. If + // not, get the calibration database parameters and set + // the response properly + set_response(m_obs); + + } // endcatch: file was an XML file + + } // 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)["srcmdl"].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 a possible typos."; + throw GException::invalid_value(G_GET_PARAMETERS, msg); + } + + + // Optionally read ahead parameters so that they get correctly + // dumped into the log file + if (read_ahead()) { + m_outfile = (*this)["outfile"].filename(); + } + + // TODO read other parameters here + + // Return + return; +} + + diff --git a/src/ctulimit/ctulimit.hpp b/src/ctulimit/ctulimit.hpp new file mode 100644 index 0000000..b100ea7 --- /dev/null +++ b/src/ctulimit/ctulimit.hpp @@ -0,0 +1,120 @@ +/*************************************************************************** + * ctulimit - Upper limit calculation tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2014 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 "00-01-00" + + +/***********************************************************************//** + * @class ctulimit + * + * @brief Upper limit calculation tool + * + * This class computes and upper limit for a given parameter. + * + * 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& ulimit(void) const; + +protected: + // Protected methods + void init_members(void); + void copy_members(const ctulimit& app); + void free_members(void); + void get_parameters(void); + + // User parameters + std::string m_infile; //!< Input file + std::string m_modelfile;//!< Name of model optimised model file + std::string m_srcname; //!< Name of source which is moved around + std::string m_parname; //!< Name of parameter upper limit should be computed + std::string m_outfile; //!< Output counts map or XML file + std::string m_algorithm; //!< algorithm of upper limit computation + + // Protected members + GObservations m_obs; //!< Observation container + + double m_ulimit; //!. * + * * + ***************************************************************************/ +/** + * @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" + " cttsmap. 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; +} -- libgit2 0.21.4