Complete ctexpcube tools and add unit test.

The ctexpcube tool is now functional, a unit test has been added.

Complete ctexpcube tools and add unit test.
The ctexpcube tool is now functional, a unit test has been added.
Jürgen Knödlseder
1 parent 8dab5cb9
Showing 7 changed files with 302 additions and 55 deletions
pyext/ctexpcube.i
pyext/ctools/ctools.i
src/ctexpcube/ctexpcube
src/ctexpcube/ctexpcube.cpp
src/ctexpcube/ctexpcube.hpp
src/ctexpcube/ctexpcube.par
test/test_ctools.sh
+/***************************************************************************
+ *                   ctexpcube - CTA exposure cube tool                    *
+ * ----------------------------------------------------------------------- *
+ *  copyright (C) 2014 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 <http://www.gnu.org/licenses/>.  *
+ *                                                                         *
+ ***************************************************************************/
+/**
+ * @file ctexpcube.i
+ * @brief CTA exposure cube tool definition
+ * @author Juergen Knoedlseder
+ */
+%{
+/* Put headers and other declarations here that are needed for compilation */
+#include "ctexpcube.hpp"
+#include "GTools.hpp"
+%}
+
+
+/***********************************************************************//**
+ * @class ctexpcube
+ *
+ * @brief CTA exposure cube tool
+ ***************************************************************************/
+class ctexpcube : public GApplication  {
+
+public:
+    // Constructors and destructors
+    ctexpcube(void);
+    explicit ctexpcube(GObservations obs);
+    ctexpcube(int argc, char *argv[]);
+    ctexpcube(const ctexpcube& app);
+    virtual ~ctexpcube(void);
+
+    // Operators
+    ctexpcube& operator=(const ctexpcube& app);
+
+    // Methods
+    void           clear(void);
+    void           execute(void);
+    void           run(void);
+    void           save(void);
+    GObservations& obs(void);
+    void           get_parameters(void);
+    void           get_obs(void);
+    void           set_response(void);
+    void           get_ebounds(void);
+    void           set_from_cntmap(const std::string& filename);
+};
+
+
+/***********************************************************************//**
+ * @brief CTA exposure map tool Python extension
+ ***************************************************************************/
+%extend ctexpcube {
+    ctexpcube copy() {
+        return (*self);
+    }
+}
 /***************************************************************************
  *                          ctools - SWIG file                             *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010-2013 by Juergen Knoedlseder                         *
+ *  copyright (C) 2010-2014 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -51,3 +51,4 @@
 %include "ctmodel.i"
 %include "ctselect.i"
 %include "ctskymap.i"
+%include "ctexpcube.i"
@@ -31,13 +31,12 @@
 #include <cstdio>
 #include "ctexpcube.hpp"
 #include "GTools.hpp"
+#include "GWcs.hpp"
  
 /* __ Method name definitions ____________________________________________ */
 #define G_GET_EBOUNDS                              "ctexpcube::get_ebounds()"
-
-
-#define G_SETUP_OBS                                    "ctexpcube::setup_obs()"
-#define G_MODEL_MAP                    "ctexpcube::model_map(GCTAObservation*)"
+#define G_GET_PARAMETERS                        "ctexpcube::get_parameters()"
+#define G_SET_FROM_CNTMAP          "ctexpcube::set_from_cntmap(std::string&)"
  
 /* __ Debug definitions __________________________________________________ */
  
@@ -276,37 +275,8 @@ void ctexpcube::run(void)
         log.header1("Generate exposure cube");
     }
  
