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