Pierrick Martin / gammalib

Commit 2ec38a7eec0664db97c5d22fe8028b345ff42e5d

Authored by Jürgen Knödlseder
1 parent 05fc86b0

Revise GCTAEdispRmf class (#2704)

Inline Side-by-side
Showing 6 changed files with 235 additions and 235 deletions
ChangeLog
  View file @2ec38a7
1   -2018-10-25
  1 +2018-10-31
2 2  
3 3 * Version 1.6.0 released
4 4 ========================
... ...
1 1 New Features and Important Changes in GammaLib 1.6.0
2 2  
3   -25 October 2018
  3 +31 October 2018
4 4  
5 5  
6 6 1. Introduction
... ...
inst/cta/include/GCTAEdispRmf.hpp
  View file @2ec38a7
... ... @@ -50,14 +50,12 @@ class GEbounds;
50 50 * The energy dispersion is defined as
51 51 *
52 52 * \f[
53   - * E_{\rm disp}(E_{\rm true}, E_{\rm reco}) =
54   - * \frac{E_{\rm disp}(E_{\rm true}, \log_{10} E_{\rm reco}, \theta)}
55   - * {\log_{10} E_{\rm reco}}
  53 + * E_{\rm disp}(E_{\rm reco} | E_{\rm true})
56 54 * \f]
57 55 *
58 56 * in units of MeV\f$^{-1}\f$ where
59   - * \f$E_{\rm reco}\f$ is the reconstructed energy in units of MeV, and
60   - * \f$E_{\rm true}\f$ is the true energy in units of MeV.
  57 + * \f$E_{\rm reco}\f$ is the reconstructed energy, and
  58 + * \f$E_{\rm true}\f$ is the true energy.
61 59 ***************************************************************************/
62 60 class GCTAEdispRmf : public GCTAEdisp {
63 61  
... ... @@ -114,49 +112,41 @@ private:
114 112 void copy_members(const GCTAEdispRmf& psf);
115 113 void free_members(void);
116 114 void set_matrix(void);
  115 + void set_max_edisp(void);
117 116 void set_cache(void) const;
118   - void set_max_edisp(void) const;
119   - void update(const double& logEsrc, const double& logEobs) const;
120   - void compute_ebounds_obs(const double& theta = 0.0,
121   - const double& phi = 0.0,
122   - const double& zenith = 0.0,
123   - const double& azimuth = 0.0) const;
124   - void compute_ebounds_src(const double& theta = 0.0,
125   - const double& phi = 0.0,
126   - const double& zenith = 0.0,
127   - const double& azimuth = 0.0) const;
  117 + void update(const GEnergy& ereco, const GEnergy& etrue) const;
  118 + void compute_ereco_bounds(void) const;
  119 + void compute_etrue_bounds(void) const;
128 120  
129 121 // Members
130 122 GFilename m_filename; //!< Name of response file
131 123 GRmf m_rmf; //!< Redistribution matrix file
132 124 GMatrixSparse m_matrix; //!< Normalised redistribution matrix
  125 + double m_max_edisp; //!< Maximum energy dispersion value for MC
133 126  
134 127 // Interpolation cache
135   - mutable GNodeArray m_etrue; //!< Array of log10(Etrue)
136   - mutable GNodeArray m_emeasured; //!< Array of log10(Emeasured)
137   - mutable double m_last_etrue; //!< Last log10(Etrue)
138   - mutable double m_last_emeasured; //!< Last log10(Emeasured)
139   - mutable int m_itrue1; //!< Index of left Etrue
140   - mutable int m_itrue2; //!< Index of right Etrue
141   - mutable int m_imeas1; //!< Index of left Emeasured
142   - mutable int m_imeas2; //!< Index of right Emeasured
143   - mutable double m_wgt1; //!< Weight of lower left node
144   - mutable double m_wgt2; //!< Weight of upper left node
145   - mutable double m_wgt3; //!< Weight of lower right node
146   - mutable double m_wgt4; //!< Weight of upper right node
147   -
148   - // Monte Carlo cache
149   - mutable double m_max_edisp;
150   - mutable double m_last_theta_obs;
151   - mutable double m_last_theta_src;
152   - mutable GEnergy m_last_etrue_ereco_bounds;
153   - mutable GEnergy m_last_ereco_etrue_bounds;
154   - mutable int m_index_obs;
155   - mutable int m_index_src;
156   - mutable bool m_ebounds_obs_computed;
157   - mutable bool m_ebounds_src_computed;
158   - mutable std::vector<GEbounds> m_ebounds_obs;
159   - mutable std::vector<GEbounds> m_ebounds_src;
  128 + mutable GNodeArray m_etrue; //!< Array of log10(Etrue)
  129 + mutable GNodeArray m_ereco; //!< Array of log10(Ereco)
  130 + mutable GEnergy m_last_etrue; //!< Last true energy
  131 + mutable GEnergy m_last_ereco; //!< Last reconstructed energy
  132 + mutable int m_itrue1; //!< Index of left Etrue
  133 + mutable int m_itrue2; //!< Index of right Etrue
  134 + mutable int m_ireco1; //!< Index of left Ereco
  135 + mutable int m_ireco2; //!< Index of right Ereco
  136 + mutable double m_wgt1; //!< Weight of lower left node
  137 + mutable double m_wgt2; //!< Weight of upper left node
  138 + mutable double m_wgt3; //!< Weight of lower right node
  139 + mutable double m_wgt4; //!< Weight of upper right node
  140 +
  141 + // Computation cache
  142 + mutable bool m_ereco_bounds_computed;
  143 + mutable bool m_etrue_bounds_computed;
  144 + mutable GEnergy m_last_etrue_bounds;
  145 + mutable GEnergy m_last_ereco_bounds;
  146 + mutable int m_index_ereco;
  147 + mutable int m_index_etrue;
  148 + mutable std::vector<GEbounds> m_ereco_bounds;
  149 + mutable std::vector<GEbounds> m_etrue_bounds;
160 150 };
161 151  
162 152  
... ...
inst/cta/src/GCTAEdisp2D.cpp
  View file @2ec38a7
... ... @@ -46,7 +46,7 @@
46 46  
47 47 /* __ Method name definitions ____________________________________________ */
48 48 #define G_READ "GCTAEdisp2D::read(GFitsTable&)"
49   -#define G_MC "GCTAEdisp2D::mc(GRan&, double&, double&, double&, double&, "\
  49 +#define G_MC "GCTAEdisp2D::mc(GRan&, GEnergy&, double&, double&, double&, "\
50 50 "double&)"
51 51 #define G_FETCH "GCTAEdisp2D::fetch()"
52 52  
... ... @@ -220,7 +220,7 @@ double GCTAEdisp2D::operator()(const GEnergy&amp; ereco,
220 220 if ((logEsrc >= m_logEsrc_min) && (logEsrc <= m_logEsrc_max) &&
221 221 (theta >= m_theta_min) && (theta <= m_theta_max)) {
222 222  
223   - // Compute migration (migra=Ereco/Etrue
  223 + // Compute migration (migra=Ereco/Etrue)
224 224 double migra = ereco.MeV() / etrue.MeV();
225 225  
226 226 // Continue only if migration is in validity range
... ... @@ -470,13 +470,13 @@ void GCTAEdisp2D::save(const GFilename&amp; filename, const bool&amp; clobber) const
470 470 * @param[in] phi Azimuth angle in camera system (radians). Not used.
471 471 * @param[in] zenith Zenith angle in Earth system (radians). Not used.
472 472 * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
473   - * @return Observed energy.
  473 + * @return Reconstructed energy.
474 474 *
475 475 * @exception GException::invalid_return_value
476 476 * No energy dispersion information found for parameters or
477   - * energy dispersion matrix is invalid.
  477 + * energy dispersion matrix is empty.
478 478 *
479   - * Draws observed energy value given a true energy @p etrue and offset
  479 + * Draws reconstructed energy value given a true energy @p etrue and offset
480 480 * angle @p theta. If no energy dispersion information is available the
481 481 * method will return the true photon energy.
482 482 ***************************************************************************/
... ...
inst/cta/src/GCTAEdispRmf.cpp
  View file @2ec38a7
... ... @@ -43,10 +43,13 @@
43 43  
44 44 /* __ Method name definitions ____________________________________________ */
45 45 #define G_LOAD "GCTAEdispRmf::load(GFilename&)"
  46 +#define G_MC "GCTAEdispRmf::mc(GRan&, GEnergy&, double&, double&, double&, "\
  47 + "double&)"
46 48  
47 49 /* __ Macros _____________________________________________________________ */
48 50  
49 51 /* __ Coding definitions _________________________________________________ */
  52 +//#define G_LINEAR_ERECO_INTERPOLATION //!< Linear interpolation for Ereco
50 53  
51 54 /* __ Debug definitions __________________________________________________ */
52 55  
... ... @@ -173,14 +176,12 @@ GCTAEdispRmf&amp; GCTAEdispRmf::operator=(const GCTAEdispRmf&amp; edisp)
173 176 * Returns the energy dispersion
174 177 *
175 178 * \f[
176   - * E_{\rm disp}(E_{\rm true}, E_{\rm reco}) =
177   - * \frac{E_{\rm disp}(E_{\rm true}, \log_{10} E_{\rm reco}, \theta)}
178   - * {\log_{10} E_{\rm reco}}
  179 + * E_{\rm disp}(E_{\rm reco} | E_{\rm true})
179 180 * \f]
180 181 *
181 182 * in units of MeV\f$^{-1}\f$ where
182   - * \f$E_{\rm reco}\f$ is the reconstructed energy in units of MeV, and
183   - * \f$E_{\rm true}\f$ is the true energy in units of MeV.
  183 + * \f$E_{\rm reco}\f$ is the reconstructed energy, and
  184 + * \f$E_{\rm true}\f$ is the true energy.
184 185 ***************************************************************************/
185 186 double GCTAEdispRmf::operator()(const GEnergy& ereco,
186 187 const GEnergy& etrue,
... ... @@ -189,29 +190,20 @@ double GCTAEdispRmf::operator()(const GEnergy&amp; ereco,
189 190 const double& zenith,
190 191 const double& azimuth) const
191 192 {
192   - // Get log10 of true and reconstructed photon energies
193   - double logEsrc = etrue.log10TeV();
194   - double logEobs = ereco.log10TeV();
195   -
196 193 // Update indexes and weighting factors for interpolation
197   - update(logEobs, logEsrc);
  194 + update(ereco, etrue);
198 195  
199 196 // Perform interpolation
200   - double edisp = m_wgt1 * m_matrix(m_itrue1, m_imeas1) +
201   - m_wgt2 * m_matrix(m_itrue1, m_imeas2) +
202   - m_wgt3 * m_matrix(m_itrue2, m_imeas1) +
203   - m_wgt4 * m_matrix(m_itrue2, m_imeas2);
  197 + double edisp = m_wgt1 * m_matrix(m_itrue1, m_ireco1) +
  198 + m_wgt2 * m_matrix(m_itrue1, m_ireco2) +
  199 + m_wgt3 * m_matrix(m_itrue2, m_ireco1) +
  200 + m_wgt4 * m_matrix(m_itrue2, m_ireco2);
204 201  
205 202 // Make sure that energy dispersion is not negative
206 203 if (edisp < 0.0) {
207 204 edisp = 0.0;
208 205 }
209 206  
210   - // If the energy dispersion is positive, return it per MeV
211   - else {
212   - edisp /= (gammalib::ln10 * ereco.MeV());
213   - }
214   -
215 207 // Return energy dispersion
216 208 return edisp;
217 209 }
... ... @@ -292,9 +284,14 @@ void GCTAEdispRmf::load(const GFilename&amp; filename)
292 284 * @param[in] phi Azimuth angle in camera system (radians). Not used.
293 285 * @param[in] zenith Zenith angle in Earth system (radians). Not used.
294 286 * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
295   - * @return Energy.
  287 + * @return Reconstructed energy.
296 288 *
297   - * Draws observed energy value from RMF matrix.
  289 + * @exception GException::invalid_return_value
  290 + * Energy dispersion matrix is empty.
  291 + *
  292 + * Draws reconstructed energy value from RMF matrix given a true energy
  293 + * @p etrue. If no energy dispersion information is available the method
  294 + * will return the true photon energy.
298 295 ***************************************************************************/
299 296 GEnergy GCTAEdispRmf::mc(GRan& ran,
300 297 const GEnergy& etrue,
... ... @@ -303,25 +300,33 @@ GEnergy GCTAEdispRmf::mc(GRan&amp; ran,
303 300 const double& zenith,
304 301 const double& azimuth) const
305 302 {
306   - // Initialise energies
307   - GEnergy ereco;
  303 + // Initialise reconstructed event energy with true photon energy
  304 + GEnergy ereco = etrue;
308 305  
309   - // Get boundaries for rejection method
  306 + // Get boundaries for reconstructed energy
310 307 GEbounds ebounds = ereco_bounds(etrue, theta, phi, zenith, azimuth);
311   - double emin = ebounds.emin().log10TeV();
312   - double emax = ebounds.emax().log10TeV();
  308 + double emin = ebounds.emin().TeV();
  309 + double emax = ebounds.emax().TeV();
  310 +
  311 + // Get maximum energy dispersion value (including some margin)
  312 + double max_edisp = 1.5 * m_max_edisp;
  313 +
  314 + // Throw an exception if maximum energy dispersion is zero
  315 + if (max_edisp <= 0.0) {
  316 + std::string msg = "Energy dispersion matrix is empty. Please provide "
  317 + "a valid energy dispersion matrix.";
  318 + throw GException::invalid_return_value(G_MC, msg);
  319 + }
313 320  
314 321 // Find energy by rejection method
315 322 double ewidth = emax - emin;
316   - double logEobs = 0.5*(emin+emax);
317   - if (m_max_edisp > 0.0) {
  323 + if (max_edisp > 0.0) {
318 324 double f = 0.0;
319 325 double ftest = 1.0;
320 326 while (ftest > f) {
321   - logEobs = emin + ewidth * ran.uniform();
322   - ereco.log10TeV(logEobs);
323   - f = operator()(ereco, etrue, theta, phi, zenith, azimuth);
324   - ftest = ran.uniform() * m_max_edisp;
  327 + ereco.TeV(emin + ewidth * ran.uniform());
  328 + f = operator()(ereco, etrue, theta, phi, zenith, azimuth);
  329 + ftest = ran.uniform() * max_edisp;
325 330 }
326 331 }
327 332  
... ... @@ -350,30 +355,29 @@ GEbounds GCTAEdispRmf::ereco_bounds(const GEnergy&amp; etrue,
350 355 const double& azimuth) const
351 356 {
352 357 // Compute only if parameters changed
353   - if (!m_ebounds_obs_computed || theta != m_last_theta_obs) {
  358 + if (!m_ereco_bounds_computed) {
354 359  
355 360 // Store last theta value
356   - m_ebounds_obs_computed = true;
357   - m_last_theta_obs = theta;
358   - m_last_etrue_ereco_bounds.log10TeV(-30.0); // force update
  361 + m_ereco_bounds_computed = true;
  362 + m_last_etrue_bounds.TeV(0.0); // force update
359 363  
360   - // Compute ebounds_obs
361   - compute_ebounds_obs(theta, phi, zenith, azimuth);
  364 + // Compute m_ereco_bounds
  365 + compute_ereco_bounds();
362 366  
363 367 }
364 368  
365 369 // Search index only if logEsrc has changed
366   - if (etrue != m_last_etrue_ereco_bounds) {
  370 + if (etrue != m_last_etrue_bounds) {
367 371  
368 372 // Store last photon energy
369   - m_last_etrue_ereco_bounds = etrue;
  373 + m_last_etrue_bounds = etrue;
370 374  
371 375 // Get true photon energy in TeV
372 376 double etrue_TeV = etrue.TeV();
373 377  
374 378 // Find right index with bisection
375 379 int low = 0;
376   - int high = m_ebounds_obs.size() - 1;
  380 + int high = m_ereco_bounds.size() - 1;
377 381 while ((high-low) > 1) {
378 382 int mid = (low+high) / 2;
379 383 double e = m_rmf.etrue().emin(mid).TeV();
... ... @@ -386,12 +390,12 @@ GEbounds GCTAEdispRmf::ereco_bounds(const GEnergy&amp; etrue,
386 390 }
387 391  
388 392 // Index found
389   - m_index_obs = low;
  393 + m_index_ereco = low;
390 394  
391 395 }
392 396  
393 397 // Return energy boundaries
394   - return (m_ebounds_obs[m_index_obs]);
  398 + return (m_ereco_bounds[m_index_ereco]);
395 399 }
396 400  
397 401  
... ... @@ -415,30 +419,29 @@ GEbounds GCTAEdispRmf::etrue_bounds(const GEnergy&amp; ereco,
415 419 const double& azimuth) const
416 420 {
417 421 // Compute only if parameters changed
418   - if (!m_ebounds_src_computed || theta != m_last_theta_src) {
  422 + if (!m_etrue_bounds_computed) {
419 423  
420   - // Store last value
421   - m_ebounds_src_computed = true;
422   - m_last_theta_src = theta;
423   - m_last_ereco_etrue_bounds.log10TeV(-30.0); // force update
  424 + // Set computation flag
  425 + m_etrue_bounds_computed = true;
  426 + m_last_ereco_bounds.TeV(0.0); // force update
424 427  
425   - // Compute ebounds_obs
426   - compute_ebounds_src(theta, phi, zenith, azimuth);
  428 + // Compute m_etrue_bounds
  429 + compute_etrue_bounds();
427 430  
428 431 }
429 432  
430 433 // Search index only if ereco has changed
431   - if (ereco != m_last_ereco_etrue_bounds) {
  434 + if (ereco != m_last_ereco_bounds) {
432 435  
433 436 // Store last value
434   - m_last_ereco_etrue_bounds = ereco;
  437 + m_last_ereco_bounds = ereco;
435 438  
436 439 // Get reconstructed energy in TeV
437 440 double ereco_TeV = ereco.TeV();
438 441  
439 442 // Find right index with bisection
440 443 int low = 0;
441   - int high = m_ebounds_src.size() - 1;
  444 + int high = m_etrue_bounds.size() - 1;
442 445 while ((high-low) > 1) {
443 446 int mid = (low+high) / 2;
444 447 double e = m_rmf.emeasured().emin(mid).TeV();
... ... @@ -451,12 +454,12 @@ GEbounds GCTAEdispRmf::etrue_bounds(const GEnergy&amp; ereco,
451 454 }
452 455  
453 456 // Index found
454   - m_index_src = low;
  457 + m_index_etrue = low;
455 458  
456 459 }
457 460  
458 461 // Return energy boundaries
459   - return (m_ebounds_src[m_index_src]);
  462 + return (m_etrue_bounds[m_index_etrue]);
460 463 }
461 464  
462 465  
... ... @@ -474,7 +477,7 @@ GEbounds GCTAEdispRmf::etrue_bounds(const GEnergy&amp; ereco,
474 477 *
475 478 * \f[
476 479 * \int_{E_{\rm reco}^{\rm min}}^{E_{\rm reco}^{\rm max}}
477   - * E_{\rm disp}(E_{\rm true}, E_{\rm reco}) \, dE_{\rm reco}
  480 + * E_{\rm disp}(E_{\rm reco} | E_{\rm true}) \, dE_{\rm reco}
478 481 * \f]
479 482 *
480 483 * where
... ... @@ -491,8 +494,13 @@ double GCTAEdispRmf::prob_erecobin(const GEnergy&amp; ereco_min,
491 494  
492 495 // Get log10 of energies
493 496 double logEsrc = etrue.log10TeV();
  497 + #if defined(G_LINEAR_ERECO_INTERPOLATION)
  498 + double logEobs_min = ereco_min.TeV();
  499 + double logEobs_max = ereco_max.TeV();
  500 + #else
494 501 double logEobs_min = ereco_min.log10TeV();
495 502 double logEobs_max = ereco_max.log10TeV();
  503 + #endif
496 504  
497 505 // Set logEsrc for node array interpolation
498 506 m_etrue.set_value(logEsrc);
... ... @@ -510,23 +518,27 @@ double GCTAEdispRmf::prob_erecobin(const GEnergy&amp; ereco_min,
510 518 double f2 = 0.0;
511 519  
512 520 // Loop over all measured energy nodes
513   - for (int i = 0; i < m_emeasured.size(); ++i) {
  521 + for (int i = 0; i < m_ereco.size(); ++i) {
514 522  
515 523 // If measured energy is below logEobs_min then skip node
516   - if (m_emeasured[i] <= logEobs_min) {
  524 + if (m_ereco[i] <= logEobs_min) {
517 525 continue;
518 526 }
519 527  
520 528 // If measured energy is above maximum migration value then use
521 529 // logEobs_max as measured energy
522   - if (m_emeasured[i] > logEobs_max) {
  530 + if (m_ereco[i] > logEobs_max) {
523 531 x2 = ereco_max;
524 532 f2 = this->operator()(x2, etrue);
525 533 }
526 534  
527 535 // ... otherwise use measured energy
528 536 else {
529   - x2.log10TeV(m_emeasured[i]);
  537 + #if defined(G_LINEAR_ERECO_INTERPOLATION)
  538 + x2.TeV(m_ereco[i]);
  539 + #else
  540 + x2.log10TeV(m_ereco[i]);
  541 + #endif
530 542 f2 = wgt_left * m_matrix(inx_left, i) +
531 543 wgt_right * m_matrix(inx_right, i);
532 544 if (f2 < 0.0) {
... ... @@ -535,14 +547,14 @@ double GCTAEdispRmf::prob_erecobin(const GEnergy&amp; ereco_min,
535 547 }
536 548  
537 549 // Compute integral
538   - prob += 0.5 * (f1 + f2) * (x2.log10TeV() - x1.log10TeV());
  550 + prob += 0.5 * (f1 + f2) * (x2.MeV() - x1.MeV());
539 551  
540 552 // Set second node as first node
541 553 x1 = x2;
542 554 f1 = f2;
543 555  
544 556 // If node energy is above migra_max then break now
545   - if (m_emeasured[i] > logEobs_max) {
  557 + if (m_ereco[i] > logEobs_max) {
546 558 break;
547 559 }
548 560  
... ... @@ -553,7 +565,7 @@ double GCTAEdispRmf::prob_erecobin(const GEnergy&amp; ereco_min,
553 565 if (x1 < ereco_max) {
554 566 x2 = ereco_max;
555 567 f2 = this->operator()(x2, etrue);
556   - prob += 0.5 * (f1 + f2) * (x2.log10TeV() - x1.log10TeV());
  568 + prob += 0.5 * (f1 + f2) * (x2.MeV() - x1.MeV());
557 569 }
558 570  
559 571 // Return probability
... ... @@ -616,33 +628,31 @@ void GCTAEdispRmf::init_members(void)
616 628 m_filename.clear();
617 629 m_rmf.clear();
618 630 m_matrix.clear();
  631 + m_max_edisp = 0.0;
619 632  
620 633 // Initialise interpolation cache
621 634 m_etrue.clear();
622   - m_emeasured.clear();
623   - m_last_etrue = 9999.0;
624   - m_last_emeasured = 9999.0;
625   - m_itrue1 = 0;
626   - m_itrue2 = 0;
627   - m_imeas1 = 0;
628   - m_imeas2 = 0;
629   - m_wgt1 = 0.0;
630   - m_wgt2 = 0.0;
631   - m_wgt3 = 0.0;
632   - m_wgt4 = 0.0;
633   -
634   - // Initialise Monte Carlo cache
635   - m_max_edisp = 0.0;
636   - m_last_theta_obs = -1.0;
637   - m_last_theta_src = -1.0;
638   - m_last_etrue_ereco_bounds.log10TeV(-30.0);
639   - m_last_ereco_etrue_bounds.log10TeV(-30.0);
640   - m_index_obs = 0;
641   - m_index_src = 0;
642   - m_ebounds_obs_computed = false;
643   - m_ebounds_src_computed = false;
644   - m_ebounds_obs.clear();
645   - m_ebounds_src.clear();
  635 + m_ereco.clear();
  636 + m_last_etrue.clear();
  637 + m_last_ereco.clear();
  638 + m_itrue1 = 0;
  639 + m_itrue2 = 0;
  640 + m_ireco1 = 0;
  641 + m_ireco2 = 0;
  642 + m_wgt1 = 0.0;
  643 + m_wgt2 = 0.0;
  644 + m_wgt3 = 0.0;
  645 + m_wgt4 = 0.0;
  646 +
  647 + // Initialise computation cache
  648 + m_ereco_bounds_computed = false;
  649 + m_etrue_bounds_computed = false;
  650 + m_index_ereco = 0;
  651 + m_index_etrue = 0;
  652 + m_last_etrue_bounds.clear();
  653 + m_last_ereco_bounds.clear();
  654 + m_ereco_bounds.clear();
  655 + m_etrue_bounds.clear();
646 656  
647 657 // Return
648 658 return;
... ... @@ -657,36 +667,34 @@ void GCTAEdispRmf::init_members(void)
657 667 void GCTAEdispRmf::copy_members(const GCTAEdispRmf& edisp)
658 668 {
659 669 // Copy members
660   - m_filename = edisp.m_filename;
661   - m_rmf = edisp.m_rmf;
662   - m_matrix = edisp.m_matrix;
  670 + m_filename = edisp.m_filename;
  671 + m_rmf = edisp.m_rmf;
  672 + m_matrix = edisp.m_matrix;
  673 + m_max_edisp = edisp.m_max_edisp;
663 674  
664 675 // Copy interpolation cache
665   - m_etrue = edisp.m_etrue;
666   - m_emeasured = edisp.m_emeasured;
667   - m_last_etrue = edisp.m_last_etrue;
668   - m_last_emeasured = edisp.m_last_emeasured;
669   - m_itrue1 = edisp.m_itrue1;
670   - m_itrue2 = edisp.m_itrue2;
671   - m_imeas1 = edisp.m_imeas1;
672   - m_imeas2 = edisp.m_imeas2;
673   - m_wgt1 = edisp.m_wgt1;
674   - m_wgt2 = edisp.m_wgt2;
675   - m_wgt3 = edisp.m_wgt3;
676   - m_wgt4 = edisp.m_wgt4;
677   -
678   - // Copy Monte Carlo cache
679   - m_max_edisp = edisp.m_max_edisp;
680   - m_last_theta_obs = edisp.m_last_theta_obs;
681   - m_last_theta_src = edisp.m_last_theta_src;
682   - m_last_etrue_ereco_bounds = edisp.m_last_etrue_ereco_bounds;
683   - m_last_ereco_etrue_bounds = edisp.m_last_ereco_etrue_bounds;
684   - m_index_obs = edisp.m_index_obs;
685   - m_index_src = edisp.m_index_src;
686   - m_ebounds_obs_computed = edisp.m_ebounds_obs_computed;
687   - m_ebounds_src_computed = edisp.m_ebounds_src_computed;
688   - m_ebounds_obs = edisp.m_ebounds_obs;
689   - m_ebounds_src = edisp.m_ebounds_src;
  676 + m_etrue = edisp.m_etrue;
  677 + m_ereco = edisp.m_ereco;
  678 + m_last_etrue = edisp.m_last_etrue;
  679 + m_last_ereco = edisp.m_last_ereco;
  680 + m_itrue1 = edisp.m_itrue1;
  681 + m_itrue2 = edisp.m_itrue2;
  682 + m_ireco1 = edisp.m_ireco1;
  683 + m_ireco2 = edisp.m_ireco2;
  684 + m_wgt1 = edisp.m_wgt1;
  685 + m_wgt2 = edisp.m_wgt2;
  686 + m_wgt3 = edisp.m_wgt3;
  687 + m_wgt4 = edisp.m_wgt4;
  688 +
  689 + // Copy computation cache
  690 + m_ereco_bounds_computed = edisp.m_ereco_bounds_computed;
  691 + m_etrue_bounds_computed = edisp.m_etrue_bounds_computed;
  692 + m_index_ereco = edisp.m_index_ereco;
  693 + m_index_etrue = edisp.m_index_etrue;
  694 + m_last_etrue_bounds = edisp.m_last_etrue_bounds;
  695 + m_last_ereco_bounds = edisp.m_last_ereco_bounds;
  696 + m_ereco_bounds = edisp.m_ereco_bounds;
  697 + m_etrue_bounds = edisp.m_etrue_bounds;
690 698  
691 699 // Return
692 700 return;
... ... @@ -715,28 +723,42 @@ void GCTAEdispRmf::set_matrix(void)
715 723 // Initialize matrix
716 724 m_matrix = GMatrixSparse(rows, columns);
717 725  
  726 + // Get reconstructued energy bin boundaries
  727 + GEbounds ebounds = m_rmf.emeasured();
  728 +
718 729 // Fill matrix elements
719   - for (int itrue = 0; itrue < rows; ++itrue) {
720   - for (int imeasured = 0; imeasured < columns; ++imeasured) {
721   - m_matrix(itrue, imeasured) = m_rmf(itrue, imeasured);
  730 + for (int ireco = 0; ireco < columns; ++ireco) {
  731 +
  732 + // Compute reconstructued energy bin width in MeV
  733 + double ewidth = ebounds.emax(ireco).MeV() - ebounds.emin(ireco).MeV();
  734 +
  735 + // If bin width is positive then fill all matrix bins by deviding
  736 + // the RMF values by the bin width. This assures that the stored
  737 + // values are per MeV of reconstructed energy.
  738 + if (ewidth > 0.0) {
  739 + for (int itrue = 0; itrue < rows; ++itrue) {
  740 + m_matrix(itrue, ireco) = m_rmf(itrue, ireco) / ewidth;
  741 + }
722 742 }
723   - }
  743 +
  744 + } // endfor: looped over all reconstructed energies
724 745  
725 746 // Initialise row sums
726   - std::vector<double> row_sums;
727   - row_sums.reserve(rows);
  747 + std::vector<double> row_sums(rows, 0.0);
728 748  
729 749 // Normalize matrix elements
730 750 for (int itrue = 0; itrue < rows; ++itrue) {
731 751  
  752 + // Get true energy
  753 + GEnergy etrue = m_rmf.etrue().elogmean(itrue);
  754 +
732 755 // Get integration boundaries
733   - GEbounds ebounds = ereco_bounds(m_rmf.etrue().elogmean(itrue), 0.0);
  756 + GEbounds ebounds = ereco_bounds(etrue);
734 757  
735 758 // Integrate contributions over energy boundaries
736 759 double sum = 0.0;
737 760 for (int i = 0; i < ebounds.size(); ++i) {
738   - sum += prob_erecobin(ebounds.emin(i), ebounds.emax(i),
739   - m_rmf.etrue().elogmean(itrue), 0.0);
  761 + sum += prob_erecobin(ebounds.emin(i), ebounds.emax(i), etrue, 0.0);
740 762 }
741 763  
742 764 // Store matrix sum
... ... @@ -748,8 +770,8 @@ void GCTAEdispRmf::set_matrix(void)
748 770 for (int itrue = 0; itrue < rows; ++itrue) {
749 771 double sum = row_sums[itrue];
750 772 if (sum > 0.0) {
751   - for (int imeasured = 0; imeasured < columns; ++imeasured) {
752   - m_matrix(itrue, imeasured) /= sum;
  773 + for (int ireco = 0; ireco < columns; ++ireco) {
  774 + m_matrix(itrue, ireco) /= sum;
753 775 }
754 776 }
755 777 }
... ... @@ -768,7 +790,7 @@ void GCTAEdispRmf::set_cache(void) const
768 790 {
769 791 // Clear node arrays
770 792 m_etrue.clear();
771   - m_emeasured.clear();
  793 + m_ereco.clear();
772 794  
773 795 // Set log10(Etrue) nodes
774 796 for (int i = 0; i < m_rmf.ntrue(); ++i) {
... ... @@ -777,7 +799,11 @@ void GCTAEdispRmf::set_cache(void) const
777 799  
778 800 // Set log10(Emeasured) nodes
779 801 for (int i = 0; i < m_rmf.nmeasured(); ++i) {
780   - m_emeasured.append(m_rmf.emeasured().elogmean(i).log10TeV());
  802 + #if defined(G_LINEAR_ERECO_INTERPOLATION)
  803 + m_ereco.append(m_rmf.emeasured().elogmean(i).TeV());
  804 + #else
  805 + m_ereco.append(m_rmf.emeasured().elogmean(i).log10TeV());
  806 + #endif
781 807 }
782 808  
783 809 // Return
... ... @@ -788,15 +814,15 @@ void GCTAEdispRmf::set_cache(void) const
788 814 /***********************************************************************//**
789 815 * @brief Set maximum energy dispersion value
790 816 ***************************************************************************/
791   -void GCTAEdispRmf::set_max_edisp(void) const
  817 +void GCTAEdispRmf::set_max_edisp(void)
792 818 {
793 819 // Initialise maximum
794 820 m_max_edisp = 0.0;
795 821  
796 822 // Loop over all response table elements
797   - for (int i = 0; i < m_matrix.rows(); ++i) {
798   - for (int k = 0; k < m_matrix.columns(); ++k) {
799   - double value = m_matrix(i,k);
  823 + for (int itrue = 0; itrue < m_matrix.rows(); ++itrue) {
  824 + for (int ireco = 0; ireco < m_matrix.columns(); ++ireco) {
  825 + double value = m_matrix(itrue, ireco);
800 826 if (value > m_max_edisp) {
801 827 m_max_edisp = value;
802 828 }
... ... @@ -811,37 +837,45 @@ void GCTAEdispRmf::set_max_edisp(void) const
811 837 /***********************************************************************//**
812 838 * @brief Update cache
813 839 *
814   - * @param[in] logEsrc Log10 of true energy (\f$\log_{10}\f$ TeV).
815   - * @param[in] logEobs Log10 of reconstructed energy (\f$\log_{10}\f$ TeV).
  840 + * @param[in] ereco Reconstructed energy.
  841 + * @param[in] etrue True energy.
816 842 *
817 843 * Updates the interpolation cache. The interpolation cache is composed
818 844 * of four indices and weights that define 4 data values of the RMF matrix
819 845 * that are used for bilinear interpolation.
820 846 ***************************************************************************/
821   -void GCTAEdispRmf::update(const double& logEsrc, const double& logEobs) const
  847 +void GCTAEdispRmf::update(const GEnergy& ereco, const GEnergy& etrue) const
822 848 {
823 849 // Update cache only of arguments have changed
824   - if (logEsrc != m_last_etrue || logEobs != m_last_emeasured) {
  850 + if (ereco != m_last_etrue || etrue != m_last_ereco) {
825 851  
826 852 // Store actual values
827   - m_last_etrue = logEsrc;
828   - m_last_emeasured = logEobs;
  853 + m_last_ereco = ereco;
  854 + m_last_etrue = etrue;
  855 +
  856 + // Get log10 of true and reconstructed photon energies
  857 + double logEsrc = etrue.log10TeV();
  858 + #if defined(G_LINEAR_ERECO_INTERPOLATION)
  859 + double logEobs = ereco.TeV();
  860 + #else
  861 + double logEobs = ereco.log10TeV();
  862 + #endif
829 863  
830 864 // Set values for node arrays
831 865 m_etrue.set_value(logEsrc);
832   - m_emeasured.set_value(logEobs);
  866 + m_ereco.set_value(logEobs);
833 867  
834 868 // Set indices for bi-linear interpolation
835 869 m_itrue1 = m_etrue.inx_left();
836 870 m_itrue2 = m_etrue.inx_right();
837   - m_imeas1 = m_emeasured.inx_left();
838   - m_imeas2 = m_emeasured.inx_right();
  871 + m_ireco1 = m_ereco.inx_left();
  872 + m_ireco2 = m_ereco.inx_right();
839 873  
840 874 // Set weighting factors for bi-linear interpolation
841   - m_wgt1 = m_etrue.wgt_left() * m_emeasured.wgt_left();
842   - m_wgt2 = m_etrue.wgt_left() * m_emeasured.wgt_right();
843   - m_wgt3 = m_etrue.wgt_right() * m_emeasured.wgt_left();
844   - m_wgt4 = m_etrue.wgt_right() * m_emeasured.wgt_right();
  875 + m_wgt1 = m_etrue.wgt_left() * m_ereco.wgt_left();
  876 + m_wgt2 = m_etrue.wgt_left() * m_ereco.wgt_right();
  877 + m_wgt3 = m_etrue.wgt_right() * m_ereco.wgt_left();
  878 + m_wgt4 = m_etrue.wgt_right() * m_ereco.wgt_right();
845 879  
846 880 } // endif: update requested
847 881  
... ... @@ -851,33 +885,21 @@ void GCTAEdispRmf::update(const double&amp; logEsrc, const double&amp; logEobs) const
851 885  
852 886  
853 887 /***********************************************************************//**
854   - * @brief Compute ebounds_obs vector
855   - *
856   - * @param[in] theta Offset angle in camera system (radians). Not used.
857   - * @param[in] phi Azimuth angle in camera system (radians). Not used.
858   - * @param[in] zenith Zenith angle in Earth system (radians). Not used.
859   - * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
  888 + * @brief Compute m_ereco_bounds vector
860 889 ***************************************************************************/
861   -void GCTAEdispRmf::compute_ebounds_obs(const double& theta,
862   - const double& phi,
863   - const double& zenith,
864   - const double& azimuth) const
  890 +void GCTAEdispRmf::compute_ereco_bounds(void) const
865 891 {
866 892 // Clear boundaries
867   - m_ebounds_obs.clear();
  893 + m_ereco_bounds.clear();
868 894  
869 895 // Loop over Etrue
870 896 for (int i = 0; i < m_rmf.ntrue(); ++i) {
871 897  
872 898 // Get true photon energy
873   - double logEsrc = m_rmf.etrue().elogmean(i).log10TeV();
874   -
875   - // Set true energy
876   - GEnergy etrue;
877   - etrue.log10TeV(logEsrc);
  899 + GEnergy etrue = m_rmf.etrue().elogmean(i);
878 900  
879 901 // Add ebounds to vector
880   - m_ebounds_obs.push_back(m_rmf.emeasured(etrue));
  902 + m_ereco_bounds.push_back(m_rmf.emeasured(etrue));
881 903  
882 904 } // endfor: looped over true photon energies
883 905  
... ... @@ -887,33 +909,21 @@ void GCTAEdispRmf::compute_ebounds_obs(const double&amp; theta,
887 909  
888 910  
889 911 /***********************************************************************//**
890   - * @brief Compute ebounds_src vector
891   - *
892   - * @param[in] theta Offset angle in camera system (radians). Not used.
893   - * @param[in] phi Azimuth angle in camera system (radians). Not used.
894   - * @param[in] zenith Zenith angle in Earth system (radians). Not used.
895   - * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
  912 + * @brief Compute m_etrue_bounds vector
896 913 ***************************************************************************/
897   -void GCTAEdispRmf::compute_ebounds_src(const double& theta,
898   - const double& phi,
899   - const double& zenith,
900   - const double& azimuth) const
  914 +void GCTAEdispRmf::compute_etrue_bounds(void) const
901 915 {
902 916 // Clear boundaries
903   - m_ebounds_src.clear();
  917 + m_etrue_bounds.clear();
904 918  
905 919 // Loop over Eobs
906 920 for (int i = 0; i < m_rmf.nmeasured(); ++i) {
907 921  
908   - // Get observed photon energy
909   - double logEobs = m_rmf.emeasured().elogmean(i).log10TeV();
910   -
911   - // Set measured energy
912   - GEnergy emeasured;
913   - emeasured.log10TeV(logEobs);
  922 + // Get reconstructed energy
  923 + GEnergy ereco = m_rmf.emeasured().elogmean(i);
914 924  
915 925 // Add ebounds to vector
916   - m_ebounds_src.push_back(m_rmf.emeasured(emeasured));
  926 + m_etrue_bounds.push_back(m_rmf.emeasured(ereco));
917 927  
918 928 } // endfor: looped over measured photon energies
919 929  
... ...
inst/cta/test/test_CTA.cpp
  View file @2ec38a7
... ... @@ -958,7 +958,7 @@ void TestGCTAResponse::test_response_edisp_RMF(void)
958 958 test_value(edisp2(GEnergy(30.0,"TeV"), GEnergy(1.0,"TeV")), 0.0);
959 959  
960 960 // Test that diagonal element is non-zero
961   - test_value(edisp2(GEnergy(30.0,"TeV"), GEnergy(30.0,"TeV")), 1.0785985662216e-07);
  961 + test_value(edisp2(GEnergy(30.0,"TeV"), GEnergy(30.0,"TeV")), 1.07978131144085e-07);
962 962  
963 963 // Test if non-diagonal element (above diagonal) is zero
964 964 test_value(edisp2(GEnergy(1.0,"TeV"), GEnergy(30.0,"TeV")), 0.0);
... ... @@ -967,13 +967,13 @@ void TestGCTAResponse::test_response_edisp_RMF(void)
967 967 GEnergy etrue(1.0, "TeV");
968 968 GEnergy emin(0.1, "TeV");
969 969 GEnergy emax(10.0, "TeV");
970   - test_value(edisp2.prob_erecobin(emin, emax, etrue, 0.0), 1.00002576227557);
971   - test_value(edisp2.prob_erecobin(etrue, emax, etrue, 0.0), 0.42065218271377);
972   - test_value(edisp2.prob_erecobin(emin, etrue, etrue, 0.0), 0.43096191010834);
  970 + test_value(edisp2.prob_erecobin(emin, emax, etrue, 0.0), 1.00105708071354);
  971 + test_value(edisp2.prob_erecobin(etrue, emax, etrue, 0.0), 0.50215659368534);
  972 + test_value(edisp2.prob_erecobin(emin, etrue, etrue, 0.0), 0.498805373808917);
973 973  
974 974 // Test mc method
975 975 GRan ran;
976   - test_value(edisp2.mc(ran, etrue).TeV(), 0.930126476703161);
  976 + test_value(edisp2.mc(ran, etrue).TeV(), 0.923744959981026);
977 977  
978 978 // Test normalisation
979 979 test_edisp_integration(edisp2, 0.5, 10.0);
... ... @@ -1006,8 +1006,8 @@ void TestGCTAResponse::test_response_edisp_2D(void)
1006 1006 test_value(edisp2.ereco_bounds(GEnergy(10.0,"TeV")).emax().TeV(), 13.2011115155942);
1007 1007  
1008 1008 // Test ebounds_src method
1009   - test_value(edisp2.etrue_bounds(GEnergy(10.0,"TeV")).emin().TeV(), 7.7813808018878);
1010   - test_value(edisp2.etrue_bounds(GEnergy(10.0,"TeV")).emax().TeV(), 15.8059295273666);
  1009 + test_value(edisp2.etrue_bounds(GEnergy(10.0,"TeV")).emin().TeV(), 7.75156710957461);
  1010 + test_value(edisp2.etrue_bounds(GEnergy(10.0,"TeV")).emax().TeV(), 15.6834030551103);
1011 1011  
1012 1012 // Test if non-diagonal element (below diagonal) is zero
1013 1013 test_value(edisp2(GEnergy(30.0,"TeV"), GEnergy(1.0,"TeV")), 0.0);
... ... @@ -1028,7 +1028,7 @@ void TestGCTAResponse::test_response_edisp_2D(void)
1028 1028  
1029 1029 // Test mc method
1030 1030 GRan ran;
1031   - test_value(edisp2.mc(ran, etrue).TeV(), 0.872387491302767);
  1031 + test_value(edisp2.mc(ran, etrue).TeV(), 0.991300452406315);
1032 1032  
1033 1033 // Test normalisation
1034 1034 test_edisp_integration(edisp2);
... ...