GCOSResponse.cpp

Go to the documentation of this file.

/***************************************************************************
 *                  GCOSResponse.cpp - COSI response 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 GCOSResponse.hpp
 * @brief COSI instrument response function class implementation
 * @author Juergen Knoedlseder
 */
 
/* __ Includes ___________________________________________________________ */
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <cstdio>            // std::fopen, std::fgets, std::fclose, etc...
#include <string>
#include <typeinfo>
#include "GException.hpp"
#include "GSource.hpp"
#include "GModelSpatialPointSource.hpp"
#include "GFits.hpp"
#include "GFitsBinTable.hpp"
#include "GFitsTableLongCol.hpp"
#include "GFitsTableStringCol.hpp"
#include "GFitsTableDoubleCol.hpp"
#include "GFitsTableUShortCol.hpp"
#include "GCOSResponse.hpp"
#include "GCOSObservation.hpp"
 
/* __ Method name definitions ____________________________________________ */
#define G_IRF           "GCOSResponse::irf(GEvent&, GSource&, GObservation&)"
#define G_NROI            "GCOSResponse::nroi(GModelSky&, GEnergy&, GTime&, "\
                                                             "GObservation&)"
#define G_EBOUNDS                           "GCOSResponse::ebounds(GEnergy&)"
#define G_LOAD                               "GCOSResponse::load(GFilename&)"
#define G_LOAD_RESPONSE_CHUNK       "GCOSResponse::load_response_chunk(int&)"
#define G_READ_BOUNDS       "GCOSResponse::read_bounds(FILE*, GXmlElement*, "\
                                  "GXmlElement*, GXmlElement*, std::string&)"
#define G_READ_EBOUNDS     "GCOSResponse::read_ebounds(FILE*, GXmlElement*, "\
                                  "GXmlElement*, GXmlElement*, std::string&)"
#define G_READ_NDARRAY     "GCOSResponse::read_ndarray(FILE*, GXmlElement*, "\
                                                "GXmlElement*, GXmlElement*)"
#define G_READ_RESPONSE_MATRIX   "GCOSResponse::read_response_matrix(FILE*, "\
                    "GXmlElement*, GXmlElement*, GXmlElement*, GXmlElement*)"
#define G_GET_RESPONSE_CHUNK_KEY "GCOSResponse::get_response_chunk_key(int&)"
 
/* __ Macros _____________________________________________________________ */
 
/* __ Coding definitions _________________________________________________ */
 
/* __ Debug definitions __________________________________________________ */
 
/* __ Constants __________________________________________________________ */
 
/* __ Using of namespaces ________________________________________________ */
using namespace gammalib::hdf5;
 
 
/*==========================================================================
 =                                                                         =
 =                       Constructors/destructors                          =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Void constructor
 *
 * Creates an empty COSI response.
 ***************************************************************************/
