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