Simon Bonnefoy / gammalib

Commit 6517c4861da56f1f19c34d17fe670327997c310c

Authored by Jürgen Knödlseder
1 parent 57acac30

Change GCTAEdisp2D::operator() and GCTAResponseIrf::edisp() interface (#2704)

Inline Side-by-side
Showing 20 changed files with 231 additions and 219 deletions
... ... @@ -142,10 +142,13 @@ The following methods have been renamed:
142 142 The arguments for the following methods have been changed:
143 143 - GCTAOnOffObservation(GCTAObservation&, GModelSpatial&, GEbounds&, GEbounds&, GSkyRegions&, GSkyRegions&, bool&)
144 144 - GCTAModelRadial::mc(GRan&)
  145 +- GCTAEdisp::operator(GEnergy&, GEnergy&, double&, double&, double&, double&)
  146 +- GCTAResponseIrf::edisp(GEnergy&, GEnergy&, double&, double&, double&, double&)
145 147 - GLATAeff(GFits&, std::string&)
146 148 - GLATEdisp(GFits&, std::string&)
147 149 - GLATPsf(GFits&, std::string&)
148 150  
  151 +
149 152 The return value of the following methods has been changed:
150 153 - none
151 154  
... ...
inst/cta/include/GCTAEdisp.hpp
  View file @6517c48
... ... @@ -55,12 +55,12 @@ public:
55 55 virtual ~GCTAEdisp(void);
56 56  
57 57 // Pure virtual operators
58   - virtual double operator()(const double& logEobs,
59   - const double& logEsrc,
60   - const double& theta = 0.0,
61   - const double& phi = 0.0,
62   - const double& zenith = 0.0,
63   - const double& azimuth = 0.0) const = 0;
  58 + virtual double operator()(const GEnergy& ereco,
  59 + const GEnergy& etrue,
  60 + const double& theta = 0.0,
  61 + const double& phi = 0.0,
  62 + const double& zenith = 0.0,
  63 + const double& azimuth = 0.0) const = 0;
64 64  
65 65 // Operators
66 66 GCTAEdisp& operator=(const GCTAEdisp& edisp);
... ...
inst/cta/include/GCTAEdisp2D.hpp
  View file @6517c48
... ... @@ -31,6 +31,7 @@
31 31 #include <string>
32 32 #include "GFilename.hpp"
33 33 #include "GFunction.hpp"
  34 +#include "GEnergy.hpp"
34 35 #include "GCTAEdisp.hpp"
35 36 #include "GCTAResponseTable.hpp"
36 37  
... ... @@ -105,12 +106,12 @@ public:
105 106  
106 107 // Operators
107 108 GCTAEdisp2D& operator=(const GCTAEdisp2D& edisp);
108   - double operator()(const double& logEobs,
109   - const double& logEsrc,
110   - const double& theta = 0.0,
111   - const double& phi = 0.0,
112   - const double& zenith = 0.0,
113   - const double& azimuth = 0.0) const;
  109 + double operator()(const GEnergy& ereco,
  110 + const GEnergy& etrue,
  111 + const double& theta = 0.0,
  112 + const double& phi = 0.0,
  113 + const double& zenith = 0.0,
  114 + const double& azimuth = 0.0) const;
114 115  
115 116 // Implemented methods
116 117 void clear(void);
... ... @@ -215,15 +216,15 @@ private:
215 216 class edisp_kern : public GFunction {
216 217 public:
217 218 edisp_kern(const GCTAEdisp2D* parent,
218   - const double& logEsrc,
  219 + const GEnergy& etrue,
219 220 const double& theta) :
220 221 m_parent(parent),
221   - m_logEsrc(logEsrc),
  222 + m_etrue(etrue),
222 223 m_theta(theta) { }
223 224 double eval(const double& logEobs);
224 225 protected:
225 226 const GCTAEdisp2D* m_parent; //!< Pointer to parent class
226   - double m_logEsrc; //!< True photon energy
  227 + GEnergy m_etrue; //!< True photon energy
227 228 double m_theta; //!< Offset angle
228 229 };
229 230  
... ...
inst/cta/include/GCTAEdispPerfTable.hpp
  View file @6517c48
... ... @@ -75,12 +75,12 @@ public:
75 75  
76 76 // Operators
77 77 GCTAEdispPerfTable& operator=(const GCTAEdispPerfTable& psf);
78   - double operator()(const double& logEobs,
79   - const double& logEsrc,
80   - const double& theta = 0.0,
81   - const double& phi = 0.0,
82   - const double& zenith = 0.0,
83   - const double& azimuth = 0.0) const;
  78 + double operator()(const GEnergy& ereco,
  79 + const GEnergy& etrue,
  80 + const double& theta = 0.0,
  81 + const double& phi = 0.0,
  82 + const double& zenith = 0.0,
  83 + const double& azimuth = 0.0) const;
84 84  
85 85 // Implemented pure virtual methods
86 86 void clear(void);
... ...
inst/cta/include/GCTAEdispRmf.hpp
  View file @6517c48
... ... @@ -70,12 +70,12 @@ public:
70 70  
71 71 // Operators
72 72 GCTAEdispRmf& operator=(const GCTAEdispRmf& edisp);
73   - double operator()(const double& logEobs,
74   - const double& logEsrc,
75   - const double& theta = 0.0,
76   - const double& phi = 0.0,
77   - const double& zenith = 0.0,
78   - const double& azimuth = 0.0) const;
  73 + double operator()(const GEnergy& ereco,
  74 + const GEnergy& etrue,
  75 + const double& theta = 0.0,
  76 + const double& phi = 0.0,
  77 + const double& zenith = 0.0,
  78 + const double& azimuth = 0.0) const;
79 79  
80 80 // Implemented pure virtual methods
81 81 void clear(void);
... ...
inst/cta/include/GCTAResponseIrf.hpp
  View file @6517c48
... ... @@ -137,12 +137,12 @@ public:
137 137 const double& zenith,
138 138 const double& azimuth,
139 139 const double& srcLogEng) const;
140   - double edisp(const GEnergy& obsEng,
  140 + double edisp(const GEnergy& ereco,
  141 + const GEnergy& etrue,
141 142 const double& theta,
142 143 const double& phi,
143 144 const double& zenith,
144   - const double& azimuth,
145   - const double& srcLogEng) const;
  145 + const double& azimuth) const;
