GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAResponse_helpers.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAResponse_helpers.cpp - CTA response helper classes *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2012-2020 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 GCTAResponse_helpers.cpp
23  * @brief CTA response helper classes implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <cmath>
32 #include "GTools.hpp"
33 #include "GMath.hpp"
34 #include "GIntegral.hpp"
35 #include "GVector.hpp"
36 #include "GCTAResponse_helpers.hpp"
37 #include "GCTAEdisp.hpp"
38 #include "GCTASupport.hpp"
42 #include "GWcs.hpp"
43 
44 /* __ Method name definitions ____________________________________________ */
45 
46 /* __ Macros _____________________________________________________________ */
47 
48 /* __ Coding definitions _________________________________________________ */
49 //#define G_USE_OBSDIR_FOR_AEFF //!< Use event offset for Aeff computation
50 
51 /* __ Debug definitions __________________________________________________ */
52 //#define G_DEBUG_INTEGRAL //!< Debug integration
53 //#define G_DEBUG_MODEL_ZERO //!< Debug check for zero model values
54 
55 /* __ Constants __________________________________________________________ */
56 const double g_kludge_radius = 1.0e-12; //!< Tiny angle (radians)
57 const double g_ellipse_kludge_radius = 1.0e-6; //!< About 0.2 arc seconds
58  // A larger radius is used
59  // as the ellipse is defined
60  // in cartesian coordinates
61  // while distances are
62  // computed in spherical
63  // trigonometry
64 
65 
66 /*==========================================================================
67  = =
68  = Helper functions =
69  = =
70  ==========================================================================*/
71 
72 /***********************************************************************//**
73  * @brief Limit omega interval
74  *
75  * @param[in] min Interval minimum (radians).
76  * @param[in] max Interval maximum (radians).
77  * @param[in] domega Half length of interval (radians).
78  * @return Vector of intervals.
79  *
80  * Limits an omega interval [@p min,@p max] to the interval specified by
81  * [-@p domega,@p domega]. This may lead to a split of [@p min,@p max] in
82  * several intervals. The method thus returns a vector of intervals that
83  * overlap with [-@p domega,@p domega]. If there is no overlap with the
84  * interval, the method returns an empty vector.
85  *
86  * The method takes care of wrap arounds. It is assumed that on input
87  * [@p min,@p max] is contained within [-2pi,+2pi].
88  ***************************************************************************/
90  const double& max,
91  const double& domega)
92 {
93  // Allocate intervals
94  cta_omega_intervals intervals;
95 
96  // Continue only if domega is smaller than pi
97  if (domega < gammalib::pi) {
98 
99  // Set limiting intervals. To take care of a possible wrap around in
100  // omega we consider also intervals that are shifted by +/-2pi
101  double omega_min = -domega;
102  double omega_max = +domega;
103  double omega_min_plus = omega_min + gammalib::twopi;
104  double omega_max_plus = omega_max + gammalib::twopi;
105  double omega_min_minus = omega_min - gammalib::twopi;
106  double omega_max_minus = omega_max - gammalib::twopi;
107 
108  // If the [min,max] interval overlaps with the unshifted
109  // [-domega,domega] interval then constrain the interval
110  if (max > omega_min && min < omega_max) {
111  double interval_min = min;
112  double interval_max = max;
113  if (interval_min < omega_min) {
114  interval_min = omega_min;
115  }
116  if (interval_max > omega_max) {
117  interval_max = omega_max;
118  }
119  intervals.push_back(std::make_pair(interval_min,interval_max));
120  }
121 
122  // If the [min,max] interval overlaps with the [-domega,domega]
123  // interval shifted by +2pi then constrain the interval using the
124  // shifted interval
125  if (max > omega_min_plus && min < omega_max_plus) {
126  double interval_min = min;
127  double interval_max = max;
128  if (interval_min < omega_min_plus) {
129  interval_min = omega_min_plus;
130  }
131  if (interval_max > omega_max_plus) {
132  interval_max = omega_max_plus;
133  }
134  intervals.push_back(std::make_pair(interval_min,interval_max));
135  }
136 
137  // If the [min,max] interval overlaps with the [-domega,domega]
138  // interval shifted by -2pi then constrain the interval using the
139  // shifted interval
140  if (max > omega_min_minus && min < omega_max_minus) {
141  double interval_min = min;
142  double interval_max = max;
143  if (interval_min < omega_min_minus) {
144  interval_min = omega_min_minus;
145  }
146  if (interval_max > omega_max_minus) {
147  interval_max = omega_max_minus;
148  }
149  intervals.push_back(std::make_pair(interval_min,interval_max));
150  }
151 
152  } // endif: interval was not the full circle
153 
154  // ... otherwise append the provided interval
155  else {
156  intervals.push_back(std::make_pair(min,max));
157  }
158 
159  // Return intervals
160  return intervals;
161 }
162 
163 
164 /***********************************************************************//**
165  * @brief Determine resolution of spatial model
166  *
167  * @param[in] model Pointer to spatial model.
168  * @return Resolution of spatial model (radians).
169  *
170  * Determine the resolution of a spatial model. So far the method only works
171  * for a spatial map or cube model holding a WCS projection. If a constant
172  * spatial model is encountered a resolution of 180 deg is returned.
173  *
174  * If the resolution of the model could not be determined, the method returns
175  * a resolution of 0.01 deg.
176  ***************************************************************************/
178 {
179  // Initialise resolution to default resolution
180  double resolution = 0.01 * gammalib::deg2rad;
181 
182  // Extract pointer to spatial map. The pointer will be NULL if no spatial
183  // map exists.
184  const GSkyMap* map = NULL;
185  const GModelSpatialDiffuseMap* pmap =
186  dynamic_cast<const GModelSpatialDiffuseMap*>(model);
187  if (pmap != NULL) {
188  map = &(pmap->map());
189  }
190  else {
191  const GModelSpatialDiffuseCube* pcube =
192  dynamic_cast<const GModelSpatialDiffuseCube*>(model);
193  if (pcube != NULL) {
194  map = &(pcube->cube());
195  }
196  else {
197  const GModelSpatialDiffuseConst* pconst =
198  dynamic_cast<const GModelSpatialDiffuseConst*>(model);
199  if (pconst != NULL) {
200  resolution = gammalib::pi;
201  }
202  }
203  }
204 
205  // If a spatial map exists then get it's resolution. This so far only
206  // works for WCS maps.
207  if (map != NULL) {
208  const GWcs* wcs = dynamic_cast<const GWcs*>(map->projection());
209  if (wcs != NULL) {
210  double dx = std::abs(wcs->cdelt(0));
211  double dy = std::abs(wcs->cdelt(1));
212  resolution = ((dx < dy) ? dx : dy) * gammalib::deg2rad;
213  }
214  }
215 
216  // Return resolution
217  return resolution;
218 }
219 
220 
221 /*==========================================================================
222  = =
223  = Helper class methods for response computation =
224  = =
225  ==========================================================================*/
226 
227 /***********************************************************************//**
228  * @brief Integration kernel for npsf() method
229  *
230  * @param[in] delta Distance from PSF centre (radians).
231  * @return Azimuthally integrated PSF.
232  *
233  * Computes
234  *
235  * \f[
236  * \int_{0}^{\phi} PSF(\delta) d\phi
237  * \f]
238  *
239  * for a given offset angle \f$\delta\f$.
240  *
241  * The azimuthal integration is performed over an arclength given by
242  * \f$\phi\f$. The method actually assumes that the PSF is azimuthally
243  * symmetric, hence it just multiplies the PSF value by the arclength times
244  * the sinus of the offset angle.
245  ***************************************************************************/
246 double cta_npsf_kern_rad_azsym::eval(const double& delta)
247 {
248  // Initialise PSF value
249  double value = 0.0;
250 
251  // Get arclength for given radius in radians
252  double phi = gammalib::roi_arclength(delta,
253  m_psf,
254  m_cospsf,
255  m_sinpsf,
256  m_roi,
257  m_cosroi);
258 
259  // If arclength is positive then compute the PSF value
260  if (phi > 0) {
261 
262  // Compute PSF value
263  value = m_rsp->psf(delta, m_theta, m_phi, m_zenith, m_azimuth, m_logE) *
264  phi * std::sin(delta);
265 
266  // Compile option: Check for NaN/Inf
267  #if defined(G_NAN_CHECK)
268  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
269  std::cout << "*** ERROR: cta_npsf_kern_rad_azsym::eval";
270  std::cout << " NaN/Inf encountered";
271  std::cout << " (value=" << value;
272  std::cout << ", delta=" << delta;
273  std::cout << ", phi=" << phi << ")";
274  std::cout << std::endl;
275  }
276  #endif
277 
278  } // endif: arclength was positive
279 
280  // Return
281  return value;
282 }
283 
284 
285 /***********************************************************************//**
286  * @brief Integration kernel for GCTAResponseIrf::nroi method
287  *
288  * @param[in] etrue True photon energy in MeV.
289  * @return Nroi.
290  ***************************************************************************/
291 double cta_nroi_kern::eval(const double& etrue)
292 {
293  // Set true energy
294  GEnergy srcEng;
295  srcEng.MeV(etrue);
296 
297  // Compute response components
298  double nroi_spatial = m_rsp->nroi(*m_model, srcEng, m_srcTime, m_obsEng, m_obsTime, *m_obs);
299  double nroi_spectral = m_model->spectral()->eval(srcEng, m_srcTime);
300  double nroi_temporal = m_model->temporal()->eval(m_srcTime);
301 
302  // Compute response
303  double nroi = nroi_spatial * nroi_spectral * nroi_temporal;
304 
305  // Return response
306  return nroi;
307 }
308 
309 
310 /*==========================================================================
311  = =
312  = Helper class methods for radial models =
313  = =
314  ==========================================================================*/
315 
316 /***********************************************************************//**
317  * @brief Kernel for radial model zenith angle integration of Irf
318  *
319  * @param[in] rho Zenith angle with respect to model centre [radians].
320  *
321  * Computes the kernel
322  *
323  * \f[
324  * K(\rho | E, t) = \sin \rho \times S_{\rm p}(\rho | E, t) \times
325  * \int_{\omega_{\rm min}}^{\omega_{\rm max}}
326  * IRF(\rho, \omega) d\omega
327  * \f]
328  *
329  * for the zenith angle integration of radial models.
330  ***************************************************************************/
331 double cta_irf_radial_kern_rho::eval(const double& rho)
332 {
333  // Initialise result
334  double irf = 0.0;
335 
336  // Continue only if rho is positive (otherwise the integral will be
337  // zero)
338  if (rho > 0.0) {
339 
340  // Compute half length of arc that lies within PSF validity circle
341  // (in radians)
342  double domega = 0.5 * gammalib::roi_arclength(rho,
343  m_zeta,
344  m_cos_zeta,
345  m_sin_zeta,
346  m_delta_max,
348 
349  // Continue only if arc length is positive
350  if (domega > 0.0) {
351 
352  // Compute omega integration range
353  double omega_min = -domega;
354  double omega_max = +domega;
355 
356  // Reduce rho by an infinite amount to avoid rounding errors
357  // at the boundary of a sharp edged model
358  double rho_kludge = rho - g_kludge_radius;
359  if (rho_kludge < 0.0) {
360  rho_kludge = 0.0;
361  }
362 
363  // Evaluate sky model
364  double model = m_model->eval(rho_kludge, m_srcEng, m_srcTime);
365 
366  // Debug: test if model is non positive
367  #if defined(G_DEBUG_MODEL_ZERO)
368  if (model <= 0.0) {
369  std::cout << "*** WARNING: cta_irf_radial_kern_rho::eval";
370  std::cout << " zero model for (rho)=(";
371  std::cout << rho*gammalib::rad2deg << ")";
372  std::cout << " rho-r_model=" << (rho-m_model.theta_max());
373  std::cout << " radians" << std::endl;
374  }
375  #endif
376 
377  // Continue only if model is positive
378  if (model > 0.0) {
379 
380  // Precompute cosine and sine terms for azimuthal
381  // integration
382  double cos_rho = std::cos(rho);
383  double sin_rho = std::sin(rho);
384  double cos_psf = cos_rho*m_cos_zeta;
385  double sin_psf = sin_rho*m_sin_zeta;
386  double cos_ph = cos_rho*m_cos_lambda;
387  double sin_ph = sin_rho*m_sin_lambda;
388 
389  // Setup integration kernel
391  m_zenith,
392  m_azimuth,
393  m_srcEng,
394  m_obsEng,
395  m_zeta,
396  m_lambda,
397  m_omega0,
398  rho,
399  cos_psf,
400  sin_psf,
401  cos_ph,
402  sin_ph);
403 
404  // Integrate over phi
405  GIntegral integral(&integrand);
406  integral.fixed_iter(m_iter);
407  irf = integral.romberg(omega_min, omega_max, m_iter) *
408  model * sin_rho;
409 
410  // Compile option: Check for NaN/Inf
411  #if defined(G_NAN_CHECK)
413  std::cout << "*** ERROR: cta_irf_radial_kern_rho";
414  std::cout << "(rho=" << rho << "):";
415  std::cout << " NaN/Inf encountered";
416  std::cout << " (irf=" << irf;
417  std::cout << ", domega=" << domega;
418  std::cout << ", model=" << model;
419  std::cout << ", sin_rho=" << sin_rho << ")";
420  std::cout << std::endl;
421  }
422  #endif
423 
424  } // endif: model was positive
425 
426  } // endif: arclength was positive
427 
428  } // endif: rho was positive
429 
430  // Return result
431  return irf;
432 }
433 
434 
435 /***********************************************************************//**
436  * @brief Kernel for radial model azimuth angle IRF integration
437  *
438  * @param[in] omega Azimuth angle (radians).
439  *
440  * Computes the kernel
441  *
442  * \f[
443  * IRF(\rho,\omega)
444  * \f]
445  *
446  * for the azimuth angle integration of radial models.
447  *
448  * From the model coordinates \f$(\rho,\omega)\f$, the method computes the
449  * angle between the true (\f$\vec{p}\f$) and observed (\f$\vec{p'}\f$)
450  * photon arrival direction using
451  *
452  * \f[\delta = \arccos(\cos \rho \cos \zeta +
453  * \sin \rho \sin \zeta \cos \omega)\f]
454  *
455  * where
456  * \f$\zeta\f$ is the angular distance between the observed photon arrival
457  * direction \f$\vec{p'}\f$ and the model centre \f$\vec{m}\f$. This angle
458  * \f$\delta\f$ is used to compute the \f$PSF(\delta)\f$ value.
459  *
460  * The method computes also the angle \f$\theta\f$ between the observed
461  * photon arrival direction \f$\vec{p'}\f$ and the camera pointing
462  * \f$\vec{d}\f$ using
463  *
464  * \f[\theta = \arccos(\cos \rho \cos \lambda +
465  * \sin \rho \sin \lambda \cos \omega_0 - \omega)\f]
466  *
467  * where
468  * \f$\lambda\f$ is the angular distance between the model centre
469  * \f$\vec{m}\f$ and the camera pointing direction \f$\vec{d}\f$.
470  * The angle \f$\theta\f$ is used in the computation of the IRF (no
471  * azimuthal dependence is so far implemented for the IRF computation).
472  ***************************************************************************/
473 double cta_irf_radial_kern_omega::eval(const double& omega)
474 {
475  // Get log10(E/TeV) of true photon energy
476  double srcLogEng = m_srcEng.log10TeV();
477 
478  // Compute PSF offset angle [radians]
479  double delta = std::acos(m_cos_psf + m_sin_psf * std::cos(omega));
480 
481  // Compute true photon offset angle in camera system [radians]
482  double offset = std::acos(m_cos_ph + m_sin_ph * std::cos(m_omega0 - omega));
483 
484  //TODO: Compute true photon azimuth angle in camera system [radians]
485  double azimuth = 0.0;
486 
487  // Evaluate IRF
488  double irf = m_rsp->aeff(offset, azimuth, m_zenith, m_azimuth, srcLogEng) *
489  m_rsp->psf(delta, offset, azimuth, m_zenith, m_azimuth, srcLogEng);
490 
491  // Optionally take energy dispersion into account
492  if (m_rsp->use_edisp() && irf > 0.0) {
493  irf *= m_rsp->edisp(m_obsEng, m_srcEng, offset, azimuth, m_zenith, m_azimuth);
494  }
495 
496  // Compile option: Check for NaN/Inf
497  #if defined(G_NAN_CHECK)
499  std::cout << "*** ERROR: cta_irf_radial_kern_omega::eval";
500  std::cout << "(omega=" << omega << "):";
501  std::cout << " NaN/Inf encountered";
502  std::cout << " (irf=" << irf;
503  std::cout << ", delta=" << delta;
504  std::cout << ", offset=" << offset;
505  std::cout << ", azimuth=" << azimuth << ")";
506  std::cout << std::endl;
507  }
508  #endif
509 
510  // Return
511  return irf;
512 }
513 
514 
515 /***********************************************************************//**
516  * @brief Kernel for zenith angle Nroi integration or radial model
517  *
518  * @param[in] rho Radial model zenith angle (radians).
519  *
520  * Computes
521  *
522  * \f[
523  * K(\rho | E, t) = \sin \rho \times S_{\rm p}(\rho | E, t) \times
524  * \int_{\omega_{\rm min}}^{\omega_{\rm max}}
525  * N_{\rm pred}(\rho,\omega) d\omega
526  * \f]
527  *
528  * The azimuth angle integration range
529  * \f$[\omega_{\rm min}, \omega_{\rm max}]\f$
530  * is limited to an arc around the vector connecting the model centre to
531  * the ROI centre. This limitation assures that the integration converges
532  * properly.
533  ***************************************************************************/
534 double cta_nroi_radial_kern_rho::eval(const double& rho)
535 {
536  // Initialise Nroi value
537  double nroi = 0.0;
538 
539  // Continue only if rho is positive
540  if (rho > 0.0) {
541 
542  // Compute half length of arc that lies within ROI+PSF radius (radians)
543  double domega = 0.5 * gammalib::roi_arclength(rho,
544  m_dist,
545  m_cos_dist,
546  m_sin_dist,
547  m_radius,
548  m_cos_radius);
549 
550  // Continue only if arc length is positive
551  if (domega > 0.0) {
552 
553  // Compute omega integration range
554  double omega_min = m_omega0 - domega;
555  double omega_max = m_omega0 + domega;
556 
557  // Reduce rho by an infinite amount to avoid rounding errors
558  // at the boundary of a sharp edged model
559  double rho_kludge = rho - g_kludge_radius;
560  if (rho_kludge < 0.0) {
561  rho_kludge = 0.0;
562  }
563 
564  // Get radial model value
565  double model = m_model->eval(rho_kludge, m_srcEng, m_srcTime);
566 
567  // Debug: test if model is non positive
568  #if defined(G_DEBUG_MODEL_ZERO)
569  if (model <= 0.0) {
570  std::cout << "*** WARNING: cta_nroi_radial_kern_rho::eval";
571  std::cout << " zero model for (rho)=(";
572  std::cout << rho*gammalib::rad2deg << ")";
573  std::cout << " rho-r_model=" << (rho-m_model.theta_max());
574  std::cout << " radians" << std::endl;
575  }
576  #endif
577 
578  // Continue only if we have a positive model value
579  if (model > 0.0) {
580 
581  // Compute sine and cosine of offset angle
582  double sin_rho = std::sin(rho);
583  double cos_rho = std::cos(rho);
584 
585  // Setup phi integration kernel
587  m_obs,
588  m_rot,
589  m_srcEng,
590  m_srcTime,
591  m_obsEng,
592  m_obsTime,
593  sin_rho,
594  cos_rho);
595 
596  // Integrate over phi
597  GIntegral integral(&integrand);
598  integral.fixed_iter(m_iter);
599  nroi = integral.romberg(omega_min, omega_max, m_iter) *
600  sin_rho * model;
601 
602  // Debug: Check for NaN
603  #if defined(G_NAN_CHECK)
604  if (gammalib::is_notanumber(nroi) || gammalib::is_infinite(nroi)) {
605  std::cout << "*** ERROR: cta_nroi_radial_kern_rho::eval";
606  std::cout << "(rho=" << rho << "):";
607  std::cout << " NaN/Inf encountered";
608  std::cout << " (nroi=" << nroi;
609  std::cout << ", model=" << model;
610  std::cout << ", omega=[" << omega_min << "," << omega_max << "]";
611  std::cout << ", sin_rho=" << sin_rho;
612  std::cout << ")" << std::endl;
613  }
614  #endif
615 
616  } // endif: model was positive
617 
618  } // endif: arc length was positive
619 
620  } // endif: rho was positive
621 
622  // Return Nroi
623  return nroi;
624 }
625 
626 
627 /***********************************************************************//**
628  * @brief Kernel for azimuth angle Nroi integration of radial model
629  *
630  * @param[in] omega Azimuth angle (radians).
631  *
632  * @todo Re-consider formula for possible simplification (dumb matrix
633  * multiplication is definitely not the fastest way to do that
634  * computation).
635  ***************************************************************************/
636 double cta_nroi_radial_kern_omega::eval(const double& omega)
637 {
638  // Compute sky direction vector in native coordinates
639  double cos_omega = std::cos(omega);
640  double sin_omega = std::sin(omega);
641  GVector native(-cos_omega*m_sin_rho, sin_omega*m_sin_rho, m_cos_rho);
642 
643  // Rotate from native into celestial system
644  GVector cel = *m_rot * native;
645 
646  // Set sky direction
647  GSkyDir srcDir;
648  srcDir.celvector(cel);
649 
650  // Set Photon
651  GPhoton photon(srcDir, m_srcEng, m_srcTime);
652 
653  // Compute point source Nroi for this sky direction
654  double nroi = m_rsp->nirf(photon, m_obsEng, m_obsTime, *m_obs);
655 
656  // Debug: Check for NaN
657  #if defined(G_NAN_CHECK)
658  if (gammalib::is_notanumber(nroi) || gammalib::is_infinite(nroi)) {
659  std::cout << "*** ERROR: cta_nroi_radial_kern_omega::eval";
660  std::cout << "(omega=" << omega << "):";
661  std::cout << " NaN/Inf encountered";
662  std::cout << " (nroi=" << nroi;
663  std::cout << ", cos_omega=" << cos_omega;
664  std::cout << ", sin_omega=" << sin_omega;
665  std::cout << ")" << std::endl;
666  }
667  #endif
668 
669  // Return Nroi
670  return nroi;
671 }
672 
673 
674 /*==========================================================================
675  = =
676  = Helper class methods for elliptical models =
677  = =
678  ==========================================================================*/
679 
680 /***********************************************************************//**
681  * @brief Kernel for elliptical model integration over model's zenith angle
682  *
683  * @param[in] rho Radial distance from model centre [radians].
684  * @return Radial IRF integration kernel.
685  *
686  * Computes
687  *
688  * \f[
689  * K(\rho | E, t) = \sin \rho \times
690  * \int_{\omega_{\rm min}}^{\omega_{\rm max}}
691  * S_{\rm p}(\rho, \omega | E, t) \, IRF(\rho, \omega)
692  * d\omega
693  * \f]
694  *
695  * where
696  * \f$\omega\f$ is the azimuth angle with respect to the model centre,
697  * counted counterclockwise from the vector connecting the model centre
698  * to the observed photon direction, and
699  * \f$\rho\f$ is the radial distance from the model centre.
700  * \f$S_{\rm p}(\rho, \omega | E, t)\f$ is the elliptical source model
701  * for a given true photon energy and photon arrival time,
702  * \f$IRF(\rho, \omega)\f$ is the instrument response function.
703  *
704  * The method performs the required coordinate transformations from the
705  * model system, spanned by \f$(\rho, \omega)\f$, to the system needed for
706  * IRF computations. Furthermore, the method limits the integration range
707  * to area where the ellipse intersects the IRF.
708  ***************************************************************************/
709 double cta_irf_elliptical_kern_rho::eval(const double& rho)
710 {
711  // Initialise result
712  double irf = 0.0;
713 
714  // Continue only if rho is positive
715  if (rho > 0.0) {
716 
717  // Compute half length of the arc (in radians) from a circle with
718  // radius rho that intersects with the point spread function, defined
719  // as a circle with maximum radius m_delta_max
720  double domega = 0.5 * gammalib::roi_arclength(rho,
721  m_rho_obs,
724  m_delta_max,
726 
727  // Continue only if arc length is positive
728  if (domega > 0.0) {
729 
730  // Precompute cosine and sine terms for azimuthal integration
731  double cos_rho = std::cos(rho);
732  double sin_rho = std::sin(rho);
733  double cos_psf = cos_rho * m_cos_rho_obs;
734  double sin_psf = sin_rho * m_sin_rho_obs;
735  double cos_ph = cos_rho * m_cos_rho_pnt;
736  double sin_ph = sin_rho * m_sin_rho_pnt;
737 
738  // Reduce rho by an infinite amount to avoid rounding errors
739  // at the boundary of a sharp edged model (e.g. an elliptical
740  // disk model)
741  double rho_kludge = rho - g_ellipse_kludge_radius;
742  if (rho_kludge < 0.0) {
743  rho_kludge = 0.0;
744  }
745 
746  // Setup integration kernel
748  m_model,
749  m_zenith,
750  m_azimuth,
751  m_srcEng,
752  m_srcTime,
753  m_obsEng,
755  m_omega_pnt,
756  rho_kludge,
757  cos_psf,
758  sin_psf,
759  cos_ph,
760  sin_ph);
761 
762  // Setup integrator
763  GIntegral integral(&integrand);
764  integral.fixed_iter(m_iter);
765 
766  // If the radius rho is not larger than the semiminor axis
767  // boundary, the circle with that radius is fully contained in
768  // the ellipse and we can just integrate over the relevant arc
769  if (rho <= m_semiminor) {
770 
771  // Compute omega integration range
772  double omega_min = -domega;
773  double omega_max = +domega;
774 
775  // Integrate over omega
776  irf = integral.romberg(omega_min, omega_max, m_iter) *
777  sin_rho;
778 
779  } // endif: circle comprised in ellipse
780 
781  // ... otherwise there are arcs that intersect with the Psf circle
782  else {
783 
784  // Compute half the arc length (in radians) of a circle of
785  // radius rho, centred on the model, that intersects with
786  // the ellipse boundary
787  double arg1 = 1.0 - (m_semiminor*m_semiminor) / (rho*rho);
788  double arg2 = 1.0 - (m_semiminor*m_semiminor) /
790  double omega_width = std::acos(std::sqrt(arg1/arg2));
791 
792  // Continue only if the arclength is positive
793  if (omega_width > 0.0) {
794 
795  // Compute azimuth angle difference between ellipse
796  // position angle and position angle of observed
797  // photon in the model system. This angle will define
798  // the reference point around which the circle arc. Make
799  // sure that omega_0 is within [-pi,pi] thus that the omega
800  // intervals are between [-2pi,2pi]
801  double omega_0 = m_posangle - m_posangle_obs;
802  if (omega_0 > gammalib::pi) {
803  omega_0 -= gammalib::pi;
804  }
805  else if (omega_0 < -gammalib::pi) {
806  omega_0 += gammalib::pi;
807  }
808 
809  // Compute azimuth angle intervals
810  double omega1_min = omega_0 - omega_width;
811  double omega1_max = omega_0 + omega_width;
812  double omega2_min = omega1_min + gammalib::pi;
813  double omega2_max = omega1_max + gammalib::pi;
814 
815  // Limit intervals to the intersection of the ellipse with
816  // the Psf circle. This may lead to a split of intervals,
817  // and we gather all these intervals in a special interval
818  // pair containers
819  cta_omega_intervals intervals1 =
820  gammalib::limit_omega(omega1_min, omega1_max, domega);
821  cta_omega_intervals intervals2 =
822  gammalib::limit_omega(omega2_min, omega2_max, domega);
823 
824  // Integrate over all intervals for omega1
825  for (int i = 0; i < intervals1.size(); ++i) {
826  double min = intervals1[i].first;
827  double max = intervals1[i].second;
828  irf += integral.romberg(min, max, m_iter) * sin_rho;
829  }
830 
831  // Integrate over all intervals for omega2
832  for (int i = 0; i < intervals2.size(); ++i) {
833  double min = intervals2[i].first;
834  double max = intervals2[i].second;
835  irf += integral.romberg(min, max, m_iter) * sin_rho;
836  }
837 
838  } // endif: arc length was positive
839 
840  } // endelse: circle was not comprised in ellipse
841 
842  // Compile option: Check for NaN/Inf
843  #if defined(G_NAN_CHECK)
845  std::cout << "*** ERROR: cta_irf_elliptical_kern_rho";
846  std::cout << "(rho=" << rho << "):";
847  std::cout << " NaN/Inf encountered";
848  std::cout << " (irf=" << irf;
849  std::cout << ", domega=" << domega;
850  std::cout << ", sin_rho=" << sin_rho << ")";
851  std::cout << std::endl;
852  }
853  #endif
854 
855  } // endif: arc length was positive
856 
857  } // endif: rho was positive
858 
859  // Return result
860  return irf;
861 }
862 
863 
864 /***********************************************************************//**
865  * @brief Kernel for elliptical model integration over model's azimuth angle
866  *
867  * @param[in] omega Azimuth angle (radians).
868  *
869  * Computes
870  *
871  * \f[
872  * S_{\rm p}(\omega | \rho, E, t) \, IRF(\omega | \rho)
873  * \f]
874  *
875  * where
876  * \f$\omega\f$ is the azimuth angle with respect to the model centre,
877  * counted counterclockwise from the vector connecting the model centre
878  * to the observed photon direction, and
879  * \f$\rho\f$ is the radial distance from the model centre.
880  *
881  * From the coordinates \f$(\rho,\omega)\f$ in the model system, the method
882  * computes the angle between the true (\f$\vec{p}\f$) and observed
883  * (\f$\vec{p'}\f$) photon arrival direction using
884  *
885  * \f[\delta = \arccos(\cos \rho \cos \zeta +
886  * \sin \rho \sin \zeta \cos \omega)\f]
887  *
888  * where
889  * \f$\zeta\f$ is the angular distance between the observed photon arrival
890  * direction \f$\vec{p'}\f$ and the model centre \f$\vec{m}\f$. This angle
891  * \f$\delta\f$ is used to compute the \f$PSF(\delta)\f$ value.
892  *
893  * The method computes also the angle \f$\theta\f$ between the observed
894  * photon arrival direction \f$\vec{p'}\f$ and the camera pointing
895  * \f$\vec{d}\f$ using
896  *
897  * \f[\theta = \arccos(\cos \rho \cos \lambda +
898  * \sin \rho \sin \lambda \cos \omega_0 - \omega)\f]
899  *
900  * where
901  * \f$\lambda\f$ is the angular distance between the model centre
902  * \f$\vec{m}\f$ and the camera pointing direction \f$\vec{d}\f$.
903  * The angle \f$\theta\f$ is used in the computation of the IRF (no
904  * azimuthal dependence is so far implemented for the IRF computation).
905  ***************************************************************************/
906 double cta_irf_elliptical_kern_omega::eval(const double& omega)
907 {
908  // Initialise IRF value
909  double irf = 0.0;
910 
911  // Compute azimuth angle in model coordinate system (radians)
912  double omega_model = omega + m_posangle_obs;
913 
914  // Evaluate sky model
915  double model = m_model->eval(m_rho, omega_model, m_srcEng, m_srcTime);
916 
917  // Debug: test if model is non positive
918  #if defined(G_DEBUG_MODEL_ZERO)
919  if (model <= 0.0) {
920  double m_semiminor_rad = m_model->semiminor() * gammalib::deg2rad;
921  double m_semimajor_rad = m_model->semimajor() * gammalib::deg2rad;
922  double diff_angle = omega_model - m_model->posangle() * gammalib::deg2rad;
923  double cosinus = std::cos(diff_angle);
924  double sinus = std::sin(diff_angle);
925  double arg1 = m_semiminor_rad * cosinus;
926  double arg2 = m_semimajor_rad * sinus;
927  double r_ellipse = m_semiminor_rad * m_semimajor_rad /
928  std::sqrt(arg1*arg1 + arg2*arg2);
929  std::cout << "*** WARNING: cta_irf_elliptical_kern_omega::eval:";
930  std::cout << " zero model for (rho,omega)=(";
931  std::cout << m_rho*gammalib::rad2deg << ",";
932  std::cout << omega*gammalib::rad2deg << ")";
933  std::cout << " semiminor=" << m_semiminor_rad;
934  std::cout << " semimajor=" << m_semimajor_rad;
935  std::cout << " posang=" << m_model->posangle() * gammalib::deg2rad;
936  std::cout << " rel_posang=" << diff_angle * gammalib::deg2rad;
937  std::cout << " rho/r_ellipse=" << (m_rho/r_ellipse);
938  std::cout << std::endl;
939  }
940  #endif
941 
942  // Continue only if model is positive
943  if (model > 0.0) {
944 
945  // Get log10(E/TeV) of true photon energy
946  double srcLogEng = m_srcEng.log10TeV();
947 
948  // Compute Psf offset angle [radians]
949  double delta = std::acos(m_cos_psf + m_sin_psf * std::cos(omega));
950 
951  // Compute true photon offset and azimuth angle in camera system
952  // [radians]
953  double theta = std::acos(m_cos_ph + m_sin_ph * std::cos(m_omega_pnt - omega));
954  double phi = 0.0; //TODO: Implement IRF Phi dependence
955 
956  // Evaluate IRF * model
957  irf = m_rsp->aeff(theta, phi, m_zenith, m_azimuth, srcLogEng) *
958  m_rsp->psf(delta, theta, phi, m_zenith, m_azimuth, srcLogEng) *
959  model;
960 
961  // Optionally take energy dispersion into account
962  if (m_rsp->use_edisp() && irf > 0.0) {
963  irf *= m_rsp->edisp(m_obsEng, m_srcEng, theta, phi, m_zenith, m_azimuth);
964  }
965 
966  // Compile option: Check for NaN/Inf
967  #if defined(G_NAN_CHECK)
969  std::cout << "*** ERROR: cta_irf_elliptical_kern_omega::eval";
970  std::cout << "(omega=" << omega << "):";
971  std::cout << " NaN/Inf encountered";
972  std::cout << " (irf=" << irf;
973  std::cout << ", model=" << model;
974  std::cout << ", delta=" << delta;
975  std::cout << ", theta=" << theta;
976  std::cout << ", phi=" << phi << ")";
977  std::cout << std::endl;
978  }
979  #endif
980 
981  } // endif: model is positive
982 
983  // Return
984  return irf;
985 }
986 
987 
988 /***********************************************************************//**
989  * @brief Kernel for zenith angle Nroi integration of elliptical model
990  *
991  * @param[in] rho Elliptical model offset angle (radians).
992  *
993  * Computes
994  *
995  * \f[
996  * K(\rho | E, t) = \sin \rho \times
997  * \int_{\omega_{\rm min}}^{\omega_{\rm max}}
998  * S_{\rm p}(\rho,\omega | E, t) \,
999  * N_{\rm pred}(\rho,\omega) d\omega
1000  * \f]
1001  *
1002  * where
1003  * \f$\omega\f$ is the azimuth angle with respect to the model centre,
1004  * counted counterclockwise from the vector connecting the model centre
1005  * to the ROI centre, and
1006  * \f$\rho\f$ is the radial distance from the model centre.
1007  * \f$S_{\rm p}(\rho, \omega | E, t)\f$ is the elliptical source model
1008  * for a given true photon energy and photon arrival time,
1009  * \f$N_{\rm pred}(\rho,\omega)\f$ is the data space integral over the
1010  * response function.
1011  *
1012  * The method performs the required coordinate transformations from the
1013  * model system, spanned by \f$(\rho, \omega)\f$, to the system needed for
1014  * Npred computations. Furthermore, the method limits the integration range
1015  * to area where the ellipse intersects the ROI.
1016  ***************************************************************************/
1017 double cta_nroi_elliptical_kern_rho::eval(const double& rho)
1018 {
1019  // Initialise Nroi value
1020  double nroi = 0.0;
1021 
1022  // Continue only if rho is positive
1023  if (rho > 0.0) {
1024 
1025  // Compute half length of arc that lies within ROI+PSF radius (radians)
1026  double domega = 0.5 * gammalib::roi_arclength(rho,
1027  m_rho_roi,
1028  m_cos_rho_roi,
1029  m_sin_rho_roi,
1030  m_radius_roi,
1032 
1033  // Continue only if arc length is positive
1034  if (domega > 0.0) {
1035 
1036  // Compute sine and cosine terms for azimuthal integration
1037  double sin_rho = std::sin(rho);
1038  double cos_rho = std::cos(rho);
1039 
1040  // Reduce rho by an infinite amount to avoid rounding errors
1041  // at the boundary of a sharp edged model (e.g. an elliptical
1042  // disk model)
1043  double rho_kludge = rho - g_ellipse_kludge_radius;
1044  if (rho_kludge < 0.0) {
1045  rho_kludge = 0.0;
1046  }
1047 
1048  // Setup phi integration kernel
1050  m_obs,
1051  m_model,
1052  m_rot,
1053  m_srcEng,
1054  m_srcTime,
1055  m_obsEng,
1056  m_obsTime,
1057  rho_kludge,
1058  sin_rho,
1059  cos_rho,
1060  m_posangle_roi);
1061 
1062  // Setup integrator
1063  GIntegral integral(&integrand);
1064  integral.fixed_iter(m_iter);
1065 
1066  // If the radius rho is not larger than the semiminor axis
1067  // boundary, the circle with that radius is fully contained in
1068  // the ellipse and we can just integrate over the relevant arc
1069  if (rho <= m_semiminor) {
1070 
1071  // Compute omega integration range
1072  double omega_min = -domega;
1073  double omega_max = +domega;
1074 
1075  // Integrate over omega
1076  nroi = integral.romberg(omega_min, omega_max, m_iter) *
1077  sin_rho;
1078 
1079  } // endif: circle comprised in ellipse
1080 
1081  // ... otherwise there are arcs that intersect with the ROI
1082  // circle
1083  else {
1084 
1085  // Compute half the arc length (in radians) of a circle of
1086  // radius rho, centred on the model, that intersects with
1087  // the ellipse boundary
1088  double arg1 = 1.0 - (m_semiminor*m_semiminor) / (rho*rho);
1089  double arg2 = 1.0 - (m_semiminor*m_semiminor) /
1091  double omega_width = std::acos(std::sqrt(arg1/arg2));
1092 
1093  // Continue only if the arclength is positive
1094  if (omega_width > 0.0) {
1095 
1096  // Compute azimuth angle difference between ellipse
1097  // position angle and position angle of ROI in the model
1098  // system. This angle will define the reference point for
1099  // the circle arcs. Make sure that omega_0 is within [-pi,pi]
1100  // thus that the omega intervals are between [-2pi,2pi]
1101  double omega_0 = m_posangle - m_posangle_roi;
1102  if (omega_0 > gammalib::pi) {
1103  omega_0 -= gammalib::pi;
1104  }
1105  else if (omega_0 < -gammalib::pi) {
1106  omega_0 += gammalib::pi;
1107  }
1108 
1109  // Compute azimuth angle intervals
1110  double omega1_min = omega_0 - omega_width;
1111  double omega1_max = omega_0 + omega_width;
1112  double omega2_min = omega1_min + gammalib::pi;
1113  double omega2_max = omega1_max + gammalib::pi;
1114 
1115  // Limit intervals to the intersection of the ellipse with
1116  // the Psf circle. This may lead to a split of intervals,
1117  // and we gather all these intervals in a special interval
1118  // pair containers
1119  cta_omega_intervals intervals1 =
1120  gammalib::limit_omega(omega1_min, omega1_max, domega);
1121  cta_omega_intervals intervals2 =
1122  gammalib::limit_omega(omega2_min, omega2_max, domega);
1123 
1124  // Integrate over all intervals for omega1
1125  for (int i = 0; i < intervals1.size(); ++i) {
1126  double min = intervals1[i].first;
1127  double max = intervals1[i].second;
1128  nroi += integral.romberg(min, max, m_iter) * sin_rho;
1129  }
1130 
1131  // Integrate over all intervals for omega2
1132  for (int i = 0; i < intervals2.size(); ++i) {
1133  double min = intervals2[i].first;
1134  double max = intervals2[i].second;
1135  nroi += integral.romberg(min, max, m_iter) * sin_rho;
1136  }
1137 
1138  } // endif: arc length was positive
1139 
1140  } // endelse: circle was not comprised in ellipse
1141 
1142  // Debug: Check for NaN
1143  #if defined(G_NAN_CHECK)
1144  if (gammalib::is_notanumber(nroi) || gammalib::is_infinite(nroi)) {
1145  std::cout << "*** ERROR: cta_nroi_elliptical_kern_rho::eval";
1146  std::cout << "(rho=" << rho << "):";
1147  std::cout << " NaN/Inf encountered";
1148  std::cout << " (nroi=" << nroi;
1149  std::cout << ", sin_rho=" << sin_rho;
1150  std::cout << ", cos_rho=" << cos_rho;
1151  std::cout << ")" << std::endl;
1152  }
1153  #endif
1154 
1155  } // endif: arc length was positive
1156 
1157  } // endif: rho was positive
1158 
1159  // Return Nroi
1160  return nroi;
1161 }
1162 
1163 
1164 /***********************************************************************//**
1165  * @brief Kernel for azimuth angle Nroi integration of elliptical model
1166  *
1167  * @param[in] omega Azimuth angle (radians).
1168  *
1169  * Computes
1170  *
1171  * \f[
1172  * S_{\rm p}(\omega | \rho, E, t) \, N_{\rm pred}(\omega | \rho)
1173  * \f]
1174  *
1175  * where
1176  * \f$\omega\f$ is the azimuth angle with respect to the model centre,
1177  * counted counterclockwise from the vector connecting the model centre
1178  * to the centre of the Region of Interest (ROI), and
1179  * \f$\rho\f$ is the radial distance from the model centre.
1180  *
1181  * @todo Npred computation goes over sky coordinates. This can maybe be
1182  * optimized to reduce the number of coordinate transformations.
1183  * @todo Check whether the Npred omega argument is the right one.
1184  ***************************************************************************/
1185 double cta_nroi_elliptical_kern_omega::eval(const double& omega)
1186 {
1187  // Initialise Nroi value
1188  double nroi = 0.0;
1189 
1190  // Compute azimuth angle in model coordinate system (radians)
1191  double omega_model = omega + m_posangle_roi;
1192 
1193  // Evaluate sky model
1194  double model = m_model->eval(m_rho, omega_model, m_srcEng, m_srcTime);
1195 
1196  // Debug: test if model is non positive
1197  #if defined(G_DEBUG_MODEL_ZERO)
1198  if (model <= 0.0) {
1199  double m_semiminor_rad = m_model->semiminor() * gammalib::deg2rad;
1200  double m_semimajor_rad = m_model->semimajor() * gammalib::deg2rad;
1201  double diff_angle = omega_model - m_model->posangle() * gammalib::deg2rad;
1202  double cosinus = std::cos(diff_angle);
1203  double sinus = std::sin(diff_angle);
1204  double arg1 = m_semiminor_rad * cosinus;
1205  double arg2 = m_semimajor_rad * sinus;
1206  double r_ellipse = m_semiminor_rad * m_semimajor_rad /
1207  std::sqrt(arg1*arg1 + arg2*arg2);
1208  std::cout << "*** WARNING: cta_nroi_elliptical_kern_omega::eval";
1209  std::cout << " zero model for (rho,omega)=(";
1210  std::cout << m_rho*gammalib::rad2deg << ",";
1211  std::cout << omega*gammalib::rad2deg << ")";
1212  std::cout << " rho-r_ellipse=" << (m_rho-r_ellipse) << " radians";
1213  std::cout << std::endl;
1214  }
1215  #endif
1216 
1217  // Continue only if model is positive
1218  if (model > 0.0) {
1219 
1220  // Compute sky direction vector in native coordinates
1221  double cos_omega = std::cos(omega_model);
1222  double sin_omega = std::sin(omega_model);
1223  GVector native(-cos_omega*m_sin_rho, sin_omega*m_sin_rho, m_cos_rho);
1224 
1225  // Rotate from native into celestial system
1226  GVector cel = *m_rot * native;
1227 
1228  // Set sky direction
1229  GSkyDir srcDir;
1230  srcDir.celvector(cel);
1231 
1232  // Set Photon
1233  GPhoton photon(srcDir, m_srcEng, m_srcTime);
1234 
1235  // Compute Nroi for this sky direction
1236  nroi = m_rsp->nirf(photon, m_obsEng, m_obsTime, *m_obs) * model;
1237 
1238  // Debug: Check for NaN
1239  #if defined(G_NAN_CHECK)
1240  if (gammalib::is_notanumber(nroi) || gammalib::is_infinite(nroi)) {
1241  std::cout << "*** ERROR: cta_nroi_elliptical_kern_omega::eval";
1242  std::cout << "(omega=" << omega << "):";
1243  std::cout << " NaN/Inf encountered";
1244  std::cout << " (nroi=" << nroi;
1245  std::cout << ", model=" << model;
1246  std::cout << ", cos_omega=" << cos_omega;
1247  std::cout << ", sin_omega=" << sin_omega;
1248  std::cout << ")" << std::endl;
1249  }
1250  #endif
1251 
1252  } // endif: sky intensity was positive
1253 
1254  // Return Nroi
1255  return nroi;
1256 }
1257 
1258 
1259 /*==========================================================================
1260  = =
1261  = Helper class methods for diffuse models =
1262  = =
1263  ==========================================================================*/
1264 
1265 /***********************************************************************//**
1266  * @brief Kernel for IRF offest angle integration of the diffuse source model
1267  *
1268  * @param[in] theta Offset angle with respect to observed photon direction
1269  * (radians).
1270  *
1271  * Computes
1272  *
1273  * \f[
1274  * K(\theta | E, t) = \sin \theta \times PSF(\theta)
1275  * \int_{0}^{2\pi}
1276  * S_{\rm p}(\theta, \phi | E, t) \,
1277  * Aeff(\theta, \phi) \,
1278  * Edisp(\theta, \phi) d\phi
1279  * \f]
1280  *
1281  * The PSF is assumed to be azimuthally symmetric, hence the PSF is computed
1282  * outside the azimuthal integration.
1283  *
1284  * Note that the integration is only performed for \f$\theta>0\f$. Otherwise
1285  * zero is returned.
1286  ***************************************************************************/
1287 double cta_irf_diffuse_kern_theta::eval(const double& theta)
1288 {
1289  // Initialise result
1290  double irf = 0.0;
1291 
1292  // Continue only if offset angle is positive
1293  if (theta > 0.0) {
1294 
1295  // Get PSF value. We can do this externally to the azimuthal
1296  // integration as the PSF is so far azimuthally symmetric. Once
1297  // we introduce asymmetries, we have to move this done into the
1298  // Phi kernel method/
1299  double psf = m_rsp->psf(theta, m_theta, m_phi, m_zenith, m_azimuth, m_srcLogEng);
1300 
1301  // Continue only if PSF is positive
1302  if (psf > 0.0) {
1303 
1304  // Precompute terms needed for the computation of the angular
1305  // distance between the true photon direction and the pointing
1306  // direction (i.e. the camera centre)
1307  double sin_theta = std::sin(theta);
1308  double cos_theta = std::cos(theta);
1309  double sin_ph = sin_theta * m_sin_eta;
1310  double cos_ph = cos_theta * m_cos_eta;
1311 
1312  // Setup kernel for azimuthal integration
1313  cta_irf_diffuse_kern_phi integrand(m_rsp,
1314  m_model,
1315  m_rot,
1316  m_zenith,
1317  m_azimuth,
1318  m_srcEng,
1319  m_srcTime,
1320  m_srcLogEng,
1321  m_obsEng,
1322  sin_theta,
1323  cos_theta,
1324  sin_ph,
1325  cos_ph);
1326 
1327  // Set number of azimuthal integration iterations
1328  int iter = gammalib::iter_phi(theta, m_resolution,
1330 
1331  // Integrate over phi
1332  GIntegral integral(&integrand);
1333  integral.fixed_iter(iter);
1334  irf = integral.romberg(0.0, gammalib::twopi) * psf * sin_theta;
1335 
1336  // Compile option: Debugging
1337  #if defined(G_DEBUG_INTEGRAL)
1338  if (!integral.isvalid()) {
1339  std::cout << "cta_irf_diffuse_kern_theta(theta=";
1340  std::cout << theta*gammalib::rad2deg << " deg) psf=";
1341  std::cout << psf << ":" << std::endl;
1342  std::cout << integral.print() << std::endl;
1343  }
1344  #endif
1345 
1346  // Compile option: Check for NaN/Inf
1347  #if defined(G_NAN_CHECK)
1349  std::cout << "*** ERROR: cta_irf_diffuse_kern_theta";
1350  std::cout << "(theta=" << theta << "):";
1351  std::cout << " NaN/Inf encountered";
1352  std::cout << " (irf=" << irf;
1353  std::cout << ", psf=" << psf;
1354  std::cout << ", sin_theta=" << sin_theta;
1355  std::cout << ", cos_theta=" << cos_theta;
1356  std::cout << ", sin_ph=" << sin_ph;
1357  std::cout << ", cos_ph=" << cos_ph;
1358  std::cout << ")";
1359  std::cout << std::endl;
1360  }
1361  #endif
1362 
1363  } // endif: PSF was positive
1364 
1365  } // endif: offset angle was positive
1366 
1367  // Return result
1368  return irf;
1369 }
1370 
1371 
1372 /***********************************************************************//**
1373  * @brief Kernel for IRF azimuth angle integration of the diffuse source model
1374  *
1375  * @param[in] phi Azimuth angle around observed photon direction (radians).
1376  *
1377  * Computes
1378  *
1379  * \f[
1380  * S_{\rm p}(\theta, \phi | E, t) \,
1381  * Aeff(\theta, \phi) \,
1382  * Edisp(\theta, \phi)
1383  * \f]
1384  *
1385  * As the coordinates \f$(\theta, \phi)\f$ are given in the reference frame
1386  * of the observed photon direction, some coordinate transformations have
1387  * to be performed.
1388  *
1389  * First, \f$(\theta, \phi)\f$ are transformed into the celestial reference
1390  * frame using the rotation matrix.
1391  *
1392  * Then, the offset angle of the true photon direction is computed in the
1393  * camera system (so far we do not compute the azimuth angle as we assume an
1394  * azimuthally symmetric response).
1395  *
1396  * @todo Optimize computation of sky direction in native coordinates
1397  * @todo Implement azimuth angle computation of true photon in camera
1398  * @todo Replace (theta,phi) by (delta,alpha)
1399  ***************************************************************************/
1400 double cta_irf_diffuse_kern_phi::eval(const double& phi)
1401 {
1402  // Initialise result
1403  double irf = 0.0;
1404 
1405  // Compute sine and cosine of azimuth angle
1406  double sin_phi = std::sin(phi);
1407  double cos_phi = std::cos(phi);
1408 
1409  // Compute sky direction vector in native coordinates
1410  m_native[0] = -cos_phi * m_sin_theta;
1411  m_native[1] = sin_phi * m_sin_theta;
1412  m_native[2] = m_cos_theta;
1413 
1414  // Rotate from native into celestial system
1415  GVector cel = *m_rot * m_native;
1416 
1417  // Set sky direction
1418  GSkyDir srcDir;
1419  srcDir.celvector(cel);
1420 
1421  // Set photon
1422  m_photon.dir(srcDir);
1425 
1426  // Get sky intensity for photon
1427  double intensity = m_model->eval(m_photon);
1428 
1429  // Continue only if sky intensity is positive
1430  if (intensity > 0.0) {
1431 
1432  // Compute true photon offset angle in camera system [radians]
1433  double offset = std::acos(m_cos_ph + m_sin_ph * cos_phi);
1434 
1435  //TODO: Compute true photon azimuth angle in camera system [radians]
1436  double azimuth = 0.0;
1437 
1438  // Evaluate model times the effective area
1439  irf = intensity *
1440  m_rsp->aeff(offset, azimuth, m_zenith, m_azimuth, m_srcLogEng);
1441 
1442  // Optionally take energy dispersion into account
1443  if (m_rsp->use_edisp() && irf > 0.0) {
1444  irf *= m_rsp->edisp(m_obsEng, m_srcEng, offset, azimuth,
1445  m_zenith, m_azimuth);
1446  }
1447 
1448  // Compile option: Check for NaN/Inf
1449  #if defined(G_NAN_CHECK)
1451  std::cout << "*** ERROR: cta_irf_diffuse_kern_phi::eval";
1452  std::cout << "(phi=" << phi << "):";
1453  std::cout << " NaN/Inf encountered";
1454  std::cout << " (irf=" << irf;
1455  std::cout << ", intensity=" << intensity;
1456  std::cout << ", offset=" << offset;
1457  std::cout << ", azimuth=" << azimuth;
1458  std::cout << ")";
1459  std::cout << std::endl;
1460  }
1461  #endif
1462 
1463  } // endif: sky intensity was positive
1464 
1465  // Return
1466  return irf;
1467 }
1468 
1469 
1470 /***********************************************************************//**
1471  * @brief Kernel for Nroi offest angle integration of diffuse model
1472  *
1473  * @param[in] theta Offset angle with respect to ROI centre (radians).
1474  *
1475  * Computes
1476  *
1477  * \f[
1478  * K(\theta | E, t) = \sin \theta \times
1479  * \int_{0}^{2\pi}
1480  * S_{\rm p}(\theta, \phi | E, t) \,
1481  * N_{\rm pred}(\theta, \phi) d\phi
1482  * \f]
1483  *
1484  * This method integrates a diffuse model for a given offset angle with
1485  * respect to the ROI centre over all azimuth angles (from 0 to 2pi). The
1486  * integration is only performed for positive offset angles, otherwise 0 is
1487  * returned.
1488  *
1489  * Integration is done using the Romberg algorithm. The integration kernel
1490  * is defined by the helper class cta_npred_diffuse_kern_phi.
1491  *
1492  * Note that the integration precision was adjusted trading-off between
1493  * computation time and computation precision. A value of 1e-4 was judged
1494  * appropriate.
1495  ***************************************************************************/
1496 double cta_nroi_diffuse_kern_theta::eval(const double& theta)
1497 {
1498  // Initialise Nroi value
1499  double nroi = 0.0;
1500 
1501  // Continue only if offset angle is positive
1502  if (theta > 0.0) {
1503 
1504  // Compute sine of offset angle
1505  double sin_theta = std::sin(theta);
1506 
1507  // Setup phi integration kernel
1508  cta_nroi_diffuse_kern_phi integrand(m_rsp,
1509  m_obs,
1510  m_model,
1511  m_rot,
1512  m_srcEng,
1513  m_srcTime,
1514  m_obsEng,
1515  m_obsTime,
1516  theta,
1517  sin_theta);
1518 
1519  // Integrate over phi
1520  GIntegral integral(&integrand);
1521  integral.fixed_iter(m_iter);
1522  nroi = integral.romberg(0.0, gammalib::twopi, m_iter) * sin_theta;
1523  #if defined(G_DEBUG_INTEGRAL)
1524  if (!integral.isvalid()) {
1525  std::cout << "cta_nroi_diffuse_kern_theta(theta=";
1526  std::cout << theta*gammalib::rad2deg << " deg):" << std::endl;
1527  std::cout << integral.print() << std::endl;
1528  }
1529  #endif
1530 
1531  // Debug: Check for NaN
1532  #if defined(G_NAN_CHECK)
1533  if (gammalib::is_notanumber(nroi) || gammalib::is_infinite(nroi)) {
1534  std::cout << "*** ERROR: cta_nroi_radial_kern_theta::eval";
1535  std::cout << "(theta=" << theta << "):";
1536  std::cout << " NaN/Inf encountered";
1537  std::cout << " (nroi=" << nroi;
1538  std::cout << ", sin_theta=" << sin_theta;
1539  std::cout << ")" << std::endl;
1540  }
1541  #endif
1542 
1543  } // endif: offset angle was positive
1544 
1545  // Return Nroi
1546  return nroi;
1547 }
1548 
1549 
1550 /***********************************************************************//**
1551  * @brief Kernel for Nroi azimuth angle integration of diffuse model
1552  *
1553  * @param[in] phi Azimuth angle with respect to ROI centre (radians).
1554  *
1555  * Computes
1556  *
1557  * \f[
1558  * S_{\rm p}(\theta, \phi | E, t) \, N_{\rm pred}(\theta, \phi)
1559  * \f]
1560  *
1561  * @todo Re-consider formula for possible simplification (dumb matrix
1562  * multiplication is definitely not the fastest way to do that
1563  * computation).
1564  ***************************************************************************/
1565 double cta_nroi_diffuse_kern_phi::eval(const double& phi)
1566 {
1567  // Initialise Nroi value
1568  double nroi = 0.0;
1569 
1570  // Compute sky direction vector in native coordinates
1571  double cos_phi = std::cos(phi);
1572  double sin_phi = std::sin(phi);
1573  GVector native(-cos_phi*m_sin_theta, sin_phi*m_sin_theta, m_cos_theta);
1574 
1575  // Rotate from native into celestial system
1576  GVector cel = *m_rot * native;
1577 
1578  // Set sky direction
1579  GSkyDir srcDir;
1580  srcDir.celvector(cel);
1581 
1582  // Set Photon
1583  GPhoton photon(srcDir, m_srcEng, m_srcTime);
1584 
1585  // Get sky intensity for this photon
1586  double intensity = m_model->eval(photon);
1587 
1588  // Continue only if sky intensity is positive
1589  if (intensity > 0.0) {
1590 
1591  // Compute Nroi for this sky direction
1592  nroi = m_rsp->nirf(photon, m_obsEng, m_obsTime, *m_obs) * intensity;
1593 
1594  // Debug: Check for NaN
1595  #if defined(G_NAN_CHECK)
1596  if (gammalib::is_notanumber(nroi) || gammalib::is_infinite(nroi)) {
1597  std::cout << "*** ERROR: cta_nroi_diffuse_kern_phi::eval";
1598  std::cout << "(phi=" << phi << "):";
1599  std::cout << " NaN/Inf encountered";
1600  std::cout << " (nroi=" << nroi;
1601  std::cout << ", intensity=" << intensity;
1602  std::cout << ", cos_phi=" << cos_phi;
1603  std::cout << ", sin_phi=" << sin_phi;
1604  std::cout << ")" << std::endl;
1605  }
1606  #endif
1607 
1608  } // endif: sky intensity was positive
1609 
1610  // Return Nroi
1611  return nroi;
1612 }
1613 
1614 
1615 /*==========================================================================
1616  = =
1617  = Helper class methods for stacked analysis =
1618  = =
1619  ==========================================================================*/
1620 
1621 /***********************************************************************//**
1622  * @brief Kernel for radial model integration over zenith angle
1623  *
1624  * @param[in] rho Radial distance from model centre (radians).
1625  * @return Integration kernel.
1626  *
1627  * Computes the integration kernel
1628  *
1629  * \f[
1630  * K(\rho | E, t) = \sin \rho \times
1631  * S_{\rm p}(\rho | E, t) \times
1632  * \int_{\omega} PSF(\rho, \omega) d\omega
1633  * \f]
1634  *
1635  * where
1636  * \f$\omega\f$ is the azimuth angle with respect to the model centre,
1637  * counted counterclockwise from the vector connecting the model centre
1638  * to the observed photon direction, and
1639  * \f$\rho\f$ is the radial distance from the model centre.
1640  * \f$S_{\rm p}(\rho | E, t)\f$ is the radial source model
1641  * for a given true photon energy and photon arrival time,
1642  * \f$PSF(\rho, \omega)\f$ is the point spread function.
1643  ***************************************************************************/
1644 double cta_psf_radial_kern_rho::eval(const double& rho)
1645 {
1646  // Initialise result
1647  double irf = 0.0;
1648 
1649  // Continue only if rho is positive
1650  if (rho > 0.0) {
1651 
1652  // Compute half length of the arc (in radians) from a circle with
1653  // radius rho that intersects with the point spread function, defined
1654  // as a circle with maximum radius m_delta_max
1655  double domega = 0.5 * gammalib::roi_arclength(rho,
1656  m_rho_obs,
1657  m_cos_rho_obs,
1658  m_sin_rho_obs,
1659  m_delta_max,
1660  m_cos_delta_max);
1661 
1662  // Continue only if arc length is positive
1663  if (domega > 0.0) {
1664 
1665  // Reduce rho by an infinite amount to avoid rounding errors
1666  // at the boundary of a sharp edged model
1667  double rho_kludge = rho - g_kludge_radius;
1668  if (rho_kludge < 0.0) {
1669  rho_kludge = 0.0;
1670  }
1671 
1672  // Evaluate sky model
1673  double model = m_model->eval(rho_kludge, m_srcEng, m_srcTime);
1674 
1675  // Debug: test if model is non positive
1676  #if defined(G_DEBUG_MODEL_ZERO)
1677  if (model <= 0.0) {
1678  std::cout << "*** WARNING: cta_psf_radial_kern_rho::eval";
1679  std::cout << " zero model for (rho)=(";
1680  std::cout << rho*gammalib::rad2deg << ")";
1681  std::cout << " rho-r_model=" << (rho-m_model->theta_max());
1682  std::cout << " radians" << std::endl;
1683  }
1684  #endif
1685 
1686  // Continue only if model is positive
1687  if (model > 0.0) {
1688 
1689  // Compute omega integration range
1690  double omega_min = -domega;
1691  double omega_max = +domega;
1692 
1693  // Precompute cosine and sine terms for azimuthal integration
1694  double cos_rho = std::cos(rho);
1695  double sin_rho = std::sin(rho);
1696  double cos_psf = cos_rho * m_cos_rho_obs;
1697  double sin_psf = sin_rho * m_sin_rho_obs;
1698 
1699  // Setup integration kernel
1700  cta_psf_radial_kern_omega integrand(m_rsp,
1701  m_model,
1702  m_srcDir,
1703  m_srcEng,
1704  m_srcTime,
1705  cos_psf,
1706  sin_psf);
1707 
1708  // Setup integrator
1709  GIntegral integral(&integrand);
1710  integral.fixed_iter(m_iter);
1711 
1712  // Integrate over omega
1713  irf = integral.romberg(omega_min, omega_max, m_iter) *
1714  model * sin_rho;
1715 
1716  // Compile option: Check for NaN/Inf
1717  #if defined(G_NAN_CHECK)
1719  std::cout << "*** ERROR: cta_psf_radial_kern_rho";
1720  std::cout << "(rho=" << rho << "):";
1721  std::cout << " NaN/Inf encountered";
1722  std::cout << " (irf=" << irf;
1723  std::cout << ", domega=" << domega;
1724  std::cout << ", model=" << model;
1725  std::cout << ", sin_rho=" << sin_rho << ")";
1726  std::cout << std::endl;
1727  }
1728  #endif
1729 
1730  } // endif: model was positive
1731 
1732  } // endif: arclength was positive
1733 
1734  } // endif: rho was positive
1735 
1736  // Return result
1737  return irf;
1738 }
1739 
1740 
1741 /***********************************************************************//**
1742  * @brief Kernel for radial model integration over azimuth angle
1743  *
1744  * @param[in] omega Azimuth angle (radians).
1745  *
1746  * Computes
1747  *
1748  * \f[
1749  * K(\omega | \rho, E, t) = PSF(\omega | \rho)
1750  * \f]
1751  *
1752  * where
1753  * \f$\omega\f$ is the azimuth angle with respect to the model centre,
1754  * counted counterclockwise from the vector connecting the model centre
1755  * to the observed photon direction, and
1756  * \f$\rho\f$ is the radial distance from the model centre.
1757  *
1758  * From the coordinates \f$(\rho,\omega)\f$ in the model system, the method
1759  * computes the angle between the true (\f$\vec{p}\f$) and observed
1760  * (\f$\vec{p'}\f$) photon arrival direction using
1761  *
1762  * \f[\delta = \arccos(\cos \rho \cos \rho_{\rm obs} +
1763  * \sin \rho \sin \rho_{\rm obs} \cos \omega)\f]
1764  *
1765  * where
1766  * \f$\rho_{\rm obs}\f$ is the angular distance between the observed photon
1767  * arrival direction \f$\vec{p'}\f$ and the model centre \f$\vec{m}\f$.
1768  * \f$\delta\f$ is used to compute the value of the point spread function.
1769  ***************************************************************************/
1770 double cta_psf_radial_kern_omega::eval(const double& omega)
1771 {
1772  // Initialise Irf value
1773  double irf = 0.0;
1774 
1775  // Compute Psf offset angle (radians)
1776  double delta = std::acos(m_cos_psf + m_sin_psf * std::cos(omega));
1777 
1778  // Evaluate Psf * model for this delta
1779  irf = m_rsp->psf()(m_srcDir, delta, m_srcEng);
1780 
1781  // Compile option: Check for NaN/Inf
1782  #if defined(G_NAN_CHECK)
1784  std::cout << "*** ERROR: cta_psf_radial_kern_omega::eval";
1785  std::cout << "(omega=" << omega << "):";
1786  std::cout << " NaN/Inf encountered";
1787  std::cout << " (irf=" << irf;
1788  std::cout << ", delta=" << delta;
1789  std::cout << ", omega=" << omega << ")";
1790  std::cout << std::endl;
1791  }
1792  #endif
1793 
1794  // Return
1795  return irf;
1796 }
1797 
1798 
1799 /***********************************************************************//**
1800  * @brief Kernel for PSF integration of radial model
1801  *
1802  * @param[in] delta PSF offset angle (radians).
1803  * @return Azimuthally integrated product between PSF and radial model.
1804  *
1805  * Computes the azimuthally integrated product of point spread function and
1806  * the radial model intensity. As the PSF is azimuthally symmetric, it is
1807  * not included in the azimuthally integration, but just multiplied on the
1808  * azimuthally integrated model. The method returns thus
1809  *
1810  * \f[
1811  * {\rm PSF}(\delta) \times
1812  * \int_0^{2\pi} {\rm M}(\delta, \phi) \sin \delta {\rm d}\phi
1813  * \f]
1814  *
1815  * where \f${\rm M}(\delta, \phi)\f$ is the radial model in the coordinate
1816  * system of the point spread function, defined by the angle \f$\delta\f$
1817  * between the true and the measured photon direction and the azimuth angle
1818  * \f$\phi\f$ around the measured photon direction.
1819  ***************************************************************************/
1820 double cta_psf_radial_kern_delta::eval(const double& delta)
1821 {
1822  // Initialise value
1823  double value = 0.0;
1824 
1825  // If we're at the Psf peak the model is zero (due to the sin(delta)
1826  // term. We thus only integrate for positive deltas.
1827  if (delta > 0.0) {
1828 
1829  // Get Psf for this delta
1830  double psf = m_rsp->psf()(m_srcDir, delta, m_srcEng);
1831 
1832  // Continue only if Psf is positive
1833  if (psf > 0.0) {
1834 
1835  // Compute half length of the arc (in radians) from a circle with
1836  // radius delta that intersects with the model, defined as a circle
1837  // with maximum radius m_theta_max
1838  double dphi = 0.5 * gammalib::roi_arclength(delta,
1839  m_delta_mod,
1842  m_theta_max,
1843  m_cos_theta_max);
1844 
1845  // Continue only if arc length is positive
1846  if (dphi > 0.0) {
1847 
1848  // Compute phi integration range
1849  double phi_min = -dphi;
1850  double phi_max = +dphi;
1851 
1852  // Precompute cosine and sine terms for azimuthal integration
1853  double sin_delta = std::sin(delta);
1854  double cos_delta = std::cos(delta);
1855  double sin_fact = sin_delta * m_sin_delta_mod;
1856  double cos_fact = cos_delta * m_cos_delta_mod;
1857 
1858  // Setup kernel for azimuthal integration of the spatial model
1859  cta_psf_radial_kern_phi integrand(m_model,
1860  m_srcEng,
1861  m_srcTime,
1862  sin_fact,
1863  cos_fact);
1864 
1865  // Setup integrator
1866  GIntegral integral(&integrand);
1867  integral.fixed_iter(m_iter);
1868 
1869  // Integrate over azimuth
1870  value = integral.romberg(phi_min, phi_max, m_iter) *
1871  psf * sin_delta;
1872 
1873  // Debug: Check for NaN
1874  #if defined(G_NAN_CHECK)
1875  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
1876  std::cout << "*** ERROR: cta_psf_radial_kern_delta::eval";
1877  std::cout << "(delta=" << delta << "):";
1878  std::cout << " NaN/Inf encountered";
1879  std::cout << " (value=" << value;
1880  std::cout << ")" << std::endl;
1881  }
1882  #endif
1883 
1884  } // endif: Psf value was positive
1885 
1886  } // endif: arc length was positive
1887 
1888  } // endif: delta was positive
1889 
1890  // Return kernel value
1891  return value;
1892 }
1893 
1894 
1895 /***********************************************************************//**
1896  * @brief Kernel for azimuthal radial model integration
1897  *
1898  * @param[in] phi Azimuth angle (radians).
1899  * @return Radial model value.
1900  *
1901  * Computes the value of the radial model at the position \f$(\delta,\phi)\f$
1902  * given in point spread function coordinates. The \f$\theta\f$ angle of the
1903  * radial model is computed using
1904  *
1905  * \f[
1906  * \theta = \arccos \left( \cos \delta \cos \zeta +
1907  * \sin \delta \sin \zeta \cos \phi \right)
1908  * \f]
1909  *
1910  * where \f$\delta\f$ is the angle between true and measured photon
1911  * direction, \f$\zeta\f$ is the angle between model centre and measured
1912  * photon direction, and \f$\phi\f$ is the azimuth angle with respect to the
1913  * measured photon direction, where \f$\phi=0\f$ corresponds to the
1914  * connecting line between model centre and measured photon direction.
1915  ***************************************************************************/
1916 double cta_psf_radial_kern_phi::eval(const double& phi)
1917 {
1918  // Compute radial model theta angle
1919  double theta = std::acos(m_cos_fact + m_sin_fact * std::cos(phi));
1920 
1921  // Reduce theta by an infinite amount to avoid rounding errors at the
1922  // boundary of a sharp edged model
1923  double theta_kluge = theta - 1.0e-12;
1924  if (theta_kluge < 0.0) {
1925  theta_kluge = 0.0;
1926  }
1927 
1928  // Get radial model value
1929  double value = m_model->eval(theta_kluge, m_srcEng, m_srcTime);
1930 
1931  // Debug: test if model is non positive
1932  #if defined(G_DEBUG_MODEL_ZERO)
1933  if (value <= 0.0) {
1934  std::cout << "*** WARNING: cta_psf_radial_kern_phi::eval";
1935  std::cout << " zero model for (phi)=(";
1936  std::cout << phi*gammalib::rad2deg << ")";
1937  std::cout << " theta-r_model=" << (theta-m_model->theta_max());
1938  std::cout << " radians" << std::endl;
1939  }
1940  #endif
1941 
1942  // Debug: Check for NaN
1943  #if defined(G_NAN_CHECK)
1944  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
1945  std::cout << "*** ERROR: cta_psf_radial_kern_phi::eval";
1946  std::cout << "(phi=" << phi << "):";
1947  std::cout << " NaN/Inf encountered";
1948  std::cout << " (value=" << value;
1949  std::cout << ")" << std::endl;
1950  }
1951  #endif
1952 
1953  // Return kernel value
1954  return value;
1955 }
1956 
1957 
1958 /***********************************************************************//**
1959  * @brief Kernel for elliptical model integration over zenith angle
1960  *
1961  * @param[in] rho Radial distance from model centre (radians).
1962  * @return Integration kernel.
1963  *
1964  * Computes the integration kernel
1965  *
1966  * \f[
1967  * K(\rho | E, t) = \sin \rho \times
1968  * \int_{\omega}
1969  * S_{\rm p}(\rho, \omega | E, t) \, PSF(\rho, \omega)
1970  * d\omega
1971  * \f]
1972  *
1973  * where
1974  * \f$\omega\f$ is the azimuth angle with respect to the model centre,
1975  * counted counterclockwise from the vector connecting the model centre
1976  * to the observed photon direction, and
1977  * \f$\rho\f$ is the radial distance from the model centre.
1978  * \f$S_{\rm p}(\rho, \omega | E, t)\f$ is the elliptical source model
1979  * for a given true photon energy and photon arrival time,
1980  * \f$PSF(\rho, \omega)\f$ is the point spread function.
1981  * The integration is over the \f$\omega\f$ values that are comprised within
1982  * the ellipse boundaries.
1983  ***************************************************************************/
1984 double cta_psf_elliptical_kern_rho::eval(const double& rho)
1985 {
1986  // Initialise result
1987  double irf = 0.0;
1988 
1989  // Continue only if rho is positive
1990  if (rho > 0.0) {
1991 
1992  // Compute half length of the arc (in radians) from a circle with
1993  // radius rho that intersects with the point spread function, defined
1994  // as a circle with maximum radius m_delta_max
1995  double domega = 0.5 * gammalib::roi_arclength(rho,
1996  m_rho_obs,
1997  m_cos_rho_obs,
1998  m_sin_rho_obs,
1999  m_delta_max,
2000  m_cos_delta_max);
2001 
2002  // Continue only if arc length is positive
2003  if (domega > 0.0) {
2004 
2005  // Precompute cosine and sine terms for azimuthal integration
2006  double cos_rho = std::cos(rho);
2007  double sin_rho = std::sin(rho);
2008  double cos_psf = cos_rho * m_cos_rho_obs;
2009  double sin_psf = sin_rho * m_sin_rho_obs;
2010 
2011  // Reduce rho by an infinite amount to avoid rounding errors
2012  // at the boundary of a sharp edged model (e.g. an elliptical
2013  // disk model)
2014  double rho_kludge = rho - g_ellipse_kludge_radius;
2015  if (rho_kludge < 0.0) {
2016  rho_kludge = 0.0;
2017  }
2018 
2019  // Setup integration kernel
2021  m_model,
2022  m_srcDir,
2023  m_srcEng,
2024  m_srcTime,
2026  rho_kludge,
2027  cos_psf,
2028  sin_psf);
2029 
2030  // Setup integrator
2031  GIntegral integral(&integrand);
2032  integral.fixed_iter(m_iter);
2033 
2034  // If the radius rho is not larger than the semiminor axis
2035  // boundary, the circle with that radius is fully contained in
2036  // the ellipse and we can just integrate over the relevant arc
2037  if (rho <= m_semiminor) {
2038 
2039  // Compute omega integration range
2040  double omega_min = -domega;
2041  double omega_max = +domega;
2042 
2043  // Integrate over omega
2044  irf = integral.romberg(omega_min, omega_max, m_iter) *
2045  sin_rho;
2046 
2047  } // endif: circle comprised in ellipse
2048 
2049  // ... otherwise there are arcs that intersect with the Psf circle
2050  else {
2051 
2052  // Compute half the arc length (in radians) of a circle of
2053  // radius rho, centred on the model, that intersects with
2054  // the ellipse boundary
2055  double arg1 = 1.0 - (m_semiminor*m_semiminor) / (rho*rho);
2056  double arg2 = 1.0 - (m_semiminor*m_semiminor) /
2058  double omega_width = std::acos(std::sqrt(arg1/arg2));
2059 
2060  // Continue only if the arclength is positive
2061  if (omega_width > 0.0) {
2062 
2063  // Compute azimuth angle difference between ellipse
2064  // position angle and position angle of observed
2065  // photon in the model system. This angle will define
2066  // the reference point around which the circle arc. Make
2067  // sure that omega_0 is within [-pi,pi] thus that the omega
2068  // intervals are between [-2pi,2pi]
2069  double omega_0 = m_posangle - m_posangle_obs;
2070  if (omega_0 > gammalib::pi) {
2071  omega_0 -= gammalib::pi;
2072  }
2073  else if (omega_0 < -gammalib::pi) {
2074  omega_0 += gammalib::pi;
2075  }
2076 
2077  // Compute azimuth angle intervals
2078  double omega1_min = omega_0 - omega_width;
2079  double omega1_max = omega_0 + omega_width;
2080  double omega2_min = omega1_min + gammalib::pi;
2081  double omega2_max = omega1_max + gammalib::pi;
2082 
2083  // Limit intervals to the intersection of the ellipse with
2084  // the Psf circle. This may lead to a split of intervals,
2085  // and we gather all these intervals in a special interval
2086  // pair containers
2087  cta_omega_intervals intervals1 =
2088  gammalib::limit_omega(omega1_min, omega1_max, domega);
2089  cta_omega_intervals intervals2 =
2090  gammalib::limit_omega(omega2_min, omega2_max, domega);
2091 
2092  // Integrate over all intervals for omega1
2093  for (int i = 0; i < intervals1.size(); ++i) {
2094  double min = intervals1[i].first;
2095  double max = intervals1[i].second;
2096  irf += integral.romberg(min, max, m_iter) * sin_rho;
2097  }
2098 
2099  // Integrate over all intervals for omega2
2100  for (int i = 0; i < intervals2.size(); ++i) {
2101  double min = intervals2[i].first;
2102  double max = intervals2[i].second;
2103  irf += integral.romberg(min, max, m_iter) * sin_rho;
2104  }
2105 
2106  } // endif: arc length was positive
2107 
2108  } // endelse: circle was not comprised in ellipse
2109 
2110  // Compile option: Check for NaN/Inf
2111  #if defined(G_NAN_CHECK)
2113  std::cout << "*** ERROR: cta_psf_elliptical_kern_rho";
2114  std::cout << "(rho=" << rho << "):";
2115  std::cout << " NaN/Inf encountered";
2116  std::cout << " (irf=" << irf;
2117  std::cout << ", domega=" << domega;
2118  std::cout << ", sin_rho=" << sin_rho << ")";
2119  std::cout << std::endl;
2120  }
2121  #endif
2122 
2123  } // endif: arc length was positive
2124 
2125  } // endif: rho was positive
2126 
2127  // Return result
2128  return irf;
2129 }
2130 
2131 
2132 /***********************************************************************//**
2133  * @brief Kernel for elliptical model integration over azimuth angle
2134  *
2135  * @param[in] omega Azimuth angle (radians).
2136  *
2137  * Computes
2138  *
2139  * \f[
2140  * K(\omega | \rho, E, t) = S_{\rm p}(\omega | \rho, E, t) \,
2141  * PSF(\omega | \rho)
2142  * \f]
2143  *
2144  * where
2145  * \f$\omega\f$ is the azimuth angle with respect to the model centre,
2146  * counted counterclockwise from the vector connecting the model centre
2147  * to the observed photon direction, and
2148  * \f$\rho\f$ is the radial distance from the model centre.
2149  *
2150  * From the coordinates \f$(\rho,\omega)\f$ in the model system, the method
2151  * computes the angle between the true (\f$\vec{p}\f$) and observed
2152  * (\f$\vec{p'}\f$) photon arrival direction using
2153  *
2154  * \f[\delta = \arccos(\cos \rho \cos \rho_{\rm obs} +
2155  * \sin \rho \sin \rho_{\rm obs} \cos \omega)\f]
2156  *
2157  * where
2158  * \f$\rho_{\rm obs}\f$ is the angular distance between the observed photon
2159  * arrival direction \f$\vec{p'}\f$ and the model centre \f$\vec{m}\f$.
2160  * \f$\delta\f$ is used to compute the value of the point spread function.
2161  ***************************************************************************/
2162 double cta_psf_elliptical_kern_omega::eval(const double& omega)
2163 {
2164  // Initialise Irf value
2165  double irf = 0.0;
2166 
2167  // Compute azimuth angle in model coordinate system (radians)
2168  double omega_model = omega + m_posangle_obs;
2169 
2170  // Evaluate sky model
2171  double model = m_model->eval(m_rho, omega_model, m_srcEng, m_srcTime);
2172 
2173  // Debug: test if model is non positive
2174  #if defined(G_DEBUG_MODEL_ZERO)
2175  if (model <= 0.0) {
2176  double semiminor_rad = m_model->semiminor() * gammalib::deg2rad;
2177  double semimajor_rad = m_model->semimajor() * gammalib::deg2rad;
2178  double diff_angle = omega_model - m_model->posangle() * gammalib::deg2rad;
2179  double cosinus = std::cos(diff_angle);
2180  double sinus = std::sin(diff_angle);
2181  double arg1 = semiminor_rad * cosinus;
2182  double arg2 = semimajor_rad * sinus;
2183  double r_ellipse = semiminor_rad * semimajor_rad /
2184  std::sqrt(arg1*arg1 + arg2*arg2);
2185  std::cout << "*** WARNING: cta_psf_elliptical_kern_omega::eval";
2186  std::cout << " zero model for (rho,omega)=";
2187  std::cout << m_rho*gammalib::rad2deg << ",";
2188  std::cout << omega*gammalib::rad2deg << ")";
2189  std::cout << " rho-r_ellipse=" << (m_rho-r_ellipse) << " radians";
2190  std::cout << std::endl;
2191  }
2192  #endif
2193 
2194  // Continue only if model is positive
2195  if (model > 0.0) {
2196 
2197  // Compute Psf offset angle (radians)
2198  double delta = std::acos(m_cos_psf + m_sin_psf * std::cos(omega));
2199 
2200  // Evaluate Psf * model for this delta
2201  irf = m_rsp->psf()(m_srcDir, delta, m_srcEng) * model;
2202 
2203  // Compile option: Check for NaN/Inf
2204  #if defined(G_NAN_CHECK)
2206  std::cout << "*** ERROR: cta_psf_elliptical_kern_omega::eval";
2207  std::cout << "(omega=" << omega << "):";
2208  std::cout << " NaN/Inf encountered";
2209  std::cout << " (irf=" << irf;
2210  std::cout << ", model=" << model;
2211  std::cout << ", delta=" << delta;
2212  std::cout << ", rho=" << m_rho;
2213  std::cout << ", omega=" << omega << ")";
2214  std::cout << std::endl;
2215  }
2216  #endif
2217 
2218  } // endif: model is positive
2219 
2220  // Return
2221  return irf;
2222 }
2223 
2224 
2225 /***********************************************************************//**
2226  * @brief Kernel for PSF integration of spatial model
2227  *
2228  * @param[in] delta PSF offset angle (radians).
2229  * @return Azimuthally integrated product between PSF and model.
2230  *
2231  * Computes the azimuthally integrated product of point spread function and
2232  * the spatial model intensity. As the PSF is azimuthally symmetric, it is
2233  * not included in the azimuthally integration, but just multiplied on the
2234  * azimuthally integrated model. The method returns thus
2235  *
2236  * \f[
2237  * {\rm PSF}(\delta) \times
2238  * \int_0^{2\pi} {\rm M}(\delta, \phi) \sin \delta {\rm d}\phi
2239  * \f]
2240  *
2241  * where \f${\rm M}(\delta, \phi)\f$ is the spatial model in the coordinate
2242  * system of the point spread function, defined by the angle \f$\delta\f$
2243  * between the true and the measured photon direction and the azimuth angle
2244  * \f$\phi\f$ around the measured photon direction.
2245  ***************************************************************************/
2246 double cta_psf_diffuse_kern_delta::eval(const double& delta)
2247 {
2248  // Initialise value
2249  double value = 0.0;
2250 
2251  // If we're at the PSF peak the model is zero (due to the sin(delta)
2252  // term. We thus only integrate for positive deltas.
2253  if (delta > 0.0) {
2254 
2255  // Get PSF for this delta
2256  value = m_rsp->psf()(m_srcDir, delta, m_srcEng);
2257 
2258  // Continue only if PSF is positive
2259  if (value > 0.0) {
2260 
2261  // Compute sine and cosine of delta
2262  double sin_delta = std::sin(delta);
2263  double cos_delta = std::cos(delta);
2264 
2265  // Setup kernel for azimuthal integration of the diffuse model
2268  sin_delta, cos_delta);
2269 
2270  // Set number of azimuthal integration iterations
2271  int iter = gammalib::iter_phi(delta, m_resolution,
2273 
2274  // Setup untegration
2275  GIntegral integral(&integrand);
2276 
2277  // Set fixed number of iterations
2278  integral.fixed_iter(iter);
2279 
2280  // Azimuthally integrate model
2281  value *= integral.romberg(0.0, gammalib::twopi) * sin_delta;
2282 
2283  } // endif: PSF value was positive
2284 
2285  } // endif: delta was positive
2286 
2287  // Debug: Check for NaN
2288  #if defined(G_NAN_CHECK)
2289  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
2290  std::cout << "*** ERROR: cta_psf_diffuse_kern_delta::eval";
2291  std::cout << "(delta=" << delta << "):";
2292  std::cout << " NaN/Inf encountered";
2293  std::cout << " (value=" << value;
2294  std::cout << ")" << std::endl;
2295  }
2296  #endif
2297 
2298  // Return kernel value
2299  return value;
2300 }
2301 
2302 
2303 /***********************************************************************//**
2304  * @brief Kernel for map integration of spatial model
2305  *
2306  * @param[in] phi Azimuth angle (radians).
2307  * @return Spatial model value.
2308  *
2309  * Computes the value of the spatial model at the position (delta,phi) given
2310  * in point spread function coordinates. The transformation from point
2311  * spread function coordinates into sky coordinates is done using a rotation
2312  * matrix that is pre-computed on entry.
2313  ***************************************************************************/
2314 double cta_psf_diffuse_kern_phi::eval(const double& phi)
2315 {
2316  // Compute sky direction vector in native coordinates
2317  double cos_phi = std::cos(phi);
2318  double sin_phi = std::sin(phi);
2319  GVector native(-cos_phi*m_sin_delta, sin_phi*m_sin_delta, m_cos_delta);
2320 
2321  // Rotate from native into celestial system
2322  GVector cel = *m_rot * native;
2323 
2324  // Set sky direction
2325  GSkyDir srcDir;
2326  srcDir.celvector(cel);
2327 
2328  // Set photon
2329  GPhoton photon(srcDir, m_srcEng, m_srcTime);
2330 
2331  // Compute map value for this sky direction
2332  double value = m_model->eval(photon);
2333 
2334  // Debug: Check for NaN
2335  #if defined(G_NAN_CHECK)
2336  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
2337  std::cout << "*** ERROR: cta_psf_diffuse_kern_phi::eval";
2338  std::cout << "(phi=" << phi << "):";
2339  std::cout << " NaN/Inf encountered";
2340  std::cout << " (value=" << value;
2341  std::cout << ")" << std::endl;
2342  }
2343  #endif
2344 
2345  // Return kernel value
2346  return value;
2347 }
double m_zeta
Distance model centre - measured photon.
const GCTAObservation * m_obs
CTA observation.
double m_cos_psf
Cosine term for PSF offset angle computation.
double eval(const double &rho)
Kernel for zenith angle Nroi integration or radial model.
CTA response helper classes definition.
Sky map class.
Definition: GSkyMap.hpp:89
std::vector< std::pair< double, double > > cta_omega_intervals
virtual double nroi(const GModelSky &model, const GEnergy &obsEng, const GTime &obsTime, const GObservation &obs) const
Return integral of event probability for a given sky model over ROI.
const GCTAResponseCube * m_rsp
CTA response.
double m_omega0
Position angle of ROI.
GEnergy m_srcEng
True photon energy.
GEnergy m_srcEng
True photon energy.
GEnergy m_srcEng
True photon energy.
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:1860
GTime m_srcTime
True photon time.
GTime m_obsTime
Observed photon arrival time.
GEnergy m_srcEng
True photon energy.
double m_theta_max
Maximum model radius.
double m_cos_ph
Cosine term for photon offset angle computation.
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:380
const GSkyProjection * projection(void) const
Returns pointer to sky projection.
Definition: GSkyMap.hpp:414
double eval(const double &rho)
Kernel for zenith angle Nroi integration of elliptical model.
virtual double eval(const GPhoton &photon, const bool &gradients=false) const =0
Kernel for radial model azimuth angle IRF integration.
double m_sin_zeta
Sine of zeta.
double m_sin_rho
Sine of offset angle.
double m_cos_rho_obs
Cosine of m_rho_obs.
double m_posangle_roi
Position angle of ROI.
double m_rho_obs
Distance of model centre from measured photon.
double m_sin_rho
Sine of offset angle.
GModelSpectral * spectral(void) const
Return spectral model component.
Definition: GModelSky.hpp:262
double m_cos_radius
Cosine of ROI+PSF radius.
double m_zenith
Zenith angle of source in Earth system.
const GModelSpatialElliptical * m_model
Elliptical model.
double m_cosroi
Cosine of ROI radius.
double m_delta_max
Maximum PSF radius.
const double pi
Definition: GMath.hpp:35
const GModelSpatialElliptical * m_model
Elliptical model.
virtual bool use_edisp(void) const
Signal if response uses energy dispersion.
const GMatrix * m_rot
Rotation matrix.
double eval(const double &rho)
Kernel for elliptical model integration over model&#39;s zenith angle.
GVector abs(const GVector &vector)
Computes absolute of vector elements.
Definition: GVector.cpp:1163
const GCTAResponseCube * m_rsp
CTA response.
const GModelSpatialElliptical * m_model
Elliptical model.
GEnergy m_srcEng
True photon energy.
double m_theta
Photon offset angle.
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1100
double eval(const double &theta)
Kernel for Nroi offest angle integration of diffuse model.
double m_radius_roi
ROI+PSF radius.
const GSkyMap & map(void) const
Get map.
double m_sin_psf
Sine term for PSF offset angle computation.
GEnergy m_obsEng
Observed photon energy.
virtual double eval(const GEnergy &srcEng, const GTime &srcTime=GTime(), const bool &gradients=false) const =0
const GMatrix * m_rot
Rotation matrix.
GModelTemporal * temporal(void) const
Return temporal model component.
Definition: GModelSky.hpp:278
const GCTAResponseIrf * m_rsp
CTA response.
double m_semiminor
Ellipse boundary semiminor axis.
double m_dist
Distance model-ROI centre.
double m_omega0
Azimuth of pointing in model system.
const GModelSpatial * m_model
Spatial model.
GTime m_srcTime
True photon arrival time.
double m_azimuth
Azimuth angle.
double m_sin_theta
Sine of offset angle.
GVector m_native
Pre-allocate memory for native coordinates.
GTime m_obsTime
Measured arrival time.
double m_cos_psf
Cosine term for PSF offset angle computation.
const GMatrix * m_rot
Rotation matrix.
GEnergy m_srcEng
True photon energy.
const GModelSpatial * m_model
Spatial model.
double m_delta_max
Maximum PSF radius.
double m_sin_rho_obs
Sine of m_rho_obs.
GPhoton m_photon
Pre-allocate memory for a photon object.
const GModelSpatialRadial * m_model
Radial model.
GEnergy m_obsEng
Observed photon energy.
double m_logE
Log10 of true photon energy (E/TeV).
double m_cos_delta_max
Cosine of maximum PSF radius.
double m_cos_ph
Cosine term in angular distance equation.
GTime m_srcTime
True photon arrival time.
double eval(const double &delta)
Integration kernel for npsf() method.
const GCTAResponseIrf * m_rsp
CTA response.
GEnergy m_srcEng
True photon energy.
int m_iter
Integration iterations.
double m_cos_eta
Cosine of angular distance between.
const GModelSpatialElliptical * m_model
Spatial model.
GEnergy m_obsEng
Observed photon energy.
double acos(const double &arg)
Computes acos by avoiding NaN due to rounding errors.
Definition: GMath.cpp:69
double log10TeV(void) const
Return log10 of energy in TeV.
Definition: GEnergy.cpp:457
GEnergy m_obsEng
Measured energy.
const GCTAResponseIrf * m_rsp
CTA response.
double m_zenith
Zenith angle.
int m_iter
Integration iterations.
GTime m_srcTime
True photon arrival time.
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
double m_resolution
Spatial map resolution.
int m_max_iter
Maximum number of Romberg iterations.
GEnergy m_obsEng
Measured event energy.
Gammalib tools definition.
GTime m_srcTime
True photon arrival time.
GIntegral class interface definition.
Definition: GIntegral.hpp:45
double m_sin_lambda
Sine of lambda.
double eval(const double &omega)
Kernel for radial model integration over azimuth angle.
const GModelSpatialElliptical * m_model
Elliptical model.
double m_lambda
Distance model centre - pointing.
int m_iter
Integration iterations.
double eval(const double &rho)
Kernel for radial model integration over zenith angle.
GTime m_srcTime
True photon arrival time.
double m_phi
Azimuth angle of source in camera system.
const GModelSpatialRadial * m_model
Radial spatial model.
double m_zenith
Pointing zenith angle.
const GMatrix * m_rot
Rotation matrix.
double min(const GVector &vector)
Computes minimum vector element.
Definition: GVector.cpp:843
Kernel for ellitpical model azimuth angle IRF integration.
double m_posangle
Ellipse boundary position angle.
double m_psf
PSF-ROI centre distance in radians.
const GCTAResponseIrf * m_rsp
CTA response.
const GCTAResponseIrf * m_rsp
CTA response.
const double g_ellipse_kludge_radius
About 0.2 arc seconds.
const GMatrix * m_rot
Rotation matrix.
GTime m_srcTime
True photon time.
GTime m_obsTime
Observed photon arrival time.
GSkyDir m_srcDir
True photon arrival direction.
const GCTAResponseIrf * m_rsp
CTA response.
GSkyDir m_srcDir
True photon sky direction.
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:184
GEnergy m_srcEng
True photon energy.
GSkyDir m_srcDir
True photon arrival direction.
double m_cos_zeta
Cosine of zeta.
double m_sin_delta_mod
Sine of m_delta_mod.
double m_cos_delta_mod
Cosine of m_delta_mod.
Class that handles photons.
Definition: GPhoton.hpp:47
const GCTAObservation * m_obs
CTA observation.
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:167
GSkyDir m_srcDir
True photon arrival direction.
GTime m_obsTime
Observed photon arrival time.
Definition of support function used by CTA classes.
double eval(const double &etrue)
Integration kernel for GCTAResponseIrf::nroi method.
double eval(const double &theta)
Kernel for IRF offest angle integration of the diffuse source model.
const GCTAPsf * psf(void) const
Return pointer to point spread function.
virtual double eval(const double &theta, const GEnergy &energy, const GTime &time, const bool &gradients=false) const =0
double m_radius
ROI+PSF radius.
double semiminor(void) const
Return semi-minor axis of ellipse.
double eval(const double &omega)
Kernel for radial model azimuth angle IRF integration.
double m_cos_lambda
Cosine of lambda.
const double deg2rad
Definition: GMath.hpp:43
GVector sqrt(const GVector &vector)
Computes square root of vector elements.
Definition: GVector.cpp:1268
const GCTAObservation * m_obs
CTA observation.
const GMatrix * m_rot
Rotation matrix.
double eval(const double &omega)
Kernel for azimuth angle Nroi integration of radial model.
double m_cos_dist
Cosine of distance model-ROI centre.
GEnergy m_srcEng
True photon energy.
GSkyDir m_srcDir
True photon sky direction.
double m_sin_ph
Sine term in angular distance equation.
double m_cos_rho_obs
Cosine of m_rho_obs.
const GCTAResponseIrf * m_rsp
CTA response function.
double m_posangle_obs
Photon position angle measured from model centre.
GEnergy m_srcEng
True photon energy.
Kernel for Psf phi angle integration used for stacked analysis.
double m_rho_roi
Distance between model and ROI centre.
GTime m_srcTime
True photon arrival time.
double m_delta_max
Maximum PSF radius.
GTime m_srcTime
True photon arrival time.
const GModelSpatial * m_model
Spatial model.
GTime m_srcTime
True photon time.
const GCTAAeff * aeff(void) const
Return pointer to effective area response.
int m_iter
Integration iterations.
double m_sin_eta
Sine of angular distance between.
const GCTACubePsf & psf(void) const
Return cube analysis point spread function.
const GCTAObservation * m_obs
CTA observation.
double m_cos_rho_roi
Cosine of m_rho_roi.
double m_sin_ph
Sine term for photon offset angle computation.
const GModelSpatial * m_model
Spatial model.
double eval(const double &omega)
Kernel for elliptical model integration over model&#39;s azimuth angle.
double m_azimuth
Azimuth angle.
double m_posangle_obs
Measured photon position angle from model centre.
double eval(const double &rho)
Kernel for radial model zenith angle integration of Irf.
void celvector(const GVector &vector)
Set sky direction from 3D vector in celestial coordinates.
Definition: GSkyDir.cpp:291
const GModelSpatialRadial * m_model
Radial model.
const GCTAResponseCube * m_rsp
Response cube.
double eval(const double &phi)
Kernel for map integration of spatial model.
const GCTAResponseIrf * m_rsp
CTA response.
double m_posangle_obs
Photon position angle measured from model centre.
double m_sin_psf
Sine term for PSF offset angle computation.
GEnergy m_srcEng
True photon energy.
Kernel for IRF azimuth angle integration of the diffuse source model.
double cdelt(const int &inx) const
Return pixel size.
Definition: GWcs.cpp:923
double nirf(const GPhoton &photon, const GEnergy &obsEng, const GTime &obsTime, const GObservation &obs) const
Return spatial integral of Instrument Response Function.
GEnergy m_obsEng
Observed photon energy.
double m_posangle
Ellipse boundary position angle.
double m_cos_theta
Cosine of offset angle.
double m_srcLogEng
True photon log energy.
double m_sin_ph
Sine term for photon offset angle computation.
double m_cos_rho
Cosine of offset angle.
GEnergy m_obsEng
Measured event energy.
double m_sin_rho_roi
Sine of m_rho_roi.
double m_semimajor
Ellipse boundary semimajor axis.
double m_cos_rho_pnt
Cosine of m_rho_pnt.
const GTime & time(void) const
Return photon time.
Definition: GPhoton.hpp:134
double m_omega_pnt
Azimuth of pointing in model system.
const GMatrix * m_rot
Rotation matrix.
double resolution(const GModelSpatial *model)
Determine resolution of spatial model.
double m_cos_fact
cos(delta)*cos(delta_mod)
Abstract CTA energy dispersion base class definition.
double m_roi
ROI radius in radians.
const GModelSpatialElliptical * m_model
Elliptical model.
virtual double eval(const GTime &srcTime, const bool &gradients=false) const =0
double m_cos_psf
Cosine term for PSF offset angle computation.
double m_sin_fact
sin(delta)*sin(delta_mod)
double m_srcLogEng
True photon log energy.
Isotropic spatial model class interface definition.
const GModelSpatialRadial * m_model
Radial spatial model.
Vector class interface definition.
GSkyDir m_srcDir
True photon arrival direction.
GEnergy m_srcEng
True photon energy.
double m_sin_rho_pnt
Sine of m_rho_pnt.
GTime m_srcTime
True arrival time.
GEnergy m_srcEng
True photon energy.
double m_cos_delta_max
Cosine of maximum PSF radius.
double m_cospsf
Cosine of PSF-ROI centre distance.
double max(const GVector &vector)
Computes maximum vector element.
Definition: GVector.cpp:872
GTime m_srcTime
True photon arrival time.
const double g_kludge_radius
Tiny angle (radians)
double m_posangle_roi
Position angle of ROI.
double m_cos_theta
Cosine of offset angle.
const GCTAResponseIrf * m_rsp
CTA response function.
double eval(const double &phi)
Kernel for azimuthal radial model integration.
GEnergy m_obsEng
Measured event energy.
double eval(const double &omega)
Kernel for elliptical model integration over azimuth angle.
GEnergy m_obsEng
Observed photon energy.
double eval(const double &rho)
Kernel for elliptical model integration over zenith angle.
const GCTAResponseIrf * m_rsp
CTA response.
GEnergy m_srcEng
True photon energy.
int m_iter
Integration iterations.
double m_phi
Photon azimuth angle.
const GCTAObservation * m_obs
CTA observation.
double m_rho_obs
Distance of model centre from measured photon.
const GCTAResponseCube * m_rsp
Response cube.
const GCTAResponseIrf * m_rsp
CTA response.
const GMatrix * m_rot
Rotation matrix.
double m_sin_dist
Sine of distance model-ROI centre.
cta_omega_intervals limit_omega(const double &min, const double &max, const double &domega)
Limit omega interval.
double m_sin_psf
Sine term for PSF offset angle computation.
Abstract world coordinate system base class.
Definition: GWcs.hpp:51
double eval(const double &delta)
Kernel for PSF integration of radial model.
int iter_phi(const double &rho, const double &resolution, const int &iter_min, const int &iter_max)
Determine number of azimuthal Romberg iterations.
GEnergy m_srcEng
True photon energy.
GEnergy m_srcEng
True photon energy.
double m_posangle_obs
Measured photon position angle from model centre.
std::string print(const GChatter &chatter=NORMAL) const
Print integral information.
Definition: GIntegral.cpp:918
double eval(const double &omega)
Kernel for azimuth angle Nroi integration of elliptical model.
GTime m_srcTime
True photon arrival time.
double m_rho_obs
Distance of model centre from measured photon.
double eval(const double &phi)
Kernel for Nroi azimuth angle integration of diffuse model.
double m_cos_theta_max
Cosine of m_theta_max.
double m_sin_rho_obs
Sine of m_rho_obs.
GEnergy m_srcEng
True photon energy.
GEnergy m_obsEng
Observed photon energy.
double eval(const double &delta)
Kernel for PSF integration of spatial model.
double m_cos_rho
Cosine of offset angle.
const GMatrix * m_rot
Rotation matrix.
double posangle(void) const
Return Position Angle of model.
const GModelSpatial * m_model
Spatial model.
double eval(const double &phi)
Kernel for IRF azimuth angle integration of the diffuse source model.
double m_cos_rho_obs
Cosine of m_rho_obs.
Spatial map cube model class interface definition.
const GObservation * m_obs
Observation.
GEnergy m_srcEng
True photon energy.
double m_semimajor
Ellipse boundary semimajor axis.
GEnergy m_obsEng
Measured event energy.
double m_cos_delta_max
Cosine of maximum PSF radius.
double m_semiminor
Ellipse boundary semiminor axis.
double m_azimuth
Azimuth angle of source in Earth system.
double m_azimuth
Azimuth angle.
double m_rho
Model zenith angle.
double m_posangle
Ellipse boundary position angle.
double m_cos_radius_roi
Cosine of m_radius_roi.
double m_theta
Offset angle of source in camera system.
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1226
double m_sin_psf
Sine term for PSF offset angle computation.
const GEnergy & energy(void) const
Return photon energy.
Definition: GPhoton.hpp:122
double m_sin_rho_obs
Sine of m_rho_obs.
virtual double theta_max(void) const =0
const GModelSpatialRadial * m_model
Radial model.
double m_cos_ph
Cosine term for photon offset angle computation.
int m_iter
Integration iterations.
void fixed_iter(const int &iter)
Set fixed number of iterations.
Definition: GIntegral.hpp:181
GTime m_obsTime
Observed photon arrival time.
Abstract spatial model base class.
double m_rho
Model zenith angle.
const GCTAResponseIrf * m_rsp
CTA response.
GTime m_obsTime
Observed photon arrival time.
const GModelSky * m_model
Sky model.
GTime m_srcTime
True photon time.
GTime m_srcTime
True photon arrival time.
const GMatrix * m_rot
Rotation matrix.
const double twopi
Definition: GMath.hpp:36
Vector class.
Definition: GVector.hpp:46
double m_sin_theta
Sine of offset angle.
virtual double eval(const double &theta, const double &posangle, const GEnergy &energy, const GTime &time, const bool &gradients=false) const =0
double m_sinpsf
Sine of PSF-ROI centre distance.
const double rad2deg
Definition: GMath.hpp:44
const GCTAResponseIrf * m_rsp
CTA response.
GEnergy m_srcEng
True photon energy.
const GCTAEdisp * edisp(void) const
Return pointer to energy dispersion.
double m_delta_mod
Distance of model from Psf.
double m_cos_delta_max
Cosine of maximum PSF radius.
const GCTAResponseIrf * m_rsp
CTA response.
Kernel for azimuth angle Nroi integration of elliptical model.
Integration class interface definition.
const GCTAResponseCube * m_rsp
CTA response.
double m_delta_max
Maximum PSF radius.
double m_omega0
Azimuth of pointing in model system.
Sky direction class.
Definition: GSkyDir.hpp:62
const GSkyMap & cube(void) const
Get map cube.
GTime m_obsTime
Observed photon arrival time.
int m_iter
Integration iterations.
const GCTAObservation * m_obs
CTA observation.
Kernel for Nroi azimuth angle integration of diffuse model.
GEnergy m_srcEng
True photon energy.
double m_azimuth
Pointing azimuth angle.
double semimajor(void) const
Return semi-major axis of ellipse.
int m_iter
Integration iterations.
const GSkyDir & dir(void) const
Return photon sky direction.
Definition: GPhoton.hpp:110
Spatial map model class interface definition.
GTime m_srcTime
True photon arrival time.
GEnergy m_obsEng
Measured event energy.
GTime m_srcTime
True photon arrival time.
Kernel for Psf phi angle integration used for stacked analysis.
const GModelSpatial * m_model
Spatial model.
const GCTAResponseCube * m_rsp
CTA response.
Abstract world coordinate system base class definition.
int m_min_iter
Minimum number of Romberg iterations.
double m_semimajor
Ellipse boundary semimajor axis.
double m_omega_pnt
Azimuth of pointing in model system.
GEnergy m_obsEng
Measured event energy.
Mathematical function definitions.
double m_cos_psf
Cosine term for PSF offset angle computation.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
Kernel for radial model azimuth angle integration.
double m_semiminor
Ellipse boundary semiminor axis.