GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAEdispRmf.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAEdispRmf.cpp - CTA RMF 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 GCTAEdispRmf.cpp
23  * @brief CTA RMF 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 <cmath>
32 #include "GTools.hpp"
33 #include "GMath.hpp"
34 #include "GFilename.hpp"
35 #include "GRan.hpp"
36 #include "GIntegral.hpp"
37 #include "GFunction.hpp"
38 #include "GEnergy.hpp"
39 #include "GRmf.hpp"
40 #include "GNodeArray.hpp"
41 #include "GMatrixSparse.hpp"
42 #include "GCTAEdispRmf.hpp"
43 
44 /* __ Method name definitions ____________________________________________ */
45 #define G_LOAD "GCTAEdispRmf::load(GFilename&)"
46 #define G_MC "GCTAEdispRmf::mc(GRan&, GEnergy&, double&, double&, double&, "\
47  "double&)"
48 
49 /* __ Macros _____________________________________________________________ */
50 
51 /* __ Coding definitions _________________________________________________ */
52 //#define G_LINEAR_ERECO_INTERPOLATION //!< Linear interpolation for Ereco
53 
54 /* __ Debug definitions __________________________________________________ */
55 
56 /* __ Constants __________________________________________________________ */
57 
58 
59 /*==========================================================================
60  = =
61  = Constructors/destructors =
62  = =
63  ==========================================================================*/
64 
65 /***********************************************************************//**
66  * @brief Void constructor
67  ***************************************************************************/
69 {
70  // Initialise class members
71  init_members();
72 
73  // Return
74  return;
75 }
76 
77 
78 /***********************************************************************//**
79  * @brief File constructor
80  *
81  * @param[in] filename Redistribution Matrix File name.
82  *
83  * Construct instance by loading the Redistribution Matrix File.
84  ***************************************************************************/
86 {
87  // Initialise class members
88  init_members();
89 
90  // Load energy dispersion from file
91  load(filename);
92 
93  // Return
94  return;
95 }
96 
97 
98 /***********************************************************************//**
99  * @brief Copy constructor
100  *
101  * @param[in] edisp Energy dispersion
102  ***************************************************************************/
104 {
105  // Initialise class members
106  init_members();
107 
108  // Copy members
109  copy_members(edisp);
110 
111  // Return
112  return;
113 }
114 
115 
116 /***********************************************************************//**
117  * @brief Destructor
118  ***************************************************************************/
120 {
121  // Free members
122  free_members();
123 
124  // Return
125  return;
126 }
127 
128 
129 /*==========================================================================
130  = =
131  = Operators =
132  = =
133  ==========================================================================*/
134 
135 /***********************************************************************//**
136  * @brief Assignment operator
137  *
138  * @param[in] edisp Energy dispersion
139  * @return Energy dispersion
140  ***************************************************************************/
142 {
143  // Execute only if object is not identical
144  if (this != &edisp) {
145 
146  // Copy base class members
147  this->GCTAEdisp::operator=(edisp);
148 
149  // Free members
150  free_members();
151 
152  // Initialise private members
153  init_members();
154 
155  // Copy members
156  copy_members(edisp);
157 
158  } // endif: object was not identical
159 
160  // Return this object
161  return *this;
162 }
163 
164 
165 /***********************************************************************//**
166  * @brief Return energy dispersion in units of MeV\f$^{-1}\f$
167  *
168  * @param[in] ereco Reconstructed photon energy.
169  * @param[in] etrue True photon energy.
170  * @param[in] theta Offset angle in camera system (radians). Not used.
171  * @param[in] phi Azimuth angle in camera system (radians). Not used.
172  * @param[in] zenith Zenith angle in Earth system (radians). Not used.
173  * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
174  * @return Energy dispersion (MeV\f$^{-1}\f$)
175  *
176  * Returns the energy dispersion
177  *
178  * \f[
179  * E_{\rm disp}(E_{\rm reco} | E_{\rm true})
180  * \f]
181  *
182  * in units of MeV\f$^{-1}\f$ where
183  * \f$E_{\rm reco}\f$ is the reconstructed energy, and
184  * \f$E_{\rm true}\f$ is the true energy.
185  ***************************************************************************/
186 double GCTAEdispRmf::operator()(const GEnergy& ereco,
187  const GEnergy& etrue,
188  const double& theta,
189  const double& phi,
190  const double& zenith,
191  const double& azimuth) const
192 {
193  // Update indexes and weighting factors for interpolation
194  update(ereco, etrue);
195 
196  // Perform interpolation
197  double edisp = m_wgt1 * m_matrix(m_itrue1, m_ireco1) +
201 
202  // Make sure that energy dispersion is not negative
203  if (edisp < 0.0) {
204  edisp = 0.0;
205  }
206 
207  // Return energy dispersion
208  return edisp;
209 }
210 
211 
212 /*==========================================================================
213  = =
214  = Public methods =
215  = =
216  ==========================================================================*/
217 
218 /***********************************************************************//**
219  * @brief Clear instance
220  *
221  * This method properly resets the object to an initial state.
222  ***************************************************************************/
224 {
225  // Free class members (base and derived classes, derived class first)
226  free_members();
227  this->GCTAEdisp::free_members();
228 
229  // Initialise members
230  this->GCTAEdisp::init_members();
231  init_members();
232 
233  // Return
234  return;
235 }
236 
237 
238 /***********************************************************************//**
239  * @brief Clone instance
240  *
241  * @return Deep copy of instance.
242  ***************************************************************************/
244 {
245  return new GCTAEdispRmf(*this);
246 }
247 
248 
249 /***********************************************************************//**
250  * @brief Load energy dispersion from RMF file
251  *
252  * @param[in] filename of RMF file.
253  *
254  * This method loads the energy dispersion information from an RMF file.
255  ***************************************************************************/
256 void GCTAEdispRmf::load(const GFilename& filename)
257 {
258  // Load RMF file
259  m_rmf.load(filename);
260 
261  // Store the filename
263 
264  // Set cache (has to come before set_matrix())
265  set_cache();
266 
267  // Set matrix
268  set_matrix();
269 
270  // Set maximum edisp value
271  set_max_edisp();
272 
273  // Return
274  return;
275 }
276 
277 
278 /***********************************************************************//**
279  * @brief Simulate energy dispersion
280  *
281  * @param[in] ran Random number generator.
282  * @param[in] etrue True photon energy.
283  * @param[in] theta Offset angle in camera system (radians). Not used.
284  * @param[in] phi Azimuth angle in camera system (radians). Not used.
285  * @param[in] zenith Zenith angle in Earth system (radians). Not used.
286  * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
287  * @return Reconstructed energy.
288  *
289  * @exception GException::invalid_return_value
290  * Energy dispersion matrix is empty.
291  *
292  * Draws reconstructed energy value from RMF matrix given a true energy
293  * @p etrue. If no energy dispersion information is available the method
294  * will return the true photon energy.
295  ***************************************************************************/
297  const GEnergy& etrue,
298  const double& theta,
299  const double& phi,
300  const double& zenith,
301  const double& azimuth) const
302 {
303  // Initialise reconstructed event energy with true photon energy
304  GEnergy ereco = etrue;
305 
306  // Get boundaries for reconstructed energy
307  GEbounds ebounds = ereco_bounds(etrue, theta, phi, zenith, azimuth);
308  double emin = ebounds.emin().TeV();
309  double emax = ebounds.emax().TeV();
310 
311  // Get maximum energy dispersion value (including some margin)
312  double max_edisp = 1.5 * m_max_edisp;
313 
314  // Throw an exception if maximum energy dispersion is zero
315  if (max_edisp <= 0.0) {
316  std::string msg = "Energy dispersion matrix is empty. Please provide "
317  "a valid energy dispersion matrix.";
319  }
320 
321  // Find energy by rejection method
322  double ewidth = emax - emin;
323  if (max_edisp > 0.0) {
324  double f = 0.0;
325  double ftest = 1.0;
326  while (ftest > f) {
327  ereco.TeV(emin + ewidth * ran.uniform());
328  f = operator()(ereco, etrue, theta, phi, zenith, azimuth);
329  ftest = ran.uniform() * max_edisp;
330  }
331  }
332 
333  // Return reconstructed photon energy
334  return ereco;
335 }
336 
337 
338 /***********************************************************************//**
339  * @brief Return observed energy interval that contains the energy dispersion.
340  *
341  * @param[in] etrue True photon energy.
342  * @param[in] theta Offset angle in camera system (radians). Not used.
343  * @param[in] phi Azimuth angle in camera system (radians). Not used.
344  * @param[in] zenith Zenith angle in Earth system (radians). Not used.
345  * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
346  * @return Reconstructed energy boundaries.
347  *
348  * Returns the band of observed energies outside of which the energy
349  * dispersion becomes negligible for a given true energy @p etrue.
350  ***************************************************************************/
352  const double& theta,
353  const double& phi,
354  const double& zenith,
355  const double& azimuth) const
356 {
357  // Compute only if parameters changed
359 
360  // Store last theta value
362  m_last_etrue_bounds.TeV(0.0); // force update
363 
364  // Compute m_ereco_bounds
366 
367  }
368 
369  // Search index only if logEsrc has changed
370  if (etrue != m_last_etrue_bounds) {
371 
372  // Store last photon energy
373  m_last_etrue_bounds = etrue;
374 
375  // Get true photon energy in TeV
376  double etrue_TeV = etrue.TeV();
377 
378  // Find right index with bisection
379  int low = 0;
380  int high = m_ereco_bounds.size() - 1;
381  while ((high-low) > 1) {
382  int mid = (low+high) / 2;
383  double e = m_rmf.etrue().emin(mid).TeV();
384  if (etrue_TeV < e) {
385  high = mid;
386  }
387  else {
388  low = mid;
389  }
390  }
391 
392  // Index found
393  m_index_ereco = low;
394 
395  }
396 
397  // Return energy boundaries
398  return (m_ereco_bounds[m_index_ereco]);
399 }
400 
401 
402 /***********************************************************************//**
403  * @brief Return true energy interval that contains the energy dispersion.
404  *
405  * @param[in] ereco Reconstructed event energy.
406  * @param[in] theta Offset angle in camera system (radians). Not used.
407  * @param[in] phi Azimuth angle in camera system (radians). Not used.
408  * @param[in] zenith Zenith angle in Earth system (radians). Not used.
409  * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
410  * @return True energy boundaries.
411  *
412  * Returns the band of true photon energies outside of which the energy
413  * dispersion becomes negligible for a given observed energy @p logEobs.
414  ***************************************************************************/
416  const double& theta,
417  const double& phi,
418  const double& zenith,
419  const double& azimuth) const
420 {
421  // Compute only if parameters changed
423 
424  // Set computation flag
426  m_last_ereco_bounds.TeV(0.0); // force update
427 
428  // Compute m_etrue_bounds
430 
431  }
432 
433  // Search index only if ereco has changed
434  if (ereco != m_last_ereco_bounds) {
435 
436  // Store last value
437  m_last_ereco_bounds = ereco;
438 
439  // Get reconstructed energy in TeV
440  double ereco_TeV = ereco.TeV();
441 
442  // Find right index with bisection
443  int low = 0;
444  int high = m_etrue_bounds.size() - 1;
445  while ((high-low) > 1) {
446  int mid = (low+high) / 2;
447  double e = m_rmf.emeasured().emin(mid).TeV();
448  if (ereco_TeV < e) {
449  high = mid;
450  }
451  else {
452  low = mid;
453  }
454  }
455 
456  // Index found
457  m_index_etrue = low;
458 
459  }
460 
461  // Return energy boundaries
462  return (m_etrue_bounds[m_index_etrue]);
463 }
464 
465 
466 /***********************************************************************//**
467  * @brief Return energy dispersion probability for reconstructed energy
468  * interval
469  *
470  * @param[in] ereco_min Minimum of reconstructed energy interval.
471  * @param[in] ereco_max Maximum of reconstructed energy interval.
472  * @param[in] etrue True energy.
473  * @param[in] theta Offset angle. Not used.
474  * @return Integrated energy dispersion probability.
475  *
476  * Computes
477  *
478  * \f[
479  * \int_{E_{\rm reco}^{\rm min}}^{E_{\rm reco}^{\rm max}}
480  * E_{\rm disp}(E_{\rm reco} | E_{\rm true}) \, dE_{\rm reco}
481  * \f]
482  *
483  * where
484  * \f$E_{\rm reco}\f$ is the reconstructed energy and
485  * \f$E_{\rm true}\f$ is the true energy.
486  ***************************************************************************/
487 double GCTAEdispRmf::prob_erecobin(const GEnergy& ereco_min,
488  const GEnergy& ereco_max,
489  const GEnergy& etrue,
490  const double& theta) const
491 {
492  // Initalize probability
493  double prob = 0.0;
494 
495  // Get log10 of energies
496  double logEsrc = etrue.log10TeV();
497  #if defined(G_LINEAR_ERECO_INTERPOLATION)
498  double logEobs_min = ereco_min.TeV();
499  double logEobs_max = ereco_max.TeV();
500  #else
501  double logEobs_min = ereco_min.log10TeV();
502  double logEobs_max = ereco_max.log10TeV();
503  #endif
504 
505  // Set logEsrc for node array interpolation
506  m_etrue.set_value(logEsrc);
507 
508  // Store indices and weights for interpolation
509  int inx_left = m_etrue.inx_left();
510  int inx_right = m_etrue.inx_right();
511  double wgt_left = m_etrue.wgt_left();
512  double wgt_right = m_etrue.wgt_right();
513 
514  // Initialise first and second node
515  GEnergy x1 = ereco_min;
516  double f1 = this->operator()(x1, etrue);
517  GEnergy x2;
518  double f2 = 0.0;
519 
520  // Loop over all measured energy nodes
521  for (int i = 0; i < m_ereco.size(); ++i) {
522 
523  // If measured energy is below logEobs_min then skip node
524  if (m_ereco[i] <= logEobs_min) {
525  continue;
526  }
527 
528  // If measured energy is above maximum migration value then use
529  // logEobs_max as measured energy
530  if (m_ereco[i] > logEobs_max) {
531  x2 = ereco_max;
532  f2 = this->operator()(x2, etrue);
533  }
534 
535  // ... otherwise use measured energy
536  else {
537  #if defined(G_LINEAR_ERECO_INTERPOLATION)
538  x2.TeV(m_ereco[i]);
539  #else
540  x2.log10TeV(m_ereco[i]);
541  #endif
542  f2 = wgt_left * m_matrix(inx_left, i) +
543  wgt_right * m_matrix(inx_right, i);
544  if (f2 < 0.0) {
545  f2 = 0.0;
546  }
547  }
548 
549  // Compute integral
550  prob += 0.5 * (f1 + f2) * (x2.MeV() - x1.MeV());
551 
552  // Set second node as first node
553  x1 = x2;
554  f1 = f2;
555 
556  // If node energy is above migra_max then break now
557  if (m_ereco[i] > logEobs_max) {
558  break;
559  }
560 
561  } // endfor: looped over all nodes
562 
563  // If last node energy is below migra_max then compute last part of
564  // integral up to emax
565  if (x1 < ereco_max) {
566  x2 = ereco_max;
567  f2 = this->operator()(x2, etrue);
568  prob += 0.5 * (f1 + f2) * (x2.MeV() - x1.MeV());
569  }
570 
571  // Return probability
572  return prob;
573 }
574 
575 
576 /***********************************************************************//**
577  * @brief Print RMF information
578  *
579  * @param[in] chatter Chattiness.
580  * @return String containing energy dispersion information.
581  ***************************************************************************/
582 std::string GCTAEdispRmf::print(const GChatter& chatter) const
583 {
584  // Initialise result string
585  std::string result;
586 
587  // Continue only if chatter is not silent
588  if (chatter != SILENT) {
589 
590  // Append header
591  result.append("=== GCTAEdispRmf ===");
592 
593  // Append energy boundary information
594  result.append("\n"+gammalib::parformat("Filename")+m_filename);
595  result.append("\n"+gammalib::parformat("Number of true energy bins"));
596  result.append(gammalib::str(m_rmf.etrue().size()));
597  result.append("\n"+gammalib::parformat("Number of measured bins"));
598  result.append(gammalib::str(m_rmf.emeasured().size()));
599  result.append("\n"+gammalib::parformat("True energy range"));
600  result.append(m_rmf.etrue().emin().print());
601  result.append(" - ");
602  result.append(m_rmf.etrue().emax().print());
603  result.append("\n"+gammalib::parformat("Measured energy range"));
604  result.append(m_rmf.emeasured().emin().print());
605  result.append(" - ");
606  result.append(m_rmf.emeasured().emax().print());
607 
608  } // endif: chatter was not silent
609 
610  // Return result
611  return result;
612 
613 }
614 
615 
616 /*==========================================================================
617  = =
618  = Private methods =
619  = =
620  ==========================================================================*/
621 
622 /***********************************************************************//**
623  * @brief Initialise class members
624  ***************************************************************************/
626 {
627  // Initialise members
628  m_filename.clear();
629  m_rmf.clear();
630  m_matrix.clear();
631  m_max_edisp = 0.0;
632 
633  // Initialise interpolation cache
634  m_etrue.clear();
635  m_ereco.clear();
638  m_itrue1 = 0;
639  m_itrue2 = 0;
640  m_ireco1 = 0;
641  m_ireco2 = 0;
642  m_wgt1 = 0.0;
643  m_wgt2 = 0.0;
644  m_wgt3 = 0.0;
645  m_wgt4 = 0.0;
646 
647  // Initialise computation cache
648  m_ereco_bounds_computed = false;
649  m_etrue_bounds_computed = false;
650  m_index_ereco = 0;
651  m_index_etrue = 0;
654  m_ereco_bounds.clear();
655  m_etrue_bounds.clear();
656 
657  // Return
658  return;
659 }
660 
661 
662 /***********************************************************************//**
663  * @brief Copy class members
664  *
665  * @param[in] edisp Energy dispersion
666  ***************************************************************************/
668 {
669  // Copy members
670  m_filename = edisp.m_filename;
671  m_rmf = edisp.m_rmf;
672  m_matrix = edisp.m_matrix;
673  m_max_edisp = edisp.m_max_edisp;
674 
675  // Copy interpolation cache
676  m_etrue = edisp.m_etrue;
677  m_ereco = edisp.m_ereco;
678  m_last_etrue = edisp.m_last_etrue;
679  m_last_ereco = edisp.m_last_ereco;
680  m_itrue1 = edisp.m_itrue1;
681  m_itrue2 = edisp.m_itrue2;
682  m_ireco1 = edisp.m_ireco1;
683  m_ireco2 = edisp.m_ireco2;
684  m_wgt1 = edisp.m_wgt1;
685  m_wgt2 = edisp.m_wgt2;
686  m_wgt3 = edisp.m_wgt3;
687  m_wgt4 = edisp.m_wgt4;
688 
689  // Copy computation cache
698 
699  // Return
700  return;
701 }
702 
703 
704 /***********************************************************************//**
705  * @brief Delete class members
706  ***************************************************************************/
708 {
709  // Return
710  return;
711 }
712 
713 
714 /***********************************************************************//**
715  * @brief Set redistribution matrix
716  ***************************************************************************/
718 {
719  // Determine matrix rows and columns
720  int rows = m_rmf.ntrue();
721  int columns = m_rmf.nmeasured();
722 
723  // Initialize matrix
724  m_matrix = GMatrixSparse(rows, columns);
725 
726  // Get reconstructued energy bin boundaries
727  GEbounds ebounds = m_rmf.emeasured();
728 
729  // Fill matrix elements
730  for (int ireco = 0; ireco < columns; ++ireco) {
731 
732  // Compute reconstructued energy bin width in MeV
733  double ewidth = ebounds.emax(ireco).MeV() - ebounds.emin(ireco).MeV();
734 
735  // If bin width is positive then fill all matrix bins by deviding
736  // the RMF values by the bin width. This assures that the stored
737  // values are per MeV of reconstructed energy.
738  if (ewidth > 0.0) {
739  for (int itrue = 0; itrue < rows; ++itrue) {
740  m_matrix(itrue, ireco) = m_rmf(itrue, ireco) / ewidth;
741  }
742  }
743 
744  } // endfor: looped over all reconstructed energies
745 
746  // Initialise row sums
747  std::vector<double> row_sums(rows, 0.0);
748 
749  // Normalize matrix elements
750  for (int itrue = 0; itrue < rows; ++itrue) {
751 
752  // Get true energy
753  GEnergy etrue = m_rmf.etrue().elogmean(itrue);
754 
755  // Get integration boundaries
756  GEbounds ebounds = ereco_bounds(etrue);
757 
758  // Integrate contributions over energy boundaries
759  double sum = 0.0;
760  for (int i = 0; i < ebounds.size(); ++i) {
761  sum += prob_erecobin(ebounds.emin(i), ebounds.emax(i), etrue, 0.0);
762  }
763 
764  // Store matrix sum
765  row_sums.push_back(sum);
766 
767  } // endfor: looped over all true energies
768 
769  // Normalize matrix elements
770  for (int itrue = 0; itrue < rows; ++itrue) {
771  double sum = row_sums[itrue];
772  if (sum > 0.0) {
773  for (int ireco = 0; ireco < columns; ++ireco) {
774  m_matrix(itrue, ireco) /= sum;
775  }
776  }
777  }
778 
779  // Return
780  return;
781 }
782 
783 
784 /***********************************************************************//**
785  * @brief Set interpolation cache
786  *
787  * Sets the interpolation cache.
788  ***************************************************************************/
789 void GCTAEdispRmf::set_cache(void) const
790 {
791  // Clear node arrays
792  m_etrue.clear();
793  m_ereco.clear();
794 
795  // Set log10(Etrue) nodes
796  for (int i = 0; i < m_rmf.ntrue(); ++i) {
798  }
799 
800  // Set log10(Emeasured) nodes
801  for (int i = 0; i < m_rmf.nmeasured(); ++i) {
802  #if defined(G_LINEAR_ERECO_INTERPOLATION)
804  #else
806  #endif
807  }
808 
809  // Return
810  return;
811 }
812 
813 
814 /***********************************************************************//**
815  * @brief Set maximum energy dispersion value
816  ***************************************************************************/
818 {
819  // Initialise maximum
820  m_max_edisp = 0.0;
821 
822  // Loop over all response table elements
823  for (int itrue = 0; itrue < m_matrix.rows(); ++itrue) {
824  for (int ireco = 0; ireco < m_matrix.columns(); ++ireco) {
825  double value = m_matrix(itrue, ireco);
826  if (value > m_max_edisp) {
827  m_max_edisp = value;
828  }
829  }
830  }
831 
832  // Return
833  return;
834 }
835 
836 
837 /***********************************************************************//**
838  * @brief Update cache
839  *
840  * @param[in] ereco Reconstructed energy.
841  * @param[in] etrue True energy.
842  *
843  * Updates the interpolation cache. The interpolation cache is composed
844  * of four indices and weights that define 4 data values of the RMF matrix
845  * that are used for bilinear interpolation.
846  ***************************************************************************/
847 void GCTAEdispRmf::update(const GEnergy& ereco, const GEnergy& etrue) const
848 {
849  // Update cache only of arguments have changed
850  if (ereco != m_last_etrue || etrue != m_last_ereco) {
851 
852  // Store actual values
853  m_last_ereco = ereco;
854  m_last_etrue = etrue;
855 
856  // Get log10 of true and reconstructed photon energies
857  double logEsrc = etrue.log10TeV();
858  #if defined(G_LINEAR_ERECO_INTERPOLATION)
859  double logEobs = ereco.TeV();
860  #else
861  double logEobs = ereco.log10TeV();
862  #endif
863 
864  // Set values for node arrays
865  m_etrue.set_value(logEsrc);
866  m_ereco.set_value(logEobs);
867 
868  // Set indices for bi-linear interpolation
873 
874  // Set weighting factors for bi-linear interpolation
879 
880  } // endif: update requested
881 
882  // Return
883  return;
884 }
885 
886 
887 /***********************************************************************//**
888  * @brief Compute m_ereco_bounds vector
889  ***************************************************************************/
891 {
892  // Clear boundaries
893  m_ereco_bounds.clear();
894 
895  // Loop over Etrue
896  for (int i = 0; i < m_rmf.ntrue(); ++i) {
897 
898  // Get true photon energy
899  GEnergy etrue = m_rmf.etrue().elogmean(i);
900 
901  // Add ebounds to vector
902  m_ereco_bounds.push_back(m_rmf.emeasured(etrue));
903 
904  } // endfor: looped over true photon energies
905 
906  // Return
907  return;
908 }
909 
910 
911 /***********************************************************************//**
912  * @brief Compute m_etrue_bounds vector
913  ***************************************************************************/
915 {
916  // Clear boundaries
917  m_etrue_bounds.clear();
918 
919  // Loop over Eobs
920  for (int i = 0; i < m_rmf.nmeasured(); ++i) {
921 
922  // Get reconstructed energy
923  GEnergy ereco = m_rmf.emeasured().elogmean(i);
924 
925  // Add ebounds to vector
926  m_etrue_bounds.push_back(m_rmf.emeasured(ereco));
927 
928  } // endfor: looped over measured photon energies
929 
930  // Return
931  return;
932 }
virtual ~GCTAEdispRmf(void)
Destructor.
int size(void) const
Return number of nodes in node array.
Definition: GNodeArray.hpp:188
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.
#define G_MC
Energy value class definition.
Random number generator class definition.
Sparse matrix class interface definition.
void clear(void)
Clear instance.
void set_max_edisp(void)
Set maximum energy dispersion value.
int size(void) const
Return number of energy boundaries.
Definition: GEbounds.hpp:157
int m_itrue2
Index of right Etrue.
void free_members(void)
Delete class members.
Definition: GCTAEdisp.cpp:164
const double & wgt_left(void) const
Returns left node weight.
Definition: GNodeArray.hpp:260
Abstract base class for the CTA energy dispersion.
Definition: GCTAEdisp.hpp:49
void clear(void)
Clear node array.
Definition: GNodeArray.cpp:235
CTA Redistribution Matrix File (RMF) energy dispersion class.
double TeV(void) const
Return energy in TeV.
Definition: GEnergy.cpp:348
double sum(const GVector &vector)
Computes vector sum.
Definition: GVector.cpp:901
double log10TeV(void) const
Return log10 of energy in TeV.
Definition: GEnergy.cpp:457
Random number generator class.
Definition: GRan.hpp:44
void init_members(void)
Initialise class members.
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
double m_wgt3
Weight of lower right node.
GCTAEdispRmf(void)
Void constructor.
int nmeasured(void) const
Return number of measured energy bins in redistribution matrix.
Definition: GRmf.hpp:207
Gammalib tools definition.
void set_value(const double &value) const
Set indices and weighting factors for interpolation.
Definition: GNodeArray.cpp:580
virtual void clear(void)
Clear matrix.
void copy_members(const GCTAEdispRmf &psf)
Copy class members.
void free_members(void)
Delete class members.
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.
const GEbounds & emeasured(void) const
Return measured energy boundaries.
Definition: GRmf.hpp:265
void compute_ereco_bounds(void) const
Compute m_ereco_bounds vector.
GMatrixSparse m_matrix
Normalised redistribution matrix.
int m_itrue1
Index of left Etrue.
GEnergy m_last_ereco
Last reconstructed energy.
Node array class interface definition.
XSPEC Redistribution Matrix File class definition.
const double & wgt_right(void) const
Returns right node weight.
Definition: GNodeArray.hpp:275
Energy boundaries container class.
Definition: GEbounds.hpp:60
Single parameter function abstract base class definition.
std::string print(const GChatter &chatter=NORMAL) const
Print energy.
Definition: GEnergy.cpp:748
GFilename m_filename
Name of response file.
Filename class.
Definition: GFilename.hpp:62
const int & rows(void) const
Return number of matrix rows.
std::vector< GEbounds > m_etrue_bounds
const GEnergy & emin(void) const
Return minimum energy of all intervals.
Definition: GEbounds.hpp:181
GNodeArray m_etrue
Array of log10(Etrue)
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
GCTAEdispRmf & operator=(const GCTAEdispRmf &edisp)
Assignment operator.
CTA RMF energy dispersion class definition.
int ntrue(void) const
Return number of true energy bins in redistribution matrix.
Definition: GRmf.hpp:193
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.
GChatter
Definition: GTypemaps.hpp:33
GEnergy elogmean(const int &index) const
Returns logarithmic mean energy for a given energy interval.
Definition: GEbounds.cpp:1140
const int & inx_left(void) const
Returns left node index.
Definition: GNodeArray.hpp:231
bool m_ereco_bounds_computed
void clear(void)
Clear object.
Definition: GRmf.cpp:315
void load(const GFilename &filename)
Load energy dispersion from RMF file.
GFilename filename(void) const
Return filename.
const int & inx_right(void) const
Returns right node index.
Definition: GNodeArray.hpp:245
GCTAEdispRmf * clone(void) const
Clone instance.
int m_ireco2
Index of right Ereco.
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 .
GCTAEdisp & operator=(const GCTAEdisp &edisp)
Assignment operator.
Definition: GCTAEdisp.cpp:106
bool m_etrue_bounds_computed
void set_cache(void) const
Set interpolation cache.
void clear(void)
Clear file name.
Definition: GFilename.cpp:188
void update(const GEnergy &ereco, const GEnergy &etrue) const
Update cache.
GEnergy m_last_etrue
Last true energy.
double m_wgt1
Weight of lower left node.
Sparse matrix class definition.
double m_wgt2
Weight of upper left node.
void compute_etrue_bounds(void) const
Compute m_etrue_bounds vector.
GEnergy m_last_etrue_bounds
void init_members(void)
Initialise class members.
Definition: GCTAEdisp.cpp:142
GNodeArray m_ereco
Array of log10(Ereco)
GEnergy m_last_ereco_bounds
const GEnergy & emax(void) const
Return maximum energy of all intervals.
Definition: GEbounds.hpp:193
const GEbounds & etrue(void) const
Return true energy boundaries.
Definition: GRmf.hpp:251
int m_ireco1
Index of left Ereco.
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
double m_wgt4
Weight of upper right node.
void set_matrix(void)
Set redistribution matrix.
GRmf m_rmf
Redistribution matrix file.
double m_max_edisp
Maximum energy dispersion value for MC.
Integration class interface definition.
std::vector< GEbounds > m_ereco_bounds
const int & columns(void) const
Return number of matrix columns.
std::string print(const GChatter &chatter=NORMAL) const
Print RMF information.
void clear(void)
Clear instance.
Definition: GEnergy.cpp:261
Filename class interface definition.
void load(const GFilename &filename)
Load Redistribution Matrix File.
Definition: GRmf.cpp:506
Mathematical function definitions.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
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 str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413