Simon Bonnefoy / gammalib

Commit e5c337473c9557302d85ed51d0e126aaa766f302

Authored by Jürgen Knödlseder
1 parent f2695a75

Adapt to new GCTAEdisp2D interface. Do not use diffuse value cache in case of energy dispersion.

Inline Side-by-side
Showing 2 changed files with 76 additions and 48 deletions
inst/cta/src/GCTAResponseIrf.cpp
  View file @e5c3374
... ... @@ -115,6 +115,7 @@
115 115 //#define G_DEBUG_IRF_RADIAL //!< Debug irf_radial method
116 116 //#define G_DEBUG_IRF_DIFFUSE //!< Debug irf_diffuse method
117 117 //#define G_DEBUG_IRF_ELLIPTICAL //!< Debug irf_elliptical method
  118 +//#define G_DEBUG_NROI_EDISP //!< Debug nroi energy dispersion computation
118 119 //#define G_DEBUG_NROI_RADIAL //!< Debug npred_radial method
119 120 //#define G_DEBUG_NROI_ELLIPTICAL //!< Debug npred_elliptical method
120 121 //#define G_DEBUG_NROI_DIFFUSE //!< Debug npred_diffuse method
... ... @@ -455,12 +456,12 @@ double GCTAResponseIrf::nroi(const GModelSky&amp; model,
455 456 // Loop over all boundaries
456 457 for (int i = 0; i < ebounds.size(); ++i) {
457 458  
458   - // Get boundaries in MeV
459   - double emin = ebounds.emin(i).MeV();
460   - double emax = ebounds.emax(i).MeV();
  459 + // Get true energy boundaries in log10 MeV
  460 + double log10_etrue_min = ebounds.emin(i).log10MeV();
  461 + double log10_etrue_max = ebounds.emax(i).log10MeV();
461 462  
462 463 // Continue only if valid
463   - if (emax > emin) {
  464 + if (log10_etrue_max > log10_etrue_min) {
464 465  
465 466 // Setup integration function
466 467 cta_nroi_kern integrand(this, &obs, &model, srcTime, obsEng, obsTime);
... ... @@ -470,14 +471,32 @@ double GCTAResponseIrf::nroi(const GModelSky&amp; model,
470 471 integral.fixed_iter(iter);
471 472  
472 473 // Do Romberg integration
473   - emin = std::log(emin);
474   - emax = std::log(emax);
475   - nroi += integral.romberg(emin, emax);
  474 + nroi += integral.romberg(log10_etrue_min, log10_etrue_max);
476 475  
477 476 } // endif: interval was valid
478 477  
479 478 } // endfor: looped over energy intervals
480 479  
  480 + // Debug energy dispersion computation
  481 + #if defined(G_DEBUG_NROI_EDISP)
  482 + apply_edisp(false);
  483 + const GEnergy& srcEng = obsEng;
  484 + double nroi_spatial = this->nroi(model, srcEng, srcTime, obsEng, obsTime, obs);
  485 + double nroi_spectral = model.spectral()->eval(srcEng, srcTime);
  486 + double nroi_temporal = model.temporal()->eval(srcTime);
  487 + double nroi_noedisp = nroi_spatial * nroi_spectral * nroi_temporal;
  488 + apply_edisp(true);
  489 + std::cout << "GCTAResponseIrf::nroi:";
  490 + std::cout << " obs_id=" << obs.id();
  491 + std::cout << " ereco=" << obsEng;
  492 + std::cout << " nroi(noedisp)=" << nroi_noedisp;
  493 + std::cout << " nroi(edisp)=" << nroi;
  494 + if (nroi_noedisp > 0.0) {
  495 + std::cout << " ratio(edisp/noedisp)=" << nroi/nroi_noedisp;
  496 + }
  497 + std::cout << std::endl;
  498 + #endif
  499 +
481 500 } // endif: energy dispersion requested
482 501  
483 502 // ... otherwise evaluate
... ... @@ -487,12 +506,12 @@ double GCTAResponseIrf::nroi(const GModelSky&amp; model,
487 506 const GEnergy& srcEng = obsEng;
488 507  
489 508 // Compute response components
490   - double npred_spatial = this->nroi(model, srcEng, srcTime, obsEng, obsTime, obs);
491   - double npred_spectral = model.spectral()->eval(srcEng, srcTime);
492   - double npred_temporal = model.temporal()->eval(srcTime);
  509 + double nroi_spatial = this->nroi(model, srcEng, srcTime, obsEng, obsTime, obs);
  510 + double nroi_spectral = model.spectral()->eval(srcEng, srcTime);
  511 + double nroi_temporal = model.temporal()->eval(srcTime);
493 512  
494 513 // Compute response
495   - nroi = npred_spatial * npred_spectral * npred_temporal;
  514 + nroi = nroi_spatial * nroi_spectral * nroi_temporal;
496 515  
497 516 } // endelse: no energy dispersion requested
498 517  
... ... @@ -1783,7 +1802,7 @@ double GCTAResponseIrf::psf_delta_max(const double&amp; theta,
1783 1802  
1784 1803  
1785 1804 /***********************************************************************//**
1786   - * @brief Return energy dispersion (in units of MeV^-1)
  1805 + * @brief Return energy dispersion (in units of MeV\f$^{-1}\f$)
1787 1806 *
1788 1807 * @param[in] obsEng Measured event energy.
1789 1808 * @param[in] theta Radial offset angle in camera (radians).
... ... @@ -1791,6 +1810,7 @@ double GCTAResponseIrf::psf_delta_max(const double&amp; theta,
1791 1810 * @param[in] zenith Zenith angle of telescope pointing (radians).
1792 1811 * @param[in] azimuth Azimuth angle of telescope pointing (radians).
1793 1812 * @param[in] srcLogEng Log10 of true photon energy (E/TeV).
  1813 + * @return Energy dispersion (MeV\f$^{-1}\f$).
1794 1814 ***************************************************************************/
1795 1815 double GCTAResponseIrf::edisp(const GEnergy& obsEng,
1796 1816 const double& theta,
... ... @@ -1799,13 +1819,11 @@ double GCTAResponseIrf::edisp(const GEnergy&amp; obsEng,
1799 1819 const double& azimuth,
1800 1820 const double& srcLogEng) const
1801 1821 {
1802   - // Compute log10 energy in TeV and linear energy in MeV
  1822 + // Compute log10 energy in TeV
1803 1823 double obsLogEng = obsEng.log10TeV();
1804   - double energy = obsEng.MeV();
1805 1824  
1806 1825 // Compute energy dispersion
1807   - double edisp = (*m_edisp)(obsLogEng, srcLogEng, theta, phi, zenith, azimuth) /
1808   - (gammalib::ln10 * energy);
  1826 + double edisp = (*m_edisp)(obsLogEng, srcLogEng, theta, phi, zenith, azimuth);
1809 1827  
1810 1828 // Return energy dispersion
1811 1829 return edisp;
... ... @@ -2707,11 +2725,12 @@ double GCTAResponseIrf::irf_diffuse(const GEvent&amp; event,
2707 2725 // Retrieve CTA pointing
2708 2726 const GCTAPointing& pnt = gammalib::cta_pnt(G_IRF_DIFFUSE, obs);
2709 2727  
2710   - // Try getting the IRF value from cache
  2728 + // Try getting the IRF value from cache. Do not do this for energy
  2729 + // dispersion since the the spectral model will change the IRF values
2711 2730 #if defined(G_USE_IRF_CACHE)
2712 2731 const GCTAEventList* list = dynamic_cast<const GCTAEventList*>(obs.events());
2713 2732 const GCTAEventAtom* atom = dynamic_cast<const GCTAEventAtom*>(&event);
2714   - if (list != NULL && atom != NULL) {
  2733 + if (list != NULL && atom != NULL && !use_edisp()) {
2715 2734 irf = list->irf_cache(source.name(), atom->index());
2716 2735 if (irf >= 0.0) {
2717 2736 has_irf = true;
... ... @@ -2826,9 +2845,11 @@ double GCTAResponseIrf::irf_diffuse(const GEvent&amp; event,
2826 2845 // Apply deadtime correction
2827 2846 irf *= obs.deadc(srcTime);
2828 2847  
2829   - // Put IRF value in cache
  2848 + // Put IRF value in cache. Only do this in no energy dispersion is
  2849 + // requested, since for energy dispersion the IRF values will
  2850 + // depend on the spectrum
2830 2851 #if defined(G_USE_IRF_CACHE)
2831   - if (list != NULL && atom != NULL) {
  2852 + if (list != NULL && atom != NULL && !use_edisp()) {
2832 2853 list->irf_cache(source.name(), atom->index(), irf);
2833 2854 }
2834 2855 #endif
... ... @@ -3315,9 +3336,11 @@ double GCTAResponseIrf::nroi_diffuse(const GModelSky&amp; model,
3315 3336 // Build unique identifier
3316 3337 std::string id = model.name() + "::" + obs.id();
3317 3338  
3318   - // Check if Nroi value is already in cache
  3339 + // Check if Nroi value is already in cache. Do not use the Npred cache
  3340 + // for energy dispersion since the actual value will depend on the
  3341 + // spectral model
3319 3342 #if defined(G_USE_NPRED_CACHE)
3320   - if (!m_npred_names.empty()) {
  3343 + if (!m_npred_names.empty() && !use_edisp()) {
3321 3344  
3322 3345 // Search for unique identifier, and if found, recover Npred value
3323 3346 // and break
... ... @@ -3407,18 +3430,23 @@ double GCTAResponseIrf::nroi_diffuse(const GModelSky&amp; model,
3407 3430 #if defined(G_DEBUG_NROI_DIFFUSE)
3408 3431 std::cout << "GCTAResponseIrf::nroi_diffuse:";
3409 3432 std::cout << " roi_psf_radius=" << roi_psf_radius;
3410   - std::cout << " nroi=" << npred;
  3433 + std::cout << " nroi=" << nroi;
  3434 + std::cout << " Etrue=" << srcEng;
  3435 + std::cout << " Ereco=" << obsEng;
3411 3436 std::cout << " id=" << id << std::endl;
3412 3437 #endif
3413 3438  
3414 3439 } // endif: offset angle range was valid
3415 3440  
3416   - // Store result in Npred cache
  3441 + // Store result in Npred cache in case that no energy dispersion
  3442 + // is requested
3417 3443 #if defined(G_USE_NPRED_CACHE)
3418   - m_npred_names.push_back(id);
3419   - m_npred_energies.push_back(srcEng);
3420   - m_npred_times.push_back(srcTime);
3421   - m_npred_values.push_back(nroi);
  3444 + if (!use_edisp()) {
  3445 + m_npred_names.push_back(id);
  3446 + m_npred_energies.push_back(srcEng);
  3447 + m_npred_times.push_back(srcTime);
  3448 + m_npred_values.push_back(nroi);
  3449 + }
3422 3450 #endif
3423 3451  
3424 3452 // Debug: Check for NaN
... ...
inst/cta/src/GCTAResponse_helpers.cpp
  View file @e5c3374
... ... @@ -133,7 +133,7 @@ cta_omega_intervals gammalib::limit_omega(const double&amp; min,
133 133 }
134 134 intervals.push_back(std::make_pair(interval_min,interval_max));
135 135 }
136   -
  136 +
137 137 // If the [min,max] interval overlaps with the [-domega,domega]
138 138 // interval shifted by -2pi then constrain the interval using the
139 139 // shifted interval
... ... @@ -339,7 +339,7 @@ double cta_npsf_kern_rad_azsym::eval(const double&amp; delta)
339 339 {
340 340 // Initialise PSF value
341 341 double value = 0.0;
342   -
  342 +
343 343 // Get arclength for given radius in radians
344 344 double phi = gammalib::cta_roi_arclength(delta,
345 345 m_psf,
... ... @@ -350,7 +350,7 @@ double cta_npsf_kern_rad_azsym::eval(const double&amp; delta)
350 350  
351 351 // If arclength is positive then compute the PSF value
352 352 if (phi > 0) {
353   -
  353 +
354 354 // Compute PSF value
355 355 value = m_rsp->psf(delta, m_theta, m_phi, m_zenith, m_azimuth, m_logE) *
356 356 phi * std::sin(delta);
... ... @@ -366,7 +366,7 @@ double cta_npsf_kern_rad_azsym::eval(const double&amp; delta)
366 366 std::cout << std::endl;
367 367 }
368 368 #endif
369   -
  369 +
370 370 } // endif: arclength was positive
371 371  
372 372 // Return
... ... @@ -375,25 +375,25 @@ double cta_npsf_kern_rad_azsym::eval(const double&amp; delta)
375 375  
376 376  
377 377 /***********************************************************************//**
378   - * @brief Integration kernel for nedisp() method
  378 + * @brief Integration kernel for GCTAResponseIrf::nroi method
379 379 *
380   - * @param[in] logEsrc Logarithm of true photon energy.
  380 + * @param[in] logEsrc Base 10 logarithm of true photon energy in MeV.
381 381 * @return Nroi.
382 382 ***************************************************************************/
383 383 double cta_nroi_kern::eval(const double& logEsrc)
384 384 {
385 385 // Set true energy
386 386 GEnergy srcEng;
387   - double expx = std::exp(logEsrc);
388   - srcEng.MeV(expx);
  387 + srcEng.log10MeV(logEsrc);
389 388  
390 389 // Compute response components
391 390 double nroi_spatial = m_rsp->nroi(*m_model, srcEng, m_srcTime, m_obsEng, m_obsTime, *m_obs);
392 391 double nroi_spectral = m_model->spectral()->eval(srcEng, m_srcTime);
393 392 double nroi_temporal = m_model->temporal()->eval(m_srcTime);
394   -
  393 +
395 394 // Compute response
396   - double nroi = nroi_spatial * nroi_spectral * nroi_temporal * expx;
  395 + double nroi = nroi_spatial * nroi_spectral * nroi_temporal *
  396 + gammalib::ln10 * m_obsEng.MeV();
397 397  
398 398 // Return response
399 399 return nroi;
... ... @@ -567,10 +567,10 @@ double cta_irf_radial_kern_omega::eval(const double&amp; omega)
567 567 {
568 568 // Compute PSF offset angle [radians]
569 569 double delta = std::acos(m_cos_psf + m_sin_psf * std::cos(omega));
570   -
  570 +
571 571 // Compute true photon offset angle in camera system [radians]
572 572 double offset = std::acos(m_cos_ph + m_sin_ph * std::cos(m_omega0 - omega));
573   -
  573 +
574 574 //TODO: Compute true photon azimuth angle in camera system [radians]
575 575 double azimuth = 0.0;
576 576  
... ... @@ -582,7 +582,7 @@ double cta_irf_radial_kern_omega::eval(const double&amp; omega)
582 582 if (m_rsp->use_edisp() && irf > 0.0) {
583 583 irf *= m_rsp->edisp(m_obsEng, offset, azimuth, m_zenith, m_azimuth, m_srcLogEng);
584 584 }
585   -
  585 +
586 586 // Compile option: Check for NaN/Inf
587 587 #if defined(G_NAN_CHECK)
588 588 if (gammalib::is_notanumber(irf) || gammalib::is_infinite(irf)) {
... ... @@ -944,7 +944,7 @@ double cta_irf_elliptical_kern_rho::eval(const double&amp; rho)
944 944 #endif
945 945  
946 946 } // endif: arc length was positive
947   -
  947 +
948 948 } // endif: rho was positive
949 949  
950 950 // Return result
... ... @@ -1035,7 +1035,7 @@ double cta_irf_elliptical_kern_omega::eval(const double&amp; omega)
1035 1035  
1036 1036 // Compute Psf offset angle [radians]
1037 1037 double delta = std::acos(m_cos_psf + m_sin_psf * std::cos(omega));
1038   -
  1038 +
1039 1039 // Compute true photon offset and azimuth angle in camera system
1040 1040 // [radians]
1041 1041 double theta = std::acos(m_cos_ph + m_sin_ph * std::cos(m_omega_pnt - omega));
... ... @@ -1242,7 +1242,7 @@ double cta_nroi_elliptical_kern_rho::eval(const double&amp; rho)
1242 1242 #endif
1243 1243  
1244 1244 } // endif: arc length was positive
1245   -
  1245 +
1246 1246 } // endif: rho was positive
1247 1247  
1248 1248 // Return Nroi
... ... @@ -1302,10 +1302,10 @@ double cta_nroi_elliptical_kern_omega::eval(const double&amp; omega)
1302 1302 std::cout << std::endl;
1303 1303 }
1304 1304 #endif
1305   -
  1305 +
1306 1306 // Continue only if model is positive
1307 1307 if (model > 0.0) {
1308   -
  1308 +
1309 1309 // Compute sky direction vector in native coordinates
1310 1310 double cos_omega = std::cos(omega_model);
1311 1311 double sin_omega = std::sin(omega_model);
... ... @@ -1339,7 +1339,7 @@ double cta_nroi_elliptical_kern_omega::eval(const double&amp; omega)
1339 1339 #endif
1340 1340  
1341 1341 } // endif: sky intensity was positive
1342   -
  1342 +
1343 1343 // Return Nroi
1344 1344 return nroi;
1345 1345 }
... ... @@ -2210,7 +2210,7 @@ double cta_psf_elliptical_kern_rho::eval(const double&amp; rho)
2210 2210 #endif
2211 2211  
2212 2212 } // endif: arc length was positive
2213   -
  2213 +
2214 2214 } // endif: rho was positive
2215 2215  
2216 2216 // Return result
... ...