GCOSEventCube.cpp

Go to the documentation of this file.

/***************************************************************************
 *            GCOSEventCube.cpp - COSI event bin container class           *
 * ----------------------------------------------------------------------- *
 *  copyright (C) 2026 by Juergen Knoedlseder                              *
 * ----------------------------------------------------------------------- *
 *                                                                         *
 *  This program is free software: you can redistribute it and/or modify   *
 *  it under the terms of the GNU General Public License as published by   *
 *  the Free Software Foundation, either version 3 of the License, or      *
 *  (at your option) any later version.                                    *
 *                                                                         *
 *  This program is distributed in the hope that it will be useful,        *
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of         *
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the          *
 *  GNU General Public License for more details.                           *
 *                                                                         *
 *  You should have received a copy of the GNU General Public License      *
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.  *
 *                                                                         *
 ***************************************************************************/
/**
 * @file GCOSEventCube.hpp
 * @brief COSI event bin container class implementation
 * @author Juergen Knoedlseder
 */
 
/* __ Includes ___________________________________________________________ */
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include "GBounds.hpp"
#include "GEbounds.hpp"
#include "GSkyMap.hpp"
#include "GCOSEventCube.hpp"
#include "GCOSEventList.hpp"
 
/* __ Method name definitions ____________________________________________ */
#define G_NAXIS                                   "GCOSEventCube::naxis(int)"
#define G_SET         "GCOSEventCube::set(std::string&, int&, std::string&, "\
                                                       "GBounds&, GEbounds&)"
#define G_FILL                          "GCOSEventCube::fill(GCOSEventList&)"
#define G_SET_SCATTER_DIRECTIONS    "GCOSEventCube::set_scatter_directions()"
#define G_SET_SCATTER_ANGLES            "GCOSEventCube::set_scatter_angles()"
#define G_SET_ENERGIES                        "GCOSEventCube::set_energies()"
#define G_SET_TIMES                              "GCOSEventCube::set_times()"
#define G_SET_BIN                              "GCOSEventCube::set_bin(int&)"
 
/* __ Macros _____________________________________________________________ */
 
/* __ Coding definitions _________________________________________________ */
 
/* __ Debug definitions __________________________________________________ */
 
 
 
/*==========================================================================
 =                                                                         =
 =                        Constructors/destructors                         =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Void constructor
 *
 * Constructs an empty COSI event cube.
 ***************************************************************************/
