Sergio Hernández Cadena / gammalib

Commit 1125319657da763967bb935bce560e85e844b4cf

Authored by Eirik Hatlen
Committed by Jürgen Knödlseder
1 parent 4e0a5d23

Included model for ebounds in GResponse (3386) 

Only relevant for GSpectralGaussian when colvoling with the energy dispersion.
The ebounds will just be one bin allowing for proper convolving of response when
using a narrow (<< eres) Gaussian
Inline Side-by-side
Showing 2 changed files with 113 additions and 969 deletions
include/GResponse.hpp
  View file @1125319
... ... @@ -33,13 +33,9 @@
33 33 #include "GFunction.hpp"
34 34 #include "GTime.hpp"
35 35 #include "GResponseCache.hpp"
36   -#include "GResponseVectorCache.hpp"
37   -#include "GFunctions.hpp"
38 36  
39 37 /* __ Forward declarations _______________________________________________ */
40   -class GVector;
41 38 class GMatrix;
42   -class GMatrixSparse;
43 39 class GEvent;
44 40 class GPhoton;
45 41 class GSource;
... ... @@ -47,7 +43,6 @@ class GSkyDir;
47 43 class GEnergy;
48 44 class GEbounds;
49 45 class GObservation;
50   -class GModelPar;
51 46 class GModelSky;
52 47 class GModelSpatialRadial;
53 48 class GModelSpatialElliptical;
... ... @@ -106,34 +101,24 @@ public:
106 101 const GEvent& event,
107 102 const GObservation& obs,
108 103 const bool& grad = true) const;
109   - virtual GVector convolve(const GModelSky& model,
110   - const GObservation& obs,
111   - GMatrixSparse* gradients = NULL) const;
112 104 virtual double irf_spatial(const GEvent& event,
113 105 const GSource& source,
114 106 const GObservation& obs) const;
115   - virtual GVector irf_spatial(const GModelSky& model,
116   - const GObservation& obs,
117   - GMatrix* gradients = NULL) const;
118 107 virtual void remove_response_cache(const std::string& name);
119 108  
  109 + virtual GEbounds ebounds_model(const GModelSky& model) const;
  110 +
120 111 protected:
121 112 // Protected methods
122   - void init_members(void);
123   - void copy_members(const GResponse& rsp);
124   - void free_members(void);
125   - double eval_prob(const GModelSky& model,
126   - const GEvent& event,
127   - const GEnergy& srcEng,
128   - const GTime& srcTime,
129   - const GObservation& obs,
130   - const bool& grad) const;
131   - GVector eval_probs(const GModelSky& model,
132   - const GObservation& obs,
133   - GMatrixSparse* gradients = NULL) const;
134   - int size_edisp_vector(const GModelSky& model,
135   - const GObservation& obs,
136   - const bool& grad) const;
  113 + void init_members(void);
  114 + void copy_members(const GResponse& rsp);
  115 + void free_members(void);
  116 + double eval_prob(const GModelSky& model,
  117 + const GEvent& event,
  118 + const GEnergy& srcEng,
  119 + const GTime& srcTime,
  120 + const GObservation& obs,
  121 + const bool& grad) const;
