Unnormalize background cube values in ctbkgcube (#2685)

Jürgen Knödlseder
1 parent 8f475e42
Showing 4 changed files with 34 additions and 9 deletions
ChangeLog
NEWS
doc/source/admin/release_history/1.6.rst
src/ctbkgcube/ctbkgcube.cpp
@@ -3,6 +3,7 @@
         * Version 1.6.0 released
         ========================
  
+        Unnormalize background cube values in ctbkgcube (#2685)
         Account for varying energy thresholds in ctbin weight computation (#2685)
         Add ctbkgcube::cntcube() method (#2685)
         Avoid ctselect error in case that RoI radius is unchanged (#2682)
@@ -61,8 +61,11 @@ ctbkgcube - CTA background cube generation
 ------------------------------------------
 Add ctbkgcube::cntcube() method so that a counts cube be specified as
 instance before running the tool. This avoids storing the counts cube on
-disk if the couts cube definition should be used to generate a background
-cube (#2685).
+disk if the counts cube definition should be used to generate a background
+cube. In addition, the background cube values are now devided by the
+weight values of the counts cube since later, during model evaluation,
+the model is multiplied by the weight values. This makes sure that the
+normalisation is correct (#2685).
  
  
 ctbutterfly - Generate butterfly diagrams
@@ -36,7 +36,8 @@ Improvements
  
 * [`2685 <https://cta-redmine.irap.omp.eu/issues/2685>`_] -
   Add ``ctbkgcube::cntcube()`` method
-  Accout for varying energy thresholds in :ref:`ctbin` weight computation
+  Account for varying energy thresholds in :ref:`ctbin` weight computation
+  Unnormalize background cube values in :ref:`ctbkgcube`
 * Add ``use_irf_bkg`` parameter to :ref:`cslightcrv`
 * [`2671 <https://cta-redmine.irap.omp.eu/issues/2671>`_] -
   Do not skip bins with zero counts in :ref:`csspec`
@@ -296,6 +296,18 @@ void ctbkgcube::run(void)
     // Fill background cube from observations
     m_background.fill(m_obs, logger);
  
+    // Un-normalize background cube. We do this since the counts cube weight is
+    // later multiplied when computing the model value for binned analysis.
+    GSkyMap* bkg = const_cast<GSkyMap*>(&m_background.cube());
+    for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
+        for (int iebin = 0; iebin < m_cube.ebins(); ++iebin) {
+            double weight = m_cube.weights()(pixel, iebin);
+            if (weight > 0.0) {
+                (*bkg)(pixel, iebin) /= weight;
+            }
+        }
+    }
+
     // Create a background model for the output background cube and append
     // that model to the input model in place of the original
     // background models
@@ -432,6 +444,7 @@ void ctbkgcube::init_members(void)
  
     // Initialise protected members
     m_background.clear();
+    m_cube.clear();
     m_bkgmdl.clear();
     m_outmdl.clear();
  
@@ -456,6 +469,7 @@ void ctbkgcube::copy_members(const ctbkgcube&amp; app)
  
     // Copy protected members
     m_background = app.m_background;
+    m_cube       = app.m_cube;
     m_bkgmdl     = app.m_bkgmdl;
     m_outmdl     = app.m_outmdl;
  
@@ -486,15 +500,21 @@ void ctbkgcube::get_parameters(void)
     // event lists, but do not accept counts cubes.
     setup_observations(m_obs, true, true, false);
  
-    // If the "incube" file name is valid then setup the background cube from
-    // the counts cube. Otherwise create a counts cube from the user
+    // If there is a valid event cube than use it for construction.
+    // Otherwise check if "incube" file name is valid then use the specified
+    // counts cube. Otherwise create a counts cube from the user
     // parameters
-    GCTAEventCube cube = (*this)["incube"].is_valid()
-                         ? GCTAEventCube((*this)["incube"].filename())
-                         : create_cube(m_obs);
+    if (m_cube.size() == 0) {
+        if ((*this)["incube"].is_valid()) {
+            m_cube = GCTAEventCube((*this)["incube"].filename());
+        }
+        else {
+            m_cube = create_cube(m_obs);
+        }
+    }
  
     // Define background cube
-    m_background = GCTACubeBackground(cube);
+    m_background = GCTACubeBackground(m_cube);
  
     // Setup models from "inmodel" parameter
     setup_models(m_obs);