Apply event bin weighting in ctmodel for stacked computation (#2721)

Jürgen Knödlseder
1 parent 1c81f76b
Showing 4 changed files with 20 additions and 2 deletions
ChangeLog
NEWS
doc/source/admin/release_history/1.6.rst
src/ctmodel/ctmodel.cpp
@@ -3,6 +3,7 @@
         * Version 1.6.0 released
         ========================
+        Apply event bin weighting in ctmodel for stacked computation (#2721)
         Add code generator (#2709)
         Add csbkgmodel script (#2694)
         Compute DETX and DETY in ctmodel if they are not available (#2694)
@@ -118,6 +118,13 @@ None
 ctmodel - CTA counts cube modelling
 -----------------------------------
+When running ctmodel on a counts cube, and specifically a stacked count cube,
+the event bin weighting is taken into account in the model evaluation. This
+is needed to properly take into account the changing thresholds. No such
+weighting is done if the model is computed from an event list, since in that
+case the model is properly computed for each observation using the run-wise
+response (#2721).
+
 The tools now computes DETX and DETY for every event bin if they are not
 available so that model that rely on the existence of these parameters can
 be handled (#2694).
@@ -21,6 +21,8 @@ In particular, this release provides:
 Bug fixes
 ---------
+* [`2721 <https://cta-redmine.irap.omp.eu/issues/2721>`_] -
+  Apply event bin weighting in :ref:`ctmodel` for stacked computation
 * [`2694 <https://cta-redmine.irap.omp.eu/issues/2694>`_] -
   Compute ``DETX`` and ``DETY`` in :ref:`ctmodel` if they are not available
 * [`2686 <https://cta-redmine.irap.omp.eu/issues/2686>`_] -
@@ -893,6 +893,9 @@ void ctmodel::fill_cube(const GCTAObservation* obs, GModels&amp; models)
     bin.ontime(obs->ontime());
     bin.weight(1.0);
+    // Flag if the CTA observation is a counts cube
+    bool obs_is_cube = (obs->eventtype() == "CountsCube");
+
     // Loop over all the spatial bins
     for (int i = 0; i < npix; ++i) {
@@ -923,11 +926,16 @@ void ctmodel::fill_cube(const GCTAObservation* obs, GModels&amp; models)
                 continue;
             }
-            // Set energy, energy width, energy index and weight of bin
+            // Set energy, energy width and energy index
             bin.energy(m_energy[iebin]);
             bin.ewidth(m_ewidth[iebin]);
             bin.ieng(iebin);
-            //bin.weight(m_cube[ibin]->weight()); // Bad residuals for H.E.S.S. Crab
+
+            // If observation is a counts cube then set the proper weight for
+            // this bin
+            if (obs_is_cube) {
+                bin.weight(m_cube[ibin]->weight());
+            }
             // Compute model value for cube bin
             double model = models.eval(bin, *obs) * bin.size();