GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
cta_helpers_response_stacked_vector.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * CTA helper classes for stacked vector response *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2020-2022 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 cta_helpers_response_stacked_vector.hpp
23  * @brief Implementation of CTA helper classes for stacked 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 "GCTAResponse_helpers.hpp"
34 #include "GModelPar.hpp"
35 #include "GModelSpatialRadial.hpp"
36 #include "GCTAResponseCube.hpp"
37 #include "GIntegrals.hpp"
38 
39 /* __ Method name definitions ____________________________________________ */
40 #define G_PSF_RADIAL_KERNS_PHI "cta_psf_radial_kerns_phi::eval(double&)"
41 
42 /* __ Macros _____________________________________________________________ */
43 
44 /* __ Coding definitions _________________________________________________ */
45 
46 /* __ Debug definitions __________________________________________________ */
47 
48 /* __ Constants __________________________________________________________ */
49 
50 
51 /*==========================================================================
52  = =
53  = Helper methods =
54  = =
55  ==========================================================================*/
56 
57 /***********************************************************************//**
58  * @brief Return size of vector kernel for PSF integration of radial model
59  *
60  * @return Size of vector kernel.
61  *
62  * The size of the vector kernel depends on whether gradient computation is
63  * requested or not. With gradient computation, the size is given by the
64  * product of the number of true energies and the number of model parameters
65  * plus one, without gradient computation the size is simply the number of
66  * true energies.
67  ***************************************************************************/
69 {
70  // Set size
71  int nengs = m_srcEngs.size();
72  int npars = m_model->size();
73  int size = (m_grad) ? nengs * (npars+1) : nengs;
74 
75  // Return size
76  return size;
77 }
78 
79 
80 /***********************************************************************//**
81  * @brief Kernel constructor for PSF integration of radial model
82  *
83  * @param[in] rsp Pointer to CTA response cube.
84  * @param[in] model Pointer to radial spatial model.
85  * @param[in] obsDir Observed sky direction.
86  * @param[in] srcEngs True energies.
87  * @param[in] zeta Angular distance between event direction of model centre (radians).
88  * @param[in] theta_max Maximum model radius (radians).
89  * @param[in] iter Number of iterations for azimuth integration.
90  * @param[in] grad Compute parameter gradients?
91  ***************************************************************************/
93  const GModelSpatialRadial* model,
94  const GSkyDir& obsDir,
95  const GEnergies& srcEngs,
96  const double& zeta,
97  const double& theta_max,
98  const int& iter,
99  const bool& grad)
100 {
101  // Store constructor arguments
102  m_rsp = rsp;
103  m_model = model;
104  m_obsDir = obsDir;
105  m_srcEngs = srcEngs;
106  m_zeta = zeta;
107  m_theta_max = theta_max;
108  m_iter = iter;
109  m_grad = grad;
110 
111  // Pre-compute trigonometric functions
112  m_cos_zeta = std::cos(zeta);
113  m_sin_zeta = std::sin(zeta);
114  m_cos_theta_max = std::cos(theta_max);
115 
116  // Set coordinate system
117  m_par_cel = (model->coordsys() == "CEL");
118 
119  // Store references to longitude and latitude parameters
120  if (m_par_cel) {
121  m_par_lon = &((*(const_cast<GModelSpatialRadial*>(model)))["RA"]);
122  m_par_lat = &((*(const_cast<GModelSpatialRadial*>(model)))["DEC"]);
123  }
124  else {
125  m_par_lon = &((*(const_cast<GModelSpatialRadial*>(model)))["GLON"]);
126  m_par_lat = &((*(const_cast<GModelSpatialRadial*>(model)))["GLAT"]);
127  }
128 
129  // Pre-compute constants needed for gradient computation. The following
130  // exceptions need to be handled
131  // * beta_0 = +/-90 deg
132  // * m_sin_zeta = 0
133  // * denom = 0
134  if (grad) {
135  double alpha_0 = m_par_cel ? model->dir().ra() : model->dir().l();
136  double beta_0 = m_par_cel ? model->dir().dec() : model->dir().b();
137  double alpha_reco = m_par_cel ? obsDir.ra() : obsDir.l();
138  double beta_reco = m_par_cel ? obsDir.dec() : obsDir.b();
139  double sin_beta_0 = std::sin(beta_0);
140  double cos_beta_0 = std::cos(beta_0);
141  double tan_beta_0 = std::tan(beta_0); // Exception: beta_0 = 90 deg
142  double sin_beta_reco = std::sin(beta_reco);
143  double cos_beta_reco = std::cos(beta_reco);
144  double sin_dalpha = std::sin(alpha_0 - alpha_reco);
145  double cos_dalpha = std::cos(alpha_0 - alpha_reco);
146  double arg = cos_beta_reco * tan_beta_0 -
147  sin_beta_reco * cos_dalpha;
148  double denom = sin_dalpha * sin_dalpha + arg * arg;
149  if (m_sin_zeta != 0.0) {
150  m_dzeta_dalpha_0 = cos_beta_0 * cos_beta_reco * sin_dalpha / m_sin_zeta;
151  m_dzeta_dbeta_0 = (sin_beta_0 * cos_beta_reco * cos_dalpha -
152  cos_beta_0 * sin_beta_reco) / m_sin_zeta;
153  }
154  else {
155  m_dzeta_dalpha_0 = cos_beta_0 * cos_beta_reco;
156  m_dzeta_dbeta_0 = 1.0;
157  }
158  if (denom != 0.0) {
159  m_dphi_dalpha_0 = (sin_beta_reco - cos_dalpha * cos_beta_reco * tan_beta_0) /
160  denom;
161  m_dphi_dbeta_0 = (sin_dalpha * cos_beta_reco * (1.0 + tan_beta_0 * tan_beta_0)) /
162  denom;
163  }
164  else {
165  m_dphi_dalpha_0 = 0.0;
166  m_dphi_dbeta_0 = 0.0;
167  }
168  }
169  else {
170  m_dzeta_dalpha_0 = 0.0;
171  m_dzeta_dbeta_0 = 0.0;
172  m_dphi_dalpha_0 = 0.0;
173  m_dphi_dbeta_0 = 0.0;
174  }
175 
176  // Return
177  return;
178 }
179 
180 
181 /***********************************************************************//**
182  * @brief Kernel for PSF integration of radial model
183  *
184  * @param[in] delta PSF offset angle (radians).
185  * @return Azimuthally integrated product between PSF and radial model
186  * values for all energies.
187  *
188  * Computes the azimuthally integrated product of point spread function and
189  * the radial model intensity. As the PSF is azimuthally symmetric, it is
190  * not included in the azimuthally integration, but just multiplied on the
191  * azimuthally integrated model. The method returns thus
192  *
193  * \f[
194  * {\rm PSF}(\delta) \times
195  * \int_0^{2\pi} {\rm M}(\delta, \phi) \sin \delta {\rm d}\phi
196  * \f]
197  *
198  * where \f${\rm M}(\delta, \phi)\f$ is the radial model in the coordinate
199  * system of the point spread function, defined by the angle \f$\delta\f$
200  * between the true and the measured photon direction and the azimuth angle
201  * \f$\phi\f$ around the measured photon direction.
202  ***************************************************************************/
204 {
205  // Determine size of the result vector
206  int npars = m_model->size();
207  int nengs = m_srcEngs.size();
208 
209  // Initialise result vector
210  GVector values(size());
211 
212  // If we're at the Psf peak the model is zero (due to the sin(delta)
213  // term. We thus only integrate for positive deltas.
214  if (delta > 0.0) {
215 
216  // Compute half length of the arc (in radians) from a circle with
217  // radius delta that intersects with the model, defined as a circle
218  // with maximum radius m_theta_max
219  double dphi = 0.5 * gammalib::roi_arclength(delta,
220  m_zeta,
221  m_cos_zeta,
222  m_sin_zeta,
223  m_theta_max,
225 
226  // Continue only if arc length is positive
227  if (dphi > 0.0) {
228 
229  // Compute phi integration range
230  double phi_min = -dphi;
231  double phi_max = +dphi;
232 
233  // Precompute cosine and sine terms for azimuthal integration
234  double sin_delta = std::sin(delta);
235  double cos_delta = std::cos(delta);
236  double sin_delta_sin_zeta = sin_delta * m_sin_zeta;
237  double sin_delta_cos_zeta = sin_delta * m_cos_zeta;
238 
239  // If gradients are requested then use vector integrator
240  if (m_grad) {
241 
242  // Precompute terms for gradient computation
243  double cos_delta_sin_zeta = cos_delta * m_sin_zeta;
244  double cos_delta_cos_zeta = cos_delta * m_cos_zeta;
245 
246  // Setup kernel for azimuthal integration of the spatial model
247  cta_psf_radial_kerns_phi integrand(this,
248  sin_delta_sin_zeta,
249  sin_delta_cos_zeta,
250  cos_delta_sin_zeta,
251  cos_delta_cos_zeta);
252 
253  // Setup integrator
254  GIntegrals integral(&integrand);
255  integral.fixed_iter(m_iter);
256 
257  // Integrate over azimuth
258  GVector irf = integral.romberg(phi_min, phi_max, m_iter) * sin_delta;
259 
260  // Extract value and gradients
261  double value = irf[0];
262 
263  // Multiply in energy dependent PSF
264  for (int i = 0; i < nengs; ++i) {
265 
266  // Get Psf for this energy. We approximate here the true sky
267  // direction by the reconstructed sky direction.
268  double psf = m_rsp->psf()(m_obsDir, delta, m_srcEngs[i]);
269 
270  // Continue only if PSF is positive
271  if (psf > 0.0) {
272 
273  // Compute value
274  values[i] = value * psf;
275 
276  // If gradient computation is requested the multiply also
277  // all factor gradients with the Psf value
278  if (m_grad) {
279  for (int k = 1, ig = i+nengs; k <= npars; ++k, ig += nengs) {
280  values[ig] = irf[k] * psf;
281  }
282  }
283 
284  } // endif: PSF was positive
285 
286  } // endfor: looped over energies
287 
288  } // endif: gradient computation was requested
289 
290  // ... otherwise use scalar integration
291  else {
292 
293  // Allocate dummy energy and time to satisfy interface
294  static const GEnergy srcEng;
295  static const GTime srcTime;
296 
297  // Setup kernel for azimuthal integration of the spatial model
299  srcEng,
300  srcTime,
301  sin_delta_sin_zeta,
302  sin_delta_cos_zeta);
303 
304  // Setup integrator
305  GIntegral integral(&integrand);
306  integral.fixed_iter(m_iter);
307 
308  // Integrate over azimuth
309  double value = integral.romberg(phi_min, phi_max, m_iter) * sin_delta;
310 
311  // Multiply in energy dependent PSF
312  for (int i = 0; i < nengs; ++i) {
313 
314  // Get Psf for this energy. We approximate here the true sky
315  // direction by the reconstructed sky direction.
316  double psf = m_rsp->psf()(m_obsDir, delta, m_srcEngs[i]);
317 
318  // If PSF is positive then compute corresponding value
319  if (psf > 0.0) {
320  values[i] = value * psf;
321  }
322 
323  } // endfor: looped over energies
324 
325  } // endelse: scalar integration requested
326 
327  } // endif: arc length was positive
328 
329  } // endif: delta was positive
330 
331  // Return kernel values
332  return values;
333 }
334 
335 
336 /***********************************************************************//**
337  * @brief Kernel for azimuthal radial model integration
338  *
339  * @param[in] phi Azimuth angle (radians).
340  * @return Vector of radial model value and parameter factor gradients.
341  *
342  * Computes the value and parameter factor gradients of the radial model at
343  * the position \f$(\delta,\phi)\f$ given in point spread function
344  * coordinates.
345  *
346  * The angle \f$\theta\f$ of the radial model is computed using
347  *
348  * \f[
349  * \theta = \arccos \left( \cos \delta \cos \zeta +
350  * \sin \delta \sin \zeta \cos \phi \right)
351  * \f]
352  *
353  * where \f$\delta\f$ is the angle between true and measured photon
354  * direction, \f$\zeta\f$ is the angle between model centre and measured
355  * photon direction, and \f$\phi\f$ is the azimuth angle with respect to the
356  * measured photon direction, where \f$\phi=0\f$ corresponds to the
357  * connecting line between model centre and measured photon direction.
358  *
359  * If gradient computation is requested, the method also returns the gradients
360  *
361  * \f[
362  * \frac{\partial f}{\partial \alpha_0}
363  * \f]
364  *
365  * and
366  *
367  * \f[
368  * \frac{\partial f}{\partial \beta_0}
369  * \f]
370  *
371  * in the second and third slot of the returned vector. The following slots
372  * contains gradients with respect to other model parameters.
373  ***************************************************************************/
375 {
376  // Allocate dummy energy and time to satisfy interface
377  static const GEnergy srcEng;
378  static const GTime srcTime;
379 
380  // Compute sin(phi) and cos(phi). The following code allows usage of the
381  // sincos() method in case that gradients are requested.
382  double sin_phi;
383  double cos_phi;
384  if (m_outer->m_grad) {
385  sin_phi = std::sin(phi);
386  cos_phi = std::cos(phi);
387  }
388  else {
389  cos_phi = std::cos(phi);
390  }
391 
392  // Compute radial model theta angle
393  double cos_theta = m_cos_delta_cos_zeta + m_sin_delta_sin_zeta * cos_phi;
394  double theta = std::acos(cos_theta);
395 
396  // Reduce theta by an infinite amount to avoid rounding errors at the
397  // boundary of a sharp edged model
398  double theta_kluge = theta - 1.0e-12;
399  if (theta_kluge < 0.0) {
400  theta_kluge = 0.0;
401  }
402 
403  // If gradients are requested then compute partial derivatives of theta
404  // and phi with respect to Right Ascension and Declination and store them
405  // into the respective factor gradients so that the
406  // GModelSpatialRadial::eval() method can use the information.
407  if (m_outer->m_grad) {
408  double g_lon = 0.0;
409  double g_lat = 0.0;
410  if (theta > 0.0) {
411 
412  // Compute sin(theta) by avoiding to use a call to sine
413  double sin_theta2 = 1.0 - cos_theta * cos_theta;
414  double sin_theta = (sin_theta2 > 0.0) ? std::sqrt(sin_theta2) : 0.0;
415 
416  // Continue only if sine is non-zero
417  if (sin_theta != 0.0) {
418  double norm = 1.0 / sin_theta;
419  double dtheta_dzeta = (m_cos_delta_sin_zeta -
420  m_sin_delta_cos_zeta * cos_phi) *
421  norm;
422  double dtheta_dphi = (m_sin_delta_sin_zeta * sin_phi) *
423  norm;
424  if (m_outer->m_zeta != 0.0) {
425  g_lon = dtheta_dzeta * m_outer->m_dzeta_dalpha_0 +
426  dtheta_dphi * m_outer->m_dphi_dalpha_0;
427  g_lat = dtheta_dzeta * m_outer->m_dzeta_dbeta_0 +
428  dtheta_dphi * m_outer->m_dphi_dbeta_0;
429  }
430  else {
431  g_lon = 0.0; //TODO: Set correct value
432  g_lat = -cos_phi;
433  }
434  }
435  }
438  }
439 
440  // Compute model value and optionally model parameter gradients
441  m_values[0] = m_outer->m_model->eval(theta_kluge, srcEng, srcTime,
442  m_outer->m_grad);
443 
444  // If gradients are requested, extract now the fully computed model
445  // parameter gradients into the kernel values vector
446  if (m_outer->m_grad) {
447  for (int i = 1, k = 0; i < m_size; ++i, ++k) {
448  m_values[i] = (*(m_outer->m_model))[k].factor_gradient();
449  }
450  }
451 
452  // Debug: Check for NaN
453  #if defined(G_NAN_CHECK)
456  std::string msg = "NaN/Inf encountered for phi="+gammalib::str(phi);
458  }
459  #endif
460 
461  // Return kernel values
462  return m_values;
463 }
CTA response helper classes definition.
GVector eval(const double &delta)
Kernel for PSF integration of radial model.
double roi_arclength(const double &rad, const double &dist, const double &cosdist, const double &sindist, const double &roi, const double &cosroi)
Returns length of circular arc within circular ROI.
Definition: GTools.cpp:2107
const double & factor_gradient(void) const
Return parameter factor gradient.
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:864
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:381
CTA cube-style response function class definition.
void warning(const std::string &origin, const std::string &message)
Emits warning.
Definition: GTools.cpp:1386
int size(void) const
Return number of parameters.
const double & b(void) const
Return galactic latitude in radians.
Definition: GSkyDir.hpp:187
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1190
Abstract radial spatial model base class interface definition.
CTA cube-style response function class.
void fixed_iter(const int &iter)
Set fixed number of iterations.
Definition: GIntegrals.hpp:179
GSkyDir m_obsDir
Reconstructed event direction.
std::string coordsys(void) const
Return coordinate system.
double acos(const double &arg)
Computes acos by avoiding NaN due to rounding errors.
Definition: GMath.cpp:69
double m_cos_delta_sin_zeta
cos(delta) * sin(zeta)
cta_psf_radial_kerns_delta * m_outer
Pointer to outer integr.
double m_sin_delta_cos_zeta
sin(delta) * cos(zeta)
Time class.
Definition: GTime.hpp:55
GIntegral class interface definition.
Definition: GIntegral.hpp:46
GVector eval(const double &phi)
Kernel for azimuthal radial model integration.
Energy container class.
Definition: GEnergies.hpp:60
const double & ra(void) const
Return Right Ascension in radians.
Definition: GSkyDir.hpp:214
bool m_par_cel
Celestial or galactic coordinates.
GModelPar * m_par_lat
Latitude parameter.
const GCTAResponseCube * m_rsp
Response cube.
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:201
Model parameter class interface definition.
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:184
virtual double eval(const double &theta, const GEnergy &energy, const GTime &time, const bool &gradients=false) const =0
double m_zeta
Distance of model from Psf.
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1358
const GSkyDir & dir(void) const
Return position of radial spatial model.
Kernel for Psf phi angle integration used for stacked analysis.
GVector tan(const GVector &vector)
Computes tangens of vector elements.
Definition: GVector.cpp:1379
cta_psf_radial_kerns_delta(const GCTAResponseCube *rsp, const GModelSpatialRadial *model, const GSkyDir &obsDir, const GEnergies &srcEngs, const double &zeta, const double &theta_max, const int &iter, const bool &grad)
Kernel constructor for PSF integration of radial model.
const GCTACubePsf & psf(void) const
Return cube analysis point spread function.
Integration class for set of functions interface definition.
int size(void) const
Return number of energies in container.
Definition: GEnergies.hpp:170
Implementation of CTA helper classes for stacked vector response.
GModelPar * m_par_lon
Longitude parameter.
double m_cos_delta_cos_zeta
cos(delta) * cos(zeta)
const double & l(void) const
Return galactic longitude in radians.
Definition: GSkyDir.hpp:160
const double & dec(void) const
Return Declination in radians.
Definition: GSkyDir.hpp:241
GEnergies m_srcEngs
True photon energies.
Kernel for radial spatial model PSF phi angle integration.
double m_sin_delta_sin_zeta
sin(delta) * sin(zeta)
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
void fixed_iter(const int &iter)
Set fixed number of iterations.
Definition: GIntegral.hpp:183
int size(void) const
Return size of vector kernel for PSF integration of radial model.
Abstract radial spatial model base class.
Vector class.
Definition: GVector.hpp:46
Integration class for set of functions.
Definition: GIntegrals.hpp:47
double m_cos_theta_max
Cosine of m_theta_max.
Sky direction class.
Definition: GSkyDir.hpp:62
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
const GModelSpatialRadial * m_model
Radial model.
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489