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