GammaLib  2.0.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAResponseIrf.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAResponseIrf.cpp - CTA instrument response function class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2010-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 GCTAResponseIrf.cpp
23  * @brief CTA response class 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 <vector>
33 #include <string>
34 #include "GFits.hpp"
35 #include "GTools.hpp"
36 #include "GMath.hpp"
37 #include "GFilename.hpp"
38 #include "GIntegral.hpp"
39 #include "GCaldb.hpp"
40 #include "GSource.hpp"
41 #include "GRan.hpp"
42 #include "GModelSky.hpp"
44 #include "GModelSpatialRadial.hpp"
49 #include "GArf.hpp"
50 #include "GCTAObservation.hpp"
51 #include "GCTAResponseIrf.hpp"
52 #include "GCTAResponse_helpers.hpp"
53 #include "GCTAPointing.hpp"
54 #include "GCTAEventAtom.hpp"
55 #include "GCTAEventList.hpp"
56 #include "GCTARoi.hpp"
57 #include "GCTAException.hpp"
58 #include "GCTASupport.hpp"
59 #include "GCTAAeff.hpp"
60 #include "GCTAAeff2D.hpp"
61 #include "GCTAAeffArf.hpp"
62 #include "GCTAAeffPerfTable.hpp"
63 #include "GCTAPsf.hpp"
64 #include "GCTAPsf2D.hpp"
65 #include "GCTAPsfVector.hpp"
66 #include "GCTAPsfPerfTable.hpp"
67 #include "GCTAPsfKing.hpp"
68 #include "GCTAPsfTable.hpp"
69 #include "GCTAEdisp.hpp"
70 #include "GCTAEdisp2D.hpp"
71 #include "GCTAEdispRmf.hpp"
72 #include "GCTAEdispPerfTable.hpp"
73 #include "GCTABackground.hpp"
75 #include "GCTABackground3D.hpp"
76 
77 /* __ Method name definitions ____________________________________________ */
78 #define G_IRF "GCTAResponseIrf::irf(GInstDir&, GEnergy&, GTime&, GSkyDir&,"\
79  " GEnergy&, GTime&, GObservation&)"
80 #define G_MC "GCTAResponseIrf::mc(double&, GPhoton&, GObservation&, GRan&)"
81 #define G_READ "GCTAResponseIrf::read(GXmlElement&)"
82 #define G_WRITE "GCTAResponseIrf::write(GXmlElement&)"
83 #define G_LOAD_AEFF "GCTAResponseIrf::load_aeff(GFilename&)"
84 #define G_LOAD_PSF "GCTAResponseIrf::load_psf(GFilename&)"
85 #define G_LOAD_EDISP "GCTAResponseIrf::load_edisp(GFilename&)"
86 #define G_LOAD_BACKGROUND "GCTAResponseIrf::load_background(GFilename&)"
87 #define G_NIRF "GCTAResponseIrf::nirf(GPhoton&, GEnergy&, GTime&,"\
88  " GObservation&)"
89 #define G_IRF_RADIAL "GCTAResponseIrf::irf_radial(GEvent&, GSource&,"\
90  " GObservation&)"
91 #define G_IRF_ELLIPTICAL "GCTAResponseIrf::irf_elliptical(GEvent&, GSource&,"\
92  " GObservation&)"
93 #define G_IRF_DIFFUSE "GCTAResponseIrf::irf_diffuse(GEvent&, GSource&,"\
94  " GObservation&)"
95 #define G_NROI_RADIAL "GCTAResponseIrf::nroi_radial(GModelSky&, GEnergy&,"\
96  " GTime&, GEnergy&, GTime&, GObservation&)"
97 #define G_NROI_ELLIPTICAL "GCTAResponseIrf::nroi_elliptical(GModelSky&,"\
98  " GEnergy&, GTime&, GEnergy&, GTime&, GObservation&)"
99 #define G_NROI_DIFFUSE "GCTAResponseIrf::nroi_diffuse(GModelSky&, GEnergy&,"\
100  " GTime&, GEnergy&, GTime&, GObservation&)"
101 #define G_AEFF "GCTAResponseIrf::aeff(double&, double&, double&, double&,"\
102  " double&)"
103 #define G_PSF "GCTAResponseIrf::psf(double&, double&, double&, double&,"\
104  " double&)"
105 #define G_PSF_DELTA_MAX "GCTAResponseIrf::psf_delta_max(double&, double&,"\
106  " double&, double&, double&)"
107 
108 /* __ Macros _____________________________________________________________ */
109 
110 /* __ Coding definitions _________________________________________________ */
111 //#define G_USE_PSF_SYSTEM //!< Do radial Irf integrations in Psf system
112 
113 /* __ Debug definitions __________________________________________________ */
114 //#define G_DEBUG_IRF_RADIAL //!< Debug irf_radial method
115 //#define G_DEBUG_IRF_DIFFUSE //!< Debug irf_diffuse method
116 //#define G_DEBUG_IRF_ELLIPTICAL //!< Debug irf_elliptical method
117 //#define G_DEBUG_NROI_RADIAL //!< Debug npred_radial method
118 //#define G_DEBUG_NROI_ELLIPTICAL //!< Debug npred_elliptical method
119 //#define G_DEBUG_NROI_DIFFUSE //!< Debug npred_diffuse method
120 //#define G_DEBUG_PRINT_AEFF //!< Debug print() Aeff method
121 //#define G_DEBUG_PRINT_PSF //!< Debug print() Psf method
122 //#define G_DEBUG_PSF_DUMMY_SIGMA //!< Debug psf_dummy_sigma method
123 
124 /* __ Constants __________________________________________________________ */
125 
126 
127 /*==========================================================================
128  = =
129  = Constructors/destructors =
130  = =
131  ==========================================================================*/
132 
133 /***********************************************************************//**
134  * @brief Void constructor
135  *
136  * Constructs void CTA response.
137  ***************************************************************************/
139 {
140  // Initialise members
141  init_members();
142 
143  // Return
144  return;
145 }
146 
147 
148 /***********************************************************************//**
149  * @brief Copy constructor
150  *
151  * @param[in] rsp CTA response.
152  *
153  * Constructs CTA response by making a deep copy of an existing object.
154  **************************************************************************/
156 {
157  // Initialise members
158  init_members();
159 
160  // Copy members
161  copy_members(rsp);
162 
163  // Return
164  return;
165 }
166 
167 
168 /***********************************************************************//**
169  * @brief XML constructor
170  *
171  * @param[in] xml XML element.
172  *
173  * Construct CTA response from XML element.
174  ***************************************************************************/
176 {
177  // Initialise members
178  init_members();
179 
180  // Read information from XML element
181  read(xml);
182 
183  // Return
184  return;
185 }
186 
187 
188 /***********************************************************************//**
189  * @brief Response constructor
190  *
191  * @param[in] rspname Response file name.
192  * @param[in] caldb Calibration database.
193  *
194  * Create instance of CTA response by specifying the response name and the
195  * calibration database. The response name can be either a response identifier
196  * or a filename (see GCTAResponseIrf::load for more information).
197  ***************************************************************************/
198 GCTAResponseIrf::GCTAResponseIrf(const std::string& rspname,
199  const GCaldb& caldb) : GCTAResponse()
200 {
201  // Initialise members
202  init_members();
203 
204  // Set calibration database
205  m_caldb = caldb;
206 
207  // Load IRF
208  load(rspname);
209 
210  // Return
211  return;
212 }
213 
214 
215 /***********************************************************************//**
216  * @brief Destructor
217  *
218  * Destroys instance of CTA response object.
219  ***************************************************************************/
221 {
222  // Free members
223  free_members();
224 
225  // Return
226  return;
227 }
228 
229 
230 /*==========================================================================
231  = =
232  = Operators =
233  = =
234  ==========================================================================*/
235 
236 /***********************************************************************//**
237  * @brief Assignment operator
238  *
239  * @param[in] rsp CTA response.
240  * @return CTA response.
241  *
242  * Assigns CTA response object to another CTA response object. The assignment
243  * performs a deep copy of all information, hence the original object from
244  * which the assignment has been performed can be destroyed after this
245  * operation without any loss of information.
246  ***************************************************************************/
248 {
249  // Execute only if object is not identical
250  if (this != &rsp) {
251 
252  // Copy base class members
253  this->GCTAResponse::operator=(rsp);
254 
255  // Free members
256  free_members();
257 
258  // Initialise members
259  init_members();
260 
261  // Copy members
262  copy_members(rsp);
263 
264  } // endif: object was not identical
265 
266  // Return this object
267  return *this;
268 }
269 
270 
271 /*==========================================================================
272  = =
273  = Public methods =
274  = =
275  ==========================================================================*/
276 
277 /***********************************************************************//**
278  * @brief Clear instance
279  *
280  * Clears CTA response object by resetting all members to an initial state.
281  * Any information that was present in the object before will be lost.
282  ***************************************************************************/
284 {
285  // Free class members (base and derived classes, derived class first)
286  free_members();
288  this->GResponse::free_members();
289 
290  // Initialise members
291  this->GResponse::init_members();
293  init_members();
294 
295  // Return
296  return;
297 }
298 
299 
300 /***********************************************************************//**
301  * @brief Clone instance
302  *
303  * @return Pointer to deep copy of CTA response.
304  *
305  * Creates a clone (deep copy) of a CTA response object.
306  ***************************************************************************/
308 {
309  return new GCTAResponseIrf(*this);
310 }
311 
312 
313 /***********************************************************************//**
314  * @brief Return value of instrument response function
315  *
316  * @param[in] event Observed event.
317  * @param[in] photon Incident photon.
318  * @param[in] obs Observation.
319  *
320  * @todo Set polar angle phi of photon in camera system
321  ***************************************************************************/
322 double GCTAResponseIrf::irf(const GEvent& event,
323  const GPhoton& photon,
324  const GObservation& obs) const
325 {
326  // Retrieve CTA pointing and instrument direction
327  const GCTAPointing& pnt = gammalib::cta_pnt(G_IRF, obs);
328  const GCTAInstDir& dir = gammalib::cta_dir(G_IRF, event);
329 
330  // Get event attributes
331  const GSkyDir& obsDir = dir.dir();
332  const GEnergy& obsEng = event.energy();
333 
334  // Get photon attributes
335  const GSkyDir& srcDir = photon.dir();
336  const GEnergy& srcEng = photon.energy();
337  const GTime& srcTime = photon.time();
338 
339  // Get pointing direction zenith angle and azimuth [radians]
340  double zenith = pnt.zenith();
341  double azimuth = pnt.azimuth();
342 
343  // Get radial offset and polar angles of true photon in camera [radians]
344  double theta = pnt.dir().dist(srcDir);
345  double phi = 0.0; //TODO: Implement Phi dependence
346 
347  // Get log10(E/TeV) of true photon energy.
348  double srcLogEng = srcEng.log10TeV();
349 
350  // Determine angular separation between true and measured photon
351  // direction in radians
352  double delta = obsDir.dist(srcDir);
353 
354  // Get maximum angular separation for which PSF is significant
355  double delta_max = psf_delta_max(theta, phi, zenith, azimuth, srcLogEng);
356 
357  // Initialise IRF value
358  double irf = 0.0;
359 
360  // Compute only if we're sufficiently close to PSF
361  if (delta <= delta_max) {
362 
363  // Get effective area component
364  irf = aeff(theta, phi, zenith, azimuth, srcLogEng);
365 
366  // Multiply-in PSF
367  if (irf > 0.0) {
368 
369  // Get PSF component
370  irf *= psf(delta, theta, phi, zenith, azimuth, srcLogEng);
371 
372  // Multiply-in energy dispersion
373  if (use_edisp() && irf > 0.0) {
374 
375  // Multiply-in energy dispersion
376  irf *= edisp(obsEng, srcEng, theta, phi, zenith, azimuth);
377 
378  } // endif: energy dispersion was available and PSF was non-zero
379 
380  // Apply deadtime correction
381  irf *= obs.deadc(srcTime);
382 
383  } // endif: Aeff was non-zero
384 
385  } // endif: we were sufficiently close to PSF
386 
387  // Compile option: Check for NaN/Inf
388  #if defined(G_NAN_CHECK)
390  std::cout << "*** ERROR: GCTAResponseIrf::irf:";
391  std::cout << " NaN/Inf encountered";
392  std::cout << " (irf=" << irf;
393  std::cout << ", theta=" << theta;
394  std::cout << ", phi=" << phi << ")";
395  std::cout << std::endl;
396  }
397  #endif
398 
399  // Return IRF value
400  return irf;
401 }
402 
403 
404 /***********************************************************************//**
405  * @brief Return integral of event probability for a given sky model over ROI
406  *
407  * @param[in] model Sky model.
408  * @param[in] obsEng Observed photon energy.
409  * @param[in] obsTime Observed photon arrival time.
410  * @param[in] obs Observation.
411  * @return Event probability.
412  *
413  * @exception GException::feature_not_implemented
414  * Method is not implemented.
415  *
416  * Computes the integral
417  *
418  * \f[
419  * N_{\rm ROI}(E',t') = \int_{\rm ROI} P(p',E',t') dp'
420  * \f]
421  *
422  * of the event probability
423  *
424  * \f[
425  * P(p',E',t') = \int \int \int
426  * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp \, dE \, dt
427  * \f]
428  *
429  * for a given sky model \f$S(p,E,t)\f$ and response function
430  * \f$R(p',E',t'|p,E,t)\f$ over the Region of Interest (ROI).
431  ***************************************************************************/
432 double GCTAResponseIrf::nroi(const GModelSky& model,
433  const GEnergy& obsEng,
434  const GTime& obsTime,
435  const GObservation& obs) const
436 {
437  // Set number of iterations for Romberg integration.
438  static const int iter = 6;
439 
440  // Initialise Nroi value
441  double nroi = 0.0;
442 
443  // No time dispersion supported
444  const GTime& srcTime = obsTime;
445 
446  // If energy dispersion is requested then integrate over the relevant
447  // true photon energies ...
448  if (use_edisp()) {
449 
450  // Retrieve true energy boundaries
451  GEbounds ebounds = edisp()->etrue_bounds(obsEng);
452 
453  // Loop over all boundaries
454  for (int i = 0; i < ebounds.size(); ++i) {
455 
456  // Get true energy boundaries in MeV
457  double etrue_min = ebounds.emin(i).MeV();
458  double etrue_max = ebounds.emax(i).MeV();
459 
460  // Continue only if valid
461  if (etrue_max > etrue_min) {
462 
463  // Setup integration function
464  cta_nroi_kern integrand(this, &obs, &model, srcTime, obsEng, obsTime);
465  GIntegral integral(&integrand);
466 
467  // Set fixed number of iterations
468  integral.fixed_iter(iter);
469 
470  // Do Romberg integration
471  nroi += integral.romberg(etrue_min, etrue_max);
472 
473  } // endif: interval was valid
474 
475  } // endfor: looped over energy intervals
476 
477  } // endif: energy dispersion requested
478 
479  // ... otherwise evaluate
480  else {
481 
482  // No energy dispersion
483  const GEnergy& srcEng = obsEng;
484 
485  // Compute response components
486  double nroi_spatial = this->nroi(model, srcEng, srcTime, obsEng, obsTime, obs);
487  double nroi_spectral = model.spectral()->eval(srcEng, srcTime);
488  double nroi_temporal = model.temporal()->eval(srcTime);
489 
490  // Compute response
491  nroi = nroi_spatial * nroi_spectral * nroi_temporal;
492 
493  } // endelse: no energy dispersion requested
494 
495  // If required, apply instrument specific model scaling
496  if (model.has_scales()) {
497  nroi *= model.scale(obs.instrument()).value();
498  }
499 
500  // Compile option: Check for NaN/Inf
501  #if defined(G_NAN_CHECK)
502  if (gammalib::is_notanumber(nroi) || gammalib::is_infinite(nroi)) {
503  std::cout << "*** ERROR: GCTAResponseIrf::nroi:";
504  std::cout << " NaN/Inf encountered";
505  std::cout << " (nroi=" << nroi;
506  std::cout << ", obsEng=" << obsEng;
507  std::cout << ", obsTime=" << obsTime;
508  std::cout << ")" << std::endl;
509  }
510  #endif
511 
512  // Return response value
513  return nroi;
514 }
515 
516 
517 /***********************************************************************//**
518  * @brief Return true energy boundaries for a specific observed energy
519  *
520  * @param[in] obsEnergy Observed Energy.
521  * @return True energy boundaries for given observed energy.
522  ***************************************************************************/
524 {
525  // Initialise an empty boundary object
527 
528  // If energy dispersion is available then set the energy boundaries
529  if (edisp() != NULL) {
530  ebounds = edisp()->etrue_bounds(obsEnergy);
531  }
532 
533  // Return energy boundaries
534  return ebounds;
535 }
536 
537 
538 /***********************************************************************//**
539  * @brief Simulate event from photon
540  *
541  * @param[in] area Simulation surface area.
542  * @param[in] photon Photon.
543  * @param[in] obs Observation.
544  * @param[in] ran Random number generator.
545  * @return Simulated event.
546  *
547  * Simulates a CTA event using the response function from an incident photon.
548  * If the event is not detected a NULL pointer is returned.
549  *
550  * The method also applies a deadtime correction using a Monte Carlo process,
551  * taking into account temporal deadtime variations. For this purpose, the
552  * method makes use of the time dependent GObservation::deadc method.
553  *
554  * @todo Set polar angle phi of photon in camera system
555  * @todo Implement energy dispersion
556  ***************************************************************************/
557 GCTAEventAtom* GCTAResponseIrf::mc(const double& area,
558  const GPhoton& photon,
559  const GObservation& obs,
560  GRan& ran) const
561 {
562  // Initialise event
563  GCTAEventAtom* event = NULL;
564 
565  // Retrieve CTA pointing
566  const GCTAPointing& pnt = gammalib::cta_pnt(G_MC, obs);
567 
568  // Get pointing direction zenith angle and azimuth [radians]
569  double zenith = pnt.zenith();
570  double azimuth = pnt.azimuth();
571 
572  // Get radial offset and polar angles of true photon in camera [radians]
573  double theta = pnt.dir().dist(photon.dir());
574  double phi = 0.0; //TODO Implement Phi dependence
575 
576  // Compute effective area for photon
577  double srcLogEng = photon.energy().log10TeV();
578  double effective_area = aeff(theta, phi, zenith, azimuth, srcLogEng);
579 
580  // Compute limiting value
581  double ulimite = effective_area / area;
582 
583  // Warning if ulimite is larger than one
584  if (ulimite > 1.0) {
585  std::string msg = "Effective area "+
586  gammalib::str(effective_area)+
587  " cm2 is larger than simulation surface area "+
588  gammalib::str(area)+
589  " cm2 for photon energy "+
590  gammalib::str(photon.energy().TeV())+
591  " TeV. Simulations are inaccurate.";
592  gammalib::warning(G_MC, msg);
593  }
594 
595  // Continue only if event is detected
596  if ((ulimite > 0.0) && (ran.uniform() <= ulimite)) {
597 
598  // Apply deadtime correction
599  double deadc = obs.deadc(photon.time());
600  if (deadc >= 1.0 || ran.uniform() <= deadc) {
601 
602  // Simulate PSF and energy dispersion. Skip event if exception
603  // occurs
604  try {
605  // Simulate offset from photon arrival direction.
606  double delta = psf()->mc(ran, srcLogEng, theta, phi, zenith, azimuth) *
608  double alpha = 360.0 * ran.uniform();
609 
610  // Rotate sky direction by offset
611  GSkyDir sky_dir = photon.dir();
612  sky_dir.rotate_deg(alpha, delta);
613 
614  // Set measured photon arrival direction in instrument direction
615  GCTAInstDir inst_dir = pnt.instdir(sky_dir);
616 
617  // Set measured photon energy
618  GEnergy energy = photon.energy();
619  if (use_edisp()) {
620  energy = edisp()->mc(ran, photon.energy(), theta, phi,
621  zenith, azimuth);
622  }
623 
624  // Allocate event
625  event = new GCTAEventAtom;
626 
627  // Set event attributes
628  event->dir(inst_dir);
629  event->energy(energy);
630  event->time(photon.time());
631 
632  } // endtry: PSF and energy dispersion simulation successful
633 
634  // ... otherwise catch exception
636  ;
637  }
638 
639  } // endif: detector was alive
640 
641  } // endif: event was detected
642 
643  // Return event
644  return event;
645 }
646 
647 
648 /***********************************************************************//**
649  * @brief Read response from XML element
650  *
651  * @param[in] xml XML element.
652  *
653  * Reads information for a CTA observation from an XML element. The
654  * calibration database and response name can be specified using
655  *
656  * <observation name="..." id="..." instrument="...">
657  * ...
658  * <parameter name="Calibration" database="..." response="..."/>
659  * </observation>
660  *
661  * If even more control is required over the response, individual file names
662  * can be specified using
663  *
664  * <observation name="..." id="..." instrument="...">
665  * ...
666  * <parameter name="EffectiveArea" file="..."/>
667  * <parameter name="PointSpreadFunction" file="..."/>
668  * <parameter name="EnergyDispersion" file="..."/>
669  * <parameter name="Background" file="..."/>
670  * </observation>
671  *
672  ***************************************************************************/
674 {
675  // First check for "Calibration" parameter
676  if (gammalib::xml_has_par(xml, "Calibration")) {
677 
678  // Get parameter
679  const GXmlElement* par = gammalib::xml_get_par(G_READ, xml, "Calibration");
680 
681  // Read database and response
682  std::string xml_caldb = gammalib::strip_whitespace(par->attribute("database"));
683  std::string xml_rspname = gammalib::strip_whitespace(par->attribute("response"));
684 
685  // Set calibration database
686  GCaldb caldb("cta", xml_caldb);
687  this->caldb(caldb);
688 
689  // Load response
690  this->load(xml_rspname);
691 
692  // If there is a sigma attribute then set sigma values of Aeff and
693  // Background
694  if (par->has_attribute("sigma")) {
695 
696  // Get sigma value
697  double sigma = gammalib::todouble(par->attribute("sigma"));
698 
699  // If we have an effective area performance table then set sigma
700  // value
701  GCTAAeffPerfTable* perf = const_cast<GCTAAeffPerfTable*>(dynamic_cast<const GCTAAeffPerfTable*>(aeff()));
702  if (perf != NULL) {
703  perf->sigma(sigma);
704  }
705 
706  // If we have a background performance table then set sigma value
707  GCTABackgroundPerfTable* bgm = const_cast<GCTABackgroundPerfTable*>(dynamic_cast<const GCTABackgroundPerfTable*>(background()));
708  if (bgm != NULL) {
709  bgm->sigma(sigma);
710  }
711 
712  } // endif: sigma attribute specified
713 
714  // Store database and response names for later writing into the XML
715  // file (we do this now since the load() method clears all
716  // GCTAResponseIrf data members)
717  //m_xml_caldb = xml_caldb;
718  //m_xml_rspname = xml_rspname;
719 
720  } // endif: "Calibration" parameter found
721 
722  // ... otherwise check for components
723  else {
724 
725  // Handle effective area
726  if (gammalib::xml_has_par(xml, "EffectiveArea")) {
727 
728  // Get parameter
729  const GXmlElement* par = gammalib::xml_get_par(G_READ, xml, "EffectiveArea");
730 
731  // Get filename
732  std::string filename = gammalib::strip_whitespace(par->attribute("file"));
733 
734  // If filename is not empty then load effective area
735  if (!filename.empty()) {
736 
737  // Load effective area
738  load_aeff(filename);
739 
740  // Optional attributes
741  double thetacut = 0.0;
742  double scale = 1.0;
743  double sigma = 3.0;
744 
745  // Optionally extract thetacut (0.0 if no thetacut)
746  if (par->has_attribute("thetacut")) {
747  thetacut = gammalib::todouble(par->attribute("thetacut"));
748  }
749 
750  // Optionally extract scale factor (1.0 if no scale)
751  if (par->has_attribute("scale")) {
752  scale = gammalib::todouble(par->attribute("scale"));
753  }
754 
755  // Optionally extract sigma (3.0 if no sigma)
756  if (par->has_attribute("sigma")) {
757  sigma = gammalib::todouble(par->attribute("sigma"));
758  }
759 
760  // If we have an ARF then set attributes
761  GCTAAeffArf* arf = const_cast<GCTAAeffArf*>(dynamic_cast<const GCTAAeffArf*>(aeff()));
762  if (arf != NULL) {
763  arf->thetacut(thetacut);
764  arf->scale(scale);
765  arf->sigma(sigma);
766  }
767 
768  // If we have a performance table then set attributes
769  GCTAAeffPerfTable* perf = const_cast<GCTAAeffPerfTable*>(dynamic_cast<const GCTAAeffPerfTable*>(aeff()));
770  if (perf != NULL) {
771  perf->sigma(sigma);
772  }
773 
774  } // endif: effective area filename was valid
775 
776  } // endif: handled effective area
777 
778  // Handle PSF
779  if (gammalib::xml_has_par(xml, "PointSpreadFunction")) {
780 
781  // Get parameter
782  const GXmlElement* par = gammalib::xml_get_par(G_READ, xml, "PointSpreadFunction");
783 
784  // Get filename
785  std::string filename = gammalib::strip_whitespace(par->attribute("file"));
786 
787  // If filename is not empty then load point spread function
788  if (!filename.empty()) {
789  load_psf(filename);
790  }
791 
792  } // endif: handled PSF
793 
794  // Handle energy dispersion
795  if (gammalib::xml_has_par(xml, "EnergyDispersion")) {
796 
797  // Get parameter
798  const GXmlElement* par = gammalib::xml_get_par(G_READ, xml, "EnergyDispersion");
799 
800  // Get filename
801  std::string filename = gammalib::strip_whitespace(par->attribute("file"));
802 
803  // If filename is not empty then load energy dispersion
804  if (!filename.empty()) {
805  load_edisp(filename);
806  }
807 
808  } // endif: handled energy dispersion
809 
810  // Handle Background
811  if (gammalib::xml_has_par(xml, "Background")) {
812 
813  // Get parameter
814  const GXmlElement* par = gammalib::xml_get_par(G_READ, xml, "Background");
815 
816  // Get filename
817  std::string filename = gammalib::strip_whitespace(par->attribute("file"));
818 
819  // If filename is not empty then load background model
820  if (!filename.empty()) {
821  load_background(filename);
822  }
823 
824  // Optional attributes
825  double sigma = 3.0;
826 
827  // Optionally extract sigma (3.0 if no sigma)
828  if (par->has_attribute("sigma")) {
829  sigma = gammalib::todouble(par->attribute("sigma"));
830  }
831 
832  // If we have a performance table then set attributes
833  GCTABackgroundPerfTable* bgm = const_cast<GCTABackgroundPerfTable*>(dynamic_cast<const GCTABackgroundPerfTable*>(background()));
834  if (bgm != NULL) {
835  bgm->sigma(sigma);
836  }
837 
838  }
839 
840  } // endelse: handled components
841 
842  // If we have an ARF then remove thetacut if necessary
843  GCTAAeffArf* arf = const_cast<GCTAAeffArf*>(dynamic_cast<const GCTAAeffArf*>(aeff()));
844  if (arf != NULL) {
845  if (arf->thetacut() > 0.0) {
846  arf->remove_thetacut(*this);
847  }
848  }
849 
850  // Return
851  return;
852 }
853 
854 
855 /***********************************************************************//**
856  * @brief Write response information into XML element
857  *
858  * @param[in] xml XML element.
859  *
860  * Writes information for a CTA response into an XML element. If the
861  * calibration database and response name had been specified, the following
862  * output is written
863  *
864  * <observation name="..." id="..." instrument="...">
865  * ...
866  * <parameter name="Calibration" database="..." response="..."/>
867  * </observation>
868  *
869  * If even more control was required over the response and individual file
870  * names were specified, the following output is written
871  *
872  * <observation name="..." id="..." instrument="...">
873  * ...
874  * <parameter name="EffectiveArea" file="..."/>
875  * <parameter name="PointSpreadFunction" file="..."/>
876  * <parameter name="EnergyDispersion" file="..."/>
877  * <parameter name="Background" file="..."/>
878  * </observation>
879  ***************************************************************************/
881 {
882  // Determine number of existing parameter nodes in XML element
883  //int npars = xml.elements("parameter");
884 
885  // If we have a calibration database and response name, then set
886  // the information ...
887  if (!m_xml_caldb.empty() || !m_xml_rspname.empty()) {
888  GXmlElement* par = gammalib::xml_need_par(G_WRITE, xml, "Calibration");
889  par->attribute("database", m_xml_caldb);
890  par->attribute("response", m_xml_rspname);
891  }
892 
893  // ... otherwise add response components if they exist
894  else {
895 
896  // Add effective area if it exists
897  if (aeff() != NULL) {
898 
899  // Get effective area filename
900  GFilename filename = aeff()->filename();
901 
902  // Continue only if filename is not empty
903  if (!(filename.is_empty())) {
904 
905  // Get pointer to effective area
906  GXmlElement* par =
907  gammalib::xml_need_par(G_WRITE, xml, "EffectiveArea");
908 
909  // Initialise attributes
910  double thetacut = 0.0;
911  double scale = 1.0;
912  double sigma = 0.0;
913 
914  // Get optional ARF attributes
915  const GCTAAeffArf* arf =
916  dynamic_cast<const GCTAAeffArf*>(aeff());
917  if (arf != NULL) {
918  thetacut = arf->thetacut();
919  scale = arf->scale();
920  sigma = arf->sigma();
921  }
922 
923  // Get optional performance table attributes
924  const GCTAAeffPerfTable* perf =
925  dynamic_cast<const GCTAAeffPerfTable*>(aeff());
926  if (perf != NULL) {
927  sigma = perf->sigma();
928  }
929 
930  // Set attributes
931  par->attribute("file", filename);
932  if (thetacut > 0.0) {
933  par->attribute("thetacut", gammalib::str(thetacut));
934  }
935  if (scale != 1.0) {
936  par->attribute("scale", gammalib::str(scale));
937  }
938  if (sigma > 0.0) {
939  par->attribute("sigma", gammalib::str(sigma));
940  }
941 
942  }
943  }
944 
945  // Add PSF if it exists
946  if (psf() != NULL) {
947  if (!(psf()->filename().is_empty())) {
948  GXmlElement* par =
949  gammalib::xml_need_par(G_WRITE, xml, "PointSpreadFunction");
950  par->attribute("file", psf()->filename());
951  }
952  }
953 
954  // Add Edisp if it exists
955  if (edisp() != NULL) {
956  if (!(edisp()->filename().is_empty())) {
957  GXmlElement* par =
958  gammalib::xml_need_par(G_WRITE, xml, "EnergyDispersion");
959  par->attribute("file", edisp()->filename());
960  }
961  }
962 
963  // Add background if it exists
964  if (background() != NULL) {
965  if (!(background()->filename().is_empty())) {
966  GXmlElement* par =
967  gammalib::xml_need_par(G_WRITE, xml, "Background");
968  par->attribute("file", background()->filename());
969  }
970  }
971 
972  } // endelse: response components added
973 
974  // Return
975  return;
976 }
977 
978 
979 /***********************************************************************//**
980  * @brief Load CTA response
981  *
982  * @param[in] rspname CTA response name.
983  *
984  * Loads the CTA response with specified name @p rspname. The method first
985  * searchs for an appropriate response in the calibration database. If no
986  * appropriate response is found, the method takes the database root path
987  * and response name to build the full path to the response file, and tries
988  * to load the response from these paths.
989  *
990  * The method sets the calibration database and response names for writing
991  * of the response information into the XML file.
992  ***************************************************************************/
993 void GCTAResponseIrf::load(const std::string& rspname)
994 {
995  // Clear instance but conserve calibration database
996  GCaldb caldb = m_caldb;
997  clear();
998  m_caldb = caldb;
999 
1000  // First attempt reading the response using the GCaldb interface
1001  std::string expr = "NAME("+rspname+")";
1002  GFilename aeffname = m_caldb.filename("","","EFF_AREA","","",expr);
1003  GFilename psfname = m_caldb.filename("","","RPSF","","",expr);
1004  GFilename edispname = m_caldb.filename("","","EDISP","","",expr);
1005  GFilename bgdname = m_caldb.filename("","","BKG","","",expr);
1006 
1007  // Kluge: if the background filename is empty it may be because we have
1008  // and old response file that used "BGD" as the name of the background
1009  // extension. So we try to get here the old name
1010  if (bgdname.is_empty()) {
1011  bgdname = m_caldb.filename("","","BGD","","",expr);
1012  }
1013 
1014  // Signal usage of calibration database
1015  bool use_caldb = true;
1016 
1017  // If filenames are empty then build filenames from CALDB root path and
1018  // response name
1019  if (aeffname.is_empty()) {
1021  use_caldb = false;
1022  }
1023  if (psfname.is_empty()) {
1025  use_caldb = false;
1026  }
1027  if (edispname.is_empty()) {
1029  use_caldb = false;
1030  }
1031  if (bgdname.is_empty()) {
1033  use_caldb = false;
1034  }
1035 
1036  // Load effective area
1037  load_aeff(aeffname);
1038 
1039  // Load point spread function
1040  load_psf(psfname);
1041 
1042  // Load energy dispersion
1043  load_edisp(edispname);
1044 
1045  // Load background
1046  load_background(bgdname);
1047 
1048  // Remove theta cut
1049  GCTAAeffArf* arf = const_cast<GCTAAeffArf*>(dynamic_cast<const GCTAAeffArf*>(m_aeff));
1050  if (arf != NULL) {
1051  arf->remove_thetacut(*this);
1052  }
1053 
1054  // Store response name
1055  m_rspname = rspname;
1056 
1057  // If the calibration database has been used then store database and
1058  // response names for later writing into the XML file
1059  if (use_caldb) {
1060  m_xml_caldb = caldb.instrument();
1062  }
1063 
1064  // Return
1065  return;
1066 }
1067 
1068 
1069 /***********************************************************************//**
1070  * @brief Load effective area
1071  *
1072  * @param[in] filename Effective area filename.
1073  *
1074  * @exception GException::file_error
1075  * File or extension not found.
1076  *
1077  * Loads the effective area from a response file.
1078  *
1079  * If the file is a FITS file, the method will either use the extension name
1080  * specified with the filename, or if no extension name is given, search for
1081  * an `EFFECTIVE AREA` or `SPECRESP` extension in the file and open the
1082  * corresponding FITS table. If a column named `SPECRESP` is found in the
1083  * table, a GCTAAeffArf object will be allocated. Otherwise, a GCTAAeff2D
1084  * object will be allocated. In both cases, the method will extract the
1085  * optional `LO_THRES` and `HI_THRES` safe energy thresholds from the FITS
1086  * file header.
1087  *
1088  * If the file is not a FITS file, it will be interpreted as a
1089  * GCTAAeffPerfTable performance table.
1090  ***************************************************************************/
1092 {
1093  // Free any existing effective area instance
1094  if (m_aeff != NULL) delete m_aeff;
1095  m_aeff = NULL;
1096 
1097  // Check for existence of file
1098  if (!filename.exists()) {
1099  std::string msg = "File \""+filename+"\" not found. Please specify "
1100  "a valid effective area response file.";
1101  throw GException::file_error(G_LOAD_AEFF, msg);
1102  }
1103 
1104  // If file is a FITS file ...
1105  if (filename.is_fits()) {
1106 
1107  // Open FITS file
1108  GFits fits(filename);
1109 
1110  // Get the extension name. If an extension name has been specified
1111  // then use this name, otherwise preset the extension name according
1112  // to the GADF specifications. If no GADF compliant name was found,
1113  // search either for "EFFECTIVE AREA" or the "SPECRESP" extension.
1114  std::string extname = gammalib::gadf_hduclas4(fits, "AEFF_2D");
1115  if (filename.has_extname()) {
1116  extname = filename.extname();
1117  }
1118  if (extname.empty()) {
1120  extname = gammalib::extname_cta_aeff2d;
1121  }
1122  else if (fits.contains(gammalib::extname_arf)) {
1123  extname = gammalib::extname_arf;
1124  }
1125  }
1126 
1127  // Continue only if extension name is not empty
1128  if (!extname.empty()) {
1129 
1130  // Get FITS table
1131  const GFitsTable& table = *fits.table(extname);
1132 
1133  // Read safe energy thresholds if available
1134  if (table.has_card("LO_THRES")) {
1135  m_lo_safe_thres = table.real("LO_THRES");
1136  }
1137  if (table.has_card("HI_THRES")) {
1138  m_hi_safe_thres = table.real("HI_THRES");
1139  }
1140 
1141  // Check for specific table column
1142  if (table.contains(gammalib::extname_arf)) {
1143 
1144  // Close file
1145  fits.close();
1146 
1147  // Allocate Aeff from file
1148  m_aeff = new GCTAAeffArf(filename);
1149 
1150  } // endif: load as GCTAAeffArf
1151 
1152  else {
1153 
1154  // Close file
1155  fits.close();
1156 
1157  // Allocate Aeff from file
1158  m_aeff = new GCTAAeff2D(filename);
1159 
1160  } // endelse: load as GCTAAeff2D
1161 
1162  } // endif: extension name is not empty
1163 
1164  // Signal that no extension was found
1165  else {
1166  std::string msg = "FITS file \""+filename+"\" does not "
1167  "contain a valid effective area table. "
1168  "Please specify a valid effective area "
1169  "response file.";
1170  throw GException::file_error(G_LOAD_AEFF, msg);
1171  }
1172 
1173  } // endif: file was FITS file
1174 
1175  // ... else try handling an ASCII performance table
1176  else {
1177 
1178  // Allocate a performance table
1179  m_aeff = new GCTAAeffPerfTable(filename);
1180 
1181  } // endif: load as GCTAAeffPerfTable
1182 
1183  // Return
1184  return;
1185 }
1186 
1187 
1188 /***********************************************************************//**
1189  * @brief Load CTA PSF vector
1190  *
1191  * @param[in] filename FITS file name.
1192  *
1193  * @exception GException::file_error
1194  * File or extension not found.
1195  *
1196  * Loads the point spead function from a response file.
1197  *
1198  * If the file is a FITS file, the method will either use the extension name
1199  * specified with the filename, or if no extension name is given, search for
1200  * one of the following extension names
1201  *
1202  * POINT SPREAD FUNCTION
1203  * PSF
1204  * PSF_2D_TABLE
1205  *
1206  * in the file and open the corresponding FITS table. If columns named `GAMMA`
1207  * and `SIGMA` are found in the table, a GCTAPsfKing object will be allocated.
1208  *
1209  * If columns named `SCALE`, `SIGMA_1`, `AMPL_2`, `SIGMA_2`, `AMPL_3` and
1210  * `SIGMA_3` are found in the table, a GCTAPsf2D object will be allocated.
1211  *
1212  * If columns named `RAD_LO`, `RAD_HI`, and `RPSF` are found in the table, a
1213  * GCTAPsfTable object will be allocated.
1214  *
1215  * Otherwise, a CTAPsfVector object will be allocated.
1216  *
1217  * If the file is not a FITS file, it will be interpreted as a
1218  * GCTAPsfPerfTable performance table.
1219  ***************************************************************************/
1221 {
1222  // Free any existing point spread function instance
1223  if (m_psf != NULL) delete m_psf;
1224  m_psf = NULL;
1225 
1226  // Check for existence of file
1227  if (!filename.exists()) {
1228  std::string msg = "File \""+filename+"\" not found. Please specify "
1229  "a valid point spread function response file.";
1230  throw GException::file_error(G_LOAD_PSF, msg);
1231  }
1232 
1233  // If file is a FITS file ...
1234  if (filename.is_fits()) {
1235 
1236  // Open FITS file
1237  GFits fits(filename);
1238 
1239  // Get the extension name. If an extension name has been specified
1240  // then use this name, otherwise preset the extension name according
1241  // to the GADF specifications. If no GADF compliant name was found,
1242  // search either for "POINT SPREAD FUNCTION" or the "PSF" extension.
1243  std::string extname = gammalib::gadf_hduclas4(fits, "PSF_TABLE");
1244  if (filename.has_extname()) {
1245  extname = filename.extname();
1246  }
1247  if (extname.empty()) {
1248  if (fits.contains("POINT SPREAD FUNCTION")) {
1249  extname = "POINT SPREAD FUNCTION";
1250  }
1251  else if (fits.contains("PSF")) {
1252  extname = "PSF";
1253  }
1254  }
1255 
1256  // Continue only if extension name is not empty
1257  if (!extname.empty()) {
1258 
1259  // Get FITS table
1260  const GFitsTable& table = *fits.table(extname);
1261 
1262  // Check for King profile specific table columns
1263  if (table.contains("GAMMA") && table.contains("SIGMA")) {
1264 
1265  // Close FITS file
1266  fits.close();
1267 
1268  // Allocate King profile PSF
1269  m_psf = new GCTAPsfKing(filename);
1270 
1271  }
1272 
1273  // ... otherwise check for Gaussian profile specific table
1274  // columns
1275  else if (table.contains("SCALE") && table.contains("SIGMA_1") &&
1276  table.contains("AMPL_2") && table.contains("SIGMA_2") &&
1277  table.contains("AMPL_3") && table.contains("SIGMA_3")) {
1278 
1279  // Close FITS file
1280  fits.close();
1281 
1282  // Allocate Gaussian profile PSF
1283  m_psf = new GCTAPsf2D(filename);
1284 
1285  }
1286 
1287  // ... otherwise check for PSF table specific table columns
1288  else if (table.contains("RAD_LO") && table.contains("RAD_HI") &&
1289  table.contains("RPSF")) {
1290 
1291  // Close FITS file
1292  fits.close();
1293 
1294  // Allocate PSF table
1295  m_psf = new GCTAPsfTable(filename);
1296 
1297  }
1298 
1299  // ... otherwise try opening as vector PSF
1300  else {
1301 
1302  // Close FITS file
1303  fits.close();
1304 
1305  // Allocate vector PSF
1306  m_psf = new GCTAPsfVector(filename);
1307 
1308  }
1309 
1310  } // endif: extension name is not empty
1311 
1312  // Signal that no extension was found
1313  else {
1314  std::string msg = "FITS file \""+filename+"\" does not "
1315  "contain a valid point spread function table. "
1316  "Please specify a valid point spread function "
1317  "response file.";
1318  throw GException::file_error(G_LOAD_PSF, msg);
1319  }
1320 
1321  } // endif: file was FITS file
1322 
1323  // ... otherwise load file as a performance table
1324  else {
1325 
1326  // Allocate a performance table
1327  m_psf = new GCTAPsfPerfTable(filename);
1328 
1329  }
1330 
1331  // Return
1332  return;
1333 }
1334 
1335 
1336 /***********************************************************************//**
1337  * @brief Load energy dispersion information
1338  *
1339  * @param[in] filename Energy dispersion file name.
1340  *
1341  * @exception GException::file_error
1342  * File or extension not found.
1343  *
1344  * Loads the energy dispersion from a response file.
1345  *
1346  * If the file is a FITS file, the method will either use the extension name
1347  * specified with the filename, or if no extension name is given, search for
1348  * an `ENERGY DISPERSION` or `MATRIX` extension in the file and open the
1349  * corresponding FITS table. If columns named "MIGRA_LO" and "MIGRA_HI" are
1350  * found in the table, a GCTAEdisp2D object will be allocated. Otherwise, a
1351  * GCTAEdispRmf object will be allocated.
1352  *
1353  * If the file is not a FITS file, it will be interpreted as a
1354  * GCTAEdispPerfTable performance table.
1355  ***************************************************************************/
1357 {
1358  // Free any existing energy dispersion instance
1359  if (m_edisp != NULL) delete m_edisp;
1360  m_edisp = NULL;
1361 
1362  // Check for existence of file
1363  if (!filename.exists()) {
1364  std::string msg = "File \""+filename+"\" not found. Please specify "
1365  "a valid energy dispersion response file.";
1367  }
1368 
1369  // If file is a FITS file ...
1370  if (filename.is_fits()) {
1371 
1372  // Open FITS file
1373  GFits fits(filename);
1374 
1375  // Get the extension name. If an extension name has been specified
1376  // then use this name, otherwise preset the extension name according
1377  // to the GADF specifications. If no GADF compliant name was found,
1378  // search either for "ENERGY DISPERSION" or the "MATRIX" extension.
1379  std::string extname = gammalib::gadf_hduclas4(fits, "EDISP_2D");
1380  if (filename.has_extname()) {
1381  extname = filename.extname();
1382  }
1383  if (extname.empty()) {
1386  }
1387  else if (fits.contains(gammalib::extname_rmf)) {
1388  extname = gammalib::extname_rmf;
1389  }
1390  }
1391 
1392  // Continue only if extension name is not empty
1393  if (!extname.empty()) {
1394 
1395  // Get FITS table
1396  const GFitsTable& table = *fits.table(extname);
1397 
1398  // Check for 2D migration matrix
1399  if (table.contains("MIGRA_LO") && table.contains("MIGRA_HI")) {
1400 
1401  // Close FITS file
1402  fits.close();
1403 
1404  // Allocate 2D migration matrix
1405  m_edisp = new GCTAEdisp2D(filename);
1406 
1407  }
1408 
1409  // ... otherwise allocate RMF
1410  else {
1411 
1412  // Close FITS file
1413  fits.close();
1414 
1415  // Allocate Gaussian profile PSF
1416  m_edisp = new GCTAEdispRmf(filename);
1417 
1418  }
1419 
1420  } // endif: extension name is not empty
1421 
1422  // Signal that no extension was found
1423  else {
1424  std::string msg = "FITS file \""+filename+"\" does not "
1425  "contain a valid energy dispersion table. "
1426  "Please specify a valid energy dispersion "
1427  "response file.";
1429  }
1430 
1431  } // endif: file was FITS file
1432 
1433  // ... otherwise load file as a performance table
1434  else {
1435 
1436  // Allocate a performance table
1437  m_edisp = new GCTAEdispPerfTable(filename);
1438 
1439  }
1440 
1441  // Return
1442  return;
1443 }
1444 
1445 
1446 /***********************************************************************//**
1447  * @brief Load background model
1448  *
1449  * @param[in] filename Background model file name.
1450  *
1451  * @exception GException::file_error
1452  * File not found.
1453  *
1454  * Loads the background model from a response file.
1455  *
1456  * If the file is a FITS file the method allocates a GCTABackground3D
1457  * response and loads the information from the FITS file.
1458  *
1459  * If the file is not a FITS file it will be interpreted as a
1460  * GCTABackgroundPerfTable performance table.
1461  ***************************************************************************/
1463 {
1464  // Free any existing background model instance
1465  if (m_background != NULL) delete m_background;
1466  m_background = NULL;
1467 
1468  // Check for existence of file
1469  if (!filename.exists()) {
1470  std::string msg = "File \""+filename+"\" not found. Please specify "
1471  "a valid background response file.";
1473  }
1474 
1475  // If file is a FITS file than load background as 3D background
1476  if (filename.is_fits()) {
1477  m_background = new GCTABackground3D(filename);
1478  }
1479 
1480  // ... otherwise load background as performance table
1481  else {
1482  m_background = new GCTABackgroundPerfTable(filename);
1483  }
1484 
1485  // Return
1486  return;
1487 }
1488 
1489 
1490 /***********************************************************************//**
1491  * @brief Set offset angle dependence (degrees)
1492  *
1493  * @param[in] sigma Offset angle dependence value (degrees).
1494  *
1495  * Set the offset angle dependence for 1D effective area functions. The
1496  * method set the sigma value in case that the effective area function
1497  * is of type GCTAAeffArf or GCTAAeffPerfTable. Otherwise, nothing will
1498  * be done.
1499  ***************************************************************************/
1500 void GCTAResponseIrf::offset_sigma(const double& sigma)
1501 {
1502  // If effective area is an ARF then set offset angle
1503  GCTAAeffArf* arf = dynamic_cast<GCTAAeffArf*>(m_aeff);
1504  if (arf != NULL) {
1505  arf->sigma(sigma);
1506  }
1507 
1508  // If effective area is a performance table then set offset angle
1509  GCTAAeffPerfTable* prf = dynamic_cast<GCTAAeffPerfTable*>(m_aeff);
1510  if (prf != NULL) {
1511  prf->sigma(sigma);
1512  }
1513 
1514  // Return
1515  return;
1516 }
1517 
1518 
1519 /***********************************************************************//**
1520  * @brief Return offset angle dependence (degrees)
1521  *
1522  * @return Offset angle dependence value (degrees).
1523  *
1524  * Return the offset angle dependence for 1D effective area functions. The
1525  * method returns the sigma value in case that the effective area function
1526  * is of type GCTAAeffArf or GCTAAeffPerfTable. Otherwise, 0.0 will be
1527  * returned.
1528  ***************************************************************************/
1530 {
1531  // Initialise value
1532  double sigma = 0.0;
1533 
1534  // If effective area is an ARF then get offset angle
1535  GCTAAeffArf* arf = dynamic_cast<GCTAAeffArf*>(m_aeff);
1536  if (arf != NULL) {
1537  sigma = arf->sigma();
1538  }
1539 
1540  // If effective area is a performance table then get offset angle
1541  GCTAAeffPerfTable* prf = dynamic_cast<GCTAAeffPerfTable*>(m_aeff);
1542  if (prf != NULL) {
1543  sigma = prf->sigma();
1544  }
1545 
1546  // Return sigma
1547  return sigma;
1548 }
1549 
1550 
1551 /***********************************************************************//**
1552  * @brief Print CTA response information
1553  *
1554  * @param[in] chatter Chattiness.
1555  * @return String containing CTA response information.
1556  ***************************************************************************/
1557 std::string GCTAResponseIrf::print(const GChatter& chatter) const
1558 {
1559  // Initialise result string
1560  std::string result;
1561 
1562  // Continue only if chatter is not silent
1563  if (chatter != SILENT) {
1564 
1565  // Append header
1566  result.append("=== GCTAResponseIrf ===");
1567 
1568  // Append response information
1569  result.append("\n"+gammalib::parformat("Caldb mission")+m_caldb.mission());
1570  result.append("\n"+gammalib::parformat("Caldb instrument")+m_caldb.instrument());
1571  result.append("\n"+gammalib::parformat("Response name")+m_rspname);
1572  result.append("\n"+gammalib::parformat("Energy dispersion"));
1573  if (use_edisp()) {
1574  result.append("Used");
1575  }
1576  else {
1577  if (apply_edisp()) {
1578  result.append("Not available");
1579  }
1580  else {
1581  result.append("Not used");
1582  }
1583  }
1584 
1585  // Append safe energy threshold information
1586  result.append("\n"+gammalib::parformat("Safe energy range"));
1587  if (m_lo_safe_thres > 0.0 && m_hi_safe_thres) {
1588  result.append(gammalib::str(m_lo_safe_thres));
1589  result.append(" - ");
1590  result.append(gammalib::str(m_hi_safe_thres));
1591  result.append(" TeV");
1592  }
1593  else if (m_lo_safe_thres > 0.0) {
1594  result.append("> ");
1595  result.append(gammalib::str(m_lo_safe_thres));
1596  result.append(" TeV");
1597  }
1598  else if (m_hi_safe_thres > 0.0) {
1599  result.append("< ");
1600  result.append(gammalib::str(m_hi_safe_thres));
1601  result.append(" TeV");
1602  }
1603  else {
1604  result.append("undefined");
1605  }
1606 
1607  // Get reduced chatter level
1608  GChatter reduced_chatter = gammalib::reduce(chatter);
1609 
1610  // Append detailed information
1611  if (chatter >= NORMAL) {
1612 
1613  // Append calibration database
1614  result.append("\n"+m_caldb.print(reduced_chatter));
1615 
1616  // Append effective area information
1617  if (m_aeff != NULL) {
1618  result.append("\n"+m_aeff->print(reduced_chatter));
1619  }
1620 
1621  // Append point spread function information
1622  if (m_psf != NULL) {
1623  result.append("\n"+m_psf->print(reduced_chatter));
1624  }
1625 
1626  // Append energy dispersion information
1627  if (m_edisp != NULL) {
1628  result.append("\n"+m_edisp->print(reduced_chatter));
1629  }
1630 
1631  // Append background information
1632  if (m_background != NULL) {
1633  result.append("\n"+m_background->print(reduced_chatter));
1634  }
1635 
1636  // Append cache information
1637  result.append("\n"+m_irf_cache.print(reduced_chatter));
1638  result.append("\n"+m_nroi_cache.print(reduced_chatter));
1639 
1640  } // endif: appended detailed information
1641 
1642  } // endif: chatter was not silent
1643 
1644  // Return result
1645  return result;
1646 }
1647 
1648 
1649 /*==========================================================================
1650  = =
1651  = Low-level CTA response methods =
1652  = =
1653  ==========================================================================*/
1654 
1655 /***********************************************************************//**
1656  * @brief Return effective area (in units of cm2)
1657  *
1658  * @param[in] theta Radial offset angle of photon in camera (radians).
1659  * @param[in] phi Polar angle of photon in camera (radians).
1660  * @param[in] zenith Zenith angle of telescope pointing (radians).
1661  * @param[in] azimuth Azimuth angle of telescope pointing (radians).
1662  * @param[in] srcLogEng Log10 of true photon energy (E/TeV).
1663  * @return Effective area in units fo cm2.
1664  *
1665  * @exception GException::invalid_value
1666  * No effective area information found.
1667  *
1668  * Returns the effective area as function of the true photon position in the
1669  * camera system and the telescope pointing direction in the Earth system.
1670  ***************************************************************************/
1671 double GCTAResponseIrf::aeff(const double& theta,
1672  const double& phi,
1673  const double& zenith,
1674  const double& azimuth,
1675  const double& srcLogEng) const
1676 {
1677  // Throw an exception if instrument response is not defined
1678  if (m_aeff == NULL) {
1679  std::string msg = "No effective area information found in response.\n"
1680  "Please make sure that the instrument response is"
1681  " properly defined.";
1682  throw GException::invalid_value(G_AEFF, msg);
1683  }
1684 
1685  // Get effective area
1686  double aeff = (*m_aeff)(srcLogEng, theta, phi, zenith, azimuth);
1687 
1688  // Return effective area
1689  return aeff;
1690 }
1691 
1692 
1693 /***********************************************************************//**
1694  * @brief Return point spread function (in units of sr^-1)
1695  *
1696  * @param[in] delta Angular separation between true and measured photon
1697  * directions (radians).
1698  * @param[in] theta Radial offset angle of photon in camera (radians).
1699  * @param[in] phi Polar angle of photon in camera (radians).
1700  * @param[in] zenith Zenith angle of telescope pointing (radians).
1701  * @param[in] azimuth Azimuth angle of telescope pointing (radians).
1702  * @param[in] srcLogEng Log10 of true photon energy (E/TeV).
1703  *
1704  * @exception GException::invalid_value
1705  * No point spread function information found.
1706  *
1707  * Returns the point spread function for a given offset angle as function
1708  * of the true photon position in the camera system and the telescope
1709  * pointing direction in the Earth system.
1710  ***************************************************************************/
1711 double GCTAResponseIrf::psf(const double& delta,
1712  const double& theta,
1713  const double& phi,
1714  const double& zenith,
1715  const double& azimuth,
1716  const double& srcLogEng) const
1717 {
1718  // Throw an exception if instrument response is not defined
1719  if (m_psf == NULL) {
1720  std::string msg = "No point spread function information found in"
1721  " response.\n"
1722  "Please make sure that the instrument response is"
1723  " properly defined.";
1724  throw GException::invalid_value(G_PSF, msg);
1725  }
1726 
1727  // Compute PSF
1728  double psf = (*m_psf)(delta, srcLogEng, theta, phi, zenith, azimuth);
1729 
1730  // Return PSF
1731  return psf;
1732 }
1733 
1734 
1735 /***********************************************************************//**
1736  * @brief Return maximum angular separation (in radians)
1737  *
1738  * @param[in] theta Radial offset angle in camera (radians).
1739  * @param[in] phi Polar angle in camera (radians).
1740  * @param[in] zenith Zenith angle of telescope pointing (radians).
1741  * @param[in] azimuth Azimuth angle of telescope pointing (radians).
1742  * @param[in] srcLogEng Log10 of true photon energy (E/TeV).
1743  *
1744  * @exception GException::invalid_value
1745  * No point spread function information found.
1746  *
1747  * This method returns the maximum angular separation between true and
1748  * measured photon directions for which the PSF is non zero as function
1749  * of the true photon position in the camera system and the telescope
1750  * pointing direction in the Earth system.
1751  ***************************************************************************/
1752 double GCTAResponseIrf::psf_delta_max(const double& theta,
1753  const double& phi,
1754  const double& zenith,
1755  const double& azimuth,
1756  const double& srcLogEng) const
1757 {
1758  // Throw an exception if instrument response is not defined
1759  if (m_psf == NULL) {
1760  std::string msg = "No point spread function information found in"
1761  " response.\n"
1762  "Please make sure that the instrument response is"
1763  " properly defined.";
1765  }
1766 
1767  // Compute PSF
1768  double delta_max = m_psf->delta_max(srcLogEng, theta, phi, zenith, azimuth);
1769 
1770  // Return PSF
1771  return delta_max;
1772 }
1773 
1774 
1775 /***********************************************************************//**
1776  * @brief Return energy dispersion (in units of MeV\f$^{-1}\f$)
1777  *
1778  * @param[in] ereco Reconstructed event energy.
1779  * @param[in] etrue True photon energy.
1780  * @param[in] theta Radial offset angle in camera (radians).
1781  * @param[in] phi Polar angle in camera (radians).
1782  * @param[in] zenith Zenith angle of telescope pointing (radians).
1783  * @param[in] azimuth Azimuth angle of telescope pointing (radians).
1784  * @return Energy dispersion (MeV\f$^{-1}\f$).
1785  ***************************************************************************/
1786 double GCTAResponseIrf::edisp(const GEnergy& ereco,
1787  const GEnergy& etrue,
1788  const double& theta,
1789  const double& phi,
1790  const double& zenith,
1791  const double& azimuth) const
1792 {
1793  // Compute energy dispersion
1794  double edisp = (*m_edisp)(ereco, etrue, theta, phi, zenith, azimuth);
1795 
1796  // Return energy dispersion
1797  return edisp;
1798 }
1799 
1800 
1801 /***********************************************************************//**
1802  * @brief Return spatial integral of sky model
1803  *
1804  * @param[in] model Sky Model.
1805  * @param[in] srcEng True photon energy.
1806  * @param[in] srcTime True photon arrival time.
1807  * @param[in] obsEng Observed event energy.
1808  * @param[in] obsTime Observed event arrival time.
1809  * @param[in] obs Observation.
1810  *
1811  * Computes the integral
1812  *
1813  * \f[
1814  * N_{\rm ROI}(E',t'|E,t) = \int_{\rm ROI} P(p',E',t'|E,t) dp'
1815  * \f]
1816  *
1817  * of
1818  *
1819  * \f[
1820  * P(p',E',t'|E,t) = \int
1821  * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp
1822  * \f]
1823  *
1824  * over the Region of Interest (ROI) for a sky model \f$S(p,E,t)\f$ and the
1825  * response function \f$R(p',E',t'|p,E,t)\f$.
1826  ***************************************************************************/
1827 double GCTAResponseIrf::nroi(const GModelSky& model,
1828  const GEnergy& srcEng,
1829  const GTime& srcTime,
1830  const GEnergy& obsEng,
1831  const GTime& obsTime,
1832  const GObservation& obs) const
1833 {
1834  // Initialise response value
1835  double nroi = 0.0;
1836 
1837  // Set response value attributes
1838  std::string name = obs.id() + "::" + model.name();
1839 
1840  // Signal if spatial model has free parameters
1841  bool has_free_pars = model.spatial()->has_free_pars();
1842 
1843  // If the spatial model component has free parameters, or the response
1844  // cache should not be used, or the cache does not contain the requested
1845  // IRF value then compute the IRF value for the spatial model.
1846  if (has_free_pars ||
1847  !m_use_nroi_cache ||
1848  !m_nroi_cache.contains(name, obsEng, srcEng, &nroi)) {
1849 
1850  // Select method depending on the spatial model type
1851  switch (model.spatial()->code()) {
1853  nroi = nroi_ptsrc(model, srcEng, srcTime, obsEng, obsTime, obs);
1854  break;
1855  case GMODEL_SPATIAL_RADIAL:
1856  nroi = nroi_radial(model, srcEng, srcTime, obsEng, obsTime, obs);
1857  break;
1859  nroi = nroi_elliptical(model, srcEng, srcTime, obsEng, obsTime, obs);
1860  break;
1862  nroi = nroi_diffuse(model, srcEng, srcTime, obsEng, obsTime, obs);
1863  break;
1865  nroi = nroi_composite(model, srcEng, srcTime, obsEng, obsTime, obs);
1866  break;
1867  default:
1868  break;
1869  }
1870 
1871  } // endif: computed spatial model
1872 
1873  // If the spatial model has no free parameters and the response cache
1874  // should be used then put the IRF value in the response cache.
1875  if (!has_free_pars && m_use_nroi_cache) {
1876  m_nroi_cache.set(name, obsEng, srcEng, nroi);
1877  }
1878 
1879  // Return response value
1880  return nroi;
1881 }
1882 
1883 
1884 /***********************************************************************//**
1885  * @brief Return spatial integral of Instrument Response Function
1886  *
1887  * @param[in] photon Photon.
1888  * @param[in] obsEng Observed event energy.
1889  * @param[in] obsTime Observed event time.
1890  * @param[in] obs Observation.
1891  *
1892  * Computes the integral of the instrument response function over the Region
1893  * of Interest
1894  *
1895  * \f[
1896  * R(E',t'|p,E,t) = \int_{\rm ROI} R(p',E',t'|p,E,t) dp'
1897  * \f]
1898  ***************************************************************************/
1899 double GCTAResponseIrf::nirf(const GPhoton& photon,
1900  const GEnergy& obsEng,
1901  const GTime& obsTime,
1902  const GObservation& obs) const
1903 {
1904  // Retrieve CTA observation, ROI and pointing
1905  const GCTAObservation& cta = gammalib::cta_obs(G_NIRF, obs);
1906  const GCTARoi& roi = gammalib::cta_event_list(G_NIRF, obs).roi();
1907  const GCTAPointing& pnt = cta.pointing();
1908 
1909  // Get photon attributes
1910  const GSkyDir& srcDir = photon.dir();
1911  const GEnergy& srcEng = photon.energy();
1912  const GTime& srcTime = photon.time();
1913 
1914  // Get pointing direction zenith angle and azimuth [radians]
1915  double zenith = pnt.zenith();
1916  double azimuth = pnt.azimuth();
1917 
1918  // Get radial offset and polar angles of true photon in camera [radians]
1919  double theta = pnt.dir().dist(srcDir);
1920  double phi = 0.0; //TODO: Implement Phi dependence
1921 
1922  // Get log10(E/TeV) of true photon energy.
1923  double srcLogEng = srcEng.log10TeV();
1924 
1925  // Get effectve area components
1926  double nroi = aeff(theta, phi, zenith, azimuth, srcLogEng);
1927 
1928  // Multiply-in PSF
1929  if (nroi > 0.0) {
1930 
1931  // Get PSF
1932  nroi *= npsf(srcDir, srcLogEng, srcTime, pnt, roi);
1933 
1934  // Multiply-in energy dispersion
1935  if (use_edisp() && nroi > 0.0) {
1936 
1937  // Multiply-in energy dispersion
1938  nroi *= edisp(obsEng, srcEng, theta, phi, zenith, azimuth);
1939 
1940  } // endif: had energy dispersion
1941 
1942  // Apply deadtime correction
1943  nroi *= obs.deadc(srcTime);
1944 
1945  } // endif: had non-zero effective area
1946 
1947  // Compile option: Check for NaN/Inf
1948  #if defined(G_NAN_CHECK)
1949  if (gammalib::is_notanumber(nroi) || gammalib::is_infinite(nroi)) {
1950  std::cout << "*** ERROR: GCTAResponseIrf::nroi:";
1951  std::cout << " NaN/Inf encountered";
1952  std::cout << " (nroi=" << nroi;
1953  std::cout << ", theta=" << theta;
1954  std::cout << ", phi=" << phi << ")";
1955  std::cout << std::endl;
1956  }
1957  #endif
1958 
1959  // Return Nroi
1960  return nroi;
1961 }
1962 
1963 
1964 /***********************************************************************//**
1965  * @brief Return result of PSF integration over ROI.
1966  *
1967  * @param[in] srcDir True photon direction.
1968  * @param[in] srcLogEng Log10 of true photon energy (E/TeV).
1969  * @param[in] srcTime True photon arrival time (not used).
1970  * @param[in] pnt CTA pointing.
1971  * @param[in] roi CTA region of interest.
1972  *
1973  * This method integrates the PSF over the circular region of interest (ROI).
1974  * Integration is done in a polar coordinate system centred on the PSF since
1975  * the PSF is assumed to be azimuthally symmetric. The polar integration is
1976  * done using the method npsf_kern_rad_azsym() that computes analytically
1977  * the arclength that is comprised within the ROI.
1978  *
1979  * Note that the integration is only performed when the PSF is spilling out
1980  * of the ROI border, otherwise the integral is simply 1. Numerical
1981  * integration is done using the standard Romberg method. The integration
1982  * boundaries are computed so that only the PSF section that falls in the ROI
1983  * is considered.
1984  *
1985  * @todo Enhance romberg() integration method for small integration regions
1986  * (see comment about kluge below)
1987  * @todo Implement phi dependence in camera system
1988  ***************************************************************************/
1989 double GCTAResponseIrf::npsf(const GSkyDir& srcDir,
1990  const double& srcLogEng,
1991  const GTime& srcTime,
1992  const GCTAPointing& pnt,
1993  const GCTARoi& roi) const
1994 {
1995  // Set number of iterations for Romberg integration.
1996  static const int iter = 6;
1997 
1998  // Declare result
1999  double value = 0.0;
2000 
2001  // Get pointing direction zenith angle and azimuth [radians]
2002  double zenith = pnt.zenith();
2003  double azimuth = pnt.azimuth();
2004 
2005  // Compute offset angle of source direction in camera system
2006  double theta = pnt.dir().dist(srcDir);
2007 
2008  // Compute azimuth angle of source direction in camera system
2009  double phi = 0.0; //TODO: Implement phi dependence
2010 
2011  // Extract relevant parameters from arguments
2012  double roi_radius = roi.radius() * gammalib::deg2rad;
2013  double roi_psf_distance = roi.centre().dir().dist(srcDir);
2014  double rmax = psf_delta_max(theta, phi, zenith, azimuth, srcLogEng);
2015 
2016  // If PSF is fully enclosed by the ROI then skip the numerical
2017  // integration and assume that the integral is 1.0
2018  if (roi_psf_distance + rmax <= roi_radius) {
2019  value = 1.0;
2020  }
2021 
2022  // ... otherwise perform numerical integration
2023  else {
2024 
2025  // Compute minimum PSF integration radius
2026  double rmin = (roi_psf_distance > roi_radius)
2027  ? roi_psf_distance - roi_radius : 0.0;
2028 
2029  // Continue only if integration range is valid
2030  if (rmax > rmin) {
2031 
2032  // Setup integration kernel
2033  cta_npsf_kern_rad_azsym integrand(this,
2034  roi_radius,
2035  roi_psf_distance,
2036  srcLogEng,
2037  theta,
2038  phi,
2039  zenith,
2040  azimuth);
2041 
2042  // Setup integration
2043  GIntegral integral(&integrand);
2044 
2045  // Set fixed number of iterations
2046  integral.fixed_iter(iter);
2047 
2048  // Radially integrate PSF. In case that the radial integration
2049  // region is small, we do the integration using a simple
2050  // trapezoidal rule. This is a kluge to prevent convergence
2051  // problems in the romberg() method for small integration intervals.
2052  // Ideally, the romberg() method should be enhanced to handle this
2053  // case automatically. The kluge threshold was fixed manually!
2054  if (rmax-rmin < 1.0e-12) {
2055  value = integral.trapzd(rmin, rmax);
2056  }
2057  else {
2058  value = integral.romberg(rmin, rmax);
2059  }
2060 
2061  // Compile option: Check for NaN/Inf
2062  #if defined(G_NAN_CHECK)
2063  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
2064  std::cout << "*** ERROR: GCTAResponseIrf::npsf:";
2065  std::cout << " NaN/Inf encountered";
2066  std::cout << " (value=" << value;
2067  std::cout << ", roi_radius=" << roi_radius * gammalib::rad2deg;
2068  std::cout << ", roi_psf_distance=";
2069  std::cout << roi_psf_distance * gammalib::rad2deg;
2070  std::cout << ", r=[" << rmin * gammalib::rad2deg;
2071  std::cout << "," << rmax * gammalib::rad2deg << "])";
2072  std::cout << std::endl;
2073  }
2074  #endif
2075 
2076  } // endif: integration range was valid
2077 
2078  } // endelse: numerical integration required
2079 
2080  // Return integrated PSF
2081  return value;
2082 }
2083 
2084 
2085 /*==========================================================================
2086  = =
2087  = Private methods =
2088  = =
2089  ==========================================================================*/
2090 
2091 /***********************************************************************//**
2092  * @brief Initialise class members
2093  ***************************************************************************/
2095 {
2096  // Initialise members
2097  m_caldb.clear();
2098  m_rspname.clear();
2099  m_aeff = NULL;
2100  m_psf = NULL;
2101  m_edisp = NULL;
2102  m_background = NULL;
2103  m_apply_edisp = false; //!< Switched off by default
2104  m_lo_safe_thres = 0.0;
2105  m_hi_safe_thres = 0.0;
2106 
2107  // XML response filenames
2108  m_xml_caldb.clear();
2109  m_xml_rspname.clear();
2110 
2111  // Return
2112  return;
2113 }
2114 
2115 
2116 /***********************************************************************//**
2117  * @brief Copy class members
2118  *
2119  * @param[in] rsp Response to be copied
2120  ***************************************************************************/
2122 {
2123  // Copy members
2124  m_caldb = rsp.m_caldb;
2125  m_rspname = rsp.m_rspname;
2129 
2130  // Copy response filenames
2131  m_xml_caldb = rsp.m_xml_caldb;
2133 
2134  // Copy nroi cache
2135 
2136  // Clone members
2137  m_aeff = (rsp.m_aeff != NULL) ? rsp.m_aeff->clone() : NULL;
2138  m_psf = (rsp.m_psf != NULL) ? rsp.m_psf->clone() : NULL;
2139  m_edisp = (rsp.m_edisp != NULL) ? rsp.m_edisp->clone() : NULL;
2140  m_background = (rsp.m_background != NULL) ? rsp.m_background->clone() : NULL;
2141 
2142  // Return
2143  return;
2144 }
2145 
2146 
2147 /***********************************************************************//**
2148  * @brief Delete class members
2149  ***************************************************************************/
2151 {
2152  // Free memory
2153  if (m_aeff != NULL) delete m_aeff;
2154  if (m_psf != NULL) delete m_psf;
2155  if (m_edisp != NULL) delete m_edisp;
2156  if (m_background != NULL) delete m_background;
2157 
2158  // Initialise pointers
2159  m_aeff = NULL;
2160  m_psf = NULL;
2161  m_edisp = NULL;
2162  m_background = NULL;
2163 
2164  // Return
2165  return;
2166 }
2167 
2168 
2169 /***********************************************************************//**
2170  * @brief Return filename with appropriate extension
2171  *
2172  * @param[in] filename File name.
2173  * @return File name.
2174  *
2175  * Checks if the specified @p filename exists, and if not, checks whether a
2176  * file with the added suffix .dat exists. Returns the file name with the
2177  * appropriate extension.
2178  ***************************************************************************/
2179 GFilename GCTAResponseIrf::irf_filename(const std::string& filename) const
2180 {
2181  // Set input filename as result filename
2182  GFilename result = filename;
2183 
2184  // If file does not exist then try a variant with extension .dat
2185  if (!result.exists()) {
2186  GFilename testname = filename + ".dat";
2187  if (testname.exists()) {
2188  result = testname;
2189  }
2190  }
2191 
2192  // Return result
2193  return result;
2194 }
2195 
2196 
2197 /***********************************************************************//**
2198  * @brief Return value of point source instrument response function
2199  *
2200  * @param[in] event Observed event.
2201  * @param[in] source Source.
2202  * @param[in] obs Observation.
2203  * @return Value of instrument response function for a point source.
2204  *
2205  * This method returns the value of the instrument response function for a
2206  * point source. The method assumes that source.model() is of type
2207  * GModelSpatialPointSource.
2208  ***************************************************************************/
2210  const GSource& source,
2211  const GObservation& obs) const
2212 {
2213  // Get point source spatial model
2214  const GModelSpatialPointSource* src =
2215  static_cast<const GModelSpatialPointSource*>(source.model());
2216 
2217  // Set Photon
2218  GPhoton photon(src->dir(), source.energy(), source.time());
2219 
2220  // Compute IRF
2221  double irf = this->irf(event, photon, obs);
2222 
2223  // Return IRF
2224  return irf;
2225 }
2226 
2227 
2228 /***********************************************************************//**
2229  * @brief Return IRF value for radial source model
2230  *
2231  * @param[in] event Observed event.
2232  * @param[in] source Source.
2233  * @param[in] obs Observation.
2234  *
2235  * @exception GCTAException::bad_model_type
2236  * Model is not a radial model.
2237  *
2238  * Integrates the product of the spatial model and the instrument response
2239  * function over the true photon arrival direction using
2240  *
2241  * \f[
2242  * \int_{\rho_{\rm min}}^{\rho_{\rm max}}
2243  * \sin \rho \times S_{\rm p}(\rho | E, t) \times
2244  * \int_{\omega_{\rm min}}^{\omega_{\rm max}}
2245  * {\rm Irf}(\rho, \omega) d\omega d\rho
2246  * \f]
2247  *
2248  * where
2249  * \f$S_{\rm p}(\rho | E, t)\f$ is the radial spatial model,
2250  * \f${\rm Irf}(\rho, \omega)\f$ is the instrument response function,
2251  * \f$\rho\f$ is the radial distance from the model centre, and
2252  * \f$\omega\f$ is the position angle with respect to the connecting line
2253  * between the model centre and the observed photon arrival direction.
2254  *
2255  * The integration is performed in the coordinate system of the source
2256  * model spanned by \f$\rho\f$ and \f$\omega\f$ which allows to benefit
2257  * from the symmetry of the source model.
2258  *
2259  * The source centre is located at \f$\vec{m}\f$, and a spherical system
2260  * is defined around this location with \f$(\omega,\rho)\f$ being the
2261  * azimuth and zenith angles, respectively. \f$\omega=0\f$ is defined
2262  * by the direction that connects the source centre \f$\vec{m}\f$ to the
2263  * measured photon direction \f$\vec{p'}\f$, and \f$\omega\f$ increases
2264  * counterclockwise.
2265  *
2266  * Note that this method approximates the true theta angle (angle between
2267  * incident photon and pointing direction) by the measured theta angle
2268  * (angle between the measured photon arrival direction and the pointing
2269  * direction). Given the slow variation of the Psf shape over the field of
2270  * view, this approximation should be fine. It helps in fact a lot in
2271  * speeding up the computations.
2272  ***************************************************************************/
2274  const GSource& source,
2275  const GObservation& obs) const
2276 {
2277  // Retrieve CTA pointing
2278  const GCTAPointing& pnt = gammalib::cta_pnt(G_IRF_RADIAL, obs);
2279  const GCTAInstDir& dir = gammalib::cta_dir(G_IRF_RADIAL, event);
2280 
2281  // Get pointer on radial model
2282  const GModelSpatialRadial* model =
2283  dynamic_cast<const GModelSpatialRadial*>(source.model());
2284  if (model == NULL) {
2286  }
2287 
2288  // Get pointer on shell model (will be NULL for other models)
2289  const GModelSpatialRadialShell* shell =
2290  dynamic_cast<const GModelSpatialRadialShell*>(model);
2291 
2292  // Get pointer on ring model (will be NULL for other models)
2293  const GModelSpatialRadialRing* ring =
2294  dynamic_cast<const GModelSpatialRadialRing*>(model);
2295 
2296  // Set number of iterations for Romberg integration.
2297  // These values have been determined after careful testing, see
2298  // https://cta-redmine.irap.omp.eu/issues/1299
2299  // https://cta-redmine.irap.omp.eu/issues/1521
2300  // https://cta-redmine.irap.omp.eu/issues/2860
2301  static const int iter_rho = 6;
2302  static const int iter_phi = 6;
2303 
2304  // Get event attributes
2305  const GSkyDir& obsDir = dir.dir();
2306  const GEnergy& obsEng = event.energy();
2307 
2308  // Get source attributes
2309  const GSkyDir& centre = model->dir();
2310  const GEnergy& srcEng = source.energy();
2311  const GTime& srcTime = source.time();
2312 
2313  // Get pointing direction zenith angle and azimuth [radians]
2314  double zenith = pnt.zenith();
2315  double azimuth = pnt.azimuth();
2316 
2317  // Determine angular distance between measured photon direction and model
2318  // centre [radians]
2319  double zeta = centre.dist(obsDir);
2320 
2321  // Determine angular distance between measured photon direction and
2322  // pointing direction [radians]
2323  double eta = pnt.dir().dist(obsDir);
2324 
2325  // Determine angular distance between model centre and pointing direction
2326  // [radians]
2327  double lambda = centre.dist(pnt.dir());
2328 
2329  // Compute azimuth angle of pointing in model system [radians]
2330  // Will be comprised in interval [0,pi]
2331  double omega0 = 0.0;
2332  double denom = std::sin(lambda) * std::sin(zeta);
2333  if (denom != 0.0) {
2334  double arg = (std::cos(eta) - std::cos(lambda) * std::cos(zeta))/denom;
2335  omega0 = gammalib::acos(arg);
2336  }
2337 
2338  // Get log10(E/TeV) of true photon energy
2339  double srcLogEng = srcEng.log10TeV();
2340 
2341  // Assign the observed theta angle (eta) as the true theta angle
2342  // between the source and the pointing directions. This is a (not
2343  // too bad) approximation which helps to speed up computations.
2344  // If we want to do this correctly, however, we would need to move
2345  // the psf_dummy_sigma down to the integration kernel, and we would
2346  // need to make sure that psf_delta_max really gives the absolute
2347  // maximum (this is certainly less critical)
2348  double theta = eta;
2349  double phi = 0.0; //TODO: Implement IRF Phi dependence
2350 
2351  // Get maximum PSF and source radius in radians.
2352  double delta_max = psf_delta_max(theta, phi, zenith, azimuth, srcLogEng);
2353  double src_max = model->theta_max();
2354 
2355  // Set radial model zenith angle range
2356  double rho_min = (zeta > delta_max) ? zeta - delta_max : 0.0;
2357  double rho_max = zeta + delta_max;
2358  if (rho_max > src_max) {
2359  rho_max = src_max;
2360  }
2361 
2362  // Initialise IRF value
2363  double irf = 0.0;
2364 
2365  // Perform zenith angle integration if interval is valid
2366  if (rho_max > rho_min) {
2367 
2368  // Setup integration kernel
2369  cta_irf_radial_kern_rho integrand(this,
2370  model,
2371  zenith,
2372  azimuth,
2373  srcEng,
2374  srcTime,
2375  obsEng,
2376  zeta,
2377  lambda,
2378  omega0,
2379  delta_max,
2380  iter_phi);
2381 
2382  // Integrate over model's zenith angle
2383  GIntegral integral(&integrand);
2384  integral.fixed_iter(iter_rho);
2385 
2386  // Setup integration boundaries
2387  std::vector<double> bounds;
2388  bounds.push_back(rho_min);
2389  bounds.push_back(rho_max);
2390 
2391  // If the integration range includes a transition between full
2392  // containment of model within Psf and partial containment, then
2393  // add a boundary at this location
2394  double transition_point = delta_max - zeta;
2395  if (transition_point > rho_min && transition_point < rho_max) {
2396  bounds.push_back(transition_point);
2397  }
2398 
2399  // If we have a shell model then add an integration boundary for the
2400  // shell radius as a function discontinuity will occur at this
2401  // location
2402  if (shell != NULL) {
2403  double shell_radius = shell->radius() * gammalib::deg2rad;
2404  if (shell_radius > rho_min && shell_radius < rho_max) {
2405  bounds.push_back(shell_radius);
2406  }
2407  }
2408 
2409  // If we have a ring model then add an integration boundary for the
2410  // ring radius as a function discontinuity will occur at this
2411  // location
2412  if (ring != NULL) {
2413  double ring_radius = ring->radius() * gammalib::deg2rad;
2414  if (ring_radius > rho_min && ring_radius < rho_max) {
2415  bounds.push_back(ring_radius);
2416  }
2417  }
2418 
2419  // Integrate kernel
2420  irf = integral.romberg(bounds, iter_rho);
2421 
2422  // Compile option: Check for NaN/Inf
2423  #if defined(G_NAN_CHECK)
2425  std::cout << "*** ERROR: GCTAResponseIrf::irf_radial:";
2426  std::cout << " NaN/Inf encountered";
2427  std::cout << " (irf=" << irf;
2428  std::cout << ", rho_min=" << rho_min;
2429  std::cout << ", rho_max=" << rho_max;
2430  std::cout << ", omega0=" << omega0 << ")";
2431  std::cout << std::endl;
2432  }
2433  #endif
2434 
2435  // Apply deadtime correction
2436  irf *= obs.deadc(srcTime);
2437 
2438  } // endif: integration interval is valid
2439 
2440  // Compile option: Show integration results
2441  #if defined(G_DEBUG_IRF_RADIAL)
2442  std::cout << "GCTAResponseIrf::irf_radial:";
2443  std::cout << " rho=[" << rho_min*gammalib::rad2deg;
2444  std::cout << ", " << rho_max*gammalib::rad2deg << " deg;";
2445  std::cout << " zeta=" << zeta*gammalib::rad2deg << " deg;";
2446  std::cout << " eta=" << eta*gammalib::rad2deg << " deg;";
2447  std::cout << " lambda=" << lambda*gammalib::rad2deg << " deg;";
2448  std::cout << " omega0=" << omega0*gammalib::rad2deg << " deg;";
2449  std::cout << " theta_max=" << src_max*gammalib::rad2deg << " deg;";
2450  std::cout << " delta_max=" << delta_max*gammalib::rad2deg << " deg;";
2451  std::cout << " obsDir=" << obsDir << ";";
2452  std::cout << " modelDir=" << model->dir() << ";";
2453  std::cout << " irf=" << irf << std::endl;
2454  #endif
2455 
2456  // Return IRF value
2457  return irf;
2458 }
2459 
2460 
2461 /***********************************************************************//**
2462  * @brief Return Irf value for elliptical source model
2463  *
2464  * @param[in] event Observed event.
2465  * @param[in] source Source.
2466  * @param[in] obs Observation.
2467  *
2468  * @exception GCTAException::bad_instdir_type
2469  * Instrument direction is not a valid CTA instrument direction.
2470  * @exception GCTAException::bad_model_type
2471  * Model is not an elliptical model.
2472  *
2473  * Integrates the product of the model and the IRF over the true photon
2474  * arrival direction using
2475  *
2476  * \f[
2477  * \int_{\rho_{\rm min}}^{\rho_{\rm max}}
2478  * \sin \rho \times
2479  * \int_{\omega}
2480  * S_{\rm p}(\rho, \omega | E, t) \, IRF(\rho, \omega) d\omega d\rho
2481  * \f]
2482  *
2483  * where
2484  * \f$S_{\rm p}(\rho, \omega | E, t)\f$ is the elliptical model,
2485  * \f$IRF(\rho, \omega)\f$ is the instrument response function,
2486  * \f$\rho\f$ is the distance from the model centre, and
2487  * \f$\omega\f$ is the position angle with respect to the connecting line
2488  * between the model centre and the observed photon arrival direction.
2489  *
2490  * The source model centre is located at \f$\vec{m}\f$, and a spherical
2491  * coordinate system is defined around this location with \f$(\rho,\omega)\f$
2492  * being the zenith and azimuth angles, respectively. The azimuth angle
2493  * \f$(\omega)\f$ is counted counterclockwise from the vector that runs from
2494  * the model centre \f$\vec{m}\f$ to the measured photon direction
2495  * \f$\vec{p'}\f$.
2496  ***************************************************************************/
2498  const GSource& source,
2499  const GObservation& obs) const
2500 {
2501  // Set number of iterations for Romberg integration.
2502  // These values have been determined after careful testing, see
2503  // https://cta-redmine.irap.omp.eu/issues/1299
2504  // https://cta-redmine.irap.omp.eu/issues/1521
2505  // https://cta-redmine.irap.omp.eu/issues/2860
2506  static const int iter_rho = 6;
2507  static const int iter_phi = 6;
2508 
2509  // Retrieve CTA pointing
2510  const GCTAPointing& pnt = gammalib::cta_pnt(G_IRF_ELLIPTICAL, obs);
2511  const GCTAInstDir& dir = gammalib::cta_dir(G_IRF_ELLIPTICAL, event);
2512 
2513  // Get pointer on elliptical model
2514  const GModelSpatialElliptical* model =
2515  dynamic_cast<const GModelSpatialElliptical*>(source.model());
2516  if (model == NULL) {
2518  }
2519 
2520  // Get event attributes (measured photon)
2521  const GSkyDir& obsDir = dir.dir();
2522  const GEnergy& obsEng = event.energy();
2523 
2524  // Get source attributes
2525  const GSkyDir& centre = model->dir();
2526  const GEnergy& srcEng = source.energy();
2527  const GTime& srcTime = source.time();
2528 
2529  // Get pointing direction zenith angle and azimuth [radians]
2530  double zenith = pnt.zenith();
2531  double azimuth = pnt.azimuth();
2532 
2533  // Determine angular distance between observed photon direction and model
2534  // centre and position angle of observed photon direction seen from the
2535  // model centre [radians]
2536  double rho_obs = centre.dist(obsDir);
2537  double posangle_obs = centre.posang(obsDir); // Celestial system
2538 
2539  // Determine angular distance between model centre and pointing direction
2540  // [radians]
2541  double rho_pnt = centre.dist(pnt.dir());
2542  double posangle_pnt = centre.posang(pnt.dir()); // Celestial system
2543 
2544  // Compute azimuth angle of pointing in model coordinate system [radians]
2545  double omega_pnt = posangle_pnt - posangle_obs;
2546 
2547  // Get log10(E/TeV) of true photon energy
2548  double srcLogEng = srcEng.log10TeV();
2549 
2550  // Get maximum PSF radius [radians]. We assign here the measured theta
2551  // angle (eta) as the true theta angle between the source and the pointing
2552  // directions. As we only use the angle to determine the maximum PSF size,
2553  // this should be sufficient.
2554  double theta = pnt.dir().dist(obsDir);
2555  double phi = 0.0; //TODO: Implement IRF Phi dependence
2556  double delta_max = psf_delta_max(theta, phi, zenith, azimuth, srcLogEng);
2557 
2558  // Get the ellipse boundary (radians). Note that these are NOT the
2559  // parameters of the ellipse but of a boundary ellipse that is used
2560  // for computing the relevant omega angle intervals for a given angle
2561  // rho. The boundary ellipse takes care of the possibility that the
2562  // semiminor axis is larger than the semimajor axis
2563  double semimajor; // Will be the larger axis
2564  double semiminor; // Will be the smaller axis
2565  double posangle; // Will be the corrected position angle
2566  double aspect_ratio; // Ratio between smaller/larger axis of model
2567  if (model->semimajor() >= model->semiminor()) {
2568  aspect_ratio = (model->semimajor() > 0.0) ?
2569  model->semiminor() / model->semimajor() : 0.0;
2570  posangle = model->posangle() * gammalib::deg2rad;
2571  }
2572  else {
2573  aspect_ratio = (model->semiminor() > 0.0) ?
2574  model->semimajor() / model->semiminor() : 0.0;
2575  posangle = model->posangle() * gammalib::deg2rad + gammalib::pihalf;
2576  }
2577  semimajor = model->theta_max();
2578  semiminor = semimajor * aspect_ratio;
2579 
2580  // Set zenith angle integration range for elliptical model
2581  double rho_min = (rho_obs > delta_max) ? rho_obs - delta_max : 0.0;
2582  double rho_max = rho_obs + delta_max;
2583  if (rho_max > semimajor) {
2584  rho_max = semimajor;
2585  }
2586 
2587  // Initialise IRF value
2588  double irf = 0.0;
2589 
2590  // Perform zenith angle integration if interval is valid
2591  if (rho_max > rho_min) {
2592 
2593  // Setup integration kernel
2594  cta_irf_elliptical_kern_rho integrand(this,
2595  model,
2596  semimajor,
2597  semiminor,
2598  posangle,
2599  zenith,
2600  azimuth,
2601  srcEng,
2602  srcTime,
2603  obsEng,
2604  rho_obs,
2605  posangle_obs,
2606  rho_pnt,
2607  omega_pnt,
2608  delta_max,
2609  iter_phi);
2610 
2611  // Integrate over model's zenith angle
2612  GIntegral integral(&integrand);
2613  integral.fixed_iter(iter_rho);
2614 
2615  // Setup integration boundaries
2616  std::vector<double> bounds;
2617  bounds.push_back(rho_min);
2618  bounds.push_back(rho_max);
2619 
2620  // If the integration range includes the semiminor boundary, then
2621  // add an integration boundary at that location
2622  if (semiminor > rho_min && semiminor < rho_max) {
2623  bounds.push_back(semiminor);
2624  }
2625 
2626  // Integrate kernel
2627  irf = integral.romberg(bounds, iter_rho);
2628 
2629  // Compile option: Check for NaN/Inf
2630  #if defined(G_NAN_CHECK)
2632  std::cout << "*** ERROR: GCTAResponseIrf::irf_elliptical:";
2633  std::cout << " NaN/Inf encountered";
2634  std::cout << " (irf=" << irf;
2635  std::cout << ", rho_min=" << rho_min;
2636  std::cout << ", rho_max=" << rho_max << ")";
2637  std::cout << std::endl;
2638  }
2639  #endif
2640 
2641  // Apply deadtime correction
2642  irf *= obs.deadc(srcTime);
2643 
2644  } // endif: integration interval is valid
2645 
2646  // Compile option: Show integration results
2647  #if defined(G_DEBUG_IRF_ELLIPTICAL)
2648  std::cout << "GCTAResponseIrf::irf_elliptical:";
2649  std::cout << " rho_min=" << rho_min;
2650  std::cout << " rho_max=" << rho_max;
2651  std::cout << " irf=" << irf << std::endl;
2652  #endif
2653 
2654  // Return IRF value
2655  return irf;
2656 }
2657 
2658 
2659 /***********************************************************************//**
2660  * @brief Return value of diffuse source instrument response function
2661  *
2662  * @param[in] event Observed event.
2663  * @param[in] source Source.
2664  * @param[in] obs Observation.
2665  *
2666  * Integrates the product of the model and the IRF over the true photon
2667  * arrival direction using
2668  *
2669  * \f[
2670  * \int_{0}^{\theta_{\rm max}}
2671  * \sin \theta \times PSF(\theta)
2672  * \int_{0}^{2\pi}
2673  * S_{\rm p}(\theta, \phi | E, t) \,
2674  * Aeff(\theta, \phi) \,
2675  * Edisp(\theta, \phi) d\phi
2676  * \f]
2677  *
2678  * where
2679  * - \f$S_{\rm p}(\theta, \phi | E, t)\f$ is the diffuse model,
2680  * - \f$PSF(\theta)\f$ is the azimuthally symmetric Point Spread Function,
2681  * - \f$Aeff(\theta, \phi)\f$ is the effective area,
2682  * - \f$Edisp(\theta, \phi)\f$ is the energy dispersion,
2683  * - \f$\theta\f$ is the distance from the PSF centre, and
2684  * - \f$\phi\f$ is the azimuth angle.
2685  *
2686  * The integration is performed in the reference of the observed arrival
2687  * direction. Integration is done first over the azimuth angle \f$\phi\f$ and
2688  * then over the offset angle \f$\theta\f$.
2689  *
2690  * The integration kernels for this method are implemented by the response
2691  * helper classes cta_irf_diffuse_kern_theta and cta_irf_diffuse_kern_phi.
2692  ***************************************************************************/
2694  const GSource& source,
2695  const GObservation& obs) const
2696 {
2697  // Set minimum and maximum number of iterations for Romberg integration.
2698  // These values have been determined after careful testing, see
2699  // https://cta-redmine.irap.omp.eu/issues/1248
2700  // https://cta-redmine.irap.omp.eu/issues/2458
2701  // https://cta-redmine.irap.omp.eu/issues/2860
2702  static const int min_iter_rho = 6;
2703  static const int min_iter_phi = 6;
2704  static const int max_iter_rho = 8;
2705  static const int max_iter_phi = 8;
2706 
2707  // Initialise IRF value
2708  double irf = 0.0;
2709 
2710  // Retrieve CTA pointing
2711  const GCTAPointing& pnt = gammalib::cta_pnt(G_IRF_DIFFUSE, obs);
2712 
2713  // Get CTA instrument direction
2714  const GCTAInstDir& dir = gammalib::cta_dir(G_IRF_DIFFUSE, event);
2715 
2716  // Get pointer on spatial model
2717  const GModelSpatial* model =
2718  dynamic_cast<const GModelSpatial*>(source.model());
2719  if (model == NULL) {
2721  }
2722 
2723  // Get resolution of spatial model
2724  double resolution = gammalib::resolution(model);
2725 
2726  // Get event attributes
2727  const GEnergy& obsEng = event.energy();
2728 
2729  // Get source attributes
2730  const GEnergy& srcEng = source.energy();
2731  const GTime& srcTime = source.time();
2732 
2733  // Get pointing direction zenith angle and azimuth [radians]
2734  double zenith = pnt.zenith();
2735  double azimuth = pnt.azimuth();
2736 
2737  // Determine angular distance between measured photon direction and
2738  // pointing direction [radians]
2739  double eta = pnt.dir().dist(dir.dir());
2740 
2741  // Get log10(E/TeV) of true photon energy
2742  double srcLogEng = srcEng.log10TeV();
2743 
2744  // Assign the observed theta angle (eta) as the true theta angle
2745  // between the source and the pointing directions. This is a (not
2746  // too bad) approximation which helps to speed up computations.
2747  // If we want to do this correctly, however, we would need to move
2748  // the psf_dummy_sigma down to the integration kernel, and we would
2749  // need to make sure that psf_delta_max really gives the absolute
2750  // maximum (this is certainly less critical)
2751  double theta = eta;
2752  double phi = 0.0; //TODO: Implement Phi dependence
2753 
2754  // Get maximum PSF radius in radians
2755  double delta_max = psf_delta_max(theta, phi, zenith, azimuth, srcLogEng);
2756 
2757  // Perform zenith angle integration if interval is valid
2758  if (delta_max > 0.0) {
2759 
2760  // Compute rotation matrix to convert from coordinates (theta,phi)
2761  // in the reference frame of the observed arrival direction into
2762  // celestial coordinates
2763  GMatrix ry;
2764  GMatrix rz;
2765  ry.eulery(dir.dir().dec_deg() - 90.0);
2766  rz.eulerz(-dir.dir().ra_deg());
2767  GMatrix rot = (ry * rz).transpose();
2768 
2769  // Setup integration kernel
2770  cta_irf_diffuse_kern_theta integrand(this,
2771  model,
2772  &rot,
2773  theta,
2774  phi,
2775  zenith,
2776  azimuth,
2777  srcEng,
2778  srcTime,
2779  srcLogEng,
2780  obsEng,
2781  eta,
2782  min_iter_phi,
2783  max_iter_phi,
2784  resolution);
2785 
2786  // Set number of radial integration iterations
2787  int iter = gammalib::iter_rho(delta_max, resolution,
2788  min_iter_rho, max_iter_rho);
2789 
2790  // Integrate over Psf delta angle
2791  GIntegral integral(&integrand);
2792  integral.fixed_iter(iter);
2793  irf = integral.romberg(0.0, delta_max);
2794 
2795  // Compile option: Check for NaN/Inf
2796  #if defined(G_NAN_CHECK)
2798  std::cout << "*** ERROR: GCTAResponseIrf::irf_diffuse:";
2799  std::cout << " NaN/Inf encountered";
2800  std::cout << " (irf=" << irf;
2801  std::cout << ", delta_max=" << delta_max << ")";
2802  std::cout << std::endl;
2803  }
2804  #endif
2805  }
2806 
2807  // Apply deadtime correction
2808  irf *= obs.deadc(srcTime);
2809 
2810  // Compile option: Show integration results
2811  #if defined(G_DEBUG_IRF_DIFFUSE)
2812  std::cout << "GCTAResponseIrf::irf_diffuse:";
2813  std::cout << " srcLogEng=" << srcLogEng;
2814  std::cout << " obsLogEng=" << obsLogEng;
2815  std::cout << " eta=" << eta;
2816  std::cout << " delta_max=" << delta_max;
2817  std::cout << " irf=" << irf << std::endl;
2818  #endif
2819 
2820  // Return IRF value
2821  return irf;
2822 }
2823 
2824 
2825 /***********************************************************************//**
2826  * @brief Return spatial integral of point source model
2827  *
2828  * @param[in] model Sky Model.
2829  * @param[in] srcEng True photon energy.
2830  * @param[in] srcTime True photon arrival time.
2831  * @param[in] obsEng Observed event energy.
2832  * @param[in] obsTime Observed event arrival time.
2833  * @param[in] obs Observation.
2834  *
2835  * Computes the integral
2836  *
2837  * \f[
2838  * N_{\rm ROI}(E',t'|E,t) = \int_{\rm ROI} P(p',E',t'|E,t) dp'
2839  * \f]
2840  *
2841  * of
2842  *
2843  * \f[
2844  * P(p',E',t'|E,t) = \int
2845  * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp
2846  * \f]
2847  *
2848  * over the Region of Interest (ROI) for a point source model \f$S(p,E,t)\f$
2849  * and the response function \f$R(p',E',t'|p,E,t)\f$.
2850  ***************************************************************************/
2852  const GEnergy& srcEng,
2853  const GTime& srcTime,
2854  const GEnergy& obsEng,
2855  const GTime& obsTime,
2856  const GObservation& obs) const
2857 {
2858  // Get point source spatial model
2859  const GModelSpatialPointSource* src =
2860  static_cast<const GModelSpatialPointSource*>(model.spatial());
2861 
2862  // Set Photon
2863  GPhoton photon(src->dir(), srcEng, srcTime);
2864 
2865  // Compute Nroi
2866  double nroi = nirf(photon, obsEng, obsTime, obs);
2867 
2868  // Return Nroi
2869  return nroi;
2870 }
2871 
2872 
2873 /***********************************************************************//**
2874  * @brief Return spatial integral of radial source model
2875  *
2876  * @param[in] model Sky Model.
2877  * @param[in] srcEng True photon energy.
2878  * @param[in] srcTime True photon arrival time.
2879  * @param[in] obsEng Observed event energy.
2880  * @param[in] obsTime Observed event arrival time.
2881  * @param[in] obs Observation.
2882  *
2883  * Computes the integral
2884  *
2885  * \f[
2886  * N_{\rm ROI}(E',t'|E,t) = \int_{\rm ROI} P(p',E',t'|E,t) dp'
2887  * \f]
2888  *
2889  * of
2890  *
2891  * \f[
2892  * P(p',E',t'|E,t) = \int
2893  * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp
2894  * \f]
2895  *
2896  * over the Region of Interest (ROI) for a radial source model \f$S(p,E,t)\f$
2897  * and the response function \f$R(p',E',t'|p,E,t)\f$.
2898  ***************************************************************************/
2900  const GEnergy& srcEng,
2901  const GTime& srcTime,
2902  const GEnergy& obsEng,
2903  const GTime& obsTime,
2904  const GObservation& obs) const
2905 {
2906  // Set number of iterations for Romberg integration.
2907  // These values have been determined after careful testing, see
2908  // https://cta-redmine.irap.omp.eu/issues/1299
2909  static const int iter_rho = 6;
2910  static const int iter_phi = 6;
2911 
2912  // Initialise Nroi value
2913  double nroi = 0.0;
2914 
2915  // Retrieve CTA observation, ROI and pointing
2916  const GCTAObservation& cta = gammalib::cta_obs(G_NROI_RADIAL, obs);
2917  const GCTARoi& roi = gammalib::cta_event_list(G_NROI_RADIAL, obs).roi();
2918  const GCTAPointing& pnt = cta.pointing();
2919 
2920  // Get pointer on radial model
2921  const GModelSpatialRadial* spatial =
2922  dynamic_cast<const GModelSpatialRadial*>(model.spatial());
2923  if (spatial == NULL) {
2925  }
2926 
2927  // Get source attributes
2928  const GSkyDir& centre = spatial->dir();
2929 
2930  // Get pointing direction zenith angle and azimuth [radians]
2931  double zenith = pnt.zenith();
2932  double azimuth = pnt.azimuth();
2933 
2934  // Get log10(E/TeV) of true photon energy
2935  double srcLogEng = srcEng.log10TeV();
2936 
2937  // Get maximum PSF radius (radians). We do this for the onaxis PSF only,
2938  // as this allows us doing this computation in the outer loop. This
2939  // should be sufficient here, unless the offaxis PSF becomes much worse
2940  // than the onaxis PSF. In this case, we may add a safety factor here
2941  // to make sure we encompass the entire PSF.
2942  double psf_max_radius = psf_delta_max(0.0, 0.0, zenith, azimuth, srcLogEng);
2943 
2944  // Extract ROI radius (radians)
2945  double roi_radius = roi.radius() * gammalib::deg2rad;
2946 
2947  // Compute distance between ROI and model centre (radians)
2948  double roi_model_distance = roi.centre().dir().dist(centre);
2949 
2950  // Compute the ROI radius plus maximum PSF radius (radians). Any photon
2951  // coming from beyond this radius will not make it in the dataspace and
2952  // thus can be neglected.
2953  double roi_psf_radius = roi_radius + psf_max_radius;
2954 
2955  // Set offset angle integration range. We take here the ROI+PSF into
2956  // account to make no integrations beyond the point where the
2957  // contribution drops to zero.
2958  double rho_min = (roi_model_distance > roi_psf_radius)
2959  ? roi_model_distance - roi_psf_radius: 0.0;
2960  double rho_max = spatial->theta_max();
2961 
2962  // Perform offset angle integration only if interval is valid
2963  if (rho_max > rho_min) {
2964 
2965  // Compute rotation matrix to convert from native model coordinates,
2966  // given by (rho,omega), into celestial coordinates.
2967  GMatrix ry;
2968  GMatrix rz;
2969  ry.eulery(spatial->dec() - 90.0);
2970  rz.eulerz(-spatial->ra());
2971  GMatrix rot = (ry * rz).transpose();
2972 
2973  // Compute position angle of ROI centre with respect to model
2974  // centre (radians)
2975  double omega0 = centre.posang(roi.centre().dir()); // Celestial system
2976 
2977  // Setup integration kernel
2978  cta_nroi_radial_kern_rho integrand(this,
2979  &cta,
2980  spatial,
2981  &rot,
2982  srcEng,
2983  srcTime,
2984  obsEng,
2985  obsTime,
2986  roi_model_distance,
2987  roi_psf_radius,
2988  omega0,
2989  iter_phi);
2990 
2991  // Integrate over model's zenith angle
2992  GIntegral integral(&integrand);
2993  integral.fixed_iter(iter_rho);
2994 
2995  // Setup integration boundaries
2996  std::vector<double> bounds;
2997  bounds.push_back(rho_min);
2998  bounds.push_back(rho_max);
2999 
3000  // If the integration range includes a transition between full
3001  // containment of model within ROI and partial containment, then
3002  // add a boundary at this location
3003  double transition_point = roi_psf_radius - roi_model_distance;
3004  if (transition_point > rho_min && transition_point < rho_max) {
3005  bounds.push_back(transition_point);
3006  }
3007 
3008  // If the integration range includes a transition between full
3009  // containment of ROI within model and partial containment, then
3010  // add a boundary at this location
3011  transition_point = roi_psf_radius + roi_model_distance;
3012  if (transition_point > rho_min && transition_point < rho_max) {
3013  bounds.push_back(transition_point);
3014  }
3015 
3016  // If we have a shell model then add an integration boundary for the
3017  // shell radius as a function discontinuity will occur at this
3018  // location
3019  const GModelSpatialRadialShell* shell = dynamic_cast<const GModelSpatialRadialShell*>(spatial);
3020  if (shell != NULL) {
3021  double shell_radius = shell->radius() * gammalib::deg2rad;
3022  if (shell_radius > rho_min && shell_radius < rho_max) {
3023  bounds.push_back(shell_radius);
3024  }
3025  }
3026 
3027  // If we have a ring model then add an integration boundary for the
3028  // ring radius as a function discontinuity will occur at this
3029  // location
3030  const GModelSpatialRadialRing* ring = dynamic_cast<const GModelSpatialRadialRing*>(spatial);
3031  if (ring != NULL) {
3032  double ring_radius = ring->radius() * gammalib::deg2rad;
3033  if (ring_radius > rho_min && ring_radius < rho_max) {
3034  bounds.push_back(ring_radius);
3035  }
3036  }
3037 
3038  // Integrate kernel
3039  nroi = integral.romberg(bounds, iter_rho);
3040 
3041  // Compile option: Show integration results
3042  #if defined(G_DEBUG_NROI_RADIAL)
3043  std::cout << "GCTAResponseIrf::nroi_radial:";
3044  std::cout << " rho_min=" << rho_min;
3045  std::cout << " rho_max=" << rho_max;
3046  std::cout << " nroi=" << nroi << std::endl;
3047  #endif
3048 
3049  } // endif: offset angle range was valid
3050 
3051  // Debug: Check for NaN
3052  #if defined(G_NAN_CHECK)
3053  if (gammalib::is_notanumber(nroi) || gammalib::is_infinite(nroi)) {
3054  std::cout << "*** ERROR: GCTAResponseIrf::nroi_radial:";
3055  std::cout << " NaN/Inf encountered";
3056  std::cout << " (nroi=" << nroi;
3057  std::cout << ", rho_min=" << rho_min;
3058  std::cout << ", rho_max=" << rho_max;
3059  std::cout << ")" << std::endl;
3060  }
3061  #endif
3062 
3063  // Return Nroi
3064  return nroi;
3065 }
3066 
3067 
3068 /***********************************************************************//**
3069  * @brief Return spatial integral of elliptical source model
3070  *
3071  * @param[in] model Sky Model.
3072  * @param[in] srcEng True photon energy.
3073  * @param[in] srcTime True photon arrival time.
3074  * @param[in] obsEng Observed event energy.
3075  * @param[in] obsTime Observed event arrival time.
3076  * @param[in] obs Observation.
3077  *
3078  * Computes the integral
3079  *
3080  * \f[
3081  * N_{\rm ROI}(E',t'|E,t) = \int_{\rm ROI} P(p',E',t'|E,t) dp'
3082  * \f]
3083  *
3084  * of
3085  *
3086  * \f[
3087  * P(p',E',t'|E,t) = \int
3088  * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp
3089  * \f]
3090  *
3091  * over the Region of Interest (ROI) for an elliptical source model
3092  * \f$S(p,E,t)\f$ and the response function \f$R(p',E',t'|p,E,t)\f$.
3093  ***************************************************************************/
3095  const GEnergy& srcEng,
3096  const GTime& srcTime,
3097  const GEnergy& obsEng,
3098  const GTime& obsTime,
3099  const GObservation& obs) const
3100 {
3101  // Set number of iterations for Romberg integration.
3102  // These values have been determined after careful testing, see
3103  // https://cta-redmine.irap.omp.eu/issues/1299
3104  static const int iter_rho = 6;
3105  static const int iter_phi = 6;
3106 
3107  // Initialise Nroi value
3108  double nroi = 0.0;
3109 
3110  // Retrieve CTA observation, ROI and pointing
3113  const GCTAPointing& pnt = cta.pointing();
3114 
3115  // Get pointer on elliptical model
3116  const GModelSpatialElliptical* spatial =
3117  dynamic_cast<const GModelSpatialElliptical*>(model.spatial());
3118  if (spatial == NULL) {
3120  }
3121 
3122  // Get source attributes
3123  const GSkyDir& centre = spatial->dir();
3124 
3125  // Get pointing direction zenith angle and azimuth [radians]
3126  double zenith = pnt.zenith();
3127  double azimuth = pnt.azimuth();
3128 
3129  // Get log10(E/TeV) of true photon energy
3130  double srcLogEng = srcEng.log10TeV();
3131 
3132  // Get maximum PSF radius (radians). We do this for the onaxis PSF only,
3133  // as this allows us doing this computation in the outer loop. This
3134  // should be sufficient here, unless the offaxis PSF becomes much worse
3135  // than the onaxis PSF. In this case, we may add a safety factor here
3136  // to make sure we encompass the entire PSF.
3137  double psf_max_radius = psf_delta_max(0.0, 0.0, zenith, azimuth, srcLogEng);
3138 
3139  // Extract ROI radius plus maximum PSF radius (radians). Any photon
3140  // coming from beyond this radius will not make it in the dataspace and
3141  // thus can be neglected.
3142  double radius_roi = roi.radius() * gammalib::deg2rad + psf_max_radius;
3143 
3144  // Compute distance between ROI and model centre (radians)
3145  double rho_roi = roi.centre().dir().dist(centre);
3146 
3147  // Get the ellipse boundary (radians). Note that these are NOT the
3148  // parameters of the ellipse but of a boundary ellipse that is used
3149  // for computing the relevant omega angle intervals for a given angle
3150  // rho. The boundary ellipse takes care of the possibility that the
3151  // semiminor axis is larger than the semimajor axis
3152  double semimajor; // Will be the larger axis
3153  double semiminor; // Will be the smaller axis
3154  double posangle; // Will be the corrected position angle
3155  double aspect_ratio; // Ratio between smaller/larger axis of model
3156  if (spatial->semimajor() >= spatial->semiminor()) {
3157  aspect_ratio = (spatial->semimajor() > 0.0) ?
3158  spatial->semiminor() / spatial->semimajor() : 0.0;
3159  posangle = spatial->posangle() * gammalib::deg2rad;
3160  }
3161  else {
3162  aspect_ratio = (spatial->semiminor() > 0.0) ?
3163  spatial->semimajor() / spatial->semiminor() : 0.0;
3164  posangle = spatial->posangle() * gammalib::deg2rad + gammalib::pihalf;
3165  }
3166  semimajor = spatial->theta_max();
3167  semiminor = semimajor * aspect_ratio;
3168 
3169  // Set offset angle integration range. We take here the ROI+PSF into
3170  // account to make no integrations beyond the point where the
3171  // contribution drops to zero.
3172  double rho_min = (rho_roi > radius_roi) ? rho_roi - radius_roi: 0.0;
3173  double rho_max = rho_roi + radius_roi;
3174  if (rho_max > semimajor) {
3175  rho_max = semimajor;
3176  }
3177 
3178  // Perform offset angle integration only if interval is valid
3179  if (rho_max > rho_min) {
3180 
3181  // Compute rotation matrix to convert from native model coordinates,
3182  // given by (rho,omega), into celestial coordinates.
3183  GMatrix ry;
3184  GMatrix rz;
3185  ry.eulery(spatial->dec() - 90.0);
3186  rz.eulerz(-spatial->ra());
3187  GMatrix rot = (ry * rz).transpose();
3188 
3189  // Compute position angle of ROI centre with respect to model
3190  // centre (radians)
3191  double posangle_roi = centre.posang(roi.centre().dir()); // Celestial system
3192 
3193  // Setup integration kernel
3194  cta_nroi_elliptical_kern_rho integrand(this,
3195  &cta,
3196  spatial,
3197  &rot,
3198  semimajor,
3199  semiminor,
3200  posangle,
3201  srcEng,
3202  srcTime,
3203  obsEng,
3204  obsTime,
3205  rho_roi,
3206  posangle_roi,
3207  radius_roi,
3208  iter_phi);
3209 
3210  // Integrate over model's zenith angle
3211  GIntegral integral(&integrand);
3212  integral.fixed_iter(iter_rho);
3213 
3214  // Setup integration boundaries
3215  std::vector<double> bounds;
3216  bounds.push_back(rho_min);
3217  bounds.push_back(rho_max);
3218 
3219  // If the integration range includes the semiminor boundary, then
3220  // add an integration boundary at that location
3221  if (semiminor > rho_min && semiminor < rho_max) {
3222  bounds.push_back(semiminor);
3223  }
3224 
3225  // Integrate kernel
3226  nroi = integral.romberg(bounds, iter_rho);
3227 
3228  // Compile option: Show integration results
3229  #if defined(G_DEBUG_NROI_ELLIPTICAL)
3230  std::cout << "GCTAResponseIrf::nroi_elliptical:";
3231  std::cout << " rho_min=" << rho_min;
3232  std::cout << " rho_max=" << rho_max;
3233  std::cout << " nroi=" << nroi << std::endl;
3234  #endif
3235 
3236  } // endif: offset angle range was valid
3237 
3238  // Debug: Check for NaN
3239  #if defined(G_NAN_CHECK)
3240  if (gammalib::is_notanumber(nroi) || gammalib::is_infinite(nroi)) {
3241  std::cout << "*** ERROR: GCTAResponseIrf::nroi_elliptical:";
3242  std::cout << " NaN/Inf encountered";
3243  std::cout << " (nroi=" << nroi;
3244  std::cout << ", rho_min=" << rho_min;
3245  std::cout << ", rho_max=" << rho_max;
3246  std::cout << ")" << std::endl;
3247  }
3248  #endif
3249 
3250  // Return Nroi
3251  return nroi;
3252 }
3253 
3254 
3255 /***********************************************************************//**
3256  * @brief Return spatial integral of diffuse source model
3257  *
3258  * @param[in] model Sky Model.
3259  * @param[in] srcEng True photon energy.
3260  * @param[in] srcTime True photon arrival time.
3261  * @param[in] obsEng Observed event energy.
3262  * @param[in] obsTime Observed event arrival time.
3263  * @param[in] obs Observation.
3264  *
3265  * Computes the integral
3266  *
3267  * \f[
3268  * N_{\rm ROI}(E',t'|E,t) = \int_{\rm ROI} P(p',E',t'|E,t) dp'
3269  * \f]
3270  *
3271  * of
3272  *
3273  * \f[
3274  * P(p',E',t'|E,t) = \int
3275  * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp
3276  * \f]
3277  *
3278  * over the Region of Interest (ROI) for a diffuse source model
3279  * \f$S(p,E,t)\f$ and the response function \f$R(p',E',t'|p,E,t)\f$.
3280  ***************************************************************************/
3282  const GEnergy& srcEng,
3283  const GTime& srcTime,
3284  const GEnergy& obsEng,
3285  const GTime& obsTime,
3286  const GObservation& obs) const
3287 {
3288  // Set number of iterations for Romberg integration.
3289  // These values have been determined after careful testing, see
3290  // https://cta-redmine.irap.omp.eu/issues/1248
3291  static const int iter_rho = 9;
3292  static const int iter_phi = 9;
3293 
3294  // Initialise Nroi value
3295  double nroi = 0.0;
3296 
3297  // Retrieve CTA observation, ROI and pointing
3299  const GCTARoi& roi = gammalib::cta_event_list(G_NROI_DIFFUSE, obs).roi();
3300  const GCTAPointing& pnt = cta.pointing();
3301 
3302  // Get pointer on spatial model
3303  const GModelSpatial* spatial =
3304  dynamic_cast<const GModelSpatial*>(model.spatial());
3305  if (spatial == NULL) {
3307  }
3308 
3309  // Get pointing direction zenith angle and azimuth [radians]
3310  double zenith = pnt.zenith();
3311  double azimuth = pnt.azimuth();
3312 
3313  // Get log10(E/TeV) of true photon energy
3314  double srcLogEng = srcEng.log10TeV();
3315 
3316  // Get maximum PSF radius (radians). We do this for the onaxis PSF only,
3317  // as this allows us doing this computation in the outer loop. This
3318  // should be sufficient here, unless the offaxis PSF becomes much worse
3319  // than the onaxis PSF. In this case, we may add a safety factor here
3320  // to make sure we encompass the entire PSF.
3321  double psf_max_radius = psf_delta_max(0.0, 0.0, zenith, azimuth, srcLogEng);
3322 
3323  // Extract ROI radius (radians)
3324  double roi_radius = roi.radius() * gammalib::deg2rad;
3325 
3326  // Compute the ROI radius plus maximum PSF radius (radians). Any photon
3327  // coming from beyond this radius will not make it in the dataspace and
3328  // thus can be neglected.
3329  double roi_psf_radius = roi_radius + psf_max_radius;
3330 
3331  // Perform offset angle integration only if interval is valid
3332  if (roi_psf_radius > 0.0) {
3333 
3334  // Compute rotation matrix to convert from native ROI coordinates,
3335  // given by (theta,phi), into celestial coordinates.
3336  GMatrix ry;
3337  GMatrix rz;
3338  ry.eulery(roi.centre().dir().dec_deg() - 90.0);
3339  rz.eulerz(-roi.centre().dir().ra_deg());
3340  GMatrix rot = (ry * rz).transpose();
3341 
3342  // Setup integration kernel
3343  cta_nroi_diffuse_kern_theta integrand(this,
3344  &cta,
3345  spatial,
3346  &rot,
3347  srcEng,
3348  srcTime,
3349  obsEng,
3350  obsTime,
3351  iter_phi);
3352 
3353  // Integrate over model's zenith angle
3354  GIntegral integral(&integrand);
3355  integral.fixed_iter(iter_rho);
3356  nroi = integral.romberg(0.0, roi_psf_radius);
3357 
3358  // Compile option: Show integration results
3359  #if defined(G_DEBUG_NROI_DIFFUSE)
3360  std::cout << "GCTAResponseIrf::nroi_diffuse:";
3361  std::cout << " roi_psf_radius=" << roi_psf_radius;
3362  std::cout << " nroi=" << nroi;
3363  std::cout << " Etrue=" << srcEng;
3364  std::cout << " Ereco=" << obsEng;
3365  std::cout << " id=" << id << std::endl;
3366  #endif
3367 
3368  } // endif: offset angle range was valid
3369 
3370  // Debug: Check for NaN
3371  #if defined(G_NAN_CHECK)
3372  if (gammalib::is_notanumber(nroi) || gammalib::is_infinite(nroi)) {
3373  std::cout << "*** ERROR: GCTAResponseIrf::nroi_diffuse:";
3374  std::cout << " NaN/Inf encountered";
3375  std::cout << " (nroi=" << nroi;
3376  std::cout << ", roi_psf_radius=" << roi_psf_radius;
3377  std::cout << ")" << std::endl;
3378  }
3379  #endif
3380 
3381  // Return Nroi
3382  return nroi;
3383 }
3384 
3385 
3386 /***********************************************************************//**
3387  * @brief Return spatial integral of composite source model
3388  *
3389  * @param[in] model Sky Model.
3390  * @param[in] srcEng True photon energy.
3391  * @param[in] srcTime True photon arrival time.
3392  * @param[in] obsEng Observed event energy.
3393  * @param[in] obsTime Observed event arrival time.
3394  * @param[in] obs Observation.
3395  *
3396  * Computes the integral
3397  *
3398  * \f[
3399  * N_{\rm ROI}(E',t'|E,t) = \int_{\rm ROI} P(p',E',t'|E,t) dp'
3400  * \f]
3401  *
3402  * of
3403  *
3404  * \f[
3405  * P(p',E',t'|E,t) = \int
3406  * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp
3407  * \f]
3408  *
3409  * over the Region of Interest (ROI) for a composite source model
3410  * \f$S(p,E,t)\f$ and the response function \f$R(p',E',t'|p,E,t)\f$.
3411  ***************************************************************************/
3413  const GEnergy& srcEng,
3414  const GTime& srcTime,
3415  const GEnergy& obsEng,
3416  const GTime& obsTime,
3417  const GObservation& obs) const
3418 {
3419  // Initialise nroi
3420  double nroi = 0.0;
3421 
3422  // Get composite model
3423  GModelSpatialComposite* comp =
3424  dynamic_cast<GModelSpatialComposite*>(model.spatial());
3425 
3426  // Loop over model components
3427  for (int i = 0; i < comp->components(); ++i) {
3428 
3429  // Create new sky model
3430  GModelSky sky = GModelSky(model);
3431 
3432  // Assign spatial part
3433  sky.spatial(comp->component(i));
3434 
3435  // Compute nroi
3436  nroi += this->nroi(sky, srcEng, srcTime, obsEng, obsTime, obs) *
3437  comp->scale(i);
3438 
3439  }
3440 
3441  // Divide by number of model components
3442  double sum = comp->sum_of_scales();
3443  if (sum > 0.0) {
3444  nroi /= sum;
3445  }
3446 
3447  // Return nroi
3448  return nroi;
3449 
3450 }
CTA background model base class definition.
void load_aeff(const GFilename &filename)
Load effective area.
GSkyDir dir(void) const
Return position of elliptical spatial model.
GCTAResponseIrf(void)
Void constructor.
void load_background(const GFilename &filename)
Load background model.
CTA response helper classes definition.
CTA 2D energy dispersion class.
Definition: GCTAEdisp2D.hpp:90
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.
#define G_NIRF
void remove_thetacut(const GCTAResponseIrf &rsp)
Remove thetacut.
CTA event list class interface definition.
virtual double mc(GRan &ran, const double &logE, const double &theta=0.0, const double &phi=0.0, const double &zenith=0.0, const double &azimuth=0.0, const bool &etrue=true) const =0
GFitsTable * table(const int &extno)
Get pointer to table HDU.
Definition: GFits.cpp:472
virtual GCTAEdisp * clone(void) const =0
Clones object.
CTA performance table effective area class.
bool contains(const int &extno) const
Check if HDU exists in FITS file.
Definition: GFits.hpp:282
Kernel for IRF offest angle integration of the diffuse source model.
bool has_card(const int &cardno) const
Check existence of header card.
Definition: GFitsHDU.hpp:233
double m_hi_safe_thres
Safe high energy threshold.
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:381
Point source spatial model class interface definition.
virtual bool apply_edisp(void) const
Signal if energy dispersion should be applied.
Radial shell model class interface definition.
const GCTAInstDir & cta_dir(const std::string &origin, const GEvent &event)
Retrieve CTA instrument direction from generic event.
CTA event atom class definition.
virtual void write(GXmlElement &xml) const
Write response information into XML element.
Abstract elliptical spatial model base class.
Random number generator class definition.
void warning(const std::string &origin, const std::string &message)
Emits warning.
Definition: GTools.cpp:1276
double irf_elliptical(const GEvent &event, const GSource &source, const GObservation &obs) const
Return Irf value for elliptical source model.
virtual std::string instrument(void) const =0
GModelSpectral * spectral(void) const
Return spectral model component.
Definition: GModelSky.hpp:266
virtual std::string print(const GChatter &chatter=NORMAL) const =0
Print content of object.
virtual bool use_edisp(void) const
Signal if response uses energy dispersion.
virtual void roi(const GRoi &roi)
Set ROI.
CTA vector point spread function class.
bool has_free_pars(void) const
Checks if the spatial model has free parameters.
const GModelSpatial * component(const int &index) const
Returns pointer to spatial component element.
int size(void) const
Return number of energy boundaries.
Definition: GEbounds.hpp:162
virtual GCTAResponseIrf & operator=(const GCTAResponseIrf &rsp)
Assignment operator.
const std::string & mission(void) const
Return mission.
Definition: GCaldb.hpp:137
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1107
GCTAAeff * m_aeff
Effective area.
virtual double eval(const GEnergy &srcEng, const GTime &srcTime=GTime(), const bool &gradients=false) const =0
#define G_LOAD_EDISP
GModelTemporal * temporal(void) const
Return temporal model component.
Definition: GModelSky.hpp:282
Abstract radial spatial model base class interface definition.
#define G_LOAD_BACKGROUND
const GCTABackground * background(void) const
Return pointer to background model.
CTA performance table background class definition.
void free_members(void)
Delete class members.
bool is_empty(void) const
Signal if filename is empty.
Definition: GFilename.hpp:160
CTA instrument response function class definition.
const GCTAPointing & cta_pnt(const std::string &origin, const GObservation &obs)
Retrieve CTA pointing from generic observation.
bool contains(const std::string &colname) const
Checks the presence of a column in table.
Definition: GFitsTable.cpp:716
Abstract interface for the event classes.
Definition: GEvent.hpp:71
GFilename filename(const std::string &detector, const std::string &filter, const std::string &codename, const std::string &date, const std::string &time, const std::string &expr)
Return calibration file name based on selection parameters.
Definition: GCaldb.cpp:507
Kernel for zenith angle Nroi integration of elliptical model.
const double & radius(void) const
Returns radius of region of interest in degrees.
Definition: GCTARoi.hpp:124
XML element node class.
Definition: GXmlElement.hpp:47
CTA pointing class interface definition.
bool has_extname(void) const
Signal if filename has an extension name.
Definition: GFilename.hpp:255
Kernel for zenith angle Nroi integration of radial model.
CTA point spread function class with a King profile.
Definition: GCTAPsfKing.hpp:54
CTA Redistribution Matrix File (RMF) energy dispersion class.
const double & zenith(void) const
Return pointing zenith angle.
const std::string & rspname(void) const
Return response name.
double TeV(void) const
Return energy in TeV.
Definition: GEnergy.cpp:348
double sum(const GVector &vector)
Computes vector sum.
Definition: GVector.cpp:908
double acos(const double &arg)
Computes acos by avoiding NaN due to rounding errors.
Definition: GMath.cpp:69
std::string extname(const std::string &defaultname="") const
Return extension name.
Definition: GFilename.cpp:385
CTA ARF effective area class definition.
double log10TeV(void) const
Return log10 of energy in TeV.
Definition: GEnergy.cpp:457
Random number generator class.
Definition: GRan.hpp:44
void init_members(void)
Initialise class members.
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
Interface for the CTA region of interest class.
Definition: GCTARoi.hpp:48
virtual GEnergy mc(GRan &ran, const GEnergy &etrue, const double &theta=0.0, const double &phi=0.0, const double &zenith=0.0, const double &azimuth=0.0) const =0
Time class.
Definition: GTime.hpp:54
Gammalib tools definition.
GIntegral class interface definition.
Definition: GIntegral.hpp:46
FITS file class.
Definition: GFits.hpp:63
FITS file class interface definition.
CTA 3D background class.
void thetacut(const double &thetacut)
Set theta cut angle.
const GCaldb & caldb(void) const
Return calibration database.
#define G_IRF_ELLIPTICAL
void pointing(const GCTAPointing &pointing)
Set CTA pointing.
void init_members(void)
Initialise class members.
Definition: GResponse.cpp:754
std::string strip_whitespace(const std::string &arg)
Strip leading and trailing whitespace from string.
Definition: GTools.cpp:75
Source class definition.
virtual double delta_max(const double &logE, const double &theta=0.0, const double &phi=0.0, const double &zenith=0.0, const double &azimuth=0.0, const bool &etrue=true) const =0
virtual GClassCode code(void) const =0
GCTABackground * m_background
Background.
void dir(const GSkyDir &dir)
Set sky direction.
void init_members(void)
Initialise class members.
#define G_NROI_RADIAL
CTA 2D point spread function class.
Definition: GCTAPsf2D.hpp:54
void eulery(const double &angle)
Set Euler rotation matrix around y axis.
Definition: GMatrix.cpp:1093
virtual double theta_max(void) const =0
#define G_PSF_DELTA_MAX
void scale(const double &scale)
Set effective area scaling factor.
double psf_delta_max(const double &theta, const double &phi, const double &zenith, const double &azimuth, const double &srcLogEng) const
Return maximum angular separation (in radians)
std::string centre(const std::string &s, const int &n, const char &c= ' ')
Centre string to achieve a length of n characters.
Definition: GTools.cpp:1008
void id(const std::string &id)
Set observation identifier.
const std::string & instrument(void) const
Return instrument.
Definition: GCaldb.hpp:149
GCTAEdisp * m_edisp
Energy dispersion.
CTA performance table effective area class definition.
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:184
double sum_of_scales(void) const
Returns sum of all model scales.
Class that handles photons.
Definition: GPhoton.hpp:47
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:167
virtual GCTAResponse & operator=(const GCTAResponse &rsp)
Assignment operator.
Definition of support function used by CTA classes.
bool m_apply_edisp
Apply energy dispersion.
double real(const std::string &keyname) const
Return card value as double precision.
Definition: GFitsHDU.hpp:423
const GCTAPsf * psf(void) const
Return pointer to point spread function.
double radius(void) const
Return shell radius.
CTA effective area base class definition.
double semiminor(void) const
Return semi-minor axis of ellipse.
const double pihalf
Definition: GMath.hpp:38
const double deg2rad
Definition: GMath.hpp:43
const std::string extname_arf
Definition: GArf.hpp:45
double m_lo_safe_thres
Safe low energy threshold.
GFilename irf_filename(const std::string &filename) const
Return filename with appropriate extension.
Calibration database class.
Definition: GCaldb.hpp:66
CTA 2D effective area class.
Definition: GCTAAeff2D.hpp:55
Energy boundaries container class.
Definition: GEbounds.hpp:60
const GXmlAttribute * attribute(const int &index) const
Return attribute.
virtual double irf(const GEvent &event, const GPhoton &photon, const GObservation &obs) const
Return value of instrument response function.
double ra(void) const
Return Right Ascencion of model centre.
CTA performance table energy dispersion class.
#define G_AEFF
virtual GEbounds ebounds(const GEnergy &obsEnergy) const
Return true energy boundaries for a specific observed energy.
const std::string & name(void) const
Return parameter name.
Definition: GModel.hpp:257
const GCTAAeff * aeff(void) const
Return pointer to effective area response.
Filename class.
Definition: GFilename.hpp:62
void rotate_deg(const double &phi, const double &theta)
Rotate sky direction by zenith and azimuth angle.
Definition: GSkyDir.cpp:389
const GEnergy & emin(void) const
Return minimum energy of all intervals.
Definition: GEbounds.hpp:186
double dec(void) const
Return Declination of model centre.
const GCTAInstDir & centre(void) const
Returns region of interest centre.
Definition: GCTARoi.hpp:110
XSPEC Auxiliary Response File class definition.
#define G_IRF
double trapzd(const double &a, const double &b, const int &n=1, double result=0.0)
Perform Trapezoidal integration.
Definition: GIntegral.cpp:585
void eulerz(const double &angle)
Set Euler rotation matrix around z axis.
Definition: GMatrix.cpp:1130
CTA 2D effective area class definition.
virtual std::string print(const GChatter &chatter=NORMAL) const =0
Print content of object.
bool has_attribute(const std::string &name) const
Check if element has a given attribute.
CTA performance table energy dispersion class definition.
GCTAInstDir instdir(const GSkyDir &skydir) const
Get instrument direction from sky direction.
GXmlElement * xml_need_par(const std::string &origin, GXmlElement &xml, const std::string &name)
Return pointer to parameter with given name in XML element.
Definition: GTools.cpp:1527
double dec(void) const
Return Declination of model centre.
bool exists(void) const
Checks whether file exists.
Definition: GFilename.cpp:223
Kernel for Nroi offest angle integration of diffuse model.
virtual std::string print(const GChatter &chatter=NORMAL) const =0
Print content of object.
double nirf(const GPhoton &photon, const GEnergy &obsEng, const GTime &obsTime, const GObservation &obs) const
Return spatial integral of Instrument Response Function.
Radial ring model class interface definition.
double uniform(void)
Returns random double precision floating value in range 0 to 1.
Definition: GRan.cpp:242
CTA pointing class.
void free_members(void)
Delete class members.
CTA RMF energy dispersion class definition.
virtual double deadc(const GTime &time=GTime()) const =0
const GCTAEventList & cta_event_list(const std::string &origin, const GObservation &obs)
Retrieve CTA event list from generic observation.
bool m_use_nroi_cache
Control usage of nroi cache.
Definition: GResponse.hpp:313
Abstract interface for FITS table.
Definition: GFitsTable.hpp:44
void load_psf(const GFilename &filename)
Load CTA PSF vector.
virtual GCTAResponseIrf * clone(void) const
Clone instance.
const GTime & time(void) const
Return photon time.
Definition: GPhoton.hpp:134
const GTime & time(void) const
Return photon arrival time.
Definition: GSource.hpp:156
double resolution(const GModelSpatial *model)
Determine resolution of spatial model.
GChatter
Definition: GTypemaps.hpp:33
Abstract CTA energy dispersion base class definition.
virtual GCTAAeff * clone(void) const =0
Clones object.
const std::string extname_cta_aeff2d
Definition: GCTAAeff2D.hpp:43
CTA 3D background class definition.
GResponseCache m_nroi_cache
Definition: GResponse.hpp:324
CTA 2D point spread function class definition.
bool is_fits(void) const
Checks whether file is a FITS file.
Definition: GFilename.cpp:313
virtual double eval(const GTime &srcTime, const bool &gradients=false) const =0
double ra(void) const
Return Right Ascencion of model centre.
const GCTAObservation & cta_obs(const std::string &origin, const GObservation &obs)
Retrieve CTA observation from generic observation.
Spatial composite model class interface definition.
CTA 2D energy dispersion class definition.
const GModelSpatial * model(void) const
Return spatial model component.
Definition: GSource.hpp:128
Abstract observation base class.
const GEnergy & energy(void) const
Return photon energy.
Definition: GSource.hpp:142
Spatial composite model.
King profile CTA point spread function class definition.
virtual GCTABackground * clone(void) const =0
Clones object.
GCaldb m_caldb
Calibration database.
double npsf(const GSkyDir &srcDir, const double &srcLogEng, const GTime &srcTime, const GCTAPointing &pnt, const GCTARoi &roi) const
Return result of PSF integration over ROI.
virtual GFilename filename(void) const =0
bool has_scales(void) const
Signals that model has scales.
Definition: GModel.hpp:350
CTA region of interest class interface definition.
CTA observation class interface definition.
CTA exception handler interface definition.
virtual std::string print(const GChatter &chatter=NORMAL) const =0
Print content of object.
double offset_sigma(void) const
Return offset angle dependence (degrees)
GModelPar & scale(const int &index)
Returns reference to scale parameter by index.
Definition: GModel.cpp:365
CTA point spread function base class definition.
void load_edisp(const GFilename &filename)
Load energy dispersion information.
void load(const std::string &rspname)
Load CTA response.
virtual GEbounds etrue_bounds(const GEnergy &ereco, const double &theta=0.0, const double &phi=0.0, const double &zenith=0.0, const double &azimuth=0.0) const =0
Calibration database class interface definition.
std::string gadf_hduclas4(const GFits &fits, const std::string &hduclas4)
Return extension name for GADF response table of given HDU class 4.
double nroi_diffuse(const GModelSky &model, const GEnergy &srcEng, const GTime &srcTime, const GEnergy &obsEng, const GTime &obsTime, const GObservation &obs) const
Return spatial integral of diffuse source model.
Point source spatial model.
CTA performance table point spread function class definition.
bool xml_has_par(const GXmlElement &xml, const std::string &name)
Checks if parameter with given name in XML element exists.
Definition: GTools.cpp:1486
CTA instrument response function class.
virtual ~GCTAResponseIrf(void)
Destructor.
CTA instrument response function class.
void free_members(void)
Delete class members.
Definition: GResponse.cpp:802
std::string m_xml_rspname
Response name in XML file.
Abstract elliptical spatial model base class interface definition.
const double & azimuth(void) const
Return pointing azimuth angle.
double posang(const GSkyDir &dir, const std::string &coordsys="CEL") const
Compute position angle between sky directions in radians.
Definition: GSkyDir.cpp:1083
#define G_IRF_DIFFUSE
int components(void) const
Return number of model components.
const std::string extname_cta_edisp2d
Definition: GCTAEdisp2D.hpp:46
#define G_MC
int iter_phi(const double &rho, const double &resolution, const int &iter_min, const int &iter_max)
Determine number of azimuthal Romberg iterations.
CTA point spread function vector class definition.
GSkyDir dir(void) const
Return position of point source.
Kernel for elliptical model zenith angle integration of IRF.
std::string m_rspname
Name of the instrument response.
double irf_radial(const GEvent &event, const GSource &source, const GObservation &obs) const
Return IRF value for radial source model.
Sky model class interface definition.
Sky model class.
Definition: GModelSky.hpp:121
GSkyDir dir(void) const
Return position of radial spatial model.
virtual void clear(void)
Clear instance.
double nroi_radial(const GModelSky &model, const GEnergy &srcEng, const GTime &srcTime, const GEnergy &obsEng, const GTime &obsTime, const GObservation &obs) const
Return spatial integral of radial source model.
double posangle(void) const
Return Position Angle of model.
#define G_NROI_ELLIPTICAL
double nroi_composite(const GModelSky &model, const GEnergy &srcEng, const GTime &srcTime, const GEnergy &obsEng, const GTime &obsTime, const GObservation &obs) const
Return spatial integral of composite source model.
std::string m_xml_caldb
Calibration database string in XML file.
CTA event atom class.
GChatter reduce(const GChatter &chatter)
Reduce chattiness by one level.
Definition: GTypemaps.hpp:65
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1233
#define G_IRF_RADIAL
const GEnergy & energy(void) const
Return photon energy.
Definition: GPhoton.hpp:122
std::string print(const GChatter &chatter=NORMAL) const
Print calibration database information.
Definition: GCaldb.cpp:640
virtual GCTAPsf * clone(void) const =0
Clones object.
CTA point spread function table class definition.
virtual double theta_max(void) const =0
double dist(const GSkyDir &dir) const
Compute angular distance between sky directions in radians.
Definition: GSkyDir.hpp:269
void fixed_iter(const int &iter)
Set fixed number of iterations.
Definition: GIntegral.hpp:183
Kernel for radial model zenith angle integration of IRF.
virtual void read(const GXmlElement &xml)
Read response from XML element.
double irf_ptsrc(const GEvent &event, const GSource &source, const GObservation &obs) const
Return value of point source instrument response function.
CTA performance table background class.
GModelSpatial * spatial(void) const
Return spatial model component.
Definition: GModelSky.hpp:250
#define G_NROI_DIFFUSE
CTA instrument direction class.
Definition: GCTAInstDir.hpp:63
const GSkyDir & dir(void) const
Return pointing sky direction.
Abstract spatial model base class.
GResponseCache m_irf_cache
Definition: GResponse.hpp:323
Abstract radial spatial model base class.
double radius(void) const
Return ring inner radius.
Generic matrix class definition.
Definition: GMatrix.hpp:79
CTA point spread function table class.
const GEnergy & emax(void) const
Return maximum energy of all intervals.
Definition: GEbounds.hpp:198
std::string filepath(const std::string &pathname, const std::string &filename)
Build file path from path name and file name.
Definition: GTools.cpp:388
CTA performance table point spread function class.
void copy_members(const GCTAResponseIrf &rsp)
Copy class members.
const double rad2deg
Definition: GMath.hpp:44
bool contains(const std::string &name, const GEnergy &ereco, const GEnergy &etrue, double *value=NULL) const
Check if cache contains a value for specific parameters.
GCTAEventAtom * mc(const double &area, const GPhoton &photon, const GObservation &obs, GRan &ran) const
Simulate event from photon.
CTA observation class.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1033
Integration kernel for npsf() method.
double nroi_elliptical(const GModelSky &model, const GEnergy &srcEng, const GTime &srcTime, const GEnergy &obsEng, const GTime &obsTime, const GObservation &obs) const
Return spatial integral of elliptical source model.
const GCTAEdisp * edisp(void) const
Return pointer to energy dispersion.
#define G_READ
Class that handles gamma-ray sources.
Definition: GSource.hpp:53
#define G_LOAD_PSF
void set(const std::string &name, const GEnergy &ereco, const GEnergy &etrue, const double &value)
Set cache value.
Integration class interface definition.
int iter_rho(const double &rho_max, const double &resolution, const int &iter_min, const int &iter_max)
Determine number of radial Romberg iterations.
Sky direction class.
Definition: GSkyDir.hpp:62
#define G_PSF
void close(void)
Close FITS file.
Definition: GFits.cpp:1314
double semimajor(void) const
Return semi-major axis of ellipse.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print CTA response information.
const GSkyDir & dir(void) const
Return photon sky direction.
Definition: GPhoton.hpp:110
double irf_diffuse(const GEvent &event, const GSource &source, const GObservation &obs) const
Return value of diffuse source instrument response function.
const GXmlElement * xml_get_par(const std::string &origin, const GXmlElement &xml, const std::string &name)
Return pointer to parameter with given name in XML element.
Definition: GTools.cpp:1579
void clear(void)
Clear calibration database.
Definition: GCaldb.cpp:194
double scale(const int &index) const
Returns scale of spatial component.
std::string rootdir(void) const
Return path to CALDB root directory.
Definition: GCaldb.cpp:259
const GCTAInstDir & dir(void) const
Return instrument direction.
Filename class interface definition.
#define G_WRITE
std::string print(const GChatter &chatter=NORMAL) const
Print response cache.
double nroi_ptsrc(const GModelSky &model, const GEnergy &srcEng, const GTime &srcTime, const GEnergy &obsEng, const GTime &obsTime, const GObservation &obs) const
Return spatial integral of point source model.
void sigma(const double &sigma)
Set sigma for offset angle dependence.
Mathematical function definitions.
void sigma(const double &sigma)
Set sigma for offset angle dependence.
CTA ARF effective area class.
Definition: GCTAAeffArf.hpp:51
double todouble(const std::string &arg)
Convert string into double precision value.
Definition: GTools.cpp:816
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
GCTAPsf * m_psf
Point spread function.
#define G_LOAD_AEFF
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:415
void sigma(const double &sigma)
Set sigma for offset angle dependence.
const std::string extname_rmf
Definition: GRmf.hpp:44