Domenico Tiziani / gammalib

Commit 233655feb51273e5495f63ab6f862c125474e95c

Authored by Jürgen Knödlseder
1 parent 19469fd7

Make use of effective area classes in response. 

The CTA response class makes now use of the GCTAAeff classes to handle the effective area response. The GCTAResponse::load_aeff() method automatically detects the response type, allocates the appropriate GCTAAeff derived class, and loads the response information. The GCTAObservation::read() and GCTAObservation::write() method also handle now transparently all 3 effective area response formats.

Only the effective area correction for ARF response files is not yet implemented. This will be done once the response factorization has also be applied for the point spread function.
Inline Side-by-side
Showing 15 changed files with 400 additions and 359 deletions
inst/cta/include/GCTAAeff.hpp
  View file @233655f
... ... @@ -52,6 +52,8 @@ public:
52 52 virtual double operator()(const double& logE,
53 53 const double& theta = 0.0,
54 54 const double& phi = 0.0,
  55 + const double& zenith = 0.0,
  56 + const double& azimuth = 0.0,
55 57 const bool& etrue = true) const = 0;
56 58  
57 59 // Operators
... ... @@ -61,6 +63,7 @@ public:
61 63 virtual void clear(void) = 0;
62 64 virtual GCTAAeff* clone(void) const = 0;
63 65 virtual void load(const std::string& filename) = 0;
  66 + virtual std::string filename(void) const = 0;
64 67 virtual std::string print(void) const = 0;
65 68  
66 69 protected:
... ...
inst/cta/include/GCTAAeff2D.hpp
  View file @233655f
... ... @@ -56,12 +56,15 @@ public:
56 56 double operator()(const double& logE,
57 57 const double& theta = 0.0,
58 58 const double& phi = 0.0,
  59 + const double& zenith = 0.0,
  60 + const double& azimuth = 0.0,
59 61 const bool& etrue = true) const;
60 62  
61 63 // Implemented pure virtual methods
62 64 void clear(void);
63 65 GCTAAeff2D* clone(void) const;
64 66 void load(const std::string& filename);
  67 + std::string filename(void) const;
65 68 std::string print(void) const;
66 69  
67 70 // Methods
... ... @@ -74,7 +77,8 @@ private:
74 77 void free_members(void);
75 78  
76 79 // Members
77   - GCTAResponseTable m_aeff; //!< Aeff response table
  80 + std::string m_filename; //!< Name of Aeff response file
  81 + GCTAResponseTable m_aeff; //!< Aeff response table
78 82 };
79 83  
80 84 #endif /* GCTAAEFF2D_HPP */
... ...
inst/cta/include/GCTAAeffArf.hpp
  View file @233655f
... ... @@ -58,18 +58,21 @@ public:
58 58 double operator()(const double& logE,
59 59 const double& theta = 0.0,
60 60 const double& phi = 0.0,
  61 + const double& zenith = 0.0,
  62 + const double& azimuth = 0.0,
61 63 const bool& etrue = true) const;
62 64  
63 65 // Implemented pure virtual methods
64 66 void clear(void);
65 67 GCTAAeffArf* clone(void) const;
66 68 void load(const std::string& filename);
  69 + std::string filename(void) const;
67 70 std::string print(void) const;
68 71  
69 72 // Methods
70 73 int size(void) const { return m_logE.size(); }
71   - void offset_sigma(const double& sigma) { m_offset_sigma=sigma; }
72   - const double& offset_sigma(void) const { return m_offset_sigma; }
  74 + void sigma(const double& sigma) { m_sigma=sigma; }
  75 + const double& sigma(void) const { return m_sigma; }
73 76 void thetacut(const double& thetacut) { m_thetacut=thetacut; }
74 77 const double& thetacut(void) const { return m_thetacut; }
75 78 void scale(const double& scale) { m_scale=scale; }
... ... @@ -83,11 +86,12 @@ private:
83 86 void free_members(void);
84 87  
85 88 // Members
86   - GNodeArray m_logE; //!< log(E) nodes for Aeff interpolation
87   - std::vector<double> m_aeff; //!< Effective area in cm2
88   - double m_offset_sigma; //!< Sigma for offset angle computation (0=none)
89   - double m_thetacut; //!< Theta cut for ARF
90   - double m_scale; //!< Scale for ARF
  89 + std::string m_filename; //!< Name of Aeff response file
  90 + GNodeArray m_logE; //!< log(E) nodes for Aeff interpolation
  91 + std::vector<double> m_aeff; //!< Effective area in cm2
  92 + double m_sigma; //!< Sigma for offset angle computation (0=none)
  93 + double m_thetacut; //!< Theta cut for ARF
  94 + double m_scale; //!< Scale for ARF
91 95 };
92 96  
93 97 #endif /* GCTAAEFFARF_HPP */
... ...
inst/cta/include/GCTAAeffPerfTable.hpp
  View file @233655f
... ... @@ -58,18 +58,21 @@ public:
58 58 double operator()(const double& logE,
59 59 const double& theta = 0.0,
60 60 const double& phi = 0.0,
  61 + const double& zenith = 0.0,
  62 + const double& azimuth = 0.0,
61 63 const bool& etrue = true) const;
62 64  
63 65 // Implemented pure virtual methods
64 66 void clear(void);
65 67 GCTAAeffPerfTable* clone(void) const;
66 68 void load(const std::string& filename);
  69 + std::string filename(void) const;
67 70 std::string print(void) const;
68 71  
69 72 // Methods
70   - int size(void) const { return m_logE.size(); }
71   - void offset_sigma(const double& sigma) { m_offset_sigma=sigma; }
72   - const double& offset_sigma(void) const { return m_offset_sigma; }
  73 + int size(void) const { return m_logE.size(); }
  74 + void sigma(const double& sigma) { m_sigma=sigma; }
  75 + const double& sigma(void) const { return m_sigma; }
73 76  
74 77 private:
75 78 // Methods
... ... @@ -78,9 +81,10 @@ private:
78 81 void free_members(void);
79 82  
80 83 // Members
81   - GNodeArray m_logE; //!< log(E) nodes for Aeff interpolation
82   - std::vector<double> m_aeff; //!< Effective area in cm2
83   - double m_offset_sigma; //!< Sigma for offset angle computation (0=none)
  84 + std::string m_filename; //!< Name of Aeff response file
  85 + GNodeArray m_logE; //!< log(E) nodes for Aeff interpolation
  86 + std::vector<double> m_aeff; //!< Effective area in cm2
  87 + double m_sigma; //!< Sigma for offset angle computation (0=none)
84 88 };
85 89  
86 90 #endif /* GCTAAEFFPERFTABLE_HPP */
... ...
inst/cta/include/GCTAResponse.hpp
  View file @233655f
... ... @@ -55,6 +55,7 @@
55 55 #include "GCTAEventAtom.hpp"
56 56 #include "GCTADir.hpp"
57 57 #include "GCTAResponseTable.hpp"
  58 +#include "GCTAAeff.hpp"
