Aligned ctbkgcube to ctexpcube and ctpsfcube by using the GCTACubeBackground::fill function

Michael Mayer
1 parent d063ccd1
Showing 3 changed files with 20 additions and 184 deletions
pyext/ctbkgcube.i
src/ctbkgcube/ctbkgcube.cpp
src/ctbkgcube/ctbkgcube.hpp
@@ -47,8 +47,7 @@ public:
     void                 clear(void);
     void                 run(void);
     void                 save(void);
-    const GCTACubeBackground& bgcube(void) const;
-    const GCTAEventCube& cube(void) const;
+    const GCTACubeBackground& cube(void) const;
     const GModels&       models(void) const;
 };
  
@@ -209,9 +209,6 @@ void ctbkgcube::run(void)
         log.cout(true);
     }
  
-    // Initialise run parameters
-    m_livetime = 0.0;
-
     // Get task parameters
     get_parameters();
  
@@ -248,9 +245,9 @@ void ctbkgcube::run(void)
  
     // Copy models from observation container and reset output model
     // container
+    GModels models_orig = m_obs.models();
     m_bkgmdl = m_obs.models();
     m_outmdl.clear();
-    m_cube_model = -1;
  
     // Remove all models that are not CTA background models from the
     // container and put all removed components in the output
@@ -309,59 +306,11 @@ void ctbkgcube::run(void)
         log.header1("Generate background 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 "+gammalib::str(i));
-                }
-            }
-
-            // Increment number of observations
-            n_observations++;
-
-            // Fill the cube
-            fill_cube(obs);
-
-        } // endif: CTA observation found
-
-    } // endfor: looped over observations
-
-    // Re-normalize cube to get units of counts/MeV/s/sr
-    if (m_livetime > 0.0) {
-    
-        // Loop over all bins in background cube and divide the content
-        // by the total livetime.
-        for (int i = 0; i < m_bkgcube.size(); ++i) {
-            GCTAEventBin* bin  = m_bkgcube[i];
-            double        rate = bin->counts() / m_livetime;
-            bin->counts(rate);
-        }
-
-    } // endif: livetime was positive
-
-    // Log results
-    if (logTerse()) {
-        log << gammalib::parformat("Number of observations");
-        log << n_observations << std::endl;
-    }
+    // Assign background models to container
+    m_obs.models(m_bkgmdl);
  
-    // Set cube background response
-    m_background.set(m_bkgcube.map(), m_bkgcube.ebounds());
+    // Fill background cube from observations
+    m_background.fill(m_obs);
  
     // Create a background model for the output background cube and append
     // that model to the input model in place of the original
@@ -391,6 +340,9 @@ void ctbkgcube::run(void)
         log << m_outmdl << std::endl;
     }
  
+    // Recover original models
+    m_obs.models(models_orig);
+
     // Return
     return;
 }
@@ -444,13 +396,10 @@ void ctbkgcube::init_members(void)
     m_outcube.clear();
     m_outmodel.clear();
     m_obs.clear();
-    m_bkgcube.clear();
     m_background.clear();
     m_bkgmdl.clear();
     m_outmdl.clear();
     m_ebounds.clear();
-    m_livetime   = 0.0;
-    m_cube_model = -1;
  
     // Return
     return;
@@ -468,13 +417,10 @@ void ctbkgcube::copy_members(const ctbkgcube&amp; app)
     m_outmodel   = app.m_outmodel;
     m_outcube    = app.m_outcube;
     m_obs        = app.m_obs;
-    m_bkgcube    = app.m_bkgcube;
     m_background = app.m_background;
     m_bkgmdl     = app.m_bkgmdl;
     m_outmdl     = app.m_outmdl;
     m_ebounds    = app.m_ebounds;
-    m_livetime   = app.m_livetime;
-    m_cube_model = app.m_cube_model;
  
     // Return
     return;
@@ -522,23 +468,18 @@ void ctbkgcube::get_parameters(void)
         (gammalib::strip_whitespace(incube) == "")) {
  
         // Create an event cube based on task parameters
-        m_bkgcube = create_cube(m_obs);
+        GCTAEventCube cube = create_cube(m_obs);
+
+        // Define background cube
+        m_background = GCTACubeBackground(cube);
  
     } // endif: filename was not valid
  
     // ... otherwise setup the background cube from the counts map
     else {
-    
-        // Load event cube from filename
-        m_bkgcube.load(incube);
-        // Loop over all bins in background cube
-        for (int i = 0; i < m_bkgcube.size(); ++i) {
  
-            // Get event bin
-            GCTAEventBin* bin = m_bkgcube[i];
-            bin->counts(0.0);
-
-        }
+        // Define background cube from file
+        m_background = GCTACubeBackground(incube);
  
     } // endelse: cube was loaded from file
  
