Fixed response computations for composite model

Michael Mayer
1 parent 6d02c035
Showing 3 changed files with 24 additions and 9 deletions
inst/cta/src/GCTAResponseCube.cpp
inst/cta/src/GCTAResponseIrf.cpp
src/model/GModelSpatialComposite.cpp
@@ -1516,13 +1516,25 @@ double GCTAResponseCube::irf_composite(const GEvent&amp;       event,
     // Initialise IRF
     double irf = 0.0;
  
+    // Get pointer to composite model
     const GModelSpatialComposite* model = dynamic_cast<const GModelSpatialComposite*>(source.model());
  
+    // Loop over model components
     for (int i=0;i<model->components();++i) {
+
+        // Get pointer to spatial component
         GModelSpatial* spat = const_cast<GModelSpatial*>(model->component(i));
+
+        // Create new GSource object
         GSource src(source.name(), spat, source.energy(), source.time());
+
+        // Compute irf value
         irf += this->irf(event, source, obs);
     }
  
+    // Divide by number of model components
+    irf /= (double)model->components();
+
+
     return irf;
 }
@@ -2882,9 +2882,12 @@ double GCTAResponseIrf::irf_composite(const GEvent&amp;       event,
         GSource src(source.name(), spat, source.energy(), source.time());
  
         // Compute irf value
-        irf = this->irf(event, source, obs);
+        irf += this->irf(event, src, obs);
     }
  
+    // Divide by number of model components
+    irf /= (double)model->components();
+
     // Return IRF value
     return irf;
 }
@@ -3516,24 +3519,24 @@ double GCTAResponseIrf::nroi_composite(const GModelSky&amp;    model,
     double nroi = 0.0;
  
     // Get composite model
-    const GModelSpatialComposite* comp = dynamic_cast<const GModelSpatialComposite*>(model.spatial());
+    GModelSpatialComposite* comp = dynamic_cast<GModelSpatialComposite*>(model.spatial());
  
     // Loop over model components
-    for (int i=0;i<comp->components();++i) {
-
-        // Retrieve pointer to spatial component
-        GModelSpatial* spat = const_cast<GModelSpatial*>(comp->component(i));
+    for (int i = 0; i < comp->components(); ++i) {
  
         // Create new sky model
         GModelSky sky = GModelSky(model);
  
         // Assign spatial part
-        sky.spatial(spat);
+        sky.spatial(comp->component(i));
  
         // Compute nroi
-        nroi = this->nroi(sky, srcEng, srcTime, obsEng, obsTime, obs);
+        nroi += this->nroi(sky, srcEng, srcTime, obsEng, obsTime, obs);
     }
  
+    // Divide by number of model components
+    nroi /= (double)comp->components();
+
     // Return nroi
     return nroi;
  
@@ -628,7 +628,7 @@ const GModelSpatial* GModelSpatialComposite::component(const std::string&amp; name) 
 void GModelSpatialComposite::init_members(void)
 {
     // Initialise model type
-    m_type = "SpatialComposite";
+    m_type = "Composite";
  
     // Initialise models vector
     m_components.clear();