58 59  
59 60 /* __ Type definitions ___________________________________________________ */
60 61 typedef std::vector<double> GCTAPsfPars;
... ... @@ -136,26 +137,25 @@ public:
136 137 const GObservation& obs) const;
137 138  
138 139 // Other Methods
139   - GCTAEventAtom* mc(const double& area, const GPhoton& photon,
140   - const GObservation& obs, GRan& ran) const;
141   - void caldb(const std::string& caldb);
142   - std::string caldb(void) const { return m_caldb; }
143   - void load(const std::string& rspname);
144   - void eps(const double& eps) { m_eps=eps; }
145   - const double& eps(void) const { return m_eps; }
146   - void offset_sigma(const double& sigma) { m_offset_sigma=sigma; }
147   - const double& offset_sigma(void) const { return m_offset_sigma; }
148   - std::string arffile(void) const { return m_arffile; }
149   - std::string rmffile(void) const { return m_rmffile; }
150   - std::string psffile(void) const { return m_psffile; }
151   - double arf_thetacut(void) const { return m_arf_thetacut; }
152   - double arf_scale(void) const { return m_arf_scale; }
153   - void arf_thetacut(const double& value) { m_arf_thetacut=value; }
154   - void arf_scale(const double& value) { m_arf_scale=value; }
155   - void load_arf(const std::string& filename);
156   - void load_psf(const std::string& filename);
157   - void read_arf(const GFitsTable* hdu);
158   - void read_psf(const GFitsTable* hdu);
  140 + GCTAEventAtom* mc(const double& area, const GPhoton& photon,
  141 + const GObservation& obs, GRan& ran) const;
  142 + void caldb(const std::string& caldb);
  143 + std::string caldb(void) const { return m_caldb; }
  144 + void load(const std::string& rspname);
  145 + void eps(const double& eps) { m_eps=eps; }
  146 + const double& eps(void) const { return m_eps; }
  147 + void offset_sigma(const double& sigma) { m_offset_sigma=sigma; }
  148 + const double& offset_sigma(void) const { return m_offset_sigma; }
  149 + std::string rmffile(void) const { return m_rmffile; }
  150 + std::string psffile(void) const { return m_psffile; }
  151 + void read_psf(const GFitsTable* hdu);
  152 +
  153 + // New factorisation
  154 + void load_aeff(const std::string& filename);
  155 + void load_psf(const std::string& filename);
  156 + void load_edisp(const std::string& filename) {}
  157 + const GCTAAeff* aeff(void) const { return m_aeff; }
  158 + void aeff(GCTAAeff* aeff) { m_aeff=aeff; }
159 159  
160 160 // Low-level response methods
161 161 double aeff(const double& theta,
... ... @@ -208,15 +208,10 @@ private:
208 208 // Private data members
209 209 std::string m_caldb; //!< Name of or path to the calibration database
210 210 std::string m_rspname; //!< Name of the instrument response
211   - std::string m_arffile; //!< Name of ARF file
212 211 std::string m_rmffile; //!< Name of RMF file
213 212 std::string m_psffile; //!< Name of PSF file
214   - double m_arf_thetacut; //!< Theta cut for ARF
215   - double m_arf_scale; //!< Scale for ARF
216 213 GNodeArray m_psf_logE; //!< log(E) nodes for PSF interpolation
217   - GNodeArray m_aeff_logE; //!< log(E) nodes for Aeff interpolation
218 214 std::vector<double> m_logE; //!< log(E) = log10(E/TeV) - bin centre
219   - std::vector<double> m_aeff; //!< Effective area in cm2 after all cuts
220 215 std::vector<double> m_r68; //!< 68% containment radius of PSF post cuts in degrees
221 216 std::vector<double> m_r80; //!< 80% containment radius of PSF post cuts in degrees
222 217 double m_eps; //!< Integration precision
... ... @@ -225,6 +220,9 @@ private:
225 220 // New PSF handling
226 221 int m_psf_version; //!< PSF version
227 222 GCTAResponseTable m_psf_table; //!< PSF response table
  223 +
  224 + // New PSF factorisation
  225 + GCTAAeff* m_aeff; //!< Effective area component
228 226 };
229 227  
230 228 #endif /* GCTARESPONSE_HPP */
... ...
inst/cta/pyext/GCTAAeff.i
  View file @233655f
... ... @@ -47,12 +47,15 @@ public:
47 47 virtual double operator()(const double& logE,
48 48 const double& theta = 0.0,
49 49 const double& phi = 0.0,
  50 + const double& zenith = 0.0,
  51 + const double& azimuth = 0.0,
50 52 const bool& etrue = true) const = 0;
51 53  
52 54 // Pure virtual methods
53   - virtual void clear(void) = 0;
54   - virtual GCTAAeff* clone(void) const = 0;
55   - virtual void load(const std::string& filename) = 0;
  55 + virtual void clear(void) = 0;
  56 + virtual GCTAAeff* clone(void) const = 0;
  57 + virtual void load(const std::string& filename) = 0;
  58 + virtual std::string filename(void) const = 0;
56 59 };
57 60  
58 61  
... ...
inst/cta/pyext/GCTAAeff2D.i
  View file @233655f
... ... @@ -49,12 +49,15 @@ public:
49 49 double operator()(const double& logE,
50 50 const double& theta = 0.0,
51 51 const double& phi = 0.0,
  52 + const double& zenith = 0.0,
  53 + const double& azimuth = 0.0,
52 54 const bool& etrue = true) const;
53 55  
54 56 // Implemented pure virtual methods
55 57 void clear(void);
56 58 GCTAAeff2D* clone(void) const;
57 59 void load(const std::string& filename);
  60 + std::string filename(void) const;
58 61  
59 62 // Methods
60 63 void read(const GFits* file);
... ...
inst/cta/pyext/GCTAAeffArf.i
  View file @233655f
... ... @@ -49,17 +49,20 @@ public:
49 49 double operator()(const double& logE,
50 50 const double& theta = 0.0,
51 51 const double& phi = 0.0,
  52 + const double& zenith = 0.0,
  53 + const double& azimuth = 0.0,
52 54 const bool& etrue = true) const;
53 55  
54 56 // Implemented pure virtual methods
55 57 void clear(void);
56 58 GCTAAeffArf* clone(void) const;
57 59 void load(const std::string& filename);
  60 + std::string filename(void) const;
58 61  
59 62 // Methods
60 63 int size(void) const;
61   - void offset_sigma(const double& sigma);
62   - const double& offset_sigma(void) const;
  64 + void sigma(const double& sigma);
  65 + const double& sigma(void) const;
63 66 void thetacut(const double& thetacut);
64 67 const double& thetacut(void) const;
65 68 void scale(const double& scale);
... ...
inst/cta/pyext/GCTAAeffPerfTable.i
  View file @233655f
... ... @@ -49,17 +49,20 @@ public:
49 49 double operator()(const double& logE,
50 50 const double& theta = 0.0,
51 51 const double& phi = 0.0,
  52 + const double& zenith = 0.0,
  53 + const double& azimuth = 0.0,
52 54 const bool& etrue = true) const;
53 55  
54 56 // Implemented pure virtual methods
55 57 void clear(void);
56 58 GCTAAeffPerfTable* clone(void) const;
57 59 void load(const std::string& filename);
  60 + std::string filename(void) const;
58 61  
59 62 // Methods
60 63 int size(void) const;
61   - void offset_sigma(const double& sigma);
62   - const double& offset_sigma(void) const;
  64 + void sigma(const double& sigma);
  65 + const double& sigma(void) const;
63 66 };
64 67  
65 68  
... ...
inst/cta/pyext/GCTAResponse.i
  View file @233655f