137 122  
138 123 // Virtual protected methods
139 124 virtual double irf_ptsrc(const GEvent& event,
... ... @@ -151,42 +136,30 @@ protected:
151 136 virtual double irf_composite(const GEvent& event,
152 137 const GSource& source,
153 138 const GObservation& obs) const;
154   - virtual GVector irf_ptsrc(const GModelSky& model,
155   - const GObservation& obs,
156   - GMatrix* gradients = NULL) const;
157   - virtual GVector irf_radial(const GModelSky& model,
158   - const GObservation& obs,
159   - GMatrix* gradients = NULL) const;
160   - virtual GVector irf_elliptical(const GModelSky& model,
161   - const GObservation& obs,
162   - GMatrix* gradients = NULL) const;
163   - virtual GVector irf_diffuse(const GModelSky& model,
164   - const GObservation& obs,
165   - GMatrix* gradients = NULL) const;
166   - virtual GVector irf_composite(const GModelSky& model,
167   - const GObservation& obs,
168   - GMatrix* gradients = NULL) const;
169 139  
170 140 // Protected classes
171   - class edisp_kerns : public GFunctions {
  141 + class edisp_kern : public GFunction {
172 142 public:
173   - edisp_kerns(const GResponse* parent,
174   - const GObservation* obs,
175   - const GModelSky* model,
176   - const GEvent* event,
177   - const GTime& srcTime,
178   - const bool& grad);
179   - int size(void) const { return m_size; }
180   - GVector eval(const double& etrue);
  143 + edisp_kern(const GResponse* parent,
  144 + const GObservation* obs,
  145 + const GModelSky* model,
  146 + const GEvent* event,
  147 + const GTime& srcTime,
  148 + const bool& grad) :
  149 + m_parent(parent),
  150 + m_obs(obs),
  151 + m_model(model),
  152 + m_event(event),
  153 + m_srcTime(srcTime),
  154 + m_grad(grad) { }
  155 + double eval(const double& etrue);
181 156 protected:
182   - const GResponse* m_parent; //!< Response
183   - const GObservation* m_obs; //!< Observation
184   - const GModelSky* m_model; //!< Sky model
185   - const GEvent* m_event; //!< Event
186   - int m_size; //!< Array of values and gradients
187   - std::vector<GModelPar*> m_pars; //!< Parameter pointers
188   - GTime m_srcTime; //!< True arrival time
189   - bool m_grad; //!< Gradient flag
  157 + const GResponse* m_parent; //!< Response
  158 + const GObservation* m_obs; //!< Observation
  159 + const GModelSky* m_model; //!< Sky model
  160 + const GEvent* m_event; //!< Event
  161 + GTime m_srcTime; //!< True arrival time
  162 + bool m_grad; //!< Gradient flag
190 163 };
191 164 class irf_radial_kern_theta : public GFunction {
192 165 public:
... ... @@ -320,9 +293,8 @@ protected:
320 293 // Cache for irf(GEvent&, GSource&, GObservation&) and
321 294 // nroi(GModelSky&, GEnergy&, GTime&, GEnergy&, GTime&, GObservation&)
322 295 // computations
323   - mutable GResponseCache m_irf_cache;
324   - mutable GResponseCache m_nroi_cache;
325   - mutable GResponseVectorCache m_irf_vector_cache;
  296 + mutable GResponseCache m_irf_cache;
  297 + mutable GResponseCache m_nroi_cache;
326 298 };
327 299  
328 300 #endif /* GRESPONSE_HPP */
... ...
src/obs/GResponse.cpp
  View file @1125319
1 1 /***************************************************************************
2 2 * GResponse.cpp - Abstract response base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2008-2021 by Juergen Knoedlseder *
  4 + * copyright (C) 2008-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -35,9 +35,7 @@
35 35 #include "GMath.hpp"
36 36 #include "GVector.hpp"
37 37 #include "GMatrix.hpp"
38   -#include "GMatrixSparse.hpp"
39 38 #include "GIntegral.hpp"
40   -#include "GIntegrals.hpp"
41 39 #include "GResponse.hpp"
42 40 #include "GEvent.hpp"
43 41 #include "GPhoton.hpp"
... ... @@ -48,6 +46,7 @@
48 46 #include "GObservation.hpp"
49 47 #include "GSkyMap.hpp"
50 48 #include "GModelSky.hpp"
  49 +#include "GModelSpectralGauss.hpp"
51 50 #include "GModelSpatialPointSource.hpp"
52 51 #include "GModelSpatialRadial.hpp"
53 52 #include "GModelSpatialElliptical.hpp"
... ... @@ -57,8 +56,6 @@
57 56 #include "GModelSpatialComposite.hpp"
58 57  
59 58 /* __ Method name definitions ____________________________________________ */
60   -#define G_CONVOLVE "GResponse::convolve(GModelSky&, GObservation&, "\
61   - "GMatrixSparse*)"
62 59 #define G_IRF_RADIAL "GResponse::irf_radial(GEvent&, GSource&,"\
63 60 " GObservation&)"
64 61 #define G_IRF_ELLIPTICAL "GResponse::irf_elliptical(GEvent&, GSource&,"\
... ... @@ -71,6 +68,7 @@
71 68 /* __ Coding definitions _________________________________________________ */
72 69  
73 70 /* __ Debug definitions __________________________________________________ */
  71 +#define G_DEBUG_CONVOLVE_EDISP //!< Debug convolve for energy dispersion
74 72  
75 73  
76 74 /*==========================================================================
... ... @@ -168,7 +166,7 @@ GResponse&amp; GResponse::operator=(const GResponse&amp; rsp)
168 166 * @param[in] model Sky model.
169 167 * @param[in] event Event.
170 168 * @param[in] obs Observation.
171   - * @param[in] grad Should model gradients be computed?
  169 + * @param[in] grad Should model gradients be computed? (default: true)
172 170 * @return Event probability.
173 171 *
174 172 * Computes the event probability
... ... @@ -205,11 +203,11 @@ double GResponse::convolve(const GModelSky&amp; model,
205 203 // Retrieve true energy boundaries
206 204 GEbounds ebounds = this->ebounds(event.energy());
207 205  
208   - // Determine size of result array
209   - int size = size_edisp_vector(model, obs, grad);
210   -
211 206 // Initialise vector array
212   - GVector array(size);
  207 + //GVector array(ebounds.size()); // need this?
  208 +
  209 + // Get Valid Ebounds
  210 + GEbounds valid_ebounds = this->ebounds_model(model);
213 211  
214 212 // Loop over all boundaries
215 213 for (int i = 0; i < ebounds.size(); ++i) {
... ... @@ -218,67 +216,38 @@ double GResponse::convolve(const GModelSky&amp; model,
218 216 double etrue_min = ebounds.emin(i).MeV();
219 217 double etrue_max = ebounds.emax(i).MeV();
220 218  
  219 +
  220 + // Limit energy bounds for model
  221 + if (valid_ebounds.size() > 0 )
  222 + {
  223 + if (valid_ebounds.emin(0).MeV() > etrue_min )
  224 + {
  225 + etrue_min = valid_ebounds.emin(0).MeV();
  226 + }
  227 + if (valid_ebounds.emax(0).MeV() < etrue_max )
  228 + {
  229 + etrue_max = valid_ebounds.emax(0).MeV();
  230 + }
  231 + }
  232 +
221 233 // Continue only if valid
222 234 if (etrue_max > etrue_min) {
223 235  
224 236 // Setup integration function
225   - edisp_kerns integrand(this, &obs, &model, &event, srcTime, grad);
226   - GIntegrals integral(&integrand);
  237 + edisp_kern integrand(this, &obs, &model, &event, srcTime, grad);
  238 + GIntegral integral(&integrand);
227 239  
228 240 // Set number of iterations
229 241 integral.fixed_iter(iter);
230 242  
231 243 // Do Romberg integration
232   - array += integral.romberg(etrue_min, etrue_max);
  244 + prob += integral.romberg(etrue_min, etrue_max);
233 245  
234 246 } // endif: interval was valid
235 247  
236 248 } // endfor: looped over intervals
237 249  
238   - // Initialise array index
239   - int index = 0;
240   -
241   - // Get probability
242   - prob = array[index++];
243   -
244   - // Set gradients
245   - if (grad) {
246   -
247   - // Set spectral gradients
248   - if (model.spectral() != NULL) {
249   - for (int i = 0; i < model.spectral()->size(); ++i) {
250   - GModelPar& par = (*(model.spectral()))[i];
251   - if (par.is_free() && par.has_grad()) {
252   - par.factor_gradient(array[index++]);
253   - }
254   - }
255   - }
256   -
257   - // Set temporal gradients
258   - if (model.temporal() != NULL) {
259   - for (int i = 0; i < model.temporal()->size(); ++i) {
260   - GModelPar& par = (*(model.temporal()))[i];
261   - if (par.is_free() && par.has_grad()) {
262   - par.factor_gradient(array[index++]);
263   - }
264   - }
265   - }
266   -
267   - // Optionally set scale gradient for instrument
268   - if (model.has_scales()) {
269   - for (int i = 0; i < model.scales(); ++i) {
270   - const GModelPar& par = model.scale(i);
271   - if (par.name() == obs.instrument()) {
272   - if (par.is_free() && par.has_grad()) {
273   - par.factor_gradient(array[index++]);
274   - }
275   - }
276   - }
277   - }
278   -
279   - } // endif: set gradients
280   -
281   - } // endif: Case A
  250 + }
282 251  
283 252 // Case B: No integration (assume no energy dispersion)
284 253 else {
... ... @@ -289,7 +258,7 @@ double GResponse::convolve(const GModelSky&amp; model,
289 258 // Evaluate probability
290 259 prob = eval_prob(model, event, srcEng, srcTime, obs, grad);
291 260  
292   - } // endelse: Case B
  261 + }
293 262  
294 263 // Compile option: Check for NaN/Inf
295 264 #if defined(G_NAN_CHECK)
... ... @@ -310,214 +279,51 @@ double GResponse::convolve(const GModelSky&amp; model,
310 279 }
311 280  
312 281  
  282 +
313 283 /***********************************************************************//**
314   - * @brief Convolve sky model with the instrument response
  284 + * @brief Validity of true energy interval
315 285 *
316   - * @param[in] model Sky model.
317   - * @param[in] obs Observation.
318   - * @param[out] gradients Pointer to matrix of gradients.
319   - * @return Vector of event probabilities.
  286 + * @param[in] GModelSky& model
  287 + * @return GEbounds
320 288 *
321   - * Computes the event probability
  289 + * Computes the valitdity of the true energy model.
  290 + * Special treatment of narrow Gaussians as approximation for
  291 + * monochromatic lines.
  292 + * For a Gaussian Model return one sigle boundry 5 sigma around the line
  293 + * energy.
  294 + * Other models return empty boundry container.
322 295 *
323   - * \f[
324   - * P(p',E',t') = \int \int \int
325   - * S(p,E,t) \times R(p',E',t'|p,E,t) \, dp \, dE \, dt
326   - * \f]
327   - *
328   - * without taking into account any time dispersion. Energy dispersion is
329   - * correctly handled by this method. If time dispersion is indeed needed,
330   - * an instrument specific method needs to be provided.
331 296 ***************************************************************************/
332   -GVector GResponse::convolve(const GModelSky& model,
333   - const GObservation& obs,
334   - GMatrixSparse* gradients) const
  297 +GEbounds GResponse::ebounds_model(const GModelSky& model) const
335 298 {
336   - // Set number of iterations for Romberg integration.
337   - static const int iter = 6;
338   -
339   - // Get number of model parameters and number of events
340   - int npars = model.size();
341   - int nevents = obs.events()->size();
342   -
343   - // Initialise gradients flag
344   - bool grad = (gradients != NULL);
345   -
346   - // Check matrix consistency
347   - if (grad) {
348   - if (gradients->columns() != npars) {
349   - std::string msg = "Number of "+gammalib::str(gradients->columns())+
350   - " columns in gradient matrix differs from number "
351   - "of "+gammalib::str(npars)+" parameters "
352   - "in model. Please provide a compatible gradient "
353   - "matrix.";
354   - throw GException::invalid_argument(G_CONVOLVE, msg);
355   - }
356   - if (gradients->rows() != nevents) {
357   - std::string msg = "Number of "+gammalib::str(gradients->rows())+
358   - " rows in gradient matrix differs from number "
359   - "of "+gammalib::str(nevents)+" events in "
360   - "observation. Please provide a compatible "
361   - "gradient matrix.";
362   - throw GException::invalid_argument(G_CONVOLVE, msg);
363   - }
364   - }
365   -
366   - // Initialise result
367   - GVector probs(nevents);
368   -
369   - // Continue only if the model has a spatial component
370   - if (model.spatial() != NULL) {
371   -
372   - // Case A: Integration
373   - if (use_edisp()) {
374   -
375   - // Initialise temporary vectors to hold gradients
376   - GVector* tmp_gradients = NULL;
377   - if (grad) {
378   - tmp_gradients = new GVector[npars];
379   - for (int i = 0; i < npars; ++i) {
380   - tmp_gradients[i] = GVector(nevents);
381   - }
382   - }
383   -
384   - // Determine size of result array
385   - int size = size_edisp_vector(model, obs, grad);
386   -
387   - // Loop over events
388   - for (int k = 0; k < nevents; ++k) {
389   -
390   - // Get reference to event
391   - const GEvent& event = *((*obs.events())[k]);
392   -
393   - // Get source time (no dispersion)
394   - GTime srcTime = event.time();
395   -
396   - // Retrieve true energy boundaries
397   - GEbounds ebounds = this->ebounds(event.energy());
398   -
399   - // Initialise vector array
400   - GVector array(size);
401   -
402   - // Loop over all boundaries
403   - for (int i = 0; i < ebounds.size(); ++i) {
404   -
405   - // Get true energy boundaries in MeV
406   - double etrue_min = ebounds.emin(i).MeV();
407   - double etrue_max = ebounds.emax(i).MeV();
408   -
409   - // Continue only if valid
410   - if (etrue_max > etrue_min) {
411   -
412   - // Setup integration function
413   - edisp_kerns integrand(this, &obs, &model, &event, srcTime, grad);
414   - GIntegrals integral(&integrand);
415   -
416   - // Set number of iterations
417   - integral.fixed_iter(iter);
418   -
419   - // Do Romberg integration
420   - array += integral.romberg(etrue_min, etrue_max);
421   -
422   - } // endif: interval was valid
423   -
424   - } // endfor: looped over intervals
425   -
426   - // Initialise array index
427   - int index = 0;
428   -
429   - // Get probability
430   - probs[k] = array[index++];
431   -
432   - // Set gradients
433   - if (grad) {
434   -
435   - // Get number of spatial, spectral and temporal parameters
436   - int n_spat = model.spatial()->size();
437   - int n_spec = model.spectral()->size();
438   - int n_temp = model.temporal()->size();
439   -
440   - // Set spectral gradients
441   - if (model.spectral() != NULL) {
442   - int offset = n_spat;
443   - for (int i = 0; i < n_spec; ++i) {
444   - const GModelPar& par = (*(model.spectral()))[i];
445   - if (par.is_free() && par.has_grad()) {
446   - tmp_gradients[offset+i][k] = array[index++];
447   - }
448   - }
449   - }
450   -
451   - // Set temporal gradients
452   - if (model.temporal() != NULL) {
453   - int offset = n_spat + n_spec;
454   - for (int i = 0; i < n_temp; ++i) {
455   - const GModelPar& par = (*(model.temporal()))[i];
456   - if (par.is_free() && par.has_grad()) {
457   - tmp_gradients[offset+i][k] = array[index++];
458   - }
459   - }
460   - }
461   -
462   - // Optionally set scale gradient for instrument
463   - if (model.has_scales()) {
464   - int offset = n_spat + n_spec + n_temp;
465   - for (int i = 0; i < model.scales(); ++i) {
466   - const GModelPar& par = model.scale(i);
467   - if (par.name() == obs.instrument()) {
468   - if (par.is_free() && par.has_grad()) {
469   - tmp_gradients[offset+i][k] = array[index++];
470   - }
471   - }
472   - }
473   - }
474   -
475   - } // endif: set gradients
476   -
477   - // Compile option: Check for NaN/Inf
478   - #if defined(G_NAN_CHECK)
479   - if (gammalib::is_notanumber(probs[k]) || gammalib::is_infinite(probs[k])) {
480   - std::cout << "*** ERROR: GResponse::convolve:";
481   - std::cout << " NaN/Inf encountered";
482   - std::cout << " (probs[" << k << "]=" << probs[k];
483   - std::cout << ", event=" << event;
484   - std::cout << ", srcTime=" << srcTime;
485   - std::cout << ")" << std::endl;
486   - }
487   - #endif
488   -
489   - } // endfor: looped over events
490   -
491   - // Post-process gradients
492   - if (grad) {
493 299  
494   - // Fill gradients into matrix
495   - for (int i = 0; i < npars; ++i) {
496   - gradients->column(i, tmp_gradients[i]);
497   - }
498   -
499   - // Delete temporal gradients
500   - delete [] tmp_gradients;
501   -
502   - } // endif: post-processed gradients
503   -
504   - } // endif: Case A
505   -
506   - // Case B: No integration (assume no energy dispersion)
507   - else {
508   -
509   - // Evaluate probability
510   - probs = eval_probs(model, obs, gradients);
511   -
512   - } // endelse: Case B
513   -
514   - } // endif: spatial component valid
515   -
516   - // Return probabilities
517   - return probs;
  300 + if(model.spectral()->type() == "Gaussian")
  301 + {
  302 + const double size = 5.0;
  303 + const GModelSpectralGauss* gauss = dynamic_cast<const GModelSpectralGauss*>(model.spectral());
  304 +
  305 +
  306 + GEnergy mu = gauss->mean();
  307 + GEnergy sigma = gauss->sigma();
  308 +
  309 + GEnergy ebound_min = mu - size * sigma;
  310 + GEnergy ebound_max = mu + size * sigma;
  311 +
  312 + GEbounds emebounds(ebound_min, ebound_max);
  313 +
  314 + return emebounds;
  315 +
  316 + }
  317 + else
  318 + {
  319 + GEbounds empty_ebounds;
  320 + return empty_ebounds;
  321 + }
  322 +
518 323 }
519 324  
520 325  
  326 +
521 327 /***********************************************************************//**
522 328 * @brief Return instrument response integrated over the spatial model
523 329 *
... ... @@ -569,7 +375,7 @@ double GResponse::irf_spatial(const GEvent&amp; event,
569 375 double irf = 0.0;
570 376  
571 377 // Set IRF value attributes
572   - std::string name = obs.id() + ":" + source.name();
  378 + std::string name = obs.id() + "::" + source.name();
573 379 const GInstDir& dir = event.dir();
574 380 const GEnergy& ereco = event.energy();
575 381 const GEnergy& etrue = source.energy();
... ... @@ -619,112 +425,6 @@ double GResponse::irf_spatial(const GEvent&amp; event,
619 425  
620 426  
621 427 /***********************************************************************//**
622   - * @brief Return instrument response vector integrated over the spatial model
623   - *
624   - * @param[in] model Sky model.
625   - * @param[in] obs Observation.
626   - * @param[in] gradients Gradients matrix.
627   - * @return Instrument response vector to a spatial model.
628   - *
629   - * Returns the instrument response to a sky model integrated over the spatial
630   - * model component for all events in a given observation. The method computes
631   - *
632   - * \f[
633   - * {\tt irf}(p', E', t') = \int_p M_{\rm S}(p | E, t) \,
634   - * R(p', E', t' | p, E, t) \, d\,p
635   - * \f]
636   - *
637   - * where
638   - * * \f$M_{\rm S}(p | E, t)\f$ is the spatial model component,
639   - * * \f$R(p', E', t' | p, E, t)\f$ is the Instrument Response Function (IRF),
640   - * * \f$p'\f$ is the measured instrument direction,
641   - * * \f$E'\f$ is the measured or reconstructed energy,
642   - * * \f$t'\f$ is the measured arrival time,
643   - * * \f$p\f$ is the true photon arrival direction,
644   - * * \f$E\f$ is the true photon energy, and
645   - * * \f$t\f$ is the true trigger time.
646   - *
647   - * The integration is done over all relevant true sky directions \f$p\f$.
648   - *
649   - * Depending on the type of the source model the method branches to the
650   - * following methods to perform the actual computations
651   - *
652   - * * irf_ptsrc() - for the handling of a point source
653   - * * irf_radial() - for radial models
654   - * * irf_elliptical() - for elliptical models
655   - * * irf_diffuse() - for diffuse models
656   - * * irf_composite() - for composite models
657   - *
658   - * The method implements a caching mechanism for spatial models that have all
659   - * parameters fixed. For those models the instrument response for a given
660   - * event and observation is only computed once and then stored in an internal
661   - * cache from which it is fetched back in case that the method is called
662   - * again for the same event and observation.
663   - ***************************************************************************/
664   -GVector GResponse::irf_spatial(const GModelSky& model,
665   - const GObservation& obs,
666   - GMatrix* gradients) const
667   -{
668   - // Get number of model parameters and number of events
669   - int npars = model.size();
670   - int nevents = obs.events()->size();
671   -
672   - // Initialise result
673   - GVector irfs(nevents);
674   -
675   - // Continue only if model has a spatial component
676   - if (model.spatial() != NULL) {
677   -
678   - // Set IRF cache identifier
679   - std::string cache_id = obs.id() + ":" + model.name();
680   -
681   - // Signal if spatial model has free parameters
682   - bool has_free_pars = model.spatial()->has_free_pars();
683   -
684   - // If the spatial model component has free parameters, or the response
685   - // cache should not be used, or the cache does not contain the requested
686   - // IRF value then compute the IRF value for the spatial model.
687   - if (has_free_pars ||
688   - !m_use_irf_cache ||
689   - !m_irf_vector_cache.contains(cache_id, &irfs)) {
690   -
691   - // Compute IRF for spatial model
692   - switch (model.spatial()->code()) {
693   - case GMODEL_SPATIAL_POINT_SOURCE:
694   - irfs = irf_ptsrc(model, obs, gradients);
695   - break;
696   - case GMODEL_SPATIAL_RADIAL:
697   - irfs = irf_radial(model, obs, gradients);
698   - break;
699   - case GMODEL_SPATIAL_ELLIPTICAL:
700   - irfs = irf_elliptical(model, obs, gradients);
701   - break;
702   - case GMODEL_SPATIAL_DIFFUSE:
703   - irfs = irf_diffuse(model, obs, gradients);
704   - break;
705   - case GMODEL_SPATIAL_COMPOSITE:
706   - irfs = irf_composite(model, obs, gradients);
707   - break;
708   - default:
709   - break;
710   - }
711   -
712   - // If the spatial model has no free parameters and the response cache
713   - // should be used then put the IRF value in the response cache.
714   - if (!has_free_pars && m_use_irf_cache) {
715   - m_irf_vector_cache.set(cache_id, irfs);
716   - }
717   -
718   - } // endif: computed spatial model
719   -
720   - } // endif: model had spatial component
721   -
722   - // Return IRF values
723   - return irfs;
724   -}
725   -
726   -
727   -/***********************************************************************//**
728 428 * @brief Remove response cache for model
729 429 *
730 430 * @param[in] name Model name.
... ... @@ -736,7 +436,6 @@ void GResponse::remove_response_cache(const std::string&amp; name)
736 436 // Remove model from response caches
737 437 m_irf_cache.remove(name);
738 438 m_nroi_cache.remove(name);
739   - m_irf_vector_cache.remove(name);
740 439  
741 440 // Return
742 441 return;
... ... @@ -766,7 +465,6 @@ void GResponse::init_members(void)
766 465 // Clear cache
767 466 m_irf_cache.clear();
768 467 m_nroi_cache.clear();
769   - m_irf_vector_cache.clear();
770 468  
771 469 // Return
772 470 return;
... ... @@ -790,7 +488,6 @@ void GResponse::copy_members(const GResponse&amp; rsp)
790 488 m_irf_diffuse_resolution = rsp.m_irf_diffuse_resolution;
791 489 m_irf_cache = rsp.m_irf_cache;
792 490 m_nroi_cache = rsp.m_nroi_cache;
793   - m_irf_vector_cache = rsp.m_irf_vector_cache;
794 491  
795 492 // Return
796 493 return;
... ... @@ -1094,6 +791,9 @@ double GResponse::irf_diffuse(const GEvent&amp; event,
1094 791 const GEnergy& srcEng = source.energy();
1095 792 const GTime& srcTime = source.time();
1096 793  
  794 + // Initialise map
  795 + //GSkyMap skymap;
  796 +
1097 797 // If model is a diffuse map model then extract sky map
1098 798 const GModelSpatialDiffuseMap* map =
1099 799 dynamic_cast<const GModelSpatialDiffuseMap*>(source.model());
... ... @@ -1257,232 +957,6 @@ double GResponse::irf_composite(const GEvent&amp; event,
1257 957  
1258 958  
1259 959 /***********************************************************************//**
1260   - * @brief Return instrument response to point source sky model
1261   - *
1262   - * @param[in] model Sky model.
1263   - * @param[in] obs Observation.
1264   - * @param[in] gradients Gradients matrix.
1265   - * @return Instrument response to point source sky model.
1266   - *
1267   - * Returns the instrument response to a point source sky model for all
1268   - * events.
1269   - ***************************************************************************/
1270   -GVector GResponse::irf_ptsrc(const GModelSky& model,
1271   - const GObservation& obs,
1272   - GMatrix* gradients) const
1273   -{
1274   - // Get number of events
1275   - int nevents = obs.events()->size();
1276   -
1277   - // Initialise result
1278   - GVector irfs(nevents);
1279   -
1280   - // Loop over events
1281   - for (int k = 0; k < nevents; ++k) {
1282   -
1283   - // Get reference to event
1284   - const GEvent& event = *((*obs.events())[k]);
1285   -
1286   - // Get source energy and time (no dispersion)
1287   - GEnergy srcEng = event.energy();
1288   - GTime srcTime = event.time();
1289   -
1290   - // Setup source
1291   - GSource source(model.name(), model.spatial(), srcEng, srcTime);
1292   -
1293   - // Get IRF value for event
1294   - irfs[k] = this->irf_ptsrc(event, source, obs);
1295   -
1296   - } // endfor: looped over events
1297   -
1298   - // Return IRF value
1299   - return irfs;
1300   -}
1301   -
1302   -
1303   -/***********************************************************************//**
1304   - * @brief Return instrument response to radial source sky model
1305   - *
1306   - * @param[in] model Sky model.
1307   - * @param[in] obs Observation.
1308   - * @param[in] gradients Gradients matrix.
1309   - * @return Instrument response to radial source sky model.
1310   - *
1311   - * Returns the instrument response to a radial source sky model for all
1312   - * events.
1313   - ***************************************************************************/
1314   -GVector GResponse::irf_radial(const GModelSky& model,
1315   - const GObservation& obs,
1316   - GMatrix* gradients) const
1317   -{
1318   - // Get number of events
1319   - int nevents = obs.events()->size();
1320   -
1321   - // Initialise result
1322   - GVector irfs(nevents);
1323   -
1324   - // Loop over events
1325   - for (int k = 0; k < nevents; ++k) {
1326   -
1327   - // Get reference to event
1328   - const GEvent& event = *((*obs.events())[k]);
1329   -
1330   - // Get source energy and time (no dispersion)
1331   - GEnergy srcEng = event.energy();
1332   - GTime srcTime = event.time();
1333   -
1334   - // Setup source
1335   - GSource source(model.name(), model.spatial(), srcEng, srcTime);
1336   -
1337   - // Get IRF value for event
1338   - irfs[k] = this->irf_radial(event, source, obs);
1339   -
1340   - } // endfor: looped over events
1341   -
1342   - // Return IRF value
1343   - return irfs;
1344   -}
1345   -
1346   -
1347   -/***********************************************************************//**
1348   - * @brief Return instrument response to ellipitical source sky model
1349   - *
1350   - * @param[in] model Sky model.
1351   - * @param[in] obs Observation.
1352   - * @param[in] gradients Gradients matrix.
1353   - * @return Instrument response to ellipitical source sky model.
1354   - *
1355   - * Returns the instrument response to a ellipitical source sky model for all
1356   - * events.
1357   - ***************************************************************************/
1358   -GVector GResponse::irf_elliptical(const GModelSky& model,
1359   - const GObservation& obs,
1360   - GMatrix* gradients) const
1361   -{
1362   - // Get number of events
1363   - int nevents = obs.events()->size();
1364   -
1365   - // Initialise result
1366   - GVector irfs(nevents);
1367   -
1368   - // Loop over events
1369   - for (int k = 0; k < nevents; ++k) {
1370   -
1371   - // Get reference to event
1372   - const GEvent& event = *((*obs.events())[k]);
1373   -
1374   - // Get source energy and time (no dispersion)
1375   - GEnergy srcEng = event.energy();
1376   - GTime srcTime = event.time();
1377   -
1378   - // Setup source
1379   - GSource source(model.name(), model.spatial(), srcEng, srcTime);
1380   -
1381   - // Get IRF value for event
1382   - irfs[k] = this->irf_elliptical(event, source, obs);
1383   -
1384   - } // endfor: looped over events
1385   -
1386   - // Return IRF value
1387   - return irfs;
1388   -}
1389   -
1390   -
1391   -/***********************************************************************//**
1392   - * @brief Return instrument response to diffuse source sky model
1393   - *
1394   - * @param[in] model Sky model.
1395   - * @param[in] obs Observation.
1396   - * @param[in] gradients Gradients matrix.
1397   - * @return Instrument response to diffuse source sky model.
1398   - *
1399   - * Returns the instrument response to a diffuse source sky model for all
1400   - * events.
1401   - ***************************************************************************/
1402   -GVector GResponse::irf_diffuse(const GModelSky& model,
1403   - const GObservation& obs,
1404   - GMatrix* gradients) const
1405   -{
1406   - // Get number of events
1407   - int nevents = obs.events()->size();
1408   -
1409   - // Initialise result
1410   - GVector irfs(nevents);
1411   -
1412   - // Loop over events
1413   - for (int k = 0; k < nevents; ++k) {
1414   -
1415   - // Get reference to event
1416   - const GEvent& event = *((*obs.events())[k]);
1417   -
1418   - // Get source energy and time (no dispersion)
1419   - GEnergy srcEng = event.energy();
1420   - GTime srcTime = event.time();
1421   -
1422   - // Setup source
1423   - GSource source(model.name(), model.spatial(), srcEng, srcTime);
1424   -
1425   - // Get IRF value for event
1426   - irfs[k] = this->irf_diffuse(event, source, obs);
1427   -
1428   - } // endfor: looped over events
1429   -
1430   - // Return IRF value
1431   - return irfs;
1432   -}
1433   -
1434   -
1435   -/***********************************************************************//**
1436   - * @brief Return instrument response to composite source sky model
1437   - *
1438   - * @param[in] model Sky model.
1439   - * @param[in] obs Observation.
1440   - * @param[in] gradients Gradients matrix.
1441   - * @return Instrument response to composite source sky model.
1442   - *
1443   - * Returns the instrument response to a composite source sky model for all
1444   - * events.
1445   - ***************************************************************************/
1446   -GVector GResponse::irf_composite(const GModelSky& model,
1447   - const GObservation& obs,
1448   - GMatrix* gradients) const
1449   -{
1450   - // Get number of events
1451   - int nevents = obs.events()->size();
1452   -
1453   - // Initialise result
1454   - GVector irfs(nevents);
1455   -
1456   - // Get pointer to composite spatial model
1457   - const GModelSpatialComposite* composite =
1458   - dynamic_cast<const GModelSpatialComposite*>(model.spatial());
1459   -
1460   - // Create copy of sky model
1461   - GModelSky sky = model;
1462   -
1463   - // Loop over model components
1464   - for (int i = 0; i < composite->components(); ++i) {
1465   -
1466   - // Set spatial component of sky model
1467   - sky.spatial(const_cast<GModelSpatial*>(composite->component(i)));
1468   -
1469   - // Compute and add IRF values
1470   - irfs += irf_spatial(sky, obs, gradients) * composite->scale(i);
1471   -
1472   - } // endfor: looped over all model components
1473   -
1474   - // Divide by number of model components
1475   - double sum = composite->sum_of_scales();
1476   - if (sum > 0.0) {
1477   - irfs /= sum;
1478   - }
1479   -
1480   - // Return IRF values
1481   - return irfs;
1482   -}
1483   -
1484   -
1485   -/***********************************************************************//**
1486 960 * @brief Convolve sky model with the instrument response
1487 961 *
1488 962 * @param[in] model Sky model.
... ... @@ -1623,346 +1097,41 @@ double GResponse::eval_prob(const GModelSky&amp; model,
1623 1097 }
1624 1098  
1625 1099  
1626   -/***********************************************************************//**
1627   - * @brief Convolve sky model with the instrument response
1628   - *
1629   - * @param[in] model Sky model.
1630   - * @param[in] obs Observation.
1631   - * @param[out] gradients Pointer to matrix of gradients.
1632   - * @return Event probabilities.
1633   - *
1634   - * Computes the event probability
1635   - *
1636   - * \f[
1637   - * P(p',E',t'|E,t) = \int S(p,E,t) \times R(p',E',t'|p,E,t) \, dp
1638   - * \f]
1639   - *
1640   - * for a all events.
1641   - ***************************************************************************/
1642   -GVector GResponse::eval_probs(const GModelSky& model,
1643   - const GObservation& obs,
1644   - GMatrixSparse* gradients) const
1645   -{
1646   - // Get number of model parameters and number of events
1647   - int npars = model.size();
1648   - int nevents = obs.events()->size();
1649   -
1650   - // Initialise gradients flag
1651   - bool grad = (gradients != NULL);
1652   -
1653   - // Initialise result
1654   - GVector probs(nevents);
1655 1100  
1656   - // Continue only if the model has a spatial component
1657   - if (model.spatial() != NULL) {
1658 1101  
1659   - // Initialise matrix to hold spatial gradients
1660   - GMatrix spat_gradients(nevents, model.spatial()->size());
1661   -
1662   - // Compute IRF value
1663   - probs = irf_spatial(model, obs, &spat_gradients);
1664   -
1665   - // Get global model scaling
1666   - double scale = (model.has_scales())
1667   - ? model.scale(obs.instrument()).value() : 1.0;
1668   -
1669   - // Initialise temporary vectors to hold gradients
1670   - GVector* tmp_gradients = NULL;
1671   - if (grad) {
1672   - tmp_gradients = new GVector[npars];
1673   - for (int i = 0; i < npars; ++i) {
1674   - tmp_gradients[i] = GVector(nevents);
1675   - }
1676   - }
1677   -
1678   - // Loop over events
1679   - for (int k = 0; k < nevents; ++k) {
1680   -
1681   - // Get probability
1682   - double spat = probs[k];
1683   -
1684   - // Continue only if spatial value is positive
1685   - if (spat > 0.0) {
1686   -
1687   - // Get reference to event
1688   - const GEvent& event = *((*obs.events())[k]);
1689   -
1690   - // Get source energy and time
1691   - GEnergy srcEng = event.energy();
1692   - GTime srcTime = event.time();
1693   -
1694   - // Evaluate spectral and temporal components
1695   - double spec = (model.spectral() != NULL)
1696   - ? model.spectral()->eval(srcEng, srcTime, grad)
1697   - : 1.0;
1698   - double temp = (model.temporal() != NULL)
1699   - ? model.temporal()->eval(srcTime, grad)
1700   - : 1.0;
1701   -
1702   - // Compute probability
1703   - probs[k] = spat * spec * temp * scale;
1704   -
1705   - // Optionally compute partial derivatives
1706   - if (grad) {
1707   -
1708   - // Get number of spatial, spectral and temporal parameters
1709   - int n_spat = model.spatial()->size();
1710   - int n_spec = model.spectral()->size();
1711   - int n_temp = model.temporal()->size();
1712   -
1713   - // Multiply factors to spatial gradients
1714   - if (model.spatial() != NULL) {
1715   - double fact = spec * temp * scale;
1716   - for (int i = 0; i < n_spat; ++i) {
1717   - tmp_gradients[i][k] = spat_gradients(k,i) * fact;
1718   - }
1719   - }
1720   -
1721   - // Multiply factors to spectral gradients
1722   - if (model.spectral() != NULL) {
1723   - int offset = n_spat;
1724   - double fact = spat * temp * scale;
1725   - for (int i = 0; i < n_spec; ++i) {
1726   - tmp_gradients[offset+i][k] =
1727   - (*(model.spectral()))[i].factor_gradient() * fact;
1728   - }
1729   - }
1730   -
1731   - // Multiply factors to temporal gradients
1732   - if (model.temporal() != NULL) {
1733   - int offset = n_spat + n_spec;
1734   - double fact = spat * spec * scale;
1735   - for (int i = 0; i < n_temp; ++i) {
1736   - tmp_gradients[offset+i][k] =
1737   - (*(model.temporal()))[i].factor_gradient() * fact;
1738   - }
1739   - }
1740   -
1741   - // Optionally set scale gradient for instrument
1742   - if (model.has_scales()) {
1743   - int offset = n_spat + n_spec + n_temp;
1744   - double fact = spat * spec * temp;
1745   - for (int i = 0; i < model.scales(); ++i) {
1746   - const GModelPar& par = model.scale(i);
1747   - double g_scale = 0.0;
1748   - if (par.name() == obs.instrument()) {
1749   - if (par.is_free()) {
1750   - tmp_gradients[offset+i][k] = par.scale() * fact;
1751   - }
1752   - }
1753   - }
1754   - }
1755   -
1756   - } // endif: computed partial derivatives
1757   -
1758   - } // endif: IRF spatial value was positive
1759   -
1760   - // Compile option: Check for NaN/Inf
1761   - #if defined(G_NAN_CHECK)
1762   - if (gammalib::is_notanumber(probs[k]) || gammalib::is_infinite(probs[k])) {
1763   - std::cout << "*** ERROR: GResponse::eval_prob:";
1764   - std::cout << " NaN/Inf encountered";
1765   - std::cout << " (probs[" << k << "]=" << probs[k];
1766   - std::cout << ", spat=" << spat;
1767   - std::cout << ")" << std::endl;
1768   - }
1769   - #endif
1770   -
1771   - } // endfor: looped over events
1772   -
1773   - // Post-process gradients
1774   - if (grad) {
1775   -
1776   - // Fill gradients into matrix
1777   - for (int i = 0; i < npars; ++i) {
1778   - gradients->column(i, tmp_gradients[i]);
1779   - }
1780   -
1781   - // Delete temporal gradients
1782   - delete [] tmp_gradients;
1783   -
1784   - } // endif: gradient post processing
1785   -
1786   - } // endif: Gamma-ray source model had a spatial component
1787   -
1788   - // Return event probabilities
1789   - return probs;
1790   -}
1791   -
1792   -
1793   -/***********************************************************************//**
1794   - * @brief Return size of vector for energy dispersion computation
1795   - *
1796   - * @param[in] model Sky model.
1797   - * @param[in] obs Observation.
1798   - * @param[out] grad Signals whether gradients should be computed.
1799   - * @return Size of vector for energy dispersion computation.
1800   - *
1801   - * Computes the size of the vector that will be computed for the computation
1802   - * of the energy dispersion.
1803   - ***************************************************************************/
1804   -int GResponse::size_edisp_vector(const GModelSky& model,
1805   - const GObservation& obs,
1806   - const bool& grad) const
1807   -{
1808   - // Initialise vector size
1809   - int size = 1;
1810   -
1811   - // Determine size of gradient part
1812   - if (grad) {
1813   - if (model.spectral() != NULL) {
1814   - for (int i = 0; i < model.spectral()->size(); ++i) {
1815   - GModelPar& par = (*(model.spectral()))[i];
1816   - if (par.is_free() && par.has_grad()) {
1817   - size++;
1818   - }
1819   - }
1820   - }
1821   - if (model.temporal() != NULL) {
1822   - for (int i = 0; i < model.temporal()->size(); ++i) {
1823   - GModelPar& par = (*(model.temporal()))[i];
1824   - if (par.is_free() && par.has_grad()) {
1825   - size++;
1826   - }
1827   - }
1828   - }
1829   - if (model.has_scales()) {
1830   - for (int i = 0; i < model.scales(); ++i) {
1831   - GModelPar& par = const_cast<GModelPar&>(model.scale(i));
1832   - if (par.name() == obs.instrument()) {
1833   - if (par.is_free() && par.has_grad()) {
1834   - size++;
1835   - }
1836   - }
1837   - }
1838   - }
1839   - }
1840   -
1841   - // Return size
1842   - return size;
1843   -}
1844 1102  
1845 1103  
1846 1104 /***********************************************************************//**
1847   - * @brief Constructor for energy dispersion integration kernels class
1848   - *
1849   - * @param[in] parent Pointer to response.
1850   - * @param[in] obs Pointer to observation.
1851   - * @param[in] model Pointer to sky model.
1852   - * @param[in] event Pointer to event.
1853   - * @param[in] srcTime True time.
1854   - * @param[in] grad Evaluate gradients?
1855   - *
1856   - * This method constructs the integration kernel needed for the energy
1857   - * dispersion computation.
1858   - ***************************************************************************/
1859   -GResponse::edisp_kerns::edisp_kerns(const GResponse* parent,
1860   - const GObservation* obs,
1861   - const GModelSky* model,
1862   - const GEvent* event,
1863   - const GTime& srcTime,
1864   - const bool& grad)
1865   -{
1866   - // Set members
1867   - m_parent = parent;
1868   - m_obs = obs;
1869   - m_model = model;
1870   - m_event = event;
1871   - m_size = 1;
1872   - m_srcTime = srcTime;
1873   - m_grad = grad;
1874   -
1875   - // If gradients are requested then put pointers to all relevant parameter
1876   - // in the parameter pointer vector
1877   - m_pars.clear();
1878   - if (m_grad) {
1879   -
1880   - // Add free spectral parameters with analytical gradients
1881   - if (m_model->spectral() != NULL) {
1882   - for (int i = 0; i < m_model->spectral()->size(); ++i) {
1883   - GModelPar& par = (*(m_model->spectral()))[i];
1884   - if (par.is_free() && par.has_grad()) {
1885   - m_pars.push_back(&par);
1886   - }
1887   - }
1888   - }
1889   -
1890   - // Add free temporal parameters with analytical gradients
1891   - if (m_model->temporal() != NULL) {
1892   - for (int i = 0; i < m_model->temporal()->size(); ++i) {
1893   - GModelPar& par = (*(m_model->temporal()))[i];
1894   - if (par.is_free() && par.has_grad()) {
1895   - m_pars.push_back(&par);
1896   - }
1897   - }
1898   - }
1899   -
1900   - // Add free scaling factors
1901   - if (m_model->has_scales()) {
1902   - for (int i = 0; i < m_model->scales(); ++i) {
1903   - GModelPar& par = const_cast<GModelPar&>(m_model->scale(i));
1904   - if (par.name() == m_obs->instrument()) {
1905   - if (par.is_free() && par.has_grad()) {
1906   - m_pars.push_back(&par);
1907   - }
1908   - }
1909   - }
1910   - }
1911   -
1912   - // Update size
1913   - m_size += m_pars.size();
1914   -
1915   - } // endif: gradients were requested
1916   -
1917   - // Return
1918   - return;
1919   -}
1920   -
1921   -
1922   -/***********************************************************************//**
1923   - * @brief Evaluate energy dispersion integration kernel
  1105 + * @brief Integration kernel for GResponse::edisp_kern() class
1924 1106 *
1925 1107 * @param[in] etrue True photon energy in MeV.
1926 1108 *
1927 1109 * This method implements the integration kernel needed for the
1928 1110 * GResponse::edisp_kern() class.
1929 1111 ***************************************************************************/
1930   -GVector GResponse::edisp_kerns::eval(const double& etrue)
  1112 +double GResponse::edisp_kern::eval(const double& etrue)
1931 1113 {
1932   - // Initialise result
1933   - GVector kernels(m_size);
1934   -
1935   - // Initialise Ndarray array index
1936   - int index = 0;
1937   -
1938 1114 // Set true energy
1939 1115 GEnergy srcEng;
1940 1116 srcEng.MeV(etrue);
1941 1117  
1942 1118 // Get function value
1943   - kernels[index++] = m_parent->eval_prob(*m_model, *m_event,
1944   - srcEng, m_srcTime, *m_obs,
1945   - m_grad);
1946   -
1947   - // If gradients are requested then extract them and put them into the
1948   - // array
1949   - if (m_grad) {
1950   - for (int i = 0; i < m_pars.size(); ++i, ++index) {
1951   - kernels[index] = m_pars[i]->factor_gradient();
1952   - }
1953   - }
  1119 + double value = m_parent->eval_prob(*m_model, *m_event,
  1120 + srcEng, m_srcTime, *m_obs, m_grad);
1954 1121  
1955 1122 // Compile option: Check for NaN
1956 1123 #if defined(G_NAN_CHECK)
1957   - if (gammalib::is_notanumber(kernels[0]) || gammalib::is_infinite(kernels[0])) {
  1124 + if (gammalib::is_notanumber(value) || gammalib::is_infinite(value)) {
1958 1125 std::cout << "*** ERROR: GResponse::edisp_kern::eval";
1959   - std::cout << "(etrue=" << etrue << "): NaN/Inf encountered";
1960   - std::cout << " (value=" << kernels[0] << ")" << std::endl;
  1126 + std::cout << "(etrue=" << etrue << "): ";
  1127 + std::cout << " NaN/Inf encountered";
  1128 + std::cout << " (value=" << value;
  1129 + std::cout << ")" << std::endl;
1961 1130 }
1962 1131 #endif
1963 1132  
1964   - // Return kernels
1965   - return kernels;
  1133 + // Return value
  1134 + return value;
1966 1135 }
1967 1136  
1968 1137  
... ... @@ -2141,3 +1310,6 @@ double GResponse::irf_elliptical_kern_phi::eval(const double&amp; phi)
2141 1310 // Return
2142 1311 return irf;
2143 1312 }
  1313 +
  1314 +
  1315 +
... ...