Johan Bregeon / gammalib

Commit 5624e7c0481d59fec8345895df18c4124c7a2ebd

Authored by knodlseder
1 parent 4ce81325

remove non-generic response methods

Inline Side-by-side
Showing 11 changed files with 359 additions and 294 deletions
include/GResponse.hpp
  View file @5624e7c
... ... @@ -12,7 +12,7 @@
12 12 ***************************************************************************/
13 13 /**
14 14 * @file GResponse.hpp
15   - * @brief Abstract response base class definition.
  15 + * @brief Abstract response base class definition
16 16 * @author J. Knodlseder
17 17 */
18 18  
... ... @@ -24,12 +24,6 @@
24 24 #include <iostream>
25 25 #include "GLog.hpp"
26 26 #include "GEvent.hpp"
27   -#include "GPointing.hpp"
28   -#include "GInstDir.hpp"
29   -#include "GRoi.hpp"
30   -#include "GGti.hpp"
31   -#include "GEbounds.hpp"
32   -#include "GSkyDir.hpp"
33 27 #include "GEnergy.hpp"
34 28 #include "GTime.hpp"
35 29  
... ... @@ -43,10 +37,14 @@ class GObservation;
43 37 *
44 38 * @brief Abstract interface for the instrument response function
45 39 *
46   - * The response function provides conversion methods between physical
47   - * parameters (such as source position, flux, ...) and the measured
48   - * instrumental parameters (such as measured energy, photon interaction,
49   - * ...).
  40 + * The response function provides conversion between physical parameters
  41 + * (such as source position, flux, ...) and the measured instrumental
  42 + * parameters (such as measured energy, photon interaction, ...).
  43 + * For a given observation, the irf method returns the instrument response
  44 + * for a given event and source model as function of the true photon energy
  45 + * and the true photon arrival time.
  46 + * The npred method returns the integral of the instrument response function
  47 + * over the dataspace. This method is only required for unbinned analysis.
50 48 ***************************************************************************/
51 49 class GResponse {
52 50  
... ... @@ -66,35 +64,21 @@ public:
66 64 // Pure virtual methods
67 65 virtual void clear(void) = 0;
68 66 virtual GResponse* clone(void) const = 0;
69   - virtual void load(const std::string& irfname) = 0;
70 67 virtual bool hasedisp(void) const = 0;
71 68 virtual bool hastdisp(void) const = 0;
72   - virtual double irf(const GInstDir& obsDir, const GEnergy& obsEng,
73   - const GTime& obsTime,
74   - const GSkyDir& srcDir, const GEnergy& srcEng,
75   - const GTime& srcTime,
76   - const GObservation& obs) const = 0;
77 69 virtual double irf(const GEvent& event, const GModel& model,
78 70 const GEnergy& srcEng, const GTime& srcTime,
79 71 const GObservation& obs) const = 0;
80   - virtual double nirf(const GSkyDir& srcDir, const GEnergy& srcEng,
81   - const GTime& srcTime,
82   - const GObservation& obs) const = 0;
  72 + virtual double npred(const GModel& model, const GEnergy& srcEng,
  73 + const GTime& srcTime,
  74 + const GObservation& obs) const = 0;
83 75 virtual std::string print(void) const = 0;
84 76  
85   - // Other methods
86   - virtual void caldb(const std::string& caldb);
87   - virtual std::string caldb(void) const { return m_caldb; }
88   -
89 77 protected:
90 78 // Protected methods
91 79 void init_members(void);
92 80 void copy_members(const GResponse& rsp);
93 81 void free_members(void);
94   -
95   - // Protected data area
96   - std::string m_caldb; //!< Name of or path to the calibration database
97   - std::string m_rspname; //!< Name of the instrument response
98 82 };
99 83  
100 84 #endif /* GRESPONSE_HPP */
... ...
inst/cta/include/GCTAResponse.hpp
  View file @5624e7c
... ... @@ -58,26 +58,33 @@ public:
58 58 // Implement pure virtual base class methods
59 59 void clear(void);
60 60 GCTAResponse* clone(void) const;
61   - void load(const std::string& rspname);
62 61 bool hasedisp(void) const { return false; }
63 62 bool hastdisp(void) const { return false; }
64   - double irf(const GInstDir& obsDir, const GEnergy& obsEng, const GTime& obsTime,
65   - const GSkyDir& srcDir, const GEnergy& srcEng, const GTime& srcTime,
66   - const GObservation& obs) const;
67 63 double irf(const GEvent& event, const GModel& model,
68 64 const GEnergy& srcEng, const GTime& srcTime,
69 65 const GObservation& obs) const;
70   - double nirf(const GSkyDir& srcDir, const GEnergy& srcEng, const GTime& srcTime,
71   - const GObservation& obs) const;
  66 + double npred(const GModel& model, const GEnergy& srcEng,
  67 + const GTime& srcTime,
  68 + const GObservation& obs) const;
72 69 std::string print(void) const;
73 70  
74 71 // Other Methods
  72 + void caldb(const std::string& caldb);
  73 + std::string caldb(void) const { return m_caldb; }
  74 + void load(const std::string& rspname);
  75 +
  76 + // Other response methods
  77 + double irf(const GInstDir& obsDir, const GEnergy& obsEng, const GTime& obsTime,
  78 + const GSkyDir& srcDir, const GEnergy& srcEng, const GTime& srcTime,
  79 + const GObservation& obs) const;
75 80 double irf(const GCTAEventAtom& event, const GModel& model,
76 81 const GEnergy& srcEng, const GTime& srcTime,
77 82 const GObservation& obs) const;
78 83 double irf(const GCTAEventBin& event, const GModel& model,
79 84 const GEnergy& srcEng, const GTime& srcTime,
80 85 const GObservation& obs) const;
  86 + double npred(const GSkyDir& srcDir, const GEnergy& srcEng, const GTime& srcTime,
  87 + const GObservation& obs) const;
