GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCOMIaq.hpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCOMIaq.hpp - COMPTEL instrument response representation *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2017 by Juergen Knoedlseder *
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 GCOMIaq.hpp
23  * @brief COMPTEL instrument response representation class interface definition
24  * @author Juergen Knoedlseder
25  */
26 
27 #ifndef GCOMIAQ_HPP
28 #define GCOMIAQ_HPP
29 
30 /* __ Includes ___________________________________________________________ */
31 #include <string>
32 #include "GMath.hpp"
33 #include "GEbounds.hpp"
34 #include "GFunction.hpp"
35 #include "GFitsImageFloat.hpp"
36 #include "GCOMD1Response.hpp"
37 #include "GCOMD2Response.hpp"
38 #include "GCOMInstChars.hpp"
39 
40 /* __ Type definitions ___________________________________________________ */
41 
42 /* __ Forward declarations _______________________________________________ */
43 class GFilename;
44 class GModelSpectral;
45 
46 
47 /***********************************************************************//**
48  * @class GCOMIaq
49  *
50  * @brief Interface for the COMPTEL instrument response representation class
51  ***************************************************************************/
52 class GCOMIaq : public GBase {
53 
54 public:
55  // Constructors and destructors
56  GCOMIaq(void);
57  GCOMIaq(const GCOMIaq& iaq);
58  GCOMIaq(const double& phigeo_max, const double& phigeo_bin_size,
59  const double& phibar_max, const double& phibar_bin_size);
60  ~GCOMIaq(void);
61 
62  // Operators
63  GCOMIaq& operator=(const GCOMIaq & iaq);
64 
65  // Methods
66  void clear(void);
67  GCOMIaq* clone(void) const;
68  std::string classname(void) const;
69  void set(const GEnergy& energy, const GEbounds& ebounds);
70  void set(const GModelSpectral& spectrum, const GEbounds& ebounds);
71  void save(const GFilename& filename, const bool& clobber = false) const;
72  std::string print(const GChatter& chatter = NORMAL) const;
73 
74 private:
75  // Private methods
76  void init_members(void);
77  void copy_members(const GCOMIaq& iaq);
78  void free_members(void);
79  void init_response(void);
80  void remove_cards(void);
81 
82  // RESPSI methods
83  void compton_kinematics(const double& energy);
84  double compute_iaq_bin(const double& etrue1,
85  const double& etrue2,
86  const double& phibar);
87  void klein_nishina(const double& energy);
88  double klein_nishina_bin(const double& energy,
89  const double& phigeo);
90  double klein_nishina_integral(const double& v,
91  const double& a);
92  void weight_iaq(const double& energy);
93  void location_smearing(const double& zenith);
94  std::vector<double> location_spread(const double& zenith) const;
95 
96  // Reponse integration kernel
97  class response_kernel : public GFunction {
98  public:
99  response_kernel(const GCOMD1Response& response_d1,
100  const GCOMD2Response& response_d2,
101  const double& etrue1,
102  const double& etrue2,
103  const double& phibar,
104  const double& etmin,
105  const double& etmax,
106  const double& e1min,
107  const double& e1max,
108  const double& e2min,
109  const double& e2max) :
110  m_response_d1(response_d1),
111  m_response_d2(response_d2),
112  m_etrue1(etrue1),
113  m_etrue2(etrue2),
114  m_cos_phibar(std::cos(phibar*gammalib::deg2rad)),
115  m_sin_phibar(std::sin(phibar*gammalib::deg2rad)),
116  m_etmin(etmin),
117  m_etmax(etmax),
118  m_e1min(e1min),
119  m_e1max(e1max),
120  m_e2min(e2min),
121  m_e2max(e2max) {}
122  double eval(const double& energy1);
123  protected:
124  const GCOMD1Response& m_response_d1; //!< Reference to D1 module response
125  const GCOMD2Response& m_response_d2; //!< Reference to D2 module response
126  double m_etrue1; //!< True D1 energy (MeV)
127  double m_etrue2; //!< True D2 energy (MeV)
128  double m_cos_phibar; //!< cos(phibar)
129  double m_sin_phibar; //!< sin(phibar)
130  double m_etmin; //!< Minimum total energy (MeV)
131  double m_etmax; //!< Maximum total energy (MeV)
132  double m_e1min; //!< Minimum D1 energy (MeV)
133  double m_e1max; //!< Maximum D1 energy (MeV)
134  double m_e2min; //!< Minimum D2 energy (MeV)
135  double m_e2max; //!< Maximum D2 energy (MeV)
136  };
137 
138  // Smearing integration kernel
139  class smearing_kernel : public GFunction {
140  public:
141  smearing_kernel(const double& phigeo,
142  const double& sigma,
143  const GNodeArray& phigeos,
144  const std::vector<double>& values) :
145  m_phigeos(phigeos),
146  m_values(values),
147  m_phigeo(phigeo),
148  m_wgt(1.0/sigma),
149  m_norm(1.0/(sigma*gammalib::sqrt_twopi)) {}
150  double eval(const double& phigeo);
151  protected:
152  const GNodeArray& m_phigeos; //!< Geometrical scatter angles
153  const std::vector<double>& m_values; //!< IAQ values
154  double m_phigeo; //!< Current phigeo value
155  double m_wgt; //!< Inverse of Gaussian std (1/deg)
156  double m_norm; //!< Gaussian normalisation
157  };
158 
159  // Private data members
160  GFitsImageFloat m_iaq; //!< Response
161  GEbounds m_ebounds; //!< Energy boundaries
162  GCOMD1Response m_response_d1; //!< D1 module response
163  GCOMD2Response m_response_d2; //!< D2 module response
164  GCOMInstChars m_ict; //!< Instrument characteristics
165  double m_phigeo_max; //!< Maximum geometrical scatter angle (deg)
166  double m_phibar_max; //!< Maximum Compton scatter angle (deg)
167  double m_phigeo_bin_size; //!< Bin size in geometrical scatter angle (deg)
168  double m_phibar_bin_size; //!< Bin size in Compton scatter angle (deg)
169  double m_phibar_resolution; //!< Bin size for oversampling (deg)
170  double m_e1min; //!< Minimum D1 energy (MeV)
171  double m_e1max; //!< Maximum D1 energy (MeV)
172  double m_e2min; //!< Minimum D2 energy (MeV)
173  double m_e2max; //!< Maximum D2 energy (MeV)
174  int m_num_energies; //!< Number of energies for continuum IAQ
175  bool m_psd_correct; //!< PSD correction usage flag
176  double m_zenith; //!< Zenith angle for location smearing (deg)
177 };
178 
179 
180 /***********************************************************************//**
181  * @brief Return class name
182  *
183  * @return String containing the class name ("GCOMIaq").
184  ***************************************************************************/
185 inline
186 std::string GCOMIaq::classname(void) const
187 {
188  return ("GCOMIaq");
189 }
190 
191 #endif /* GCOMIAQ_HPP */
const double sqrt_twopi
Definition: GMath.hpp:54
COMPTEL Instrument Characteristics interface definition.
double m_phigeo_bin_size
Bin size in geometrical scatter angle (deg)
Definition: GCOMIaq.hpp:167
std::string classname(void) const
Return class name.
Definition: GCOMIaq.hpp:186
void compton_kinematics(const double &energy)
Generate IAQ matrix based on Compton kinematics.
Definition: GCOMIaq.cpp:681
double m_norm
Gaussian normalisation.
Definition: GCOMIaq.hpp:156
Node array class.
Definition: GNodeArray.hpp:60
double m_wgt
Inverse of Gaussian std (1/deg)
Definition: GCOMIaq.hpp:155
double m_phibar_bin_size
Bin size in Compton scatter angle (deg)
Definition: GCOMIaq.hpp:168
Abstract spectral model base class.
double eval(const double &phigeo)
Computes product of D1 and D2 responses at energy E1 and phibar.
Definition: GCOMIaq.cpp:1463
Interface for the COMPTEL instrument response representation class.
Definition: GCOMIaq.hpp:52
void init_response(void)
Initialise COMPTEL response function and instrument characteristics.
Definition: GCOMIaq.cpp:621
double m_e1min
Minimum D1 energy (MeV)
Definition: GCOMIaq.hpp:132
double m_phigeo
Current phigeo value.
Definition: GCOMIaq.hpp:154
Single precision FITS image class definition.
double m_etrue1
True D1 energy (MeV)
Definition: GCOMIaq.hpp:126
void save(const GFilename &filename, const bool &clobber=false) const
Save COMPTEL instrument response representation into FITS file.
Definition: GCOMIaq.cpp:421
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1100
COMPTEL D2 module response class interface definition.
double m_e1max
Maximum D1 energy (MeV)
Definition: GCOMIaq.hpp:133
smearing_kernel(const double &phigeo, const double &sigma, const GNodeArray &phigeos, const std::vector< double > &values)
Definition: GCOMIaq.hpp:141
double m_phigeo_max
Maximum geometrical scatter angle (deg)
Definition: GCOMIaq.hpp:165
void set(const GEnergy &energy, const GEbounds &ebounds)
Set mono-energetic IAQ.
Definition: GCOMIaq.cpp:253
double m_e1max
Maximum D1 energy (MeV)
Definition: GCOMIaq.hpp:171
double m_sin_phibar
sin(phibar)
Definition: GCOMIaq.hpp:129
GCOMIaq(void)
Void constructor.
Definition: GCOMIaq.cpp:78
std::vector< double > location_spread(const double &zenith) const
Compute location spread vector.
Definition: GCOMIaq.cpp:1247
Interface for the COMPTEL Instrument Characteristics class.
double eval(const double &energy1)
Computes product of D1 and D2 responses at energy E1 and phibar.
Definition: GCOMIaq.cpp:1366
double m_phibar_resolution
Bin size for oversampling (deg)
Definition: GCOMIaq.hpp:169
double m_phibar_max
Maximum Compton scatter angle (deg)
Definition: GCOMIaq.hpp:166
void klein_nishina(const double &energy)
Multiply IAQ matrix by the Klein-Nishina formula.
Definition: GCOMIaq.cpp:866
void location_smearing(const double &zenith)
Perform location smearing.
Definition: GCOMIaq.cpp:1154
double m_e2max
Maximum D2 energy (MeV)
Definition: GCOMIaq.hpp:135
double m_etmin
Minimum total energy (MeV)
Definition: GCOMIaq.hpp:130
double m_e2max
Maximum D2 energy (MeV)
Definition: GCOMIaq.hpp:173
const double deg2rad
Definition: GMath.hpp:43
Energy boundaries container class.
Definition: GEbounds.hpp:60
Single parameter function abstract base class definition.
Interface for the COMPTEL D2 module response class.
Interface for the COMPTEL D1 module response class.
double compute_iaq_bin(const double &etrue1, const double &etrue2, const double &phibar)
Computes the IAQ for one bin.
Definition: GCOMIaq.cpp:791
Filename class.
Definition: GFilename.hpp:62
double klein_nishina_bin(const double &energy, const double &phigeo)
Computes Klein-Nishina probability for one bin.
Definition: GCOMIaq.cpp:935
Interface class for all GammaLib classes.
Definition: GBase.hpp:52
GEbounds m_ebounds
Energy boundaries.
Definition: GCOMIaq.hpp:161
GCOMD1Response m_response_d1
D1 module response.
Definition: GCOMIaq.hpp:162
bool m_psd_correct
PSD correction usage flag.
Definition: GCOMIaq.hpp:175
GChatter
Definition: GTypemaps.hpp:33
response_kernel(const GCOMD1Response &response_d1, const GCOMD2Response &response_d2, const double &etrue1, const double &etrue2, const double &phibar, const double &etmin, const double &etmax, const double &e1min, const double &e1max, const double &e2min, const double &e2max)
Definition: GCOMIaq.hpp:99
int m_num_energies
Number of energies for continuum IAQ.
Definition: GCOMIaq.hpp:174
GCOMD2Response m_response_d2
D2 module response.
Definition: GCOMIaq.hpp:163
double m_etrue2
True D2 energy (MeV)
Definition: GCOMIaq.hpp:127
double m_e1min
Minimum D1 energy (MeV)
Definition: GCOMIaq.hpp:170
Single precision FITS image class.
COMPTEL D1 module response class interface definition.
GCOMIaq & operator=(const GCOMIaq &iaq)
Assignment operator.
Definition: GCOMIaq.cpp:183
GCOMIaq * clone(void) const
Clone instance.
Definition: GCOMIaq.cpp:236
double m_cos_phibar
cos(phibar)
Definition: GCOMIaq.hpp:128
GCOMInstChars m_ict
Instrument characteristics.
Definition: GCOMIaq.hpp:164
double m_e2min
Minimum D2 energy (MeV)
Definition: GCOMIaq.hpp:134
const GNodeArray & m_phigeos
Geometrical scatter angles.
Definition: GCOMIaq.hpp:152
double m_zenith
Zenith angle for location smearing (deg)
Definition: GCOMIaq.hpp:176
double m_e2min
Minimum D2 energy (MeV)
Definition: GCOMIaq.hpp:172
void free_members(void)
Delete class members.
Definition: GCOMIaq.cpp:611
const GCOMD1Response & m_response_d1
Reference to D1 module response.
Definition: GCOMIaq.hpp:124
Single parameter function abstract base class.
Definition: GFunction.hpp:44
void weight_iaq(const double &energy)
Weight IAQ matrix.
Definition: GCOMIaq.cpp:1036
GFitsImageFloat m_iaq
Response.
Definition: GCOMIaq.hpp:160
std::string print(const GChatter &chatter=NORMAL) const
Print COMPTEL instrument response representation information.
Definition: GCOMIaq.cpp:464
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1226
Energy boundaries class interface definition.
void init_members(void)
Initialise class members.
Definition: GCOMIaq.cpp:547
void remove_cards(void)
Remove any extra header cards.
Definition: GCOMIaq.cpp:639
~GCOMIaq(void)
Destructor.
Definition: GCOMIaq.cpp:161
const GCOMD2Response & m_response_d2
Reference to D2 module response.
Definition: GCOMIaq.hpp:125
double m_etmax
Maximum total energy (MeV)
Definition: GCOMIaq.hpp:131
void copy_members(const GCOMIaq &iaq)
Copy class members.
Definition: GCOMIaq.cpp:580
void clear(void)
Clear instance.
Definition: GCOMIaq.cpp:217
double klein_nishina_integral(const double &v, const double &a)
Computes Klein-Nishina integral.
Definition: GCOMIaq.cpp:1003
Mathematical function definitions.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
const std::vector< double > & m_values
IAQ values.
Definition: GCOMIaq.hpp:153