Simon Bonnefoy / gammalib

Commit b3e3958fc5e3c04332acc1d1c2a228e2bc167a83

Authored by Jürgen Knödlseder
1 parent 00c9c201

Rename GCTAEdisp::ebounds_obs() and GCTAEdisp::ebounds_src() to GCTAEdisp::ereco… 

…_bounds() and GCTAEdisp::etrue_bounds() (#2704)
Inline Side-by-side
Showing 16 changed files with 216 additions and 195 deletions
ChangeLog
  View file @b3e3958
... ... @@ -3,6 +3,7 @@
3 3 * Version 1.6.0 released
4 4 ========================
5 5  
  6 + Revise GCTAEdisp class interface (#2704)
6 7 GCTAEdisp returns now energy dispersion values per MeV (#2704)
7 8 Protect GCTAInstDir from access to invalid information (#2694)
8 9 Implement general CTA background model (#2694)
... ...
... ... @@ -138,6 +138,8 @@ The following methods have been removed:
138 138  
139 139 The following methods have been renamed:
140 140 - GPhoton::mc_id() from GPhoton::mcid()
  141 +- GCTAEdisp::ereco_bounds() from GCTAEdisp::ebounds_obs()
  142 +- GCTAEdisp::etrue_bounds() from GCTAEdisp::ebounds_src()
141 143  
142 144 The arguments for the following methods have been changed:
143 145 - GCTAOnOffObservation(GCTAObservation&, GModelSpatial&, GEbounds&, GEbounds&, GSkyRegions&, GSkyRegions&, bool&)
... ... @@ -310,6 +312,10 @@ assurance the coherence between the various classes. In particular, no re-
310 312 normalization is needed anymore to digest the H.E.S.S. energy dispersion, and
311 313 also the CTA energy dispersion is handled correctly (#2704).
312 314  
  315 +The GCTAEdisp class now uses GEnergy objects instead of log10(E) to deal with
  316 +true and reconstructed energies. This avoids computational overhead when going
  317 +from one unit to another (#2704).
  318 +
313 319 The GCTAInstDir class now assures that only valid information can be accessed.
314 320 Since the class can hold a sky direction and detector coordinates, flags were
315 321 added that signal what information actually exists in the class. If a non-
... ...
doc/source/admin/release_history/1.6.rst
  View file @b3e3958
... ... @@ -46,7 +46,8 @@ Improvements
46 46 ------------
47 47  
48 48 * [`2704 <https://cta-redmine.irap.omp.eu/issues/2704>`_] -
49   - ``GCTAEdisp`` returns now energy dispersion values per MeV
  49 + ``GCTAEdisp`` returns now energy dispersion values per MeV.
  50 + The ``GCTAEdisp`` interface was revised.
50 51 * [`2694 <https://cta-redmine.irap.omp.eu/issues/2694>`_] -
51 52 Protect ``GCTAInstDir`` from access to invalid information
52 53 * [`2677 <https://cta-redmine.irap.omp.eu/issues/2677>`_] -
... ...
inst/cta/include/GCTAEdisp.hpp
  View file @b3e3958
... ... @@ -77,16 +77,16 @@ public:
77 77 const double& phi = 0.0,
78 78 const double& zenith = 0.0,
79 79 const double& azimuth = 0.0) const = 0;
80   - virtual GEbounds ebounds_obs(const double& logEsrc,
81   - const double& theta = 0.0,
82   - const double& phi = 0.0,
83   - const double& zenith = 0.0,
84   - const double& azimuth = 0.0) const = 0;
85   - virtual GEbounds ebounds_src(const double& logEobs,
86   - const double& theta = 0.0,
87   - const double& phi = 0.0,
88   - const double& zenith = 0.0,
89   - const double& azimuth = 0.0) const = 0;
  80 + virtual GEbounds ereco_bounds(const GEnergy& etrue,
  81 + const double& theta = 0.0,
  82 + const double& phi = 0.0,
  83 + const double& zenith = 0.0,
  84 + const double& azimuth = 0.0) const = 0;
  85 + virtual GEbounds etrue_bounds(const GEnergy& ereco,
  86 + const double& theta = 0.0,
  87 + const double& phi = 0.0,
  88 + const double& zenith = 0.0,
  89 + const double& azimuth = 0.0) const = 0;
90 90 virtual double prob_erecobin(const GEnergy& ereco_min,
91 91 const GEnergy& ereco_max,
92 92 const GEnergy& etrue,
... ...
inst/cta/include/GCTAEdisp2D.hpp
  View file @b3e3958
... ... @@ -125,16 +125,16 @@ public:
125 125 const double& phi = 0.0,
126 126 const double& zenith = 0.0,
127 127 const double& azimuth = 0.0) const;
128   - GEbounds ebounds_obs(const double& logEsrc,
129   - const double& theta = 0.0,
130   - const double& phi = 0.0,
131   - const double& zenith = 0.0,
132   - const double& azimuth = 0.0) const;
133   - GEbounds ebounds_src(const double& logEobs,
134   - const double& theta = 0.0,
135   - const double& phi = 0.0,
136   - const double& zenith = 0.0,
137   - const double& azimuth = 0.0) const;
  128 + GEbounds ereco_bounds(const GEnergy& etrue,
  129 + const double& theta = 0.0,
  130 + const double& phi = 0.0,
  131 + const double& zenith = 0.0,
  132 + const double& azimuth = 0.0) const;
  133 + GEbounds etrue_bounds(const GEnergy& ereco,
  134 + const double& theta = 0.0,
  135 + const double& phi = 0.0,
  136 + const double& zenith = 0.0,
  137 + const double& azimuth = 0.0) const;
138 138 double prob_erecobin(const GEnergy& ereco_min,
139 139 const GEnergy& ereco_max,
140 140 const GEnergy& etrue,
... ... @@ -248,8 +248,8 @@ private:
248 248 mutable bool m_ebounds_src_computed;
249 249 mutable double m_last_theta_obs;
250 250 mutable double m_last_theta_src;
251   - mutable double m_last_logEsrc;
252   - mutable double m_last_logEobs;
  251 + mutable GEnergy m_last_etrue;
  252 + mutable GEnergy m_last_ereco;
253 253 mutable int m_index_obs;
254 254 mutable int m_index_src;
255 255 mutable double m_max_edisp;
... ...
inst/cta/include/GCTAEdispPerfTable.hpp
  View file @b3e3958
... ... @@ -94,16 +94,16 @@ public:
94 94 const double& phi = 0.0,
95 95 const double& zenith = 0.0,
96 96 const double& azimuth = 0.0) const;
97   - GEbounds ebounds_obs(const double& logEsrc,
98   - const double& theta = 0.0,
99   - const double& phi = 0.0,
100   - const double& zenith = 0.0,
101   - const double& azimuth = 0.0) const;
102   - GEbounds ebounds_src(const double& logEobs,
103   - const double& theta = 0.0,
104   - const double& phi = 0.0,
105   - const double& zenith = 0.0,
106   - const double& azimuth = 0.0) const;
  97 + GEbounds ereco_bounds(const GEnergy& etrue,
  98 + const double& theta = 0.0,
  99 + const double& phi = 0.0,
  100 + const double& zenith = 0.0,
  101 + const double& azimuth = 0.0) const;
  102 + GEbounds etrue_bounds(const GEnergy& ereco,
  103 + const double& theta = 0.0,
  104 + const double& phi = 0.0,
  105 + const double& zenith = 0.0,
  106 + const double& azimuth = 0.0) const;
107 107 double prob_erecobin(const GEnergy& ereco_min,
108 108 const GEnergy& ereco_max,
109 109 const GEnergy& etrue,
... ...
inst/cta/include/GCTAEdispRmf.hpp
  View file @b3e3958
... ... @@ -89,16 +89,16 @@ public:
89 89 const double& phi = 0.0,
90 90 const double& zenith = 0.0,
91 91 const double& azimuth = 0.0) const;
92   - GEbounds ebounds_obs(const double& logEsrc,
93   - const double& theta = 0.0,
94   - const double& phi = 0.0,
95   - const double& zenith = 0.0,
96   - const double& azimuth = 0.0) const;
97   - GEbounds ebounds_src(const double& logEobs,
98   - const double& theta = 0.0,
99   - const double& phi = 0.0,
100   - const double& zenith = 0.0,
101   - const double& azimuth = 0.0) const;
  92 + GEbounds ereco_bounds(const GEnergy& etrue,
  93 + const double& theta = 0.0,
  94 + const double& phi = 0.0,
  95 + const double& zenith = 0.0,
  96 + const double& azimuth = 0.0) const;
  97 + GEbounds etrue_bounds(const GEnergy& ereco,
  98 + const double& theta = 0.0,
  99 + const double& phi = 0.0,
  100 + const double& zenith = 0.0,
  101 + const double& azimuth = 0.0) const;
102 102 double prob_erecobin(const GEnergy& ereco_min,
103 103 const GEnergy& ereco_max,
104 104 const GEnergy& etrue,
... ... @@ -149,8 +149,8 @@ private:
149 149 mutable double m_max_edisp;
150 150 mutable double m_last_theta_obs;
151 151 mutable double m_last_theta_src;
152   - mutable double m_last_logEsrc;
153   - mutable double m_last_logEobs;
  152 + mutable GEnergy m_last_etrue_ereco_bounds;
  153 + mutable GEnergy m_last_ereco_etrue_bounds;
154 154 mutable int m_index_obs;
155 155 mutable int m_index_src;
156 156 mutable bool m_ebounds_obs_computed;
... ...
inst/cta/pyext/GCTAEdisp.i
  View file @b3e3958
... ... @@ -62,16 +62,16 @@ public:
62 62 const double& phi = 0.0,
63 63 const double& zenith = 0.0,
64 64 const double& azimuth = 0.0) const = 0;
65   - virtual GEbounds ebounds_obs(const double& logEsrc,
66   - const double& theta = 0.0,
67   - const double& phi = 0.0,
68   - const double& zenith = 0.0,
69   - const double& azimuth = 0.0) const = 0;
70   - virtual GEbounds ebounds_src(const double& logEobs,
71   - const double& theta = 0.0,
72   - const double& phi = 0.0,
73   - const double& zenith = 0.0,
74   - const double& azimuth = 0.0) const = 0;
  65 + virtual GEbounds ereco_bounds(const GEnergy& etrue,
  66 + const double& theta = 0.0,
  67 + const double& phi = 0.0,
  68 + const double& zenith = 0.0,
  69 + const double& azimuth = 0.0) const = 0;
  70 + virtual GEbounds etrue_bounds(const GEnergy& ereco,
  71 + const double& theta = 0.0,
  72 + const double& phi = 0.0,
  73 + const double& zenith = 0.0,
  74 + const double& azimuth = 0.0) const = 0;
75 75 virtual double prob_erecobin(const GEnergy& ereco_min,
76 76 const GEnergy& ereco_max,
77 77 const GEnergy& etrue,
... ...
inst/cta/pyext/GCTAEdisp2D.i
  View file @b3e3958
... ... @@ -63,16 +63,16 @@ public:
63 63 const double& phi = 0.0,
64 64 const double& zenith = 0.0,
65 65 const double& azimuth = 0.0) const;
66   - GEbounds ebounds_obs(const double& logEsrc,
67   - const double& theta = 0.0,
68   - const double& phi = 0.0,
69   - const double& zenith = 0.0,
70   - const double& azimuth = 0.0) const;
71   - GEbounds ebounds_src(const double& logEobs,
72   - const double& theta = 0.0,
73   - const double& phi = 0.0,
74   - const double& zenith = 0.0,
75   - const double& azimuth = 0.0) const;
  66 + GEbounds ereco_bounds(const GEnergy& etrue,
  67 + const double& theta = 0.0,
  68 + const double& phi = 0.0,
  69 + const double& zenith = 0.0,
  70 + const double& azimuth = 0.0) const;
  71 + GEbounds etrue_bounds(const GEnergy& ereco,
  72 + const double& theta = 0.0,
  73 + const double& phi = 0.0,
  74 + const double& zenith = 0.0,
  75 + const double& azimuth = 0.0) const;
76 76 double prob_erecobin(const GEnergy& ereco_min,
77 77 const GEnergy& ereco_max,
78 78 const GEnergy& etrue,
... ...
inst/cta/pyext/GCTAEdispPerfTable.i
  View file @b3e3958
... ... @@ -63,16 +63,16 @@ public:
63 63 const double& phi = 0.0,
64 64 const double& zenith = 0.0,
65 65 const double& azimuth = 0.0) const;
66   - GEbounds ebounds_obs(const double& logEsrc,
67   - const double& theta = 0.0,
68   - const double& phi = 0.0,
69   - const double& zenith = 0.0,
70   - const double& azimuth = 0.0) const;
71   - GEbounds ebounds_src(const double& logEobs,
72   - const double& theta = 0.0,
73   - const double& phi = 0.0,
74   - const double& zenith = 0.0,
75   - const double& azimuth = 0.0) const;
  66 + GEbounds ereco_bounds(const GEnergy& etrue,
  67 + const double& theta = 0.0,
  68 + const double& phi = 0.0,
  69 + const double& zenith = 0.0,
  70 + const double& azimuth = 0.0) const;
  71 + GEbounds etrue_bounds(const GEnergy& ereco,
  72 + const double& theta = 0.0,
  73 + const double& phi = 0.0,
  74 + const double& zenith = 0.0,
  75 + const double& azimuth = 0.0) const;
76 76 double prob_erecobin(const GEnergy& ereco_min,
77 77 const GEnergy& ereco_max,
78 78 const GEnergy& etrue,
... ...
inst/cta/pyext/GCTAEdispRmf.i
  View file @b3e3958
... ... @@ -63,16 +63,16 @@ public:
63 63 const double& phi = 0.0,
64 64 const double& zenith = 0.0,
65 65 const double& azimuth = 0.0) const;
66   - GEbounds ebounds_obs(const double& logEsrc,
67   - const double& theta = 0.0,
68   - const double& phi = 0.0,
69   - const double& zenith = 0.0,
70   - const double& azimuth = 0.0) const;
71   - GEbounds ebounds_src(const double& logEobs,
72   - const double& theta = 0.0,
73   - const double& phi = 0.0,
74   - const double& zenith = 0.0,
75   - const double& azimuth = 0.0) const;
  66 + GEbounds ereco_bounds(const GEnergy& etrue,
  67 + const double& theta = 0.0,
  68 + const double& phi = 0.0,
  69 + const double& zenith = 0.0,
  70 + const double& azimuth = 0.0) const;
  71 + GEbounds etrue_bounds(const GEnergy& ereco,
  72 + const double& theta = 0.0,
  73 + const double& phi = 0.0,
  74 + const double& zenith = 0.0,
  75 + const double& azimuth = 0.0) const;
76 76 double prob_erecobin(const GEnergy& ereco_min,
77 77 const GEnergy& ereco_max,
78 78 const GEnergy& etrue,
... ...
inst/cta/src/GCTAEdisp2D.cpp
  View file @b3e3958
... ... @@ -496,7 +496,7 @@ GEnergy GCTAEdisp2D::mc(GRan&amp; ran,
496 496 GEnergy ereco;
497 497  
498 498 // Get boundaries for observed energy
499   - GEbounds ebounds = ebounds_obs(etrue.log10TeV(), theta, phi, zenith, azimuth);
  499 + GEbounds ebounds = ereco_bounds(etrue, theta, phi, zenith, azimuth);
500 500 double emin = ebounds.emin().log10TeV();
501 501 double emax = ebounds.emax().log10TeV();
502 502  
... ... @@ -576,7 +576,7 @@ GEnergy GCTAEdisp2D::mc(GRan&amp; ran,
576 576 /***********************************************************************//**
577 577 * @brief Return observed energy interval that contains the energy dispersion.
578 578 *
579   - * @param[in] logEsrc Log10 of the true photon energy (\f$\log_{10}\f$ TeV).
  579 + * @param[in] etrue True photon energy.
580 580 * @param[in] theta Offset angle in camera system (radians).
581 581 * @param[in] phi Azimuth angle in camera system (radians). Not used.
582 582 * @param[in] zenith Zenith angle in Earth system (radians). Not used.
... ... @@ -587,11 +587,11 @@ GEnergy GCTAEdisp2D::mc(GRan&amp; ran,
587 587 * dispersion becomes negligible for a given true energy @p logEsrc and
588 588 * offset angle @p theta.
589 589 ***************************************************************************/
590   -GEbounds GCTAEdisp2D::ebounds_obs(const double& logEsrc,
591   - const double& theta,
592   - const double& phi,
593   - const double& zenith,
594   - const double& azimuth) const
  590 +GEbounds GCTAEdisp2D::ereco_bounds(const GEnergy& etrue,
  591 + const double& theta,
  592 + const double& phi,
  593 + const double& zenith,
  594 + const double& azimuth) const
595 595 {
596 596 // Make sure that energy dispersion is online
597 597 fetch();
... ... @@ -602,18 +602,21 @@ GEbounds GCTAEdisp2D::ebounds_obs(const double&amp; logEsrc,
602 602 // Set computation flag
603 603 m_ebounds_obs_computed = true;
604 604 m_last_theta_obs = theta;
605   - m_last_logEsrc = -30.0; // force update
  605 + m_last_etrue.log10TeV(-30.0); // force update
606 606  
607 607 // Compute ebounds_obs
608 608 compute_ebounds_obs(theta, phi, zenith, azimuth);
609 609  
610 610 }
611 611  
612   - // Search index only if logEsrc has changed
613   - if (logEsrc != m_last_logEsrc) {
  612 + // Search index only if true photon energy has changed
  613 + if (etrue != m_last_etrue) {
  614 +
  615 + // Store last true photon energy
  616 + m_last_etrue = etrue;
614 617  
615   - // Store last log(Esrc)
616   - m_last_logEsrc = logEsrc;
  618 + // Get true energy in TeV
  619 + double etrue_TeV = etrue.TeV();
617 620  
618 621 // Find right index with bisection
619 622 int low = 0;
... ... @@ -621,7 +624,7 @@ GEbounds GCTAEdisp2D::ebounds_obs(const double&amp; logEsrc,
621 624 while ((high-low) > 1) {
622 625 int mid = (low+high) / 2;
623 626 double e = m_edisp.axis_lo(m_inx_etrue, mid);
624   - if (logEsrc < std::log10(e)) {
  627 + if (etrue_TeV < e) {
625 628 high = mid;
626 629 }
627 630 else {
... ... @@ -642,7 +645,7 @@ GEbounds GCTAEdisp2D::ebounds_obs(const double&amp; logEsrc,
642 645 /***********************************************************************//**
643 646 * @brief Return true energy interval that contains the energy dispersion.
644 647 *
645   - * @param[in] logEobs Log10 of the observed event energy (\f$\log_{10}\f$ TeV).
  648 + * @param[in] ereco Reconstructed event energy.
646 649 * @param[in] theta Offset angle in camera system (radians).
647 650 * @param[in] phi Azimuth angle in camera system (radians). Not used.
648 651 * @param[in] zenith Zenith angle in Earth system (radians). Not used.
... ... @@ -653,11 +656,11 @@ GEbounds GCTAEdisp2D::ebounds_obs(const double&amp; logEsrc,
653 656 * dispersion becomes negligible for a given observed energy @p logEobs and
654 657 * offset angle @p theta.
655 658 ***************************************************************************/
656   -GEbounds GCTAEdisp2D::ebounds_src(const double& logEobs,
657   - const double& theta,
658   - const double& phi,
659   - const double& zenith,
660   - const double& azimuth) const
  659 +GEbounds GCTAEdisp2D::etrue_bounds(const GEnergy& ereco,
  660 + const double& theta,
  661 + const double& phi,
  662 + const double& zenith,
  663 + const double& azimuth) const
661 664 {
662 665 // Make sure that energy dispersion is online
663 666 fetch();
... ... @@ -668,18 +671,21 @@ GEbounds GCTAEdisp2D::ebounds_src(const double&amp; logEobs,
668 671 // Set computation flag
669 672 m_ebounds_src_computed = true;
670 673 m_last_theta_src = theta;
671   - m_last_logEobs = -30.0; // force update
  674 + m_last_ereco.log10TeV(-30.0); // force update
672 675  
673 676 // Compute ebounds_src
674 677 compute_ebounds_src(theta, phi, zenith, azimuth);
675 678  
676 679 }
677 680  
678   - // Search index only if logEobs has changed
679   - if (logEobs != m_last_logEobs) {
  681 + // Search index only if reconstructed event energy has changed
  682 + if (ereco != m_last_ereco) {
  683 +
  684 + // Store reconstructed event energy
  685 + m_last_ereco = ereco;
680 686  
681   - // Store observed energy
682   - m_last_logEobs = logEobs;
  687 + // Get reconstructed event energy in TeV
  688 + double ereco_TeV = ereco.TeV();
683 689  
684 690 // Find right index with bisection
685 691 int low = 0;
... ... @@ -687,7 +693,7 @@ GEbounds GCTAEdisp2D::ebounds_src(const double&amp; logEobs,
687 693 while ((high-low) > 1) {
688 694 int mid = (low+high) / 2;
689 695 double e = m_edisp.axis_lo(m_inx_etrue, mid);
690   - if (logEobs < std::log10(e)) {
  696 + if (ereco_TeV < e) {
691 697 high = mid;
692 698 }
693 699 else {
... ... @@ -1046,8 +1052,8 @@ void GCTAEdisp2D::init_members(void)
1046 1052 m_ebounds_src_computed = false;
1047 1053 m_last_theta_obs = -1.0;
1048 1054 m_last_theta_src = -1.0;
1049   - m_last_logEsrc = -30.0;
1050   - m_last_logEobs = -30.0;
  1055 + m_last_etrue.log10TeV(-30.0);
  1056 + m_last_ereco.log10TeV(-30.0);
1051 1057 m_index_obs = 0;
1052 1058 m_index_src = 0;
1053 1059 m_max_edisp = 0.0;
... ... @@ -1086,8 +1092,8 @@ void GCTAEdisp2D::copy_members(const GCTAEdisp2D&amp; edisp)
1086 1092 m_ebounds_src_computed = edisp.m_ebounds_src_computed;
1087 1093 m_last_theta_obs = edisp.m_last_theta_obs;
1088 1094 m_last_theta_src = edisp.m_last_theta_src;
1089   - m_last_logEsrc = edisp.m_last_logEsrc;
1090   - m_last_logEobs = edisp.m_last_logEobs;
  1095 + m_last_etrue = edisp.m_last_etrue;
  1096 + m_last_ereco = edisp.m_last_ereco;
1091 1097 m_index_obs = edisp.m_index_obs;
1092 1098 m_index_src = edisp.m_index_src;
1093 1099 m_max_edisp = edisp.m_max_edisp;
... ... @@ -1532,7 +1538,7 @@ void GCTAEdisp2D::normalize_table(void)
1532 1538 double sum = 0.0;
1533 1539  
1534 1540 // Get integration boundaries
1535   - GEbounds ebounds = ebounds_obs(logEsrc, theta);
  1541 + GEbounds ebounds = ereco_bounds(etrue, theta);
1536 1542  
1537 1543 // Loop over all energy intervals
1538 1544 for (int i = 0; i < ebounds.size(); ++i) {
... ...
inst/cta/src/GCTAEdispPerfTable.cpp
  View file @b3e3958
... ... @@ -360,7 +360,7 @@ GEnergy GCTAEdispPerfTable::mc(GRan&amp; ran,
360 360 /***********************************************************************//**
361 361 * @brief Return observed energy interval that contains the energy dispersion.
362 362 *
363   - * @param[in] logEsrc Log10 of the true photon energy (\f$\log_{10}\f$ TeV).
  363 + * @param[in] etrue True photon energy.
364 364 * @param[in] theta Offset angle in camera system (radians). Not used.
365 365 * @param[in] phi Azimuth angle in camera system (radians). Not used.
366 366 * @param[in] zenith Zenith angle in Earth system (radians). Not used.
... ... @@ -372,23 +372,26 @@ GEnergy GCTAEdispPerfTable::mc(GRan&amp; ran,
372 372 * band is set to \f$\pm 5 \times \sigma\f$, where \f$\sigma\f$ is the
373 373 * Gaussian width of the energy dispersion.
374 374 ***************************************************************************/
375   -GEbounds GCTAEdispPerfTable::ebounds_obs(const double& logEsrc,
376   - const double& theta,
377   - const double& phi,
378   - const double& zenith,
379   - const double& azimuth) const
  375 +GEbounds GCTAEdispPerfTable::ereco_bounds(const GEnergy& etrue,
  376 + const double& theta,
  377 + const double& phi,
  378 + const double& zenith,
  379 + const double& azimuth) const
380 380 {
381 381 // Set energy band constant
382 382 const double number_of_sigmas = 5.0;
383 383  
  384 + // Get log10 of true photon energy in TeV
  385 + double etrue_log10TeV = etrue.log10TeV();
  386 +
384 387 // Get energy dispersion sigma
385   - double sigma = m_logE.interpolate(logEsrc, m_sigma);
  388 + double sigma = m_logE.interpolate(etrue_log10TeV, m_sigma);
386 389  
387 390 // Compute energy boundaries
388 391 GEnergy emin;
389 392 GEnergy emax;
390   - emin.log10TeV(logEsrc - number_of_sigmas * sigma);
391   - emax.log10TeV(logEsrc + number_of_sigmas * sigma);
  393 + emin.log10TeV(etrue_log10TeV - number_of_sigmas * sigma);
  394 + emax.log10TeV(etrue_log10TeV + number_of_sigmas * sigma);
392 395  
393 396 // Return energy boundaries
394 397 return (GEbounds(emin, emax));
... ... @@ -398,7 +401,7 @@ GEbounds GCTAEdispPerfTable::ebounds_obs(const double&amp; logEsrc,
398 401 /***********************************************************************//**
399 402 * @brief Return true energy interval that contains the energy dispersion.
400 403 *
401   - * @param[in] logEobs Log10 of the observed event energy (\f$\log_{10}\f$ TeV).
  404 + * @param[in] ereco Reconstructed event energy.
402 405 * @param[in] theta Offset angle in camera system (radians). Not used.
403 406 * @param[in] phi Azimuth angle in camera system (radians). Not used.
404 407 * @param[in] zenith Zenith angle in Earth system (radians). Not used.
... ... @@ -410,23 +413,26 @@ GEbounds GCTAEdispPerfTable::ebounds_obs(const double&amp; logEsrc,
410 413 * This band is set to \f$\pm 5 \times \sigma\f$, where \f$\sigma\f$ is the
411 414 * Gaussian width of the energy dispersion.
412 415 ***************************************************************************/
413   -GEbounds GCTAEdispPerfTable::ebounds_src(const double& logEobs,
414   - const double& theta,
415   - const double& phi,
416   - const double& zenith,
417   - const double& azimuth) const
  416 +GEbounds GCTAEdispPerfTable::etrue_bounds(const GEnergy& ereco,
  417 + const double& theta,
  418 + const double& phi,
  419 + const double& zenith,
  420 + const double& azimuth) const
418 421 {
419 422 // Set energy band constant
420 423 const double number_of_sigmas = 5.0;
421 424  
  425 + // Get log10 of reconstructed event energy in TeV
  426 + double ereco_log10TeV = ereco.log10TeV();
  427 +
422 428 // Get energy dispersion sigma
423   - double sigma = m_logE.interpolate(logEobs, m_sigma);
  429 + double sigma = m_logE.interpolate(ereco_log10TeV, m_sigma);
424 430  
425 431 // Compute energy boundaries
426 432 GEnergy emin;
427 433 GEnergy emax;
428   - emin.log10TeV(logEobs - number_of_sigmas * sigma);
429   - emax.log10TeV(logEobs + number_of_sigmas * sigma);
  434 + emin.log10TeV(ereco_log10TeV - number_of_sigmas * sigma);
  435 + emax.log10TeV(ereco_log10TeV + number_of_sigmas * sigma);
430 436  
431 437 // Return energy boundaries
432 438 return (GEbounds(emin, emax));
... ...
inst/cta/src/GCTAEdispRmf.cpp
  View file @b3e3958
... ... @@ -307,7 +307,7 @@ GEnergy GCTAEdispRmf::mc(GRan&amp; ran,
307 307 GEnergy ereco;
308 308  
309 309 // Get boundaries for rejection method
310   - GEbounds ebounds = ebounds_obs(etrue.log10TeV(), theta, phi, zenith, azimuth);
  310 + GEbounds ebounds = ereco_bounds(etrue, theta, phi, zenith, azimuth);
311 311 double emin = ebounds.emin().log10TeV();
312 312 double emax = ebounds.emax().log10TeV();
313 313  
... ... @@ -333,7 +333,7 @@ GEnergy GCTAEdispRmf::mc(GRan&amp; ran,
333 333 /***********************************************************************//**
334 334 * @brief Return observed energy interval that contains the energy dispersion.
335 335 *
336   - * @param[in] logEsrc Log10 of the true photon energy (\f$\log_{10}\f$ TeV).
  336 + * @param[in] etrue True photon energy.
337 337 * @param[in] theta Offset angle in camera system (radians). Not used.
338 338 * @param[in] phi Azimuth angle in camera system (radians). Not used.
339 339 * @param[in] zenith Zenith angle in Earth system (radians). Not used.
... ... @@ -341,13 +341,13 @@ GEnergy GCTAEdispRmf::mc(GRan&amp; ran,
341 341 * @return Reconstructed energy boundaries.
342 342 *
343 343 * Returns the band of observed energies outside of which the energy
344   - * dispersion becomes negligible for a given true energy @p logEsrc.
  344 + * dispersion becomes negligible for a given true energy @p etrue.
345 345 ***************************************************************************/
346   -GEbounds GCTAEdispRmf::ebounds_obs(const double& logEsrc,
347   - const double& theta,
348   - const double& phi,
349   - const double& zenith,
350   - const double& azimuth) const
  346 +GEbounds GCTAEdispRmf::ereco_bounds(const GEnergy& etrue,
  347 + const double& theta,
  348 + const double& phi,
  349 + const double& zenith,
  350 + const double& azimuth) const
351 351 {
352 352 // Compute only if parameters changed
353 353 if (!m_ebounds_obs_computed || theta != m_last_theta_obs) {
... ... @@ -355,6 +355,7 @@ GEbounds GCTAEdispRmf::ebounds_obs(const double&amp; logEsrc,
355 355 // Store last theta value
356 356 m_ebounds_obs_computed = true;
357 357 m_last_theta_obs = theta;
  358 + m_last_etrue_ereco_bounds.log10TeV(-30.0); // force update
358 359  
359 360 // Compute ebounds_obs
360 361 compute_ebounds_obs(theta, phi, zenith, azimuth);
... ... @@ -362,18 +363,21 @@ GEbounds GCTAEdispRmf::ebounds_obs(const double&amp; logEsrc,
362 363 }
363 364  
364 365 // Search index only if logEsrc has changed
365   - if (logEsrc != m_last_logEsrc) {
  366 + if (etrue != m_last_etrue_ereco_bounds) {
366 367  
367 368 // Store last photon energy
368   - m_last_logEsrc = logEsrc;
  369 + m_last_etrue_ereco_bounds = etrue;
  370 +
  371 + // Get true photon energy in TeV
  372 + double etrue_TeV = etrue.TeV();
369 373  
370 374 // Find right index with bisection
371 375 int low = 0;
372 376 int high = m_ebounds_obs.size() - 1;
373 377 while ((high-low) > 1) {
374 378 int mid = (low+high) / 2;
375   - double e = m_rmf.etrue().emin(mid).log10TeV();
376   - if (logEsrc < e) {
  379 + double e = m_rmf.etrue().emin(mid).TeV();
  380 + if (etrue_TeV < e) {
377 381 high = mid;
378 382 }
379 383 else {
... ... @@ -394,7 +398,7 @@ GEbounds GCTAEdispRmf::ebounds_obs(const double&amp; logEsrc,
394 398 /***********************************************************************//**
395 399 * @brief Return true energy interval that contains the energy dispersion.
396 400 *
397   - * @param[in] logEobs Log10 of the observed event energy (\f$\log_{10}\f$ TeV).
  401 + * @param[in] ereco Reconstructed event energy.
398 402 * @param[in] theta Offset angle in camera system (radians). Not used.
399 403 * @param[in] phi Azimuth angle in camera system (radians). Not used.
400 404 * @param[in] zenith Zenith angle in Earth system (radians). Not used.
... ... @@ -404,11 +408,11 @@ GEbounds GCTAEdispRmf::ebounds_obs(const double&amp; logEsrc,
404 408 * Returns the band of true photon energies outside of which the energy
405 409 * dispersion becomes negligible for a given observed energy @p logEobs.
406 410 ***************************************************************************/
407   -GEbounds GCTAEdispRmf::ebounds_src(const double& logEobs,
408   - const double& theta,
409   - const double& phi,
410   - const double& zenith,
411   - const double& azimuth) const
  411 +GEbounds GCTAEdispRmf::etrue_bounds(const GEnergy& ereco,
  412 + const double& theta,
  413 + const double& phi,
  414 + const double& zenith,
  415 + const double& azimuth) const
412 416 {
413 417 // Compute only if parameters changed
414 418 if (!m_ebounds_src_computed || theta != m_last_theta_src) {
... ... @@ -416,25 +420,29 @@ GEbounds GCTAEdispRmf::ebounds_src(const double&amp; logEobs,
416 420 // Store last value
417 421 m_ebounds_src_computed = true;
418 422 m_last_theta_src = theta;
  423 + m_last_ereco_etrue_bounds.log10TeV(-30.0); // force update
419 424  
420 425 // Compute ebounds_obs
421 426 compute_ebounds_src(theta, phi, zenith, azimuth);
422 427  
423 428 }
424 429  
425   - // Search index only if logEobs has changed
426   - if (logEobs != m_last_logEobs) {
  430 + // Search index only if ereco has changed
  431 + if (ereco != m_last_ereco_etrue_bounds) {
427 432  
428 433 // Store last value
429   - m_last_logEobs = logEobs;
  434 + m_last_ereco_etrue_bounds = ereco;
  435 +
  436 + // Get reconstructed energy in TeV
  437 + double ereco_TeV = ereco.TeV();
430 438  
431 439 // Find right index with bisection
432 440 int low = 0;
433 441 int high = m_ebounds_src.size() - 1;
434 442 while ((high-low) > 1) {
435 443 int mid = (low+high) / 2;
436   - double e = m_rmf.emeasured().emin(mid).log10TeV();
437   - if (logEobs < e) {
  444 + double e = m_rmf.emeasured().emin(mid).TeV();
  445 + if (ereco_TeV < e) {
438 446 high = mid;
439 447 }
440 448 else {
... ... @@ -627,8 +635,8 @@ void GCTAEdispRmf::init_members(void)
627 635 m_max_edisp = 0.0;
628 636 m_last_theta_obs = -1.0;
629 637 m_last_theta_src = -1.0;
630   - m_last_logEsrc = -30.0;
631   - m_last_logEobs = -30.0;
  638 + m_last_etrue_ereco_bounds.log10TeV(-30.0);
  639 + m_last_ereco_etrue_bounds.log10TeV(-30.0);
632 640 m_index_obs = 0;
633 641 m_index_src = 0;
634 642 m_ebounds_obs_computed = false;
... ... @@ -668,17 +676,17 @@ void GCTAEdispRmf::copy_members(const GCTAEdispRmf&amp; edisp)
668 676 m_wgt4 = edisp.m_wgt4;
669 677  
670 678 // Copy Monte Carlo cache
671   - m_max_edisp = edisp.m_max_edisp;
672   - m_last_theta_obs = edisp.m_last_theta_obs;
673   - m_last_theta_src = edisp.m_last_theta_src;
674   - m_last_logEsrc = edisp.m_last_logEsrc;
675   - m_last_logEobs = edisp.m_last_logEobs;
676   - m_index_obs = edisp.m_index_obs;
677   - m_index_src = edisp.m_index_src;
678   - m_ebounds_obs_computed = edisp.m_ebounds_obs_computed;
679   - m_ebounds_src_computed = edisp.m_ebounds_src_computed;
680   - m_ebounds_obs = edisp.m_ebounds_obs;
681   - m_ebounds_src = edisp.m_ebounds_src;
  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;
682 690  
683 691 // Return
684 692 return;
... ... @@ -721,11 +729,8 @@ void GCTAEdispRmf::set_matrix(void)
721 729 // Normalize matrix elements
722 730 for (int itrue = 0; itrue < rows; ++itrue) {
723 731  
724   - // Get true photon energy
725   - double logEsrc = m_rmf.etrue().elogmean(itrue).log10TeV();
726   -
727 732 // Get integration boundaries
728   - GEbounds ebounds = ebounds_obs(logEsrc, 0.0);
  733 + GEbounds ebounds = ereco_bounds(m_rmf.etrue().elogmean(itrue), 0.0);
729 734  
730 735 // Integrate contributions over energy boundaries
731 736 double sum = 0.0;
... ...
inst/cta/src/GCTAResponseIrf.cpp
  View file @b3e3958
... ... @@ -450,8 +450,7 @@ double GCTAResponseIrf::nroi(const GModelSky&amp; model,
450 450 if (use_edisp()) {
451 451  
452 452 // Retrieve true energy boundaries
453   - double obsLogEng = obsEng.log10TeV();
454   - GEbounds ebounds = edisp()->ebounds_src(obsLogEng);
  453 + GEbounds ebounds = edisp()->etrue_bounds(obsEng);
455 454  
456 455 // Loop over all boundaries
457 456 for (int i = 0; i < ebounds.size(); ++i) {
... ... @@ -550,8 +549,7 @@ GEbounds GCTAResponseIrf::ebounds(const GEnergy&amp; obsEnergy) const
550 549  
551 550 // If energy dispersion is available then set the energy boundaries
552 551 if (edisp() != NULL) {
553   - double obsLogEng = obsEnergy.log10TeV();
554   - ebounds = edisp()->ebounds_src(obsLogEng); // Requires TeV
  552 + ebounds = edisp()->etrue_bounds(obsEnergy);
555 553 }
556 554  
557 555 // Return energy boundaries
... ...
inst/cta/test/test_CTA.cpp
  View file @b3e3958
... ... @@ -890,12 +890,12 @@ void TestGCTAResponse::test_response_edisp_PerfTable(void)
890 890 "Test filename after loading perfromance table");
891 891  
892 892 // Test ebounds_obs method
893   - test_value(edisp2.ebounds_obs(-1.0).emin().GeV(), 31.6795127349936);
894   - test_value(edisp2.ebounds_obs(-1.0).emax().GeV(), 315.661420794325);
  893 + test_value(edisp2.ereco_bounds(GEnergy(0.1,"TeV")).emin().GeV(), 31.6795127349936);
  894 + test_value(edisp2.ereco_bounds(GEnergy(0.1,"TeV")).emax().GeV(), 315.661420794325);
895 895  
896 896 // Test ebounds_src method
897   - test_value(edisp2.ebounds_src(-1.0).emin().GeV(), 31.6795127349936);
898   - test_value(edisp2.ebounds_src(-1.0).emax().GeV(), 315.661420794325);
  897 + test_value(edisp2.etrue_bounds(GEnergy(0.1,"TeV")).emin().GeV(), 31.6795127349936);
  898 + test_value(edisp2.etrue_bounds(GEnergy(0.1,"TeV")).emax().GeV(), 315.661420794325);
899 899  
900 900 // Test if non-diagonal element (below diagonal) is zero
901 901 test_value(edisp2(GEnergy(30.0,"TeV"), GEnergy(10.0,"TeV")), 0.0);
... ... @@ -947,12 +947,12 @@ void TestGCTAResponse::test_response_edisp_RMF(void)
947 947 "Test filename after loading RMF file");
948 948  
949 949 // Test ebounds_obs method
950   - test_value(edisp2.ebounds_obs(1.0).emin().TeV(), 5.75439929962158);
951   - test_value(edisp2.ebounds_obs(1.0).emax().TeV(), 15.8489322662354);
  950 + test_value(edisp2.ereco_bounds(GEnergy(10.0,"TeV")).emin().TeV(), 5.75439929962158);
  951 + test_value(edisp2.ereco_bounds(GEnergy(10.0,"TeV")).emax().TeV(), 15.8489322662354);
952 952  
953 953 // Test ebounds_src method
954   - test_value(edisp2.ebounds_src(1.0).emin().TeV(), 5.75439929962158);
955   - test_value(edisp2.ebounds_src(1.0).emax().TeV(), 15.8489322662354);
  954 + test_value(edisp2.etrue_bounds(GEnergy(10.0,"TeV")).emin().TeV(), 5.75439929962158);
  955 + test_value(edisp2.etrue_bounds(GEnergy(10.0,"TeV")).emax().TeV(), 15.8489322662354);
956 956  
957 957 // Test if non-diagonal element (below diagonal) is zero
958 958 test_value(edisp2(GEnergy(30.0,"TeV"), GEnergy(1.0,"TeV")), 0.0);
... ... @@ -1002,12 +1002,12 @@ void TestGCTAResponse::test_response_edisp_2D(void)
1002 1002 "Test filename after loading response table file");
1003 1003  
1004 1004 // Test ebounds_obs method
1005   - test_value(edisp2.ebounds_obs(1.0).emin().TeV(), 6.82816152379073);
1006   - test_value(edisp2.ebounds_obs(1.0).emax().TeV(), 13.2011115155942);
  1005 + test_value(edisp2.ereco_bounds(GEnergy(10.0,"TeV")).emin().TeV(), 6.82816152379073);
  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.ebounds_src(1.0).emin().TeV(), 7.7813808018878);
1010   - test_value(edisp2.ebounds_src(1.0).emax().TeV(), 15.8059295273666);
  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);
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);
... ... @@ -1427,10 +1427,9 @@ void TestGCTAResponse::test_response_edisp_integration(const GCTAResponseIrf&amp; rs
1427 1427  
1428 1428 // Compute log10 of true energy
1429 1429 GEnergy etrue(e_src,"TeV");
1430   - double log10_e_src = etrue.log10TeV();
1431 1430  
1432   - // Retrieve boundaries in observed energy
1433   - GEbounds ebounds = rsp.edisp()->ebounds_obs(log10_e_src);
  1431 + // Retrieve observed energy boundaries
  1432 + GEbounds ebounds = rsp.edisp()->ereco_bounds(etrue);
1434 1433 double emin = ebounds.emin().MeV();
1435 1434 double emax = ebounds.emax().MeV();
1436 1435  
... ... @@ -1480,10 +1479,9 @@ void TestGCTAResponse::test_edisp_integration(const GCTAEdisp&amp; edisp,
1480 1479  
1481 1480 // Compute log10 of true energy
1482 1481 GEnergy etrue(e_src,"TeV");
1483   - double log10_e_src = etrue.log10TeV();
1484 1482  
1485   - // Retrieve boundaries in observed energy
1486   - GEbounds ebounds = edisp.ebounds_obs(log10_e_src);
  1483 + // Retrieve observed energy boundaries
  1484 + GEbounds ebounds = edisp.ereco_bounds(etrue);
1487 1485  
1488 1486 // Skip this energy if the boundaries are empty
1489 1487 if (ebounds.is_empty()) {
... ...