81 88 double live(const GSkyDir& srcDir, const GEnergy& srcEng, const GTime& srcTime,
82 89 const GPointing& pnt) const;
83 90 double aeff(const GSkyDir& srcDir, const GEnergy& srcEng, const GTime& srcTime,
... ... @@ -97,8 +104,6 @@ public:
97 104 const GPointing& pnt, const GEbounds& ebds) const;
98 105 double ntdisp(const GSkyDir& srcDir, const GEnergy& srcEng, const GTime& srcTime,
99 106 const GPointing& pnt, const GGti& gti) const;
100   -
101   - // Other Methods
102 107 double psf(const double& theta, const double& sigma) const;
103 108 double psf_sigma(const GEnergy& srcEng) const;
104 109 double npsf(const double& psf, const double& radroi, const double& sigma) const;
... ... @@ -135,11 +140,13 @@ private:
135 140 };
136 141  
137 142 // Private data members
138   - GNodeArray m_nodes; //!< log(E) nodes for interpolation
139   - std::vector<double> m_logE; //!< log(E) = log10(E/TeV) - bin centre
140   - std::vector<double> m_aeff; //!< Effective area in square metres after all cuts
141   - std::vector<double> m_r68; //!< 68% containment radius of PSF post cuts in degrees
142   - std::vector<double> m_r80; //!< 80% containment radius of PSF post cuts in degrees
  143 + std::string m_caldb; //!< Name of or path to the calibration database
  144 + std::string m_rspname; //!< Name of the instrument response
  145 + GNodeArray m_nodes; //!< log(E) nodes for interpolation
  146 + std::vector<double> m_logE; //!< log(E) = log10(E/TeV) - bin centre
  147 + std::vector<double> m_aeff; //!< Effective area in square metres after all cuts
  148 + std::vector<double> m_r68; //!< 68% containment radius of PSF post cuts in degrees
  149 + std::vector<double> m_r80; //!< 80% containment radius of PSF post cuts in degrees
143 150  
144 151 };
145 152  
... ...
inst/cta/pyext/GCTAResponse.i
  View file @5624e7c
... ... @@ -37,25 +37,32 @@ public:
37 37 // Implement pure virtual base class methods
38 38 void clear(void);
39 39 GCTAResponse* clone(void) const;
40   - void load(const std::string& rspname);
41 40 bool hasedisp(void) const { return false; }
42 41 bool hastdisp(void) const { return false; }
43   - double irf(const GInstDir& obsDir, const GEnergy& obsEng, const GTime& obsTime,
44   - const GSkyDir& srcDir, const GEnergy& srcEng, const GTime& srcTime,
45   - const GObservation& obs) const;
46 42 double irf(const GEvent& event, const GModel& model,
47 43 const GEnergy& srcEng, const GTime& srcTime,
48 44 const GObservation& obs) const;
49   - double nirf(const GSkyDir& srcDir, const GEnergy& srcEng, const GTime& srcTime,
50   - const GObservation& obs) const;
  45 + double npred(const GModel& model, const GEnergy& srcEng,
  46 + const GTime& srcTime,
  47 + const GObservation& obs) const;
51 48  
52 49 // Other Methods
  50 + void caldb(const std::string& caldb);
  51 + std::string caldb(void) const { return m_caldb; }
  52 + void load(const std::string& rspname);
  53 +
  54 + // Other response methods
  55 + double irf(const GInstDir& obsDir, const GEnergy& obsEng, const GTime& obsTime,
  56 + const GSkyDir& srcDir, const GEnergy& srcEng, const GTime& srcTime,
  57 + const GObservation& obs) const;
53 58 double irf(const GCTAEventAtom& event, const GModel& model,
54 59 const GEnergy& srcEng, const GTime& srcTime,
55 60 const GObservation& obs) const;
56 61 double irf(const GCTAEventBin& event, const GModel& model,
57 62 const GEnergy& srcEng, const GTime& srcTime,
58 63 const GObservation& obs) const;
  64 + double npred(const GSkyDir& srcDir, const GEnergy& srcEng, const GTime& srcTime,
  65 + const GObservation& obs) const;
59 66 double live(const GSkyDir& srcDir, const GEnergy& srcEng, const GTime& srcTime,
60 67 const GPointing& pnt) const;
61 68 double aeff(const GSkyDir& srcDir, const GEnergy& srcEng, const GTime& srcTime,
... ... @@ -75,8 +82,6 @@ public:
75 82 const GPointing& pnt, const GEbounds& ebds) const;
76 83 double ntdisp(const GSkyDir& srcDir, const GEnergy& srcEng, const GTime& srcTime,
77 84 const GPointing& pnt, const GGti& gti) const;
78   -
79   - // Other Methods
80 85 double psf(const double& theta, const double& sigma) const;
81 86 double psf_sigma(const GEnergy& srcEng) const;
82 87 double npsf(const double& psf, const double& radroi, const double& sigma) const;
... ...
inst/cta/src/GCTAResponse.cpp
  View file @5624e7c
... ... @@ -33,6 +33,13 @@
33 33 #include "GIntegrand.hpp"
34 34  
35 35 /* __ Method name definitions ____________________________________________ */
  36 +#define G_CALDB "GCTAResponse::caldb(std::string&)"
  37 +#define G_IRF_ATOM "GCTAResponse::irf(GCTAEventAtom&, GModel&, GEnergy&,"\
  38 + " GTime&, GObservation&)"
  39 +#define G_IRF_BIN "GCTAResponse::irf(GCTAEventBin&, GModel&, GEnergy&,"\
  40 + " GTime&, GObservation&)"
  41 +#define G_NPRED "GCTAResponse::npred(GModel&, GEnergy&, GTime&,"\
  42 + " GObservation&)"
