diff --git a/configure.ac b/configure.ac index da38597..b9e546d 100644 --- a/configure.ac +++ b/configure.ac @@ -957,6 +957,7 @@ AC_CONFIG_FILES([Makefile src/ctselect/Makefile src/ctexpcube/Makefile src/ctpsfcube/Makefile + src/ctedispcube/Makefile src/ctbkgcube/Makefile src/ctcubemask/Makefile src/cttsmap/Makefile diff --git a/src/Makefile.am b/src/Makefile.am index 25b1a9b..dddc604 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -14,6 +14,7 @@ SUBDIRS = support \ ctskymap \ ctexpcube \ ctpsfcube \ + ctedispcube \ ctbkgcube \ ctcubemask \ cttsmap \ @@ -40,6 +41,7 @@ libctools_la_LIBADD = support/libsupport.la \ ctskymap/libctskymap.la \ ctexpcube/libctexpcube.la \ ctpsfcube/libctpsfcube.la \ + ctedispcube/libctedispcube.la \ ctbkgcube/libctbkgcube.la \ ctcubemask/libctcubemask.la \ cttsmap/libcttsmap.la \ diff --git a/src/ctedispcube/Makefile.am b/src/ctedispcube/Makefile.am new file mode 100644 index 0000000..bf905d9 --- /dev/null +++ b/src/ctedispcube/Makefile.am @@ -0,0 +1,25 @@ +## Process this file with automake to produce Makefile.in + +# Set include directories +INCLUDES = -I$(top_srcdir)/src/support + +# Program name +bin_PROGRAMS = ctedispcube + +# Sources +sources = ctedispcube.cpp + +# Libtool library +noinst_LTLIBRARIES = libctedispcube.la +libctedispcube_la_SOURCES = $(sources) +libctedispcube_la_LDFLAGS = @OPENMP_LIB_LDFLAGS@ + +# Program definition +ctedispcube_SOURCES = main.cpp +ctedispcube_LDFLAGS = $(LIB_LINK_OPTION) @LDFLAGS@ +ctedispcube_CPPFLAGS = @CPPFLAGS@ +ctedispcube_LDADD = libctedispcube.la $(top_srcdir)/src/support/libsupport.la + +# Parameter files +pfilesdir = $(prefix)/syspfiles +pfiles_DATA = ctedispcube.par diff --git a/src/ctedispcube/ctedispcube.cpp b/src/ctedispcube/ctedispcube.cpp new file mode 100644 index 0000000..c993e0d --- /dev/null +++ b/src/ctedispcube/ctedispcube.cpp @@ -0,0 +1,414 @@ +/*************************************************************************** + * ctedispcube - EDISP cube generation tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2016 by Maria Haupt * + * ----------------------------------------------------------------------- * + * * + * 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 ctedispcube.cpp + * @brief EDISP cube generation tool implementation + * @author Maria Haupt + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#include "ctedispcube.hpp" +#include "GTools.hpp" +#include "GWcs.hpp" + +/* __ Method name definitions ____________________________________________ */ +#define G_GET_PARAMETERS "ctedispcube::get_parameters()" + +/* __ Debug definitions __________________________________________________ */ + +/* __ Coding definitions _________________________________________________ */ + + +/*========================================================================== + = = + = Constructors/destructors = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Void constructor + ***************************************************************************/ +ctedispcube::ctedispcube(void) : ctool(CTEDISPCUBE_NAME, CTEDISPCUBE_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. + ***************************************************************************/ +ctedispcube::ctedispcube(const GObservations& obs) : + ctool(CTEDISPCUBE_NAME, CTEDISPCUBE_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. + ***************************************************************************/ +ctedispcube::ctedispcube(int argc, char *argv[]) : + ctool(CTEDISPCUBE_NAME, CTEDISPCUBE_VERSION, argc, argv) +{ + // Initialise members + init_members(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy constructor + * + * @param[in] app Application. + ***************************************************************************/ +ctedispcube::ctedispcube(const ctedispcube& app) : ctool(app) +{ + // Initialise members + init_members(); + + // Copy members + copy_members(app); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Destructor + ***************************************************************************/ +ctedispcube::~ctedispcube(void) +{ + // Free members + free_members(); + + // Return + return; +} + + +/*========================================================================== + = = + = Operators = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Assignment operator + * + * @param[in] app Application. + * @return Application. + ***************************************************************************/ +ctedispcube& ctedispcube::operator=(const ctedispcube& 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 ctedispcube::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 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 an EDISP cube from the CTA observations. + ***************************************************************************/ +void ctedispcube::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; + } + + // Set energy dispersion flag for all CTA observations and save old + // values in save_edisp vector +// std::vector save_edisp; +// save_edisp.assign(m_obs.size(), false); +// for (int i = 0; i < m_obs.size(); ++i) { +// GCTAObservation* obs = dynamic_cast(m_obs[i]); +// if (obs != NULL) { +// save_edisp[i] = obs->response()->apply_edisp(); +// obs->response()->apply_edisp(m_apply_edisp); +// } +// } + + // 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("Generate EDISP cube"); + } + + // Fill EDISP + m_edispcube.fill(m_obs, &log); + + // Restore energy dispersion flag for all CTA observations +// for (int i = 0; i < m_obs.size(); ++i) { +// GCTAObservation* obs = dynamic_cast(m_obs[i]); +// if (obs != NULL) { +// obs->response()->apply_edisp(save_edisp[i]); +// } +// } + + // Return + return; +} + + +/***********************************************************************//** + * @brief Save EDISP cube + ***************************************************************************/ +void ctedispcube::save(void) +{ + // Write header + if (logTerse()) { + log << std::endl; + log.header1("Save EDISP cube"); + } + + // Get output filename + m_outcube = (*this)["outcube"].filename(); + + // Save only if filename is non-empty + if (!m_outcube.is_empty()) { + + // Log filename + if (logTerse()) { + log << "Save EDISP cube into file \""+m_outcube+"\"."; + log << std::endl; + } + + // Save EDISP cube + m_edispcube.save(m_outcube, clobber()); + + } + + // Return + return; +} + + +/*========================================================================== + = = + = Private methods = + = = + ==========================================================================*/ + +/***********************************************************************//** + * @brief Initialise class members + ***************************************************************************/ +void ctedispcube::init_members(void) +{ + // Initialise members + m_outcube.clear(); +// m_apply_edisp = false; + + // Initialise protected members + m_obs.clear(); + m_edispcube.clear(); + + // Return + return; +} + + +/***********************************************************************//** + * @brief Copy class members + * + * @param[in] app Application. + ***************************************************************************/ +void ctedispcube::copy_members(const ctedispcube& app) +{ + // Copy attributes + m_outcube = app.m_outcube; +// m_apply_edisp = app.m_apply_edisp; + + // Copy protected members + m_obs = app.m_obs; + m_edispcube = app.m_edispcube; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Delete class members + ***************************************************************************/ +void ctedispcube::free_members(void) +{ + // 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 ctedispcube::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); + + // Throw exception if counts cube is given + require_inobs_nocube(G_GET_PARAMETERS); + + // Build observation container + m_obs = get_observations(); + + } // endif: there was no observation in the container + + // Get the incube filename + std::string incube = (*this)["incube"].filename(); + + // Get additional binning parameters + double amax = (*this)["amax"].real(); + int anumbins = (*this)["anumbins"].integer(); + + // Check for filename validity + if ((incube == "NONE") || (gammalib::strip_whitespace(incube) == "")) { + + // Create an event cube based on task parameters + GCTAEventCube cube = create_cube(m_obs); + + // Define edisp cube + m_edispcube = GCTACubeEdisp(cube, amax, anumbins); + + } + + // ... otherwise setup the exposure cube from the counts map + else { + + // Load event cube from filename + GCTAEventCube cube(incube); + + // Define edisp cube + m_edispcube = GCTACubeEdisp(cube, amax, anumbins); + + } // endelse: cube loaded from file + + // Read energy dispersion flag +// m_apply_edisp = (*this)["edisp"].boolean(); + + // Read output filename (if needed) + if (read_ahead()) { + m_outcube = (*this)["outcube"].filename(); + } + + // Return + return; +} diff --git a/src/ctedispcube/ctedispcube.hpp b/src/ctedispcube/ctedispcube.hpp new file mode 100644 index 0000000..da663c2 --- /dev/null +++ b/src/ctedispcube/ctedispcube.hpp @@ -0,0 +1,92 @@ +/*************************************************************************** + * ctedispcube - EDISP cube generation tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2014-2016 by Chia-Chun Lu * + * ----------------------------------------------------------------------- * + * * + * 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 ctedispcube.hpp + * @brief EDISP cube generation tool definition + * @author Chia-Chun Lu + */ + +#ifndef CTEDISPCUBE_HPP +#define CTEDISPCUBE_HPP + +/* __ Includes ___________________________________________________________ */ +#include "GammaLib.hpp" +#include "GCTALib.hpp" +#include "ctool.hpp" + +/* __Definitions _________________________________________________________ */ +#define CTEDISPCUBE_NAME "ctedispcube" +#define CTEDISPCUBE_VERSION "1.1.0" + + +/***********************************************************************//** + * @class ctedispcube + * + * @brief EDISP cube generation tool + ***************************************************************************/ +class ctedispcube : public ctool { + +public: + // Constructors and destructors + ctedispcube(void); + explicit ctedispcube(const GObservations& obs); + ctedispcube(int argc, char *argv[]); + ctedispcube(const ctedispcube& app); + virtual ~ctedispcube(void); + + // Operators + ctedispcube& operator=(const ctedispcube& app); + + // Methods + void clear(void); + void run(void); + void save(void); + const GCTACubeEdisp& edispcube(void) const; + +protected: + // Protected methods + void init_members(void); + void copy_members(const ctedispcube& app); + void free_members(void); + void get_parameters(void); + + // User parameters + GFilename m_outcube; //!< Output EDISP cube file name + bool m_apply_edisp; //!< Apply energy dispersion? + + // Protected members + GObservations m_obs; //!< Observation container + GCTACubeEdisp m_edispcube; //!< EDISP cube +}; + + +/***********************************************************************//** + * @brief Return EDISP cube + * + * @return EDISP cube + ***************************************************************************/ +inline +const GCTACubeEdisp& ctedispcube::edispcube(void) const +{ + return (m_edispcube); +} + +#endif /* CTEDISPCUBE_HPP */ diff --git a/src/ctedispcube/ctedispcube.par b/src/ctedispcube/ctedispcube.par new file mode 100644 index 0000000..ab8961c --- /dev/null +++ b/src/ctedispcube/ctedispcube.par @@ -0,0 +1,61 @@ +############################################################################# +# # +# CTA Science Analysis Tools # +# # +# ------------------------------------------------------------------------- # +# # +# File: ctpsfcube.par # +# Version: 1.0.0 # +# # +# Author: Chia-Chun Lu # +# c.lu@mpi-hd.mpg.de # +# MPIK # +# # +# Purpose: Parameter file for the ctpsfcube application # +# # +############################################################################# +# +# Observational parameters +#=================== +inobs, f, a, "NONE",,, "Input event list or observation definition XML file" +incube, f, a, "NONE",,, "Input counts cube file to extract PSF cube definition" +caldb, s, a, "prod2",,, "Calibration database" +irf, s, a, "South_0.5h",,, "Instrument response function" +edisp, b, h, no,,, "Apply energy dispersion?" +outcube, f, a, "psfcube.fits",,, "Output PSF cube file" + +# +# Energy binning parameters +#========================== +ebinalg, s, h, "LOG", FILE|LIN|LOG,, "Algorithm for defining energy bins" +emin, r, a, 0.1,,, "Lower energy limit (TeV)" +emax, r, a, 100.0,,, "Upper energy limit (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 +#=========================== +usepnt, b, h, no,,, "Use pointing instead of xref/yref parameters?" +nxpix, i, a, 10,,, "Size of the X axis in pixels" +nypix, i, a, 10,,, "Size of the Y axis in pixels" +binsz, r, a, 1.0,,, "Image scale (in degrees/pixel)" +coordsys, s, a, "CEL", CEL|GAL,,"Coordinate system (CEL - celestial, GAL - galactic)" +xref, r, a, 83.63,0,360, "First coordinate of image center in degrees (RA or galactic l)" +yref, r, a, 22.01,-90,90, "Second coordinate of image center in degrees (DEC or galactic b)" +proj, s, a, "CAR", AIT|AZP|CAR|MER|MOL|STG|TAN,, "Projection method" + +# +# Angular separation binning parameters +#====================================== +amax, r, h, 0.3,,, "Upper bound of angular separation between true and measued photon direction (degrees)" +anumbins, i, h, 200,,, "Number of angular separation bins" + +# +# 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, "ctpsfcube.log",,, "Log filename" diff --git a/src/ctedispcube/main.cpp b/src/ctedispcube/main.cpp new file mode 100644 index 0000000..9304f2c --- /dev/null +++ b/src/ctedispcube/main.cpp @@ -0,0 +1,79 @@ +/*************************************************************************** + * ctedispcube - CTA EDISP cube tool * + * ----------------------------------------------------------------------- * + * copyright (C) 2016 by Maria Haupt * + * ----------------------------------------------------------------------- * + * * + * 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 ctedispcube/main.cpp + * @brief CTA EDISP cube tool main code + * @author Maria Haupt + */ + +/* __ Includes ___________________________________________________________ */ +#ifdef HAVE_CONFIG_H +#include +#endif +#include +#include "ctedispcube.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 ctedispcube application. It allocates a + * ctedispcube object and runs the application. + ***************************************************************************/ +int main (int argc, char *argv[]) +{ + // Initialise return code + int rc = 1; + + // Create instance of application + ctedispcube 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" + " ctedispcube. 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; +}