GSPIResponse.hpp

Go to the documentation of this file.

/***************************************************************************
 *              GSPIResponse.hpp - INTEGRAL/SPI response class             *
 * ----------------------------------------------------------------------- *
 *  copyright (C) 2020-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 GSPIResponse.hpp
 * @brief INTEGRAL/SPI instrument response function class definition
 * @author Juergen Knoedlseder
 */
 
#ifndef GSPIRESPONSE_HPP
#define GSPIRESPONSE_HPP
 
/* __ Includes ___________________________________________________________ */
#include <vector>
#include "GNodeArray.hpp"
#include "GFilename.hpp"
#include "GVector.hpp"
#include "GMatrix.hpp"
#include "GEbounds.hpp"
#include "GSkyDir.hpp"
#include "GSkyMap.hpp"
#include "GResponse.hpp"
 
/* __ Forward declaration ________________________________________________ */
class GEvent;
class GPhoton;
class GEnergy;
class GTime;
class GSource;
class GObservation;
class GModelSky;
class GModelSpatialDiffuseMap;
class GEbounds;
class GSPIObservation;
class GSPIEventCube;
class GSPIEventBin;
 
/* __ Type definitions ___________________________________________________ */
typedef struct {
    int               nevents;     // Number of events
    int               nirf;        // Dimension of IRF
    int               nx;          // Number of IRF pixels in x direction
    int               ny;          // Number of IRF pixels in y direction
    int               nxy;         // Number of IRF pixels
    int               ndet;        // Number of detectors
    int               nreg;        // Number of IRF regions
    int               ireg;        // IRF region (0 = photopeak)
    std::vector<int>  detids;      // Vector of detector IDs
    std::vector<bool> zero;        // Vector of IRF pixel usage
    GVector           zeniths;     // Zenith angles for all IRF pixels
    GVector           sin_zeniths; // Sine of zenith angles for all IRF pixels
    GVector           cos_zeniths; // Cosine of zenith angles for all IRF pixels
    GVector           azimuths;    // Azimuth angles for all IRF pixels
    GVector           solidangle;  // Solid angles for all IRF pixels
    GMatrix           rot;         // Current rotation matrix
    GSkyDir           dir;         // Current sky direction of IRF pixel
} GSPIResponseIrf;
 
/* __ Constants __________________________________________________________ */
 
 
/***********************************************************************//**
 * @class GSPIResponse
 *
 * @brief INTEGRAL/SPI instrument response function class
 *
 * The INTEGRAL/SPI instrument response function class defines the function
 * that translates from physical quantities to measured quantities.
 *
 * @todo Complete the class description.
 ***************************************************************************/
class GSPIResponse : public GResponse {
 
    // Friends
    friend class spi_radial_kerns_rho;
    friend class spi_radial_kerns_omega;
    friend class spi_elliptical_kerns_rho;
    friend class spi_elliptical_kerns_omega;
 
public:
    // Constructors and destructors
    GSPIResponse(void);
    GSPIResponse(const GSPIResponse& rsp);
    explicit GSPIResponse(const GFilename& rspname);
    virtual ~GSPIResponse(void);
 
    // Operators
    virtual GSPIResponse& operator=(const GSPIResponse & rsp);
 
    // Implement pure virtual base class methods
    virtual void          clear(void);
    virtual GSPIResponse* clone(void) const;
    virtual std::string   classname(void) const;
    virtual bool          use_edisp(void) const;
    virtual bool          use_tdisp(void) const;
    virtual double        irf(const GEvent&       event,
                              const GPhoton&      photon,
                              const GObservation& obs) const;
    virtual double        nroi(const GModelSky&    model,
                               const GEnergy&      obsEng,
                               const GTime&        obsTime,
                               const GObservation& obs) const;
    virtual GEbounds      ebounds(const GEnergy& obsEnergy) const;
    virtual std::string   print(const GChatter& chatter = NORMAL) const;
 
