Domenico Tiziani / gammalib

Commit 9a6488ba7f798b692f2de1ab39d7223b2b4b8df4

Authored by knodlseder
1 parent 015e0426

update of generic response interface

Inline Side-by-side
Showing 10 changed files with 855 additions and 406 deletions
include/GResponse.hpp
  View file @9a6488b
@@ -24,9 +24,15 @@ @@ -24,9 +24,15 @@
24 #include <iostream> 24 #include <iostream>
25 #include "GLog.hpp" 25 #include "GLog.hpp"
26 #include "GEvent.hpp" 26 #include "GEvent.hpp"
  27 +#include "GSkyDir.hpp"
27 #include "GEnergy.hpp" 28 #include "GEnergy.hpp"
28 #include "GTime.hpp" 29 #include "GTime.hpp"
29 #include "GModelSky.hpp" 30 #include "GModelSky.hpp"
  31 +#include "GModelPointSource.hpp"
  32 +#include "GModelExtendedSource.hpp"
  33 +#include "GModelDiffuseSource.hpp"
  34 +#include "GIntegrand.hpp"
  35 +#include "GMatrix.hpp"
30 36
31 /* __ Forward declarations _______________________________________________ */ 37 /* __ Forward declarations _______________________________________________ */
32 class GObservation; 38 class GObservation;
@@ -66,19 +72,125 @@ public: @@ -66,19 +72,125 @@ public:
66 virtual GResponse* clone(void) const = 0; 72 virtual GResponse* clone(void) const = 0;
67 virtual bool hasedisp(void) const = 0; 73 virtual bool hasedisp(void) const = 0;
68 virtual bool hastdisp(void) const = 0; 74 virtual bool hastdisp(void) const = 0;
69 - virtual double irf(const GEvent& event, const GModelSky& model,  
70 - const GEnergy& srcEng, const GTime& srcTime, 75 + virtual double irf(const GInstDir& obsDir,
  76 + const GEnergy& obsEng,
  77 + const GTime& obsTime,
  78 + const GSkyDir& srcDir,
  79 + const GEnergy& srcEng,
  80 + const GTime& srcTime,
71 const GObservation& obs) const = 0; 81 const GObservation& obs) const = 0;
72 - virtual double npred(const GModelSky& model, const GEnergy& srcEng,  
73 - const GTime& srcTime, 82 + virtual double npred(const GSkyDir& srcDir,
  83 + const GEnergy& srcEng,
  84 + const GTime& srcTime,
74 const GObservation& obs) const = 0; 85 const GObservation& obs) const = 0;
75 virtual std::string print(void) const = 0; 86 virtual std::string print(void) const = 0;
76 87
  88 + // Virtual methods
  89 + virtual double irf(const GEvent& event,
  90 + const GModelSky& model,
  91 + const GEnergy& srcEng,
  92 + const GTime& srcTime,
  93 + const GObservation& obs) const;
  94 + virtual double irf_ptsrc(const GInstDir& obsDir,
  95 + const GEnergy& obsEng,
  96 + const GTime& obsTime,
  97 + const GModelPointSource& model,
  98 + const GEnergy& srcEng,
  99 + const GTime& srcTime,
  100 + const GObservation& obs) const;
  101 + virtual double irf_extended(const GInstDir& obsDir,
  102 + const GEnergy& obsEng,
  103 + const GTime& obsTime,
  104 + const GModelExtendedSource& model,
  105 + const GEnergy& srcEng,
  106 + const GTime& srcTime,
  107 + const GObservation& obs) const;
  108 + virtual double irf_diffuse(const GInstDir& obsDir,
  109 + const GEnergy& obsEng,
  110 + const GTime& obsTime,
  111 + const GModelDiffuseSource& model,
  112 + const GEnergy& srcEng,
  113 + const GTime& srcTime,
  114 + const GObservation& obs) const;
  115 + virtual double npred(const GModelSky& model,
  116 + const GEnergy& srcEng,
  117 + const GTime& srcTime,
  118 + const GObservation& obs) const;
  119 + virtual double npred_ptsrc(const GModelPointSource& model,
  120 + const GEnergy& srcEng,
  121 + const GTime& srcTime,
  122 + const GObservation& obs) const;
  123 + virtual double npred_extended(const GModelExtendedSource& model,
  124 + const GEnergy& srcEng,
  125 + const GTime& srcTime,
  126 + const GObservation& obs) const;
  127 + virtual double npred_diffuse(const GModelDiffuseSource& model,
  128 + const GEnergy& srcEng,
  129 + const GTime& srcTime,
  130 + const GObservation& obs) const;
  131 +
77 protected: 132 protected:
78 // Protected methods 133 // Protected methods
79 void init_members(void); 134 void init_members(void);
80 void copy_members(const GResponse& rsp); 135 void copy_members(const GResponse& rsp);
81 void free_members(void); 136 void free_members(void);
  137 +
  138 + // Npred theta integration kernel
  139 + class npred_kern_theta : public GIntegrand {
  140 + public:
  141 + npred_kern_theta(const GResponse* rsp,
  142 + const GModelRadial* radial,
  143 + const GEnergy* srcEng,
  144 + const GTime* srcTime,
  145 + const GObservation* obs,
  146 + const GMatrix* rot) :
  147 + m_rsp(rsp),
  148 + m_radial(radial),
  149 + m_srcEng(srcEng),
  150 + m_srcTime(srcTime),
  151 + m_obs(obs),
  152 + m_rot(rot) { return; }
  153 + double eval(double theta);
  154 + protected:
  155 + const GResponse* m_rsp; //!< Pointer to response
  156 + const GModelRadial* m_radial; //!< Pointer to radial spatial model
  157 + const GEnergy* m_srcEng; //!< Pointer to true photon energy
  158 + const GTime* m_srcTime; //!< Pointer to true photon arrival time
  159 + const GObservation* m_obs; //!< Pointer to observation
  160 + const GMatrix* m_rot; //!< Rotation matrix
  161 + };
  162 +
  163 + // Npred phi integration kernel
  164 + class npred_kern_phi : public GIntegrand {
  165 + public:
  166 + npred_kern_phi(const GResponse* rsp,
  167 + const GEnergy* srcEng,
  168 + const GTime* srcTime,
  169 + const GObservation* obs,
  170 + const GMatrix* rot,
  171 + double theta,
  172 + double sin_theta) :
  173 + m_rsp(rsp),
  174 + m_srcEng(srcEng),
  175 + m_srcTime(srcTime),
  176 + m_obs(obs),
  177 + m_rot(rot),
  178 + m_theta(theta),
  179 + m_cos_theta(std::cos(theta)),
  180 + m_sin_theta(sin_theta) { return; }
  181 + double eval(double phi);
  182 + protected:
  183 + const GResponse* m_rsp; //!< Pointer to response
  184 + const GModelRadial* m_radial; //!< Pointer to radial spatial model
  185 + const GEnergy* m_srcEng; //!< Pointer to true photon energy
  186 + const GTime* m_srcTime; //!< Pointer to true photon arrival time
  187 + const GObservation* m_obs; //!< Pointer to observation
  188 + const GMatrix* m_rot; //!< Rotation matrix
  189 + double m_theta; //!< Offset angle (radians)
  190 + double m_cos_theta; //!< cosine of offset angle
  191 + double m_sin_theta; //!< Sine of offset angle
  192 + };
  193 +
82 }; 194 };
83 195
84 #endif /* GRESPONSE_HPP */ 196 #endif /* GRESPONSE_HPP */
inst/cta/include/GCTAResponse.hpp
  View file @9a6488b
@@ -12,7 +12,7 @@ @@ -12,7 +12,7 @@
12 ***************************************************************************/ 12 ***************************************************************************/
13 /** 13 /**
14 * @file GCTAResponse.hpp 14 * @file GCTAResponse.hpp
15 - * @brief GCTAResponse class interface definition. 15 + * @brief CTA instrument response function class interface definition
16 * @author J. Knodlseder 16 * @author J. Knodlseder
17 */ 17 */
18 18
@@ -73,13 +73,35 @@ public: @@ -73,13 +73,35 @@ public:
73 virtual GCTAResponse* clone(void) const; 73 virtual GCTAResponse* clone(void) const;
74 virtual bool hasedisp(void) const { return false; } 74 virtual bool hasedisp(void) const { return false; }
75 virtual bool hastdisp(void) const { return false; } 75 virtual bool hastdisp(void) const { return false; }
76 - virtual double irf(const GEvent& event, const GModelSky& model,  
77 - const GEnergy& srcEng, const GTime& srcTime, 76 + virtual double irf(const GInstDir& obsDir,
  77 + const GEnergy& obsEng,
  78 + const GTime& obsTime,
  79 + const GSkyDir& srcDir,
  80 + const GEnergy& srcEng,
  81 + const GTime& srcTime,
78 const GObservation& obs) const; 82 const GObservation& obs) const;
79 - virtual double npred(const GModelSky& model, const GEnergy& srcEng,  
80 - const GTime& srcTime, const GObservation& obs) const; 83 + virtual double npred(const GSkyDir& srcDir,
  84 + const GEnergy& srcEng,
  85 + const GTime& srcTime,
  86 + const GObservation& obs) const;
81 virtual std::string print(void) const; 87 virtual std::string print(void) const;
82 88
  89 + // Overload virtual base class methods
  90 + virtual double irf_extended(const GInstDir& obsDir,
  91 + const GEnergy& obsEng,
  92 + const GTime& obsTime,
  93 + const GModelExtendedSource& model,
  94 + const GEnergy& srcEng,
  95 + const GTime& srcTime,
  96 + const GObservation& obs) const;
  97 + virtual double irf_diffuse(const GInstDir& obsDir,
  98 + const GEnergy& obsEng,
  99 + const GTime& obsTime,
  100 + const GModelDiffuseSource& model,
  101 + const GEnergy& srcEng,
  102 + const GTime& srcTime,
  103 + const GObservation& obs) const;
  104 +
