Pierrick Martin / gammalib

Commit d511521a2b0ab0d45af1341da75d042dee83bd68

Authored by Jürgen Knödlseder
1 parent 52ba4ad2

Implement GResponse::irf(const GEvent& event, const GSource& source, const GObservation& obs) method 

The method is now implemented at the base class level and branches to the model specific IRF computation methods that are also implemented at the GResponse level. This will allow implementing generic model dependent IRF methods at the base class level, making all spatial model types automatically accessible to all instrument modules. These methods can always be overloaded later in case that more efficient code needs to be generated. But at least, all spatial models work a priori out of the box.
Inline Side-by-side
Showing 14 changed files with 347 additions and 296 deletions
include/GResponse.hpp
  View file @d511521
1 1 /***************************************************************************
2 2 * GResponse.hpp - Abstract response base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2008-2018 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 *
... ... @@ -32,6 +32,7 @@
32 32 #include "GBase.hpp"
33 33 #include "GFunction.hpp"
34 34 #include "GTime.hpp"
  35 +#include "GResponseCache.hpp"
35 36  
36 37 /* __ Forward declarations _______________________________________________ */
37 38 class GEvent;
... ... @@ -84,9 +85,6 @@ public:
84 85 virtual double irf(const GEvent& event,
85 86 const GPhoton& photon,
86 87 const GObservation& obs) const = 0;
87   - virtual double irf(const GEvent& event,
88   - const GSource& source,
89   - const GObservation& obs) const = 0;
90 88 virtual double nroi(const GModelSky& model,
91 89 const GEnergy& obsEng,
92 90 const GTime& obsTime,
... ... @@ -99,6 +97,9 @@ public:
99 97 const GEvent& event,
100 98 const GObservation& obs,
101 99 const bool& grad = true) const;
  100 + virtual double irf(const GEvent& event,
  101 + const GSource& source,
  102 + const GObservation& obs) const;
102 103  
103 104 protected:
104 105 // Protected methods
... ... @@ -112,6 +113,23 @@ protected:
112 113 const GObservation& obs,
113 114 const bool& grad) const;
114 115  
  116 + // Virtual protected methods
  117 + virtual double irf_ptsrc(const GEvent& event,
  118 + const GSource& source,
  119 + const GObservation& obs) const;
  120 + virtual double irf_radial(const GEvent& event,
  121 + const GSource& source,
  122 + const GObservation& obs) const;
  123 + virtual double irf_elliptical(const GEvent& event,
  124 + const GSource& source,
  125 + const GObservation& obs) const;
  126 + virtual double irf_diffuse(const GEvent& event,
  127 + const GSource& source,
  128 + const GObservation& obs) const;
  129 + virtual double irf_composite(const GEvent& event,
  130 + const GSource& source,
  131 + const GObservation& obs) const;
  132 +
115 133 // Protected classes
116 134 class edisp_kern : public GFunction {
117 135 public:
... ... @@ -136,6 +154,16 @@ protected:
136 154 GTime m_srcTime; //!< True arrival time
137 155 bool m_grad; //!< Gradient flag
138 156 };
  157 +
  158 + // Protected members
  159 + bool m_use_irf_cache; //!< Control usage of irf cache
  160 + bool m_use_nroi_cache; //!< Control usage of nroi cache
  161 +
  162 + // Cache for irf(GEvent&, GSource&, GObservation&) and
  163 + // nroi(GModelSky&, GEnergy&, GTime&, GEnergy&, GTime&, GObservation&)
  164 + // computations
  165 + mutable GResponseCache m_irf_cache;
  166 + mutable GResponseCache m_nroi_cache;
139 167 };
140 168  
141 169 #endif /* GRESPONSE_HPP */
... ...
include/GResponseCache.hpp
  View file @d511521
... ... @@ -69,6 +69,7 @@ public:
69 69 std::string classname(void) const;
70 70 int size(void) const;
71 71 bool is_empty(void) const;
  72 + void energy_scale(const GEnergy& energy);
72 73 void set(const std::string& name,
73 74 const GEnergy& ereco,
74 75 const GEnergy& etrue,
... ... @@ -102,7 +103,8 @@ protected:
102 103 const GEnergy& etrue) const;
103 104  
104 105 // Protected members
105   - GResponseCacheName m_cache;
  106 + GResponseCacheName m_cache; //<! Response cache
  107 + double m_energy_scale; //<! Energy scale 1/MeV
