Merge branch '1439-change-usage-of-binned-observations-in-GCTACube-classes' into integration

Jürgen Knödlseder
2 parents ac71a9ff 0658d867
Showing 11 changed files with 278 additions and 196 deletions
inst/cta/include/GCTACubeBackground.hpp
inst/cta/include/GCTACubeExposure.hpp
inst/cta/include/GCTACubePsf.hpp
inst/cta/include/GCTAObservation.hpp
inst/cta/pyext/GCTACubeBackground.i
inst/cta/pyext/GCTACubeExposure.i
inst/cta/pyext/GCTACubePsf.i
inst/cta/pyext/GCTAObservation.i
inst/cta/src/GCTACubeBackground.cpp
inst/cta/src/GCTACubeExposure.cpp
inst/cta/src/GCTACubePsf.cpp
@@ -65,7 +65,7 @@ public:
     void                clear(void);
     GCTACubeBackground* clone(void) const;
     std::string         classname(void) const;
-    void                fill(const GObservations& obs);
+    void                fill(const GObservations& obs, GLog* log = NULL);
     double              integral(const double& logE) const;
     void                read(const GFits& fits);
     void                write(GFits& file) const;
@@ -74,7 +74,7 @@ public:
     GCTACubeExposure*  clone(void) const;
     std::string        classname(void) const;
     void               set(const GCTAObservation& obs);
-    void               fill(const GObservations& obs);
+    void               fill(const GObservations& obs, GLog* log = NULL);
     const GSkymap&     cube(void) const;
     const GEbounds&    ebounds(void) const;
     const GGti&        gti(void) const;
@@ -82,7 +82,7 @@ public:
     GCTACubePsf*       clone(void) const;
     std::string        classname(void) const;
     void               set(const GCTAObservation& obs);
-    void               fill(const GObservations& obs);
+    void               fill(const GObservations& obs, GLog* log = NULL);
     const GSkymap&     map(void) const;
     const GEbounds&    ebounds(void) const;
     const GNodeArray&  deltas(void) const;
@@ -124,6 +124,7 @@ public:
     void                deadc(const double& deadc);
     void                eventfile(const std::string& filename);
     const std::string&  eventfile(void) const;
+    const std::string&  eventtype(void) const;
     void                dispose_events(void);
     const double&       lo_user_thres(void) const;
     const double&       hi_user_thres(void) const;
@@ -413,6 +414,18 @@ const std::string&amp; GCTAObservation::eventfile(void) const
 /***********************************************************************//**
+ * @brief Return type of events
+ *
+ * @return Type of events.
+ ***************************************************************************/
+inline
+const std::string& GCTAObservation::eventtype(void) const
+{
+    return m_eventtype;
+}
+
+
+/***********************************************************************//**
  * @brief Return user low energy threshold
  *
  * @return User low energy threshold.
@@ -52,7 +52,7 @@ public:
     void                clear(void);
     GCTACubeBackground* clone(void) const;
     std::string         classname(void) const;
-    void                fill(const GObservations& obs);
+    void                fill(const GObservations& obs, GLog* log = NULL);
     double              integral(const double& logE) const;
     void                read(const GFits& fits);
     void                write(GFits& file) const;
@@ -65,7 +65,7 @@ public:
     GCTACubeExposure* clone(void) const;
     std::string       classname(void) const;
     void              set(const GCTAObservation& obs);
-    void              fill(const GObservations& obs);
+    void              fill(const GObservations& obs, GLog* log = NULL);
     const GSkymap&    cube(void) const;
     const GEbounds&   ebounds(void) const;
     const GGti&       gti(void) const;
@@ -74,7 +74,7 @@ public:
     GCTACubePsf*       clone(void) const;
     std::string        classname(void) const;
     void               set(const GCTAObservation& obs);
-    void               fill(const GObservations& obs);
+    void               fill(const GObservations& obs, GLog* log = NULL);
     const GSkymap&     map(void) const;
     const GEbounds&    ebounds(void) const;
     const GNodeArray&  deltas(void) const;
@@ -96,6 +96,7 @@ public:
     void                deadc(const double& deadc);
     void                eventfile(const std::string& filename);
     const std::string&  eventfile(void) const;
+    const std::string&  eventtype(void) const;
     void                dispose_events(void);
     const double&       lo_user_thres(void) const;
     const double&       hi_user_thres(void) const;
@@ -40,6 +40,7 @@
 #include "GCTARoi.hpp"
 #include "GCTAInstDir.hpp"
 #include "GCTACubeBackground.hpp"
+#include "GLog.hpp"
 /* __ Method name definitions ____________________________________________ */
 #define G_READ                             "GCTACubeBackground::read(GFits&)"
