GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAEdisp2D.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAEdisp2D.cpp - CTA 2D energy dispersion class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2015-2019 by Florent Forest *
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 GCTAEdisp2D.cpp
23  * @brief CTA 2D energy dispersion class implementation
24  * @author Florent Forest
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <cmath>
32 #include <vector>
33 #include "GTools.hpp"
34 #include "GMath.hpp"
35 #include "GNdarray.hpp"
36 #include "GFft.hpp"
37 #include "GException.hpp"
38 #include "GIntegral.hpp"
39 #include "GFilename.hpp"
40 #include "GRan.hpp"
41 #include "GFits.hpp"
42 #include "GFitsTable.hpp"
43 #include "GFitsBinTable.hpp"
44 #include "GCTAEdisp2D.hpp"
45 #include "GCTASupport.hpp"
46 
47 /* __ Method name definitions ____________________________________________ */
48 #define G_READ "GCTAEdisp2D::read(GFitsTable&)"
49 #define G_MC "GCTAEdisp2D::mc(GRan&, GEnergy&, double&, double&, double&, "\
50  "double&)"
51 #define G_FETCH "GCTAEdisp2D::fetch()"
52 
53 /* __ Macros _____________________________________________________________ */
54 
55 /* __ Coding definitions _________________________________________________ */
56 #define G_SMOOTH_EDISP_KLUDGE //!< Get rid of noise in energy dispersion
57 //#define G_SMOOTH_EDISP_SAVE_TEST //!< Save test matrix
58 
59 /* __ Debug definitions __________________________________________________ */
60 //#define G_DEBUG_HESS_RENORMALIZATION
61 
62 /* __ Constants __________________________________________________________ */
63 
64 
65 /*==========================================================================
66  = =
67  = Constructors/destructors =
68  = =
69  ==========================================================================*/
70 
71 /***********************************************************************//**
72  * @brief Void constructor
73  *
74  * Constructs empty energy dispersion.
75  ***************************************************************************/
77 {
78  // Initialise class members
79  init_members();
80 
81  // Return
82  return;
83 }
84 
85 
86 /***********************************************************************//**
87  * @brief File constructor
88  *
89  * @param[in] filename FITS file name.
90  *
91  * Constructs energy dispersion from a FITS file.
92  ***************************************************************************/
94 {
95 
96  // Initialise class members
97  init_members();
98 
99  // Load energy dispersion from file
100  load(filename);
101 
102  // Return
103  return;
104 }
105 
106 
107 /***********************************************************************//**
108  * @brief Copy constructor
109  *
110  * @param[in] edisp Energy dispersion.
111  *
112  * Constructs energy dispersion by copying from another energy dispersion.
113  ***************************************************************************/
115 {
116  // Initialise class members
117  init_members();
118 
119  // Copy members
120  copy_members(edisp);
121 
122  // Return
123  return;
124 }
125 
126 
127 /***********************************************************************//**
128  * @brief Destructor
129  *
130  * Destructs energy dispersion.
131  ***************************************************************************/
133 {
134  // Free members
135  free_members();
136 
137  // Return
138  return;
139 }
140 
141 
142 /*==========================================================================
143  = =
144  = Operators =
145  = =
146  ==========================================================================*/
147 
148 /***********************************************************************//**
149  * @brief Assignment operator
150  *
151  * @param[in] edisp Energy dispersion.
152  * @return Energy dispersion.
153  *
154  * Assigns energy dispersion.
155  ***************************************************************************/
157 {
158  // Execute only if object is not identical
159  if (this != &edisp) {
160 
161  // Copy base class members
162  this->GCTAEdisp::operator=(edisp);
163 
164  // Free members
165  free_members();
166 
167  // Initialise private members
168  init_members();
169 
170  // Copy members
171  copy_members(edisp);
172 
173  } // endif: object was not identical
174 
175  // Return this object
176  return *this;
177 }
178 
179 
180 /***********************************************************************//**
181  * @brief Return energy dispersion in units of MeV\f$^{-1}\f$
182  *
183  * @param[in] ereco Reconstructed photon energy.
184  * @param[in] etrue True photon energy.
185  * @param[in] theta Offset angle in camera system (radians).
186  * @param[in] phi Azimuth angle in camera system (radians). Not used.
187  * @param[in] zenith Zenith angle in Earth system (radians). Not used.
188  * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
189  * @return Energy dispersion (MeV\f$^{-1}\f$)
190  *
191  * Returns the energy dispersion
192  *
193  * \f[
194  * E_{\rm disp}(E_{\rm reco} | E_{\rm true}, \theta) =
195  * \frac{E_{\rm disp}(m | E_{\rm true}, \theta)}{E_{\rm true}}
196  * \f]
197  *
198  * in units of MeV\f$^{-1}\f$ where
199  * \f$E_{\rm reco}\f$ is the reconstructed energy,
200  * \f$E_{\rm true}\f$ is the true energy, and
201  * \f$\theta\f$ is the offset angle.
202  ***************************************************************************/
203 double GCTAEdisp2D::operator()(const GEnergy& ereco,
204  const GEnergy& etrue,
205  const double& theta,
206  const double& phi,
207  const double& zenith,
208  const double& azimuth) const
209 {
210  // Make sure that energy dispersion is online
211  fetch();
212 
213  // Initalize edisp
214  double edisp = 0.0;
215 
216  // Get log10 of true photon energy
217  double logEsrc = etrue.log10TeV();
218 
219  // Continue only if logEsrc and theta are in validity range
220  if ((logEsrc >= m_logEsrc_min) && (logEsrc <= m_logEsrc_max) &&
221  (theta >= m_theta_min) && (theta <= m_theta_max)) {
222 
223  // Compute migration (migra=Ereco/Etrue)
224  double migra = ereco.MeV() / etrue.MeV();
225 
226  // Continue only if migration is in validity range
227  if ((migra >= m_migra_min) && (migra <= m_migra_max)) {
228 
229  // Setup argument vector
230  double arg[3];
231  arg[m_inx_etrue] = logEsrc;
232  arg[m_inx_migra] = migra;
233  arg[m_inx_theta] = theta;
234 
235  // Compute edisp
236  edisp = m_edisp(m_inx_matrix, arg[0], arg[1], arg[2]);
237 
238  // Make sure that energy dispersion is non-negative
239  if (edisp < 0.0) {
240  edisp = 0.0;
241  }
242 
243  // If energy dispersion is positive then divide the energy
244  // dispersion value by true photon energy to get the energy
245  // dispersion per MeV
246  else {
247  edisp /= etrue.MeV();
248  }
249 
250  } // endif: migration was in valid range
251 
252  } // endif: true energy and offset angle was in valid range
253 
254  // Return
255  return edisp;
256 }
257 
258 
259 /*==========================================================================
260  = =
261  = Public methods =
262  = =
263  ==========================================================================*/
264 
265 /***********************************************************************//**
266  * @brief Clear energy dispersion
267  *
268  * Clears energy dispersion.
269  ***************************************************************************/
271 {
272  // Free class members
273  free_members();
274  this->GCTAEdisp::free_members();
275 
276  // Initialise members
277  this->GCTAEdisp::init_members();
278  init_members();
279 
280  // Return
281  return;
282 }
283 
284 
285 /***********************************************************************//**
286  * @brief Clone energy dispersion
287  *
288  * @return Deep copy of energy dispersion.
289  *
290  * Returns a pointer to a deep copy of the point spread function.
291  ***************************************************************************/
293 {
294  return new GCTAEdisp2D(*this);
295 }
296 
297 
298 /***********************************************************************//**
299  * @brief Assign response table
300  *
301  * @param[in] table Response table.
302  *
303  * Assigns the response table for an energy dispersion.
304  ***************************************************************************/
306 {
307  // Assign response table
308  m_edisp = table;
309 
310  // Set table
311  set_table();
312 
313  // Return
314  return;
315 }
316 
317 
318 /***********************************************************************//**
319  * @brief Read energy dispersion from FITS table
320  *
321  * @param[in] table FITS table.
322  *
323  * @exception GException::invalid_value
324  * Response table is not three-dimensional.
325  *
326  * Reads the energy dispersion form the FITS @p table. The following column
327  * names are mandatory:
328  *
329  * ENERG_LO - True energy lower bin boundaries (alternative name: ETRUE_LO)
330  * ENERG_HI - True energy upper bin boundaries (alternative name: ETRUE_HI)
331  * MIGRA_LO - Migration lower bin boundaries
332  * MIGRA_HI - Migration upper bin boundaries
333  * THETA_LO - Offset angle lower bin boundaries
334  * THETA_HI - Offset angle upper bin boundaries
335  * MATRIX - Migration matrix
336  *
337  * The data are stored in the m_edisp member. The energy axis will be set to
338  * \f$log_{10} E_{\rm true}\f$, the offset angle axis to radians.
339  *
340  * The method assures that the energy dispersion is properly normalised.
341  ***************************************************************************/
342 void GCTAEdisp2D::read(const GFitsTable& table)
343 {
344  // Clear response table
345  m_edisp.clear();
346 
347  // Read energy dispersion table
348  m_edisp.read(table);
349 
350  // Throw an exception if the table is not two-dimensional
351  if (m_edisp.axes() != 3) {
352  std::string msg = "Expected three-dimensional energy dispersion "
353  "response table but found "+
355  " dimensions. Please specify a three-dimensional "
356  "energy dispersion.";
357  throw GException::invalid_value(G_READ, msg);
358  }
359 
360  // Set table
361  set_table();
362 
363  // Return
364  return;
365 }
366 
367 
368 /***********************************************************************//**
369  * @brief Write energy dispersion into FITS binary table
370  *
371  * @param[in] table FITS binary table.
372  *
373  * Writes the energy dispersion into the FITS binary @p table.
374  *
375  * @todo Add keywords.
376  ***************************************************************************/
378 {
379  // Make sure that energy dispersion is online
380  fetch();
381 
382  // Write background table
383  m_edisp.write(table);
384 
385  // Return
386  return;
387 }
388 
389 
390 /***********************************************************************//**
391  * @brief Load energy dispersion from FITS file
392  *
393  * @param[in] filename FITS file name.
394  *
395  * Loads the energy dispersion from a FITS file.
396  *
397  * The method does not actually load the FITS file, which will only be loaded
398  * on request. Only the filename is stored. See the fetch() method for the
399  * actually loading of the energy dispersion.
400  ***************************************************************************/
401 void GCTAEdisp2D::load(const GFilename& filename)
402 {
403  // Clear object
404  clear();
405 
406  // Store filename
408 
409  // Return
410  return;
411 }
412 
413 
414 /***********************************************************************//**
415  * @brief Save energy dispersion table into FITS file
416  *
417  * @param[in] filename FITS file name.
418  * @param[in] clobber Overwrite existing file?
419  *
420  * Saves energy dispersion into a FITS file. If a file with the given
421  * @p filename does not yet exist it will be created, otherwise the method
422  * opens the existing file. The method will create a (or replace an existing)
423  * energy dispersion extension. The extension name can be specified as part
424  * of the @p filename, or if no extension name is given, is assumed to be
425  * `ENERGY DISPERSION`.
426  *
427  * An existing file will only be modified if the @p clobber flag is set to
428  * true.
429  ***************************************************************************/
430 void GCTAEdisp2D::save(const GFilename& filename, const bool& clobber) const
431 {
432  // Get extension name
433  std::string extname = filename.extname(gammalib::extname_cta_edisp2d);
434 
435  // Open or create FITS file (without extension name since the requested
436  // extension may not yet exist in the file)
437  GFits fits(filename.url(), true);
438 
439  // Remove extension if it exists already
440  if (fits.contains(extname)) {
441  fits.remove(extname);
442  }
443 
444  // Create binary table
446 
447  // Write the background table
448  write(table);
449 
450  // Set binary table extension name
451  table.extname(extname);
452 
453  // Append table to FITS file
454  fits.append(table);
455 
456  // Save to file
457  fits.save(clobber);
458 
459  // Return
460  return;
461 }
462 
463 
464 /***********************************************************************//**
465  * @brief Simulate energy dispersion
466  *
467  * @param[in] ran Random number generator.
468  * @param[in] etrue True photon energy.
469  * @param[in] theta Offset angle in camera system (radians).
470  * @param[in] phi Azimuth angle in camera system (radians). Not used.
471  * @param[in] zenith Zenith angle in Earth system (radians). Not used.
472  * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
473  * @return Reconstructed energy.
474  *
475  * @exception GException::invalid_return_value
476  * No energy dispersion information found for parameters or
477  * energy dispersion matrix is empty.
478  *
479  * Draws reconstructed energy value given a true energy @p etrue and offset
480  * angle @p theta. If no energy dispersion information is available the
481  * method will return the true photon energy.
482  ***************************************************************************/
484  const GEnergy& etrue,
485  const double& theta,
486  const double& phi,
487  const double& zenith,
488  const double& azimuth) const
489 {
490  // Make sure that energy dispersion is online
491  fetch();
492 
493  // Initialise reconstructed event energy with true photon energy
494  GEnergy ereco = etrue;
495 
496  // Get boundaries for reconstructed energy
497  GEbounds ebounds = ereco_bounds(etrue, theta, phi, zenith, azimuth);
498  double emin = ebounds.emin().TeV();
499  double emax = ebounds.emax().TeV();
500 
501  // Throw an exception if minimum energy is equal or larger than maximum
502  // energy (this means that no energy dispersion information is available
503  // for that energy)
504  if (emin >= emax) {
505  std::string msg = "No valid energy dispersion information available "
506  "for true photon energy "+etrue.print()+", "
507  "offset angle "+
508  gammalib::str(theta*gammalib::rad2deg)+" deg "
509  "and azimuth angle "+
510  gammalib::str(phi*gammalib::rad2deg)+" deg. Either "
511  "provide an energy dispersion matrix covering these "
512  "parameters or restrict the parameter space.";
514  }
515 
516  // Get maximum energy dispersion value (including some margin)
517  double max_edisp = 1.5 * get_max_edisp(etrue, theta);
518 
519  // Throw an exception if maximum energy dispersion is zero
520  if (max_edisp <= 0.0) {
521  std::string msg = "Energy dispersion matrix is empty. Please provide "
522  "a valid energy dispersion matrix.";
524  }
525 
526  // Initialise rejection method
527  double ewidth = emax - emin;
528  double f = 0.0;
529  double ftest = 1.0;
530  int zeros = 0;
531  const int max_subsequent_zeros = 10;
532 
533  // Find energy by rejection method
534  while (ftest > f) {
535 
536  // Draw random observed energy and evaluate energy dispersion matrix
537  // until a non-zero energy dispersion value is found. Subsequent zero
538  // energy dispersion values may indicate that there is no valid energy
539  // dispersion information available. After a limited number of
540  // subsequent zeros the loop is terminated with an exception
541  int zeros = 0;
542  do {
543  ereco.TeV(emin + ewidth * ran.uniform());
544  f = operator()(ereco, etrue, theta, phi, zenith, azimuth);
545  if (f == 0.0) {
546  zeros++;
547  if (zeros > max_subsequent_zeros) {
548  std::string msg = "No valid energy dispersion information "
549  "available for true photon energy "+
550  etrue.print()+", offset angle "+
552  " deg and azimuth angle "+
553  gammalib::str(phi*gammalib::rad2deg)+
554  " deg. Either provide an energy "
555  "dispersion matrix covering these "
556  "parameters or restrict the parameter "
557  "space.";
559  }
560  }
561  } while (f == 0);
562 
563  // Get uniform random value between zero and the maximum energy
564  // dispersion value. The loop is quit if the random number is smaller
565  // than the energy dispersion value
566  ftest = ran.uniform() * max_edisp;
567 
568  } // endwhile:
569 
570  // Return energy
571  return ereco;
572 }
573 
574 
575 /***********************************************************************//**
576  * @brief Return observed energy interval that contains the energy dispersion.
577  *
578  * @param[in] etrue True photon energy.
579  * @param[in] theta Offset angle in camera system (radians).
580  * @param[in] phi Azimuth angle in camera system (radians). Not used.
581  * @param[in] zenith Zenith angle in Earth system (radians). Not used.
582  * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
583  * @return Observed energy boundaries.
584  *
585  * Returns the band of observed energies outside of which the energy
586  * dispersion becomes negligible for a given true energy @p etrue and
587  * offset angle @p theta.
588  ***************************************************************************/
590  const double& theta,
591  const double& phi,
592  const double& zenith,
593  const double& azimuth) const
594 {
595  // Make sure that energy dispersion is online
596  fetch();
597 
598  // Compute only if parameters changed
599  if (!m_ereco_bounds_computed || theta != m_last_theta_ereco) {
600 
601  // Set computation flag
604  m_last_etrue.TeV(0.0); // force update
605 
606  // Compute reconstructed energy bounds
607  compute_ereco_bounds(theta, phi, zenith, azimuth);
608 
609  }
610 
611  // Search index only if true photon energy has changed
612  if (etrue != m_last_etrue) {
613 
614  // Store last true photon energy
616 
617  // Get true energy in TeV
618  double etrue_TeV = etrue.TeV();
619 
620  // Find right index with bisection
621  int low = 0;
622  int high = m_ereco_bounds.size() - 1;
623  while ((high-low) > 1) {
624  int mid = (low+high) / 2;
625  double e = m_edisp.axis_lo(m_inx_etrue, mid);
626  if (etrue_TeV < e) {
627  high = mid;
628  }
629  else {
630  low = mid;
631  }
632  }
633 
634  // Index found
635  m_index_ereco = low;
636 
637  }
638 
639  // Return energy boundaries
641 }
642 
643 
644 /***********************************************************************//**
645  * @brief Return true energy interval that contains the energy dispersion.
646  *
647  * @param[in] ereco Reconstructed event energy.
648  * @param[in] theta Offset angle in camera system (radians).
649  * @param[in] phi Azimuth angle in camera system (radians). Not used.
650  * @param[in] zenith Zenith angle in Earth system (radians). Not used.
651  * @param[in] azimuth Azimuth angle in Earth system (radians). Not used.
652  * @return True energy boundaries.
653  *
654  * Returns the band of true photon energies outside of which the energy
655  * dispersion becomes negligible for a given observed energy @p ereco and
656  * offset angle @p theta.
657  ***************************************************************************/
659  const double& theta,
660  const double& phi,
661  const double& zenith,
662  const double& azimuth) const
663 {
664  // Make sure that energy dispersion is online
665  fetch();
666 
667  // Compute true energy boundaries only if parameters changed
668  if (!m_etrue_bounds_computed || theta != m_last_theta_etrue) {
669 
670  // Set computation flag
673  m_last_ereco.TeV(0.0); // force update
674 
675  // Compute true energy boundaries
676  compute_etrue_bounds(theta, phi, zenith, azimuth);
677 
678  }
679 
680  // Search index only if reconstructed event energy has changed
681  if (ereco != m_last_ereco) {
682 
683  // Store reconstructed event energy
684  m_last_ereco = ereco;
685 
686  // Get reconstructed event energy in TeV
687  double ereco_TeV = ereco.TeV();
688 
689  // Find right index with bisection
690  int low = 0;
691  int high = m_etrue_bounds.size() - 1;
692  while ((high-low) > 1) {
693  int mid = (low+high) / 2;
694  double e = m_edisp.axis_lo(m_inx_etrue, mid);
695  if (ereco_TeV < e) {
696  high = mid;
697  }
698  else {
699  low = mid;
700  }
701  }
702 
703  // Index found
704  m_index_etrue = low;
705 
706  }
707 
708  // Return energy boundaries
710 }
711 
712 
713 /***********************************************************************//**
714  * @brief Fetch energy dispersion
715  *
716  * @exception GException::file_error
717  * File not found.
718  * Unable to load energy dispersion.
719  *
720  * Fetches the energy dispersion by reading it from a FITS file.
721  *
722  * If the filename contains no extension name the method scans the `HDUCLASS`
723  * keywords of all extensions and loads the energy dispersion from the first
724  * extension for which `HDUCLAS4=EDISP_2D`.
725  *
726  * Otherwise, the background will be loaded from the `ENERGY DISPERSION`
727  * extension.
728  *
729  * This method does nothing if the energy dispersion is already loaded, if
730  * there is nothing to fetch, or if the m_filename member is empty.
731  *
732  * The method is thread save. The method checks whether the file from which
733  * the energy dispersion should be loaded actually exists.
734  ***************************************************************************/
735 void GCTAEdisp2D::fetch(void) const
736 {
737  // Continue only if energy dispersion has not yet been fetched
738  if (!m_fetched) {
739 
740  // Continue only if the file name is not empty
741  if (!m_filename.is_empty()) {
742 
743  // Throw an exception if the file does not exist
744  if (!m_filename.exists()) {
745  std::string msg = "File \""+m_filename+"\" not found. Cannot "
746  "fetch energy dispersion. Maybe the file has "
747  "been deleted in the meantime.";
749  }
750 
751  // Signal that energy dispersion will be fetched (has to come
752  // before reading since the ebounds_obs() and ebounds_src() methods
753  // will otherwise call the fetch() method recursively.
754  m_fetched = true;
755 
756  // Initialise exception flag
757  bool has_exception = false;
758 
759  // Load energy dispersion. Catch any exception. Put the code into
760  // a critical zone as it might be called from within a parallelized
761  // thread.
762  #pragma omp critical(GCTAEdisp2D_fetch)
763  {
764  try {
765 
766 
767  // Open FITS file
768  GFits fits(m_filename);
769 
770  // Get the default extension name. If no GADF compliant name
771  // was found then set the default extension name to
772  // "ENERGY DISPERSION".
773  std::string extname = gammalib::gadf_hduclas4(fits, "EDISP_2D");
774  if (extname.empty()) {
776  }
777 
778  // Get energy dispersion table
779  const GFitsTable& table = *fits.table(m_filename.extname(extname));
780 
781  // Read energy dispersion from table
782  const_cast<GCTAEdisp2D*>(this)->read(table);
783 
784  // Close FITS file
785  fits.close();
786 
787  }
788  catch (...) {
789  has_exception = true;
790  }
791  }
792 
793  // Throw an exception if an exception has occured
794  if (has_exception) {
795  std::string msg = "Unable to load energy dispersion from "
796  "file \""+m_filename+"\".";
797  throw GException::file_error(G_FETCH, msg);
798  }
799 
800  } // endif: filename was not empty
801 
802  } // endif: energy dispersion had not yet been fetched
803 
804  // Return
805  return;
806 }
807 
808 
809 /***********************************************************************//**
810  * @brief Return energy dispersion probability for reconstructed energy
811  * interval
812  *
813  * @param[in] ereco_min Minimum of reconstructed energy interval.
814  * @param[in] ereco_max Maximum of reconstructed energy interval.
815  * @param[in] etrue True energy.
816  * @param[in] theta Offset angle (radians).
817  * @return Integrated energy dispersion probability.
818  *
819  * Computes
820  *
821  * \f[
822  * \int_{E_{\rm reco}^{\rm min}}^{E_{\rm reco}^{\rm max}}
823  * E_{\rm disp}(E_{\rm reco} | E_{\rm true}, \theta) \, dE_{\rm reco}
824  * \f]
825  *
826  * where
827  * \f$E_{\rm reco}\f$ is the reconstructed energy,
828  * \f$E_{\rm true}\f$ is the true energy and
829  * \f$\theta\f$ is the offset angle.
830  *
831  * The method takes into account that the energy dispersion data are stored
832  * in a matrix and uses the trapezoidal rule for integration of the tabulated
833  * data. This assures the fastest possible integration.
834  ***************************************************************************/
835 double GCTAEdisp2D::prob_erecobin(const GEnergy& ereco_min,
836  const GEnergy& ereco_max,
837  const GEnergy& etrue,
838  const double& theta) const
839 {
840  // Make sure that energy dispersion is online
841  fetch();
842 
843  // Initalize probability
844  double prob = 0.0;
845 
846  // Get log10 of true energy
847  double logEsrc = etrue.log10TeV();
848 
849  // Get migration limits
850  double migra_min = ereco_min / etrue;
851  double migra_max = ereco_max / etrue;
852 
853  // Continue only if logEsrc and theta are in validity range and migration
854  // interval overlaps with response table
855  if ((logEsrc >= m_logEsrc_min) && (logEsrc <= m_logEsrc_max) &&
856  (theta >= m_theta_min) && (theta <= m_theta_max) &&
857  (migra_max > m_migra_min) && (migra_min < m_migra_max)) {
858 
859  // Constrain migration limits to response table
860  if (migra_min < m_migra_min) {
861  migra_min = m_migra_min;
862  }
863  if (migra_max > m_migra_max) {
864  migra_max = m_migra_max;
865  }
866 
867  // Retrieve references to node arrays
868  const GNodeArray& etrue_nodes = m_edisp.axis_nodes(m_inx_etrue);
869  const GNodeArray& theta_nodes = m_edisp.axis_nodes(m_inx_theta);
870  const GNodeArray& migra_nodes = m_edisp.axis_nodes(m_inx_migra);
871 
872  // Set logEsrc and theta for node array interpolation
873  etrue_nodes.set_value(logEsrc);
874  theta_nodes.set_value(theta);
875 
876  // Compute base indices
877  int base_ll = table_index(etrue_nodes.inx_left(), 0, theta_nodes.inx_left());
878  int base_lr = table_index(etrue_nodes.inx_left(), 0, theta_nodes.inx_right());
879  int base_rl = table_index(etrue_nodes.inx_right(), 0, theta_nodes.inx_left());
880  int base_rr = table_index(etrue_nodes.inx_right(), 0, theta_nodes.inx_right());
881 
882  // Get migration stride
883  int stride = table_stride(m_inx_migra);
884 
885  // Initialise first and second node
886  double x1 = migra_min;
887  double f1 = table_value(base_ll, base_lr, base_rl, base_rr,
888  etrue_nodes.wgt_left(),
889  etrue_nodes.wgt_right(),
890  theta_nodes.wgt_left(),
891  theta_nodes.wgt_right(), x1);
892  double x2 = 0.0;
893  double f2 = 0.0;
894 
895  // Loop over all migration nodes
896  for (int i = 0, offset = 0; i < migra_nodes.size(); ++i, offset += stride) {
897 
898  // If migration value is below migra_min then skip node
899  if (migra_nodes[i] <= migra_min) {
900  continue;
901  }
902 
903  // If migration value is above maximum migration value then use
904  // migra_max as migration value
905  if (migra_nodes[i] > migra_max) {
906  x2 = migra_max;
907  f2 = table_value(base_ll, base_lr, base_rl, base_rr,
908  etrue_nodes.wgt_left(),
909  etrue_nodes.wgt_right(),
910  theta_nodes.wgt_left(),
911  theta_nodes.wgt_right(), x2);
912  }
913 
914  // ... otherwise use migration value
915  else {
916  x2 = migra_nodes[i];
917  f2 = table_value(base_ll, base_lr, base_rl, base_rr,
918  etrue_nodes.wgt_left(),
919  etrue_nodes.wgt_right(),
920  theta_nodes.wgt_left(),
921  theta_nodes.wgt_right(), offset);
922  }
923 
924  // Compute integral
925  prob += 0.5 * (f1 + f2) * (x2 - x1);
926 
927  // Set second node as first node
928  x1 = x2;
929  f1 = f2;
930 
931  // If node energy is above migra_max then break now
932  if (migra_nodes[i] > migra_max) {
933  break;
934  }
935 
936  } // endfor: looped over all nodes
937 
938  // If last node energy is below migra_max then compute last part of
939  // integral up to emax
940  if (x1 < migra_max) {
941  x2 = migra_max;
942  f2 = table_value(base_ll, base_lr, base_rl, base_rr,
943  etrue_nodes.wgt_left(),
944  etrue_nodes.wgt_right(),
945  theta_nodes.wgt_left(),
946  theta_nodes.wgt_right(), x2);
947  prob += 0.5 * (f1 + f2) * (x2 - x1);
948  }
949 
950  } // endif: logEsrc, theta and migration range were valid
951 
952  // Return probability
953  return prob;
954 }
955 
956 
957 /***********************************************************************//**
958  * @brief Print energy dispersion information
959  *
960  * @param[in] chatter Chattiness.
961  * @return String containing energy dispersion information.
962  ***************************************************************************/
963 std::string GCTAEdisp2D::print(const GChatter& chatter) const
964 {
965  // Initialise result string
966  std::string result;
967 
968  // Continue only if chatter is not silent
969  if (chatter != SILENT) {
970 
971  // Append header
972  result.append("=== GCTAEdisp2D ===");
973  result.append("\n"+gammalib::parformat("Filename")+m_filename);
974 
975  // Make sure that energy dispersion is online
976  fetch();
977 
978  // Initialise information
979  int nebins = 0;
980  int nmigrabins = 0;
981  int nthetabins = 0;
982  double emin = 0.0;
983  double emax = 0.0;
984  double mmin = 0.0;
985  double mmax = 0.0;
986  double omin = 0.0;
987  double omax = 0.0;
988 
989  // Extract information if there are axes in the response table
990  if (m_edisp.axes() > 0) {
991  nebins = m_edisp.axis_bins(m_inx_etrue);
992  nmigrabins = m_edisp.axis_bins(m_inx_migra);
993  nthetabins = m_edisp.axis_bins(m_inx_theta);
994  emin = m_edisp.axis_lo(m_inx_etrue,0);
995  emax = m_edisp.axis_hi(m_inx_etrue,nebins-1);
996  mmin = m_edisp.axis_lo(m_inx_migra,0);
997  mmax = m_edisp.axis_hi(m_inx_migra,nmigrabins-1);
998  omin = m_edisp.axis_lo(m_inx_theta,0);
999  omax = m_edisp.axis_hi(m_inx_theta,nthetabins-1);
1000  }
1001 
1002  // Append information
1003  result.append("\n"+gammalib::parformat("Number of energy bins") +
1004  gammalib::str(nebins));
1005  result.append("\n"+gammalib::parformat("Number of migration bins") +
1006  gammalib::str(nmigrabins));
1007  result.append("\n"+gammalib::parformat("Number of offset bins") +
1008  gammalib::str(nthetabins));
1009  result.append("\n"+gammalib::parformat("Log10(Energy) range"));
1010  result.append(gammalib::str(emin)+" - "+gammalib::str(emax)+" TeV");
1011  result.append("\n"+gammalib::parformat("Migration range"));
1012  result.append(gammalib::str(mmin)+" - "+gammalib::str(mmax));
1013  result.append("\n"+gammalib::parformat("Offset angle range"));
1014  result.append(gammalib::str(omin)+" - "+gammalib::str(omax)+" deg");
1015 
1016  } // endif: chatter was not silent
1017 
1018  // Return result
1019  return result;
1020 }
1021 
1022 
1023 /*==========================================================================
1024  = =
1025  = Private methods =
1026  = =
1027  ==========================================================================*/
1028 
1029 /***********************************************************************//**
1030  * @brief Initialise class members
1031  ***************************************************************************/
1033 {
1034  // Initialise members
1035  m_filename.clear();
1036  m_edisp.clear();
1037  m_fetched = false;
1038  m_inx_etrue = 0;
1039  m_inx_migra = 1;
1040  m_inx_theta = 2;
1041  m_inx_matrix = 0;
1042  m_logEsrc_min = 0.0;
1043  m_logEsrc_max = 0.0;
1044  m_migra_min = 0.0;
1045  m_migra_max = 0.0;
1046  m_theta_min = 0.0;
1047  m_theta_max = 0.0;
1048  m_max_edisp.clear();
1049 
1050  // Initialise computation cache
1051  m_ereco_bounds_computed = false;
1052  m_etrue_bounds_computed = false;
1053  m_index_ereco = 0;
1054  m_index_etrue = 0;
1055  m_last_theta_ereco = -1.0;
1056  m_last_theta_etrue = -1.0;
1057  m_last_etrue.TeV(0.0);
1058  m_last_ereco.TeV(0.0);
1059  m_ereco_bounds.clear();
1060  m_etrue_bounds.clear();
1061 
1062  // Return
1063  return;
1064 }
1065 
1066 
1067 /***********************************************************************//**
1068  * @brief Copy class members
1069  *
1070  * @param[in] edisp Energy dispersion.
1071  ***************************************************************************/
1073 {
1074  // Copy members
1075  m_filename = edisp.m_filename;
1076  m_edisp = edisp.m_edisp;
1077  m_fetched = edisp.m_fetched;
1078  m_inx_etrue = edisp.m_inx_etrue;
1079  m_inx_migra = edisp.m_inx_migra;
1080  m_inx_theta = edisp.m_inx_theta;
1081  m_inx_matrix = edisp.m_inx_matrix;
1082  m_logEsrc_min = edisp.m_logEsrc_min;
1083  m_logEsrc_max = edisp.m_logEsrc_max;
1084  m_migra_min = edisp.m_migra_min;
1085  m_migra_max = edisp.m_migra_max;
1086  m_theta_min = edisp.m_theta_min;
1087  m_theta_max = edisp.m_theta_max;
1088  m_max_edisp = edisp.m_max_edisp;
1089 
1090  // Copy computation cache
1093  m_index_ereco = edisp.m_index_ereco;
1094  m_index_etrue = edisp.m_index_etrue;
1097  m_last_etrue = edisp.m_last_etrue;
1098  m_last_ereco = edisp.m_last_ereco;
1101 
1102  // Return
1103  return;
1104 }
1105 
1106 
1107 /***********************************************************************//**
1108  * @brief Delete class members
1109  ***************************************************************************/
1111 {
1112  // Return
1113  return;
1114 }
1115 
1116 
1117 /***********************************************************************//**
1118  * @brief Get true energy
1119  *
1120  * @param[in] ietrue True energy index.
1121  * @return True energy.
1122  ***************************************************************************/
1123 GEnergy GCTAEdisp2D::etrue(const int& ietrue) const
1124 {
1125  // Compute true energy in TeV
1126  double etrue_TeV = std::sqrt(m_edisp.axis_lo(m_inx_etrue, ietrue) *
1127  m_edisp.axis_hi(m_inx_etrue, ietrue));
1128 
1129  // Set true energy
1130  GEnergy etrue;
1131  etrue.TeV(etrue_TeV);
1132 
1133  // Return true energy
1134  return etrue;
1135 }
1136 
1137 
1138 /***********************************************************************//**
1139  * @brief Get migration
1140  *
1141  * @param[in] imigra Migration index.
1142  * @return Migration.
1143  ***************************************************************************/
1144 double GCTAEdisp2D::migra(const int& imigra) const
1145 {
1146  // Compute migration
1147  double migra = 0.5 * (m_edisp.axis_hi(m_inx_migra, imigra) +
1148  m_edisp.axis_lo(m_inx_migra, imigra));
1149 
1150  // Return migration
1151  return migra;
1152 }
1153 
1154 
1155 /***********************************************************************//**
1156  * @brief Get offset angle in radiaus
1157  *
1158  * @param[in] itheta Offset angle index.
1159  * @return Offset angle (radians).
1160  ***************************************************************************/
1161 double GCTAEdisp2D::theta(const int& itheta) const
1162 {
1163  // Compute offset angle in radians
1164  double theta = 0.5 * (m_edisp.axis_lo(m_inx_theta, itheta) +
1165  m_edisp.axis_hi(m_inx_theta, itheta)) *
1167 
1168  // Return offset angle
1169  return theta;
1170 }
1171 
1172 
1173 /***********************************************************************//**
1174  * @brief Compute vector of reconstructed energy bounds
1175  *
1176  * @param[in] theta Offset angle (radians).
1177  * @param[in] phi Azimuth angle (radians).
1178  * @param[in] zenith Zenith angle (radians).
1179  * @param[in] azimuth Azimuth angle (radians).
1180  *
1181  * Computes for all true energies the energy boundaries of the reconstructed
1182  * energies covered by valid migration matrix elements. Only matrix elements
1183  * with values >= 1.0e-12 are considered as valid elements. In case that no
1184  * matrix elements are found of a given true energy, the interval of observed
1185  * energies will be set to [0,0] (i.e. an empty interval).
1186  ***************************************************************************/
1187 void GCTAEdisp2D::compute_ereco_bounds(const double& theta,
1188  const double& phi,
1189  const double& zenith,
1190  const double& azimuth) const
1191 {
1192  // Clear ebounds_obs vector
1193  m_ereco_bounds.clear();
1194 
1195  // Set epsilon
1196  const double eps = 1.0e-12;
1197 
1198  // Determine number of true energy and migration bins
1199  int netrue = m_edisp.axis_bins(m_inx_etrue);
1200  int nmigra = m_edisp.axis_bins(m_inx_migra);
1201 
1202  // Loop over true photon energy
1203  for (int ietrue = 0; ietrue < netrue; ++ietrue) {
1204 
1205  // Compute true photon energy
1206  GEnergy etrue = this->etrue(ietrue);
1207 
1208  // Initialise results
1209  double ereco_min = 0.0;
1210  double ereco_max = 0.0;
1211  bool found_min = false;
1212  bool found_max = false;
1213 
1214  // Find minimum boundary
1215  for (int imigra = 0; imigra < nmigra; ++imigra) {
1216 
1217  // Compute migra value
1218  double migra = this->migra(imigra);
1219 
1220  // Get matrix term
1221  double arg[3];
1222  arg[m_inx_etrue] = etrue.log10TeV();
1223  arg[m_inx_migra] = migra;
1224  arg[m_inx_theta] = theta;
1225  double edisp = m_edisp(m_inx_matrix, arg[0], arg[1], arg[2]);
1226 
1227  // Find first non-negligible matrix term
1228  if (edisp >= eps) {
1229  found_min = true;
1230  ereco_min = migra * etrue.TeV();
1231  break;
1232  }
1233 
1234  } // endfor: find minimum boundary
1235 
1236  // Find maximum boundary
1237  for (int imigra = nmigra-1; imigra >= 0; imigra--) {
1238 
1239  // Compute migra value
1240  double migra = this->migra(imigra);
1241 
1242  // Get matrix term
1243  double arg[3];
1244  arg[m_inx_etrue] = etrue.log10TeV();
1245  arg[m_inx_migra] = migra;
1246  arg[m_inx_theta] = theta;
1247  double edisp = m_edisp(m_inx_matrix, arg[0], arg[1], arg[2]);
1248 
1249  // Find first non-negligible matrix term
1250  if (edisp >= eps) {
1251  found_max = true;
1252  ereco_max = migra * etrue.TeV();
1253  break;
1254  }
1255 
1256  } // endfor: find maximum boundary
1257 
1258  // If we did not find boundaries then set the interval to
1259  // a zero interval for safety
1260  if (!found_min || !found_max) {
1261  ereco_min = 0.0;
1262  ereco_max = 0.0;
1263  }
1264 
1265  // Set energy boundaries
1266  GEnergy emin;
1267  GEnergy emax;
1268  emin.TeV(ereco_min);
1269  emax.TeV(ereco_max);
1270 
1271  // Add energy boundaries to vector
1272  m_ereco_bounds.push_back(GEbounds(emin, emax));
1273 
1274  } // endfor: looped over logEsrc
1275 
1276  // Return
1277  return;
1278 }
1279 
1280 
1281 /***********************************************************************//**
1282  * @brief Compute vector of true energy boundaries
1283  *
1284  * @param[in] theta Offset angle (radians).
1285  * @param[in] phi Azimuth angle (radians).
1286  * @param[in] zenith Zenith angle (radians).
1287  * @param[in] azimuth Azimuth angle (radians).
1288  *
1289  * Computes for all observed energies the energy boundaries of the true
1290  * energies covered by valid migration matrix elements. Only matrix elements
1291  * with values >= 1.0e-12 are considered as valid elements. In case that no
1292  * matrix elements are found of a given observed energy, the interval of true
1293  * energies will be set to [0,0] (i.e. an empty interval).
1294  ***************************************************************************/
1295 void GCTAEdisp2D::compute_etrue_bounds(const double& theta,
1296  const double& phi,
1297  const double& zenith,
1298  const double& azimuth) const
1299 {
1300  // Initialise energies
1301  GEnergy etrue;
1302  GEnergy ereco;
1303 
1304  // Clear ebounds_src vector
1305  m_etrue_bounds.clear();
1306 
1307  // Set epsilon
1308  const double eps = 1.0e-12;
1309 
1310  // Determine number of true energy bins
1311  int nereco = m_edisp.axis_bins(m_inx_etrue);
1312  int nmigras = m_edisp.axis_bins(m_inx_migra);
1313 
1314  // Loop over reconstructed energy
1315  for (int iereco = 0; iereco < nereco; ++iereco) {
1316 
1317  // Set ereco (we use the true energy axis for the computation)
1318  GEnergy ereco = this->etrue(iereco);
1319 
1320  // Initialise results
1321  GEnergy etrue_min;
1322  GEnergy etrue_max;
1323  bool found_min = false;
1324  bool found_max = false;
1325 
1326  // Find minimum boundary
1327  for (int imigra = 0; imigra < nmigras; ++imigra) {
1328 
1329  // Compute migra value
1330  double migra = this->migra(imigra);
1331 
1332  // Fall through if migration is zero
1333  if (migra == 0.0) {
1334  continue;
1335  }
1336 
1337  // Compute true energy
1338  etrue = ereco / migra;
1339 
1340  // Get energy dispersion
1341  double edisp = operator()(ereco, etrue, theta);
1342 
1343  // Find first non-negligible matrix term
1344  if (edisp >= eps) {
1345  found_max = true;
1346  etrue_max = etrue;
1347  break;
1348  }
1349 
1350  } // endfor: find minimum boundary
1351 
1352  // Find maximum boundary
1353  for (int imigra = nmigras-1; imigra >= 0; imigra--) {
1354 
1355  // Compute migra value
1356  double migra = this->migra(imigra);
1357 
1358  // Fall through if migration is zero
1359  if (migra == 0.0) {
1360  continue;
1361  }
1362 
1363  // Compute true energy
1364  etrue = ereco / migra;
1365 
1366  // Get energy dispersion
1367  double edisp = operator()(ereco, etrue, theta);
1368 
1369  // Find first non-negligible matrix term
1370  if (edisp >= eps) {
1371  found_min = true;
1372  etrue_min = etrue;
1373  break;
1374  }
1375 
1376  } // endfor: find maximum boundary
1377 
1378  // If we did not find boundaries then set the interval to a zero
1379  // interval for safety
1380  if (!found_min || !found_max) {
1381  etrue_min.clear();
1382  etrue_max.clear();
1383  }
1384 
1385  // Add energy boundaries to vector
1386  m_etrue_bounds.push_back(GEbounds(etrue_min, etrue_max));
1387 
1388  } // endfor : loopend over observed energy
1389 
1390  // Return
1391  return;
1392 }
1393 
1394 
1395 /***********************************************************************//**
1396  * @brief Set table
1397  *
1398  * After assigning or loading a response table, performs all necessary steps
1399  * to set the table.
1400  *
1401  * For CTA, the method applies a kludge that smoothes the energy dispersion
1402  * table to get rid of numerical fluctuations.
1403  *
1404  * Sets the data members
1405  * m_inx_etrue
1406  * m_inx_migra
1407  * m_inx_theta
1408  * m_inx_matrix
1409  * m_logEsrc_min
1410  * m_logEsrc_max
1411  * m_migra_min
1412  * m_migra_max
1413  * m_theta_min
1414  * m_theta_max
1415  * m_max_edisp
1416  ***************************************************************************/
1418 {
1419  // Set table indices
1420  if (m_edisp.has_axis("ENERG")) {
1421  m_inx_etrue = m_edisp.axis("ENERG");
1422  }
1423  else {
1424  m_inx_etrue = m_edisp.axis("ETRUE"); // Old name, should not be used
1425  }
1426  m_inx_migra = m_edisp.axis("MIGRA");
1427  m_inx_theta = m_edisp.axis("THETA");
1428  m_inx_matrix = m_edisp.table("MATRIX");
1429 
1430  // Set true energy axis to logarithmic scale
1432 
1433  // Set offset angle axis to radians
1435 
1436  // Set table boundaries
1437  set_boundaries();
1438 
1439  // Smooth energy dispersion table (kludge only for CTA)
1440  #if defined(G_SMOOTH_EDISP_KLUDGE)
1441  if (gammalib::strip_whitespace(m_edisp.telescope()) == "CTA") {
1442  smooth_table();
1443  }
1444  #endif
1445 
1446  // Normalize energy dispersion table
1447  normalize_table();
1448 
1449  // Set maximum energy dispersion value array
1450  set_max_edisp();
1451 
1452  // Save test
1453  #if defined(G_SMOOTH_EDISP_SAVE_TEST)
1454  this->save("test_edisp.fits", true);
1455  #endif
1456 
1457  // Return
1458  return;
1459 }
1460 
1461 
1462 /***********************************************************************//**
1463  * @brief Set energy dispersion boundaries
1464  *
1465  * Sets the data members m_logEsrc_min, m_logEsrc_max, m_migra_min,
1466  * m_migra_max, m_theta_min and m_theta_max that define the validity range
1467  * of the energy dispersion.
1468  ***************************************************************************/
1470 {
1471  // Get number of bins
1472  int netrue = m_edisp.axis_bins(m_inx_etrue);
1473  int nmigra = m_edisp.axis_bins(m_inx_migra);
1474  int ntheta = m_edisp.axis_bins(m_inx_theta);
1475 
1476  // Compute minimum and maximum logEsrc boundaries
1479 
1480  // Compute minimum and maximum migration boundaries
1482  m_migra_max = m_edisp.axis_hi(m_inx_migra, nmigra-1);
1483 
1484  // Compute minimum and maximum theta boundaries
1487 
1488  // Return
1489  return;
1490 }
1491 
1492 
1493 /***********************************************************************//**
1494  * @brief Set array of maximum energy dispersion values
1495  *
1496  * Sets an internal array of maximum energy dispersion values for each given
1497  * true photon energy and offset angle in the response table. This array will
1498  * be used for Monte Carlo simulations.
1499  ***************************************************************************/
1501 {
1502  // Get number of bins
1503  int netrue = m_edisp.axis_bins(m_inx_etrue);
1504  int nmigra = m_edisp.axis_bins(m_inx_migra);
1505  int ntheta = m_edisp.axis_bins(m_inx_theta);
1506 
1507  // Initialise maximum energy dispersion
1508  m_max_edisp.assign(netrue*ntheta, 0.0);
1509 
1510  // Loop over offset angle
1511  for (int itheta = 0, index = 0; itheta < ntheta; ++itheta) {
1512 
1513  // Loop over true photon energy
1514  for (int ietrue = 0; ietrue < netrue; ++ietrue, ++index) {
1515 
1516  // Get true photon energy
1517  GEnergy etrue = this->etrue(ietrue);
1518 
1519  // Get offset and stride
1520  int offset = table_index(ietrue, 0, itheta);
1521  int stride = table_stride(m_inx_migra);
1522 
1523  // Find maximum
1524  double maximum = 0.0;
1525  for (int imigra = 0, i = offset; imigra < nmigra; ++imigra, i += stride) {
1526  double value = m_edisp(m_inx_matrix, i);
1527  if (value > maximum) {
1528  maximum = value;
1529  }
1530  }
1531 
1532  // Set maximum
1533  m_max_edisp[index] = maximum / etrue.MeV();
1534 
1535  } // endfor: looped over true photon energy
1536 
1537  } // endfor: looped over offset angle
1538 
1539  // Return
1540  return;
1541 }
1542 
1543 
1544 /***********************************************************************//**
1545  * @brief Get maximum energy dispersion value
1546  *
1547  * @param[in] etrue True photon energy.
1548  * @param[in] theta Offset angle (radians).
1549  * @return Maximum energy dispersion value.
1550  *
1551  * For a given true energy @p etrue and offset angle @p theta, returns the
1552  * maximum energy dispersion value that may occur. This method makes use
1553  * of the internal array pre-computed by set_max_edisp().
1554  *
1555  * If set_max_edisp() was not called, the method returns 0.
1556  ***************************************************************************/
1558  const double& theta) const
1559 {
1560  // Initialise maximum energy dispersion
1561  double max_edisp = 0.0;
1562 
1563  // Continue only if maximum energy dispersion array is set
1564  if (m_max_edisp.size() > 0) {
1565 
1566  // Get number of true energy bins
1567  int netrue = m_edisp.axis_bins(m_inx_etrue);
1568 
1569  // Get references to axis nodes
1570  const GNodeArray& etrue_nodes = m_edisp.axis_nodes(m_inx_etrue);
1571  const GNodeArray& theta_nodes = m_edisp.axis_nodes(m_inx_theta);
1572 
1573  // Set values
1574  etrue_nodes.set_value(etrue.log10TeV());
1575  theta_nodes.set_value(theta);
1576 
1577  // Get indices that bound the true energy and offset angle
1578  int inx_left = theta_nodes.inx_left() * netrue;
1579  int inx_right = theta_nodes.inx_right() * netrue;
1580  int index_1 = inx_left + etrue_nodes.inx_left();
1581  int index_2 = inx_left + etrue_nodes.inx_right();
1582  int index_3 = inx_right + etrue_nodes.inx_left();
1583  int index_4 = inx_right + etrue_nodes.inx_right();
1584 
1585  // Get maximum energy dispersion
1586  max_edisp = m_max_edisp[index_1];
1587  double max_edisp_2 = m_max_edisp[index_2];
1588  double max_edisp_3 = m_max_edisp[index_3];
1589  double max_edisp_4 = m_max_edisp[index_4];
1590  if (max_edisp_2 > max_edisp) {
1591  max_edisp = max_edisp_2;
1592  }
1593  if (max_edisp_3 > max_edisp) {
1594  max_edisp = max_edisp_3;
1595  }
1596  if (max_edisp_4 > max_edisp) {
1597  max_edisp = max_edisp_4;
1598  }
1599 
1600  } // endif: Maximum energy dispersion array was set
1601 
1602  // Return maximum energy dispersion
1603  return max_edisp;
1604 }
1605 
1606 
1607 /***********************************************************************//**
1608  * @brief Normalize energy dispersion table
1609  *
1610  * Normalize the energy dispersion table using
1611  *
1612  * \f[
1613  * \int_{E_{\rm reco}^{\rm min}}^{E_{\rm reco}^{\rm max}}
1614  * E_{\rm disp}(E_{\rm reco} | E_{\rm true}, \theta) \, dE_{\rm reco} = 1
1615  * \f]
1616  *
1617  * where
1618  * \f$E_{\rm reco}_{\rm min}\f$ and \f$E_{\rm reco}_{\rm max}\f$ is the
1619  * minimum and maximum migration value for a given true energy as returned
1620  * by the ebounds_obs() method, and
1621  * \f$E_{\rm disp}(E_{\rm true}, E_{\rm reco}, \theta)\f$ is the energy
1622  * dispersion returned by the operator(), given in units of \f$MeV^{-1}\f$.
1623  *
1624  * The normalisation is performed for each \f$E_{\rm true}\f$ and
1625  * \f$\theta\f$ bin using a Romberg integration method. Since the
1626  * normalisation may affect the energy boundaries
1627  * \f$E_{\rm reco}^{\rm min}\f$ and \f$E_{\rm reco}^{\rm max}\f$, two
1628  * normalisation passes are performed to stabilize the result.
1629  ***************************************************************************/
1631 {
1632  // Get axes dimensions
1633  int etrue_size = m_edisp.axis_bins(m_inx_etrue);
1634  int migra_size = m_edisp.axis_bins(m_inx_migra);
1635  int theta_size = m_edisp.axis_bins(m_inx_theta);
1636 
1637  // Now normalize the migration matrix vectors. We do this by integrating
1638  // for each true energy and offset angle the measured energy. We perform
1639  // here two passes as the ebounds_obs() method depends on the content of
1640  // the matrix, and by having two passes we are sure that we get a
1641  // consistent energy interval.
1642  for (int pass = 0; pass < 2; ++pass) {
1643 
1644  // Initialise sums
1645  std::vector<double> sums;
1646  sums.reserve(theta_size*etrue_size);
1647 
1648  // Loop over offset angle.
1649  for (int i_theta = 0; i_theta < theta_size; ++i_theta) {
1650 
1651  // Get offset angle (in radians)
1652  double theta = this->theta(i_theta);
1653 
1654  // Loop over true photon energy
1655  for (int i_etrue = 0; i_etrue < etrue_size; ++i_etrue) {
1656 
1657  // Get true photon energy
1658  GEnergy etrue = this->etrue(i_etrue);
1659 
1660  // Initialise integration
1661  double sum = 0.0;
1662 
1663  // Get integration boundaries
1664  GEbounds ebounds = ereco_bounds(etrue, theta);
1665 
1666  // Loop over all energy intervals
1667  for (int i = 0; i < ebounds.size(); ++i) {
1668 
1669  // Get energy boundaries in log10 of energy in MeV
1670  double emin = ebounds.emin(i).log10MeV();
1671  double emax = ebounds.emax(i).log10MeV();
1672 
1673  // Setup integration function
1674  edisp_ereco_kern integrand(this, etrue, theta);
1675  GIntegral integral(&integrand);
1676 
1677  // Set integration precision
1678  integral.eps(1.0e-6);
1679 
1680  // Do Romberg integration
1681  sum += integral.romberg(emin, emax);
1682 
1683  } // endfor: looped over energy intervals
1684 
1685  // Store sum
1686  sums.push_back(sum);
1687 
1688  } // endfor: looped over Etrue
1689 
1690  } // endfor: looped over theta
1691 
1692  // Normalize vectors of migration matrix for all etrue and theta.
1693  for (int i_theta = 0, inx = 0; i_theta < theta_size; ++i_theta) {
1694  for (int i_etrue = 0; i_etrue < etrue_size; ++i_etrue, ++inx) {
1695  double sum = sums[inx];
1696  if (sum > 0.0) {
1697  int offset = table_index(i_etrue, 0, i_theta);
1698  int stride = table_stride(m_inx_migra);
1699  for (int k = 0, i = offset; k < migra_size; ++k, i += stride) {
1700  m_edisp(m_inx_matrix,i) /= sum;
1701  }
1702  }
1703  }
1704  }
1705 
1706  } // endfor: made 2 passes
1707 
1708  // Return
1709  return;
1710 }
1711 
1712 
1713 /***********************************************************************//**
1714  * @brief Return index of response table element
1715  *
1716  * @param[in] ietrue True energy index.
1717  * @param[in] imigra Migration index.
1718  * @param[in] itheta Offset index.
1719  * @return Index of response table element.
1720  ***************************************************************************/
1721 int GCTAEdisp2D::table_index(const int& ietrue,
1722  const int& imigra,
1723  const int& itheta) const
1724 {
1725  // Set index vector
1726  int inx[3];
1727  inx[m_inx_etrue] = ietrue;
1728  inx[m_inx_migra] = imigra;
1729  inx[m_inx_theta] = itheta;
1730 
1731  // Compute index
1732  int index = inx[0] + (inx[1] + inx[2] * m_edisp.axis_bins(1)) *
1733  m_edisp.axis_bins(0);
1734 
1735  // Return index
1736  return index;
1737 }
1738 
1739 
1740 /***********************************************************************//**
1741  * @brief Return stride of response table axis
1742  *
1743  * @param[in] axis Response table axis.
1744  * @return Stride of response table axis.
1745  ***************************************************************************/
1746 int GCTAEdisp2D::table_stride(const int& axis) const
1747 {
1748  // Initialise stride
1749  int stride = 1;
1750 
1751  // Multiply in higher
1752  if (axis > 0) {
1753  stride *= m_edisp.axis_bins(0);
1754  }
1755  if (axis > 1) {
1756  stride *= m_edisp.axis_bins(1);
1757  }
1758 
1759  // Return stride
1760  return stride;
1761 }
1762 
1763 
1764 /***********************************************************************//**
1765  * @brief Return bi-linearly interpolate table value for given migration bin
1766  *
1767  * @param[in] base_ll Base index for left true energy and left offset angle.
1768  * @param[in] base_lr Base index for left true energy and right offset angle.
1769  * @param[in] base_rl Base index for right true energy and left offset angle.
1770  * @param[in] base_rr Base index for right true energy and right offset angle.
1771  * @param[in] wgt_el Weighting for left true energy.
1772  * @param[in] wgt_er Weighting for right true energy.
1773  * @param[in] wgt_tl Weighting for left offset angle.
1774  * @param[in] wgt_tr Weighting for right offset angle.
1775  * @param[in] offset Offset of migration bin with respect to base indices.
1776  * @return Bi-linearly interpolate table value.
1777  ***************************************************************************/
1778 double GCTAEdisp2D::table_value(const int& base_ll,
1779  const int& base_lr,
1780  const int& base_rl,
1781  const int& base_rr,
1782  const double& wgt_el,
1783  const double& wgt_er,
1784  const double& wgt_tl,
1785  const double& wgt_tr,
1786  const int& offset) const
1787 {
1788  // Compute table element indices
1789  int inx_ll = base_ll + offset;
1790  int inx_lr = base_lr + offset;
1791  int inx_rl = base_rl + offset;
1792  int inx_rr = base_rr + offset;
1793 
1794  // Get
1795  double value = wgt_el * (m_edisp(m_inx_matrix, inx_ll) * wgt_tl +
1796  m_edisp(m_inx_matrix, inx_lr) * wgt_tr) +
1797  wgt_er * (m_edisp(m_inx_matrix, inx_rl) * wgt_tl +
1798  m_edisp(m_inx_matrix, inx_rr) * wgt_tr);
1799 
1800  // Return value
1801  return value;
1802 }
1803 
1804 
1805 /***********************************************************************//**
1806  * @brief Return bi-linearly interpolate table value for given migration value
1807  *
1808  * @param[in] base_ll Base index for left true energy and left offset angle.
1809  * @param[in] base_lr Base index for left true energy and right offset angle.
1810  * @param[in] base_rl Base index for right true energy and left offset angle.
1811  * @param[in] base_rr Base index for right true energy and right offset angle.
1812  * @param[in] wgt_el Weighting for left true energy.
1813  * @param[in] wgt_er Weighting for right true energy.
1814  * @param[in] wgt_tl Weighting for left offset angle.
1815  * @param[in] wgt_tr Weighting for right offset angle.
1816  * @param[in] migra Migration value.
1817  * @return Bi-linearly interpolate table value.
1818  ***************************************************************************/
1819 double GCTAEdisp2D::table_value(const int& base_ll,
1820  const int& base_lr,
1821  const int& base_rl,
1822  const int& base_rr,
1823  const double& wgt_el,
1824  const double& wgt_er,
1825  const double& wgt_tl,
1826  const double& wgt_tr,
1827  const double& migra) const
1828 {
1829  // Retrieve references to migration node array
1830  const GNodeArray& migra_nodes = m_edisp.axis_nodes(m_inx_migra);
1831 
1832  // Set migration value for node array interpolation
1833  migra_nodes.set_value(migra);
1834 
1835  // Get migration stride
1836  int stride = table_stride(m_inx_migra);
1837 
1838  // Get left value
1839  double value_left = table_value(base_ll, base_lr, base_rl, base_rr,
1840  wgt_el, wgt_er, wgt_tl, wgt_tr,
1841  migra_nodes.inx_left() * stride);
1842 
1843  // Get right value
1844  double value_right = table_value(base_ll, base_lr, base_rl, base_rr,
1845  wgt_el, wgt_er, wgt_tl, wgt_tr,
1846  migra_nodes.inx_right() * stride);
1847 
1848  // Interpolate result
1849  double value = migra_nodes.wgt_left() * value_left +
1850  migra_nodes.wgt_right() * value_right;
1851 
1852  // Make sure that interpolation result is not negative
1853  if (value < 0.0) {
1854  value = 0.0;
1855  }
1856 
1857  // Return
1858  return value;
1859 }
1860 
1861 
1862 /***********************************************************************//**
1863  * @brief Smoothed energy dispersion table
1864  *
1865  * This method implements the kludge for CTA that reduces the statistical
1866  * noise in the energy dispersion matrix.
1867  ***************************************************************************/
1869 {
1870  // Get axes dimensions
1871  int netrue = m_edisp.axis_bins(m_inx_etrue);
1872  int nmigra = m_edisp.axis_bins(m_inx_migra);
1873  int ntheta = m_edisp.axis_bins(m_inx_theta);
1874  int npix = netrue * nmigra;
1875 
1876  // Loop over all offset angles
1877  for (int itheta = 0; itheta < ntheta; ++itheta) {
1878 
1879  // Compute means, rms and total as function of etrue
1880  GNdarray mean(netrue);
1881  GNdarray rms(netrue);
1882  GNdarray total(netrue);
1883  get_moments(itheta, &mean, &rms, &total);
1884 
1885  // Smooth all three
1886  mean = smooth_array(mean, 30, 30, 0.5);
1887  rms = smooth_array(rms, 30, 30, 0.03);
1888  total = smooth_array(total, 30, 30, 0.0);
1889 
1890  // Make sure that total is not negative
1891  for (int ietrue = 0; ietrue < netrue; ++ietrue) {
1892  if (total(ietrue) < 0.0) {
1893  total(ietrue) = 0.0;
1894  }
1895  }
1896 
1897  // Replace matrix by Gaussians
1898  for (int ietrue = 0; ietrue < netrue; ++ietrue) {
1899 
1900  // Continue only if there is information
1901  if (total(ietrue) > 0.0) {
1902 
1903  // Get Gaussian
1904  GNdarray gauss = gaussian_array(mean(ietrue), rms(ietrue), total(ietrue));
1905 
1906  // Compute base index
1907  int ibase = itheta * npix + ietrue;
1908  if ((mean(ietrue) > 0.0) && (rms(ietrue))) {
1909  for (int imigra = 0, i = ibase; imigra < nmigra; ++imigra, i += netrue) {
1910  m_edisp(m_inx_matrix, i) = gauss(imigra);
1911  }
1912  }
1913  else {
1914  for (int imigra = 0, i = ibase; imigra < nmigra; ++imigra, i += netrue) {
1915  m_edisp(m_inx_matrix, i) = 0.0;
1916  }
1917  }
1918 
1919  } // endif: there was information
1920 
1921  // ... otherwise reset the matrix to zero
1922  else {
1923  int ibase = itheta * npix + ietrue;
1924  for (int imigra = 0, i = ibase; imigra < nmigra; ++imigra, i += netrue) {
1925  m_edisp(m_inx_matrix, i) = 0.0;
1926  }
1927  }
1928 
1929  } // endfor: looped over true energies
1930 
1931  } // endfor: looped over offset angles
1932 
1933  // Return
1934  return;
1935 }
1936 
1937 
1938 /***********************************************************************//**
1939  * @brief Smoothed array
1940  *
1941  * @param[in] array Array.
1942  * @param[in] nstart Number of nodes used for regression at first array value.
1943  * @param[in] nstop Number of nodes used for regression at last array value.
1944  * @param[in] limit Limit for array element exclusion.
1945  * @return Smoothed array.
1946  *
1947  * Returns a smoothed array that is computed by locally adjusting a linear
1948  * regression law to the array elements.
1949  ***************************************************************************/
1951  const int& nstart,
1952  const int& nstop,
1953  const double& limit) const
1954 {
1955  // Initialise smoothed array
1956  GNdarray smoothed_array(array.size());
1957 
1958  // Compute node step size
1959  double nstep = double(nstop - nstart)/double(array.size());
1960 
1961  // Loop over all array elements and computed the smoothed value
1962  for (int i = 0; i < array.size(); ++i) {
1963  int nodes = nstart + int(i*nstep);
1964  smoothed_array(i) = smoothed_array_value(i, array, nodes, limit);
1965  }
1966 
1967  // Return smoothed array
1968  return smoothed_array;
1969 }
1970 
1971 
1972 /***********************************************************************//**
1973  * @brief Get smoothed array value
1974  *
1975  * @param[in] inx Index of array element.
1976  * @param[in] array Array.
1977  * @param[in] nodes Number of nodes used for regression.
1978  * @param[in] limit Limit for array element exclusion.
1979  * @return Smoothed array value.
1980  *
1981  * Returns a smoothed array value that is derived from adjusting a linear
1982  * law using regression to @p nodes adjacent array elements. Array elements
1983  * with values below @p limit are excluded from the regression.
1984  ***************************************************************************/
1986  const GNdarray& array,
1987  const int& nodes,
1988  const double& limit) const
1989 {
1990  // Initialise variables
1991  double mean_x = 0.0;
1992  double mean_y = 0.0;
1993  int ileft = inx - 1;
1994  int iright = inx + 1;
1995  int nleft = 0;
1996  int nright = 0;
1997 
1998  // Allocate vector of elements used for regression
1999  std::vector<double> x_vals;
2000  std::vector<double> y_vals;
2001 
2002  // Add nodes on the left of the element of interest
2003  while (nleft < nodes) {
2004  if (ileft < 0) {
2005  break;
2006  }
2007  if (array(ileft) > limit) {
2008  double x = double(ileft);
2009  double y = array(ileft);
2010  x_vals.push_back(x);
2011  y_vals.push_back(y);
2012  mean_x += x;
2013  mean_y += y;
2014  nleft++;
2015  }
2016  ileft--;
2017  }
2018 
2019  // Add remaining nodes on the right of the element of interest
2020  while (nright < 2*nodes-nleft) {
2021  if (iright >= array.size()) {
2022  break;
2023  }
2024  if (array(iright) > limit) {
2025  double x = double(iright);
2026  double y = array(iright);
2027  x_vals.push_back(x);
2028  y_vals.push_back(y);
2029  mean_x += x;
2030  mean_y += y;
2031  nright++;
2032  }
2033  iright++;
2034  }
2035 
2036  // Add remaining nodes on the left if all right nodes were not exhausted
2037  while (nleft < 2*nodes-nright) {
2038  if (ileft < 0) {
2039  break;
2040  }
2041  if (array(ileft) > limit) {
2042  double x = double(ileft);
2043  double y = array(ileft);
2044  x_vals.push_back(x);
2045  y_vals.push_back(y);
2046  mean_x += x;
2047  mean_y += y;
2048  nleft++;
2049  }
2050  ileft--;
2051  }
2052 
2053  // Compute mean x and y values
2054  double total = double(nleft+nright);
2055  if (total > 0.0) {
2056  mean_x /= total;
2057  mean_y /= total;
2058  }
2059 
2060  // Compute regression line slope
2061  double beta_nom = 0.0;
2062  double beta_denom = 0.0;
2063  for (int i = 0; i < x_vals.size(); ++i) {
2064  double x = x_vals[i];
2065  double y = y_vals[i];
2066  beta_nom += (x - mean_x) * (y - mean_y);
2067  beta_denom += (x - mean_x) * (x - mean_x);
2068  }
2069  double beta = beta_nom / beta_denom;
2070 
2071  // Compute regression line offset
2072  double alpha = mean_y - beta * mean_x;
2073 
2074  // Compute value
2075  double value = alpha + beta * double(inx);
2076 
2077  // Return value
2078  return value;
2079 }
2080 
2081 
2082 /***********************************************************************//**
2083  * @brief Compute moments
2084  *
2085  * @param[in] itheta Offset angle index.
2086  * @param[out] mean Pointer to mean array.
2087  * @param[out] rms Pointer to rms array.
2088  * @param[out] total Pointer to total array.
2089  *
2090  * Computes the mean and root mean square values as function of true energy
2091  * for a given offset angle.
2092  ***************************************************************************/
2093 void GCTAEdisp2D::get_moments(const int& itheta,
2094  GNdarray* mean,
2095  GNdarray* rms,
2096  GNdarray* total) const
2097 {
2098  // Get axes dimensions
2099  int netrue = m_edisp.axis_bins(m_inx_etrue);
2100  int nmigra = m_edisp.axis_bins(m_inx_migra);
2101  int npix = netrue * nmigra;
2102 
2103  // Loop over all true energies
2104  for (int ietrue = 0; ietrue < netrue; ++ietrue) {
2105 
2106  // Compute base index
2107  int ibase = itheta * npix + ietrue;
2108 
2109  // Extract array
2110  GNdarray array(nmigra);
2111  for (int imigra = 0, i = ibase; imigra < nmigra; ++imigra, i += netrue) {
2112  array(imigra) = m_edisp(m_inx_matrix, i);
2113  }
2114 
2115  // Store sum
2116  (*total)(ietrue) = sum(array);
2117 
2118  // Get Gaussian mean and rms
2119  double mean_value = 0.0;
2120  double rms_value = 0.0;
2121  get_mean_rms(array, &mean_value, &rms_value);
2122 
2123  // Store mean and rms
2124  (*mean)(ietrue) = mean_value;
2125  (*rms)(ietrue) = rms_value;
2126 
2127  } // endfor: looped over all true energies
2128 
2129  // Return
2130  return;
2131 }
2132 
2133 
2134 /***********************************************************************//**
2135  * @brief Compute mean and root mean square of migration array
2136  *
2137  * @param[in] array Energy dispersion array.
2138  * @param[out] mean Pointer to mean migration value.
2139  * @param[out] rms Pointer to root mean square value.
2140  *
2141  * Computes the mean and the root mean square of the migration array. If the
2142  * migration array is empty the mean and root mean square values are set to
2143  * zero.
2144  *
2145  * The method does not check whether the @p mean and @p rms pointers are
2146  * valid.
2147  ***************************************************************************/
2149  double* mean,
2150  double* rms) const
2151 {
2152  // Initialise mean and rms
2153  *mean = 0.0;
2154  *rms = 0.0;
2155 
2156  // Get reference to migration values
2157  const GNodeArray& migras = m_edisp.axis_nodes(m_inx_migra);
2158  int nmigra = migras.size();
2159 
2160  // Pre-compute values for mean and rms computation
2161  double sum = 0.0;
2162  for (int imigra = 0; imigra < nmigra; ++imigra) {
2163  double migra = migras[imigra];
2164  double weight = migra * array(imigra);
2165  *mean += migra * array(imigra);
2166  *rms += weight * migra;
2167  sum += array(imigra);
2168  }
2169 
2170  // If the array is not empty then compute now mean and standard deviation
2171  if (sum > 0.0) {
2172  *mean /= sum;
2173  *rms /= sum;
2174  *rms -= (*mean) * (*mean);
2175  *rms = std::sqrt(*rms);
2176  }
2177 
2178  // Return
2179  return;
2180 }
2181 
2182 
2183 /***********************************************************************//**
2184  * @brief Return Gaussian approximation of energy dispersion array
2185  *
2186  * @param[in] mean Gaussian mean.
2187  * @param[in] rms Gaussian rms.
2188  * @param[in] total Gaussian total.
2189  * @return Gaussian approximation of energy dispersion array.
2190  *
2191  * Returns a Gaussian approximation of the energy dispersion array by
2192  * computing the mean migration value and its root mean square and by using
2193  * these values as the centre and the width of a Gaussian function. The
2194  * Gaussian function is normalized so that the sum of the output array is
2195  * unity.
2196  ***************************************************************************/
2198  const double& rms,
2199  const double& total) const
2200 {
2201  // Initialise empty Gaussian array
2202  int nmigra = m_edisp.axis_bins(m_inx_migra);
2203  GNdarray gaussian(nmigra);
2204 
2205  // If the rms is valid then compute Gaussian
2206  if (rms > 0.0) {
2207 
2208  // Get reference to migration values
2209  const GNodeArray& migras = m_edisp.axis_nodes(m_inx_migra);
2210  int nmigra = migras.size();
2211 
2212  // Compute Gaussian
2213  double total_gauss = 0.0;
2214  for (int imigra = 0; imigra < nmigra; ++imigra) {
2215  double arg = (migras[imigra] - mean) / rms;
2216  double value = std::exp(-0.5*arg*arg);
2217  gaussian(imigra) = value;
2218  total_gauss += value;
2219  }
2220 
2221  // Normalise Gaussian
2222  if (total_gauss > 0.0) {
2223  gaussian *= total / total_gauss;
2224  }
2225 
2226  } // endif: computed Gaussian
2227 
2228  // Return Gaussian array
2229  return gaussian;
2230 }
2231 
2232 
2233 /***********************************************************************//**
2234  * @brief Integration kernel for GCTAEdisp2D::edisp_ereco_kern class
2235  *
2236  * @param[in] log10Ereco Log10 of reconstructed energy (\f$\log_{10}\f$ MeV).
2237  * @return Energy dispersion (\f$(\log_{10}\f$ MeV\f$)^{-1}\f$).
2238  *
2239  * This method implements the function
2240  *
2241  * \f[
2242  * E_{\rm disp}(\log_{10} E_{\rm reco} | E_{\rm true}, \theta) =
2243  * E_{\rm disp}(m | E_{\rm true}, \theta) \times
2244  * \frac{\log 10 \times E_{\rm reco}}{E_{\rm true}}
2245  * \f]
2246  *
2247  * which is the integration kernel needed for the
2248  * GCTAEdisp2D::edisp_ereco_kern class. The class is used by
2249  * GCTAEdisp2D::normalize_table() to integrate the energy dispersion
2250  * information using
2251  *
2252  * \f[
2253  * \int_{\log_{10} E_{\rm reco}^{\rm min}}^{\log_{10} E_{\rm reco}^{\rm max}}
2254  * E_{\rm disp}(\log_{10} E_{\rm reco} | E_{\rm true}, \theta) \,
2255  * d(\log_{10} E_{\rm reco})
2256  * \f]
2257  ***************************************************************************/
2258 double GCTAEdisp2D::edisp_ereco_kern::eval(const double& log10Ereco)
2259 {
2260  // Set reconstructued event energy
2261  GEnergy ereco;
2262  ereco.log10MeV(log10Ereco);
2263 
2264  // Get function value
2265  double value = m_parent->operator()(ereco, m_etrue, m_theta);
2266 
2267  // Normalize energy dispersion so that the units are per log10 MeV
2268  value *= gammalib::ln10 * ereco.MeV();
2269 
2270  // Compile option: Check for NaN
2271  #if defined(G_NAN_CHECK)
2272  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
2273  std::cout << "*** ERROR: GCTAEdisp2D::edisp_ereco_kern::eval";
2274  std::cout << "(log10Ereco=" << log10Ereco << "): ";
2275  std::cout << " NaN/Inf encountered";
2276  std::cout << " (value=" << value;
2277  std::cout << ")" << std::endl;
2278  }
2279  #endif
2280 
2281  // Return value
2282  return value;
2283 }
int m_inx_etrue
True energy index.
const GNodeArray & axis_nodes(const int &axis) const
Return axis nodes.
void normalize_table(void)
Normalize energy dispersion table.
int size(void) const
Return number of nodes in node array.
Definition: GNodeArray.hpp:188
CTA 2D energy dispersion class.
Definition: GCTAEdisp2D.hpp:90
std::string print(const GChatter &chatter=NORMAL) const
Print energy dispersion information.
GEnergy m_etrue
True photon energy.
GFilename filename(void) const
Return filename.
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition: GFits.cpp:472
Node array class.
Definition: GNodeArray.hpp:60
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:380
GEnergy etrue(const int &ietrue) const
Get true energy.
Random number generator class definition.
const GCTAEdisp2D * m_parent
Pointer to parent class.
double smoothed_array_value(const int &inx, const GNdarray &array, const int &nodes, const double &limit) const
Get smoothed array value.
double theta(const int &itheta) const
Get offset angle in radiaus.
int size(void) const
Return number of energy boundaries.
Definition: GEbounds.hpp:157
GCTAEdisp2D(void)
Void constructor.
Definition: GCTAEdisp2D.cpp:76
void set_max_edisp(void)
Set array of maximum energy dispersion values.
bool is_empty(void) const
Signal if filename is empty.
Definition: GFilename.hpp:160
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
double m_last_theta_etrue
bool m_ereco_bounds_computed
double TeV(void) const
Return energy in TeV.
Definition: GEnergy.cpp:348
double sum(const GVector &vector)
Computes vector sum.
Definition: GVector.cpp:901
std::string extname(const std::string &defaultname="") const
Return extension name.
Definition: GFilename.cpp:385
double log10TeV(void) const
Return log10 of energy in TeV.
Definition: GEnergy.cpp:457
double m_last_theta_ereco
Random number generator class.
Definition: GRan.hpp:44
double migra(const int &imigra) const
Get migration.
void copy_members(const GCTAEdisp2D &edisp)
Copy class members.
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
virtual ~GCTAEdisp2D(void)
Destructor.
GCTAResponseTable m_edisp
Edisp response table.
Gammalib tools definition.
void fetch(void) const
Fetch energy dispersion.
GIntegral class interface definition.
Definition: GIntegral.hpp:45
FITS file class.
Definition: GFits.hpp:63
void axis_radians(const int &axis)
Set nodes for a radians axis.
FITS file class interface definition.
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 set_value(const double &value) const
Set indices and weighting factors for interpolation.
Definition: GNodeArray.cpp:580
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.
std::string strip_whitespace(const std::string &arg)
Strip leading and trailing whitespace from string.
Definition: GTools.cpp:73
bool m_etrue_bounds_computed
void read(const GFitsTable &table)
Read response table from FITS table HDU.
const double ln10
Definition: GMath.hpp:46
double m_theta_max
Maximum theta (radians)
double log10MeV(void) const
Return log10 of energy in MeV.
Definition: GEnergy.cpp:423
bool m_fetched
Signals that Edisp has been fetched.
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:184
void write(GFitsBinTable &table) const
Write energy dispersion into FITS binary table.
const double & axis_lo(const int &axis, const int &bin) const
Return lower bin boundary for bin in axis.
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:167
Definition of support function used by CTA classes.
void compute_ereco_bounds(const double &theta=0.0, const double &phi=0.0, const double &zenith=0.0, const double &azimuth=0.0) const
Compute vector of reconstructed energy bounds.
void write(GFitsTable &table) const
Write response table into FITS table HDU.
int axis(const std::string &name) const
Determine index of an axis.
const double deg2rad
Definition: GMath.hpp:43
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1268
#define G_FETCH
Definition: GCTAEdisp2D.cpp:51
const double & wgt_right(void) const
Returns right node weight.
Definition: GNodeArray.hpp:275
void load(const GFilename &filename)
Load energy dispersion from FITS file.
Energy boundaries container class.
Definition: GEbounds.hpp:60
void remove(const int &extno)
Remove HDU from FITS file.
Definition: GFits.cpp:848
void set_table(void)
Set table.
std::string print(const GChatter &chatter=NORMAL) const
Print energy.
Definition: GEnergy.cpp:748
double eval(const double &log10Ereco)
Integration kernel for GCTAEdisp2D::edisp_ereco_kern class.
int table_stride(const int &axis) const
Return stride of response table axis.
void clear(void)
Clear response table.
Filename class.
Definition: GFilename.hpp:62
void smooth_table(void)
Smoothed energy dispersion table.
N-dimensional array class interface definition.
const GEnergy & emin(void) const
Return minimum energy of all intervals.
Definition: GEbounds.hpp:181
double m_logEsrc_min
Minimum logE (log10(E/TeV))
void get_mean_rms(const GNdarray &array, double *mean, double *rms) const
Compute mean and root mean square of migration array.
bool exists(void) const
Checks whether file exists.
Definition: GFilename.cpp:223
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
int m_inx_matrix
Matrix.
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
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.
int size(void) const
Return number of elements in array.
Definition: GNdarray.hpp:293
GChatter
Definition: GTypemaps.hpp:33
#define G_READ
Definition: GCTAEdisp2D.cpp:48
const int & inx_left(void) const
Returns left node index.
Definition: GNodeArray.hpp:231
CTA 2D energy dispersion class definition.
void save(const GFilename &filename, const bool &clobber=false) const
Save energy dispersion table into FITS file.
N-dimensional array class.
Definition: GNdarray.hpp:44
#define G_MC
Definition: GCTAEdisp2D.cpp:49
const std::string & extname(void) const
Return extension name.
Definition: GFitsHDU.hpp:162
void axis_log10(const int &axis)
Set nodes for a logarithmic (base 10) axis.
double get_max_edisp(const GEnergy &etrue, const double &theta) const
Get maximum energy dispersion value.
GNdarray smooth_array(const GNdarray &array, const int &nstart, const int &nstop, const double &limit) const
Smoothed array.
int table(const std::string &name) const
Determine index of table.
void compute_etrue_bounds(const double &theta=0.0, const double &phi=0.0, const double &zenith=0.0, const double &azimuth=0.0) const
Compute vector of true energy boundaries.
const GCTAResponseTable & table(void) const
Return response table.
std::string gadf_hduclas4(const GFits &fits, const std::string &hduclas4)
Return extension name for GADF response table of given HDU class 4.
double m_logEsrc_max
Maximum logE (log10(E/TeV))
const int & inx_right(void) const
Returns right node index.
Definition: GNodeArray.hpp:245
void free_members(void)
Delete class members.
Fast Fourier transformation class interface definition.
std::vector< GEbounds > m_ereco_bounds
GCTAEdisp & operator=(const GCTAEdisp &edisp)
Assignment operator.
Definition: GCTAEdisp.cpp:106
void clear(void)
Clear energy dispersion.
double m_theta_min
Minimum theta (radians)
std::vector< double > m_max_edisp
Maximum Edisp.
const std::string & telescope(void) const
Return telescope string.
GCTAEdisp2D * clone(void) const
Clone energy dispersion.
const std::string extname_cta_edisp2d
Definition: GCTAEdisp2D.hpp:46
void init_members(void)
Initialise class members.
void eps(const double &eps)
Set relative precision.
Definition: GIntegral.hpp:206
int m_inx_theta
Theta index.
GFilename m_filename
Name of Edisp response file.
std::string url(void) const
Return Uniform Resource Locator (URL)
Definition: GFilename.hpp:189
int axis_bins(const int &axis) const
Return number bins in an axis.
void clear(void)
Clear file name.
Definition: GFilename.cpp:188
double table_value(const int &base_ll, const int &base_lr, const int &base_rl, const int &base_rr, const double &wgt_el, const double &wgt_er, const double &wgt_tl, const double &wgt_tr, const int &offset) const
Return bi-linearly interpolate table value for given migration bin.
int table_index(const int &ietrue, const int &imigra, const int &itheta) const
Return index of response table element.
int m_inx_migra
Migration index.
GNdarray gaussian_array(const double &mean, const double &rms, const double &total) const
Return Gaussian approximation of energy dispersion array.
FITS binary table class.
bool has_axis(const std::string &name) const
Check whether an axis exists.
const double & axis_hi(const int &axis, const int &bin) const
Return upper bin boundary for bin in axis.
Exception handler interface definition.
FITS binary table class definition.
GVector exp(const GVector &vector)
Computes exponential of vector elements.
Definition: GVector.cpp:1142
CTA response table class.
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 .
void init_members(void)
Initialise class members.
Definition: GCTAEdisp.cpp:142
double m_migra_min
Minimum migration.
const GEnergy & emax(void) const
Return maximum energy of all intervals.
Definition: GEbounds.hpp:193
void get_moments(const int &itheta, GNdarray *mean, GNdarray *rms, GNdarray *total) const
Compute moments.
const double rad2deg
Definition: GMath.hpp:44
GEnergy m_last_etrue
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
Integration class interface definition.
GEnergy m_last_ereco
void close(void)
Close FITS file.
Definition: GFits.cpp:1314
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.
const int & axes(void) const
Return number of axes of the tables.
GCTAEdisp2D & operator=(const GCTAEdisp2D &edisp)
Assignment operator.
double m_theta
Offset angle.
void clear(void)
Clear instance.
Definition: GEnergy.cpp:261
void read(const GFitsTable &table)
Read energy dispersion from FITS table.
double m_migra_max
Maximum migration.
Filename class interface definition.
Mathematical function definitions.
GVector log10(const GVector &vector)
Computes base10 logarithm of vector elements.
Definition: GVector.cpp:1205
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
std::vector< GEbounds > m_etrue_bounds
void set_boundaries(void)
Set energy dispersion boundaries.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:413
FITS table abstract base class interface definition.