From f048b0ed5a5ea50c751f03e9af93fa7db7660422 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Mon, 15 Oct 2012 23:47:57 +0200 Subject: [PATCH] Create ctmodel tool. --- src/ctmodel/Makefile.am | 21 +++++++++++++++++++++ src/ctmodel/ctmodel.cpp | 854 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/ctmodel/ctmodel.hpp | 115 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/ctmodel/ctmodel.par | 57 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++ src/ctmodel/main.cpp | 79 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 5 files changed, 1126 insertions(+), 0 deletions(-) create mode 100644 src/ctmodel/Makefile.am create mode 100644 src/ctmodel/ctmodel.cpp create mode 100644 src/ctmodel/ctmodel.hpp create mode 100644 src/ctmodel/ctmodel.par create mode 100644 src/ctmodel/main.cpp diff --git a/src/ctmodel/Makefile.am b/src/ctmodel/Makefile.am new file mode 100644 index 0000000..a968db2 --- /dev/null +++ b/src/ctmodel/Makefile.am @@ -0,0 +1,21 @@ +## Process this file with automake to produce Makefile.in + +# Program name +bin_PROGRAMS = ctmodel + +# Sources +sources = ctmodel.cpp + +# Libtool library +noinst_LTLIBRARIES = libctmodel.la +libctmodel_la_SOURCES = $(sources) + +# Program definition +ctmodel_SOURCES = main.cpp +ctmodel_LDFLAGS = $(LIB_LINK_OPTION) @LDFLAGS@ +ctmodel_CPPFLAGS = @CPPFLAGS@ +ctmodel_LDADD = libctmodel.la + +# Parameter files +pfilesdir = $(prefix)/syspfiles +pfiles_DATA = ctmodel.par diff --git a/src/ctmodel/ctmodel.cpp b/src/ctmodel/ctmodel.cpp new file mode 100644 index 0000000..0ce9a7a --- /dev/null +++ b/src/ctmodel/ctmodel.cpp @@ -0,0 +1,854 @@ +/*************************************************************************** + * ctmodel - CTA counts model tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2012 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * 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 ctmodel.cpp + * @brief CTA counts model tool implementation + * @author Juergen Knoedlseder + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#include "ctmodel.hpp" +#include "GTools.hpp" + +/* __ Method name definitions ____________________________________________ */ +#define G_SETUP_OBS "ctmodel::setup_obs()" +#define G_MODEL_MAP "ctmodel::model_map(GCTAObservation*)" + +/* __ Debug definitions __________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + ***************************************************************************/ +ctmodel::ctmodel(void) : GApplication(CTMODEL_NAME, CTMODEL_VERSION) +{ + // Initialise members + init_members(); + + // Write header into logger + log_header(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Observations constructor + * + * This method creates an instance of the class by copying an existing + * observations container. + ***************************************************************************/ +ctmodel::ctmodel(GObservations obs) : GApplication(CTMODEL_NAME, CTMODEL_VERSION) +{ + // Initialise members + init_members(); + + // Set observations + m_obs = obs; + + // Write header into logger + log_header(); + + // Return + return; +} + + + +/***********************************************************************//** + * @brief Command line constructor + * + * @param[in] argc Number of arguments in command line. + * @param[in] argv Array of command line arguments. + ***************************************************************************/ +ctmodel::ctmodel(int argc, char *argv[]) : + GApplication(CTMODEL_NAME, CTMODEL_VERSION, argc, argv) +{ + // Initialise members + init_members(); + + // Write header into logger + log_header(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] app Application. + ***************************************************************************/ +ctmodel::ctmodel(const ctmodel& app) : GApplication(app) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(app); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +ctmodel::~ctmodel(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] app ctmodel application. + * @return Returns ctmodel application. + ***************************************************************************/ +ctmodel& ctmodel::operator= (const ctmodel& app) +{ + // Execute only if object is not identical + if (this != &app) { + + // Copy base class members + this->GApplication::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 ctmodel::clear(void) +{ + // Free members + free_members(); + this->GApplication::free_members(); + + // Initialise members + this->GApplication::init_members(); + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Execute application + * + * This is the main execution method of the ctmodel class. It is invoked + * when the executable is called from command line. The method generates + * the model maps and saves the results. + ***************************************************************************/ +void ctmodel::execute(void) +{ + // Signal that some parameters should be read ahead + m_read_ahead = true; + + // Create the model map(s) + run(); + + // Save the model map(s) into FITS file + save(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Generate the model map(s) + * + * This method reads the task parameters from the parfile, sets up the + * observation container, loops over all CTA observations in the container + * and generates a model map for each CTA observation. + ***************************************************************************/ +void ctmodel::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(); + + // Setup observation container + setup_obs(); + + // 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; + if (m_obs.size() > 1) { + log.header1("Generate model maps"); + } + else { + log.header1("Generate model map"); + } + } + + // Initialise observation counter + int n_observations = 0; + + // Loop over all observations in the container + for (int i = 0; i < m_obs.size(); ++i) { + + // Initialise event input and output filenames + m_infiles.push_back(""); + + // Get CTA observation + GCTAObservation* obs = dynamic_cast(&m_obs[i]); + + // Continue only if observation is a CTA observation + if (obs != NULL) { + + // Write header for observation + if (logTerse()) { + if (obs->name().length() > 1) { + log.header3("Observation "+obs->name()); + } + else { + log.header3("Observation"); + } + } + + // Increment number of observations + n_observations++; + + // Save event file name (for possible saving) + m_infiles[i] = obs->eventfile(); + + // Generate model map + model_map(obs, m_obs.models()); + + } // endif: CTA observation found + + } // endfor: looped over observations + + // If more than a single observation has been handled then make sure + // that an XML file will be used for storage + if (n_observations > 1) { + m_use_xml = true; + } + + // Write observation(s) into logger + if (logTerse()) { + log << std::endl; + if (m_obs.size() > 1) { + log.header1("Observations after model map generation"); + } + else { + log.header1("Observation after model map generation"); + } + log << m_obs << std::endl; + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Save model map(s) + * + * This method saves the model map(s) into FITS file(s). There are two + * modes, depending on the m_use_xml flag. + * + * If m_use_xml is true, all model map(s) will be saved into FITS files, + * where the output filenames are constructued from the input filenames by + * prepending the m_prefix string to name. Any path information will be + * stripped form the input name, hence event files will be written into the + * local working directory (unless some path information is present + * in the prefix). In addition, an XML file will be created that gathers + * the filename information for the model map(s). If an XML file was present + * on input, all metadata information will be copied from this input file. + * + * If m_use_xml is false, the model map will be saved into a FITS file. + ***************************************************************************/ +void ctmodel::save(void) +{ + // Write header + if (logTerse()) { + log << std::endl; + if (m_obs.size() > 1) { + log.header1("Save observations"); + } + else { + log.header1("Save observation"); + } + } + + // Case A: Save model map(s) and XML metadata information + if (m_use_xml) { + save_xml(); + } + + // Case B: Save model map as FITS file + else { + save_fits(); + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Get application parameters + * + * Get all task parameters from parameter file or (if required) by querying + * the user. The parameters are read in the correct order. + ***************************************************************************/ +void ctmodel::get_parameters(void) +{ + // If we do not have any observations in the container then get all + // information to build (at least) one + if (m_obs.size() == 0) { + + // Read general parameters. Only read the output parameters if + // required, and also only read the model XML filename if no + // model is yet in the container. + m_infile = (*this)["infile"].filename(); + if (m_read_ahead) { + m_outfile = (*this)["outfile"].filename(); + m_prefix = (*this)["prefix"].string(); + } + m_caldb = (*this)["caldb"].string(); + m_irf = (*this)["irf"].string(); + if (m_obs.models().size() == 0) { + m_srcmdl = (*this)["srcmdl"].filename(); + } + + // If there is no input filename then read all parameters that + // are required to build a model map from scratch + if ((m_infile == "NONE") || (strip_whitespace(m_infile) == "")) { + m_tmin = (*this)["tmin"].real(); + m_tmax = (*this)["tmax"].real(); + m_emin = (*this)["emin"].real(); + m_emax = (*this)["emax"].real(); + m_enumbins = (*this)["enumbins"].integer(); + m_proj = (*this)["proj"].string(); + m_coordsys = (*this)["coordsys"].string(); + m_xref = (*this)["xref"].real(); + m_yref = (*this)["yref"].real(); + m_binsz = (*this)["binsz"].real(); + m_nxpix = (*this)["nxpix"].integer(); + m_nypix = (*this)["nypix"].integer(); + } + } + + // ... otherwise, read only the parameters that are required for the + // model generation. + else { + if (m_read_ahead) { + m_outfile = (*this)["outfile"].filename(); + m_prefix = (*this)["prefix"].string(); + } + m_caldb = (*this)["caldb"].string(); + m_irf = (*this)["irf"].string(); + if (m_obs.models().size() == 0) { + m_srcmdl = (*this)["srcmdl"].filename(); + } + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Setup observation container + * + * @exception GException::no_cube + * No event cube found in CTA observation. + * + * This method sets up the observation container for processing. There are + * two cases: + * + * If there are no observations in the actual observation container, the + * method will check in "infile" parameter. If this parameter is "NONE" or + * empty, the task parameters will be used to construct a model map. + * Otherwise, the method first tries to interpret the "infile" parameter as + * a counts map, and attemps loading of the file in an event cube. If this + * fails, the method tries to interpret the "infile" parameter as an + * observation definition XML file. If this also fails, an exception will + * be thrown. + * + * If observations exist already in the observation container, the method + * will simply keep them. + * + * Test if all CTA observations contain counts maps. + * + * Finally, if no models exist so far in the observation container, the + * models will be loaded from the model XML file. + ***************************************************************************/ +void ctmodel::setup_obs(void) +{ + // If there are no observations in the container then try to build some + if (m_obs.size() == 0) { + + // If no input filename has been specified, then create a model map + // from the task parameters + if ((m_infile == "NONE") || (strip_whitespace(m_infile) == "")) { + + // Setup energy range covered by model + GEnergy emin; + GEnergy emax; + GEbounds ebds; + emin.TeV(m_emin); + emax.TeV(m_emax); + ebds.setlog(emin, emax, m_enumbins); + + // Setup time interval covered by model + GGti gti; + GTime tmin; + GTime tmax; + tmin.met(m_tmin); + tmax.met(m_tmax); + gti.append(tmin, tmax); + + // Setup skymap + GSkymap map = GSkymap(m_proj, m_coordsys, + m_xref, m_yref, m_binsz, m_binsz, + m_nxpix, m_nypix, m_enumbins); + + // Create model cube from sky map + GCTAEventCube cube(map, ebds, gti); + + // Allocate CTA observation + GCTAObservation obs; + + // Set event cube in observation + obs.events(&cube); + + // Append CTA observation to container + m_obs.append(obs); + + // Signal that no XML file should be used for storage + m_use_xml = false; + + } // endif: created model map from task parameters + + // ... otherwise try to load information from the file + else { + + // First try to open the file as a counts map + try { + + // Allocate CTA observation + GCTAObservation obs; + + // Load counts map in CTA observation + obs.load_binned(m_infile); + + // Append CTA observation to container + m_obs.append(obs); + + // Signal that no XML file should be used for storage + m_use_xml = false; + + } + + // ... otherwise try to open as XML file + catch (GException::fits_open_error &e) { + + // Load observations from XML file. This will throw + // an exception if it fails. + m_obs.load(m_infile); + + // Signal that XML file should be used for storage + m_use_xml = true; + + } + + } // endelse: loaded information from input file + + } // endif: there was no observation in the container + + // If there are no models associated with the observations then + // load now the model definition from the XML file + if (m_obs.models().size() == 0) { + m_obs.models(GModels(m_srcmdl)); + } + + // Check if all CTA observations contain an event cube and setup response + // for all observations + for (int i = 0; i < m_obs.size(); ++i) { + + // Is this observation a CTA observation? + GCTAObservation* obs = dynamic_cast(&m_obs[i]); + + // Yes ... + if (obs != NULL) { + + // Throw an exception if this observation does not contain + // an event cube + if (dynamic_cast(obs->events()) == NULL) { + throw GException::no_cube(G_SETUP_OBS); + } + + // Set response + obs->response(m_irf, m_caldb); + + } // endif: observation was a CTA observation + } // endfor: looped over all observations + + // Return + return; +} + + +/***********************************************************************//** + * @brief Generate model map + * + * @param[in] obs CTA observation pointer. + * @param[in] models Model container. + * + * @exception GException::no_cube + * No event cube found in CTA observation. + ***************************************************************************/ +void ctmodel::model_map(GCTAObservation* obs, const GModels& models) +{ + // Continue only if observation pointer is valid + if (obs != NULL) { + + // Get event cube pointer + GCTAEventCube* cube = + const_cast(dynamic_cast(obs->events())); + + // Throw an exception if the observation does not hold and event + // cube + if (cube == NULL) { + throw GException::no_cube(G_MODEL_MAP); + } + + // Initialise statistics + double sum = 0.0; + + // Loop over all events in counts map + for (int i = 0; i < cube->size(); ++i) { + + // Get event bin + GCTAEventBin* bin = (*cube)[i]; + + // Compute model value for event bin + double model = + models.eval(*(const_cast(bin)), *obs) * + bin->size(); + + // Store value + bin->counts(model); + + // Sum all events + sum += model; + } + + // Log results + if (logTerse()) { + log << parformat("Model events in cube"); + log << sum << std::endl; + } + + } // endif: observation pointer was not valid + + // Return + return; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void ctmodel::init_members(void) +{ + // Initialise members + m_infile.clear(); + m_outfile.clear(); + m_prefix.clear(); + m_caldb.clear(); + m_irf.clear(); + m_srcmdl.clear(); + m_proj.clear(); + m_coordsys.clear(); + m_tmin = 0.0; + m_tmax = 0.0; + m_deadc = 1.0; + m_emin = 0.0; + m_emax = 0.0; + m_enumbins = 0; + m_xref = 0.0; + m_yref = 0.0; + m_binsz = 0.0; + m_nxpix = 0; + m_nypix = 0; + + // Initialise protected members + m_obs.clear(); + m_infiles.clear(); + m_use_xml = false; + m_read_ahead = false; + + // Set logger properties + log.date(true); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] app Application. + ***************************************************************************/ +void ctmodel::copy_members(const ctmodel& app) +{ + // Copy attributes + m_infile = app.m_infile; + m_outfile = app.m_outfile; + m_prefix = app.m_prefix; + m_caldb = app.m_caldb; + m_irf = app.m_irf; + m_srcmdl = app.m_srcmdl; + m_proj = app.m_proj; + m_coordsys = app.m_coordsys; + m_tmin = app.m_tmin; + m_tmax = app.m_tmax; + m_deadc = app.m_deadc; + m_emin = app.m_emin; + m_emax = app.m_emax; + m_enumbins = app.m_enumbins; + m_xref = app.m_xref; + m_yref = app.m_yref; + m_binsz = app.m_binsz; + m_nxpix = app.m_nxpix; + m_nypix = app.m_nypix; + + // Copy protected members + m_obs = app.m_obs; + m_infiles = app.m_infiles; + m_use_xml = app.m_use_xml; + m_read_ahead = app.m_read_ahead; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void ctmodel::free_members(void) +{ + // Write separator into logger + if (logTerse()) { + log << std::endl; + } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Set output file name. + * + * @param[in] filename Input file name. + * + * This method converts an input filename into an output filename by + * prepending the prefix specified by m_prefix to the input filename. Any + * path will be stripped from the input filename. + ***************************************************************************/ +std::string ctmodel::set_outfile_name(const std::string& filename) const +{ + // Split input filename into path elements + std::vector elements = split(filename, "/"); + + // The last path element is the filename + std::string outname = m_prefix + elements[elements.size()-1]; + + // Return output filename + return outname; +} + + +/***********************************************************************//** + * @brief Save model map in FITS format. + * + * Save the model map as a FITS file. The filename of the FITS file is + * specified by the m_outfile member. + ***************************************************************************/ +void ctmodel::save_fits(void) +{ + // Get output filename + m_outfile = (*this)["outfile"].filename(); + + // Get CTA observation from observation container + GCTAObservation* obs = dynamic_cast(&m_obs[0]); + + // Save model map + save_model_map(obs, m_outfile); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Save model map(s) in XML format. + * + * Save the model map(s) into FITS files and write the file path information + * into a XML file. The filename of the XML file is specified by the + * m_outfile member, the filename(s) of the model map(s) are built by + * prepending the prefix given by the m_prefix member to the input model + * map(s) filenames. Any path present in the input filename will be stripped, + * i.e. the model map(s) will be written in the local working directory + * (unless a path is specified in the m_prefix member). + ***************************************************************************/ +void ctmodel::save_xml(void) +{ + // Get output filename and prefix + m_outfile = (*this)["outfile"].filename(); + m_prefix = (*this)["prefix"].string(); + + // Loop over all observation in the container + for (int i = 0; i < m_obs.size(); ++i) { + + // Get CTA observation + GCTAObservation* obs = dynamic_cast(&m_obs[i]); + + // Handle only CTA observations + if (obs != NULL) { + + // Set event output file name + std::string outfile = set_outfile_name(m_infiles[i]); + + // Store output file name in observation + obs->eventfile(outfile); + + // Save event list + save_model_map(obs, outfile); + + } // endif: observation was a CTA observations + + } // endfor: looped over observations + + // Save observations in XML file + m_obs.save(m_outfile); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Save a single model map into a FITS file + * + * @param[in] obs Pointer to CTA observation. + * @param[in] outfile Output file name. + * + * This method saves a single model map into a FITS file. The method does + * nothing if the observation pointer is not valid. + ***************************************************************************/ +void ctmodel::save_model_map(const GCTAObservation* obs, + const std::string& outfile) const +{ + // Save only if observation is valid + if (obs != NULL) { + + // Save observation into FITS file + obs->save(outfile, clobber()); + + } // endif: observation was valid + + // Return + return; +} diff --git a/src/ctmodel/ctmodel.hpp b/src/ctmodel/ctmodel.hpp new file mode 100644 index 0000000..89a69be --- /dev/null +++ b/src/ctmodel/ctmodel.hpp @@ -0,0 +1,115 @@ +/*************************************************************************** + * ctmodel - CTA counts model tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2012 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * 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 ctmodel.hpp + * @brief CTA counts model tool definition + * @author Juergen Knoedlseder + */ + +#ifndef CTMODEL_HPP +#define CTMODEL_HPP + +/* __ Includes ___________________________________________________________ */ +#include "GammaLib.hpp" +#include "GCTALib.hpp" + +/* __Definitions _________________________________________________________ */ +#define CTMODEL_NAME "ctmodel" +#define CTMODEL_VERSION "00-01-00" + + +/***********************************************************************//** + * @class ctmodel + * + * @brief CTA counts model tool interface defintion + * + * This class creates counts model map(s). The definition of the counts model + * can be taken from a predefined observation container, from the counts maps + * found in an observation definition XML file, from an individual counts + * map, or from the parameter definition. + * + * Results are stored in an observation container that can be written to disk + * in form of FITS files (model maps) and an updated observation definition + * XML file. + ***************************************************************************/ +class ctmodel : public GApplication { + +public: + // Constructors and destructors + ctmodel(void); + explicit ctmodel(GObservations obs); + ctmodel(int argc, char *argv[]); + ctmodel(const ctmodel& app); + virtual ~ctmodel(void); + + // Operators + ctmodel& operator= (const ctmodel& app); + + // Methods + void clear(void); + void execute(void); + void run(void); + void save(void); + GObservations& obs(void) { return m_obs; } + void get_parameters(void); + void setup_obs(void); + void model_map(GCTAObservation* obs, const GModels& models); + +protected: + // Protected methods + void init_members(void); + void copy_members(const ctmodel& app); + void free_members(void); + std::string set_outfile_name(const std::string& filename) const; + void save_fits(void); + void save_xml(void); + void save_model_map(const GCTAObservation* obs, + const std::string& outfile) const; + + // User parameters + std::string m_infile; //!< Input counts map or XML file + std::string m_outfile; //!< Output model map or XML file + std::string m_prefix; //!< Prefix for multiple model maps + std::string m_caldb; //!< Calibration database + std::string m_irf; //!< Instrument response functions + std::string m_srcmdl; //!< Source model + double m_tmin; //!< Start time + double m_tmax; //!< Stop time + double m_deadc; //!< Deadtime correction factor + double m_emin; //!< Lower energy + double m_emax; //!< Upper energy + int m_enumbins; //!< Number of energy bins + std::string m_proj; //!< WCS projection + std::string m_coordsys; //!< Coordinate system + double m_xref; //!< Longitude reference coordinate + double m_yref; //!< Latitude reference coordinate + double m_binsz; //!< Pixel size + int m_nxpix; //!< Number of pixels in longitude + int m_nypix; //!< Number of pixels in latitude + + // Protected members + GObservations m_obs; //!< Observation container + std::vector m_infiles; //!< Input event filenames + bool m_use_xml; //!< Use XML file instead of FITS file + bool m_read_ahead; //!< Read ahead parameters +}; + +#endif /* CTMODEL_HPP */ diff --git a/src/ctmodel/ctmodel.par b/src/ctmodel/ctmodel.par new file mode 100644 index 0000000..cab13e1 --- /dev/null +++ b/src/ctmodel/ctmodel.par @@ -0,0 +1,57 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: ctmodel.par # +# Version: 00-01-00 # +# # +# Author: Juergen Knoedlseder # +# Jurgen.Knodlseder@irap.omp.eu # +# IRAP # +# # +# Purpose: Parameter file for the ctmodel application # +# # +############################################################################# +# +# General parameters +#=================== +infile, f, a, "NONE",,, "Input counts map or observation definition file" +outfile, f, a, "modmap.fits",,, "Output counts map or observation definition file" +prefix, s, h, "modmap_",,, "Prefix for counts maps in observation definition file" +caldb, s, a, "test/irf",,, "Calibration database" +irf, s, a, "kb_E_50h_v3",,, "Instrument response function" +srcmdl, f, a, "test/data/crab.xml",,, "Source model" +# +# Time range parameters +#====================== +tmin, r, a, 0.0, 0, , "Start time (MET in s)" +tmax, r, a, 1800.0, 0, , "End time (MET in s)" +deadc, r, h, 0.95, 0, 1, "Average deadtime correction factor" +# +# Energy binning parameters +#========================== +ebinalg, s, h, "LOG", FILE|LIN|LOG,,"Algorithm for defining energy bins" +emin, r, a, 0.1,,,"Start value for first energy bin in TeV" +emax, r, a, 100.0,,,"Stop value for last energy bin in TeV" +enumbins, i, a, 20,,,"Number of energy bins" +ebinfile, f, h, "NONE",,,"Name of the file containing the energy bin definition" +# +# Spatial binning parameters +#=========================== +nxpix, i, a, 200,,, "Size of the X axis in pixels" +nypix, i, a, 200,,, "Size of the Y axis in pixels" +binsz, r, a, 0.02,,, "Image scale (in degrees/pixel)" +coordsys, s, a, "CEL", CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)" +xref, r, a, 83.63,,, "First coordinate of image center in degrees (RA or galactic l)" +yref, r, a, 22.01,,, "Second coordinate of image center in degrees (DEC or galactic b)" +axisrot, r, h, 0.0,,, "Rotation angle of image axis, in degrees" +proj, s, a, "CAR",,, "Projection method e.g. AIT|ARC|CAR|GLS|MER|NCP|SIN|STG|TAN" +# +# Standard parameters +#==================== +chatter, i, h, 1,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" diff --git a/src/ctmodel/main.cpp b/src/ctmodel/main.cpp new file mode 100644 index 0000000..7240904 --- /dev/null +++ b/src/ctmodel/main.cpp @@ -0,0 +1,79 @@ +/*************************************************************************** + * ctmodel - CTA counts model tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2012 by Juergen Knoedlseder * + * ----------------------------------------------------------------------- * + * * + * 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 main.cpp + * @brief CTA counts model tool main code + * @author Juergen Knoedlseder + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#include "ctmodel.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 gtbin application. It allocates a + * ctmodel object and runs the application. + ***************************************************************************/ +int main (int argc, char *argv[]) +{ + // Initialise return code + int rc = 1; + + // Create instance of application + ctmodel 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" + " ctmodel. 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