106 108 };
107 109  
108 110  
... ... @@ -133,6 +135,21 @@ bool GResponseCache::is_empty(void) const
133 135  
134 136  
135 137 /***********************************************************************//**
  138 + * @brief Set minimum energy of instrument
  139 + *
  140 + * @param[in] energy Energy.
  141 + *
  142 + * Sets minimum energy of instrument for hash scaling.
  143 + ***************************************************************************/
  144 +inline
  145 +void GResponseCache::energy_scale(const GEnergy& energy)
  146 +{
  147 + m_energy_scale = 1.0 / energy.MeV();
  148 + return;
  149 +}
  150 +
  151 +
  152 +/***********************************************************************//**
136 153 * @brief Remove cache for source
137 154 *
138 155 * @param[in] name Source name.
... ...
inst/cta/include/GCTAInstDir.hpp
  View file @d511521
... ... @@ -128,7 +128,7 @@ std::string GCTAInstDir::classname(void) const
128 128 inline
129 129 double GCTAInstDir::hash(void) const
130 130 {
131   - double hash = m_dir.ra_deg() + m_dir.dec_deg() * 1.0e2;
  131 + double hash = m_dir.ra_deg() * 1.0e3 + m_dir.dec_deg();
132 132 return hash;
133 133 }
134 134  
... ...
inst/cta/include/GCTAResponse.hpp
  View file @d511521
1 1 /***************************************************************************
2 2 * GCTAResponse.hpp - CTA response abstract base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2019 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -86,45 +86,24 @@ public:
86 86 virtual std::string print(const GChatter& chatter = NORMAL) const = 0;
87 87  
88 88 // Implemeted methods
89   - virtual double irf(const GEvent& event,
90   - const GSource& source,
91   - const GObservation& obs) const;
92 89 virtual void remove_response_cache(const std::string& name);
93 90  
  91 + // Make sure the base class methods are visible. This is required here
  92 + // since there are different overloaded irf() methods, and not making
  93 + // the base class methods available would hide the
  94 + //
  95 + // double GResponse::irf(const GEvent& event,
  96 + // const GSource& source,
  97 + // const GObservation& obs) const
  98 + //
  99 + // method.
  100 + using GResponse::irf;
  101 +
94 102 protected:
95 103 // Protected methods
96 104 void init_members(void);
97 105 void copy_members(const GCTAResponse& rsp);
98 106 void free_members(void);
99   -
100   - // Pure virtual protected methods
101   - virtual double irf_ptsrc(const GEvent& event,
102   - const GSource& source,
103   - const GObservation& obs) const = 0;
104   - virtual double irf_radial(const GEvent& event,
105   - const GSource& source,
106   - const GObservation& obs) const = 0;
107   - virtual double irf_elliptical(const GEvent& event,
108   - const GSource& source,
109   - const GObservation& obs) const = 0;
110   - virtual double irf_diffuse(const GEvent& event,
111   - const GSource& source,
112   - const GObservation& obs) const = 0;
113   -
114   - // Implemented protected methods
115   - virtual double irf_composite(const GEvent& event,
116   - const GSource& source,
117   - const GObservation& obs) const;
118   -
119   - // Protected members
120   - bool m_use_irf_cache; //!< Control usage of irf cache
121   - bool m_use_nroi_cache; //!< Control usage of nroi cache
122   -
123   - // Cache for irf(GEvent&, GSource&, GObservation&) and
124   - // nroi(GModelSky&, GEnergy&, GTime&, GEnergy&, GTime&, GObservation&)
125   - // computations
126   - mutable GCTAResponseCache m_irf_cache;
127   - mutable GCTAResponseCache m_nroi_cache;
128 107 };
129 108  
130 109 #endif /* GCTARESPONSE_HPP */
... ...
inst/cta/src/GCTAResponse.cpp
  View file @d511521