@@ -263,6 +264,7 @@ GCTACubeBackground* GCTACubeBackground::clone(void) const
  * @brief Fill background cube from observation container
  *
  * @param[in] obs Observation container.
+ * @param[in] log Pointer to logger (optional).
  *
  * @exception GException::invalid_value
  *            No event list found in CTA observations.
@@ -271,7 +273,7 @@ GCTACubeBackground* GCTACubeBackground::clone(void) const
  * for all CTA observations in an observation container. The cube pixel
  * values are computed as the sum over the background rates.
  ***************************************************************************/
-void GCTACubeBackground::fill(const GObservations& obs)
+void GCTACubeBackground::fill(const GObservations& obs, GLog* log)
 {
     // Clear background cube
     m_cube = 0.0;
@@ -288,49 +290,63 @@ void GCTACubeBackground::fill(const GObservations&amp; obs)
         // Get observation and continue only if it is a CTA observation
         const GCTAObservation* cta = dynamic_cast<const GCTAObservation*>(obs[i]);
-        if (cta != NULL) {
+        // Skip observation if it's not CTA
+        if (cta == NULL) {
+            if (log != NULL) {
+                *log << "Warning: Skipping "+obs[i]->instrument()+" observation ";
+                *log << "\"" << obs[i]->name() << "\"" <<std::endl;
+            }
+            continue;
+        }
-            // Extract region of interest from CTA observation
-            GCTARoi roi = cta->roi();
+        // Skip observation if we don't have an unbinned observation
+        if (cta->eventtype() != "EventList") {
+            if (log != NULL) {
+                *log << "Warning: Skipping binned observation ";
+                *log << "\"" << cta->name() << "\"" <<std::endl;
+            }
+            continue;
+        }
-            // Set GTI of actual observations as the GTI of the event cube
-            eventcube.gti(cta->gti());
+        // Extract region of interest from CTA observation
+        GCTARoi roi = cta->roi();
-            // Get observation livetime
-            double livetime = cta->livetime();
+        // Set GTI of actual observations as the GTI of the event cube
+        eventcube.gti(cta->gti());
-            // Loop over all bins in background cube
-            for (int i = 0; i < eventcube.size(); ++i) {
+        // Get observation livetime
+        double livetime = cta->livetime();
-                // Get event bin
-                GCTAEventBin* bin = eventcube[i];
+        // Loop over all bins in background cube
+        for (int i = 0; i < eventcube.size(); ++i) {
-                // Continue only if binned in contained in ROI
-                if (roi.contains(*bin)) {
+            // Get event bin
+            GCTAEventBin* bin = eventcube[i];
-                    // Compute model value for event bin. The model value is
-                    // given in counts/MeV/s/sr.
-                    double model = obs.models().eval(*bin, *cta);
+            // Continue only if binned in contained in ROI
+            if (roi.contains(*bin)) {
-                    // Multiply by livetime to get the correct weighting for
-                    // each observation. We divide by the total livetime later
-                    // to get the background model in units of counts/MeV/s/sr.
-                    model *= livetime;
+                // Compute model value for event bin. The model value is
+                // given in counts/MeV/s/sr.
+                double model = obs.models().eval(*bin, *cta);
-                    // Add existing number of counts
-                    model += bin->counts();
+                // Multiply by livetime to get the correct weighting for
+                // each observation. We divide by the total livetime later
+                // to get the background model in units of counts/MeV/s/sr.
+                model *= livetime;
-                    // Store cumulated value (units: counts/MeV/sr)
-                    bin->counts(model);
+                // Add existing number of counts
+                model += bin->counts();
-                } // endif: bin was contained in RoI
+                // Store cumulated value (units: counts/MeV/sr)
+                bin->counts(model);
-            } // endfor: looped over all bins
+            } // endif: bin was contained in RoI
-            // Accumulate livetime
-            total_livetime += livetime;
+        } // endfor: looped over all bins
-        } // endif: cta observation was vaild
+        // Accumulate livetime
+        total_livetime += livetime;
     } // endfor: looped over all observations
