GCOMResponse.cpp

Go to the documentation of this file.

/***************************************************************************
 *                GCOMResponse.cpp - COMPTEL Response class                *
 * ----------------------------------------------------------------------- *
 *  copyright (C) 2012-2024 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 GCOMResponse.cpp
 * @brief COMPTEL response class implementation
 * @author Juergen Knoedlseder
 */
 
/* __ Includes ___________________________________________________________ */
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <string>
#include <typeinfo>
#include "GTools.hpp"
#include "GMath.hpp"
#include "GVector.hpp"
#include "GMatrix.hpp"
#include "GIntegrals.hpp"
#include "GFits.hpp"
#include "GCaldb.hpp"
#include "GEvent.hpp"
#include "GPhoton.hpp"
#include "GSource.hpp"
#include "GEnergy.hpp"
#include "GTime.hpp"
#include "GObservation.hpp"
#include "GFitsImage.hpp"
#include "GFitsImageFloat.hpp"
#include "GFilename.hpp"
#include "GModelSky.hpp"
#include "GModelSpatialPointSource.hpp"
#include "GModelSpatialRadial.hpp"
#include "GModelSpatialElliptical.hpp"
#include "GModelSpatialDiffuse.hpp"
#include "GCOMResponse.hpp"
#include "GCOMObservation.hpp"
#include "GCOMEventCube.hpp"
#include "GCOMEventBin.hpp"
#include "GCOMInstDir.hpp"
#include "com_helpers_response_vector.hpp"
 
 
/* __ Method name definitions ____________________________________________ */
#define G_IRF           "GCOMResponse::irf(GEvent&, GPhoton&, GObservation&)"
#define G_IRF_SPATIAL         "GCOMResponse::irf_spatial(GEvent&, GSource&, "\
                                                             "GObservation&)"
#define G_NROI            "GCOMResponse::nroi(GModelSky&, GEnergy&, GTime&, "\
                                                             "GObservation&)"
#define G_EBOUNDS                           "GCOMResponse::ebounds(GEnergy&)"
#define G_BACKPROJECT   "GCOMResponse::backproject(GObservation&, GEvents*, "\
                                                                  "GSkyMap*)"
#define G_IRF_PTSRC     "GCOMResponse::irf_ptsrc(GModelSky&, GObservation&, "\
                                                                  "GMatrix*)"
#define G_IRF_RADIAL   "GCOMResponse::irf_radial(GModelSky&, GObservation&, "\
                                                                  "GMatrix*)"
#define G_IRF_ELLIPTICAL          "GCOMResponse::irf_elliptical(GModelSky&, "\
                                                  " GObservation&, GMatrix*)"
#define G_IRF_DIFFUSE "GCOMResponse::irf_diffuse(GModelSky&, GObservation&, "\
                                                                  "GMatrix*)"
 
/* __ Macros _____________________________________________________________ */
 
/* __ Coding definitions _________________________________________________ */
#define G_IRF_RADIAL_METHOD     2 //!< 0=integration, 1=direct, 2=mix
#define G_IRF_ELLIPTICAL_METHOD 2 //!< 0=integration, 1=direct, 2=mix
#define G_IRF_DIFFUSE_DIRECT      //!< Use direct computation for diffuse response
 
/* __ Debug definitions __________________________________________________ */
//#define G_DEBUG_COMPUTE_FAQ    //!< Debug FAQ computation
//#define G_DEBUG_IRF_TIMING     //!< Time response computation
 
/* __ Constants __________________________________________________________ */
 
 
/*==========================================================================
 =                                                                         =
 =                       Constructors/destructors                          =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Void constructor
 *
 * Creates an empty COMPTEL response.
 ***************************************************************************/