1 1 /***************************************************************************
2 2 * GCTAResponse.cpp - CTA Response class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2019 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2020 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -152,73 +152,6 @@ GCTAResponse&amp; GCTAResponse::operator=(const GCTAResponse&amp; rsp)
152 152 ==========================================================================*/
153 153  
154 154 /***********************************************************************//**
155   - * @brief Return instrument response
156   - *
157   - * @param[in] event Event.
158   - * @param[in] source Source.
159   - * @param[in] obs Observation.
160   - * @return Instrument response.
161   - *
162   - * Returns the instrument response for a given event, source and observation.
163   - ***************************************************************************/
164   -double GCTAResponse::irf(const GEvent& event,
165   - const GSource& source,
166   - const GObservation& obs) const
167   -{
168   - // Initialise IRF value
169   - double irf = 0.0;
170   -
171   - // Set IRF value attributes
172   - std::string name = obs.id() + "::" + source.name();
173   - const GCTAInstDir& dir = gammalib::cta_dir(G_IRF, event);
174   - const GEnergy& ereco = event.energy();
175   - const GEnergy& etrue = source.energy();
176   -
177   - // Signal if spatial model has free parameters
178   - bool has_free_pars = source.model()->has_free_pars();
179   -
180   - // If the spatial model component has free parameters, or the response
181   - // cache should not be used, or the cache does not contain the requested
182   - // IRF value then compute the IRF value for the spatial model.
183   - if (has_free_pars ||
184   - !m_use_irf_cache ||
185   - !m_irf_cache.contains(name, dir, ereco, etrue, &irf)) {
186   -
187   - // Compute IRF for spatial model
188   - switch (source.model()->code()) {
189   - case GMODEL_SPATIAL_POINT_SOURCE:
190   - irf = irf_ptsrc(event, source, obs);
191   - break;
192   - case GMODEL_SPATIAL_RADIAL:
193   - irf = irf_radial(event, source, obs);
194   - break;
195   - case GMODEL_SPATIAL_ELLIPTICAL:
196   - irf = irf_elliptical(event, source, obs);
197   - break;
198   - case GMODEL_SPATIAL_DIFFUSE:
199   - irf = irf_diffuse(event, source, obs);
200   - break;
201   - case GMODEL_SPATIAL_COMPOSITE:
202   - irf = irf_composite(event, source, obs);
203   - break;
204   - default:
205   - break;
206   - }
207   -
208   - } // endif: computed spatial model
209   -
210   - // If the spatial model has no free parameters and the response cache
211   - // should be used then put the IRF value in the response cache.
212   - if (!has_free_pars && m_use_irf_cache) {
213   - m_irf_cache.set(name, dir, ereco, etrue, irf);
214   - }
215   -
216   - // Return IRF value
217   - return irf;
218   -}
219   -
220   -
221   -/***********************************************************************//**
222 155 * @brief Remove response cache for model
223 156 *
224 157 * @param[in] name Model name.
... ... @@ -250,8 +183,10 @@ void GCTAResponse::init_members(void)
250 183 // Initialize members
251 184 m_use_irf_cache = true; //!< Switched on by default
252 185 m_use_nroi_cache = true; //!< Switched on by default
253   - m_irf_cache.clear();
254   - m_nroi_cache.clear();
  186 +
  187 + // Set energy scale for response cache to TeV
  188 + m_irf_cache.energy_scale(GEnergy(1.0, "TeV"));
  189 + m_nroi_cache.energy_scale(GEnergy(1.0, "TeV"));
255 190  
256 191 // Return
257 192 return;
... ... @@ -265,12 +200,6 @@ void GCTAResponse::init_members(void)
265 200 ***************************************************************************/
266 201 void GCTAResponse::copy_members(const GCTAResponse& rsp)
267 202 {
268   - // Copy members
269   - m_use_irf_cache = rsp.m_use_irf_cache;
270   - m_use_nroi_cache = rsp.m_use_nroi_cache;
271   - m_irf_cache = rsp.m_irf_cache;
272   - m_nroi_cache = rsp.m_nroi_cache;
273   -
274 203 // Return
275 204 return;
276 205 }
... ... @@ -284,49 +213,3 @@ void GCTAResponse::free_members(void)
284 213 // Return
285 214 return;
286 215 }
287   -
288   -
289   -/***********************************************************************//**
290   - * @brief Return instrument response to composite source
291   - *
292   - * @param[in] event Observed event.
293   - * @param[in] source Source.
294   - * @param[in] obs Observation.
295   - * @return Instrument response to composite source.
296   - *
297   - * Returns the instrument response to a specified composite source.
298   - ***************************************************************************/
299   -double GCTAResponse::irf_composite(const GEvent& event,
300   - const GSource& source,
301   - const GObservation& obs) const
302   -{
303   - // Initialise IRF
304   - double irf = 0.0;
305   -
306   - // Get pointer to composite model
307   - const GModelSpatialComposite* model =
308   - dynamic_cast<const GModelSpatialComposite*>(source.model());
309   -
310   - // Loop over model components
311   - for (int i = 0; i < model->components(); ++i) {
312   -
313   - // Get pointer to spatial component
314   - GModelSpatial* spat = const_cast<GModelSpatial*>(model->component(i));
315   -
316   - // Create new GSource object
317   - GSource src(source.name(), spat, source.energy(), source.time());
318   -
319   - // Compute irf value
320   - irf += this->irf(event, src, obs) * model->scale(i);
321   -
322   - }
323   -
324   - // Divide by number of model components
325   - double sum = model->sum_of_scales();
326   - if (sum > 0.0) {
327   - irf /= sum;
328   - }
329   -
330   - // Return IRF value
331   - return irf;
332   -}
... ...
inst/spi/include/GSPIResponse.hpp
  View file @d511521