@@ -34,6 +34,7 @@
 #include "GCTAEventList.hpp"
 #include "GMath.hpp"
 #include "GTools.hpp"
+#include "GLog.hpp"
 /* __ Method name definitions ____________________________________________ */
@@ -298,55 +299,61 @@ GCTACubeExposure* GCTACubeExposure::clone(void) const
  ***************************************************************************/
 void GCTACubeExposure::set(const GCTAObservation& obs)
 {
-    // Clear GTIs, reset livetime and exposure cube pixels
-    m_gti.clear();
-    m_livetime = 0.0;
-    m_cube     = 0.0;
-    // Extract region of interest from CTA observation
-    GCTARoi roi = obs.roi();
+    // Only continue if we have an unbinned observation
+    if (obs.eventtype() == "EventList") {
-    // Get references on CTA response and pointing direction
-    const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(obs.response());
-    const GSkyDir&         pnt = obs.pointing().dir();
+        // Clear GTIs, reset livetime and exposure cube pixels
+        m_gti.clear();
+        m_livetime = 0.0;
+        m_cube     = 0.0;
-    // Continue only if response is valid
-    if (rsp != NULL) {
+        // Extract region of interest from CTA observation
+        GCTARoi roi = obs.roi();
-        // Loop over all pixels in sky map
-        for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
+        // Get references on CTA response and pointing direction
+        const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(obs.response());
+        const GSkyDir&         pnt = obs.pointing().dir();
-            // Get pixel sky direction
-            GSkyDir dir = m_cube.inx2dir(pixel);
-            
-            // Continue only if pixel is within RoI
-            if (roi.centre().dir().dist_deg(dir) <= roi.radius()) {
+        // Continue only if response is valid
+        if (rsp != NULL) {
-                // Compute theta angle with respect to pointing direction
-                // in radians
-                double theta = pnt.dist(dir);
-    
-                // Loop over all exposure cube energy bins
-                for (int iebin = 0; iebin < m_ebounds.size(); ++iebin){
+            // Loop over all pixels in sky map
+            for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
-                    // Get logE/TeV
-                    double logE = m_ebounds.elogmean(iebin).log10TeV();
+                // Get pixel sky direction
+                GSkyDir dir = m_cube.inx2dir(pixel);
-                    // Set exposure cube (effective area * lifetime)
-                    m_cube(pixel, iebin) = rsp->aeff(theta, 0.0, 0.0, 0.0, logE) *
-                                           obs.livetime();
+                // Continue only if pixel is within RoI
+                if (roi.centre().dir().dist_deg(dir) <= roi.radius()) {
-                } // endfor: looped over energy bins
+                    // Compute theta angle with respect to pointing direction
+                    // in radians
+                    double theta = pnt.dist(dir);
-            } // endif: pixel was within RoI
+                    // Loop over all exposure cube energy bins
+                    for (int iebin = 0; iebin < m_ebounds.size(); ++iebin){
-        } // endfor: looped over all pixels
+                        // Get logE/TeV
+                        double logE = m_ebounds.elogmean(iebin).log10TeV();
-        // Append GTIs and increment livetime
-        m_gti.extend(obs.gti());
-        m_livetime += obs.livetime();
+                        // Set exposure cube (effective area * lifetime)
+                        m_cube(pixel, iebin) = rsp->aeff(theta, 0.0, 0.0, 0.0, logE) *
+                                               obs.livetime();
+
+                    } // endfor: looped over energy bins
+
+                } // endif: pixel was within RoI
+    
+            } // endfor: looped over all pixels
-    } // endif: response was valid
+            // Append GTIs and increment livetime
+            m_gti.extend(obs.gti());
+            m_livetime += obs.livetime();
+
+        } // endif: response was valid
+
+    } // endif: observation was unbinned
     // Return
     return;
@@ -357,6 +364,7 @@ void GCTACubeExposure::set(const GCTAObservation&amp; obs)
  * @brief Fill exposure cube from observation container
  *
  * @param[in] obs Observation container.
+ * @param[in] log Pointer to logger (optional).
  *
  * @exception GException::invalid_value
  *            No event list found in CTA observations.
@@ -365,7 +373,7 @@ void GCTACubeExposure::set(const GCTAObservation&amp; obs)
  * an observation container. The cube pixel values are computed as the sum
  * over the products of the effective area and the livetime.
  ***************************************************************************/
