ctmodel now derives from ctool. Add ctool::set_response() method that sets the r…

…esponse for all CTA observation in the container that have not yet a response using the task parameters.

ctmodel now derives from ctool. Add ctool::set_response() method that sets the r…
…esponse for all CTA observation in the container that have not yet a response using the task parameters.
Jürgen Knödlseder
1 parent 6ec18d45
Showing 9 changed files with 68 additions and 201 deletions
pyext/ctmodel.i
src/ctexpcube/ctexpcube.cpp
src/ctexpcube/ctexpcube.hpp
src/ctlike/ctlike.cpp
src/ctmodel/Makefile.am
src/ctmodel/ctmodel.cpp
src/ctmodel/ctmodel.hpp
src/support/ctool.cpp
src/support/ctool.hpp
 /***************************************************************************
- *                ctmodel - CTA model cube generation tool                 *
+ *                  ctmodel - Model cube generation tool                   *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2012-2014 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
@@ -20,23 +20,21 @@
  ***************************************************************************/
 /**
  * @file ctmodel.i
- * @brief CTA model cube tool definition
+ * @brief Model cube generation tool definition
  * @author Juergen Knoedlseder
  */
 %{
 /* Put headers and other declarations here that are needed for compilation */
 #include "ctmodel.hpp"
-#include "GTools.hpp"
 %}
 /***********************************************************************//**
  * @class ctmodel
  *
- * @brief CTA counts model tool interface defintion
+ * @brief Model cube generation tool
  ***************************************************************************/
