From c238a196d2e5259c529643eee2963543d33d169e Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Wed, 29 Oct 2014 16:12:49 +0100 Subject: [PATCH] Correct transition point for shell model (was in degrees but the method expected radians; this led to convergence problems with the shell model). --- inst/cta/include/GCTAResponseCube.hpp | 2 +- inst/cta/src/GCTAResponseCube.cpp | 131 +++++++++++++++++++++-------------------------------------------------------------------------------------------------------------- inst/cta/src/GCTAResponseIrf.cpp | 183 ++------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- inst/cta/src/GCTAResponse_helpers.cpp | 569 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ inst/cta/src/GCTAResponse_helpers.hpp | 357 ++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- 5 files changed, 145 insertions(+), 1097 deletions(-) diff --git a/inst/cta/include/GCTAResponseCube.hpp b/inst/cta/include/GCTAResponseCube.hpp index 1736f81..c8b5a89 100644 --- a/inst/cta/include/GCTAResponseCube.hpp +++ b/inst/cta/include/GCTAResponseCube.hpp @@ -111,7 +111,7 @@ private: void free_members(void); int cache_index(const std::string& name) const; double psf_radial(const GModelSpatialRadial* model, - const double& zeta, + const double& rho_obs, const GSkyDir& obsDir, const GEnergy& srcEng, const GTime& srcTime) const; diff --git a/inst/cta/src/GCTAResponseCube.cpp b/inst/cta/src/GCTAResponseCube.cpp index 61db959..a97a7ea 100644 --- a/inst/cta/src/GCTAResponseCube.cpp +++ b/inst/cta/src/GCTAResponseCube.cpp @@ -37,6 +37,7 @@ #include "GCTASourceCubeDiffuse.hpp" #include "GModelSpatialPointSource.hpp" #include "GModelSpatialRadial.hpp" +#include "GModelSpatialRadialShell.hpp" #include "GModelSpatialElliptical.hpp" #include "GModelSpatialDiffuse.hpp" #include "GPhoton.hpp" @@ -63,8 +64,6 @@ /* __ Macros _____________________________________________________________ */ /* __ Coding definitions _________________________________________________ */ -#define G_NEW_PSF_RADIAL -#define G_NEW_PSF_ELLIPTICAL /* __ Debug definitions __________________________________________________ */ @@ -464,11 +463,11 @@ double GCTAResponseCube::irf_radial(const GEvent& event, // Compute angle between model centre and measured photon direction // (radians) - double zeta = model->dir().dist(obsDir); + double rho_obs = model->dir().dist(obsDir); // Continue only if we're sufficiently close to the model centre to get a // non-zero response - if (zeta <= model->theta_max()+psf().delta_max()) { + if (rho_obs <= model->theta_max()+psf().delta_max()) { // Get exposure irf = exposure()(obsDir, obsEng); @@ -476,7 +475,7 @@ double GCTAResponseCube::irf_radial(const GEvent& event, // Continue only if exposure is positive if (irf > 0.0) { irf /= obs.livetime(); //!< Recover effective area from exposure - irf *= psf_radial(model, zeta, obsDir, obsEng, obsTime); + irf *= psf_radial(model, rho_obs, obsDir, obsEng, obsTime); irf *= obs.deadc(obsTime); } @@ -535,8 +534,8 @@ double GCTAResponseCube::irf_elliptical(const GEvent& event, // Compute angle between model centre and measured photon direction and // position angle (radians) - double rho_obs = model->dir().dist(obsDir); // was zeta - double posangle_obs = model->dir().posang(obsDir); // was omega + double rho_obs = model->dir().dist(obsDir); + double posangle_obs = model->dir().posang(obsDir); // Continue only if we're sufficiently close to the model centre to get a // non-zero response @@ -927,8 +926,8 @@ int GCTAResponseCube::cache_index(const std::string& name) const * @brief Integrate PSF over radial model * * @param[in] model Radial model. - * @param[in] zeta Angle between model centre and measured photon direction - * (radians). + * @param[in] rho_obs Angle between model centre and measured photon direction + * (radians). * @param[in] obsDir Observed event direction. * @param[in] srcEng True photon energy. * @param[in] srcTime True photon arrival time. @@ -949,9 +948,8 @@ int GCTAResponseCube::cache_index(const std::string& name) const * \f${\rm PSF}(\delta)\f$ is the azimuthally symmetric point spread * function. ***************************************************************************/ -#if defined(G_NEW_PSF_ELLIPTICAL) double GCTAResponseCube::psf_radial(const GModelSpatialRadial* model, - const double& rho_obs, // was zeta + const double& rho_obs, const GSkyDir& obsDir, const GEnergy& srcEng, const GTime& srcTime) const @@ -968,7 +966,7 @@ double GCTAResponseCube::psf_radial(const GModelSpatialRadial* model, // Get maximum PSF radius (radians) double delta_max = psf().delta_max(); - // Set zenith angle integration range for elliptical model (radians) + // Set zenith angle integration range for radial model (radians) double rho_min = (rho_obs > delta_max) ? rho_obs - delta_max : 0.0; double rho_max = rho_obs + delta_max; double src_max = model->theta_max(); @@ -1016,6 +1014,17 @@ double GCTAResponseCube::psf_radial(const GModelSpatialRadial* model, bounds.push_back(transition_point); } + // If we have a shell model then add an integration boundary for the + // shell radius as a function discontinuity will occur at this + // location + const GModelSpatialRadialShell* shell = dynamic_cast(model); + if (shell != NULL) { + double shell_radius = shell->radius() * gammalib::deg2rad; + if (shell_radius > rho_min && shell_radius < rho_max) { + bounds.push_back(shell_radius); + } + } + // Integrate kernel irf = integral.romberg(bounds, iter_rho); @@ -1036,49 +1045,6 @@ double GCTAResponseCube::psf_radial(const GModelSpatialRadial* model, // Return PSF return irf; } -#else -double GCTAResponseCube::psf_radial(const GModelSpatialRadial* model, - const double& zeta, - const GSkyDir& obsDir, - const GEnergy& srcEng, - const GTime& srcTime) const -{ - // Set integration precision and Romberg order. These values have been - // determined after careful testing, see - // https://cta-redmine.irap.omp.eu/issues/1291 - const double eps_delta = 1.0e-2; - const int order_delta = 5; - const double eps_phi = 1.0e-2; - const int order_phi = 5; - - // Initialise value - double value = 0.0; - - // Compute sin and cosine of angle between model centre and measured - // photon direction - double sin_zeta = std::sin(zeta); - double cos_zeta = std::cos(zeta); - - // Get offset angle integration interval in radians - double delta_min = 0.0; - double delta_max = psf().delta_max(); - - // Setup integration kernel. We take here the observed photon arrival - // direction as the true photon arrival direction because the PSF does - // not vary significantly over a small region. - cta_psf_radial_kern_delta integrand(this, model, obsDir, srcEng, srcTime, - sin_zeta, cos_zeta, - eps_phi, order_phi); - - // Integrate over PSF delta angle - GIntegral integral(&integrand); - integral.eps(eps_delta); - value = integral.romberg(delta_min, delta_max, order_delta); - - // Return PSF - return value; -} -#endif /***********************************************************************//** @@ -1124,7 +1090,6 @@ double GCTAResponseCube::psf_radial(const GModelSpatialRadial* model, * \f${\rm PSF}(\delta)\f$ is the azimuthally symmetric point spread * function. ***************************************************************************/ -#if defined(G_NEW_PSF_ELLIPTICAL) double GCTAResponseCube::psf_elliptical(const GModelSpatialElliptical* model, const double& rho_obs, const double& posangle_obs, @@ -1201,16 +1166,6 @@ double GCTAResponseCube::psf_elliptical(const GModelSpatialElliptical* model, bounds.push_back(rho_min); bounds.push_back(rho_max); - // If the integration range includes a transition between full - // containment and partial containment, then add a boundary at - // this location - /* - double transition_point = delta_max - rho_obs; - if (transition_point > rho_min && transition_point < rho_max) { - bounds.push_back(transition_point); - } - */ - // If the integration range includes the semiminor boundary, then // add an integration boundary at that location if (semiminor > rho_min && semiminor < rho_max) { @@ -1237,47 +1192,3 @@ double GCTAResponseCube::psf_elliptical(const GModelSpatialElliptical* model, // Return PSF return irf; } -#else -double GCTAResponseCube::psf_elliptical(const GModelSpatialElliptical* model, - const double& rho_obs, - const double& posangle_obs, - const GSkyDir& obsDir, - const GEnergy& srcEng, - const GTime& srcTime) const -{ - // Set integration precision and Romberg order. These values have been - // determined after careful testing, see - // https://cta-redmine.irap.omp.eu/issues/1291 - const double eps_delta = 1.0e-2; - const int order_delta = 5; - const double eps_phi = 1.0e-2; - const int order_phi = 5; - - // Initialise value - double value = 0.0; - - // Compute sin and cosine of angle between model centre and measured - // photon direction - double sin_rho_obs = std::sin(rho_obs); - double cos_rho_obs = std::cos(rho_obs); - - // Get offset angle integration interval in radians - double delta_min = 0.0; - double delta_max = psf().delta_max(); - - // Setup integration kernel. We take here the observed photon arrival - // direction as the true photon arrival direction because the PSF does - // not vary significantly over a small region. - cta_psf_elliptical_kern_delta integrand(this, model, obsDir, srcEng, srcTime, - sin_rho_obs, cos_rho_obs, posangle_obs, - eps_phi, order_phi); - - // Integrate over PSF delta angle - GIntegral integral(&integrand); - integral.eps(eps_delta); - value = integral.romberg(delta_min, delta_max, order_delta); - - // Return PSF - return value; -} -#endif diff --git a/inst/cta/src/GCTAResponseIrf.cpp b/inst/cta/src/GCTAResponseIrf.cpp index 22fe4c1..c3f1be2 100644 --- a/inst/cta/src/GCTAResponseIrf.cpp +++ b/inst/cta/src/GCTAResponseIrf.cpp @@ -1282,174 +1282,6 @@ std::string GCTAResponseIrf::print(const GChatter& chatter) const * @exception GCTAException::bad_model_type * Model is not a radial model. * - * Integrates the product of the model and the IRF over the true photon - * arrival direction using - * - * \f[ - * \int_{\rho_{\rm min}}^{\rho_{\rm max}} - * \sin \rho \times S_{\rm p}(\rho | E, t) \times - * \int_{\omega_{\rm min}}^{\omega_{\rm max}} - * IRF(\rho, \omega) d\omega d\rho - * \f] - * - * where - * - \f$S_{\rm p}(\rho | E, t)\f$ is the radial model, - * - \f$IRF(\rho, \omega)\f$ is the IRF - * - \f$\rho\f$ is the distance from the model centre, and - * - \f$\omega\f$ is the azimuth angle is the position angle with respect to - * the connecting line between the model centre and the observed photon - * arrival direction. - * - * The integration is performed in the coordinate system of the source - * model spanned by \f$\rho\f$ and \f$\omega\f$ which allows to benefit - * from the symmetry of the source model. - * - * The source centre is located at \f$\vec{m}\f$, and a spherical system - * is defined around this location with \f$(\omega,\rho)\f$ being the - * azimuth and zenith angles, respectively. \f$\omega=0\f$ is defined - * by the direction that connects the source centre \f$\vec{m}\f$ to the - * measured photon direction \f$\vec{p'}\f$, and \f$\omega\f$ increases - * counterclockwise. - * - * Note that this method approximates the true theta angle (angle between - * incident photon and pointing direction) by the measured theta angle - * (angle between the measured photon arrival direction and the pointing - * direction). Given the slow variation of the PSF shape over the field of - * view, this approximation should be fine. It helps in fact a lot in - * speeding up the computations. - ***************************************************************************/ -#if defined(G_USE_PSF_SYSTEM) -double GCTAResponseIrf::irf_radial(const GEvent& event, - const GSource& source, - const GObservation& obs) const -{ - // Set integration precision and Romberg order. These values have been - // determined after careful testing, see - // https://cta-redmine.irap.omp.eu/issues/1299 - static const int iter_delta = 5; - static const int iter_phi = 5; - - // Set a tiny angle (in radians) to be removed from the model boundary - // to avoid requesting values beyond the boundary due to round-off - // errors (this could produce a discontinuity that the integrator does - // not like) - static const double tiny = 1.0e-9; - - // Initialise Irf - double irf = 0.0; - - // Retrieve CTA pointing - const GCTAPointing& pnt = retrieve_pnt(G_IRF_RADIAL, obs); - const GCTAInstDir& dir = retrieve_dir(G_IRF_RADIAL, event); - - // Get pointer to radial model - const GModelSpatialRadial* model = static_cast(source.model()); - - // Get event attributes - const GSkyDir& obsDir = dir.dir(); - const GEnergy& obsEng = event.energy(); - - // Get source attributes - const GEnergy& srcEng = source.energy(); - const GTime& srcTime = source.time(); - double srcLogEng = srcEng.log10TeV(); - - // Angle between measured photon and model centre (radians) - double zeta = obsDir.dist(model->dir()); - - // Angle between measured photon and pointing (radians) - double obsOffset = obsDir.dist(pnt.dir()); - - // Compute position angle of pointing in Psf system (radians). - // The angle Phi=0 of the Psf system is given by the connecting line - // between the measured photon direction and the model centre - double phi0 = obsDir.posang(pnt.dir()) - obsDir.posang(model->dir()); - - // Get Psf delta integration interval [delta_min, delta_max] (radians) - double theta_max = model->theta_max() - tiny; - double delta_min = (zeta > theta_max) ? zeta - theta_max : 0.0; - double delta_max = psf_delta_max(obsOffset, 0.0, - pnt.zenith(), pnt.azimuth(), - srcLogEng); - double distance = zeta + theta_max; - if (distance < delta_max) { - delta_max = distance; - } - - // Continue only if interval has positive length - if (delta_max > delta_min) { - - // Setup integration kernel. - cta_irf_radial_kern_delta integrand(this, model, pnt, - srcEng, srcLogEng, srcTime, - obsOffset, obsEng, - zeta, phi0, theta_max, - iter_phi); - - // Setup integration over Psf delta angle - GIntegral integral(&integrand); - integral.fixed_iter(iter_delta); - - // Setup integration boundaries - std::vector bounds; - bounds.push_back(delta_min); - bounds.push_back(delta_max); - - // If the integration range includes a transition between full - // containment and partial containment, then add a boundary at - // this location - double transition_point = theta_max - zeta; - if (transition_point > delta_min && transition_point < delta_max) { - bounds.push_back(transition_point); - } - - // Integrate kernel - irf = integral.romberg(bounds, iter_delta); - - // Compile option: Check for NaN/Inf - #if defined(G_NAN_CHECK) - if (gammalib::is_notanumber(irf) || gammalib::is_infinite(irf)) { - std::cout << "*** ERROR: GCTAResponseIrf::irf_radial:"; - std::cout << " NaN/Inf encountered"; - std::cout << " (irf=" << irf; - std::cout << ", delta_min=" << delta_min; - std::cout << ", delta_max=" << delta_max << ")"; - std::cout << std::endl; - } - #endif - - // Apply deadtime correction - irf *= obs.deadc(srcTime); - - } // endif: delta interval has positive length - - // Compile option: Show integration results - #if defined(G_DEBUG_IRF_RADIAL) - std::cout << "GCTAResponseIrf::irf_radial:"; - std::cout << " delta=[" << delta_min*gammalib::rad2deg; - std::cout << ", " << delta_max*gammalib::rad2deg << " deg;"; - std::cout << " zeta=" << zeta*gammalib::rad2deg << " deg;"; - std::cout << " phi0=" << phi0*gammalib::rad2deg << " deg;"; - std::cout << " theta_max=" << theta_max*gammalib::rad2deg << " deg;"; - std::cout << " obsDir=" << obsDir << ";"; - std::cout << " modelDir=" << model->dir() << ";"; - std::cout << " irf=" << irf << std::endl; - #endif - - // Return Irf value - return irf; -} -#else -/***********************************************************************//** - * @brief Return IRF value for radial source model - * - * @param[in] event Observed event. - * @param[in] source Source. - * @param[in] obs Observation. - * - * @exception GCTAException::bad_model_type - * Model is not a radial model. - * * Integrates the product of the model and the Irf over the true photon * arrival direction using * @@ -1617,7 +1449,7 @@ double GCTAResponseIrf::irf_radial(const GEvent& event, // location const GModelSpatialRadialShell* shell = dynamic_cast(model); if (shell != NULL) { - double shell_radius = shell->radius(); + double shell_radius = shell->radius() * gammalib::deg2rad; if (shell_radius > rho_min && shell_radius < rho_max) { bounds.push_back(shell_radius); } @@ -1663,7 +1495,6 @@ double GCTAResponseIrf::irf_radial(const GEvent& event, // Return IRF value return irf; } -#endif /***********************************************************************//** @@ -1824,16 +1655,6 @@ double GCTAResponseIrf::irf_elliptical(const GEvent& event, bounds.push_back(rho_min); bounds.push_back(rho_max); - // If the integration range includes a transition between full - // containment and partial containment, then add a boundary at - // this location - /* - double transition_point = delta_max - rho_obs; - if (transition_point > rho_min && transition_point < rho_max) { - bounds.push_back(transition_point); - } - */ - // If the integration range includes the semiminor boundary, then // add an integration boundary at that location if (semiminor > rho_min && semiminor < rho_max) { @@ -2224,7 +2045,7 @@ double GCTAResponseIrf::npred_radial(const GSource& source, // location const GModelSpatialRadialShell* shell = dynamic_cast(model); if (shell != NULL) { - double shell_radius = shell->radius(); + double shell_radius = shell->radius() * gammalib::deg2rad; if (shell_radius > rho_min && shell_radius < rho_max) { bounds.push_back(shell_radius); } diff --git a/inst/cta/src/GCTAResponse_helpers.cpp b/inst/cta/src/GCTAResponse_helpers.cpp index 905fab7..a912a16 100644 --- a/inst/cta/src/GCTAResponse_helpers.cpp +++ b/inst/cta/src/GCTAResponse_helpers.cpp @@ -283,8 +283,15 @@ double cta_irf_radial_kern_rho::eval(const double& rho) double omega_min = -domega; double omega_max = +domega; + // Reduce rho by an infinite amount to avoid rounding errors + // at the boundary of a sharp edged model + double rho_kluge = rho - 1.0e-12; + if (rho_kluge < 0.0) { + rho_kluge = 0.0; + } + // Evaluate sky model - double model = m_model.eval(rho, m_srcEng, m_srcTime); + double model = m_model.eval(rho_kluge, m_srcEng, m_srcTime); // Debug: test if model is non positive #if defined(G_DEBUG_MODEL_ZERO) @@ -292,7 +299,7 @@ double cta_irf_radial_kern_rho::eval(const double& rho) std::cout << "*** WARNING: cta_irf_radial_kern_rho::eval"; std::cout << " zero model for (rho)=("; std::cout << rho*gammalib::rad2deg << ")"; - std::cout << " rho-r_model" << (rho-m_model.theta_max()); + std::cout << " rho-r_model=" << (rho-m_model.theta_max()); std::cout << " radians" << std::endl; } #endif @@ -474,8 +481,15 @@ double cta_npred_radial_kern_rho::eval(const double& rho) double omega_min = m_omega0 - domega; double omega_max = m_omega0 + domega; + // Reduce rho by an infinite amount to avoid rounding errors + // at the boundary of a sharp edged model + double rho_kluge = rho - 1.0e-12; + if (rho_kluge < 0.0) { + rho_kluge = 0.0; + } + // Get radial model value - double model = m_model.eval(rho, m_srcEng, m_srcTime); + double model = m_model.eval(rho_kluge, m_srcEng, m_srcTime); // Debug: test if model is non positive #if defined(G_DEBUG_MODEL_ZERO) @@ -483,7 +497,7 @@ double cta_npred_radial_kern_rho::eval(const double& rho) std::cout << "*** WARNING: cta_npred_radial_kern_rho::eval"; std::cout << " zero model for (rho)=("; std::cout << rho*gammalib::rad2deg << ")"; - std::cout << " rho-r_model" << (rho-m_model.theta_max()); + std::cout << " rho-r_model=" << (rho-m_model.theta_max()); std::cout << " radians" << std::endl; } #endif @@ -1500,237 +1514,11 @@ double cta_npred_diffuse_kern_phi::eval(const double& phi) } -/***********************************************************************//** - * @brief Kernel for PSF integration of spatial model - * - * @param[in] delta PSF offset angle (radians). - * @return Azimuthally integrated product between PSF and model. - * - * Computes the azimuthally integrated product of point spread function and - * the spatial model intensity. As the PSF is azimuthally symmetric, it is - * not included in the azimuthally integration, but just multiplied on the - * azimuthally integrated model. The method returns thus - * - * \f[ - * {\rm PSF}(\delta) \times - * \int_0^{2\pi} {\rm M}(\delta, \phi) \sin \delta {\rm d}\phi - * \f] - * - * where \f${\rm M}(\delta, \phi)\f$ is the spatial model in the coordinate - * system of the point spread function, defined by the angle \f$\delta\f$ - * between the true and the measured photon direction and the azimuth angle - * \f$\phi\f$ around the measured photon direction. - ***************************************************************************/ -double cta_psf_diffuse_kern_delta::eval(const double& delta) -{ - // Initialise value - double value = 0.0; - - // If we're at the PSF peak the model is zero (due to the sin(delta) - // term. We thus only integrate for positive deltas. - if (delta > 0.0) { - - // Get PSF for this delta - value = m_rsp->psf()(m_srcDir, delta, m_srcEng); - - // Continue only if PSF is positive - if (value > 0.0) { - - // Compute sine and cosine of delta - double sin_delta = std::sin(delta); - double cos_delta = std::cos(delta); - - // Setup kernel for azimuthal integration of the diffuse model - cta_psf_diffuse_kern_phi integrand(m_model, m_srcEng, m_srcTime, - m_rot, - sin_delta, cos_delta); - - // Azimuthally integrate model - GIntegral integral(&integrand); - integral.eps(m_eps); - value *= integral.romberg(0.0, gammalib::twopi, m_order) * sin_delta; - - } // endif: PSF value was positive - - } // endif: delta was positive - - // Debug: Check for NaN - #if defined(G_NAN_CHECK) - if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: cta_psf_diffuse_kern_delta::eval"; - std::cout << "(delta=" << delta << "):"; - std::cout << " NaN/Inf encountered"; - std::cout << " (value=" << value; - std::cout << ")" << std::endl; - } - #endif - - // Return kernel value - return value; -} - - -/***********************************************************************//** - * @brief Kernel for map integration of spatial model - * - * @param[in] phi Azimuth angle (radians). - * @return Spatial model value. - * - * Computes the value of the spatial model at the position (delta,phi) given - * in point spread function coordinates. The transformation from point - * spread function coordinates into sky coordinates is done using a rotation - * matrix that is pre-computed on entry. - ***************************************************************************/ -double cta_psf_diffuse_kern_phi::eval(const double& phi) -{ - // Compute sky direction vector in native coordinates - double cos_phi = std::cos(phi); - double sin_phi = std::sin(phi); - GVector native(-cos_phi*m_sin_delta, sin_phi*m_sin_delta, m_cos_delta); - - // Rotate from native into celestial system - GVector cel = m_rot * native; - - // Set sky direction - GSkyDir srcDir; - srcDir.celvector(cel); - - // Set photon - GPhoton photon(srcDir, m_srcEng, m_srcTime); - - // Compute map value this sky direction - double value = m_model->eval(photon); - - // Debug: Check for NaN - #if defined(G_NAN_CHECK) - if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: cta_psf_diffuse_kern_phi::eval"; - std::cout << "(phi=" << phi << "):"; - std::cout << " NaN/Inf encountered"; - std::cout << " (value=" << value; - std::cout << ")" << std::endl; - } - #endif - - // Return kernel value - return value; -} - - -/***********************************************************************//** - * @brief Kernel for PSF integration of radial model - * - * @param[in] delta PSF offset angle (radians). - * @return Azimuthally integrated product between PSF and radial model. - * - * Computes the azimuthally integrated product of point spread function and - * the radial model intensity. As the PSF is azimuthally symmetric, it is - * not included in the azimuthally integration, but just multiplied on the - * azimuthally integrated model. The method returns thus - * - * \f[ - * {\rm PSF}(\delta) \times - * \int_0^{2\pi} {\rm M}(\delta, \phi) \sin \delta {\rm d}\phi - * \f] - * - * where \f${\rm M}(\delta, \phi)\f$ is the radial model in the coordinate - * system of the point spread function, defined by the angle \f$\delta\f$ - * between the true and the measured photon direction and the azimuth angle - * \f$\phi\f$ around the measured photon direction. - ***************************************************************************/ -double cta_psf_radial_kern_delta::eval(const double& delta) -{ - // Initialise value - double value = 0.0; - - // If we're at the PSF peak the model is zero (due to the sin(delta) - // term. We thus only integrate for positive deltas. - if (delta > 0.0) { - - // Get PSF for this delta - value = m_rsp->psf()(m_srcDir, delta, m_srcEng); - - // Continue only if PSF is positive - if (value > 0.0) { - - // Compute sin(delta), sin(delta)*sin(zeta), cos(delta)*cos(zeta) - double sin_delta = std::sin(delta); - double sin_fact = sin_delta * m_sin_zeta; - double cos_fact = std::cos(delta)* m_cos_zeta; - - // Setup kernel for azimuthal integration of the spatial model - cta_psf_radial_kern_phi integrand(m_model, m_srcEng, m_srcTime, - sin_fact, cos_fact); - - // Azimuthally integrate model - GIntegral integral(&integrand); - integral.eps(m_eps); - value *= integral.romberg(0.0, gammalib::twopi, m_order) * sin_delta; - - } // endif: PSF value was positive - - } // endif: delta was positive - - // Debug: Check for NaN - #if defined(G_NAN_CHECK) - if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: cta_psf_radial_kern_delta::eval"; - std::cout << "(delta=" << delta << "):"; - std::cout << " NaN/Inf encountered"; - std::cout << " (value=" << value; - std::cout << ")" << std::endl; - } - #endif - - // Return kernel value - return value; -} - - -/***********************************************************************//** - * @brief Kernel for azimuthal radial model integration - * - * @param[in] phi Azimuth angle (radians). - * @return Radial model value. - * - * Computes the value of the radial model at the position \f$(\delta,\phi)\f$ - * given in point spread function coordinates. The \f$\theta\f$ angle of the - * radial model is computed using - * - * \f[ - * \theta = \arccos \left( \cos \delta \cos \zeta + - * \sin \delta \sin \zeta \cos \phi - * \f] - * - * where \f$\delta\f$ is the angle between true and measured photon - * direction, \f$\zeta\f$ is the angle between model centre and measured - * photon direction, and \f$\phi\f$ is the azimuth angle with respect to the - * measured photon direction, where \f$\phi=0\f$ corresponds to the - * connecting line between model centre and measured photon direction. - ***************************************************************************/ -double cta_psf_radial_kern_phi::eval(const double& phi) -{ - // Compute radial model theta angle - double theta = std::acos(m_cos_fact + m_sin_fact * std::cos(phi)); - - // Get radial model value - double value = m_model->eval(theta, m_srcEng, m_srcTime); - - // Debug: Check for NaN - #if defined(G_NAN_CHECK) - if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: cta_psf_radial_kern_phi::eval"; - std::cout << "(phi=" << phi << "):"; - std::cout << " NaN/Inf encountered"; - std::cout << " (value=" << value; - std::cout << ")" << std::endl; - } - #endif - - // Return kernel value - return value; -} - +/*========================================================================== + = = + = Helper class methods for stacked analysis = + = = + ==========================================================================*/ /***********************************************************************//** * @brief Kernel for radial model integration over zenith angle @@ -1777,8 +1565,7 @@ double cta_psf_radial_kern_rho::eval(const double& rho) if (domega > 0.0) { // Reduce rho by an infinite amount to avoid rounding errors - // at the boundary of a sharp edged model (e.g. an elliptical - // disk model) + // at the boundary of a sharp edged model double rho_kluge = rho - 1.0e-12; if (rho_kluge < 0.0) { rho_kluge = 0.0; @@ -1793,7 +1580,7 @@ double cta_psf_radial_kern_rho::eval(const double& rho) std::cout << "*** WARNING: cta_psf_radial_kern_rho::eval"; std::cout << " zero model for (rho)=("; std::cout << rho*gammalib::rad2deg << ")"; - std::cout << " rho-r_model" << (rho-m_model->theta_max()); + std::cout << " rho-r_model=" << (rho-m_model->theta_max()); std::cout << " radians" << std::endl; } #endif @@ -1912,133 +1699,6 @@ double cta_psf_radial_kern_omega::eval(const double& omega) /***********************************************************************//** - * @brief Kernel for PSF integration of elliptical model - * - * @param[in] delta PSF offset angle (radians). - * @return Azimuthally integrated product between PSF and elliptical model. - * - * Computes the azimuthally integrated product of point spread function and - * the elliptical model intensity. As the PSF is azimuthally symmetric, it is - * not included in the azimuthally integration, but just multiplied on the - * azimuthally integrated model. The method returns thus - * - * \f[ - * {\rm PSF}(\delta) \times - * \int_0^{2\pi} {\rm M}(\delta, \phi) \sin \delta {\rm d}\phi - * \f] - * - * where \f${\rm M}(\delta, \phi)\f$ is the elliptical model in the coordinate - * system of the point spread function, defined by the angle \f$\delta\f$ - * between the true and the measured photon direction and the azimuth angle - * \f$\phi\f$ around the measured photon direction. - ***************************************************************************/ -double cta_psf_elliptical_kern_delta::eval(const double& delta) -{ - // Initialise value - double value = 0.0; - - // If we're at the PSF peak the model is zero (due to the sin(delta) - // term. We thus only integrate for positive deltas. - if (delta > 0.0) { - - // Get PSF for this delta - value = m_rsp->psf()(m_srcDir, delta, m_srcEng); - - // Continue only if PSF is positive - if (value > 0.0) { - - // Compute sin(delta), sin(delta)*sin(zeta), cos(delta)*cos(zeta) - double sin_delta = std::sin(delta); - double sin_fact = sin_delta * m_sin_zeta; - double cos_fact = std::cos(delta)* m_cos_zeta; - - // Setup kernel for azimuthal integration of the elliptical model - cta_psf_elliptical_kern_phi integrand(m_model, m_srcEng, m_srcTime, - m_omega, - sin_delta, sin_fact, cos_fact); - - // Azimuthally integrate model - GIntegral integral(&integrand); - integral.eps(m_eps); - value *= integral.romberg(0.0, gammalib::twopi, m_order) * sin_delta; - - } // endif: PSF value was positive - - } // endif: delta was positive - - // Debug: Check for NaN - #if defined(G_NAN_CHECK) - if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: cta_psf_elliptical_kern_delta::eval"; - std::cout << "(delta=" << delta << "):"; - std::cout << " NaN/Inf encountered"; - std::cout << " (value=" << value; - std::cout << ")" << std::endl; - } - #endif - - // Return kernel value - return value; -} - - -/***********************************************************************//** - * @brief Kernel for azimuthal elliptical model integration - * - * @param[in] phi Azimuth angle (radians). - * @return Radial model value. - * - * Computes the value of the elliptical model at the position \f$(\delta,\phi)\f$ - * given in point spread function coordinates. The \f$\theta\f$ angle of the - * elliptical model is computed using - * - * \f[ - * \theta = \arccos \left( \cos \delta \cos \zeta + - * \sin \delta \sin \zeta \cos \phi \right) - * \f] - * - * where \f$\delta\f$ is the angle between true and measured photon - * direction, \f$\zeta\f$ is the angle between model centre and measured - * photon direction, and \f$\phi\f$ is the azimuth angle with respect to the - * measured photon direction, where \f$\phi=0\f$ corresponds to the - * connecting line between model centre and measured photon direction. - * - * The position angle of the elliptical model is computed using - * - * \f[ - * {\rm posangle} = \arcsin \left( \frac{\sin \phi \sin \delta} - * {\sin \theta} \right) - * \f] - ***************************************************************************/ -double cta_psf_elliptical_kern_phi::eval(const double& phi) -{ - // Compute radial model theta angle - double theta = std::acos(m_cos_fact + m_sin_fact * std::cos(phi)); - - // Compute radial model position angle - double posangle = std::asin((std::sin(phi) * m_sin_delta)/std::sin(theta)) + - m_omega; - - // Get radial model value - double value = m_model->eval(theta, posangle, m_srcEng, m_srcTime); - - // Debug: Check for NaN - #if defined(G_NAN_CHECK) - if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: cta_psf_elliptical_kern_phi::eval"; - std::cout << "(phi=" << phi << "):"; - std::cout << " NaN/Inf encountered"; - std::cout << " (value=" << value; - std::cout << ")" << std::endl; - } - #endif - - // Return kernel value - return value; -} - - -/***********************************************************************//** * @brief Kernel for elliptical model integration over zenith angle * * @param[in] rho Radial distance from model centre (radians). @@ -2305,112 +1965,64 @@ double cta_psf_elliptical_kern_omega::eval(const double& omega) } -/*========================================================================== - = = - = Helper class methods for response computation testing = - = = - = WARNING: These helper classes are for testing purposes only and not for = - = production. These classes may be removed at any moment without = - = any warning. = - ==========================================================================*/ - - /***********************************************************************//** - * @brief Kernel for integration of radial model in Psf system + * @brief Kernel for PSF integration of spatial model * - * @param[in] delta Psf offset angle (radians). - * @return Azimuthally integrated product between Psf and radial model. + * @param[in] delta PSF offset angle (radians). + * @return Azimuthally integrated product between PSF and model. * * Computes the azimuthally integrated product of point spread function and - * the radial model intensity. As the Psf is azimuthally symmetric, it is + * the spatial model intensity. As the PSF is azimuthally symmetric, it is * not included in the azimuthally integration, but just multiplied on the * azimuthally integrated model. The method returns thus * * \f[ - * {\rm Psf}(\delta) \times + * {\rm PSF}(\delta) \times * \int_0^{2\pi} {\rm M}(\delta, \phi) \sin \delta {\rm d}\phi * \f] * - * where \f${\rm M}(\delta, \phi)\f$ is the radial model in the coordinate + * where \f${\rm M}(\delta, \phi)\f$ is the spatial model in the coordinate * system of the point spread function, defined by the angle \f$\delta\f$ * between the true and the measured photon direction and the azimuth angle * \f$\phi\f$ around the measured photon direction. ***************************************************************************/ -double cta_irf_radial_kern_delta::eval(const double& delta) +double cta_psf_diffuse_kern_delta::eval(const double& delta) { // Initialise value double value = 0.0; - // If we're at the Psf peak the model is zero (due to the sin(delta) + // If we're at the PSF peak the model is zero (due to the sin(delta) // term. We thus only integrate for positive deltas. if (delta > 0.0) { - // Get Psf for this delta. We use here the event offset instead of - // the true photon offset because the Psf varies little over the - // dimension of the Psf, hence it should basically be invariant - // over such small scales - value = m_rsp->psf(delta, m_obsOffset, 0.0, m_pnt.zenith(), m_pnt.azimuth(), m_srcLogEng); - - // Compile option: Get Aeff now using the event offset instead of - // computing Aeff later in the azimuthal integration kernel where - // the true photon offset is available. This is an approximation and - // assumes that Aeff varies little over the dimension of the Psf. - // This option is only for testing purposes and should not be used - // for production. - #if defined(G_USE_OBSDIR_FOR_AEFF) - value *= m_rsp->aeff(m_obsOffset, 0.0, m_pnt.zenith(), m_pnt.azimuth(), m_srcLogEng); - #endif + // Get PSF for this delta + value = m_rsp->psf()(m_srcDir, delta, m_srcEng); - // Continue only if Psf is positive + // Continue only if PSF is positive if (value > 0.0) { - // Compute half length of arc that lies within model (radians) - double dphi = 0.5 * gammalib::cta_roi_arclength(delta, - m_zeta, - m_cos_zeta, - m_sin_zeta, - m_theta_max, - m_cos_theta_max); - - // Continue only if arc length is positive - if (dphi > 0.0) { - - // Compute azimuth integration range [phi_min,phi_max]. - double phi_min = -dphi; - double phi_max = +dphi; - - // Pre-compute terms for angular distance computation in - // azimuth integration kernel. - double sin_delta = std::sin(delta); - double cos_delta = std::cos(delta); - double sin_fact_model = sin_delta * m_sin_zeta; - double cos_fact_model = cos_delta * m_cos_zeta; - double sin_fact_inst = sin_delta * m_sin_obs_offset; - double cos_fact_inst = cos_delta * m_cos_obs_offset; - - // Setup kernel for azimuthal integration of the spatial model - cta_irf_radial_kern_phi integrand(m_rsp, m_model, m_pnt, - m_srcEng, m_srcLogEng, m_srcTime, - m_obsEng, - m_phi0, - sin_fact_model, cos_fact_model, - sin_fact_inst, cos_fact_inst); - - // Integrate kernel - GIntegral integral(&integrand); - integral.fixed_iter(m_iter); - value *= integral.romberg(phi_min, phi_max, m_iter) * sin_delta; + // Compute sine and cosine of delta + double sin_delta = std::sin(delta); + double cos_delta = std::cos(delta); + + // Setup kernel for azimuthal integration of the diffuse model + cta_psf_diffuse_kern_phi integrand(m_model, m_srcEng, m_srcTime, + m_rot, + sin_delta, cos_delta); - } // endif: arc length was positive + // Azimuthally integrate model + GIntegral integral(&integrand); + integral.eps(m_eps); + value *= integral.romberg(0.0, gammalib::twopi, m_order) * sin_delta; - } // endif: Psf value was positive + } // endif: PSF value was positive } // endif: delta was positive // Debug: Check for NaN #if defined(G_NAN_CHECK) if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: cta_irf_radial_kern_delta::eval"; + std::cout << "*** ERROR: cta_psf_diffuse_kern_delta::eval"; std::cout << "(delta=" << delta << "):"; std::cout << " NaN/Inf encountered"; std::cout << " (value=" << value; @@ -2424,79 +2036,40 @@ double cta_irf_radial_kern_delta::eval(const double& delta) /***********************************************************************//** - * @brief Kernel for azimuthal radial model integration + * @brief Kernel for map integration of spatial model * * @param[in] phi Azimuth angle (radians). - * @return Radial model value. - * - * Computes the value of the radial model at the position \f$(\delta,\phi)\f$ - * given in point spread function coordinates. The \f$\theta\f$ angle of the - * radial model is computed using - * - * \f[ - * \theta = \arccos \left( \cos \delta \cos \zeta + - * \sin \delta \sin \zeta \cos \phi - * \f] + * @return Spatial model value. * - * where \f$\delta\f$ is the angle between true and measured photon - * direction, \f$\zeta\f$ is the angle between model centre and measured - * photon direction, and \f$\phi\f$ is the azimuth angle with respect to the - * measured photon direction, where \f$\phi=0\f$ corresponds to the - * connecting line between model centre and measured photon direction. + * Computes the value of the spatial model at the position (delta,phi) given + * in point spread function coordinates. The transformation from point + * spread function coordinates into sky coordinates is done using a rotation + * matrix that is pre-computed on entry. ***************************************************************************/ -double cta_irf_radial_kern_phi::eval(const double& phi) +double cta_psf_diffuse_kern_phi::eval(const double& phi) { - // Compute radial model theta angle - double theta = gammalib::acos(m_cos_fact_model + - m_sin_fact_model * std::cos(phi)); - - // Get radial model value - double value = m_model->eval(theta, m_srcEng, m_srcTime); - - // For debugging: Check boundary violation - if (value == 0.0) { - std::cout << "Model 0 at theta=" << theta*gammalib::rad2deg; - std::cout << " deg; phi="; - std::cout << phi*gammalib::rad2deg << " deg"; - std::cout << std::endl; - throw; - } - - // If model is positive the multiply it with the effective area and - // optionally fold in energy dispersion - #if defined(G_USE_OBSDIR_FOR_AEFF) - if (m_rsp->use_edisp() && value > 0.0) { - - // Compute offset angle of true photon in camera system (radians) - double offset = gammalib::acos(m_cos_fact_inst + - m_sin_fact_inst * std::cos(m_phi0 - phi)); - - // Take energy dispersion into account - value *= m_rsp->edisp(m_obsEng, offset, 0.0, m_pnt.zenith(), m_pnt.azimuth(), m_srcLogEng); + // Compute sky direction vector in native coordinates + double cos_phi = std::cos(phi); + double sin_phi = std::sin(phi); + GVector native(-cos_phi*m_sin_delta, sin_phi*m_sin_delta, m_cos_delta); - } // endif: model value was positive - #else - if (value > 0.0) { + // Rotate from native into celestial system + GVector cel = m_rot * native; - // Compute offset angle of true photon in camera system (radians) - double offset = gammalib::acos(m_cos_fact_inst + - m_sin_fact_inst * std::cos(m_phi0 - phi)); + // Set sky direction + GSkyDir srcDir; + srcDir.celvector(cel); - // Multiply with effective area - value *= m_rsp->aeff(offset, 0.0, m_pnt.zenith(), m_pnt.azimuth(), m_srcLogEng); + // Set photon + GPhoton photon(srcDir, m_srcEng, m_srcTime); - // Optionally take energy dispersion into account - if (m_rsp->use_edisp() && value > 0.0) { - value *= m_rsp->edisp(m_obsEng, offset, 0.0, m_pnt.zenith(), m_pnt.azimuth(), m_srcLogEng); - } + // Compute map value this sky direction + double value = m_model->eval(photon); - } // endif: model value was positive - #endif - // Debug: Check for NaN #if defined(G_NAN_CHECK) if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { - std::cout << "*** ERROR: cta_irf_radial_kern_phi::eval"; + std::cout << "*** ERROR: cta_psf_diffuse_kern_phi::eval"; std::cout << "(phi=" << phi << "):"; std::cout << " NaN/Inf encountered"; std::cout << " (value=" << value; diff --git a/inst/cta/src/GCTAResponse_helpers.hpp b/inst/cta/src/GCTAResponse_helpers.hpp index 18175ed..49d677f 100644 --- a/inst/cta/src/GCTAResponse_helpers.hpp +++ b/inst/cta/src/GCTAResponse_helpers.hpp @@ -982,142 +982,6 @@ protected: /***********************************************************************//** - * @class cta_psf_diffuse_kern_delta - * - * @brief Kernel for Psf delta angle integration used for stacked analysis - ***************************************************************************/ -class cta_psf_diffuse_kern_delta : public GFunction { -public: - cta_psf_diffuse_kern_delta(const GCTAResponseCube* rsp, - const GModelSpatial* model, - const GSkyDir& srcDir, - const GEnergy& srcEng, - const GTime& srcTime, - const GMatrix& rot, - const double& eps, - const int& order) : - m_rsp(rsp), - m_model(model), - m_srcDir(srcDir), - m_srcEng(srcEng), - m_srcTime(srcTime), - m_rot(rot), - m_eps(eps), - m_order (order), - m_psf_max(rsp->psf()(srcDir, 0.0, srcEng)) { } - double eval(const double& delta); -protected: - const GCTAResponseCube* m_rsp; //!< Response cube - const GModelSpatial* m_model; //!< Spatial model - const GSkyDir& m_srcDir; //!< True photon arrival direction - const GEnergy& m_srcEng; //!< True photon energy - const GTime& m_srcTime; //!< True photon arrival time - const GMatrix& m_rot; //!< Rotation matrix - const double& m_eps; //!< Integration precision - const int& m_order; //!< Romberg order - double m_psf_max; //!< Maximum PSF value -}; - - -/***********************************************************************//** - * @class cta_psf_diffuse_kern_phi - * - * @brief Kernel for Psf phi angle integration used for stacked analysis - ***************************************************************************/ -class cta_psf_diffuse_kern_phi : public GFunction { -public: - cta_psf_diffuse_kern_phi(const GModelSpatial* model, - const GEnergy& srcEng, - const GTime& srcTime, - const GMatrix& rot, - const double& sin_delta, - const double& cos_delta) : - m_model(model), - m_srcEng(srcEng), - m_srcTime(srcTime), - m_rot(rot), - m_sin_delta(sin_delta), - m_cos_delta(cos_delta) { } - double eval(const double& phi); -protected: - const GModelSpatial* m_model; //!< Spatial model - const GEnergy& m_srcEng; //!< True photon energy - const GTime& m_srcTime; //!< True photon arrival time - const GMatrix& m_rot; //!< Rotation matrix - const double& m_sin_delta; //!< sin(delta) - const double& m_cos_delta; //!< cos(delta) -}; - - -/***********************************************************************//** - * @class cta_psf_radial_kern_delta - * - * @brief Kernel for Psf delta angle integration used for stacked analysis - ***************************************************************************/ -class cta_psf_radial_kern_delta : public GFunction { -public: - cta_psf_radial_kern_delta(const GCTAResponseCube* rsp, - const GModelSpatialRadial* model, - const GSkyDir& srcDir, - const GEnergy& srcEng, - const GTime& srcTime, - const double& sin_zeta, - const double& cos_zeta, - const double& eps, - const int& order) : - m_rsp(rsp), - m_model(model), - m_srcDir(srcDir), - m_srcEng(srcEng), - m_srcTime(srcTime), - m_sin_zeta(sin_zeta), - m_cos_zeta(cos_zeta), - m_eps(eps), - m_order (order), - m_psf_max(rsp->psf()(srcDir, 0.0, srcEng)) { } - double eval(const double& delta); -protected: - const GCTAResponseCube* m_rsp; //!< Response cube - const GModelSpatialRadial* m_model; //!< Radial model - const GSkyDir& m_srcDir; //!< True photon arrival direction - const GEnergy& m_srcEng; //!< True photon energy - const GTime& m_srcTime; //!< True photon arrival time - const double& m_sin_zeta; //!< sin(zeta) - const double& m_cos_zeta; //!< cos(zeta) - const double& m_eps; //!< Integration precision - const int& m_order; //!< Romberg order - double m_psf_max; //!< Maximum PSF value -}; - - -/***********************************************************************//** - * @class cta_psf_radial_kern_phi - * - * @brief Kernel for Psf phi angle integration used for stacked analysis - ***************************************************************************/ -class cta_psf_radial_kern_phi : public GFunction { -public: - cta_psf_radial_kern_phi(const GModelSpatialRadial* model, - const GEnergy& srcEng, - const GTime& srcTime, - const double& sin_fact, - const double& cos_fact) : - m_model(model), - m_srcEng(srcEng), - m_srcTime(srcTime), - m_sin_fact(sin_fact), - m_cos_fact(cos_fact) { } - double eval(const double& phi); -protected: - const GModelSpatialRadial* m_model; //!< Radial model - const GEnergy& m_srcEng; //!< True photon energy - const GTime& m_srcTime; //!< True photon arrival time - const double& m_sin_fact; //!< sin(delta)*sin(zeta) - const double& m_cos_fact; //!< cos(delta)*cos(zeta) -}; - - -/***********************************************************************//** * @class cta_irf_radial_kern_rho * * @brief Kernel for radial model zenith angle integration @@ -1227,83 +1091,6 @@ public: /***********************************************************************//** - * @class cta_psf_elliptical_kern_delta - * - * @brief Kernel for Psf delta angle integration used for stacked analysis - ***************************************************************************/ -class cta_psf_elliptical_kern_delta : public GFunction { -public: - cta_psf_elliptical_kern_delta(const GCTAResponseCube* rsp, - const GModelSpatialElliptical* model, - const GSkyDir& srcDir, - const GEnergy& srcEng, - const GTime& srcTime, - const double& sin_zeta, - const double& cos_zeta, - const double& omega, - const double& eps, - const int& order) : - m_rsp(rsp), - m_model(model), - m_srcDir(srcDir), - m_srcEng(srcEng), - m_srcTime(srcTime), - m_sin_zeta(sin_zeta), - m_cos_zeta(cos_zeta), - m_omega(omega), - m_eps(eps), - m_order (order), - m_psf_max(rsp->psf()(srcDir, 0.0, srcEng)) { } - double eval(const double& delta); -protected: - const GCTAResponseCube* m_rsp; //!< Response cube - const GModelSpatialElliptical* m_model; //!< Elliptical model - const GSkyDir& m_srcDir; //!< True photon arrival direction - const GEnergy& m_srcEng; //!< True photon energy - const GTime& m_srcTime; //!< True photon arrival time - const double& m_sin_zeta; //!< sin(zeta) - const double& m_cos_zeta; //!< cos(zeta) - const double& m_omega; //!< Position angle of photon w/r model - const double& m_eps; //!< Integration precision - const int& m_order; //!< Romberg order - double m_psf_max; //!< Maximum PSF value -}; - - -/***********************************************************************//** - * @class cta_psf_elliptical_kern_delta - * - * @brief Kernel for Psf phi angle integration used for stacked analysis - ***************************************************************************/ -class cta_psf_elliptical_kern_phi : public GFunction { -public: - cta_psf_elliptical_kern_phi(const GModelSpatialElliptical* model, - const GEnergy& srcEng, - const GTime& srcTime, - const double& omega, - const double& sin_delta, - const double& sin_fact, - const double& cos_fact) : - m_model(model), - m_srcEng(srcEng), - m_srcTime(srcTime), - m_omega(omega), - m_sin_delta(sin_delta), - m_sin_fact(sin_fact), - m_cos_fact(cos_fact) { } - double eval(const double& phi); -protected: - const GModelSpatialElliptical* m_model; //!< Elliptical model - const GEnergy& m_srcEng; //!< True photon energy - const GTime& m_srcTime; //!< True photon arrival time - const double& m_omega; //!< Position angle of photon w/r model - const double& m_sin_delta; //!< sin(delta) - const double& m_sin_fact; //!< sin(delta)*sin(zeta) - const double& m_cos_fact; //!< cos(delta)*cos(zeta) -}; - - -/***********************************************************************//** * @class cta_irf_elliptical_kern_rho * * @brief Kernel for elliptical model zenith angle integration @@ -1433,114 +1220,70 @@ public: /***********************************************************************//** - * @class cta_irf_radial_kern_delta - * - * @brief Test kernel for Irf delta angle integration + * @class cta_psf_diffuse_kern_delta * - * This class is used for Psf system based integration. It is only here for - * testing and this class may be removed in the future. + * @brief Kernel for Psf delta angle integration used for stacked analysis ***************************************************************************/ -class cta_irf_radial_kern_delta : public GFunction { +class cta_psf_diffuse_kern_delta : public GFunction { public: - cta_irf_radial_kern_delta(const GCTAResponseIrf* rsp, - const GModelSpatialRadial* model, - const GCTAPointing& pnt, - const GEnergy& srcEng, - const double& srcLogEng, - const GTime& srcTime, - const double& obsOffset, - const GEnergy& obsEng, - const double& zeta, - const double& phi0, - const double& theta_max, - const int& iter) : - m_rsp(rsp), - m_model(model), - m_pnt(pnt), - m_srcEng(srcEng), - m_srcLogEng(srcLogEng), - m_srcTime(srcTime), - m_obsOffset(obsOffset), - m_sin_obs_offset(std::sin(obsOffset)), - m_cos_obs_offset(std::cos(obsOffset)), - m_obsEng(obsEng), - m_zeta(zeta), - m_phi0(phi0), - m_sin_zeta(std::sin(zeta)), - m_cos_zeta(std::cos(zeta)), - m_theta_max(theta_max), - m_cos_theta_max(std::cos(theta_max)), - m_iter(iter) { } + cta_psf_diffuse_kern_delta(const GCTAResponseCube* rsp, + const GModelSpatial* model, + const GSkyDir& srcDir, + const GEnergy& srcEng, + const GTime& srcTime, + const GMatrix& rot, + const double& eps, + const int& order) : + m_rsp(rsp), + m_model(model), + m_srcDir(srcDir), + m_srcEng(srcEng), + m_srcTime(srcTime), + m_rot(rot), + m_eps(eps), + m_order (order), + m_psf_max(rsp->psf()(srcDir, 0.0, srcEng)) { } double eval(const double& delta); protected: - const GCTAResponseIrf* m_rsp; //!< Response - const GModelSpatialRadial* m_model; //!< Radial model - const GCTAPointing& m_pnt; //!< CTA pointing - const GEnergy& m_srcEng; //!< True photon energy - const double& m_srcLogEng; //!< Log10(true photon energy) - const GTime& m_srcTime; //!< True photon arrival time - const double& m_obsOffset; //!< Observed photon offset angle - const double m_sin_obs_offset; - const double m_cos_obs_offset; - const GEnergy& m_obsEng; //!< Measured event energy - const double& m_zeta; //!< zeta - const double& m_phi0; //!< phi0 - const double m_sin_zeta; //!< sin(zeta) - const double m_cos_zeta; //!< cos(zeta) - const double& m_theta_max; //!< Model radius - const double m_cos_theta_max; - const int& m_iter; //!< Integration iterations + const GCTAResponseCube* m_rsp; //!< Response cube + const GModelSpatial* m_model; //!< Spatial model + const GSkyDir& m_srcDir; //!< True photon arrival direction + const GEnergy& m_srcEng; //!< True photon energy + const GTime& m_srcTime; //!< True photon arrival time + const GMatrix& m_rot; //!< Rotation matrix + const double& m_eps; //!< Integration precision + const int& m_order; //!< Romberg order + double m_psf_max; //!< Maximum PSF value }; /***********************************************************************//** - * @class cta_irf_radial_kern_phi - * - * @brief Test kernel for Irf phi angle integration + * @class cta_psf_diffuse_kern_phi * - * This class is used for Psf system based integration. It is only here for - * testing and this class may be removed in the future. + * @brief Kernel for Psf phi angle integration used for stacked analysis ***************************************************************************/ -class cta_irf_radial_kern_phi : public GFunction { +class cta_psf_diffuse_kern_phi : public GFunction { public: - cta_irf_radial_kern_phi(const GCTAResponseIrf* rsp, - const GModelSpatialRadial* model, - const GCTAPointing& pnt, - const GEnergy& srcEng, - const double& srcLogEng, - const GTime& srcTime, - const GEnergy& obsEng, - const double& phi0, - const double& sin_fact_model, - const double& cos_fact_model, - const double& sin_fact_inst, - const double& cos_fact_inst) : - m_rsp(rsp), - m_model(model), - m_pnt(pnt), - m_srcEng(srcEng), - m_srcLogEng(srcLogEng), - m_srcTime(srcTime), - m_obsEng(obsEng), - m_phi0(phi0), - m_sin_fact_model(sin_fact_model), - m_cos_fact_model(cos_fact_model), - m_sin_fact_inst(sin_fact_inst), - m_cos_fact_inst(cos_fact_inst) { } + cta_psf_diffuse_kern_phi(const GModelSpatial* model, + const GEnergy& srcEng, + const GTime& srcTime, + const GMatrix& rot, + const double& sin_delta, + const double& cos_delta) : + m_model(model), + m_srcEng(srcEng), + m_srcTime(srcTime), + m_rot(rot), + m_sin_delta(sin_delta), + m_cos_delta(cos_delta) { } double eval(const double& phi); protected: - const GCTAResponseIrf* m_rsp; //!< Response - const GModelSpatialRadial* m_model; //!< Radial model - const GCTAPointing& m_pnt; //!< CTA pointing - const GEnergy& m_srcEng; //!< True photon energy - const double& m_srcLogEng; //!< Log10(true photon energy) - const GTime& m_srcTime; //!< True photon arrival time - const GEnergy& m_obsEng; //!< Measured event energy - const double& m_phi0; //!< phi0 - const double& m_sin_fact_model; //!< sin(delta)*sin(zeta) - const double& m_cos_fact_model; //!< cos(delta)*cos(zeta) - const double& m_sin_fact_inst; //!< sin(delta)*sin(obsOffset) - const double& m_cos_fact_inst; //!< cos(delta)*cos(obsOffset) + const GModelSpatial* m_model; //!< Spatial model + const GEnergy& m_srcEng; //!< True photon energy + const GTime& m_srcTime; //!< True photon arrival time + const GMatrix& m_rot; //!< Rotation matrix + const double& m_sin_delta; //!< sin(delta) + const double& m_cos_delta; //!< cos(delta) }; #endif /* GCTARESPONSE_HELPERS_HPP */ -- libgit2 0.21.4