@@ -557,8 +498,8 @@ void ctbkgcube::get_parameters(void)
         m_obs.models(models);
     }
  
-    // Get energy definition
-    m_ebounds = m_bkgcube.ebounds();
+//    // Get energy definition
+//    m_ebounds = m_background.ebounds();
  
     // Read output filenames (if needed)
     if (read_ahead()) {
@@ -569,92 +510,3 @@ void ctbkgcube::get_parameters(void)
     // Return
     return;
 }
-
-/***********************************************************************//**
- * @brief Generate background cube
- *
- * @param[in] obs Pointer to CTA observation.
- *
- * Fills the background cube with the model value for a CTA observation.
- * This method will not need to load the event file as it gets the event
- * ROI and GTI directly from the observation.
- ***************************************************************************/
-void ctbkgcube::fill_cube(GCTAObservation* obs)
-{
-
-    // Continue only if observation pointer is valid
-    if (obs != NULL) {
-
-        // Initialise statistics
-        double sum            = 0.0;
-        int    n_bins_inside  = 0;
-        int    n_bins_outside = 0;
-
-        // Extract region of interest from CTA observation
-        GCTARoi roi = obs->roi();
-
-        // Set GTI of actual observations as the GTI of the event cube
-        m_bkgcube.gti(obs->gti());
-
-        // Get observation livetime
-        double livetime = obs->livetime();
-
-        // Loop over all bins in background cube
-        for (int i = 0; i < m_bkgcube.size(); ++i) {
-
-            // Get event bin
-            GCTAEventBin* bin = m_bkgcube[i];
-
-            // Continue only if binned in contained in ROI
-            if (roi.contains(*bin)) {
-            
-                // Compute model value for event bin. The model value is
-                // given in counts/MeV/s/sr.
-                double model = m_bkgmdl.eval(*bin, *obs);
-
-                // Compute number events
-                sum += model * bin->size();
-
-                // Multiply by livetime to get the correct weighting for
-                // each observation. We divide by the total livetime later
-                // to get the background model in units of counts/MeV/s/sr.
-                model *= livetime;
-
-                // Add existing number of counts
-                model += bin->counts();
-
-                // Store cumulated value (units: counts/MeV/sr)
-                bin->counts(model);
-
-                // Increment contained bin counter
-                n_bins_inside++;
-                
-            } // endif: bin was contained in RoI
-            
-            // ... otherwise count bin as an outsider
-            else {
-                n_bins_outside++;
-            }
-
-        } // endfor: looped over all bins
-
-        // Accumulate livetime
-        m_livetime += livetime;
-
-        // Log results
-        if (logTerse()) {
-            log << gammalib::parformat("Bins within RoI");
-            log << n_bins_inside << std::endl;
-            log << gammalib::parformat("Bins outside RoI");
-            log << n_bins_outside << std::endl;
-            log << gammalib::parformat("Background events in cube");
-            log << sum << std::endl;
-            log << gammalib::parformat("Cube livetime");
-            log << livetime << " sec" << std::endl;
-        }
-
-    } // endif: observation pointer was not valid
-
-    // Return
-    return;
-}
@@ -59,8 +59,7 @@ public:
     void                 clear(void);
     void                 run(void);
     void                 save(void);
-    const GCTACubeBackground& bgcube(void) const;
-    const GCTAEventCube& cube(void) const;
+    const GCTACubeBackground& cube(void) const;
     const GModels&       models(void) const;
  
 protected:
@@ -75,27 +74,13 @@ protected:
     std::string   m_outcube;     //!< Filename of output cube
     std::string   m_outmodel;    //!< Filename of output XML model
     GObservations m_obs;         //!< Observation container
-    GCTAEventCube m_bkgcube;     //!< Background event cube
     GCTACubeBackground m_background;     //!< Background cube response
     GModels       m_bkgmdl;      //!< CTA background models
     GModels       m_outmdl;      //!< Output models
     GEbounds      m_ebounds;     //!< Energy boundaries
-    double        m_livetime;    //!< Total livetime for normalization
-    int           m_cube_model;  //!< Cube model index
-};
  
+};
  
-/***********************************************************************//**
- * @brief Return event cube containing background rate
- *
- * @return Background event cube.
- ***************************************************************************/
-inline
-const GCTAEventCube& ctbkgcube::cube(void) const
-{
-    // Return background cube
-    return (m_bkgcube);
-}
  
 /***********************************************************************//**
  * @brief Return background response cube containing background rate
@@ -103,7 +88,7 @@ const GCTAEventCube&amp; ctbkgcube::cube(void) const
  * @return Background response cube.
  ***************************************************************************/
 inline
-const GCTACubeBackground& ctbkgcube::bgcube(void) const
+const GCTACubeBackground& ctbkgcube::cube(void) const
 {
     // Return background cube
     return (m_background);