    // Other Methods
    void             rspname(const GFilename& rspname);
    const GFilename& rspname(void) const;
    bool             is_precomputed(void) const;
    const double&    energy_keV(void) const;
    const double&    dlogE(void) const;
    const double&    gamma(void) const;
    void             set(const GSPIObservation& obs,
                         const GEnergy&         energy = GEnergy());
    double           irf_value(const GSkyDir&      srcDir,
                               const GSPIEventBin& bin,
                               const int&          ireg) const;
    double           zenith(const int& ipt, const GSkyDir& dir) const;
    double           azimuth(const int& ipt, const GSkyDir& dir) const;
    void             read(const GFits& fits);
    void             write(GFits& fits) const;
    void             load(const GFilename& filename);
    void             save(const GFilename& filename,
                          const bool&      clobber = false) const;
 
private:
    // Private methods
    void    init_members(void);
    void    copy_members(const GSPIResponse& rsp);
    void    free_members(void);
    void    read_detids(const GFits& fits);
    void    read_energies(const GFits& fits);
    void    write_detids(GFits& fits) const;
    void    write_energies(GFits& fits) const;
    void    load_irfs(const int& region);
    GSkyMap load_irf(const GFilename& rspname) const;
    GSkyMap compute_irf(const double& emin, const double& emax) const;
    void    set_wcs(const GFitsImage* image);
    void    set_detids(const GSPIEventCube* cube);
    void    set_cache(const GSPIEventCube* cube);
    void    irf_vector(const GSPIEventCube* cube,
                       const double&        xpix,
                       const double&        ypix,
                       const int&           ireg,
                       const double*        livetimes,
                       GVector*             irf) const;
    int     irf_detid(const int& detid) const;
    double  irf_weight(const double& beta,
                       const double& emin,
                       const double& emax) const;
    void    irf_diffuse_map(const GModelSpatialDiffuseMap* diffuse,
                            const GObservation&            obs,
                            GVector*                       irfs,
                            GMatrix*                       gradients = NULL) const;
    void    irf_init(const GObservation& obs, GSPIResponseIrf* irf) const;
    void    irf_rot_matrix(const int&       ipt,
                           GSPIResponseIrf* irf) const;
    void    irf_skydir(const int&       ipt,
                       const int&       inx,
                       GSPIResponseIrf* irf) const;
    int     irf_map(const int&             idet,
                    const int&             ieng,
                    const GSPIResponseIrf* irf) const;
 
    // Overloaded virtual base class methods
    virtual GVector irf_ptsrc(const GModelSky&    model,
                              const GObservation& obs,
                              GMatrix*            gradients = NULL) const;
    virtual GVector irf_radial(const GModelSky&    model,
                               const GObservation& obs,
                               GMatrix*            gradients = NULL) const;
    virtual GVector irf_elliptical(const GModelSky&    model,
                                   const GObservation& obs,
                                   GMatrix*            gradients = NULL) const;
    virtual GVector irf_diffuse(const GModelSky&    model,
                                const GObservation& obs,
                                GMatrix*            gradients = NULL) const;
 
    // Private data members
    GFilename            m_rspname;    //!< File name of response group
    std::vector<int>     m_detids;     //!< Vector of detector IDs
    GNodeArray           m_energies;   //!< Node array of IRF energies
    GEbounds             m_ebounds;    //!< Energy bounaries of IRF
    GSkyMap              m_irfs;       //!< IRFs stored as sky map
    double               m_energy_keV; //!< IRF line energy (optional)
    double               m_dlogE;      //!< Logarithmic energy step for IRF band
    double               m_gamma;      //!< Power-law spectral index for IRF band
 
    // Private cache
    std::vector<GSkyDir> m_spix;         //!< SPI pointing direction
    std::vector<double>  m_posang;       //!< Position angle of Y axis (CEL, radians)
    mutable bool         m_has_wcs;      //!< Has WCS information
    mutable double       m_wcs_xmin;     //!< Minimum X value (radians)
    mutable double       m_wcs_ymin;     //!< Minimum Y value (radians)
    mutable double       m_wcs_xmax;     //!< Maximum X value (radians)
    mutable double       m_wcs_ymax;     //!< Maximum Y value (radians)
    mutable double       m_wcs_xbin;     //!< X value bin size (radians)
    mutable double       m_wcs_ybin;     //!< Y value bin size (radians)
    mutable double       m_wcs_xpix_max; //!< Maximum X pixel index
    mutable double       m_wcs_ypix_max; //!< Maximum Y pixel index
    mutable double       m_max_zenith;   //!< Maximum zenith angle (radians)
};
 
 
/***********************************************************************//**
 * @brief Return class name
 *
 * @return String containing the class name ("GSPIResponse").
 ***************************************************************************/
inline
std::string GSPIResponse::classname(void) const
{
    return ("GSPIResponse");
}
 
 
/***********************************************************************//**
 * @brief Signal if energy dispersion will be used
 *
 * @return False.
 *
 * @todo Implement method as needed.
 ***************************************************************************/
