GammaLib  2.0.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
GResponse.cpp
Go to the documentation of this file.
1 /***************************************************************************
2  * GResponse.cpp - Abstract response base class *
3  * ----------------------------------------------------------------------- *
4  * copyright (C) 2008-2021 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 GResponse.hpp
23  * @brief Abstract response base class interface definition
24  * @author Juergen Knoedlseder
25  */
26 
27 /* __ Includes ___________________________________________________________ */
28 #ifdef HAVE_CONFIG_H
29 #include <config.h>
30 #endif
31 #include <string>
32 #include <unistd.h> // access() function
33 #include "GException.hpp"
34 #include "GTools.hpp"
35 #include "GMath.hpp"
36 #include "GVector.hpp"
37 #include "GMatrix.hpp"
38 #include "GMatrixSparse.hpp"
39 #include "GIntegral.hpp"
40 #include "GIntegrals.hpp"
41 #include "GResponse.hpp"
42 #include "GEvent.hpp"
43 #include "GPhoton.hpp"
44 #include "GEnergy.hpp"
45 #include "GTime.hpp"
46 #include "GSource.hpp"
47 #include "GEbounds.hpp"
48 #include "GObservation.hpp"
49 #include "GSkyMap.hpp"
50 #include "GModelSky.hpp"
52 #include "GModelSpatialRadial.hpp"
54 #include "GModelSpatialDiffuse.hpp"
58 #include "GModelSpectralGauss.hpp"
59 
60 /* __ Method name definitions ____________________________________________ */
61 #define G_CONVOLVE "GResponse::convolve(GModelSky&, GObservation&, "\
62  "GMatrixSparse*)"
63 #define G_IRF_RADIAL "GResponse::irf_radial(GEvent&, GSource&,"\
64  " GObservation&)"
65 #define G_IRF_ELLIPTICAL "GResponse::irf_elliptical(GEvent&, GSource&,"\
66  " GObservation&)"
67 #define G_IRF_DIFFUSE "GResponse::irf_diffuse(GEvent&, GSource&,"\
68  " GObservation&)"
69 
70 /* __ Macros _____________________________________________________________ */
71 
72 /* __ Coding definitions _________________________________________________ */
73 
74 /* __ Debug definitions __________________________________________________ */
75 
76 
77 /*==========================================================================
78  = =
79  = Constructors/destructors =
80  = =
81  ==========================================================================*/
82 
83 /***********************************************************************//**
84  * @brief Void constructor
85  ***************************************************************************/
87 {
88  // Initialise class members for clean destruction
89  init_members();
90 
91  // Return
92  return;
93 }
94 
95 
96 /***********************************************************************//**
97  * @brief Copy constructor
98  *
99  * @param[in] rsp Response.
100  ***************************************************************************/
102 {
103  // Initialise class members for clean destruction
104  init_members();
105 
106  // Copy members
107  copy_members(rsp);
108 
109  // Return
110  return;
111 }
112 
113 
114 /***********************************************************************//**
115  * @brief Destructor
116  ***************************************************************************/
118 {
119  // Free members
120  free_members();
121 
122  // Return
123  return;
124 }
125 
126 
127 /*==========================================================================
128  = =
129  = Operators =
130  = =
131  ==========================================================================*/
132 
133 /***********************************************************************//**
134  * @brief Assignment operator
135  *
136  * @param[in] rsp Response.
137  * @return Response.
138  ***************************************************************************/
140 {
141  // Execute only if object is not identical
142  if (this != &rsp) {
143 
144  // Free members
145  free_members();
146 
147  // Initialise private members for clean destruction
148  init_members();
149 
150  // Copy members
151  copy_members(rsp);
152 
153  } // endif: object was not identical
154 
155  // Return this object
156  return *this;
157 }
158 
159 
160 /*==========================================================================
161  = =
162  = Public methods =
163  = =
164  ==========================================================================*/
165 
166 /***********************************************************************//**
167  * @brief Convolve sky model with the instrument response
168  *
169  * @param[in] model Sky model.
170  * @param[in] event Event.
171  * @param[in] obs Observation.
172  * @param[in] grad Should model gradients be computed?
173  * @return Event probability.
174  *
175  * Computes the event probability
176  *
177  * \f[
178  * P(p',E',t') = \int \int \int
179  * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp \, dE \, dt
180  * \f]
181  *
182  * without taking into account any time dispersion. Energy dispersion is
183  * correctly handled by this method. If time dispersion is indeed needed,
184  * an instrument specific method needs to be provided.
185  ***************************************************************************/
186 double GResponse::convolve(const GModelSky& model,
187  const GEvent& event,
188  const GObservation& obs,
189  const bool& grad) const
190 {
191  // Set number of iterations for Romberg integration.
192  static const int iter = 6;
193 
194  // Initialise result
195  double prob = 0.0;
196 
197  // Continue only if the model has a spatial component
198  if (model.spatial() != NULL) {
199 
200  // Get source time (no dispersion)
201  GTime srcTime = event.time();
202 
203  // Case A: Integration
204  if (use_edisp()) {
205 
206  // Retrieve true energy boundaries
207  GEbounds ebounds = this->ebounds(event.energy());
208 
209  // Retrieve model energy boundaries
210  GEbounds ebounds_model = this->ebounds_model(model);
211 
212  // Determine size of result array
213  int size = size_edisp_vector(model, obs, grad);
214 
215  // Initialise vector array
216  GVector array(size);
217 
218  // Loop over all boundaries
219  for (int i = 0; i < ebounds.size(); ++i) {
220 
221  // Get true energy boundaries in MeV
222  double etrue_min = ebounds.emin(i).MeV();
223  double etrue_max = ebounds.emax(i).MeV();
224 
225  // Limit energy boundaries
226  if (ebounds_model.size() > 0) {
227  if (ebounds_model.emin(0).MeV() > etrue_min) {
228  etrue_min = ebounds_model.emin(0).MeV();
229  }
230  if (ebounds_model.emax(0).MeV() < etrue_max) {
231  etrue_max = ebounds_model.emax(0).MeV();
232  }
233  }
234 
235  // Continue only if valid
236  if (etrue_max > etrue_min) {
237 
238  // Setup integration function
239  edisp_kerns integrand(this, &obs, &model, &event, srcTime, grad);
240  GIntegrals integral(&integrand);
241 
242  // Set number of iterations
243  integral.fixed_iter(iter);
244 
245  // Do Romberg integration
246  array += integral.romberg(etrue_min, etrue_max);
247 
248  } // endif: interval was valid
249 
250  } // endfor: looped over intervals
251 
252  // Initialise array index
253  int index = 0;
254 
255  // Get probability
256  prob = array[index++];
257 
258  // Set gradients
259  if (grad) {
260 
261  // Set spectral gradients
262  if (model.spectral() != NULL) {
263  for (int i = 0; i < model.spectral()->size(); ++i) {
264  GModelPar& par = (*(model.spectral()))[i];
265  if (par.is_free() && par.has_grad()) {
266  par.factor_gradient(array[index++]);
267  }
268  }
269  }
270 
271  // Set temporal gradients
272  if (model.temporal() != NULL) {
273  for (int i = 0; i < model.temporal()->size(); ++i) {
274  GModelPar& par = (*(model.temporal()))[i];
275  if (par.is_free() && par.has_grad()) {
276  par.factor_gradient(array[index++]);
277  }
278  }
279  }
280 
281  // Optionally set scale gradient for instrument
282  if (model.has_scales()) {
283  for (int i = 0; i < model.scales(); ++i) {
284  const GModelPar& par = model.scale(i);
285  if (par.name() == obs.instrument()) {
286  if (par.is_free() && par.has_grad()) {
287  par.factor_gradient(array[index++]);
288  }
289  }
290  }
291  }
292 
293  } // endif: set gradients
294 
295  } // endif: Case A
296 
297  // Case B: No integration (assume no energy dispersion)
298  else {
299 
300  // Get source energy (no dispersion)
301  GEnergy srcEng = event.energy();
302 
303  // Evaluate probability
304  prob = eval_prob(model, event, srcEng, srcTime, obs, grad);
305 
306  } // endelse: Case B
307 
308  // Compile option: Check for NaN/Inf
309  #if defined(G_NAN_CHECK)
310  if (gammalib::is_notanumber(prob) || gammalib::is_infinite(prob)) {
311  std::cout << "*** ERROR: GResponse::convolve:";
312  std::cout << " NaN/Inf encountered";
313  std::cout << " (prob=" << prob;
314  std::cout << ", event=" << event;
315  std::cout << ", srcTime=" << srcTime;
316  std::cout << ")" << std::endl;
317  }
318  #endif
319 
320  } // endif: spatial component valid
321 
322  // Return probability
323  return prob;
324 }
325 
326 
327 /***********************************************************************//**
328  * @brief Convolve sky model with the instrument response
329  *
330  * @param[in] model Sky model.
331  * @param[in] obs Observation.
332  * @param[out] gradients Pointer to matrix of gradients.
333  * @return Vector of event probabilities.
334  *
335  * Computes the event probability
336  *
337  * \f[
338  * P(p',E',t') = \int \int \int
339  * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp \, dE \, dt
340  * \f]
341  *
342  * without taking into account any time dispersion. Energy dispersion is
343  * correctly handled by this method. If time dispersion is indeed needed,
344  * an instrument specific method needs to be provided.
345  ***************************************************************************/
347  const GObservation& obs,
348  GMatrixSparse* gradients) const
349 {
350  // Set number of iterations for Romberg integration.
351  static const int iter = 6;
352 
353  // Get number of model parameters and number of events
354  int npars = model.size();
355  int nevents = obs.events()->size();
356 
357  // Initialise gradients flag
358  bool grad = (gradients != NULL);
359 
360  // Check matrix consistency
361  if (grad) {
362 
363  // Check number of columns
364  int ncolumns = (nevents > 0) ? npars : 0;
365  if (gradients->columns() != ncolumns) {
366  std::string msg = "Number of "+gammalib::str(gradients->columns())+" "
367  "columns in gradient matrix differs from expected "
368  "number of "+gammalib::str(ncolumns)+". Please "
369  "specify compatible arguments.";
371  }
372 
373  // Check number of rows
374  int nrows = (ncolumns > 0) ? nevents : 0;
375  if (gradients->rows() != nrows) {
376  std::string msg = "Number of "+gammalib::str(gradients->rows())+" "
377  "rows in gradient matrix differs from expected "
378  "number of "+gammalib::str(nrows)+". Please "
379  "specify compatible arguments.";
381  }
382 
383  } // endif: gadient was requested
384 
385  // Initialise result
386  GVector probs(nevents);
387 
388  // Continue only if the model has a spatial component
389  if (model.spatial() != NULL) {
390 
391  // Case A: Integration
392  if (use_edisp()) {
393 
394  // Initialise temporary vectors to hold gradients
395  GVector* tmp_gradients = NULL;
396  if (grad) {
397  tmp_gradients = new GVector[npars];
398  for (int i = 0; i < npars; ++i) {
399  tmp_gradients[i] = GVector(nevents);
400  }
401  }
402 
403  // Determine size of result array
404  int size = size_edisp_vector(model, obs, grad);
405 
406  // Get model energy boundaries
407  GEbounds ebounds_model = this->ebounds_model(model);
408 
409  // Loop over events
410  for (int k = 0; k < nevents; ++k) {
411 
412  // Get reference to event
413  const GEvent& event = *((*obs.events())[k]);
414 
415  // Get source time (no dispersion)
416  GTime srcTime = event.time();
417 
418  // Retrieve true energy boundaries
419  GEbounds ebounds = this->ebounds(event.energy());
420 
421  // Initialise vector array
422  GVector array(size);
423 
424  // Loop over all boundaries
425  for (int i = 0; i < ebounds.size(); ++i) {
426 
427  // Get true energy boundaries in MeV
428  double etrue_min = ebounds.emin(i).MeV();
429  double etrue_max = ebounds.emax(i).MeV();
430 
431  // Limit energy boundaries
432  if (ebounds_model.size() > 0) {
433  if (ebounds_model.emin(0).MeV() > etrue_min) {
434  etrue_min = ebounds_model.emin(0).MeV();
435  }
436  if (ebounds_model.emax(0).MeV() < etrue_max) {
437  etrue_max = ebounds_model.emax(0).MeV();
438  }
439  }
440 
441  // Continue only if valid
442  if (etrue_max > etrue_min) {
443 
444  // Setup integration function
445  edisp_kerns integrand(this, &obs, &model, &event, srcTime, grad);
446  GIntegrals integral(&integrand);
447 
448  // Set number of iterations
449  integral.fixed_iter(iter);
450 
451  // Do Romberg integration
452  array += integral.romberg(etrue_min, etrue_max);
453 
454  } // endif: interval was valid
455 
456  } // endfor: looped over intervals
457 
458  // Initialise array index
459  int index = 0;
460 
461  // Get probability
462  probs[k] = array[index++];
463 
464  // Set gradients
465  if (grad) {
466 
467  // Get number of spatial, spectral and temporal parameters
468  int n_spat = model.spatial()->size();
469  int n_spec = model.spectral()->size();
470  int n_temp = model.temporal()->size();
471 
472  // Set spectral gradients
473  if (model.spectral() != NULL) {
474  int offset = n_spat;
475  for (int i = 0; i < n_spec; ++i) {
476  const GModelPar& par = (*(model.spectral()))[i];
477  if (par.is_free() && par.has_grad()) {
478  tmp_gradients[offset+i][k] = array[index++];
479  }
480  }
481  }
482 
483  // Set temporal gradients
484  if (model.temporal() != NULL) {
485  int offset = n_spat + n_spec;
486  for (int i = 0; i < n_temp; ++i) {
487  const GModelPar& par = (*(model.temporal()))[i];
488  if (par.is_free() && par.has_grad()) {
489  tmp_gradients[offset+i][k] = array[index++];
490  }
491  }
492  }
493 
494  // Optionally set scale gradient for instrument
495  if (model.has_scales()) {
496  int offset = n_spat + n_spec + n_temp;
497  for (int i = 0; i < model.scales(); ++i) {
498  const GModelPar& par = model.scale(i);
499  if (par.name() == obs.instrument()) {
500  if (par.is_free() && par.has_grad()) {
501  tmp_gradients[offset+i][k] = array[index++];
502  }
503  }
504  }
505  }
506 
507  } // endif: set gradients
508 
509  // Compile option: Check for NaN/Inf
510  #if defined(G_NAN_CHECK)
511  if (gammalib::is_notanumber(probs[k]) || gammalib::is_infinite(probs[k])) {
512  std::cout << "*** ERROR: GResponse::convolve:";
513  std::cout << " NaN/Inf encountered";
514  std::cout << " (probs[" << k << "]=" << probs[k];
515  std::cout << ", event=" << event;
516  std::cout << ", srcTime=" << srcTime;
517  std::cout << ")" << std::endl;
518  }
519  #endif
520 
521  } // endfor: looped over events
522 
523  // Post-process gradients
524  if (grad) {
525 
526  // Fill gradients into matrix
527  for (int i = 0; i < npars; ++i) {
528  gradients->column(i, tmp_gradients[i]);
529  }
530 
531  // Delete temporal gradients
532  delete [] tmp_gradients;
533 
534  } // endif: post-processed gradients
535 
536  } // endif: Case A
537 
538  // Case B: No integration (assume no energy dispersion)
539  else {
540 
541  // Evaluate probability
542  probs = eval_probs(model, obs, gradients);
543 
544  } // endelse: Case B
545 
546  } // endif: spatial component valid
547 
548  // Return probabilities
549  return probs;
550 }
551 
552 
553 /***********************************************************************//**
554  * @brief Return instrument response integrated over the spatial model
555  *
556  * @param[in] event Event.
557  * @param[in] source Source.
558  * @param[in] obs Observation.
559  * @return Instrument response to a spatial model.
560  *
561  * Returns the instrument response for a given event, source and observation
562  * integrated over the spatial model component. The method computes
563  *
564  * \f[
565  * {\tt irf}(p', E', t') = \int_p M_{\rm S}(p | E, t) \,
566  * R(p', E', t' | p, E, t) \, d\,p
567  * \f]
568  *
569  * where
570  * * \f$M_{\rm S}(p | E, t)\f$ is the spatial model component,
571  * * \f$R(p', E', t' | p, E, t)\f$ is the Instrument Response Function (IRF),
572  * * \f$p'\f$ is the measured instrument direction,
573  * * \f$E'\f$ is the measured or reconstructed energy,
574  * * \f$t'\f$ is the measured arrival time,
575  * * \f$p\f$ is the true photon arrival direction,
576  * * \f$E\f$ is the true photon energy, and
577  * * \f$t\f$ is the true trigger time.
578  *
579  * The integration is done over all relevant true sky directions \f$p\f$.
580  *
581  * Depending on the type of the source model the method branches to the
582  * following methods to perform the actual computations
583  *
584  * * irf_ptsrc() - for the handling of a point source
585  * * irf_radial() - for radial models
586  * * irf_elliptical() - for elliptical models
587  * * irf_diffuse() - for diffuse models
588  * * irf_composite() - for composite models
589  *
590  * The method implements a caching mechanism for spatial models that have all
591  * parameters fixed. For those models the instrument response for a given
592  * event and observation is only computed once and then stored in an internal
593  * cache from which it is fetched back in case that the method is called
594  * again for the same event and observation.
595  ***************************************************************************/
596 double GResponse::irf_spatial(const GEvent& event,
597  const GSource& source,
598  const GObservation& obs) const
599 {
600  // Initialise IRF value
601  double irf = 0.0;
602 
603  // Set IRF value attributes
604  std::string name = obs.id() + ":" + source.name();
605  const GInstDir& dir = event.dir();
606  const GEnergy& ereco = event.energy();
607  const GEnergy& etrue = source.energy();
608 
609  // Signal if spatial model has free parameters
610  bool has_free_pars = source.model()->has_free_pars();
611 
612  // If the spatial model component has free parameters, or the response
613  // cache should not be used, or the cache does not contain the requested
614  // IRF value then compute the IRF value for the spatial model.
615  if (has_free_pars ||
616  !m_use_irf_cache ||
617  !m_irf_cache.contains(name, dir, ereco, etrue, &irf)) {
618 
619  // Compute IRF for spatial model
620  switch (source.model()->code()) {
622  irf = irf_ptsrc(event, source, obs);
623  break;
625  irf = irf_radial(event, source, obs);
626  break;
628  irf = irf_elliptical(event, source, obs);
629  break;
631  irf = irf_diffuse(event, source, obs);
632  break;
634  irf = irf_composite(event, source, obs);
635  break;
636  default:
637  break;
638  }
639 
640  } // endif: computed spatial model
641 
642  // If the spatial model has no free parameters and the response cache
643  // should be used then put the IRF value in the response cache.
644  if (!has_free_pars && m_use_irf_cache) {
645  m_irf_cache.set(name, dir, ereco, etrue, irf);
646  }
647 
648  // Return IRF value
649  return irf;
650 }
651 
652 
653 /***********************************************************************//**
654  * @brief Return instrument response vector integrated over the spatial model
655  *
656  * @param[in] model Sky model.
657  * @param[in] obs Observation.
658  * @param[in] gradients Gradients matrix.
659  * @return Instrument response vector to a spatial model.
660  *
661  * Returns the instrument response to a sky model integrated over the spatial
662  * model component for all events in a given observation. The method computes
663  *
664  * \f[
665  * {\tt irf}(p', E', t') = \int_p M_{\rm S}(p | E, t) \,
666  * R(p', E', t' | p, E, t) \, d\,p
667  * \f]
668  *
669  * where
670  * * \f$M_{\rm S}(p | E, t)\f$ is the spatial model component,
671  * * \f$R(p', E', t' | p, E, t)\f$ is the Instrument Response Function (IRF),
672  * * \f$p'\f$ is the measured instrument direction,
673  * * \f$E'\f$ is the measured or reconstructed energy,
674  * * \f$t'\f$ is the measured arrival time,
675  * * \f$p\f$ is the true photon arrival direction,
676  * * \f$E\f$ is the true photon energy, and
677  * * \f$t\f$ is the true trigger time.
678  *
679  * The integration is done over all relevant true sky directions \f$p\f$.
680  *
681  * Depending on the type of the source model the method branches to the
682  * following methods to perform the actual computations
683  *
684  * * irf_ptsrc() - for the handling of a point source
685  * * irf_radial() - for radial models
686  * * irf_elliptical() - for elliptical models
687  * * irf_diffuse() - for diffuse models
688  * * irf_composite() - for composite models
689  *
690  * The method implements a caching mechanism for spatial models that have all
691  * parameters fixed. For those models the instrument response for a given
692  * event and observation is only computed once and then stored in an internal
693  * cache from which it is fetched back in case that the method is called
694  * again for the same event and observation.
695  ***************************************************************************/
697  const GObservation& obs,
698  GMatrix* gradients) const
699 {
700  // Get number of model parameters and number of events
701  int npars = model.size();
702  int nevents = obs.events()->size();
703 
704  // Initialise result
705  GVector irfs(nevents);
706 
707  // Continue only if model has a spatial component
708  if (model.spatial() != NULL) {
709 
710  // Set IRF cache identifier
711  std::string cache_id = obs.id() + ":" + model.name();
712 
713  // Signal if spatial model has free parameters
714  bool has_free_pars = model.spatial()->has_free_pars();
715 
716  // If the spatial model component has free parameters, or the response
717  // cache should not be used, or the cache does not contain the requested
718  // IRF value then compute the IRF value for the spatial model.
719  if (has_free_pars ||
720  !m_use_irf_cache ||
721  !m_irf_vector_cache.contains(cache_id, &irfs)) {
722 
723  // Compute IRF for spatial model
724  switch (model.spatial()->code()) {
726  irfs = irf_ptsrc(model, obs, gradients);
727  break;
729  irfs = irf_radial(model, obs, gradients);
730  break;
732  irfs = irf_elliptical(model, obs, gradients);
733  break;
735  irfs = irf_diffuse(model, obs, gradients);
736  break;
738  irfs = irf_composite(model, obs, gradients);
739  break;
740  default:
741  break;
742  }
743 
744  // If the spatial model has no free parameters and the response cache
745  // should be used then put the IRF value in the response cache.
746  if (!has_free_pars && m_use_irf_cache) {
747  m_irf_vector_cache.set(cache_id, irfs);
748  }
749 
750  } // endif: computed spatial model
751 
752  } // endif: model had spatial component
753 
754  // Return IRF values
755  return irfs;
756 }
757 
758 
759 /***********************************************************************//**
760  * @brief Remove response cache for model
761  *
762  * @param[in] name Model name.
763  *
764  * Remove response cache for model @p name from response cache.
765  ***************************************************************************/
766 void GResponse::remove_response_cache(const std::string& name)
767 {
768  // Remove model from response caches
769  m_irf_cache.remove(name);
770  m_nroi_cache.remove(name);
772 
773  // Return
774  return;
775 }
776 
777 
778 /*==========================================================================
779  = =
780  = Protected methods =
781  = =
782  ==========================================================================*/
783 
784 /***********************************************************************//**
785  * @brief Initialise class members
786  ***************************************************************************/
788 {
789  // Initialize members
790  m_use_irf_cache = true; //!< Switched on by default
791  m_use_nroi_cache = true; //!< Switched on by default
796  m_irf_diffuse_resolution = 0.1; //!< Angular resolution (deg)
797 
798  // Clear cache
799  m_irf_cache.clear();
802 
803  // Return
804  return;
805 }
806 
807 
808 /***********************************************************************//**
809  * @brief Copy class members
810  *
811  * @param[in] rsp Response.
812  ***************************************************************************/
814 {
815  // Copy members
823  m_irf_cache = rsp.m_irf_cache;
826 
827  // Return
828  return;
829 }
830 
831 
832 /***********************************************************************//**
833  * @brief Delete class members
834  ***************************************************************************/
836 {
837  // Return
838  return;
839 }
840 
841 
842 /***********************************************************************//**
843  * @brief Return instrument response to point source
844  *
845  * @param[in] event Observed event.
846  * @param[in] source Source.
847  * @param[in] obs Observation.
848  * @return Instrument response to point source.
849  *
850  * Returns the instrument response to a point source.
851  ***************************************************************************/
852 double GResponse::irf_ptsrc(const GEvent& event,
853  const GSource& source,
854  const GObservation& obs) const
855 {
856  // Get pointer to point source model
857  const GModelSpatialPointSource* model =
858  static_cast<const GModelSpatialPointSource*>(source.model());
859 
860  // Setup photon
861  GPhoton photon(model->dir(), source.energy(), source.time());
862 
863  // Get IRF
864  double irf = this->irf(event, photon, obs);
865 
866  // Return IRF value
867  return irf;
868 }
869 
870 
871 /***********************************************************************//**
872  * @brief Return instrument response to radial source
873  *
874  * @param[in] event Observed event.
875  * @param[in] source Source.
876  * @param[in] obs Observation.
877  * @return Instrument response to radial source.
878  *
879  * Returns the instrument response to a radial source for a given event,
880  * source and observation. The method computes
881  *
882  * \f[
883  * {\tt irf}(p', E', t') = \int_0^{\theta_{\rm max}}
884  * M_{\rm S}(\theta | E, t) \, sin \, \theta
885  * \int_0^{2\pi}
886  * R(p', E', t' | \theta, \phi, E, t) \,
887  * d\,\phi \, d\,\theta
888  * \f]
889  *
890  * where
891  * * \f$M_{\rm S}(\theta | E, t)\f$ is the radial model component,
892  * * \f$R(p', E', t' | \theta, \phi, E, t)\f$ is the Instrument Response
893  * Function (IRF),
894  * * \f$p'\f$ is the measured instrument direction,
895  * * \f$E'\f$ is the measured or reconstructed energy,
896  * * \f$t'\f$ is the measured arrival time,
897  * * \f$\theta\f$ is the radial distance from the model centre,
898  * * \f$\phi\f$ is the azimuthal angle around the model centre,
899  * * \f$E\f$ is the true photon energy, and
900  * * \f$t\f$ is the true trigger time.
901  *
902  * The azimuth integration is done over \f$2\pi\f$ using the
903  * irf_radial_kern_phi::eval() method.
904  * The radial integration is done out to a maximum angle
905  * \f$\theta_{\rm max}\f$, given by GModelSpatialRadial::theta_max(), by the
906  * irf_radial_kern_theta::eval() method.
907  ***************************************************************************/
908 double GResponse::irf_radial(const GEvent& event,
909  const GSource& source,
910  const GObservation& obs) const
911 {
912  // Initialise IRF
913  double irf = 0.0;
914 
915  // Get pointer to radial model
916  const GModelSpatialRadial* model =
917  static_cast<const GModelSpatialRadial*>(source.model());
918 
919  // Get source attributes
920  const GEnergy& srcEng = source.energy();
921  const GTime& srcTime = source.time();
922 
923  // Set radial integration range (radians)
924  double theta_min = 0.0;
925  double theta_max = model->theta_max() - 1.0e-12; // Kludge to stay inside
926  // the boundary of a
927  // sharp edged model
928 
929  // Integrate if interval is valid
930  if (theta_max > theta_min) {
931 
932  // Allocate Euler and rotation matrices
933  GMatrix ry;
934  GMatrix rz;
935 
936  // Set up rotation matrix to rotate from native model coordinates to
937  // celestial coordinates
938  ry.eulery( model->dir().dec_deg() - 90.0);
939  rz.eulerz(-model->dir().ra_deg());
940  GMatrix rot = (ry * rz).transpose();
941 
942  // Setup integration kernel
943  irf_radial_kern_theta integrand(this,
944  &event,
945  &obs,
946  model,
947  &rot,
948  &srcEng,
949  &srcTime,
951 
952  // Integrate over model's radial coordinate
953  GIntegral integral(&integrand);
955 
956  // Integrate kernel
957  irf = integral.romberg(theta_min, theta_max, m_irf_radial_iter_theta);
958 
959  // Compile option: Check for NaN/Inf
960  #if defined(G_NAN_CHECK)
962  std::cout << "*** ERROR: GResponse::irf_radial:";
963  std::cout << " NaN/Inf encountered";
964  std::cout << " (irf=" << irf;
965  std::cout << ", theta_min=" << theta_min;
966  std::cout << ", theta_max=" << theta_max << ")";
967  std::cout << std::endl;
968  }
969  #endif
970 
971  } // endif: integration interval was valid
972 
973  // Return IRF value
974  return irf;
975 }
976 
977 
978 /***********************************************************************//**
979  * @brief Return instrument response to elliptical source
980  *
981  * @param[in] event Observed event.
982  * @param[in] source Source.
983  * @param[in] obs Observation.
984  * @return Instrument response to elliptical source.
985  *
986  * Returns the instrument response to an elliptical source for a given event,
987  * source and observation. The method computes
988  *
989  * \f[
990  * {\tt irf}(p', E', t') = \int_0^{\theta_{\rm max}} \int_0^{2\pi}
991  * M_{\rm S}(\theta, \phi | E, t) \,
992  * R(p', E', t' | \theta, \phi, E, t) \,
993  * sin \, \theta \,
994  * d\,\phi \, d\,\theta
995  * \f]
996  *
997  * where
998  * * \f$M_{\rm S}(\theta, \phi | E, t)\f$ is the elliptical model component,
999  * * \f$R(p', E', t' | \theta, \phi, E, t)\f$ is the Instrument Response
1000  * Function (IRF),
1001  * * \f$p'\f$ is the measured instrument direction,
1002  * * \f$E'\f$ is the measured or reconstructed energy,
1003  * * \f$t'\f$ is the measured arrival time,
1004  * * \f$\theta\f$ is the radial distance from the model centre,
1005  * * \f$\phi\f$ is the azimuthal angle around the model centre,
1006  * * \f$E\f$ is the true photon energy, and
1007  * * \f$t\f$ is the true trigger time.
1008  *
1009  * The azimuth integration is done over \f$2\pi\f$ using the
1010  * irf_elliptical_kern_phi::eval() method.
1011  * The radial integration is done out to a maximum angle
1012  * \f$\theta_{\rm max}\f$, given by GModelSpatialElliptical::theta_max(), by
1013  * the irf_elliptical_kern_theta::eval() method.
1014  ***************************************************************************/
1015 double GResponse::irf_elliptical(const GEvent& event,
1016  const GSource& source,
1017  const GObservation& obs) const
1018 {
1019  // Initialise IRF
1020  double irf = 0.0;
1021 
1022  // Get pointer to elliptical model
1023  const GModelSpatialElliptical* model =
1024  static_cast<const GModelSpatialElliptical*>(source.model());
1025 
1026  // Get source attributes
1027  const GEnergy& srcEng = source.energy();
1028  const GTime& srcTime = source.time();
1029 
1030  // Set radial integration range (radians)
1031  double theta_min = 0.0;
1032  double theta_max = model->theta_max() - 1.0e-6; // Kludge to stay inside
1033  // the boundary of a
1034  // sharp edged model
1035 
1036  // Integrate if interval is valid
1037  if (theta_max > theta_min) {
1038 
1039  // Allocate Euler and rotation matrices
1040  GMatrix ry;
1041  GMatrix rz;
1042 
1043  // Set up rotation matrix to rotate from native model coordinates to
1044  // celestial coordinates
1045  ry.eulery( model->dir().dec_deg() - 90.0);
1046  rz.eulerz(-model->dir().ra_deg());
1047  GMatrix rot = (ry * rz).transpose();
1048 
1049  // Setup integration kernel
1050  irf_elliptical_kern_theta integrand(this,
1051  &event,
1052  &obs,
1053  model,
1054  &rot,
1055  &srcEng,
1056  &srcTime,
1058 
1059  // Integrate over model's radial coordinate
1060  GIntegral integral(&integrand);
1062 
1063  // Integrate kernel
1064  irf = integral.romberg(theta_min, theta_max, m_irf_elliptical_iter_theta);
1065 
1066  // Compile option: Check for NaN/Inf
1067  #if defined(G_NAN_CHECK)
1069  std::cout << "*** ERROR: GResponse::irf_elliptical:";
1070  std::cout << " NaN/Inf encountered";
1071  std::cout << " (irf=" << irf;
1072  std::cout << ", theta_min=" << theta_min;
1073  std::cout << ", theta_max=" << theta_max << ")";
1074  std::cout << std::endl;
1075  }
1076  #endif
1077 
1078  } // endif: integration interval was valid
1079 
1080  // Return IRF value
1081  return irf;
1082 }
1083 
1084 
1085 /***********************************************************************//**
1086  * @brief Return instrument response to diffuse source
1087  *
1088  * @param[in] event Observed event.
1089  * @param[in] source Source.
1090  * @param[in] obs Observation.
1091  * @return Instrument response to diffuse source.
1092  *
1093  * Returns the instrument response to a diffuse source for a given event,
1094  * source and observation. The method computes
1095  *
1096  * \f[
1097  * {\tt irf}(p', E', t') = \sum_p M_{\rm S}(p | E, t) \, \Delta(p)
1098  * R(p', E', t' | p, E, t)
1099  * \f]
1100  *
1101  * where
1102  * * \f$M_{\rm S}(p | E, t)\f$ is the diffuse model component,
1103  * * \f$\Delta(p)\f$ is the solid angle of the map pixel
1104  * * \f$R(p', E', t' | p, E, t)\f$ is the Instrument Response Function (IRF),
1105  * * \f$p'\f$ is the measured instrument direction,
1106  * * \f$E'\f$ is the measured or reconstructed energy,
1107  * * \f$t'\f$ is the measured arrival time,
1108  * * \f$p\f$ is the true photon arrival direction,
1109  * * \f$E\f$ is the true photon energy, and
1110  * * \f$t\f$ is the true trigger time.
1111  *
1112  * The method simply adds up all sky map pixels multiplied with the IRF for
1113  * all diffuse model pixels. For models that do not contain sky map pixels
1114  * (such as the GModelSpatialDiffuseConst model) the method allocates an
1115  * all-sky sky map with an angular resolution that is specified by the
1116  * m_irf_diffuse_resolution member.
1117  ***************************************************************************/
1118 double GResponse::irf_diffuse(const GEvent& event,
1119  const GSource& source,
1120  const GObservation& obs) const
1121 {
1122  // Initialise IRF
1123  double irf = 0.0;
1124 
1125  // Get source attributes
1126  const GEnergy& srcEng = source.energy();
1127  const GTime& srcTime = source.time();
1128 
1129  // If model is a diffuse map model then extract sky map
1130  const GModelSpatialDiffuseMap* map =
1131  dynamic_cast<const GModelSpatialDiffuseMap*>(source.model());
1132  if (map != NULL) {
1133 
1134  // Get reference to sky map
1135  const GSkyMap& skymap = map->map();
1136 
1137  // Loop over all sky map pixels
1138  for (int i = 0; i < skymap.npix(); ++i) {
1139 
1140  // Get map value
1141  double value = skymap(i);
1142 
1143  // Go to next pixel if map value is not positive
1144  if (value <= 0.0) {
1145  continue;
1146  }
1147 
1148  // Get sky direction
1149  GSkyDir srcDir = skymap.inx2dir(i);
1150 
1151  // Setup photon
1152  GPhoton photon(srcDir, srcEng, srcTime);
1153 
1154  // Add IRF multiplied by flux in map pixel. Since the map value
1155  // is per steradian we have to multiply with the solid angle of
1156  // the map pixel
1157  irf += this->irf(event, photon, obs) * value * skymap.solidangle(i);
1158 
1159  } // endfor: looped over all sky map pixels
1160 
1161  } // endif: model was diffuse map
1162 
1163  // ... otherwise if model is a cube map model then extract sky map
1164  else {
1165  const GModelSpatialDiffuseCube* cube =
1166  dynamic_cast<const GModelSpatialDiffuseCube*>(source.model());
1167  if (cube != NULL) {
1168 
1169  // Get reference to sky cube
1170  const GSkyMap& skymap = cube->cube();
1171 
1172  // Loop over all sky map pixels
1173  for (int i = 0; i < skymap.npix(); ++i) {
1174 
1175  // Get sky direction
1176  GSkyDir srcDir = skymap.inx2dir(i);
1177 
1178  // Setup photon
1179  GPhoton photon(srcDir, srcEng, srcTime);
1180 
1181  // Get map value
1182  double value = cube->eval(photon);
1183 
1184  // Go to next pixel if map value is not positive
1185  if (value <= 0.0) {
1186  continue;
1187  }
1188 
1189  // Add IRF multiplied by flux in map pixel. Since the map
1190  // value is per steradian we have to multiply with the solid
1191  // angle of the map pixel
1192  irf += this->irf(event, photon, obs) * value * skymap.solidangle(i);
1193 
1194  } // endfor: looped over all sky map pixels
1195 
1196  } // endelse: model was diffuse cube
1197 
1198  // ... otherwise allocate all sky map with specified resolution
1199  else {
1200 
1201  // Allocate sky map
1202  int nx = int(360.0/m_irf_diffuse_resolution + 0.5);
1203  int ny = int(180.0/m_irf_diffuse_resolution + 0.5);
1204  double dx = 360.0 / double(nx);
1205  double dy = 180.0 / double(ny);
1206  GSkyMap skymap = GSkyMap("CAR", "CEL", 0.0, 0.0, -dx, dy, nx, ny);
1207 
1208  // Get pointer to diffuse model
1209  const GModelSpatialDiffuse* model =
1210  static_cast<const GModelSpatialDiffuse*>(source.model());
1211 
1212  // Loop over all sky map pixels
1213  for (int i = 0; i < skymap.npix(); ++i) {
1214 
1215  // Get sky direction
1216  GSkyDir srcDir = skymap.inx2dir(i);
1217 
1218  // Setup photon
1219  GPhoton photon(srcDir, srcEng, srcTime);
1220 
1221  // Get model value
1222  double value = model->eval(photon);
1223 
1224  // Go to next pixel if map value is not positive
1225  if (value <= 0.0) {
1226  continue;
1227  }
1228 
1229  // Add IRF multiplied by flux in map pixel. Since the map
1230  // value is per steradian we have to multiply with the solid
1231  // angle of the map pixel
1232  irf += this->irf(event, photon, obs) * value * skymap.solidangle(i);
1233 
1234  } // endfor: looped over sky map pixels
1235 
1236  } // endelse: model was not a cube
1237 
1238  } // endelse: model was not a map
1239 
1240  // Return IRF value
1241  return irf;
1242 }
1243 
1244 
1245 /***********************************************************************//**
1246  * @brief Return instrument response to composite source
1247  *
1248  * @param[in] event Observed event.
1249  * @param[in] source Source.
1250  * @param[in] obs Observation.
1251  * @return Instrument response to composite source.
1252  *
1253  * Returns the instrument response to a specified composite source.
1254  ***************************************************************************/
1255 double GResponse::irf_composite(const GEvent& event,
1256  const GSource& source,
1257  const GObservation& obs) const
1258 {
1259  // Initialise IRF
1260  double irf = 0.0;
1261 
1262  // Get pointer to composite model
1263  const GModelSpatialComposite* model =
1264  dynamic_cast<const GModelSpatialComposite*>(source.model());
1265 
1266  // Loop over model components
1267  for (int i = 0; i < model->components(); ++i) {
1268 
1269  // Get pointer to spatial component
1270  GModelSpatial* spat = const_cast<GModelSpatial*>(model->component(i));
1271 
1272  // Create new GSource object
1273  GSource src(source.name(), spat, source.energy(), source.time());
1274 
1275  // Compute irf value
1276  irf += irf_spatial(event, src, obs) * model->scale(i);
1277 
1278  }
1279 
1280  // Divide by number of model components
1281  double sum = model->sum_of_scales();
1282  if (sum > 0.0) {
1283  irf /= sum;
1284  }
1285 
1286  // Return IRF value
1287  return irf;
1288 }
1289 
1290 
1291 /***********************************************************************//**
1292  * @brief Return instrument response to point source sky model
1293  *
1294  * @param[in] model Sky model.
1295  * @param[in] obs Observation.
1296  * @param[in] gradients Gradients matrix.
1297  * @return Instrument response to point source sky model.
1298  *
1299  * Returns the instrument response to a point source sky model for all
1300  * events.
1301  ***************************************************************************/
1303  const GObservation& obs,
1304  GMatrix* gradients) const
1305 {
1306  // Get number of events
1307  int nevents = obs.events()->size();
1308 
1309  // Initialise result
1310  GVector irfs(nevents);
1311 
1312  // Loop over events
1313  for (int k = 0; k < nevents; ++k) {
1314 
1315  // Get reference to event
1316  const GEvent& event = *((*obs.events())[k]);
1317 
1318  // Get source energy and time (no dispersion)
1319  GEnergy srcEng = event.energy();
1320  GTime srcTime = event.time();
1321 
1322  // Setup source
1323  GSource source(model.name(), model.spatial(), srcEng, srcTime);
1324 
1325  // Get IRF value for event
1326  irfs[k] = this->irf_ptsrc(event, source, obs);
1327 
1328  } // endfor: looped over events
1329 
1330  // Return IRF value
1331  return irfs;
1332 }
1333 
1334 
1335 /***********************************************************************//**
1336  * @brief Return instrument response to radial source sky model
1337  *
1338  * @param[in] model Sky model.
1339  * @param[in] obs Observation.
1340  * @param[in] gradients Gradients matrix.
1341  * @return Instrument response to radial source sky model.
1342  *
1343  * Returns the instrument response to a radial source sky model for all
1344  * events.
1345  ***************************************************************************/
1347  const GObservation& obs,
1348  GMatrix* gradients) const
1349 {
1350  // Get number of events
1351  int nevents = obs.events()->size();
1352 
1353  // Initialise result
1354  GVector irfs(nevents);
1355 
1356  // Loop over events
1357  for (int k = 0; k < nevents; ++k) {
1358 
1359  // Get reference to event
1360  const GEvent& event = *((*obs.events())[k]);
1361 
1362  // Get source energy and time (no dispersion)
1363  GEnergy srcEng = event.energy();
1364  GTime srcTime = event.time();
1365 
1366  // Setup source
1367  GSource source(model.name(), model.spatial(), srcEng, srcTime);
1368 
1369  // Get IRF value for event
1370  irfs[k] = this->irf_radial(event, source, obs);
1371 
1372  } // endfor: looped over events
1373 
1374  // Return IRF value
1375  return irfs;
1376 }
1377 
1378 
1379 /***********************************************************************//**
1380  * @brief Return instrument response to ellipitical source sky model
1381  *
1382  * @param[in] model Sky model.
1383  * @param[in] obs Observation.
1384  * @param[in] gradients Gradients matrix.
1385  * @return Instrument response to ellipitical source sky model.
1386  *
1387  * Returns the instrument response to a ellipitical source sky model for all
1388  * events.
1389  ***************************************************************************/
1391  const GObservation& obs,
1392  GMatrix* gradients) const
1393 {
1394  // Get number of events
1395  int nevents = obs.events()->size();
1396 
1397  // Initialise result
1398  GVector irfs(nevents);
1399 
1400  // Loop over events
1401  for (int k = 0; k < nevents; ++k) {
1402 
1403  // Get reference to event
1404  const GEvent& event = *((*obs.events())[k]);
1405 
1406  // Get source energy and time (no dispersion)
1407  GEnergy srcEng = event.energy();
1408  GTime srcTime = event.time();
1409 
1410  // Setup source
1411  GSource source(model.name(), model.spatial(), srcEng, srcTime);
1412 
1413  // Get IRF value for event
1414  irfs[k] = this->irf_elliptical(event, source, obs);
1415 
1416  } // endfor: looped over events
1417 
1418  // Return IRF value
1419  return irfs;
1420 }
1421 
1422 
1423 /***********************************************************************//**
1424  * @brief Return instrument response to diffuse source sky model
1425  *
1426  * @param[in] model Sky model.
1427  * @param[in] obs Observation.
1428  * @param[in] gradients Gradients matrix.
1429  * @return Instrument response to diffuse source sky model.
1430  *
1431  * Returns the instrument response to a diffuse source sky model for all
1432  * events.
1433  ***************************************************************************/
1435  const GObservation& obs,
1436  GMatrix* gradients) const
1437 {
1438  // Get number of events
1439  int nevents = obs.events()->size();
1440 
1441  // Initialise result
1442  GVector irfs(nevents);
1443 
1444  // Loop over events
1445  for (int k = 0; k < nevents; ++k) {
1446 
1447  // Get reference to event
1448  const GEvent& event = *((*obs.events())[k]);
1449 
1450  // Get source energy and time (no dispersion)
1451  GEnergy srcEng = event.energy();
1452  GTime srcTime = event.time();
1453 
1454  // Setup source
1455  GSource source(model.name(), model.spatial(), srcEng, srcTime);
1456 
1457  // Get IRF value for event
1458  irfs[k] = this->irf_diffuse(event, source, obs);
1459 
1460  } // endfor: looped over events
1461 
1462  // Return IRF value
1463  return irfs;
1464 }
1465 
1466 
1467 /***********************************************************************//**
1468  * @brief Return instrument response to composite source sky model
1469  *
1470  * @param[in] model Sky model.
1471  * @param[in] obs Observation.
1472  * @param[in] gradients Gradients matrix.
1473  * @return Instrument response to composite source sky model.
1474  *
1475  * Returns the instrument response to a composite source sky model for all
1476  * events.
1477  ***************************************************************************/
1479  const GObservation& obs,
1480  GMatrix* gradients) const
1481 {
1482  // Get number of events
1483  int nevents = obs.events()->size();
1484 
1485  // Initialise result
1486  GVector irfs(nevents);
1487 
1488  // Get pointer to composite spatial model
1489  const GModelSpatialComposite* composite =
1490  dynamic_cast<const GModelSpatialComposite*>(model.spatial());
1491 
1492  // Create copy of sky model
1493  GModelSky sky = model;
1494 
1495  // Loop over model components
1496  for (int i = 0; i < composite->components(); ++i) {
1497 
1498  // Set spatial component of sky model
1499  sky.spatial(const_cast<GModelSpatial*>(composite->component(i)));
1500 
1501  // Compute and add IRF values
1502  irfs += irf_spatial(sky, obs, gradients) * composite->scale(i);
1503 
1504  } // endfor: looped over all model components
1505 
1506  // Divide by number of model components
1507  double sum = composite->sum_of_scales();
1508  if (sum > 0.0) {
1509  irfs /= sum;
1510  }
1511 
1512  // Return IRF values
1513  return irfs;
1514 }
1515 
1516 
1517 /***********************************************************************//**
1518  * @brief Convolve sky model with the instrument response
1519  *
1520  * @param[in] model Sky model.
1521  * @param[in] event Event.
1522  * @param[in] srcEng Source energy.
1523  * @param[in] srcTime Source time.
1524  * @param[in] obs Observation.
1525  * @param[in] grad Should model gradients be computed? (default: true)
1526  * @return Event probability.
1527  *
1528  * Computes the event probability
1529  *
1530  * \f[
1531  * P(p',E',t'|E,t) = \int S(p,E,t) \times R(p',E',t'|p,E,t) \, dp
1532  * \f]
1533  *
1534  * for a given true energy \f$E\f$ and time \f$t\f$.
1535  ***************************************************************************/
1536 double GResponse::eval_prob(const GModelSky& model,
1537  const GEvent& event,
1538  const GEnergy& srcEng,
1539  const GTime& srcTime,
1540  const GObservation& obs,
1541  const bool& grad) const
1542 {
1543  // Initialise result
1544  double prob = 0.0;
1545 
1546  // Continue only if the model has a spatial component
1547  if (model.spatial() != NULL) {
1548 
1549  // Set source
1550  GSource source(model.name(), model.spatial(), srcEng, srcTime);
1551 
1552  // Compute IRF value
1553  double irf = irf_spatial(event, source, obs);
1554 
1555  // Continue only if IRF value is positive
1556  if (irf > 0.0) {
1557 
1558  // Optionally get scaling
1559  double scale = (model.has_scales())
1560  ? model.scale(obs.instrument()).value() : 1.0;
1561 
1562  // Evaluate spectral and temporal components
1563  double spec = (model.spectral() != NULL)
1564  ? model.spectral()->eval(srcEng, srcTime, grad) : 1.0;
1565  double temp = (model.temporal() != NULL)
1566  ? model.temporal()->eval(srcTime, grad) : 1.0;
1567 
1568  // Compute probability
1569  prob = spec * temp * irf * scale;
1570 
1571  // Optionally compute partial derivatives
1572  if (grad) {
1573 
1574  // Multiply factors to spectral gradients
1575  if (model.spectral() != NULL) {
1576  double fact = temp * irf * scale;
1577  if (fact != 1.0) {
1578  for (int i = 0; i < model.spectral()->size(); ++i) {
1579  (*model.spectral())[i].factor_gradient((*model.spectral())[i].factor_gradient() * fact);
1580  }
1581  }
1582  }
1583 
1584  // Multiply factors to temporal gradients
1585  if (model.temporal() != NULL) {
1586  double fact = spec * irf * scale;
1587  if (fact != 1.0) {
1588  for (int i = 0; i < model.temporal()->size(); ++i) {
1589  (*model.temporal())[i].factor_gradient((*model.temporal())[i].factor_gradient() * fact);
1590  }
1591  }
1592  }
1593 
1594  // Optionally compute scale gradient for instrument
1595  if (model.has_scales()) {
1596  for (int i = 0; i < model.scales(); ++i) {
1597  double g_scale = 0.0;
1598  if (model.scale(i).name() == obs.instrument()) {
1599  if (model.scale(i).is_free()) {
1600  g_scale = model.scale(i).scale() * spec * temp * irf;
1601  }
1602  }
1603  model.scale(i).factor_gradient(g_scale);
1604  }
1605  }
1606 
1607  } // endif: computed partial derivatives
1608 
1609  // Compile option: Check for NaN/Inf
1610  #if defined(G_NAN_CHECK)
1611  if (gammalib::is_notanumber(prob) || gammalib::is_infinite(prob)) {
1612  std::cout << "*** ERROR: GResponse::eval_prob:";
1613  std::cout << " NaN/Inf encountered";
1614  std::cout << " (prob=" << prob;
1615  std::cout << ", spec=" << spec;
1616  std::cout << ", temp=" << temp;
1617  std::cout << ", irf=" << irf;
1618  std::cout << ")" << std::endl;
1619  }
1620  #endif
1621 
1622  } // endif: IRF value was positive
1623 
1624  // ... otherwise if gradient computation is requested then set the
1625  // spectral and temporal gradients to zero
1626  else if (grad) {
1627 
1628  // Reset spectral gradients
1629  if (model.spectral() != NULL) {
1630  for (int i = 0; i < model.spectral()->size(); ++i) {
1631  (*model.spectral())[i].factor_gradient(0.0);
1632  }
1633  }
1634 
1635  // Reset temporal gradients
1636  if (model.temporal() != NULL) {
1637  for (int i = 0; i < model.temporal()->size(); ++i) {
1638  (*model.temporal())[i].factor_gradient(0.0);
1639  }
1640  }
1641 
1642  // Optionally reset scale gradients
1643  if (model.has_scales()) {
1644  for (int i = 0; i < model.scales(); ++i) {
1645  model.scale(i).factor_gradient(0.0);
1646  }
1647  }
1648 
1649  } // endelse: IRF value was not positive
1650 
1651  } // endif: Gamma-ray source model had a spatial component
1652 
1653  // Return event probability
1654  return prob;
1655 }
1656 
1657 
1658 /***********************************************************************//**
1659  * @brief Convolve sky model with the instrument response
1660  *
1661  * @param[in] model Sky model.
1662  * @param[in] obs Observation.
1663  * @param[out] gradients Pointer to matrix of gradients.
1664  * @return Event probabilities.
1665  *
1666  * Computes the event probability
1667  *
1668  * \f[
1669  * P(p',E',t'|E,t) = \int S(p,E,t) \times R(p',E',t'|p,E,t) \, dp
1670  * \f]
1671  *
1672  * for a all events.
1673  ***************************************************************************/
1675  const GObservation& obs,
1676  GMatrixSparse* gradients) const
1677 {
1678  // Get number of model parameters and number of events
1679  int npars = model.size();
1680  int nevents = obs.events()->size();
1681 
1682  // Initialise gradients flag
1683  bool grad = (gradients != NULL);
1684 
1685  // Initialise result
1686  GVector probs(nevents);
1687 
1688  // Continue only if the model has a spatial component
1689  if (model.spatial() != NULL) {
1690 
1691  // Initialise matrix to hold spatial gradients
1692  GMatrix spat_gradients(nevents, model.spatial()->size());
1693 
1694  // Compute IRF value
1695  probs = irf_spatial(model, obs, &spat_gradients);
1696 
1697  // Get global model scaling
1698  double scale = (model.has_scales())
1699  ? model.scale(obs.instrument()).value() : 1.0;
1700 
1701  // Initialise temporary vectors to hold gradients
1702  GVector* tmp_gradients = NULL;
1703  if (grad) {
1704  tmp_gradients = new GVector[npars];
1705  for (int i = 0; i < npars; ++i) {
1706  tmp_gradients[i] = GVector(nevents);
1707  }
1708  }
1709 
1710  // Loop over events
1711  for (int k = 0; k < nevents; ++k) {
1712 
1713  // Get probability
1714  double spat = probs[k];
1715 
1716  // Continue only if spatial value is positive
1717  if (spat > 0.0) {
1718 
1719  // Get reference to event
1720  const GEvent& event = *((*obs.events())[k]);
1721 
1722  // Get source energy and time
1723  GEnergy srcEng = event.energy();
1724  GTime srcTime = event.time();
1725 
1726  // Evaluate spectral and temporal components
1727  double spec = (model.spectral() != NULL)
1728  ? model.spectral()->eval(srcEng, srcTime, grad)
1729  : 1.0;
1730  double temp = (model.temporal() != NULL)
1731  ? model.temporal()->eval(srcTime, grad)
1732  : 1.0;
1733 
1734  // Compute probability
1735  probs[k] = spat * spec * temp * scale;
1736 
1737  // Optionally compute partial derivatives
1738  if (grad) {
1739 
1740  // Get number of spatial, spectral and temporal parameters
1741  int n_spat = model.spatial()->size();
1742  int n_spec = model.spectral()->size();
1743  int n_temp = model.temporal()->size();
1744 
1745  // Multiply factors to spatial gradients
1746  if (model.spatial() != NULL) {
1747  double fact = spec * temp * scale;
1748  for (int i = 0; i < n_spat; ++i) {
1749  tmp_gradients[i][k] = spat_gradients(k,i) * fact;
1750  }
1751  }
1752 
1753  // Multiply factors to spectral gradients
1754  if (model.spectral() != NULL) {
1755  int offset = n_spat;
1756  double fact = spat * temp * scale;
1757  for (int i = 0; i < n_spec; ++i) {
1758  tmp_gradients[offset+i][k] =
1759  (*(model.spectral()))[i].factor_gradient() * fact;
1760  }
1761  }
1762 
1763  // Multiply factors to temporal gradients
1764  if (model.temporal() != NULL) {
1765  int offset = n_spat + n_spec;
1766  double fact = spat * spec * scale;
1767  for (int i = 0; i < n_temp; ++i) {
1768  tmp_gradients[offset+i][k] =
1769  (*(model.temporal()))[i].factor_gradient() * fact;
1770  }
1771  }
1772 
1773  // Optionally set scale gradient for instrument
1774  if (model.has_scales()) {
1775  int offset = n_spat + n_spec + n_temp;
1776  double fact = spat * spec * temp;
1777  for (int i = 0; i < model.scales(); ++i) {
1778  const GModelPar& par = model.scale(i);
1779  double g_scale = 0.0;
1780  if (par.name() == obs.instrument()) {
1781  if (par.is_free()) {
1782  tmp_gradients[offset+i][k] = par.scale() * fact;
1783  }
1784  }
1785  }
1786  }
1787 
1788  } // endif: computed partial derivatives
1789 
1790  } // endif: IRF spatial value was positive
1791 
1792  // Compile option: Check for NaN/Inf
1793  #if defined(G_NAN_CHECK)
1794  if (gammalib::is_notanumber(probs[k]) || gammalib::is_infinite(probs[k])) {
1795  std::cout << "*** ERROR: GResponse::eval_prob:";
1796  std::cout << " NaN/Inf encountered";
1797  std::cout << " (probs[" << k << "]=" << probs[k];
1798  std::cout << ", spat=" << spat;
1799  std::cout << ")" << std::endl;
1800  }
1801  #endif
1802 
1803  } // endfor: looped over events
1804 
1805  // Post-process gradients
1806  if (grad) {
1807 
1808  // Fill gradients into matrix
1809  for (int i = 0; i < npars; ++i) {
1810  gradients->column(i, tmp_gradients[i]);
1811  }
1812 
1813  // Delete temporal gradients
1814  delete [] tmp_gradients;
1815 
1816  } // endif: gradient post processing
1817 
1818  } // endif: Gamma-ray source model had a spatial component
1819 
1820  // Return event probabilities
1821  return probs;
1822 }
1823 
1824 
1825 /***********************************************************************//**
1826  * @brief Return size of vector for energy dispersion computation
1827  *
1828  * @param[in] model Sky model.
1829  * @param[in] obs Observation.
1830  * @param[out] grad Signals whether gradients should be computed.
1831  * @return Size of vector for energy dispersion computation.
1832  *
1833  * Computes the size of the vector that will be computed for the computation
1834  * of the energy dispersion.
1835  ***************************************************************************/
1837  const GObservation& obs,
1838  const bool& grad) const
1839 {
1840  // Initialise vector size
1841  int size = 1;
1842 
1843  // Determine size of gradient part
1844  if (grad) {
1845  if (model.spectral() != NULL) {
1846  for (int i = 0; i < model.spectral()->size(); ++i) {
1847  GModelPar& par = (*(model.spectral()))[i];
1848  if (par.is_free() && par.has_grad()) {
1849  size++;
1850  }
1851  }
1852  }
1853  if (model.temporal() != NULL) {
1854  for (int i = 0; i < model.temporal()->size(); ++i) {
1855  GModelPar& par = (*(model.temporal()))[i];
1856  if (par.is_free() && par.has_grad()) {
1857  size++;
1858  }
1859  }
1860  }
1861  if (model.has_scales()) {
1862  for (int i = 0; i < model.scales(); ++i) {
1863  GModelPar& par = const_cast<GModelPar&>(model.scale(i));
1864  if (par.name() == obs.instrument()) {
1865  if (par.is_free() && par.has_grad()) {
1866  size++;
1867  }
1868  }
1869  }
1870  }
1871  }
1872 
1873  // Return size
1874  return size;
1875 }
1876 
1877 
1878 /***********************************************************************//**
1879  * @brief Constructor for energy dispersion integration kernels class
1880  *
1881  * @param[in] parent Pointer to response.
1882  * @param[in] obs Pointer to observation.
1883  * @param[in] model Pointer to sky model.
1884  * @param[in] event Pointer to event.
1885  * @param[in] srcTime True time.
1886  * @param[in] grad Evaluate gradients?
1887  *
1888  * This method constructs the integration kernel needed for the energy
1889  * dispersion computation.
1890  ***************************************************************************/
1892  const GObservation* obs,
1893  const GModelSky* model,
1894  const GEvent* event,
1895  const GTime& srcTime,
1896  const bool& grad)
1897 {
1898  // Set members
1899  m_parent = parent;
1900  m_obs = obs;
1901  m_model = model;
1902  m_event = event;
1903  m_size = 1;
1904  m_srcTime = srcTime;
1905  m_grad = grad;
1906 
1907  // If gradients are requested then put pointers to all relevant parameter
1908  // in the parameter pointer vector
1909  m_pars.clear();
1910  if (m_grad) {
1911 
1912  // Add free spectral parameters with analytical gradients
1913  if (m_model->spectral() != NULL) {
1914  for (int i = 0; i < m_model->spectral()->size(); ++i) {
1915  GModelPar& par = (*(m_model->spectral()))[i];
1916  if (par.is_free() && par.has_grad()) {
1917  m_pars.push_back(&par);
1918  }
1919  }
1920  }
1921 
1922  // Add free temporal parameters with analytical gradients
1923  if (m_model->temporal() != NULL) {
1924  for (int i = 0; i < m_model->temporal()->size(); ++i) {
1925  GModelPar& par = (*(m_model->temporal()))[i];
1926  if (par.is_free() && par.has_grad()) {
1927  m_pars.push_back(&par);
1928  }
1929  }
1930  }
1931 
1932  // Add free scaling factors
1933  if (m_model->has_scales()) {
1934  for (int i = 0; i < m_model->scales(); ++i) {
1935  GModelPar& par = const_cast<GModelPar&>(m_model->scale(i));
1936  if (par.name() == m_obs->instrument()) {
1937  if (par.is_free() && par.has_grad()) {
1938  m_pars.push_back(&par);
1939  }
1940  }
1941  }
1942  }
1943 
1944  // Update size
1945  m_size += m_pars.size();
1946 
1947  } // endif: gradients were requested
1948 
1949  // Return
1950  return;
1951 }
1952 
1953 
1954 /***********************************************************************//**
1955  * @brief Return true energy intervals for sky model
1956  *
1957  * @param[in] model Sky model.
1958  * @return True energy intervals.
1959  *
1960  * Returns the true energy intervals for a sky model. For all spectral
1961  * models other than the GModelSpectralGauss model the method will return
1962  * an empty energy boundaries structure. For the GModelSpectralGauss model
1963  * the method will return the interval [mean - 5 * sigma, mean + 5 * sigma].
1964  ***************************************************************************/
1966 {
1967  // Initialise empty energy boundaries
1968  GEbounds ebounds;
1969 
1970  // If the spectral model is a Gaussian line then return
1971  const GModelSpectralGauss* gauss =
1972  dynamic_cast<const GModelSpectralGauss*>(model.spectral());
1973  if (gauss != NULL) {
1974 
1975  // Define width of energy interval
1976  const double width = 5.0;
1977 
1978  // Compute energy boundaries
1979  GEnergy mean = gauss->mean();
1980  GEnergy sigma = gauss->sigma();
1981  GEnergy ebound_min = mean - width * sigma;
1982  GEnergy ebound_max = mean + width * sigma;
1983 
1984  // Make sure that energies are not negative positive
1985  if (ebound_min.MeV() < 0.0) {
1986  ebound_min.MeV(0.0);
1987  }
1988  if (ebound_max.MeV() < 0.0) {
1989  ebound_max.MeV(0.0);
1990  }
1991 
1992  // If energy boundaries define a positive interval then append them
1993  if (ebound_max > ebound_min) {
1994  ebounds.append(ebound_min, ebound_max);
1995  }
1996 
1997  } // endif: spectral model was a Gaussian
1998 
1999  // Return energy boudnaries
2000  return ebounds;
2001 }
2002 
2003 
2004 /***********************************************************************//**
2005  * @brief Evaluate energy dispersion integration kernel
2006  *
2007  * @param[in] etrue True photon energy in MeV.
2008  *
2009  * This method implements the integration kernel needed for the
2010  * GResponse::edisp_kern() class.
2011  ***************************************************************************/
2013 {
2014  // Initialise result
2015  GVector kernels(m_size);
2016 
2017  // Initialise Ndarray array index
2018  int index = 0;
2019 
2020  // Set true energy
2021  GEnergy srcEng;
2022  srcEng.MeV(etrue);
2023 
2024  // Get function value
2025  kernels[index++] = m_parent->eval_prob(*m_model, *m_event,
2026  srcEng, m_srcTime, *m_obs,
2027  m_grad);
2028 
2029  // If gradients are requested then extract them and put them into the
2030  // array
2031  if (m_grad) {
2032  for (int i = 0; i < m_pars.size(); ++i, ++index) {
2033  kernels[index] = m_pars[i]->factor_gradient();
2034  }
2035  }
2036 
2037  // Compile option: Check for NaN
2038  #if defined(G_NAN_CHECK)
2039  if (gammalib::is_notanumber(kernels[0]) || gammalib::is_infinite(kernels[0])) {
2040  std::cout << "*** ERROR: GResponse::edisp_kern::eval";
2041  std::cout << "(etrue=" << etrue << "): NaN/Inf encountered";
2042  std::cout << " (value=" << kernels[0] << ")" << std::endl;
2043  }
2044  #endif
2045 
2046  // Return kernels
2047  return kernels;
2048 }
2049 
2050 
2051 /***********************************************************************//**
2052  * @brief Zenith angle integration kernel for radial model
2053  *
2054  * @param[in] theta Zenith angle (radians).
2055  * @return IRF value.
2056  *
2057  * Integrated the IRF multiplied by the model value for a given zenith angle
2058  * over all azimuth angles.
2059  ***************************************************************************/
2060 double GResponse::irf_radial_kern_theta::eval(const double& theta)
2061 {
2062  // Initialise result
2063  double irf = 0.0;
2064 
2065  // Continue only for positive zenith angles (otherwise the integral will
2066  // be zero)
2067  if (theta > 0.0) {
2068 
2069  // Evaluate sky model
2070  double model = m_model->eval(theta, *m_srcEng, *m_srcTime);
2071 
2072  // Continue only if model is positive
2073  if (model > 0.0) {
2074 
2075  // Set azimuthal integration range
2076  double phi_min = 0.0;
2077  double phi_max = gammalib::twopi;
2078 
2079  // Setup integration kernel
2080  irf_radial_kern_phi integrand(m_rsp,
2081  m_event,
2082  m_obs,
2083  m_rot,
2084  theta,
2085  m_srcEng,
2086  m_srcTime);
2087 
2088  // Setup integration
2089  GIntegral integral(&integrand);
2090  integral.fixed_iter(m_iter_phi);
2091 
2092  // Integrate over phi
2093  irf = integral.romberg(phi_min, phi_max, m_iter_phi) *
2094  model * std::sin(theta);
2095 
2096  } // endif: model was positive
2097 
2098  } // endif: theta was positive
2099 
2100  // Return result
2101  return irf;
2102 }
2103 
2104 
2105 /***********************************************************************//**
2106  * @brief Azimuth angle integration kernel for radial model
2107  *
2108  * @param[in] phi Azimuth angle (radians).
2109  * @return IRF value.
2110  *
2111  * Computes the IRF at a given zenith and azimuth angle with respect to a
2112  * specified centre.
2113  ***************************************************************************/
2114 double GResponse::irf_radial_kern_phi::eval(const double& phi)
2115 {
2116  // Set up native coordinate vector
2117  double cos_phi = std::cos(phi);
2118  double sin_phi = std::sin(phi);
2119  GVector native(-cos_phi*m_sin_theta, sin_phi*m_sin_theta, m_cos_theta);
2120 
2121  // Rotate vector into celestial coordinates
2122  GVector vector = (*m_rot) * native;
2123 
2124  // Convert vector into sky direction
2125  GSkyDir dir(vector);
2126 
2127  // Setup photon
2128  GPhoton photon(dir, *m_srcEng, *m_srcTime);
2129 
2130  // Evaluate IRF for photon
2131  double irf = m_rsp->irf(*m_event, photon, *m_obs);
2132 
2133  // Return
2134  return irf;
2135 }
2136 
2137 
2138 /***********************************************************************//**
2139  * @brief Zenith angle integration kernel for elliptical model
2140  *
2141  * @param[in] theta Zenith angle (radians).
2142  * @return IRF value.
2143  *
2144  * Integrated the IRF multiplied by the model value for a given zenith angle
2145  * over all azimuth angles.
2146  ***************************************************************************/
2148 {
2149  // Initialise result
2150  double irf = 0.0;
2151 
2152  // Continue only for positive zenith angles (otherwise the integral will
2153  // be zero)
2154  if (theta > 0.0) {
2155 
2156  // Set azimuthal integration range
2157  double phi_min = 0.0;
2158  double phi_max = gammalib::twopi;
2159 
2160  // Setup integration kernel
2161  irf_elliptical_kern_phi integrand(m_rsp,
2162  m_event,
2163  m_obs,
2164  m_model,
2165  m_rot,
2166  theta,
2167  m_srcEng,
2168  m_srcTime);
2169 
2170  // Setup integration
2171  GIntegral integral(&integrand);
2172  integral.fixed_iter(m_iter_phi);
2173 
2174  // Integrate over phi
2175  irf = integral.romberg(phi_min, phi_max, m_iter_phi) * std::sin(theta);
2176 
2177  } // endif: theta was positive
2178 
2179  // Return result
2180  return irf;
2181 }
2182 
2183 
2184 /***********************************************************************//**
2185  * @brief Azimuth angle integration kernel for elliptical model
2186  *
2187  * @param[in] phi Azimuth angle (radians).
2188  * @return IRF value.
2189  *
2190  * Computes the IRF at a given zenith and azimuth angle with respect to a
2191  * specified centre.
2192  ***************************************************************************/
2194 {
2195  // Initialise result
2196  double irf = 0.0;
2197 
2198  // Evaluate sky model
2199  double model = m_model->eval(m_theta, phi, *m_srcEng, *m_srcTime);
2200 
2201  // Continue only if model is positive
2202  if (model > 0.0) {
2203 
2204  // Set up native coordinate vector
2205  double cos_phi = std::cos(phi);
2206  double sin_phi = std::sin(phi);
2207  GVector native(-cos_phi*m_sin_theta, sin_phi*m_sin_theta, m_cos_theta);
2208 
2209  // Rotate vector into celestial coordinates
2210  GVector vector = (*m_rot) * native;
2211 
2212  // Convert vector into sky direction
2213  GSkyDir dir(vector);
2214 
2215  // Setup photon
2216  GPhoton photon(dir, *m_srcEng, *m_srcTime);
2217 
2218  // Evaluate IRF for photon
2219  irf = m_rsp->irf(*m_event, photon, *m_obs) * model;
2220 
2221  } // endif: model was positive
2222 
2223  // Return
2224  return irf;
2225 }
double m_irf_diffuse_resolution
Angular resolution for diffuse model.
Definition: GResponse.hpp:319
double eval(const double &phi)
Azimuth angle integration kernel for radial model.
Definition: GResponse.cpp:2114
Sky map class.
Definition: GSkyMap.hpp:89
int scales(void) const
Return number of scale parameters in model.
Definition: GModel.hpp:247
bool m_use_irf_cache
Control usage of irf cache.
Definition: GResponse.hpp:313
bool m_grad
Gradient flag.
Definition: GResponse.hpp:190
const GObservation * m_obs
Observation.
Definition: GResponse.hpp:184
const double & factor_gradient(void) const
Return parameter factor gradient.
double dec_deg(void) const
Returns Declination in degrees.
Definition: GSkyDir.hpp:256
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:381
Point source spatial model class interface definition.
Energy value class definition.
Abstract elliptical spatial model base class.
const std::string & name(void) const
Return parameter name.
Sparse matrix class interface definition.
virtual std::string instrument(void) const =0
GModelSpectral * spectral(void) const
Return spectral model component.
Definition: GModelSky.hpp:311
GResponse(void)
Void constructor.
Definition: GResponse.cpp:86
const GSkyDir & dir(void) const
Return position of point source.
int size(void) const
Return number of parameters.
bool has_free_pars(void) const
Checks if the spatial model has free parameters.
int size(void) const
Return number of parameters.
const GModelSpatial * component(const int &index) const
Returns pointer to spatial component element.
virtual double irf(const GEvent &event, const GPhoton &photon, const GObservation &obs) const =0
int size(void) const
Return number of energy boundaries.
Definition: GEbounds.hpp:163
void append(const GEnergy &emin, const GEnergy &emax)
Append energy interval.
Definition: GEbounds.cpp:301
GVector cos(const GVector &vector)
Computes cosine of vector elements.
Definition: GVector.cpp:1190
const GSkyMap & map(void) const
Get map.
Generic matrix class definition.
virtual double eval(const GEnergy &srcEng, const GTime &srcTime=GTime(), const bool &gradients=false) const =0
GModelTemporal * temporal(void) const
Return temporal model component.
Definition: GModelSky.hpp:327
Abstract radial spatial model base class interface definition.
virtual void remove_response_cache(const std::string &name)
Remove response cache for model.
Definition: GResponse.cpp:766
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
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 set(const std::string &cache_id, const GVector &vector)
Set cache value.
double sum(const GVector &vector)
Computes vector sum.
Definition: GVector.cpp:944
std::vector< GModelPar * > m_pars
Parameter pointers.
Definition: GResponse.hpp:188
double eval(const double &phi)
Azimuth angle integration kernel for elliptical model.
Definition: GResponse.cpp:2193
int m_irf_elliptical_iter_theta
Elliptical model integration theta iterations.
Definition: GResponse.hpp:317
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
double eval(const double &phi)
Zenith angle integration kernel for radial model.
Definition: GResponse.cpp:2060
Time class.
Definition: GTime.hpp:55
Gammalib tools definition.
GIntegral class interface definition.
Definition: GIntegral.hpp:46
const std::string & name(void) const
Return model name.
Definition: GSource.hpp:114
void init_members(void)
Initialise class members.
Definition: GResponse.cpp:787
const GEvent * m_event
Event.
Definition: GResponse.hpp:186
#define G_CONVOLVE
Definition: GResponse.cpp:61
Sky map class definition.
Source class definition.
virtual GClassCode code(void) const =0
bool is_free(void) const
Signal if parameter is free.
virtual double irf_spatial(const GEvent &event, const GSource &source, const GObservation &obs) const
Return instrument response integrated over the spatial model.
Definition: GResponse.cpp:596
void eulery(const double &angle)
Set Euler rotation matrix around y axis.
Definition: GMatrix.cpp:1157
virtual double theta_max(void) const =0
GResponseVectorCache m_irf_vector_cache
Definition: GResponse.hpp:326
void id(const std::string &id)
Set observation identifier.
Gaussian spectral model class.
GEbounds ebounds_model(const GModelSky &model) const
Return true energy intervals for sky model.
Definition: GResponse.cpp:1965
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
int size(void) const
Return number of parameters in model.
Definition: GModel.hpp:233
double sum_of_scales(void) const
Returns sum of all model scales.
Class that handles photons.
Definition: GPhoton.hpp:47
virtual double irf_elliptical(const GEvent &event, const GSource &source, const GObservation &obs) const
Return instrument response to elliptical source.
Definition: GResponse.cpp:1015
const double & scale(void) const
Return parameter scale.
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:184
virtual double irf_diffuse(const GEvent &event, const GSource &source, const GObservation &obs) const
Return instrument response to diffuse source.
Definition: GResponse.cpp:1118
Model parameter class.
Definition: GModelPar.hpp:87
void remove(const std::string &name)
Remove cache for source.
int m_irf_elliptical_iter_phi
Elliptical model integration phi iterations.
Definition: GResponse.hpp:318
bool has_grad(void) const
Signal if parameter gradient is computed analytically.
double eval(const double &phi)
Zenith angle integration kernel for elliptical model.
Definition: GResponse.cpp:2147
virtual double eval(const GPhoton &photon, const bool &gradients=false) const =0
const GSkyDir & dir(void) const
Return position of radial spatial model.
Energy boundaries container class.
Definition: GEbounds.hpp:60
GSkyDir inx2dir(const int &index) const
Returns sky direction of pixel.
Definition: GSkyMap.cpp:1331
const std::string & name(void) const
Return parameter name.
Definition: GModel.hpp:261
const int & rows(void) const
Return number of matrix rows.
const GEnergy & emin(void) const
Return minimum energy of all intervals.
Definition: GEbounds.hpp:187
void eulerz(const double &angle)
Set Euler rotation matrix around z axis.
Definition: GMatrix.cpp:1194
Integration class for set of functions interface definition.
int size(void) const
Return number of parameters.
const GResponse * m_parent
Response.
Definition: GResponse.hpp:183
virtual const GEnergy & energy(void) const =0
Abstract event base class definition.
GEnergy mean(void) const
Return mean energy.
void remove(const std::string &cache_id)
Remove cache.
double ra_deg(void) const
Returns Right Ascension in degrees.
Definition: GSkyDir.hpp:229
bool m_use_nroi_cache
Control usage of nroi cache.
Definition: GResponse.hpp:314
virtual ~GResponse(void)
Destructor.
Definition: GResponse.cpp:117
const GTime & time(void) const
Return photon arrival time.
Definition: GSource.hpp:156
virtual GEbounds ebounds(const GEnergy &obsEng) const =0
GResponseCache m_nroi_cache
Definition: GResponse.hpp:325
virtual double eval(const GTime &srcTime, const bool &gradients=false) const =0
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
Spatial composite model.
Vector class interface definition.
bool has_scales(void) const
Signals that model has scales.
Definition: GModel.hpp:355
Photon class definition.
Abstract observation base class interface definition.
GModelPar & scale(const int &index)
Returns reference to scale parameter by index.
Definition: GModel.cpp:369
Point source spatial model.
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.
virtual bool use_edisp(void) const =0
Abstract response base class definition.
int m_irf_radial_iter_phi
Radial model integration phi iterations.
Definition: GResponse.hpp:316
virtual GResponse & operator=(const GResponse &rsp)
Assignment operator.
Definition: GResponse.cpp:139
virtual double irf_radial(const GEvent &event, const GSource &source, const GObservation &obs) const
Return instrument response to radial source.
Definition: GResponse.cpp:908
int components(void) const
Return number of model components.
edisp_kerns(const GResponse *parent, const GObservation *obs, const GModelSky *model, const GEvent *event, const GTime &srcTime, const bool &grad)
Constructor for energy dispersion integration kernels class.
Definition: GResponse.cpp:1891
void clear(void)
Clear response vector cache.
Sky model class interface definition.
double eval_prob(const GModelSky &model, const GEvent &event, const GEnergy &srcEng, const GTime &srcTime, const GObservation &obs, const bool &grad) const
Convolve sky model with the instrument response.
Definition: GResponse.cpp:1536
Sky model class.
Definition: GModelSky.hpp:122
virtual GVector column(const int &column) const
Extract column as vector from matrix.
Spatial map cube model class interface definition.
Sparse matrix class definition.
virtual double eval(const GPhoton &photon, const bool &gradients=false) const
Evaluate function.
void copy_members(const GResponse &rsp)
Copy class members.
Definition: GResponse.cpp:813
bool contains(const std::string &cache_id, GVector *irfs=NULL) const
Check if cache contains a value for specific parameters.
GVector sin(const GVector &vector)
Computes sine of vector elements.
Definition: GVector.cpp:1316
virtual GEvents * events(void)
Return events.
int m_irf_radial_iter_theta
Radial model integration theta iterations.
Definition: GResponse.hpp:315
Energy boundaries class interface definition.
int size_edisp_vector(const GModelSky &model, const GObservation &obs, const bool &grad) const
Return size of vector for energy dispersion computation.
Definition: GResponse.cpp:1836
GVector romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegrals.cpp:210
Exception handler interface definition.
virtual double theta_max(void) const =0
int m_size
Array of values and gradients.
Definition: GResponse.hpp:187
const GModelSky * m_model
Sky model.
Definition: GResponse.hpp:185
void fixed_iter(const int &iter)
Set fixed number of iterations.
Definition: GIntegral.hpp:183
GModelSpatial * spatial(void) const
Return spatial model component.
Definition: GModelSky.hpp:295
virtual double irf_ptsrc(const GEvent &event, const GSource &source, const GObservation &obs) const
Return instrument response to point source.
Definition: GResponse.cpp:852
Abstract spatial model base class.
GResponseCache m_irf_cache
Definition: GResponse.hpp:324
Abstract radial spatial model base class.
const double twopi
Definition: GMath.hpp:36
Generic matrix class definition.
Definition: GMatrix.hpp:79
Vector class.
Definition: GVector.hpp:46
const GEnergy & emax(void) const
Return maximum energy of all intervals.
Definition: GEbounds.hpp:199
void clear(void)
Clear response cache.
Abstract diffuse spatial model base class.
double solidangle(const int &index) const
Returns solid angle of pixel.
Definition: GSkyMap.cpp:1858
GVector eval(const double &etrue)
Evaluate energy dispersion integration kernel.
Definition: GResponse.cpp:2012
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
Abstract instrument response base class.
Definition: GResponse.hpp:77
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.
const int & npix(void) const
Returns number of pixels.
Definition: GSkyMap.hpp:345
virtual double convolve(const GModelSky &model, const GEvent &event, const GObservation &obs, const bool &grad=true) const
Convolve sky model with the instrument response.
Definition: GResponse.cpp:186
Integration class interface definition.
Sky direction class.
Definition: GSkyDir.hpp:62
const GSkyMap & cube(void) const
Get map cube.
const int & columns(void) const
Return number of matrix columns.
Time class interface definition.
GTime m_srcTime
True arrival time.
Definition: GResponse.hpp:189
Spatial map model class interface definition.
double scale(const int &index) const
Returns scale of spatial component.
Gaussian spectral model class interface definition.
GVector eval_probs(const GModelSky &model, const GObservation &obs, GMatrixSparse *gradients=NULL) const
Convolve sky model with the instrument response.
Definition: GResponse.cpp:1674
GEnergy sigma(void) const
Return energy width.
Mathematical function definitions.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48
virtual int size(void) const =0
std::string str(const unsigned short int &value)
Convert unsigned short integer value into string.
Definition: GTools.cpp:489