... ... @@ -85,26 +85,20 @@ public:
85 85 const GObservation& obs) const;
86 86  
87 87 // Other Methods
88   - GCTAEventAtom* mc(const double& area, const GPhoton& photon,
89   - const GObservation& obs, GRan& ran) const;
90   - void caldb(const std::string& caldb);
91   - std::string caldb(void) const;
92   - void load(const std::string& rspname);
93   - void eps(const double& eps);
94   - const double& eps(void) const;
95   - void offset_sigma(const double& sigma);
96   - const double& offset_sigma(void) const;
97   - std::string arffile(void) const;
98   - std::string rmffile(void) const;
99   - std::string psffile(void) const;
100   - double arf_thetacut(void) const;
101   - double arf_scale(void) const;
102   - void arf_thetacut(const double& value);
103   - void arf_scale(const double& value);
104   - void load_arf(const std::string& filename);
105   - void load_psf(const std::string& filename);
106   - void read_arf(const GFitsTable* hdu);
107   - void read_psf(const GFitsTable* hdu);
  88 + GCTAEventAtom* mc(const double& area, const GPhoton& photon,
  89 + const GObservation& obs, GRan& ran) const;
  90 + void caldb(const std::string& caldb);
  91 + std::string caldb(void) const;
  92 + void load(const std::string& rspname);
  93 + void eps(const double& eps);
  94 + const double& eps(void) const;
  95 + void offset_sigma(const double& sigma);
  96 + const double& offset_sigma(void) const;
  97 + void load_aeff(const std::string& filename);
  98 + void load_psf(const std::string& filename);
  99 + void load_edisp(const std::string& filename);
  100 + const GCTAAeff* aeff(void) const;
  101 + void aeff(GCTAAeff* aeff);
