GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAEdisp2D.hpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAEdisp2D.hpp - CTA 2D energy dispersion class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2015-2018 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.hpp
23  * @brief CTA 2D energy dispersion class definition
24  * @author Florent Forest
25  */
26 
27 #ifndef GCTAEDISP2D_HPP
28 #define GCTAEDISP2D_HPP
29 
30 /* __ Includes ___________________________________________________________ */
31 #include <string>
32 #include "GFilename.hpp"
33 #include "GFunction.hpp"
34 #include "GEnergy.hpp"
35 #include "GCTAEdisp.hpp"
36 #include "GCTAResponseTable.hpp"
37 
38 /* __ Forward declarations _______________________________________________ */
39 class GFft;
40 class GRan;
41 class GFits;
42 class GFitsBinTable;
43 
44 /* __ Constants __________________________________________________________ */
45 namespace gammalib {
46  const std::string extname_cta_edisp2d = "ENERGY DISPERSION";
47 }
48 
49 
50 /***********************************************************************//**
51  * @class GCTAEdisp2D
52  *
53  * @brief CTA 2D energy dispersion class
54  *
55  * This class implements the energy dispersion for the CTA 2D response. The
56  * CTA 2D energy dispersion is in fact a 3-dimensional function, where the
57  * energy dispersion
58  *
59  * \f[
60  * E_{\rm disp}(E_{\rm reco} | E_{\rm true}, \theta)
61  * \f]
62  *
63  * is given as function of true energy \f$E_{\rm true}\f$, reconstructed
64  * energy \f$E_{\rm reco}\f$ and offset angle \f$\theta\f$. The
65  * GCTAEdisp2D::operator() returns the value of the function
66  * \f$E_{\rm disp}(E_{\rm reco} | E_{\rm true}, \theta)\f$ in units of
67  * MeV\f$^{-1}\f$ with
68  *
69  * \f[
70  * \int_{E_{\rm reco}^{\rm min}}^{E_{\rm reco}^{\rm max}}
71  * E_{\rm disp}(E_{\rm reco} | E_{\rm true}, \theta) \, dE_{\rm reco} = 1
72  * \f]
73  *
74  * The energy dispersion is stored internally as a 3-dimensional response
75  * table, spanned by true energy \f$E_{\rm true}\f$, migration
76  * \f$m=E_{\rm reco}/E_{\rm true}\f$ and offset angle \f$\theta\f$, with
77  *
78  * \f[
79  * E_{\rm disp}(E_{\rm reco} | E_{\rm true}, \theta) =
80  * \frac{E_{\rm disp}(m | E_{\rm true}, \theta)}{E_{\rm true}}
81  * \f]
82  *
83  * and
84  *
85  * \f[
86  * \int_{m_{\rm min}}^{m_{\rm max}}
87  * E_{\rm disp}(m | E_{\rm true}, \theta) \, dm = 1
88  * \f]
89  ***************************************************************************/
90 class GCTAEdisp2D : public GCTAEdisp {
91 
92 public:
93  // Constructors and destructors
94  GCTAEdisp2D(void);
95  explicit GCTAEdisp2D(const GFilename& filename);
96  GCTAEdisp2D(const GCTAEdisp2D& edisp);
97  virtual ~GCTAEdisp2D(void);
98 
99  // Operators
100  GCTAEdisp2D& operator=(const GCTAEdisp2D& edisp);
101  double operator()(const GEnergy& ereco,
102  const GEnergy& etrue,
103  const double& theta = 0.0,
104  const double& phi = 0.0,
105  const double& zenith = 0.0,
106  const double& azimuth = 0.0) const;
107 
108  // Implemented methods
109  void clear(void);
110  GCTAEdisp2D* clone(void) const;
111  std::string classname(void) const;
112  void load(const GFilename& filename);
113  GFilename filename(void) const;
114  GEnergy mc(GRan& ran,
115  const GEnergy& etrue,
116  const double& theta = 0.0,
117  const double& phi = 0.0,
118  const double& zenith = 0.0,
119  const double& azimuth = 0.0) const;
120  GEbounds ereco_bounds(const GEnergy& etrue,
121  const double& theta = 0.0,
122  const double& phi = 0.0,
123  const double& zenith = 0.0,
124  const double& azimuth = 0.0) const;
125  GEbounds etrue_bounds(const GEnergy& ereco,
126  const double& theta = 0.0,
127  const double& phi = 0.0,
128  const double& zenith = 0.0,
129  const double& azimuth = 0.0) const;
130  double prob_erecobin(const GEnergy& ereco_min,
131  const GEnergy& ereco_max,
132  const GEnergy& etrue,
133  const double& theta) const;
134  std::string print(const GChatter& chatter = NORMAL) const;
135 
136  // Methods
137  void fetch(void) const;
138  const GCTAResponseTable& table(void) const;
139  void table(const GCTAResponseTable& table);
140  void read(const GFitsTable& table);
141  void write(GFitsBinTable& table) const;
142  void save(const GFilename& filename,
143  const bool& clobber = false) const;
144 
145 private:
146  // Methods
147  void init_members(void);
148  void copy_members(const GCTAEdisp2D& edisp);
149  void free_members(void);
150  GEnergy etrue(const int& ietrue) const;
151  double migra(const int& imigra) const;
152  double theta(const int& itheta) const;
153  void compute_ereco_bounds(const double& theta = 0.0,
154  const double& phi = 0.0,
155  const double& zenith = 0.0,
156  const double& azimuth = 0.0) const;
157  void compute_etrue_bounds(const double& theta = 0.0,
158  const double& phi = 0.0,
159  const double& zenith = 0.0,
160  const double& azimuth = 0.0) const;
161  void set_table(void);
162  void set_boundaries(void);
163  void set_max_edisp(void);
164  double get_max_edisp(const GEnergy& etrue, const double& theta) const;
165  void normalize_table(void);
166  int table_index(const int& ietrue,
167  const int& imigra,
168  const int& itheta) const;
169  int table_stride(const int& axis) const;
170  double table_value(const int& base_ll,
171  const int& base_lr,
172  const int& base_rl,
173  const int& base_rr,
174  const double& wgt_el,
175  const double& wgt_er,
176  const double& wgt_tl,
177  const double& wgt_tr,
178  const int& offset) const;
179  double table_value(const int& base_ll,
180  const int& base_lr,
181  const int& base_rl,
182  const int& base_rr,
183  const double& wgt_el,
184  const double& wgt_er,
185  const double& wgt_tl,
186  const double& wgt_tr,
187  const double& migra) const;
188 
189  // Kludge
190  void smooth_table(void);
191  GNdarray smooth_array(const GNdarray& array,
192  const int& nstart,
193  const int& nstop,
194  const double& limit) const;
195  double smoothed_array_value(const int& inx,
196  const GNdarray& array,
197  const int& nodes,
198  const double& limit) const;
199  void get_moments(const int& itheta,
200  GNdarray* mean,
201  GNdarray* rms,
202  GNdarray* total) const;
203  void get_mean_rms(const GNdarray& array, double* mean, double* rms) const;
204  GNdarray gaussian_array(const double& mean,
205  const double& rms,
206  const double& total) const;
207 
208  // Protected classes
209  class edisp_ereco_kern : public GFunction {
210  public:
212  const GEnergy& etrue,
213  const double& theta) :
214  m_parent(parent),
215  m_etrue(etrue),
216  m_theta(theta) { }
217  double eval(const double& log10Ereco);
218  protected:
219  const GCTAEdisp2D* m_parent; //!< Pointer to parent class
220  GEnergy m_etrue; //!< True photon energy
221  double m_theta; //!< Offset angle
222  };
223 
224  // Members
225  mutable GFilename m_filename; //!< Name of Edisp response file
226  GCTAResponseTable m_edisp; //!< Edisp response table
227  mutable bool m_fetched; //!< Signals that Edisp has been fetched
228  int m_inx_etrue; //!< True energy index
229  int m_inx_migra; //!< Migration index
230  int m_inx_theta; //!< Theta index
231  int m_inx_matrix; //!< Matrix
232  double m_logEsrc_min; //!< Minimum logE (log10(E/TeV))
233  double m_logEsrc_max; //!< Maximum logE (log10(E/TeV))
234  double m_migra_min; //!< Minimum migration
235  double m_migra_max; //!< Maximum migration
236  double m_theta_min; //!< Minimum theta (radians)
237  double m_theta_max; //!< Maximum theta (radians)
238  std::vector<double> m_max_edisp; //!< Maximum Edisp
239 
240  // Computation cache
243  mutable double m_last_theta_ereco;
244  mutable double m_last_theta_etrue;
247  mutable int m_index_ereco;
248  mutable int m_index_etrue;
249  mutable std::vector<GEbounds> m_ereco_bounds;
250  mutable std::vector<GEbounds> m_etrue_bounds;
251 };
252 
253 
254 /***********************************************************************//**
255  * @brief Return class name
256  *
257  * @return String containing the class name ("GCTAEdisp2D").
258  ***************************************************************************/
259 inline
260 std::string GCTAEdisp2D::classname(void) const
261 {
262  return ("GCTAEdisp2D");
263 }
264 
265 
266 /***********************************************************************//**
267  * @brief Return filename
268  *
269  * @return Name of FITS file from which energy dispersion was loaded.
270  ***************************************************************************/
271 inline
273 {
274  // Return filename
275  return (m_filename);
276 }
277 
278 
279 /***********************************************************************//**
280  * @brief Return response table
281  *
282  * @return Reference to response table.
283  ***************************************************************************/
284 inline
286 {
287  fetch();
288  return (m_edisp);
289 }
290 
291 #endif /* GCTAEDISP2D_HPP */
int m_inx_etrue
True energy index.
void normalize_table(void)
Normalize energy dispersion table.
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.
GEnergy etrue(const int &ietrue) const
Get true energy.
Energy value 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.
std::string classname(void) const
Return class name.
double theta(const int &itheta) const
Get offset angle in radiaus.
GCTAEdisp2D(void)
Void constructor.
Definition: GCTAEdisp2D.cpp:76
void set_max_edisp(void)
Set array of maximum energy dispersion values.
Abstract base class for the CTA energy dispersion.
Definition: GCTAEdisp.hpp:49
double m_last_theta_etrue
edisp_ereco_kern(const GCTAEdisp2D *parent, const GEnergy &etrue, const double &theta)
bool m_ereco_bounds_computed
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.
virtual ~GCTAEdisp2D(void)
Destructor.
GCTAResponseTable m_edisp
Edisp response table.
void fetch(void) const
Fetch energy dispersion.
FITS file class.
Definition: GFits.hpp:63
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.
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.
bool m_etrue_bounds_computed
double m_theta_max
Maximum theta (radians)
bool m_fetched
Signals that Edisp has been fetched.
void write(GFitsBinTable &table) const
Write energy dispersion into FITS binary table.
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.
CTA response table class definition.
void load(const GFilename &filename)
Load energy dispersion from FITS file.
Energy boundaries container class.
Definition: GEbounds.hpp:60
Single parameter function abstract base class definition.
void set_table(void)
Set table.
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.
Filename class.
Definition: GFilename.hpp:62
void smooth_table(void)
Smoothed energy dispersion table.
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.
Fast Fourier Transformation class.
Definition: GFft.hpp:57
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.
GChatter
Definition: GTypemaps.hpp:33
Abstract CTA energy dispersion base 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
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.
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.
double m_logEsrc_max
Maximum logE (log10(E/TeV))
void free_members(void)
Delete class members.
std::vector< GEbounds > m_ereco_bounds
void clear(void)
Clear energy dispersion.
double m_theta_min
Minimum theta (radians)
std::vector< double > m_max_edisp
Maximum Edisp.
GCTAEdisp2D * clone(void) const
Clone energy dispersion.
const std::string extname_cta_edisp2d
Definition: GCTAEdisp2D.hpp:46
void init_members(void)
Initialise class members.
int m_inx_theta
Theta index.
GFilename m_filename
Name of Edisp response file.
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.
Single parameter function abstract base class.
Definition: GFunction.hpp:44
FITS binary table class.
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 .
double m_migra_min
Minimum migration.
void get_moments(const int &itheta, GNdarray *mean, GNdarray *rms, GNdarray *total) const
Compute moments.
GEnergy m_last_etrue
GEnergy m_last_ereco
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.
GCTAEdisp2D & operator=(const GCTAEdisp2D &edisp)
Assignment operator.
double m_theta
Offset angle.
void read(const GFitsTable &table)
Read energy dispersion from FITS table.
double m_migra_max
Maximum migration.
Filename class interface definition.
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.