GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAEdispPerfTable.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAEdispPerfTable.cpp - CTA performance table energy dispersion class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2014-2021 by Christoph Deil & Ellis Owen *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file GCTAEdispPerfTable.cpp
23  * @brief CTA performance table energy dispersion class implementation
24  * @author Christoph Deil & Ellis Owen
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <cstdio> // std::fopen, std::fgets, and std::fclose
32 #include <cmath>
33 #include "GException.hpp"
34 #include "GTools.hpp"
35 #include "GMath.hpp"
36 #include "GRan.hpp"
37 #include "GCTAEdispPerfTable.hpp"
38 
39 /* __ Method name definitions ____________________________________________ */
40 #define G_LOAD "GCTAEdispPerfTable::load(GFilename&)"
41 
42 /* __ Macros _____________________________________________________________ */
43 
44 /* __ Coding definitions _________________________________________________ */
45 
46 /* __ Debug definitions __________________________________________________ */
47 
48 /* __ Constants __________________________________________________________ */
49 
50 
51 /*==========================================================================
52  = =
53  = Constructors/destructors =
54  = =
55  ==========================================================================*/
56 
57 /***********************************************************************//**
58  * @brief Void constructor
59  ***************************************************************************/
61 {
62  // Initialise class members
63  init_members();
64 
65  // Return
66  return;
67 }
68 
69 
70 /***********************************************************************//**
71  * @brief File constructor
72  *
73  * @param[in] filename Performance table file name.
74  *
75  * Construct instance by loading the energy dispersion information from
76  * an ASCII performance table.
77  ***************************************************************************/
79  GCTAEdisp()
80 {
81  // Initialise class members
82  init_members();
83 
84  // Load energy dispersion from file
85  load(filename);
86 
87  // Return
88  return;
89 }
90 
91 
92 /***********************************************************************//**
93  * @brief Copy constructor
94  *
95  * @param[in] edisp Energy dispersion
96  ***************************************************************************/
98  GCTAEdisp(edisp)
99 {
100  // Initialise class members
101  init_members();
102 
103  // Copy members
104  copy_members(edisp);
105 
106  // Return
107  return;
108 }
109 
110 
111 /***********************************************************************//**
112  * @brief Destructor
113  ***************************************************************************/
115 {
116  // Free members
117  free_members();
118 
119  // Return
120  return;
121 }
122 
123 
124 /*==========================================================================
125  = =
126  = Operators =
127  = =
128  ==========================================================================*/
129 
130 /***********************************************************************//**
131  * @brief Assignment operator
132  *
133  * @param[in] edisp Energy dispersion
134  * @return Energy dispersion
135  ***************************************************************************/
137 {
138  // Execute only if object is not identical
139  if (this != &edisp) {
140 
141  // Copy base class members
142  this->GCTAEdisp::operator=(edisp);
143 
144  // Free members
145  free_members();
146 
147  // Initialise private members
148  init_members();
149 
150  // Copy members
151  copy_members(edisp);
152 
153  } // endif: object was not identical
154 
155  // Return this object
156  return *this;
157 }
158 
159 
160 /***********************************************************************//**
161  * @brief Return energy dispersion in units of MeV\f$^{-1}\f$
162  *
163  * @param[in] ereco Reconstructed photon energy.
164  * @param[in] etrue True photon energy.
165  * @param[in] theta Offset angle in camera system (radians). Not used.
166  * @param[in] phi Azimuth angle in camera system (radians). Not used.
167  * @param[in] zenith Zenith angle in Earth system (radians). Not used.
168  * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
169  * @return Energy dispersion (MeV\f$^{-1}\f$)
170  *
171  * Returns the energy dispersion
172  *
173  * \f[
174  * E_{\rm disp}(E_{\rm reco} | E_{\rm true}) =
175  * \frac{1}{\sqrt{2\pi}\sigma(E_{\rm true})}
176  * \exp \left(\frac{-(\log_{10} E_{\rm reco} - \log_{10} E_{\rm true})^2}
177  * {2 \sigma(E_{\rm true})^2} \right) \times
178  * \frac{1}{\log_{10} E_{\rm reco}}
179  * \f]
180  *
181  * in units of MeV\f$^{-1}\f$ where
182  * \f$E_{\rm reco}\f$ is the reconstructed energy,
183  * \f$E_{\rm true}\f$ is the true energy, and
184  * \f$\sigma(E_{\rm true})\f$ is the standard deviation of the energy
185  * dispersion that depends on the true photon energy.
186  ***************************************************************************/
188  const GEnergy& etrue,
189  const double& theta,
190  const double& phi,
191  const double& zenith,
192  const double& azimuth) const
193 {
194  // Get log10 of true and reconstructed photon energies
195  double logEsrc = etrue.log10TeV();
196  double logEobs = ereco.log10TeV();
197 
198  // Update the parameter cache
199  update(logEsrc);
200 
201  // Compute energy dispersion value
202  double delta = logEobs - logEsrc;
203  double edisp = m_par_scale * std::exp(m_par_width * delta * delta);
204 
205  // Compute energy dispersion per MeV
206  edisp /= (gammalib::ln10 * ereco.MeV());
207 
208  // Return energy dispersion
209  return edisp;
210 }
211 
212 
213 /*==========================================================================
214  = =
215  = Public methods =
216  = =
217  ==========================================================================*/
218 
219 /***********************************************************************//**
220  * @brief Clear instance
221  *
222  * This method properly resets the object to an initial state.
223  ***************************************************************************/
225 {
226  // Free class members (base and derived classes, derived class first)
227  free_members();
228  this->GCTAEdisp::free_members();
229 
230  // Initialise members
231  this->GCTAEdisp::init_members();
232  init_members();
233 
234  // Return
235  return;
236 }
237 
238 
239 /***********************************************************************//**
240  * @brief Clone instance
241  *
242  * @return Deep copy of instance.
243  ***************************************************************************/
245 {
246  return new GCTAEdispPerfTable(*this);
247 }
248 
249 
250 /***********************************************************************//**
251  * @brief Load energy dispersion from performance table
252  *
253  * @param[in] filename Performance table file name.
254  *
255  * @exception GException::file_error
256  * File could not be opened for read access.
257  *
258  * This method loads the energy dispersion information from an ASCII
259  * performance table. The energy resolution is stored in the 5th column
260  * of the performance table as RMS(ln(Eest/Etrue)). The method converts
261  * this internally to a sigma value by multiplying the stored values by
262  * 1/ln(10).
263  ***************************************************************************/
264 void GCTAEdispPerfTable::load(const GFilename& filename)
265 {
266 
267  // Clear arrays
268  m_logE.clear();
269  m_sigma.clear();
270 
271  // Allocate line buffer
272  const int n = 1000;
273  char line[n];
274 
275  // Open performance table readonly
276  FILE* fptr = std::fopen(filename.url().c_str(), "r");
277  if (fptr == NULL) {
278  std::string msg = "Energy dispersion file \""+filename.url()+
279  "\" not found or readable. Please specify a valid "
280  "and readable energy dispersion file.";
281  throw GException::file_error(G_LOAD, msg);
282  }
283 
284  // Read lines
285  while (std::fgets(line, n, fptr) != NULL) {
286 
287  // Split line in elements. Strip empty elements from vector.
288  std::vector<std::string> elements = gammalib::split(line, " ");
289  for (int i = elements.size()-1; i >= 0; i--) {
290  if (gammalib::strip_whitespace(elements[i]).length() == 0) {
291  elements.erase(elements.begin()+i);
292  }
293  }
294 
295  // Skip header
296  if (elements[0].find("log(E)") != std::string::npos) {
297  continue;
298  }
299 
300  // Break loop if end of data table has been reached
301  if (elements[0].find("----------") != std::string::npos) {
302  break;
303  }
304 
305  // Push elements in node array and vector
306  m_logE.append(gammalib::todouble(elements[0]));
307 
308  // The energy resolution is stored as RMS(ln(Eest/Etrue))
309  m_sigma.push_back(gammalib::todouble(elements[4]) * gammalib::inv_ln10);
310 
311  } // endwhile: looped over lines
312 
313  // Close file
314  std::fclose(fptr);
315 
316  // Store filename
318 
319  // Return
320  return;
321 }
322 
323 
324 /***********************************************************************//**
325  * @brief Simulate energy dispersion
326  *
327  * @param[in] ran Random number generator.
328  * @param[in] etrue True photon energy.
329  * @param[in] theta Offset angle in camera system (radians). Not used.
330  * @param[in] phi Azimuth angle in camera system (radians). Not used.
331  * @param[in] zenith Zenith angle in Earth system (radians). Not used.
332  * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
333  * @return Energy.
334  *
335  * Draws observed energy value from a normal distribution of width
336  * m_par_sigma around @p logE.
337  ***************************************************************************/
339  const GEnergy& etrue,
340  const double& theta,
341  const double& phi,
342  const double& zenith,
343  const double& azimuth) const
344 {
345  // Get log10 TeV of true photon energy
346  double logEsrc = etrue.log10TeV();
347 
348  // Update the parameter cache
349  update(logEsrc);
350 
351  // Draw log observed energy in TeV
352  double logEobs = m_par_sigma * ran.normal() + logEsrc;
353 
354  // Set energy
355  GEnergy energy;
356  energy.log10TeV(logEobs);
357 
358  // Return energy
359  return energy;
360 }
361 
362 
363 /***********************************************************************//**
364  * @brief Return observed energy interval that contains the energy dispersion.
365  *
366  * @param[in] etrue True photon energy.
367  * @param[in] theta Offset angle in camera system (radians). Not used.
368  * @param[in] phi Azimuth angle in camera system (radians). Not used.
369  * @param[in] zenith Zenith angle in Earth system (radians). Not used.
370  * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
371  * @return Boundaries of reconstruced energies.
372  *
373  * Returns the band of observed energies outside of which the energy
374  * dispersion becomes negligible for a given true energy @p logEsrc. This
375  * band is set to \f$\pm 5 \times \sigma\f$, where \f$\sigma\f$ is the
376  * Gaussian width of the energy dispersion.
377  ***************************************************************************/
379  const double& theta,
380  const double& phi,
381  const double& zenith,
382  const double& azimuth) const
383 {
384  // Set energy band constant
385  const double number_of_sigmas = 5.0;
386 
387  // Get log10 of true photon energy in TeV
388  double etrue_log10TeV = etrue.log10TeV();
389 
390  // Get energy dispersion sigma
391  double sigma = m_logE.interpolate(etrue_log10TeV, m_sigma);
392 
393  // Compute energy boundaries
394  GEnergy emin;
395  GEnergy emax;
396  emin.log10TeV(etrue_log10TeV - number_of_sigmas * sigma);
397  emax.log10TeV(etrue_log10TeV + number_of_sigmas * sigma);
398 
399  // Return energy boundaries
400  return (GEbounds(emin, emax));
401 }
402 
403 
404 /***********************************************************************//**
405  * @brief Return true energy interval that contains the energy dispersion.
406  *
407  * @param[in] ereco Reconstructed event energy.
408  * @param[in] theta Offset angle in camera system (radians). Not used.
409  * @param[in] phi Azimuth angle in camera system (radians). Not used.
410  * @param[in] zenith Zenith angle in Earth system (radians). Not used.
411  * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
412  * @return Boundaries of true energies.
413  *
414  * Returns the band of true photon energies outside of which the energy
415  * dispersion becomes negligible for a given observed energy @p logEobs.
416  * This band is set to \f$\pm 5 \times \sigma\f$, where \f$\sigma\f$ is the
417  * Gaussian width of the energy dispersion.
418  ***************************************************************************/
420  const double& theta,
421  const double& phi,
422  const double& zenith,
423  const double& azimuth) const
424 {
425  // Set energy band constant
426  const double number_of_sigmas = 5.0;
427 
428  // Get log10 of reconstructed event energy in TeV
429  double ereco_log10TeV = ereco.log10TeV();
430 
431  // Get energy dispersion sigma
432  double sigma = m_logE.interpolate(ereco_log10TeV, m_sigma);
433 
434  // Compute energy boundaries
435  GEnergy emin;
436  GEnergy emax;
437  emin.log10TeV(ereco_log10TeV - number_of_sigmas * sigma);
438  emax.log10TeV(ereco_log10TeV + number_of_sigmas * sigma);
439 
440  // Return energy boundaries
441  return (GEbounds(emin, emax));
442 }
443 
444 
445 /***********************************************************************//**
446  * @brief Return energy dispersion probability for reconstructed energy
447  * interval
448  *
449  * @param[in] ereco_min Minimum of reconstructed energy interval.
450  * @param[in] ereco_max Maximum of reconstructed energy interval.
451  * @param[in] etrue True energy.
452  * @param[in] theta Offset angle (radians). Not used.
453  * @return Integrated energy dispersion probability.
454  *
455  * Computes
456  *
457  * \f[
458  * \int_{E_{\rm reco}^{\rm min}}^{E_{\rm reco}^{\rm max}}
459  * E_{\rm disp}(E_{\rm reco} | E_{\rm true}) \, dE_{\rm reco}
460  * \f]
461  *
462  * where
463  * \f$E_{\rm reco}\f$ is the reconstructed energy and
464  * \f$E_{\rm true}\f$ is the true energy.
465  ***************************************************************************/
467  const GEnergy& ereco_max,
468  const GEnergy& etrue,
469  const double& theta) const
470 {
471  // Update the parameter cache
472  update(etrue.log10TeV());
473 
474  // Get normalized energy boundaries in log10 energy
475  double xmin = (ereco_min.log10TeV() - etrue.log10TeV()) / m_par_sigma;
476  double xmax = (ereco_max.log10TeV() - etrue.log10TeV()) / m_par_sigma;
477 
478  // Compute fraction of probability within the energy boundaries
479  double prob = gammalib::gauss_integral(xmin, xmax);
480 
481  // Return
482  return prob;
483 }
484 
485 
486 /***********************************************************************//**
487  * @brief Print energy dispersion information
488  *
489  * @param[in] chatter Chattiness.
490  * @return String containing energy dispersion information.
491  ***************************************************************************/
492 std::string GCTAEdispPerfTable::print(const GChatter& chatter) const
493 {
494  // Initialise result string
495  std::string result;
496 
497  // Continue only if chatter is not silent
498  if (chatter != SILENT) {
499 
500  // Compute energy boundaries in TeV
501  int num = m_logE.size();
502  double emin = std::pow(10.0, m_logE[0]);
503  double emax = std::pow(10.0, m_logE[num-1]);
504 
505  // Append header
506  result.append("=== GCTAEdispPerfTable ===");
507 
508  // Append information
509  result.append("\n"+gammalib::parformat("Filename")+m_filename);
510  result.append("\n"+gammalib::parformat("Number of energy bins") +
511  gammalib::str(num));
512  result.append("\n"+gammalib::parformat("Energy range"));
513  result.append(gammalib::str(emin)+" - "+gammalib::str(emax)+" TeV");
514 
515  // Append detailed information
516  GChatter reduced_chatter = gammalib::reduce(chatter);
517  if (reduced_chatter > SILENT) {
518  for(int i=0; i < num; ++i) {
519  double sigma = m_sigma[i];
520  double logE=m_logE[i];
521  result.append("\n"+gammalib::str(logE)+" "+gammalib::str(sigma));
522  }
523  }
524 
525  } // endif: chatter was not silent
526 
527  // Return result
528  return result;
529 }
530 
531 
532 /*==========================================================================
533  = =
534  = Private methods =
535  = =
536  ==========================================================================*/
537 
538 /***********************************************************************//**
539  * @brief Initialise class members
540  ***************************************************************************/
542 {
543  // Initialise members
544  m_filename.clear();
545  m_logE.clear();
546  m_sigma.clear();
547  m_par_logE = -1.0e30;
548  m_par_scale = 1.0;
549  m_par_sigma = 0.0;
550  m_par_width = 0.0;
551 
552  // Return
553  return;
554 }
555 
556 
557 /***********************************************************************//**
558  * @brief Copy class members
559  *
560  * @param[in] edisp Energy dispersion
561  ***************************************************************************/
563 {
564  // Copy members
565  m_filename = edisp.m_filename;
566  m_logE = edisp.m_logE;
567  m_sigma = edisp.m_sigma;
568  m_par_logE = edisp.m_par_logE;
569  m_par_scale = edisp.m_par_scale;
570  m_par_sigma = edisp.m_par_sigma;
571  m_par_width = edisp.m_par_width;
572 
573  // Return
574  return;
575 }
576 
577 
578 /***********************************************************************//**
579  * @brief Delete class members
580  ***************************************************************************/
582 {
583  // Return
584  return;
585 }
586 
587 
588 /***********************************************************************//**
589  * @brief Update energy dispersion parameter cache
590  *
591  * @param[in] logE Log10 of the true photon energy (\f$\log_{10}\f$ TeV).
592  *
593  * This method updates the energy dispersion parameter cache. As the
594  * performance table energy dispersion only depends on true photon energy,
595  * the only parameter on which the cache values depend is the true photon
596  * energy.
597  ***************************************************************************/
598 void GCTAEdispPerfTable::update(const double& logE) const
599 {
600  // Only compute energy dispersion parameters if arguments have changed
601  if (logE != m_par_logE) {
602 
603  // Save energy
604  m_par_logE = logE;
605 
606  // Determine Gaussian sigma and pre-compute Gaussian parameters
609  m_par_width = -0.5 / (m_par_sigma * m_par_sigma);
610 
611  }
612 
613  // Return
614  return;
615 }
void free_members(void)
Delete class members.
int size(void) const
Return number of nodes in node array.
Definition: GNodeArray.hpp:192
std::string print(const GChatter &chatter=NORMAL) const
Print energy dispersion information.
Random number generator class definition.
std::vector< double > m_sigma
Sigma value (rms) of energy resolution.
double gauss_integral(const double &x1, const double &x2)
Returns the integral of a Gaussian function.
Definition: GMath.cpp:605
void free_members(void)
Delete class members.
Definition: GCTAEdisp.cpp:164
Abstract base class for the CTA energy dispersion.
Definition: GCTAEdisp.hpp:49
double m_par_width
Gaussian width parameter.
void clear(void)
Clear node array.
Definition: GNodeArray.cpp:235
double log10TeV(void) const
Return log10 of energy in TeV.
Definition: GEnergy.cpp:457
GEbounds ereco_bounds(const GEnergy &etrue, const double &theta=0.0, const double &phi=0.0, const double &zenith=0.0, const double &azimuth=0.0) const
Return observed energy interval that contains the energy dispersion.
Random number generator class.
Definition: GRan.hpp:44
std::vector< std::string > split(const std::string &s, const std::string &sep)
Split string.
Definition: GTools.cpp:983
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
Gammalib tools definition.
GCTAEdispPerfTable & operator=(const GCTAEdispPerfTable &psf)
Assignment operator.
std::string strip_whitespace(const std::string &arg)
Strip leading and trailing whitespace from string.
Definition: GTools.cpp:80
GCTAEdispPerfTable * clone(void) const
Clone instance.
const double ln10
Definition: GMath.hpp:46
GEnergy mc(GRan &ran, const GEnergy &etrue, const double &theta=0.0, const double &phi=0.0, const double &zenith=0.0, const double &azimuth=0.0) const
Simulate energy dispersion.
void copy_members(const GCTAEdispPerfTable &psf)
Copy class members.
void clear(void)
Clear instance.
double m_par_logE
Energy for which precomputation is done.
Energy boundaries container class.
Definition: GEbounds.hpp:60
GFilename filename(void) const
Return filename.
CTA performance table energy dispersion class.
virtual ~GCTAEdispPerfTable(void)
Destructor.
Filename class.
Definition: GFilename.hpp:62
CTA performance table energy dispersion class definition.
double m_par_sigma
Gaussian sigma.
double operator()(const GEnergy &ereco, const GEnergy &etrue, const double &theta=0.0, const double &phi=0.0, const double &zenith=0.0, const double &azimuth=0.0) const
Return energy dispersion in units of MeV .
double interpolate(const double &value, const std::vector< double > &vector) const
Interpolate value.
Definition: GNodeArray.cpp:536
GChatter
Definition: GTypemaps.hpp:33
double normal(void)
Returns normal deviates.
Definition: GRan.cpp:257
double m_par_scale
Gaussian normalization.
GCTAEdisp & operator=(const GCTAEdisp &edisp)
Assignment operator.
Definition: GCTAEdisp.cpp:106
void load(const GFilename &filename)
Load energy dispersion from performance table.
GCTAEdispPerfTable(void)
Void constructor.
std::string url(void) const
Return Uniform Resource Locator (URL)
Definition: GFilename.hpp:189
void clear(void)
Clear file name.
Definition: GFilename.cpp:188
#define G_LOAD
GNodeArray m_logE
log(E) nodes for interpolation
GChatter reduce(const GChatter &chatter)
Reduce chattiness by one level.
Definition: GTypemaps.hpp:65
const double inv_sqrt2pi
Definition: GMath.hpp:41
GVector pow(const GVector &vector, const double &power)
Computes tanh of vector elements.
Definition: GVector.cpp:1422
Exception handler interface definition.
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1232
GEbounds etrue_bounds(const GEnergy &ereco, const double &theta=0.0, const double &phi=0.0, const double &zenith=0.0, const double &azimuth=0.0) const
Return true energy interval that contains the energy dispersion.
void init_members(void)
Initialise class members.
Definition: GCTAEdisp.cpp:142
double prob_erecobin(const GEnergy &ereco_min, const GEnergy &ereco_max, const GEnergy &etrue, const double &theta) const
Return energy dispersion probability for reconstructed energy interval.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1143
void append(const double &node)
Append one node to array.
Definition: GNodeArray.cpp:351
GFilename m_filename
Name of response file.
const double inv_ln10
Definition: GMath.hpp:48
void update(const double &logE) const
Update energy dispersion parameter cache.
void init_members(void)
Initialise class members.
Mathematical function definitions.
double todouble(const std::string &arg)
Convert string into double precision value.
Definition: GTools.cpp:926
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489