-void GCTACubeExposure::fill(const GObservations& obs)
+void GCTACubeExposure::fill(const GObservations& obs, GLog* log)
 {
     // Clear GTIs, reset livetime and exposure cube pixels
     m_gti.clear();
@@ -377,55 +385,74 @@ void GCTACubeExposure::fill(const GObservations&amp; obs)
         // Get observation and continue only if it is a CTA observation
         const GCTAObservation* cta = dynamic_cast<const GCTAObservation*>(obs[i]);
-        if (cta != NULL) {
-            // Extract region of interest from CTA observation
-            GCTARoi roi = cta->roi();
+        // Skip observation if it's not CTA
+        if (cta == NULL) {
+            if (log != NULL) {
+                *log << "Warning: Skipping "+obs[i]->instrument()+" observation ";
+                *log << "\"" << obs[i]->name() << "\"" <<std::endl;
+            }
+            continue;
+        }
+
+        // Skip observation if we don't have an unbinned observation
+        if (cta->eventtype() != "EventList") {
+            if (log != NULL) {
+                *log << "Warning: Skipping binned observation ";
+                *log << "\"" << cta->name() << "\"" <<std::endl;
+            }
+            continue;
+        }
+
+        // Extract region of interest from CTA observation
+        GCTARoi roi = cta->roi();
-            // Get references on CTA response and pointing direction
-            const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(cta->response());
-            const GSkyDir&         pnt = cta->pointing().dir();
+        // Get references on CTA response and pointing direction
+        const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(cta->response());
+        const GSkyDir&         pnt = cta->pointing().dir();
-            // Continue only if response is valid
-            if (rsp != NULL) {
+        // Skip observation if we don't have an unbinned observation
+        if (rsp == NULL) {
+            if (log != NULL) {
+                *log << "Warning: Observation \"" << cta->name();
+                *log << "\" contains no IRF response." <<std::endl;
+            }
+            continue;
+        }
-                // Loop over all pixels in sky map
-                for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
+        // Loop over all pixels in sky map
+        for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
-                    // Get pixel sky direction
-                    GSkyDir dir = m_cube.inx2dir(pixel);
+            // Get pixel sky direction
+            GSkyDir dir = m_cube.inx2dir(pixel);
-                    // Continue only if pixel is within RoI
-                    if (roi.centre().dir().dist_deg(dir) <= roi.radius()) {
+            // Continue only if pixel is within RoI
+            if (roi.centre().dir().dist_deg(dir) <= roi.radius()) {
-                        // Compute theta angle with respect to pointing
-                        // direction in radians
-                        double theta = pnt.dist(dir);
-    
-                        // Loop over all exposure cube energy bins
-                        for (int iebin = 0; iebin < m_ebounds.size(); ++iebin){
+                // Compute theta angle with respect to pointing
+                // direction in radians
+                double theta = pnt.dist(dir);
+
+                // 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();
+                    // 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();
-                        } // endfor: looped over energy bins
+                } // endfor: looped over energy bins
-                    } // endif: pixel within RoI
+            } // endif: pixel within RoI
-                } // endfor: looped over all pixels
+        } // endfor: looped over all pixels
-                // Append GTIs and increment livetime
-                m_gti.extend(cta->gti());
-                m_livetime += cta->livetime();
+        // Append GTIs and increment livetime
+        m_gti.extend(cta->gti());
+        m_livetime += cta->livetime();
-            } // endif: response was valid
-    
-        } // endif: observation was a CTA observation
-
     } // endfor: looped over observations
     // Return
@@ -34,6 +34,7 @@
 #include "GCTAEventList.hpp"
 #include "GMath.hpp"
 #include "GTools.hpp"
+#include "GLog.hpp"
 /* __ Method name definitions ____________________________________________ */
 #define G_SET                            "GCTACubePsf::set(GCTAObservation&)"
@@ -357,64 +358,69 @@ void GCTACubePsf::set(const GCTAObservation&amp; obs)
     // Clear PSF cube
     clear_cube();