108 102  
109 103 // Low-level response methods
110 104 double aeff(const double& theta,
... ...
inst/cta/src/GCTAAeff2D.cpp
  View file @233655f
... ... @@ -28,7 +28,6 @@
28 28 #ifdef HAVE_CONFIG_H
29 29 #include <config.h>
30 30 #endif
31   -//#include <cmath>
32 31 #include "GTools.hpp"
33 32 #include "GFitsTable.hpp"
34 33 #include "GCTAAeff2D.hpp"
... ... @@ -155,8 +154,10 @@ GCTAAeff2D&amp; GCTAAeff2D::operator= (const GCTAAeff2D&amp; aeff)
155 154 * @brief Return effective area in units of cm2
156 155 *
157 156 * @param[in] logE Log10 of the true photon energy (TeV).
158   - * @param[in] theta Offset angle (rad). Defaults to 0.0.
159   - * @param[in] phi Azimuth angle (rad). Not used in this method.
  157 + * @param[in] theta Offset angle in camera system (rad). Defaults to 0.0.
  158 + * @param[in] phi Azimuth angle in camera system (rad). Not used in this method.
  159 + * @param[in] zenith Zenith angle in Earth system (rad). Not used in this method.
  160 + * @param[in] azimuth Azimuth angle in Earth system (rad). Not used in this method.
160 161 * @param[in] etrue Use true energy (true/false). Defaults to true.
161 162 *
162 163 * Returns the effective area in units of cm2 for a given energy and
... ... @@ -172,6 +173,8 @@ GCTAAeff2D&amp; GCTAAeff2D::operator= (const GCTAAeff2D&amp; aeff)
172 173 double GCTAAeff2D::operator()(const double& logE,
173 174 const double& theta,
174 175 const double& phi,
  176 + const double& zenith,
  177 + const double& azimuth,
175 178 const bool& etrue) const
176 179 {
177 180 // Set parameter index
... ... @@ -242,12 +245,30 @@ void GCTAAeff2D::load(const std::string&amp; filename)
242 245 // Read effective area from file
243 246 read(&fits);
244 247  
  248 + // Close FITS file
  249 + fits.close();
  250 +
  251 + // Store filename
  252 + m_filename = filename;
  253 +
245 254 // Return
246 255 return;
247 256 }
248 257  
249 258  
250 259 /***********************************************************************//**
  260 + * @brief Return filename
  261 + *
  262 + * @return Returns filename from which effective area was loaded
  263 + ***************************************************************************/
  264 +std::string GCTAAeff2D::filename(void) const
  265 +{
  266 + // Return filename
  267 + return m_filename;
  268 +}
  269 +
  270 +
  271 +/***********************************************************************//**
251 272 * @brief Read effective area from FITS file
252 273 *
253 274 * @param[in] fits FITS file pointer.
... ... @@ -258,8 +279,8 @@ void GCTAAeff2D::load(const std::string&amp; filename)
258 279 ***************************************************************************/
259 280 void GCTAAeff2D::read(const GFits* fits)
260 281 {
261   - // Clear instance
262   - clear();
  282 + // Clear response table
  283 + m_aeff.clear();
263 284  
264 285 // Get pointer to effective area HDU
265 286 GFitsTable* hdu = fits->table("EFFECTIVE AREA");
... ... @@ -302,6 +323,7 @@ std::string GCTAAeff2D::print(void) const
302 323  
303 324 // Append header
304 325 result.append("=== GCTAAeff2D ===");
  326 + result.append("\n"+parformat("Filename")+m_filename);
305 327 result.append("\n"+parformat("Number of energy bins")+str(m_aeff.axis(0)));
306 328 result.append("\n"+parformat("Number of offset bins")+str(m_aeff.axis(1)));
307 329 result.append("\n"+parformat("Log10(Energy) range"));
... ... @@ -326,6 +348,7 @@ std::string GCTAAeff2D::print(void) const
326 348 void GCTAAeff2D::init_members(void)
327 349 {
328 350 // Initialise members
  351 + m_filename.clear();
329 352 m_aeff.clear();
330 353  
331 354 // Return
... ... @@ -341,7 +364,8 @@ void GCTAAeff2D::init_members(void)
341 364 void GCTAAeff2D::copy_members(const GCTAAeff2D& aeff)
342 365 {
343 366 // Copy members
344   - m_aeff = aeff.m_aeff;
  367 + m_filename = aeff.m_filename;
  368 + m_aeff = aeff.m_aeff;
345 369  
346 370 // Return
347 371 return;
... ...
inst/cta/src/GCTAAeffArf.cpp
  View file @233655f
... ... @@ -157,8 +157,10 @@ GCTAAeffArf&amp; GCTAAeffArf::operator= (const GCTAAeffArf&amp; aeff)
157 157 * @brief Return effective area in units of cm2
158 158 *
159 159 * @param[in] logE Log10 of the true photon energy (TeV).
160   - * @param[in] theta Offset angle (rad). Defaults to 0.0.
161   - * @param[in] phi Azimuth angle (rad). Not used.
  160 + * @param[in] theta Offset angle in camera system (rad). Defaults to 0.0.
  161 + * @param[in] phi Azimuth angle in camera system (rad). Not used in this method.
  162 + * @param[in] zenith Zenith angle in Earth system (rad). Not used in this method.
  163 + * @param[in] azimuth Azimuth angle in Earth system (rad). Not used in this method.
162 164 * @param[in] etrue Use true energy (true/false). Not used.
163 165 *
164 166 * Returns the effective area in units of cm2 for a given energy and
... ... @@ -167,9 +169,11 @@ GCTAAeffArf&amp; GCTAAeffArf::operator= (const GCTAAeffArf&amp; aeff)
167 169 * becomes negative.
168 170 ***************************************************************************/
169 171 double GCTAAeffArf::operator()(const double& logE,
170   - const double& theta,
171   - const double& phi,
172   - const bool& etrue) const
  172 + const double& theta,
  173 + const double& phi,
  174 + const double& zenith,
  175 + const double& azimuth,
  176 + const bool& etrue) const
173 177 {
174 178 // Get effective area value in cm2
175 179 double aeff = m_logE.interpolate(logE, m_aeff);
... ... @@ -180,9 +184,9 @@ double GCTAAeffArf::operator()(const double&amp; logE,
180 184 }
181 185  
182 186 // Optionally add in Gaussian offset angle dependence
183   - if (m_offset_sigma != 0.0) {
  187 + if (m_sigma != 0.0) {
184 188 double offset = theta * rad2deg;
185   - double arg = offset * offset / m_offset_sigma;
  189 + double arg = offset * offset / m_sigma;
186 190 double scale = exp(-0.5 * arg * arg);
187 191 aeff *= scale;
188 192 }
... ... @@ -242,10 +246,6 @@ GCTAAeffArf* GCTAAeffArf::clone(void) const
242 246 ***************************************************************************/
243 247 void GCTAAeffArf::load(const std::string& filename)
244 248 {
245   - // Clear arrays
246   - m_logE.clear();
247   - m_aeff.clear();
248   -
249 249 // Open ARF FITS file
250 250 GFits file(filename);
251 251  
... ... @@ -258,12 +258,27 @@ void GCTAAeffArf::load(const std::string&amp; filename)
258 258 // Close ARF FITS file
259 259 file.close();
260 260  
  261 + // Store filename
  262 + m_filename = filename;
  263 +
261 264 // Return
262 265 return;
263 266 }
264 267  
265 268  
266 269 /***********************************************************************//**
  270 + * @brief Return filename
  271 + *
  272 + * @return Returns filename from which effective area was loaded
  273 + ***************************************************************************/
  274 +std::string GCTAAeffArf::filename(void) const
  275 +{
  276 + // Return filename
  277 + return m_filename;
  278 +}
  279 +
  280 +
  281 +/***********************************************************************//**
267 282 * @brief Read CTA ARF vector
268 283 *
269 284 * @param[in] hdu FITS table pointer.
... ... @@ -388,7 +403,7 @@ void GCTAAeffArf::read_arf(const GFitsTable* hdu)
388 403 } // endfor: looped over nodes
389 404  
390 405 // Disable offset angle dependence
391   - m_offset_sigma = 0.0;
  406 + m_sigma = 0.0;
392 407  
393 408 // Return
394 409 return;
... ... @@ -411,16 +426,17 @@ std::string GCTAAeffArf::print(void) const
411 426  
412 427 // Append information
413 428 result.append("=== GCTAAeffArf ===");
  429 + result.append("\n"+parformat("Filename")+m_filename);
414 430 result.append("\n"+parformat("Number of energy bins")+str(size()));
415 431 result.append("\n"+parformat("Log10(Energy) range"));
416 432 result.append(str(emin)+" - "+str(emax)+" TeV");
417 433  
418 434 // Append offset angle dependence
419   - if (m_offset_sigma == 0) {
  435 + if (m_sigma == 0) {
420 436 result.append("\n"+parformat("Offset angle dependence")+"none");
421 437 }
422 438 else {
423   - std::string txt = "Fixed sigma="+str(m_offset_sigma);
  439 + std::string txt = "Fixed sigma="+str(m_sigma);
424 440 result.append("\n"+parformat("Offset angle dependence")+txt);
425 441 }
426 442  
... ... @@ -441,11 +457,12 @@ std::string GCTAAeffArf::print(void) const
441 457 void GCTAAeffArf::init_members(void)
442 458 {
443 459 // Initialise members
  460 + m_filename.clear();
444 461 m_logE.clear();
445 462 m_aeff.clear();
446   - m_offset_sigma = 3.0;
447   - m_thetacut = 0.0;
448   - m_scale = 1.0;
  463 + m_sigma = 0.0;
  464 + m_thetacut = 0.0;
  465 + m_scale = 1.0;
449 466  
450 467 // Return
451 468 return;
... ... @@ -460,11 +477,12 @@ void GCTAAeffArf::init_members(void)
460 477 void GCTAAeffArf::copy_members(const GCTAAeffArf& aeff)
461 478 {
462 479 // Copy members
463   - m_logE = aeff.m_logE;
464   - m_aeff = aeff.m_aeff;
465   - m_offset_sigma = aeff.m_offset_sigma;
466   - m_thetacut = aeff.m_thetacut;
467   - m_scale = aeff.m_scale;
  480 + m_filename = aeff.m_filename;
  481 + m_logE = aeff.m_logE;
  482 + m_aeff = aeff.m_aeff;
  483 + m_sigma = aeff.m_sigma;
  484 + m_thetacut = aeff.m_thetacut;
  485 + m_scale = aeff.m_scale;
468 486  
469 487 // Return
470 488 return;
... ...
inst/cta/src/GCTAAeffPerfTable.cpp
  View file @233655f
... ... @@ -157,8 +157,10 @@ GCTAAeffPerfTable&amp; GCTAAeffPerfTable::operator= (const GCTAAeffPerfTable&amp; aeff)
157 157 * @brief Return effective area in units of cm2
158 158 *
159 159 * @param[in] logE Log10 of the true photon energy (TeV).
160   - * @param[in] theta Offset angle (rad). Defaults to 0.0.
161   - * @param[in] phi Azimuth angle (rad). Not used.
  160 + * @param[in] theta Offset angle in camera system (rad). Defaults to 0.0.
  161 + * @param[in] phi Azimuth angle in camera system (rad). Not used in this method.
  162 + * @param[in] zenith Zenith angle in Earth system (rad). Not used in this method.
  163 + * @param[in] azimuth Azimuth angle in Earth system (rad). Not used in this method.
162 164 * @param[in] etrue Use true energy (true/false). Not used.
163 165 *
164 166 * Returns the effective area in units of cm2 for a given energy and
... ... @@ -169,6 +171,8 @@ GCTAAeffPerfTable&amp; GCTAAeffPerfTable::operator= (const GCTAAeffPerfTable&amp; aeff)
169 171 double GCTAAeffPerfTable::operator()(const double& logE,
170 172 const double& theta,
171 173 const double& phi,
  174 + const double& zenith,
  175 + const double& azimuth,
172 176 const bool& etrue) const
173 177 {
174 178 // Get effective area value in cm2
... ... @@ -180,9 +184,9 @@ double GCTAAeffPerfTable::operator()(const double&amp; logE,
180 184 }
181 185  
182 186 // Optionally add in Gaussian offset angle dependence
183   - if (m_offset_sigma != 0.0) {
  187 + if (m_sigma != 0.0) {
184 188 double offset = theta * rad2deg;
185   - double arg = offset * offset / m_offset_sigma;
  189 + double arg = offset * offset / m_sigma;
186 190 double scale = exp(-0.5 * arg * arg);
187 191 aeff *= scale;
188 192 }
... ... @@ -283,12 +287,27 @@ void GCTAAeffPerfTable::load(const std::string&amp; filename)
283 287  
284 288 } // endwhile: looped over lines
285 289  
  290 + // Store filename
  291 + m_filename = filename;
  292 +
286 293 // Return
287 294 return;
288 295 }
289 296  
290 297  
291 298 /***********************************************************************//**
  299 + * @brief Return filename
  300 + *
  301 + * @return Returns filename from which effective area was loaded
  302 + ***************************************************************************/
  303 +std::string GCTAAeffPerfTable::filename(void) const
  304 +{
  305 + // Return filename
  306 + return m_filename;
  307 +}
  308 +
  309 +
  310 +/***********************************************************************//**
292 311 * @brief Print effective area information
293 312 *
294 313 * @return Content of effective area instance.
... ... @@ -304,16 +323,17 @@ std::string GCTAAeffPerfTable::print(void) const
304 323  
305 324 // Append information
306 325 result.append("=== GCTAAeffPerfTable ===");
  326 + result.append("\n"+parformat("Filename")+m_filename);
307 327 result.append("\n"+parformat("Number of energy bins")+str(size()));
308 328 result.append("\n"+parformat("Log10(Energy) range"));
309 329 result.append(str(emin)+" - "+str(emax)+" TeV");
310 330  
311 331 // Append offset angle dependence
312   - if (m_offset_sigma == 0) {
  332 + if (m_sigma == 0) {
313 333 result.append("\n"+parformat("Offset angle dependence")+"none");
314 334 }
315 335 else {
316   - std::string txt = "Fixed sigma="+str(m_offset_sigma);
  336 + std::string txt = "Fixed sigma="+str(m_sigma);
317 337 result.append("\n"+parformat("Offset angle dependence")+txt);
318 338 }
319 339  
... ... @@ -334,9 +354,10 @@ std::string GCTAAeffPerfTable::print(void) const
334 354 void GCTAAeffPerfTable::init_members(void)
335 355 {
336 356 // Initialise members
  357 + m_filename.clear();
337 358 m_logE.clear();
338 359 m_aeff.clear();
339   - m_offset_sigma = 3.0;
  360 + m_sigma = 3.0;
340 361  
341 362 // Return
342 363 return;
... ... @@ -351,9 +372,10 @@ void GCTAAeffPerfTable::init_members(void)
351 372 void GCTAAeffPerfTable::copy_members(const GCTAAeffPerfTable& aeff)
352 373 {
353 374 // Copy members
354   - m_logE = aeff.m_logE;
355   - m_aeff = aeff.m_aeff;
356   - m_offset_sigma = aeff.m_offset_sigma;
  375 + m_filename = aeff.m_filename;
  376 + m_logE = aeff.m_logE;
  377 + m_aeff = aeff.m_aeff;
  378 + m_sigma = aeff.m_sigma;
357 379  
358 380 // Return
359 381 return;
... ...
inst/cta/src/GCTAObservation.cpp
  View file @233655f
... ... @@ -29,17 +29,21 @@
29 29 #include <config.h>
30 30 #endif
31 31 #include "GObservationRegistry.hpp"
32   -#include "GCTAException.hpp"
33   -#include "GCTAObservation.hpp"
34   -#include "GCTAEventList.hpp"
35   -#include "GCTAEventCube.hpp"
36   -#include "GCTARoi.hpp"
37 32 #include "GException.hpp"
38 33 #include "GFits.hpp"
39 34 #include "GTools.hpp"
40 35 #include "GModelSpatialPtsrc.hpp"
41 36 #include "GIntegral.hpp"
42 37 #include "GIntegrand.hpp"
  38 +#include "GCTAException.hpp"
  39 +#include "GCTAObservation.hpp"
  40 +#include "GCTAEventList.hpp"
  41 +#include "GCTAEventCube.hpp"
  42 +#include "GCTARoi.hpp"
  43 +#include "GCTAAeff.hpp"
  44 +#include "GCTAAeff2D.hpp"
  45 +#include "GCTAAeffArf.hpp"
  46 +#include "GCTAAeffPerfTable.hpp"
43 47  
44 48 /* __ Globals ____________________________________________________________ */
45 49 const GCTAObservation g_obs_cta_seed("CTA");
... ... @@ -339,18 +343,18 @@ void GCTAObservation::pointing(const GCTAPointing&amp; pointing)
339 343 * for a binned observation.
340 344 *
341 345 * @todo Still supports old ARF, PSF and RMF parameter names.
  346 + * @todo PSF correction for ARF not yet implemented.
342 347 ***************************************************************************/
343 348 void GCTAObservation::read(const GXmlElement& xml)
344 349 {
345 350 // Clear observation
346 351 clear();
347 352  
348   - // Initialise response parameters
349   - std::string aeff;
350   - std::string psf;
351   - std::string rmf;
352   - double thetacut = 0.0;
353   - double scale = 1.0;
  353 + // Delete old response function
  354 + if (m_response != NULL) delete m_response;
  355 +
  356 + // Allocate new CTA response function
  357 + m_response = new GCTAResponse;
354 358  
355 359 // Extract instrument name
356 360 m_instrument = xml.attribute("instrument");
... ... @@ -401,20 +405,53 @@ void GCTAObservation::read(const GXmlElement&amp; xml)
401 405 else if ((par->attribute("name") == "EffectiveArea") ||
402 406 (par->attribute("name") == "ARF")) {
403 407  
404   - // Store filename
405   - aeff = par->attribute("file");
  408 + // Get filename
  409 + std::string filename = par->attribute("file");
406 410  
407   - // Optionally extract thetacut (0.0 if no thetacut)
408   - std::string s_thetacut = par->attribute("thetacut");
409   - if (s_thetacut.length() > 0) {
410   - thetacut = todouble(s_thetacut);
411   - }
412   -
413   - // Optionally extract scale factor (1.0 if no scale)
414   - std::string s_scale = par->attribute("scale");
415   - if (s_scale.length() > 0) {
416   - scale = todouble(s_scale);
417   - }
  411 + // If filename is not empty then load effective area
  412 + if (strip_whitespace(filename).length() > 0) {
  413 +
  414 + // Load effective area
  415 + m_response->load_aeff(filename);
  416 +
  417 + // Optional attributes
  418 + double thetacut = 0.0;
  419 + double scale = 1.0;
  420 + double sigma = 0.0;
  421 +
  422 + // Optionally extract thetacut (0.0 if no thetacut)
  423 + std::string s_thetacut = par->attribute("thetacut");
  424 + if (s_thetacut.length() > 0) {
  425 + thetacut = todouble(s_thetacut);
  426 + }
  427 +
  428 + // Optionally extract scale factor (1.0 if no scale)
  429 + std::string s_scale = par->attribute("scale");
  430 + if (s_scale.length() > 0) {
  431 + scale = todouble(s_scale);
  432 + }
  433 +
  434 + // Optionally extract sigma (0.0 if no sigma)
  435 + std::string s_sigma = par->attribute("sigma");
  436 + if (s_sigma.length() > 0) {
  437 + sigma = todouble(s_sigma);
  438 + }
  439 +
  440 + // If we have an ARF then set attributes
  441 + GCTAAeffArf* arf = const_cast<GCTAAeffArf*>(dynamic_cast<const GCTAAeffArf*>(m_response->aeff()));
  442 + if (arf != NULL) {
  443 + arf->thetacut(thetacut);
  444 + arf->scale(scale);
  445 + arf->sigma(sigma);
  446 + }
  447 +
  448 + // If we have a performance table then set attributes
  449 + GCTAAeffPerfTable* perf = const_cast<GCTAAeffPerfTable*>(dynamic_cast<const GCTAAeffPerfTable*>(m_response->aeff()));
  450 + if (perf != NULL) {
  451 + perf->sigma(sigma);
  452 + }
  453 +
  454 + } // endif: effective area filename was valid
418 455  
419 456 // Increase number of parameters
420 457 npar[1]++;
... ... @@ -423,7 +460,16 @@ void GCTAObservation::read(const GXmlElement&amp; xml)
423 460 // Handle PSF
424 461 else if ((par->attribute("name") == "PointSpreadFunction") ||
425 462 (par->attribute("name") == "PSF")) {
426   - psf = par->attribute("file");
  463 +
  464 + // Get filename
  465 + std::string filename = par->attribute("file");
  466 +
  467 + // If filename is not empty then load point spread function
  468 + if (strip_whitespace(filename).length() > 0) {
  469 + m_response->load_psf(filename);
  470 + }
  471 +
  472 + // Increase number of parameters
427 473 npar[2]++;
428 474 }
429 475  
... ... @@ -431,7 +477,16 @@ void GCTAObservation::read(const GXmlElement&amp; xml)
431 477 // Handle RMF
432 478 else if ((par->attribute("name") == "EnergyDispersion") ||
433 479 (par->attribute("name") == "RMF")) {
434   - rmf = par->attribute("file");
  480 +
  481 + // Get filename
  482 + std::string filename = par->attribute("file");
  483 +
  484 + // If filename is not empty then load energy dispersion
  485 + if (strip_whitespace(filename).length() > 0) {
  486 + m_response->load_edisp(filename);
  487 + }
  488 +
  489 + // Increase number of parameters
435 490 npar[3]++;
436 491 }
437 492  
... ... @@ -444,22 +499,14 @@ void GCTAObservation::read(const GXmlElement&amp; xml)
444 499 ", \"PointSpreadFunction\" and \"EnergyDispersion\" parameters.");
445 500 }
446 501  
447   - // Delete old response function
448   - if (m_response != NULL) delete m_response;
449   -
450   - // Allocate new CTA response function
451   - m_response = new GCTAResponse;
452   -
453   - // Load PSF (need this before Aeff to handle thetacut)
454   - if (strip_whitespace(psf).length() > 0) {
455   - m_response->load_psf(psf);
456   - }
457   -
458   - // Load effective area
459   - if (strip_whitespace(aeff).length() > 0) {
460   - m_response->arf_thetacut(thetacut);
461   - m_response->arf_scale(scale);
462   - m_response->load_arf(aeff);
  502 + // If we have an ARF then apply effective area correction if necessary
  503 + GCTAAeffArf* arf = const_cast<GCTAAeffArf*>(dynamic_cast<const GCTAAeffArf*>(m_response->aeff()));
  504 + if (arf != NULL) {
  505 + if (arf->thetacut() > 0.0) {
  506 + // TODO
  507 + // Loop over all energies in ARF vector, determine PSF integral,
  508 + // and perform correction
  509 + }
463 510 }
464 511  
465 512 // Return
... ... @@ -554,23 +601,52 @@ void GCTAObservation::write(GXmlElement&amp; xml) const
554 601 npar[0]++;
555 602 }
556 603  
557   - // Handle ARF
  604 + // Handle effective area
558 605 else if (par->attribute("name") == "EffectiveArea") {
559   - std::string aeff = "";
  606 +
  607 + // Initialise attributes
  608 + std::string filename = "";
560 609 double thetacut = 0.0;
561 610 double scale = 1.0;
  611 + double sigma = 0.0;
  612 +
  613 + // Continue only if response and effective area are valid
562 614 if (m_response != NULL) {
563   - aeff = m_response->arffile();
564   - thetacut = m_response->arf_thetacut();
565   - scale = m_response->arf_scale();
566   - }
567   - par->attribute("file", aeff);
  615 + if (m_response->aeff() != NULL) {
  616 +
  617 + // Get filename
  618 + filename = m_response->aeff()->filename();
  619 +
  620 + // Get optional ARF attributes
  621 + const GCTAAeffArf* arf =
  622 + dynamic_cast<const GCTAAeffArf*>(m_response->aeff());
  623 + if (arf != NULL) {
  624 + thetacut = arf->thetacut();
  625 + scale = arf->scale();
  626 + sigma = arf->sigma();
  627 + }
  628 +
  629 + // Get optional performance table attributes
  630 + const GCTAAeffPerfTable* perf =
  631 + dynamic_cast<const GCTAAeffPerfTable*>(m_response->aeff());
  632 + if (perf != NULL) {
  633 + sigma = arf->sigma();
  634 + }
  635 +
  636 + } // endif: effective area was valid
  637 + } // endif: response was valid
  638 +
  639 + // Set attributes
  640 + par->attribute("file", filename);
568 641 if (thetacut > 0.0) {
569 642 par->attribute("thetacut", str(thetacut));
570 643 }
571 644 if (scale != 1.0) {
572 645 par->attribute("scale", str(scale));
573 646 }
  647 + if (sigma > 0.0) {
  648 + par->attribute("sigma", str(sigma));
  649 + }
574 650 npar[1]++;
575 651 }
576 652  
... ...
inst/cta/src/GCTAResponse.cpp
  View file @233655f
... ... @@ -47,6 +47,9 @@
47 47 #include "GCTAException.hpp"
48 48 #include "GCTASupport.hpp"
49 49 #include "GCTADir.hpp"
  50 +#include "GCTAAeff2D.hpp"
  51 +#include "GCTAAeffArf.hpp"
  52 +#include "GCTAAeffPerfTable.hpp"
50 53  
51 54 /* __ Method name definitions ____________________________________________ */
52 55 #define G_CALDB "GCTAResponse::caldb(std::string&)"
... ... @@ -587,6 +590,9 @@ void GCTAResponse::load(const std::string&amp; irfname)
587 590 // Read performance table
588 591 read_performance_table(filename);
589 592  
  593 + // Read effective area
  594 + m_aeff = new GCTAAeffPerfTable(filename);
  595 +
590 596 // Store response name
591 597 m_rspname = irfname;
592 598  
... ... @@ -599,29 +605,68 @@ void GCTAResponse::load(const std::string&amp; irfname)
599 605  
600 606  
601 607 /***********************************************************************//**
602   - * @brief Load CTA ARF vector
  608 + * @brief Load effective area
603 609 *
604   - * @param[in] filename FITS file name.
  610 + * @param[in] filename Effective area filename.
  611 + *
  612 + * This method allocates an effective area instance and load the effective
  613 + * area information from a response file. The following response file formats
  614 + * are supported:
  615 + *
  616 + * (1) A CTA performance table. This is an ASCII file which specifies the
  617 + * on-axis effective area as function of energy.
  618 + *
  619 + * (2) A ARF FITS file. This is a FITS file which stores the effective area
  620 + * in a vector.
  621 + *
  622 + * (3) A CTA response table. This is a FITS file which specifies the
  623 + * effective area as function of energy and offset angle.
  624 + *
  625 + * This method examines the file, and depending on the detected format,
  626 + * allocates the appropriate effective area class and loads the data.
605 627 *
606   - * This method loads a CTA ARF vector from a FITS file. See the
607   - * GCTAResponse::read_arf method for more information on ARF vector reading.
  628 + * First, the method checks whether the file is a FITS file or not. If the
  629 + * file is not a FITS file, it is assumed that the file is an ASCII
  630 + * performance table. If the file is a FITS file, the number of rows found
  631 + * in the table is used to distinguish between an ARF (multiple rows) and
  632 + * a CTA response table (single row).
  633 + *
  634 + * @todo Implement a method that checks if a file is a FITS file instead
  635 + * of using try-catch.
608 636 ***************************************************************************/
609   -void GCTAResponse::load_arf(const std::string& filename)
  637 +void GCTAResponse::load_aeff(const std::string& filename)
610 638 {
611   - // Open ARF FITS file
612   - GFits file(filename);
613   -
614   - // Get ARF table
615   - GFitsTable* table = file.table("SPECRESP");
  639 + // Free any existing effective area instance
  640 + if (m_aeff != NULL) delete m_aeff;
  641 + m_aeff = NULL;
  642 +
  643 + // Try opening the file as a FITS file
  644 + try {
  645 +
  646 + // Open FITS file
  647 + GFits file(filename);
  648 +
  649 + // If file contains an "EFFECTIVE AREA" extension then load it
  650 + // as CTA response table
  651 + if (file.hashdu("EFFECTIVE AREA")) {
  652 + file.close();
  653 + m_aeff = new GCTAAeff2D(filename);
  654 + }
616 655  
617   - // Read ARF
618   - read_arf(table);
  656 + // ... else if file contains a "SPECRESP" extension then load it
  657 + // as ARF
  658 + else if (file.hashdu("SPECRESP")) {
  659 + file.close();
  660 + m_aeff = new GCTAAeffArf(filename);
  661 + }
619 662  
620   - // Close ARF FITS file
621   - file.close();
  663 + }
622 664  
623   - // Store filename
624   - m_arffile = filename;
  665 + // If FITS file opening failed then assume that we have a performance
  666 + // table
  667 + catch (GException::fits_open_error &e) {
  668 + m_aeff = new GCTAAeffPerfTable(filename);
  669 + }
625 670  
626 671 // Return
627 672 return;
... ... @@ -697,137 +742,6 @@ void GCTAResponse::load_psf(const std::string&amp; filename)
697 742  
698 743  
699 744 /***********************************************************************//**
700   - * @brief Read CTA ARF vector
701   - *
702   - * @param[in] hdu FITS table pointer.
703   - *
704   - * This method reads a CTA ARF vector from the FITS HDU. Note that the
705   - * energies are converted to TeV and the effective area is converted to cm2.
706   - * Conversion is done based on the units provided for the energy and
707   - * effective area columns. Units that are recognized are 'keV', 'MeV', 'GeV',
708   - * 'TeV', 'm^2', 'm2', 'cm^2' and 'cm^2' (case independent).
709   - *
710   - * @todo Assign appropriate theta angle for PSF. So far we use onaxis.
711   - * For appropriate theta angle assignment, we would need this
712   - * information in the response header.
713   - ***************************************************************************/
714   -void GCTAResponse::read_arf(const GFitsTable* hdu)
715   -{
716   - // Clear arrays
717   - m_aeff_logE.clear();
718   - m_aeff.clear();
719   -
720   - // Get pointers to table columns
721   - const GFitsTableCol* energy_lo = &(*hdu)["ENERG_LO"];
722   - const GFitsTableCol* energy_hi = &(*hdu)["ENERG_HI"];
723   - const GFitsTableCol* specresp = &(*hdu)["SPECRESP"];
724   -
725   - // Determine unit conversion factors (default: TeV and cm^2)
726   - std::string u_energy_lo = tolower(strip_whitespace(energy_lo->unit()));
727   - std::string u_energy_hi = tolower(strip_whitespace(energy_hi->unit()));
728   - std::string u_specresp = tolower(strip_whitespace(specresp->unit()));
729   - double c_energy_lo = 1.0;
730   - double c_energy_hi = 1.0;
731   - double c_specresp = 1.0;
732   - if (u_energy_lo == "kev") {
733   - c_energy_lo = 1.0e-9;
734   - }
735   - else if (u_energy_lo == "mev") {
736   - c_energy_lo = 1.0e-6;
737   - }
738   - else if (u_energy_lo == "gev") {
739   - c_energy_lo = 1.0e-3;
740   - }
741   - if (u_energy_hi == "kev") {
742   - c_energy_hi = 1.0e-9;
743   - }
744   - else if (u_energy_hi == "mev") {
745   - c_energy_hi = 1.0e-6;
746   - }
747   - else if (u_energy_hi == "gev") {
748   - c_energy_hi = 1.0e-3;
749   - }
750   - if (u_specresp == "m^2" || u_specresp == "m2") {
751   - c_specresp = 10000.0;
752   - }
753   -
754   - // Extract number of energy bins
755   - int num = energy_lo->length();
756   -
757   - // Set nodes
758   - for (int i = 0; i < num; ++i) {
759   -
760   - // Compute log10 mean energy in TeV
761   - double e_min = energy_lo->real(i) * c_energy_lo;
762   - double e_max = energy_hi->real(i) * c_energy_hi;
763   - double logE = 0.5 * (log10(e_min) + log10(e_max));
764   -
765   - // Initialise scale factor
766   - double scale = m_arf_scale;
767   -
768   - // Optionally compute scaling factor from thetacut. This is done
769   - // by computing the containment fraction for the specified thetacut.
770   - if (m_arf_thetacut > 0.0) {
771   -
772   - // Get PSF parameters for node energy and theta angle
773   - //TODO: Implement theta angle computation
774   - GCTAPsfPars pars = psf_dummy_sigma(logE, 0.0);
775   -
776   - // Get maximum integration radius
777   - double rmax = m_arf_thetacut * deg2rad;
778   -
779   - // Setup integration kernel
780   - cta_npsf_kern_rad_azsym integrand(this,
781   - rmax,
782   - 0.0,
783   - pars);
784   -
785   - // Setup integration
786   - GIntegral integral(&integrand);
787   - integral.eps(m_eps);
788   -
789   - // Perform integration
790   - double fraction = integral.romb(0.0, rmax);
791   -
792   - // Update scale factor
793   - if (fraction > 0.0) {
794   - scale /= fraction;
795   - #if defined(G_DEBUG_READ_ARF)
796   - std::cout << "GCTAResponse::read_arf:";
797   - std::cout << " e_min=" << e_min;
798   - std::cout << " e_max=" << e_max;
799   - std::cout << " logE=" << logE;
800   - std::cout << " scale=" << scale;
801   - std::cout << " fraction=" << fraction;
802   - std::cout << std::endl;
803   - #endif
804   - }
805   - else {
806   - std::cout << "WARNING: GCTAResponse::read_arf:";
807   - std::cout << " Non-positive integral occured in";
808   - std::cout << " PSF integration in GCTAResponse::read_arf.";
809   - std::cout << std::endl;
810   - }
811   - }
812   -
813   - // Compute effective area in cm2
814   - double aeff = specresp->real(i) * c_specresp * scale;
815   -
816   - // Store log10 mean energy and effective area value
817   - m_aeff_logE.append(logE);
818   - m_aeff.push_back(aeff);
819   -
820   - } // endfor: looped over nodes
821   -
822   - // Disable offset angle dependence
823   - m_offset_sigma = 0.0;
824   -
825   - // Return
826   - return;
827   -}
828   -
829   -
830   -/***********************************************************************//**
831 745 * @brief Read CTA PSF vector
832 746 *
833 747 * @param[in] hdu FITS table pointer.
... ... @@ -925,25 +839,13 @@ std::string GCTAResponse::print(void) const
925 839 result.append("=== GCTAResponse ===");
926 840 result.append("\n"+parformat("Calibration database")+m_caldb);
927 841 result.append("\n"+parformat("Response name")+m_rspname);
928   - result.append("\n"+parformat("ARF file name")+m_arffile);
929 842 result.append("\n"+parformat("RMF file name")+m_rmffile);
930 843 result.append("\n"+parformat("PSF file name")+m_psffile);
931   - result.append("\n"+parformat("ARF theta angle cut"));
932   - if (m_arf_thetacut > 0) {
933   - result.append(str(m_arf_thetacut));
934   - }
935   - else {
936   - result.append("none");
937   - }
938   - result.append("\n"+parformat("ARF scaling")+str(m_arf_scale));
939   - if (m_offset_sigma == 0) {
940   - result.append("\n"+parformat("Offset angle dependence")+"none");
941   - }
942   - else {
943   - std::string txt = "Fixed sigma="+str(m_offset_sigma);
944   - result.append("\n"+parformat("Offset angle dependence")+txt);
  844 +
  845 + // Append effective area information
  846 + if (m_aeff != NULL) {
  847 + result.append("\n"+m_aeff->print());
945 848 }
946   - result.append("\n"+parformat("Effective area nodes")+str(m_aeff_logE.size()));
947 849  
948 850 // Append PSF information
949 851 result.append("\n"+parformat("PSF version"));
... ... @@ -967,15 +869,6 @@ std::string GCTAResponse::print(void) const
967 869 result.append("\n"+parformat("PSF nodes")+str(m_psf_logE.size()));
968 870 }
969 871  
970   - // Debug option: Plot Aeff
971   - #if defined(G_DEBUG_PRINT_AEFF)
972   - result.append("\n"+parformat("Effective area"));
973   - for (int i = 0; i < m_aeff_logE.size(); ++i) {
974   - result.append("\n"+parformat("logE="+str(m_aeff_logE[i])));
975   - result.append("Aeff="+str(m_aeff.at(i))+" m2");
976   - }
977   - #endif
978   -
979 872 // Debug option: Plot PSF
980 873 #if defined(G_DEBUG_PRINT_PSF)
981 874 result.append("\n"+parformat("Point spread function"));
... ... @@ -1547,14 +1440,10 @@ double GCTAResponse::npred_diffuse(const GModelDiffuseSource&amp; model,
1547 1440 * @param[in] zenith Zenith angle of telescope pointing (radians).
1548 1441 * @param[in] azimuth Azimuth angle of telescope pointing (radians).
1549 1442 * @param[in] srcLogEng Log10 of true photon energy (E/TeV).
  1443 + * @return Effective area in units fo cm2.
1550 1444 *
1551   - * A simple offset angle dependence has been implemented using a Gaussian
1552   - * if theta^2 (similar to the radial acceptance model implemented for the
1553   - * background). The width of the Gaussian is assumed independent of energy
1554   - * and given by the fixed parameter m_offset_sigma.
1555   - *
1556   - * @todo Implement offset angle dependence from response file.
1557   - * @todo So far the parameters phi, zenith, and azimuth are not used.
  1445 + * Returns the effective area as function of position in the camera system
  1446 + * and in the earth system.
1558 1447 ***************************************************************************/
1559 1448 double GCTAResponse::aeff(const double& theta,
1560 1449 const double& phi,
... ... @@ -1562,19 +1451,10 @@ double GCTAResponse::aeff(const double&amp; theta,
1562 1451 const double& azimuth,
1563 1452 const double& srcLogEng) const
1564 1453 {
1565   - // Interpolate effective area using node array
1566   - double aeff = m_aeff_logE.interpolate(srcLogEng, m_aeff);
1567   - if (aeff < 0) {
1568   - aeff = 0.0;
1569   - }
1570   -
1571   - // Optionally add in offset angle dependence
1572   - if (m_offset_sigma != 0.0) {
1573   - double offset = theta * rad2deg;
1574   - double arg = offset * offset / m_offset_sigma;
1575   - double scale = exp(-0.5 * arg * arg);
1576   - aeff *= scale;
1577   - }
  1454 + // Get effective area
  1455 + double aeff = (m_aeff != NULL)
  1456 + ? (*m_aeff)(srcLogEng, theta, phi, zenith, azimuth)
  1457 + : 0.0;
1578 1458  
1579 1459 // Return effective area
1580 1460 return aeff;
... ... @@ -2043,23 +1923,21 @@ void GCTAResponse::init_members(void)
2043 1923 // Initialise members
2044 1924 m_caldb.clear();
2045 1925 m_rspname.clear();
2046   - m_arffile.clear();
2047 1926 m_rmffile.clear();
2048 1927 m_psffile.clear();
2049 1928 m_psf_logE.clear();
2050   - m_aeff_logE.clear();
2051 1929 m_logE.clear();
2052   - m_aeff.clear();
2053 1930 m_r68.clear();
2054 1931 m_r80.clear();
2055 1932 m_eps = 1.0e-5; // Precision for Romberg integration
2056 1933 m_offset_sigma = 3.0; // Default value for now ...
2057   - m_arf_thetacut = 0.0; // Default: no thetacut
2058   - m_arf_scale = 1.0; // Default: ARF is unscaled
2059 1934  
2060 1935 // Initialise PSF information
2061 1936 m_psf_version = -99; // Unknown PSF version
2062 1937 m_psf_table.clear();
  1938 +
  1939 + // New factorisation
  1940 + m_aeff = NULL;
2063 1941  
2064 1942 // Return
2065 1943 return;
... ... @@ -2076,24 +1954,22 @@ void GCTAResponse::copy_members(const GCTAResponse&amp; rsp)
2076 1954 // Copy attributes
2077 1955 m_caldb = rsp.m_caldb;
2078 1956 m_rspname = rsp.m_rspname;
2079   - m_arffile = rsp.m_arffile;
2080 1957 m_rmffile = rsp.m_rmffile;
2081 1958 m_psffile = rsp.m_psffile;
2082 1959 m_psf_logE = rsp.m_psf_logE;
2083   - m_aeff_logE = rsp.m_aeff_logE;
2084 1960 m_logE = rsp.m_logE;
2085   - m_aeff = rsp.m_aeff;
2086 1961 m_r68 = rsp.m_r68;
2087 1962 m_r80 = rsp.m_r80;
2088 1963 m_eps = rsp.m_eps;
2089 1964 m_offset_sigma = rsp.m_offset_sigma;
2090   - m_arf_thetacut = rsp.m_arf_thetacut;
2091   - m_arf_scale = rsp.m_arf_scale;
2092 1965  
2093 1966 // Copy PSF information
2094 1967 m_psf_version = rsp.m_psf_version;
2095 1968 m_psf_table = rsp.m_psf_table;
2096 1969  
  1970 + // New factorisation
  1971 + m_aeff = (rsp.m_aeff != NULL) ? rsp.m_aeff->clone() : NULL;
  1972 +
2097 1973 // Return
2098 1974 return;
2099 1975 }
... ... @@ -2104,6 +1980,12 @@ void GCTAResponse::copy_members(const GCTAResponse&amp; rsp)
2104 1980 ***************************************************************************/
2105 1981 void GCTAResponse::free_members(void)
2106 1982 {
  1983 + // Free memory
  1984 + if (m_aeff != NULL) delete m_aeff;
  1985 +
  1986 + // Initialise pointers
  1987 + m_aeff = NULL;
  1988 +
2107 1989 // Return
2108 1990 return;
2109 1991 }
... ... @@ -2125,9 +2007,9 @@ void GCTAResponse::read_performance_table(const std::string&amp; filename)
2125 2007 {
2126 2008 // Clear arrays
2127 2009 m_psf_logE.clear();
2128   - m_aeff_logE.clear();
  2010 + //m_aeff_logE.clear();
2129 2011 m_logE.clear();
2130   - m_aeff.clear();
  2012 + //m_aeff.clear();
2131 2013 m_r68.clear();
2132 2014 m_r80.clear();
2133 2015  
... ... @@ -2164,7 +2046,7 @@ void GCTAResponse::read_performance_table(const std::string&amp; filename)
2164 2046  
2165 2047 // Push elements in vectors
2166 2048 m_logE.push_back(todouble(elements[0]));
2167   - m_aeff.push_back(todouble(elements[1])*10000.0);
  2049 + //m_aeff.push_back(todouble(elements[1])*10000.0);
2168 2050 m_r68.push_back(todouble(elements[2]));
2169 2051 m_r80.push_back(todouble(elements[3]));
2170 2052  
... ... @@ -2175,7 +2057,7 @@ void GCTAResponse::read_performance_table(const std::string&amp; filename)
2175 2057 if (num > 0) {
2176 2058 for (int i = 0; i < num; ++i) {
2177 2059 m_psf_logE.append(m_logE.at(i));
2178   - m_aeff_logE.append(m_logE.at(i));
  2060 + //m_aeff_logE.append(m_logE.at(i));
2179 2061 }
2180 2062 }
2181 2063  
... ...