146 146 double nroi(const GModelSky& model,
147 147 const GEnergy& srcEng,
148 148 const GTime& srcTime,
... ...
inst/cta/pyext/GCTAEdisp.i
  View file @6517c48
... ... @@ -43,12 +43,12 @@ public:
43 43 virtual ~GCTAEdisp(void);
44 44  
45 45 // Pure virtual operators
46   - virtual double operator()(const double& logEobs,
47   - const double& logEsrc,
48   - const double& theta = 0.0,
49   - const double& phi = 0.0,
50   - const double& zenith = 0.0,
51   - const double& azimuth = 0.0) const = 0;
  46 + virtual double operator()(const GEnergy& ereco,
  47 + const GEnergy& etrue,
  48 + const double& theta = 0.0,
  49 + const double& phi = 0.0,
  50 + const double& zenith = 0.0,
  51 + const double& azimuth = 0.0) const = 0;
52 52  
53 53 // Pure virtual methods
54 54 virtual void clear(void) = 0;
... ...
inst/cta/pyext/GCTAEdisp2D.i
  View file @6517c48
... ... @@ -44,12 +44,12 @@ public:
44 44 virtual ~GCTAEdisp2D(void);
45 45  
46 46 // Operators
47   - double operator()(const double& logEobs,
48   - const double& logEsrc,
49   - const double& theta = 0.0,
50   - const double& phi = 0.0,
51   - const double& zenith = 0.0,
52   - const double& azimuth = 0.0) const;
  47 + double operator()(const GEnergy& ereco,
  48 + const GEnergy& etrue,
  49 + const double& theta = 0.0,
  50 + const double& phi = 0.0,
  51 + const double& zenith = 0.0,
  52 + const double& azimuth = 0.0) const;
53 53  
54 54 // Implemented methods
55 55 void clear(void);
... ...
inst/cta/pyext/GCTAEdispPerfTable.i
  View file @6517c48
... ... @@ -44,12 +44,12 @@ public:
44 44 virtual ~GCTAEdispPerfTable(void);
45 45  
46 46 // Operators
47   - double operator()(const double& logEobs,
48   - const double& logEsrc,
49   - const double& theta = 0.0,
50   - const double& phi = 0.0,
51   - const double& zenith = 0.0,
52   - const double& azimuth = 0.0) const;
  47 + double operator()(const GEnergy& ereco,
  48 + const GEnergy& etrue,
  49 + const double& theta = 0.0,
  50 + const double& phi = 0.0,
  51 + const double& zenith = 0.0,
  52 + const double& azimuth = 0.0) const;
53 53  
54 54 // Implemented pure virtual methods
55 55 void clear(void);
... ...
inst/cta/pyext/GCTAEdispRmf.i
  View file @6517c48
... ... @@ -44,12 +44,12 @@ public:
44 44 virtual ~GCTAEdispRmf(void);
45 45  
46 46 // Operators
47   - double operator()(const double& logEobs,
48   - const double& logEsrc,
49   - const double& theta = 0.0,
50   - const double& phi = 0.0,
51   - const double& zenith = 0.0,
52   - const double& azimuth = 0.0) const;
  47 + double operator()(const GEnergy& ereco,
  48 + const GEnergy& etrue,
  49 + const double& theta = 0.0,
  50 + const double& phi = 0.0,
  51 + const double& zenith = 0.0,
  52 + const double& azimuth = 0.0) const;
53 53  
54 54 // Implemented pure virtual methods
55 55 void clear(void);
... ...
inst/cta/pyext/GCTAResponseIrf.i
  View file @6517c48
... ... @@ -105,12 +105,12 @@ public:
105 105 const double& zenith,
106 106 const double& azimuth,
107 107 const double& srcLogEng) const;
108   - double edisp(const GEnergy& obsEng,
  108 + double edisp(const GEnergy& ereco,
  109 + const GEnergy& etrue,
109 110 const double& theta,
110 111 const double& phi,
111 112 const double& zenith,
112   - const double& azimuth,
113   - const double& srcLogEng) const;
  113 + const double& azimuth) const;
