Johan Bregeon / gammalib

Commit edad52f329b8bf86541a501684ed490cfe5fd4de

Authored by Jürgen Knödlseder
1 parent 6c547174

Fix Npred integration precisions to 1e-5. Do not use m_eps anymore. 

Fix all Npred integration precision to 1e-5 as this is needed for accurate results. Do not use m_eps for the integration precision anymore as each integrator may be optimized individually.
Inline Side-by-side
Showing 2 changed files with 9 additions and 5 deletions
inst/cta/src/GCTAResponse.cpp
  View file @edad52f
@@ -985,7 +985,7 @@ double GCTAResponse::irf_radial(const GEvent& event, @@ -985,7 +985,7 @@ double GCTAResponse::irf_radial(const GEvent& event,
985 985
986 // Integrate over zenith angle 986 // Integrate over zenith angle
987 GIntegral integral(&integrand); 987 GIntegral integral(&integrand);
988 - integral.eps(m_eps); 988 + integral.eps(1.0e-5);
989 irf = integral.romb(rho_min, rho_max); 989 irf = integral.romb(rho_min, rho_max);
990 990
991 // Compile option: Check for NaN/Inf 991 // Compile option: Check for NaN/Inf
@@ -1153,7 +1153,7 @@ double GCTAResponse::irf_elliptical(const GEvent& event, @@ -1153,7 +1153,7 @@ double GCTAResponse::irf_elliptical(const GEvent& event,
1153 1153
1154 // Integrate over zenith angle 1154 // Integrate over zenith angle
1155 GIntegral integral(&integrand); 1155 GIntegral integral(&integrand);
1156 - integral.eps(m_eps); 1156 + integral.eps(1.0e-5);
1157 irf = integral.romb(rho_min, rho_max); 1157 irf = integral.romb(rho_min, rho_max);
1158 1158
1159 // Compile option: Check for NaN/Inf 1159 // Compile option: Check for NaN/Inf
@@ -1496,6 +1496,7 @@ double GCTAResponse::npred_radial(const GSource& source, @@ -1496,6 +1496,7 @@ double GCTAResponse::npred_radial(const GSource& source,
1496 1496
1497 // Integrate over theta 1497 // Integrate over theta
1498 GIntegral integral(&integrand); 1498 GIntegral integral(&integrand);
  1499 + integral.eps(1.0e-5);
1499 npred = integral.romb(rho_min, rho_max); 1500 npred = integral.romb(rho_min, rho_max);
1500 1501
1501 // Compile option: Show integration results 1502 // Compile option: Show integration results
@@ -1656,6 +1657,7 @@ double GCTAResponse::npred_elliptical(const GSource& source, @@ -1656,6 +1657,7 @@ double GCTAResponse::npred_elliptical(const GSource& source,
1656 1657
1657 // Integrate over theta 1658 // Integrate over theta
1658 GIntegral integral(&integrand); 1659 GIntegral integral(&integrand);
  1660 + integral.eps(1.0e-5);
1659 npred = integral.romb(rho_min, rho_max); 1661 npred = integral.romb(rho_min, rho_max);
1660 1662
1661 // Compile option: Show integration results 1663 // Compile option: Show integration results
@@ -1822,7 +1824,7 @@ double GCTAResponse::npred_diffuse(const GSource& source, @@ -1822,7 +1824,7 @@ double GCTAResponse::npred_diffuse(const GSource& source,
1822 1824
1823 // Integrate over theta 1825 // Integrate over theta
1824 GIntegral integral(&integrand); 1826 GIntegral integral(&integrand);
1825 - integral.eps(1.0e-4); 1827 + integral.eps(1.0e-5);
1826 npred = integral.romb(0.0, roi_psf_radius); 1828 npred = integral.romb(0.0, roi_psf_radius);
1827 1829
1828 // Compile option: Show integration results 1830 // Compile option: Show integration results
@@ -2068,7 +2070,7 @@ double GCTAResponse::npsf(const GSkyDir& srcDir, @@ -2068,7 +2070,7 @@ double GCTAResponse::npsf(const GSkyDir& srcDir,
2068 2070
2069 // Setup integration 2071 // Setup integration
2070 GIntegral integral(&integrand); 2072 GIntegral integral(&integrand);
2071 - integral.eps(m_eps); 2073 + integral.eps(1.0e-5);
2072 2074
2073 // Radially integrate PSF. In case that the radial integration 2075 // Radially integrate PSF. In case that the radial integration
2074 // region is small, we do the integration using a simple 2076 // region is small, we do the integration using a simple
inst/cta/src/GCTAResponse_helpers.cpp
  View file @edad52f
@@ -339,6 +339,7 @@ double cta_npred_radial_kern_rho::eval(const double& rho) @@ -339,6 +339,7 @@ double cta_npred_radial_kern_rho::eval(const double& rho)
339 339
340 // Integrate over phi 340 // Integrate over phi
341 GIntegral integral(&integrand); 341 GIntegral integral(&integrand);
  342 + integral.eps(1.0e-5);
342 npred = integral.romb(omega_min, omega_max) * sin_rho * model; 343 npred = integral.romb(omega_min, omega_max) * sin_rho * model;
343 344
344 // Debug: Check for NaN 345 // Debug: Check for NaN
@@ -643,6 +644,7 @@ double cta_npred_elliptical_kern_rho::eval(const double& rho) @@ -643,6 +644,7 @@ double cta_npred_elliptical_kern_rho::eval(const double& rho)
643 644
644 // Integrate over phi 645 // Integrate over phi
645 GIntegral integral(&integrand); 646 GIntegral integral(&integrand);
  647 + integral.eps(1.0e-5);
646 npred = integral.romb(omega_min, omega_max) * sin_rho; 648 npred = integral.romb(omega_min, omega_max) * sin_rho;
647 649
648 // Debug: Check for NaN 650 // Debug: Check for NaN
@@ -966,7 +968,7 @@ double cta_npred_diffuse_kern_theta::eval(const double& theta) @@ -966,7 +968,7 @@ double cta_npred_diffuse_kern_theta::eval(const double& theta)
966 968
967 // Integrate over phi 969 // Integrate over phi
968 GIntegral integral(&integrand); 970 GIntegral integral(&integrand);
969 - integral.eps(1.0e-4); 971 + integral.eps(1.0e-5);
970 npred = integral.romb(0.0, gammalib::twopi) * sin_theta; 972 npred = integral.romb(0.0, gammalib::twopi) * sin_theta;
971 973
972 // Debug: Check for NaN 974 // Debug: Check for NaN