Michael Mayer / gammalib

Commit fa87da8876ad9b2ff26ffd058b62e4fc715aa969

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

Replace fixed integration precision by fixed number of iterations.

Inline Side-by-side
Showing 10 changed files with 80 additions and 41 deletions
inst/cta/include/GCTAModelIrfBackground.hpp
  View file @fa87da8
... ... @@ -93,13 +93,16 @@ private:
93 93 class npred_roi_kern_theta : public GFunction {
94 94 public:
95 95 npred_roi_kern_theta(const GCTABackground* bgd,
96   - const double& logE) :
  96 + const double& logE,
  97 + const int& iter) :
97 98 m_bgd(bgd),
98   - m_logE(logE) { }
  99 + m_logE(logE),
  100 + m_iter(iter) { }
99 101 double eval(const double& theta);
100 102 protected:
101 103 const GCTABackground* m_bgd; //!< Pointer to background
102 104 const double& m_logE; //!< Log10 of energy
  105 + const int& m_iter; //!< Romberg iterations
103 106 };
104 107  
105 108 // ROI integration kernel over phi
... ...
inst/cta/src/GCTAAeffArf.cpp
  View file @fa87da8
... ... @@ -364,6 +364,9 @@ void GCTAAeffArf::read(const GFitsTable&amp; hdu)
364 364 ***************************************************************************/
365 365 void GCTAAeffArf::remove_thetacut(const GCTAResponseIrf& rsp)
366 366 {
  367 + // Set number of iterations for Romberg integration.
  368 + static const int iter = 6;
  369 +
367 370 // Continue only if thetacut value has been set
368 371 if (m_thetacut > 0.0) {
369 372  
... ... @@ -385,7 +388,9 @@ void GCTAAeffArf::remove_thetacut(const GCTAResponseIrf&amp; rsp)
385 388  
386 389 // Setup integration
387 390 GIntegral integral(&integrand);
388   - integral.eps(1.0e-5);
  391 +
  392 + // Set fixed number of iterations
  393 + integral.fixed_iter(iter);
389 394  
390 395 // Perform integration
391 396 double fraction = integral.romberg(0.0, rmax);
... ...
inst/cta/src/GCTACubeSourceDiffuse.cpp
  View file @fa87da8
... ... @@ -359,13 +359,9 @@ double GCTACubeSourceDiffuse::psf(const GCTAResponseCube* rsp,
359 359 const GEnergy& srcEng,
360 360 const GTime& srcTime) const
361 361 {
362   - // Set integration precision and Romberg order. These values have been
363   - // determined after careful testing, see
364   - // https://cta-redmine.irap.omp.eu/issues/1291
365   - const double eps_delta = 1.0e-2;
366   - const int order_delta = 5;
367   - const double eps_phi = 1.0e-2;
368   - const int order_phi = 5;
  362 + // Set number of iterations for Romberg integration.
  363 + static const int iter_delta = 5;
  364 + static const int iter_phi = 5;
369 365  
370 366 // Initialise PSF
371 367 double psf = 0.0;
... ... @@ -387,12 +383,16 @@ double GCTACubeSourceDiffuse::psf(const GCTAResponseCube* rsp,
387 383 // direction as the true photon arrival direction because the PSF does
388 384 // not vary significantly over a small region.
389 385 cta_psf_diffuse_kern_delta integrand(rsp, model, obsDir, srcEng, srcTime,
390   - rot, eps_phi, order_phi);
  386 + rot, iter_phi);
391 387  
392   - // Integrate over PSF delta angle
  388 + // Setup integration
393 389 GIntegral integral(&integrand);
394   - integral.eps(eps_delta);
395   - psf = integral.romberg(delta_min, delta_max, order_delta);
  390 +
  391 + // Set fixed number of iterations
  392 + integral.fixed_iter(iter_delta);
  393 +
  394 + // Integrate over PSF delta angle
  395 + psf = integral.romberg(delta_min, delta_max);
396 396  
397 397 // Return PSF
398 398 return psf;
... ...
inst/cta/src/GCTAModelIrfBackground.cpp
  View file @fa87da8
... ... @@ -65,8 +65,6 @@ const GModelRegistry g_cta_inst_background_registry(&amp;g_cta_inst_backgrou
65 65 //#define G_DEBUG_NPRED
66 66  
67 67 /* __ Constants __________________________________________________________ */
68   -const double g_cta_inst_background_npred_theta_eps = 1.0e-4;
69   -const double g_cta_inst_background_npred_phi_eps = 1.0e-4;
70 68  
71 69  
72 70 /*==========================================================================
... ... @@ -448,6 +446,10 @@ double GCTAModelIrfBackground::npred(const GEnergy&amp; obsEng,
448 446 const GTime& obsTime,
449 447 const GObservation& obs) const
450 448 {
  449 + // Set number of iterations for Romberg integration.
  450 + static const int iter_theta = 6;
  451 + static const int iter_phi = 6;
  452 +
451 453 // Initialise result
452 454 double npred = 0.0;
453 455 bool has_npred = false;
... ... @@ -525,11 +527,15 @@ double GCTAModelIrfBackground::npred(const GEnergy&amp; obsEng,
525 527 double logE = obsEng.log10TeV();
526 528  
527 529 // Setup integration function
528   - GCTAModelIrfBackground::npred_roi_kern_theta integrand(bgd, logE);
  530 + GCTAModelIrfBackground::npred_roi_kern_theta integrand(bgd,
  531 + logE,
  532 + iter_phi);
529 533  
530   - // Setup integrator
  534 + // Setup integration
531 535 GIntegral integral(&integrand);
532   - integral.eps(g_cta_inst_background_npred_theta_eps);
  536 +
  537 + // Set fixed number of iterations
  538 + integral.fixed_iter(iter_theta);
533 539  
534 540 // Spatially integrate radial component
535 541 npred = integral.romberg(0.0, roi_radius);
... ... @@ -1208,9 +1214,13 @@ double GCTAModelIrfBackground::npred_roi_kern_theta::eval(const double&amp; theta)
1208 1214 m_logE,
1209 1215 theta);
1210 1216  
1211   - // Integrate over phi
  1217 + // Setup integration
1212 1218 GIntegral integral(&integrand);
1213   - integral.eps(g_cta_inst_background_npred_phi_eps);
  1219 +
  1220 + // Set fixed number of iterations
  1221 + integral.fixed_iter(m_iter);
  1222 +
  1223 + // Integrate over phi
1214 1224 value = integral.romberg(0.0, gammalib::twopi) * std::sin(theta);
1215 1225  
1216 1226 // Debug: Check for NaN
... ...
inst/cta/src/GCTAResponseIrf.cpp
  View file @fa87da8
... ... @@ -470,6 +470,9 @@ double GCTAResponseIrf::nroi(const GModelSky&amp; model,
470 470 const GTime& obsTime,
471 471 const GObservation& obs) const
472 472 {
  473 + // Set number of iterations for Romberg integration.
  474 + static const int iter = 4;
  475 +
473 476 // Initialise Nroi value
474 477 double nroi = 0.0;
475 478  
... ... @@ -498,8 +501,8 @@ double GCTAResponseIrf::nroi(const GModelSky&amp; model,
498 501 cta_nroi_kern integrand(model, *this, srcTime, obsEng, obsTime, obs);
499 502 GIntegral integral(&integrand);
500 503  
501   - // Set integration precision
502   - integral.eps(1.0e-3);
  504 + // Set fixed number of iterations
  505 + integral.fixed_iter(iter);
503 506  
504 507 // Do Romberg integration
505 508 emin = std::log(emin);
... ... @@ -1811,6 +1814,9 @@ double GCTAResponseIrf::npsf(const GSkyDir&amp; srcDir,
1811 1814 const GCTAPointing& pnt,
1812 1815 const GCTARoi& roi) const
1813 1816 {
  1817 + // Set number of iterations for Romberg integration.
  1818 + static const int iter = 6;
  1819 +
1814 1820 // Declare result
1815 1821 double value = 0.0;
1816 1822  
... ... @@ -1857,7 +1863,9 @@ double GCTAResponseIrf::npsf(const GSkyDir&amp; srcDir,
1857 1863  
1858 1864 // Setup integration
1859 1865 GIntegral integral(&integrand);
1860   - integral.eps(1.0e-5);
  1866 +
  1867 + // Set fixed number of iterations
  1868 + integral.fixed_iter(iter);
1861 1869  
1862 1870 // Radially integrate PSF. In case that the radial integration
1863 1871 // region is small, we do the integration using a simple
... ...
inst/cta/src/GCTAResponse_helpers.cpp
  View file @fa87da8
... ... @@ -2185,10 +2185,14 @@ double cta_psf_diffuse_kern_delta::eval(const double&amp; delta)
2185 2185 m_rot,
2186 2186 sin_delta, cos_delta);
2187 2187  
2188   - // Azimuthally integrate model
  2188 + // Setup untegration
2189 2189 GIntegral integral(&integrand);
2190   - integral.eps(m_eps);
2191   - value *= integral.romberg(0.0, gammalib::twopi, m_order) * sin_delta;
  2190 +
  2191 + // Set fixed number of iterations
  2192 + integral.fixed_iter(m_iter);
  2193 +
  2194 + // Azimuthally integrate model
  2195 + value *= integral.romberg(0.0, gammalib::twopi) * sin_delta;
2192 2196  
2193 2197 } // endif: PSF value was positive
2194 2198  
... ...
inst/cta/src/GCTAResponse_helpers.hpp
  View file @fa87da8
... ... @@ -1375,16 +1375,14 @@ public:
1375 1375 const GEnergy& srcEng,
1376 1376 const GTime& srcTime,
1377 1377 const GMatrix& rot,
1378   - const double& eps,
1379   - const int& order) :
  1378 + const int& iter) :
1380 1379 m_rsp(rsp),
1381 1380 m_model(model),
1382 1381 m_srcDir(srcDir),
1383 1382 m_srcEng(srcEng),
1384 1383 m_srcTime(srcTime),
1385 1384 m_rot(rot),
1386   - m_eps(eps),
1387   - m_order (order),
  1385 + m_iter (iter),
1388 1386 m_psf_max(rsp->psf()(srcDir, 0.0, srcEng)) { }
1389 1387 double eval(const double& delta);
1390 1388 protected:
... ... @@ -1394,8 +1392,7 @@ protected:
1394 1392 const GEnergy& m_srcEng; //!< True photon energy
1395 1393 const GTime& m_srcTime; //!< True photon arrival time
1396 1394 const GMatrix& m_rot; //!< Rotation matrix
1397   - const double& m_eps; //!< Integration precision
1398   - const int& m_order; //!< Romberg order
  1395 + const int& m_iter; //!< Romberg iterations
1399 1396 double m_psf_max; //!< Maximum PSF value
1400 1397 };
1401 1398  
... ...
src/model/GModelSpectralLogParabola.cpp
  View file @fa87da8
... ... @@ -408,6 +408,9 @@ double GModelSpectralLogParabola::eval_gradients(const GEnergy&amp; srcEng,
408 408 double GModelSpectralLogParabola::flux(const GEnergy& emin,
409 409 const GEnergy& emax) const
410 410 {
  411 + // Set number of iterations for Romberg integration.
  412 + static const int iter = 9;
  413 +
411 414 // Initialise flux
412 415 double flux = 0.0;
413 416  
... ... @@ -420,8 +423,8 @@ double GModelSpectralLogParabola::flux(const GEnergy&amp; emin,
420 423 // Initialise integral class with function
421 424 GIntegral integral(&kernel);
422 425  
423   - // Set integration precision
424   - integral.eps(1.0e-8);
  426 + // Set number of iterations
  427 + integral.fixed_iter(iter);
425 428  
426 429 // Calculate integral between emin and emax
427 430 flux = integral.romberg(emin.MeV(), emax.MeV());
... ... @@ -454,6 +457,9 @@ double GModelSpectralLogParabola::flux(const GEnergy&amp; emin,
454 457 double GModelSpectralLogParabola::eflux(const GEnergy& emin,
455 458 const GEnergy& emax) const
456 459 {
  460 + // Set number of iterations for Romberg integration.
  461 + static const int iter = 9;
  462 +
457 463 // Initialise flux
458 464 double eflux = 0.0;
459 465  
... ... @@ -466,8 +472,8 @@ double GModelSpectralLogParabola::eflux(const GEnergy&amp; emin,
466 472 // Initialise integral class with function
467 473 GIntegral integral(&kernel);
468 474  
469   - // Set integration precision
470   - integral.eps(1.0e-8);
  475 + // Set number of iterations
  476 + integral.fixed_iter(iter);
471 477  
472 478 // Calculate integral between emin and emax
473 479 eflux = integral.romberg(emin.MeV(), emax.MeV());
... ...
src/obs/GObservation.cpp
  View file @fa87da8
... ... @@ -1409,6 +1409,9 @@ double GObservation::npred_kern::eval(const double&amp; x)
1409 1409 double GObservation::npred_spec(const GModel& model,
1410 1410 const GTime& obsTime) const
1411 1411 {
  1412 + // Set number of iterations for Romberg integration.
  1413 + static const int iter = 8;
  1414 +
1412 1415 // Initialise result
1413 1416 double result = 0.0;
1414 1417  
... ... @@ -1429,8 +1432,8 @@ double GObservation::npred_spec(const GModel&amp; model,
1429 1432 GObservation::npred_spec_kern integrand(this, &model, &obsTime);
1430 1433 GIntegral integral(&integrand);
1431 1434  
1432   - // Set integration precision
1433   - integral.eps(1.0e-6); // Needed for fluctuating bgd. model !!!
  1435 + // Set number of iterations
  1436 + integral.fixed_iter(iter);
1434 1437  
1435 1438 // Do Romberg integration
1436 1439 #if defined(G_LN_ENERGY_INT)
... ...
src/obs/GResponse.cpp
  View file @fa87da8
... ... @@ -174,6 +174,9 @@ double GResponse::convolve(const GModelSky&amp; model,
174 174 const GObservation& obs,
175 175 const bool& grad) const
176 176 {
  177 + // Set number of iterations for Romberg integration.
  178 + static const int iter = 4;
  179 +
177 180 // Initialise result
178 181 double prob = 0.0;
179 182  
... ... @@ -203,8 +206,8 @@ double GResponse::convolve(const GModelSky&amp; model,
203 206 edisp_kern integrand(this, model, event, srcTime, obs, grad);
204 207 GIntegral integral(&integrand);
205 208  
206   - // Set integration precision
207   - integral.eps(1.0e-3);
  209 + // Set number of iterations
  210 + integral.fixed_iter(iter);
208 211  
209 212 // Do Romberg integration
210 213 emin = std::log(emin);
... ...