From f522be3ac0b17b470e5a56a3e9875ce320c54fbb Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Thu, 30 Oct 2014 00:06:12 +0100 Subject: [PATCH] Finalize response optimization. --- inst/cta/src/GCTAResponseCube.cpp | 285 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++------------------------------------------------------------------------------------------------------- inst/cta/src/GCTAResponseIrf.cpp | 30 ++++++++++-------------------- inst/cta/src/GCTAResponse_helpers.cpp | 165 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--- inst/cta/src/GCTAResponse_helpers.hpp | 69 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ inst/cta/test/test_irf_radial.py | 7 ++++--- inst/cta/test/test_likelihood_profile.py | 12 ++++++------ inst/cta/test/test_npred_elliptical.py | 15 +++++++-------- 7 files changed, 440 insertions(+), 143 deletions(-) diff --git a/inst/cta/src/GCTAResponseCube.cpp b/inst/cta/src/GCTAResponseCube.cpp index a97a7ea..b4cd689 100644 --- a/inst/cta/src/GCTAResponseCube.cpp +++ b/inst/cta/src/GCTAResponseCube.cpp @@ -64,6 +64,7 @@ /* __ Macros _____________________________________________________________ */ /* __ Coding definitions _________________________________________________ */ +#define G_RADIAL_PSF_BASED //!< Use Psf-based integration for radial model /* __ Debug definitions __________________________________________________ */ @@ -237,11 +238,12 @@ GCTAResponseCube* GCTAResponseCube::clone(void) const /***********************************************************************//** - * @brief Return value of instrument response function + * @brief Return instrument response * * @param[in] event Observed event. * @param[in] photon Incident photon. * @param[in] obs Observation. + * @return Instrument response. ***************************************************************************/ double GCTAResponseCube::irf(const GEvent& event, const GPhoton& photon, @@ -308,14 +310,14 @@ double GCTAResponseCube::irf(const GEvent& event, /***********************************************************************//** - * @brief Return value of instrument response function + * @brief Return instrument response * * @param[in] event Event. * @param[in] source Source. * @param[in] obs Observation. + * @return Instrument response. * - * Returns the instrument response function for a given event, source and - * observation. + * Returns the instrument response for a given event, source and observation. ***************************************************************************/ double GCTAResponseCube::irf(const GEvent& event, const GSource& source, @@ -348,15 +350,14 @@ double GCTAResponseCube::irf(const GEvent& event, /***********************************************************************//** - * @brief Return value of point source instrument response function + * @brief Return instrument response to point source * * @param[in] event Observed event. * @param[in] source Source. * @param[in] obs Observation. - * @return Value of instrument response function for a point source. + * @return Instrument response to point source. * - * This method returns the value of the instrument response function for a - * point source. + * Returns the instrument response to a specified point source. ***************************************************************************/ double GCTAResponseCube::irf_ptsrc(const GEvent& event, const GSource& source, @@ -426,15 +427,14 @@ double GCTAResponseCube::irf_ptsrc(const GEvent& event, /***********************************************************************//** - * @brief Return value of radial source instrument response function + * @brief Return instrument response to radial source * * @param[in] event Observed event. * @param[in] source Source. * @param[in] obs Observation. - * @return Value of instrument response function for a radial source. + * @return Instrument response to radial source. * - * This method returns the value of the instrument response function for a - * radial source. Response values are computed on the fly. + * Returns the instrument response to a specified radial source. ***************************************************************************/ double GCTAResponseCube::irf_radial(const GEvent& event, const GSource& source, @@ -497,15 +497,14 @@ double GCTAResponseCube::irf_radial(const GEvent& event, /***********************************************************************//** - * @brief Return value of elliptical source instrument response function + * @brief Return instrument response to elliptical source * * @param[in] event Observed event. * @param[in] source Source. * @param[in] obs Observation. - * @return Value of instrument response function for an elliptical source. + * @return Instrument response to elliptical source. * - * This method returns the value of the instrument response function for an - * elliptical source. Response values are computed on the fly. + * Returns the instrument response to a specified elliptical source. ***************************************************************************/ double GCTAResponseCube::irf_elliptical(const GEvent& event, const GSource& source, @@ -569,21 +568,20 @@ double GCTAResponseCube::irf_elliptical(const GEvent& event, /***********************************************************************//** - * @brief Return value of extended source instrument response function + * @brief Return instrument response to diffuse source * * @param[in] event Observed event. * @param[in] source Source. * @param[in] obs Observation. - * @return Value of instrument response function for an extended source. + * @return Instrument response to diffuse source. * - * This method returns the value of the instrument response function for an - * extended source. It uses a pre-computation cache to store the IRF for - * the spatial model component. + * Returns the instrument response to a specified diffuse source. * - * The pre-computation cache is initialised if no cache has yet been - * allocated, or if at the beginning of a scan over the events, the model - * parameters have changed. The beginning of a scan is defined by an event - * bin index of 0. + * The method uses a pre-computation cache to store the instrument response + * for the spatial model component. The pre-computation cache is initialised + * if no cache has yet been allocated, or if at the beginning of a scan over + * the events, the model parameters have changed. The beginning of a scan is + * defined by an event bin index of 0. ***************************************************************************/ double GCTAResponseCube::irf_diffuse(const GEvent& event, const GSource& source, @@ -652,9 +650,13 @@ double GCTAResponseCube::irf_diffuse(const GEvent& event, * * @param[in] photon Incident photon. * @param[in] obs Observation. + * @return Spatial integral of point spread function. * * @exception GException::feature_not_implemented * Method not implemented. + * + * This method is a dummy method that is required for a class derived from + * GResponse but that needs not to be implemented for a binned response. ***************************************************************************/ double GCTAResponseCube::npred(const GPhoton& photon, const GObservation& obs) const @@ -669,21 +671,21 @@ double GCTAResponseCube::npred(const GPhoton& photon, /***********************************************************************//** - * @brief Read response + * @brief Read response information from Xml element * * @param[in] xml XML element. * * @exception GException::xml_invalid_parnames * Invalid parameter names found in XML element. * - * Reads information for a CTA observation from an XML element. The exposure - * and PSF cubes are specified using + * Reads response information from an Xml element. The Exposure and Psf cubes + * are specified using * - * - * ... - * - * - * + * + * ... + * + * + * * ***************************************************************************/ void GCTAResponseCube::read(const GXmlElement& xml) @@ -738,18 +740,18 @@ void GCTAResponseCube::read(const GXmlElement& xml) /***********************************************************************//** - * @brief Write response + * @brief Write response information into Xml element * * @param[in] xml XML element. * - * Writes information for a CTA observation into an XML element. The exposure - * and PSF cubes are specified using + * Writes response information into an Xml element. The Exposure and Psf + * cubes are specified using * - * - * ... - * - * - * + * + * ... + * + * + * * ***************************************************************************/ void GCTAResponseCube::write(GXmlElement& xml) const @@ -774,10 +776,10 @@ void GCTAResponseCube::write(GXmlElement& xml) const /***********************************************************************//** - * @brief Print CTA response information + * @brief Print response information * * @param[in] chatter Chattiness (defaults to NORMAL). - * @return String containing CTA response information. + * @return String containing response information. ***************************************************************************/ std::string GCTAResponseCube::print(const GChatter& chatter) const { @@ -819,20 +821,6 @@ std::string GCTAResponseCube::print(const GChatter& chatter) const /*========================================================================== = = - = Model type dependent CTA response methods = - = = - ==========================================================================*/ - - -/*========================================================================== - = = - = Low-level CTA response methods = - = = - ==========================================================================*/ - - -/*========================================================================== - = = = Private methods = = = ==========================================================================*/ @@ -858,7 +846,7 @@ void GCTAResponseCube::init_members(void) /***********************************************************************//** * @brief Copy class members * - * @param[in] rsp Response to be copied + * @param[in] rsp Response. ***************************************************************************/ void GCTAResponseCube::copy_members(const GCTAResponseCube& rsp) { @@ -894,10 +882,12 @@ void GCTAResponseCube::free_members(void) /***********************************************************************//** - * @brief Determines the cache index for a given model name + * @brief Determine cache index for model * * @param[in] name Model name. * @return Cache index (-1 if model has not been found). + * + * Determines the cache index for a given model @p name. ***************************************************************************/ int GCTAResponseCube::cache_index(const std::string& name) const { @@ -922,8 +912,114 @@ int GCTAResponseCube::cache_index(const std::string& name) const } +#if defined(G_RADIAL_PSF_BASED) /***********************************************************************//** - * @brief Integrate PSF over radial model + * @brief Integrate Psf over radial model + * + * @param[in] model Radial model. + * @param[in] delta_mod 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. + * + * Integrates the product of the spatial model and the point spread + * function over the true photon arrival direction using + * + * \f[ + * \int_{\delta_{\rm min}}^{\delta_{\rm max}} + * \sin \delta \times {\rm Psf}(\delta) \times + * \int_{\phi_{\rm min}}^{\phi_{\rm max}} + * S_{\rm p}(\delta, \phi | E, t) d\phi d\delta + * \f] + * + * where + * \f$S_{\rm p}(\delta, \phi | E, t)\f$ is the radial spatial model, + * \f${\rm Psf}(\delta)\f$ is the point spread function, + * \f$\delta\f$ is angular distance between the true and the measured + * photon direction, and + * \f$\phi\f$ is the position angle around the observed photon direction + * measured counterclockwise from the connecting line between the model + * centre and the observed photon arrival direction. + ***************************************************************************/ +double GCTAResponseCube::psf_radial(const GModelSpatialRadial* model, + const double& delta_mod, + const GSkyDir& obsDir, + const GEnergy& srcEng, + const GTime& srcTime) const +{ + // Set number of iterations for Romberg integration. + // These values have been determined after careful testing, see + // https://cta-redmine.irap.omp.eu/issues/1291 + static const int iter_delta = 5; + static const int iter_phi = 6; + + // Initialise value + double value = 0.0; + + // Get maximum Psf radius (radians) + double psf_max = psf().delta_max(); + + // Get maximum model radius (radians) + double theta_max = model->theta_max(); + + // Set offset angle integration range (radians) + double delta_min = (delta_mod > theta_max) ? delta_mod - theta_max : 0.0; + double delta_max = delta_mod + theta_max; + if (delta_max > psf_max) { + delta_max = psf_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, + delta_mod, + theta_max, + iter_phi); + + // Integrate over model's zenith 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 of Psf within model and partial containment, then + // add a boundary at this location + double transition_point = theta_max - delta_mod; + if (transition_point > delta_min && transition_point < delta_max) { + bounds.push_back(transition_point); + } + + // Integrate kernel + value = integral.romberg(bounds, iter_delta); + + // Compile option: Check for NaN/Inf + #if defined(G_NAN_CHECK) + if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) { + std::cout << "*** ERROR: GCTAResponseCube::psf_radial:"; + std::cout << " NaN/Inf encountered"; + std::cout << " (value=" << value; + std::cout << ", delta_min=" << delta_min; + std::cout << ", delta_max=" << delta_max << ")"; + std::cout << std::endl; + } + #endif + + // Return integral + return value; +} +#else +/***********************************************************************//** + * @brief Integrate Psf over radial model * * @param[in] model Radial model. * @param[in] rho_obs Angle between model centre and measured photon direction @@ -932,21 +1028,23 @@ int GCTAResponseCube::cache_index(const std::string& name) const * @param[in] srcEng True photon energy. * @param[in] srcTime True photon arrival time. * - * Computes the integral + * Integrates the product of the spatial model and the point spread + * function over the true photon arrival direction using * * \f[ - * \int_0^{\delta_{\rm max}} - * {\rm PSF}(\delta) \times - * \int_0^{2\pi} {\rm M}(\delta, \phi) \sin \delta - * {\rm d}\phi {\rm d}\delta + * \int_{\rho_{\rm min}}^{\rho_{\rm max}} + * \sin \rho \times S_{\rm p}(\rho | E, t) \times + * \int_{\omega_{\rm min}}^{\omega_{\rm max}} + * {\rm Psf}(\rho, \omega) d\omega d\rho * \f] * - * where \f${\rm M}(\delta, \phi)\f$ is a 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. - * \f${\rm PSF}(\delta)\f$ is the azimuthally symmetric point spread - * function. + * where + * \f$S_{\rm p}(\rho | E, t)\f$ is the radial spatial model, + * \f${\rm Psf}(\rho, \omega)\f$ is the point spread function, + * \f$\rho\f$ is the radial distance from the model centre, and + * \f$\omega\f$ is the position angle around the model centre measured + * counterclockwise from the connecting line between the model centre and + * the observed photon arrival direction. ***************************************************************************/ double GCTAResponseCube::psf_radial(const GModelSpatialRadial* model, const double& rho_obs, @@ -1006,14 +1104,6 @@ double GCTAResponseCube::psf_radial(const GModelSpatialRadial* model, bounds.push_back(transition_point); } - // If the integration range includes a transition between full - // containment of Psf within model and partial containment, then - // add a boundary at this location - transition_point = src_max - delta_max; - if (transition_point > rho_min && transition_point < rho_max) { - 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 @@ -1045,10 +1135,11 @@ double GCTAResponseCube::psf_radial(const GModelSpatialRadial* model, // Return PSF return irf; } +#endif /***********************************************************************//** - * @brief Integrate PSF over elliptical model + * @brief Integrate Psf over elliptical model * * @param[in] model Elliptical model. * @param[in] rho_obs Angle between model centre and measured photon direction @@ -1059,36 +1150,24 @@ double GCTAResponseCube::psf_radial(const GModelSpatialRadial* model, * @param[in] srcEng True photon energy. * @param[in] srcTime True photon arrival time. * - * Computes the integral + * Integrates the product of the spatial model and the point spread + * function over the true photon arrival direction using * * \f[ * \int_{\rho_{\rm min}}^{\rho_{\rm max}} - * \sin \rho \times + * \sin \rho \times * \int_{\omega} - * S_{\rm p}(\rho, \omega | E, t) \, PSF(\rho, \omega) d\omega d\rho + * S_{\rm p}(\rho, \omega | E, t) \times + * {\rm Psf}(\rho, \omega) d\omega d\rho * \f] * * where - * \f$S_{\rm p}(\rho, \omega | E, t)\f$ is the elliptical model, - * \f$PSF(\rho, \omega)\f$ is the point spread function, - * \f$\rho\f$ is the distance from the model centre, and - * \f$\omega\f$ is the position angle with respect to the connecting line - * between the model centre and the observed photon arrival direction. - ***************************************************************************/ -/* - * \f[ - * \int_0^{\delta_{\rm max}} - * {\rm PSF}(\delta) \times - * \int_0^{2\pi} {\rm M}(\delta, \phi) \sin \delta - * {\rm d}\phi {\rm d}\delta - * \f] - * - * where \f${\rm M}(\delta, \phi)\f$ is an 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. - * \f${\rm PSF}(\delta)\f$ is the azimuthally symmetric point spread - * function. + * \f$S_{\rm p}(\rho, \omega | E, t)\f$ is the elliptical spatial model, + * \f${\rm Psf}(\rho, \omega)\f$ is the point spread function, + * \f$\rho\f$ is the radial distance from the model centre, and + * \f$\omega\f$ is the position angle around the model centre measured + * counterclockwise from the connecting line between the model centre and + * the observed photon arrival direction. ***************************************************************************/ double GCTAResponseCube::psf_elliptical(const GModelSpatialElliptical* model, const double& rho_obs, diff --git a/inst/cta/src/GCTAResponseIrf.cpp b/inst/cta/src/GCTAResponseIrf.cpp index c3f1be2..ffc2063 100644 --- a/inst/cta/src/GCTAResponseIrf.cpp +++ b/inst/cta/src/GCTAResponseIrf.cpp @@ -1268,7 +1268,7 @@ std::string GCTAResponseIrf::print(const GChatter& chatter) const /*========================================================================== = = - = Model type dependent CTA response methods = + = Model type dependent CTA response methods = = = ==========================================================================*/ @@ -1282,24 +1282,22 @@ 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 + * Integrates the product of the spatial model and the instrument response + * function 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 + * {\rm 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. + * \f$S_{\rm p}(\rho | E, t)\f$ is the radial spatial model, + * \f${\rm Irf}(\rho, \omega)\f$ is the instrument response function, + * \f$\rho\f$ is the radial distance from the model centre, and + * \f$\omega\f$ 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 @@ -1436,14 +1434,6 @@ double GCTAResponseIrf::irf_radial(const GEvent& event, bounds.push_back(transition_point); } - // If the integration range includes a transition between full - // containment of Psf within model and partial containment, then - // add a boundary at this location - transition_point = src_max - delta_max; - if (transition_point > rho_min && transition_point < rho_max) { - 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 @@ -1498,7 +1488,7 @@ double GCTAResponseIrf::irf_radial(const GEvent& event, /***********************************************************************//** - * @brief Return IRF value for elliptical source model + * @brief Return Irf value for elliptical source model * * @param[in] event Observed event. * @param[in] source Source. diff --git a/inst/cta/src/GCTAResponse_helpers.cpp b/inst/cta/src/GCTAResponse_helpers.cpp index a912a16..54486c0 100644 --- a/inst/cta/src/GCTAResponse_helpers.cpp +++ b/inst/cta/src/GCTAResponse_helpers.cpp @@ -692,7 +692,7 @@ double cta_irf_elliptical_kern_rho::eval(const double& rho) // If the radius rho is not larger than the semiminor axis // boundary, the circle with that radius is fully contained in // the ellipse and we can just integrate over the relevant arc - if (rho < m_semiminor) { + if (rho <= m_semiminor) { // Compute omega integration range double omega_min = -domega; @@ -984,7 +984,7 @@ double cta_npred_elliptical_kern_rho::eval(const double& rho) // If the radius rho is not larger than the semiminor axis // boundary, the circle with that radius is fully contained in // the ellipse and we can just integrate over the relevant arc - if (rho < m_semiminor) { + if (rho <= m_semiminor) { // Compute omega integration range double omega_min = -domega; @@ -1699,6 +1699,165 @@ double cta_psf_radial_kern_omega::eval(const double& omega) /***********************************************************************//** + * @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 + double psf = m_rsp->psf()(m_srcDir, delta, m_srcEng); + + // Continue only if Psf is positive + if (psf > 0.0) { + + // Compute half length of the arc (in radians) from a circle with + // radius delta that intersects with the model, defined as a circle + // with maximum radius m_theta_max + double dphi = 0.5 * gammalib::cta_roi_arclength(delta, + m_delta_mod, + m_cos_delta_mod, + m_sin_delta_mod, + m_theta_max, + m_cos_theta_max); + + // Continue only if arc length is positive + if (dphi > 0.0) { + + // Compute phi integration range + double phi_min = -dphi; + double phi_max = +dphi; + + // Precompute cosine and sine terms for azimuthal integration + double sin_delta = std::sin(delta); + double cos_delta = std::cos(delta); + double sin_fact = sin_delta * m_sin_delta_mod; + double cos_fact = cos_delta * m_cos_delta_mod; + + // 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); + + // Setup integrator + GIntegral integral(&integrand); + integral.fixed_iter(m_iter); + + // Integrate over azimuth + value = integral.romberg(phi_min, phi_max, m_iter) * + psf * sin_delta; + + // 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 + + } // endif: Psf value was positive + + } // endif: arc length was positive + + } // endif: delta was positive + + // 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 \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. + ***************************************************************************/ +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)); + + // Reduce theta by an infinite amount to avoid rounding errors at the + // boundary of a sharp edged model + double theta_kluge = theta - 1.0e-12; + if (theta_kluge < 0.0) { + theta_kluge = 0.0; + } + + // Get radial model value + double value = m_model->eval(theta_kluge, m_srcEng, m_srcTime); + + // Debug: test if model is non positive + #if defined(G_DEBUG_MODEL_ZERO) + if (value <= 0.0) { + std::cout << "*** WARNING: cta_psf_radial_kern_phi::eval"; + std::cout << " zero model for (phi)=("; + std::cout << phi*gammalib::rad2deg << ")"; + std::cout << " theta-r_model=" << (theta-m_model->theta_max()); + std::cout << " radians" << std::endl; + } + #endif + + // 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; +} + + +/***********************************************************************//** * @brief Kernel for elliptical model integration over zenith angle * * @param[in] rho Radial distance from model centre (radians). @@ -1777,7 +1936,7 @@ double cta_psf_elliptical_kern_rho::eval(const double& rho) // If the radius rho is not larger than the semiminor axis // boundary, the circle with that radius is fully contained in // the ellipse and we can just integrate over the relevant arc - if (rho < m_semiminor) { + if (rho <= m_semiminor) { // Compute omega integration range double omega_min = -domega; diff --git a/inst/cta/src/GCTAResponse_helpers.hpp b/inst/cta/src/GCTAResponse_helpers.hpp index 49d677f..6af453d 100644 --- a/inst/cta/src/GCTAResponse_helpers.hpp +++ b/inst/cta/src/GCTAResponse_helpers.hpp @@ -1091,6 +1091,75 @@ public: /***********************************************************************//** + * @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& delta_mod, + const double& theta_max, + const int& iter) : + m_rsp(rsp), + m_model(model), + m_srcDir(srcDir), + m_srcEng(srcEng), + m_srcTime(srcTime), + m_delta_mod(delta_mod), + m_cos_delta_mod(std::cos(delta_mod)), + m_sin_delta_mod(std::sin(delta_mod)), + m_theta_max(theta_max), + m_cos_theta_max(std::cos(theta_max)), + m_iter(iter) { } + 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_delta_mod; //!< Distance of model from Psf + double m_cos_delta_mod; //!< Cosine of m_delta_mod + double m_sin_delta_mod; //!< Sine of m_delta_mod + const double& m_theta_max; //!< Maximum model radius + double m_cos_theta_max; //!< Cosine of m_theta_max + const int& m_iter; //!< Integration iterations +}; + + +/***********************************************************************//** + * @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(delta_mod) + const double& m_cos_fact; //!< cos(delta)*cos(delta_mod) +}; + + +/***********************************************************************//** * @class cta_irf_elliptical_kern_rho * * @brief Kernel for elliptical model zenith angle integration diff --git a/inst/cta/test/test_irf_radial.py b/inst/cta/test/test_irf_radial.py index 62090a3..a581733 100755 --- a/inst/cta/test/test_irf_radial.py +++ b/inst/cta/test/test_irf_radial.py @@ -210,18 +210,19 @@ if __name__ == '__main__': sys.exit() # Set Psf size - delta_max = 0.6 + delta_max = 0.8 # Set model parameters centre = gammalib.GSkyDir() centre.radec_deg(0.0, 0.0) - radius = 0.8 + radius = 1.0 # Set observed photon direction obsDir = gammalib.GSkyDir() - obsDir.radec_deg(0.4, -0.3) + #obsDir.radec_deg(0.4, -0.3) #obsDir.radec_deg(-0.4, -0.3) #obsDir.radec_deg(0.4, 0.4) + obsDir.radec_deg(0.0, 0.3) # Show elliptical response geometry show_response(centre, radius, obsDir, delta_max) diff --git a/inst/cta/test/test_likelihood_profile.py b/inst/cta/test/test_likelihood_profile.py index 672290b..363a7b9 100755 --- a/inst/cta/test/test_likelihood_profile.py +++ b/inst/cta/test/test_likelihood_profile.py @@ -82,7 +82,7 @@ def show_likelihood_profile(result): # ========================== # # Compute likelihood profile # # ========================== # -def compute_likelihood_profile(obs, parname="RA", scale=1.0e-4, steps=10): +def compute_likelihood_profile(obs, parname="Radius", scale=1.0e-5, steps=10): """ Compute the log-likelihood profile for a given parameter. @@ -93,12 +93,12 @@ def compute_likelihood_profile(obs, parname="RA", scale=1.0e-4, steps=10): """ # Free parameter for which profile is to be determined obs.models()[0][parname].free() - obs.models()[0]["DEC"].free() # Allocate Levenberg-Marquardt optimizer log = gammalib.GLog() log.cout(True) opt = gammalib.GOptimizerLM(log) + opt.eps(5.0e-3) # Optimize model parameters obs.optimize(opt) @@ -256,8 +256,8 @@ if __name__ == '__main__': # Set test #test = "point" - test = "gauss" - #test = "disk" + #test = "gauss" + test = "disk" #test = "ldisk" #test = "shell" #test = "ellipse" @@ -309,10 +309,10 @@ if __name__ == '__main__': cntref = 337.5 # +/- 1.8 (from ctobssim) # Perform unbinned computation - original_likelihood_profile(model, events, irf, gammalib.GCaldb(caldb), cntref) + #original_likelihood_profile(model, events, irf, gammalib.GCaldb(caldb), cntref) # Perform binned computation #original_likelihood_profile(model, cntmap, irf, gammalib.GCaldb(caldb), cntref) # Perform stacked computation - #stacked_likelihood_profile(model, cntmap, expcube, psfcube, cntref) + stacked_likelihood_profile(model, cntmap, expcube, psfcube, cntref) diff --git a/inst/cta/test/test_npred_elliptical.py b/inst/cta/test/test_npred_elliptical.py index 4f35236..dc54bc8 100755 --- a/inst/cta/test/test_npred_elliptical.py +++ b/inst/cta/test/test_npred_elliptical.py @@ -279,7 +279,7 @@ def show_response(centre, semimajor, semiminor, posang, roiDir, pntDir, roi_radi if (domega > 0.0): # Full containment within ellipse - if (rho < bound_minor): + if (rho <= bound_minor): # Compute phi integration range omega_min = -domega @@ -400,17 +400,16 @@ if __name__ == '__main__': # Set model parameters centre = gammalib.GSkyDir() centre.radec_deg(0.0, 0.0) - #semimajor = 1.0 - #semiminor = 0.50 - semiminor = 1.0 - semimajor = 0.50 - posang = 310.0 + semimajor = 1.0 + semiminor = 0.50 + posang = 30.0 # Set ROI direction roiDir = gammalib.GSkyDir() - #roiDir.radec_deg(-0.4, -0.3) - roiDir.radec_deg(+0.4, -0.3) + roiDir.radec_deg(-0.4, -0.3) + #roiDir.radec_deg(+0.4, -0.3) #roiDir.radec_deg(+0.4, +0.3) + #roiDir.radec_deg(0.0, 0.1) # Set pointing direction pntDir = gammalib.GSkyDir() -- libgit2 0.21.4