Enable ctbin to produce counts cubes for a joint binned analysis (#2714)

Jürgen Knödlseder
1 parent cf3f1ec9
Showing 8 changed files with 448 additions and 313 deletions
ChangeLog
NEWS
doc/source/admin/release_history/1.6.rst
pyext/ctbin.i
src/ctbin/ctbin.cpp
src/ctbin/ctbin.hpp
src/ctbin/ctbin.par
test/test_ctbin.py
-2018-11-03
+2018-11-07
  
         * Version 1.6.0 released
         ========================
  
+        Enable ctbin to produce counts cubes for a joint binned analysis (#2714)
         Add code generator (#2709)
         Add csbkgmodel script (#2694)
         Compute DETX and DETY in ctmodel if they are not available (#2694)
@@ -55,6 +55,14 @@ of the filename is now done by the GApplicationPar class.
  
 ctbin - CTA event binning
 -------------------------
+Added the option to produce an individual counts cube for each event list in
+the input observation definition XML file. The option is enabled by setting
+the new hidden parameter "stack" to "no". The cube() method now has an optional
+index argument, providing access to the i'th counts cube that was produced.
+The cubes() method specifies the number of counts cubes that were produced.
+If no counts cube was produced, the cube() method now throws an exception
+(#2714).
+
 ctbin now accouts for varying energy thresholds when computing the weight
 array (#2685).
  
@@ -41,6 +41,8 @@ Bug fixes
 Improvements
 ------------
  
+* [`2714 <https://cta-redmine.irap.omp.eu/issues/2714>`_] -
+  Enable :ref:`ctbin` to produce counts cubes for a joint binned analysis
 * [`2685 <https://cta-redmine.irap.omp.eu/issues/2685>`_] -
   Add ``ctbkgcube::cntcube()`` method
   Account for varying energy thresholds in :ref:`ctbin` weight computation
 /***************************************************************************
  *                        ctbin - Event binning tool                       *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010-2016 by Juergen Knoedlseder                         *
+ *  copyright (C) 2010-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -46,10 +46,10 @@ public:
     // Methods
     void                 clear(void);
     void                 run(void);
+    int                  cubes(void) const;
+    const GCTAEventCube& cube(const int& index = 0) const;
     void                 save(void);
     void                 publish(const std::string& name = "");
-    const GCTAEventCube& cube(void) const;
-    const GCTAEventCube& cube(const size_t& index) const;
 };
  
  
@@ -38,6 +38,7 @@
 #endif
  
 /* __ Method name definitions ____________________________________________ */
+#define G_CUBE                                            "ctbin::cube(int&)"
 #define G_GET_PARAMETERS                            "ctbin::get_parameters()"
 #define G_FILL_CUBE                      "ctbin::fill_cube(GCTAObservation*)"
 #define G_SET_WEIGHTS                  "ctbin::set_weights(GCTAObservation*)"
@@ -234,6 +235,9 @@ void ctbin::run(void)
     // Get task parameters
     get_parameters();
  
+    // Initialise sky direction cache for stacking
+    init_sky_dir_cache();
+
     // Write input observation container into logger
     log_observations(NORMAL, m_obs, "Input observation");
  
@@ -255,33 +259,33 @@ void ctbin::run(void)
         log_string(NORMAL, msg);
     }
  
-    // Initialise cubes and ids if unstacked
-    if (!m_stack) {
-	m_cubes = std::vector<GCTAEventCube>(obs_list.size());
-	m_ids = std::vector<std::string>(obs_list.size());
+    // Set number of relevant observations
+    int nobs = obs_list.size();
+
+    // Initialise arrays
+    if (nobs > 0) {
+        m_counts  = std::vector<GSkyMap>(nobs);
+        m_weights = std::vector<GSkyMap>(nobs);
+        m_cubes   = std::vector<GCTAEventCube>(nobs);
+        m_ids     = std::vector<std::string>(nobs);
+    }
+    else {
+        m_counts.clear();
+        m_weights.clear();
+        m_cubes.clear();
+        m_ids.clear();
     }
  
     // Write header into logger
-    log_header1(TERSE, gammalib::number("Bin observation", m_obs.size()));
+    log_header1(TERSE, gammalib::number("Bin observation", nobs));
  
     // Loop over all unbinned CTA observations in the container
     #pragma omp parallel for
-    for (int i = 0; i < obs_list.size(); ++i) {
+    for (int i = 0; i < nobs; ++i) {
  
         // Get pointer to observation
         GCTAObservation* obs = obs_list[i];
  
-	// Make sure that the observation holds a CTA event list. If this
-	// is not the case then throw an exception.
-	const GCTAEventList* events = dynamic_cast<const GCTAEventList*>
-	    (obs->events());
-	if (events == NULL) {
-	    std::string msg = "CTA Observation does not contain an event "
-		"list. An event list is needed to fill the "
-		"counts cube.";
-	    throw GException::invalid_value(G_FILL_CUBE, msg);
-	}
-
         // Fill the cube
         GSkyMap counts = fill_cube(obs);
  
@@ -291,88 +295,118 @@ void ctbin::run(void)
         // Dispose events to free memory
         obs->dispose_events();
  
+        // Post process cubes
         #pragma omp critical(ctbin_fill_cube)
-	if (m_stack) {
-	    
-	    
-	    m_counts += counts;
-
-	    for (int iebin = 0; iebin < m_ebounds.size(); ++iebin) {
-		for (int pixel = 0; pixel < m_counts.npix(); ++pixel) {
-		    if (weights(pixel, iebin) == 1.0) {
-			m_weights(pixel, iebin) = 1.0;
-		    }
-		}
-	    }
-	
-	    
-	    // If GTI is empty then set its time reference from the observation.
-	    // From then on we keep that time reference
+        {
+            // If stacking of cubes is requested then add up counts cube
+            // and set weights in the first slot of the m_counts and
+            // m_weights arrays
+            if (m_stack) {
+                if (i == 0) {
+                    m_counts[0]  = counts;
+                    m_weights[0] = weights;
+                }
+                else {
+                    for (int iebin = 0; iebin < m_ebounds.size(); ++iebin) {
+                        for (int pixel = 0; pixel < counts.npix(); ++pixel) {
+                            m_counts[0](pixel, iebin) += counts(pixel, iebin);
+                            if (weights(pixel, iebin) == 1.0) {
+                                m_weights[0](pixel, iebin) = 1.0;
+                            }
+                        }
+                    }
+                }
+            }
+
+            // ... otherwise store counts, weights and observation ID
+        	else {
+                m_counts[i]  = counts;
+                m_weights[i] = weights;
+                m_ids[i]     = obs->id();
+            }
+
+        } // end: omp critical section
+
+    } // endfor: looped over observations
+
+    // Combine GTIs of all observations and compute total ontime and livetime
+    for (int i = 0; i < nobs; ++i) {
+
+        // Get pointer to observation
+        GCTAObservation* obs = obs_list[i];
+
+        // Get Good Time Intervals for observation
+        GGti gti = obs->gti();
+
+        // If resulting GTI is not filled then copy over the reference from
+        // the observation and use it for the resulting GTI
 	    if (m_gti.is_empty()) {
-		m_gti.reference(events->gti().reference());
+            m_gti.reference(gti.reference());
 	    }
-	    
+
 	    // Append GTIs
-	    m_gti.extend(events->gti());
-	    
+	    m_gti.extend(gti);
+
 	    // Update ontime and livetime
 	    m_ontime   += obs->ontime();
 	    m_livetime += obs->livetime();
-	    
-	} else { // don't stack cubes
  
-	    // Construct a cube from this observation
-	    GCTAEventCube cube(counts, weights, m_ebounds, events->gti());
-	    m_cubes[i] = cube;
-	    obs->events(cube);
+    } // endfor: looped over all observations
+
+    // Post process data for a stacked cube
+    if (m_stack && nobs > 0) {
+
+        // Compute livetime fraction per energy bin so that varying energy
+        // thresholds are correctly taken into account in the weighting
+        // array. Computational speed is not optimum, we have to see later
+        // whether this is an issue.
+        for (int iebin = 0; iebin < m_ebounds.size(); ++iebin) {
+            double livetime_ebin = 0.0;
+            for (int i = 0; i < obs_list.size(); ++i) {
+                std::vector<bool> usage =
+                     cube_layer_usage(m_ebounds, obs_list[i]->ebounds());
+                if (usage[iebin]) {
+                    livetime_ebin += obs_list[i]->livetime();
+                }
+            }
+            if (m_livetime > 0.0) {
+                double fraction = livetime_ebin / m_livetime;
+                for (int pixel = 0; pixel < m_counts[0].npix(); ++pixel) {
+                    if (m_weights[0](pixel, iebin) > 0.0) {
+                        m_weights[0](pixel, iebin) *= fraction;
+                    }
+                }
+            }
+        }
  
-	    // Store id of the observation
-	    m_ids[i] = obs->id();
-	    
-	}
+        // Set a single cube in the observation container
+        obs_cube_stacked();
  
+    } // endif: stacked cube post processing
  
-    } // endfor: looped over observations
+    // Post process data for non-stacked cubes
+    else {
  
-    if (m_stack) {
-	// Compute total livetime in all observations
-	double livetime = 0.0;
-	for (int i = 0; i < obs_list.size(); ++i) {
-	    livetime += obs_list[i]->livetime();
-	}
-	
-	// Compute livetime fraction per energy bin so that varying energy
-	// thresholds are correctly taken into account in the weighting
-	// array. Computational speed is not optimum, we have to see later
-	// whether this is an issue.
-	for (int iebin = 0; iebin < m_ebounds.size(); ++iebin) {
-	    double livetime_ebin = 0.0;
-	    for (int i = 0; i < obs_list.size(); ++i) {
-		const GCTAEventList* events = dynamic_cast<const GCTAEventList*>
-		    (obs_list[i]->events());
-		std::vector<bool> usage  = cube_layer_usage(m_ebounds,
-							    events->ebounds());
-		if (usage[iebin]) {
-		    livetime_ebin += obs_list[i]->livetime();
-		}
-	    }
-	    if (livetime > 0.0) {
-		double fraction = livetime_ebin / livetime;
-		for (int pixel = 0; pixel < m_counts.npix(); ++pixel) {
-		    if (m_weights(pixel, iebin) > 0.0) {
-			m_weights(pixel, iebin) *= fraction;
-		    }
-		}
-	    }
-	}
+        // Loop over all relevant observations
+        for (int i = 0; i < nobs; ++i) {
  
-	// Build event cube (needs to come before obs_cube() since this method
-	// relies on correct setting of m_cube)
-	m_cube = GCTAEventCube(m_counts, m_weights, m_ebounds, m_gti);
+            // Get pointer to observation
+            GCTAObservation* obs = obs_list[i];
  
-	// Set a single cube in the observation container
-	obs_cube();
-    }
+            // Build counts cube
+            m_cubes[i] = GCTAEventCube(m_counts[i], m_weights[i], m_ebounds,
+                                       obs->gti());
+
+            // Attach event cube in the correct slot of the observation
+            // container
+            obs->events(m_cubes[i]);
+
+            // Store ID of the observation for saving of XML file
+            m_ids[i] = obs->id();
+
+        } // endfor: looped over relevant observations
+
+    } // endelse: non-stacked cube post processing
  
     // Write resulting observation container into logger
     log_observations(NORMAL, m_obs, "Binned observation");
@@ -388,72 +422,108 @@ void ctbin::run(void)
  
  
 /***********************************************************************//**
+ * @brief Return event cube at index
+ *
+ * @param[in] index Cube index (0,...,cubes()-1).
+ * @return Reference to event cube.
+ *
+ * @exception GException::out_of_range
+ *            Event cube index out of range.
+ *
+ * Returns a reference to the event cube at the given index.
+ ***************************************************************************/
+const GCTAEventCube& ctbin::cube(const int& index) const
+{
+    // Compile option: raise an exception if index is out of range
+    if (index < 0 || index >= cubes()) {
+        throw GException::out_of_range(G_CUBE, "Cube index", index, cubes());
+    }
+
+    // Return cube
+    return m_cubes[index];
+}
+
+
+/***********************************************************************//**
  * @brief Save counts cube
  *
- * Saves the counts cube into a FITS file.
+ * Saves the counts cube. If the counts cube was stacked single FITS file
+ * will be written. Otherwise each generated counts cube will be written
+ * into a FITS file, and an observation definition XML file will be written
+ * out that contains an observation for each counts cube.
  ***************************************************************************/
 void ctbin::save(void)
 {
     // Write header
     log_header1(TERSE, "Save counts cube");
  
-    if (!m_stack) {
-
-	// Loop over all observations
-	for (size_t i = 0; i < m_obs.size(); ++i) {
-	    
-	    // Loop over all ids of observations that were binned
-	    for (size_t j = 0; j < m_ids.size(); ++j) {
-	    
-		// Check if ids are equal
-		if (m_ids[j] == m_obs[i]->id()) {
-		    
-		    // Construct filename
-		    std::string filename(m_ids[j]);
-		    filename.insert(0, m_prefix);
-		    filename.append(".fits");
-
-		    // Set filename for counts cube 
-		    GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
-		    obs->eventfile(filename);
-
-		    // Save countscube
-		    obs->save(filename, clobber());
-		}
-	    }
-	}
-	if ((*this)["outcube"].is_valid()) {
-	    // Get counts cube filename
-	    m_outcube = (*this)["outcube"].filename();
-	    m_obs.save(m_outcube);
-	}
-	
-    } else { // else: stack
-
-	// Save only if filename is valid and if there is at least one
-	// observation
-	if ((*this)["outcube"].is_valid() && m_obs.size() > 0) {
-	    
-	    // Get counts cube filename
-	    m_outcube = (*this)["outcube"].filename();
-	    
-	    // Get CTA observation from observation container
-	    GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]);
-	    
-	    // Save only if observation is valid
-	    if (obs != NULL) {
-		
-		// Log counts cube file name
-		log_value(NORMAL, "Counts cube file", m_outcube.url());
-		
-		// Save cube
-		obs->save(m_outcube, clobber());
-		
-	    } // endif: observation was valid
-	    
-	} // endif: outcube file was valid
-	
-    } // endelse: stack 
+    // If counts cube was stacked then save first observation
+    if (m_stack) {
+
+        // Save only if filename is valid and if there is at least one
+        // observation
+        if ((*this)["outcube"].is_valid() && m_obs.size() > 0) {
+
+            // Get counts cube filename
+            m_outcube = (*this)["outcube"].filename();
+
+            // Get CTA observation from observation container
+            GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[0]);
+
+            // Save only if observation is valid
+            if (obs != NULL) {
+
+                // Log counts cube file name
+                log_value(NORMAL, "Counts cube file", m_outcube.url());
+
+                // Save cube
+                obs->save(m_outcube, clobber());
+
+            } // endif: observation was valid
+
+        } // endif: outcube file was valid
+
+    } // endif: counts cube was stacked
+
+    // ... otherwise write out all relevant counts cubes and the observation
+    // definition XML file
+    else {
+
+        // Loop over all observations
+        for (int i = 0; i < m_obs.size(); ++i) {
+
+            // Loop over all ids of observations that were binned
+            for (int j = 0; j < m_ids.size(); ++j) {
+
+                // Skip if ID differes
+                if (m_ids[j] != m_obs[i]->id()) {
+                    continue;
+                }
+
+                // Construct filename
+                std::string filename(m_ids[j]);
+                filename.insert(0, m_prefix);
+                filename.append(".fits");
+
+                // Set filename for counts cube
+                GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
+                obs->eventfile(filename);
+
+                // Save countscube
+                obs->save(filename, clobber());
+
+            } // endfor: looped over all IDs
+
+        } // endfor: looped over all observations
+
+        // Write XML file
+        if ((*this)["outcube"].is_valid()) {
+            m_outcube = (*this)["outcube"].filename();
+            m_obs.save(m_outcube);
+        }
+
+    } // endelse: counts cube were not stacked
+
     // Return
     return;
 }
@@ -478,8 +548,10 @@ void ctbin::publish(const std::string&amp; name)
     // Write counts cube name into logger
     log_value(NORMAL, "Counts cube name", user_name);
  
-    // Publish counts cube
-    m_counts.publish(user_name);
+    // Publish first counts cube if it exists
+    if (cubes() > 0) {
+        m_counts[0].publish(user_name);
+    }
  
     // Return
     return;
@@ -505,15 +577,14 @@ void ctbin::init_members(void)
     m_chatter = static_cast<GChatter>(2);
  
     // Initialise protected members
+    m_cubes.clear();
     m_counts.clear();
     m_weights.clear();
+    m_ids.clear();
     m_ebounds.clear();
     m_gti.clear();
-    m_cube.clear();
     m_ontime   = 0.0;
     m_livetime = 0.0;
-    m_cubes.clear();
-    m_ids.clear();
  
     // Initialise cache members
     m_dirs.clear();
@@ -538,15 +609,14 @@ void ctbin::copy_members(const ctbin&amp; app)
     m_chatter = app.m_chatter;
  
     // Copy protected members
+    m_cubes    = app.m_cubes;
     m_counts   = app.m_counts;
     m_weights  = app.m_weights;
+    m_ids      = app.m_ids;
     m_ebounds  = app.m_ebounds;
     m_gti      = app.m_gti;
-    m_cube     = app.m_cube;
     m_ontime   = app.m_ontime;
     m_livetime = app.m_livetime;
-    m_cubes = app.m_cubes;
-    m_ids = app.m_ids;
  
     // Copy cache members
     m_dirs = app.m_dirs;
@@ -581,14 +651,8 @@ void ctbin::get_parameters(void)
     // information and do not accept counts cubes.
     setup_observations(m_obs, false, true, false);
  
-    // Create an event cube based on task parameters
-    GCTAEventCube cube = create_cube(m_obs);
-
-    // Get the skymap from the cube and initialise all counts cube bins and
-    // weights to zero
-    m_counts  = cube.counts();
-    m_counts  = 0.0;
-    m_weights = m_counts;
+    // Create an event cube to query task parameters
+    GCTAEventCube cube = this->ctool::create_cube(m_obs);
  
     // Get energy boundaries
     m_ebounds = cube.ebounds();
@@ -608,13 +672,6 @@ void ctbin::get_parameters(void)
         (*this)["outcube"].query();
     }
  
-    // 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());
-    }
-
     // Set number of OpenMP threads
     #ifdef _OPENMP
     int nthreads = (*this)["nthreads"].integer();
@@ -632,6 +689,95 @@ void ctbin::get_parameters(void)
  
  
 /***********************************************************************//**
+ * @brief Initialise sky direction cache for cube stack
+ *
+ * Initialises a cache that holds the sky directions of each counts cube
+ * pixel for a stacked analysis. This will speed up the computations later.
+ *
+ * If no stacking is requested, the method does nothing.
+ ***************************************************************************/
+void ctbin::init_sky_dir_cache(void)
+{
+    // Continue only is stacking is requested
+    if (m_stack) {
+
+        // Create a sky map from task parameters
+        GSkyMap map = this->ctool::create_map(m_obs);
+
+        // Initialise sky direction
+        m_dirs.resize(map.npix());
+
+        // Compute sky direction for each pixel
+        for (int pixel = 0; pixel < map.npix(); ++pixel) {
+            m_dirs[pixel] = map.inx2dir(pixel);
+        }
+
+    } // endif: stacking was requested
+
+    // Return
+    return;
+}
+
+
+/***********************************************************************//**
+ * @brief Create counts cube sky map for current observation
+ *
+ * @param[in] obs CTA observation.
+ *
+ * Creates a counts cube sky map for current observations.
+ ***************************************************************************/
+GSkyMap ctbin::create_cube(const GCTAObservation* obs)
+{
+    // Read coordinate system and projection
+    std::string coordsys = (*this)["coordsys"].string();
+    std::string proj     = (*this)["proj"].string();
+
+    // Read sky map reference
+    double xref   = 0.0;
+    double yref   = 0.0;
+    bool   usepnt = (*this)["usepnt"].boolean();
+    if (!usepnt) {
+        xref = (*this)["xref"].real();
+        yref = (*this)["yref"].real();
+    }
+
+    // Read sky map pixel size and number of bins
+    double binsz = (*this)["binsz"].real();
+    int    nxpix = (*this)["nxpix"].integer();
+    int    nypix = (*this)["nypix"].integer();
+
+    // Read energy boundaries
+    GEbounds ebounds = create_ebounds();
+
+    // If requested, get pointing from observations
+    if (usepnt) {
+
+        // Get pointing direction for observation
+        GSkyDir pnt = obs->pointing().dir();
+
+        // Set xref/yref based on the coordinate system
+        if (gammalib::toupper(coordsys) == "GAL") {
+            xref = pnt.l_deg();
+            yref = pnt.b_deg();
+        }
+        else {
+            xref = pnt.ra_deg();
+            yref = pnt.dec_deg();
+        }
+
+    } // endif: got pointing from observation
+
+    // Initialise counts cube
+    GSkyMap cube = GSkyMap(proj, coordsys, xref, yref, -binsz, binsz,
+                           nxpix, nypix, ebounds.size());
+
+
+    // Return cube
+    return cube;
+}
+
+
+/***********************************************************************//**
  * @brief Fill events into counts cube
  *
  * @param[in] obs CTA observation.
@@ -641,10 +787,10 @@ void ctbin::get_parameters(void)
  *
  * Fills the events from an event list into the counts cube.
  ***************************************************************************/
-GSkyMap ctbin::fill_cube(GCTAObservation* obs)
+GSkyMap ctbin::fill_cube(const GCTAObservation* obs)
 {
-    // Make sure that the observation holds a CTA event list. If this
-    // is not the case then throw an exception.
+    // Make sure that the observation holds a CTA event list. If this is
+    // not the case then throw an exception.
     const GCTAEventList* events = dynamic_cast<const GCTAEventList*>
                                   (obs->events());
     if (events == NULL) {
@@ -655,7 +801,7 @@ GSkyMap ctbin::fill_cube(GCTAObservation* obs)
     }
  
     // Get the RoI
-    const GCTARoi& roi = events->roi();
+    const GCTARoi& roi = obs->roi();
  
     // Check for RoI sanity
     if (!roi.is_valid()) {
@@ -667,7 +813,7 @@ GSkyMap ctbin::fill_cube(GCTAObservation* obs)
     }
  
     // Get counts cube usage flags
-    std::vector<bool> usage = cube_layer_usage(m_ebounds, events->ebounds());
+    std::vector<bool> usage = cube_layer_usage(m_ebounds, obs->ebounds());
  
     // Initialise binning statistics
     int num_outside_roi  = 0;
@@ -676,9 +822,12 @@ GSkyMap ctbin::fill_cube(GCTAObservation* obs)
     int num_outside_ebds = 0;
     int num_in_map       = 0;
  
-    // Initialise new empty counts map
-    GSkyMap counts(m_counts);
-    counts *= 0.;
+    // Create counts cube
+    GSkyMap counts = create_cube(obs);
+
+    // Extract counts cube boundaries
+    double pixel_x_max = counts.nx() - 0.5;
+    double pixel_y_max = counts.ny() - 0.5;
  
     // Fill counts sky map
     for (int i = 0; i < events->size(); ++i) {
@@ -699,7 +848,7 @@ GSkyMap ctbin::fill_cube(GCTAObservation* obs)
         // Determine sky pixel
         GSkyPixel pixel;
         try {
-            pixel = m_counts.dir2pix(dir);
+            pixel = counts.dir2pix(dir);
         }
         catch (std::exception &e) {
             num_invalid_wcs++;
@@ -708,8 +857,8 @@ GSkyMap ctbin::fill_cube(GCTAObservation* obs)
  
         // Skip event if corresponding counts cube pixel is outside the
         // counts cube map range
-        if (pixel.x() < -0.5 || pixel.x() > (m_counts.nx()-0.5) ||
-            pixel.y() < -0.5 || pixel.y() > (m_counts.ny()-0.5)) {
+        if (pixel.x() < -0.5 || pixel.x() > pixel_x_max ||
+            pixel.y() < -0.5 || pixel.y() > pixel_y_max) {
             num_outside_map++;
             continue;
         }
@@ -726,7 +875,6 @@ GSkyMap ctbin::fill_cube(GCTAObservation* obs)
         }
  
         // Fill event in skymap
-        //#pragma omp critical(ctbin_fill_cube)
         counts(pixel, iebin) += 1.0;
  
         // Increment number of maps
@@ -734,23 +882,6 @@ GSkyMap ctbin::fill_cube(GCTAObservation* obs)
  
     } // endfor: looped over all events
  
-    // Update time information
-    /*#pragma omp critical(ctbin_fill_cube)
-    {
-        // If GTI is empty then set its time reference from the observation.
-        // From then on we keep that time reference
-        if (m_gti.is_empty()) {
-            m_gti.reference(events->gti().reference());
-        }
-
-        // Append GTIs
-        m_gti.extend(events->gti());
-
-        // Update ontime and livetime
-        m_ontime   += obs->ontime();
-        m_livetime += obs->livetime();
-	}*/
-
     // Log filling results
     #pragma omp critical(ctbin_fill_cube)
     {
@@ -779,21 +910,10 @@ GSkyMap ctbin::fill_cube(GCTAObservation* obs)
  * Sets the counts cube weights for all bins that are considered for the
  * specific observation to unity.
  ***************************************************************************/
-GSkyMap ctbin::set_weights(GCTAObservation* obs)
+GSkyMap ctbin::set_weights(const GCTAObservation* obs)
 {
-    // Make sure that the observation holds a CTA event list. If this
-    // is not the case then throw an exception.
-    const GCTAEventList* events = dynamic_cast<const GCTAEventList*>
-                                  (obs->events());
-    if (events == NULL) {
-        std::string msg = "CTA Observation does not contain an event "
-                          "list. An event list is needed to fill the "
-                          "counts cube.";
-        throw GException::invalid_value(G_SET_WEIGHTS, msg);
-    }
-
     // Get the RoI
-    const GCTARoi& roi = events->roi();
+    const GCTARoi& roi = obs->roi();
  
     // Check for RoI sanity
     if (!roi.is_valid()) {
@@ -804,19 +924,32 @@ GSkyMap ctbin::set_weights(GCTAObservation* obs)
         throw GException::invalid_value(G_SET_WEIGHTS, msg);
     }
  
+    // Get the RoI centre and radius in radians
+    GSkyDir roi_centre = roi.centre().dir();
+    double  roi_radius = roi.radius() * gammalib::deg2rad;
+
     // Get counts cube usage flags
-    std::vector<bool> usage = cube_layer_usage(m_ebounds, events->ebounds());
+    std::vector<bool> usage = cube_layer_usage(m_ebounds, obs->ebounds());
+
+    // Create weights cube
+    GSkyMap weights = create_cube(obs);
  
-    // Initialise new empty weights map
-    GSkyMap weights(m_weights);
-    weights *= 0.;
+    // Set sky direction cache availability flag
+    bool cache = (m_stack && m_dirs.size() == weights.npix());
  
     // Loop over all pixels in counts cube
-    for (int pixel = 0; pixel < m_counts.npix(); ++pixel) {
+    for (int pixel = 0; pixel < weights.npix(); ++pixel) {
  
         // Skip pixel if it is outside the RoI
-        if (roi.centre().dir().dist_deg(m_dirs[pixel]) > roi.radius()) {
-            continue;
+        if (cache) {
+            if (m_dirs[pixel].dist(roi_centre) > roi_radius) {
+                continue;
+            }
+        }
+        else {
+            if (roi_centre.dist(weights.inx2dir(pixel)) > roi_radius) {
+                continue;
+            }
         }
  
         // Loop over all energy layers of counts cube
@@ -828,7 +961,6 @@ GSkyMap ctbin::set_weights(GCTAObservation* obs)
             }
  
             // Signal that bin was filled
-            //#pragma omp critical(ctbin_set_weights)
             weights(pixel, iebin) = 1.0;
  
         } // endfor: looped over energy layers of counts cube
@@ -851,8 +983,11 @@ GSkyMap ctbin::set_weights(GCTAObservation* obs)
  * information in case that a single CTA observation is provided to support
  * binned analysis.
  ***************************************************************************/
-void ctbin::obs_cube(void)
+void ctbin::obs_cube_stacked(void)
 {
+    // Set event cube in first slot
+    m_cubes[0] = GCTAEventCube(m_counts[0], m_weights[0], m_ebounds, m_gti);
+
     // If we have only a single CTA observation in the container, then
     // keep that observation and just attach the event cube to it. Reset
     // the filename, otherwise we still will have the old event filename
@@ -865,9 +1000,9 @@ void ctbin::obs_cube(void)
  
             // Change the event type if we had an unbinned observation
             if (obs->eventtype() == "EventList") {
-		
+
                 // Assign cube to the observation
-		obs->events(this->cube());
+                obs->events(m_cubes[0]);
  
             }
  
@@ -877,7 +1012,8 @@ void ctbin::obs_cube(void)
             else {
  
                 // Create empty counts cube
-                GCTAEventCube cube(m_counts, m_weights, m_ebounds, obs->gti());
+                GCTAEventCube cube(m_counts[0], m_weights[0], m_ebounds,
+                                   obs->gti());
  
                 // Assign empty cube to observation
                 obs->events(cube);
@@ -899,49 +1035,49 @@ void ctbin::obs_cube(void)
         // Allocate observation container
         GObservations container;
  
-	// Allocate CTA observation.
-	GCTAObservation obs;
-	
-	// Attach event cube to CTA observation
-	obs.events(this->cube());
-	
-	// Compute average pointing direction for all CTA event lists
-	double ra     = 0.0;
-	double dec    = 0.0;
-	double number = 0.0;
-	for (int i = 0; i < m_obs.size(); ++i) {
-	    GCTAObservation* cta = dynamic_cast<GCTAObservation*>(m_obs[i]);
-	    if ((cta != NULL) && (cta->eventtype() == "EventList")) {
-		ra     += cta->pointing().dir().ra();
-		dec    += cta->pointing().dir().dec();
-		number += 1.0;
-	    }
-	}
-	if (number > 0.0) {
-	    ra  /= number;
-	    dec /= number;
-	}
-	GSkyDir dir;
-	dir.radec(ra, dec);
-	GCTAPointing pointing(dir);
-	
-	// Compute deadtime correction
-	double deadc = (m_ontime > 0.0) ? m_livetime / m_ontime : 0.0;
-	
-	// Set CTA observation attributes
-	obs.pointing(pointing);
-	obs.ra_obj(dir.ra_deg());   //!< Dummy
-	obs.dec_obj(dir.dec_deg()); //!< Dummy
-	obs.ontime(m_ontime);
-	obs.livetime(m_livetime);
-	obs.deadc(deadc);
-	
-	// Append CTA observation
-	container.append(obs);
+        // Allocate CTA observation.
+        GCTAObservation obs;
+
+        // Attach event cube to CTA observation
+        obs.events(m_cubes[0]);
+
+        // Compute average pointing direction for all CTA event lists
+        double ra     = 0.0;
+        double dec    = 0.0;
+        double number = 0.0;
+        for (int i = 0; i < m_obs.size(); ++i) {
+            GCTAObservation* cta = dynamic_cast<GCTAObservation*>(m_obs[i]);
+            if ((cta != NULL) && (cta->eventtype() == "EventList")) {
+                ra     += cta->pointing().dir().ra();
+                dec    += cta->pointing().dir().dec();
+                number += 1.0;
+            }
+        }
+        if (number > 0.0) {
+            ra  /= number;
+            dec /= number;
+        }
+        GSkyDir dir;
+        dir.radec(ra, dec);
+        GCTAPointing pointing(dir);
+
+        // Compute deadtime correction
+        double deadc = (m_ontime > 0.0) ? m_livetime / m_ontime : 0.0;
+
+        // Set CTA observation attributes
+        obs.pointing(pointing);
+        obs.ra_obj(dir.ra_deg());   //!< Dummy
+        obs.dec_obj(dir.dec_deg()); //!< Dummy
+        obs.ontime(m_ontime);
+        obs.livetime(m_livetime);
+        obs.deadc(deadc);
  
         // Set models in observation container
         container.models(m_obs.models());
  
+        // Append CTA observation
+        container.append(obs);
+
         // Copy over all remaining non-CTA observations
         for (int i = 0; i < m_obs.size(); ++i) {
             GCTAObservation* obs = dynamic_cast<GCTAObservation*>(m_obs[i]);
 /***************************************************************************
  *                        ctbin - Event binning tool                       *
  * ----------------------------------------------------------------------- *
- *  copyright (C) 2010-2017 by Juergen Knoedlseder                         *
+ *  copyright (C) 2010-2018 by Juergen Knoedlseder                         *
  * ----------------------------------------------------------------------- *
  *                                                                         *
  *  This program is free software: you can redistribute it and/or modify   *
@@ -81,10 +81,10 @@ public:
     // Methods
     void                 clear(void);
     void                 run(void);
+    int                  cubes(void) const;
+    const GCTAEventCube& cube(const int& index = 0) const;
     void                 save(void);
     void                 publish(const std::string& name = "");
-    const GCTAEventCube& cube(void) const;
-    const GCTAEventCube& cube(const size_t& index) const;
  
 protected:
     // Protected methods
@@ -92,60 +92,46 @@ protected:
     void    copy_members(const ctbin& app);
     void    free_members(void);
     void    get_parameters(void);
-    GSkyMap fill_cube(GCTAObservation* obs);
-    GSkyMap set_weights(GCTAObservation* obs);
-    void    obs_cube(void);
+    void    init_sky_dir_cache(void);
+    GSkyMap create_cube(const GCTAObservation* obs);
+    GSkyMap fill_cube(const GCTAObservation* obs);
+    GSkyMap set_weights(const GCTAObservation* obs);
+    void    obs_cube_stacked(void);
  
     // User parameters
-    GFilename     m_outcube;  //!< Output counts map file name
-    bool          m_usepnt;   //!< Use pointing instead of xref/yref parameters
-    bool          m_stack;    //!< Output one stacked cube or multiple cubes
-    std::string   m_prefix;   //!< Prefix for output path of multiple cubes
-    bool          m_publish;  //!< Publish counts cube?
-    GChatter      m_chatter;  //!< Chattiness
+    GFilename   m_outcube;  //!< Output counts map file name
+    bool        m_usepnt;   //!< Use pointing instead of xref/yref parameters
+    bool        m_stack;    //!< Output one stacked cube or multiple cubes
+    std::string m_prefix;   //!< Prefix for output path of multiple cubes
+    bool        m_publish;  //!< Publish counts cube?
+    GChatter    m_chatter;  //!< Chattiness
  
     // Protected members
-    GSkyMap       m_counts;   //!< Event cube counts
-    GSkyMap       m_weights;  //!< Event cube weights
-    GEbounds      m_ebounds;  //!< Energy boundaries
-    GGti          m_gti;      //!< Good time intervals
-    GCTAEventCube m_cube;     //!< Events cube (for cube() method)
-    double        m_ontime;   //!< Total ontime
-    double        m_livetime; //!< Total livetime
-    std::vector<GCTAEventCube> m_cubes; //!< Event cubes (for cube() method)
-    std::vector<std::string> m_ids;     //!< Observation ids
+    std::vector<GCTAEventCube> m_cubes;    //!< Event cubes
+    std::vector<GSkyMap>       m_counts;   //!< List of event cube counts
+    std::vector<GSkyMap>       m_weights;  //!< List of event cube weights
+    std::vector<std::string>   m_ids;      //!< List of observation ids
+    GEbounds                   m_ebounds;  //!< Energy boundaries
+    GGti                       m_gti;      //!< Stacked Good time intervals
+    double                     m_ontime;   //!< Total ontime
+    double                     m_livetime; //!< Total livetime
  
     // Cache members
-    std::vector<GSkyDir> m_dirs; //!< Cached GSkyDir for each pixel in m_counts
+    std::vector<GSkyDir>       m_dirs;     //!< Cached GSkyDir for all pixels
 };
  
  
 /***********************************************************************//**
- * @brief Return event cube
+ * @brief Return number of event cubes
  *
- * @return Reference to event cube
+ * @return Number of event cubes.
  *
- * Returns a reference to the event cube.
+ * Returns number of event cubes.
  ***************************************************************************/
 inline
-const GCTAEventCube& ctbin::cube(void) const
+int ctbin::cubes(void) const
 {
-    return m_cube;
-}
-
-
-/***********************************************************************//**
- * @brief Return event cube at index
- *
- * @param[in] index Index of cube
- * @return Reference to event cube
- *
- * Returns a reference to the event cube at the given index.
- ***************************************************************************/
-inline
-const GCTAEventCube& ctbin::cube(const size_t& index) const
-{
-    return m_cubes[index];
+    return m_cubes.size();
 }
  
 #endif /* CTBIN_HPP */
@@ -18,7 +18,7 @@
 #===================
 inobs,   f, a, events.fits,,, "Input event list or observation definition XML file"
 outcube, f, a, cntcube.fits,,, "Output counts cube file"
-stack,   b, h, yes,,, "Stack events to one cube?"
+stack,   b, h, yes,,, "Stack events into one cube?"
 prefix,  s, h, binned_,,, "Prefix for counts cubes in observation definition XML file"
  
 #
@@ -165,7 +165,7 @@ class Test(test):
         # cube
         self.test_value(binning.obs().size(), 0,
              'Check that empty ctbin has an empty observation container')
-        self.test_value(binning.cube().size(), 0,
+        self.test_value(binning.cubes(), 0,
              'Check that empty ctbin has an empty counts cube')
  
         # Check that saving does not nothing
@@ -220,11 +220,13 @@ class Test(test):
         cpy_bin.logFileOpen()
         cpy_bin.run()
  
-        # Check content of observation container and counts cube. Now an empty
-        # event cube is expected since on input the observation is binned, and
-        # any binned observation will be skipped.
-        self._check_observation(cpy_bin, 0)
-        self._check_cube(cpy_bin.cube(), 0)
+        # Check content of observation container and number of counts cubes.
+        # There should be a single binned observation in the observation
+        # container, which is the one that was produced in the run before.
+        # Since ctbin finds no unbinned observation in the container, the
+        # number of cubes should be zero.
+        self._check_observation(cpy_bin, 245)
+        self.test_value(cpy_bin.cubes(), 0, 'Check that there are no counts cubes')
  
         # Save counts cube
         binning.save()
@@ -240,7 +242,7 @@ class Test(test):
         # counts cube
         self.test_value(binning.obs().size(), 0,
              'Check that empty ctbin has an empty observation container')
-        self.test_value(binning.cube().size(), 0,
+        self.test_value(binning.cubes(), 0,
              'Check that empty ctbin has an empty counts cube')
  
         # Prepare observation container for stacking of events into a