83 // Other Methods 105 // Other Methods
84 GCTAEventAtom* mc(const double& area, const GPhoton& photon, 106 GCTAEventAtom* mc(const double& area, const GPhoton& photon,
85 const GPointing& pnt, GRan& ran) const; 107 const GPointing& pnt, GRan& ran) const;
@@ -89,41 +111,6 @@ public: @@ -89,41 +111,6 @@ public:
89 void eps(const double& eps) { m_eps=eps; } 111 void eps(const double& eps) { m_eps=eps; }
90 const double& eps(void) const { return m_eps; } 112 const double& eps(void) const { return m_eps; }
91 113
92 - // Model type dependent response methods  
93 - double irf_ptsrc(const GCTAInstDir& obsDir,  
94 - const GEnergy& obsEng,  
95 - const GTime& obsTime,  
96 - const GModelPointSource& model,  
97 - const GEnergy& srcEng,  
98 - const GTime& srcTime,  
99 - const GCTAObservation& obs) const;  
100 - double irf_extended(const GCTAInstDir& obsDir,  
101 - const GEnergy& obsEng,  
102 - const GTime& obsTime,  
103 - const GModelExtendedSource& model,  
104 - const GEnergy& srcEng,  
105 - const GTime& srcTime,  
106 - const GCTAObservation& obs) const;  
107 - double irf_diffuse(const GCTAInstDir& obsDir,  
108 - const GEnergy& obsEng,  
109 - const GTime& obsTime,  
110 - const GModelDiffuseSource& model,  
111 - const GEnergy& srcEng,  
112 - const GTime& srcTime,  
113 - const GCTAObservation& obs) const;  
114 - double npred_ptsrc(const GModelPointSource& model,  
115 - const GEnergy& srcEng,  
116 - const GTime& srcTime,  
117 - const GCTAObservation& obs) const;  
118 - double npred_extended(const GModelExtendedSource& model,  
119 - const GEnergy& srcEng,  
120 - const GTime& srcTime,  
121 - const GCTAObservation& obs) const;  
122 - double npred_diffuse(const GModelDiffuseSource& model,  
123 - const GEnergy& srcEng,  
124 - const GTime& srcTime,  
125 - const GCTAObservation& obs) const;  
126 -  
127 // Low-level response methods 114 // Low-level response methods
128 double aeff(const double& theta, 115 double aeff(const double& theta,
129 const double& phi, 116 const double& phi,
inst/cta/pyext/GCTAResponse.i
  View file @9a6488b
@@ -27,6 +27,7 @@ @@ -27,6 +27,7 @@
27 * @brief CTA instrument response function 27 * @brief CTA instrument response function
28 ***************************************************************************/ 28 ***************************************************************************/
29 class GCTAResponse : public GResponse { 29 class GCTAResponse : public GResponse {
  30 +
30 public: 31 public:
31 // Constructors and destructors 32 // Constructors and destructors
32 GCTAResponse(void); 33 GCTAResponse(void);
@@ -39,11 +40,33 @@ public: @@ -39,11 +40,33 @@ public:
39 virtual GCTAResponse* clone(void) const; 40 virtual GCTAResponse* clone(void) const;
40 virtual bool hasedisp(void) const; 41 virtual bool hasedisp(void) const;
41 virtual bool hastdisp(void) const; 42 virtual bool hastdisp(void) const;
42 - virtual double irf(const GEvent& event, const GModelSky& model,  
43 - const GEnergy& srcEng, const GTime& srcTime, 43 + virtual double irf(const GInstDir& obsDir,
  44 + const GEnergy& obsEng,
  45 + const GTime& obsTime,
  46 + const GSkyDir& srcDir,
  47 + const GEnergy& srcEng,
  48 + const GTime& srcTime,
44 const GObservation& obs) const; 49 const GObservation& obs) const;
45 - virtual double npred(const GModelSky& model, const GEnergy& srcEng,  
46 - const GTime& srcTime, const GObservation& obs) const; 50 + virtual double npred(const GSkyDir& srcDir,
  51 + const GEnergy& srcEng,
  52 + const GTime& srcTime,
  53 + const GObservation& obs) const;
  54 +
  55 + // Overload virtual base class methods
  56 + virtual double irf_extended(const GInstDir& obsDir,
  57 + const GEnergy& obsEng,
  58 + const GTime& obsTime,
  59 + const GModelExtendedSource& model,
  60 + const GEnergy& srcEng,
  61 + const GTime& srcTime,
  62 + const GObservation& obs) const;
  63 + virtual double irf_diffuse(const GInstDir& obsDir,
  64 + const GEnergy& obsEng,
  65 + const GTime& obsTime,
  66 + const GModelDiffuseSource& model,
  67 + const GEnergy& srcEng,
  68 + const GTime& srcTime,
  69 + const GObservation& obs) const;
47 70
48 // Other Methods 71 // Other Methods
49 GCTAEventAtom* mc(const double& area, const GPhoton& photon, 72 GCTAEventAtom* mc(const double& area, const GPhoton& photon,
@@ -54,41 +77,6 @@ public: @@ -54,41 +77,6 @@ public:
54 void eps(const double& eps); 77 void eps(const double& eps);
55 const double& eps(void) const; 78 const double& eps(void) const;
56 79
57 - // Model type dependent response methods  
58 - double irf_ptsrc(const GCTAInstDir& obsDir,  
59 - const GEnergy& obsEng,  
60 - const GTime& obsTime,  
61 - const GModelPointSource& model,  
62 - const GEnergy& srcEng,  
63 - const GTime& srcTime,  
64 - const GCTAObservation& obs) const;  
65 - double irf_extended(const GCTAInstDir& obsDir,  
66 - const GEnergy& obsEng,  
67 - const GTime& obsTime,  
68 - const GModelExtendedSource& model,  
69 - const GEnergy& srcEng,  
70 - const GTime& srcTime,  
71 - const GCTAObservation& obs) const;  
72 - double irf_diffuse(const GCTAInstDir& obsDir,  
73 - const GEnergy& obsEng,  
74 - const GTime& obsTime,  
75 - const GModelDiffuseSource& model,  
76 - const GEnergy& srcEng,  
77 - const GTime& srcTime,  
78 - const GCTAObservation& obs) const;  
79 - double npred_ptsrc(const GModelPointSource& model,  
80 - const GEnergy& srcEng,  
81 - const GTime& srcTime,  
82 - const GCTAObservation& obs) const;  
83 - double npred_extended(const GModelExtendedSource& model,  
84 - const GEnergy& srcEng,  
85 - const GTime& srcTime,  
86 - const GCTAObservation& obs) const;  
87 - double npred_diffuse(const GModelDiffuseSource& model,  
88 - const GEnergy& srcEng,  
89 - const GTime& srcTime,  
90 - const GCTAObservation& obs) const;  
91 -  
92 // Low-level response methods 80 // Low-level response methods
93 double aeff(const double& theta, 81 double aeff(const double& theta,
94 const double& phi, 82 const double& phi,
@@ -126,6 +114,7 @@ public: @@ -126,6 +114,7 @@ public:
126 // Analytical PSF implementation 114 // Analytical PSF implementation
127 double psf_dummy(const double& delta, const double& sigma) const; 115 double psf_dummy(const double& delta, const double& sigma) const;
128 double psf_dummy_sigma(const double& srcLogEng) const; 116 double psf_dummy_sigma(const double& srcLogEng) const;
  117 + double psf_dummy_max(const double& sigma) const;
129 }; 118 };
130 119
131 120
inst/cta/src/GCTAResponse.cpp
  View file @9a6488b
@@ -40,20 +40,14 @@ @@ -40,20 +40,14 @@
40 40
41 /* __ Method name definitions ____________________________________________ */ 41 /* __ Method name definitions ____________________________________________ */
42 #define G_CALDB "GCTAResponse::caldb(std::string&)" 42 #define G_CALDB "GCTAResponse::caldb(std::string&)"
43 -#define G_IRF "GCTAResponse::irf(GEvent&, GModelSky&, GEnergy&, GTime&,"\  
44 - " GObservation&)"  
45 -#define G_NPRED "GCTAResponse::npred(GModelSky&, GEnergy&, GTime&,"\ 43 +#define G_IRF "GCTAResponse::irf(GInstDir&, GEnergy&, GTime&, GSkyDir&,"\
  44 + " GEnergy&, GTime&, GObservation&)"
  45 +#define G_NPRED "GCTAResponse::npred(GSkyDir&, GEnergy&, GTime&,"\
46 " GObservation&)" 46 " GObservation&)"
47 #define G_IRF_EXTENDED "GCTAResponse::irf_extended(GInstDir&, GEnergy&,"\ 47 #define G_IRF_EXTENDED "GCTAResponse::irf_extended(GInstDir&, GEnergy&,"\
48 " GTime&, GModelExtendedSource&, GEnergy&, GTime&, GObservation&)" 48 " GTime&, GModelExtendedSource&, GEnergy&, GTime&, GObservation&)"
49 #define G_IRF_DIFFUSE "GCTAResponse::irf_diffuse(GCTAInstDir&, GEnergy&,"\ 49 #define G_IRF_DIFFUSE "GCTAResponse::irf_diffuse(GCTAInstDir&, GEnergy&,"\
50 " GTime&, GModelDiffuseSource&, GEnergy&, GTime&, GCTAObservation&)" 50 " GTime&, GModelDiffuseSource&, GEnergy&, GTime&, GCTAObservation&)"
51 -#define G_NPRED_PTSRC "GCTAResponse::npred_ptsrc("\  
52 - "GModelPointSource&, GEnergy&, GTime&, GCTAObservation&)"  
53 -#define G_NPRED_EXTENDED "GCTAResponse::npred_extended("\  
54 - "GModelExtendedSource&, GEnergy&, GTime&, GCTAObservation&)"  
55 -#define G_NPRED_DIFFUSE "GCTAResponse::npred_diffuse("\  
56 - "GModelDiffuseSource&, GEnergy&, GTime&, GCTAObservation&)"  
57 #define G_READ "GCTAResponse::read_performance_table(std::string&)" 51 #define G_READ "GCTAResponse::read_performance_table(std::string&)"
58 52
59 /* __ Macros _____________________________________________________________ */ 53 /* __ Macros _____________________________________________________________ */
@@ -213,53 +207,70 @@ GCTAResponse* GCTAResponse::clone(void) const @@ -213,53 +207,70 @@ GCTAResponse* GCTAResponse::clone(void) const
213 /***********************************************************************//** 207 /***********************************************************************//**
214 * @brief Return value of instrument response function 208 * @brief Return value of instrument response function
215 * 209 *
216 - * @param[in] event Event.  
217 - * @param[in] model Source model. 210 + * @param[in] obsDir Observed photon direction.
  211 + * @param[in] obsEng Observed energy of photon.
  212 + * @param[in] obsTime Observed photon arrival time.
  213 + * @param[in] srcDir True photon arrival direction.
218 * @param[in] srcEng True energy of photon. 214 * @param[in] srcEng True energy of photon.
219 * @param[in] srcTime True photon arrival time. 215 * @param[in] srcTime True photon arrival time.
220 * @param[in] obs Observation. 216 * @param[in] obs Observation.
221 * 217 *
222 - * @exception GCTAException::bad_observation_type  
223 - * Observation is not a CTA observation.  
224 * @exception GCTAException::bad_instdir_type 218 * @exception GCTAException::bad_instdir_type
225 - * Instrument direction is not a CTA instrument direction.  
226 - *  
227 - * @todo Throw exception if model pointer is invalid. 219 + * Instrument direction is not a valid CTA instrument direction.
228 ***************************************************************************/ 220 ***************************************************************************/
229 -double GCTAResponse::irf(const GEvent& event,  
230 - const GModelSky& model, 221 +double GCTAResponse::irf(const GInstDir& obsDir,
  222 + const GEnergy& obsEng,
  223 + const GTime& obsTime,
  224 + const GSkyDir& srcDir,
231 const GEnergy& srcEng, 225 const GEnergy& srcEng,
232 const GTime& srcTime, 226 const GTime& srcTime,
233 const GObservation& obs) const 227 const GObservation& obs) const
234 { 228 {
  229 + // Get pointer on CTA instrument direction
  230 + const GCTAInstDir* dir = dynamic_cast<const GCTAInstDir*>(&obsDir);
  231 + if (dir == NULL)
  232 + throw GCTAException::bad_instdir_type(G_IRF);
  233 +
  234 + // Get pointing direction zenith angle and azimuth
  235 + double zenith = 0.0;
  236 + double azimuth = 0.0;
  237 +
  238 + // Get radial offset and polar angles in camera
  239 + double theta = 0.0;
  240 + double phi = 0.0;
  241 +
  242 + // Get log10(E/TeV) of true photon energy.
  243 + double srcLogEng = srcEng.log10TeV();
  244 +
  245 + // Determine angular separation between true and measured photon
  246 + // direction in radians
  247 + double delta = dir->dist(srcDir);
  248 +
  249 + // Get maximum angular separation for which PSF is significant
  250 + double delta_max = psf_delta_max(theta, phi, zenith, azimuth, srcLogEng);
  251 +
235 // Initialise IRF value 252 // Initialise IRF value
236 double irf = 0.0; 253 double irf = 0.0;
237 254
238 - // Get pointer on CTA observation  
239 - const GCTAObservation* ctaobs = dynamic_cast<const GCTAObservation*>(&obs);  
240 - if (ctaobs == NULL)  
241 - throw GCTAException::bad_observation_type(G_IRF); 255 + // Compute only if we're sufficiently close to PSF
  256 + if (delta <= delta_max) {
242 257
243 - // Get pointer on CTA instrument direction  
244 - const GCTAInstDir* ctadir = dynamic_cast<const GCTAInstDir*>(&event.dir());  
245 - if (ctadir == NULL)  
246 - throw GCTAException::bad_instdir_type(G_IRF); 258 + // Get point source IRF
  259 + irf = aeff(theta, phi, zenith, azimuth, srcLogEng);
  260 + irf *= psf(delta, theta, phi, zenith, azimuth, srcLogEng);
247 261
248 - // Get model pointers  
249 - const GModelPointSource* ptsrc = dynamic_cast<const GModelPointSource*>(&model);  
250 - const GModelExtendedSource* extsrc = dynamic_cast<const GModelExtendedSource*>(&model);  
251 - const GModelDiffuseSource* difsrc = dynamic_cast<const GModelDiffuseSource*>(&model);  
252 -  
253 - // Call model dependent method  
254 - if (ptsrc != NULL)  
255 - irf = irf_ptsrc(*ctadir, event.energy(), event.time(),  
256 - *ptsrc, srcEng, srcTime, *ctaobs);  
257 - else if (extsrc != NULL)  
258 - irf = irf_extended(*ctadir, event.energy(), event.time(),  
259 - *extsrc, srcEng, srcTime, *ctaobs);  
260 - else if (difsrc != NULL)  
261 - irf = irf_diffuse(*ctadir, event.energy(), event.time(),  
262 - *difsrc, srcEng, srcTime, *ctaobs); 262 + // Multiply-in energy dispersion
  263 + if (hasedisp()) {
  264 +
  265 + // Get log10(E/TeV) of measured photon energy.
  266 + double obsLogEng = obsEng.log10TeV();
  267 +
  268 + // Multiply-in energy dispersion
  269 + irf *= edisp(obsLogEng, theta, phi, zenith, azimuth, srcLogEng);
  270 +
  271 + } // endif: energy dispersion was available
  272 +
  273 + } // endif: we were sufficiently close to PSF
263 274
264 // Return IRF value 275 // Return IRF value
265 return irf; 276 return irf;
@@ -267,45 +278,62 @@ double GCTAResponse::irf(const GEvent&amp; event, @@ -267,45 +278,62 @@ double GCTAResponse::irf(const GEvent&amp; event,
267 278
268 279
269 /***********************************************************************//** 280 /***********************************************************************//**
270 - * @brief Return spatial integral of instrument response function 281 + * @brief Return spatial integral of point spread function
271 * 282 *
272 - * @param[in] model Source model. 283 + * @param[in] srcDir True photon arrival direction.
273 * @param[in] srcEng True energy of photon. 284 * @param[in] srcEng True energy of photon.
274 * @param[in] srcTime True photon arrival time. 285 * @param[in] srcTime True photon arrival time.
275 * @param[in] obs Observation. 286 * @param[in] obs Observation.
276 * 287 *
277 * @exception GCTAException::bad_observation_type 288 * @exception GCTAException::bad_observation_type
278 * Observation is not a CTA observation. 289 * Observation is not a CTA observation.
279 - *  
280 - * @todo Throw exception if model pointer is invalid. 290 + * @exception GException::no_list
  291 + * Observation does not contain a valid CTA event list.
281 ***************************************************************************/ 292 ***************************************************************************/
282 -double GCTAResponse::npred(const GModelSky& model, 293 +double GCTAResponse::npred(const GSkyDir& srcDir,
283 const GEnergy& srcEng, 294 const GEnergy& srcEng,
284 const GTime& srcTime, 295 const GTime& srcTime,
285 const GObservation& obs) const 296 const GObservation& obs) const
286 { 297 {
287 - // Initialise Npred value  
288 - double npred = 0.0;  
289 -  
290 // Get pointer on CTA observation 298 // Get pointer on CTA observation
291 const GCTAObservation* ctaobs = dynamic_cast<const GCTAObservation*>(&obs); 299 const GCTAObservation* ctaobs = dynamic_cast<const GCTAObservation*>(&obs);
292 if (ctaobs == NULL) 300 if (ctaobs == NULL)
293 throw GCTAException::bad_observation_type(G_NPRED); 301 throw GCTAException::bad_observation_type(G_NPRED);
294 302
295 - // Get model pointers  
296 - const GModelPointSource* ptsrc = dynamic_cast<const GModelPointSource*>(&model);  
297 - const GModelExtendedSource* extsrc = dynamic_cast<const GModelExtendedSource*>(&model);  
298 - const GModelDiffuseSource* difsrc = dynamic_cast<const GModelDiffuseSource*>(&model); 303 + // Get pointer on CTA events list
  304 + const GCTAEventList* events = dynamic_cast<const GCTAEventList*>(ctaobs->events());
  305 + if (events == NULL)
  306 + throw GException::no_list(G_NPRED);
299 307
300 - // Call model dependent method  
301 - if (ptsrc != NULL)  
302 - npred = npred_ptsrc(*ptsrc, srcEng, srcTime, *ctaobs);  
303 - else if (extsrc != NULL)  
304 - npred = npred_extended(*extsrc, srcEng, srcTime, *ctaobs);  
305 - else if (difsrc != NULL)  
306 - npred = npred_diffuse(*difsrc, srcEng, srcTime, *ctaobs); 308 + // Get pointer on CTA pointing
  309 + const GCTAPointing *pnt = ctaobs->pointing(srcTime);
307 310
308 - // Return response value 311 + // Get pointing direction zenith angle and azimuth
  312 + // @todo Not yet implemented (and not yet needed)
  313 + double zenith = 0.0;
  314 + double azimuth = 0.0;
  315 +
  316 + // Get radial offset and polar angles in camera
  317 + // @todo Not yet implemented (and not yet needed)
  318 + double theta = 0.0;
  319 + double phi = 0.0;
  320 +
  321 + // Get log10(E/TeV) of true photon energy.
  322 + double srcLogEng = srcEng.log10TeV();
  323 +
  324 + // Get IRF components
  325 + double npred = aeff(theta, phi, zenith, azimuth, srcLogEng);
  326 + npred *= npsf(srcDir, srcLogEng, srcTime, *pnt, events->roi());
  327 +
  328 + // Multiply-in energy dispersion
  329 + if (hasedisp()) {
  330 +
  331 + // Multiply-in energy dispersion
  332 + npred *= nedisp(srcDir, srcEng, srcTime, *pnt, events->ebounds());
  333 +
  334 + } // endif: had energy dispersion
  335 +
  336 + // Return Npred
309 return npred; 337 return npred;
310 } 338 }
311 339
@@ -470,81 +498,6 @@ std::string GCTAResponse::print(void) const @@ -470,81 +498,6 @@ std::string GCTAResponse::print(void) const
470 ==========================================================================*/ 498 ==========================================================================*/
471 499
472 /***********************************************************************//** 500 /***********************************************************************//**
473 - * @brief Return value of point source instrument response function  
474 - *  
475 - * @param[in] obsDir Observed photon direction.  
476 - * @param[in] obsEng Observed energy of photon.  
477 - * @param[in] obsTime Observed photon arrival time.  
478 - * @param[in] model Point source model.  
479 - * @param[in] srcEng True energy of photon.  
480 - * @param[in] srcTime True photon arrival time.  
481 - * @param[in] obs CTA Observation.  
482 - *  
483 - * @todo The telescope zenith and azimuth angles as well as the radial offset  
484 - * and polar angles in the camera are not yet computed correctly (as  
485 - * the actual IRF implement does not need these values).  
486 - ***************************************************************************/  
487 -double GCTAResponse::irf_ptsrc(const GCTAInstDir& obsDir,  
488 - const GEnergy& obsEng,  
489 - const GTime& obsTime,  
490 - const GModelPointSource& model,  
491 - const GEnergy& srcEng,  
492 - const GTime& srcTime,  
493 - const GCTAObservation& obs) const  
494 -{  
495 - // Get point source location  
496 - GSkyDir srcDir = model.dir();  
497 -  
498 - // Get CTA pointing  
499 - const GCTAPointing* pnt = obs.pointing(srcTime);  
500 -  
501 - // Get pointing direction zenith angle and azimuth  
502 - double zenith = 0.0;  
503 - double azimuth = 0.0;  
504 -  
505 - // Get radial offset and polar angles in camera  
506 - double theta = 0.0;  
507 - double phi = 0.0;  
508 -  
509 - // Get log10(E/TeV) of true photon energy.  
510 - double srcLogEng = srcEng.log10TeV();  
511 -  
512 - // Determine angular separation between true and measured photon  
513 - // direction in radians  
514 - double delta = obsDir.dist(srcDir);  
515 -  
516 - // Get maximum angular separation for which PSF is significant  
517 - double delta_max = psf_delta_max(theta, phi, zenith, azimuth, srcLogEng);  
518 -  
519 - // Initialise IRF value  
520 - double irf = 0.0;  
521 -  
522 - // Compute only if we're sufficiently close to PSF  
523 - if (delta <= delta_max) {  
524 -  
525 - // Get point source IRF  
526 - irf = aeff(theta, phi, zenith, azimuth, srcLogEng);  
527 - irf *= psf(delta, theta, phi, zenith, azimuth, srcLogEng);  
528 -  
529 - // Multiply-in energy dispersion  
530 - if (hasedisp()) {  
531 -  
532 - // Get log10(E/TeV) of measured photon energy.  
533 - double obsLogEng = obsEng.log10TeV();  
534 -  
535 - // Multiply-in energy dispersion  
536 - irf *= edisp(obsLogEng, theta, phi, zenith, azimuth, srcLogEng);  
537 -  
538 - } // endif: energy dispersion was available  
539 -  
540 - } // endif: we were sufficiently close to PSF  
541 -  
542 - // Return IRF value  
543 - return irf;  
544 -}  
545 -  
546 -  
547 -/***********************************************************************//**  
548 * @brief Return value of extended source instrument response function 501 * @brief Return value of extended source instrument response function
549 * 502 *
550 * @param[in] obsDir Observed photon direction. 503 * @param[in] obsDir Observed photon direction.
@@ -553,7 +506,12 @@ double GCTAResponse::irf_ptsrc(const GCTAInstDir&amp; obsDir, @@ -553,7 +506,12 @@ double GCTAResponse::irf_ptsrc(const GCTAInstDir&amp; obsDir,
553 * @param[in] model Extended source model. 506 * @param[in] model Extended source model.
554 * @param[in] srcEng True energy of photon. 507 * @param[in] srcEng True energy of photon.
555 * @param[in] srcTime True photon arrival time. 508 * @param[in] srcTime True photon arrival time.
556 - * @param[in] obs CTA Observation. 509 + * @param[in] obs Observation.
  510 + *
  511 + * @exception GCTAException::bad_observation_type
  512 + * Observation is not a valid CTA observation
  513 + * @exception GCTAException::bad_instdir_type
  514 + * Instrument direction is not a valid CTA instrument direction
557 * 515 *
558 * @todo The telescope zenith and azimuth angles as well as the radial offset 516 * @todo The telescope zenith and azimuth angles as well as the radial offset
559 * and polar angles in the camera are not yet computed correctly (as 517 * and polar angles in the camera are not yet computed correctly (as
@@ -562,16 +520,26 @@ double GCTAResponse::irf_ptsrc(const GCTAInstDir&amp; obsDir, @@ -562,16 +520,26 @@ double GCTAResponse::irf_ptsrc(const GCTAInstDir&amp; obsDir,
562 * @todo The theta limits could still improved as only the theta range is 520 * @todo The theta limits could still improved as only the theta range is
563 * relevant for which the source and the PSF circle overlap. 521 * relevant for which the source and the PSF circle overlap.
564 ***************************************************************************/ 522 ***************************************************************************/
565 -double GCTAResponse::irf_extended(const GCTAInstDir& obsDir, 523 +double GCTAResponse::irf_extended(const GInstDir& obsDir,
566 const GEnergy& obsEng, 524 const GEnergy& obsEng,
567 const GTime& obsTime, 525 const GTime& obsTime,
568 const GModelExtendedSource& model, 526 const GModelExtendedSource& model,
569 const GEnergy& srcEng, 527 const GEnergy& srcEng,
570 const GTime& srcTime, 528 const GTime& srcTime,
571 - const GCTAObservation& obs) const 529 + const GObservation& obs) const
572 { 530 {
573 - // Get CTA pointing  
574 - const GCTAPointing* pnt = obs.pointing(srcTime); 531 + // Get pointer on CTA observation
  532 + const GCTAObservation* ctaobs = dynamic_cast<const GCTAObservation*>(&obs);
  533 + if (ctaobs == NULL)
  534 + throw GCTAException::bad_observation_type(G_IRF_EXTENDED);
  535 +
  536 + // Get pointer on CTA instrument direction
  537 + const GCTAInstDir* dir = dynamic_cast<const GCTAInstDir*>(&obsDir);
  538 + if (dir == NULL)
  539 + throw GCTAException::bad_instdir_type(G_IRF_EXTENDED);
  540 +
  541 + // Get pointer on CTA pointing
  542 + const GCTAPointing *pnt = ctaobs->pointing(srcTime);
575 543
576 // Get pointing direction zenith angle and azimuth 544 // Get pointing direction zenith angle and azimuth
577 double zenith = 0.0; 545 double zenith = 0.0;
@@ -589,7 +557,7 @@ double GCTAResponse::irf_extended(const GCTAInstDir&amp; obsDir, @@ -589,7 +557,7 @@ double GCTAResponse::irf_extended(const GCTAInstDir&amp; obsDir,
589 double sigma = psf_dummy_sigma(srcLogEng); 557 double sigma = psf_dummy_sigma(srcLogEng);
590 558
591 // Determine camera direction of measured photon 559 // Determine camera direction of measured photon
592 - GCTADir obsCam(obsDir, *pnt); 560 + GCTADir obsCam(*dir, *pnt);
593 561
594 // Determine camera direction of source model centre 562 // Determine camera direction of source model centre
595 GCTADir srcCam(model.dir(), *pnt); 563 GCTADir srcCam(model.dir(), *pnt);
@@ -667,13 +635,13 @@ double GCTAResponse::irf_extended(const GCTAInstDir&amp; obsDir, @@ -667,13 +635,13 @@ double GCTAResponse::irf_extended(const GCTAInstDir&amp; obsDir,
667 * 635 *
668 * @todo Implement method. 636 * @todo Implement method.
669 ***************************************************************************/ 637 ***************************************************************************/
670 -double GCTAResponse::irf_diffuse(const GCTAInstDir& obsDir, 638 +double GCTAResponse::irf_diffuse(const GInstDir& obsDir,
671 const GEnergy& obsEng, 639 const GEnergy& obsEng,
672 const GTime& obsTime, 640 const GTime& obsTime,
673 const GModelDiffuseSource& model, 641 const GModelDiffuseSource& model,
674 const GEnergy& srcEng, 642 const GEnergy& srcEng,
675 const GTime& srcTime, 643 const GTime& srcTime,
676 - const GCTAObservation& obs) const 644 + const GObservation& obs) const
677 { 645 {
678 // Feature not yet implemented 646 // Feature not yet implemented
679 throw GException::feature_not_implemented(G_IRF_DIFFUSE, 647 throw GException::feature_not_implemented(G_IRF_DIFFUSE,
@@ -684,113 +652,6 @@ double GCTAResponse::irf_diffuse(const GCTAInstDir&amp; obsDir, @@ -684,113 +652,6 @@ double GCTAResponse::irf_diffuse(const GCTAInstDir&amp; obsDir,
684 } 652 }
685 653
686 654
687 -/***********************************************************************//**  
688 - * @brief Return spatial integral of point source model  
689 - *  
690 - * @param[in] model Point source model.  
691 - * @param[in] srcEng True energy of photon.  
692 - * @param[in] srcTime True photon arrival time.  
693 - * @param[in] obs Observation.  
694 - *  
695 - * @exception GException::no_list  
696 - * Observation does not contain an event list.  
697 - ***************************************************************************/  
698 -double GCTAResponse::npred_ptsrc(const GModelPointSource& model,  
699 - const GEnergy& srcEng,  
700 - const GTime& srcTime,  
701 - const GCTAObservation& obs) const  
702 -{  
703 - // Get point source location  
704 - GSkyDir srcDir = model.dir();  
705 -  
706 - // Get pointer on CTA events list  
707 - const GCTAEventList* ptr = dynamic_cast<const GCTAEventList*>(obs.events());  
708 - if (ptr == NULL)  
709 - throw GException::no_list(G_NPRED_PTSRC);  
710 -  
711 - // Get pointing direction zenith angle and azimuth  
712 - const GCTAPointing *pnt = obs.pointing(srcTime);  
713 - double zenith = 0.0;  
714 - double azimuth = 0.0;  
715 -  
716 - // Get radial offset and polar angles in camera  
717 - double theta = 0.0;  
718 - double phi = 0.0;  
719 -  
720 - // Get log10(E/TeV) of true photon energy.  
721 - double srcLogEng = srcEng.log10TeV();  
722 -  
723 - // Get IRF components  
724 - double nirf = aeff(theta, phi, zenith, azimuth, srcLogEng);  
725 - nirf *= npsf(srcDir, srcLogEng, srcTime, *pnt, ptr->roi());  
726 -  
727 - // Multiply-in energy dispersion  
728 - if (hasedisp()) {  
729 -  
730 - // Multiply-in energy dispersion  
731 - nirf *= nedisp(srcDir, srcEng, srcTime, *pnt, ptr->ebounds());  
732 -  
733 - } // endif: had energy dispersion  
734 -  
735 - // Return integrated IRF value  
736 - return nirf;  
737 -}  
738 -  
739 -  
740 -/***********************************************************************//**  
741 - * @brief Return spatial integral of extended source model  
742 - *  
743 - * @param[in] model Extended source model.  
744 - * @param[in] srcEng True energy of photon.  
745 - * @param[in] srcTime True photon arrival time.  
746 - * @param[in] obs Observation.  
747 - *  
748 - * @exception GException::feature_not_implemented  
749 - * Method not yet implemented.  
750 - *  
751 - * @todo Implement method.  
752 - ***************************************************************************/  
753 -double GCTAResponse::npred_extended(const GModelExtendedSource& model,  
754 - const GEnergy& srcEng,  
755 - const GTime& srcTime,  
756 - const GCTAObservation& obs) const  
757 -{  
758 - // Feature not yet implemented  
759 - throw GException::feature_not_implemented(G_NPRED_EXTENDED,  
760 - "Method for extended source not yet implemented.");  
761 -  
762 - // Return integrated IRF value  
763 - return 0.0;  
764 -}  
765 -  
766 -  
767 -/***********************************************************************//**  
768 - * @brief Return spatial integral of diffuse source model  
769 - *  
770 - * @param[in] model Diffuse source model.  
771 - * @param[in] srcEng True energy of photon.  
772 - * @param[in] srcTime True photon arrival time.  
773 - * @param[in] obs Observation.  
774 - *  
775 - * @exception GException::feature_not_implemented  
776 - * Method not yet implemented.  
777 - *  
778 - * @todo Implement method.  
779 - ***************************************************************************/  
780 -double GCTAResponse::npred_diffuse(const GModelDiffuseSource& model,  
781 - const GEnergy& srcEng,  
782 - const GTime& srcTime,  
783 - const GCTAObservation& obs) const  
784 -{  
785 - // Feature not yet implemented  
786 - throw GException::feature_not_implemented(G_NPRED_DIFFUSE,  
787 - "Method for extended source not yet implemented.");  
788 -  
789 - // Return integrated IRF value  
790 - return 0.0;  
791 -}  
792 -  
793 -  
794 /*========================================================================== 655 /*==========================================================================
795 = = 656 = =
796 = Low-level CTA response methods = 657 = Low-level CTA response methods =
inst/lat/include/GLATResponse.hpp
  View file @9a6488b
@@ -32,8 +32,6 @@ @@ -32,8 +32,6 @@
32 #include "GModel.hpp" 32 #include "GModel.hpp"
33 #include "GObservation.hpp" 33 #include "GObservation.hpp"
34 #include "GResponse.hpp" 34 #include "GResponse.hpp"
35 -#include "GEnergy.hpp"  
36 -#include "GTime.hpp"  
37 35
38 36
39 /***********************************************************************//** 37 /***********************************************************************//**
@@ -52,18 +50,28 @@ public: @@ -52,18 +50,28 @@ public:
52 // Operators 50 // Operators
53 GLATResponse& operator= (const GLATResponse & rsp); 51 GLATResponse& operator= (const GLATResponse & rsp);
54 52
55 - // Implement pure virtual base class methods  
56 - void clear(void);  
57 - GLATResponse* clone(void) const;  
58 - bool hasedisp(void) const { return false; }  
59 - bool hastdisp(void) const { return false; }  
60 - double irf(const GEvent& event, const GModelSky& model,  
61 - const GEnergy& srcEng, const GTime& srcTime,  
62 - const GObservation& obs) const;  
63 - double npred(const GModelSky& model, const GEnergy& srcEng,  
64 - const GTime& srcTime,  
65 - const GObservation& obs) const;  
66 - std::string print(void) const; 53 + // Implemented pure virtual methods
  54 + virtual void clear(void);
  55 + virtual GLATResponse* clone(void) const;
  56 + virtual bool hasedisp(void) const { return false; }
  57 + virtual bool hastdisp(void) const { return false; }
  58 + virtual double irf(const GInstDir& obsDir,
  59 + const GEnergy& obsEng,
  60 + const GTime& obsTime,
  61 + const GSkyDir& srcDir,
  62 + const GEnergy& srcEng,
  63 + const GTime& srcTime,
  64 + const GObservation& obs) const;
  65 + virtual double npred(const GSkyDir& srcDir,
  66 + const GEnergy& srcEng,
  67 + const GTime& srcTime,
  68 + const GObservation& obs) const;
  69 + virtual std::string print(void) const;
  70 +
  71 + // Implemented virtual methods
  72 + virtual double irf(const GEvent& event, const GModelSky& model,
  73 + const GEnergy& srcEng, const GTime& srcTime,
  74 + const GObservation& obs) const;
67 75
68 // Other Methods 76 // Other Methods
69 void caldb(const std::string& caldb); 77 void caldb(const std::string& caldb);
inst/lat/pyext/GLATResponse.i
  View file @9a6488b
1 /*************************************************************************** 1 /***************************************************************************
2 - * GLATResponse.i - Fermi LAT Response class * 2 + * GLATResponse.i - Fermi-LAT response class *
3 * ----------------------------------------------------------------------- * 3 * ----------------------------------------------------------------------- *
4 * copyright (C) 2008-2011 by Jurgen Knodlseder * 4 * copyright (C) 2008-2011 by Jurgen Knodlseder *
5 * ----------------------------------------------------------------------- * 5 * ----------------------------------------------------------------------- *
@@ -12,7 +12,7 @@ @@ -12,7 +12,7 @@
12 ***************************************************************************/ 12 ***************************************************************************/
13 /** 13 /**
14 * @file GLATResponse.i 14 * @file GLATResponse.i
15 - * @brief Fermi LAT Response class python bindings 15 + * @brief Fermi-LAT Response class Python interface definition
16 * @author J. Knodlseder 16 * @author J. Knodlseder
17 */ 17 */
18 %{ 18 %{
@@ -24,26 +24,37 @@ @@ -24,26 +24,37 @@
24 /***********************************************************************//** 24 /***********************************************************************//**
25 * @class GLATResponse 25 * @class GLATResponse
26 * 26 *
27 - * @brief Fermi LAT Response class python bindings 27 + * @brief Fermi-LAT response class
28 ***************************************************************************/ 28 ***************************************************************************/
29 class GLATResponse : public GResponse { 29 class GLATResponse : public GResponse {
  30 +
30 public: 31 public:
31 // Constructors and destructors 32 // Constructors and destructors
32 GLATResponse(void); 33 GLATResponse(void);
33 GLATResponse(const GLATResponse& rsp); 34 GLATResponse(const GLATResponse& rsp);
34 virtual ~GLATResponse(void); 35 virtual ~GLATResponse(void);
35 36
36 - // Implement pure virtual base class methods  
37 - void clear(void);  
38 - GLATResponse* clone(void) const;  
39 - bool hasedisp(void) const;  
40 - bool hastdisp(void) const;  
41 - double irf(const GEvent& event, const GModelSky& model,  
42 - const GEnergy& srcEng, const GTime& srcTime,  
43 - const GObservation& obs) const;  
44 - double npred(const GModelSky& model, const GEnergy& srcEng,  
45 - const GTime& srcTime,  
46 - const GObservation& obs) const; 37 + // Implemented pure virtual methods
  38 + virtual void clear(void);
  39 + virtual GLATResponse* clone(void) const;
  40 + virtual bool hasedisp(void) const;
  41 + virtual bool hastdisp(void) const;
  42 + virtual double irf(const GInstDir& obsDir,
  43 + const GEnergy& obsEng,
  44 + const GTime& obsTime,
  45 + const GSkyDir& srcDir,
  46 + const GEnergy& srcEng,
  47 + const GTime& srcTime,
  48 + const GObservation& obs) const;
  49 + virtual double npred(const GSkyDir& srcDir,
  50 + const GEnergy& srcEng,
  51 + const GTime& srcTime,
  52 + const GObservation& obs) const;
  53 +
  54 + // Implemented virtual methods
  55 + virtual double irf(const GEvent& event, const GModelSky& model,
  56 + const GEnergy& srcEng, const GTime& srcTime,
  57 + const GObservation& obs) const;
47 58
48 // Other Methods 59 // Other Methods
49 void caldb(const std::string& caldb); 60 void caldb(const std::string& caldb);
inst/lat/src/GLATResponse.cpp
  View file @9a6488b
1 /*************************************************************************** 1 /***************************************************************************
2 - * GLATResponse.cpp - Fermi LAT Response class * 2 + * GLATResponse.cpp - Fermi-LAT response class *
3 * ----------------------------------------------------------------------- * 3 * ----------------------------------------------------------------------- *
4 * copyright (C) 2008-2011 by Jurgen Knodlseder * 4 * copyright (C) 2008-2011 by Jurgen Knodlseder *
5 * ----------------------------------------------------------------------- * 5 * ----------------------------------------------------------------------- *
@@ -12,7 +12,7 @@ @@ -12,7 +12,7 @@
12 ***************************************************************************/ 12 ***************************************************************************/
13 /** 13 /**
14 * @file GLATResponse.cpp 14 * @file GLATResponse.cpp
15 - * @brief Fermi LAT Response class implementation. 15 + * @brief Fermi-LAT response class implementation
16 * @author J. Knodlseder 16 * @author J. Knodlseder
17 */ 17 */
18 18
@@ -23,6 +23,7 @@ @@ -23,6 +23,7 @@
23 #include <unistd.h> // access() function 23 #include <unistd.h> // access() function
24 #include <cstdlib> // std::getenv() function 24 #include <cstdlib> // std::getenv() function
25 #include <string> 25 #include <string>
  26 +#include "GLATInstDir.hpp"
26 #include "GLATResponse.hpp" 27 #include "GLATResponse.hpp"
27 #include "GLATObservation.hpp" 28 #include "GLATObservation.hpp"
28 #include "GLATEventAtom.hpp" 29 #include "GLATEventAtom.hpp"
@@ -37,6 +38,8 @@ @@ -37,6 +38,8 @@
37 /* __ Method name definitions ____________________________________________ */ 38 /* __ Method name definitions ____________________________________________ */
38 #define G_CALDB "GLATResponse::caldb(std::string&)" 39 #define G_CALDB "GLATResponse::caldb(std::string&)"
39 #define G_LOAD "GLATResponse::load(std::string&)" 40 #define G_LOAD "GLATResponse::load(std::string&)"
  41 +#define G_IRF "GLATResponse::irf(GInstDir&, GEnergy&, GTime&, GSkyDir&,"\
  42 + " GEnergy&, GTime&, GObservation&)"
40 #define G_AEFF "GLATResponse::aeff(int&)" 43 #define G_AEFF "GLATResponse::aeff(int&)"
41 #define G_PSF "GLATResponse::psf(int&)" 44 #define G_PSF "GLATResponse::psf(int&)"
42 #define G_EDISP "GLATResponse::edisp(int&)" 45 #define G_EDISP "GLATResponse::edisp(int&)"
@@ -44,7 +47,7 @@ @@ -44,7 +47,7 @@
44 "GTime&, GObservation&)" 47 "GTime&, GObservation&)"
45 #define G_IRF_BIN "GLATResponse::irf(GLATEventBin&, GModel&, GEnergy&,"\ 48 #define G_IRF_BIN "GLATResponse::irf(GLATEventBin&, GModel&, GEnergy&,"\
46 "GTime&, GObservation&)" 49 "GTime&, GObservation&)"
47 -#define G_NPRED "GLATResponse::npred(GModel&, GEnergy&, GTime&,"\ 50 +#define G_NPRED "GLATResponse::npred(GSkyDir&, GEnergy&, GTime&,"\
48 " GObservation&)" 51 " GObservation&)"
49 52
50 /* __ Macros _____________________________________________________________ */ 53 /* __ Macros _____________________________________________________________ */
@@ -69,7 +72,7 @@ @@ -69,7 +72,7 @@
69 ***************************************************************************/ 72 ***************************************************************************/
70 GLATResponse::GLATResponse(void) : GResponse() 73 GLATResponse::GLATResponse(void) : GResponse()
71 { 74 {
72 - // Initialise class members for clean destruction 75 + // Initialise members
73 init_members(); 76 init_members();
74 77
75 // Return 78 // Return
@@ -84,7 +87,7 @@ GLATResponse::GLATResponse(void) : GResponse() @@ -84,7 +87,7 @@ GLATResponse::GLATResponse(void) : GResponse()
84 ***************************************************************************/ 87 ***************************************************************************/
85 GLATResponse::GLATResponse(const GLATResponse& rsp) : GResponse(rsp) 88 GLATResponse::GLATResponse(const GLATResponse& rsp) : GResponse(rsp)
86 { 89 {
87 - // Initialise class members for clean destruction 90 + // Initialise members
88 init_members(); 91 init_members();
89 92
90 // Copy members 93 // Copy members
@@ -130,7 +133,7 @@ GLATResponse&amp; GLATResponse::operator= (const GLATResponse&amp; rsp) @@ -130,7 +133,7 @@ GLATResponse&amp; GLATResponse::operator= (const GLATResponse&amp; rsp)
130 // Free members 133 // Free members
131 free_members(); 134 free_members();
132 135
133 - // Initialise private members for clean destruction 136 + // Initialise members
134 init_members(); 137 init_members();
135 138
136 // Copy members 139 // Copy members
@@ -198,7 +201,7 @@ void GLATResponse::caldb(const std::string&amp; caldb) @@ -198,7 +201,7 @@ void GLATResponse::caldb(const std::string&amp; caldb)
198 // Check if calibration database directory is accessible 201 // Check if calibration database directory is accessible
199 if (access(caldb.c_str(), R_OK) != 0) 202 if (access(caldb.c_str(), R_OK) != 0)
200 throw GException::caldb_not_found(G_CALDB, caldb); 203 throw GException::caldb_not_found(G_CALDB, caldb);
201 - 204 +
202 // Store the path to the calibration database 205 // Store the path to the calibration database
203 m_caldb = caldb; 206 m_caldb = caldb;
204 207
@@ -208,6 +211,78 @@ void GLATResponse::caldb(const std::string&amp; caldb) @@ -208,6 +211,78 @@ void GLATResponse::caldb(const std::string&amp; caldb)
208 211
209 212
210 /***********************************************************************//** 213 /***********************************************************************//**
  214 + * @brief Return value of point source IRF
  215 + *
  216 + * @param[in] obsDir Measured photon arrival direction.
  217 + * @param[in] obsEng Measured photon energy.
  218 + * @param[in] obsTime Measured photon arrival time.
  219 + * @param[in] srcDir True photon arrival direction.
  220 + * @param[in] srcEng True energy of photon.
  221 + * @param[in] srcTime True photon arrival time.
  222 + * @param[in] obs Observation.
  223 + *
  224 + * @exception GLATException::bad_instdir_type
  225 + * Instrument direction is not a valid LAT instrument direction.
  226 + *
  227 + * @todo The IRF value is not devided by ontime of the event, but it is
  228 + * already time integrated.
  229 + ***************************************************************************/
  230 +double GLATResponse::irf(const GInstDir& obsDir,
  231 + const GEnergy& obsEng,
  232 + const GTime& obsTime,
  233 + const GSkyDir& srcDir,
  234 + const GEnergy& srcEng,
  235 + const GTime& srcTime,
  236 + const GObservation& obs) const
  237 +{
  238 + // Get LAT instrument direction
  239 + const GLATInstDir* dir = dynamic_cast<const GLATInstDir*>(&obsDir);
  240 + if (dir == NULL)
  241 + throw GLATException::bad_instdir_type(G_IRF);
  242 +
  243 + // Search for mean PSF
  244 + int ipsf = -1;
  245 + for (int i = 0; i < m_ptsrc.size(); ++i) {
  246 + if (m_ptsrc[i]->dir() == srcDir) {
  247 + ipsf = i;
  248 + break;
  249 + }
  250 + }
  251 +
  252 + // If mean PSF has not been found then create it now
  253 + if (ipsf == -1) {
  254 +
  255 + // Allocate new mean PSF
  256 + GLATMeanPsf* psf = new GLATMeanPsf(srcDir, static_cast<const GLATObservation&>(obs));
  257 +
  258 + // Set source name
  259 + std::string name = "SRC("+str(srcDir.ra_deg())+","+str(srcDir.dec_deg())+")";
  260 + psf->name(name);
  261 +
  262 + // Push mean PSF on stack
  263 + ((GLATResponse*)this)->m_ptsrc.push_back(psf);
  264 +
  265 + // Set index of mean PSF
  266 + ipsf = m_ptsrc.size()-1;
  267 +
  268 + // Debug option: dump mean PSF
  269 + #if G_DUMP_MEAN_PSF
  270 + std::cout << "Added new mean PSF \""+name+"\"" << std::endl;
  271 + std::cout << *psf << std::endl;
  272 + #endif
  273 +
  274 + } // endif: created new mean PSF
  275 +
  276 + // Get IRF value
  277 + double offset = dir->dist_deg(srcDir);
  278 + double irf = (*m_ptsrc[ipsf])(offset, srcEng.log10MeV());
  279 +
  280 + // Return IRF value
  281 + return irf;
  282 +}
  283 +
  284 +
  285 +/***********************************************************************//**
211 * @brief Return value of model instrument response function 286 * @brief Return value of model instrument response function
212 * 287 *
213 * @param[in] event Event. 288 * @param[in] event Event.
@@ -265,7 +340,7 @@ double GLATResponse::irf(const GLATEventAtom&amp; event, const GModelSky&amp; model, @@ -265,7 +340,7 @@ double GLATResponse::irf(const GLATEventAtom&amp; event, const GModelSky&amp; model,
265 340
266 /***********************************************************************//** 341 /***********************************************************************//**
267 * @brief Return value of model IRF for event bin 342 * @brief Return value of model IRF for event bin
268 - * 343 + *
269 * @param[in] event Event bin. 344 * @param[in] event Event bin.
270 * @param[in] model Source model. 345 * @param[in] model Source model.
271 * @param[in] srcEng True energy of photon. 346 * @param[in] srcEng True energy of photon.
@@ -404,14 +479,14 @@ double GLATResponse::irf(const GLATEventBin&amp; event, const GModelSky&amp; model, @@ -404,14 +479,14 @@ double GLATResponse::irf(const GLATEventBin&amp; event, const GModelSky&amp; model,
404 /***********************************************************************//** 479 /***********************************************************************//**
405 * @brief Return integral of instrument response function. 480 * @brief Return integral of instrument response function.
406 * 481 *
407 - * @param[in] model Source model. 482 + * @param[in] srcDir True photon arrival direction.
408 * @param[in] srcEng True energy of photon. 483 * @param[in] srcEng True energy of photon.
409 * @param[in] srcTime True photon arrival time. 484 * @param[in] srcTime True photon arrival time.
410 * @param[in] obs Observation. 485 * @param[in] obs Observation.
411 * 486 *
412 * @todo Not yet implemented. 487 * @todo Not yet implemented.
413 ***************************************************************************/ 488 ***************************************************************************/
414 -double GLATResponse::npred(const GModelSky& model, const GEnergy& srcEng, 489 +double GLATResponse::npred(const GSkyDir& srcDir, const GEnergy& srcEng,
415 const GTime& srcTime, 490 const GTime& srcTime,
416 const GObservation& obs) const 491 const GObservation& obs) const
417 { 492 {
inst/mwl/include/GMWLResponse.hpp
  View file @9a6488b
@@ -12,7 +12,7 @@ @@ -12,7 +12,7 @@
12 ***************************************************************************/ 12 ***************************************************************************/
13 /** 13 /**
14 * @file GMWLResponse.hpp 14 * @file GMWLResponse.hpp
15 - * @brief GMWLResponse class interface definition. 15 + * @brief Multi-wavelength response class interface definition
16 * @author J. Knodlseder 16 * @author J. Knodlseder
17 */ 17 */
18 18
@@ -26,11 +26,14 @@ @@ -26,11 +26,14 @@
26 /***********************************************************************//** 26 /***********************************************************************//**
27 * @class GMWLResponse 27 * @class GMWLResponse
28 * 28 *
29 - * @brief Interface for the Multi-wavelength response class 29 + * @brief Multi-wavelength response class
30 * 30 *
31 - * The Multi-wavelength response class is needed to evaluate the response  
32 - * to a source model. Since the multi-wavelength instrument classes work  
33 - * directly in photon space, the response is by definition unity. 31 + * This class implements a dummy response class for multi-wavelength
  32 + * observations. Since the multi-wavelength instrument classes handles data
  33 + * that a provided in photon space, no instrument response is in fact needed.
  34 + * The dummy response implemented by this class provides a simple diagonal
  35 + * response matrix that allows integration of multi-wavelength observations
  36 + * using the standard instrument specific interface.
34 ***************************************************************************/ 37 ***************************************************************************/
35 class GMWLResponse : public GResponse { 38 class GMWLResponse : public GResponse {
36 39
@@ -41,22 +44,25 @@ public: @@ -41,22 +44,25 @@ public:
41 virtual ~GMWLResponse(void); 44 virtual ~GMWLResponse(void);
42 45
43 // Operators 46 // Operators
44 - GMWLResponse& operator= (const GMWLResponse & rsp); 47 + virtual GMWLResponse& operator= (const GMWLResponse & rsp);
45 48
46 - // Reponse function computation methods  
47 - double irf(const GEvent& event, const GModelSky& model,  
48 - const GEnergy& srcEng, const GTime& srcTime,  
49 - const GObservation& obs) const { return 1.0; }  
50 - double npred(const GModelSky& model, const GEnergy& srcEng,  
51 - const GTime& srcTime,  
52 - const GObservation& obs) const { return 1.0; }  
53 -  
54 - // Pure virtual base class methods  
55 - void clear(void);  
56 - GMWLResponse* clone(void) const;  
57 - bool hasedisp(void) const { return false; }  
58 - bool hastdisp(void) const { return false; }  
59 - std::string print(void) const; 49 + // Implemented pure virtual methods
  50 + virtual void clear(void);
  51 + virtual GMWLResponse* clone(void) const;
  52 + virtual bool hasedisp(void) const { return false; }
  53 + virtual bool hastdisp(void) const { return false; }
  54 + virtual double irf(const GInstDir& obsDir,
  55 + const GEnergy& obsEng,
  56 + const GTime& obsTime,
  57 + const GSkyDir& srcDir,
  58 + const GEnergy& srcEng,
  59 + const GTime& srcTime,
  60 + const GObservation& obs) const { return 1.0; }
  61 + virtual double npred(const GSkyDir& srcDir,
  62 + const GEnergy& srcEng,
  63 + const GTime& srcTime,
  64 + const GObservation& obs) const { return 1.0; }
  65 + virtual std::string print(void) const;
60 66
61 protected: 67 protected:
62 // Protected methods 68 // Protected methods
pyext/GResponse.i
  View file @9a6488b
@@ -12,7 +12,7 @@ @@ -12,7 +12,7 @@
12 ***************************************************************************/ 12 ***************************************************************************/
13 /** 13 /**
14 * @file GResponse.i 14 * @file GResponse.i
15 - * @brief Abstract response base class python bindings 15 + * @brief Abstract response base class Python interface definition
16 * @author J. Knodlseder 16 * @author J. Knodlseder
17 */ 17 */
18 %{ 18 %{
@@ -25,9 +25,10 @@ @@ -25,9 +25,10 @@
25 /***********************************************************************//** 25 /***********************************************************************//**
26 * @class GResponse 26 * @class GResponse
27 * 27 *
28 - * @brief Abstract response base class python bindings 28 + * @brief Abstract response base class
29 ***************************************************************************/ 29 ***************************************************************************/
30 class GResponse { 30 class GResponse {
  31 +
31 public: 32 public:
32 // Constructors and destructors 33 // Constructors and destructors
33 GResponse(void); 34 GResponse(void);
@@ -35,16 +36,65 @@ public: @@ -35,16 +36,65 @@ public:
35 virtual ~GResponse(void); 36 virtual ~GResponse(void);
36 37
37 // Pure virtual methods 38 // Pure virtual methods
38 - virtual void clear(void) = 0;  
39 - virtual GResponse* clone(void) const = 0;  
40 - virtual bool hasedisp(void) const = 0;  
41 - virtual bool hastdisp(void) const = 0;  
42 - virtual double irf(const GEvent& event, const GModelSky& model,  
43 - const GEnergy& srcEng, const GTime& srcTime,  
44 - const GObservation& obs) const = 0;  
45 - virtual double npred(const GModelSky& model, const GEnergy& srcEng,  
46 - const GTime& srcTime,  
47 - const GObservation& obs) const = 0; 39 + virtual void clear(void) = 0;
  40 + virtual GResponse* clone(void) const = 0;
  41 + virtual bool hasedisp(void) const = 0;
  42 + virtual bool hastdisp(void) const = 0;
  43 + virtual double irf(const GInstDir& obsDir,
  44 + const GEnergy& obsEng,
  45 + const GTime& obsTime,
  46 + const GSkyDir& srcDir,
  47 + const GEnergy& srcEng,
  48 + const GTime& srcTime,
  49 + const GObservation& obs) const = 0;
  50 + virtual double npred(const GSkyDir& srcDir,
  51 + const GEnergy& srcEng,
  52 + const GTime& srcTime,
  53 + const GObservation& obs) const = 0;
  54 +
  55 + // Virtual methods
  56 + virtual double irf(const GEvent& event,
  57 + const GModelSky& model,
  58 + const GEnergy& srcEng,
  59 + const GTime& srcTime,
  60 + const GObservation& obs) const;
  61 + virtual double irf_ptsrc(const GInstDir& obsDir,
  62 + const GEnergy& obsEng,
  63 + const GTime& obsTime,
  64 + const GModelPointSource& model,
  65 + const GEnergy& srcEng,
  66 + const GTime& srcTime,
  67 + const GObservation& obs) const;
  68 + virtual double irf_extended(const GInstDir& obsDir,
  69 + const GEnergy& obsEng,
  70 + const GTime& obsTime,
  71 + const GModelExtendedSource& model,
  72 + const GEnergy& srcEng,
  73 + const GTime& srcTime,
  74 + const GObservation& obs) const;
  75 + virtual double irf_diffuse(const GInstDir& obsDir,
  76 + const GEnergy& obsEng,
  77 + const GTime& obsTime,
  78 + const GModelDiffuseSource& model,
  79 + const GEnergy& srcEng,
  80 + const GTime& srcTime,
  81 + const GObservation& obs) const;
  82 + virtual double npred(const GModelSky& model,
  83 + const GEnergy& srcEng,
  84 + const GTime& srcTime,
  85 + const GObservation& obs) const;
  86 + virtual double npred_ptsrc(const GModelPointSource& model,
  87 + const GEnergy& srcEng,
  88 + const GTime& srcTime,
  89 + const GObservation& obs) const;
  90 + virtual double npred_extended(const GModelExtendedSource& model,
  91 + const GEnergy& srcEng,
  92 + const GTime& srcTime,
  93 + const GObservation& obs) const;
  94 + virtual double npred_diffuse(const GModelDiffuseSource& model,
  95 + const GEnergy& srcEng,
  96 + const GTime& srcTime,
  97 + const GObservation& obs) const;
48 }; 98 };
49 99
50 100
src/obs/GResponse.cpp
  View file @9a6488b
@@ -18,16 +18,29 @@ @@ -18,16 +18,29 @@
18 #include <string> 18 #include <string>
19 #include <unistd.h> // access() function 19 #include <unistd.h> // access() function
20 #include "GResponse.hpp" 20 #include "GResponse.hpp"
  21 +#include "GObservation.hpp"
  22 +#include "GIntegral.hpp"
  23 +#include "GVector.hpp"
  24 +#include "GSkyDir.hpp"
21 #include "GException.hpp" 25 #include "GException.hpp"
  26 +#include "GTools.hpp"
22 27
23 /* __ Method name definitions ____________________________________________ */ 28 /* __ Method name definitions ____________________________________________ */
  29 +#define G_IRF_EXTENDED "GResponse::irf_extended(GInstDir&, GEnergy&,"\
  30 + " GTime&, GModelExtendedSource&, GEnergy&, GTime&, GObservation&)"
  31 +#define G_IRF_DIFFUSE "GResponse::irf_diffuse(GInstDir&, GEnergy&,"\
  32 + " GTime&, GModelDiffuseSource&, GEnergy&, GTime&, GObservation&)"
  33 +#define G_NPRED_EXTENDED "GResponse::npred_extended(GModelExtendedSource&,"\
  34 + " GEnergy&, GTime&, GObservation&)"
  35 +#define G_NPRED_DIFFUSE "GResponse::npred_diffuse(GModelDiffuseSource&,"\
  36 + " GEnergy&, GTime&, GObservation&)"
24 37
25 /* __ Macros _____________________________________________________________ */ 38 /* __ Macros _____________________________________________________________ */
26 39
27 /* __ Coding definitions _________________________________________________ */ 40 /* __ Coding definitions _________________________________________________ */
28 41
29 /* __ Debug definitions __________________________________________________ */ 42 /* __ Debug definitions __________________________________________________ */
30 - 43 +#define G_DEBUG_NPRED_EXTENDED 0 //!< Debug npred_extended
31 44
32 45
33 /*========================================================================== 46 /*==========================================================================
@@ -118,10 +131,286 @@ GResponse&amp; GResponse::operator= (const GResponse&amp; rsp) @@ -118,10 +131,286 @@ GResponse&amp; GResponse::operator= (const GResponse&amp; rsp)
118 = = 131 = =
119 ==========================================================================*/ 132 ==========================================================================*/
120 133
  134 +/***********************************************************************//**
  135 + * @brief Return value of instrument response function
  136 + *
  137 + * @param[in] event Event.
  138 + * @param[in] model Source model.
  139 + * @param[in] srcEng True energy of photon.
  140 + * @param[in] srcTime True photon arrival time.
  141 + * @param[in] obs Observation.
  142 + *
  143 + * @todo Throw exception if model pointer is invalid.
  144 + ***************************************************************************/
  145 +double GResponse::irf(const GEvent& event,
  146 + const GModelSky& model,
  147 + const GEnergy& srcEng,
  148 + const GTime& srcTime,
  149 + const GObservation& obs) const
  150 +{
  151 + // Initialise IRF value
  152 + double irf = 0.0;
  153 +
  154 + // Get model pointers
  155 + const GModelPointSource* ptsrc = dynamic_cast<const GModelPointSource*>(&model);
  156 + const GModelExtendedSource* extsrc = dynamic_cast<const GModelExtendedSource*>(&model);
  157 + const GModelDiffuseSource* difsrc = dynamic_cast<const GModelDiffuseSource*>(&model);
  158 +
  159 + // Call model dependent method
  160 + if (ptsrc != NULL)
  161 + irf = irf_ptsrc(event.dir(), event.energy(), event.time(),
  162 + *ptsrc, srcEng, srcTime, obs);
  163 + else if (extsrc != NULL)
  164 + irf = irf_extended(event.dir(), event.energy(), event.time(),
  165 + *extsrc, srcEng, srcTime, obs);
  166 + else if (difsrc != NULL)
  167 + irf = irf_diffuse(event.dir(), event.energy(), event.time(),
  168 + *difsrc, srcEng, srcTime, obs);
  169 +
  170 + // Return IRF value
  171 + return irf;
  172 +}
  173 +
  174 +
  175 +/***********************************************************************//**
  176 + * @brief Return value of point source instrument response function
  177 + *
  178 + * @param[in] obsDir Observed photon direction.
  179 + * @param[in] obsEng Observed energy of photon.
  180 + * @param[in] obsTime Observed photon arrival time.
  181 + * @param[in] model Point source model.
  182 + * @param[in] srcEng True energy of photon.
  183 + * @param[in] srcTime True photon arrival time.
  184 + * @param[in] obs Observation.
  185 + ***************************************************************************/
  186 +double GResponse::irf_ptsrc(const GInstDir& obsDir,
  187 + const GEnergy& obsEng,
  188 + const GTime& obsTime,
  189 + const GModelPointSource& model,
  190 + const GEnergy& srcEng,
  191 + const GTime& srcTime,
  192 + const GObservation& obs) const
  193 +{
  194 + // Get point source location
  195 + GSkyDir srcDir = model.dir();
  196 +
  197 + // Compute IRF
  198 + double irf = this->irf(obsDir, obsEng, obsTime, srcDir, srcEng, srcTime, obs);
  199 +
  200 + // Return IRF
  201 + return irf;
  202 +}
  203 +
  204 +
  205 +/***********************************************************************//**
  206 + * @brief Return value of extended source instrument response function
  207 + *
  208 + * @param[in] obsDir Observed photon direction.
  209 + * @param[in] obsEng Observed energy of photon.
  210 + * @param[in] obsTime Observed photon arrival time.
  211 + * @param[in] model Exteded source model.
  212 + * @param[in] srcEng True energy of photon.
  213 + * @param[in] srcTime True photon arrival time.
  214 + * @param[in] obs Observation.
  215 + *
  216 + * @todo Generic method is not yet implemented.
  217 + ***************************************************************************/
  218 +double GResponse::irf_extended(const GInstDir& obsDir,
  219 + const GEnergy& obsEng,
  220 + const GTime& obsTime,
  221 + const GModelExtendedSource& model,
  222 + const GEnergy& srcEng,
  223 + const GTime& srcTime,
  224 + const GObservation& obs) const
  225 +{
  226 + // Feature not yet implemented
  227 + throw GException::feature_not_implemented(G_IRF_EXTENDED,
  228 + "Extended IRF not yet implemented.");
  229 +
  230 + // Return IRF
  231 + return 0.0;
  232 +}
  233 +
  234 +
  235 +/***********************************************************************//**
  236 + * @brief Return value of diffuse source instrument response function
  237 + *
  238 + * @param[in] obsDir Observed photon direction.
  239 + * @param[in] obsEng Observed energy of photon.
  240 + * @param[in] obsTime Observed photon arrival time.
  241 + * @param[in] model Diffuse source model.
  242 + * @param[in] srcEng True energy of photon.
  243 + * @param[in] srcTime True photon arrival time.
  244 + * @param[in] obs Observation.
  245 + *
  246 + * @todo Generic method is not yet implemented.
  247 + ***************************************************************************/
  248 +double GResponse::irf_diffuse(const GInstDir& obsDir,
  249 + const GEnergy& obsEng,
  250 + const GTime& obsTime,
  251 + const GModelDiffuseSource& model,
  252 + const GEnergy& srcEng,
  253 + const GTime& srcTime,
  254 + const GObservation& obs) const
  255 +{
  256 + // Feature not yet implemented
  257 + throw GException::feature_not_implemented(G_IRF_DIFFUSE,
  258 + "Diffuse IRF not yet implemented.");
  259 +
  260 + // Return IRF
  261 + return 0.0;
  262 +}
  263 +
  264 +
  265 +/***********************************************************************//**
  266 + * @brief Return spatial integral of instrument response function
  267 + *
  268 + * @param[in] model Source model.
  269 + * @param[in] srcEng True energy of photon.
  270 + * @param[in] srcTime True photon arrival time.
  271 + * @param[in] obs Observation.
  272 + *
  273 + * @todo Throw exception if model pointer is invalid.
  274 + ***************************************************************************/
  275 +double GResponse::npred(const GModelSky& model,
  276 + const GEnergy& srcEng,
  277 + const GTime& srcTime,
  278 + const GObservation& obs) const
  279 +{
  280 + // Initialise Npred value
  281 + double npred = 0.0;
  282 +
  283 + // Get model pointers
  284 + const GModelPointSource* ptsrc = dynamic_cast<const GModelPointSource*>(&model);
  285 + const GModelExtendedSource* extsrc = dynamic_cast<const GModelExtendedSource*>(&model);
  286 + const GModelDiffuseSource* difsrc = dynamic_cast<const GModelDiffuseSource*>(&model);
  287 +
  288 + // Call model dependent method
  289 + if (ptsrc != NULL)
  290 + npred = npred_ptsrc(*ptsrc, srcEng, srcTime, obs);
  291 + else if (extsrc != NULL)
  292 + npred = npred_extended(*extsrc, srcEng, srcTime, obs);
  293 + else if (difsrc != NULL)
  294 + npred = npred_diffuse(*difsrc, srcEng, srcTime, obs);
  295 +
  296 + // Return response value
  297 + return npred;
  298 +}
  299 +
  300 +
  301 +/***********************************************************************//**
  302 + * @brief Return spatial integral of point source model
  303 + *
  304 + * @param[in] model Point source model.
  305 + * @param[in] srcEng True energy of photon.
  306 + * @param[in] srcTime True photon arrival time.
  307 + * @param[in] obs Observation.
  308 + ***************************************************************************/
  309 +double GResponse::npred_ptsrc(const GModelPointSource& model,
  310 + const GEnergy& srcEng,
  311 + const GTime& srcTime,
  312 + const GObservation& obs) const
  313 +{
  314 + // Get point source location
  315 + GSkyDir srcDir = model.dir();
  316 +
  317 + // Compute Npred
  318 + double npred = this->npred(srcDir, srcEng, srcTime, obs);
  319 +
  320 + // Return Npred
  321 + return npred;
  322 +}
  323 +
  324 +
  325 +/***********************************************************************//**
  326 + * @brief Return spatial integral of extended source model
  327 + *
  328 + * @param[in] model Extended source model.
  329 + * @param[in] srcEng True energy of photon.
  330 + * @param[in] srcTime True photon arrival time.
  331 + * @param[in] obs Observation.
  332 + ***************************************************************************/
  333 +double GResponse::npred_extended(const GModelExtendedSource& model,
  334 + const GEnergy& srcEng,
  335 + const GTime& srcTime,
  336 + const GObservation& obs) const
  337 +{
  338 + // Initialise Npred value
  339 + double npred = 0.0;
  340 +
  341 + // Compute rotation matrix to convert from native coordinates given by
  342 + // (theta,phi) into celestial coordinates.
  343 + GMatrix ry;
  344 + GMatrix rz;
  345 + GMatrix rot;
  346 + ry.eulery(model.radial()->dec() - 90.0);
  347 + rz.eulerz(-model.radial()->ra());
  348 + rot = transpose(ry * rz);
  349 +
  350 + // Set offset angle integration range
  351 + double theta_min = 0.0;
  352 + double theta_max = model.radial()->theta_max();
  353 +
  354 + // Perform offset angle integration if interval is valid
  355 + if (theta_max > theta_min) {
  356 +
  357 + // Setup integration kernel
  358 + GResponse::npred_kern_theta integrand(this,
  359 + model.radial(),
  360 + &srcEng,
  361 + &srcTime,
  362 + &obs,
  363 + &rot);
  364 +
  365 + // Integrate over theta
  366 + GIntegral integral(&integrand);
  367 + npred = integral.romb(theta_min, theta_max);
  368 +
  369 + // Compile option: Show integration results
  370 + #if G_DEBUG_NPRED_EXTENDED
  371 + std::cout << "npred_extended:";
  372 + std::cout << " theta_min=" << theta_min;
  373 + std::cout << " theta_max=" << theta_max;
  374 + std::cout << " npred=" << npred << std::endl;
  375 + #endif
  376 +
  377 + } // endif: offset angle range was valid
  378 +
  379 + // Return Npred
  380 + return npred;
  381 +}
  382 +
  383 +
  384 +/***********************************************************************//**
  385 + * @brief Return spatial integral of diffuse source model
  386 + *
  387 + * @param[in] model Diffuse source model.
  388 + * @param[in] srcEng True energy of photon.
  389 + * @param[in] srcTime True photon arrival time.
  390 + * @param[in] obs Observation.
  391 + *
  392 + * @exception GException::feature_not_implemented
  393 + * Method not yet implemented.
  394 + *
  395 + * @todo Implement method.
  396 + ***************************************************************************/
  397 +double GResponse::npred_diffuse(const GModelDiffuseSource& model,
  398 + const GEnergy& srcEng,
  399 + const GTime& srcTime,
  400 + const GObservation& obs) const
  401 +{
  402 + // Feature not yet implemented
  403 + throw GException::feature_not_implemented(G_NPRED_DIFFUSE,
  404 + "Method for extended source not yet implemented.");
  405 +
  406 + // Return Npred
  407 + return 0.0;
  408 +}
  409 +
121 410
122 /*========================================================================== 411 /*==========================================================================
123 = = 412 = =
124 - = Private methods = 413 + = Protected methods =
125 = = 414 = =
126 ==========================================================================*/ 415 ==========================================================================*/
127 416
@@ -157,6 +446,67 @@ void GResponse::free_members(void) @@ -157,6 +446,67 @@ void GResponse::free_members(void)
157 } 446 }
158 447
159 448
  449 +/***********************************************************************//**
  450 + * @brief Kernel for offset angle Npred integration
  451 + *
  452 + * @param[in] theta Radial model offset angle (radians).
  453 + ***************************************************************************/
  454 +double GResponse::npred_kern_theta::eval(double theta)
  455 +{
  456 + // Initialise Npred value
  457 + double npred = 0.0;
  458 +
  459 + // Get radial model value
  460 + double model = m_radial->eval(theta);
  461 +
  462 + // Compute sine of offset angle
  463 + double sin_theta = std::sin(theta);
  464 +
  465 + // Setup phi integration kernel
  466 + GResponse::npred_kern_phi integrand(m_rsp,
  467 + m_srcEng,
  468 + m_srcTime,
  469 + m_obs,
  470 + m_rot,
  471 + theta,
  472 + sin_theta);
  473 +
  474 + // Integrate over phi
  475 + GIntegral integral(&integrand);
  476 + npred = integral.romb(0.0, twopi) * sin_theta * model;
  477 +
  478 + // Return Npred
  479 + return npred;
  480 +}
  481 +
  482 +
  483 +/***********************************************************************//**
  484 + * @brief Kernel for azimuth angle Npred integration
  485 + *
  486 + * @param[in] phi Azimuth angle (radians).
  487 + ***************************************************************************/
  488 +double GResponse::npred_kern_phi::eval(double phi)
  489 +{
  490 + // Compute sky direction vector in native coordinates
  491 + double cos_phi = std::cos(phi);
  492 + double sin_phi = std::sin(phi);
  493 + GVector native(-cos_phi*m_sin_theta, sin_phi*m_sin_theta, m_cos_theta);
  494 +
  495 + // Rotate from native into celestial system
  496 + GVector cel = *m_rot * native;
  497 +
  498 + // Set sky direction
  499 + GSkyDir srcDir;
  500 + srcDir.celvector(cel);
  501 +
  502 + // Compute Npred for this sky direction
  503 + double npred = m_rsp->npred(srcDir, *m_srcEng, *m_srcTime, *m_obs);
  504 +
  505 + // Return Npred
  506 + return npred;
  507 +}
  508 +
  509 +
160 /*========================================================================== 510 /*==========================================================================
161 = = 511 = =
162 = Friends = 512 = Friends =