inline
bool GSPIResponse::use_edisp(void) const
{
    return false;
}
 
 
/***********************************************************************//**
 * @brief Signal if time dispersion will be used
 *
 * @return False.
 *
 * @todo Implement method as needed.
 ***************************************************************************/
inline
bool GSPIResponse::use_tdisp(void) const
{
    return false;
}
 
 
/***********************************************************************//**
 * @brief Set response name
 *
 * @param[in] rspname Response group file name.
 *
 * Sets the response group file name.
 ***************************************************************************/
inline
void GSPIResponse::rspname(const GFilename& rspname)
{
    m_rspname = rspname;
    return;
}
 
 
/***********************************************************************//**
 * @brief Get response group file name
 *
 * @return Response group file name.
 *
 * Returns the response group file name.
 ***************************************************************************/
inline
const GFilename& GSPIResponse::rspname(void) const
{
    return m_rspname;
}
 
 
/***********************************************************************//**
 * @brief Signals if response is precomputed
 *
 * @return True if response is precomputed.
 *
 * Signals if the response was precomputed.
 ***************************************************************************/
inline
bool GSPIResponse::is_precomputed(void) const
{
    return (!m_ebounds.is_empty());
}
 
 
/***********************************************************************//**
 * @brief Return line IRF energy in keV
 *
 * @return Line IRF energy (keV).
 *
 * Returns the energy in keV for a line IRF. If the IRF is a continuum IRF
 * the method returns 0.
 ***************************************************************************/
inline
const double& GSPIResponse::energy_keV(void) const
{
    return (m_energy_keV);
}
 
 
/***********************************************************************//**
 * @brief Return logarithmic step size for continuum IRFs
 *
 * @return Logarithmic step size for continuum IRFs.
 *
 * Returns the logarithmic step size for the computation of continuum IRFs.
 ***************************************************************************/
inline
const double& GSPIResponse::dlogE(void) const
{
    return (m_dlogE);
}
 
 
/***********************************************************************//**
 * @brief Return power-law index for continuum IRFs
 *
 * @return Power-law index for continuum IRFs.
 *
 * Returns the power-law index for the computation of continuum IRFs.
 ***************************************************************************/
inline
const double& GSPIResponse::gamma(void) const
{
    return (m_gamma);
}
 
 
/***********************************************************************//**
 * @brief Return zenith angle of sky direction for pointing in radians
 *
 * @param[in] ipt Pointing index.
 * @param[in] dir Sky direction.
 * @return Zenith angle (radians).
 *
 * Returns zenith angle of sky direction for pointing in radians.
 ***************************************************************************/
inline
double GSPIResponse::zenith(const int& ipt, const GSkyDir& dir) const
{
    return (m_spix[ipt].dist(dir));
}
 
 
/***********************************************************************//**
 * @brief Return azimuth angle of sky direction for pointing in radians
 *
 * @param[in] ipt Pointing index.
 * @param[in] dir Sky direction.
 * @return Azimuth angle (radians).
 *
 * Returns azimuth angle of sky direction for pointing in radians.
 ***************************************************************************/
inline
double GSPIResponse::azimuth(const int& ipt, const GSkyDir& dir) const
{
    double azimuth = m_posang[ipt] - m_spix[ipt].posang(dir); // Celestial system
    if (azimuth < 0.0) {
        azimuth += gammalib::twopi;
    }
    return (azimuth);
}
 
 
/***********************************************************************//**
 * @brief Return map index for detector and energy index
 *
 * @param[in] idet Detector index.
 * @param[in] ieng Energy index.
 * @param[in] irf Stucture for projection computation of response.
 * @return Map index.
 *
 * Returns the map index for a specific detector index @p idet and energy
 * index @p ieng.
 *
 * The client needs to make sure that the indices are comprised within the
 * valid ranges and that the GSPIResponse::irf_init() method was called
 * prior to invoking this method.
 ***************************************************************************/
inline
int GSPIResponse::irf_map(const int&             idet,
                          const int&             ieng,
                          const GSPIResponseIrf* irf) const
{
    // Get IRF map index
    int map = irf->detids[idet] + (irf->ireg + ieng * irf->nreg) * irf->ndet;
 
    // Return map index
    return map;
}
 
#endif /* GSPIRESPONSE_HPP */