-    // Initialise observation counter
-    int n_observations = 0;
-
-    // Loop over all observations in the container
-    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) {
-
-            // 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++;
-
-            // Generate exposure cube
-            //model_map(obs, m_obs.models());
-
-        } // endif: CTA observation found
-
-    } // endfor: looped over observations
+    // Fill exposure
+    m_expcube.fill(m_obs);
  
     // Restore energy dispersion flag for all CTA observations
     for (int i = 0; i < m_obs.size(); ++i) {
@@ -349,24 +319,27 @@ void ctexpcube::save(void)
  *
  * Get all task parameters from parameter file or (if required) by querying
  * the user. The parameters are read in the correct order.
+ *
+ * @todo Setup exposure cube from counts map
  ***************************************************************************/
 void ctexpcube::get_parameters(void)
 {
     // If we do not have any observations in the container then get an
     // input file name or observation descriptor file
     if (m_obs.size() == 0) {
-        m_infile = (*this)["infile"].filename();
+        get_obs();
     }
  
-    // Read response parameters
-    m_caldb       = (*this)["caldb"].string();
-    m_irf         = (*this)["irf"].string();
+    // Make sure that response is set
+    set_response();
+
+    // Read energy dispersion flag
     m_apply_edisp = (*this)["edisp"].boolean();
  
     // If no counts map is specified then setup the exposure cube from
     // the user parameters
-    m_cntmap = (*this)["cntmap"].filename();
-    if ((m_cntmap == "NONE") || (gammalib::strip_whitespace(m_cntmap) == "")) {
+    std::string cntmap = (*this)["cntmap"].filename();
+    if ((cntmap == "NONE") || (gammalib::strip_whitespace(cntmap) == "")) {
  
         // Get user parameters for counts map definition
         std::string wcs      = (*this)["proj"].string();
@@ -388,6 +361,10 @@ void ctexpcube::get_parameters(void)
  
     // ... otherwise setup the exposure cube from the counts map
     else {
+    
+        // Set exposure cube from counts map
+        set_from_cntmap(cntmap);
+    
     }
  
     // Read output filename (if needed)
@@ -401,8 +378,98 @@ void ctexpcube::get_parameters(void)
  
  
 /***********************************************************************//**
+ * @brief Get observation definition
+ *
+ * Get observation definition from the user parameters.
+ ***************************************************************************/
+void ctexpcube::get_obs(void)
+{
+    // Get input filename
+    std::string filename = (*this)["infile"].filename();
+
+    // Try first to open as FITS file
+    try {
+
+        // Allocate CTA observation
+        GCTAObservation obs;
+        
+        // Load input file in CTA observation
+        obs.load(filename);
+
+        // Append CTA observation to container
+        m_obs.append(obs);
+
+            
+    }
+        
+    // ... otherwise try to open as XML file
+    catch (GException::fits_open_error &e) {
+
+        // Load observations from XML file
+        m_obs.load(filename);
+
+    }
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @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) {
+
+            // Set response if it isn't set already
+            if (obs->response().aeff() == NULL) {
+
+                // Load response information
+                std::string database = (*this)["caldb"].string();
+                std::string irf      = (*this)["irf"].string();
+
+                // 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(database)) {
+                    caldb.rootdir(database);
+                }
+                else {
+                    caldb.open("cta", database);
+                }
+
+                // Set reponse
+            	obs->response(irf, caldb);
+
+            } // endif: observation already has a response
+
+        } // endif: observation was a CTA observation
+
+    } // endfor: looped over all observations
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
  * @brief Get energy boundaries from parameters
  *
+ * @exception GException::invalid_value
+ *            Invalid extension name encountered.
+ *
  * Get the energy boundaries from the user parameters.
  ***************************************************************************/
 void ctexpcube::get_ebounds(void)
@@ -467,6 +534,74 @@ void ctexpcube::get_ebounds(void)
 }
  
  
+/***********************************************************************//**
+ * @brief Set exposure cube definition from counts map
+ *
+ * @exception GException::invalid_value
+ *            Invalid counts map projection or invalid events encountered.
+ *
+ * Set exposure cube definition from counts map.
+ ***************************************************************************/
+void ctexpcube::set_from_cntmap(const std::string& filename)
+{
+    // Allocate CTA observation
+    GCTAObservation obs;
+        
+    // Load counts map in CTA observation
+    obs.load(filename);
+
+    // Set exposure cube from counts map
+    const GCTAEventCube* cube = dynamic_cast<const GCTAEventCube*>(obs.events());
+
+    // Continue only if cube is valid
+    if (cube != NULL) {
+
+        // Get sky map projection
+        const GWcs* wcs = dynamic_cast<const GWcs*>(cube->map().projection());
+        
+        // Continue only if projection is valid
+        if (wcs != NULL) {
+            
+            // Get user parameters for counts map definition
+            std::string proj     = wcs->code();
+            std::string coordsys = wcs->coordsys();
+            double      xref     = wcs->crval(0);
+            double      yref     = wcs->crval(1);
+            double      dx       = wcs->cdelt(0);
+            double      dy       = wcs->cdelt(1);
+            int         nx       = cube->map().nx();
+            int         ny       = cube->map().ny();
+
+            // Get energy definition
+            m_ebounds = cube->ebounds();
+
+            // Define exposure cube
+            m_expcube = GCTAExposure(proj, coordsys, xref, yref,
+                                     dx, dy, nx, ny,
+                                     m_ebounds);
+        
+        } // endif: WCS projection was valid
+
+        // ... projection is not of WCS type
+        else {
+            std::string msg = "Counts map project is not of WCS type.";
+            throw GException::invalid_value(G_SET_FROM_CNTMAP, msg);
+        }
+
+    } // endif: observation contained an events cube
+
+    // ... there is not events cube
+    else {
+        std::string msg = "No events cube found in file \""
+                          ""+filename+"\".";
+        throw GException::invalid_value(G_SET_FROM_CNTMAP, msg);
+    }
+
+    // Return
+    return;
+}
+
+
 /*==========================================================================
  =                                                                         =
  =                             Private methods                             =
@@ -479,10 +614,7 @@ void ctexpcube::get_ebounds(void)
 void ctexpcube::init_members(void)
 {
     // Initialise members
-    m_infile.clear();
     m_outfile.clear();
-    m_caldb.clear();
-    m_irf.clear();
     m_apply_edisp = false;
  
     // Initialise protected members
@@ -507,10 +639,7 @@ void ctexpcube::init_members(void)
 void ctexpcube::copy_members(const ctexpcube& app)
 {
     // Copy attributes
-    m_infile      = app.m_infile;
     m_outfile     = app.m_outfile;
-    m_caldb       = app.m_caldb;
-    m_irf         = app.m_irf;
     m_apply_edisp = app.m_apply_edisp;
  
     // Copy protected members
@@ -61,7 +61,10 @@ public:
     void           save(void);
     GObservations& obs(void) { return m_obs; }
     void           get_parameters(void);
+    void           get_obs(void);
+    void           set_response(void);
     void           get_ebounds(void);
+    void           set_from_cntmap(const std::string& filename);
  
 protected:
     // Protected methods
@@ -70,11 +73,7 @@ protected:
     void           free_members(void);
  
     // User parameters
-    std::string   m_infile;      //!< Input counts map or XML file
-    std::string   m_cntmap;      //!< Counts map for exposure cube definition
     std::string   m_outfile;     //!< Output exposure cube file
-    std::string   m_caldb;       //!< Calibration database
-    std::string   m_irf;         //!< Instrument response functions
     bool          m_apply_edisp; //!< Apply energy dispersion?
  
     // Protected members
@@ -17,7 +17,7 @@
 #
 # General parameters
 #===================
-infile,   f, a, "NONE",,, "Input event list, counts map or observation definition file"
+infile,   f, a, "NONE",,, "Observation definition file"
 cntmap,   f, a, "NONE",,, "Counts map for exposure cube definition"
 outfile,  f, a, "expcube.fits",,, "Output exposure cube file"
 caldb,    s, a, "dummy",,, "Calibration database"
@@ -2,7 +2,7 @@
 # ==========================================================================
 # This script tests all ctools
 #
-# Copyright (C) 2011-2012 Juergen Knoedlseder
+# Copyright (C) 2011-2014 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
@@ -40,6 +40,7 @@ ctmodel=../src/ctmodel/ctmodel
 ctobssim=../src/ctobssim/ctobssim
 ctselect=../src/ctselect/ctselect
 ctskymap=../src/ctskymap/ctskymap
+ctexpcube=../src/ctexpcube/ctexpcube
  
  
 #
@@ -340,4 +341,49 @@ else
   exit 1
 fi
 $ECHO " ok"
-#sleep 100
+
+
+#
+# Test ctexpcube
+# ==============
+$ECHO -n "Test ctexpcube: "
+#
+# Run 1
+$ctexpcube infile="data/crab_events.fits" \
+           cntmap="NONE" \
+           outfile="expcube1.fits" \
+           caldb="irf" \
+           irf="cta_dummy_irf" \
+           emin=0.1 \
+           emax=100.0 \
+           enumbins=20 \
+           nxpix=200 \
+           nypix=200 \
+           binsz=0.02 \
+           coordsys="CEL" \
+           xref=83.63 \
+           yref=22.01 \
+           proj="CAR"
+$ECHO -n "."
+if [ -s "expcube1.fits" ]
+then
+  $ECHO -n "."
+else
+  $ECHO " expcube1.fits file is not valid"
+  exit 1
+fi
+#
+# Run 2
+$ctexpcube infile="data/crab_events.fits" \
+           cntmap="data/crab_cntmap.fits" \
+           outfile="expcube2.fits" \
+           caldb="irf" \
+           irf="cta_dummy_irf"
+$ECHO -n "."
+if [ -s "expcube2.fits" ]
+then
+  $ECHO -n "."
+else
+  $ECHO " expcube2.fits file is not valid"
+  exit 1
+fi