GCOSEventCube::GCOSEventCube(void) : GEventCube()
{
    // Initialise members
    init_members();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Cube constructor
 *
 * @param[in] coords Event cube coordinate system (CEL or GAL).
 * @param[in] nside Event cube Nside parameter.
 * @param[in] order Event cube pixel ordering (RING or NEST).
 * @param[in] phibounds Phi layer boundaries.
 * @param[in] ebounds Event cube energy boundaries.
 * @param[in] gti Good Time intervals.
 *
 * Constructs a HealPix event cube by specifying the coordinate system, the
 * HealPix Nside parameter and pixel ordering scheme, the event cube Phi
 * layer boundaries, the event cube energy boundaries and the event cube
 * Good Time Invervals.
 ***************************************************************************/
GCOSEventCube::GCOSEventCube(const std::string& coords,
                             const int&         nside,
                             const std::string& order,
                             const GBounds&     phibounds,
                             const GEbounds&    ebounds,
                             const GGti&        gti) : GEventCube()
{
    // Initialise members
    init_members();
 
    // Set event cube
    this->set(coords, nside, order, phibounds, ebounds, gti);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Load constructor
 *
 * @param[in] filename COSI event cube FITS filename.
 *
 * Construct a COSI event cube by loading it from a FITS file.
 ***************************************************************************/
GCOSEventCube::GCOSEventCube(const GFilename& filename) : GEventCube()
{
    // Initialise members
    init_members();
 
    // Load event cube
    load(filename);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Copy constructor
 *
 * @param[in] cube COSI event cube.
 ***************************************************************************/
GCOSEventCube::GCOSEventCube(const GCOSEventCube& cube) : GEventCube(cube)
{
    // Initialise members
    init_members();
 
    // Copy members
    copy_members(cube);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Destructor
 ***************************************************************************/
GCOSEventCube::~GCOSEventCube(void)
{
    // Free members
    free_members();
 
    // Return
    return;
}
 
 
/*==========================================================================
 =                                                                         =
 =                               Operators                                 =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Assignment operator
 *
 * @param[in] cube COSI event cube.
 * @return COSI event cube.
 ***************************************************************************/
GCOSEventCube& GCOSEventCube::operator=(const GCOSEventCube& cube)
{
    // Execute only if object is not identical
    if (this != &cube) {
 
        // Copy base class members
        this->GEventCube::operator=(cube);
 
        // Free members
        free_members();
 
        // Initialise members
        init_members();
 
        // Copy members
        copy_members(cube);
 
    } // endif: object was not identical
 
    // Return this object
    return *this;
}
 
 
/***********************************************************************//**
 * @brief Event bin access operator
 *
 * @param[in] index Event index [0,...,size()-1].
 * @return Pointer to event bin.
 *
 * Returns pointer to an event bin. Note that the returned pointer is in
 * fact always the same, but the method sets the pointers within the
 * event bin so that they point to the appropriate information.
 ***************************************************************************/
GCOSEventBin* GCOSEventCube::operator[](const int& index)
{
    // Set event bin
    set_bin(index);
 
    // Return pointer
    return (&m_bin);
}
 
 
/***********************************************************************//**
 * @brief Event bin access operator (const version)
 *
 * @param[in] index Event index [0,...,size()-1].
 * @return Const pointer to event bin.
 *
 * Returns pointer to an event bin. Note that the returned pointer is in
 * fact always the same, but the method sets the pointers within the
 * event bin so that they point to the appropriate information.
 ***************************************************************************/
const GCOSEventBin* GCOSEventCube::operator[](const int& index) const
{
    // Set event bin (circumvent const correctness)
    const_cast<GCOSEventCube*>(this)->set_bin(index);
 
    // Return pointer
    return (&m_bin);
}
 
 
/*==========================================================================
 =                                                                         =
 =                             Public methods                              =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Clear COSI event cube
 *
 * Clears COSI event cube by resetting all class members to an
 * initial state. Any information that was present before will be lost.
 ***************************************************************************/
void GCOSEventCube::clear(void)
{
    // Free class members (base and derived classes, derived class first)
    free_members();
    this->GEventCube::free_members();
    this->GEvents::free_members();
 
    // Initialise members
    this->GEvents::init_members();
    this->GEventCube::init_members();
    init_members();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Clone event cube
 *
 * @return Pointer to deep copy of COSI event cube.
 ***************************************************************************/
GCOSEventCube* GCOSEventCube::clone(void) const
{
    return new GCOSEventCube(*this);
}
 
 
/***********************************************************************//**
 * @brief Return number of bins in event cube
 *
 * @return Number of bins in event cube.
 ***************************************************************************/
int GCOSEventCube::size(void) const
{
    // Compute number of bins
    int nbins = m_cube.npix() * m_cube.nmaps();
 
    // Return number of bins
    return nbins;
}
 
 
/***********************************************************************//**
 * @brief Return dimension of event cube
 *
 * @return Number of dimensions in event cube.
 *
 * The number of dimensions of the event cube counts one for the sky map
 * pixels, one for the phibar axis and one for the energy axis. So typically
 * an event cube has 3 dimensions.
 ***************************************************************************/
int GCOSEventCube::dim(void) const
{
    // Compute dimension
    int dim = (m_cube.is_empty()) ? 0 : m_cube.ndim() + 1;
 
    // Return dimension
    return dim;
}
 
 
/***********************************************************************//**
 * @brief Return number of bins in axis
 *
 * @param[in] axis Axis [0,...,dim()-1].
 * @return Number of bins in axis.
 *
 * @exception GException::out_of_range
 *            Axis is out of range.
 *
 * Returns the number of bins along a given event cube @p axis.
 * If @p axis=0 the number if HealPix pixels is returned.
 * If @p axis=1 the number of Phi laters is returned.
 * If @p axis=2 the number of energy bins is returned.
 ***************************************************************************/
int GCOSEventCube::naxis(const int& axis) const
{
    // Optionally check if the axis is valid
    #if defined(G_RANGE_CHECK)
    if (axis < 0 || axis >= dim()) {
        throw GException::out_of_range(G_NAXIS, "COSI event cube axis",
                                       axis, dim());
    }
    #endif
 
    // Set result
    int naxis = (axis == 0) ? m_cube.npix() : m_cube.shape()[axis-1];
 
    // Return result
    return naxis;
}
 
 
/***********************************************************************//**
 * @brief Load COSI event cube from FITS file
 *
 * @param[in] filename FITS file name.
 *
 * The method clears the object before loading, thus any events residing in
 * the object before loading will be lost.
 ***************************************************************************/
void GCOSEventCube::load(const GFilename& filename)
{
    // Open FITS file
    GFits fits(filename);
 
    // Read event cube from FITS file
    read(fits);
 
    // Close FITS file
    fits.close();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Save COSI event cube into FITS file
 *
 * @param[in] filename FITS file name.
 * @param[in] clobber Overwrite existing FITS file?
 *
 * Saves the COSI event cube into a FITS file.
 ***************************************************************************/
void GCOSEventCube::save(const GFilename& filename, const bool& clobber) const
{
    // Create empty FITS file
    GFits fits;
 
    // Write event cube into FITS file
    write(fits);
 
    // Save FITS file
    fits.saveto(filename, clobber);
 
    // Close FITS file
    fits.close();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Read COSI event cube from FITS file
 *
 * @param[in] fits FITS file.
 *
 * Reads a COSI event cube from a FITS file.
 ***************************************************************************/
void GCOSEventCube::read(const GFits& fits)
{
    // Clear object
    clear();
 
    // Get event cube HDU and Phi boundaries table
    const GFitsHDU&   hdu       = *fits["EVENT CUBE"];
    const GFitsTable& phibounds = *fits.table("PHI BOUNDARIES");
    const GFitsTable& ebounds   = *fits.table(gammalib::extname_ebounds);
    const GFitsTable& gti       = *fits.table(gammalib::extname_gti);
 
    // Read event cube
    m_cube.read(hdu);
 
    // Read Phi boundaries
    m_phibounds.read(phibounds);
 
    // Read energy boundaries
    m_ebounds.read(ebounds);
 
    // Read Good Time Intervals
    m_gti.read(gti);
 
    // Read event selection statistics
    if (hdu.has_card("NEVCHK")) {
        m_num_events_checked = hdu.integer("NEVCHK");
    }
    if (hdu.has_card("NEVUSE")) {
        m_num_events_used = hdu.integer("NEVUSE");
    }
    if (hdu.has_card("NEVGTI")) {
        m_num_outside_gti = hdu.integer("NEVGTI");
    }
    if (hdu.has_card("NEVEBD")) {
        m_num_outside_ebounds = hdu.integer("NEVEBD");
    }
    if (hdu.has_card("NEVPBD")) {
        m_num_outside_phibounds = hdu.integer("NEVPBD");
    }
 
    // Set event cube shape
    m_cube.shape(m_phibounds.size(), m_ebounds.size());
 
    // Set scatter directions
    set_scatter_directions();
 
    // Set scatter angles
    set_scatter_angles();
 
    // Set energies
    set_energies();
 
    // Set times
    set_times();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Write COSI event cube into FITS file.
 *
 * @param[in] fits FITS file.
 *
 * Writes the COSI event cube into a FITS file. The method does nothing if
 * size() is not positive.
 ***************************************************************************/
void GCOSEventCube::write(GFits& fits) const
{
    // Continue only if size is positive
    if (size() > 0) {
 
        // Write event cube
        m_cube.write(fits, "EVENT CUBE");
 
        // Write Phi boundaries
        m_phibounds.write(fits, "PHI BOUNDARIES");
 
        // Write energy boundaries
        m_ebounds.write(fits);
 
        // Write Good Time Intervals
        m_gti.write(fits);
 
        // Optionally write event selection statistics into event cube
        if (m_num_events_checked > 0) {
 
            // Get event cube HDU
            GFitsHDU* hdu = fits["EVENT CUBE"];
 
            // Derive number of rejected events
            int num_events_rejected = m_num_events_checked - m_num_events_used;
 
            // Write event selection statistics
            hdu->card("NEVCHK", m_num_events_checked,    "Number of checked events");
            hdu->card("NEVUSE", m_num_events_used,       "Number of used events");
            hdu->card("NEVREJ", num_events_rejected,     "Number of rejected events");
            hdu->card("NEVGTI", m_num_outside_gti,       "Number of events outside GTI");
            hdu->card("NEVEBD", m_num_outside_ebounds,   "Number of events outside energy boundaries");
            hdu->card("NEVPBD", m_num_outside_phibounds, "Number of events outside Phi boundaries");
 
        }
 
    } // endif: size was positive
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Return number of events in cube
 *
 * @return Number of events in event cube.
 *
 * This method returns the number of events in the event cube.
 ***************************************************************************/
double GCOSEventCube::number(void) const
{
    // Initialise number of events
    double number = 0.0;
 
    // Get pointer to event cube bins
    const double *ptr = m_cube.pixels();
 
    // Loop over event cube and sum events
    for (int i = 0; i < size(); ++i, ++ptr) {
        number += *ptr;
    }
 
    // Return
    return number;
}
 
 
/***********************************************************************//**
 * @brief Set event cube
 *
 * @param[in] coords Event cube coordinate system (CEL or GAL).
 * @param[in] nside Event cube Nside parameter.
 * @param[in] order Event cube pixel ordering (RING or NEST).
 * @param[in] phibounds Phi layer boundaries.
 * @param[in] ebounds Event cube energy boundaries.
 * @param[in] gti Good Time intervals.
 *
 * @exception GException::invalid_argument
 *            Invalid units specified for Phi layer boundaries.
 *
 * Sets a HealPix event cube by specifying the coordinate system, the
 * HealPix Nside parameter and pixel ordering scheme, the event cube Phi
 * layer boundaries, the event cube energy boundaries and the event cube
 * Good Time Intervals.
 ***************************************************************************/
void GCOSEventCube::set(const std::string& coords,
                        const int&         nside,
                        const std::string& order,
                        const GBounds&     phibounds,
                        const GEbounds&    ebounds,
                        const GGti&        gti)
{
    // Clear event cube
    clear();
 
    // Compute number of maps in event cube
    int nmaps = phibounds.size() * ebounds.size();
 
    // Allocate event cube
    m_cube = GSkyMap(coords, nside, order, nmaps);
 
    // Set event cube shape
    m_cube.shape(phibounds.size(), ebounds.size());
 
    // Make sure that Phi boundaries are in degrees
    std::string unit = gammalib::tolower(phibounds.unit());
    if ((unit != "") && ((unit != "deg") || (unit != "degrees"))) {
        std::string msg = "Invalid unit \""+phibounds.unit()+"\"specified for "
                          "Phi boundaries. Please specify Phi boundaries in "
                          "units of \"deg\".";
        throw GException::invalid_argument(G_SET, msg);
    }
 
    // Store Phi boundaries and set units to deg
    m_phibounds = phibounds;
    m_phibounds.unit("deg");
 
    // Store energy boundaries
    m_ebounds = ebounds;
 
    // Store Good Time intervals
    m_gti = gti;
 
    // Set scatter directions
    set_scatter_directions();
 
    // Set scatter angles
    set_scatter_angles();
 
    // Set energies
    set_energies();
 
    // Set times
    set_times();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Fill event cube with events from event list
 *
 * @param[in] list COSI event list.
 *
 * @exception GException::invalid_value
 *            No Good Time Intervals defined.
 *            No energy boundaries defined.
 *            No scatter angle boundaries defined.
 *
 * Fills events from an event list into the event cube. Events that existed
 * in the event cube will be preserved so that events from multiple event
 * lists can be filled into the event cube.
 *
 * The method requires that Good Time Intervals, energy boundaries and
 * scatter angle boundaries are defined. If boundaries are not defined an
 * exception will be thrown.
 *
 * Only events that fall within the Good Time Intervals, energy boundaries
 * and scatter angle boundaries will be filled into the cube, other events
 * will be ignored.
 *
 * The method assumes that the event list and the Good Time Intervals are
 * sorted by increasing time. This allows for a fast check of event times
 * against Good Time Intervals.
 *
 * The method updates the event filling statistics.
 ***************************************************************************/
void GCOSEventCube::fill(const GCOSEventList& list)
{
    // Check that boundaries are defined
    if (gti().is_empty()) {
        std::string msg = "No Good Time Intervals were defined for the event "
                          "cube. Please set the Good Time Intervals before "
                          "calling the method.";
        throw GException::invalid_value(G_FILL, msg);
    }
    if (ebounds().is_empty()) {
        std::string msg = "No energy boundaries were defined for the event "
                          "cube. Please set the energy boundaries before "
                          "calling the method.";
        throw GException::invalid_value(G_FILL, msg);
    }
    if (phibounds().is_empty()) {
        std::string msg = "No scatter angle boundaries were defined for the "
                          "event cube. Please set the scatter angle "
                          "boundaries before calling the method.";
        throw GException::invalid_value(G_FILL, msg);
    }
 
    // Initialise Good Time Interval index
    int igti = 0;
 
    // Loop over events
    for (int i = 0; i < list.size(); ++i) {
 
        // Get event and attributes
        const GCOSEventAtom* event = list[i];
        const GTime&         time  = event->time();
 
        // Update selection statistics
        m_num_events_checked++;
 
        // If event is before current Good Time Interval then skip event
        if (time < gti().tstart(igti)) {
            m_num_outside_gti++;
            continue;
        }
 
        // If event time is after current Good Time Interval then search
        // for the next Good Time Interval that has a stop time after the
        // event time. If such an interval is found but the event time is
        // not in the interval the drop the event. If no such interval was
        // found then exit the loop.
        if (time >= gti().tstop(igti)) {
            bool found = false;                       // Signal that no GTI was found
            igti++;                                   // Move to next GTI
            while (igti < gti().size()) {             // Search GTIs until exhaustion
                if (time < gti().tstop(igti)) {       // If event is before stop...
                    if (time >= gti().tstart(igti)) { // then check if after start...
                        found = true;                 // and if so we found a new GTI
                    }
                    break;                            // Break since event is before GTI stop
                }
                igti++;                               // Move to next GTI
            }
            if (!found) {                             // No GTI found
                if (igti < gti().size()) {            // If GTIs are not exhausted
                    m_num_outside_gti++;              // the event is dropped
                }
                else {                                // otherwise we exit the event loop
                    m_num_events_checked += list.size() - i - 1;
                    break;
                }
            }
        }
 
        // Get energy boundary index
        int ieng = m_ebounds.index(event->energy());
 
        // Skip event if it does not fall into energy boundaries
        if (ieng == -1) {
            m_num_outside_ebounds++;
            continue;
        }
 
        // Get Phi boundary index
        int iphi = m_phibounds.index(event->dir().phi());
 
        // Skip event if it does not fall into Phi boundaries
        if (iphi == -1) {
            m_num_outside_phibounds++;
            continue;
        }
 
        // If we survived until now the event will be used
        m_num_events_used++;
 
        // Compute map index
        int imap = iphi + ieng * m_phibounds.size();
 
        // Get skymap pixel index
        int ipix = m_cube.dir2inx(event->dir().dir());
 
        // Increment event cube
        m_cube(ipix, imap) += 1.0;
 
    } // endfor: looped over events
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Print COSI event cube information
 *
 * @param[in] chatter Chattiness.
 * @return String containing COSI event cube information.
 ***************************************************************************/
std::string GCOSEventCube::print(const GChatter& chatter) const
{
    // Initialise result string
    std::string result;
 
    // Continue only if chatter is not silent
    if (chatter != SILENT) {
 
        // Append header
        result.append("=== GCOSEventCube ===");
 
        // Append information
        result.append("\n"+gammalib::parformat("Number of events"));
        result.append(gammalib::str(number()));
        result.append("\n"+gammalib::parformat("Number of elements"));
        result.append(gammalib::str(size()));
        result.append("\n"+gammalib::parformat("Number of pixels"));
        result.append(gammalib::str(m_npix));
        result.append("\n"+gammalib::parformat("Number of Phi layers"));
        result.append(gammalib::str(m_phibounds.size()));
        result.append("\n"+gammalib::parformat("Number of energy bins"));
        result.append(gammalib::str(m_ebounds.size()));
        result.append("\n"+gammalib::parformat("Phi boundaries"));
        result.append(gammalib::str(m_phibounds.min())+" - ");
        result.append(gammalib::str(m_phibounds.max())+" ");
        result.append(m_phibounds.unit());
        result.append("\n"+gammalib::parformat("Energy boundaries"));
        result.append(gammalib::str(emin().keV())+" - ");
        result.append(gammalib::str(emax().keV())+" keV");
        result.append("\n"+gammalib::parformat("Ontime"));
        result.append(gammalib::str(m_ontime)+" s");
        result.append("\n"+gammalib::parformat("Mean time"));
        result.append(m_time.print(gammalib::reduce(chatter)));
 
        // Append GTI interval
        result.append("\n"+gammalib::parformat("Good time intervals"));
        result.append(gammalib::str(gti().size()));
        result.append("\n"+gammalib::parformat("Time interval"));
        if (gti().size() > 0) {
            result.append(gammalib::str(tstart().mjd()));
            result.append(" - ");
            result.append(gammalib::str(tstop().mjd())+" Modified Julian days");
        }
        else {
            result.append("not defined");
        }
        result.append("\n"+gammalib::parformat("Date interval"));
        if (gti().size() > 0) {
            result.append(tstart().utc());
            result.append(" - ");
            result.append(tstop().utc());
        }
        else {
            result.append("not defined");
        }
 
    } // endif: chatter was not silent
 
    // Return result
    return result;
}
 
 
/*==========================================================================
 =                                                                         =
 =                             Private methods                             =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Initialise class members
 ***************************************************************************/
void GCOSEventCube::init_members(void)
{
    // Initialise members
    m_cube.clear();
    m_phibounds.clear();
    m_bin.clear();
    m_dir.clear();
    m_time.clear();
    m_ontime     = 0.0;
    m_solidangle = 0.0;
    m_npix       = 0;
    m_nphi       = 0;
    m_dirs.clear();
    m_phis.clear();
    m_phiwidths.clear();
    m_energies.clear();
    m_ewidths.clear();
 
    // Initialise filling selection statistics
    m_num_events_checked    = 0;
    m_num_events_used       = 0;
    m_num_outside_gti       = 0;
    m_num_outside_ebounds   = 0;
    m_num_outside_phibounds = 0;
 
    // Prepare event bin
    init_bin();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Copy class members
 *
 * @param[in] cube COSI event cube.
 *
 * This method copies the class members from another event cube in the actual
 * object. It also prepares the event bin member that will be returned in
 * case of an operator access to the class.
 ***************************************************************************/
void GCOSEventCube::copy_members(const GCOSEventCube& cube)
{
    // Copy members. Note that the event bin is copied as it will be
    // initialised by the init_bin() method. The event bin serves just as a
    // container of pointers, hence we do not want to copy over the pointers
    // from the original class.
    m_cube       = cube.m_cube;
    m_phibounds  = cube.m_phibounds;
    m_dir        = cube.m_dir;
    m_time       = cube.m_time;
    m_ontime     = cube.m_ontime;
    m_solidangle = cube.m_solidangle;
    m_npix       = cube.m_npix;
    m_nphi       = cube.m_nphi;
    m_dirs       = cube.m_dirs;
    m_phis       = cube.m_phis;
    m_phiwidths  = cube.m_phiwidths;
    m_energies   = cube.m_energies;
    m_ewidths    = cube.m_ewidths;
 
    // Copy filling selection statistics
    m_num_events_checked    = cube.m_num_events_checked;
    m_num_events_used       = cube.m_num_events_used;
    m_num_outside_gti       = cube.m_num_outside_gti;
    m_num_outside_ebounds   = cube.m_num_outside_ebounds;
    m_num_outside_phibounds = cube.m_num_outside_phibounds;
 
    // Prepare event bin
    init_bin();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Delete class members
 ***************************************************************************/
void GCOSEventCube::free_members(void)
{
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Set scatter directions of all events cube pixels
 *
 * @exception GException::invalid_value
 *            No sky pixels have been defined.
 *
 * Computes the scatter directions of all events cube pixels. This method
 * sets the m_dirs, m_npix and m_solidangle members of the class.
 ***************************************************************************/
void GCOSEventCube::set_scatter_directions(void)
{
    // Clear scatter directions
    m_dirs.clear();
 
    // Store number of pixels
    m_npix = m_cube.npix();
 
    // Throw an exception if we have no sky pixels
    if (m_npix < 1) {
        std::string msg = "No sky pixels were defined for event cube. "
                          "Please define the sky pixels of the event "
                          "cube.";
        throw GException::invalid_value(G_SET_SCATTER_DIRECTIONS, msg);
    }
 
    // Reserve pixel scatter direction
    m_dirs.reserve(m_npix);
 
    // Set pixel scatter direction
    for (int i = 0; i < m_npix; ++i) {
        m_dirs[i] = m_cube.inx2dir(i);
    }
 
    // Store solid angle of pixels (HealPix map solid angles are constant)
    m_solidangle = m_cube.solidangle(0);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Set Compton scatter angles of event cube
 *
 * @exception GException::invalid_value
 *            No scatter angles have been defined.
 *
 * Computes the Compton scatter angles for all Phibar layers. This method
 * sets the m_phis, m_phiwidths and m_nphi members of the class.
 ***************************************************************************/
void GCOSEventCube::set_scatter_angles(void)
{
    // Clear scatter angles and widths
    m_phis.clear();
    m_phiwidths.clear();
 
    // Store number of Phi boundaries
    m_nphi = m_phibounds.size();
 
    // Throw an exception if we have no scatter angles
    if (m_nphi < 1) {
        std::string msg = "No scatter angles were defined for event cube. "
                          "Please define the scatter angles of the event "
                          "cube.";
        throw GException::invalid_value(G_SET_SCATTER_ANGLES, msg);
    }
 
    // Reserve space for scatter angles and widths
    m_phis.reserve(m_nphi);
    m_phiwidths.reserve(m_nphi);
 
    // Set mean scatter angles and scatter angle bin widths
    for (int i = 0; i < m_nphi; ++i) {
        m_phis[i]      = m_phibounds.mean(i);
        m_phiwidths[i] = m_phibounds.max(i) - m_phibounds.min(i);
    }
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Set log mean energies and energy widths of energy bins
 *
 * @exception GException::invalid_value
 *            No energy boundaries defined.
 *
 * Computes the long mean energy and energy bin width of the event cube.
 ***************************************************************************/
void GCOSEventCube::set_energies(void)
{
    // Clear energies and energy widths
    m_energies.clear();
    m_ewidths.clear();
 
    // Throw an exception if energy boundaries are empty
    if (m_ebounds.size() < 1) {
        std::string msg = "No energy boundaries were defined for event cube. "
                          "Please define the energy boundaries of the event "
                          "cube.";
        throw GException::invalid_value(G_SET_ENERGIES, msg);
    }
 
    // Reserve energies and energy widths
    m_energies.reserve(m_ebounds.size());
    m_ewidths.reserve(m_ebounds.size());
 
    // Loop over energy bins
    for (int i = 0; i < m_ebounds.size(); ++i) {
        m_energies[i] = m_ebounds.elogmean(i);
        m_ewidths[i]  = m_ebounds.emax(i) - m_ebounds.emin(i);
    }
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Set mean event time and ontime of event cube
 *
 * @exception GException::invalid_value
 *            No Good Time Intervals found.
 *
 * Computes the mean time of the event cube by taking the mean between start
 * and stop time. Computes also the ontime by summing up of all good time
 * intervals.
 ***************************************************************************/
void GCOSEventCube::set_times(void)
{
    // Throw an exception if GTI is empty
    if (m_gti.size() < 1) {
        std::string msg = "No Good Time Intervals were defined for event "
                          "cube. Please define the Good Time Intervals of "
                          "the event cube.";
        throw GException::invalid_value(G_SET_TIMES, msg);
    }
 
    // Compute mean time
    m_time = m_gti.tstart() + 0.5 * (m_gti.tstop() - m_gti.tstart());
 
    // Set ontime
    m_ontime = m_gti.ontime();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Initialise event bin
 *
 * Initialises the event bin. The event bin is cleared and all fixed pointers
 * are set. Only the m_counts, m_phiwidth, m_energy and m_ewidth members of
 * the event bin will be set to NULL as those will change for individual bins
 * and set by the set_bin() method. The set_bin() method is called before any
 * event bin access.
 ***************************************************************************/
void GCOSEventCube::init_bin(void)
{
    // Prepare event bin
    m_bin.free_members();
    m_bin.m_counts     = NULL;           //!< Will be set by set_bin method
    m_bin.m_dir        = &m_dir;         //!< Content will be set by set_bin method
    m_bin.m_solidangle = &m_solidangle;  //!< Fixed content
    m_bin.m_time       = &m_time;        //!< Fixed content
    m_bin.m_ontime     = &m_ontime;      //!< Fixed content
    m_bin.m_phiwidth   = NULL;           //!< Will be set by set_bin method
    m_bin.m_energy     = NULL;           //!< Will be set by set_bin method
    m_bin.m_ewidth     = NULL;           //!< Will be set by set_bin method
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Set event bin
 *
 * @param[in] index Event index [0,...,size()[.
 *
 * @exception GException::out_of_range
 *            Event index is outside valid range.
 *
 * Sets the pointers of the event bin to the event cube cell that corresponds
 * to the specified @p index. The event bin is in fact physically stored in
 * the event cube, and only a single event bin is indeed allocated. This
 * method sets up the pointers in the event bin so that a client can easily
 * access the information of individual bins as if they were stored in an
 * array.
 *
 * The method updates the following members of the event bin: m_index,
 * m_counts, m_phiwidth, m_energy and m_ewidth.
 *
 * The method also updates the event bin direction m_dir.
 ***************************************************************************/
void GCOSEventCube::set_bin(const int& index)
{
    // Optionally check if the index is valid
    #if defined(G_RANGE_CHECK)
    if (index < 0 || index >= size()) {
        throw GException::out_of_range(G_SET_BIN, "Event index", index, size());
    }
    #endif
 
    // Get pixel, phi layer and energy bin indices.
    int ipix = index % m_npix;
    int imap = index / m_npix;
    int iphi = imap  % m_nphi;
    int ieng = imap  / m_nphi;
 
    // Set event bin member
    m_bin.m_index    = index;
    m_bin.m_counts   = &((const_cast<double*>(m_cube.pixels()))[index]);
    m_bin.m_phiwidth = &(m_phiwidths[iphi]);
    m_bin.m_energy   = &(m_energies[ieng]);
    m_bin.m_ewidth   = &(m_ewidths[ieng]);
 
    // Set instrument direction member
    m_dir.dir(m_dirs[ipix]);
    m_dir.phi(m_phis[iphi]);
 
    // Return
    return;
}