-    // Extract region of interest from CTA observation
-    const GCTAEventList* list = dynamic_cast<const GCTAEventList*>(obs.events());
-    if (list == NULL) {
-        std::string msg = "CTA Observation does not contain an event "
-                          "list. Event list information is needed to "
-                          "retrieve the Region of Interest for each "
-                          "CTA observation.";
-        throw GException::invalid_value(G_SET, msg);
-    }
-    const GCTARoi& roi = list->roi();
+    // Only continue if we have an unbinned observation
+    if (obs.eventtype() == "EventList") {
+
+        // Extract region of interest from CTA observation
+        const GCTAEventList* list = dynamic_cast<const GCTAEventList*>(obs.events());
+        if (list == NULL) {
+            std::string msg = "CTA Observation does not contain an event "
+                              "list. Event list information is needed to "
+                              "retrieve the Region of Interest for each "
+                              "CTA observation.";
+            throw GException::invalid_value(G_SET, msg);
+        }
+        const GCTARoi& roi = list->roi();
-    // Get references on CTA response and pointing direction
-    const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(obs.response());
-    const GSkyDir&         pnt = obs.pointing().dir();
+        // Get references on CTA response and pointing direction
+        const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(obs.response());
+        const GSkyDir&         pnt = obs.pointing().dir();
-    // Continue only if response is valid
-    if (rsp != NULL) {
+        // Continue only if response is valid
+        if (rsp != NULL) {
-        // Loop over all pixels in sky map
-        for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
+            // Loop over all pixels in sky map
+            for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
+    
+                // Get pixel sky direction
+                GSkyDir dir = m_cube.inx2dir(pixel);
+                
+                // Continue only if pixel is within RoI
+                if (roi.centre().dir().dist_deg(dir) <= roi.radius()) {
-            // Get pixel sky direction
-            GSkyDir dir = m_cube.inx2dir(pixel);
-            
-            // Continue only if pixel is within RoI
-            if (roi.centre().dir().dist_deg(dir) <= roi.radius()) {
+                    // Compute theta angle with respect to pointing direction
+                    // in radians
+                    double  theta = pnt.dist(dir);
-                // Compute theta angle with respect to pointing direction
-                // in radians
-                double  theta = pnt.dist(dir);
-    
-                // Loop over all exposure cube energy bins
-                for (int iebin = 0; iebin < m_ebounds.size(); ++iebin){
+                    // 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();
+                        // Get logE/TeV
+                        double logE = m_ebounds.elogmean(iebin).log10TeV();
-                    // Loop over delta values
-                    for (int idelta = 0; idelta < m_deltas.size(); ++idelta) {
+                        // Loop over delta values
+                        for (int idelta = 0; idelta < m_deltas.size(); ++idelta) {
-                        // Compute delta in radians
-                        double delta = m_deltas[idelta] * gammalib::deg2rad;
+                            // Compute delta in radians
+                            double delta = m_deltas[idelta] * gammalib::deg2rad;
-                        // Set map index
-                        int imap = offset(idelta, iebin);
-                
-                        // Set PSF cube
-                        m_cube(pixel, imap) = rsp->psf(delta, theta, 0.0, 0.0, 0.0, logE);
+                            // Set map index
+                            int imap = offset(idelta, iebin);
-                    } // endfor: looped over delta bins
+                            // Set PSF cube
+                            m_cube(pixel, imap) = rsp->psf(delta, theta, 0.0, 0.0, 0.0, logE);
-                } // endfor: looped over energy bins
+                        } // endfor: looped over delta bins
-            } // endif: pixel was within RoI
+                    } // endfor: looped over energy bins
-        } // endfor: looped over all pixels
+                } // endif: pixel was within RoI
-    } // endif: response was valid
+            } // endfor: looped over all pixels
+
+        } // endif: response was valid
+
+    } // endif: observation was unbinned
     // Compile option: guarantee smooth Psf
     #if defined(G_SMOOTH_PSF)
@@ -431,7 +437,7 @@ void GCTACubePsf::set(const GCTAObservation&amp; obs)
  *
  * @param[in] obs Observation container.
  ***************************************************************************/
