GammaLib  2.1.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-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 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 "GPhoton.hpp"
35 #include "GSource.hpp"
36 #include "GEvent.hpp"
37 #include "GSkyDir.hpp"
38 #include "GEnergy.hpp"
39 #include "GEnergies.hpp"
40 #include "GTime.hpp"
41 #include "GIntegral.hpp"
42 #include "GIntegrals.hpp"
43 #include "GObservation.hpp"
44 #include "GNdarray.hpp"
45 #include "GMatrixSparse.hpp"
47 #include "GModelSpatialRadial.hpp"
51 #include "GModelSpatialDiffuse.hpp"
52 #include "GCTAResponseCube.hpp"
53 #include "GCTAResponse_helpers.hpp"
54 #include "GCTAInstDir.hpp"
55 #include "GCTAEventBin.hpp"
56 #include "GCTASupport.hpp"
57 #include "GCTAEventCube.hpp" // Kludge
58 #include "GCTAEventBin.hpp" // Kludge
60 
61 /* __ Method name definitions ____________________________________________ */
62 #define G_IRF "GCTAResponseCube::irf(GEvent&, GPhoton& GObservation&)"
63 #define G_IRF_PTSRC "GCTAResponseCube::irf_ptsrc(GEvent&, GSource&,"\
64  " GObservation&)"
65 #define G_IRF_RADIAL "GCTAResponseCube::irf_radial(GEvent&, GSource&,"\
66  " GObservation&)"
67 #define G_IRF_DIFFUSE "GCTAResponseCube::irf_diffuse(GEvent&, GSource&,"\
68  " GObservation&)"
69 #define G_NROI "GCTAResponseCube::nroi(GModelSky&, GEnergy&, GTime&, "\
70  "GObservation&)"
71 #define G_EBOUNDS "GCTAResponseCube::ebounds(GEnergy&)"
72 #define G_READ "GCTAResponseCube::read(GXmlElement&)"
73 #define G_WRITE "GCTAResponseCube::write(GXmlElement&)"
74 #define G_IRF_RADIAL2 "GCTAResponseCube::irf_radial(GModelSky&, "\
75  "GObservation&, GMatrixSparse*)"
76 
77 /* __ Macros _____________________________________________________________ */
78 
79 /* __ Coding definitions _________________________________________________ */
80 
81 /* __ Debug definitions __________________________________________________ */
82 
83 /* __ Constants __________________________________________________________ */
84 
85 
86 /*==========================================================================
87  = =
88  = Constructors/destructors =
89  = =
90  ==========================================================================*/
91 
92 /***********************************************************************//**
93  * @brief Void constructor
94  *
95  * Constructs void CTA response.
96  ***************************************************************************/
98 {
99  // Initialise members
100  init_members();
101 
102  // Return
103  return;
104 }
105 
106 
107 /***********************************************************************//**
108  * @brief Copy constructor
109  *
110  * @param[in] rsp CTA response.
111  *
112  * Constructs CTA cube response by making a deep copy of an existing
113  * object.
114  **************************************************************************/
116  GCTAResponse(rsp)
117 {
118  // Initialise members
119  init_members();
120 
121  // Copy members
122  copy_members(rsp);
123 
124  // Return
125  return;
126 }
127 
128 
129 /***********************************************************************//**
130  * @brief XML constructor
131  *
132  * @param[in] xml XML element.
133  *
134  * Construct CTA response from XML element.
135  ***************************************************************************/
137 {
138  // Initialise members
139  init_members();
140 
141  // Read information from XML element
142  read(xml);
143 
144  // Return
145  return;
146 }
147 
148 
149 /***********************************************************************//**
150  * @brief Response constructor
151  *
152  * @param[in] exposure CTA cube analysis exposure.
153  * @param[in] psf CTA cube analysis point spread function.
154  * @param[in] background CTA cube background response.
155  *
156  * Constructs CTA cube analysis response from a cube analysis exposure,
157  * a point spread function cube and a background cube.
158  **************************************************************************/
160  const GCTACubePsf& psf,
161  const GCTACubeBackground& background) :
162  GCTAResponse()
163 {
164  // Initialise members
165  init_members();
166 
167  // Set members
169  m_psf = psf;
171 
172  // Signal that no energy dispersion was given
173  m_has_edisp = false;
174 
175  // Return
176  return;
177 }
178 
179 
180 /***********************************************************************//**
181  * @brief Response constructor
182  *
183  * @param[in] exposure CTA cube analysis exposure.
184  * @param[in] psf CTA cube analysis point spread function.
185  * @param[in] edisp CTA cube energy dispersion response.
186  * @param[in] background CTA cube background response.
187  *
188  * Constructs CTA cube analysis response from a cube analysis exposure,
189  * a point spread function cube, an energy dispersion cube and a background cube.
190  **************************************************************************/
192  const GCTACubePsf& psf,
193  const GCTACubeEdisp& edisp,
194  const GCTACubeBackground& background) :
195  GCTAResponse()
196 {
197  // Initialise members
198  init_members();
199 
200  // Set members
202  m_psf = psf;
203  m_edisp = edisp;
205 
206  // Signal that edisp is available
207  m_has_edisp = true;
208 
209  // Return
210  return;
211 }
212 
213 
214 /***********************************************************************//**
215  * @brief Destructor
216  *
217  * Destroys instance of CTA response object.
218  ***************************************************************************/
220 {
221  // Free members
222  free_members();
223 
224  // Return
225  return;
226 }
227 
228 
229 /*==========================================================================
230  = =
231  = Operators =
232  = =
233  ==========================================================================*/
234 
235 /***********************************************************************//**
236  * @brief Assignment operator
237  *
238  * @param[in] rsp CTA response.
239  * @return CTA response.
240  *
241  * Assigns CTA response object to another CTA response object. The assignment
242  * performs a deep copy of all information, hence the original object from
243  * which the assignment has been performed can be destroyed after this
244  * operation without any loss of information.
245  ***************************************************************************/
247 {
248  // Execute only if object is not identical
249  if (this != &rsp) {
250 
251  // Copy base class members
252  this->GCTAResponse::operator=(rsp);
253 
254  // Free members
255  free_members();
256 
257  // Initialise members
258  init_members();
259 
260  // Copy members
261  copy_members(rsp);
262 
263  } // endif: object was not identical
264 
265  // Return this object
266  return *this;
267 }
268 
269 
270 /*==========================================================================
271  = =
272  = Public methods =
273  = =
274  ==========================================================================*/
275 
276 /***********************************************************************//**
277  * @brief Clear instance
278  *
279  * Clears CTA response object by resetting all members to an initial state.
280  * Any information that was present in the object before will be lost.
281  ***************************************************************************/
283 {
284  // Free class members (base and derived classes, derived class first)
285  free_members();
287  this->GResponse::free_members();
288 
289  // Initialise members
290  this->GResponse::init_members();
292  init_members();
293 
294  // Return
295  return;
296 }
297 
298 
299 /***********************************************************************//**
300  * @brief Clone instance
301  *
302  * @return Pointer to deep copy of CTA response.
303  *
304  * Creates a clone (deep copy) of a CTA response object.
305  ***************************************************************************/
307 {
308  return new GCTAResponseCube(*this);
309 }
310 
311 
312 /***********************************************************************//**
313  * @brief Return instrument response
314  *
315  * @param[in] event Observed event.
316  * @param[in] photon Incident photon.
317  * @param[in] obs Observation (not used).
318  * @return Instrument response.
319  ***************************************************************************/
320 double GCTAResponseCube::irf(const GEvent& event,
321  const GPhoton& photon,
322  const GObservation& obs) const
323 {
324  // Retrieve event instrument direction
325  const GCTAInstDir& dir = gammalib::cta_dir(G_IRF, event);
326 
327  // Get event attributes
328  const GSkyDir& obsDir = dir.dir();
329  const GEnergy& obsEng = event.energy();
330 
331  // Get photon attributes
332  const GSkyDir& srcDir = photon.dir();
333  const GEnergy& srcEng = photon.energy();
334  //const GTime& srcTime = photon.time();
335 
336  // Determine angular separation between true and measured photon
337  // direction in radians
338  double delta = obsDir.dist(srcDir);
339 
340  // Get maximum angular separation for PSF (in radians)
341  double delta_max = psf().delta_max();
342 
343  // Initialise IRF value
344  double irf = 0.0;
345 
346  // Get livetime (in seconds)
347  double livetime = exposure().livetime();
348 
349  // Continue only if livetime is >0 and if we're sufficiently close
350  // to the PSF
351  if ((livetime > 0.0) && (delta <= delta_max)) {
352 
353  // Get exposure
354  irf = exposure()(srcDir, srcEng);
355 
356  // Multiply-in PSF
357  if (irf > 0.0) {
358 
359  // Get PSF component
360  irf *= psf()(srcDir, delta, srcEng);
361 
362  // Multiply-in energy dispersion
363  if (use_edisp() && irf > 0.0) {
364  irf *= edisp()(obsEng, srcEng, srcDir);
365  }
366 
367  // Divide by livetime
368  irf /= livetime;
369 
370  // Apply deadtime correction
371  irf *= exposure().deadc();
372 
373  } // endif: Aeff was non-zero
374 
375  } // endif: we were sufficiently close to PSF and livetime was >0
376 
377  // Compile option: Check for NaN/Inf
378  #if defined(G_NAN_CHECK)
380  std::cout << "*** ERROR: GCTAResponseCube::irf:";
381  std::cout << " NaN/Inf encountered";
382  std::cout << " irf=" << irf;
383  std::cout << std::endl;
384  }
385  #endif
386 
387  // Return IRF value
388  return irf;
389 }
390 
391 
392 /***********************************************************************//**
393  * @brief Return integral of event probability for a given sky model over ROI
394  *
395  * @param[in] model Sky model.
396  * @param[in] obsEng Observed photon energy.
397  * @param[in] obsTime Observed photon arrival time.
398  * @param[in] obs Observation.
399  * @return 0.0
400  *
401  * @exception GException::feature_not_implemented
402  * Method is not implemented.
403  *
404  * Computes the integral
405  *
406  * \f[
407  * N_{\rm ROI}(E',t') = \int_{\rm ROI} P(p',E',t') dp'
408  * \f]
409  *
410  * of the event probability
411  *
412  * \f[
413  * P(p',E',t') = \int \int \int
414  * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp \, dE \, dt
415  * \f]
416  *
417  * for a given sky model \f$S(p,E,t)\f$ and response function
418  * \f$R(p',E',t'|p,E,t)\f$ over the Region of Interest (ROI).
419  *
420  * @todo Implement method (is maybe not really needed)
421  ***************************************************************************/
422 double GCTAResponseCube::nroi(const GModelSky& model,
423  const GEnergy& obsEng,
424  const GTime& obsTime,
425  const GObservation& obs) const
426 {
427  // Method is not implemented
428  std::string msg = "Spatial integration of sky model over the data space "
429  "is not implemented.";
431 
432  // Return Npred
433  return (0.0);
434 }
435 
436 
437 /***********************************************************************//**
438  * @brief Return true energy boundaries for a specific observed energy
439  *
440  * @param[in] obsEng Observed photon energy.
441  * @return Boundaries in true energy.
442  ***************************************************************************/
444 {
445  // Get energy boundaries from energy dispersion
446  GEbounds ebounds = m_edisp.ebounds(obsEng);
447 
448  // Return energy boundaries
449  return ebounds;
450 }
451 
452 
453 /***********************************************************************//**
454  * @brief Return instrument response integrated over the spatial model
455  *
456  * @param[in] event Event.
457  * @param[in] source Source.
458  * @param[in] obs Observation.
459  * @return Instrument response to a spatial model.
460  *
461  * Returns the instrument response for a given event, source and observation
462  * integrated over the spatial model component. The method computes
463  *
464  * \f[
465  * {\tt irf}(p', E', t') = \int_p M_{\rm S}(p | E, t) \,
466  * R(p', E', t' | p, E, t) \, d\,p
467  * \f]
468  *
469  * where
470  * * \f$M_{\rm S}(p | E, t)\f$ is the spatial model component,
471  * * \f$R(p', E', t' | p, E, t)\f$ is the Instrument Response Function (IRF),
472  * * \f$p'\f$ is the measured instrument direction,
473  * * \f$E'\f$ is the measured or reconstructed energy,
474  * * \f$t'\f$ is the measured arrival time,
475  * * \f$p\f$ is the true photon arrival direction,
476  * * \f$E\f$ is the true photon energy, and
477  * * \f$t\f$ is the true trigger time.
478  *
479  * The integration is done over all relevant true sky directions \f$p\f$.
480  *
481  * Depending on the type of the source model the method branches to the
482  * following methods to perform the actual computations
483  *
484  * irf_ptsrc() - for the handling of a point source
485  * irf_radial() - for radial models
486  * irf_elliptical() - for elliptical models
487  * irf_diffuse() - for diffuse models
488  * irf_composite() - for composite models
489  *
490  * The method implements a caching mechanism for spatial models that have all
491  * parameters fixed. For those models the instrument response for a given
492  * event and observation is only computed once and then stored in an internal
493  * cache from which it is fetched back in case that the method is called
494  * again for the same event and observation.
495  ***************************************************************************/
497  const GSource& source,
498  const GObservation& obs) const
499 {
500  // Initialise IRF value
501  double irf = 0.0;
502 
503  // Set IRF value attributes
504  std::string name = obs.id() + "::" + source.name();
505  const GInstDir& dir = event.dir();
506  const GEnergy& ereco = event.energy();
507  const GEnergy& etrue = source.energy();
508 
509  // Signal if spatial model has free parameters
510  bool has_free_pars = source.model()->has_free_pars();
511 
512  // Kludge: the IRF response cache should be used, the model has no
513  // free parameters and there is no energy dispersion
514  if (m_use_irf_cache && !has_free_pars && !use_edisp()) {
515 
516  // Build unique cache name
517  std::string name = obs.id() + "::" + source.name();
518 
519  // Get index in cache, returns -1 if name is not found in cache
520  int index = -1;
521  for (int i = 0; i < m_cache_names.size(); ++i) {
522  if (m_cache_names[i] == name) {
523  index = i;
524  break;
525  }
526  }
527 
528  // If index was not found then allocate a new cache map
529  if (index == -1) {
530 
531  // Get pointer to event cube
532  const GCTAEventCube* cube =
533  static_cast<const GCTAEventCube*>(obs.events());
534 
535  // Allocate cache
536  GNdarray cache(cube->nx()*cube->ny(), cube->ebins());
537 
538  // Initialise all cache values with -1 (not set)
539  for (int i = 0; i < cache.size(); ++i) {
540  cache(i) = -1.0;
541  }
542 
543  // Insert cache
544  m_cache_names.push_back(name);
545  m_cache_values.push_back(cache);
546 
547  // Set index
548  index = m_cache_names.size()-1;
549 
550  } // endif: allocated new cache
551 
552  // Get reference to CTA event bin
553  const GCTAEventBin& bin = static_cast<const GCTAEventBin&>(event);
554 
555  // Get cache value
556  irf = m_cache_values[index](bin.ipix(), bin.ieng());
557 
558  // If cache value is not valid then copute IRF
559  if (irf < 0.0) {
560 
561  // Compute IRF for spatial model
562  switch (source.model()->code()) {
564  irf = irf_ptsrc(event, source, obs);
565  break;
567  irf = irf_radial(event, source, obs);
568  break;
570  irf = irf_elliptical(event, source, obs);
571  break;
573  irf = irf_diffuse(event, source, obs);
574  break;
576  irf = irf_composite(event, source, obs);
577  break;
578  default:
579  break;
580  }
581 
582  // Set cache value
583  m_cache_values[index](bin.ipix(), bin.ieng()) = irf;
584 
585  } // endif: computed IRF
586 
587  } // endif: kludge
588 
589  // ... otherwise use release 1.7 response cache
590  else {
591 
592  // If the spatial model component has free parameters, or the response
593  // cache should not be used, or the cache does not contain the requested
594  // IRF value then compute the IRF value for the spatial model.
595  if (has_free_pars ||
596  !m_use_irf_cache ||
597  !m_irf_cache.contains(name, dir, ereco, etrue, &irf)) {
598 
599  // Compute IRF for spatial model
600  switch (source.model()->code()) {
602  irf = irf_ptsrc(event, source, obs);
603  break;
605  irf = irf_radial(event, source, obs);
606  break;
608  irf = irf_elliptical(event, source, obs);
609  break;
611  irf = irf_diffuse(event, source, obs);
612  break;
614  irf = irf_composite(event, source, obs);
615  break;
616  default:
617  break;
618  }
619 
620  } // endif: computed spatial model
621 
622  // If the spatial model has no free parameters and the response cache
623  // should be used then put the IRF value in the response cache.
624  if (!has_free_pars && m_use_irf_cache) {
625  m_irf_cache.set(name, dir, ereco, etrue, irf);
626  }
627 
628  } // endelse: used release 1.7 response cache
629 
630  // Return IRF value
631  return irf;
632 }
633 
634 
635 /***********************************************************************//**
636  * @brief Read response information from XML element
637  *
638  * @param[in] xml XML element.
639  *
640  * Reads response information from an XML element. The Exposure, Psf,
641  * background cubes, and optionally the energy dispersion cube, are specified
642  * using
643  *
644  * <observation name="..." id="..." instrument="...">
645  * ...
646  * <parameter name="ExposureCube" file="..."/>
647  * <parameter name="PsfCube" file="..."/>
648  * <parameter name="EdispCube" file="..."/>
649  * <parameter name="BkgCube" file="..."/>
650  * </observation>
651  *
652  ***************************************************************************/
654 {
655  // Clear response
656  clear();
657 
658  // Get exposure cube information and load cube
659  const GXmlElement* exppar = gammalib::xml_get_par(G_READ, xml, "ExposureCube");
660  std::string expname = gammalib::xml_file_expand(xml,
661  exppar->attribute("file"));
662 
663  // Get PSF cube information and load cube
664  const GXmlElement* psfpar = gammalib::xml_get_par(G_READ, xml, "PsfCube");
665  std::string psfname = gammalib::xml_file_expand(xml,
666  psfpar->attribute("file"));
667 
668  // Get background cube information and load cube
669  const GXmlElement* bkgpar = gammalib::xml_get_par(G_READ, xml, "BkgCube");
670  std::string bkgname = gammalib::xml_file_expand(xml,
671  bkgpar->attribute("file"));
672 
673  // Load cubes
674  m_exposure.load(expname);
675  m_psf.load(psfname);
676  m_background.load(bkgname);
677 
678  // Optionally load energy dispersion cube
679  if (gammalib::xml_has_par(xml, "EdispCube")) {
680  const GXmlElement* edisppar = gammalib::xml_get_par(G_READ, xml, "EdispCube");
681  std::string edispname = gammalib::xml_file_expand(xml,
682  edisppar->attribute("file"));
683  m_edisp.load(edispname);
684  m_has_edisp = true;
685  }
686 
687  // Return
688  return;
689 }
690 
691 
692 /***********************************************************************//**
693  * @brief Write response information into XML element
694  *
695  * @param[in] xml XML element.
696  *
697  * Writes response information into an XML element. The Exposure, Psf
698  * and background cubes, and optionally the energy dispersion cube, are
699  * specified using
700  *
701  * <observation name="..." id="..." instrument="...">
702  * ...
703  * <parameter name="ExposureCube" file="..."/>
704  * <parameter name="PsfCube" file="..."/>
705  * <parameter name="EdispCube" file="..."/>
706  * <parameter name="BkgCube" file="..."/>
707  * </observation>
708  *
709  ***************************************************************************/
711 {
712  // Add exposure cube filename
713  std::string filename = gammalib::xml_file_reduce(xml, m_exposure.filename());
714  if (!(filename.empty())) {
715  GXmlElement* par = gammalib::xml_need_par(G_WRITE, xml, "ExposureCube");
716  par->attribute("file", filename);
717  }
718 
719  // Add PSF cube filename
720  filename = gammalib::xml_file_reduce(xml, m_psf.filename());
721  if (!(filename.empty())) {
722  GXmlElement* par = gammalib::xml_need_par(G_WRITE, xml, "PsfCube");
723  par->attribute("file", filename);
724  }
725 
726  // Optionally add energy dispersions cube filename
727  if (m_has_edisp) {
728  filename = gammalib::xml_file_reduce(xml, m_edisp.filename());
729  if (!(filename.empty())) {
730  GXmlElement* par = gammalib::xml_need_par(G_WRITE, xml, "EdispCube");
731  par->attribute("file", filename);
732  }
733  }
734 
735  // Add background cube filename
737  if (!(filename.empty())) {
738  GXmlElement* par = gammalib::xml_need_par(G_WRITE, xml, "BkgCube");
739  par->attribute("file", filename);
740  }
741 
742  // Return
743  return;
744 }
745 
746 
747 /***********************************************************************//**
748  * @brief Print response information
749  *
750  * @param[in] chatter Chattiness.
751  * @return String containing response information.
752  ***************************************************************************/
753 std::string GCTAResponseCube::print(const GChatter& chatter) const
754 {
755  // Initialise result string
756  std::string result;
757 
758  // Continue only if chatter is not silent
759  if (chatter != SILENT) {
760 
761  // Append header
762  result.append("=== GCTAResponseCube ===");
763 
764  // Append response information
765  result.append("\n"+gammalib::parformat("Energy dispersion"));
766  if (use_edisp()) {
767  result.append("Used");
768  }
769  else {
770  if (apply_edisp()) {
771  result.append("Not available");
772  }
773  else {
774  result.append("Not used");
775  }
776  }
777 
778  // Get reduced chatter level
779  GChatter reduced_chatter = gammalib::reduce(chatter);
780 
781  // Append detailed information
782  if (chatter >= NORMAL) {
783 
784  // Append exposure cube information
785  result.append("\n"+m_exposure.print(reduced_chatter));
786 
787  // Append point spread function information
788  result.append("\n"+m_psf.print(reduced_chatter));
789 
790  // Optionally append energy dispersion information
791  if (m_has_edisp) {
792  result.append("\n"+m_edisp.print(reduced_chatter));
793  }
794 
795  // Append background information
796  result.append("\n"+m_background.print(reduced_chatter));
797 
798  // Append cache information
799  result.append("\n"+m_irf_cache.print(reduced_chatter));
800 
801  } // endif: appended detailed information
802 
803  } // endif: chatter was not silent
804 
805  // Return result
806  return result;
807 }
808 
809 
810 /*==========================================================================
811  = =
812  = Private methods =
813  = =
814  ==========================================================================*/
815 
816 /***********************************************************************//**
817  * @brief Initialise class members
818  ***************************************************************************/
820 {
821  // Initialise members
822  m_exposure.clear();
823  m_psf.clear();
824  m_edisp.clear();
826  m_apply_edisp = false;
827  m_has_edisp = false;
828 
829  // Kludge: Initialise cube response cache
830  m_cache_names.clear();
831  m_cache_values.clear();
832 
833  // Return
834  return;
835 }
836 
837 
838 /***********************************************************************//**
839  * @brief Copy class members
840  *
841  * @param[in] rsp Response.
842  ***************************************************************************/
844 {
845  // Copy members
846  m_exposure = rsp.m_exposure;
847  m_psf = rsp.m_psf;
848  m_edisp = rsp.m_edisp;
851  m_has_edisp = rsp.m_has_edisp;
852 
853  // Kludge: Copy cube response cache
856 
857  // Return
858  return;
859 }
860 
861 
862 /***********************************************************************//**
863  * @brief Delete class members
864  ***************************************************************************/
866 {
867  // Return
868  return;
869 }
870 
871 
872 /***********************************************************************//**
873  * @brief Integrate Psf over radial model
874  *
875  * @param[in] model Radial model.
876  * @param[in] delta_mod Angle between model centre and measured photon
877  * direction (radians).
878  * @param[in] obsDir Observed event direction.
879  * @param[in] srcEng True photon energy.
880  * @param[in] srcTime True photon arrival time.
881  *
882  * Integrates the product of the spatial model and the point spread
883  * function over the true photon arrival direction using
884  *
885  * \f[
886  * \int_{\delta_{\rm min}}^{\delta_{\rm max}}
887  * \sin \delta \times {\rm Psf}(\delta) \times
888  * \int_{\phi_{\rm min}}^{\phi_{\rm max}}
889  * S_{\rm p}(\delta, \phi | E, t) d\phi d\delta
890  * \f]
891  *
892  * where
893  * \f$S_{\rm p}(\delta, \phi | E, t)\f$ is the radial spatial model,
894  * \f${\rm Psf}(\delta)\f$ is the point spread function,
895  * \f$\delta\f$ is angular distance between the true and the measured
896  * photon direction, and
897  * \f$\phi\f$ is the position angle around the observed photon direction
898  * measured counterclockwise from the connecting line between the model
899  * centre and the observed photon arrival direction.
900  ***************************************************************************/
902  const double& delta_mod,
903  const GSkyDir& obsDir,
904  const GEnergy& srcEng,
905  const GTime& srcTime) const
906 {
907  // Set number of iterations for Romberg integration.
908  // These values have been determined after careful testing, see
909  // https://cta-redmine.irap.omp.eu/issues/1291
910  static const int iter_delta = 5;
911  static const int iter_phi = 6;
912 
913  // Initialise value
914  double value = 0.0;
915 
916  // Get maximum Psf radius (radians)
917  double psf_max = psf().delta_max();
918 
919  // Get maximum model radius (radians)
920  double theta_max = model->theta_max();
921 
922  // Set offset angle integration range (radians)
923  double delta_min = (delta_mod > theta_max) ? delta_mod - theta_max : 0.0;
924  double delta_max = delta_mod + theta_max;
925  if (delta_max > psf_max) {
926  delta_max = psf_max;
927  }
928 
929  // Setup integration kernel. We take here the observed photon arrival
930  // direction as the true photon arrival direction because the PSF does
931  // not vary significantly over a small region.
932  cta_psf_radial_kern_delta integrand(this,
933  model,
934  obsDir,
935  srcEng,
936  srcTime,
937  delta_mod,
938  theta_max,
939  iter_phi);
940 
941  // Integrate over model's zenith angle
942  GIntegral integral(&integrand);
943  integral.fixed_iter(iter_delta);
944 
945  // Setup integration boundaries
946  std::vector<double> bounds;
947  bounds.push_back(delta_min);
948  bounds.push_back(delta_max);
949 
950  // If the integration range includes a transition between full
951  // containment of Psf within model and partial containment, then
952  // add a boundary at this location
953  double transition_point = theta_max - delta_mod;
954  if (transition_point > delta_min && transition_point < delta_max) {
955  bounds.push_back(transition_point);
956  }
957 
958  // Integrate kernel
959  value = integral.romberg(bounds, iter_delta);
960 
961  // Compile option: Check for NaN/Inf
962  #if defined(G_NAN_CHECK)
963  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
964  std::cout << "*** ERROR: GCTAResponseCube::psf_radial:";
965  std::cout << " NaN/Inf encountered";
966  std::cout << " (value=" << value;
967  std::cout << ", delta_min=" << delta_min;
968  std::cout << ", delta_max=" << delta_max << ")";
969  std::cout << std::endl;
970  }
971  #endif
972 
973  // Return integral
974  return value;
975 }
976 
977 
978 /***********************************************************************//**
979  * @brief Integrate Psf over elliptical model
980  *
981  * @param[in] model Elliptical model.
982  * @param[in] rho_obs Angle between model centre and measured photon direction
983  * (radians).
984  * @param[in] posangle_obs Position angle of measured photon direction with
985  * respect to model centre (radians).
986  * @param[in] obsDir Observed event direction.
987  * @param[in] srcEng True photon energy.
988  * @param[in] srcTime True photon arrival time.
989  *
990  * Integrates the product of the spatial model and the point spread
991  * function over the true photon arrival direction using
992  *
993  * \f[
994  * \int_{\rho_{\rm min}}^{\rho_{\rm max}}
995  * \sin \rho \times
996  * \int_{\omega}
997  * S_{\rm p}(\rho, \omega | E, t) \times
998  * {\rm Psf}(\rho, \omega) d\omega d\rho
999  * \f]
1000  *
1001  * where
1002  * \f$S_{\rm p}(\rho, \omega | E, t)\f$ is the elliptical spatial model,
1003  * \f${\rm Psf}(\rho, \omega)\f$ is the point spread function,
1004  * \f$\rho\f$ is the radial distance from the model centre, and
1005  * \f$\omega\f$ is the position angle around the model centre measured
1006  * counterclockwise from the connecting line between the model centre and
1007  * the observed photon arrival direction.
1008  ***************************************************************************/
1010  const double& rho_obs,
1011  const double& posangle_obs,
1012  const GSkyDir& obsDir,
1013  const GEnergy& srcEng,
1014  const GTime& srcTime) const
1015 {
1016  // Set number of iterations for Romberg integration.
1017  // These values have been determined after careful testing, see
1018  // https://cta-redmine.irap.omp.eu/issues/1299
1019  static const int iter_rho = 5;
1020  static const int iter_phi = 5;
1021 
1022  // Initialise value
1023  double irf = 0.0;
1024 
1025  // Get maximum PSF radius (radians)
1026  double delta_max = psf().delta_max();
1027 
1028  // Get the ellipse boundary (radians). Note that these are NOT the
1029  // parameters of the ellipse but of a boundary ellipse that is used
1030  // for computing the relevant omega angle intervals for a given angle
1031  // rho. The boundary ellipse takes care of the possibility that the
1032  // semiminor axis is larger than the semimajor axis
1033  double semimajor; // Will be the larger axis
1034  double semiminor; // Will be the smaller axis
1035  double posangle; // Will be the corrected position angle
1036  double aspect_ratio; // Ratio between smaller/larger axis of model
1037  if (model->semimajor() >= model->semiminor()) {
1038  aspect_ratio = (model->semimajor() > 0.0) ?
1039  model->semiminor() / model->semimajor() : 0.0;
1040  posangle = model->posangle() * gammalib::deg2rad;
1041  }
1042  else {
1043  aspect_ratio = (model->semiminor() > 0.0) ?
1044  model->semimajor() / model->semiminor() : 0.0;
1045  posangle = model->posangle() * gammalib::deg2rad + gammalib::pihalf;
1046  }
1047  semimajor = model->theta_max();
1048  semiminor = semimajor * aspect_ratio;
1049 
1050  // Set zenith angle integration range for elliptical model
1051  double rho_min = (rho_obs > delta_max) ? rho_obs - delta_max : 0.0;
1052  double rho_max = rho_obs + delta_max;
1053  if (rho_max > semimajor) {
1054  rho_max = semimajor;
1055  }
1056 
1057  // Perform zenith angle integration if interval is valid
1058  if (rho_max > rho_min) {
1059 
1060  // Setup integration kernel. We take here the observed photon arrival
1061  // direction as the true photon arrival direction because the PSF does
1062  // not vary significantly over a small region.
1063  cta_psf_elliptical_kern_rho integrand(this,
1064  model,
1065  semimajor,
1066  semiminor,
1067  posangle,
1068  obsDir,
1069  srcEng,
1070  srcTime,
1071  rho_obs,
1072  posangle_obs,
1073  delta_max,
1074  iter_phi);
1075 
1076  // Integrate over model's zenith angle
1077  GIntegral integral(&integrand);
1078  integral.fixed_iter(iter_rho);
1079 
1080  // Setup integration boundaries
1081  std::vector<double> bounds;
1082  bounds.push_back(rho_min);
1083  bounds.push_back(rho_max);
1084 
1085  // Kluge: add this transition point as this allows to fit the test
1086  // case without any stalls. Not clear why this is the case, maybe
1087  // simply because the rho integral gets cut down into one more
1088  // sub-interval which may increase precision and smoothed the
1089  // likelihood contour
1090  double transition_point = delta_max - rho_obs;
1091  if (transition_point > rho_min && transition_point < rho_max) {
1092  bounds.push_back(transition_point);
1093  }
1094 
1095  // If the integration range includes the semiminor boundary, then
1096  // add an integration boundary at that location
1097  if (semiminor > rho_min && semiminor < rho_max) {
1098  bounds.push_back(semiminor);
1099  }
1100 
1101  // Integrate kernel
1102  irf = integral.romberg(bounds, iter_rho);
1103 
1104  // Compile option: Check for NaN/Inf
1105  #if defined(G_NAN_CHECK)
1107  std::cout << "*** ERROR: GCTAResponseCube::psf_elliptical:";
1108  std::cout << " NaN/Inf encountered";
1109  std::cout << " (irf=" << irf;
1110  std::cout << ", rho_min=" << rho_min;
1111  std::cout << ", rho_max=" << rho_max << ")";
1112  std::cout << std::endl;
1113  }
1114  #endif
1115 
1116  } // endif: integration interval is valid
1117 
1118  // Return PSF
1119  return irf;
1120 }
1121 
1122 
1123 /***********************************************************************//**
1124  * @brief Integrate PSF over diffuse model
1125  *
1126  * @param[in] model Diffuse spatial model.
1127  * @param[in] obsDir Observed event direction.
1128  * @param[in] srcEng True photon energy.
1129  * @param[in] srcTime True photon arrival time.
1130  *
1131  * Computes the integral
1132  *
1133  * \f[
1134  * \int_0^{\delta_{\rm max}}
1135  * {\rm PSF}(\delta) \times
1136  * \int_0^{2\pi} {\rm Map}(\delta, \phi) \sin \delta
1137  * {\rm d}\phi {\rm d}\delta
1138  * \f]
1139  *
1140  * where \f${\rm Map}(\delta, \phi)\f$ is the diffuse map in the coordinate
1141  * system of the point spread function, defined by the angle \f$\delta\f$
1142  * between the true and the measured photon direction and the azimuth angle
1143  * \f$\phi\f$ around the measured photon direction.
1144  * \f${\rm PSF}(\delta)\f$ is the azimuthally symmetric point spread
1145  * function.
1146  ***************************************************************************/
1148  const GSkyDir& obsDir,
1149  const GEnergy& srcEng,
1150  const GTime& srcTime) const
1151 {
1152  // Set minimum and maximum number of iterations for Romberg integration.
1153  // These values have been determined after careful testing, see
1154  // https://cta-redmine.irap.omp.eu/issues/2458
1155  static const int min_iter_delta = 5;
1156  static const int min_iter_phi = 5;
1157  static const int max_iter_delta = 8;
1158  static const int max_iter_phi = 8;
1159 
1160  // Compute rotation matrix to convert from PSF centred coordinate system
1161  // spanned by delta and phi into the reference frame of the observed
1162  // arrival direction given in Right Ascension and Declination.
1163  GMatrix ry;
1164  GMatrix rz;
1165  ry.eulery(obsDir.dec_deg() - 90.0);
1166  rz.eulerz(-obsDir.ra_deg());
1167  GMatrix rot = (ry * rz).transpose();
1168 
1169  // Get offset angle integration interval in radians
1170  double delta_min = 0.0;
1171  double delta_max = 1.1 * this->psf().delta_max();
1172 
1173  // Get resolution of spatial model
1174  double resolution = gammalib::resolution(model);
1175 
1176 
1177  // Setup integration kernel. We take here the observed photon arrival
1178  // direction as the true photon arrival direction because the PSF does
1179  // not vary significantly over a small region.
1180  cta_psf_diffuse_kern_delta integrand(this, model, &rot,
1181  obsDir, srcEng, srcTime,
1182  min_iter_phi, max_iter_phi,
1183  resolution);
1184 
1185  // Set number of radial integration iterations
1186  int iter = gammalib::iter_rho(delta_max, resolution,
1187  min_iter_delta, max_iter_delta);
1188 
1189  // Setup integration
1190  GIntegral integral(&integrand);
1191 
1192  // Set fixed number of iterations
1193  integral.fixed_iter(iter);
1194 
1195  // Integrate over PSF delta angle
1196  double psf = integral.romberg(delta_min, delta_max);
1197 
1198  // Return PSF
1199  return psf;
1200 }
1201 
1202 
1203 /***********************************************************************//**
1204  * @brief Return instrument response to point source
1205  *
1206  * @param[in] event Observed event.
1207  * @param[in] source Source.
1208  * @param[in] obs Observation (not used).
1209  * @return Instrument response to point source.
1210  *
1211  * Returns the instrument response to a specified point source.
1212  ***************************************************************************/
1214  const GSource& source,
1215  const GObservation& obs) const
1216 {
1217  // Initialise IRF
1218  double irf = 0.0;
1219 
1220  // Get pointer to model source model
1221  const GModelSpatialPointSource* ptsrc =
1222  static_cast<const GModelSpatialPointSource*>(source.model());
1223 
1224  // Get point source direction
1225  GSkyDir srcDir = ptsrc->dir();
1226 
1227  // Get energy of source model
1228  GEnergy srcEng = source.energy();
1229 
1230  // Get pointer on CTA event bin
1231  if (!event.is_bin()) {
1232  std::string msg = "The current event is not a CTA event bin. "
1233  "This method only works on binned CTA data. Please "
1234  "make sure that a CTA observation containing binned "
1235  "CTA data is provided.";
1237  }
1238  const GCTAEventBin* bin = static_cast<const GCTAEventBin*>(&event);
1239 
1240  // Determine angular separation between true and measured photon
1241  // direction in radians
1242  double delta = bin->dir().dir().dist(srcDir);
1243 
1244  // Get maximum angular separation for PSF (in radians)
1245  double delta_max = psf().delta_max();
1246 
1247  // Get livetime (in seconds)
1248  double livetime = exposure().livetime();
1249 
1250  // Continue only if livetime is >0 and if we're sufficiently close
1251  // to the PSF
1252  if ((livetime > 0.0) && (delta <= delta_max)) {
1253 
1254  // Get exposure
1255  irf = exposure()(srcDir, srcEng);
1256 
1257  // Multiply-in PSF
1258  if (irf > 0.0) {
1259 
1260  // Recover effective area from exposure
1261  irf /= livetime;
1262 
1263  // Get PSF component
1264  irf *= psf()(srcDir, delta, srcEng);
1265 
1266  // Multiply-in energy dispersion
1267  if (use_edisp() && irf > 0.0) {
1268  irf *= edisp()(bin->energy(), srcEng, srcDir);
1269  }
1270 
1271  // Apply deadtime correction
1272  irf *= exposure().deadc();
1273 
1274  } // endif: exposure was non-zero
1275 
1276  } // endif: we were sufficiently close to PSF and livetime >0
1277 
1278  // Compile option: Check for NaN/Inf
1279  #if defined(G_NAN_CHECK)
1281  std::cout << "*** ERROR: GCTAResponseCube::irf_ptsrc:";
1282  std::cout << " NaN/Inf encountered";
1283  std::cout << " irf=" << irf;
1284  std::cout << std::endl;
1285  }
1286  #endif
1287 
1288  // Return IRF value
1289  return irf;
1290 }
1291 
1292 
1293 /***********************************************************************//**
1294  * @brief Return instrument response to radial source
1295  *
1296  * @param[in] event Observed event.
1297  * @param[in] source Source.
1298  * @param[in] obs Observation (not used).
1299  * @return Instrument response to radial source.
1300  *
1301  * Returns the instrument response to a specified radial source.
1302  *
1303  * @todo Correct assumptions in exposure computation, PSF computation and
1304  * energy dispersion computation.
1305  ***************************************************************************/
1307  const GSource& source,
1308  const GObservation& obs) const
1309 {
1310  // Initialise IRF
1311  double irf = 0.0;
1312 
1313  // Get pointer to CTA event bin
1314  if (!event.is_bin()) {
1315  std::string msg = "The current event is not a CTA event bin. "
1316  "This method only works on binned CTA data. Please "
1317  "make sure that a CTA observation containing binned "
1318  "CTA data is provided.";
1320  }
1321  const GCTAEventBin* bin = static_cast<const GCTAEventBin*>(&event);
1322 
1323  // Get event attribute references
1324  const GSkyDir& obsDir = bin->dir().dir();
1325  const GEnergy& obsEng = bin->energy();
1326  const GTime& obsTime = bin->time();
1327 
1328  // Get energy of source model
1329  GEnergy srcEng = source.energy();
1330 
1331  // Get pointer to radial model
1332  const GModelSpatialRadial* model =
1333  static_cast<const GModelSpatialRadial*>(source.model());
1334 
1335  // Compute angle between model centre and measured photon direction
1336  // (radians)
1337  double rho_obs = model->dir().dist(obsDir);
1338 
1339  // Get livetime (in seconds)
1340  double livetime = exposure().livetime();
1341 
1342  // Continue only if livetime is >0 and if we're sufficiently close to
1343  // the model centre to get a non-zero response
1344  if ((livetime > 0.0) && (rho_obs <= model->theta_max()+psf().delta_max())) {
1345 
1346  // Get exposure at the observed event direction.
1347  //
1348  // The current code assumes that the exposure at the observed and
1349  // true event direction does not vary significantly. In other words,
1350  // the code assumes that the exposure is constant over the size of
1351  // the PSF.
1352  irf = exposure()(obsDir, srcEng);
1353 
1354  // Continue only if exposure is positive
1355  if (irf > 0.0) {
1356 
1357  // Recover effective area from exposure
1358  irf /= livetime;
1359 
1360  // Get PSF component
1361  irf *= psf_radial(model, rho_obs, obsDir, srcEng, obsTime);
1362 
1363  // Multiply-in energy dispersion
1364  //
1365  // The current code assumes that the energy dispersion at the
1366  // observed and true event direction does not vary
1367  // significantly. In other words, the code assumes that the
1368  // energy dispersion is constant over the size of the PSF.
1369  if (use_edisp() && irf > 0.0) {
1370  irf *= edisp()(bin->energy(), srcEng, obsDir);
1371  }
1372 
1373  // Apply deadtime correction
1374  irf *= exposure().deadc();
1375 
1376  } // endif: exposure was positive
1377 
1378  } // endif: we were sufficiently close and livetime >0
1379 
1380  // Compile option: Check for NaN/Inf
1381  #if defined(G_NAN_CHECK)
1383  std::cout << "*** ERROR: GCTAResponseCube::irf_radial:";
1384  std::cout << " NaN/Inf encountered";
1385  std::cout << " irf=" << irf;
1386  std::cout << std::endl;
1387  }
1388  #endif
1389 
1390  // Return IRF value
1391  return irf;
1392 }
1393 
1394 
1395 /***********************************************************************//**
1396  * @brief Return instrument response to elliptical source
1397  *
1398  * @param[in] event Observed event.
1399  * @param[in] source Source.
1400  * @param[in] obs Observation (not used).
1401  * @return Instrument response to elliptical source.
1402  *
1403  * Returns the instrument response to a specified elliptical source.
1404  *
1405  * @todo Correct assumptions in exposure computation, PSF computation and
1406  * energy dispersion computation.
1407  ***************************************************************************/
1409  const GSource& source,
1410  const GObservation& obs) const
1411 {
1412  // Initialise IRF
1413  double irf = 0.0;
1414 
1415  // Get pointer to CTA event bin
1416  if (!event.is_bin()) {
1417  std::string msg = "The current event is not a CTA event bin. "
1418  "This method only works on binned CTA data. Please "
1419  "make sure that a CTA observation containing binned "
1420  "CTA data is provided.";
1422  }
1423  const GCTAEventBin* bin = static_cast<const GCTAEventBin*>(&event);
1424 
1425  // Get event attribute references
1426  const GSkyDir& obsDir = bin->dir().dir();
1427  const GEnergy& obsEng = bin->energy();
1428  const GTime& obsTime = bin->time();
1429 
1430  // Get energy of source model
1431  GEnergy srcEng = source.energy();
1432 
1433  // Get pointer to elliptical model
1434  const GModelSpatialElliptical* model =
1435  static_cast<const GModelSpatialElliptical*>(source.model());
1436 
1437  // Compute angle between model centre and measured photon direction and
1438  // position angle (radians)
1439  double rho_obs = model->dir().dist(obsDir);
1440  double posangle_obs = model->dir().posang(obsDir); // Celestial
1441 
1442  // Get livetime (in seconds)
1443  double livetime = exposure().livetime();
1444 
1445  // Continue only if livetime is >0 and if we're sufficiently close to
1446  // the model centre to get a non-zero response
1447  if ((livetime > 0.0) && (rho_obs <= model->theta_max()+psf().delta_max())) {
1448 
1449  // Get exposure
1450  //
1451  // The current code assumes that the exposure at the observed and
1452  // true event direction does not vary significantly. In other words,
1453  // the code assumes that the exposure is constant over the size of
1454  // the PSF.
1455  irf = exposure()(obsDir, srcEng);
1456 
1457  // Continue only if exposure is positive
1458  if (irf > 0.0) {
1459 
1460  // Recover effective area from exposure
1461  irf /= livetime;
1462 
1463  // Get PSF component
1464  irf *= psf_elliptical(model, rho_obs, posangle_obs, obsDir, srcEng, obsTime);
1465 
1466  // Multiply-in energy dispersion
1467  //
1468  // The current code assumes that the energy dispersion at the
1469  // observed and true event direction does not vary
1470  // significantly. In other words, the code assumes that the
1471  // energy dispersion is constant over the size of the PSF.
1472  if (use_edisp() && irf > 0.0) {
1473  irf *= edisp()(bin->energy(), srcEng, obsDir);
1474  }
1475 
1476  // Apply deadtime correction
1477  irf *= exposure().deadc();
1478 
1479  } // endif: exposure was positive
1480 
1481  } // endif: we were sufficiently close and livetime >0
1482 
1483  // Compile option: Check for NaN/Inf
1484  #if defined(G_NAN_CHECK)
1486  std::cout << "*** ERROR: GCTAResponseCube::irf_elliptical:";
1487  std::cout << " NaN/Inf encountered";
1488  std::cout << " irf=" << irf;
1489  std::cout << std::endl;
1490  }
1491  #endif
1492 
1493  // Return IRF value
1494  return irf;
1495 }
1496 
1497 
1498 /***********************************************************************//**
1499  * @brief Return instrument response to diffuse source
1500  *
1501  * @param[in] event Observed event.
1502  * @param[in] source Source.
1503  * @param[in] obs Observation.
1504  * @return Instrument response to diffuse source.
1505  *
1506  * Returns the instrument response to a specified diffuse source.
1507  *
1508  * The method uses a pre-computation cache to store the instrument response
1509  * for the spatial model component. The pre-computation cache is initialised
1510  * if no cache has yet been allocated, or if at the beginning of a scan over
1511  * the events, the model parameters have changed. The beginning of a scan is
1512  * defined by an event bin index of 0.
1513  ***************************************************************************/
1515  const GSource& source,
1516  const GObservation& obs) const
1517 {
1518  // Initialise IRF
1519  double irf = 0.0;
1520 
1521  // Get pointer to CTA event bin
1522  if (!event.is_bin()) {
1523  std::string msg = "The current event is not a CTA event bin. "
1524  "This method only works on binned CTA data. Please "
1525  "make sure that a CTA observation containing binned "
1526  "CTA data is provided.";
1528  }
1529  const GCTAEventBin* bin = static_cast<const GCTAEventBin*>(&event);
1530 
1531  // Get event attribute references
1532  const GSkyDir& obsDir = bin->dir().dir();
1533  const GEnergy& obsEng = bin->energy();
1534  const GTime& obsTime = bin->time();
1535 
1536  // Get energy of source model
1537  GEnergy srcEng = source.energy();
1538 
1539  // Get pointer to elliptical model
1540  const GModelSpatialDiffuse* model =
1541  static_cast<const GModelSpatialDiffuse*>(source.model());
1542 
1543  // Get livetime (in seconds) and deadtime correction factor
1544  double livetime = exposure().livetime();
1545  double deadc = exposure().deadc();
1546 
1547  // Get Psf radius (in degrees)
1548  double delta_max = psf().delta_max() * gammalib::rad2deg;
1549 
1550  // Continue only if livetime is >0 and model contains reconstructed
1551  // sky direction
1552  if (livetime > 0.0 && model->contains(obsDir, delta_max)) {
1553 
1554  // Get exposure
1555  //
1556  // The current code assumes that the exposure at the observed and
1557  // true event direction does not vary significantly. In other words,
1558  // the code assumes that the exposure is constant over the size of
1559  // the PSF.
1560  irf = exposure()(obsDir, srcEng);
1561 
1562  // Continue only if exposure is positive
1563  if (irf > 0.0) {
1564 
1565  // Recover effective area from exposure
1566  irf /= livetime;
1567 
1568  // Compute product of PSF and diffuse map, integrated over the
1569  // relevant PSF area. We assume no energy dispersion and thus
1570  // compute the product using the observed energy.
1571  irf *= psf_diffuse(model, obsDir, srcEng, obsTime);
1572 
1573  // Multiply-in energy dispersion
1574  //
1575  // The current code assumes that the energy dispersion at the
1576  // observed and true event direction does not vary
1577  // significantly. In other words, the code assumes that the
1578  // energy dispersion is constant over the size of the PSF.
1579  if (use_edisp() && irf > 0.0) {
1580  irf *= edisp()(bin->energy(), srcEng, obsDir);
1581  }
1582 
1583  // Apply deadtime correction
1584  irf *= deadc;
1585 
1586  } // endif: exposure was positive
1587 
1588  } // endif: livetime was positive
1589 
1590  // Compile option: Check for NaN/Inf
1591  #if defined(G_NAN_CHECK)
1593  std::cout << "*** ERROR: GCTAResponseCube::irf_diffuse:";
1594  std::cout << " NaN/Inf encountered";
1595  std::cout << " irf=" << irf;
1596  std::cout << std::endl;
1597  }
1598  #endif
1599 
1600  // Return IRF value
1601  return irf;
1602 }
1603 
1604 
1605 /*==========================================================================
1606  = =
1607  = Vector response methods =
1608  = =
1609  ==========================================================================*/
1610 
1611 /***********************************************************************//**
1612  * @brief Return instrument response to radial source
1613  *
1614  * @param[in] model Sky model.
1615  * @param[in] obs Observation.
1616  * @param[out] gradients Optional spatial model gradients for all events.
1617  * @return Instrument response to radial source for all events in observation.
1618  *
1619  * Returns the instrument response to a specified radial source.
1620  *
1621  * @p gradients is an optional matrix where the number of rows corresponds
1622  * to the number of events in the observation and the number of columns
1623  * corresponds to the number of spatial model parameters.
1624  ***************************************************************************/
1626  const GObservation& obs,
1627  GMatrix* gradients) const
1628 {
1629  // Get number of events
1630  int nevents = obs.events()->size();
1631 
1632  // Initialise result
1633  GVector irfs(nevents);
1634 
1635  // Get livetime (in seconds)
1636  double livetime = exposure().livetime();
1637 
1638  // Continue only if livetime is positive
1639  if (livetime > 0.0) {
1640 
1641  // Set gradients flag
1642  bool grad = (gradients != NULL);
1643 
1644  // Get pointer to radial model
1645  const GModelSpatialRadial* radial =
1646  static_cast<const GModelSpatialRadial*>(model.spatial());
1647 
1648  // Get CTA event cube
1650 
1651  // Get number of radial model parameters and CTA event cube dimensions
1652  int npars = radial->size();
1653  int ndirs = cube.npix();
1654  int nengs = cube.ebins();
1655 
1656  // Check matrix consistency
1657  if (grad) {
1658  if (gradients->rows() != nevents) {
1659  std::string msg = "Number of "+gammalib::str(gradients->rows())+
1660  " rows in gradient matrix differs from number "
1661  "of "+gammalib::str(nevents)+" events in "
1662  "observation. Please provide a compatible "
1663  "gradient matrix.";
1665  }
1666  if (gradients->columns() != npars) {
1667  std::string msg = "Number of "+gammalib::str(gradients->columns())+
1668  " columns in gradient matrix differs from number "
1669  "of "+gammalib::str(npars)+" parameters in "
1670  "model. Please provide a compatible "
1671  "gradient matrix.";
1673  }
1674  }
1675 
1676  // Setup energies container
1677  GEnergies srcEngs;
1678  for (int ieng = 0; ieng < nengs; ++ieng) {
1679  srcEngs.append(cube.energy(ieng));
1680  }
1681 
1682  // If requested, setup vectors of gradients
1683  GVector* gradient = NULL;
1684  if (grad) {
1685 
1686  // Allocate vectors to hold the gradients
1687  gradient = new GVector[npars];
1688  for (int ipar = 0; ipar < npars; ++ipar) {
1689  gradient[ipar] = GVector(nevents);
1690  }
1691 
1692  // Signal all parameters that will have analytical gradients.
1693  // These are all parameters that are free for which the model
1694  // provides analytical gradients.
1695  for (int ipar = 0; ipar < npars; ++ipar) {
1696  const GModelPar& par = (*radial)[ipar];
1697  if (par.is_free() && par.has_grad()) {
1698  obs.computed_gradient(model, par);
1699  }
1700  }
1701 
1702  } // endif: setup of gradient computation
1703 
1704  // Setup cube time
1705  GTime srcTime = cube.time();
1706 
1707  // Set IRF normalisation
1708  double norm = exposure().deadc() / livetime;
1709 
1710  // Loop over event directions
1711  for (int idir = 0; idir < ndirs; ++idir) {
1712 
1713  // Get event direction
1714  GSkyDir obsDir = cube.counts().inx2dir(idir);
1715 
1716  // Compute angle between model centre and measured photon direction
1717  // (radians)
1718  double zeta = radial->dir().dist(obsDir);
1719 
1720  // Continue only if we're sufficiently close to the model centre to
1721  // get a non-zero response
1722  if (zeta <= radial->theta_max()+psf().delta_max()) {
1723 
1724  // Get radially integrated PSF
1725  GVector psf = psf_radial(radial, zeta, obsDir, srcEngs, srcTime, grad);
1726 
1727  // Loop over true energies
1728  for (int ieng = 0, index = idir; ieng < nengs; ++ieng, index += ndirs) {
1729 
1730  // Get effective area at the observed direction. This
1731  // assumes that the exposure at the observed and true event
1732  // direction does not vary significantly. In other words,
1733  // is is assumed that the exposure is constant over the
1734  // size of the PSF.
1735  double aeff = norm * exposure()(obsDir, srcEngs[ieng]);
1736 
1737  // Compute IRF value
1738  irfs[index] = aeff * psf[ieng];
1739 
1740  // Optionally compute gradients
1741  if (grad) {
1742  for (int ipar = 0, ipsf = ieng+nengs; ipar < npars;
1743  ++ipar, ipsf += nengs) {
1744  gradient[ipar][index] = aeff * psf[ipsf];
1745  }
1746  }
1747 
1748  } // endfor: looped over energies
1749 
1750  } // endif: direction close enough to model centre
1751 
1752  } // endfor: looped over event directions
1753 
1754  // If gradients were requested then insert vectors into the gradient
1755  // matrix
1756  if (grad) {
1757  for (int ipar = 0; ipar < npars; ++ipar) {
1758  gradients->column(ipar, gradient[ipar]);
1759  }
1760  }
1761 
1762  // If needed, free vectors of gradients
1763  if (gradient != NULL) {
1764  delete [] gradient;
1765  }
1766 
1767  } // endif: livetime was positive
1768 
1769  // Return IRF value
1770  return irfs;
1771 }
1772 
1773 
1774 /***********************************************************************//**
1775  * @brief Integrate Psf over radial model
1776  *
1777  * @param[in] model Radial model.
1778  * @param[in] zeta Angle between model centre and event direction (radians).
1779  * @param[in] obsDir Observed event direction.
1780  * @param[in] srcEngs True photon energies.
1781  * @param[in] srcTime True photon arrival time.
1782  * @param[in] grad Compute gradients?
1783  ***************************************************************************/
1785  const double& zeta,
1786  const GSkyDir& obsDir,
1787  const GEnergies& srcEngs,
1788  const GTime& srcTime,
1789  const bool& grad) const
1790 {
1791  // Set number of iterations for Romberg integration.
1792  // These values have been determined after careful testing, see
1793  // https://cta-redmine.irap.omp.eu/issues/1291
1794  static const int iter_delta = 5;
1795  static const int iter_phi = 6;
1796 
1797  // Get maximum Psf radius (radians)
1798  double psf_max = psf().delta_max();
1799 
1800  // Get maximum model radius (radians)
1801  double theta_max = model->theta_max();
1802 
1803  // Set offset angle integration range (radians)
1804  double delta_min = (zeta > theta_max) ? zeta - theta_max : 0.0;
1805  double delta_max = zeta + theta_max;
1806  if (delta_max > psf_max) {
1807  delta_max = psf_max;
1808  }
1809 
1810  // Setup integration kernel. We take here the observed photon arrival
1811  // direction as the true photon arrival direction because the PSF does
1812  // not vary significantly over a small region.
1813  cta_psf_radial_kerns_delta integrand(this,
1814  model,
1815  obsDir,
1816  srcEngs,
1817  zeta,
1818  theta_max,
1819  iter_phi,
1820  grad);
1821 
1822  // Integrate over model's zenith angle
1823  GIntegrals integral(&integrand);
1824  integral.fixed_iter(iter_delta);
1825 
1826  // Setup integration boundaries
1827  std::vector<double> bounds;
1828  bounds.push_back(delta_min);
1829  bounds.push_back(delta_max);
1830 
1831  // If the integration range includes a transition between full
1832  // containment of Psf within model and partial containment, then
1833  // add a boundary at this location
1834  double transition_point = theta_max - zeta;
1835  if (transition_point > delta_min && transition_point < delta_max) {
1836  bounds.push_back(transition_point);
1837  }
1838 
1839  // Integrate kernel
1840  GVector values = integral.romberg(bounds, iter_delta);
1841 
1842  // Compile option: Check for NaN/Inf
1843  #if defined(G_NAN_CHECK)
1844  for (int i = 0; i < srcEngs.size(); ++i) {
1845  if (gammalib::is_notanumber(values[i]) ||
1846  gammalib::is_infinite(values[i])) {
1847  std::cout << "*** ERROR: GCTAResponseCube::psf_radial:";
1848  std::cout << " NaN/Inf encountered";
1849  std::cout << " (value=" << values[i];
1850  std::cout << ", delta_min=" << delta_min;
1851  std::cout << ", delta_max=" << delta_max << ")";
1852  std::cout << std::endl;
1853  }
1854  }
1855  #endif
1856 
1857  // Return integrals
1858  return values;
1859 }
const GFilename & filename(void) const
Return edisp cube filename.
#define G_IRF_DIFFUSE
void copy_members(const GCTAResponseCube &rsp)
Copy class members.
const int & ieng(void) const
Return the energy layer index.
std::string print(const GChatter &chatter=NORMAL) const
Print energy dispersion cube information.
CTA response helper classes definition.
CTA exposure cube class.
bool m_use_irf_cache
Control usage of irf cache.
Definition: GResponse.hpp:313
double irf_radial(const GEvent &event, const GSource &source, const GObservation &obs) const
Return instrument response to radial source.
void computed_gradient(const GModel &model, const GModelPar &par) const
Signals that an analytical gradient was computed for a model parameter.
double norm(const GVector &vector)
Computes vector norm.
Definition: GVector.cpp:864
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:382
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.
const GSkyDir & dir(void) const
Return position of point source.
int size(void) const
Return number of parameters.
void clear(void)
Clear instance.
const GTime & time(void) const
Return event cube mean time.
bool has_free_pars(void) const
Checks if the spatial model has free parameters.
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.
#define G_IRF_RADIAL2
void clear(void)
Clear instance.
Abstract interface for the event classes.
Definition: GEvent.hpp:71
void fixed_iter(const int &iter)
Set fixed number of iterations.
Definition: GIntegrals.hpp:179
XML element node class.
Definition: GXmlElement.hpp:48
std::vector< std::string > m_cache_names
Model names.
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
const GCTAEventCube & cta_event_cube(const std::string &origin, const GObservation &obs)
Retrieve CTA event cube from generic observation.
const GSkyDir & dir(void) const
Return position of elliptical spatial model.
virtual double irf_composite(const GEvent &event, const GSource &source, const GObservation &obs) const
Return instrument response to composite source.
Definition: GResponse.cpp:1255
void counts(const GSkyMap &counts)
Set event cube counts from sky map.
virtual GVector column(const int &column) const
Extract column as vector from matrix.
Definition: GMatrix.cpp:718
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:55
Gammalib tools definition.
virtual double irf_spatial(const GEvent &event, const GSource &source, const GObservation &obs) const
Return instrument response integrated over the spatial model.
GIntegral class interface definition.
Definition: GIntegral.hpp:46
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.
const std::string & name(void) const
Return model name.
Definition: GSource.hpp:114
void init_members(void)
Initialise class members.
Definition: GResponse.cpp:787
GFilename xml_file_reduce(const GXmlElement &xml, const std::string &filename)
Reduce file name provided for writing as XML attribute.
Definition: GTools.cpp:1946
Source class definition.
virtual GClassCode code(void) const =0
void dir(const GSkyDir &dir)
Set sky direction.
void init_members(void)
Initialise class members.
bool is_free(void) const
Signal if parameter is free.
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:1157
virtual double theta_max(void) const =0
Energy container class.
Definition: GEnergies.hpp:60
std::string print(const GChatter &chatter=NORMAL) const
Print exposure cube information.
void id(const std::string &id)
Set observation identifier.
int ny(void) const
Return number of bins in Y direction.
Abstract instrument direction base class.
Definition: GInstDir.hpp:51
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:201
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:184
virtual GCTAResponse & operator=(const GCTAResponse &rsp)
Assignment operator.
Definition of support function used by CTA classes.
Model parameter class.
Definition: GModelPar.hpp:87
std::vector< GNdarray > m_cache_values
Cached values.
Kernel for radial spatial model PSF delta angle integration.
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
CTA event bin container class.
void clear(void)
Clear background.
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
int npix(void) const
Return number of pixels in one energy bins of the event cube.
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.
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.
const int & rows(void) const
Return number of matrix rows.
N-dimensional array class interface definition.
Energy container class definition.
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:1194
Integration class for set of functions interface definition.
virtual const GTime & time(void) const
Return time of event bin.
int nx(void) const
Return number of bins in X direction.
int size(void) const
Return number of energies in container.
Definition: GEnergies.hpp:170
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:1637
Implementation of CTA helper classes for stacked vector response.
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:229
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.
const GEnergy & energy(const int &index) const
Return energy of cube layer.
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
N-dimensional array class.
Definition: GNdarray.hpp:44
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:1596
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:835
Abstract elliptical spatial model base class interface definition.
Abstract diffuse spatial model base class interface definition.
double posang(const GSkyDir &dir, const std::string &coordsys="CEL") const
Compute position angle between sky directions in radians.
Definition: GSkyDir.cpp:1151
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 for cube analysis.
Definition: GCTACubePsf.hpp:62
Sky model class.
Definition: GModelSky.hpp:122
std::string print(const GChatter &chatter=NORMAL) const
Print PSF cube information.
double posangle(void) const
Return Position Angle of model.
#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
Sparse matrix class definition.
double delta_max(void) const
Return maximum delta value in radians.
CTA event bin container class interface definition.
const double & livetime(void) const
Return livetime.
virtual GEvents * events(void)
Return events.
virtual void read(const GXmlElement &xml)
Read response information from XML element.
const GEnergy & energy(void) const
Return photon energy.
Definition: GPhoton.hpp:122
GEnergy & append(const GEnergy &energy)
Append energy to container.
Definition: GEnergies.cpp:305
virtual double irf(const GEvent &event, const GPhoton &photon, const GObservation &obs) const
Return instrument response.
GVector romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegrals.cpp:210
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:189
GModelSpatial * spatial(void) const
Return spatial model component.
Definition: GModelSky.hpp:295
CTA instrument direction class.
Definition: GCTAInstDir.hpp:63
Abstract spatial model base class.
GResponseCache m_irf_cache
Definition: GResponse.hpp:324
Abstract radial spatial model base class.
GCTAResponseCube(void)
Void constructor.
Generic matrix class definition.
Definition: GMatrix.hpp:79
Vector class.
Definition: GVector.hpp:46
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.
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.
Integration class for set of functions.
Definition: GIntegrals.hpp:47
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:1143
void clear(void)
Clear energy dispersion cube.
Class that handles gamma-ray sources.
Definition: GSource.hpp:53
void set(const std::string &name, const GEnergy &ereco, const GEnergy &etrue, const double &value)
Set cache value.
int ebins(void) const
Return number of energy bins in the event cube.
Integration class interface definition.
const GFilename & filename(void) const
Return exposure cube filename.
const int & ipix(void) const
Return the spatial pixel index.
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
const int & columns(void) const
Return number of matrix columns.
#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:1689
std::string print(const GChatter &chatter=NORMAL) const
Print response cache.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
virtual int size(void) const =0
GFilename xml_file_expand(const GXmlElement &xml, const std::string &filename)
Expand file name provided as XML attribute for loading.
Definition: GTools.cpp:1889
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489