Merge branch '1513-bining-tools-consider-ebounds' into integration

Jürgen Knödlseder
2 parents 84c8f22c 0dae16ce
Showing 7 changed files with 76 additions and 5 deletions
include/GEbounds.hpp
inst/cta/src/GCTACubeBackground.cpp
inst/cta/src/GCTACubeExposure.cpp
inst/cta/src/GCTACubePsf.cpp
pyext/GEbounds.i
src/obs/GEbounds.cpp
test/test_GObservation.cpp
@@ -103,6 +103,7 @@ public:
     GEnergy        elogmean(const int& index) const;
     GEnergy        ewidth(const int& index) const;
     bool           contains(const GEnergy& eng) const;
+    bool           contains(const GEnergy& emin, const GEnergy& emax) const;
     std::string    print(const GChatter& chatter = NORMAL) const;
  
  
@@ -327,6 +327,9 @@ void GCTACubeBackground::fill(const GObservations&amp; obs, GLog* log)
         // Extract region of interest from CTA observation
         GCTARoi roi = cta->roi();
  
+        // Extract energy boundaries from CTA observation
+        GEbounds obs_ebounds = cta->ebounds();
+
         // Check for RoI sanity
         if (!roi.is_valid()) {
             std::string msg = "No RoI information found in input observation "
@@ -347,6 +350,15 @@ void GCTACubeBackground::fill(const GObservations&amp; obs, GLog* log)
             // Get event bin
             GCTAEventBin* bin = eventcube[i];
  
+            // Determine energy bin index. Skip if bin is not fully contained
+            // within the energy boundaries of the observation
+            int index = m_ebounds.index(bin->energy());
+            if (index == -1 || 
+                !obs_ebounds.contains(m_ebounds.emin(index),
+                                      m_ebounds.emax(index))) {
+                continue;
+            }
+
             // Continue only if binned in contained in ROI
             if (roi.contains(*bin)) {
  
@@ -423,6 +423,9 @@ void GCTACubeExposure::fill(const GObservations&amp; obs, GLog* log)
         // Extract region of interest from CTA observation
         GCTARoi roi = cta->roi();
  
+        // Extract energy boundaries from CTA observation
+        GEbounds obs_ebounds = cta->ebounds();
+
         // Check for RoI sanity
         if (!roi.is_valid()) {
             std::string msg = "No RoI information found in input observation "
@@ -473,12 +476,19 @@ void GCTACubeExposure::fill(const GObservations&amp; obs, GLog* log)
                 // Loop over all exposure cube energy bins
                 for (int iebin = 0; iebin < m_ebounds.size(); ++iebin){
  
-                    // Get logE/TeV
-                    double logE = m_ebounds.elogmean(iebin).log10TeV();
+                    // Only add exposure to pixel if energy bin is inside 
+                    // observation energy boundaries
+                    if (obs_ebounds.contains(m_ebounds.emin(iebin),
+                                             m_ebounds.emax(iebin))) {
+
+                        // Get logE/TeV
+                        double logE = m_ebounds.elogmean(iebin).log10TeV();
+
+                        // Add to exposure cube (effective area * livetime)
+                        m_cube(pixel, iebin) += rsp->aeff(theta, 0.0, 0.0, 0.0, logE) *
+                                                cta->livetime();
  
-                    // Add to exposure cube (effective area * livetime)
-                    m_cube(pixel, iebin) += rsp->aeff(theta, 0.0, 0.0, 0.0, logE) *
-                                            cta->livetime();
+                    } // endif: energy bin was inside observation energy boundaries
  
                 } // endfor: looped over energy bins
  
@@ -487,6 +487,9 @@ void GCTACubePsf::fill(const GObservations&amp; obs, GLog* log)
         // Extract region of interest from CTA observation
         GCTARoi roi = cta->roi();
  
+        // Extract energy boundaries from CTA observation
+        GEbounds obs_ebounds = cta->ebounds();
+
         // Check for RoI sanity
         if (!roi.is_valid()) {
             std::string msg = "No RoI information found in input observation "
@@ -537,6 +540,12 @@ void GCTACubePsf::fill(const GObservations&amp; obs, GLog* log)
                 // Loop over all energy bins
                 for (int iebin = 0; iebin < m_ebounds.size(); ++iebin) {
  
+                    // Skip if pixel is not within observation energy boundaries
+                    if (!obs_ebounds.contains(m_ebounds.emin(iebin),
+                                              m_ebounds.emax(iebin))) {
+                        continue;
+                    }
+
                     // Get logE/TeV
                     double logE = m_ebounds.elogmean(iebin).log10TeV();
  
@@ -84,6 +84,7 @@ public:
     GEnergy        elogmean(const int& index) const;
     GEnergy        ewidth(const int& index) const;
     bool           contains(const GEnergy& eng) const;
+    bool           contains(const GEnergy& emin, const GEnergy& emax) const;
 };
  
  
@@ -995,6 +995,34 @@ bool GEbounds::contains(const GEnergy&amp; eng) const
  
  
 /***********************************************************************//**
+ * @brief Checks whether energy boundaries contain and energy bin
+ *
+ * @param[in] emin Minimum energy of bin to be checked.
+ * @param[in] emax Maximum energy of bin to be checked.
+ * @return True if energy bin [emin, emax] is fully contained inside energy
+ * the boundaries, false otherwiese
+ *
+ * Checks if the energy interval [@p emin, @p emax ] falls is fully contained
+ * by energy boundaries
+ *
+ * @todo This method is so far only correct for contiguous energy boundaries.
+ ***************************************************************************/
+bool GEbounds::contains(const GEnergy& emin, const GEnergy& emax) const
+{
+    // Initialise result
+    bool contained = false;
+
+    // Check if emin, emax is fully contained within this ebounds
+    if (emin >= m_emin && emax <= m_emax) {
+        contained = true;
+    }
+
+    // Return result
+    return contained;
+}
+
+
+/***********************************************************************//**
  * @brief Print energy boundaries
  *
  * @param[in] chatter Chattiness (defaults to NORMAL).
@@ -242,6 +242,16 @@ void TestGObservation::test_ebounds(void)
         test_try_failure(e);
     }
  
+    // Check energy containment
+    GEbounds containment(1, GEnergy(100.0, "MeV"), GEnergy(1000.0, "MeV"));
+    test_assert(containment.contains(GEnergy(200.0, "MeV")), "Energy 200 MeV should be contained.");
+    test_assert(!containment.contains(GEnergy(10.0, "MeV")), "Energy 10 MeV should not be contained.");
+    test_assert(!containment.contains(GEnergy(2000.0, "MeV")), "Energy 2000 MeV should not be contained.");
+    test_assert(containment.contains(GEnergy(200.0, "MeV"), GEnergy(800.0, "MeV")), "Energy bin [200,800] MeV should be contained.");
+    test_assert(!containment.contains(GEnergy(80.0, "MeV"), GEnergy(800.0, "MeV")), "Energy bin [80,800] MeV should not be contained.");
+    test_assert(!containment.contains(GEnergy(200.0, "MeV"), GEnergy(2000.0, "MeV")), "Energy bin [200,2000] MeV should not be contained.");
+    test_assert(!containment.contains(GEnergy(80.0, "MeV"), GEnergy(2000.0, "MeV")), "Energy bin [80,2000] MeV should not be contained.");
+
     // Return
     return;
 }