spi_helpers_response_vector.cpp

Go to the documentation of this file.

/***************************************************************************
 *                  SPI helper classes for vector response                 *
 * ----------------------------------------------------------------------- *
 *  copyright (C) 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 spi_helpers_response_vector.cpp
 * @brief Implementation of SPI helper classes for vector response
 * @author Juergen Knoedlseder
 */
 
/* __ Includes ___________________________________________________________ */
#ifdef HAVE_CONFIG_H
#include <config.h>
#endif
#include <cmath>
#include "GVector.hpp"
#include "GMatrix.hpp"
#include "GIntegrals.hpp"
#include "GModelSpatialRadial.hpp"
#include "GModelSpatialElliptical.hpp"
#include "GSPIResponse.hpp"
#include "spi_helpers_response_vector.hpp"
 
/* __ Method name definitions ____________________________________________ */
 
/* __ Macros _____________________________________________________________ */
 
/* __ Coding definitions _________________________________________________ */
 
/* __ Debug definitions __________________________________________________ */
 
/* __ Constants __________________________________________________________ */
 
 
/*==========================================================================
 =                                                                         =
 =                             Helper methods                              =
 =                                                                         =
 ==========================================================================*/
 
/***********************************************************************//**
 * @brief Kernel for radial integration of radial models
 *
 * @param[in] rho Rho angle (radians).
 * @return Azimuthally integrated radial model.
 ***************************************************************************/
GVector spi_radial_kerns_rho::eval(const double& rho)
{
    // Set static model energy and time
    static const GEnergy energy;
    static const GTime   time;
 
    // Initialise kernel vector values to zero
    m_irfs = 0.0;
 
    // Continue only if rho is positive
    if (rho > 0.0) {
 
        // Compute radial model value
        double model = m_radial->eval(rho, energy, time);
 
        // Continue only if model is positive
        if (model > 0.0) {
 
            // Precompute cosine and sine terms for azimuthal integration
            double sin_rho = std::sin(rho);
            double cos_rho = std::cos(rho);
 
            // Setup azimuthal integration kernel
            spi_radial_kerns_omega integrands(m_cube,
                                              m_rsp,
                                              m_ipt,
                                              m_livetimes,
                                              m_rot,
                                              sin_rho,
                                              cos_rho,
                                              m_irfs);
 
            // Setup integrator
            GIntegrals integral(&integrands);
            integral.fixed_iter(m_iter);
 
            // Integrate over Omega angle
            m_irfs = integral.romberg(0.0, gammalib::twopi, m_iter) * model * sin_rho;
 
        } // endif: radial model was positive
 
    } // endif: phigeo was positive
 
    // Return kernel values
    return m_irfs;
}
 
 
/***********************************************************************//**
 * @brief Kernel for azimuthal integration of radial models
 *
 * @param[in] omega Omega angle (radians).
 * @return Kernel value for radial model.
 ***************************************************************************/
GVector spi_radial_kerns_omega::eval(const double& omega)
{
    // Initialise kernel vector values to zero
    m_irfs = 0.0;
 
    // Set up native coordinate vector
    double  sin_omega = std::sin(omega);
    double  cos_omega = std::cos(omega);
    GVector native(-cos_omega*m_sin_rho, sin_omega*m_sin_rho, m_cos_rho);
 
    // Rotate vector into celestial coordinates
    GVector vector = m_rot * native;
 
    // Convert vector into sky direction
    GSkyDir dir(vector);
 
    // Convert sky direction to zenith angle (radians)
    double zenith = m_rsp->zenith(m_ipt, dir);
 
    // Continue only if zenith angle is smaller than the maximum zenith angle.
    if (zenith < m_rsp->m_max_zenith) {
 
        // Convert sky direction to azimuth angle (radians)
        double azimuth = m_rsp->azimuth(m_ipt, dir);
 
        // Compute pixel
        double xpix = (zenith * std::cos(azimuth) - m_rsp->m_wcs_xmin) / m_rsp->m_wcs_xbin;
        double ypix = (zenith * std::sin(azimuth) - m_rsp->m_wcs_ymin) / m_rsp->m_wcs_ybin;
 
        // Continue only if pixel is within IRF
        if (xpix > 0.0 && xpix < m_rsp->m_wcs_xpix_max &&
            ypix > 0.0 && ypix < m_rsp->m_wcs_ypix_max) {
 
            // Compute IRF values for this pointing (photo peak only)
            m_rsp->irf_vector(m_cube, xpix, ypix, 0, m_livetimes, &m_irfs);
 
        } // endif: pixel was within IRF
 
    } // endif: zenith angle was within validity range
 
    // Return kernel values
    return m_irfs;
}
 
 
/***********************************************************************//**
 * @brief Kernel for radial integration of elliptical models
 *
 * @param[in] rho Rho angle (radians).
 * @return Azimuthally integrated elliptical model.
 ***************************************************************************/
