diff --git a/inst/cta/include/GCTAModelIrfBackground.hpp b/inst/cta/include/GCTAModelIrfBackground.hpp index fd185f7..dce660c 100644 --- a/inst/cta/include/GCTAModelIrfBackground.hpp +++ b/inst/cta/include/GCTAModelIrfBackground.hpp @@ -93,13 +93,16 @@ private: class npred_roi_kern_theta : public GFunction { public: npred_roi_kern_theta(const GCTABackground* bgd, - const double& logE) : + const double& logE, + const int& iter) : m_bgd(bgd), - m_logE(logE) { } + m_logE(logE), + m_iter(iter) { } double eval(const double& theta); protected: const GCTABackground* m_bgd; //!< Pointer to background const double& m_logE; //!< Log10 of energy + const int& m_iter; //!< Romberg iterations }; // ROI integration kernel over phi diff --git a/inst/cta/src/GCTAAeffArf.cpp b/inst/cta/src/GCTAAeffArf.cpp index f16ae4f..013d11f 100644 --- a/inst/cta/src/GCTAAeffArf.cpp +++ b/inst/cta/src/GCTAAeffArf.cpp @@ -364,6 +364,9 @@ void GCTAAeffArf::read(const GFitsTable& hdu) ***************************************************************************/ void GCTAAeffArf::remove_thetacut(const GCTAResponseIrf& rsp) { + // Set number of iterations for Romberg integration. + static const int iter = 6; + // Continue only if thetacut value has been set if (m_thetacut > 0.0) { @@ -385,7 +388,9 @@ void GCTAAeffArf::remove_thetacut(const GCTAResponseIrf& rsp) // Setup integration GIntegral integral(&integrand); - integral.eps(1.0e-5); + + // Set fixed number of iterations + integral.fixed_iter(iter); // Perform integration double fraction = integral.romberg(0.0, rmax); diff --git a/inst/cta/src/GCTACubeSourceDiffuse.cpp b/inst/cta/src/GCTACubeSourceDiffuse.cpp index da22233..cadceb0 100644 --- a/inst/cta/src/GCTACubeSourceDiffuse.cpp +++ b/inst/cta/src/GCTACubeSourceDiffuse.cpp @@ -359,13 +359,9 @@ double GCTACubeSourceDiffuse::psf(const GCTAResponseCube* rsp, 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; + // Set number of iterations for Romberg integration. + static const int iter_delta = 5; + static const int iter_phi = 5; // Initialise PSF double psf = 0.0; @@ -387,12 +383,16 @@ double GCTACubeSourceDiffuse::psf(const GCTAResponseCube* rsp, // direction as the true photon arrival direction because the PSF does // not vary significantly over a small region. cta_psf_diffuse_kern_delta integrand(rsp, model, obsDir, srcEng, srcTime, - rot, eps_phi, order_phi); + rot, iter_phi); - // Integrate over PSF delta angle + // Setup integration GIntegral integral(&integrand); - integral.eps(eps_delta); - psf = integral.romberg(delta_min, delta_max, order_delta); + + // Set fixed number of iterations + integral.fixed_iter(iter_delta); + + // Integrate over PSF delta angle + psf = integral.romberg(delta_min, delta_max); // Return PSF return psf; diff --git a/inst/cta/src/GCTAModelIrfBackground.cpp b/inst/cta/src/GCTAModelIrfBackground.cpp index e1e0f75..1ad96ec 100644 --- a/inst/cta/src/GCTAModelIrfBackground.cpp +++ b/inst/cta/src/GCTAModelIrfBackground.cpp @@ -65,8 +65,6 @@ const GModelRegistry g_cta_inst_background_registry(&g_cta_inst_backgrou //#define G_DEBUG_NPRED /* __ Constants __________________________________________________________ */ -const double g_cta_inst_background_npred_theta_eps = 1.0e-4; -const double g_cta_inst_background_npred_phi_eps = 1.0e-4; /*========================================================================== @@ -448,6 +446,10 @@ double GCTAModelIrfBackground::npred(const GEnergy& obsEng, const GTime& obsTime, const GObservation& obs) const { + // Set number of iterations for Romberg integration. + static const int iter_theta = 6; + static const int iter_phi = 6; + // Initialise result double npred = 0.0; bool has_npred = false; @@ -525,11 +527,15 @@ double GCTAModelIrfBackground::npred(const GEnergy& obsEng, double logE = obsEng.log10TeV(); // Setup integration function - GCTAModelIrfBackground::npred_roi_kern_theta integrand(bgd, logE); + GCTAModelIrfBackground::npred_roi_kern_theta integrand(bgd, + logE, + iter_phi); - // Setup integrator + // Setup integration GIntegral integral(&integrand); - integral.eps(g_cta_inst_background_npred_theta_eps); + + // Set fixed number of iterations + integral.fixed_iter(iter_theta); // Spatially integrate radial component npred = integral.romberg(0.0, roi_radius); @@ -1208,9 +1214,13 @@ double GCTAModelIrfBackground::npred_roi_kern_theta::eval(const double& theta) m_logE, theta); - // Integrate over phi + // Setup integration GIntegral integral(&integrand); - integral.eps(g_cta_inst_background_npred_phi_eps); + + // Set fixed number of iterations + integral.fixed_iter(m_iter); + + // Integrate over phi value = integral.romberg(0.0, gammalib::twopi) * std::sin(theta); // Debug: Check for NaN diff --git a/inst/cta/src/GCTAResponseIrf.cpp b/inst/cta/src/GCTAResponseIrf.cpp index 6255ac7..8c6b19d 100644 --- a/inst/cta/src/GCTAResponseIrf.cpp +++ b/inst/cta/src/GCTAResponseIrf.cpp @@ -470,6 +470,9 @@ double GCTAResponseIrf::nroi(const GModelSky& model, const GTime& obsTime, const GObservation& obs) const { + // Set number of iterations for Romberg integration. + static const int iter = 4; + // Initialise Nroi value double nroi = 0.0; @@ -498,8 +501,8 @@ double GCTAResponseIrf::nroi(const GModelSky& model, cta_nroi_kern integrand(model, *this, srcTime, obsEng, obsTime, obs); GIntegral integral(&integrand); - // Set integration precision - integral.eps(1.0e-3); + // Set fixed number of iterations + integral.fixed_iter(iter); // Do Romberg integration emin = std::log(emin); @@ -1811,6 +1814,9 @@ double GCTAResponseIrf::npsf(const GSkyDir& srcDir, const GCTAPointing& pnt, const GCTARoi& roi) const { + // Set number of iterations for Romberg integration. + static const int iter = 6; + // Declare result double value = 0.0; @@ -1857,7 +1863,9 @@ double GCTAResponseIrf::npsf(const GSkyDir& srcDir, // Setup integration GIntegral integral(&integrand); - integral.eps(1.0e-5); + + // Set fixed number of iterations + integral.fixed_iter(iter); // Radially integrate PSF. In case that the radial integration // region is small, we do the integration using a simple diff --git a/inst/cta/src/GCTAResponse_helpers.cpp b/inst/cta/src/GCTAResponse_helpers.cpp index b9793d8..1fa0af2 100644 --- a/inst/cta/src/GCTAResponse_helpers.cpp +++ b/inst/cta/src/GCTAResponse_helpers.cpp @@ -2185,10 +2185,14 @@ double cta_psf_diffuse_kern_delta::eval(const double& delta) m_rot, sin_delta, cos_delta); - // Azimuthally integrate model + // Setup untegration GIntegral integral(&integrand); - integral.eps(m_eps); - value *= integral.romberg(0.0, gammalib::twopi, m_order) * sin_delta; + + // Set fixed number of iterations + integral.fixed_iter(m_iter); + + // Azimuthally integrate model + value *= integral.romberg(0.0, gammalib::twopi) * sin_delta; } // endif: PSF value was positive diff --git a/inst/cta/src/GCTAResponse_helpers.hpp b/inst/cta/src/GCTAResponse_helpers.hpp index c1dc147..30d48f5 100644 --- a/inst/cta/src/GCTAResponse_helpers.hpp +++ b/inst/cta/src/GCTAResponse_helpers.hpp @@ -1375,16 +1375,14 @@ public: const GEnergy& srcEng, const GTime& srcTime, const GMatrix& rot, - const double& eps, - const int& order) : + const int& iter) : 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_iter (iter), m_psf_max(rsp->psf()(srcDir, 0.0, srcEng)) { } double eval(const double& delta); protected: @@ -1394,8 +1392,7 @@ protected: 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 + const int& m_iter; //!< Romberg iterations double m_psf_max; //!< Maximum PSF value }; diff --git a/src/model/GModelSpectralLogParabola.cpp b/src/model/GModelSpectralLogParabola.cpp index 4407479..8884881 100644 --- a/src/model/GModelSpectralLogParabola.cpp +++ b/src/model/GModelSpectralLogParabola.cpp @@ -408,6 +408,9 @@ double GModelSpectralLogParabola::eval_gradients(const GEnergy& srcEng, double GModelSpectralLogParabola::flux(const GEnergy& emin, const GEnergy& emax) const { + // Set number of iterations for Romberg integration. + static const int iter = 9; + // Initialise flux double flux = 0.0; @@ -420,8 +423,8 @@ double GModelSpectralLogParabola::flux(const GEnergy& emin, // Initialise integral class with function GIntegral integral(&kernel); - // Set integration precision - integral.eps(1.0e-8); + // Set number of iterations + integral.fixed_iter(iter); // Calculate integral between emin and emax flux = integral.romberg(emin.MeV(), emax.MeV()); @@ -454,6 +457,9 @@ double GModelSpectralLogParabola::flux(const GEnergy& emin, double GModelSpectralLogParabola::eflux(const GEnergy& emin, const GEnergy& emax) const { + // Set number of iterations for Romberg integration. + static const int iter = 9; + // Initialise flux double eflux = 0.0; @@ -466,8 +472,8 @@ double GModelSpectralLogParabola::eflux(const GEnergy& emin, // Initialise integral class with function GIntegral integral(&kernel); - // Set integration precision - integral.eps(1.0e-8); + // Set number of iterations + integral.fixed_iter(iter); // Calculate integral between emin and emax eflux = integral.romberg(emin.MeV(), emax.MeV()); diff --git a/src/obs/GObservation.cpp b/src/obs/GObservation.cpp index 704f9c5..8fadc41 100644 --- a/src/obs/GObservation.cpp +++ b/src/obs/GObservation.cpp @@ -1409,6 +1409,9 @@ double GObservation::npred_kern::eval(const double& x) double GObservation::npred_spec(const GModel& model, const GTime& obsTime) const { + // Set number of iterations for Romberg integration. + static const int iter = 8; + // Initialise result double result = 0.0; @@ -1429,8 +1432,8 @@ double GObservation::npred_spec(const GModel& model, GObservation::npred_spec_kern integrand(this, &model, &obsTime); GIntegral integral(&integrand); - // Set integration precision - integral.eps(1.0e-6); // Needed for fluctuating bgd. model !!! + // Set number of iterations + integral.fixed_iter(iter); // Do Romberg integration #if defined(G_LN_ENERGY_INT) diff --git a/src/obs/GResponse.cpp b/src/obs/GResponse.cpp index a040c04..d6e5bdd 100644 --- a/src/obs/GResponse.cpp +++ b/src/obs/GResponse.cpp @@ -174,6 +174,9 @@ double GResponse::convolve(const GModelSky& model, const GObservation& obs, const bool& grad) const { + // Set number of iterations for Romberg integration. + static const int iter = 4; + // Initialise result double prob = 0.0; @@ -203,8 +206,8 @@ double GResponse::convolve(const GModelSky& model, edisp_kern integrand(this, model, event, srcTime, obs, grad); GIntegral integral(&integrand); - // Set integration precision - integral.eps(1.0e-3); + // Set number of iterations + integral.fixed_iter(iter); // Do Romberg integration emin = std::log(emin);