Re-enable ctmodel model computation for stacked analysis.

Jürgen Knödlseder
1 parent ac17570c
Showing 5 changed files with 138 additions and 6 deletions
ChangeLog
NEWS
src/ctmodel/ctmodel.cpp
src/ctmodel/ctmodel.hpp
src/ctmodel/ctmodel.par
@@ -3,6 +3,8 @@
         * Bug fix version 1.0.1 released
         ================================
  
+        Re-enable ctmodel computation for stacked analysis (#1609)
+
  
 2015-12-18  Juergen Knoedlseder  <jurgen.knodlseder@irap.omp.eu>
  
@@ -59,7 +59,8 @@ None
  
 11. ctmodel - CTA counts cube modelling
 ---------------------------------------
-None
+The get_obs() method, that has been accidently removed, has been put back
+which allows computation of a model cube for a stacked analysis.
  
  
 12. ctexpcube - CTA exposure cube generation
@@ -411,7 +411,7 @@ void ctmodel::init_members(void)
     m_gti.clear();
     m_has_cube    = false;
     m_append_cube = false;
-    m_binned = false;
+    m_binned      = false;
  
     // Return
     return;
@@ -435,7 +435,7 @@ void ctmodel::copy_members(const ctmodel&amp; app)
     m_gti         = app.m_gti;
     m_has_cube    = app.m_has_cube;
     m_append_cube = app.m_append_cube;
-    m_binned    = app.m_binned;
+    m_binned      = app.m_binned;
  
     // Return
     return;
@@ -466,7 +466,7 @@ void ctmodel::get_parameters(void)
     // If there are no observations in container then load them via user
     // parameters
     if (m_obs.size() == 0) {
-        m_obs = get_observations();
+        get_obs();
     }
  
     // Check if we got excactly one binned CTA observation
@@ -574,6 +574,134 @@ void ctmodel::get_parameters(void)
  
  
 /***********************************************************************//**
+ * @brief Get observation container
+ *
+ * Get an observation container according to the user parameters. The method
+ * supports loading of a individual FITS file or an observation definition
+ * file in XML format.
+ *
+ * If the input filename is empty, the method checks for the existence of the
+ * "expcube", "psfcube" and "bkgcube" parameters. If file names have been
+ * specified, the method loads the files and creates a dummy events cube that
+ * is appended to the observation container.
+ *
+ * If no file names are specified for the "expcube", "psfcube" or "bkgcube"
+ * parameters, the method reads the necessary parameters to build a CTA
+ * observation from scratch.
+ ***************************************************************************/
+void ctmodel::get_obs(void)
+{
+    // Get the filename from the input parameters
+    std::string filename = (*this)["inobs"].filename();
+
+    // If no observation definition file has been specified then read all
+    // parameters that are necessary to create an observation from scratch
+    if ((filename == "NONE") || (gammalib::strip_whitespace(filename) == "")) {
+
+        // Get response cube filenames
+        std::string expcube = (*this)["expcube"].filename();
+        std::string psfcube = (*this)["psfcube"].filename();
+        std::string bkgcube = (*this)["bkgcube"].filename();
+
+        // If the filenames are valid then build an observation from cube
+        // response information
+        if ((expcube != "NONE") && (psfcube != "NONE") && (bkgcube != "NONE") &&
+            (gammalib::strip_whitespace(expcube) != "") &&
+            (gammalib::strip_whitespace(psfcube) != "") &&
+            (gammalib::strip_whitespace(bkgcube) != "")) {
+
+            // Get exposure, PSF and background cubes
+            GCTACubeExposure   exposure(expcube);
+            GCTACubePsf        psf(psfcube);
+            GCTACubeBackground background(bkgcube);
+
+            // Create energy boundaries
+            GEbounds ebounds = create_ebounds();
+
+            // Create dummy sky map cube
+            GSkyMap map("CAR","GAL",0.0,0.0,1.0,1.0,1,1,ebounds.size());
+            m_append_cube = true;
+
+            // Create event cube
+            GCTAEventCube cube(map, ebounds, exposure.gti());
+
+            // Create CTA observation
+            GCTAObservation cta;
+            cta.events(cube);
+            cta.response(exposure, psf, background);
+
+            // Append observation to container
+            m_obs.append(cta);
+
+        } // endif: cube response information was available
+
+        // ... otherwise build an observation from IRF response information
+        else {
+
+            // Create CTA observation
+            GCTAObservation cta = create_cta_obs();
+
+            // Set response
+            set_obs_response(&cta);
+
+            // Append observation to container
+            m_obs.append(cta);
+            
+        }
+
+    } // endif: filename was "NONE" or ""
+
+    // ... otherwise we have a file name
+    else {
+
+        // If file is a FITS file then create an empty CTA observation
+        // and load file into observation
+        if (gammalib::is_fits(filename)) {
+
+            // Allocate empty CTA observation
+            GCTAObservation cta;
+
+            // Load data
+            cta.load(filename);
+
+            // Set response
+            set_obs_response(&cta);
+
+            // Append observation to container
+            m_obs.append(cta);
+
+            // Signal that no XML file should be used for storage
+            m_use_xml = false;
+
+        }
+
+        // ... otherwise load file into observation container
+        else {
+
+            // Load observations from XML file
+            m_obs.load(filename);
+
+            // For all observations that have no response, set the response
+            // from the task parameters
+            set_response(m_obs);
+
+            // Set observation boundary parameters (emin, emax, rad)
+            set_obs_bounds(m_obs);
+
+            // Signal that XML file should be used for storage
+            m_use_xml = true;
+
+        } // endelse: file was an XML file
+
+    }
+
+    // Return
+    return;
+
+}
+
+
+/***********************************************************************//**
  * @brief Fill model into model cube
  *
  * @param[in] obs CTA observation.
@@ -34,7 +34,7 @@
  
 /* __Definitions _________________________________________________________ */
 #define CTMODEL_NAME    "ctmodel"
-#define CTMODEL_VERSION "1.0.0"
+#define CTMODEL_VERSION "1.0.1"
  
  
 /***********************************************************************//**
@@ -79,6 +79,7 @@ protected:
     void copy_members(const ctmodel& app);
     void free_members(void);
     void get_parameters(void);
+    void get_obs(void);
     void fill_cube(const GCTAObservation* obs);
     bool has_cube(void) const;
  
@@ -5,7 +5,7 @@
 # ------------------------------------------------------------------------- #
 #                                                                           #
 #  File:      ctmodel.par                                                   #
-#  Version:   1.0.0                                                         #
+#  Version:   1.0.1                                                         #
 #                                                                           #
 #  Author:    Juergen Knoedlseder                                           #
 #             Jurgen.Knodlseder@irap.omp.eu                                 #