GammaLib  1.7.0.dev
 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-2018 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 "GTools.hpp"
34 #include "GMath.hpp"
35 #include "GRan.hpp"
36 #include "GCTAEdispPerfTable.hpp"
37 #include "GCTAException.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 GCTAExceptionHandler::file_open_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  throw GCTAException::file_open_error(G_LOAD, filename.url());
279  }
280 
281  // Read lines
282  while (std::fgets(line, n, fptr) != NULL) {
283 
284  // Split line in elements. Strip empty elements from vector.
285  std::vector<std::string> elements = gammalib::split(line, " ");
286  for (int i = elements.size()-1; i >= 0; i--) {
287  if (gammalib::strip_whitespace(elements[i]).length() == 0) {
288  elements.erase(elements.begin()+i);
289  }
290  }
291 
292  // Skip header
293  if (elements[0].find("log(E)") != std::string::npos) {
294  continue;
295  }
296 
297  // Break loop if end of data table has been reached
298  if (elements[0].find("----------") != std::string::npos) {
299  break;
300  }
301 
302  // Push elements in node array and vector
303  m_logE.append(gammalib::todouble(elements[0]));
304 
305  // The energy resolution is stored as RMS(ln(Eest/Etrue))
306  m_sigma.push_back(gammalib::todouble(elements[4]) * gammalib::inv_ln10);
307 
308  } // endwhile: looped over lines
309 
310  // Close file
311  std::fclose(fptr);
312 
313  // Store filename
315 
316  // Return
317  return;
318 }
319 
320 
321 /***********************************************************************//**
322  * @brief Simulate energy dispersion
323  *
324  * @param[in] ran Random number generator.
325  * @param[in] etrue True photon energy.
326  * @param[in] theta Offset angle in camera system (radians). Not used.
327  * @param[in] phi Azimuth angle in camera system (radians). Not used.
328  * @param[in] zenith Zenith angle in Earth system (radians). Not used.
329  * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
330  * @return Energy.
331  *
332  * Draws observed energy value from a normal distribution of width
333  * m_par_sigma around @p logE.
334  ***************************************************************************/
336  const GEnergy& etrue,
337  const double& theta,
338  const double& phi,
339  const double& zenith,
340  const double& azimuth) const
341 {
342  // Get log10 TeV of true photon energy
343  double logEsrc = etrue.log10TeV();
344 
345  // Update the parameter cache
346  update(logEsrc);
347 
348  // Draw log observed energy in TeV
349  double logEobs = m_par_sigma * ran.normal() + logEsrc;
350 
351  // Set energy
352  GEnergy energy;
353  energy.log10TeV(logEobs);
354 
355  // Return energy
356  return energy;
357 }
358 
359 
360 /***********************************************************************//**
361  * @brief Return observed energy interval that contains the energy dispersion.
362  *
363  * @param[in] etrue True photon energy.
364  * @param[in] theta Offset angle in camera system (radians). Not used.
365  * @param[in] phi Azimuth angle in camera system (radians). Not used.
366  * @param[in] zenith Zenith angle in Earth system (radians). Not used.
367  * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
368  * @return Boundaries of reconstruced energies.
369  *
370  * Returns the band of observed energies outside of which the energy
371  * dispersion becomes negligible for a given true energy @p logEsrc. This
372  * band is set to \f$\pm 5 \times \sigma\f$, where \f$\sigma\f$ is the
373  * Gaussian width of the energy dispersion.
374  ***************************************************************************/
376  const double& theta,
377  const double& phi,
378  const double& zenith,
379  const double& azimuth) const
380 {
381  // Set energy band constant
382  const double number_of_sigmas = 5.0;
383 
384  // Get log10 of true photon energy in TeV
385  double etrue_log10TeV = etrue.log10TeV();
386 
387  // Get energy dispersion sigma
388  double sigma = m_logE.interpolate(etrue_log10TeV, m_sigma);
389 
390  // Compute energy boundaries
391  GEnergy emin;
392  GEnergy emax;
393  emin.log10TeV(etrue_log10TeV - number_of_sigmas * sigma);
394  emax.log10TeV(etrue_log10TeV + number_of_sigmas * sigma);
395 
396  // Return energy boundaries
397  return (GEbounds(emin, emax));
398 }
399 
400 
401 /***********************************************************************//**
402  * @brief Return true energy interval that contains the energy dispersion.
403  *
404  * @param[in] ereco Reconstructed event energy.
405  * @param[in] theta Offset angle in camera system (radians). Not used.
406  * @param[in] phi Azimuth angle in camera system (radians). Not used.
407  * @param[in] zenith Zenith angle in Earth system (radians). Not used.
408  * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
409  * @return Boundaries of true energies.
410  *
411  * Returns the band of true photon energies outside of which the energy
412  * dispersion becomes negligible for a given observed energy @p logEobs.
413  * This band is set to \f$\pm 5 \times \sigma\f$, where \f$\sigma\f$ is the
414  * Gaussian width of the energy dispersion.
415  ***************************************************************************/
417  const double& theta,
418  const double& phi,
419  const double& zenith,
420  const double& azimuth) const
421 {
422  // Set energy band constant
423  const double number_of_sigmas = 5.0;
424 
425  // Get log10 of reconstructed event energy in TeV
426  double ereco_log10TeV = ereco.log10TeV();
427 
428  // Get energy dispersion sigma
429  double sigma = m_logE.interpolate(ereco_log10TeV, m_sigma);
430 
431  // Compute energy boundaries
432  GEnergy emin;
433  GEnergy emax;
434  emin.log10TeV(ereco_log10TeV - number_of_sigmas * sigma);
435  emax.log10TeV(ereco_log10TeV + number_of_sigmas * sigma);
436 
437  // Return energy boundaries
438  return (GEbounds(emin, emax));
439 }
440 
441 
442 /***********************************************************************//**
443  * @brief Return energy dispersion probability for reconstructed energy
444  * interval
445  *
446  * @param[in] ereco_min Minimum of reconstructed energy interval.
447  * @param[in] ereco_max Maximum of reconstructed energy interval.
448  * @param[in] etrue True energy.
449  * @param[in] theta Offset angle (radians). Not used.
450  * @return Integrated energy dispersion probability.
451  *
452  * Computes
453  *
454  * \f[
455  * \int_{E_{\rm reco}^{\rm min}}^{E_{\rm reco}^{\rm max}}
456  * E_{\rm disp}(E_{\rm reco} | E_{\rm true}) \, dE_{\rm reco}
457  * \f]
458  *
459  * where
460  * \f$E_{\rm reco}\f$ is the reconstructed energy and
461  * \f$E_{\rm true}\f$ is the true energy.
462  ***************************************************************************/
464  const GEnergy& ereco_max,
465  const GEnergy& etrue,
466  const double& theta) const
467 {
468  // Update the parameter cache
469  update(etrue.log10TeV());
470 
471  // Get normalized energy boundaries in log10 energy
472  double xmin = (ereco_min.log10TeV() - etrue.log10TeV()) / m_par_sigma;
473  double xmax = (ereco_max.log10TeV() - etrue.log10TeV()) / m_par_sigma;
474 
475  // Compute fraction of probability within the energy boundaries
476  double prob = gammalib::gauss_integral(xmin, xmax);
477 
478  // Return
479  return prob;
480 }
481 
482 
483 /***********************************************************************//**
484  * @brief Print energy dispersion information
485  *
486  * @param[in] chatter Chattiness.
487  * @return String containing energy dispersion information.
488  ***************************************************************************/
489 std::string GCTAEdispPerfTable::print(const GChatter& chatter) const
490 {
491  // Initialise result string
492  std::string result;
493 
494  // Continue only if chatter is not silent
495  if (chatter != SILENT) {
496 
497  // Compute energy boundaries in TeV
498  int num = m_logE.size();
499  double emin = std::pow(10.0, m_logE[0]);
500  double emax = std::pow(10.0, m_logE[num-1]);
501 
502  // Append header
503  result.append("=== GCTAEdispPerfTable ===");
504 
505  // Append information
506  result.append("\n"+gammalib::parformat("Filename")+m_filename);
507  result.append("\n"+gammalib::parformat("Number of energy bins") +
508  gammalib::str(num));
509  result.append("\n"+gammalib::parformat("Log10(Energy) range"));
510  result.append(gammalib::str(emin)+" - "+gammalib::str(emax)+" TeV");
511 
512  // Append detailed information
513  GChatter reduced_chatter = gammalib::reduce(chatter);
514  if (reduced_chatter > SILENT) {
515  for(int i=0; i < num; ++i) {
516  double sigma = m_sigma[i];
517  double logE=m_logE[i];
518  result.append("\n"+gammalib::str(logE)+" "+gammalib::str(sigma));
519  }
520  }
521 
522  } // endif: chatter was not silent
523 
524  // Return result
525  return result;
526 }
527 
528 
529 /*==========================================================================
530  = =
531  = Private methods =
532  = =
533  ==========================================================================*/
534 
535 /***********************************************************************//**
536  * @brief Initialise class members
537  ***************************************************************************/
539 {
540  // Initialise members
541  m_filename.clear();
542  m_logE.clear();
543  m_sigma.clear();
544  m_par_logE = -1.0e30;
545  m_par_scale = 1.0;
546  m_par_sigma = 0.0;
547  m_par_width = 0.0;
548 
549  // Return
550  return;
551 }
552 
553 
554 /***********************************************************************//**
555  * @brief Copy class members
556  *
557  * @param[in] edisp Energy dispersion
558  ***************************************************************************/
560 {
561  // Copy members
562  m_filename = edisp.m_filename;
563  m_logE = edisp.m_logE;
564  m_sigma = edisp.m_sigma;
565  m_par_logE = edisp.m_par_logE;
566  m_par_scale = edisp.m_par_scale;
567  m_par_sigma = edisp.m_par_sigma;
568  m_par_width = edisp.m_par_width;
569 
570  // Return
571  return;
572 }
573 
574 
575 /***********************************************************************//**
576  * @brief Delete class members
577  ***************************************************************************/
579 {
580  // Return
581  return;
582 }
583 
584 
585 /***********************************************************************//**
586  * @brief Update energy dispersion parameter cache
587  *
588  * @param[in] logE Log10 of the true photon energy (\f$\log_{10}\f$ TeV).
589  *
590  * This method updates the energy dispersion parameter cache. As the
591  * performance table energy dispersion only depends on true photon energy,
592  * the only parameter on which the cache values depend is the true photon
593  * energy.
594  ***************************************************************************/
595 void GCTAEdispPerfTable::update(const double& logE) const
596 {
597  // Only compute energy dispersion parameters if arguments have changed
598  if (logE != m_par_logE) {
599 
600  // Save energy
601  m_par_logE = logE;
602 
603  // Determine Gaussian sigma and pre-compute Gaussian parameters
606  m_par_width = -0.5 / (m_par_sigma * m_par_sigma);
607 
608  }
609 
610  // Return
611  return;
612 }
void free_members(void)
Delete class members.
int size(void) const
Return number of nodes in node array.
Definition: GNodeArray.hpp:188
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:862
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:73
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
CTA exception handler interface definition.
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:1332
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1142
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:1022
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:55
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:805
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:413