Add sky position caching for each bin (#1970)

Removed dynamic scheduling as it introduces a small amount of overhead (it was only done so that the observations would be analyzed in numerical order) Added caching of sky directions for each pixel to prevent the need to re-convert the positions from index for each observation. Fill cached sky positions with RA,Dec values to prevent need to convert between galactic->celestial every time the distance is compute a distance from an observation ROI.

Add sky position caching for each bin (#1970)
Removed dynamic scheduling as it introduces a small amount of overhead (it was only done so that the observations would be analyzed in numerical order) Added caching of sky directions for each pixel to prevent the need to re-convert the positions from index for each observation. Fill cached sky positions with RA,Dec values to prevent need to convert between galactic->celestial every time the distance is compute a distance from an observation ROI.
Josh Cardenzana
1 parent 803c48d3
Showing 2 changed files with 22 additions and 15 deletions
src/ctbin/ctbin.cpp
src/ctbin/ctbin.hpp
@@ -248,7 +248,7 @@ void ctbin::run(void)
     log_header1(TERSE, gammalib::number("Bin observation", m_obs.size()));
     // Loop over all unbinned CTA observations in the container
-    #pragma omp parallel for schedule(dynamic)
+    #pragma omp parallel for
     for (int oindx=0; oindx<obs_list.size(); ++oindx) {
         GCTAObservation* obs = obs_list[oindx];
@@ -371,6 +371,7 @@ void ctbin::init_members(void)
     m_ebounds.clear();
     m_gti.clear();
     m_cube.clear();
+    m_dirs.clear();
     m_ontime   = 0.0;
     m_livetime = 0.0;
@@ -401,6 +402,7 @@ void ctbin::copy_members(const ctbin&amp; app)
     m_ebounds   = app.m_ebounds;
     m_gti       = app.m_gti;
     m_cube      = app.m_cube;
+    m_dirs      = app.m_dirs;
     m_ontime    = app.m_ontime;
     m_livetime  = app.m_livetime;
@@ -443,6 +445,13 @@ void ctbin::get_parameters(void)
     m_counts  = 0.0;
     m_weights = m_counts;
+    // Cache sky directions in m_counts
+    m_dirs.resize(m_counts.npix());
+    for (int pix=0; pix<m_counts.npix(); pix++) {
+        m_dirs[pix] = m_counts.inx2dir(pix);
+        m_dirs[pix].radec(m_dirs[pix].ra(), m_dirs[pix].dec());
+    }
+    
     // Get energy boundaries
     m_ebounds = cube.ebounds();
@@ -570,17 +579,17 @@ void ctbin::fill_cube(GCTAObservation* obs)
         // Update ontime and livetime
         m_ontime   += obs->ontime();
         m_livetime += obs->livetime();
-        
-        // Log filling results
-        log_header3(TERSE, get_obs_header(obs));
-        log_value(NORMAL, "Events in list", obs->events()->size());
-        log_value(NORMAL, "Events in cube", num_in_map);
-        log_value(NORMAL, "Events outside RoI", num_outside_roi);
-        log_value(NORMAL, "Events with invalid WCS", num_invalid_wcs);
-        log_value(NORMAL, "Events outside cube area", num_outside_map);
-        log_value(NORMAL, "Events outside energy bins", num_outside_ebds);
     }
+    // Log filling results
+    log_header3(TERSE, get_obs_header(obs));
+    log_value(NORMAL, "Events in list", obs->events()->size());
+    log_value(NORMAL, "Events in cube", num_in_map);
+    log_value(NORMAL, "Events outside RoI", num_outside_roi);
+    log_value(NORMAL, "Events with invalid WCS", num_invalid_wcs);
+    log_value(NORMAL, "Events outside cube area", num_outside_map);
+    log_value(NORMAL, "Events outside energy bins", num_outside_ebds);
+    
     // Return
     return;
 }
@@ -627,12 +636,9 @@ void ctbin::set_weights(GCTAObservation* obs)
     // Loop over all pixels in counts cube
     for (int pixel = 0; pixel < m_counts.npix(); ++pixel) {
-
-        // Get pixel sky direction
-        GSkyDir dir = m_counts.inx2dir(pixel);
-
+        
         // Skip pixel if it is outside the RoI
-        if (roi.centre().dir().dist_deg(dir) > roi.radius()) {
+        if (roi.centre().dir().dist_deg(m_dirs[pixel]) > roi.radius()) {
             continue;
         }
@@ -102,6 +102,7 @@ protected:
     GCTAEventCube m_cube;     //!< Events cube (for cube() method)
     double        m_ontime;   //!< Total ontime
     double        m_livetime; //!< Total livetime
+    std::vector<GSkyDir> m_dirs;  //!< Cached GSkyDir for each pixel in m_counts
 };