114 114 double npsf(const GSkyDir& srcDir,
115 115 const double& srcLogEng,
116 116 const GTime& srcTime,
... ...
inst/cta/src/GCTACubeEdisp.cpp
  View file @6517c48
1 1 /***************************************************************************
2 2 * GCTACubeEdisp.cpp - CTA cube analysis energy dispersion class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2016-2017 by Michael Mayer *
  4 + * copyright (C) 2016-2018 by Michael Mayer *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -994,9 +994,9 @@ void GCTACubeEdisp::fill_cube(const GCTAObservation&amp; obs,
994 994  
995 995 // Add energy dispersion cube value
996 996 m_cube(pixel, imap) += rsp->edisp(Eobs,
  997 + m_energies[iebin],
997 998 theta, 0.0,
998   - 0.0, 0.0,
999   - logEsrc) * weight;
  999 + 0.0, 0.0) * weight;
1000 1000  
1001 1001 } // endfor: looped over migration bins
1002 1002  
... ...
inst/cta/src/GCTAEdisp2D.cpp
  View file @6517c48
... ... @@ -180,8 +180,8 @@ GCTAEdisp2D&amp; GCTAEdisp2D::operator=(const GCTAEdisp2D&amp; edisp)
180 180 /***********************************************************************//**
181 181 * @brief Return energy dispersion in units of MeV\f$^{-1}\f$
182 182 *
183   - * @param[in] logEobs Log10 of the observed energy (\f$\log_{10}\f$ TeV).
184   - * @param[in] logEsrc Log10 of the true photon energy (\f$\log_{10}\f$ TeV).
  183 + * @param[in] ereco Reconstructed photon energy.
  184 + * @param[in] etrue True photon energy.
185 185 * @param[in] theta Offset angle in camera system (radians).
186 186 * @param[in] phi Azimuth angle in camera system (radians). Not used.
187 187 * @param[in] zenith Zenith angle in Earth system (radians). Not used.
... ... @@ -201,12 +201,12 @@ GCTAEdisp2D&amp; GCTAEdisp2D::operator=(const GCTAEdisp2D&amp; edisp)
201 201 * \f$\log_{10} E_{\rm true}\f$ is the base 10 logarithm of the true energy
202 202 * in TeV and \f$\theta\f$ is the offset angle in radians.
203 203 ***************************************************************************/
204   -double GCTAEdisp2D::operator()(const double& logEobs,
205   - const double& logEsrc,
206   - const double& theta,
207   - const double& phi,
208   - const double& zenith,
209   - const double& azimuth) const
  204 +double GCTAEdisp2D::operator()(const GEnergy& ereco,
  205 + const GEnergy& etrue,
  206 + const double& theta,
  207 + const double& phi,
  208 + const double& zenith,
  209 + const double& azimuth) const
210 210 {
211 211 // Make sure that energy dispersion is online
212 212 fetch();
... ... @@ -214,20 +214,24 @@ double GCTAEdisp2D::operator()(const double&amp; logEobs,
214 214 // Initalize edisp
215 215 double edisp = 0.0;
216 216  
  217 + // Get log10 of true photon energy
  218 + double logEsrc = etrue.log10TeV();
  219 +
217 220 // Continue only if logEsrc and theta are in validity range
218 221 if ((logEsrc >= m_logEsrc_min) && (logEsrc <= m_logEsrc_max) &&
219 222 (theta >= m_theta_min) && (theta <= m_theta_max)) {
220 223  
221   - // Compute Eobs/Esrc
222   - double EobsOverEsrc = std::exp((logEobs-logEsrc) * gammalib::ln10);
  224 + // Compute migration (migra=Ereco/Etrue
  225 + //double EobsOverEsrc = std::exp((logEobs-logEsrc) * gammalib::ln10);
  226 + double migra = ereco.MeV() / etrue.MeV();
223 227  
224 228 // Continue only if migration is in validity range
225   - if ((EobsOverEsrc >= m_migra_min) && (EobsOverEsrc <= m_migra_max)) {
  229 + if ((migra >= m_migra_min) && (migra <= m_migra_max)) {
226 230  
227 231 // Setup argument vector
228 232 double arg[3];
229 233 arg[m_inx_etrue] = logEsrc;
230   - arg[m_inx_migra] = EobsOverEsrc;
  234 + arg[m_inx_migra] = migra;
231 235 arg[m_inx_theta] = theta;
232 236  
233 237 // Compute edisp
... ... @@ -238,17 +242,11 @@ double GCTAEdisp2D::operator()(const double&amp; logEobs,
238 242 edisp = 0.0;
239 243 }
240 244  
241   - // If energy dispersion is positive then compute true energy
242   - // in MeV and divide the energy dispersion value by this
243   - // quantity to get the energy dispersion per MeV
  245 + // If energy dispersion is positive then divide the energy
  246 + // dispersion value by true photon energy to get the energy
  247 + // dispersion per MeV
244 248 else {
245   -
246   - // Compute true energy
247   - double etrue = std::pow(10.0, logEsrc+6.0);
248   -
249   - // Compute energy dispersion per MeV
250   - edisp /= etrue;
251   -
  249 + edisp /= etrue.MeV();
252 250 }
253 251  
254 252 } // endif: migration was in valid range
... ... @@ -494,6 +492,11 @@ GEnergy GCTAEdisp2D::mc(GRan&amp; ran,
494 492 // Make sure that energy dispersion is online
495 493 fetch();
496 494  
  495 + // Initialise energies
  496 + GEnergy etrue;
  497 + GEnergy ereco;
  498 + etrue.log10TeV(logEsrc);
  499 +
497 500 // Get boundaries for observed energy
498 501 GEbounds ebounds = ebounds_obs(logEsrc, theta, phi, zenith, azimuth);
499 502 double emin = ebounds.emin().log10TeV();
... ... @@ -541,7 +544,8 @@ GEnergy GCTAEdisp2D::mc(GRan&amp; ran,
541 544 int zeros = 0;
542 545 do {
543 546 logEobs = emin + ewidth * ran.uniform();
544   - f = operator()(logEobs, logEsrc, theta, phi, zenith, azimuth);
  547 + ereco.log10TeV(logEobs);
  548 + f = operator()(ereco, etrue, theta, phi, zenith, azimuth);
545 549 if (f == 0.0) {
546 550 zeros++;
547 551 if (zeros > max_subsequent_zeros) {
... ... @@ -568,12 +572,8 @@ GEnergy GCTAEdisp2D::mc(GRan&amp; ran,
568 572  
569 573 } // endwhile:
570 574  
571   - // Set observed energy
572   - GEnergy energy;
573   - energy.log10TeV(logEobs);
574   -
575 575 // Return energy
576   - return energy;
  576 + return ereco;
577 577 }
578 578  
579 579  
... ... @@ -1244,6 +1244,10 @@ void GCTAEdisp2D::compute_ebounds_src(const double&amp; theta,
1244 1244 const double& zenith,
1245 1245 const double& azimuth) const
1246 1246 {
  1247 + // Initialise energies
  1248 + GEnergy etrue;
  1249 + GEnergy ereco;
  1250 +
1247 1251 // Clear ebounds_src vector
1248 1252 m_ebounds_src.clear();
1249 1253  
... ... @@ -1260,13 +1264,13 @@ void GCTAEdisp2D::compute_ebounds_src(const double&amp; theta,
1260 1264 // Set Eobs
1261 1265 double Eobs = std::sqrt(m_edisp.axis_hi(m_inx_etrue,iobs) *
1262 1266 m_edisp.axis_lo(m_inx_etrue,iobs));
1263   - double logEobs = std::log10(Eobs);
  1267 + ereco.TeV(Eobs);
1264 1268  
1265 1269 // Initialise results
1266   - double logEsrcMin = 0.0;
1267   - double logEsrcMax = 0.0;
1268   - bool minFound = false;
1269   - bool maxFound = false;
  1270 + GEnergy etrue_min;
  1271 + GEnergy etrue_max;
  1272 + bool minFound = false;
  1273 + bool maxFound = false;
1270 1274  
1271 1275 // Find minimum boundary
1272 1276 for (int imigra = 0; imigra < n_migras; ++imigra) {
... ... @@ -1281,15 +1285,15 @@ void GCTAEdisp2D::compute_ebounds_src(const double&amp; theta,
1281 1285 }
1282 1286  
1283 1287 // Compute log10 of true energy
1284   - double logEsrc = logEobs - std::log10(migra);
  1288 + etrue = ereco / migra;
1285 1289  
1286 1290 // Get matrix term
1287   - double edisp = operator()(logEobs, logEsrc, theta);
  1291 + double edisp = operator()(ereco, etrue, theta);
1288 1292  
1289 1293 // Find first non-negligible matrix term
1290 1294 if (edisp >= eps) {
1291 1295 maxFound = true;
1292   - logEsrcMax = logEsrc;
  1296 + etrue_max = etrue;
1293 1297 break;
1294 1298 }
1295 1299  
... ... @@ -1308,15 +1312,15 @@ void GCTAEdisp2D::compute_ebounds_src(const double&amp; theta,
1308 1312 }
1309 1313  
1310 1314 // Compute log10 of true energy
1311   - double logEsrc = logEobs - std::log10(migra);
  1315 + etrue = ereco / migra;
1312 1316  
1313 1317 // Get matrix term
1314   - double edisp = operator()(logEobs, logEsrc, theta);
  1318 + double edisp = operator()(ereco, etrue, theta);
1315 1319  
1316 1320 // Find first non-negligible matrix term
1317 1321 if (edisp >= eps) {
1318 1322 minFound = true;
1319   - logEsrcMin = logEsrc;
  1323 + etrue_min = etrue;
1320 1324 break;
1321 1325 }
1322 1326  
... ... @@ -1325,18 +1329,12 @@ void GCTAEdisp2D::compute_ebounds_src(const double&amp; theta,
1325 1329 // If we did not find boundaries then set the interval to
1326 1330 // a zero interval for safety
1327 1331 if (!minFound || !maxFound) {
1328   - logEsrcMin = 0.0;
1329   - logEsrcMax = 0.0;
  1332 + etrue_min.clear();
  1333 + etrue_max.clear();
1330 1334 }
1331 1335  
1332   - // Set energy boundaries
1333   - GEnergy emin;
1334   - GEnergy emax;
1335   - emin.log10TeV(logEsrcMin);
1336   - emax.log10TeV(logEsrcMax);
1337   -
1338 1336 // Add energy boundaries to vector
1339   - m_ebounds_src.push_back(GEbounds(emin, emax));
  1337 + m_ebounds_src.push_back(GEbounds(etrue_min, etrue_max));
1340 1338  
1341 1339 } // endfor : loopend over observed energy
1342 1340  
... ... @@ -1530,6 +1528,10 @@ void GCTAEdisp2D::normalize_table(void)
1530 1528 double emax = std::log10(m_edisp.axis_hi(m_inx_etrue,i_etrue));
1531 1529 double logEsrc = 0.5*(emin+emax);
1532 1530  
  1531 + // Set true energy
  1532 + GEnergy etrue;
  1533 + etrue.log10TeV(logEsrc);
  1534 +
1533 1535 // Initialise integration
1534 1536 double sum = 0.0;
1535 1537  
... ... @@ -1544,7 +1546,7 @@ void GCTAEdisp2D::normalize_table(void)
1544 1546 emax = ebounds.emax(i).log10MeV();
1545 1547  
1546 1548 // Setup integration function
1547   - edisp_kern integrand(this, logEsrc, theta);
  1549 + edisp_kern integrand(this, etrue, theta);
1548 1550 GIntegral integral(&integrand);
1549 1551  
1550 1552 // Set integration precision
... ... @@ -2129,16 +2131,16 @@ GNdarray GCTAEdisp2D::gaussian_array(const double&amp; mean,
2129 2131 ***************************************************************************/
2130 2132 double GCTAEdisp2D::edisp_kern::eval(const double& logEobs)
2131 2133 {
  2134 + // Set reconstructued energy
  2135 + GEnergy ereco;
  2136 + ereco.log10MeV(logEobs);
  2137 +
2132 2138 // Get function value. Note that logEobs is in log10 MeV, while we need
2133 2139 // the value in log10 TeV.
2134   - double value = m_parent->operator()(logEobs-6.0, m_logEsrc, m_theta);
2135   -
2136   - // Compute reconstructed energy in MeV. Note that logEobs is in
2137   - // log10 MeV!
2138   - double ereco = std::pow(10.0, logEobs);
  2140 + double value = m_parent->operator()(ereco, m_etrue, m_theta);
2139 2141  
2140 2142 // Normalize energy dispersion so that the units are per log10 MeV
2141   - value *= gammalib::ln10 * ereco;
  2143 + value *= gammalib::ln10 * ereco.MeV();
2142 2144  
2143 2145 // Compile option: Check for NaN
2144 2146 #if defined(G_NAN_CHECK)
... ...
inst/cta/src/GCTAEdispPerfTable.cpp
  View file @6517c48
... ... @@ -160,8 +160,8 @@ GCTAEdispPerfTable&amp; GCTAEdispPerfTable::operator=(const GCTAEdispPerfTable&amp; edis
160 160 /***********************************************************************//**
161 161 * @brief Return energy dispersion in units of MeV\f$^{-1}\f$
162 162 *
163   - * @param[in] logEobs Log10 of the observed photon energy (\f$\log_{10}\f$ TeV).
164   - * @param[in] logEsrc Log10 of the true photon energy (\f$\log_{10}\f$ TeV).
  163 + * @param[in] ereco Reconstructed photon energy.
  164 + * @param[in] etrue True photon energy.
165 165 * @param[in] theta Offset angle in camera system (radians). Not used.
166 166 * @param[in] phi Azimuth angle in camera system (radians). Not used.
167 167 * @param[in] zenith Zenith angle in Earth system (radians). Not used.
... ... @@ -184,13 +184,17 @@ GCTAEdispPerfTable&amp; GCTAEdispPerfTable::operator=(const GCTAEdispPerfTable&amp; edis
184 184 * \f$\sigma(E_{\rm true})\f$ is the standard deviation of the energy
185 185 * dispersion that depends on the true photon energy.
186 186 ***************************************************************************/
187   -double GCTAEdispPerfTable::operator()(const double& logEobs,
188   - const double& logEsrc,
189   - const double& theta,
190   - const double& phi,
191   - const double& zenith,
192   - const double& azimuth) const
  187 +double GCTAEdispPerfTable::operator()(const GEnergy& ereco,
  188 + const GEnergy& etrue,
  189 + const double& theta,
  190 + const double& phi,
  191 + const double& zenith,
  192 + const double& azimuth) const
193 193 {
  194 + // Get log10 of true and reconstructed photon energies
  195 + double logEsrc = etrue.log10TeV();
  196 + double logEobs = ereco.log10TeV();
  197 +
194 198 // Update the parameter cache
195 199 update(logEsrc);
196 200  
... ... @@ -198,11 +202,8 @@ double GCTAEdispPerfTable::operator()(const double&amp; logEobs,
198 202 double delta = logEobs - logEsrc;
199 203 double edisp = m_par_scale * std::exp(m_par_width * delta * delta);
200 204  
201   - // Compute reconstructed energy in MeV
202   - double ereco = std::pow(10.0, logEobs+6.0);
203   -
204 205 // Compute energy dispersion per MeV
205   - edisp /= (gammalib::ln10 * ereco);
  206 + edisp /= (gammalib::ln10 * ereco.MeV());
206 207  
207 208 // Return energy dispersion
208 209 return edisp;
... ...
inst/cta/src/GCTAEdispRmf.cpp
  View file @6517c48
... ... @@ -162,8 +162,8 @@ GCTAEdispRmf&amp; GCTAEdispRmf::operator=(const GCTAEdispRmf&amp; edisp)
162 162 /***********************************************************************//**
163 163 * @brief Return energy dispersion in units of MeV\f$^{-1}\f$
164 164 *
165   - * @param[in] logEobs Log10 of the observed photon energy (\f$\log_{10}\f$ TeV).
166   - * @param[in] logEsrc Log10 of the true photon energy (\f$\log_{10}\f$ TeV).
  165 + * @param[in] ereco Reconstructed photon energy.
  166 + * @param[in] etrue True photon energy.
167 167 * @param[in] theta Offset angle in camera system (radians). Not used.
168 168 * @param[in] phi Azimuth angle in camera system (radians). Not used.
169 169 * @param[in] zenith Zenith angle in Earth system (radians). Not used.
... ... @@ -182,13 +182,17 @@ GCTAEdispRmf&amp; GCTAEdispRmf::operator=(const GCTAEdispRmf&amp; edisp)
182 182 * \f$E_{\rm reco}\f$ is the reconstructed energy in units of MeV, and
183 183 * \f$E_{\rm true}\f$ is the true energy in units of MeV.
184 184 ***************************************************************************/
185   -double GCTAEdispRmf::operator()(const double& logEobs,
186   - const double& logEsrc,
187   - const double& theta,
188   - const double& phi,
189   - const double& zenith,
190   - const double& azimuth) const
  185 +double GCTAEdispRmf::operator()(const GEnergy& ereco,
  186 + const GEnergy& etrue,
  187 + const double& theta,
  188 + const double& phi,
  189 + const double& zenith,
  190 + const double& azimuth) const
191 191 {
  192 + // Get log10 of true and reconstructed photon energies
  193 + double logEsrc = etrue.log10TeV();
  194 + double logEobs = ereco.log10TeV();
  195 +
192 196 // Update indexes and weighting factors for interpolation
193 197 update(logEobs, logEsrc);
194 198  
... ... @@ -205,7 +209,7 @@ double GCTAEdispRmf::operator()(const double&amp; logEobs,
205 209  
206 210 // If the energy dispersion is positive, return it per MeV
207 211 else {
208   - edisp /= (gammalib::ln10 * std::pow(10.0, logEobs+6.0));
  212 + edisp /= (gammalib::ln10 * ereco.MeV());
209 213 }
210 214  
211 215 // Return energy dispersion
... ... @@ -299,6 +303,11 @@ GEnergy GCTAEdispRmf::mc(GRan&amp; ran,
299 303 const double& zenith,
300 304 const double& azimuth) const
301 305 {
  306 + // Initialise energies
  307 + GEnergy etrue;
  308 + GEnergy ereco;
  309 + etrue.log10TeV(logEsrc);
  310 +
302 311 // Get boundaries for rejection method
303 312 GEbounds ebounds = ebounds_obs(logEsrc, theta, phi, zenith, azimuth);
304 313 double emin = ebounds.emin().log10TeV();
... ... @@ -312,17 +321,14 @@ GEnergy GCTAEdispRmf::mc(GRan&amp; ran,
312 321 double ftest = 1.0;
313 322 while (ftest > f) {
314 323 logEobs = emin + ewidth * ran.uniform();
315   - f = operator()(logEobs, logEsrc, theta, phi, zenith, azimuth);
  324 + ereco.log10TeV(logEobs);
  325 + f = operator()(ereco, etrue, theta, phi, zenith, azimuth);
316 326 ftest = ran.uniform() * m_max_edisp;
317 327 }
318 328 }
319 329  
320   - // Set observed energy
321   - GEnergy energy;
322   - energy.log10TeV(logEobs);
323   -
324   - // Return energy
325   - return energy;
  330 + // Return reconstructed photon energy
  331 + return ereco;
326 332 }
327 333  
328 334  
... ... @@ -492,10 +498,10 @@ double GCTAEdispRmf::prob_erecobin(const GEnergy&amp; ereco_min,
492 498 double wgt_right = m_etrue.wgt_right();
493 499  
494 500 // Initialise first and second node
495   - double x1 = logEobs_min;
496   - double f1 = this->operator()(x1, logEsrc);
497   - double x2 = 0.0;
498   - double f2 = 0.0;
  501 + GEnergy x1 = ereco_min;
  502 + double f1 = this->operator()(x1, etrue);
  503 + GEnergy x2;
  504 + double f2 = 0.0;
499 505  
500 506 // Loop over all measured energy nodes
501 507 for (int i = 0; i < m_emeasured.size(); ++i) {
... ... @@ -508,13 +514,13 @@ double GCTAEdispRmf::prob_erecobin(const GEnergy&amp; ereco_min,
508 514 // If measured energy is above maximum migration value then use
509 515 // logEobs_max as measured energy
510 516 if (m_emeasured[i] > logEobs_max) {
511   - x2 = logEobs_max;
512   - f2 = this->operator()(x2, logEsrc);
  517 + x2 = ereco_max;
  518 + f2 = this->operator()(x2, etrue);
513 519 }
514 520  
515 521 // ... otherwise use measured energy
516 522 else {
517   - x2 = m_emeasured[i];
  523 + x2.log10TeV(m_emeasured[i]);
518 524 f2 = wgt_left * m_matrix(inx_left, i) +
519 525 wgt_right * m_matrix(inx_right, i);
520 526 if (f2 < 0.0) {
... ... @@ -523,7 +529,7 @@ double GCTAEdispRmf::prob_erecobin(const GEnergy&amp; ereco_min,
523 529 }
524 530  
525 531 // Compute integral
526   - prob += 0.5 * (f1 + f2) * (x2 - x1);
  532 + prob += 0.5 * (f1 + f2) * (x2.log10TeV() - x1.log10TeV());
527 533  
528 534 // Set second node as first node
529 535 x1 = x2;
... ... @@ -538,10 +544,10 @@ double GCTAEdispRmf::prob_erecobin(const GEnergy&amp; ereco_min,
538 544  
539 545 // If last node energy is below migra_max then compute last part of
540 546 // integral up to emax
541   - if (x1 < logEobs_max) {
542   - x2 = logEobs_max;
543   - f2 = this->operator()(x2, logEsrc);
544   - prob += 0.5 * (f1 + f2) * (x2 - x1);
  547 + if (x1 < ereco_max) {
  548 + x2 = ereco_max;
  549 + f2 = this->operator()(x2, etrue);
  550 + prob += 0.5 * (f1 + f2) * (x2.log10TeV() - x1.log10TeV());
545 551 }
546 552  
547 553 // Return probability
... ...
inst/cta/src/GCTAResponseIrf.cpp
  View file @6517c48
... ... @@ -375,7 +375,7 @@ double GCTAResponseIrf::irf(const GEvent&amp; event,
375 375 if (use_edisp() && irf > 0.0) {
376 376  
377 377 // Multiply-in energy dispersion
378   - irf *= edisp(obsEng, theta, phi, zenith, azimuth, srcLogEng);
  378 + irf *= edisp(obsEng, srcEng, theta, phi, zenith, azimuth);
379 379  
380 380 } // endif: energy dispersion was available and PSF was non-zero
381 381  
... ... @@ -1804,26 +1804,23 @@ double GCTAResponseIrf::psf_delta_max(const double&amp; theta,
1804 1804 /***********************************************************************//**
1805 1805 * @brief Return energy dispersion (in units of MeV\f$^{-1}\f$)
1806 1806 *
1807   - * @param[in] obsEng Measured event energy.
  1807 + * @param[in] ereco Reconstructed event energy.
  1808 + * @param[in] etrue True photon energy.
1808 1809 * @param[in] theta Radial offset angle in camera (radians).
1809 1810 * @param[in] phi Polar angle in camera (radians).
1810 1811 * @param[in] zenith Zenith angle of telescope pointing (radians).
1811 1812 * @param[in] azimuth Azimuth angle of telescope pointing (radians).
1812   - * @param[in] srcLogEng Log10 of true photon energy (E/TeV).
1813 1813 * @return Energy dispersion (MeV\f$^{-1}\f$).
1814 1814 ***************************************************************************/
1815   -double GCTAResponseIrf::edisp(const GEnergy& obsEng,
  1815 +double GCTAResponseIrf::edisp(const GEnergy& ereco,
  1816 + const GEnergy& etrue,
1816 1817 const double& theta,
1817 1818 const double& phi,
1818 1819 const double& zenith,
1819   - const double& azimuth,
1820   - const double& srcLogEng) const
  1820 + const double& azimuth) const
1821 1821 {
1822   - // Compute log10 energy in TeV
1823   - double obsLogEng = obsEng.log10TeV();
1824   -
1825 1822 // Compute energy dispersion
1826   - double edisp = (*m_edisp)(obsLogEng, srcLogEng, theta, phi, zenith, azimuth);
  1823 + double edisp = (*m_edisp)(ereco, etrue, theta, phi, zenith, azimuth);
1827 1824  
1828 1825 // Return energy dispersion
1829 1826 return edisp;
... ... @@ -1946,7 +1943,7 @@ double GCTAResponseIrf::nirf(const GPhoton&amp; photon,
1946 1943 if (use_edisp() && nroi > 0.0) {
1947 1944  
1948 1945 // Multiply-in energy dispersion
1949   - nroi *= edisp(obsEng, theta, phi, zenith, azimuth, srcLogEng);
  1946 + nroi *= edisp(obsEng, srcEng, theta, phi, zenith, azimuth);
1950 1947  
1951 1948 } // endif: had energy dispersion
1952 1949  
... ... @@ -2396,7 +2393,6 @@ double GCTAResponseIrf::irf_radial(const GEvent&amp; event,
2396 2393 azimuth,
2397 2394 srcEng,
2398 2395 srcTime,
2399   - srcLogEng,
2400 2396 obsEng,
2401 2397 zeta,
2402 2398 lambda,
... ... @@ -2614,7 +2610,6 @@ double GCTAResponseIrf::irf_elliptical(const GEvent&amp; event,
2614 2610 azimuth,
2615 2611 srcEng,
2616 2612 srcTime,
2617   - srcLogEng,
2618 2613 obsEng,
2619 2614 rho_obs,
2620 2615 posangle_obs,
... ...
inst/cta/src/GCTAResponse_helpers.cpp
  View file @6517c48
... ... @@ -483,7 +483,7 @@ double cta_irf_radial_kern_rho::eval(const double&amp; rho)
483 483 cta_irf_radial_kern_omega integrand(m_rsp,
484 484 m_zenith,
485 485 m_azimuth,
486   - m_srcLogEng,
  486 + m_srcEng,
487 487 m_obsEng,
488 488 m_zeta,
489 489 m_lambda,
... ... @@ -565,6 +565,9 @@ double cta_irf_radial_kern_rho::eval(const double&amp; rho)
565 565 ***************************************************************************/
566 566 double cta_irf_radial_kern_omega::eval(const double& omega)
567 567 {
  568 + // Get log10(E/TeV) of true photon energy
  569 + double srcLogEng = m_srcEng.log10TeV();
  570 +
568 571 // Compute PSF offset angle [radians]
569 572 double delta = std::acos(m_cos_psf + m_sin_psf * std::cos(omega));
570 573  
... ... @@ -575,12 +578,12 @@ double cta_irf_radial_kern_omega::eval(const double&amp; omega)
575 578 double azimuth = 0.0;
576 579  
577 580 // Evaluate IRF
578   - double irf = m_rsp->aeff(offset, azimuth, m_zenith, m_azimuth, m_srcLogEng) *
579   - m_rsp->psf(delta, offset, azimuth, m_zenith, m_azimuth, m_srcLogEng);
  581 + double irf = m_rsp->aeff(offset, azimuth, m_zenith, m_azimuth, srcLogEng) *
  582 + m_rsp->psf(delta, offset, azimuth, m_zenith, m_azimuth, srcLogEng);
580 583  
581 584 // Optionally take energy dispersion into account
582 585 if (m_rsp->use_edisp() && irf > 0.0) {
583   - irf *= m_rsp->edisp(m_obsEng, offset, azimuth, m_zenith, m_azimuth, m_srcLogEng);
  586 + irf *= m_rsp->edisp(m_obsEng, m_srcEng, offset, azimuth, m_zenith, m_azimuth);
584 587 }
585 588  
586 589 // Compile option: Check for NaN/Inf
... ... @@ -840,7 +843,6 @@ double cta_irf_elliptical_kern_rho::eval(const double&amp; rho)
840 843 m_azimuth,
841 844 m_srcEng,
842 845 m_srcTime,
843   - m_srcLogEng,
844 846 m_obsEng,
845 847 m_posangle_obs,
846 848 m_omega_pnt,
... ... @@ -1033,6 +1035,9 @@ double cta_irf_elliptical_kern_omega::eval(const double&amp; omega)
1033 1035 // Continue only if model is positive
1034 1036 if (model > 0.0) {
1035 1037  
  1038 + // Get log10(E/TeV) of true photon energy
  1039 + double srcLogEng = m_srcEng.log10TeV();
  1040 +
1036 1041 // Compute Psf offset angle [radians]
1037 1042 double delta = std::acos(m_cos_psf + m_sin_psf * std::cos(omega));
1038 1043  
... ... @@ -1042,14 +1047,13 @@ double cta_irf_elliptical_kern_omega::eval(const double&amp; omega)
1042 1047 double phi = 0.0; //TODO: Implement IRF Phi dependence
1043 1048  
1044 1049 // Evaluate IRF * model
1045   - irf = m_rsp->aeff(theta, phi, m_zenith, m_azimuth, m_srcLogEng) *
1046   - m_rsp->psf(delta, theta, phi, m_zenith, m_azimuth, m_srcLogEng) *
  1050 + irf = m_rsp->aeff(theta, phi, m_zenith, m_azimuth, srcLogEng) *
  1051 + m_rsp->psf(delta, theta, phi, m_zenith, m_azimuth, srcLogEng) *
1047 1052 model;
1048 1053  
1049 1054 // Optionally take energy dispersion into account
1050 1055 if (m_rsp->use_edisp() && irf > 0.0) {
1051   - irf *= m_rsp->edisp(m_obsEng, theta, phi,
1052   - m_zenith, m_azimuth, m_srcLogEng);
  1056 + irf *= m_rsp->edisp(m_obsEng, m_srcEng, theta, phi, m_zenith, m_azimuth);
1053 1057 }
1054 1058  
1055 1059 // Compile option: Check for NaN/Inf
... ... @@ -1530,8 +1534,8 @@ double cta_irf_diffuse_kern_phi::eval(const double&amp; phi)
1530 1534  
1531 1535 // Optionally take energy dispersion into account
1532 1536 if (m_rsp->use_edisp() && irf > 0.0) {
1533   - irf *= m_rsp->edisp(m_obsEng, offset, azimuth,
1534   - m_zenith, m_azimuth, m_srcLogEng);
  1537 + irf *= m_rsp->edisp(m_obsEng, m_srcEng, offset, azimuth,
  1538 + m_zenith, m_azimuth);
1535 1539 }
1536 1540  
1537 1541 // Compile option: Check for NaN/Inf
... ...
inst/cta/src/GCTAResponse_helpers.hpp
  View file @6517c48
... ... @@ -189,7 +189,6 @@ public:
189 189 const double& azimuth,
190 190 const GEnergy& srcEng,
191 191 const GTime& srcTime,
192   - const double& srcLogEng,
193 192 const GEnergy& obsEng,
194 193 const double& zeta,
195 194 const double& lambda,
... ... @@ -202,7 +201,6 @@ public:
202 201 m_azimuth(azimuth),
203 202 m_srcEng(srcEng),
204 203 m_srcTime(srcTime),
205   - m_srcLogEng(srcLogEng),
206 204 m_obsEng(obsEng),
207 205 m_zeta(zeta),
208 206 m_cos_zeta(std::cos(zeta)),
... ... @@ -222,7 +220,6 @@ protected:
222 220 double m_azimuth; //!< Azimuth angle
223 221 GEnergy m_srcEng; //!< True photon energy
224 222 GTime m_srcTime; //!< True photon time
225   - double m_srcLogEng; //!< True photon log10 energy
226 223 GEnergy m_obsEng; //!< Measured event energy
227 224 double m_zeta; //!< Distance model centre - measured photon
228 225 double m_cos_zeta; //!< Cosine of zeta
... ... @@ -260,7 +257,7 @@ public:
260 257 cta_irf_radial_kern_omega(const GCTAResponseIrf* rsp,
261 258 const double& zenith,
262 259 const double& azimuth,
263   - const double& srcLogEng,
  260 + const GEnergy& srcEng,
264 261 const GEnergy& obsEng,
265 262 const double& zeta,
266 263 const double& lambda,
... ... @@ -273,7 +270,7 @@ public:
273 270 m_rsp(rsp),
274 271 m_zenith(zenith),
275 272 m_azimuth(azimuth),
276   - m_srcLogEng(srcLogEng),
  273 + m_srcEng(srcEng),
277 274 m_obsEng(obsEng),
278 275 m_zeta(zeta),
279 276 m_lambda(lambda),
... ... @@ -288,7 +285,7 @@ protected:
288 285 const GCTAResponseIrf* m_rsp; //!< CTA response
289 286 double m_zenith; //!< Zenith angle
290 287 double m_azimuth; //!< Azimuth angle
291   - double m_srcLogEng; //!< True photon energy
  288 + GEnergy m_srcEng; //!< True photon energy
292 289 GEnergy m_obsEng; //!< Measured event energy
293 290 double m_zeta; //!< Distance model centre - measured photon
294 291 double m_lambda; //!< Distance model centre - pointing
... ... @@ -474,7 +471,6 @@ public:
474 471 const double& azimuth,
475 472 const GEnergy& srcEng,
476 473 const GTime& srcTime,
477   - const double& srcLogEng,
478 474 const GEnergy& obsEng,
479 475 const double& rho_obs,
480 476 const double& posangle_obs,
... ... @@ -491,7 +487,6 @@ public:
491 487 m_azimuth(azimuth),
492 488 m_srcEng(srcEng),
493 489 m_srcTime(srcTime),
494   - m_srcLogEng(srcLogEng),
495 490 m_obsEng(obsEng),
496 491 m_rho_obs(rho_obs),
497 492 m_cos_rho_obs(std::cos(rho_obs)),
... ... @@ -515,7 +510,6 @@ public:
515 510 double m_azimuth; //!< Azimuth angle
516 511 GEnergy m_srcEng; //!< True photon energy
517 512 GTime m_srcTime; //!< True photon time
518   - double m_srcLogEng; //!< True photon log energy
519 513 GEnergy m_obsEng; //!< Measured event energy
520 514 double m_rho_obs; //!< Distance of model centre from measured photon
521 515 double m_cos_rho_obs; //!< Cosine of m_rho_obs
... ... @@ -558,7 +552,6 @@ public:
558 552 const double& azimuth,
559 553 const GEnergy& srcEng,
560 554 const GTime& srcTime,
561   - const double& srcLogEng,
562 555 const GEnergy& obsEng,
563 556 const double& posangle_obs,
564 557 const double& omega_pnt,
... ... @@ -573,7 +566,6 @@ public:
573 566 m_azimuth(azimuth),
574 567 m_srcEng(srcEng),
575 568 m_srcTime(srcTime),
576   - m_srcLogEng(srcLogEng),
577 569 m_obsEng(obsEng),
578 570 m_posangle_obs(posangle_obs),
579 571 m_omega_pnt(omega_pnt),
... ... @@ -590,7 +582,6 @@ public:
590 582 double m_azimuth; //!< Azimuth angle
591 583 GEnergy m_srcEng; //!< True photon energy
592 584 GTime m_srcTime; //!< True photon time
593   - double m_srcLogEng; //!< True photon log energy
594 585 GEnergy m_obsEng; //!< Measured event energy
595 586 double m_posangle_obs; //!< Measured photon position angle from model centre
596 587 double m_omega_pnt; //!< Azimuth of pointing in model system
... ... @@ -1102,7 +1093,7 @@ public:
1102 1093  
1103 1094  
1104 1095 /***********************************************************************//**
1105   - * @class cta_irf_radial_kern_omega
  1096 + * @class cta_psf_radial_kern_omega
1106 1097 *
1107 1098 * @brief Kernel for radial model azimuth angle integration
1108 1099 *
... ...
inst/cta/test/test_CTA.cpp
  View file @6517c48
... ... @@ -898,13 +898,13 @@ void TestGCTAResponse::test_response_edisp_PerfTable(void)
898 898 test_value(edisp2.ebounds_src(-1.0).emax().GeV(), 315.661420794325);
899 899  
900 900 // Test if non-diagonal element (below diagonal) is zero
901   - test_value(edisp2(std::log10(30.0), std::log10(1.0)), 0.0);
  901 + test_value(edisp2(GEnergy(30.0,"TeV"), GEnergy(10.0,"TeV")), 0.0);
902 902  
903 903 // Test that diagonal element is non-zero
904   - test_value(edisp2(std::log10(30.0),std::log10(30.0)), 2.12840930049511e-07);
  904 + test_value(edisp2(GEnergy(30.0,"TeV"), GEnergy(30.0,"TeV")), 2.12840930049511e-07);
905 905  
906 906 // Test if non-diagonal element (above diagonal) is zero
907   - test_value(edisp2(std::log10(1.0), std::log10(30.0)), 0.0);
  907 + test_value(edisp2(GEnergy(10.0,"TeV"), GEnergy(30.0,"TeV")), 0.0);
908 908  
909 909  
910 910 // Test prob_erecobin method
... ... @@ -955,13 +955,13 @@ void TestGCTAResponse::test_response_edisp_RMF(void)
955 955 test_value(edisp2.ebounds_src(1.0).emax().TeV(), 15.8489322662354);
956 956  
957 957 // Test if non-diagonal element (below diagonal) is zero
958   - test_value(edisp2(std::log10(30.0), std::log10(1.0)), 0.0);
  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(std::log10(30.0),std::log10(30.0)), 1.0785985662216e-07);
  961 + test_value(edisp2(GEnergy(30.0,"TeV"), GEnergy(30.0,"TeV")), 1.0785985662216e-07);
962 962  
963 963 // Test if non-diagonal element (above diagonal) is zero
964   - test_value(edisp2(std::log10(1.0), std::log10(30.0)), 0.0);
  964 + test_value(edisp2(GEnergy(1.0,"TeV"), GEnergy(30.0,"TeV")), 0.0);
965 965  
966 966 // Test prob_erecobin method
967 967 GEnergy etrue(1.0, "TeV");
... ... @@ -1010,13 +1010,13 @@ void TestGCTAResponse::test_response_edisp_2D(void)
1010 1010 test_value(edisp2.ebounds_src(1.0).emax().TeV(), 15.8059295273666);
1011 1011  
1012 1012 // Test if non-diagonal element (below diagonal) is zero
1013   - test_value(edisp2(std::log10(30.0), std::log10(1.0)), 0.0);
  1013 + test_value(edisp2(GEnergy(30.0,"TeV"), GEnergy(1.0,"TeV")), 0.0);
1014 1014  
1015 1015 // Test that diagonal element is non-zero
1016   - test_value(edisp2(std::log10(30.0),std::log10(30.0)), 1.74200329848547e-07);
  1016 + test_value(edisp2(GEnergy(30.0,"TeV"), GEnergy(30.0,"TeV")), 1.74200329848547e-07);
1017 1017  
1018 1018 // Test if non-diagonal element (above diagonal) is zero
1019   - test_value(edisp2(std::log10(1.0), std::log10(30.0)), 0.0);
  1019 + test_value(edisp2(GEnergy(1.0,"TeV"), GEnergy(30.0,"TeV")), 0.0);
1020 1020  
1021 1021 // Test prob_erecobin method
1022 1022 GEnergy etrue(1.0, "TeV");
... ... @@ -1426,7 +1426,8 @@ void TestGCTAResponse::test_response_edisp_integration(const GCTAResponseIrf&amp; rs
1426 1426 for (double e_src = e_src_min; e_src < e_src_max; e_src *= 2.0) {
1427 1427  
1428 1428 // Compute log10 of true energy
1429   - double log10_e_src = std::log10(e_src);
  1429 + GEnergy etrue(e_src,"TeV");
  1430 + double log10_e_src = etrue.log10TeV();
1430 1431  
1431 1432 // Retrieve boundaries in observed energy
1432 1433 GEbounds ebounds = rsp.edisp()->ebounds_obs(log10_e_src);
... ... @@ -1441,17 +1442,16 @@ void TestGCTAResponse::test_response_edisp_integration(const GCTAResponseIrf&amp; rs
1441 1442 double E_obs = emin;
1442 1443 for (int i = 0; i < nE; ++i, E_obs += dE) {
1443 1444 GEnergy ereco(E_obs,"MeV");
1444   - double dp1 = (*rsp.edisp())(ereco.log10TeV(), log10_e_src);
1445   - double dp2 = rsp.edisp(ereco, 0.0, 0.0, 0.0, 0.0, log10_e_src);
  1445 + double dp1 = (*rsp.edisp())(ereco, etrue);
  1446 + double dp2 = rsp.edisp(ereco, etrue, 0.0, 0.0, 0.0, 0.0);
1446 1447 sum1 += dp1 * dE;
1447 1448 sum2 += dp2 * dE;
1448 1449 }
1449   - GEnergy eng(e_src, "TeV");
1450 1450 test_value(sum1, 1.0, 0.001, "Energy dispersion integration using "+
1451 1451 rsp.edisp()->classname()+"::operator() "
1452   - "for "+eng.print());
  1452 + "for "+etrue.print());
1453 1453 test_value(sum2, 1.0, 0.001, "Energy dispersion integration using "
1454   - "GCTAResponseIrf::edisp() for "+eng.print());
  1454 + "GCTAResponseIrf::edisp() for "+etrue.print());
1455 1455 }
1456 1456  
1457 1457 } // endif: energy dispersion was available
... ... @@ -1479,7 +1479,8 @@ void TestGCTAResponse::test_edisp_integration(const GCTAEdisp&amp; edisp,
1479 1479 for (double e_src = e_src_min; e_src <= e_src_max; e_src *= 2.0) {
1480 1480  
1481 1481 // Compute log10 of true energy
1482   - double log10_e_src = std::log10(e_src);
  1482 + GEnergy etrue(e_src,"TeV");
  1483 + double log10_e_src = etrue.log10TeV();
1483 1484  
1484 1485 // Retrieve boundaries in observed energy
1485 1486 GEbounds ebounds = edisp.ebounds_obs(log10_e_src);
... ... @@ -1498,13 +1499,12 @@ void TestGCTAResponse::test_edisp_integration(const GCTAEdisp&amp; edisp,
1498 1499 double E_obs = emin;
1499 1500 for (int i = 0; i < nE; ++i, E_obs += dE) {
1500 1501 GEnergy ereco(E_obs,"MeV");
1501   - double dp = edisp(ereco.log10TeV(), log10_e_src);
  1502 + double dp = edisp(ereco, etrue);
1502 1503 sum += dp * dE;
1503 1504 }
1504   - GEnergy eng(e_src, "TeV");
1505 1505 test_value(sum, 1.0, 0.001, "Energy dispersion integration using "+
1506 1506 edisp.classname()+"::operator() for "+
1507   - eng.print());
  1507 + etrue.print());
1508 1508 }
1509 1509  
1510 1510 // Return
... ...
inst/cta/test/test_CTA.py
  View file @6517c48
... ... @@ -237,25 +237,34 @@ class Test(gammalib.GPythonTestSuite):
237 237 """
238 238 Test GCTAEdisp classes
239 239 """
  240 + # Set some energies
  241 + eng0 = gammalib.GEnergy(1.0, 'TeV')
  242 + eng01 = gammalib.GEnergy(1.2589254, 'TeV')
  243 + eng001 = gammalib.GEnergy(1.023293, 'TeV')
  244 + eng0001 = gammalib.GEnergy(1.0023052, 'TeV')
  245 + eng1001 = gammalib.GEnergy(10.023052, 'TeV')
  246 + eng1 = gammalib.GEnergy(10.0, 'TeV')
  247 + eng30 = gammalib.GEnergy(30.0, 'TeV')
  248 +
240 249 # Test GCTAEdispRmf file constructor
241 250 filename = self._caldb + '/dc1/rmf.fits'
242 251 edisp = gammalib.GCTAEdispRmf(filename)
243 252  
244 253 # Test energy dispersion values
245   - self.test_value(edisp(math.log10(30),math.log10(1)), 0.0, 1.0e-9)
246   - self.test_value(edisp(math.log10(1),math.log10(30)), 0.0, 1.0e-9)
  254 + self.test_value(edisp(eng30, eng0), 0.0, 1.0e-9)
  255 + self.test_value(edisp(eng0, eng30), 0.0, 1.0e-9)
247 256  
248 257 # Test GCTAEdispPerfTable file constructor
249 258 filename = self._caldb + '/cta_dummy_irf.dat'
250 259 edisp = gammalib.GCTAEdispPerfTable(filename)
251 260  
252 261 # Test energy dispersion values
253   - self.test_value(edisp(0.0, 0.0), 4.3386871e-06, 1.0e-6)
254   - self.test_value(edisp(0.001, 0.0), 4.3273514e-06, 1.0e-6)
255   - self.test_value(edisp(0.01, 0.0), 4.1090489e-06, 1.0e-6)
256   - self.test_value(edisp(0.1, 0.0), 1.4984969e-07, 1.0e-6)
257   - self.test_value(edisp(1.0, 1.0), 7.8454726e-07, 1.0e-6)
258   - self.test_value(edisp(1.001, 1.0, 0.0), 7.8194077e-07, 1.0e-6)
  262 + self.test_value(edisp(eng0, eng0), 4.3386871e-06, 1.0e-6)
  263 + self.test_value(edisp(eng0001, eng0), 4.3273514e-06, 1.0e-6)
  264 + self.test_value(edisp(eng001, eng0), 4.1090489e-06, 1.0e-6)
  265 + self.test_value(edisp(eng01, eng0), 1.4984969e-07, 1.0e-6)
  266 + self.test_value(edisp(eng1, eng1), 7.8454726e-07, 1.0e-6)
  267 + self.test_value(edisp(eng1001, eng1, 0.0), 7.8194077e-07, 1.0e-6)
259 268  
260 269 # Test GCTAResponseIrf file constructor
261 270 db = gammalib.GCaldb(self._caldb)
... ...