36 43 #define G_READ "GCTAResponse::read_performance_table(std::string&)"
37 44  
38 45 /* __ Macros _____________________________________________________________ */
... ... @@ -51,7 +58,7 @@
51 58 ==========================================================================*/
52 59  
53 60 /***********************************************************************//**
54   - * @brief Constructor
  61 + * @brief Void constructor
55 62 ***************************************************************************/
56 63 GCTAResponse::GCTAResponse(void) : GResponse()
57 64 {
... ... @@ -66,7 +73,7 @@ GCTAResponse::GCTAResponse(void) : GResponse()
66 73 /***********************************************************************//**
67 74 * @brief Copy constructor
68 75 *
69   - * @param[in] rsp Response to be copied
  76 + * @param[in] rsp Response.
70 77 ***************************************************************************/
71 78 GCTAResponse::GCTAResponse(const GCTAResponse& rsp) : GResponse(rsp)
72 79 {
... ... @@ -103,7 +110,7 @@ GCTAResponse::~GCTAResponse(void)
103 110 /***********************************************************************//**
104 111 * @brief Assignment operator
105 112 *
106   - * @param[in] rsp Response to be assigned
  113 + * @param[in] rsp Response.
107 114 ***************************************************************************/
108 115 GCTAResponse& GCTAResponse::operator= (const GCTAResponse& rsp)
109 116 {
... ... @@ -136,36 +143,93 @@ GCTAResponse&amp; GCTAResponse::operator= (const GCTAResponse&amp; rsp)
136 143 ==========================================================================*/
137 144  
138 145 /***********************************************************************//**
139   - * @brief Return value of point source instrument response function.
  146 + * @brief Clear instance
  147 +***************************************************************************/
  148 +void GCTAResponse::clear(void)
  149 +{
  150 + // Free class members (base and derived classes, derived class first)
  151 + free_members();
  152 + this->GResponse::free_members();
  153 +
  154 + // Initialise members
  155 + this->GResponse::init_members();
  156 + init_members();
  157 +
  158 + // Return
  159 + return;
  160 +}
  161 +
  162 +
  163 +/***********************************************************************//**
  164 + * @brief Clone instance
  165 +***************************************************************************/
  166 +GCTAResponse* GCTAResponse::clone(void) const
  167 +{
  168 + return new GCTAResponse(*this);
  169 +}
  170 +
  171 +
  172 +/***********************************************************************//**
  173 + * @brief Set path to the calibration database
140 174 *
141   - * @param[in] obsDir Observed photon direction.
142   - * @param[in] obsEng Observed energy of photon.
143   - * @param[in] obsTime Observed photon arrival time.
144   - * @param[in] srcDir True photon direction.
145   - * @param[in] srcEng True energy of photon.
146   - * @param[in] srcTime True photon arrival time.
147   - * @param[in] obs Observations.
  175 + * @param[in] caldb Path to calibration database
  176 + *
  177 + * @exception GException::caldb_not_found
  178 + * Calibration database repository not found.
148 179 *
149   - * This method implements the default and complete instrument response
150   - * function (IRF).
  180 + * This default method simply checks if the calibration database directory
  181 + * exists. If the directory exists, the path will be stored. No checking is
  182 + * implemented that checks for the consistency of the calibration database.
  183 + *
  184 + * @todo Implement a GCalDB class that handles any calibration database
  185 + * issues. GCalDB may be an abstract class for which instrument
  186 + * specific methods are implement to handle any instrument specific
  187 + * IRF database issues.
151 188 ***************************************************************************/
152   -double GCTAResponse::irf(const GInstDir& obsDir, const GEnergy& obsEng,
153   - const GTime& obsTime,
154   - const GSkyDir& srcDir, const GEnergy& srcEng,
155   - const GTime& srcTime, const GObservation& obs) const
  189 +void GCTAResponse::caldb(const std::string& caldb)
156 190 {
157   - // Get pointing
158   - const GPointing *pnt = obs.pointing(srcTime);
  191 + // Check if calibration database directory is accessible
  192 + if (access(caldb.c_str(), R_OK) != 0)
  193 + throw GException::caldb_not_found(G_CALDB, caldb);
  194 +
  195 + // Store the path to the calibration database
  196 + m_caldb = caldb;
159 197  
160   - // Get point source IRF components
161   - double rsp = live(srcDir, srcEng, srcTime, *pnt);
162   - rsp *= aeff(srcDir, srcEng, srcTime, *pnt);
163   - rsp *= psf(obsDir, srcDir, srcEng, srcTime, *pnt);
164   - rsp *= edisp(obsEng, srcDir, srcEng, srcTime, *pnt);
165   - rsp *= tdisp(obsTime, srcDir, srcEng, srcTime, *pnt);
  198 + // Return
  199 + return;
  200 +}
166 201  
167   - // Return IRF value
168   - return rsp;
  202 +
  203 +/***********************************************************************//**
  204 + * @brief Load CTA response.
  205 + *
  206 + * @param[in] irfname Name of CTA response (without any file extension).
  207 + *
  208 + * The actually dummy version of the CTA response loads a CTA performance
  209 + * table given in ASCII format into memory.
  210 + ***************************************************************************/
  211 +void GCTAResponse::load(const std::string& irfname)
  212 +{
  213 + // Save calibration database name
  214 + std::string caldb = m_caldb;
  215 +
  216 + // Clear instance
  217 + clear();
  218 +
  219 + // Restore calibration database name
  220 + m_caldb = caldb;
  221 +
  222 + // Build filename
  223 + std::string filename = m_caldb + "/" + irfname + ".dat";
  224 +
  225 + // Read performance table
  226 + read_performance_table(filename);
  227 +
  228 + // Store response name
  229 + m_rspname = irfname;
  230 +
  231 + // Return
  232 + return;
169 233 }
170 234  
171 235  
... ... @@ -204,6 +268,9 @@ double GCTAResponse::irf(const GEvent&amp; event, const GModel&amp; model,
204 268 * @param[in] srcEng True energy of photon.
205 269 * @param[in] srcTime True photon arrival time.
206 270 * @param[in] obs Observation.
  271 + *
  272 + * @exception GException::feature_not_implemented
  273 + * Diffuse IRF not yet implemented.
207 274 ***************************************************************************/
208 275 double GCTAResponse::irf(const GCTAEventAtom& event, const GModel& model,
209 276 const GEnergy& srcEng, const GTime& srcTime,
... ... @@ -224,6 +291,15 @@ double GCTAResponse::irf(const GCTAEventAtom&amp; event, const GModel&amp; model,
224 291  
225 292 } // endif: model was point source
226 293  
  294 + // ... otherwise return diffuse IRF
  295 + else {
  296 +
  297 + // Feature not yet implemented
  298 + throw GException::feature_not_implemented(G_IRF_ATOM,
  299 + "Diffuse IRF not yet implemented.");
  300 +
  301 + } // endelse: model was not a point source
  302 +
227 303 // Return IRF value
228 304 return rsp;
229 305 }
... ... @@ -237,6 +313,9 @@ double GCTAResponse::irf(const GCTAEventAtom&amp; event, const GModel&amp; model,
237 313 * @param[in] srcEng True energy of photon.
238 314 * @param[in] srcTime True photon arrival time.
239 315 * @param[in] obs Observation.
  316 + *
  317 + * @exception GException::feature_not_implemented
  318 + * Diffuse IRF not yet implemented.
240 319 ***************************************************************************/
241 320 double GCTAResponse::irf(const GCTAEventBin& event, const GModel& model,
242 321 const GEnergy& srcEng, const GTime& srcTime,
... ... @@ -257,28 +336,106 @@ double GCTAResponse::irf(const GCTAEventBin&amp; event, const GModel&amp; model,
257 336  
258 337 } // endif: model was point source
259 338  
  339 + // ... otherwise return diffuse IRF
  340 + else {
  341 +
  342 + // Feature not yet implemented
  343 + throw GException::feature_not_implemented(G_IRF_BIN,
  344 + "Diffuse IRF not yet implemented.");
  345 +
  346 + } // endelse: model was not a point source
  347 +
260 348 // Return IRF value
261 349 return rsp;
262 350 }
263 351  
264 352  
265 353 /***********************************************************************//**
266   - * @brief Return integral of instrument response function.
  354 + * @brief Return value of point source instrument response function.
267 355 *
  356 + * @param[in] obsDir Observed photon direction.
  357 + * @param[in] obsEng Observed energy of photon.
  358 + * @param[in] obsTime Observed photon arrival time.
268 359 * @param[in] srcDir True photon direction.
269 360 * @param[in] srcEng True energy of photon.
270 361 * @param[in] srcTime True photon arrival time.
  362 + * @param[in] obs Observations.
  363 + ***************************************************************************/
  364 +double GCTAResponse::irf(const GInstDir& obsDir, const GEnergy& obsEng,
  365 + const GTime& obsTime,
  366 + const GSkyDir& srcDir, const GEnergy& srcEng,
  367 + const GTime& srcTime, const GObservation& obs) const
  368 +{
  369 + // Get pointing
  370 + const GPointing *pnt = obs.pointing(srcTime);
  371 +
  372 + // Get point source IRF components
  373 + double rsp = live(srcDir, srcEng, srcTime, *pnt);
  374 + rsp *= aeff(srcDir, srcEng, srcTime, *pnt);
  375 + rsp *= psf(obsDir, srcDir, srcEng, srcTime, *pnt);
  376 + rsp *= edisp(obsEng, srcDir, srcEng, srcTime, *pnt);
  377 + rsp *= tdisp(obsTime, srcDir, srcEng, srcTime, *pnt);
  378 +
  379 + // Return IRF value
  380 + return rsp;
  381 +}
  382 +
  383 +
  384 +/***********************************************************************//**
  385 + * @brief Return integral of instrument response function.
  386 + *
  387 + * @param[in] model Source model.
  388 + * @param[in] srcEng True energy of photon.
  389 + * @param[in] srcTime True photon arrival time.
271 390 * @param[in] obs Observation.
272 391 * @param[in] roi Region of interest of data selection.
273 392 * @param[in] ebds Energy boundaries of data selection.
274 393 * @param[in] gti Good Time Intervals of data selection.
  394 + ***************************************************************************/
  395 +double GCTAResponse::npred(const GModel& model, const GEnergy& srcEng,
  396 + const GTime& srcTime, const GObservation& obs) const
  397 +{
  398 + // Initialise response value
  399 + double rsp = 0.0;
  400 +
  401 + // If model is a point source then return the point source IRF
  402 + if (model.spatial()->isptsource()) {
  403 +
  404 + // Get point source location
  405 + GSkyDir srcDir = static_cast<GModelSpatialPtsrc*>(model.spatial())->dir();
  406 +
  407 + // Compute IRF
  408 + rsp = npred(srcDir, srcEng, srcTime, obs);
  409 +
  410 + } // endif: model was point source
  411 +
  412 + // ... otherwise return diffuse IRF
  413 + else {
  414 +
  415 + // Feature not yet implemented
  416 + throw GException::feature_not_implemented(G_NPRED,
  417 + "Diffuse IRF not yet implemented.");
  418 +
  419 + } // endelse: model was not a point source
  420 +
  421 + // Return response value
  422 + return rsp;
  423 +}
  424 +
  425 +
  426 +/***********************************************************************//**
  427 + * @brief Return integral of instrument response function.
275 428 *
276   - * This method implements the default and complete integral of the instrument
277   - * response function (IRF). It may be overwritted by a specific method in the
278   - * derived class that drops response terms that are not used.
  429 + * @param[in] srcDir True photon direction.
  430 + * @param[in] srcEng True energy of photon.
  431 + * @param[in] srcTime True photon arrival time.
  432 + * @param[in] obs Observation.
  433 + * @param[in] roi Region of interest of data selection.
  434 + * @param[in] ebds Energy boundaries of data selection.
  435 + * @param[in] gti Good Time Intervals of data selection.
279 436 ***************************************************************************/
280   -double GCTAResponse::nirf(const GSkyDir& srcDir, const GEnergy& srcEng,
281   - const GTime& srcTime, const GObservation& obs) const
  437 +double GCTAResponse::npred(const GSkyDir& srcDir, const GEnergy& srcEng,
  438 + const GTime& srcTime, const GObservation& obs) const
282 439 {
283 440 // Get pointers
284 441 const GPointing *pnt = obs.pointing(srcTime);
... ... @@ -509,61 +666,6 @@ double GCTAResponse::ntdisp(const GSkyDir&amp; srcDir, const GEnergy&amp; srcEng,
509 666  
510 667  
511 668 /***********************************************************************//**
512   - * @brief Load CTA response.
513   - *
514   - * @param[in] irfname Name of CTA response (without any file extension).
515   - *
516   - * The actually dummy version of the CTA response loads a CTA performance
517   - * table given in ASCII format into memory.
518   - ***************************************************************************/
519   -void GCTAResponse::load(const std::string& irfname)
520   -{
521   - // Initialise response members
522   - free_members();
523   - init_members();
524   -
525   - // Build filename
526   - std::string filename = m_caldb + "/" + irfname + ".dat";
527   -
528   - // Read performance table
529   - read_performance_table(filename);
530   -
531   - // Store response name
532   - m_rspname = irfname;
533   -
534   - // Return
535   - return;
536   -}
537   -
538   -
539   -/***********************************************************************//**
540   - * @brief Clear instance
541   -***************************************************************************/
542   -void GCTAResponse::clear(void)
543   -{
544   - // Free class members (base and derived classes, derived class first)
545   - free_members();
546   - this->GResponse::free_members();
547   -
548   - // Initialise members
549   - this->GResponse::init_members();
550   - init_members();
551   -
552   - // Return
553   - return;
554   -}
555   -
556   -
557   -/***********************************************************************//**
558   - * @brief Clone instance
559   -***************************************************************************/
560   -GCTAResponse* GCTAResponse::clone(void) const
561   -{
562   - return new GCTAResponse(*this);
563   -}
564   -
565   -
566   -/***********************************************************************//**
567 669 * @brief Print CTA response information
568 670 ***************************************************************************/
569 671 std::string GCTAResponse::print(void) const
... ... @@ -786,6 +888,8 @@ double GCTAResponse::npsf_kern_azsym(const double&amp; rad,
786 888 void GCTAResponse::init_members(void)
787 889 {
788 890 // Initialise members
  891 + m_caldb.clear();
  892 + m_rspname.clear();
789 893 m_logE.clear();
790 894 m_aeff.clear();
791 895 m_r68.clear();
... ... @@ -804,11 +908,13 @@ void GCTAResponse::init_members(void)
804 908 void GCTAResponse::copy_members(const GCTAResponse& rsp)
805 909 {
806 910 // Copy attributes
807   - m_nodes = rsp.m_nodes;
808   - m_logE = rsp.m_logE;
809   - m_aeff = rsp.m_aeff;
810   - m_r68 = rsp.m_r68;
811   - m_r80 = rsp.m_r80;
  911 + m_caldb = rsp.m_caldb;
  912 + m_rspname = rsp.m_rspname;
  913 + m_nodes = rsp.m_nodes;
  914 + m_logE = rsp.m_logE;
  915 + m_aeff = rsp.m_aeff;
  916 + m_r68 = rsp.m_r68;
  917 + m_r80 = rsp.m_r80;
812 918  
813 919 // Return
814 920 return;
... ...
inst/lat/include/GLATResponse.hpp
  View file @5624e7c
1 1 /***************************************************************************
2   - * GLATResponse.hpp - GLAST LAT Response class *
  2 + * GLATResponse.hpp - Fermi LAT Response class *
3 3 * ----------------------------------------------------------------------- *
4 4 * copyright (C) 2008-2010 by Jurgen Knodlseder *
5 5 * ----------------------------------------------------------------------- *
... ... @@ -12,7 +12,7 @@
12 12 ***************************************************************************/
13 13 /**
14 14 * @file GLATResponse.hpp
15   - * @brief GLATResponse class interface definition.
  15 + * @brief Fermi LAT Response class definition
16 16 * @author J. Knodlseder
17 17 */
18 18  
... ... @@ -21,9 +21,9 @@
21 21  
22 22 /* __ Includes ___________________________________________________________ */
23 23 #include <vector>
  24 +#include <string>
24 25 #include "GLATEventAtom.hpp"
25 26 #include "GLATEventBin.hpp"
26   -#include "GLATResponseTable.hpp"
27 27 #include "GLATAeff.hpp"
28 28 #include "GLATPsf.hpp"
29 29 #include "GLATEdisp.hpp"
... ... @@ -32,17 +32,8 @@
32 32 #include "GModel.hpp"
33 33 #include "GObservation.hpp"
34 34 #include "GResponse.hpp"
35   -#include "GPointing.hpp"
36   -#include "GInstDir.hpp"
37   -#include "GRoi.hpp"
38   -#include "GGti.hpp"
39   -#include "GEbounds.hpp"
40   -#include "GSkyDir.hpp"
41 35 #include "GEnergy.hpp"
42 36 #include "GTime.hpp"
43   -#include "GVector.hpp"
44   -#include "GFits.hpp"
45   -#include "GFitsTable.hpp"
46 37  
47 38  
48 39 /***********************************************************************//**
... ... @@ -64,34 +55,33 @@ public:
64 55 // Implement pure virtual base class methods
65 56 void clear(void);
66 57 GLATResponse* clone(void) const;
67   - void load(const std::string& rspname);
68 58 bool hasedisp(void) const { return false; }
69 59 bool hastdisp(void) const { return false; }
70   - double irf(const GInstDir& obsDir, const GEnergy& obsEng,
71   - const GTime& obsTime,
72   - const GSkyDir& srcDir, const GEnergy& srcEng,
73   - const GTime& srcTime,
74   - const GObservation& obs) const;
75 60 double irf(const GEvent& event, const GModel& model,
76 61 const GEnergy& srcEng, const GTime& srcTime,
77 62 const GObservation& obs) const;
78   - double nirf(const GSkyDir& srcDir, const GEnergy& srcEng,
79   - const GTime& srcTime,
80   - const GObservation& obs) const;
  63 + double npred(const GModel& model, const GEnergy& srcEng,
  64 + const GTime& srcTime,
  65 + const GObservation& obs) const;
81 66 std::string print(void) const;
82 67  
83 68 // Other Methods
84   - int size(void) const { return m_aeff.size(); }
85   - GLATAeff* aeff(const int& index) const;
86   - GLATPsf* psf(const int& index) const;
87   - GLATEdisp* edisp(const int& index) const;
88   - void save(const std::string& rspname) const;
89   - double irf(const GLATEventAtom& event, const GModel& model,
90   - const GEnergy& srcEng, const GTime& srcTime,
91   - const GObservation& obs) const;
92   - double irf(const GLATEventBin& event, const GModel& model,
93   - const GEnergy& srcEng, const GTime& srcTime,
94   - const GObservation& obs) const;
  69 + void caldb(const std::string& caldb);
  70 + std::string caldb(void) const { return m_caldb; }
  71 + void load(const std::string& rspname);
  72 + int size(void) const { return m_aeff.size(); }
  73 + GLATAeff* aeff(const int& index) const;
  74 + GLATPsf* psf(const int& index) const;
  75 + GLATEdisp* edisp(const int& index) const;
  76 + void save(const std::string& rspname) const;
  77 +
  78 + // Reponse methods
  79 + double irf(const GLATEventAtom& event, const GModel& model,
  80 + const GEnergy& srcEng, const GTime& srcTime,
  81 + const GObservation& obs) const;
  82 + double irf(const GLATEventBin& event, const GModel& model,
  83 + const GEnergy& srcEng, const GTime& srcTime,
  84 + const GObservation& obs) const;
95 85  
96 86 private:
97 87 // Private methods
... ... @@ -100,12 +90,14 @@ private:
100 90 void free_members(void);
101 91  
102 92 // Private members
103   - bool m_hasfront; //!< Front IRF loaded
104   - bool m_hasback; //!< Back IRF loaded
105   - std::vector<GLATAeff*> m_aeff; //!< Effective areas
106   - std::vector<GLATPsf*> m_psf; //!< Point spread functions
107   - std::vector<GLATEdisp*> m_edisp; //!< Energy dispersions
108   - std::vector<GLATMeanPsf*> m_ptsrc; //!< Mean PSFs for point sources
  93 + std::string m_caldb; //!< Name of or path to the calibration database
  94 + std::string m_rspname; //!< Name of the instrument response
  95 + bool m_hasfront; //!< Front IRF loaded?
  96 + bool m_hasback; //!< Back IRF loaded?
  97 + std::vector<GLATAeff*> m_aeff; //!< Effective areas
  98 + std::vector<GLATPsf*> m_psf; //!< Point spread functions
  99 + std::vector<GLATEdisp*> m_edisp; //!< Energy dispersions
  100 + std::vector<GLATMeanPsf*> m_ptsrc; //!< Mean PSFs for point sources
109 101 };
110 102  
111 103 #endif /* GLATRESPONSE_HPP */
... ...
inst/lat/pyext/GLATResponse.i
  View file @5624e7c
... ... @@ -36,33 +36,32 @@ public:
36 36 // Implement pure virtual base class methods
37 37 void clear(void);
38 38 GLATResponse* clone(void) const;
39   - void load(const std::string& rspname);
40 39 bool hasedisp(void) const { return false; }
41 40 bool hastdisp(void) const { return false; }
42   - double irf(const GInstDir& obsDir, const GEnergy& obsEng,
43   - const GTime& obsTime,
44   - const GSkyDir& srcDir, const GEnergy& srcEng,
45   - const GTime& srcTime,
46   - const GObservation& obs) const;
47 41 double irf(const GEvent& event, const GModel& model,
48 42 const GEnergy& srcEng, const GTime& srcTime,
49 43 const GObservation& obs) const;
50   - double nirf(const GSkyDir& srcDir, const GEnergy& srcEng,
51   - const GTime& srcTime,
52   - const GObservation& obs) const;
  44 + double npred(const GModel& model, const GEnergy& srcEng,
  45 + const GTime& srcTime,
  46 + const GObservation& obs) const;
53 47  
54 48 // Other Methods
55   - int size(void) const { return m_aeff.size(); }
56   - GLATAeff* aeff(const int& index) const;
57   - GLATPsf* psf(const int& index) const;
58   - GLATEdisp* edisp(const int& index) const;
59   - void save(const std::string& rspname) const;
60   - double irf(const GLATEventAtom& event, const GModel& model,
61   - const GEnergy& srcEng, const GTime& srcTime,
62   - const GObservation& obs) const;
63   - double irf(const GLATEventBin& event, const GModel& model,
64   - const GEnergy& srcEng, const GTime& srcTime,
65   - const GObservation& obs) const;
  49 + void caldb(const std::string& caldb);
  50 + std::string caldb(void) const { return m_caldb; }
  51 + void load(const std::string& rspname);
  52 + int size(void) const { return m_aeff.size(); }
  53 + GLATAeff* aeff(const int& index) const;
  54 + GLATPsf* psf(const int& index) const;
  55 + GLATEdisp* edisp(const int& index) const;
  56 + void save(const std::string& rspname) const;
  57 +
  58 + // Reponse methods
  59 + double irf(const GLATEventAtom& event, const GModel& model,
  60 + const GEnergy& srcEng, const GTime& srcTime,
  61 + const GObservation& obs) const;
  62 + double irf(const GLATEventBin& event, const GModel& model,
  63 + const GEnergy& srcEng, const GTime& srcTime,
  64 + const GObservation& obs) const;
66 65 };
67 66  
68 67  
... ...
inst/lat/src/GLATResponse.cpp
  View file @5624e7c
... ... @@ -12,7 +12,7 @@
12 12 ***************************************************************************/
13 13 /**
14 14 * @file GLATResponse.cpp
15   - * @brief GLATResponse class implementation.
  15 + * @brief Fermi LAT Response class implementation.
16 16 * @author J. Knodlseder
17 17 */
18 18  
... ... @@ -27,19 +27,23 @@
27 27 #include "GLATEventBin.hpp"
28 28 #include "GLATEventCube.hpp"
29 29 #include "GLATException.hpp"
  30 +#include "GException.hpp"
30 31 #include "GModelSpatialPtsrc.hpp"
31   -#include "GVector.hpp"
32 32 #include "GFits.hpp"
33   -#include "GFitsTable.hpp"
34 33 #include "GTools.hpp"
35 34  
36 35 /* __ Method name definitions ____________________________________________ */
37   -#define G_LOAD "load(std::string&)"
  36 +#define G_CALDB "GLATResponse::caldb(std::string&)"
  37 +#define G_LOAD "GLATResponse::load(std::string&)"
38 38 #define G_AEFF "GLATResponse::aeff(int&)"
39 39 #define G_PSF "GLATResponse::psf(int&)"
40 40 #define G_EDISP "GLATResponse::edisp(int&)"
  41 +#define G_IRF_ATOM "GLATResponse::irf(GLATEventAtom&, GModel&, GEnergy&,"\
  42 + "GTime&, GObservation&)"
41 43 #define G_IRF_BIN "GLATResponse::irf(GLATEventBin&, GModel&, GEnergy&,"\
42 44 "GTime&, GObservation&)"
  45 +#define G_NPRED "GLATResponse::npred(GModel&, GEnergy&, GTime&,"\
  46 + " GObservation&)"
43 47  
44 48 /* __ Macros _____________________________________________________________ */
45 49  
... ... @@ -170,28 +174,33 @@ GLATResponse* GLATResponse::clone(void) const
170 174  
171 175  
172 176 /***********************************************************************//**
173   - * @brief Return value of point source instrument response function
  177 + * @brief Set path to the calibration database
174 178 *
175   - * @param[in] obsDir Observed photon direction.
176   - * @param[in] obsEng Observed energy of photon.
177   - * @param[in] obsTime Observed photon arrival time.
178   - * @param[in] srcDir True photon direction.
179   - * @param[in] srcEng True energy of photon.
180   - * @param[in] srcTime True photon arrival time.
181   - * @param[in] obs Observation.
  179 + * @param[in] caldb Path to calibration database
182 180 *
183   - * @todo Not yet implemented.
  181 + * @exception GException::caldb_not_found
  182 + * Calibration database repository not found.
  183 + *
  184 + * This default method simply checks if the calibration database directory
  185 + * exists. If the directory exists, the path will be stored. No checking is
  186 + * implemented that checks for the consistency of the calibration database.
  187 + *
  188 + * @todo Implement a GCalDB class that handles any calibration database
  189 + * issues. GCalDB may be an abstract class for which instrument
  190 + * specific methods are implement to handle any instrument specific
  191 + * IRF database issues.
184 192 ***************************************************************************/
185   -double GLATResponse::irf(const GInstDir& obsDir, const GEnergy& obsEng,
186   - const GTime& obsTime,
187   - const GSkyDir& srcDir, const GEnergy& srcEng,
188   - const GTime& srcTime, const GObservation& obs) const
  193 +void GLATResponse::caldb(const std::string& caldb)
189 194 {
190   - // Initialise
191   - double rsp = 0.0;
  195 + // Check if calibration database directory is accessible
  196 + if (access(caldb.c_str(), R_OK) != 0)
  197 + throw GException::caldb_not_found(G_CALDB, caldb);
192 198  
193   - // Return IRF value
194   - return rsp;
  199 + // Store the path to the calibration database
  200 + m_caldb = caldb;
  201 +
  202 + // Return
  203 + return;
195 204 }
196 205  
197 206  
... ... @@ -243,8 +252,8 @@ double GLATResponse::irf(const GLATEventAtom&amp; event, const GModel&amp; model,
243 252 // Initialise IRF with "no response"
244 253 double irf = 0.0;
245 254  
246   - // Determine index for diffuse response
247   - //int inx = event.diffinx(model.name());
  255 + // Dump warning that integration is not yet implemented
  256 + throw GException::feature_not_implemented(G_IRF_ATOM);
248 257  
249 258 // Return IRF value
250 259 return irf;
... ... @@ -385,22 +394,25 @@ double GLATResponse::irf(const GLATEventBin&amp; event, const GModel&amp; model,
385 394 /***********************************************************************//**
386 395 * @brief Return integral of instrument response function.
387 396 *
388   - * @param[in] srcDir True photon direction.
  397 + * @param[in] model Source model.
389 398 * @param[in] srcEng True energy of photon.
390 399 * @param[in] srcTime True photon arrival time.
391 400 * @param[in] obs Observation.
392 401 *
393 402 * @todo Not yet implemented.
394 403 ***************************************************************************/
395   -double GLATResponse::nirf(const GSkyDir& srcDir, const GEnergy& srcEng,
396   - const GTime& srcTime,
397   - const GObservation& obs) const
  404 +double GLATResponse::npred(const GModel& model, const GEnergy& srcEng,
  405 + const GTime& srcTime,
  406 + const GObservation& obs) const
398 407 {
399 408 // Initialise
400   - double nirf = 0.0;
  409 + double npred = 0.0;
  410 +
  411 + // Notify that method is not yet implemented
  412 + throw GException::feature_not_implemented(G_NPRED);
401 413  
402 414 // Return integrated IRF value
403   - return nirf;
  415 + return npred;
404 416 }
405 417  
406 418  
... ... @@ -422,9 +434,14 @@ double GLATResponse::nirf(const GSkyDir&amp; srcDir, const GEnergy&amp; srcEng,
422 434 ***************************************************************************/
423 435 void GLATResponse::load(const std::string& rspname)
424 436 {
425   - // Free members (cannot call clear() here since this will overwrite the
426   - // calibration database path m_caldb)
427   - free_members();
  437 + // Save calibration database name
  438 + std::string caldb = m_caldb;
  439 +
  440 + // Clear instance
  441 + clear();
  442 +
  443 + // Restore calibration database name
  444 + m_caldb = caldb;
428 445  
429 446 // Determine response types to be loaded
430 447 std::vector<std::string> array = split(rspname, "::");
... ... @@ -615,6 +632,8 @@ std::string GLATResponse::print(void) const
615 632 void GLATResponse::init_members(void)
616 633 {
617 634 // Initialise members
  635 + m_caldb.clear();
  636 + m_rspname.clear();
618 637 m_hasfront = false;
619 638 m_hasback = false;
620 639 m_aeff.clear();
... ... @@ -640,6 +659,8 @@ void GLATResponse::init_members(void)
640 659 void GLATResponse::copy_members(const GLATResponse& rsp)
641 660 {
642 661 // Copy members
  662 + m_caldb = rsp.m_caldb;
  663 + m_rspname = rsp.m_rspname;
643 664 m_hasfront = rsp.m_hasfront;
644 665 m_hasback = rsp.m_hasback;
645 666 m_aeff = rsp.m_aeff;
... ...
inst/mwl/include/GMWLResponse.hpp
  View file @5624e7c
... ... @@ -44,19 +44,16 @@ public:
44 44 GMWLResponse& operator= (const GMWLResponse & rsp);
45 45  
46 46 // Reponse function computation methods
47   - double irf(const GInstDir& obsDir, const GEnergy& obsEng, const GTime& obsTime,
48   - const GSkyDir& srcDir, const GEnergy& srcEng, const GTime& srcTime,
49   - const GObservation& obs) const { return 1.0; }
50 47 double irf(const GEvent& event, const GModel& model,
51 48 const GEnergy& srcEng, const GTime& srcTime,
52 49 const GObservation& obs) const { return 1.0; }
53   - double nirf(const GSkyDir& srcDir, const GEnergy& srcEng, const GTime& srcTime,
54   - const GObservation& obs) const { return 1.0; }
  50 + double npred(const GModel& model, const GEnergy& srcEng,
  51 + const GTime& srcTime,
  52 + const GObservation& obs) const { return 1.0; }
55 53  
56 54 // Pure virtual base class methods
57 55 void clear(void);
58 56 GMWLResponse* clone(void) const;
59   - void load(const std::string& irfname) { return; }
60 57 bool hasedisp(void) const { return false; }
61 58 bool hastdisp(void) const { return false; }
62 59 std::string print(void) const;
... ...
inst/mwl/src/GMWLResponse.cpp
  View file @5624e7c
... ... @@ -160,8 +160,8 @@ std::string GMWLResponse::print(void) const
160 160  
161 161 // Append header
162 162 result.append("=== GMWLResponse ===\n");
163   - result.append(parformat("Calibration database")+m_caldb+"\n");
164   - result.append(parformat("Response functions")+m_rspname);
  163 + //result.append(parformat("Calibration database")+m_caldb+"\n");
  164 + //result.append(parformat("Response functions")+m_rspname);
165 165  
166 166 // Return result
167 167 return result;
... ...
pyext/GResponse.i
  View file @5624e7c
1 1 /***************************************************************************
2   - * GResponse.i - Response abstract base class SWIG interface *
  2 + * GResponse.i - Abstract response base class *
3 3 * ----------------------------------------------------------------------- *
4 4 * copyright (C) 2008-2010 by Jurgen Knodlseder *
5 5 * ----------------------------------------------------------------------- *
... ... @@ -12,7 +12,7 @@
12 12 ***************************************************************************/
13 13 /**
14 14 * @file GResponse.i
15   - * @brief GResponse class SWIG file.
  15 + * @brief Abstract response base class python bindings
16 16 * @author J. Knodlseder
17 17 */
18 18 %{
... ... @@ -24,7 +24,7 @@
24 24 /***********************************************************************//**
25 25 * @class GResponse
26 26 *
27   - * @brief Abstract SWIG interface for the instrument response function.
  27 + * @brief Abstract response base class python bindings
28 28 ***************************************************************************/
29 29 class GResponse {
30 30 public:
... ... @@ -36,24 +36,14 @@ public:
36 36 // Pure virtual methods
37 37 virtual void clear(void) = 0;
38 38 virtual GResponse* clone(void) const = 0;
39   - virtual void load(const std::string& irfname) = 0;
40 39 virtual bool hasedisp(void) const = 0;
41 40 virtual bool hastdisp(void) const = 0;
42   - virtual double irf(const GInstDir& obsDir, const GEnergy& obsEng,
43   - const GTime& obsTime,
44   - const GSkyDir& srcDir, const GEnergy& srcEng,
45   - const GTime& srcTime,
46   - const GObservation& obs) const = 0;
47 41 virtual double irf(const GEvent& event, const GModel& model,
48 42 const GEnergy& srcEng, const GTime& srcTime,
49 43 const GObservation& obs) const = 0;
50   - virtual double nirf(const GSkyDir& srcDir, const GEnergy& srcEng,
51   - const GTime& srcTime,
52   - const GObservation& obs) const = 0;
53   -
54   - // Other methods
55   - virtual void caldb(const std::string& caldb);
56   - virtual std::string caldb(void) const { return m_caldb; }
  44 + virtual double npred(const GModel& model, const GEnergy& srcEng,
  45 + const GTime& srcTime,
  46 + const GObservation& obs) const = 0;
57 47 };
58 48  
59 49  
... ...
src/obs/GResponse.cpp
  View file @5624e7c
... ... @@ -21,7 +21,6 @@
21 21 #include "GException.hpp"
22 22  
23 23 /* __ Method name definitions ____________________________________________ */
24   -#define G_CALDB "GResponse::caldb(std::string&)"
25 24  
26 25 /* __ Macros _____________________________________________________________ */
27 26  
... ... @@ -53,7 +52,7 @@ GResponse::GResponse(void)
53 52 /***********************************************************************//**
54 53 * @brief Copy constructor
55 54 *
56   - * @param[in] rsp Response from which the instance should be built.
  55 + * @param[in] rsp Response.
57 56 ***************************************************************************/
58 57 GResponse::GResponse(const GResponse& rsp)
59 58 {
... ... @@ -90,7 +89,7 @@ GResponse::~GResponse(void)
90 89 /***********************************************************************//**
91 90 * @brief Assignment operator
92 91 *
93   - * @param[in] rsp Response which should be assigned.
  92 + * @param[in] rsp Response.
94 93 ***************************************************************************/
95 94 GResponse& GResponse::operator= (const GResponse& rsp)
96 95 {
... ... @@ -119,33 +118,6 @@ GResponse&amp; GResponse::operator= (const GResponse&amp; rsp)
119 118 = =
120 119 ==========================================================================*/
121 120  
122   -/***********************************************************************//**
123   - * @brief Set path to the calibration database
124   - *
125   - * @param[in] caldb Path to calibration database
126   - *
127   - * This default method simply checks if the calibration database directory
128   - * exists. If the directory exists, the path will be stored. No checking is
129   - * implemented that checks for the consistency of the calibration database.
130   - *
131   - * @todo Implement a GCalDB class that handles any calibration database
132   - * issues. GCalDB may be an abstract class for which instrument
133   - * specific methods are implement to handle any instrument specific
134   - * IRF database issues.
135   - ***************************************************************************/
136   -void GResponse::caldb(const std::string& caldb)
137   -{
138   - // Check if calibration database directory is accessible
139   - if (access(caldb.c_str(), R_OK) != 0)
140   - throw GException::caldb_not_found(G_CALDB, caldb);
141   -
142   - // Store the path to the calibration database
143   - m_caldb = caldb;
144   -
145   - // Return
146   - return;
147   -}
148   -
149 121  
150 122 /*==========================================================================
151 123 = =
... ... @@ -158,10 +130,6 @@ void GResponse::caldb(const std::string&amp; caldb)
158 130 ***************************************************************************/
159 131 void GResponse::init_members(void)
160 132 {
161   - // Initialise members
162   - m_caldb.clear();
163   - m_rspname.clear();
164   -
165 133 // Return
166 134 return;
167 135 }
... ... @@ -174,10 +142,6 @@ void GResponse::init_members(void)
174 142 ***************************************************************************/
175 143 void GResponse::copy_members(const GResponse& rsp)
176 144 {
177   - // Copy attributes
178   - m_caldb = rsp.m_caldb;
179   - m_rspname = rsp.m_rspname;
180   -
181 145 // Return
182 146 return;
183 147 }
... ...