Justine Devin / gammalib

Commit 6c547174abc39f740b7859e337f889f45ac5b67c

Authored by Jürgen Knödlseder
1 parent 8a054f18

Put back original integration precision of 1e-4 and remove double deadtime correction. 

Deadtime correction was applied in npred() and npred_radial(), npred_elliptical() and npred_diffuse(). The latter three methods all call npred() at the end, which led to a double correction of the deadtime. This has been corrected now.
Inline Side-by-side
Showing 1 changed file with 14 additions and 13 deletions
inst/cta/src/GCTAResponse.cpp
  View file @6c54717
... ... @@ -64,8 +64,7 @@
64 64 " GEnergy&, GTime&, GObservation&)"
65 65 #define G_NPRED "GCTAResponse::npred(GSkyDir&, GEnergy&, GTime&,"\
66 66 " GObservation&)"
67   -#define G_MC "GCTAResponse::mc(double&,GPhoton&,GPointing&,GRan&)"
68   -
  67 +#define G_MC "GCTAResponse::mc(double&, GPhoton&, GObservation&, GRan&)"
69 68 #define G_IRF_RADIAL "GCTAResponse::irf_radial(GEvent&, GSource&,"\
70 69 " GObservation&)"
71 70 #define G_IRF_ELLIPTICAL "GCTAResponse::irf_elliptical(GEvent&, GSource&,"\
... ... @@ -471,6 +470,18 @@ GCTAEventAtom* GCTAResponse::mc(const double& area, const GPhoton& photon,
471 470 // Compute limiting value
472 471 double ulimite = effective_area / area;
473 472  
  473 + // Warning if ulimite is larger than one
  474 + if (ulimite > 1.0) {
  475 + std::string msg = "Effective area "+
  476 + gammalib::str(effective_area)+
  477 + " cm2 is larger than simulation surface area "+
  478 + gammalib::str(area)+
  479 + " cm2 for photon energy "+
  480 + gammalib::str(photon.energy().TeV())+
  481 + " TeV. Simulations are inaccurate.";
  482 + gammalib::warning(G_MC, msg);
  483 + }
  484 +
474 485 // Continue only if event is detected
475 486 if (ran.uniform() <= ulimite) {
476 487  
... ... @@ -1487,9 +1498,6 @@ double GCTAResponse::npred_radial(const GSource&amp; source,
1487 1498 GIntegral integral(&integrand);
1488 1499 npred = integral.romb(rho_min, rho_max);
1489 1500  
1490   - // Apply deadtime correction
1491   - npred *= obs.deadc(srcTime);
1492   -
1493 1501 // Compile option: Show integration results
1494 1502 #if defined(G_DEBUG_NPRED_RADIAL)
1495 1503 std::cout << "GCTAResponse::npred_radial:";
... ... @@ -1650,9 +1658,6 @@ double GCTAResponse::npred_elliptical(const GSource&amp; source,
1650 1658 GIntegral integral(&integrand);
1651 1659 npred = integral.romb(rho_min, rho_max);
1652 1660  
1653   - // Apply deadtime correction
1654   - npred *= obs.deadc(srcTime);
1655   -
1656 1661 // Compile option: Show integration results
1657 1662 #if defined(G_DEBUG_NPRED_ELLIPTICAL)
1658 1663 std::cout << "GCTAResponse::npred_elliptical:";
... ... @@ -1817,13 +1822,9 @@ double GCTAResponse::npred_diffuse(const GSource&amp; source,
1817 1822  
1818 1823 // Integrate over theta
1819 1824 GIntegral integral(&integrand);
1820   - //integral.eps(1.0e-4);
1821   - integral.eps(1.0e-5);
  1825 + integral.eps(1.0e-4);
1822 1826 npred = integral.romb(0.0, roi_psf_radius);
1823 1827  
1824   - // Apply deadtime correction
1825   - npred *= obs.deadc(srcTime);
1826   -
1827 1828 // Compile option: Show integration results
1828 1829 #if defined(G_DEBUG_NPRED_DIFFUSE)
1829 1830 std::cout << "GCTAResponse::npred_diffuse:";
... ...