-class ctmodel : public GApplication  {
-
+class ctmodel : public ctool {
 public:
     // Constructors and destructors
     ctmodel(void);
@@ -56,7 +54,7 @@ public:
 /***********************************************************************//**
- * @brief CTA counts model tool Python extension
+ * @brief Model cube generation tool Python extension
  ***************************************************************************/
 %extend ctmodel {
     ctmodel copy() {
@@ -394,8 +394,9 @@ void ctexpcube::get_parameters(void)
         get_obs();
     }
-    // Make sure that response is set
-    set_response();
+    // For all observations that have no response, set the response
+    // from the task parameters
+    set_response(m_obs);
     // Read energy dispersion flag
     m_apply_edisp = (*this)["edisp"].boolean();
@@ -480,37 +481,6 @@ void ctexpcube::get_obs(void)
 /***********************************************************************//**
- * @brief Set observation response
- *
- * Set response for all observations that so have no response.
- ***************************************************************************/
-void ctexpcube::set_response(void)
-{
-    // Loop over all observations
-    for (int i = 0; i < m_obs.size(); ++i) {
-
-        // Is this observation a CTA observation?
-        GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
-
-        // Yes ...
-        if (obs != NULL) {
-
-            // If the observation has no response, then set it from
-            // the task parameters
-            if (!obs->has_response()) {
-                set_obs_response(obs);
-            }
-
-        } // endif: observation was a CTA observation
-
-    } // endfor: looped over all observations
-
-    // Return
-    return;
-}
-
-
-/***********************************************************************//**
  * @brief Set exposure cube definition from counts map
  *
  * @exception GException::invalid_value
@@ -69,7 +69,6 @@ protected:
     void free_members(void);
     void get_parameters(void);
     void get_obs(void);
-    void set_response(void);
     void set_from_cntmap(const std::string& filename);
     // User parameters
@@ -415,26 +415,9 @@ void ctlike::get_parameters(void)
             // 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
-            for (int i = 0; i < m_obs.size(); ++i) {
-
-                // Get CTA observation
-                GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
-
-                // Continue only if observation is a CTA observation
-                if (obs != NULL) {
-
-                    // If response is not valid then set response from
-                    // task parameters
-                    if (!obs->has_response()) {
-                        set_obs_response(obs);
-                    }
-
-                } // endif: observation was a CTA observation
-
-            } // endfor: looped over observations
+            // For all observations that have no response, set the response
+            // from the task parameters
+            set_response(m_obs);
         } // endcatch: file was an XML file
 ## Process this file with automake to produce Makefile.in
+# Set include directories
+INCLUDES = -I$(top_srcdir)/src/support
+
 # Program name
 bin_PROGRAMS = ctmodel
@@ -14,7 +17,7 @@ libctmodel_la_SOURCES = $(sources)
 ctmodel_SOURCES = main.cpp
 ctmodel_LDFLAGS = $(LIB_LINK_OPTION) @LDFLAGS@
 ctmodel_CPPFLAGS = @CPPFLAGS@
-ctmodel_LDADD = libctmodel.la
+ctmodel_LDADD = libctmodel.la $(top_srcdir)/src/support/libsupport.la
 # Parameter files
 pfilesdir = $(prefix)/syspfiles
 /***************************************************************************
- *                ctmodel - CTA model cube generation tool                 *
+ *                  ctmodel - Model cube generation tool                   *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2012-2014 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
@@ -20,7 +20,7 @@
  ***************************************************************************/
 /**
  * @file ctmodel.cpp
- * @brief CTA model cube tool implementation
+ * @brief Model cube generation tool implementation
  * @author Juergen Knoedlseder
  */
@@ -33,7 +33,6 @@
 #include "GTools.hpp"
 /* __ Method name definitions ____________________________________________ */
-#define G_GET_EBOUNDS                                "ctmodel::get_ebounds()"
 /* __ Debug definitions __________________________________________________ */
@@ -49,14 +48,11 @@
 /***********************************************************************//**
  * @brief Void constructor
  ***************************************************************************/
-ctmodel::ctmodel(void) : GApplication(CTMODEL_NAME, CTMODEL_VERSION)
+ctmodel::ctmodel(void) : ctool(CTMODEL_NAME, CTMODEL_VERSION)
 {
     // Initialise members
     init_members();
-    // Write header into logger
-    log_header();
-
     // Return
     return;
 }
@@ -70,8 +66,7 @@ ctmodel::ctmodel(void) : GApplication(CTMODEL_NAME, CTMODEL_VERSION)
  * This method creates an instance of the class by copying an existing
  * observations container.
  ***************************************************************************/
-ctmodel::ctmodel(const GObservations& obs) :
-         GApplication(CTMODEL_NAME, CTMODEL_VERSION)
+ctmodel::ctmodel(const GObservations& obs) : ctool(CTMODEL_NAME, CTMODEL_VERSION)
 {
     // Initialise members
     init_members();
@@ -79,9 +74,6 @@ ctmodel::ctmodel(const GObservations&amp; obs) :
     // Set observations
     m_obs = obs;
-    // Write header into logger
-    log_header();
-
     // Return
     return;
 }
@@ -95,14 +87,11 @@ ctmodel::ctmodel(const GObservations&amp; obs) :
  * @param[in] argv Array of command line arguments.
  ***************************************************************************/
 ctmodel::ctmodel(int argc, char *argv[]) :
-         GApplication(CTMODEL_NAME, CTMODEL_VERSION, argc, argv)
+         ctool(CTMODEL_NAME, CTMODEL_VERSION, argc, argv)
 {
     // Initialise members
     init_members();
-    // Write header into logger
-    log_header();
-
     // Return
     return;
 }
@@ -113,7 +102,7 @@ ctmodel::ctmodel(int argc, char *argv[]) :
  *
  * @param[in] app Application.
  ***************************************************************************/
-ctmodel::ctmodel(const ctmodel& app) : GApplication(app)
+ctmodel::ctmodel(const ctmodel& app) : ctool(app)
 {
     // Initialise members
     init_members();
@@ -151,13 +140,13 @@ ctmodel::~ctmodel(void)
  * @param[in] app ctmodel application.
  * @return Returns ctmodel application.
  ***************************************************************************/
-ctmodel& ctmodel::operator= (const ctmodel& app)
+ctmodel& ctmodel::operator=(const ctmodel& app)
 {
     // Execute only if object is not identical
     if (this != &app) {
         // Copy base class members
-        this->GApplication::operator=(app);
+        this->ctool::operator=(app);
         // Free members
         free_members();
@@ -188,10 +177,12 @@ void ctmodel::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
@@ -454,11 +445,6 @@ void ctmodel::copy_members(const ctmodel&amp; app)
  ***************************************************************************/
 void ctmodel::free_members(void)
 {
-    // Write separator into logger
-    if (logTerse()) {
-        log << std::endl;
-    }
-
     // Return
     return;
 }
@@ -654,38 +640,9 @@ void ctmodel::setup_obs(void)
         m_obs.models(GModels(m_srcmdl));
     }
-    // Setup response for all observations
-    for (int i = 0; i < m_obs.size(); ++i) {
-
-        // Is this observation a CTA observation?
-        GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
-
-        // Yes ...
-        if (obs != NULL) {
-
-            // Set response if we don't have one
-            if (!obs->has_response()) {
-
-                // Set calibration database. If specified parameter is a
-                // directory then use this as the pathname to the calibration
-                // database. Otherwise interpret this as the instrument name,
-                // the mission being "cta"
-                GCaldb caldb;
-                if (gammalib::dir_exists(m_caldb)) {
-                    caldb.rootdir(m_caldb);
-                }
-                else {
-                    caldb.open("cta", m_caldb);
-                }
-
-                // Set reponse
-            	obs->response(m_irf, caldb);
-
-            } // endif: observation already has a response
-
-        } // endif: observation was a CTA observation
-
-    } // endfor: looped over all observations
+    // For all observations that have no response, set the response
+    // from the task parameters
+    set_response(m_obs);
     // Return
     return;
@@ -716,7 +673,7 @@ void ctmodel::init_cube(void)
         gti.append(GTime(0.0), GTime(0.1234));
         // Set energy boundaries
-        get_ebounds();
+        m_ebounds = get_ebounds();
         // Setup skymap
         GSkymap map = GSkymap(m_proj, m_coordsys,
@@ -824,80 +781,3 @@ void ctmodel::fill_cube(const GCTAObservation* obs)
     // Return
     return;
 }
-
-
-/***********************************************************************//**
- * @brief Get the energy boundaries
- *
- * @exception GException::invalid_value
- *            No valid energy boundary extension found.
- *
- * Get the energy boundaries according to the user parameters. The method
- * supports loading of energy boundary information from the EBOUNDS or
- * ENERGYBINS extension, or setting energy boundaries using a linear
- * or logarithmical spacing.
- *
- * @todo Ultimately this method should go in a support library as it is
- * used by several ctools.
- ***************************************************************************/
-void ctmodel::get_ebounds(void)
-{
-    // Initialse energy boundary information
-    m_ebounds.clear();
-
-    // Check whether energy binning information should be read from a FITS
-    // file ...
-    if (m_ebinalg == "FILE") {
-
-        // Open energy boundary file using the EBOUNDS or ENERGYBINS
-        // extension. Throw an exception if opening fails.
-        GFits file(m_ebinfile);
-        if (file.contains("EBOUNDS")) {
-            file.close();
-            m_ebounds.load(m_ebinfile,"EBOUNDS");
-        }
-        else if (file.contains("ENERGYBINS")) {
-            file.close();
-            m_ebounds.load(m_ebinfile,"ENERGYBINS");
-        }
-        else {
-            file.close();
-            std::string msg = "No extension with name \"EBOUNDS\" or"
-                              " \"ENERGYBINS\" found in FITS file"
-                              " \""+m_ebinfile+"\".\n"
-                              "An \"EBOUNDS\" or \"ENERGYBINS\" extension"
-                              " is required if the parameter \"ebinalg\""
-                              " is set to \"FILE\".";
-            throw GException::invalid_value(G_GET_EBOUNDS, msg);
-        }
-        
-        // Set enumbins parameter to number of ebounds
-        m_enumbins = m_ebounds.size();
-
-    } // endif: ebinalg was "FILE"
-
-    // ... otherwise use a linear or a logarithmically-spaced energy binning
-    else {
-
-        // Initialise log mode for ebinning
-        bool log = true;
-
-        // check if algorithm is linear
-        if (m_ebinalg == "LIN") {
-            log = false;
-        }
-
-        // todo: should we also check if m_ebinalg is "LOG"
-        // and throw an exception if neither LIN/LOG/FILE
-        // is given?
-
-        // Setup energy range covered by data
-        GEnergy  emin(m_emin, "TeV");
-        GEnergy  emax(m_emax, "TeV");
-        m_ebounds = GEbounds(m_enumbins, emin, emax, log);
-
-    } //endif: ebinalg was not "FILE"
-
-    // Return
-    return;
-}
 /***************************************************************************
- *                ctmodel - CTA model cube generation tool                 *
+ *                  ctmodel - Model cube generation tool                   *
  * ----------------------------------------------------------------------- *
  *  copyright (C) 2012-2014 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
@@ -20,7 +20,7 @@
  ***************************************************************************/
 /**
  * @file ctmodel.hpp
- * @brief CTA model cube tool definition
+ * @brief Model cube generation tool definition
  * @author Juergen Knoedlseder
  */
@@ -30,6 +30,7 @@
 /* __ Includes ___________________________________________________________ */
 #include "GammaLib.hpp"
 #include "GCTALib.hpp"
+#include "ctool.hpp"
 /* __Definitions _________________________________________________________ */
 #define CTMODEL_NAME    "ctmodel"
@@ -39,7 +40,7 @@
 /***********************************************************************//**
  * @class ctmodel
  *
- * @brief CTA counts model tool interface defintion
+ * @brief Model cube generation tool
  *
  * 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
@@ -50,7 +51,7 @@
  * in form of FITS files (model maps) and an updated observation definition
  * XML file.
  ***************************************************************************/
-class ctmodel : public GApplication  {
+class ctmodel : public ctool {
 public:
     // Constructors and destructors
@@ -80,7 +81,6 @@ protected:
     void setup_obs(void);
     void init_cube(void);
     void fill_cube(const GCTAObservation* obs);
-    void get_ebounds(void);
     // User parameters
     std::string m_infile;      //!< Input counts cube
@@ -303,6 +303,39 @@ GEbounds ctool::get_ebounds(void)
 /***********************************************************************//**
+ * @brief Set response for all CTA observations in container
+ *
+ * @param[in,out] obs Observation container
+ *
+ * Set the response for a CTA observations in the container that so far have
+ * no response using the "database" and "irf" task parameters.
+ ***************************************************************************/
+void ctool::set_response(GObservations& obs)
+{
+    // Setup response for all observations
+    for (int i = 0; i < obs.size(); ++i) {
+
+        // Is this observation a CTA observation?
+        GCTAObservation* cta = dynamic_cast<GCTAObservation*>(obs[i]);
+
+        // Yes ...
+        if (cta != NULL) {
+
+            // Set response if we don't have one
+            if (!cta->has_response()) {
+                set_obs_response(cta);
+            }
+
+        } // endif: observation was a CTA observation
+
+    } // endfor: looped over all observations
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Set response for CTA observation
  *
  * @param[in,out] obs CTA observation
@@ -60,6 +60,7 @@ protected:
     void     copy_members(const ctool& app);
     void     free_members(void);
     GEbounds get_ebounds(void);
+    void     set_response(GObservations& obs);
     void     set_obs_response(GCTAObservation* obs);
 };