... ... @@ -81,9 +81,6 @@ public:
81 81 virtual double irf(const GEvent& event,
82 82 const GPhoton& photon,
83 83 const GObservation& obs) const;
84   - virtual double irf(const GEvent& event,
85   - const GSource& source,
86   - const GObservation& obs) const;
87 84 virtual double nroi(const GModelSky& model,
88 85 const GEnergy& obsEng,
89 86 const GTime& obsTime,
... ... @@ -92,9 +89,9 @@ public:
92 89 virtual std::string print(const GChatter& chatter = NORMAL) const;
93 90  
94 91 // Other Methods
95   - double irf(const GSkyDir& srcDir,
96   - const GSPIEventBin& bin,
97   - const int& ireg) const;
  92 + double irf_value(const GSkyDir& srcDir,
  93 + const GSPIEventBin& bin,
  94 + const int& ireg) const;
98 95 void set(const GSPIObservation& obs, const GEnergy& energy = GEnergy());
99 96 void rspname(const GFilename& rspname);
100 97 double zenith(const int& ipt, const GSkyDir& dir) const;
... ...
inst/spi/pyext/GSPIResponse.i
  View file @d511521
... ... @@ -52,9 +52,6 @@ public:
52 52 virtual double irf(const GEvent& event,
53 53 const GPhoton& photon,
54 54 const GObservation& obs) const;
55   - virtual double irf(const GEvent& event,
56   - const GSource& source,
57   - const GObservation& obs) const;
58 55 virtual double nroi(const GModelSky& model,
59 56 const GEnergy& obsEng,
60 57 const GTime& obsTime,
... ... @@ -62,9 +59,9 @@ public:
62 59 virtual GEbounds ebounds(const GEnergy& obsEnergy) const;
63 60  
64 61 // Other Methods
65   - double irf(const GSkyDir& srcDir,
66   - const GSPIEventBin& bin,
67   - const int& ireg) const;
  62 + double irf_value(const GSkyDir& srcDir,
  63 + const GSPIEventBin& bin,
  64 + const int& ireg) const;
68 65 void set(const GSPIObservation& obs, const GEnergy& energy = GEnergy());
69 66 void rspname(const GFilename& rspname);
70 67 double zenith(const int& ipt, const GSkyDir& dir) const;
... ...
inst/spi/src/GSPIResponse.cpp
  View file @d511521
... ... @@ -288,7 +288,7 @@ double GSPIResponse::irf(const GEvent&amp; event,
288 288 if (cube->ebounds().index(photon.energy()) == bin->iebin()) {
289 289  
290 290 // Get IRF value for photo peak
291   - irf = this->irf(photon.dir(), *bin, 0);
  291 + irf = irf_value(photon.dir(), *bin, 0);
292 292  
293 293 // Compile option: Check for NaN/Inf
294 294 #if defined(G_NAN_CHECK)
... ... @@ -323,9 +323,9 @@ double GSPIResponse::irf(const GEvent&amp; event,
323 323 * event bin. The value of the IRF is bilinearly interpolated from the
324 324 * pre-computed IRFs cube that is stored in the class.
325 325 ***************************************************************************/
326   -double GSPIResponse::irf(const GSkyDir& srcDir,
327   - const GSPIEventBin& bin,
328   - const int& ireg) const
  326 +double GSPIResponse::irf_value(const GSkyDir& srcDir,
  327 + const GSPIEventBin& bin,
  328 + const int& ireg) const
329 329 {
330 330 // Initialise IRF value
331 331 double irf = 0.0;
... ... @@ -403,56 +403,6 @@ double GSPIResponse::irf(const GSkyDir&amp; srcDir,
403 403  
404 404  
405 405 /***********************************************************************//**
406   - * @brief Return value of INTEGRAL/SPI instrument response for a source
407   - *
408   - * @param[in] event Event.
409   - * @param[in] source Source.
410   - * @param[in] obs Observation.
411   - * @return Instrument response $\f(cm^2 sr^{-1})$\f
412   - *
413   - * @exception GException::feature_not_implemented
414   - * Source model is not supported.
415   - *
416   - * Returns the instrument response for a given event, source and observation.
417   - *
418   - * @todo Implement whatever is needed here.
419   - ***************************************************************************/
420   -double GSPIResponse::irf(const GEvent& event,
421   - const GSource& source,
422   - const GObservation& obs) const
423   -{
424   - // Initialise IRF value
425   - double irf = 0.0;
426   -
427   - // Select IRF depending on the spatial model type
428   - switch (source.model()->code()) {
429   - case GMODEL_SPATIAL_POINT_SOURCE:
430   - {
431   - const GModelSpatialPointSource* src =
432   - static_cast<const GModelSpatialPointSource*>(source.model());
433   - GPhoton photon(src->dir(), source.energy(), source.time());
434   - irf = this->irf(event, photon, obs);
435   - }
436   - break;
437   - case GMODEL_SPATIAL_RADIAL:
438   - case GMODEL_SPATIAL_ELLIPTICAL:
439   - case GMODEL_SPATIAL_DIFFUSE:
440   - {
441   - std::string msg = "Response computation not yet implemented for "
442   - "spatial model type \""+source.model()->type()+"\".";
443   - throw GException::feature_not_implemented(G_IRF, msg);
444   - }
445   - break;
446   - default:
447   - break;
448   - }
449   -
450   - // Return IRF value
451   - return irf;
452   -}
453   -
454   -
455   -/***********************************************************************//**
456 406 * @brief Return integral of event probability for a given sky model over ROI
457 407 *
458 408 * @param[in] model Sky model.
... ... @@ -777,6 +727,10 @@ void GSPIResponse::init_members(void)
777 727 m_wcs_ypix_max = 0.0;
778 728 m_max_zenith = 180.0 * gammalib::deg2rad;
779 729  
  730 + // Set energy scale for response cache to MeV
  731 + //m_irf_cache.energy_scale(GEnergy(1.0, "MeV"));
  732 + //m_nroi_cache.energy_scale(GEnergy(1.0, "MeV"));
  733 +
780 734 // Return
781 735 return;
782 736 }
... ...
pyext/GResponse.i
  View file @d511521
1 1 /***************************************************************************
2 2 * GResponse.i - Abstract response base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2008-2015 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 *
... ... @@ -51,9 +51,6 @@ public:
51 51 virtual double irf(const GEvent& event,
52 52 const GPhoton& photon,
53 53 const GObservation& obs) const = 0;
54   - virtual double irf(const GEvent& event,
55   - const GSource& source,
56   - const GObservation& obs) const = 0;
57 54 virtual double nroi(const GModelSky& model,
58 55 const GEnergy& obsEng,
59 56 const GTime& obsTime,
... ... @@ -65,6 +62,9 @@ public:
65 62 const GEvent& event,
66 63 const GObservation& obs,
67 64 const bool& grad = true) const;
  65 + virtual double irf(const GEvent& event,
  66 + const GSource& source,
  67 + const GObservation& obs) const;
68 68 };
69 69  
70 70  
... ...
src/obs/GResponse.cpp
  View file @d511521
1 1 /***************************************************************************
2 2 * GResponse.cpp - Abstract response base class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2008-2019 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 *
... ... @@ -39,10 +39,12 @@
39 39 #include "GPhoton.hpp"
40 40 #include "GEnergy.hpp"
41 41 #include "GTime.hpp"
42   -#include "GSource.hpp" // will become obsolete
43   -#include "GEbounds.hpp" // will become obsolete
  42 +#include "GSource.hpp"
  43 +#include "GEbounds.hpp"
44 44 #include "GObservation.hpp"
45 45 #include "GModelSky.hpp"
  46 +#include "GModelSpatialPointSource.hpp"
  47 +#include "GModelSpatialComposite.hpp"
46 48  
47 49 /* __ Method name definitions ____________________________________________ */
48 50 #define G_IRF_RADIAL "GResponse::irf_radial(GEvent&, GSource&,"\
... ... @@ -248,6 +250,73 @@ double GResponse::convolve(const GModelSky&amp; model,
248 250 }
249 251  
250 252  
  253 +/***********************************************************************//**
  254 + * @brief Return instrument response
  255 + *
  256 + * @param[in] event Event.
  257 + * @param[in] source Source.
  258 + * @param[in] obs Observation.
  259 + * @return Instrument response.
  260 + *
  261 + * Returns the instrument response for a given event, source and observation.
  262 + ***************************************************************************/
  263 +double GResponse::irf(const GEvent& event,
  264 + const GSource& source,
  265 + const GObservation& obs) const
  266 +{
  267 + // Initialise IRF value
  268 + double irf = 0.0;
  269 +
  270 + // Set IRF value attributes
  271 + std::string name = obs.id() + "::" + source.name();
  272 + const GInstDir& dir = event.dir();
  273 + const GEnergy& ereco = event.energy();
  274 + const GEnergy& etrue = source.energy();
  275 +
  276 + // Signal if spatial model has free parameters
  277 + bool has_free_pars = source.model()->has_free_pars();
  278 +
  279 + // If the spatial model component has free parameters, or the response
  280 + // cache should not be used, or the cache does not contain the requested
  281 + // IRF value then compute the IRF value for the spatial model.
  282 + if (has_free_pars ||
  283 + !m_use_irf_cache ||
  284 + !m_irf_cache.contains(name, dir, ereco, etrue, &irf)) {
  285 +
  286 + // Compute IRF for spatial model
  287 + switch (source.model()->code()) {
  288 + case GMODEL_SPATIAL_POINT_SOURCE:
  289 + irf = irf_ptsrc(event, source, obs);
  290 + break;
  291 + case GMODEL_SPATIAL_RADIAL:
  292 + irf = irf_radial(event, source, obs);
  293 + break;
  294 + case GMODEL_SPATIAL_ELLIPTICAL:
  295 + irf = irf_elliptical(event, source, obs);
  296 + break;
  297 + case GMODEL_SPATIAL_DIFFUSE:
  298 + irf = irf_diffuse(event, source, obs);
  299 + break;
  300 + case GMODEL_SPATIAL_COMPOSITE:
  301 + irf = irf_composite(event, source, obs);
  302 + break;
  303 + default:
  304 + break;
  305 + }
  306 +
  307 + } // endif: computed spatial model
  308 +
  309 + // If the spatial model has no free parameters and the response cache
  310 + // should be used then put the IRF value in the response cache.
  311 + if (!has_free_pars && m_use_irf_cache) {
  312 + m_irf_cache.set(name, dir, ereco, etrue, irf);
  313 + }
  314 +
  315 + // Return IRF value
  316 + return irf;
  317 +}
  318 +
  319 +
251 320 /*==========================================================================
252 321 = =
253 322 = Protected methods =
... ... @@ -259,6 +328,12 @@ double GResponse::convolve(const GModelSky&amp; model,
259 328 ***************************************************************************/
260 329 void GResponse::init_members(void)
261 330 {
  331 + // Initialize members
  332 + m_use_irf_cache = false; //!< Switched off by default
  333 + m_use_nroi_cache = false; //!< Switched off by default
  334 + m_irf_cache.clear();
  335 + m_nroi_cache.clear();
  336 +
262 337 // Return
263 338 return;
264 339 }
... ... @@ -271,6 +346,12 @@ void GResponse::init_members(void)
271 346 ***************************************************************************/
272 347 void GResponse::copy_members(const GResponse& rsp)
273 348 {
  349 + // Copy members
  350 + m_use_irf_cache = rsp.m_use_irf_cache;
  351 + m_use_nroi_cache = rsp.m_use_nroi_cache;
  352 + m_irf_cache = rsp.m_irf_cache;
  353 + m_nroi_cache = rsp.m_nroi_cache;
  354 +
274 355 // Return
275 356 return;
276 357 }
... ... @@ -287,6 +368,168 @@ void GResponse::free_members(void)
287 368  
288 369  
289 370 /***********************************************************************//**
  371 + * @brief Return instrument response to point source
  372 + *
  373 + * @param[in] event Observed event.
  374 + * @param[in] source Source.
  375 + * @param[in] obs Observation.
  376 + * @return Instrument response to point source.
  377 + *
  378 + * Returns the instrument response to a point source.
  379 + ***************************************************************************/
  380 +double GResponse::irf_ptsrc(const GEvent& event,
  381 + const GSource& source,
  382 + const GObservation& obs) const
  383 +{
  384 + // Get pointer to point source model
  385 + const GModelSpatialPointSource* src =
  386 + static_cast<const GModelSpatialPointSource*>(source.model());
  387 +
  388 + // Setup photon
  389 + GPhoton photon(src->dir(), source.energy(), source.time());
  390 +
  391 + // Get IRF
  392 + double irf = this->irf(event, photon, obs);
  393 +
  394 + // Return IRF value
  395 + return irf;
  396 +}
  397 +
  398 +
  399 +/***********************************************************************//**
  400 + * @brief Return instrument response to radial source
  401 + *
  402 + * @param[in] event Observed event.
  403 + * @param[in] source Source.
  404 + * @param[in] obs Observation.
  405 + * @return Instrument response to radial source.
  406 + *
  407 + * Returns the instrument response to a radial source.
  408 + *
  409 + * @todo Implement method
  410 + ***************************************************************************/
  411 +double GResponse::irf_radial(const GEvent& event,
  412 + const GSource& source,
  413 + const GObservation& obs) const
  414 +{
  415 + // Initialise IRF
  416 + double irf = 0.0;
  417 +
  418 + // Throw exception
  419 + std::string msg = "Response computation not yet implemented for spatial "
  420 + "model type \""+source.model()->type()+"\".";
  421 + throw GException::feature_not_implemented(G_IRF_RADIAL, msg);
  422 +
  423 + // Return IRF value
  424 + return irf;
  425 +}
  426 +
  427 +
  428 +/***********************************************************************//**
  429 + * @brief Return instrument response to elliptical source
  430 + *
  431 + * @param[in] event Observed event.
  432 + * @param[in] source Source.
  433 + * @param[in] obs Observation.
  434 + * @return Instrument response to elliptical source.
  435 + *
  436 + * Returns the instrument response to a elliptical source.
  437 + *
  438 + * @todo Implement method
  439 + ***************************************************************************/
  440 +double GResponse::irf_elliptical(const GEvent& event,
  441 + const GSource& source,
  442 + const GObservation& obs) const
  443 +{
  444 + // Initialise IRF
  445 + double irf = 0.0;
  446 +
  447 + // Throw exception
  448 + std::string msg = "Response computation not yet implemented for spatial "
  449 + "model type \""+source.model()->type()+"\".";
  450 + throw GException::feature_not_implemented(G_IRF_ELLIPTICAL, msg);
  451 +
  452 + // Return IRF value
  453 + return irf;
  454 +}
  455 +
  456 +
  457 +/***********************************************************************//**
  458 + * @brief Return instrument response to diffuse source
  459 + *
  460 + * @param[in] event Observed event.
  461 + * @param[in] source Source.
  462 + * @param[in] obs Observation.
  463 + * @return Instrument response to diffuse source.
  464 + *
  465 + * Returns the instrument response to a diffuse source.
  466 + *
  467 + * @todo Implement method
  468 + ***************************************************************************/
  469 +double GResponse::irf_diffuse(const GEvent& event,
  470 + const GSource& source,
  471 + const GObservation& obs) const
  472 +{
  473 + // Initialise IRF
  474 + double irf = 0.0;
  475 +
  476 + // Throw exception
  477 + std::string msg = "Response computation not yet implemented for spatial "
  478 + "model type \""+source.model()->type()+"\".";
  479 + throw GException::feature_not_implemented(G_IRF_DIFFUSE, msg);
  480 +
  481 + // Return IRF value
  482 + return irf;
  483 +}
  484 +
  485 +
  486 +/***********************************************************************//**
  487 + * @brief Return instrument response to composite source
  488 + *
  489 + * @param[in] event Observed event.
  490 + * @param[in] source Source.
  491 + * @param[in] obs Observation.
  492 + * @return Instrument response to composite source.
  493 + *
  494 + * Returns the instrument response to a specified composite source.
  495 + ***************************************************************************/
  496 +double GResponse::irf_composite(const GEvent& event,
  497 + const GSource& source,
  498 + const GObservation& obs) const
  499 +{
  500 + // Initialise IRF
  501 + double irf = 0.0;
  502 +
  503 + // Get pointer to composite model
  504 + const GModelSpatialComposite* model =
  505 + dynamic_cast<const GModelSpatialComposite*>(source.model());
  506 +
  507 + // Loop over model components
  508 + for (int i = 0; i < model->components(); ++i) {
  509 +
  510 + // Get pointer to spatial component
  511 + GModelSpatial* spat = const_cast<GModelSpatial*>(model->component(i));
  512 +
  513 + // Create new GSource object
  514 + GSource src(source.name(), spat, source.energy(), source.time());
  515 +
  516 + // Compute irf value
  517 + irf += this->irf(event, src, obs) * model->scale(i);
  518 +
  519 + }
  520 +
  521 + // Divide by number of model components
  522 + double sum = model->sum_of_scales();
  523 + if (sum > 0.0) {
  524 + irf /= sum;
  525 + }
  526 +
  527 + // Return IRF value
  528 + return irf;
  529 +}
  530 +
  531 +
  532 +/***********************************************************************//**
290 533 * @brief Convolve sky model with the instrument response
291 534 *
292 535 * @param[in] model Sky model.
... ...
src/obs/GResponseCache.cpp
  View file @d511521
... ... @@ -389,6 +389,7 @@ void GResponseCache::init_members(void)
389 389 {
390 390 // Initialise members
391 391 m_cache.clear();
  392 + m_energy_scale = 1.0;
392 393  
393 394 // Return
394 395 return;
... ... @@ -403,7 +404,8 @@ void GResponseCache::init_members(void)
403 404 void GResponseCache::copy_members(const GResponseCache& cache)
404 405 {
405 406 // Copy members
406   - m_cache = cache.m_cache;
  407 + m_cache = cache.m_cache;
  408 + m_energy_scale = cache.m_energy_scale;
407 409  
408 410 // Return
409 411 return;
... ... @@ -437,12 +439,15 @@ void GResponseCache::free_members(void)
437 439 * where
438 440 * - \f$E_{\rm reco}\f$ is the reconstructued energy in TeV, and
439 441 * - \f$E_{\rm true}\f$ is the true energy in TeV.
  442 + *
  443 + * @todo Verify unique encoding for all instruments!!!
440 444 ***************************************************************************/
441 445 double GResponseCache::encode(const GEnergy& ereco,
442 446 const GEnergy& etrue) const
443 447 {
444 448 // Construct encoded instrument direction and reconstructued energy
445   - double encoded = ereco.MeV() + etrue.MeV() * 1.0e2;
  449 + double encoded = (ereco.MeV() * 1.0e2 + etrue.MeV()) * m_energy_scale;
  450 + //double encoded = ereco.TeV() * 1.0e2 + etrue.TeV();
446 451  
447 452 // Return encoded value
448 453 return encoded;
... ... @@ -463,7 +468,7 @@ double GResponseCache::encode(const GEnergy&amp; ereco,
463 468 *
464 469 * \f[
465 470 * {\tt encode} = {\tt dir.hash()} + E_{\rm reco} \times 10^6 +
466   - E_{\rm true} \times 10^8
  471 + * E_{\rm true} \times 10^8
467 472 * \f]
468 473 *
469 474 * where
... ... @@ -472,14 +477,18 @@ double GResponseCache::encode(const GEnergy&amp; ereco,
472 477 * - \f$E_{\rm reco}\f$ is the reconstructued energy in TeV, and
473 478 * - \f$E_{\rm true}\f$ is the true energy in TeV.
474 479 *
475   - * @todo Use instrument direction encode method
  480 + * @todo Verify unique encoding for all instruments!!!
476 481 ***************************************************************************/
477 482 double GResponseCache::encode(const GInstDir& dir,
478 483 const GEnergy& ereco,
479 484 const GEnergy& etrue) const
480 485 {
481 486 // Construct encoded instrument direction and reconstructued energy
482   - double encoded = dir.hash() + ereco.MeV() * 1.0e6 + etrue.MeV() * 1.0e8;
  487 + double encoded = dir.hash() * 1.0e5 +
  488 + (ereco.TeV() * 1.0e2 + etrue.TeV()) * m_energy_scale;
  489 + //double encoded = dir.hash() * 1.0e5 +
  490 + // ereco.TeV() * 1.0e2 +
  491 + // etrue.TeV();
483 492  
484 493 // Return encoded value
485 494 return encoded;
... ...
src/template/GXXXResponse.cpp
  View file @d511521
... ... @@ -219,7 +219,7 @@ double GXXXResponse::irf(const GEvent&amp; event,
219 219 throw GException::invalid_argument(G_IRF, msg);
220 220 }
221 221  
222   - // TODO: Now comes the teal magic, this is your job !!!
  222 + // TODO: Now comes the real magic, this is your job !!!
223 223  
224 224 // Initialise IRF
225 225 double irf = 0.0;
... ... @@ -244,56 +244,6 @@ double GXXXResponse::irf(const GEvent&amp; event,
244 244  
245 245  
246 246 /***********************************************************************//**
247   - * @brief Return value of [INSTRUMENT] instrument response for a source
248   - *
249   - * @param[in] event Event.
250   - * @param[in] source Source.
251   - * @param[in] obs Observation.
252   - * @return Instrument response $\f(cm^2 sr^{-1})$\f
253   - *
254   - * @exception GException::feature_not_implemented
255   - * Source model is not supported.
256   - *
257   - * Returns the instrument response for a given event, source and observation.
258   - *
259   - * @todo Implement whatever is needed here.
260   - ***************************************************************************/
261   -double GXXXResponse::irf(const GEvent& event,
262   - const GSource& source,
263   - const GObservation& obs) const
264   -{
265   - // Initialise IRF value
266   - double irf = 0.0;
267   -
268   - // Select IRF depending on the spatial model type
269   - switch (source.model()->code()) {
270   - case GMODEL_SPATIAL_POINT_SOURCE:
271   - {
272   - const GModelSpatialPointSource* src =
273   - static_cast<const GModelSpatialPointSource*>(source.model());
274   - GPhoton photon(src->dir(), source.energy(), source.time());
275   - irf = this->irf(event, photon, obs);
276   - }
277   - break;
278   - case GMODEL_SPATIAL_RADIAL:
279   - case GMODEL_SPATIAL_ELLIPTICAL:
280   - case GMODEL_SPATIAL_DIFFUSE:
281   - {
282   - std::string msg = "Response computation not yet implemented for "
283   - "spatial model type \""+source.model()->type()+"\".";
284   - throw GException::feature_not_implemented(G_IRF, msg);
285   - }
286   - break;
287   - default:
288   - break;
289   - }
290   -
291   - // Return IRF value
292   - return irf;
293   -}
294   -
295   -
296   -/***********************************************************************//**
297 247 * @brief Return integral of event probability for a given sky model over ROI
298 248 *
299 249 * @param[in] model Sky model.
... ...
src/template/GXXXResponse.hpp
  View file @d511521
... ... @@ -73,9 +73,6 @@ public:
73 73 virtual double irf(const GEvent& event,
74 74 const GPhoton& photon,
75 75 const GObservation& obs) const;
76   - virtual double irf(const GEvent& event,
77   - const GSource& source,
78   - const GObservation& obs) const;
79 76 virtual double nroi(const GModelSky& model,
80 77 const GEnergy& obsEng,
81 78 const GTime& obsTime,
... ...
src/template/GXXXResponse.i
  View file @d511521
... ... @@ -51,9 +51,6 @@ public:
51 51 virtual double irf(const GEvent& event,
52 52 const GPhoton& photon,
53 53 const GObservation& obs) const;
54   - virtual double irf(const GEvent& event,
55   - const GSource& source,
56   - const GObservation& obs) const;
57 54 virtual double nroi(const GModelSky& model,
58 55 const GEnergy& obsEng,
59 56 const GTime& obsTime,
... ...