GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
com_helpers_response_vector.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * COMPTEL helper classes for vector response *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2021 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 com_helpers_response_vector.cpp
23  * @brief Implementation of COMPTEL helper classes for vector response
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <cmath>
33 #include "GCOMResponse.hpp"
34 #include "GIntegral.hpp"
35 #include "GIntegrals.hpp"
36 
37 /* __ Method name definitions ____________________________________________ */
38 
39 /* __ Macros _____________________________________________________________ */
40 
41 /* __ Coding definitions _________________________________________________ */
42 
43 /* __ Debug definitions __________________________________________________ */
44 
45 /* __ Constants __________________________________________________________ */
46 
47 
48 /*==========================================================================
49  = =
50  = Helper methods =
51  = =
52  ==========================================================================*/
53 
54 /***********************************************************************//**
55  * @brief Kernel for radial integration of radial models
56  *
57  * @param[in] rho Rho angle (radians).
58  * @return Azimuthally integrated radial model.
59  ***************************************************************************/
61 {
62  // Initialise kernel values
63  m_irfs = 0.0;
64 
65  // Continue only if rho is positive
66  if (rho > 0.0) {
67 
68  // Compute radial model value
69  double model = m_model.eval(rho, m_bin->energy(), m_bin->time());
70 
71  // Continue only if model is positive
72  if (model > 0.0) {
73 
74  // Precompute cosine and sine terms for azimuthal integration
75  double sin_rho = std::sin(rho);
76  double cos_rho = std::cos(rho);
77 
78  // Setup azimuthal integration kernel
79  com_radial_kerns_omega integrands(m_iaq,
80  m_irfs,
81  m_bin,
82  m_rot,
83  m_drx,
87  sin_rho,
88  cos_rho);
89 
90  // Setup integrator
91  GIntegrals integral(&integrands);
92  integral.fixed_iter(m_iter);
93 
94  // Integrate over Omega angle
95  m_irfs = integral.romberg(0.0, gammalib::twopi, m_iter) *
96  model * sin_rho;
97 
98  } // endif: radial model was positive
99 
100  } // endif: phigeo was positive
101 
102  // Return kernel values
103  return m_irfs;
104 }
105 
106 
107 /***********************************************************************//**
108  * @brief Kernel for azimuthal integration of radial models
109  *
110  * @param[in] omega Omega angle (radians).
111  * @return Kernel value for radial model.
112  ***************************************************************************/
114 {
115  // Initialise kernel values
116  m_irfs = 0.0;
117 
118  // Compute sine and cosine of azimuth angle
119  double sin_omega = std::sin(omega);
120  double cos_omega = std::cos(omega);
121 
122  // Get sky direction
123  GVector native(-cos_omega*m_sin_rho, sin_omega*m_sin_rho, m_cos_rho);
124  GVector dir = m_rot * native;
125  GSkyDir skyDir;
126  skyDir.celvector(dir);
127 
128  // Compute Phigeo of current position
129  double phigeo = m_bin->dir().dir().dist(skyDir);
130 
131  // Precompute interpolated IAQ vector for Phigeo angle
132  double phirat = phigeo / m_phigeo_bin_size; // 0.5 at bin centre
133  int iphigeo = int(phirat); // index into which Phigeo falls
134  double eps = phirat - iphigeo - 0.5; // 0.0 at bin centre [-0.5, 0.5[
135 
136  // Continue only if Phigeo index is valid
137  if (iphigeo < m_phigeo_bins) {
138 
139  // Get DRX value
140  double intensity = m_drx(skyDir);
141 
142  // Multiply result with IAQ for all Phibar layers
143  for (int iphibar = 0; iphibar < m_phibar_bins; ++iphibar) {
144 
145  // Initialise IAQ
146  double iaq = 0.0;
147 
148  // Get IAQ index
149  int i = iphibar * m_phigeo_bins + iphigeo;
150 
151  // Get interpolated IAQ value
152  if (eps < 0.0) { // interpolate towards left
153  if (iphigeo > 0) {
154  iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
155  }
156  else {
157  iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
158  }
159  }
160  else { // interpolate towards right
161  if (iphigeo < m_phigeo_bins-1) {
162  iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
163  }
164  else {
165  iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
166  }
167  }
168 
169  // If interpolated IAQ value is positive then compute IRF value
170  if (iaq > 0.0) {
171  m_irfs[iphibar] = intensity * iaq;
172  }
173 
174  } // endfor: looped over Phibar
175 
176  } // endif: Phigeo index was valid
177 
178  // Return kernel values
179  return m_irfs;
180 }
181 
182 
183 /***********************************************************************//**
184  * @brief Kernel for radial integration of elliptical models
185  *
186  * @param[in] rho Rho angle (radians).
187  * @return Azimuthally integrated elliptical model.
188  ***************************************************************************/
190 {
191  // Initialise kernel values
192  m_irfs = 0.0;
193 
194  // Continue only if rho is positive
195  if (rho > 0.0) {
196 
197  // Precompute cosine and sine terms for azimuthal integration
198  double sin_rho = std::sin(rho);
199  double cos_rho = std::cos(rho);
200 
201  // Setup azimuthal integration kernel
203  m_model,
204  m_irfs,
205  m_bin,
206  m_rot,
207  m_drx,
211  sin_rho,
212  cos_rho);
213 
214  // Setup integrator
215  GIntegrals integral(&integrands);
216  integral.fixed_iter(m_iter);
217 
218  // Integrate over Omega angle
219  m_irfs = integral.romberg(0.0, gammalib::twopi, m_iter) * sin_rho;
220 
221  } // endif: phigeo was positive
222 
223  // Return kernel values
224  return m_irfs;
225 }
226 
227 
228 /***********************************************************************//**
229  * @brief Kernel for azimuthal integration of elliptical models
230  *
231  * @param[in] omega Omega angle (radians).
232  * @return Kernel value for elliptical model.
233  ***************************************************************************/
235 {
236  // Initialise kernel values
237  m_irfs = 0.0;
238 
239  // Compute sine and cosine of azimuth angle
240  double sin_omega = std::sin(omega);
241  double cos_omega = std::cos(omega);
242 
243  // Get sky direction
244  GVector native(-cos_omega*m_sin_rho, sin_omega*m_sin_rho, m_cos_rho);
245  GVector dir = m_rot * native;
246  GSkyDir skyDir;
247  skyDir.celvector(dir);
248 
249  // Compute Phigeo of current position
250  double phigeo = m_bin->dir().dir().dist(skyDir);
251 
252  // Precompute interpolated IAQ vector for Phigeo angle
253  double phirat = phigeo / m_phigeo_bin_size; // 0.5 at bin centre
254  int iphigeo = int(phirat); // index into which Phigeo falls
255  double eps = phirat - iphigeo - 0.5; // 0.0 at bin centre [-0.5, 0.5[
256 
257  // Continue only if Phigeo index is valid
258  if (iphigeo < m_phigeo_bins) {
259 
260  // Set photon
261  GPhoton photon(skyDir, m_bin->energy(), m_bin->time());
262 
263  // Get model sky intensity for photon (unit: sr^-1)
264  double intensity = m_model.spatial()->eval(photon);
265 
266  // Continue only if intensity is positive
267  if (intensity > 0.0) {
268 
269  // Multiply-in DRX value
270  intensity *= m_drx(skyDir);
271 
272  // Multiply result with IAQ for all Phibar layers
273  for (int iphibar = 0; iphibar < m_phibar_bins; ++iphibar) {
274 
275  // Initialise IAQ
276  double iaq = 0.0;
277 
278  // Get IAQ index
279  int i = iphibar * m_phigeo_bins + iphigeo;
280 
281  // Get interpolated IAQ value
282  if (eps < 0.0) { // interpolate towards left
283  if (iphigeo > 0) {
284  iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
285  }
286  else {
287  iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
288  }
289  }
290  else { // interpolate towards right
291  if (iphigeo < m_phigeo_bins-1) {
292  iaq = (1.0 - eps) * m_iaq[i] + eps * m_iaq[i+1];
293  }
294  else {
295  iaq = (1.0 + eps) * m_iaq[i] - eps * m_iaq[i-1];
296  }
297  }
298 
299  // If interpolated IAQ value is positive then compute IRF value
300  if (iaq > 0.0) {
301  m_irfs[iphibar] = intensity * iaq;
302  }
303 
304  } // endfor: looped over Phibar
305 
306  } // endif: intensity was valid
307 
308  } // endif: Phigeo index was valid
309 
310  // Return kernel values
311  return m_irfs;
312 }
const GModelSpatialRadial & m_model
Radial spatial model.
const int & m_phigeo_bins
Number of phigeo bins.
Kernel for azimuth angle integration of elliptical models.
const int & m_phibar_bins
Number of phibar bins.
virtual double eval(const GPhoton &photon, const bool &gradients=false) const =0
const double & m_phigeo_bin_size
Phigeo bin size.
const double & m_phigeo_bin_size
Phigeo bin size.
const int & m_phibar_bins
Number of phibar bins.
Kernel for azimuth angle integration of radial models.
const double & m_cos_rho
Cosine of Rho.
const GCOMEventBin * m_bin
Event bin.
const GMatrix & m_rot
Rotation matrix.
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1190
const GModelSky & m_model
Sky model.
const double & m_phigeo_bin_size
Phigeo bin size.
const GCOMEventBin * m_bin
Event bin.
GVector eval(const double &omega)
Kernel for azimuthal integration of radial models.
GVector & m_irfs
IRF vector to update.
void fixed_iter(const int &iter)
Set fixed number of iterations.
Definition: GIntegrals.hpp:179
const GCOMEventBin * m_bin
Event bin.
const int & m_iter
Number of omega iterations.
const int & m_phibar_bins
Number of phibar bins.
virtual const GTime & time(void) const
Return time of event bin.
GVector eval(const double &omega)
Kernel for azimuthal integration of elliptical models.
const GModelSky & m_model
Sky model.
Class that handles photons.
Definition: GPhoton.hpp:47
virtual double eval(const double &theta, const GEnergy &energy, const GTime &time, const bool &gradients=false) const =0
const GMatrix & m_rot
Rotation matrix.
const int & m_phigeo_bins
Number of phigeo bins.
const int & m_phibar_bins
Number of phibar bins.
COMPTEL instrument response function class interface definition.
GVector eval(const double &phigeo)
Kernel for radial integration of elliptical models.
const GMatrix & m_rot
Rotation matrix.
void dir(const GSkyDir &dir)
Set event scatter direction.
Integration class for set of functions interface definition.
const std::vector< double > & m_iaq
IAQ vector.
void celvector(const GVector &vector)
Set sky direction from 3D vector in celestial coordinates.
Definition: GSkyDir.cpp:313
const double & m_cos_rho
Cosine of Rho.
GVector & m_irfs
IRF vector to update.
const double & m_sin_rho
Sine of Rho.
const double & m_phigeo_bin_size
Phigeo bin size.
GVector & m_irfs
IRF vector to update.
GVector & m_irfs
IRF vector to update.
Defintion of COMPTEL helper classes for vector response.
const int & m_iter
Number of omega iterations.
GVector eval(const double &phigeo)
Kernel for radial integration of radial models.
const std::vector< double > & m_iaq
IAQ vector.
virtual const GCOMInstDir & dir(void) const
Return instrument direction of event bin.
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1316
GVector romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegrals.cpp:210
GModelSpatial * spatial(void) const
Return spatial model component.
Definition: GModelSky.hpp:295
const int & m_phigeo_bins
Number of phigeo bins.
const std::vector< double > & m_iaq
IAQ vector.
const GCOMEventBin * m_bin
Event bin.
const double twopi
Definition: GMath.hpp:36
const GMatrix & m_rot
Rotation matrix.
Vector class.
Definition: GVector.hpp:46
const int & m_phigeo_bins
Number of phigeo bins.
Integration class for set of functions.
Definition: GIntegrals.hpp:47
Integration class interface definition.
Sky direction class.
Definition: GSkyDir.hpp:62
const std::vector< double > & m_iaq
IAQ vector.
virtual const GEnergy & energy(void) const
Return energy of event bin.
const double & m_sin_rho
Sine of Rho.