Updated ctmodel (#2189)

Implemented a change to remove all models that dont overlap the current observation when filling the cube. This speeds up the computation when filling from a large number of models.

Updated ctmodel (#2189)
Implemented a change to remove all models that dont overlap the current observation when filling the cube. This speeds up the computation when filling from a large number of models.
Josh Cardenzana
1 parent ed0d4b4b
Showing 2 changed files with 50 additions and 2 deletions
src/ctmodel/ctmodel.cpp
src/ctmodel/ctmodel.hpp
@@ -820,7 +820,8 @@ void ctmodel::fill_cube(const GCTAObservation* obs)
  
     // Extract copy of models from observation container
     GModels models(m_obs.models());
-
+    trim_models(models, roi);
+    
     // Get pointer to event cube pixels
     double* pixels = const_cast<double*>(m_cube.counts().pixels());
  
@@ -905,3 +906,49 @@ void ctmodel::fill_cube(const GCTAObservation* obs)
     // Return
     return;
 }
+
+
+/***********************************************************************//**
+ * @brief Remove the models falling outside a defined region of interest
+ *
+ * @param[out] models           Model container to be trimmed
+ * @param[in]  roi              Observation region of interest
+ ***************************************************************************/
+void ctmodel::trim_models(GModels& models, const GCTARoi& roi)
+{
+    // Do nothing if roi is not valid. This will be the case for binned data
+    // sets as they have the default roi radius of 0. In this case we would
+    // remove ALL models, which is obviously incorrect. So we keep all models
+    // when filling based on a counts cube.
+    if (!roi.is_valid()) {
+        return;
+    }
+    
+    // Remove all models that dont overlap with the region of interest. Note
+    // that an extra factor of 0.5 is used since point sources have regions
+    // of radius 0, which is a problem when they fall just barely outside the
+    // ROI. This ensures point sources are appropriately included.
+    GSkyRegionCircle obsreg(
+        roi.centre().dir().ra_deg(),
+        roi.centre().dir().dec_deg(),
+        roi.radius()+0.5);
+    
+    // Loop over the models in the passed container
+    for (int mod=0; mod<models.size(); mod++) {
+        GModelSky* model = dynamic_cast<GModelSky*>(models.at(mod));
+        
+        if (model == NULL) {
+            // Model is most likely a background model and should be kept
+            continue;
+        } else if (!model->spatial()->region()->overlaps(obsreg)) {
+            // Model represents a source, does not overlap the observation
+            // and should be removed.
+            models.remove(mod);
+            
+            // Decrement the index variable to prevent skipping the next model
+            mod--;
+        }
+    }
+    
+    return;
+}
@@ -82,7 +82,8 @@ protected:
     void extract_cube_properties(void);
     void fill_cube(const GCTAObservation* obs);
     bool has_cube(void) const;
-
+    void trim_models(GModels& models, const GCTARoi& roi);
+    
     // User parameters
     GFilename m_outcube;      //!< Output model cube
     bool      m_apply_edisp;  //!< Apply energy dispersion?