GCOMResponse::GCOMResponse(void) : GResponse()
{
    // Initialise members
    init_members();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Copy constructor
 *
 * @param[in] rsp COM response.
 **************************************************************************/
GCOMResponse::GCOMResponse(const GCOMResponse& rsp) : GResponse(rsp)
{
    // Initialise members
    init_members();
 
    // Copy members
    copy_members(rsp);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Response constructor
 *
 * @param[in] caldb Calibration database.
 * @param[in] iaqname IAQ file name.
 *
 * Create COMPTEL response by loading an IAQ file from a calibration
 * database.
 ***************************************************************************/
GCOMResponse::GCOMResponse(const GCaldb&      caldb,
                           const std::string& iaqname) : GResponse()
{
    // Initialise members
    init_members();
 
    // Set calibration database
    this->caldb(caldb);
 
    // Load IRF
    this->load(iaqname);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Destructor
 *
 * Destroys instance of COMPTEL response object.
 ***************************************************************************/
GCOMResponse::~GCOMResponse(void)
{
    // Free members
    free_members();
 
    // Return
    return;
}
 
 
/*==========================================================================
 =                                                                         =
 =                                Operators                                =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Assignment operator
 *
 * @param[in] rsp COMPTEL response.
 * @return COMPTEL response.
 *
 * Assigns COMPTEL response object to another COMPTEL response object. The
 * assignment performs a deep copy of all information, hence the original
 * object from which the assignment has been performed can be destroyed after
 * this operation without any loss of information.
 ***************************************************************************/
GCOMResponse& GCOMResponse::operator=(const GCOMResponse& 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 COMPTEL response object by resetting all members to an initial
 * state. Any information that was present in the object before will be lost.
 ***************************************************************************/
void GCOMResponse::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 COMPTEL response.
 ***************************************************************************/
GCOMResponse* GCOMResponse::clone(void) const
{
    return new GCOMResponse(*this);
}
 
 
/***********************************************************************//**
 * @brief Return value of instrument response function
 *
 * @param[in] event Observed event.
 * @param[in] photon Incident photon.
 * @param[in] obs Observation.
 * @return Instrument response function (\f$cm^2 sr^{-1}\f$)
 *
 * @exception GException::invalid_argument
 *            Observation is not a COMPTEL observation.
 *            Event is not a COMPTEL event bin.
 * @exception GException::invalid_value
 *            Response not initialised with a valid IAQ
 *
 * 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[
 *    {\tt IRF} = \frac{{\tt IAQ} \times {\tt DRG} \times {\tt DRX}}
 *                     {T \times {\tt TOFCOR}} \times {\tt PHASECOR}
 * \f]
 *
 * where
 * - \f${\tt IRF}\f$ is the instrument response function (\f$cm^2 sr^{-1}\f$),
 * - \f${\tt IAQ}\f$ is the COMPTEL response matrix (\f$sr^{-1}\f$),
 * - \f${\tt DRG}\f$ is the geometry factor (probability),
 * - \f${\tt DRX}\f$ is the exposure (\f$cm^2 s\f$),
 * - \f$T\f$ is the ontime (\f$s\f$), and
 * - \f${\tt TOFCOR}\f$ is a correction that accounts for the Time of Flight
 *   selection window.
 * - \f${\tt PHASECOR}\f$ is a correction that accounts for pulsar phase
 *   selection.
 *
 * The observed photon direction is spanned by the three values \f$\Chi\f$,
 * \f$\Psi\f$, and \f$\bar{\varphi})\f$. \f$\Chi\f$ and \f$\Psi\f$ is the
 * scatter direction of the event, given in sky coordinates.
 * \f$\bar{\varphi}\f$ is the Compton scatter angle, computed from the
 * energy deposits in the two detector planes.
 ***************************************************************************/
double GCOMResponse::irf(const GEvent&       event,
                         const GPhoton&      photon,
                         const GObservation& obs) const
{
    // Extract COMPTEL observation
    const GCOMObservation* observation = dynamic_cast<const GCOMObservation*>(&obs);
    if (observation == NULL) {
        std::string cls = std::string(typeid(&obs).name());
        std::string msg = "Observation of type \""+cls+"\" is not a COMPTEL "
                          "observations. Please specify a COMPTEL observation "
                          "as argument.";
        throw GException::invalid_argument(G_IRF, msg);
    }
 
    // Extract COMPTEL event cube
    const GCOMEventCube* cube = dynamic_cast<const GCOMEventCube*>(observation->events());
    if (cube == NULL) {
        std::string cls = std::string(typeid(&cube).name());
        std::string msg = "Event cube of type \""+cls+"\" is  not a COMPTEL "
                          "event cube. Please specify a COMPTEL event cube "
                          "as argument.";
        throw GException::invalid_argument(G_IRF, msg);
    }
 
    // Extract COMPTEL event bin
    const GCOMEventBin* bin = dynamic_cast<const GCOMEventBin*>(&event);
    if (bin == NULL) {
        std::string cls = std::string(typeid(&event).name());
        std::string msg = "Event of type \""+cls+"\" is  not a COMPTEL event. "
                          "Please specify a COMPTEL event as argument.";
        throw GException::invalid_argument(G_IRF, msg);
    }
 
    // Throw an exception if COMPTEL response is not set or if
    if (m_iaq.empty()) {
        std::string msg = "COMPTEL response is empty. Please initialise the "
                          "response with an \"IAQ\".";
        throw GException::invalid_value(G_IRF, msg);
    }
    else if (m_phigeo_bin_size == 0.0) {
        std::string msg = "COMPTEL response has a zero Phigeo bin size. "
                          "Please initialise the response with a valid "
                          "\"IAQ\".";
        throw GException::invalid_value(G_IRF, msg);
    }
 
    // Extract event parameters
    const GCOMInstDir& obsDir = bin->dir();
 
    // Extract photon parameters
    const GSkyDir& srcDir  = photon.dir();
    const GTime&   srcTime = photon.time();
 
    // Compute angle between true photon arrival direction and scatter
    // direction (Chi,Psi)
    double phigeo = srcDir.dist_deg(obsDir.dir());
 
    // Compute scatter angle index
    int iphibar = int(obsDir.phibar() / m_phibar_bin_size);
 
    // Initialise IRF
    double iaq = 0.0;
    double irf = 0.0;
 
    // Extract IAQ value by linear inter/extrapolation in Phigeo
    if (iphibar < m_phibar_bins) {
        double phirat  = phigeo / m_phigeo_bin_size; // 0.5 at bin centre
        int    iphigeo = int(phirat);                // index into which Phigeo falls
        double eps     = phirat - iphigeo - 0.5;     // 0.0 at bin centre [-0.5, 0.5[
        if (iphigeo < m_phigeo_bins) {
            int i = iphibar * m_phigeo_bins + iphigeo;
            if (eps < 0.0) { // interpolate towards left
                if (iphigeo > 0) {
                    iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
                }
                else {
                    iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
                }
            }
            else {           // interpolate towards right
                if (iphigeo < m_phigeo_bins-1) {
                    iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
                }
                else {
                    iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
                }
            }
        }
    }
 
    // Continue only if IAQ is positive
    if (iaq > 0.0) {
 
        // Get DRG value (unit: probability)
        double drg = observation->drg().map()(obsDir.dir(), iphibar);
 
        // Get DRX value (unit: cm^2 sec)
        double drx = observation->drx().map()(srcDir);
 
        // Get ontime (unit: s)
        double ontime = observation->ontime();
 
        // Get ToF correction
        double tofcor = cube->dre().tof_correction();
 
        // Get pulsar phase correction
        double phasecor = cube->dre().phase_correction();
 
        // Compute IRF value
        irf = iaq * drg * drx / (ontime * tofcor) * phasecor;
 
        // Apply deadtime correction
        irf *= obs.deadc(srcTime);
 
        // Make sure that IRF is positive
        if (irf < 0.0) {
            irf = 0.0;
        }
 
        // Compile option: Check for NaN/Inf
        #if defined(G_NAN_CHECK)
        if (gammalib::is_notanumber(irf) || gammalib::is_infinite(irf)) {
            std::cout << "*** ERROR: GCOMResponse::irf:";
            std::cout << " NaN/Inf encountered";
            std::cout << " (irf=" << irf;
            std::cout << ", iaq=" << iaq;
            std::cout << ", drg=" << drg;
            std::cout << ", drx=" << drx;
            std::cout << ")";
            std::cout << std::endl;
        }
        #endif
 
    } // endif: IAQ was positive
 
    // 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 GCOMResponse::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.
 ***************************************************************************/
GEbounds GCOMResponse::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 Load COMPTEL response.
 *
 * @param[in] rspname COMPTEL response name.
 *
 * Loads the COMPTEL response with specified name @p rspname.
 *
 * The method first attempts to interpret @p rspname as a file name and to
 * load the corresponding response.
 *
 * If @p rspname is not a FITS file the method searches for an appropriate
 * response in the calibration database. If no appropriate response is found,
 * the method takes the CALDB root path and response name to build the full
 * path to the response file, and tries to load the response from these
 * paths.
 *
 * If also this fails an exception is thrown.
 ***************************************************************************/
void GCOMResponse::load(const std::string& rspname)
{
    // Clear instance but conserve calibration database
    GCaldb caldb = m_caldb;
    clear();
    m_caldb = caldb;
 
    // Save response name
    m_rspname = rspname;
 
    // Interpret response name as a FITS file name
    GFilename filename(rspname);
 
    // If the filename does not exist the try getting the response from the
    // calibration database
    if (!filename.is_fits()) {
 
        // Get GCaldb response
        filename = m_caldb.filename("","","IAQ","","",rspname);
 
        // If filename is empty then build filename from CALDB root path
        // and response name
        if (filename.is_empty()) {
            filename = gammalib::filepath(m_caldb.rootdir(), rspname);
            if (!filename.exists()) {
                GFilename testname = filename + ".fits";
                if (testname.exists()) {
                    filename = testname;
                }
            }
        }
 
    } // endif: response name is not a FITS file
 
    // Open FITS file
    GFits fits(filename);
 
    // Get IAQ image
    const GFitsImage& iaq = *fits.image(0);
 
    // Read IAQ
    read(iaq);
 
    // Close IAQ FITS file
    fits.close();
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Read COMPTEL response from FITS image.
 *
 * @param[in] image FITS image.
 *
 * Read the COMPTEL response from IAQ FITS file and convert the IAQ values
 * into a probability per steradian.
 *
 * The IAQ values are divided by the solid angle of a Phigeo bin which is
 * given by
 *
 * \f{eqnarray*}{
 *    \Omega & = & 2 \pi \left[
 *      \left(
 *        1 - \cos \left( \varphi_{\rm geo} +
 *                        \frac{1}{2} \Delta \varphi_{\rm geo} \right)
 *      \right) -
 *      \left(
 *        1 - \cos \left( \varphi_{\rm geo} -
 *                        \frac{1}{2} \Delta \varphi_{\rm geo} \right)
 *      \right) \right] \\
 *     &=& 2 \pi \left[
 *        \cos \left( \varphi_{\rm geo} -
 *                    \frac{1}{2} \Delta \varphi_{\rm geo} \right) -
 *        \cos \left( \varphi_{\rm geo} +
 *                    \frac{1}{2} \Delta \varphi_{\rm geo} \right)
 *      \right] \\
 *     &=& 4 \pi \sin \left( \varphi_{\rm geo} \right)
 *             \sin \left( \frac{1}{2} \Delta \varphi_{\rm geo} \right)
 * \f}
 ***************************************************************************/
void GCOMResponse::read(const GFitsImage& image)
{
    // Continue only if there are two image axis
    if (image.naxis() == 2) {
 
        // Store IAQ dimensions
        m_phigeo_bins = image.naxes(0);
        m_phibar_bins = image.naxes(1);
 
        // Store IAQ axes definitions
        m_phigeo_ref_value = image.real("CRVAL1");
        m_phigeo_ref_pixel = image.real("CRPIX1");
        m_phigeo_bin_size  = image.real("CDELT1");
        m_phibar_ref_value = image.real("CRVAL2");
        m_phibar_ref_pixel = image.real("CRPIX2");
        m_phibar_bin_size  = image.real("CDELT2");
 
        // Get axes minima (values of first bin)
        m_phigeo_min = m_phigeo_ref_value + (1.0-m_phigeo_ref_pixel) *
                       m_phigeo_bin_size;
        m_phibar_min = m_phibar_ref_value + (1.0-m_phibar_ref_pixel) *
                       m_phibar_bin_size;
 
        // Compute IAQ size. Continue only if size is positive
        int size = m_phigeo_bins * m_phibar_bins;
        if (size > 0) {
 
            // Allocate memory for IAQ
            m_iaq.assign(size, 0.0);
 
            // Copy over IAQ values
            for (int i = 0; i < size; ++i) {
                m_iaq[i] = image.pixel(i);
            }
 
        } // endif: size was positive
 
        // Precompute variable for conversion
        double phigeo_min = m_phigeo_min      * gammalib::deg2rad;
        double phigeo_bin = m_phigeo_bin_size * gammalib::deg2rad;
        double omega0     = gammalib::fourpi  * std::sin(0.5 * phigeo_bin);
 
        // Convert IAQ matrix from probability per Phigeo bin into a
        // probability per steradian
        for (int iphigeo = 0; iphigeo < m_phigeo_bins; ++iphigeo) {
            double phigeo = iphigeo * phigeo_bin + phigeo_min;
            double omega  = omega0  * std::sin(phigeo);
            for (int iphibar = 0; iphibar < m_phibar_bins; ++iphibar) {
                m_iaq[iphigeo+iphibar*m_phigeo_bins] /= omega;
            }
        }
 
        // Compute FAQ
        compute_faq();
 
    } // endif: image had two axes
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Write COMPTEL response into FITS image.
 *
 * @param[in] image FITS image.
 *
 * Writes the COMPTEL response into an IAQ FITS file.
 ***************************************************************************/
void GCOMResponse::write(GFitsImageFloat& image) const
{
    // Continue only if response is not empty
    if (m_phigeo_bins > 0 && m_phibar_bins > 0) {
 
        // Initialise image
        image = GFitsImageFloat(m_phigeo_bins, m_phibar_bins);
 
        // Convert IAQ matrix from probability per steradian into
        //probability per Phigeo bin
        double omega0 = gammalib::fourpi *
                        std::sin(0.5 * m_phigeo_bin_size * gammalib::deg2rad);
        for (int iphigeo = 0; iphigeo < m_phigeo_bins; ++iphigeo) {
            double phigeo = iphigeo * m_phigeo_bin_size + m_phigeo_min;
            double omega  = omega0 * std::sin(phigeo * gammalib::deg2rad);
            for (int iphibar = 0; iphibar < m_phibar_bins; ++iphibar) {
                int i          = iphigeo + iphibar*m_phigeo_bins;
                image(iphigeo, iphibar) = m_iaq[i] * omega;
            }
        }
 
        // Set header keywords
        image.card("CTYPE1", "Phigeo", "Geometrical scatter angle");
        image.card("CRVAL1", m_phigeo_ref_value,
                   "[deg] Geometrical scatter angle for reference bin");
        image.card("CDELT1", m_phigeo_bin_size,
                   "[deg] Geometrical scatter angle bin size");
        image.card("CRPIX1", m_phigeo_ref_pixel,
                   "Reference bin of geometrical scatter angle");
        image.card("CTYPE2", "Phibar", "Compton scatter angle");
        image.card("CRVAL2", m_phibar_ref_value,
                   "[deg] Compton scatter angle for reference bin");
        image.card("CDELT2", m_phibar_bin_size, "[deg] Compton scatter angle bin size");
        image.card("CRPIX2", m_phibar_ref_pixel,
                   "Reference bin of Compton scatter angle");
        image.card("BUNIT", "Probability per bin", "Unit of bins");
        image.card("TELESCOP", "GRO", "Name of telescope");
        image.card("INSTRUME", "COMPTEL", "Name of instrument");
        image.card("ORIGIN", "GammaLib", "Origin of FITS file");
        image.card("OBSERVER", "Unknown", "Observer that created FITS file");
 
    } // endif: response was not empty
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Load response cache.
 *
 * @param[in] filename Response cache filename.
 *
 * Loads response cache from FITS file.
 ***************************************************************************/
void GCOMResponse::load_cache(const GFilename& filename)
{
    // Load response vector cache
    m_irf_vector_cache.load(filename);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Save response cache.
 *
 * @param[in] filename Response cache filename.
 *
 * Saves response cache from FITS file.
 ***************************************************************************/
void GCOMResponse::save_cache(const GFilename& filename) const
{
    // Save response vector cache
    m_irf_vector_cache.save(filename);
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Backproject events using instrument response function to sky map
 *
 * @param[in] obs Observation.
 * @param[in] events Events.
 * @param[in,out] map Sky map.
 *
 * @exception GException::invalid_argument
 *            Events pointer invalid.
 *            Sky map pointer invalid.
 *            Observation is not a COMPTEL observation.
 *            Events are not a COMPTEL event cube.
 * @exception GException::invalid_value
 *            Response not initialised with a valid IAQ
 *            Incompatible IAQ encountered
 *
 * Backprojects events using the instrument response function to sky map
 * using
 *
 * \f[
 *    f_j = \sum_i n_i R_{ij}
 * \f]
 *
 * where
 * - \f$f_j\f$ is the sky map pixel \f$j\f$,
 * - \f$n_i\f$ is the event bin \f$i\f$, and
 * - \f$R_{ij}\f$ is the response function for event bin \f$i\f$ and
 *   sky map pixel \f$j\f$.
 *
 * If the sky map contains at least two maps, the second map will contain
 * on output the normalisation map
 *
 * \f[
 *    f_j = \sum_i R_{ij}
 * \f]
 ***************************************************************************/
void GCOMResponse::backproject(const GObservation& obs,
                               const GEvents*      events,
                               GSkyMap*            map) const
{
    // Throw an exception if the events pointer is invalid
    if (events == NULL) {
        std::string msg = "Invalid events pointer. Please specify a "
                          "valid COMPTEL events pointer as argument.";
        throw GException::invalid_argument(G_BACKPROJECT, msg);
    }
 
    // Throw an exception if the sky map pointer is invalid
    if (map == NULL) {
        std::string msg = "Invalid sky map pointer. Please specify a "
                          "pointer to a valid sky map as argument.";
        throw GException::invalid_argument(G_BACKPROJECT, msg);
    }
 
    // Extract COMPTEL observation
    const GCOMObservation* obs_ptr = dynamic_cast<const GCOMObservation*>(&obs);
    if (obs_ptr == NULL) {
        std::string cls = std::string(typeid(&obs).name());
        std::string msg = "Observation of type \""+cls+"\" is not a COMPTEL "
                          "observations. Please specify a COMPTEL observation "
                          "as argument.";
        throw GException::invalid_argument(G_BACKPROJECT, msg);
    }
 
    // Extract COMPTEL event cube
    const GCOMEventCube* cube = dynamic_cast<const GCOMEventCube*>(events);
    if (cube == NULL) {
        std::string cls = std::string(typeid(*events).name());
        std::string msg = "Events of type \""+cls+"\" is not a COMPTEL event "
                          "cube. Please specify an event cube.";
        throw GException::invalid_argument(G_BACKPROJECT, msg);
    }
 
    // Throw an exception if COMPTEL response is not set or if
    if (m_iaq.empty()) {
        std::string msg = "COMPTEL response is empty. Please initialise the "
                          "response with an \"IAQ\".";
        throw GException::invalid_value(G_BACKPROJECT, msg);
    }
    else if (m_phigeo_bin_size == 0.0) {
        std::string msg = "COMPTEL response has a zero Phigeo bin size. "
                          "Please initialise the response with a valid "
                          "\"IAQ\".";
        throw GException::invalid_value(G_BACKPROJECT, msg);
    }
 
    // Get const reference to DRI
    const GCOMDri& dri = cube->dre();
 
    // Get number of Chi/Psi pixels, Phibar layers and event bins
    int npix    = dri.nchi() * dri.npsi();
    int nphibar = dri.nphibar();
    int nevents = dri.size();
 
    // Throw an exception if the number of Phibar bins does not match the
    // response
    if (nphibar != m_phibar_bins) {
        std::string msg = "DRI has "+gammalib::str(nphibar)+" Phibar layers "
                          "but IAQ has "+gammalib::str(m_phibar_bins)+" Phibar "
                          "bins. Please specify a compatible IAQ.";
        throw GException::invalid_value(G_BACKPROJECT, msg);
    }
 
    // Initialise some variables
    double         phigeo_min = m_phigeo_min      * gammalib::deg2rad;
    double         phigeo_bin = m_phigeo_bin_size * gammalib::deg2rad;
    double         omega0     = 2.0 * std::sin(0.5 * phigeo_bin);
    const GSkyMap& drx        = obs_ptr->drx().map();
    const double*  drg        = obs_ptr->drg().map().pixels();
 
    // Compute IAQ normalisation (1/s): DEADC / ONTIME (s)
    double iaq_norm = obs_ptr->deadc() /
                      (obs_ptr->ontime() * cube->dre().tof_correction()) *
                      cube->dre().phase_correction();
 
    // Pre-compute sky directions for all pixels and DRX values
    std::vector<GSkyDir> dirs;
    std::vector<double>  drxs;
    dirs.reserve(map->npix());
    drxs.reserve(map->npix());
    for (int isky = 0; isky < map->npix(); ++isky) {
 
        // Compute sky direction
        GSkyDir skyDir = map->inx2dir(isky);
 
        // Push back sky direction
        dirs.push_back(skyDir);
 
        // Push back DRX value
        if (drx.contains(skyDir)) {
            drxs.push_back(drx(skyDir));
        }
        else {
            drxs.push_back(0.0);
        }
 
    } // endfor: precomputation
 
    // Signal presence of normalisation map
    bool norm = (map->nmaps() > 1);
 
    // Loop over Chi and Psi pixels
    for (int ipix = 0; ipix < npix; ++ipix) {
 
        // Get sky direction to (Chi,Psi)
        GSkyDir chipsi = dri.map().inx2dir(ipix);
 
        // Loop over sky map pixels
        for (int isky = 0; isky < map->npix(); ++isky) {
 
            // Get sky direction to skymap pixel
            GSkyDir& skyDir = dirs[isky];
 
            // Get distance between (Chi,Psi) and sky direction
            double phigeo = chipsi.dist_deg(skyDir);
 
            // Compute interpolation factors
            double phirat  = phigeo / m_phigeo_bin_size; // 0.5 at bin centre
            int    iphigeo = int(phirat);                // index into which Phigeo falls
            double eps     = phirat - iphigeo - 0.5;     // 0.0 at bin centre [-0.5, 0.5[
 
            // Fall through if Phigeo is outside range
            if (iphigeo >= m_phigeo_bins) {
                continue;
            }
 
            // Get DRX value (units: cm^2 s)
            double drx_value = drxs[isky];
 
            // Loop over Phibar
            for (int iphibar = 0; iphibar < nphibar; ++iphibar) {
 
                // Get DRI index
                int idri = ipix + iphibar * npix;
 
                // Fall through if Chi and Psi pixel is empty
                if (dri[idri] == 0.0) {
                    continue;
                }
 
                // Get IAQ index
                int i = iphibar * m_phigeo_bins + iphigeo;
 
                // Initialise IAQ
                double iaq = 0.0;
 
                // Compute IAQ
                if (eps < 0.0) { // interpolate towards left
                    if (iphigeo > 0) {
                        iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
                    }
                    else {
                        iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
                    }
                }
                else {           // interpolate towards right
                    if (iphigeo < m_phigeo_bins-1) {
                        iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
                    }
                    else {
                        iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
                    }
                }
 
                // Fall through if IAQ is not positive
                if (iaq <= 0.0) {
                    continue;
                }
 
                // Compute IRF value (units: photons -> events)
                double irf = iaq * drg[idri] * iaq_norm * drx_value;
 
                // Add IRF value if it is positive
                if (irf > 0.0) {
                    (*map)(isky) += irf * dri[idri];
                    if (norm) {
                        (*map)(isky,1) += irf;
                    }
                }
 
            } // endif: looped over Phibar
 
        } // endfor: looped over sky map pixels
 
    } // endfor: looped over Chi and Psi pixels
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Print COMPTEL response information
 *
 * @param[in] chatter Chattiness.
 * @return String containing COMPTEL response information.
 ***************************************************************************/
std::string GCOMResponse::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("=== GCOMResponse ===");
 
        // Append response name
        result.append("\n"+gammalib::parformat("Response name")+m_rspname);
 
        // EXPLICIT: Append detailed information
        if (chatter >= EXPLICIT) {
 
            // Append information
            result.append("\n"+gammalib::parformat("Number of Phigeo bins"));
            result.append(gammalib::str(m_phigeo_bins));
            result.append("\n"+gammalib::parformat("Number of Phibar bins"));
            result.append(gammalib::str(m_phibar_bins));
            result.append("\n"+gammalib::parformat("Phigeo reference value"));
            result.append(gammalib::str(m_phigeo_ref_value)+" deg");
            result.append("\n"+gammalib::parformat("Phigeo reference pixel"));
            result.append(gammalib::str(m_phigeo_ref_pixel));
            result.append("\n"+gammalib::parformat("Phigeo bin size"));
            result.append(gammalib::str(m_phigeo_bin_size)+" deg");
            result.append("\n"+gammalib::parformat("Phigeo first bin value"));
            result.append(gammalib::str(m_phigeo_min)+" deg");
            result.append("\n"+gammalib::parformat("Phibar reference value"));
            result.append(gammalib::str(m_phibar_ref_value)+" deg");
            result.append("\n"+gammalib::parformat("Phibar reference pixel"));
            result.append(gammalib::str(m_phibar_ref_pixel));
            result.append("\n"+gammalib::parformat("Phibar bin size"));
            result.append(gammalib::str(m_phibar_bin_size)+" deg");
            result.append("\n"+gammalib::parformat("Phibar first bin value"));
            result.append(gammalib::str(m_phibar_min)+" deg");
 
        }
 
        // VERBOSE: Append calibration database
        if (chatter == VERBOSE) {
            result.append("\n"+m_caldb.print(chatter));
        }
 
    } // endif: chatter was not silent
 
    // Return result
    return result;
}
 
 
/*==========================================================================
 =                                                                         =
 =                             Private methods                             =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Initialise class members
 ***************************************************************************/
void GCOMResponse::init_members(void)
{
    // Initialise members
    m_caldb.clear();
    m_rspname.clear();
    m_iaq.clear();
    m_faq.clear();
    m_faq_zero.clear();
    m_faq_native.clear();
    m_faq_solidangle.clear();
    m_phigeo_bins      = 0;
    m_phibar_bins      = 0;
    m_faq_bins         = 0;
    m_phigeo_ref_value = 0.0;
    m_phigeo_ref_pixel = 0.0;
    m_phigeo_bin_size  = 0.0;
    m_phigeo_min       = 0.0;
    m_phibar_ref_value = 0.0;
    m_phibar_ref_pixel = 0.0;
    m_phibar_bin_size  = 0.0;
    m_phibar_min       = 0.0;
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Copy class members
 *
 * @param[in] rsp COMPTEL response.
 ***************************************************************************/
void GCOMResponse::copy_members(const GCOMResponse& rsp)
{
    // Copy attributes
    m_caldb            = rsp.m_caldb;
    m_rspname          = rsp.m_rspname;
    m_iaq              = rsp.m_iaq;
    m_faq              = rsp.m_faq;
    m_faq_zero         = rsp.m_faq_zero;
    m_faq_native       = rsp.m_faq_native;
    m_faq_solidangle   = rsp.m_faq_solidangle;
    m_phigeo_bins      = rsp.m_phigeo_bins;
    m_phibar_bins      = rsp.m_phibar_bins;
    m_faq_bins         = rsp.m_faq_bins;
    m_phigeo_ref_value = rsp.m_phigeo_ref_value;
    m_phigeo_ref_pixel = rsp.m_phigeo_ref_pixel;
    m_phigeo_bin_size  = rsp.m_phigeo_bin_size;
    m_phigeo_min       = rsp.m_phigeo_min;
    m_phibar_ref_value = rsp.m_phibar_ref_value;
    m_phibar_ref_pixel = rsp.m_phibar_ref_pixel;
    m_phibar_bin_size  = rsp.m_phibar_bin_size;
    m_phibar_min       = rsp.m_phibar_min;
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Delete class members
 ***************************************************************************/
void GCOMResponse::free_members(void)
{
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Compute FAQ
 *
 * Computes FAQ from IAQ that is needed for direction projection computation
 ***************************************************************************/
void GCOMResponse::compute_faq(void)
{
    // Initialise FAQ members
    m_faq.clear();
    m_faq_zero.clear();
    m_faq_solidangle.clear();
    m_faq_bins = 0;
 
    // Continue only if there are phigeo and phibar bins
    if ((m_phigeo_bins > 0) && (m_phibar_bins > 0)) {
 
        // Determine number of FAQ bins
        int bins = 2 * m_phigeo_bins - 1;
 
        // Initialise vectors
        m_faq_bins       = bins * bins;
        m_faq_zero       = std::vector<bool>(m_faq_bins, true);
        m_faq_native     = std::vector<GVector>(m_faq_bins);
        m_faq_solidangle = GVector(m_faq_bins);
 
        // Initialise FAQ skymap in celestial coordinates
        m_faq = GSkyMap("ARC", "CEL", 0.0, 0.0,
                        m_phigeo_bin_size, m_phigeo_bin_size,
                        bins, bins, m_phibar_bins);
 
        // Loop over Psi
        for (int ipsi = 0, ifaq = 0; ipsi < bins; ++ipsi) {
 
            // Compute Psi
            double psi     = (ipsi - m_phigeo_bins + 1) * m_phigeo_bin_size;
            double psi_rad = psi * gammalib::deg2rad;
            double psi2    = psi * psi;
 
            // Loop over Chi
            for (int ichi = 0; ichi < bins; ++ichi, ++ifaq) {
 
                // Compute Chi
                double chi     = (ichi - m_phigeo_bins + 1) * m_phigeo_bin_size;
                double chi_rad = chi * gammalib::deg2rad;
 
                // Compute phigeo (degrees)
                double phigeo = std::sqrt(chi*chi + psi2);
 
                // Compute native FAQ vector
                double zenith      = phigeo * gammalib::deg2rad;
                double cos_zenith  = std::cos(zenith);
                double sin_zenith  = std::sin(zenith);
                double azimuth     = 0.0 - std::atan2(psi_rad, chi_rad);
                double cos_phi     = std::cos(azimuth);
                double sin_phi     = std::sin(azimuth);
                m_faq_native[ifaq] = GVector(-cos_phi * sin_zenith,
                                              sin_phi * sin_zenith,
                                              cos_zenith);
 
                // Store solid angle of pixel
                m_faq_solidangle[ifaq] = m_faq.solidangle(ifaq);
 
                // Compute interpolation factors
                double phirat  = phigeo / m_phigeo_bin_size; // 0.5 at bin centre
                int    iphigeo = int(phirat);                // index into which Phigeo falls
                double eps     = phirat - iphigeo - 0.5;     // 0.0 at bin centre [-0.5, 0.5[
 
                // If Phigeo is in range then compute the IRF value for all
                // Phibar layers
                if (iphigeo < m_phigeo_bins) {
 
                    // Loop over Phibar
                    for (int iphibar = 0; iphibar < m_phibar_bins; ++iphibar) {
 
                        // Get IAQ index
                        int i = iphibar * m_phigeo_bins + iphigeo;
 
                        // Initialise FAQ
                        double faq = 0.0;
 
                        // Compute IAQ
                        if (eps < 0.0) { // interpolate towards left
                            if (iphigeo > 0) {
                                faq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
                            }
                            else {
                                faq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
                            }
                        }
                        else {           // interpolate towards right
                            if (iphigeo < m_phigeo_bins-1) {
                                faq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
                            }
                            else {
                                faq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
                            }
                        }
 
                        // Continue only if FAQ is positive
                        if (faq > 0.0) {
 
                            // Normalise FAQ value
                            //faq *= m_phigeo_bin_size/(gammalib::twopi * phigeo);
                            //faq *= m_phigeo_bin_size/(gammalib::twopi * phigeo)/m_faq.solidangle(ifaq);
                            //faq /= m_faq.solidangle(ifaq);
 
                            // Store FAQ value
                            m_faq(ifaq, iphibar) = faq;
 
                            // Signal that FAQ bin is not zero
                            m_faq_zero[ifaq] = false;
 
                        } // endif: faq was positive
 
                    } // endfor: looped over Phibar
 
                } // endif: Phigeo was in valid range
 
            } // endfor: looped over Psi
 
        } // endfor: looped over Chi
 
        // Debug FAQ computation
        #if defined(G_DEBUG_COMPUTE_FAQ)
        std::cout << "GCOMResponse::compute_faq" << std::endl;
        m_faq.save("faq.fits", true);
        for (int iphibar = 0; iphibar < m_phibar_bins; ++iphibar) {
            double iaq_sum = 0.0;
            double faq_sum = 0.0;
            for (int iphigeo = 0; iphigeo < m_phigeo_bins; ++iphigeo) {
                int i = iphibar * m_phigeo_bins + iphigeo;
                iaq_sum += m_iaq[i];
            }
            for (int ipix = 0; ipix < m_faq.npix(); ++ipix) {
                faq_sum += m_faq(ipix, iphibar);
            }
            std::cout << "iphibar=" << iphibar;
            std::cout << " iaq=" << iaq_sum;
            std::cout << " faq=" << faq_sum;
            std::cout << " ratio=" << faq_sum/iaq_sum << std::endl;
        }
        #endif
 
    } // endif: there were phigeo and phibar bins
 
    // Return
    return;
}
 
 
/***********************************************************************//**
 * @brief Return instrument response to point source
 *
 * @param[in] model Sky model.
 * @param[in] obs Observation.
 * @param[out] gradients Gradients matrix.
 * @return Instrument response to point source for all events in
 *         observation (\f$cm^2\f$).
 *
 * @exception GException::invalid_argument
 *            Observation is not a COMPTEL observation.
 *            Event is not a COMPTEL event bin.
 * @exception GException::invalid_value
 *            Response not initialised with a valid IAQ
 *            Incompatible IAQ encountered
 *
 * Returns the instrument response to a point source for all events in the
 * observations.
 *
 * @p gradients is an optional matrix where the number of rows corresponds
 * to the number of events in the observation and the number of columns
 * corresponds to the number of spatial model parameters. Since for point
 * sources no gradients are computed, the method does not alter the
 * content of @p gradients.
 ***************************************************************************/
GVector GCOMResponse::irf_ptsrc(const GModelSky&    model,
                                const GObservation& obs,
                                GMatrix*            gradients) const
{
    // Extract COMPTEL observation
    const GCOMObservation* obs_ptr = dynamic_cast<const GCOMObservation*>(&obs);
    if (obs_ptr == NULL) {
        std::string cls = std::string(typeid(&obs).name());
        std::string msg = "Observation of type \""+cls+"\" is not a COMPTEL "
                          "observations. Please specify a COMPTEL observation "
                          "as argument.";
        throw GException::invalid_argument(G_IRF_PTSRC, msg);
    }
 
    // Extract COMPTEL event cube
    const GCOMEventCube* cube = dynamic_cast<const GCOMEventCube*>(obs_ptr->events());
    if (cube == NULL) {
        std::string msg = "Observation \""+obs.name()+"\" ("+obs.id()+") does "
                          "not contain a COMPTEL event cube. Please specify "
                          "a COMPTEL observation containing and event cube.";
        throw GException::invalid_argument(G_IRF_PTSRC, msg);
    }
 
    // Throw an exception if COMPTEL response is not set or if
    if (m_iaq.empty()) {
        std::string msg = "COMPTEL response is empty. Please initialise the "
                          "response with an \"IAQ\".";
        throw GException::invalid_value(G_IRF_PTSRC, msg);
    }
    else if (m_phigeo_bin_size == 0.0) {
        std::string msg = "COMPTEL response has a zero Phigeo bin size. "
                          "Please initialise the response with a valid "
                          "\"IAQ\".";
        throw GException::invalid_value(G_IRF_PTSRC, msg);
    }
 
    // Get number of Chi/Psi pixels, Phibar layers and event bins
    int npix    = cube->naxis(0) * cube->naxis(1);
    int nphibar = cube->naxis(2);
    int nevents = cube->size();
 
    // Throw an exception if the number of Phibar bins does not match the
    // response
    if (nphibar != m_phibar_bins) {
        std::string msg = "DRE has "+gammalib::str(nphibar)+" Phibar layers "
                          "but IAQ has "+gammalib::str(m_phibar_bins)+" Phibar "
                          "bins. Please specify a compatible IAQ.";
        throw GException::invalid_value(G_IRF_PTSRC, msg);
    }
 
    // Initialise result
    GVector irfs(nevents);
 
    // Get point source direction
    GSkyDir srcDir =
    static_cast<const GModelSpatialPointSource*>(model.spatial())->dir();
 
    // Get IAQ normalisation (cm2): DRX (cm2 s) * DEADC / ONTIME (s)
    double iaq_norm = obs_ptr->drx().map()(srcDir) * obs_ptr->deadc() /
                      (obs_ptr->ontime() * cube->dre().tof_correction()) *
                      cube->dre().phase_correction();
 
    // Get pointer to DRG pixels
    const double* drg = obs_ptr->drg().map().pixels();
 
    // Loop over Chi and Psi
    for (int ipix = 0; ipix < npix; ++ipix) {
 
        // Get pointer to event bin
        const GCOMEventBin* bin = (*cube)[ipix];
 
        // Get reference to instrument direction
        const GCOMInstDir& obsDir = bin->dir();
 
        // Get reference to Phigeo sky direction
        const GSkyDir& phigeoDir = obsDir.dir();
 
        // Compute angle between true photon arrival direction and scatter
        // direction (Chi,Psi)
        double phigeo = srcDir.dist_deg(phigeoDir);
 
        // Compute interpolation factors
        double phirat  = phigeo / m_phigeo_bin_size; // 0.5 at bin centre
        int    iphigeo = int(phirat);                // index into which Phigeo falls
        double eps     = phirat - iphigeo - 0.5;     // 0.0 at bin centre [-0.5, 0.5[
 
        // If Phigeo is in range then compute the IRF value for all
        // Phibar layers
        if (iphigeo < m_phigeo_bins) {
 
            // Loop over Phibar
            for (int iphibar = 0; iphibar < nphibar; ++iphibar) {
 
                // Get IAQ index
                int i = iphibar * m_phigeo_bins + iphigeo;
 
                // Initialise IAQ
                double iaq = 0.0;
 
                // Compute IAQ
                if (eps < 0.0) { // interpolate towards left
                    if (iphigeo > 0) {
                        iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
                    }
                    else {
                        iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
                    }
                }
                else {           // interpolate towards right
                    if (iphigeo < m_phigeo_bins-1) {
                        iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
                    }
                    else {
                        iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
                    }
                }
 
                // Continue only if IAQ is positive
                if (iaq > 0.0) {
 
                    // Get DRI index
                    int idri = ipix + iphibar * npix;
 
                    // Compute IRF value
                    double irf = iaq * drg[idri] * iaq_norm;
 
                    // Make sure that IRF is positive
                    if (irf < 0.0) {
                        irf = 0.0;
                    }
 
                    // Store IRF value
                    irfs[idri] = irf;
 
                } // endif: IAQ was positive
 
            } // endfor: looped over Phibar
 
        } // endif: Phigeo was in valid range
 
    } // endfor: looped over Chi and Psi
 
    // Return IRF vector
    return irfs;
}
 
 
/***********************************************************************//**
 * @brief Return instrument response to radial source
 *
 * @param[in] model Sky model.
 * @param[in] obs Observation.
 * @param[out] gradients Gradients matrix.
 * @return Instrument response to radial source for all events in
 *         observation (\f$cm^2\f$).
 *
 * @exception GException::invalid_value
 *            Response not initialised
 *            Incompatible IAQ encountered
 *
 * Returns the instrument response to a radial source sky model for all
 * events. The methods works on response vectors that are computed per
 * pointing. The method assumes that the spatial model component is of
 * type GModelSpatialRadial, that the observations is of type
 * GCOMObservation and that the events are of type GCOMEventCube.
 * No checking of these assumptions is performed by the method, and not
 * respecting the assumptions will results in a segfault.
 *
 * The method toggles between two convolution methods that are used depending
 * on the extent of the radial source. Extent means here the value of the
 * radial parameter, which is the radius for a disk or the sigma parameter
 * for a Gaussian.
 *
 * For extents smaller than 10 degrees a numerical integration is performed
 * in the system of the radial model that is fast but that may become
 * inaccurate for larger extents as the model sampling will exceed the
 * resolution of the data space. The number of iterations in the zenith and
 * azimuth angle was fixed to 6.
 *
 * For extents larger than 15 degrees a direct convolution is performed in
 * the system of the instrument response function, that is slower, but that
 * samples fully the resolution of the data space. The direct convolution is
 * however not accurate for angles smaller than about 1 degree, hence it is
 * not appropriate to compute for the response for the determination of
 * upper limits.
 *
 * For intermediate extents (in the interval 10 - 15 degrees), both
 * convolutions are computed and combined according to a weight that varies
 * linearly with source extent. In principle the convolved results should be
 * identical in that intermediate range, but computing both and their
 * weighted combination assures that no steps occur when varying the source
 * extent, which is required for seamless log-likelihood fits.
 *
 * For details about the choice of the convolution method and its parameters
 * see https://gitlab.in2p3.fr/gammalib/gammalib/-/issues/14.
 *
 * @p gradients is an optional matrix where the number of rows corresponds
 * to the number of events in the observation and the number of columns
 * corresponds to the number of spatial model parameters. Since for
 * radial sources no gradients are computed, the method does not alter the
 * content of @p gradients.
 ***************************************************************************/
GVector GCOMResponse::irf_radial(const GModelSky&    model,
                                 const GObservation& obs,
                                 GMatrix*            gradients) const
{
    // Set constants
    static const int    iter_rho           = 6;    //!< 6 iterations
    static const int    iter_omega         = 6;    //!< 6 iterations
    static const double extent_integration = 10.0; //!< Integration for extent < 10 deg
    static const double extent_direct      = 15.0; //!< Direct for extent > 15 deg
 
    // Get pointers
    const GModelSpatialRadial* radial  = static_cast<const GModelSpatialRadial*>(model.spatial());
    const GCOMObservation*     obs_ptr = static_cast<const GCOMObservation*>(&obs);
    const GCOMEventCube*       cube    = static_cast<const GCOMEventCube*>(obs.events());
 
    // Set weights dependent on radial method
    #if G_IRF_RADIAL_METHOD == 0
    double wgt_integration = 1.0;
    double wgt_direct      = 0.0;
    #elif G_IRF_RADIAL_METHOD == 1
    double wgt_integration = 0.0;
    double wgt_direct      = 1.0;
    #else
    double extent          = (*radial)[2].value(); // Assume 3rd parameter is extent
    double wgt_direct      = (extent-extent_integration)/(extent_direct-extent_integration);
    double wgt_integration = 1.0 - wgt_direct;
    if (wgt_direct < 0.0) {
        wgt_direct      = 0.0;
        wgt_integration = 1.0;
    }
    else if (wgt_direct > 1.0) {
        wgt_direct      = 1.0;
        wgt_integration = 0.0;
    }
    #endif
 
    // Timing option: store entry time and log parameters
    #if defined(G_DEBUG_IRF_TIMING)
    double cstart = gammalib::get_current_clock();
    std::cout << "GCOMResponse::irf_radial(";
    std::cout << "iter_rho=" << iter_rho;
    std::cout << ",iter_omega=" << iter_omega;
    #if G_IRF_RADIAL_METHOD == 2
    std::cout << ",extent=" << extent;
    std::cout << ",extent_integration=" << extent_integration;
    std::cout << ",extent_direct=" << extent_direct;
    #endif
    std::cout << ",wgt_direct=" << wgt_direct;
    std::cout << ",wgt_integration=" << wgt_integration << ")" << std::endl;
    #endif
 
    // Get number of Chi/Psi pixels, Phibar layers and event bins
    int npix    = cube->naxis(0) * cube->naxis(1);
    int nphibar = cube->naxis(2);
    int nevents = cube->size();
 
    // Throw an exception if the number of Phibar bins does not match the
    // response
    if (nphibar != m_phibar_bins) {
        std::string msg = "DRE has "+gammalib::str(nphibar)+" Phibar layers "
                          "but IAQ has "+gammalib::str(m_phibar_bins)+" Phibar "
                          "bins. Please specify a compatible IAQ.";
        throw GException::invalid_value(G_IRF_RADIAL, msg);
    }
 
    // Initialise result
    GVector irfs(nevents);
 
    // Initialise some variables
    double         phigeo_bin = m_phigeo_bin_size * gammalib::deg2rad;
    double         phigeo_min = m_phigeo_min      * gammalib::deg2rad - 0.5 * phigeo_bin;
    double         phigeo_max = phigeo_min + phigeo_bin * m_phigeo_bins;
    const GSkyMap& drx        = obs_ptr->drx().map();
    const double*  drg        = obs_ptr->drg().map().pixels();
 
    // Compute IAQ normalisation (1/s): DEADC / ONTIME (s)
    double iaq_norm = obs_ptr->deadc() /
                      (obs_ptr->ontime() * cube->dre().tof_correction()) *
                      cube->dre().phase_correction();
 
    // If weight for direct computation is non-zero then perform direct
    // computation and add response according to the weight to the IRF
    // vector
    if (wgt_direct > 0.0) {
 
        // Throw an exception if COMPTEL FAQ response is not set
        if (m_faq.is_empty()) {
            std::string msg = "COMPTEL response is empty. Please initialise the "
                              "response with an \"IAQ\".";
            throw GException::invalid_value(G_IRF_RADIAL, msg);
        }
 
        // Compute weigthed IAQ normalisation
        double norm = iaq_norm * wgt_direct;
 
        // Loop over Chi and Psi bins
        for (int ipix = 0; ipix < npix; ++ipix) {
 
            // Get pointer to event bin
            const GCOMEventBin* bin = (*cube)[ipix];
 
            // Get reference to instrument direction
            const GCOMInstDir& obsDir = bin->dir();
 
            // Get reference to scatter sky direction
            const GSkyDir& scatter = obsDir.dir();
 
            // Compute distance between scatter direction and centre of radial
            // model in radians
            double centre_distance = scatter.dist(radial->dir());
 
            // Compute response in case that radial model overlaps with
            // COMPTEL field of view
            if (centre_distance < (phigeo_max + radial->theta_max())) {
 
                // Setup rotation matrix for current scatter direction to
                // transform from FAQ system to the celestial coordinate
                // system
                GMatrix ry;
                GMatrix rz;
                ry.eulery(scatter.dec_deg() - 90.0);
                rz.eulerz(-scatter.ra_deg());
                GMatrix rot = (ry * rz).transpose();
 
                // Loop over FAQ pixels
                for (int inx = 0; inx < m_faq_bins; ++inx) {
 
                    // Skip zero FAQ pixels
                    if (m_faq_zero[inx]) {
                        continue;
                    }
 
                    // Get sky direction for FAQ pixel
                    GVector vector = rot * m_faq_native[inx];
                    GSkyDir dir(vector);
 
                    // Compute offset angle to centre of radial model
                    double theta = radial->dir().dist(dir);
 
                    // Continue only if radial offset angle value is valid
                    if (theta < radial->theta_max()) {
 
                        // Compute radial model value
                        double model = radial->eval(theta, bin->energy(), bin->time());
 
                        // Continue only if model value is positive
                        if (model > 0.0) {
 
                            // Multiply model value by solid angle of FAQ pixel
                            model *= m_faq_solidangle[inx];
 
                            // Continue only if sky direction is within DRX
                            if (drx.contains(dir)) {
 
                                // Multiply model value DRX value and norm
                                model *= drx(dir) * norm;
 
                                // Loop over all phibar layers
                                for (int iphibar = 0, idri = ipix; iphibar < nphibar; ++iphibar, idri += npix) {
                                    double faq = m_faq(inx, iphibar);
                                    if (faq > 0.0) {
                                        irfs[idri] += model * faq * drg[idri];
                                    }
                                }
 
                            } // endif: sky direction within DRX
 
                        } // endif: model value was positive
 
                    } // endif: radial offset angle was valid
 
                } // endfor: looped over FAQ pixels
 
            } // endif: radial model overlaped with COMPTEL field of view
 
        } // endfor: looped over Chi and Psi bins
 
    } // endif: direct response computation requested
 
    // If weight for numerical integration is non-zero then perform numerical
    // integration and add response according to the weight to the IRF
    // vector
    if (wgt_integration > 0.0) {
 
        // Throw an exception if COMPTEL response is not set or if
        if (m_iaq.empty()) {
            std::string msg = "COMPTEL response is empty. Please initialise the "
                              "response with an \"IAQ\".";
            throw GException::invalid_value(G_IRF_RADIAL, msg);
        }
        else if (m_phigeo_bin_size == 0.0) {
            std::string msg = "COMPTEL response has a zero Phigeo bin size. "
                              "Please initialise the response with a valid "
                              "\"IAQ\".";
            throw GException::invalid_value(G_IRF_RADIAL, msg);
        }
 
        // Compute weigthed IAQ normalisation
        double norm = iaq_norm * wgt_integration;
 
        // Get maximum model radius, reduced by a small amount to avoid
        // rounding problems at the boundary of a sharp edged model
        double theta_max = radial->theta_max() - 1.0e-12;
 
        // Setup rotation matrix for conversion towards sky coordinates
        GMatrix ry;
        GMatrix rz;
        ry.eulery(radial->dir().dec_deg() - 90.0);
        rz.eulerz(-radial->dir().ra_deg());
        GMatrix rot = (ry * rz).transpose();
 
        // Loop over Chi and Psi
        for (int ipix = 0; ipix < npix; ++ipix) {
 
            // Get pointer to event bin
            const GCOMEventBin* bin = (*cube)[ipix];
 
            // Get reference to instrument direction
            const GCOMInstDir& obsDir = bin->dir();
 
            // Get reference to Phigeo sky direction
            const GSkyDir& phigeoDir = obsDir.dir();
 
            // Compute angle between model centre and Phigeo sky direction (radians)
            // and position angle of model with respect to Phigeo sky direction (radians)
            double zeta = phigeoDir.dist(radial->dir());
 
            // Set radial model zenith angle range
            double rho_min = (zeta > phigeo_max) ? zeta - phigeo_max : 0.0;
            double rho_max = zeta + phigeo_max;
            if (rho_min > theta_max) {
                rho_min = theta_max;
            }
            if (rho_max > theta_max) {
                rho_max = theta_max;
            }
 
            // Perform model zenith angle integration if interval is valid
            if (rho_max > rho_min) {
 
                // Initialise IRF vector
                GVector irf(nphibar);
 
                // Setup integration kernel
                com_radial_kerns_rho integrands(m_iaq,
                                                *radial,
                                                irf,
                                                bin,
                                                rot,
                                                drx,
                                                phigeo_bin,
                                                m_phigeo_bins,
                                                nphibar,
                                                iter_omega);
 
                // Setup integrator
                GIntegrals integral(&integrands);
                integral.fixed_iter(iter_rho);
 
                // Integrate over Phigeo
                irf = integral.romberg(rho_min, rho_max, iter_rho);
 
                // Add IRF to result
                for (int iphibar = 0, idri = ipix; iphibar < nphibar; ++iphibar, idri += npix) {
                    irfs[idri] += irf[iphibar] * drg[idri] * norm;
                }
 
            } // endif: Phigeo angle interval was valid
 
        } // endfor: looped over Chi and Psi pixels
 
    } // endif: numerical integration computation requested
 
    // Dump CPU consumption
    #if defined(G_DEBUG_IRF_TIMING)
    double celapse = gammalib::get_current_clock() - cstart;
    std::cout << "GCOMResponse::irf_radial consumed " << celapse;
    std::cout  << " seconds of CPU time." << std::endl;
    #endif
 
    // Return IRF vector
    return irfs;
}
 
 
/***********************************************************************//**
 * @brief Return instrument response to elliptical source
 *
 * @param[in] model Sky model.
 * @param[in] obs Observation.
 * @param[out] gradients Gradients matrix.
 * @return Instrument response to elliptical source for all events in
 *         observation (\f$cm^2\f$).
 *
 * @exception GException::invalid_value
 *            Response not initialised with a valid IAQ
 *            Incompatible IAQ encountered
 *
 * Returns the instrument response to an elliptical source sky model for all
 * events. The methods works on response vectors that are computed per
 * pointing. The method assumes that the spatial model component is of
 * type GModelSpatialElliptical, that the observations is of type
 * GCOMObservation and that the events are of type GCOMEventCube.
 * No checking of these assumptions is performed by the method, and not
 * respecting the assumptions will results in a segfault.
 *
 * The method toggles between two convolution methods that are used depending
 * on the extent of the elliptical source. Extent means here the larger of
 * the minor and major "radius" parameters, which is the radius for a disk
 * or the sigma parameter for a Gaussian.
 *
 * For extents smaller than 10 degrees a numerical integration is performed
 * in a radial system around the elliptical model that is fast but that may
 * become inaccurate for larger extents as the model sampling will exceed the
 * resolution of the data space. The number of iterations in the zenith and
 * azimuth angle was fixed to 6.
 *
 * For extents larger than 15 degrees a direct convolution is performed in
 * the system of the instrument response function, that is slower, but that
 * samples fully the resolution of the data space. The direct convolution is
 * however not accurate for angles smaller than about 1 degree, hence it is
 * not appropriate to compute for the response for the determination of
 * upper limits.
 *
 * For intermediate extents (in the interval 10 - 15 degrees), both
 * convolutions are computed and combined according to a weight that varies
 * linearly with source extent. In principle the convolved results should be
 * identical in that intermediate range, but computing both and their
 * weighted combination assures that no steps occur when varying the source
 * extent, which is required for seamless log-likelihood fits.
 *
 * For details about the choice of the convolution method and its parameters
 * see https://gitlab.in2p3.fr/gammalib/gammalib/-/issues/14.
 *
 * @p gradients is an optional matrix where the number of rows corresponds
 * to the number of events in the observation and the number of columns
 * corresponds to the number of spatial model parameters. Since for
 * elliptical sources no gradients are computed, the method does not alter
 * the content of @p gradients.
 ***************************************************************************/
GVector GCOMResponse::irf_elliptical(const GModelSky&    model,
                                     const GObservation& obs,
                                     GMatrix*            gradients) const
{
    // Set constants
    static const int    iter_rho           = 6;    //!< 6 iterations
    static const int    iter_omega         = 6;    //!< 6 iterations
    static const double extent_integration = 10.0; //!< Integration for extent < 10 deg
    static const double extent_direct      = 15.0; //!< Direct for extent > 15 deg
 
    // Get pointers
    const GModelSpatialElliptical* elliptical = static_cast<const GModelSpatialElliptical*>(model.spatial());
    const GCOMObservation*         obs_ptr    = static_cast<const GCOMObservation*>(&obs);
    const GCOMEventCube*           cube       = static_cast<const GCOMEventCube*>(obs.events());
 
    // Set weights dependent on elliptical method
    #if G_IRF_ELLIPTICAL_METHOD == 0
    double wgt_integration = 1.0;
    double wgt_direct      = 0.0;
    #elif G_IRF_ELLIPTICAL_METHOD == 1
    double wgt_integration = 0.0;
    double wgt_direct      = 1.0;
    #else
    double extent_maj      = (*elliptical)[3].value(); // Assume 4th parameter is major radius
    double extent_min      = (*elliptical)[4].value(); // Assume 5th parameter is minor radius
    double extent          = (extent_maj > extent_min) ? extent_maj : extent_min;
    double wgt_direct      = (extent-extent_integration)/(extent_direct-extent_integration);
    double wgt_integration = 1.0 - wgt_direct;
    if (wgt_direct < 0.0) {
        wgt_direct      = 0.0;
        wgt_integration = 1.0;
    }
    else if (wgt_direct > 1.0) {
        wgt_direct      = 1.0;
        wgt_integration = 0.0;
    }
    #endif
 
    // Timing option: store entry time and log parameters
    #if defined(G_DEBUG_IRF_TIMING)
    double cstart = gammalib::get_current_clock();
    std::cout << "GCOMResponse::irf_elliptical(";
    std::cout << "iter_rho=" << iter_rho;
    std::cout << ",iter_omega=" << iter_omega;
    #if G_IRF_ELLIPTICAL_METHOD == 2
    std::cout << ",extent=" << extent;
    std::cout << ",extent_integration=" << extent_integration;
    std::cout << ",extent_direct=" << extent_direct;
    #endif
    std::cout << ",wgt_direct=" << wgt_direct;
    std::cout << ",wgt_integration=" << wgt_integration << ")" << std::endl;
    #endif
 
    // Get number of Chi/Psi pixels, Phibar layers and event bins
    int npix    = cube->naxis(0) * cube->naxis(1);
    int nphibar = cube->naxis(2);
    int nevents = cube->size();
 
    // Throw an exception if the number of Phibar bins does not match the
    // response
    if (nphibar != m_phibar_bins) {
        std::string msg = "DRE has "+gammalib::str(nphibar)+" Phibar layers "
                          "but IAQ has "+gammalib::str(m_phibar_bins)+" Phibar "
                          "bins. Please specify a compatible IAQ.";
        throw GException::invalid_value(G_IRF_ELLIPTICAL, msg);
    }
 
    // Initialise result
    GVector irfs(nevents);
 
    // Initialise some variables
    double         phigeo_bin = m_phigeo_bin_size * gammalib::deg2rad;
    double         phigeo_min = m_phigeo_min      * gammalib::deg2rad - 0.5 * phigeo_bin;
    double         phigeo_max = phigeo_min + phigeo_bin * m_phigeo_bins;
    const GSkyMap& drx        = obs_ptr->drx().map();
    const double*  drg        = obs_ptr->drg().map().pixels();
 
    // Compute IAQ normalisation (1/s): DEADC / ONTIME (s)
    double iaq_norm = obs_ptr->deadc() /
                      (obs_ptr->ontime() * cube->dre().tof_correction()) *
                      cube->dre().phase_correction();
 
    // If weight for direct computation is non-zero then perform direct
    // computation and add response according to the weight to the IRF
    // vector
    if (wgt_direct > 0.0) {
 
        // Throw an exception if COMPTEL FAQ response is not set
        if (m_faq.is_empty()) {
            std::string msg = "COMPTEL response is empty. Please initialise the "
                              "response with an \"IAQ\".";
            throw GException::invalid_value(G_IRF_ELLIPTICAL, msg);
        }
 
        // Compute weigthed IAQ normalisation
        double norm = iaq_norm * wgt_direct;
 
        // Loop over Chi and Psi bins
        for (int ipix = 0; ipix < npix; ++ipix) {
 
            // Get pointer to event bin
            const GCOMEventBin* bin = (*cube)[ipix];
 
            // Get reference to instrument direction
            const GCOMInstDir& obsDir = bin->dir();
 
            // Get reference to scatter sky direction
            const GSkyDir& scatter = obsDir.dir();
 
            // Compute distance between scatter direction and centre of radial
            // model in radians
            double centre_distance = scatter.dist(elliptical->dir());
 
            // Compute response in case that elliptical model overlaps with COMPTEL
            // field of view
            if (centre_distance < (phigeo_max + elliptical->theta_max())) {
 
                // Setup rotation matrix for current scatter direction to
                // transform from FAQ system to the celestial coordinate
                // system
                GMatrix ry;
                GMatrix rz;
                ry.eulery(scatter.dec_deg() - 90.0);
                rz.eulerz(-scatter.ra_deg());
                GMatrix rot = (ry * rz).transpose();
 
                // Loop over FAQ pixels
                for (int inx = 0; inx < m_faq_bins; ++inx) {
 
                    // Skip zero FAQ pixels
                    if (m_faq_zero[inx]) {
                        continue;
                    }
 
                    // Get sky direction for FAQ pixel
                    GVector vector = rot * m_faq_native[inx];
                    GSkyDir dir(vector);
 
                    // Compute offset angle to centre of elliptical model
                    double theta = elliptical->dir().dist(dir);
 
                    // Continue only if offset angle value is valid
                    if (theta < elliptical->theta_max()) {
 
                        // Compute position angle of elliptical model in the
                        // coordinate system of the elliptical model
                        double posang = elliptical->dir().posang(dir, elliptical->coordsys());
 
                        // Compute elliptical model value
                        double model = elliptical->eval(theta, posang, bin->energy(), bin->time());
 
                        // Continue only if model value is positive
                        if (model > 0.0) {
 
                            // Multiply model value by solid angle of FAQ pixel
                            model *= m_faq_solidangle[inx];
 
                            // Continue only if sky direction is within DRX
                            if (drx.contains(dir)) {
 
                                // Multiply model value by DRX value and IAQ norm
                                model *= drx(dir) * norm;
 
                                // Loop over all phibar layers
                                for (int iphibar = 0, idri = ipix; iphibar < nphibar; ++iphibar, idri+=npix) {
                                    double faq = m_faq(inx, iphibar);
                                    if (faq > 0.0) {
                                        irfs[idri] += model * faq * drg[idri];
                                    }
                                }
 
                            } // endif: sky direction within DRX
 
                        } // endif: model value was positive
 
                    } // endif: offset angle was valid
 
                } // endfor: looped over FAQ pixels
 
            } // endif: elliptical model overlaped with COMPTEL field of view
 
        } // endfor: looped over Chi and Psi bins
 
    } // endif: direct response computation requested
 
    // If weight for numerical integration is non-zero then perform numerical
    // integration and add response according to the weight to the IRF
    // vector
    if (wgt_integration > 0.0) {
 
        // Throw an exception if COMPTEL response is not set or if
        if (m_iaq.empty()) {
            std::string msg = "COMPTEL response is empty. Please initialise the "
                              "response with an \"IAQ\".";
            throw GException::invalid_value(G_IRF_ELLIPTICAL, msg);
        }
        else if (m_phigeo_bin_size == 0.0) {
            std::string msg = "COMPTEL response has a zero Phigeo bin size. "
                              "Please initialise the response with a valid "
                              "\"IAQ\".";
            throw GException::invalid_value(G_IRF_ELLIPTICAL, msg);
        }
 
        // Compute weigthed IAQ normalisation
        double norm = iaq_norm * wgt_integration;
 
        // Get maximum model radius, reduced by a small amount to avoid
        // rounding problems at the boundary of a sharp edged model
        double theta_max = elliptical->theta_max() - 1.0e-12;
 
        // Setup rotation matrix for conversion towards sky coordinates
        GMatrix ry;
        GMatrix rz;
        ry.eulery(elliptical->dir().dec_deg() - 90.0);
        rz.eulerz(-elliptical->dir().ra_deg());
        GMatrix rot = (ry * rz).transpose();
 
        // Loop over Chi and Psi
        for (int ipix = 0; ipix < npix; ++ipix) {
 
            // Get pointer to event bin
            const GCOMEventBin* bin = (*cube)[ipix];
 
            // Get reference to instrument direction
            const GCOMInstDir& obsDir = bin->dir();
 
            // Get reference to Phigeo sky direction
            const GSkyDir& phigeoDir = obsDir.dir();
 
            // Compute angle between model centre and Phigeo sky direction (radians)
            // and position angle of model with respect to Phigeo sky direction (radians)
            double zeta = phigeoDir.dist(elliptical->dir());
 
            // Set radial model zenith angle range
            double rho_min = (zeta > phigeo_max) ? zeta - phigeo_max : 0.0;
            double rho_max = zeta + phigeo_max;
            if (rho_min > theta_max) {
                rho_min = theta_max;
            }
            if (rho_max > theta_max) {
                rho_max = theta_max;
            }
 
            // Perform model zenith angle integration if interval is valid
            if (rho_max > rho_min) {
 
                // Initialise IRF vector
                GVector irf(nphibar);
 
                // Setup integration kernel
                com_elliptical_kerns_rho integrands(m_iaq,
                                                    model,
                                                    irf,
                                                    bin,
                                                    rot,
                                                    drx,
                                                    phigeo_bin,
                                                    m_phigeo_bins,
                                                    nphibar,
                                                    iter_omega);
 
                // Setup integrator
                GIntegrals integral(&integrands);
                integral.fixed_iter(iter_rho);
 
                // Integrate over Phigeo
                irf = integral.romberg(rho_min, rho_max, iter_rho);
 
                // Add IRF to result
                for (int iphibar = 0, idri = ipix; iphibar < nphibar; ++iphibar, idri+=npix) {
                    irfs[idri] += irf[iphibar] * drg[idri] * norm;
                }
 
            } // endif: Phigeo angle interval was valid
 
        } // endfor: looped over Chi and Psi pixels
 
    } // endif: numerical integration computation requested
 
    // Dump CPU consumption
    #if defined(G_DEBUG_IRF_TIMING)
    double celapse = gammalib::get_current_clock() - cstart;
    std::cout << "GCOMResponse::irf_elliptical consumed " << celapse;
    std::cout  << " seconds of CPU time." << std::endl;
    #endif
 
    // Return IRF vector
    return irfs;
}
 
 
/***********************************************************************//**
 * @brief Return instrument response to diffuse source
 *
 * @param[in] model Sky model.
 * @param[in] obs Observation.
 * @param[out] gradients Gradients matrix.
 * @return Instrument response to diffuse source for all events in
 *         observation (\f$cm^2\f$).
 *
 * @exception GException::invalid_value
 *            Response not initialised with a valid IAQ
 *            Incompatible IAQ encountered
 *
 * Returns the instrument response to a diffuse source sky model for all
 * events. The methods works on response vectors that are computed per
 * pointing. The method assumes that the spatial model component is of
 * type GModelSpatialDiffuse, that the observations is of type
 * GCOMObservation and that the events are of type GCOMEventCube.
 * No checking of these assumptions is performed by the method, and not
 * respecting the assumptions will results in a segfault.
 *
 * The computation is done using a direct convolution in the system of the
 * instrument response function, sampling fully the resolution of the data
 * space.
 *
 * For details about the choice of the convolution method and its parameters
 * see https://gitlab.in2p3.fr/gammalib/gammalib/-/issues/14.
 *
 * @p gradients is an optional matrix where the number of rows corresponds
 * to the number of events in the observation and the number of columns
 * corresponds to the number of spatial model parameters. Since for diffuse
 * sources no gradients are computed, the method does not alter the
 * content of @p gradients.
 ***************************************************************************/
GVector GCOMResponse::irf_diffuse(const GModelSky&    model,
                                  const GObservation& obs,
                                  GMatrix*            gradients) const
{
    // Store entry time
    #if defined(G_DEBUG_IRF_TIMING)
    double cstart = gammalib::get_current_clock();
    #endif
 
    // Get pointers
    const GModelSpatialDiffuse* diffuse = static_cast<const GModelSpatialDiffuse*>(model.spatial());
    const GCOMObservation*      obs_ptr = static_cast<const GCOMObservation*>(&obs);
    const GCOMEventCube*        cube    = static_cast<const GCOMEventCube*>(obs.events());
 
    // Get number of Chi/Psi pixels, Phibar layers and event bins
    int npix    = cube->naxis(0) * cube->naxis(1);
    int nphibar = cube->naxis(2);
    int nevents = cube->size();
 
    // Throw an exception if the number of Phibar bins does not match the
    // response
    if (nphibar != m_phibar_bins) {
        std::string msg = "DRE has "+gammalib::str(nphibar)+" Phibar layers "
                          "but IAQ has "+gammalib::str(m_phibar_bins)+" Phibar "
                          "bins. Please specify a compatible IAQ.";
        throw GException::invalid_value(G_IRF_DIFFUSE, msg);
    }
 
    // Initialise result
    GVector irfs(nevents);
 
    // Initialise some variables
    double         phigeo_bin = m_phigeo_bin_size * gammalib::deg2rad;
    double         phigeo_min = m_phigeo_min      * gammalib::deg2rad - 0.5 * phigeo_bin;
    double         phigeo_max = phigeo_min + phigeo_bin * m_phigeo_bins;
    const GSkyMap& drx        = obs_ptr->drx().map();
    const double*  drg        = obs_ptr->drg().map().pixels();
 
    // Compute IAQ normalisation (1/s): DEADC / ONTIME (s)
    double iaq_norm = obs_ptr->deadc() /
                      (obs_ptr->ontime() * cube->dre().tof_correction()) *
                      cube->dre().phase_correction();
 
    #if defined(G_IRF_DIFFUSE_DIRECT)
    #if defined(G_DEBUG_IRF_TIMING)
    std::cout << "GCOMResponse::irf_diffuse - direct computation" << std::endl;
    #endif
 
    // Throw an exception if COMPTEL FAQ response is not set
    if (m_faq.is_empty()) {
        std::string msg = "COMPTEL response is empty. Please initialise the "
                          "response with an \"IAQ\".";
        throw GException::invalid_value(G_IRF_DIFFUSE, msg);
    }
 
    // Loop over Chi and Psi bins
    for (int ipix = 0; ipix < npix; ++ipix) {
 
        // Get pointer to event bin
        const GCOMEventBin* bin = (*cube)[ipix];
 
        // Get reference to instrument direction
        const GCOMInstDir& obsDir = bin->dir();
 
        // Get reference to scatter sky direction
        const GSkyDir& scatter = obsDir.dir();
 
        // transform from FAQ system to the celestial coordinate system
        GMatrix ry;
        GMatrix rz;
        ry.eulery(scatter.dec_deg() - 90.0);
        rz.eulerz(-scatter.ra_deg());
        GMatrix rot = (ry * rz).transpose();
 
        // Loop over FAQ pixels
        for (int inx = 0; inx < m_faq_bins; ++inx) {
 
            // Skip zero FAQ pixels
            if (m_faq_zero[inx]) {
                continue;
            }
 
            // Get sky direction for FAQ pixel
            GVector vector = rot * m_faq_native[inx];
            GSkyDir dir(vector);
 
            // Setup photon
            GPhoton photon(dir, bin->energy(), bin->time());
 
            // Compute model value
            double model = diffuse->eval(photon);
 
            // Continue only if model value is positive
            if (model > 0.0) {
 
                // Multiply model value by solid angle of FAQ pixel
                model *= m_faq_solidangle[inx];
 
                // Continue only if sky direction is within DRX
                if (drx.contains(dir)) {
 
                    // Multiply model value by DRX value and IAQ norm
                    model *= drx(dir) * iaq_norm;
 
                    // Loop over all phibar layers
                    for (int iphibar = 0, idri = ipix; iphibar < nphibar; ++iphibar, idri+=npix) {
                        double faq = m_faq(inx, iphibar);
                        if (faq > 0.0) {
                            irfs[idri] += model * faq * drg[idri];
                        }
                    }
 
                } // endif: sky direction within DRX
 
            } // endif: model value was positive
 
        } // endfor: looped over FAQ pixels
 
    } // endfor: looped over Chi and Psi bins
    #else
    #if defined(G_DEBUG_IRF_TIMING)
    std::cout << "GCOMResponse::irf_diffuse - numerical integration" << std::endl;
    #endif
 
    // Throw an exception if COMPTEL response is not set or if
    if (m_iaq.empty()) {
        std::string msg = "COMPTEL response is empty. Please initialise the "
                          "response with an \"IAQ\".";
        throw GException::invalid_value(G_IRF_DIFFUSE, msg);
    }
    else if (m_phigeo_bin_size == 0.0) {
        std::string msg = "COMPTEL response has a zero Phigeo bin size. "
                          "Please initialise the response with a valid "
                          "\"IAQ\".";
        throw GException::invalid_value(G_IRF_DIFFUSE, msg);
    }
 
    // Initilise some variables
    double omega0 = 2.0 * std::sin(0.5 * phigeo_bin);
 
    // Loop over Chi and Psi
    for (int ipix = 0; ipix < npix; ++ipix) {
 
        // Get pointer to event bin
        const GCOMEventBin* bin = (*cube)[ipix];
 
        // Get reference to instrument direction
        const GCOMInstDir& obsDir = bin->dir();
 
        // Get reference to Phigeo sky direction
        const GSkyDir& phigeoDir = obsDir.dir();
 
        // Loop over Phigeo
        for (int iphigeo = 0; iphigeo < m_phigeo_bins; ++iphigeo) {
 
            // Determine Phigeo value in radians. Note that phigeo_min is the value
            // for bin zero.
            double phigeo     = phigeo_min + double(iphigeo) * phigeo_bin;
            double sin_phigeo = std::sin(phigeo);
 
            // Determine number of azimuthal integration steps and step size
            // in radians
            double length = gammalib::twopi * sin_phigeo;
            int    naz    = int(length / phigeo_bin + 0.5);
            if (naz < 2) {
                naz = 2;
            }
            double az_step = gammalib::twopi / double(naz);
 
            // Computes solid angle of integration bin multiplied by Jaccobian
            double omega = omega0 * az_step * sin_phigeo;
 
            // Loop over azimuth angle
            double az = 0.0;
            for (int iaz = 0; iaz < naz; ++iaz, az += az_step) {
 
                // Get sky direction
                GSkyDir skyDir = phigeoDir;
                skyDir.rotate(az, phigeo);
 
                // Fall through if sky direction is not contained in DRX
                if (!drx.contains(skyDir)) {
                    continue;
                }
 
                // Set photon
                GPhoton photon(skyDir, bin->energy(), bin->time());
 
                // Get model sky intensity for photon (unit: sr^-1)
                double intensity = model.spatial()->eval(photon);
 
                // Fall through if intensity is zero
                if (intensity == 0.0) {
                    continue;
                }
 
                // Multiply intensity by DRX value (unit: cm^2 s sr^-1)
                intensity *= drx(skyDir);
 
                // Multiply intensity by solid angle (unit: cm^2 s)
                intensity *= omega;
 
                // Loop over Phibar
                for (int iphibar = 0; iphibar < nphibar; ++iphibar) {
 
                    // Get IAQ index
                    int i = iphibar * m_phigeo_bins + iphigeo;
 
                    // Get IAQ value
                    double iaq = m_iaq[i];
 
                    // Fall through if IAQ is not positive
                    if (iaq <= 0.0) {
                        continue;
                    }
 
                    // Get DRI index
                    int idri = ipix + iphibar * npix;
 
                    // Compute IRF value (unit: cm^2)
                    double irf = iaq * drg[idri] * iaq_norm * intensity;
 
                    // Add IRF value if it is positive
                    if (irf > 0.0) {
                        irfs[idri] += irf;
                    }
 
                } // endfor: looped over Phibar
 
            } // endfor: looped over azimuth angle around locus of IAQ
 
        } // endfor: looped over Phigeo angles
 
    } // endfor: looped over Chi and Psi pixels
    #endif
 
    // Dump CPU consumption
    #if defined(G_DEBUG_IRF_TIMING)
    double celapse = gammalib::get_current_clock() - cstart;
    std::cout << "GCOMResponse::irf_diffuse consumed " << celapse;
    std::cout  << " seconds of CPU time." << std::endl;
    #endif
 
    // Return IRF vector
    return irfs;
}