GammaLib  1.7.0.dev
 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-2019 by Juergen Knoedlseder *
5  * ----------------------------------------------------------------------- *
6  * *
7  * This program is free software: you can redistribute it and/or modify *
8  * it under the terms of the GNU General Public License as published by *
9  * the Free Software Foundation, either version 3 of the License, or *
10  * (at your option) any later version. *
11  * *
12  * This program is distributed in the hope that it will be useful, *
13  * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15  * GNU General Public License for more details. *
16  * *
17  * You should have received a copy of the GNU General Public License *
18  * along with this program. If not, see <http://www.gnu.org/licenses/>. *
19  * *
20  ***************************************************************************/
21 /**
22  * @file 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 "GTools.hpp"
34 #include "GMath.hpp"
35 #include "GException.hpp"
36 #include "GIntegral.hpp"
37 #include "GResponse.hpp"
38 #include "GEvent.hpp"
39 #include "GPhoton.hpp"
40 #include "GEnergy.hpp"
41 #include "GTime.hpp"
42 #include "GSource.hpp" // will become obsolete
43 #include "GEbounds.hpp" // will become obsolete
44 #include "GObservation.hpp"
45 #include "GModelSky.hpp"
46 
47 /* __ Method name definitions ____________________________________________ */
48 #define G_IRF_RADIAL "GResponse::irf_radial(GEvent&, GSource&,"\
49  " GObservation&)"
50 #define G_IRF_ELLIPTICAL "GResponse::irf_elliptical(GEvent&, GSource&,"\
51  " GObservation&)"
52 #define G_IRF_DIFFUSE "GResponse::irf_diffuse(GEvent&, GSource&,"\
53  " GObservation&)"
54 
55 /* __ Macros _____________________________________________________________ */
56 
57 /* __ Coding definitions _________________________________________________ */
58 
59 /* __ Debug definitions __________________________________________________ */
60 #define G_DEBUG_CONVOLVE_EDISP //!< Debug convolve for energy dispersion
61 
62 
63 /*==========================================================================
64  = =
65  = Constructors/destructors =
66  = =
67  ==========================================================================*/
68 
69 /***********************************************************************//**
70  * @brief Void constructor
71  ***************************************************************************/
73 {
74  // Initialise class members for clean destruction
75  init_members();
76 
77  // Return
78  return;
79 }
80 
81 
82 /***********************************************************************//**
83  * @brief Copy constructor
84  *
85  * @param[in] rsp Response.
86  ***************************************************************************/
88 {
89  // Initialise class members for clean destruction
90  init_members();
91 
92  // Copy members
93  copy_members(rsp);
94 
95  // Return
96  return;
97 }
98 
99 
100 /***********************************************************************//**
101  * @brief Destructor
102  ***************************************************************************/
104 {
105  // Free members
106  free_members();
107 
108  // Return
109  return;
110 }
111 
112 
113 /*==========================================================================
114  = =
115  = Operators =
116  = =
117  ==========================================================================*/
118 
119 /***********************************************************************//**
120  * @brief Assignment operator
121  *
122  * @param[in] rsp Response.
123  * @return Response.
124  ***************************************************************************/
126 {
127  // Execute only if object is not identical
128  if (this != &rsp) {
129 
130  // Free members
131  free_members();
132 
133  // Initialise private members for clean destruction
134  init_members();
135 
136  // Copy members
137  copy_members(rsp);
138 
139  } // endif: object was not identical
140 
141  // Return this object
142  return *this;
143 }
144 
145 
146 /*==========================================================================
147  = =
148  = Public methods =
149  = =
150  ==========================================================================*/
151 
152 /***********************************************************************//**
153  * @brief Convolve sky model with the instrument response
154  *
155  * @param[in] model Sky model.
156  * @param[in] event Event.
157  * @param[in] obs Observation.
158  * @param[in] grad Should model gradients be computed? (default: true)
159  * @return Event probability.
160  *
161  * Computes the event probability
162  *
163  * \f[
164  * P(p',E',t') = \int \int \int
165  * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp \, dE \, dt
166  * \f]
167  *
168  * without taking into account any time dispersion. Energy dispersion is
169  * correctly handled by this method. If time dispersion is indeed needed,
170  * an instrument specific method needs to be provided.
171  ***************************************************************************/
172 double GResponse::convolve(const GModelSky& model,
173  const GEvent& event,
174  const GObservation& obs,
175  const bool& grad) const
176 {
177  // Set number of iterations for Romberg integration.
178  static const int iter = 6;
179 
180  // Initialise result
181  double prob = 0.0;
182 
183  // Continue only if the model has a spatial component
184  if (model.spatial() != NULL) {
185 
186  // Get source time (no dispersion)
187  GTime srcTime = event.time();
188 
189  // Case A: Integration
190  if (use_edisp()) {
191 
192  // Retrieve true energy boundaries
193  GEbounds ebounds = this->ebounds(event.energy());
194 
195  // Loop over all boundaries
196  for (int i = 0; i < ebounds.size(); ++i) {
197 
198  // Get true energy boundaries in MeV
199  double etrue_min = ebounds.emin(i).MeV();
200  double etrue_max = ebounds.emax(i).MeV();
201 
202  // Continue only if valid
203  if (etrue_max > etrue_min) {
204 
205  // Setup integration function
206  edisp_kern integrand(this, &obs, &model, &event, srcTime, grad);
207  GIntegral integral(&integrand);
208 
209  // Set number of iterations
210  integral.fixed_iter(iter);
211 
212  // Do Romberg integration
213  prob += integral.romberg(etrue_min, etrue_max);
214 
215  } // endif: interval was valid
216 
217  } // endfor: looped over intervals
218 
219  }
220 
221  // Case B: No integration (assume no energy dispersion)
222  else {
223 
224  // Get source energy (no dispersion)
225  GEnergy srcEng = event.energy();
226 
227  // Evaluate probability
228  prob = eval_prob(model, event, srcEng, srcTime, obs, grad);
229 
230  }
231 
232  // Compile option: Check for NaN/Inf
233  #if defined(G_NAN_CHECK)
234  if (gammalib::is_notanumber(prob) || gammalib::is_infinite(prob)) {
235  std::cout << "*** ERROR: GResponse::convolve:";
236  std::cout << " NaN/Inf encountered";
237  std::cout << " (prob=" << prob;
238  std::cout << ", event=" << event;
239  std::cout << ", srcTime=" << srcTime;
240  std::cout << ")" << std::endl;
241  }
242  #endif
243 
244  } // endif: spatial component valid
245 
246  // Return probability
247  return prob;
248 }
249 
250 
251 /*==========================================================================
252  = =
253  = Protected methods =
254  = =
255  ==========================================================================*/
256 
257 /***********************************************************************//**
258  * @brief Initialise class members
259  ***************************************************************************/
261 {
262  // Return
263  return;
264 }
265 
266 
267 /***********************************************************************//**
268  * @brief Copy class members
269  *
270  * @param[in] rsp Response.
271  ***************************************************************************/
273 {
274  // Return
275  return;
276 }
277 
278 
279 /***********************************************************************//**
280  * @brief Delete class members
281  ***************************************************************************/
283 {
284  // Return
285  return;
286 }
287 
288 
289 /***********************************************************************//**
290  * @brief Convolve sky model with the instrument response
291  *
292  * @param[in] model Sky model.
293  * @param[in] event Event.
294  * @param[in] srcEng Source energy.
295  * @param[in] srcTime Source time.
296  * @param[in] obs Observation.
297  * @param[in] grad Should model gradients be computed? (default: true)
298  * @return Event probability.
299  *
300  * Computes the event probability
301  *
302  * \f[
303  * P(p',E',t'|E,t) = \int S(p,E,t) \times R(p',E',t'|p,E,t) \, dp
304  * \f]
305  *
306  * for a given true energy \f$E\f$ and time \f$t\f$.
307  ***************************************************************************/
308 double GResponse::eval_prob(const GModelSky& model,
309  const GEvent& event,
310  const GEnergy& srcEng,
311  const GTime& srcTime,
312  const GObservation& obs,
313  const bool& grad) const
314 {
315  // Initialise result
316  double prob = 0.0;
317 
318  // Continue only if the model has a spatial component
319  if (model.spatial() != NULL) {
320 
321  // Set source
322  GSource source(model.name(), model.spatial(), srcEng, srcTime);
323 
324  // Compute IRF value
325  double irf = this->irf(event, source, obs);
326 
327  // Continue only if IRF value is positive
328  if (irf > 0.0) {
329 
330  // Optionally get scaling
331  double scale = (model.has_scales())
332  ? model.scale(obs.instrument()).value() : 1.0;
333 
334  // Evaluate spectral and temporal components
335  double spec = (model.spectral() != NULL)
336  ? model.spectral()->eval(srcEng, srcTime, grad) : 1.0;
337  double temp = (model.temporal() != NULL)
338  ? model.temporal()->eval(srcTime, grad) : 1.0;
339 
340  // Compute probability
341  prob = spec * temp * irf * scale;
342 
343  // Optionally compute partial derivatives
344  if (grad) {
345 
346  // Multiply factors to spectral gradients
347  if (model.spectral() != NULL) {
348  double fact = temp * irf * scale;
349  if (fact != 1.0) {
350  for (int i = 0; i < model.spectral()->size(); ++i) {
351  (*model.spectral())[i].factor_gradient((*model.spectral())[i].factor_gradient() * fact);
352  }
353  }
354  }
355 
356  // Multiply factors to temporal gradients
357  if (model.temporal() != NULL) {
358  double fact = spec * irf * scale;
359  if (fact != 1.0) {
360  for (int i = 0; i < model.temporal()->size(); ++i) {
361  (*model.temporal())[i].factor_gradient((*model.temporal())[i].factor_gradient() * fact);
362  }
363  }
364  }
365 
366  // Optionally compute scale gradient for instrument
367  if (model.has_scales()) {
368  for (int i = 0; i < model.scales(); ++i) {
369  double g_scale = 0.0;
370  if (model.scale(i).name() == obs.instrument()) {
371  if (model.scale(i).is_free()) {
372  g_scale = model.scale(i).scale() * spec * temp * irf;
373  }
374  }
375  model.scale(i).factor_gradient(g_scale);
376  }
377  }
378 
379  } // endif: computed partial derivatives
380 
381  // Compile option: Check for NaN/Inf
382  #if defined(G_NAN_CHECK)
383  if (gammalib::is_notanumber(prob) || gammalib::is_infinite(prob)) {
384  std::cout << "*** ERROR: GResponse::eval_prob:";
385  std::cout << " NaN/Inf encountered";
386  std::cout << " (prob=" << prob;
387  std::cout << ", spec=" << spec;
388  std::cout << ", temp=" << temp;
389  std::cout << ", irf=" << irf;
390  std::cout << ")" << std::endl;
391  }
392  #endif
393 
394  } // endif: IRF value was positive
395 
396  // ... otherwise if gradient computation is requested then set the
397  // spectral and temporal gradients to zero
398  else if (grad) {
399 
400  // Reset spectral gradients
401  if (model.spectral() != NULL) {
402  for (int i = 0; i < model.spectral()->size(); ++i) {
403  (*model.spectral())[i].factor_gradient(0.0);
404  }
405  }
406 
407  // Reset temporal gradients
408  if (model.temporal() != NULL) {
409  for (int i = 0; i < model.temporal()->size(); ++i) {
410  (*model.temporal())[i].factor_gradient(0.0);
411  }
412  }
413 
414  // Optionally reset scale gradients
415  if (model.has_scales()) {
416  for (int i = 0; i < model.scales(); ++i) {
417  model.scale(i).factor_gradient(0.0);
418  }
419  }
420 
421  } // endelse: IRF value was not positive
422 
423  } // endif: Gamma-ray source model had a spatial component
424 
425  // Return event probability
426  return prob;
427 }
428 
429 
430 /***********************************************************************//**
431  * @brief Integration kernel for GResponse::edisp_kern() class
432  *
433  * @param[in] etrue True photon energy in MeV.
434  *
435  * This method implements the integration kernel needed for the
436  * GResponse::edisp_kern() class.
437  ***************************************************************************/
438 double GResponse::edisp_kern::eval(const double& etrue)
439 {
440  // Set true energy
441  GEnergy srcEng;
442  srcEng.MeV(etrue);
443 
444  // Get function value
445  double value = m_parent->eval_prob(*m_model, *m_event,
446  srcEng, m_srcTime, *m_obs, m_grad);
447 
448  // Compile option: Check for NaN
449  #if defined(G_NAN_CHECK)
450  if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
451  std::cout << "*** ERROR: GResponse::edisp_kern::eval";
452  std::cout << "(etrue=" << etrue << "): ";
453  std::cout << " NaN/Inf encountered";
454  std::cout << " (value=" << value;
455  std::cout << ")" << std::endl;
456  }
457  #endif
458 
459  // Return value
460  return value;
461 }
int scales(void) const
Return number of scale parameters in model.
Definition: GModel.hpp:231
const double & factor_gradient(void) const
Return parameter gradient factor.
double romberg(std::vector< double > bounds, const int &order=5)
Perform Romberg integration.
Definition: GIntegral.cpp:380
Energy value class definition.
const std::string & name(void) const
Return parameter name.
virtual std::string instrument(void) const =0
GModelSpectral * spectral(void) const
Return spectral model component.
Definition: GModelSky.hpp:262
GResponse(void)
Void constructor.
Definition: GResponse.cpp:72
int size(void) const
Return number of parameters.
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:157
const GResponse * m_parent
Response.
Definition: GResponse.hpp:132
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:278
Abstract interface for the event classes.
Definition: GEvent.hpp:71
bool m_grad
Gradient flag.
Definition: GResponse.hpp:137
double MeV(void) const
Return energy in MeV.
Definition: GEnergy.cpp:321
Time class.
Definition: GTime.hpp:54
Gammalib tools definition.
GIntegral class interface definition.
Definition: GIntegral.hpp:45
void init_members(void)
Initialise class members.
Definition: GResponse.cpp:260
Source class definition.
bool is_free(void) const
Signal if parameter is free.
bool is_notanumber(const double &x)
Signal if argument is not a number.
Definition: GTools.hpp:184
const double & scale(void) const
Return parameter scale.
bool is_infinite(const double &x)
Signal if argument is infinite.
Definition: GTools.hpp:167
Energy boundaries container class.
Definition: GEbounds.hpp:60
const std::string & name(void) const
Return parameter name.
Definition: GModel.hpp:245
const GModelSky * m_model
Sky model.
Definition: GResponse.hpp:134
const GEnergy & emin(void) const
Return minimum energy of all intervals.
Definition: GEbounds.hpp:181
int size(void) const
Return number of parameters.
const GEvent * m_event
Event.
Definition: GResponse.hpp:135
virtual const GEnergy & energy(void) const =0
Abstract event base class definition.
virtual ~GResponse(void)
Destructor.
Definition: GResponse.cpp:103
virtual GEbounds ebounds(const GEnergy &obsEng) const =0
virtual double eval(const GTime &srcTime, const bool &gradients=false) const =0
Abstract observation base class.
bool has_scales(void) const
Signals that model has scales.
Definition: GModel.hpp:338
Photon class definition.
Abstract observation base class interface definition.
GModelPar & scale(const int &index)
Returns reference to scale parameter by index.
Definition: GModel.cpp:363
void free_members(void)
Delete class members.
Definition: GResponse.cpp:282
virtual bool use_edisp(void) const =0
Abstract response base class definition.
virtual GResponse & operator=(const GResponse &rsp)
Assignment operator.
Definition: GResponse.cpp:125
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:308
Sky model class.
Definition: GModelSky.hpp:120
const GObservation * m_obs
Observation.
Definition: GResponse.hpp:133
GTime m_srcTime
True arrival time.
Definition: GResponse.hpp:136
void copy_members(const GResponse &rsp)
Copy class members.
Definition: GResponse.cpp:272
Energy boundaries class interface definition.
Exception handler interface definition.
void fixed_iter(const int &iter)
Set fixed number of iterations.
Definition: GIntegral.hpp:181
GModelSpatial * spatial(void) const
Return spatial model component.
Definition: GModelSky.hpp:246
const GEnergy & emax(void) const
Return maximum energy of all intervals.
Definition: GEbounds.hpp:193
Abstract instrument response base class.
Definition: GResponse.hpp:67
Class that handles gamma-ray sources.
Definition: GSource.hpp:53
double eval(const double &etrue)
Integration kernel for GResponse::edisp_kern() class.
Definition: GResponse.cpp:438
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:172
Integration class interface definition.
Time class interface definition.
Mathematical function definitions.
Class that handles energies in a unit independent way.
Definition: GEnergy.hpp:48