-void GCTACubePsf::fill(const GObservations& obs)
+void GCTACubePsf::fill(const GObservations& obs, GLog* log)
 {
     // Clear PSF cube
     clear_cube();
@@ -444,68 +450,87 @@ void GCTACubePsf::fill(const GObservations&amp; obs)
         // Get observation and continue only if it is a CTA observation
         const GCTAObservation* cta = dynamic_cast<const GCTAObservation*>(obs[i]);
-        if (cta != NULL) {
-            // Extract region of interest from CTA observation
-            GCTARoi roi = cta->roi();
+        // Skip observation if it's not CTA
+        if (cta == NULL) {
+            if (log != NULL) {
+                *log << "Warning: Skipping "+obs[i]->instrument()+" observation ";
+                *log << "\"" << obs[i]->name() << "\"" <<std::endl;
+            }
+            continue;
+        }
-            // Get references on CTA response and pointing direction
-            const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(cta->response());
-            const GSkyDir&         pnt = cta->pointing().dir();
+        // Skip observation if we don't have an unbinned observation
+        if (cta->eventtype() != "EventList") {
+            if (log != NULL) {
+                *log << "Warning: Skipping binned observation ";
+                *log << "\"" << cta->name() << "\"" <<std::endl;
+            }
+            continue;
+        }
-            // Continue only if response is valid
-            if (rsp != NULL) {
+        // Extract region of interest from CTA observation
+        GCTARoi roi = cta->roi();
-                // Loop over all pixels in sky map
-                for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
+        // Get references on CTA response and pointing direction
+        const GCTAResponseIrf* rsp = dynamic_cast<const GCTAResponseIrf*>(cta->response());
+        const GSkyDir&         pnt = cta->pointing().dir();
-                    // Get pixel sky direction
-                    GSkyDir dir = m_cube.inx2dir(pixel);
-                    
-                    // Continue only if pixel is within RoI
-                    if (roi.centre().dir().dist_deg(dir) <= roi.radius()) {
+        // Skip observation if we don't have an unbinned observation
+        if (rsp == NULL) {
+            if (log != NULL) {
+                *log << "Warning: Observation \"" << cta->name();
+                *log << "\" contains no IRF response." <<std::endl;
+            }
+            continue;
+        }
-                        // Compute theta angle with respect to pointing
-                        // direction in radians
-                        double theta = pnt.dist(dir);
-    
-                        // Loop over all energy bins
-                        for (int iebin = 0; iebin < m_ebounds.size(); ++iebin) {
+        // Loop over all pixels in sky map
+        for (int pixel = 0; pixel < m_cube.npix(); ++pixel) {
-                            // Get logE/TeV
-                            double logE = m_ebounds.elogmean(iebin).log10TeV();
+            // Get pixel sky direction
+            GSkyDir dir = m_cube.inx2dir(pixel);
+                    
+            // Continue only if pixel is within RoI
+            if (roi.centre().dir().dist_deg(dir) <= roi.radius()) {
-                            // Compute exposure weight
-                            double weight = rsp->aeff(theta, 0.0, 0.0, 0.0, logE) *
-                                            cta->livetime();
+                // Compute theta angle with respect to pointing
+                // direction in radians
+                double theta = pnt.dist(dir);
-                            // Accumulate weights
-                            exposure(pixel, iebin) += weight;
+                // Loop over all energy bins
+                for (int iebin = 0; iebin < m_ebounds.size(); ++iebin) {
-                            // Loop over delta values
-                            for (int idelta = 0; idelta < m_deltas.size(); ++idelta) {
+                    // Get logE/TeV
+                    double logE = m_ebounds.elogmean(iebin).log10TeV();
-                                // Compute delta in radians
-                                double delta = m_deltas[idelta] * gammalib::deg2rad;
+                    // Compute exposure weight
+                    double weight = rsp->aeff(theta, 0.0, 0.0, 0.0, logE) *
+                                    cta->livetime();
-                                // Set map index
-                                int imap = offset(idelta, iebin);
+                    // Accumulate weights
+                    exposure(pixel, iebin) += weight;
-                                // Add on PSF cube
-                                m_cube(pixel, imap) +=
-                                    rsp->psf(delta, theta, 0.0, 0.0, 0.0, logE) * weight;
+                    // Loop over delta values
+                    for (int idelta = 0; idelta < m_deltas.size(); ++idelta) {
-                            } // endfor: looped over delta bins
+                        // Compute delta in radians
+                        double delta = m_deltas[idelta] * gammalib::deg2rad;
-                        } // endfor: looped over energy bins
+                        // Set map index
+                        int imap = offset(idelta, iebin);
+
+                        // Add on PSF cube
+                        m_cube(pixel, imap) +=
+                            rsp->psf(delta, theta, 0.0, 0.0, 0.0, logE) * weight;
-                    } // endif: pixel was within RoI
+                    } // endfor: looped over delta bins
-                } // endfor: looped over all pixels
+                } // endfor: looped over energy bins
-            } // endif: response was valid
+            } // endif: pixel was within RoI
-        } // endif: observation was a CTA observation
+        } // endfor: looped over all pixels
     } // endfor: looped over observations