new ctool ctedispcube computing energy dispersion cube (#1724)

Maria Haupt
1 parent 30daf3a5
Showing 7 changed files with 674 additions and 0 deletions
configure.ac
src/Makefile.am
src/ctedispcube/Makefile.am
src/ctedispcube/ctedispcube.cpp
src/ctedispcube/ctedispcube.hpp
src/ctedispcube/ctedispcube.par
src/ctedispcube/main.cpp
@@ -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
@@ -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 \
+## 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
+/***************************************************************************
+ *                  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 <http://www.gnu.org/licenses/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file ctedispcube.cpp
+ * @brief EDISP cube generation tool implementation
+ * @author Maria Haupt
+ */
+
+/* __ Includes ___________________________________________________________ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include <cstdio>
+#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<bool> save_edisp;
+//    save_edisp.assign(m_obs.size(), false);
+//    for (int i = 0; i < m_obs.size(); ++i) {
+//        GCTAObservation* obs = dynamic_cast<GCTAObservation*>(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<GCTAObservation*>(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;
+}
+/***************************************************************************
+ *                  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 <http://www.gnu.org/licenses/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @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 */
+#############################################################################
+#                                                                           #
+#                       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"
+/***************************************************************************
+ *                    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 <http://www.gnu.org/licenses/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file ctedispcube/main.cpp
+ * @brief CTA EDISP cube tool main code
+ * @author Maria Haupt
+ */
+
+/* __ Includes ___________________________________________________________ */
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+#include <cstdio>
+#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;
+}