GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCOMD2Response.hpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCOMD2Response.hpp - COMPTEL D2 module response class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2017-2018 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 GCOMD2Response.hpp
23  * @brief COMPTEL D2 module response class interface definition
24  * @author Juergen Knoedlseder
25  */
26 
27 #ifndef GCOMD2RESPONSE_HPP
28 #define GCOMD2RESPONSE_HPP
29 
30 /* __ Includes ___________________________________________________________ */
31 #include <vector>
32 #include <string>
33 #include "GCaldb.hpp"
34 #include "GNodeArray.hpp"
35 #include "GFunction.hpp"
36 
37 /* __ Type definitions ___________________________________________________ */
38 
39 /* __ Forward declarations _______________________________________________ */
40 class GFitsTable;
41 class GFitsBinTable;
42 
43 
44 /***********************************************************************//**
45  * @class GCOMD2Response
46  *
47  * @brief Interface for the COMPTEL D2 module response class
48  ***************************************************************************/
49 class GCOMD2Response : public GBase {
50 
51 public:
52  // Constructors and destructors
53  GCOMD2Response(void);
54  GCOMD2Response(const GCOMD2Response& rsp);
55  GCOMD2Response(const GCaldb& caldb, const std::string& sdbname);
56  ~GCOMD2Response(void);
57 
58  // Operators
60  double operator()(const double& etrue, const double& ereco) const;
61 
62  // Methods
63  void clear(void);
64  GCOMD2Response* clone(void) const;
65  std::string classname(void) const;
66  void caldb(const GCaldb& caldb);
67  const GCaldb& caldb(void) const;
68  void load(const std::string& sdbname);
69  void read(const GFitsTable& table);
70  void write(GFitsBinTable& table);
71  double position(const double& etrue) const;
72  double sigma(const double& etrue) const;
73  double amplitude(const double& etrue) const;
74  double escape1(const double& etrue) const;
75  double escape2(const double& etrue) const;
76  double comptontail(const double& etrue) const;
77  double background(const double& etrue) const;
78  double emin(const double& etrue) const;
79  double ewidth(const double& etrue) const;
80  double emax(const double& etrue) const;
81  double emin(void) const;
82  double emax(void) const;
83  std::string print(const GChatter& chatter = NORMAL) const;
84 
85 private:
86  // Private methods
87  void init_members(void);
88  void copy_members(const GCOMD2Response& rsp);
89  void free_members(void);
90  void update_cache(const double& etrue) const;
91  void update_response_vector(const double& etrue) const;
92 
93  // Klein-Nishina Gaussian integration kernel
94  class kn_gauss_kernel : public GFunction {
95  public:
96  kn_gauss_kernel(const double& ereco,
97  const double& e0,
98  const double& ec,
99  const double& sigma) :
100  m_ereco(ereco),
101  m_e0(e0),
102  m_ec(ec),
103  m_wgt(1.0/sigma) {}
104  double eval(const double& e);
105  protected:
106  double m_ereco; //< Reconstructed energy (MeV)
107  double m_e0; //!< Incident energy (MeV)
108  double m_ec; //!< Compton edge energy (MeV)
109  double m_wgt; //!< Inverse of Gaussian standard deviation (1/MeV)
110  };
111 
112  // Background Gaussian integration kernel
113  class bkg_gauss_kernel : public GFunction {
114  public:
115  bkg_gauss_kernel(const double& ereco,
116  const double& e0,
117  const double& sigma) :
118  m_ereco(ereco),
119  m_e0(e0),
120  m_wgt(1.0/sigma) {}
121  double eval(const double& e);
122  protected:
123  double m_ereco; //< Reconstructed energy (MeV)
124  double m_e0; //!< Incident energy (MeV)
125  double m_wgt; //!< Inverse of Gaussian standard deviation (1/MeV)
126  };
127 
128  // Private data members
129  GCaldb m_caldb; //!< Calibration database
130  GNodeArray m_energies; //!< Input energies
131  std::vector<double> m_positions; //!< Photo peak position in MeV
132  std::vector<double> m_sigmas; //!< Photo peak width in MeV
133  std::vector<double> m_amplitudes; //!< Photo peak amplitude
134  std::vector<double> m_escapes1; //!< Amplitude of first escape peak
135  std::vector<double> m_escapes2; //!< Amplitude of second escape peak
136  std::vector<double> m_tails; //!< Amplitude of Compton tail
137  std::vector<double> m_backgrounds; //!< Amplitude of Compton background
138  std::vector<double> m_emins; //!< Lower energy threshold of D2
139  std::vector<double> m_ewidths; //!< Lower energy threshold width of D2
140  std::vector<double> m_emaxs; //!< Upper energy limit of D2
141 
142  // Pre-computation cache
143  mutable double m_energy; //!< Incident total energy (MeV)
144  mutable double m_position; //!< Position of photo peak (MeV)
145  mutable double m_sigma; //!< Width of photo peak (MeV)
146  mutable double m_amplitude; //!< Amplitude of photo peak
147  mutable double m_escape1; //!< Amplitude of first escape peak
148  mutable double m_escape2; //!< Amplitude of second escape peak
149  mutable double m_tail; //!< Amplitude of Compton tail
150  mutable double m_background; //!> Amplitude of Compton background
151  mutable double m_emin; //!< Lower energy threshold of D2 (MeV)
152  mutable double m_ewidth; //!< Lower energy threshold width of D2 (MeV)
153  mutable double m_emax; //!< Upper energy limit of D2 (MeV)
154  mutable double m_pos_escape1; //!< Position of first escape peak (MeV)
155  mutable double m_pos_escape2; //!< Position of second escape peak (MeV)
156  mutable double m_wgt_photo; //!< Inverse of width of photo peak (1/MeV)
157  mutable double m_wgt_escape1; //!< Inverse of width of first escape peak (1/MeV)
158  mutable double m_wgt_escape2; //!< Inverse of width of first escape peak (1/MeV)
159  mutable double m_compton_edge; //!< Position of Compton edge (MeV)
160 
161  // Pre-computed response vector
162  mutable double m_rsp_etrue; //!< True energy of response vector
163  mutable GNodeArray m_rsp_energies; //!< Response vector energies
164  mutable std::vector<double> m_rsp_values; //!< Response vector values
165 
166 };
167 
168 
169 /***********************************************************************//**
170  * @brief Return class name
171  *
172  * @return String containing the class name ("GCOMD2Response").
173  ***************************************************************************/
174 inline
175 std::string GCOMD2Response::classname(void) const
176 {
177  return ("GCOMD2Response");
178 }
179 
180 
181 /***********************************************************************//**
182  * @brief Return calibration database
183  *
184  * @return Calibration database.
185  ***************************************************************************/
186 inline
187 const GCaldb& GCOMD2Response::caldb(void) const
188 {
189  return m_caldb;
190 }
191 
192 
193 /***********************************************************************//**
194  * @brief Set calibration database
195  *
196  * @param[in] caldb Calibration database.
197  *
198  * Sets the calibration database containing the COMPTEL D2 module response.
199  ***************************************************************************/
200 inline
201 void GCOMD2Response::caldb(const GCaldb& caldb)
202 {
203  m_caldb = caldb;
204  return;
205 }
206 
207 
208 /***********************************************************************//**
209  * @brief Return photo peak position
210  *
211  * @param[in] etrue True energy (MeV).
212  * @return Photo peak position (MeV).
213  ***************************************************************************/
214 inline
215 double GCOMD2Response::position(const double& etrue) const
216 {
217  update_cache(etrue);
218  return (m_position);
219 }
220 
221 
222 /***********************************************************************//**
223  * @brief Return photo peak standard deviation
224  *
225  * @param[in] etrue True energy (MeV).
226  * @return Photo peak standard deviation (MeV).
227  ***************************************************************************/
228 inline
229 double GCOMD2Response::sigma(const double& etrue) const
230 {
231  update_cache(etrue);
232  return (m_sigma);
233 }
234 
235 
236 /***********************************************************************//**
237  * @brief Return photo peak amplitude
238  *
239  * @param[in] etrue True energy (MeV).
240  * @return Photo peak amplitude.
241  ***************************************************************************/
242 inline
243 double GCOMD2Response::amplitude(const double& etrue) const
244 {
245  update_cache(etrue);
246  return (m_amplitude);
247 }
248 
249 
250 /***********************************************************************//**
251  * @brief Return first escape peak amplitude
252  *
253  * @param[in] etrue True energy (MeV).
254  * @return First escape peak amplitude.
255  ***************************************************************************/
256 inline
257 double GCOMD2Response::escape1(const double& etrue) const
258 {
259  update_cache(etrue);
260  return (m_escape1);
261 }
262 
263 
264 /***********************************************************************//**
265  * @brief Return second escape peak amplitude
266  *
267  * @param[in] etrue True energy (MeV).
268  * @return Second escape peak amplitude.
269  ***************************************************************************/
270 inline
271 double GCOMD2Response::escape2(const double& etrue) const
272 {
273  update_cache(etrue);
274  return (m_escape2);
275 }
276 
277 
278 /***********************************************************************//**
279  * @brief Return Compton tail amplitude
280  *
281  * @param[in] etrue True energy (MeV).
282  * @return Compton tail amplitude.
283  ***************************************************************************/
284 inline
285 double GCOMD2Response::comptontail(const double& etrue) const
286 {
287  update_cache(etrue);
288  return (m_tail);
289 }
290 
291 
292 /***********************************************************************//**
293  * @brief Return background amplitude
294  *
295  * @param[in] etrue True energy (MeV).
296  * @return Background amplitude.
297  ***************************************************************************/
298 inline
299 double GCOMD2Response::background(const double& etrue) const
300 {
301  update_cache(etrue);
302  return (m_background);
303 }
304 
305 
306 /***********************************************************************//**
307  * @brief Return minimum energy
308  *
309  * @param[in] etrue True energy (MeV).
310  * @return Minimum energy (MeV).
311  ***************************************************************************/
312 inline
313 double GCOMD2Response::emin(const double& etrue) const
314 {
315  update_cache(etrue);
316  return (m_emin);
317 }
318 
319 
320 /***********************************************************************//**
321  * @brief Return energy threshold width
322  *
323  * @param[in] etrue True energy (MeV).
324  * @return Energy threshold width (MeV).
325  ***************************************************************************/
326 inline
327 double GCOMD2Response::ewidth(const double& etrue) const
328 {
329  update_cache(etrue);
330  return (m_ewidth);
331 }
332 
333 
334 /***********************************************************************//**
335  * @brief Return maximum energy
336  *
337  * @param[in] etrue True energy (MeV).
338  * @return Maximum energy (MeV).
339  ***************************************************************************/
340 inline
341 double GCOMD2Response::emax(const double& etrue) const
342 {
343  update_cache(etrue);
344  return (m_emax);
345 }
346 
347 
348 /***********************************************************************//**
349  * @brief Return minimum D2 input energy (MeV)
350  *
351  * @return Minimum energy D2 input energy (MeV).
352  *
353  * Returns the minimum D2 input energy (MeV). In case that no information
354  * has been read from a SDB file so far, the method returns 0.
355  ***************************************************************************/
356 inline
357 double GCOMD2Response::emin(void) const
358 {
359  double emin = (m_energies.size() > 0) ? m_energies[0] : 0.0;
360  return (emin);
361 }
362 
363 
364 /***********************************************************************//**
365  * @brief Return maximum D2 input energy (MeV)
366  *
367  * @return Maximum energy D2 input energy (MeV).
368  *
369  * Returns the maximum D2 input energy (MeV). In case that no information
370  * has been read from a SDB file so far, the method returns 0.
371  ***************************************************************************/
372 inline
373 double GCOMD2Response::emax(void) const
374 {
375  double emax = (m_energies.size() > 0) ? m_energies[m_energies.size()-1] : 0.0;
376  return (emax);
377 }
378 
379 #endif /* GCOMD2RESPONSE_HPP */
GCOMD2Response & operator=(const GCOMD2Response &rsp)
Assignment operator.
GCOMD2Response(void)
Void constructor.
double background(const double &etrue) const
Return background amplitude.
void clear(void)
Clear instance.
int size(void) const
Return number of nodes in node array.
Definition: GNodeArray.hpp:188
void copy_members(const GCOMD2Response &rsp)
Copy class members.
void read(const GFitsTable &table)
Read COMPTEL D2 module response.
double m_wgt_escape2
Inverse of width of first escape peak (1/MeV)
std::vector< double > m_emins
Lower energy threshold of D2.
Node array class.
Definition: GNodeArray.hpp:60
void update_cache(const double &etrue) const
Update computation cache.
double m_sigma
Width of photo peak (MeV)
std::vector< double > m_rsp_values
Response vector values.
double emax(void) const
Return maximum D2 input energy (MeV)
double m_pos_escape2
Position of second escape peak (MeV)
std::vector< double > m_escapes2
Amplitude of second escape peak.
kn_gauss_kernel(const double &ereco, const double &e0, const double &ec, const double &sigma)
void write(GFitsBinTable &table)
Write COMPTEL D2 module response.
std::vector< double > m_amplitudes
Photo peak amplitude.
double m_wgt_photo
Inverse of width of photo peak (1/MeV)
double m_escape2
Amplitude of second escape peak.
double sigma(const double &etrue) const
Return photo peak standard deviation.
double m_e0
Incident energy (MeV)
double m_escape1
Amplitude of first escape peak.
std::string print(const GChatter &chatter=NORMAL) const
Print COMPTEL D2 module response information.
double ewidth(const double &etrue) const
Return energy threshold width.
void init_members(void)
Initialise class members.
double m_position
Position of photo peak (MeV)
void free_members(void)
Delete class members.
double m_amplitude
Amplitude of photo peak.
Node array class interface definition.
Calibration database class.
Definition: GCaldb.hpp:66
double eval(const double &e)
Computes modified Klein-Nishina cross section multiplied with Gaussian.
Single parameter function abstract base class definition.
std::string classname(void) const
Return class name.
Interface for the COMPTEL D2 module response class.
Interface class for all GammaLib classes.
Definition: GBase.hpp:52
double m_tail
Amplitude of Compton tail.
void load(const std::string &sdbname)
Load COMPTEL D2 module response.
double m_ec
Compton edge energy (MeV)
std::vector< double > m_sigmas
Photo peak width in MeV.
double m_emax
Upper energy limit of D2 (MeV)
double m_energy
Incident total energy (MeV)
double m_emin
Amplitude of Compton background
~GCOMD2Response(void)
Destructor.
std::vector< double > m_escapes1
Amplitude of first escape peak.
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
std::vector< double > m_emaxs
Upper energy limit of D2.
double m_ewidth
Lower energy threshold width of D2 (MeV)
double emin(void) const
Return minimum D2 input energy (MeV)
GChatter
Definition: GTypemaps.hpp:33
const GCaldb & caldb(void) const
Return calibration database.
double m_wgt_escape1
Inverse of width of first escape peak (1/MeV)
double m_wgt
Inverse of Gaussian standard deviation (1/MeV)
GCOMD2Response * clone(void) const
Clone instance.
double m_pos_escape1
Position of first escape peak (MeV)
double m_wgt
Inverse of Gaussian standard deviation (1/MeV)
Calibration database class interface definition.
GNodeArray m_energies
Input energies.
std::vector< double > m_ewidths
Lower energy threshold width of D2.
std::vector< double > m_backgrounds
Amplitude of Compton background.
double m_compton_edge
Position of Compton edge (MeV)
void update_response_vector(const double &etrue) const
Update response vector.
Single parameter function abstract base class.
Definition: GFunction.hpp:44
FITS binary table class.
GNodeArray m_rsp_energies
Response vector energies.
std::vector< double > m_positions
Photo peak position in MeV.
double escape2(const double &etrue) const
Return second escape peak amplitude.
GCaldb m_caldb
Calibration database.
std::vector< double > m_tails
Amplitude of Compton tail.
double escape1(const double &etrue) const
Return first escape peak amplitude.
double eval(const double &e)
Computes Compton background multiplied with Gaussian.
double position(const double &etrue) const
Return photo peak position.
double m_e0
Incident energy (MeV)
double amplitude(const double &etrue) const
Return photo peak amplitude.
double comptontail(const double &etrue) const
Return Compton tail amplitude.
bkg_gauss_kernel(const double &ereco, const double &e0, const double &sigma)
double m_rsp_etrue
True energy of response vector.
double operator()(const double &etrue, const double &ereco) const
D2 module response evaluation operator.