GammaLib  1.7.0.dev
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GCTAResponseCube.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GCTAResponseCube.cpp - CTA cube analysis response function class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2014-2018 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 GCTAResponseCube.cpp
23  * @brief CTA cube-style response function class implementation
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <cmath>
32 #include <string>
33 #include "GTools.hpp"
34 #include "GCTAResponseCube.hpp"
35 #include "GCTAResponse_helpers.hpp"
36 #include "GCTAInstDir.hpp"
37 #include "GCTAEventBin.hpp"
38 #include "GCTASupport.hpp"
40 #include "GModelSpatialRadial.hpp"
44 #include "GModelSpatialDiffuse.hpp"
45 #include "GPhoton.hpp"
46 #include "GSource.hpp"
47 #include "GEvent.hpp"
48 #include "GSkyDir.hpp"
49 #include "GEnergy.hpp"
50 #include "GTime.hpp"
51 #include "GIntegral.hpp"
52 #include "GObservation.hpp"
53 
54 /* __ Method name definitions ____________________________________________ */
55 #define G_IRF "GCTAResponseCube::irf(GEvent&, GPhoton& GObservation&)"
56 #define G_IRF_PTSRC "GCTAResponseCube::irf_ptsrc(GEvent&, GSource&,"\
57  " GObservation&)"
58 #define G_IRF_RADIAL "GCTAResponseCube::irf_radial(GEvent&, GSource&,"\
59  " GObservation&)"
60 #define G_IRF_DIFFUSE "GCTAResponseCube::irf_diffuse(GEvent&, GSource&,"\
61  " GObservation&)"
62 #define G_NROI "GCTAResponseCube::nroi(GModelSky&, GEnergy&, GTime&, "\
63  "GObservation&)"
64 #define G_EBOUNDS "GCTAResponseCube::ebounds(GEnergy&)"
65 #define G_READ "GCTAResponseCube::read(GXmlElement&)"
66 #define G_WRITE "GCTAResponseCube::write(GXmlElement&)"
67 
68 /* __ Macros _____________________________________________________________ */
69 
70 /* __ Coding definitions _________________________________________________ */
71 #define G_RADIAL_PSF_BASED //!< Use Psf-based integration for radial model
72 
73 /* __ Debug definitions __________________________________________________ */
74 
75 /* __ Constants __________________________________________________________ */
76 
77 
78 /*==========================================================================
79  = =
80  = Constructors/destructors =
81  = =
82  ==========================================================================*/
83 
84 /***********************************************************************//**
85  * @brief Void constructor
86  *
87  * Constructs void CTA response.
88  ***************************************************************************/
90 {
91  // Initialise members
92  init_members();
93 
94  // Return
95  return;
96 }
97 
98 
99 /***********************************************************************//**
100  * @brief Copy constructor
101  *
102  * @param[in] rsp CTA response.
103  *
104  * Constructs CTA cube response by making a deep copy of an existing
105  * object.
106  **************************************************************************/
108  GCTAResponse(rsp)
109 {
110  // Initialise members
111  init_members();
112 
113  // Copy members
114  copy_members(rsp);
115 
116  // Return
117  return;
118 }
119 
120 
121 /***********************************************************************//**
122  * @brief XML constructor
123  *
124  * @param[in] xml XML element.
125  *
126  * Construct CTA response from XML element.
127  ***************************************************************************/
129 {
130  // Initialise members
131  init_members();
132 
133  // Read information from XML element
134  read(xml);
135 
136  // Return
137  return;
138 }
139 
140 
141 /***********************************************************************//**
142  * @brief Response constructor
143  *
144  * @param[in] exposure CTA cube analysis exposure.
145  * @param[in] psf CTA cube analysis point spread function.
146  * @param[in] background CTA cube background response.
147  *
148  * Constructs CTA cube analysis response from a cube analysis exposure,
149  * a point spread function cube and a background cube.
150  **************************************************************************/
152  const GCTACubePsf& psf,
153  const GCTACubeBackground& background) :
154  GCTAResponse()
155 {
156  // Initialise members
157  init_members();
158 
159  // Set members
161  m_psf = psf;
163 
164  // Signal that no energy dispersion was given
165  m_has_edisp = false;
166 
167  // Return
168  return;
169 }
170 
171 
172 /***********************************************************************//**
173  * @brief Response constructor
174  *
175  * @param[in] exposure CTA cube analysis exposure.
176  * @param[in] psf CTA cube analysis point spread function.
177  * @param[in] edisp CTA cube energy dispersion response.
178  * @param[in] background CTA cube background response.
179  *
180  * Constructs CTA cube analysis response from a cube analysis exposure,
181  * a point spread function cube, an energy dispersion cube and a background cube.
182  **************************************************************************/
184  const GCTACubePsf& psf,
185  const GCTACubeEdisp& edisp,
186  const GCTACubeBackground& background) :
187  GCTAResponse()
188 {
189  // Initialise members
190  init_members();
191 
192  // Set members
194  m_psf = psf;
195  m_edisp = edisp;
197 
198  // Signal that edisp is available
199  m_has_edisp = true;
200 
201  // Return
202  return;
203 }
204 
205 
206 /***********************************************************************//**
207  * @brief Destructor
208  *
209  * Destroys instance of CTA response object.
210  ***************************************************************************/
212 {
213  // Free members
214  free_members();
215 
216  // Return
217  return;
218 }
219 
220 
221 /*==========================================================================
222  = =
223  = Operators =
224  = =
225  ==========================================================================*/
226 
227 /***********************************************************************//**
228  * @brief Assignment operator
229  *
230  * @param[in] rsp CTA response.
231  * @return CTA response.
232  *
233  * Assigns CTA response object to another CTA response object. The assignment
234  * performs a deep copy of all information, hence the original object from
235  * which the assignment has been performed can be destroyed after this
236  * operation without any loss of information.
237  ***************************************************************************/
239 {
240  // Execute only if object is not identical
241  if (this != &rsp) {
242 
243  // Copy base class members
244  this->GCTAResponse::operator=(rsp);
245 
246  // Free members
247  free_members();
248 
249  // Initialise members
250  init_members();
251 
252  // Copy members
253  copy_members(rsp);
254 
255  } // endif: object was not identical
256 
257  // Return this object
258  return *this;
259 }
260 
261 
262 /*==========================================================================
263  = =
264  = Public methods =
265  = =
266  ==========================================================================*/
267 
268 /***********************************************************************//**
269  * @brief Clear instance
270  *
271  * Clears CTA response object by resetting all members to an initial state.
272  * Any information that was present in the object before will be lost.
273  ***************************************************************************/
275 {
276  // Free class members (base and derived classes, derived class first)
277  free_members();
279  this->GResponse::free_members();
280 
281  // Initialise members
282  this->GResponse::init_members();
284  init_members();
285 
286  // Return
287  return;
288 }
289 
290 
291 /***********************************************************************//**
292  * @brief Clone instance
293  *
294  * @return Pointer to deep copy of CTA response.
295  *
296  * Creates a clone (deep copy) of a CTA response object.
297  ***************************************************************************/
299 {
300  return new GCTAResponseCube(*this);
301 }
302 
303 
304 /***********************************************************************//**
305  * @brief Return instrument response
306  *
307  * @param[in] event Observed event.
308  * @param[in] photon Incident photon.
309  * @param[in] obs Observation (not used).
310  * @return Instrument response.
311  ***************************************************************************/
312 double GCTAResponseCube::irf(const GEvent& event,
313  const GPhoton& photon,
314  const GObservation& obs) const
315 {
316  // Retrieve event instrument direction
317  const GCTAInstDir& dir = gammalib::cta_dir(G_IRF, event);
318 
319  // Get event attributes
320  const GSkyDir& obsDir = dir.dir();
321  const GEnergy& obsEng = event.energy();
322 
323  // Get photon attributes
324  const GSkyDir& srcDir = photon.dir();
325  const GEnergy& srcEng = photon.energy();
326  //const GTime& srcTime = photon.time();
327 
328  // Determine angular separation between true and measured photon
329  // direction in radians
330  double delta = obsDir.dist(srcDir);
331 
332  // Get maximum angular separation for PSF (in radians)
333  double delta_max = psf().delta_max();
334 
335  // Initialise IRF value
336  double irf = 0.0;
337 
338  // Get livetime (in seconds)
339  double livetime = exposure().livetime();
340 
341  // Continue only if livetime is >0 and if we're sufficiently close
342  // to the PSF
343  if ((livetime > 0.0) && (delta <= delta_max)) {
344 
345  // Get exposure
346  irf = exposure()(srcDir, srcEng);
347 
348  // Multiply-in PSF
349  if (irf > 0.0) {
350 
351  // Get PSF component
352  irf *= psf()(srcDir, delta, srcEng);
353 
354  // Multiply-in energy dispersion
355  if (use_edisp() && irf > 0.0) {
356  irf *= edisp()(obsEng, srcEng, srcDir);
357  }
358 
359  // Divide by livetime
360  irf /= livetime;
361 
362  // Apply deadtime correction
363  irf *= exposure().deadc();
364 
365  } // endif: Aeff was non-zero
366 
367  } // endif: we were sufficiently close to PSF and livetime was >0
368 
369  // Compile option: Check for NaN/Inf
370  #if defined(G_NAN_CHECK)
372  std::cout << "*** ERROR: GCTAResponseCube::irf:";
373  std::cout << " NaN/Inf encountered";
374  std::cout << " irf=" << irf;
375  std::cout << std::endl;
376  }
377  #endif
378 
379  // Return IRF value
380  return irf;
381 }
382 
383 
384 /***********************************************************************//**
385  * @brief Return integral of event probability for a given sky model over ROI
386  *
387  * @param[in] model Sky model.
388  * @param[in] obsEng Observed photon energy.
389  * @param[in] obsTime Observed photon arrival time.
390  * @param[in] obs Observation.
391  * @return 0.0
392  *
393  * @exception GException::feature_not_implemented
394  * Method is not implemented.
395  *
396  * Computes the integral
397  *
398  * \f[
399  * N_{\rm ROI}(E',t') = \int_{\rm ROI} P(p',E',t') dp'
400  * \f]
401  *
402  * of the event probability
403  *
404  * \f[
405  * P(p',E',t') = \int \int \int
406  * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp \, dE \, dt
407  * \f]
408  *
409  * for a given sky model \f$S(p,E,t)\f$ and response function
410  * \f$R(p',E',t'|p,E,t)\f$ over the Region of Interest (ROI).
411  *
412  * @todo Implement method (is maybe not really needed)
413  ***************************************************************************/
414 double GCTAResponseCube::nroi(const GModelSky& model,
415  const GEnergy& obsEng,
416  const GTime& obsTime,
417  const GObservation& obs) const
418 {
419  // Method is not implemented
420  std::string msg = "Spatial integration of sky model over the data space "
421  "is not implemented.";
423 
424  // Return Npred
425  return (0.0);
426 }
427 
428 
429 /***********************************************************************//**
430  * @brief Return true energy boundaries for a specific observed energy
431  *
432  * @param[in] obsEng Observed photon energy.
433  * @return Boundaries in true energy.
434  ***************************************************************************/
436 {
437  // Get energy boundaries from energy dispersion
438  GEbounds ebounds = m_edisp.ebounds(obsEng);
439 
440  // Return energy boundaries
441  return ebounds;
442 }
443 
444 
445 /***********************************************************************//**
446  * @brief Read response information from XML element
447  *
448  * @param[in] xml XML element.
449  *
450  * Reads response information from an XML element. The Exposure, Psf,
451  * background cubes, and optionally the energy dispersion cube, are specified
452  * using
453  *
454  * <observation name="..." id="..." instrument="...">
455  * ...
456  * <parameter name="ExposureCube" file="..."/>
457  * <parameter name="PsfCube" file="..."/>
458  * <parameter name="EdispCube" file="..."/>
459  * <parameter name="BkgCube" file="..."/>
460  * </observation>
461  *
462  ***************************************************************************/
464 {
465  // Clear response
466  clear();
467 
468  // Get exposure cube information and load cube
469  const GXmlElement* exppar = gammalib::xml_get_par(G_READ, xml, "ExposureCube");
470  std::string expname = gammalib::xml_file_expand(xml,
471  exppar->attribute("file"));
472 
473  // Get PSF cube information and load cube
474  const GXmlElement* psfpar = gammalib::xml_get_par(G_READ, xml, "PsfCube");
475  std::string psfname = gammalib::xml_file_expand(xml,
476  psfpar->attribute("file"));
477 
478  // Get background cube information and load cube
479  const GXmlElement* bkgpar = gammalib::xml_get_par(G_READ, xml, "BkgCube");
480  std::string bkgname = gammalib::xml_file_expand(xml,
481  bkgpar->attribute("file"));
482 
483  // Load cubes
484  m_exposure.load(expname);
485  m_psf.load(psfname);
486  m_background.load(bkgname);
487 
488  // Optionally load energy dispersion cube
489  if (gammalib::xml_has_par(xml, "EdispCube")) {
490  const GXmlElement* edisppar = gammalib::xml_get_par(G_READ, xml, "EdispCube");
491  std::string edispname = gammalib::xml_file_expand(xml,
492  edisppar->attribute("file"));
493  m_edisp.load(edispname);
494  m_has_edisp = true;
495  }
496 
497  // Return
498  return;
499 }
500 
501 
502 /***********************************************************************//**
503  * @brief Write response information into XML element
504  *
505  * @param[in] xml XML element.
506  *
507  * Writes response information into an XML element. The Exposure, Psf
508  * and background cubes, and optionally the energy dispersion cube, are
509  * specified using
510  *
511  * <observation name="..." id="..." instrument="...">
512  * ...
513  * <parameter name="ExposureCube" file="..."/>
514  * <parameter name="PsfCube" file="..."/>
515  * <parameter name="EdispCube" file="..."/>
516  * <parameter name="BkgCube" file="..."/>
517  * </observation>
518  *
519  ***************************************************************************/
521 {
522  // Add exposure cube filename
523  std::string filename = gammalib::xml_file_reduce(xml, m_exposure.filename());
524  if (!(filename.empty())) {
525  GXmlElement* par = gammalib::xml_need_par(G_WRITE, xml, "ExposureCube");
526  par->attribute("file", filename);
527  }
528 
529  // Add PSF cube filename
530  filename = gammalib::xml_file_reduce(xml, m_psf.filename());
531  if (!(filename.empty())) {
532  GXmlElement* par = gammalib::xml_need_par(G_WRITE, xml, "PsfCube");
533  par->attribute("file", filename);
534  }
535 
536  // Optionally add energy dispersions cube filename
537  if (m_has_edisp) {
538  filename = gammalib::xml_file_reduce(xml, m_edisp.filename());
539  if (!(filename.empty())) {
540  GXmlElement* par = gammalib::xml_need_par(G_WRITE, xml, "EdispCube");
541  par->attribute("file", filename);
542  }
543  }
544 
545  // Add background cube filename
547  if (!(filename.empty())) {
548  GXmlElement* par = gammalib::xml_need_par(G_WRITE, xml, "BkgCube");
549  par->attribute("file", filename);
550  }
551 
552  // Return
553  return;
554 }
555 
556 
557 /***********************************************************************//**
558  * @brief Print response information
559  *
560  * @param[in] chatter Chattiness.
561  * @return String containing response information.
562  ***************************************************************************/
563 std::string GCTAResponseCube::print(const GChatter& chatter) const
564 {
565  // Initialise result string
566  std::string result;
567 
568  // Continue only if chatter is not silent
569  if (chatter != SILENT) {
570 
571  // Append header
572  result.append("=== GCTAResponseCube ===");
573 
574  // Append response information
575  result.append("\n"+gammalib::parformat("Energy dispersion"));
576  if (use_edisp()) {
577  result.append("Used");
578  }
579  else {
580  if (apply_edisp()) {
581  result.append("Not available");
582  }
583  else {
584  result.append("Not used");
585  }
586  }
587 
588  // Get reduced chatter level
589  GChatter reduced_chatter = gammalib::reduce(chatter);
590 
591  // Append detailed information
592  if (chatter >= NORMAL) {
593 
594  // Append exposure cube information
595  result.append("\n"+m_exposure.print(reduced_chatter));
596 
597  // Append point spread function information
598  result.append("\n"+m_psf.print(reduced_chatter));
599 
600  // Optionally append energy dispersion information
601  if (m_has_edisp) {
602  result.append("\n"+m_edisp.print(reduced_chatter));
603  }
604 
605  // Append background information
606  result.append("\n"+m_background.print(reduced_chatter));
607 
608  // Append cache information
609  result.append("\n"+m_irf_cache.print(reduced_chatter));
610 
611  } // endif: appended detailed information
612 
613  } // endif: chatter was not silent
614 
615  // Return result
616  return result;
617 }
618 
619 
620 /*==========================================================================
621  = =
622  = Private methods =
623  = =
624  ==========================================================================*/
625 
626 /***********************************************************************//**
627  * @brief Initialise class members
628  ***************************************************************************/
630 {
631  // Initialise members
632  m_exposure.clear();
633  m_psf.clear();
634  m_edisp.clear();
636  m_apply_edisp = false;
637  m_has_edisp = false;
638 
639  // Initialise cache
640  m_cache.clear();
641 
642  // Return
643  return;
644 }
645 
646 
647 /***********************************************************************//**
648  * @brief Copy class members
649  *
650  * @param[in] rsp Response.
651  ***************************************************************************/
653 {
654  // Copy members
655  m_exposure = rsp.m_exposure;
656  m_psf = rsp.m_psf;
657  m_edisp = rsp.m_edisp;
660  m_has_edisp = rsp.m_has_edisp;
661 
662  // Copy cache
663  for (int i = 0; i < rsp.m_cache.size(); ++i) {
664  m_cache.push_back((rsp.m_cache[i]->clone()));
665  }
666 
667  // Return
668  return;
669 }
670 
671 
672 /***********************************************************************//**
673  * @brief Delete class members
674  ***************************************************************************/
676 {
677  // Free cache
678  for (int i = 0; i < m_cache.size(); ++i) {
679  if (m_cache[i] != NULL) delete m_cache[i];
680  m_cache[i] = NULL;
681  }
682 
683  // Return
684  return;
685 }
686 
687 
688 /***********************************************************************//**
689  * @brief Determine cache index for model
690  *
691  * @param[in] name Model name.
692  * @return Cache index (-1 if model has not been found).
693  *
694  * Determines the cache index for a given model @p name.
695  ***************************************************************************/
696 int GCTAResponseCube::cache_index(const std::string& name) const
697 {
698  // Initialise index
699  int index = -1;
700 
701  // Continue only if there are models in cache
702  if (!m_cache.empty()) {
703 
704  // Search for model name
705  for (int i = 0; i < m_cache.size(); ++i) {
706  if (m_cache[i]->name() == name) {
707  index = i;
708  break;
709  }
710  }
711 
712  } // endif: there were models in cache
713 
714  // Return index
715  return index;
716 }
717 
718 
719 #if defined(G_RADIAL_PSF_BASED)
720 /***********************************************************************//**
721  * @brief Integrate Psf over radial model
722  *
723  * @param[in] model Radial model.
724  * @param[in] delta_mod Angle between model centre and measured photon
725  * direction (radians).
726  * @param[in] obsDir Observed event direction.
727  * @param[in] srcEng True photon energy.
728  * @param[in] srcTime True photon arrival time.
729  *
730  * Integrates the product of the spatial model and the point spread
731  * function over the true photon arrival direction using
732  *
733  * \f[
734  * \int_{\delta_{\rm min}}^{\delta_{\rm max}}
735  * \sin \delta \times {\rm Psf}(\delta) \times
736  * \int_{\phi_{\rm min}}^{\phi_{\rm max}}
737  * S_{\rm p}(\delta, \phi | E, t) d\phi d\delta
738  * \f]
739  *
740  * where
741  * \f$S_{\rm p}(\delta, \phi | E, t)\f$ is the radial spatial model,
742  * \f${\rm Psf}(\delta)\f$ is the point spread function,
743  * \f$\delta\f$ is angular distance between the true and the measured
744  * photon direction, and
745  * \f$\phi\f$ is the position angle around the observed photon direction
746  * measured counterclockwise from the connecting line between the model
747  * centre and the observed photon arrival direction.
748  ***************************************************************************/
750  const double& delta_mod,
751  const GSkyDir& obsDir,
752  const GEnergy& srcEng,
753  const GTime& srcTime) const
754 {
755  // Set number of iterations for Romberg integration.
756  // These values have been determined after careful testing, see
757  // https://cta-redmine.irap.omp.eu/issues/1291
758  static const int iter_delta = 5;
759  static const int iter_phi = 6;
760 
761  // Initialise value
762  double value = 0.0;
763 
764  // Get maximum Psf radius (radians)
765  double psf_max = psf().delta_max();
766 
767  // Get maximum model radius (radians)
768  double theta_max = model->theta_max();
769 
770  // Set offset angle integration range (radians)
771  double delta_min = (delta_mod > theta_max) ? delta_mod - theta_max : 0.0;
772  double delta_max = delta_mod + theta_max;
773  if (delta_max > psf_max) {
774  delta_max = psf_max;
775  }
776 
777  // Setup integration kernel. We take here the observed photon arrival
778  // direction as the true photon arrival direction because the PSF does
779  // not vary significantly over a small region.
780  cta_psf_radial_kern_delta integrand(this,
781  model,
782  obsDir,
783  srcEng,
784  srcTime,
785  delta_mod,
786  theta_max,
787  iter_phi);
788 
789  // Integrate over model's zenith angle
790  GIntegral integral(&integrand);
791  integral.fixed_iter(iter_delta);
792 
793  // Setup integration boundaries
794  std::vector<double> bounds;
795  bounds.push_back(delta_min);
796  bounds.push_back(delta_max);
797 
798  // If the integration range includes a transition between full
799  // containment of Psf within model and partial containment, then
800  // add a boundary at this location
801  double transition_point = theta_max - delta_mod;
802  if (transition_point > delta_min && transition_point < delta_max) {
803  bounds.push_back(transition_point);
804  }
805 
806  // Integrate kernel
807  value = integral.romberg(bounds, iter_delta);
808 
809  // Compile option: Check for NaN/Inf
810  #if defined(G_NAN_CHECK)
811  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
812  std::cout << "*** ERROR: GCTAResponseCube::psf_radial:";
813  std::cout << " NaN/Inf encountered";
814  std::cout << " (value=" << value;
815  std::cout << ", delta_min=" << delta_min;
816  std::cout << ", delta_max=" << delta_max << ")";
817  std::cout << std::endl;
818  }
819  #endif
820 
821  // Return integral
822  return value;
823 }
824 #else
825 /***********************************************************************//**
826  * @brief Integrate Psf over radial model
827  *
828  * @param[in] model Radial model.
829  * @param[in] rho_obs Angle between model centre and measured photon direction
830  * (radians).
831  * @param[in] obsDir Observed event direction.
832  * @param[in] srcEng True photon energy.
833  * @param[in] srcTime True photon arrival time.
834  *
835  * Integrates the product of the spatial model and the point spread
836  * function over the true photon arrival direction using
837  *
838  * \f[
839  * \int_{\rho_{\rm min}}^{\rho_{\rm max}}
840  * \sin \rho \times S_{\rm p}(\rho | E, t) \times
841  * \int_{\omega_{\rm min}}^{\omega_{\rm max}}
842  * {\rm Psf}(\rho, \omega) d\omega d\rho
843  * \f]
844  *
845  * where
846  * \f$S_{\rm p}(\rho | E, t)\f$ is the radial spatial model,
847  * \f${\rm Psf}(\rho, \omega)\f$ is the point spread function,
848  * \f$\rho\f$ is the radial distance from the model centre, and
849  * \f$\omega\f$ is the position angle around the model centre measured
850  * counterclockwise from the connecting line between the model centre and
851  * the observed photon arrival direction.
852  ***************************************************************************/
854  const double& rho_obs,
855  const GSkyDir& obsDir,
856  const GEnergy& srcEng,
857  const GTime& srcTime) const
858 {
859  // Set number of iterations for Romberg integration.
860  // These values have been determined after careful testing, see
861  // https://cta-redmine.irap.omp.eu/issues/1299
862  static const int iter_rho = 5;
863  static const int iter_phi = 5;
864 
865  // Initialise value
866  double irf = 0.0;
867 
868  // Get maximum PSF radius (radians)
869  double delta_max = psf().delta_max();
870 
871  // Set zenith angle integration range for radial model (radians)
872  double rho_min = (rho_obs > delta_max) ? rho_obs - delta_max : 0.0;
873  double rho_max = rho_obs + delta_max;
874  double src_max = model->theta_max();
875  if (rho_max > src_max) {
876  rho_max = src_max;
877  }
878 
879  // Perform zenith angle integration if interval is valid
880  if (rho_max > rho_min) {
881 
882  // Setup integration kernel. We take here the observed photon arrival
883  // direction as the true photon arrival direction because the PSF does
884  // not vary significantly over a small region.
885  cta_psf_radial_kern_rho integrand(this,
886  model,
887  obsDir,
888  srcEng,
889  srcTime,
890  rho_obs,
891  delta_max,
892  iter_phi);
893 
894  // Integrate over model's zenith angle
895  GIntegral integral(&integrand);
896  integral.fixed_iter(iter_rho);
897 
898  // Setup integration boundaries
899  std::vector<double> bounds;
900  bounds.push_back(rho_min);
901  bounds.push_back(rho_max);
902 
903  // If the integration range includes a transition between full
904  // containment of model within Psf and partial containment, then
905  // add a boundary at this location
906  double transition_point = delta_max - rho_obs;
907  if (transition_point > rho_min && transition_point < rho_max) {
908  bounds.push_back(transition_point);
909  }
910 
911  // If we have a shell model then add an integration boundary for the
912  // shell radius as a function discontinuity will occur at this
913  // location
914  const GModelSpatialRadialShell* shell = dynamic_cast<const GModelSpatialRadialShell*>(model);
915  if (shell != NULL) {
916  double shell_radius = shell->radius() * gammalib::deg2rad;
917  if (shell_radius > rho_min && shell_radius < rho_max) {
918  bounds.push_back(shell_radius);
919  }
920  }
921 
922  // Integrate kernel
923  irf = integral.romberg(bounds, iter_rho);
924 
925  // Compile option: Check for NaN/Inf
926  #if defined(G_NAN_CHECK)
928  std::cout << "*** ERROR: GCTAResponseCube::psf_radial:";
929  std::cout << " NaN/Inf encountered";
930  std::cout << " (irf=" << irf;
931  std::cout << ", rho_min=" << rho_min;
932  std::cout << ", rho_max=" << rho_max << ")";
933  std::cout << std::endl;
934  }
935  #endif
936 
937  } // endif: integration interval is valid
938 
939  // Return PSF
940  return irf;
941 }
942 #endif
943 
944 
945 /***********************************************************************//**
946  * @brief Integrate Psf over elliptical model
947  *
948  * @param[in] model Elliptical model.
949  * @param[in] rho_obs Angle between model centre and measured photon direction
950  * (radians).
951  * @param[in] posangle_obs Position angle of measured photon direction with
952  * respect to model centre (radians).
953  * @param[in] obsDir Observed event direction.
954  * @param[in] srcEng True photon energy.
955  * @param[in] srcTime True photon arrival time.
956  *
957  * Integrates the product of the spatial model and the point spread
958  * function over the true photon arrival direction using
959  *
960  * \f[
961  * \int_{\rho_{\rm min}}^{\rho_{\rm max}}
962  * \sin \rho \times
963  * \int_{\omega}
964  * S_{\rm p}(\rho, \omega | E, t) \times
965  * {\rm Psf}(\rho, \omega) d\omega d\rho
966  * \f]
967  *
968  * where
969  * \f$S_{\rm p}(\rho, \omega | E, t)\f$ is the elliptical spatial model,
970  * \f${\rm Psf}(\rho, \omega)\f$ is the point spread function,
971  * \f$\rho\f$ is the radial distance from the model centre, and
972  * \f$\omega\f$ is the position angle around the model centre measured
973  * counterclockwise from the connecting line between the model centre and
974  * the observed photon arrival direction.
975  ***************************************************************************/
977  const double& rho_obs,
978  const double& posangle_obs,
979  const GSkyDir& obsDir,
980  const GEnergy& srcEng,
981  const GTime& srcTime) const
982 {
983  // Set number of iterations for Romberg integration.
984  // These values have been determined after careful testing, see
985  // https://cta-redmine.irap.omp.eu/issues/1299
986  static const int iter_rho = 5;
987  static const int iter_phi = 5;
988 
989  // Initialise value
990  double irf = 0.0;
991 
992  // Get maximum PSF radius (radians)
993  double delta_max = psf().delta_max();
994 
995  // Get the ellipse boundary (radians). Note that these are NOT the
996  // parameters of the ellipse but of a boundary ellipse that is used
997  // for computing the relevant omega angle intervals for a given angle
998  // rho. The boundary ellipse takes care of the possibility that the
999  // semiminor axis is larger than the semimajor axis
1000  double semimajor; // Will be the larger axis
1001  double semiminor; // Will be the smaller axis
1002  double posangle; // Will be the corrected position angle
1003  double aspect_ratio; // Ratio between smaller/larger axis of model
1004  if (model->semimajor() >= model->semiminor()) {
1005  aspect_ratio = (model->semimajor() > 0.0) ?
1006  model->semiminor() / model->semimajor() : 0.0;
1007  posangle = model->posangle() * gammalib::deg2rad;
1008  }
1009  else {
1010  aspect_ratio = (model->semiminor() > 0.0) ?
1011  model->semimajor() / model->semiminor() : 0.0;
1012  posangle = model->posangle() * gammalib::deg2rad + gammalib::pihalf;
1013  }
1014  semimajor = model->theta_max();
1015  semiminor = semimajor * aspect_ratio;
1016 
1017  // Set zenith angle integration range for elliptical model
1018  double rho_min = (rho_obs > delta_max) ? rho_obs - delta_max : 0.0;
1019  double rho_max = rho_obs + delta_max;
1020  if (rho_max > semimajor) {
1021  rho_max = semimajor;
1022  }
1023 
1024  // Perform zenith angle integration if interval is valid
1025  if (rho_max > rho_min) {
1026 
1027  // Setup integration kernel. We take here the observed photon arrival
1028  // direction as the true photon arrival direction because the PSF does
1029  // not vary significantly over a small region.
1030  cta_psf_elliptical_kern_rho integrand(this,
1031  model,
1032  semimajor,
1033  semiminor,
1034  posangle,
1035  obsDir,
1036  srcEng,
1037  srcTime,
1038  rho_obs,
1039  posangle_obs,
1040  delta_max,
1041  iter_phi);
1042 
1043  // Integrate over model's zenith angle
1044  GIntegral integral(&integrand);
1045  integral.fixed_iter(iter_rho);
1046 
1047  // Setup integration boundaries
1048  std::vector<double> bounds;
1049  bounds.push_back(rho_min);
1050  bounds.push_back(rho_max);
1051 
1052  // Kluge: add this transition point as this allows to fit the test
1053  // case without any stalls. Not clear why this is the case, maybe
1054  // simply because the rho integral gets cut down into one more
1055  // sub-interval which may increase precision and smoothed the
1056  // likelihood contour
1057  double transition_point = delta_max - rho_obs;
1058  if (transition_point > rho_min && transition_point < rho_max) {
1059  bounds.push_back(transition_point);
1060  }
1061 
1062  // If the integration range includes the semiminor boundary, then
1063  // add an integration boundary at that location
1064  if (semiminor > rho_min && semiminor < rho_max) {
1065  bounds.push_back(semiminor);
1066  }
1067 
1068  // Integrate kernel
1069  irf = integral.romberg(bounds, iter_rho);
1070 
1071  // Compile option: Check for NaN/Inf
1072  #if defined(G_NAN_CHECK)
1074  std::cout << "*** ERROR: GCTAResponseCube::psf_elliptical:";
1075  std::cout << " NaN/Inf encountered";
1076  std::cout << " (irf=" << irf;
1077  std::cout << ", rho_min=" << rho_min;
1078  std::cout << ", rho_max=" << rho_max << ")";
1079  std::cout << std::endl;
1080  }
1081  #endif
1082 
1083  } // endif: integration interval is valid
1084 
1085  // Return PSF
1086  return irf;
1087 }
1088 
1089 
1090 /***********************************************************************//**
1091  * @brief Integrate PSF over diffuse model
1092  *
1093  * @param[in] model Diffuse spatial model.
1094  * @param[in] obsDir Observed event direction.
1095  * @param[in] srcEng True photon energy.
1096  * @param[in] srcTime True photon arrival time.
1097  *
1098  * Computes the integral
1099  *
1100  * \f[
1101  * \int_0^{\delta_{\rm max}}
1102  * {\rm PSF}(\delta) \times
1103  * \int_0^{2\pi} {\rm Map}(\delta, \phi) \sin \delta
1104  * {\rm d}\phi {\rm d}\delta
1105  * \f]
1106  *
1107  * where \f${\rm Map}(\delta, \phi)\f$ is the diffuse map in the coordinate
1108  * system of the point spread function, defined by the angle \f$\delta\f$
1109  * between the true and the measured photon direction and the azimuth angle
1110  * \f$\phi\f$ around the measured photon direction.
1111  * \f${\rm PSF}(\delta)\f$ is the azimuthally symmetric point spread
1112  * function.
1113  ***************************************************************************/
1115  const GSkyDir& obsDir,
1116  const GEnergy& srcEng,
1117  const GTime& srcTime) const
1118 {
1119  // Set minimum and maximum number of iterations for Romberg integration.
1120  // These values have been determined after careful testing, see
1121  // https://cta-redmine.irap.omp.eu/issues/2458
1122  static const int min_iter_delta = 5;
1123  static const int min_iter_phi = 5;
1124  static const int max_iter_delta = 8;
1125  static const int max_iter_phi = 8;
1126 
1127  // Initialise PSF
1128  double psf = 0.0;
1129 
1130  // Compute rotation matrix to convert from PSF centred coordinate system
1131  // spanned by delta and phi into the reference frame of the observed
1132  // arrival direction given in Right Ascension and Declination.
1133  GMatrix ry;
1134  GMatrix rz;
1135  ry.eulery(obsDir.dec_deg() - 90.0);
1136  rz.eulerz(-obsDir.ra_deg());
1137  GMatrix rot = (ry * rz).transpose();
1138 
1139  // Get offset angle integration interval in radians
1140  double delta_min = 0.0;
1141  double delta_max = 1.1 * this->psf().delta_max();
1142 
1143  // Get resolution of spatial model
1144  double resolution = gammalib::resolution(model);
1145 
1146 
1147  // Setup integration kernel. We take here the observed photon arrival
1148  // direction as the true photon arrival direction because the PSF does
1149  // not vary significantly over a small region.
1150  cta_psf_diffuse_kern_delta integrand(this, model, &rot,
1151  obsDir, srcEng, srcTime,
1152  min_iter_phi, max_iter_phi,
1153  resolution);
1154 
1155  // Set number of radial integration iterations
1156  int iter = gammalib::iter_rho(delta_max, resolution,
1157  min_iter_delta, max_iter_delta);
1158 
1159  // Setup integration
1160  GIntegral integral(&integrand);
1161 
1162  // Set fixed number of iterations
1163  integral.fixed_iter(iter);
1164 
1165  // Integrate over PSF delta angle
1166  psf = integral.romberg(delta_min, delta_max);
1167 
1168  // Return PSF
1169  return psf;
1170 }
1171 
1172 
1173 /***********************************************************************//**
1174  * @brief Return instrument response to point source
1175  *
1176  * @param[in] event Observed event.
1177  * @param[in] source Source.
1178  * @param[in] obs Observation (not used).
1179  * @return Instrument response to point source.
1180  *
1181  * Returns the instrument response to a specified point source.
1182  ***************************************************************************/
1184  const GSource& source,
1185  const GObservation& obs) const
1186 {
1187  // Initialise IRF
1188  double irf = 0.0;
1189 
1190  // Get pointer to model source model
1191  const GModelSpatialPointSource* ptsrc =
1192  static_cast<const GModelSpatialPointSource*>(source.model());
1193 
1194  // Get point source direction
1195  GSkyDir srcDir = ptsrc->dir();
1196 
1197  // Get energy of source model
1198  GEnergy srcEng = source.energy();
1199 
1200  // Get pointer on CTA event bin
1201  if (!event.is_bin()) {
1202  std::string msg = "The current event is not a CTA event bin. "
1203  "This method only works on binned CTA data. Please "
1204  "make sure that a CTA observation containing binned "
1205  "CTA data is provided.";
1207  }
1208  const GCTAEventBin* bin = static_cast<const GCTAEventBin*>(&event);
1209 
1210  // Determine angular separation between true and measured photon
1211  // direction in radians
1212  double delta = bin->dir().dir().dist(srcDir);
1213 
1214  // Get maximum angular separation for PSF (in radians)
1215  double delta_max = psf().delta_max();
1216 
1217  // Get livetime (in seconds)
1218  double livetime = exposure().livetime();
1219 
1220  // Continue only if livetime is >0 and if we're sufficiently close
1221  // to the PSF
1222  if ((livetime > 0.0) && (delta <= delta_max)) {
1223 
1224  // Get exposure
1225  irf = exposure()(srcDir, srcEng);
1226 
1227  // Multiply-in PSF
1228  if (irf > 0.0) {
1229 
1230  // Recover effective area from exposure
1231  irf /= livetime;
1232 
1233  // Get PSF component
1234  irf *= psf()(srcDir, delta, srcEng);
1235 
1236  // Multiply-in energy dispersion
1237  if (use_edisp() && irf > 0.0) {
1238  irf *= edisp()(bin->energy(), srcEng, srcDir);
1239  }
1240 
1241  // Apply deadtime correction
1242  irf *= exposure().deadc();
1243 
1244  } // endif: exposure was non-zero
1245 
1246  } // endif: we were sufficiently close to PSF and livetime >0
1247 
1248  // Compile option: Check for NaN/Inf
1249  #if defined(G_NAN_CHECK)
1251  std::cout << "*** ERROR: GCTAResponseCube::irf_ptsrc:";
1252  std::cout << " NaN/Inf encountered";
1253  std::cout << " irf=" << irf;
1254  std::cout << std::endl;
1255  }
1256  #endif
1257 
1258  // Return IRF value
1259  return irf;
1260 }
1261 
1262 
1263 /***********************************************************************//**
1264  * @brief Return instrument response to radial source
1265  *
1266  * @param[in] event Observed event.
1267  * @param[in] source Source.
1268  * @param[in] obs Observation (not used).
1269  * @return Instrument response to radial source.
1270  *
1271  * Returns the instrument response to a specified radial source.
1272  *
1273  * @todo Correct assumptions in exposure computation, PSF computation and
1274  * energy dispersion computation.
1275  ***************************************************************************/
1277  const GSource& source,
1278  const GObservation& obs) const
1279 {
1280  // Initialise IRF
1281  double irf = 0.0;
1282 
1283  // Get pointer to CTA event bin
1284  if (!event.is_bin()) {
1285  std::string msg = "The current event is not a CTA event bin. "
1286  "This method only works on binned CTA data. Please "
1287  "make sure that a CTA observation containing binned "
1288  "CTA data is provided.";
1290  }
1291  const GCTAEventBin* bin = static_cast<const GCTAEventBin*>(&event);
1292 
1293  // Get event attribute references
1294  const GSkyDir& obsDir = bin->dir().dir();
1295  const GEnergy& obsEng = bin->energy();
1296  const GTime& obsTime = bin->time();
1297 
1298  // Get energy of source model
1299  GEnergy srcEng = source.energy();
1300 
1301  // Get pointer to radial model
1302  const GModelSpatialRadial* model =
1303  static_cast<const GModelSpatialRadial*>(source.model());
1304 
1305  // Compute angle between model centre and measured photon direction
1306  // (radians)
1307  double rho_obs = model->dir().dist(obsDir);
1308 
1309  // Get livetime (in seconds)
1310  double livetime = exposure().livetime();
1311 
1312  // Continue only if livetime is >0 and if we're sufficiently close to
1313  // the model centre to get a non-zero response
1314  if ((livetime > 0.0) && (rho_obs <= model->theta_max()+psf().delta_max())) {
1315 
1316  // Get exposure at the observed event direction.
1317  //
1318  // The current code assumes that the exposure at the observed and
1319  // true event direction does not vary significantly. In other words,
1320  // the code assumes that the exposure is constant over the size of
1321  // the PSF.
1322  irf = exposure()(obsDir, srcEng);
1323 
1324  // Continue only if exposure is positive
1325  if (irf > 0.0) {
1326 
1327  // Recover effective area from exposure
1328  irf /= livetime;
1329 
1330  // Get PSF component
1331  irf *= psf_radial(model, rho_obs, obsDir, srcEng, obsTime);
1332 
1333  // Multiply-in energy dispersion
1334  //
1335  // The current code assumes that the energy dispersion at the
1336  // observed and true event direction does not vary
1337  // significantly. In other words, the code assumes that the
1338  // energy dispersion is constant over the size of the PSF.
1339  if (use_edisp() && irf > 0.0) {
1340  irf *= edisp()(bin->energy(), srcEng, obsDir);
1341  }
1342 
1343  // Apply deadtime correction
1344  irf *= exposure().deadc();
1345 
1346  } // endif: exposure was positive
1347 
1348  } // endif: we were sufficiently close and livetime >0
1349 
1350  // Compile option: Check for NaN/Inf
1351  #if defined(G_NAN_CHECK)
1353  std::cout << "*** ERROR: GCTAResponseCube::irf_radial:";
1354  std::cout << " NaN/Inf encountered";
1355  std::cout << " irf=" << irf;
1356  std::cout << std::endl;
1357  }
1358  #endif
1359 
1360  // Return IRF value
1361  return irf;
1362 }
1363 
1364 
1365 /***********************************************************************//**
1366  * @brief Return instrument response to elliptical source
1367  *
1368  * @param[in] event Observed event.
1369  * @param[in] source Source.
1370  * @param[in] obs Observation (not used).
1371  * @return Instrument response to elliptical source.
1372  *
1373  * Returns the instrument response to a specified elliptical source.
1374  *
1375  * @todo Correct assumptions in exposure computation, PSF computation and
1376  * energy dispersion computation.
1377  ***************************************************************************/
1379  const GSource& source,
1380  const GObservation& obs) const
1381 {
1382  // Initialise IRF
1383  double irf = 0.0;
1384 
1385  // Get pointer to CTA event bin
1386  if (!event.is_bin()) {
1387  std::string msg = "The current event is not a CTA event bin. "
1388  "This method only works on binned CTA data. Please "
1389  "make sure that a CTA observation containing binned "
1390  "CTA data is provided.";
1392  }
1393  const GCTAEventBin* bin = static_cast<const GCTAEventBin*>(&event);
1394 
1395  // Get event attribute references
1396  const GSkyDir& obsDir = bin->dir().dir();
1397  const GEnergy& obsEng = bin->energy();
1398  const GTime& obsTime = bin->time();
1399 
1400  // Get energy of source model
1401  GEnergy srcEng = source.energy();
1402 
1403  // Get pointer to elliptical model
1404  const GModelSpatialElliptical* model =
1405  static_cast<const GModelSpatialElliptical*>(source.model());
1406 
1407  // Compute angle between model centre and measured photon direction and
1408  // position angle (radians)
1409  double rho_obs = model->dir().dist(obsDir);
1410  double posangle_obs = model->dir().posang(obsDir);
1411 
1412  // Get livetime (in seconds)
1413  double livetime = exposure().livetime();
1414 
1415  // Continue only if livetime is >0 and if we're sufficiently close to
1416  // the model centre to get a non-zero response
1417  if ((livetime > 0.0) && (rho_obs <= model->theta_max()+psf().delta_max())) {
1418 
1419  // Get exposure
1420  //
1421  // The current code assumes that the exposure at the observed and
1422  // true event direction does not vary significantly. In other words,
1423  // the code assumes that the exposure is constant over the size of
1424  // the PSF.
1425  irf = exposure()(obsDir, srcEng);
1426 
1427  // Continue only if exposure is positive
1428  if (irf > 0.0) {
1429 
1430  // Recover effective area from exposure
1431  irf /= livetime;
1432 
1433  // Get PSF component
1434  irf *= psf_elliptical(model, rho_obs, posangle_obs, obsDir, srcEng, obsTime);
1435 
1436  // Multiply-in energy dispersion
1437  //
1438  // The current code assumes that the energy dispersion at the
1439  // observed and true event direction does not vary
1440  // significantly. In other words, the code assumes that the
1441  // energy dispersion is constant over the size of the PSF.
1442  if (use_edisp() && irf > 0.0) {
1443  irf *= edisp()(bin->energy(), srcEng, obsDir);
1444  }
1445 
1446  // Apply deadtime correction
1447  irf *= exposure().deadc();
1448 
1449  } // endif: exposure was positive
1450 
1451  } // endif: we were sufficiently close and livetime >0
1452 
1453  // Compile option: Check for NaN/Inf
1454  #if defined(G_NAN_CHECK)
1456  std::cout << "*** ERROR: GCTAResponseCube::irf_elliptical:";
1457  std::cout << " NaN/Inf encountered";
1458  std::cout << " irf=" << irf;
1459  std::cout << std::endl;
1460  }
1461  #endif
1462 
1463  // Return IRF value
1464  return irf;
1465 }
1466 
1467 
1468 /***********************************************************************//**
1469  * @brief Return instrument response to diffuse source
1470  *
1471  * @param[in] event Observed event.
1472  * @param[in] source Source.
1473  * @param[in] obs Observation.
1474  * @return Instrument response to diffuse source.
1475  *
1476  * Returns the instrument response to a specified diffuse source.
1477  *
1478  * The method uses a pre-computation cache to store the instrument response
1479  * for the spatial model component. The pre-computation cache is initialised
1480  * if no cache has yet been allocated, or if at the beginning of a scan over
1481  * the events, the model parameters have changed. The beginning of a scan is
1482  * defined by an event bin index of 0.
1483  ***************************************************************************/
1485  const GSource& source,
1486  const GObservation& obs) const
1487 {
1488  // Initialise IRF
1489  double irf = 0.0;
1490 
1491  // Get pointer to CTA event bin
1492  if (!event.is_bin()) {
1493  std::string msg = "The current event is not a CTA event bin. "
1494  "This method only works on binned CTA data. Please "
1495  "make sure that a CTA observation containing binned "
1496  "CTA data is provided.";
1498  }
1499  const GCTAEventBin* bin = static_cast<const GCTAEventBin*>(&event);
1500 
1501  // Get event attribute references
1502  const GSkyDir& obsDir = bin->dir().dir();
1503  const GEnergy& obsEng = bin->energy();
1504  const GTime& obsTime = bin->time();
1505 
1506  // Get energy of source model
1507  GEnergy srcEng = source.energy();
1508 
1509  // Get pointer to elliptical model
1510  const GModelSpatialDiffuse* model =
1511  static_cast<const GModelSpatialDiffuse*>(source.model());
1512 
1513  // Get pointer on CTA response cube
1514  //const GCTAResponseCube* rsp = gammalib::cta_rsp_cube(G_IRF_DIFFUSE, obs);
1515 
1516  // Get livetime (in seconds) and deadtime correction factor
1517  double livetime = exposure().livetime();
1518  double deadc = exposure().deadc();
1519 
1520  // Get Psf radius (in degrees)
1521  double delta_max = psf().delta_max() * gammalib::rad2deg;
1522 
1523  // Continue only if livetime is >0 and model contains reconstructed
1524  // sky direction
1525  if (livetime > 0.0 && model->contains(obsDir, delta_max)) {
1526 
1527  // Get exposure
1528  //
1529  // The current code assumes that the exposure at the observed and
1530  // true event direction does not vary significantly. In other words,
1531  // the code assumes that the exposure is constant over the size of
1532  // the PSF.
1533  irf = exposure()(obsDir, srcEng);
1534 
1535  // Continue only if exposure is positive
1536  if (irf > 0.0) {
1537 
1538  // Recover effective area from exposure
1539  irf /= livetime;
1540 
1541  // Compute product of PSF and diffuse map, integrated over the
1542  // relevant PSF area. We assume no energy dispersion and thus
1543  // compute the product using the observed energy.
1544  irf *= psf_diffuse(model, obsDir, srcEng, obsTime);
1545 
1546  // Multiply-in energy dispersion
1547  //
1548  // The current code assumes that the energy dispersion at the
1549  // observed and true event direction does not vary
1550  // significantly. In other words, the code assumes that the
1551  // energy dispersion is constant over the size of the PSF.
1552  if (use_edisp() && irf > 0.0) {
1553  irf *= edisp()(bin->energy(), srcEng, obsDir);
1554  }
1555 
1556  // Apply deadtime correction
1557  irf *= exposure().deadc();
1558 
1559  } // endif: exposure was positive
1560 
1561  } // endif: livetime was positive
1562 
1563  // Compile option: Check for NaN/Inf
1564  #if defined(G_NAN_CHECK)
1566  std::cout << "*** ERROR: GCTAResponseCube::irf_diffuse:";
1567  std::cout << " NaN/Inf encountered";
1568  std::cout << " irf=" << irf;
1569  std::cout << std::endl;
1570  }
1571  #endif
1572 
1573  // Return IRF value
1574  return irf;
1575 }
const GFilename & filename(void) const
Return edisp cube filename.
#define G_IRF_DIFFUSE
GSkyDir dir(void) const
Return position of elliptical spatial model.
void copy_members(const GCTAResponseCube &rsp)
Copy class members.
std::string print(const GChatter &chatter=NORMAL) const
Print energy dispersion cube information.
CTA response helper classes definition.
CTA exposure cube class.
double irf_radial(const GEvent &event, const GSource &source, const GObservation &obs) const
Return instrument response to radial source.
double dec_deg(void) const
Returns Declination in degrees.
Definition: GSkyDir.hpp:250
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:380
Point source spatial model class interface definition.
GCTACubePsf m_psf
Mean point spread function.
virtual bool use_edisp(void) const
Signal if response uses energy dispersion.
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 cube-style response function class definition.
Energy value class definition.
Abstract elliptical spatial model base class.
GCTACubeExposure m_exposure
Exposure cube.
virtual const GCTAInstDir & dir(void) const
Return instrument direction of event bin.
void load(const GFilename &filename)
Load PSF cube from FITS file.
GCTAResponseCache m_irf_cache
void clear(void)
Clear instance.
Sky direction class interface definition.
const GFilename & filename(void) const
Return background cube filename.
#define G_IRF_RADIAL
GCTAEventBin class interface definition.
const GCTACubeBackground & background(void) const
Return cube analysis background cube.
CTA cube background class.
Abstract radial spatial model base class interface definition.
CTA event bin class interface definition.
CTA cube-style response function class.
void clear(void)
Clear instance.
Abstract interface for the event classes.
Definition: GEvent.hpp:71
XML element node class.
Definition: GXmlElement.hpp:47
double psf_radial(const GModelSpatialRadial *model, const double &rho_obs, const GSkyDir &obsDir, const GEnergy &srcEng, const GTime &srcTime) const
Integrate Psf over radial model.
#define G_NROI
std::string print(const GChatter &chatter=NORMAL) const
Print background information.
virtual GCTAResponseCube & operator=(const GCTAResponseCube &rsp)
Assignment operator.
void load(const GFilename &filename)
Load energy dispersion cube from FITS file.
Time class.
Definition: GTime.hpp:54
Gammalib tools definition.
GIntegral class interface definition.
Definition: GIntegral.hpp:45
virtual const GEnergy & energy(void) const
Return energy of event bin.
double irf_diffuse(const GEvent &event, const GSource &source, const GObservation &obs) const
Return instrument response to diffuse source.
virtual void write(GXmlElement &xml) const
Write response information into XML element.
void init_members(void)
Initialise class members.
Definition: GResponse.cpp:260
std::vector< GCTACubeSource * > m_cache
Response cache.
GFilename xml_file_reduce(const GXmlElement &xml, const std::string &filename)
Reduce file name provided for writing as XML attribute.
Definition: GTools.cpp:1713
Source class definition.
void dir(const GSkyDir &dir)
Set sky direction.
void init_members(void)
Initialise class members.
virtual std::string print(const GChatter &chatter=NORMAL) const
Print response information.
void eulery(const double &angle)
Set Euler rotation matrix around y axis.
Definition: GMatrix.cpp:1093
virtual double theta_max(void) const =0
std::string print(const GChatter &chatter=NORMAL) const
Print exposure cube information.
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:184
void load(const GFilename &filename)
Load exposure cube from FITS file.
void init_members(void)
Initialise class members.
Class that handles photons.
Definition: GPhoton.hpp:47
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:167
virtual GCTAResponse & operator=(const GCTAResponse &rsp)
Assignment operator.
Definition of support function used by CTA classes.
double radius(void) const
Return shell radius.
const GCTACubeExposure & exposure(void) const
Return exposure cube.
double semiminor(void) const
Return semi-minor axis of ellipse.
virtual bool is_bin(void) const =0
const double pihalf
Definition: GMath.hpp:38
const double deg2rad
Definition: GMath.hpp:43
void clear(void)
Clear background.
Energy boundaries container class.
Definition: GEbounds.hpp:60
const GXmlAttribute * attribute(const int &index) const
Return attribute.
GEbounds ebounds(const GEnergy &obsEng) const
Return boundaries in true energy.
virtual ~GCTAResponseCube(void)
Destructor.
const GCTACubePsf & psf(void) const
Return cube analysis point spread function.
Kernel for Psf delta angle integration used for stacked analysis.
GCTACubeBackground m_background
Background cube.
void eulerz(const double &angle)
Set Euler rotation matrix around z axis.
Definition: GMatrix.cpp:1130
virtual const GTime & time(void) const
Return time of event bin.
bool m_has_edisp
Flag to indicate if energy.
GXmlElement * xml_need_par(const std::string &origin, GXmlElement &xml, const std::string &name)
Return pointer to parameter with given name in XML element.
Definition: GTools.cpp:1513
double psf_elliptical(const GModelSpatialElliptical *model, const double &rho_obs, const double &posangle_obs, const GSkyDir &obsDir, const GEnergy &srcEng, const GTime &srcTime) const
Integrate Psf over elliptical model.
const GFilename & filename(void) const
Return exposure cube filename.
Abstract event base class definition.
void free_members(void)
Delete class members.
const GCTACubeEdisp & edisp(void) const
Return cube analysis energy dispersion cube.
void free_members(void)
Delete class members.
virtual GCTAResponseCube * clone(void) const
Clone instance.
CTA instrument direction class interface definition.
#define G_READ
double ra_deg(void) const
Returns Right Ascension in degrees.
Definition: GSkyDir.hpp:223
virtual GEbounds ebounds(const GEnergy &obsEng) const
Return true energy boundaries for a specific observed energy.
double resolution(const GModelSpatial *model)
Determine resolution of spatial model.
CTA energy dispersion for stacked analysis.
#define G_IRF
GChatter
Definition: GTypemaps.hpp:33
double deadc(void) const
Return deadtime correction.
virtual void clear(void)
Clear instance.
Spatial composite model class interface 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
double psf_diffuse(const GModelSpatial *model, const GSkyDir &obsDir, const GEnergy &srcEng, const GTime &srcTime) const
Integrate PSF over diffuse model.
virtual bool contains(const GSkyDir &dir, const double &margin=0.0) const =0
Photon class definition.
Abstract observation base class interface definition.
double irf_ptsrc(const GEvent &event, const GSource &source, const GObservation &obs) const
Return instrument response to point source.
Point source spatial model.
void load(const GFilename &filename)
Load background cube from FITS file.
bool xml_has_par(const GXmlElement &xml, const std::string &name)
Checks if parameter with given name in XML element exists.
Definition: GTools.cpp:1475
Kernel for Psf delta angle integration used for stacked analysis.
CTA instrument response function class.
void free_members(void)
Delete class members.
Definition: GResponse.cpp:282
Abstract elliptical spatial model base class interface definition.
Abstract diffuse spatial model base class interface definition.
int iter_phi(const double &rho, const double &resolution, const int &iter_min, const int &iter_max)
Determine number of azimuthal Romberg iterations.
GSkyDir dir(void) const
Return position of point source.
CTA point spread function for cube analysis.
Definition: GCTACubePsf.hpp:62
Sky model class.
Definition: GModelSky.hpp:120
GSkyDir dir(void) const
Return position of radial spatial model.
std::string print(const GChatter &chatter=NORMAL) const
Print PSF cube information.
double posangle(void) const
Return Position Angle of model.
std::string print(const GChatter &chatter=NORMAL) const
Print CTA response cache.
#define G_IRF_PTSRC
bool m_apply_edisp
Apply energy dispersion.
GCTACubeEdisp m_edisp
Energy dispersion cube.
GChatter reduce(const GChatter &chatter)
Reduce chattiness by one level.
Definition: GTypemaps.hpp:65
double delta_max(void) const
Return maximum delta value in radians.
const double & livetime(void) const
Return livetime.
virtual void read(const GXmlElement &xml)
Read response information from XML element.
const GEnergy & energy(void) const
Return photon energy.
Definition: GPhoton.hpp:122
virtual double irf(const GEvent &event, const GPhoton &photon, const GObservation &obs) const
Return instrument response.
virtual double theta_max(void) const =0
double dist(const GSkyDir &dir) const
Compute angular distance between sky directions in radians.
Definition: GSkyDir.hpp:265
void fixed_iter(const int &iter)
Set fixed number of iterations.
Definition: GIntegral.hpp:181
CTA instrument direction class.
Definition: GCTAInstDir.hpp:59
Abstract spatial model base class.
Abstract radial spatial model base class.
GCTAResponseCube(void)
Void constructor.
Generic matrix class definition.
Definition: GMatrix.hpp:79
Abstract diffuse spatial model base class.
const double rad2deg
Definition: GMath.hpp:44
double irf_elliptical(const GEvent &event, const GSource &source, const GObservation &obs) const
Return instrument response to elliptical source.
virtual bool apply_edisp(void) const
Signal if energy dispersion should be applied.
std::string parformat(const std::string &s, const int &indent=0)
Convert string in parameter format.
Definition: GTools.cpp:1022
void clear(void)
Clear energy dispersion cube.
double posang(const GSkyDir &dir) const
Compute position angle between sky directions in radians.
Definition: GSkyDir.cpp:1022
Class that handles gamma-ray sources.
Definition: GSource.hpp:53
int cache_index(const std::string &name) const
Determine cache index for model.
Integration class interface definition.
const GFilename & filename(void) const
Return exposure cube filename.
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_WRITE
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.
double semimajor(void) const
Return semi-major axis of ellipse.
const GSkyDir & dir(void) const
Return photon sky direction.
Definition: GPhoton.hpp:110
Time class interface definition.
const GXmlElement * xml_get_par(const std::string &origin, const GXmlElement &xml, const std::string &name)
Return pointer to parameter with given name in XML element.
Definition: GTools.cpp:1562
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
GFilename xml_file_expand(const GXmlElement &xml, const std::string &filename)
Expand file name provided as XML attribute for loading.
Definition: GTools.cpp:1683