GCOSResponse::GCOSResponse(void) : GResponse()
{
    // Initialise members
    init_members();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Load constructor
 *
 * @param[in] filename COSI response file.
 *
 * Create COSI response by loading data from a response file.
 ***************************************************************************/
GCOSResponse::GCOSResponse(const GFilename& filename)
{
    // Initialise members
    init_members();
 
    // Load file
    load(filename);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Copy constructor
 *
 * @param[in] rsp COSI response.
 **************************************************************************/
GCOSResponse::GCOSResponse(const GCOSResponse& rsp) : GResponse(rsp)
{
    // Initialise members
    init_members();
 
    // Copy members
    copy_members(rsp);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Destructor
 ***************************************************************************/
GCOSResponse::~GCOSResponse(void)
{
    // Free members
    free_members();
 
    // Return
    return;
}
 
 
/*==========================================================================
 =                                                                         =
 =                                Operators                                =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Assignment operator
 *
 * @param[in] rsp COSI response.
 * @return COSI response.
 *
 * Assign COSI response to this object. The assignment performs
 * a deep copy of all information, hence the original object from which the
 * assignment was made can be destroyed after this operation without any loss
 * of information.
 ***************************************************************************/
GCOSResponse& GCOSResponse::operator=(const GCOSResponse& rsp)
{
    // Execute only if object is not identical
    if (this != &rsp) {
 
        // Copy base class members
        this->GResponse::operator=(rsp);
 
        // Free members
        free_members();
 
        // Initialise members
        init_members();
 
        // Copy members
        copy_members(rsp);
 
    } // endif: object was not identical
 
    // Return this object
    return *this;
}
 
 
/*==========================================================================
 =                                                                         =
 =                             Public methods                              =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Clear instance
 *
 * Clears COSI response by resetting all class members to an initial
 * state. Any information that was present before will be lost.
 ***************************************************************************/
void GCOSResponse::clear(void)
{
    // Free class members (base and derived classes, derived class first)
    free_members();
    this->GResponse::free_members();
 
    // Initialise members
    this->GResponse::init_members();
    init_members();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Clone instance
 *
 * @return Pointer to deep copy of COSI response.
 ***************************************************************************/
GCOSResponse* GCOSResponse::clone(void) const
{
    return new GCOSResponse(*this);
}
 
 
/***********************************************************************//**
 * @brief Return value of COSI instrument response for a photon
 *
 * @param[in] event Observed event.
 * @param[in] photon Incident photon.
 * @param[in] obs Observation.
 * @return Instrument response $\f(cm^2 sr^{-1})$\f
 *
 * @exception GException::invalid_argument
 *            Observation is not a COSI observation.
 *
 * Returns the instrument response function for a given observed photon
 * direction as function of the assumed true photon direction. The result
 * is given by
 *
 * \f[
 *    R(p'|p) = 
 * \f]
 *
 * @todo Write down formula
 * @todo Describe in detail how the response is computed.
 * @todo Implement method.
 ***************************************************************************/
double GCOSResponse::irf(const GEvent&       event,
                         const GPhoton&      photon,
                         const GObservation& obs) const
{
    // Extract COSI observation
    const GCOSObservation* observation = dynamic_cast<const GCOSObservation*>(&obs);
    if (observation == NULL) {
        std::string cls = std::string(typeid(&obs).name());
        std::string msg = "Observation of type \""+cls+"\" is not a COSI "
                          "observation. Please specify a COSI observation "
                          "as argument.";
        throw GException::invalid_argument(G_IRF, msg);
    }
 
    // Extract COSI event bin
    const GCOSEventBin* bin = dynamic_cast<const GCOSEventBin*>(&event);
    if (bin == NULL) {
        std::string cls = std::string(typeid(&event).name());
        std::string msg = "Event of type \""+cls+"\" is not a COSI event. "
                          "Please specify a COSI event as argument.";
        throw GException::invalid_argument(G_IRF, msg);
    }
 
    // Initialise IRF
    double irf = 0.0;
 
    // Get response indices for true quantities (no interpolation so far)
    int itdir = int(src_dirs().dir2pix(photon.dir()).index() + 0.5);
    int iteng = src_ebounds().index(photon.energy());
    int itpol = 0;
 
    // Continue only if indices are valid
    if ((itdir >= 0) && (iteng >= 0)) {
 
        // Get response chunk index
        int index = get_response_chunk(itdir, iteng, itpol);
 
        // Continue only if the response chunk exists
        if (index != -1) {
 
            // Get response indices for measured quantities (no interpolation so far)
            int imdir = int(obs_dirs().dir2pix(bin->dir().dir()).index() + 0.5);
            int imphi = obs_phis().index(bin->dir().phi());
            int imeng = obs_ebounds().index(bin->energy());
            if ((imdir >= 0) && (imphi >= 0) && (imeng >= 0)) {
 
                // Get response element index (no interpolation so far)
                int i = (imphi + imeng * obs_phis().size()) * obs_dirs().npix() + imdir;
 
                // Get response element
                irf = m_rsp_chunks[index][i];
 
                // Multiply with effective area
                irf *= eff_areas()(itdir);
 
                // Multiply with livetime weight
                irf *= (observation->spacecraft().livetime(photon.dir()) /
                        observation->spacecraft().livetime());
 
            } // endif: measured response indices were valid
 
        } // endif: response chunk existed
 
    } // endif: true response indices were valid
 
    // Compile option: Check for NaN/Inf
    #if defined(G_NAN_CHECK)
    if (gammalib::is_notanumber(irf) || gammalib::is_infinite(irf)) {
        std::cout << "*** ERROR: GCOSResponse::irf:";
        std::cout << " NaN/Inf encountered";
        std::cout << " (irf=" << irf;
        std::cout << ")";
        std::cout << std::endl;
    }
    #endif
 
    // Return IRF value
    return irf;
}
 
 
/***********************************************************************//**
 * @brief Return integral of event probability for a given sky model over ROI
 *
 * @param[in] model Sky model.
 * @param[in] obsEng Observed photon energy.
 * @param[in] obsTime Observed photon arrival time.
 * @param[in] obs Observation.
 * @return 0.0
 *
 * @exception GException::feature_not_implemented
 *            Method is not implemented.
 ***************************************************************************/
double GCOSResponse::nroi(const GModelSky&    model,
                          const GEnergy&      obsEng,
                          const GTime&        obsTime,
                          const GObservation& obs) const
{
    // Method is not implemented
    std::string msg = "Spatial integration of sky model over the data space "
                      "is not implemented.";
    throw GException::feature_not_implemented(G_NROI, msg);
 
    // Return Npred
    return (0.0);
}
 
 
/***********************************************************************//**
 * @brief Return true energy boundaries for a specific observed energy
 *
 * @param[in] obsEnergy Observed Energy.
 * @return True energy boundaries for given observed energy.
 *
 * @exception GException::feature_not_implemented
 *            Method is not implemented.
 *
 * @todo Implement this method if you need energy dispersion.
 ***************************************************************************/
GEbounds GCOSResponse::ebounds(const GEnergy& obsEnergy) const
{
    // Initialise an empty boundary object
    GEbounds ebounds;
 
    // Throw an exception
    std::string msg = "Energy dispersion not implemented.";
    throw GException::feature_not_implemented(G_EBOUNDS, msg);
 
    // Return energy boundaries
    return ebounds;
}
 
 
/***********************************************************************//**
 * @brief Get response chunk
 *
 * @param[in] idir True photon direction index [0,...,src_dirs.npix()-1].
 * @param[in] ieng True photon energy index [0,...,src_ebounds.size()-1].
 * @param[in] ipol True polarisation angle index [0,...,src_polarisations.size()-1].
 * @return Response chunk index
 
 * @exception GException::out_of_range
 *            Indices out of range
 *
 * Return index to a response chunk for the specified true photon direction
 * @p idir, true photon energy @p ieng, and true polarisation angle @p ipol.
 * If the corresponding response chunk has not yet been loaded it will be
 * loaded by the method.
 *
 * If no valid response chunk exists for the specified parameters, the method
 * will return -1.
 ***************************************************************************/
int GCOSResponse::get_response_chunk(const int& idir,
                                     const int& ieng,
                                     const int& ipol) const
{
    // Initialise index
    int index = -1;
 
    // Continue only if data layout pointer is valid
    if (m_rsp_chunk_datalayout != NULL) {
 
        // Get chunk dimensions
        int ndirs = m_src_dirs.npix();
        int nengs = m_src_ebounds.size();
        int npols = m_src_pols.size();
 
        // Check validity of indices
        if (idir < 0 || idir >= ndirs) {
            throw GException::out_of_range(G_LOAD_RESPONSE_CHUNK,
                                           "True photon direction index",
                                           idir, ndirs);
        }
        if (ieng < 0 || ieng >= nengs) {
            throw GException::out_of_range(G_LOAD_RESPONSE_CHUNK,
                                           "True photon energy index",
                                           ieng, nengs);
        }
        if ((npols > 0) && (ipol < 0 || ipol >= npols)) {
            throw GException::out_of_range(G_LOAD_RESPONSE_CHUNK,
                                           "True polarisation angle index",
                                           ipol, npols);
        }
 
        // Get response chunk index
        index = get_response_chunk_index(idir, ieng, ipol);
 
        // Continue only if the index was not already loaded
        if (index == -1) {
 
            // Set storage index
            index = flat_index(idir, ieng, ipol);
 
            // Get key for chunk (method checks for key validity)
            const GXmlElement* key = get_response_chunk_key(idir, ieng, ipol);
 
            // Expand environment variables in response file name
            std::string fname = gammalib::expand_env(m_rspfile);
 
            // Open HDF5 file (read-only)
            FILE* fptr = std::fopen(fname.c_str(), "rb");
            if (fptr == NULL) {
                std::string msg = "Unable to open file \""+fname+"\" for "
                                  "read access. Please load a readable "
                                  "HDF5 file before calling that method.";
                throw GException::file_error(G_LOAD_RESPONSE_CHUNK, msg);
            }
 
            // Read data (handles data filtering)
            std::string data = fread_data_chunk(fptr, key,
                                                      m_rsp_chunk_dataspace,
                                                      m_rsp_chunk_datatype,
                                                      m_rsp_chunk_datalayout,
                                                      m_rsp_chunk_datafilter);
 
            // Close file
            std::fclose(fptr);
 
            // Get number of data elements and size
            int length   = data.length();
            int size     = gammalib::toint(m_rsp_chunk_datatype->attribute("size"));
            int elements = length / size;
 
            // Initialise vector to hold data
            std::vector<uint16_t> chunk(elements, 0);
 
            // Initialise data pointer
            const char* ptr = data.c_str();
 
            // Convert data into vector (no checking)
            for (int i = 0; i < elements; ++i, ptr += size) {
                int value = gammalib::hdf5::data_to_int(ptr, m_rsp_chunk_datatype);
                chunk[i]  = value;
            }
 
            // Push back index and data
            m_rsp_chunk_indices.push_back(index);
            m_rsp_chunks.push_back(chunk);
 
        } // endif: index not already loaded
 
    } // endif: data layout pointer was valid
 
    // Return index
    return index;
}
 
 
/***********************************************************************//**
 * @brief Load COSI response
 *
 * @param[in] filename COSI response file.
 *
 * Load the COSI response either a HDF5 file or a FITS file.
 *
 * If a HDF5 file is specified, the method loads the metadata from the file
 * but does not actually load any response data. Response data are loaded
 * upon request from the information stored in the HDF5 metadata.
 *
 * If a FITS file is specified, the actual response chunks are loaded from
 * the FITS file, including all response attributes that are needed by the
 * class. In that case no HDF5 information is provided and loading of new
 * response chunks upon request is not possible.
 ***************************************************************************/
void GCOSResponse::load(const GFilename& filename)
{
    // Clear instance
    clear();
 
    // If file is a FITS file then read response from FITS file
    if (filename.is_fits()) {
 
        // Open FITS file
        GFits fits(filename);
 
        // Read response from FITS file
        read(fits);
 
        // Close FITS file
        fits.close();
 
    } // endif: file was a FITS file
 
    // ... otherwise load the response from a HDF5 file
    else {
 
        // Expand environment variables in filename
        std::string fname = gammalib::expand_env(filename.url());
 
        // Open HDF5 file (read-only)
        FILE* fptr = std::fopen(fname.c_str(), "rb");
        if (fptr == NULL) {
            std::string msg = "Unable to open file \""+fname+"\" for read access. "
                              "Please specify a readable file.";
            throw GException::file_error(G_LOAD, msg);
        }
 
        // Read HDF5 file
        m_hdf5.read(fptr);
 
        // Read response
        read_response(fptr);
 
        // Close file
        std::fclose(fptr);
 
    } // endelse: file was not a FITS file
 
    // Store filename
    m_rspfile = filename;
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Save response
 *
 * @param[in] filename FITS file name.
 * @param[in] clobber Overwrite existing file?
 *
 * Saves the response chunks into a FITS file. If the file exists already and
 * the @p clobber parameter is true, the method will overwrite the content of
 * the existing file. Otherwise, an exception is thrown.
 ***************************************************************************/
void GCOSResponse::save(const GFilename& filename, const bool& clobber) const
{
    // Initialise FITS file
    GFits fits;
 
    // Write response into FITS file
    write(fits);
 
    // Save FITS file
    fits.saveto(filename, clobber);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Read response from FITS file
 *
 * @param[in] fits FITS file.
 *
 * Reads the response from a FITS file. This method overwrites any existing
 * response information yet it does not modify the response file name or
 * any HDF5 metadata. In that was it may be used to load response chunks
 * that were saved on disk without altering the other response information.
 ***************************************************************************/
void GCOSResponse::read(const GFits& fits)
{
    // Read effective area values
    const GFitsTable* eff_area = fits.table("EFFECTIVE AREAS");
    if (!eff_area->is_empty()) {
        const GFitsTableCol* col_eff_area = (*eff_area)["Aeff"];
        int                  size         = eff_area->nrows();
        m_eff_areas                       = GNdarray(size);
        for (int i = 0; i < size; ++i) {
            m_eff_areas(i) = col_eff_area->real(i);
        }
    }
 
    // Read true energy boundaries
    const GFitsTable* src_ebounds = fits.table("TRUE ENERGY BOUNDARIES");
    if (!src_ebounds->is_empty()) {
        m_src_ebounds.read(*src_ebounds);
    }
 
    // Read true polarisation boundaries
    const GFitsTable* src_pol = fits.table("TRUE POLARISATION BOUNDARIES");
    if (!src_pol->is_empty()) {
        m_src_pols.read(*src_pol);
        m_src_pols.unit("deg");
    }
 
    // Read measured Phi boundaries
    const GFitsTable* obs_phi = fits.table("MEASURED PHI BOUNDARIES");
    if (!obs_phi->is_empty()) {
        m_obs_phis.read(*obs_phi);
        m_obs_phis.unit("deg");
    }
 
    // Read measured energy boundaries
    const GFitsTable* obs_ebounds = fits.table("MEASURED ENERGY BOUNDARIES");
    if (!obs_ebounds->is_empty()) {
        m_obs_ebounds.read(*obs_ebounds);
    }
 
    // Read measured sky projection
    const GFitsTable* obs_healpix = fits.table("MEASURED SKY DIRECTIONS");
    m_obs_dirs.read(*obs_healpix);
 
    // Read true sky projection, chunk indices, and response chunks
    const GFitsTable* src_healpix = fits.table("TRUE SKY DIRECTIONS");
    m_src_dirs.read(*src_healpix);
    if (!src_healpix->is_empty()) {
        const GFitsTableCol* col_index   = (*src_healpix)["Index"];
        const GFitsTableCol* col_name    = (*src_healpix)["Name"];
        int                  chunks      = src_healpix->nrows();
        m_rsp_chunk_indices.resize(chunks);
        m_rsp_chunks.resize(chunks);
        for (int i = 0; i < chunks; ++i) {
 
            // Read chunk index
            m_rsp_chunk_indices[i] = col_index->integer(i);
 
            // Read extension name for chunk
            std::string extname = col_name->string(i);
 
            // Get chunk binary table and data column
            const GFitsTable* chunk = fits.table(extname);
            if (!chunk->is_empty()) {
 
                const GFitsTableCol* col_chunk = (*chunk)["Response"];
 
                // Get number of pixels and columns in column
                int npix  = col_chunk->nrows();
                int ncols = col_chunk->number();
 
                // Allocate storage for pixels
                m_rsp_chunks[i].resize(npix*ncols);
 
                // Read response chunk
                for (int col = 0, inx = 0; col < ncols; ++col) {
                    for (int ipix = 0; ipix < npix; ++ipix, ++inx) {
                        m_rsp_chunks[i][inx] = (uint16_t)(col_chunk->integer(ipix, col));
                    }
                }
 
            } // endif: chunk table was not empty
 
        } // endfor: looped over chunks
    } // endif: table was not empty
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Write response into FITS file
 *
 * @param[in] file FITS file.
 *
 * Writes the response into a FITS file. The method appends the following
 * extensions to the FITS file:
 * - EFFECTIVE AREAS
 * - TRUE SKY DIRECTIONS
 * - TRUE ENERGY BOUNDARIES
 * - TRUE POLARISATION BOUNDARIES
 * - MEASURED PHI BOUNDARIES
 * - MEASURED ENERGY BOUNDARIES
 * - CHUNKi
 ***************************************************************************/
void GCOSResponse::write(GFits& fits) const
{
    // Determine number of response chunks
    int chunks = m_rsp_chunk_indices.size();
 
    // Write effective area values into binary table and append table to
    // FITS file
    GFitsBinTable eff_area;
    eff_area.extname("EFFECTIVE AREAS");
    int size = m_eff_areas.size();
    if (size > 0) {
        GFitsTableDoubleCol col_eff_area("Aeff", size);
        //col_eff_area.unit("unit");
        for (int i = 0; i < size; ++i) {
            col_eff_area(i) = m_eff_areas(i);
        }
        eff_area.append(col_eff_area);
    }
    fits.append(eff_area);
 
    // Write true sky projection, chunk indices and FITS extension names
    // into binary table and append table to FITS file
    GFitsBinTable src_healpix;
    m_src_dirs.write(src_healpix);
    src_healpix.extname("TRUE SKY DIRECTIONS");
    if (chunks > 0) {
        GFitsTableLongCol   col_index("Index", chunks);
        GFitsTableStringCol col_name("Name", chunks, 10);
        for (int i = 0; i < chunks; ++i) {
            col_index(i) = m_rsp_chunk_indices[i];
            col_name(i)  = "CHUNK" + gammalib::str(i);
        }
        src_healpix.append(col_index);
        src_healpix.append(col_name);
    }
    fits.append(src_healpix);
 
    // Write true energy boundaries into FITS file
    m_src_ebounds.write(fits, "TRUE ENERGY BOUNDARIES");
 
    // Write true polarisation boundaries into binary table and append
    // table to FITS file
    m_src_pols.write(fits, "TRUE POLARISATION BOUNDARIES");
 
    // Write measured Phi boundaries into binary table and append table
    // to FITS file
    m_obs_phis.write(fits, "MEASURED PHI BOUNDARIES");
 
    // Write measured energy boundaries into FITS file
    m_obs_ebounds.write(fits, "MEASURED ENERGY BOUNDARIES");
 
    // Write measured sky projection. The projection is written in the header
    // of an empty binary table.
    GFitsBinTable obs_healpix;
    m_obs_dirs.write(obs_healpix);
    obs_healpix.extname("MEASURED SKY DIRECTIONS");
    fits.append(obs_healpix);
 
    // Write response chunks into FITS file
    for (int i = 0; i < chunks; ++i) {
        GFitsBinTable chunk;
        m_obs_dirs.write(chunk);
        chunk.extname("CHUNK" + gammalib::str(i));
        int npix     = m_obs_dirs.npix();
        int elements = m_rsp_chunks[i].size();
        int ncols    = elements / npix;
        GFitsTableUShortCol col_chunk("Response", npix, ncols);
        for (int col = 0, inx = 0; col < ncols; ++col) {
            for (int ipix = 0; ipix < npix; ++ipix, ++inx) {
                col_chunk(ipix, col) = m_rsp_chunks[i][inx];
            }
        }
        chunk.append(col_chunk);
        fits.append(chunk);
    }
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Print COSI response information
 *
 * @param[in] chatter Chattiness.
 * @return String containing COSI response information.
 ***************************************************************************/
std::string GCOSResponse::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("=== GCOSResponse ===");
 
        // Append response information
        result.append("\n"+gammalib::parformat("Response file")+m_rspfile.url());
        result.append("\n"+gammalib::parformat("True directions"));
        if (m_src_dirs.nside() == 0) {
            result.append("not defined");
        }
        else {
            result.append(gammalib::str(m_src_dirs.npix()));
            result.append("\n"+gammalib::parformat("True Healpix nside"));
            result.append(gammalib::str(m_src_dirs.nside()));
            result.append("\n"+gammalib::parformat("True Healpix ordering"));
            result.append(m_src_dirs.ordering());
        }
        result.append("\n"+gammalib::parformat("True energies"));
        if (m_src_ebounds.size() == 0) {
            result.append("not defined");
        }
        else {
            result.append(gammalib::str(m_src_ebounds.emin().keV()));
            result.append(" - ");
            result.append(gammalib::str(m_src_ebounds.emax().keV()));
            result.append(" keV");
            result.append("\n"+gammalib::parformat("True energy bins"));
            result.append(gammalib::str(m_src_ebounds.size()));
        }
        result.append("\n"+gammalib::parformat("True polarisation angles"));
        if (m_src_pols.size() == 0) {
            result.append("not defined");
        }
        else {
            result.append(gammalib::str(m_src_pols.min()));
            result.append(" - ");
            result.append(gammalib::str(m_src_pols.max()));
            result.append(" deg");
            result.append("\n"+gammalib::parformat("True polar. angle bins"));
            result.append(gammalib::str(m_src_pols.size()));
        }
        result.append("\n"+gammalib::parformat("Measured directions"));
        if (m_obs_dirs.nside() == 0) {
            result.append("not defined");
        }
        else {
            result.append(gammalib::str(m_obs_dirs.npix()));
            result.append("\n"+gammalib::parformat("Measured Healpix nside"));
            result.append(gammalib::str(m_obs_dirs.nside()));
            result.append("\n"+gammalib::parformat("Measured Healpix ordering"));
            result.append(m_obs_dirs.ordering());
        }
        result.append("\n"+gammalib::parformat("Measured scatter angles"));
        if (m_obs_phis.size() == 0) {
            result.append("not defined");
        }
        else {
            result.append(gammalib::str(m_obs_phis.min()));
            result.append(" - ");
            result.append(gammalib::str(m_obs_phis.max()));
            result.append(" deg");
            result.append("\n"+gammalib::parformat("Measured scatt. angle bins"));
            result.append(gammalib::str(m_obs_phis.size()));
        }
        result.append("\n"+gammalib::parformat("Measured energies"));
        if (m_obs_ebounds.size() == 0) {
            result.append("not defined");
        }
        else {
            result.append(gammalib::str(m_obs_ebounds.emin().keV()));
            result.append(" - ");
            result.append(gammalib::str(m_obs_ebounds.emax().keV()));
            result.append(" keV");
            result.append("\n"+gammalib::parformat("Measured energy bins"));
            result.append(gammalib::str(m_obs_ebounds.size()));
        }
        result.append("\n"+gammalib::parformat("Effective areas"));
        if (m_eff_areas.size() == 0) {
            result.append("not defined");
        }
        else {
            result.append(gammalib::str(min(m_eff_areas)));
            result.append(" - ");
            result.append(gammalib::str(max(m_eff_areas)));
            result.append("\n"+gammalib::parformat("Effective area bins"));
            result.append(gammalib::str(m_eff_areas.size()));
        }
        result.append("\n"+gammalib::parformat("Response chunks loaded"));
        result.append(gammalib::str(m_rsp_chunk_indices.size()));
 
        // If chatter is VERBOSE and HDF5 file metadata are not empty then
        // show the HDF5 file metadata
        if ((chatter == VERBOSE) && (!m_hdf5.is_empty())) {
 
            // Get "DRM" entry and header
            const GXmlElement* drm        = m_hdf5.xml_hdf5_entry("DRM");
            const GXmlElement* drm_header = drm->element("header");
 
            // Get "AXES" entry and table and "AXIS_DESCRIPTIONS" entry and header
            const GXmlElement* axes             = m_hdf5.xml_hdf5_entry("AXES");
            const GXmlElement* axes_table       = axes->element("group")->element("table");
            const GXmlElement* axis_desc        = m_hdf5.xml_hdf5_entry("AXIS_DESCRIPTIONS");
            const GXmlElement* axis_desc_header = axis_desc->element("header");
 
            // Get "COUNTS" entry and header
            const GXmlElement* counts        = m_hdf5.xml_hdf5_entry("COUNTS");
            const GXmlElement* counts_header = counts->element("header");
 
            // Get "DRM" information
            std::string unit = xml_msg_attribute(drm_header, "UNIT", "value");
 
            // Append "DRM" information
            result.append("\n"+gammalib::parformat("Response unit")+unit);
 
            // Get number of axes
            int naxes = axes_table->elements("entry");
 
            // Get number of axis description messages
            int idesc = 0;
            int ndesc = axis_desc_header->elements("message");
 
            // Loop over axes
            for (int i = 0; i < naxes; ++i) {
 
                // Get table entry
                const GXmlElement* entry        = axes_table->element("entry", i);
                const GXmlElement* entry_header = entry->element("header");
 
                // Get header information
                std::string label    = xml_msg_attribute(entry_header, "label", "value");
                std::string unit     = xml_msg_attribute(entry_header, "unit", "value");
                std::string scale    = xml_msg_attribute(entry_header, "scale", "value");
                std::string nside    = xml_msg_attribute(entry_header, "nside", "value");
                std::string scheme   = xml_msg_attribute(entry_header, "scheme", "value");
                std::string coordsys = xml_msg_attribute(entry_header, "coordsys", "value");
 
                // Get corresponding axis description message
                std::string desc = xml_msg_attribute(axis_desc_header, label, "value");
 
                // Get corresponding axis dimension
                const GXmlElement* dataspace = xml_msg_dataspace(counts_header);
                std::string size = dataspace->attribute("size"+gammalib::str(i));
 
                // Append available axis information
                result.append("\n"+gammalib::parformat("Axis"+entry->attribute("name"))+label);
                result.append("\n"+gammalib::parformat(" Size")+size);
                if (desc.length() > 0) {
                    result.append("\n"+gammalib::parformat(" Description")+desc);
                }
                result.append("\n"+gammalib::parformat(" Size")+size);
                if (unit.length() > 0) {
                    result.append("\n"+gammalib::parformat(" Units")+unit);
                }
                if (scale.length() > 0) {
                    result.append("\n"+gammalib::parformat(" Scale")+scale);
                }
                if (nside.length() > 0) {
                    result.append("\n"+gammalib::parformat(" Nside")+nside);
                }
                if (scheme.length() > 0) {
                    result.append("\n"+gammalib::parformat(" Scheme")+scheme);
                }
                if (coordsys.length() > 0) {
                    result.append("\n"+gammalib::parformat(" Coordsys")+coordsys);
                }
 
            } // endfor: looped over axes
 
        } // endif: response file was not empty
 
    } // endif: chatter was not silent
 
    // Return result
    return result;
}
 
 
/*==========================================================================
 =                                                                         =
 =                             Private methods                             =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Initialise class members
 ***************************************************************************/
void GCOSResponse::init_members(void)
{
    // Initialise members
    m_rspfile.clear();
    m_hdf5.clear();
    m_src_dirs.clear();
    m_src_ebounds.clear();
    m_src_pols.clear();
    m_obs_dirs.clear();
    m_obs_phis.clear();
    m_obs_ebounds.clear();
    m_eff_areas.clear();
 
    // Initialise members for response chunk handling
    m_rsp_chunk_dataspace  = NULL;
    m_rsp_chunk_datatype   = NULL;
    m_rsp_chunk_datalayout = NULL;
    m_rsp_chunk_datafilter = NULL;
    m_rsp_chunk_indices.clear();
    m_rsp_chunks.clear();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Copy class members
 *
 * @param[in] rsp COSI response function.
 ***************************************************************************/
void GCOSResponse::copy_members(const GCOSResponse& rsp)
{
    // Copy members
    m_rspfile           = rsp.m_rspfile;
    m_hdf5              = rsp.m_hdf5;
    m_src_dirs          = rsp.m_src_dirs;
    m_src_ebounds       = rsp.m_src_ebounds;
    m_src_pols          = rsp.m_src_pols;
    m_obs_dirs          = rsp.m_obs_dirs;
    m_obs_phis          = rsp.m_obs_phis;
    m_obs_ebounds       = rsp.m_obs_ebounds;
    m_eff_areas         = rsp.m_eff_areas;
 
    // Copy members for response chunk handling
    m_rsp_chunk_indices = rsp.m_rsp_chunk_indices;
    m_rsp_chunks        = rsp.m_rsp_chunks;
 
    // Set pointers for response chunk handling
    set_response_chunk_pointers();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Delete class members
 ***************************************************************************/
void GCOSResponse::free_members(void)
{
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Read response from HDF5 file
 *
 * @param[in] fptr File pointer to HDF5 file.
 *
 * Read response data and information using the metadata that is provided
 * by the m_hdf5 member.
 *
 * This method does nothing if the m_hdf5 member is empty (which means that
 * no HDF5 metadata are available).
 ***************************************************************************/
void GCOSResponse::read_response(FILE* fptr)
{
    // Continue only if HDF5 metadata are not empty
    if (!m_hdf5.is_empty()) {
 
        // Initialise XML element pointers
        const GXmlElement* dataspace  = NULL;
        const GXmlElement* datatype   = NULL;
        const GXmlElement* datalayout = NULL;
        const GXmlElement* datafilter = NULL;
 
        // Get "AXES" entry and table
        const GXmlElement* axes       = m_hdf5.xml_hdf5_entry("AXES");
        const GXmlElement* axes_table = axes->element("group")->element("table");
 
        // Get number of axes
        int naxes = axes_table->elements("entry");
 
        // Loop over axes
        for (int i = 0; i < naxes; ++i) {
 
            // Get table entry
            const GXmlElement* entry        = axes_table->element("entry", i);
            const GXmlElement* entry_header = entry->element("header");
 
            // Get dataspace, datatype and datalayout
            dataspace  = xml_msg_dataspace(entry_header);
            datatype   = xml_msg_datatype(entry_header);
            datalayout = xml_msg_datalayout(entry_header);
 
            // Get header information
            std::string label    = xml_msg_attribute(entry_header, "label", "value");
            std::string unit     = xml_msg_attribute(entry_header, "unit", "value");
            std::string scale    = xml_msg_attribute(entry_header, "scale", "value");
            std::string nside    = xml_msg_attribute(entry_header, "nside", "value");
            std::string scheme   = xml_msg_attribute(entry_header, "scheme", "value");
            std::string coordsys = xml_msg_attribute(entry_header, "coordsys", "value");
 
            // Handle axes
            if (label == "Ei") {  // True energies
                m_src_ebounds = read_ebounds(fptr, dataspace, datatype, datalayout, unit);
            }
            else if (label == "Em") {  // Measured energies
                m_obs_ebounds = read_ebounds(fptr, dataspace, datatype, datalayout, unit);
            }
            else if (label == "Phi") {  // Measured energies
                m_obs_phis = read_bounds(fptr, dataspace, datatype, datalayout);
                m_obs_phis.unit("deg");
            }
            else if (label == "NuLambda") {  // True sky directions
                m_src_dirs = GHealpix(gammalib::toint(nside), scheme, "CEL");
            }
            else if (label == "PsiChi") {  // True sky directions
                m_obs_dirs = GHealpix(gammalib::toint(nside), scheme, "CEL");
            }
            else if (label == "Pol") {  // True polarisation angles
                m_src_pols = read_bounds(fptr, dataspace, datatype, datalayout);
                m_src_pols.unit("deg");
            }
 
        } // endif: looped over axes
 
        // Get "EFF_AREA" entry and header
        const GXmlElement* eff_area        = m_hdf5.xml_hdf5_entry("EFF_AREA");
        const GXmlElement* eff_area_header = eff_area->element("header");
 
        // Get dataspace, datatype and datalayout
        dataspace  = xml_msg_dataspace(eff_area_header);
        datatype   = xml_msg_datatype(eff_area_header);
        datalayout = xml_msg_datalayout(eff_area_header);
 
        // Read effective area vector
        m_eff_areas = read_ndarray(fptr, dataspace, datatype, datalayout);
 
    } // endif: HDF5 metadata were not empty
 
    // Set pointers for response chunk handling
    set_response_chunk_pointers();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Read boundaries from HDF5 file
 *
 * @param[in] fptr File pointer to HDF5 file.
 * @param[in] dataspace Pointer to "Dataspace" XML element.
 * @param[in] datatype Pointer to "Datatype" XML element.
 * @param[in] datalayout Pointer to "DataStorageLayout" XML element.
 * @return Boundaries.
 *
 * @exception GException::invalid_argument
 *            Invalid pointer to XML element
 *
 * Reads boundaries from HDF5 file and returns them as a GBounds object.
 *
 * The method expects that the @p dataspace, @p datatype and @p datalayout
 * pointers are not NULL. An exception is thrown if a NULL pointer is
 * encountered.
 ***************************************************************************/
GBounds GCOSResponse::read_bounds(FILE*              fptr,
                                  const GXmlElement* dataspace,
                                  const GXmlElement* datatype,
                                  const GXmlElement* datalayout)
{
    // Initialise boundaries
    GBounds bounds;
 
    // Check pointers
    if (dataspace == NULL) {
        std::string msg = "Invalid dataspace pointer specified. No dataspace "
                          "information was found in the COSI response file. "
                          "Please specify a valid COSI response file.";
        throw GException::invalid_argument(G_READ_BOUNDS, msg);
    }
    if (datatype == NULL) {
        std::string msg = "Invalid datatype pointer specified. No datatype "
                          "information was found in the COSI response file. "
                          "Please specify a valid COSI response file.";
        throw GException::invalid_argument(G_READ_BOUNDS, msg);
    }
    if (datalayout == NULL) {
        std::string msg = "Invalid datalayout pointer specified. No datalayout "
                          "information was found in the COSI response file. "
                          "Please specify a valid COSI response file.";
        throw GException::invalid_argument(G_READ_BOUNDS, msg);
    }
 
    // Read data
    std::string data = fread_data(fptr, dataspace, datatype, datalayout);
 
    // Get number of data elements and their size
    int elements = gammalib::toint(dataspace->attribute("elements"));
    int size     = gammalib::toint(datatype->attribute("size"));
 
    // Loop over data elements
    double min;
    const char* ptr = data.c_str();
    for (int i = 0; i < elements; ++i, ptr += size) {
        double max = data_to_double(ptr, datatype);
        if (i > 0) {
            bounds.append(min, max);
        }
        min = max;
    }
 
    // Return boundaries
    return bounds;
}
 
 
/***********************************************************************//**
 * @brief Read energy boundaries from HDF5 file
 *
 * @param[in] fptr File pointer to HDF5 file.
 * @param[in] dataspace Pointer to "Dataspace" XML element.
 * @param[in] datatype Pointer to "Datatype" XML element.
 * @param[in] datalayout Pointer to "DataStorageLayout" XML element.
 * @param[in] unit Energy unit string.
 * @return Energy boundaries
 *
 * @exception GException::invalid_argument
 *            Invalid pointer to XML element
 *
 * Reads energy boundaries from HDF5 file. If the energy @p unit string is
 * empty the energies are assumed to be given in keV.
 *
 * The method expects that the @p dataspace, @p datatype and @p datalayout
 * pointers are not NULL. An exception is thrown if a NULL pointer is
 * encountered.
 ***************************************************************************/
GEbounds GCOSResponse::read_ebounds(FILE*              fptr,
                                    const GXmlElement* dataspace,
                                    const GXmlElement* datatype,
                                    const GXmlElement* datalayout,
                                    const std::string& unit)
{
    // Initialise energy boundaries
    GEbounds ebounds;
 
    // Check pointers
    if (dataspace == NULL) {
        std::string msg = "Invalid dataspace pointer specified. No dataspace "
                          "information was found in the COSI response file. "
                          "Please specify a valid COSI response file.";
        throw GException::invalid_argument(G_READ_EBOUNDS, msg);
    }
    if (datatype == NULL) {
        std::string msg = "Invalid datatype pointer specified. No datatype "
                          "information was found in the COSI response file. "
                          "Please specify a valid COSI response file.";
        throw GException::invalid_argument(G_READ_EBOUNDS, msg);
    }
    if (datalayout == NULL) {
        std::string msg = "Invalid datalayout pointer specified. No datalayout "
                          "information was found in the COSI response file. "
                          "Please specify a valid COSI response file.";
        throw GException::invalid_argument(G_READ_EBOUNDS, msg);
    }
 
    // Read data
    std::string data = fread_data(fptr, dataspace, datatype, datalayout);
 
    // Get number of data elements and their size
    int elements = gammalib::toint(dataspace->attribute("elements"));
    int size     = gammalib::toint(datatype->attribute("size"));
 
    // Set unit string
    std::string str_unit = (unit.empty()) ? "keV" : unit;
 
    // Loop over data elements
    GEnergy     emin;
    const char* ptr = data.c_str();
    for (int i = 0; i < elements; ++i, ptr += size) {
        double  value = data_to_double(ptr, datatype);
        GEnergy emax(value, unit);
        if (i > 0) {
            ebounds.append(emin, emax);
        }
        emin = emax;
    }
 
    // Return energy boundaries
    return ebounds;
}
 
 
/***********************************************************************//**
 * @brief Read 1-dimensional array from HDF5 file
 *
 * @param[in] fptr File pointer to HDF5 file.
 * @param[in] dataspace Pointer to "Dataspace" XML element.
 * @param[in] datatype Pointer to "Datatype" XML element.
 * @param[in] datalayout Pointer to "DataStorageLayout" XML element.
 * @return Array.
 *
 * @exception GException::invalid_argument
 *            Invalid pointer to XML element
 *
 * Reads 1-dimensional array from HDF5 file and return it as GNdarray object.
 *
 * The method expects that the @p dataspace, @p datatype and @p datalayout
 * pointers are not NULL. An exception is thrown if a NULL pointer is
 * encountered.
 ***************************************************************************/
GNdarray GCOSResponse::read_ndarray(FILE*              fptr,
                                    const GXmlElement* dataspace,
                                    const GXmlElement* datatype,
                                    const GXmlElement* datalayout)
{
    // Initialise array
    GNdarray array;
 
    // Check pointers
    if (dataspace == NULL) {
        std::string msg = "Invalid dataspace pointer specified. No dataspace "
                          "information was found in the COSI response file. "
                          "Please specify a valid COSI response file.";
        throw GException::invalid_argument(G_READ_NDARRAY, msg);
    }
    if (datatype == NULL) {
        std::string msg = "Invalid datatype pointer specified. No datatype "
                          "information was found in the COSI response file. "
                          "Please specify a valid COSI response file.";
        throw GException::invalid_argument(G_READ_NDARRAY, msg);
    }
    if (datalayout == NULL) {
        std::string msg = "Invalid datalayout pointer specified. No datalayout "
                          "information was found in the COSI response file. "
                          "Please specify a valid COSI response file.";
        throw GException::invalid_argument(G_READ_NDARRAY, msg);
    }
 
    // Read data
    std::string data = fread_data(fptr, dataspace, datatype, datalayout);
 
    // Get number of data elements and their size
    int elements = gammalib::toint(dataspace->attribute("elements"));
    int size     = gammalib::toint(datatype->attribute("size"));
 
    // Allocate array
    array = GNdarray(elements);
 
    // Store effective area vector
    const char* ptr = data.c_str();
    for (int i = 0; i < elements; ++i, ptr += size) {
        array(i) = data_to_double(ptr, datatype);
    }
 
    // Return array
    return array;
}
 
 
/***********************************************************************//**
 * @brief Read response matrix from HDF5 file
 *
 * @param[in] fptr File pointer to HDF5 file.
 * @param[in] dataspace Pointer to "Dataspace" XML element.
 * @param[in] datatype Pointer to "Datatype" XML element.
 * @param[in] datalayout Pointer to "DataStorageLayout" XML element.
 * @param[in] datafilter Pointer to "DataStorageFilterPipeline" XML element.
 *
 * @exception GException::invalid_argument
 *            Invalid pointer to XML element
 *
 * @todo To be implemented (if we really need that method)
 ***************************************************************************/
void GCOSResponse::read_response_matrix(FILE*              fptr,
                                        const GXmlElement* dataspace,
                                        const GXmlElement* datatype,
                                        const GXmlElement* datalayout,
                                        const GXmlElement* datafilter)
{
    // Check pointers
    if (dataspace == NULL) {
        std::string msg = "Invalid dataspace pointer provided for response "
                          "matrix. No dataspace information was provided "
                          "in the COSI response file. Please specify a "
                          "valid COSI response file.";
        throw GException::invalid_argument(G_READ_RESPONSE_MATRIX, msg);
    }
    if (datatype == NULL) {
        std::string msg = "Invalid datatype pointer provided for response "
                          "matrix. No datatype information was provided "
                          "in the COSI response file. Please specify a "
                          "valid COSI response file.";
        throw GException::invalid_argument(G_READ_RESPONSE_MATRIX, msg);
    }
    if (datalayout == NULL) {
        std::string msg = "Invalid datalayout pointer provided for response "
                          "matrix. No datalayout information was provided "
                          "in the COSI response file. Please specify a "
                          "valid COSI response file.";
        throw GException::invalid_argument(G_READ_RESPONSE_MATRIX, msg);
    }
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Set XML element pointers for response chunk handling
 *
 * Sets pointers to the "Dataspace", "Datatype", "DataStorageLayout" and
 * "DataStorageFilterPipeline" XML elements. If no HDF5 metadata are
 * available all pointers are set to NULL.
 ***************************************************************************/
void GCOSResponse::set_response_chunk_pointers(void)
{
    // Initialise pointers
    m_rsp_chunk_dataspace  = NULL;
    m_rsp_chunk_datatype   = NULL;
    m_rsp_chunk_datalayout = NULL;
    m_rsp_chunk_datafilter = NULL;
 
    // Continue only if HDF5 metadata are not empty
    if (!m_hdf5.is_empty()) {
 
        // Get "COUNTS" entry and header
        const GXmlElement* counts        = m_hdf5.xml_hdf5_entry("COUNTS");
        const GXmlElement* counts_header = counts->element("header");
 
        // Get dataspace, datatype, datalayout and datafilter
        m_rsp_chunk_dataspace  = xml_msg_dataspace(counts_header);
        m_rsp_chunk_datatype   = xml_msg_datatype(counts_header);
        m_rsp_chunk_datalayout = xml_msg_datalayout(counts_header);
        m_rsp_chunk_datafilter = xml_msg_datafilter(counts_header);
 
    } // endif: HDF5 metadata were not empty
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Get pointer to key XML element for specified chunk index
 *
 * @param[in] idir True photon direction index (0,...,src_dirs.npix()-1).
 * @param[in] ieng True photon energy index (0,...,src_ebounds.size()-1).
 * @param[in] ipol True polarisation angle index (0,...,src_polarisations.size()-2).
 * @return Pointer to key XML element.
 *
 * @exception GException::runtime_error
 *            Key XML element not found.
 *            Key XML element with wrong index found.
 *
 * Returns key XML element within the data layout B-tree that has an "index0"
 * attribute that corresponds to the specified @p index. The method throws
 * exceptions if no key element was found or if the key element has the
 * wrong index.
 *
 * This method assumes that the data layout pointer m_rsp_chunk_datalayout is
 * valid.
 ***************************************************************************/
const GXmlElement* GCOSResponse::get_response_chunk_key(const int& idir,
                                                        const int& ieng,
                                                        const int& ipol) const
{
    // Initialise pointer to ket element
    const GXmlElement* key = NULL;
 
    // Get pointer to B-tree
    const GXmlElement* btree = m_rsp_chunk_datalayout->element("btree");
 
    // Convert 3D indices into flat 1D index
    int index = flat_index(idir, ieng, ipol);
 
    // Search B-tree structure
    while (true) {
 
        // Get number of keys in B-tree
        int nkeys = btree->elements("key");
 
        // If there is more than one key then start bi-section search
        if (nkeys > 0) {
 
            // Initialise key range
            int ikey_min = 0;
            int ikey_max = nkeys;
 
            // Set initial key as mid-point
            int ikey = nkeys / 2;
 
            // Do bi-section search
            while (ikey_min != ikey_max) {
 
                // Get key
                key = btree->element("key", ikey);
 
                // Get key index
                int idir_key  = gammalib::toint(key->attribute("index0"));
                int ieng_key  = gammalib::toint(key->attribute("index1"));
                int ipol_key  = (m_src_pols.size() == 0) ?
                                0 : gammalib::toint(key->attribute("index2"));
                int index_key = flat_index(idir_key, ieng_key, ipol_key);
 
                // If index_key is larger than index the chunk resides before
                if (index_key > index) {
                    ikey_max = ikey;
                    ikey     = ikey_min + (ikey_max-ikey_min) / 2;
                    if (ikey == ikey_max) {
                        break;
                    }
                }
 
                // If index_key is smaller than index the chunk resides after
                else if (index_key < index) {
                    ikey_min = ikey;
                    ikey     = ikey_min + (ikey_max-ikey_min) / 2;
                    if (ikey == ikey_min) {
                        break;
                    }
                }
 
                // ... otherwise the index was found and we break
                else {
                    break;
                }
 
                // If interval is zero then key was found and we break
                if (ikey_min == ikey_max) {
                    break;
                }
 
            } // endwhile: did bi-section search
 
            // If key does not contain another B-tree we are done
            if (key->elements("btree") == 0) {
                break;
            }
 
            // ... otherwise get the B-tree and continue
            else {
                btree = key->element("btree");
            }
 
        } // endif: there were keys in element
 
    } // endwhile: searched B-tree
 
    // Check that key was found
    if (key == NULL) {
        std::string msg = "Response chunk with index "+gammalib::str(index)+
                          " not found. This should never happen if the HDF5 "
                          "file is correct.";
        throw GException::runtime_error(G_GET_RESPONSE_CHUNK_KEY, msg);
    }
 
    // Check that correct chunk was found
    int idir_key  = gammalib::toint(key->attribute("index0"));
    int ieng_key  = gammalib::toint(key->attribute("index1"));
    int ipol_key  = (m_src_pols.size() == 0) ?
                    0 : gammalib::toint(key->attribute("index2"));
    int index_key = flat_index(idir_key, ieng_key, ipol_key);
    if (index_key != index) {
        std::string msg = "Bad response chunk with index " +
                          gammalib::str(index_key) +
                          " found while index " +
                          gammalib::str(index)+" was requested. This "
                          "should never happen, the code is not doing "
                          "what it is expected to do.";
        throw GException::runtime_error(G_GET_RESPONSE_CHUNK_KEY, msg);
    }
 
    // Return key
    return key;
}
 
 
/***********************************************************************//**
 * @brief Get response chunk vector index for specified chunk index
 *
 * @param[in] idir True photon direction index (0,...,src_dirs.npix()-1).
 * @param[in] ieng True photon energy index (0,...,src_ebounds.size()-1).
 * @param[in] ipol True polarisation angle index (0,...,src_polarisations.size()-2).
 * @return Response chunk vector index (-1 if chunk index was not found).
 *
 * Returns the response vector index for a specified chunk @p index. If the
 * chunk index does not exist in the response vectors a value of -1 is
 * returned.
 ***************************************************************************/
int GCOSResponse::get_response_chunk_index(const int& idir,
                                           const int& ieng,
                                           const int& ipol) const
{
    // Initialise index with invalid value
    int inx = -1;
 
    // Get number of response chunk indices
    int n = m_rsp_chunk_indices.size();
 
    // Continue only if there are chunk indices
    if (n > 0) {
 
        // Determine storage index
        int index = flat_index(idir, ieng, ipol);
 
        // Search specified chunk index
        for (int i = 0; i < n; ++i) {
            if (m_rsp_chunk_indices[i] == index) {
                inx = i;
                break;
            }
        }
 
    } // endif: there were chunk indices
 
    // Return index
    return inx;
}
 
 
/***********************************************************************//**
 * @brief Convert 3D index into flat 1D index
 *
 * @param[in] idir True photon direction index (0,...,src_dirs.npix()-1).
 * @param[in] ieng True photon energy index (0,...,src_ebounds.size()-1).
 * @param[in] ipol True polarisation angle index (0,...,src_polarisations.size()-1).
 * @return Flat index.
 *
 * Returns the flattened index of a specified true photon direction index
 * @p idir, true photon energy index @p ieng, and true polarisation angle
 * index @p ipol. If no source polarisations exists, the value of @p ipol
 * is ignored.
 *
 * The method does not check the validity of the input indices, this needs
 * to be done by the client.
 ***************************************************************************/
int GCOSResponse::flat_index(const int& idir,
                             const int& ieng,
                             const int& ipol) const
{
    // Determine number of energy boundaries and polarisation boundaries
    int nengs = m_src_ebounds.size();
    int npols = m_src_pols.size();
 
    // Determine flat index
    int index = (npols == 0) ? idir * nengs + ieng
                             : (idir * nengs + ieng) * npols + ipol;
 
    // Return flat index
    return index;
}