GVector spi_elliptical_kerns_rho::eval(const double& rho)
{
    // Initialise kernel vector values to zero
    m_irfs = 0.0;
 
    // Continue only if rho is positive
    if (rho > 0.0) {
 
        // Precompute cosine and sine terms for azimuthal integration
        double sin_rho = std::sin(rho);
        double cos_rho = std::cos(rho);
 
        // Setup azimuthal integration kernel
        spi_elliptical_kerns_omega integrands(m_cube,
                                              m_rsp,
                                              m_elliptical,
                                              m_ipt,
                                              m_livetimes,
                                              m_rot,
                                              rho,
                                              sin_rho,
                                              cos_rho,
                                              m_irfs);
 
        // Setup integrator
        GIntegrals integral(&integrands);
        integral.fixed_iter(m_iter);
 
        // Integrate over Omega angle
        m_irfs = integral.romberg(0.0, gammalib::twopi, m_iter) * sin_rho;
 
    } // endif: phigeo was positive
 
    // Return kernel values
    return m_irfs;
}
 
 
/***********************************************************************//**
 * @brief Kernel for azimuthal integration of elliptical models
 *
 * @param[in] omega Omega angle (radians).
 * @return Kernel value for elliptical model.
 ***************************************************************************/
GVector spi_elliptical_kerns_omega::eval(const double& omega)
{
    // Set static model energy and time
    static const GEnergy energy;
    static const GTime   time;
 
    // Initialise kernel vector values to zero
    m_irfs = 0.0;
 
    // Set up native coordinate vector
    double  sin_omega = std::sin(omega);
    double  cos_omega = std::cos(omega);
    GVector native(-cos_omega*m_sin_rho, sin_omega*m_sin_rho, m_cos_rho);
 
    // Rotate vector into celestial coordinates
    GVector vector = m_rot * native;
 
    // Convert vector into sky direction
    GSkyDir dir(vector);
 
    // Convert sky direction to zenith angle (radians)
    double zenith = m_rsp->zenith(m_ipt, dir);
 
    // Continue only if zenith angle is smaller than the maximum zenith angle.
    if (zenith < m_rsp->m_max_zenith) {
 
        // Convert sky direction to azimuth angle (radians)
        double azimuth = m_rsp->azimuth(m_ipt, dir);
 
        // Compute pixel
        double xpix = (zenith * std::cos(azimuth) - m_rsp->m_wcs_xmin) / m_rsp->m_wcs_xbin;
        double ypix = (zenith * std::sin(azimuth) - m_rsp->m_wcs_ymin) / m_rsp->m_wcs_ybin;
 
        // Continue only if pixel is within IRF
        if (xpix > 0.0 && xpix < m_rsp->m_wcs_xpix_max &&
            ypix > 0.0 && ypix < m_rsp->m_wcs_ypix_max) {
 
            // Compute radial model value
            double model = m_elliptical->eval(m_rho, omega, energy, time);
 
            // Continue only if model is positive
            if (model > 0.0) {
 
                // Compute IRF values for this pointing (photo peak only)
                m_rsp->irf_vector(m_cube, xpix, ypix, 0, m_livetimes, &m_irfs);
 
                // Multiply with model value
                m_irfs *= model;
 
            } // endif: model was positive
 
        } // endif: pixel was within IRF
 
    } // endif: zenith angle was within validity range
 
    // Return kernel values
    return m_irfs;
}