Samuel Flis / gammalib

Commit 63e6870b722829f42acd1144528e533c5afc3476

Authored by Jürgen Knödlseder
1 parent ca6799cf

Add GCTAInstDir theta() and phi() methods (#2308)

Inline Side-by-side
Showing 9 changed files with 215 additions and 141 deletions
... ... @@ -21,7 +21,10 @@ The following classes have been renamed:
21 21 - none
22 22  
23 23 The following methods have been added:
24   -- gammalib::plaw_integral
  24 +- gammalib::atan2(double&, double&)
  25 +- gammalib::plaw_integral(double&, double&, double&, double&)
  26 +- GCTAInstDir::theta()
  27 +- GCTAInstDir::phi()
25 28  
26 29 The following methods have been removed:
27 30 - none
... ... @@ -68,8 +71,11 @@ None
68 71  
69 72 9. Numerics module
70 73 ------------------
71   -The gammalib::plaw_integral() function has been added to compute the area under
72   -a power law between two boundaries.
  74 +The gammalib::atan2() function was added to compute the arc tangens of a ratio
  75 +in radians (#2308).
  76 +
  77 +The gammalib::plaw_integral() function was added to compute the area under a
  78 +power law between two boundaries (#2308).
73 79  
74 80  
75 81 10. Observation module
... ... @@ -124,7 +130,8 @@ None
124 130  
125 131 20. CTA interface
126 132 -----------------
127   -None
  133 +A GCTAInstDir::theta() and GCTAInstDir::phi() methods were added that return the
  134 +distance to the camera centre in radians and the azimuth angle, respectively (#2308).
128 135  
129 136  
130 137 21. COMPTEL interface
... ...
include/GMath.hpp
  View file @63e6870
... ... @@ -61,6 +61,7 @@ namespace gammalib {
61 61  
62 62 // Trigonometric functions
63 63 double acos(const double& arg);
  64 + double atan2(const double& y, const double& x);
64 65 double cosd(const double& angle);
65 66 double sind(const double& angle);
66 67 double tand(const double& angle);
... ...
inst/cta/include/GCTAInstDir.hpp
  View file @63e6870
1 1 /***************************************************************************
2 2 * GCTAInstDir.hpp - CTA instrument direction class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2014 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -38,8 +38,9 @@
38 38 *
39 39 * @brief CTA instrument direction class.
40 40 *
41   - * The CTA instrument direction is an encapsulation of a sky direction
42   - * as CTA is an imaging device.
  41 + * The CTA instrument direction defines the measured (or reconstructed)
  42 + * direction of an event. The instrument direction is given in sky and
  43 + * detector coordinates.
43 44 ***************************************************************************/
44 45 class GCTAInstDir : public GInstDir {
45 46  
... ... @@ -67,6 +68,8 @@ public:
67 68 void dety(const double &y);
68 69 const double& detx(void) const;
69 70 const double& dety(void) const;
  71 + double theta(void) const;
  72 + double phi(void) const;
70 73  
71 74 protected:
72 75 // Protected methods
... ... @@ -94,6 +97,21 @@ std::string GCTAInstDir::classname(void) const
94 97  
95 98  
96 99 /***********************************************************************//**
  100 + * @brief Set sky direction
  101 + *
  102 + * @param[in] dir Sky direction.
  103 + *
  104 + * Set the sky direction.
  105 + ***************************************************************************/
  106 +inline
  107 +void GCTAInstDir::dir(const GSkyDir& dir)
  108 +{
  109 + m_dir = dir;
  110 + return;
  111 +}
  112 +
  113 +
  114 +/***********************************************************************//**
97 115 * @brief Return reference to sky direction
98 116 *
99 117 * @return Reference to sky direction.
... ... @@ -122,16 +140,31 @@ const GSkyDir& GCTAInstDir::dir(void) const
122 140  
123 141  
124 142 /***********************************************************************//**
125   - * @brief Set sky direction
  143 + * @brief Set DETX coordinate (in radians)
126 144 *
127   - * @param[in] dir Sky direction.
  145 + * @param[in] x DETX coordinate (in radians).
128 146 *
129   - * Set the sky direction.
  147 + * Set DETX coordinate.
130 148 ***************************************************************************/
131 149 inline
132   -void GCTAInstDir::dir(const GSkyDir& dir)
  150 +void GCTAInstDir::detx(const double &x)
133 151 {
134   - m_dir = dir;
  152 + m_detx = x;
  153 + return;
  154 +}
  155 +
  156 +
  157 +/***********************************************************************//**
  158 + * @brief Set DETY coordinate (in radians)
  159 + *
  160 + * @param[in] y DETY coordinate (in radians).
  161 + *
  162 + * Set DETY coordinate.
  163 + ***************************************************************************/
  164 +inline
  165 +void GCTAInstDir::dety(const double &y)
  166 +{
  167 + m_dety = y;
135 168 return;
136 169 }
137 170  
... ... @@ -165,32 +198,41 @@ const double& GCTAInstDir::dety(void) const
165 198  
166 199  
167 200 /***********************************************************************//**
168   - * @brief Set DETX coordinate (in radians)
  201 + * @brief Return offset angle (in radians)
169 202 *
170   - * @param[in] x DETX coordinate (in radians).
  203 + * @return Offset angle (in radians).
171 204 *
172   - * Set DETX coordinate.
  205 + * Returns the offset angle from the camera centre (in radians). The offset
  206 + * angle \f$\theta\f$ is computed using
  207 + *
  208 + * \f[\theta = \sqrt{{\rm DETX}^2 + {\rm DETY}^2}\f]
  209 + *
  210 + * where \f${\rm DETX}\f$ and \f${\rm DETY}\f$ are the detector coordinates.
173 211 ***************************************************************************/
174 212 inline
175   -void GCTAInstDir::detx(const double &x)
  213 +double GCTAInstDir::theta(void) const
176 214 {
177   - m_detx = x;
178   - return;
  215 + return (std::sqrt(m_detx*m_detx + m_dety*m_dety));
179 216 }
180 217  
181 218  
182 219 /***********************************************************************//**
183   - * @brief Set DETY coordinate (in radians)
  220 + * @brief Return azimuth angle (in radians)
184 221 *
185   - * @param[in] y DETY coordinate (in radians).
  222 + * @return Azimuth angle (in radians).
186 223 *
187   - * Set DETY coordinate.
  224 + * Returns the azimuth angle from the camera centre (in radians). The
  225 + * azimuth angle \f$\phi\f$ is computed using
  226 + *
  227 + * \f[\phi = \arctan \left( \frac{{\rm DETY}}{{\rm DETX}} \right)\f]
  228 + *
  229 + * where \f${\rm DETX}\f$ and \f${\rm DETY}\f$ are the detector coordinates.
  230 + * The \c atan2(DETY,DETX) function is used for the computation.
188 231 ***************************************************************************/
189 232 inline
190   -void GCTAInstDir::dety(const double &y)
  233 +double GCTAInstDir::phi(void) const
191 234 {
192   - m_dety = y;
193   - return;
  235 + return (gammalib::atan2(m_dety, m_detx));
194 236 }
195 237  
196 238 #endif /* GCTAINSTDIR_HPP */
... ...
inst/cta/pyext/GCTAInstDir.i
  View file @63e6870
1 1 /***************************************************************************
2 2 * GCTAInstDir.i - CTA instrument direction class *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2010-2015 by Juergen Knoedlseder *
  4 + * copyright (C) 2010-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -55,6 +55,8 @@ public:
55 55 void dety(const double &y);
56 56 const double& detx(void) const;
57 57 const double& dety(void) const;
  58 + double theta(void) const;
  59 + double phi(void) const;
58 60 };
59 61  
60 62  
... ...
inst/cta/src/GCTAModelAeffBackground.cpp
  View file @63e6870
... ... @@ -318,18 +318,14 @@ double GCTAModelAeffBackground::eval(const GEvent& event,
318 318 throw GException::invalid_argument(G_EVAL, msg);
319 319 }
320 320  
321   - // Set DETX and DETY in instrument direction
  321 + // Set instrument direction
322 322 GCTAInstDir inst_dir = cta->pointing().instdir(dir->dir());
323 323  
324   - // Set theta and phi from instrument coordinates
325   - double theta = std::sqrt(inst_dir.detx() * inst_dir.detx() +
326   - inst_dir.dety() * inst_dir.dety());
327   - double phi = gammalib::atan2d(inst_dir.dety(), inst_dir.detx()) *
328   - gammalib::deg2rad;
329   -
330 324 // Evaluate function
331 325 double logE = event.energy().log10TeV();
332   - double spat = (*aeff)(logE, theta, phi,
  326 + double spat = (*aeff)(logE,
  327 + inst_dir.theta(),
  328 + inst_dir.phi(),
333 329 cta->pointing().zenith(),
334 330 cta->pointing().azimuth(), false);
335 331 double spec = (spectral() != NULL)
... ...
inst/cta/test/test_CTA.cpp
  View file @63e6870
... ... @@ -86,6 +86,29 @@ const std::string cta_onoff_onreg = datadir+"/onoff_region_on.reg";
86 86  
87 87  
88 88 /***********************************************************************//**
  89 + * @brief Set miscellaneous CTA classes test methods
  90 + ***************************************************************************/
  91 +void TestGCTA::set(void)
  92 +{
  93 + // Set test name
  94 + name("Miscellaneous CTA classes");
  95 +
  96 + // Append GCTAInstDir tests to test suite
  97 + append(static_cast<pfunction>(&TestGCTA::test_instdir),
  98 + "Test GCTAInstDir class");
  99 +
  100 + // Append GCTAPointing tests to test suite
  101 + append(static_cast<pfunction>(&TestGCTA::test_pointing_load_table),
  102 + "Test loading of table into GCTAPointing instance");
  103 + append(static_cast<pfunction>(&TestGCTA::test_pointing_interpolate_altaz),
  104 + "Test alt/az interpolation given a time");
  105 +
  106 + // Return
  107 + return;
  108 +}
  109 +
  110 +
  111 +/***********************************************************************//**
89 112 * @brief Set CTA response test methods
90 113 ***************************************************************************/
91 114 void TestGCTAResponse::set(void)
... ... @@ -133,18 +156,6 @@ void TestGCTAResponse::set(void)
133 156  
134 157  
135 158 /***********************************************************************//**
136   - * @brief Clone test suite
137   - *
138   - * @return Pointer to deep copy of test suite.
139   - ***************************************************************************/
140   -TestGCTAResponse* TestGCTAResponse::clone(void) const
141   -{
142   - // Clone test suite
143   - return new TestGCTAResponse(*this);
144   -}
145   -
146   -
147   -/***********************************************************************//**
148 159 * @brief Set CTA model test methods
149 160 ***************************************************************************/
150 161 void TestGCTAModel::set(void)
... ... @@ -166,18 +177,6 @@ void TestGCTAModel::set(void)
166 177  
167 178  
168 179 /***********************************************************************//**
169   - * @brief Clone test suite
170   - *
171   - * @return Pointer to deep copy of test suite.
172   - ***************************************************************************/
173   -TestGCTAModel* TestGCTAModel::clone(void) const
174   -{
175   - // Clone test suite
176   - return new TestGCTAModel(*this);
177   -}
178   -
179   -
180   -/***********************************************************************//**
181 180 * @brief Set CTA observation test methods
182 181 ***************************************************************************/
183 182 void TestGCTAObservation::set(void)
... ... @@ -205,18 +204,6 @@ void TestGCTAObservation::set(void)
205 204  
206 205  
207 206 /***********************************************************************//**
208   - * @brief Clone test suite
209   - *
210   - * @return Pointer to deep copy of test suite.
211   - ***************************************************************************/
212   -TestGCTAObservation* TestGCTAObservation::clone(void) const
213   -{
214   - // Clone test suite
215   - return new TestGCTAObservation(*this);
216   -}
217   -
218   -
219   -/***********************************************************************//**
220 207 * @brief Set CTA optimizer test methods
221 208 ***************************************************************************/
222 209 void TestGCTAOptimize::set(void)
... ... @@ -242,33 +229,46 @@ void TestGCTAOptimize::set(void)
242 229  
243 230  
244 231 /***********************************************************************//**
245   - * @brief Clone test suite
246   - *
247   - * @return Pointer to deep copy of test suite.
  232 + * @brief Test GCTAInstDir class
248 233 ***************************************************************************/
249   -TestGCTAOptimize* TestGCTAOptimize::clone(void) const
  234 +void TestGCTA::test_instdir(void)
250 235 {
251   - // Clone test suite
252   - return new TestGCTAOptimize(*this);
253   -}
254   -
255   -
256   -
257   -
258   -/***********************************************************************//**
259   - * @brief Set CTA pointing test methods
260   - ***************************************************************************/
261   -void TestGCTAPointing::set(void)
262   -{
263   - // Set test name
264   - name("GCTAPointing");
265   -
266   - // Append tests to test suite
267   - append(static_cast<pfunction>(&TestGCTAPointing::test_load_table),
268   - "Test load pointing from table");
269   -
270   - append(static_cast<pfunction>(&TestGCTAPointing::test_interpolate_altaz),
271   - "Test alt/az interpolation given a time");
  236 + // Check content of empty instance
  237 + GCTAInstDir instdir1;
  238 + test_value(instdir1.dir().ra(), 0.0, "Right Ascension of empty instance");
  239 + test_value(instdir1.dir().dec(), 0.0, "Declination of empty instance");
  240 + test_value(instdir1.detx(), 0.0, "DETX of empty instance");
  241 + test_value(instdir1.dety(), 0.0, "DETY of empty instance");
  242 + test_value(instdir1.theta(), 0.0, "Theta of empty instance");
  243 + test_value(instdir1.phi(), 0.0, "Phi of empty instance");
  244 +
  245 + // Set content
  246 + double ra = 83.6331 * gammalib::deg2rad;
  247 + double dec = 22.0145 * gammalib::deg2rad;
  248 + double detx = 0.5 * gammalib::deg2rad;
  249 + double dety = 1.5 * gammalib::deg2rad;
  250 + GSkyDir dir;
  251 + dir.radec(ra, dec);
  252 + instdir1.dir(dir);
  253 + instdir1.detx(detx);
  254 + instdir1.dety(dety);
  255 +
  256 + // Check content of filled instance
  257 + test_value(instdir1.dir().ra(), ra, "Right Ascension of filled instance");
  258 + test_value(instdir1.dir().dec(), dec, "Declination of filled instance");
  259 + test_value(instdir1.detx(), detx, "DETX of filled instance");
  260 + test_value(instdir1.dety(), dety, "DETY of filled instance");
  261 + test_value(instdir1.theta(), 0.02759607852, "Theta of filled instance");
  262 + test_value(instdir1.phi(), 1.24904577, "Phi of filled instance");
  263 +
  264 + // Check copying of instance
  265 + GCTAInstDir instdir2(instdir1);
  266 + test_value(instdir2.dir().ra(), ra, "Right Ascension of copied instance");
  267 + test_value(instdir2.dir().dec(), dec, "Declination of copied instance");
  268 + test_value(instdir2.detx(), detx, "DETX of copied instance");
  269 + test_value(instdir2.dety(), dety, "DETY of copied instance");
  270 + test_value(instdir2.theta(), 0.02759607852, "Theta of copied instance");
  271 + test_value(instdir2.phi(), 1.24904577, "Phi of copied instance");
272 272  
273 273 // Return
274 274 return;
... ... @@ -276,23 +276,14 @@ void TestGCTAPointing::set(void)
276 276  
277 277  
278 278 /***********************************************************************//**
279   - * @brief Clone test suite
280   - *
281   - * @return Pointer to deep copy of test suite.
282   - ***************************************************************************/
283   -TestGCTAPointing* TestGCTAPointing::clone(void) const
284   -{
285   - // Clone test suite
286   - return new TestGCTAPointing(*this);
287   -}
288   -
289   -
290   -/***********************************************************************//**
291 279 * @brief Test ability to load a CTA pointing table
292 280 ***************************************************************************/
293   -void TestGCTAPointing::test_load_table(void)
  281 +void TestGCTA::test_pointing_load_table(void)
294 282 {
  283 + // Allocate classes
295 284 GCTAPointing pnt;
  285 +
  286 + // Load pointing table
296 287 pnt.load(cta_point_table);
297 288  
298 289 // Return
... ... @@ -303,12 +294,13 @@ void TestGCTAPointing::test_load_table(void)
303 294 /***********************************************************************//**
304 295 * @brief Test interpolation of alt/az pointing dir as a function of time
305 296 ***************************************************************************/
306   -void TestGCTAPointing::test_interpolate_altaz(void)
  297 +void TestGCTA::test_pointing_interpolate_altaz(void)
307 298 {
308   -
  299 + // Allocate classes
309 300 GCTAObservation run;
  301 + GCTAPointing pnt;
310 302  
311   - GCTAPointing pnt;
  303 + // Load pointing table
312 304 pnt.load(cta_point_table);
313 305  
314 306 // Test an out-of bounds time
... ... @@ -405,7 +397,6 @@ void TestGCTAResponse::test_response_aeff(void)
405 397 for (int i = 0; i < 30; ++i) {
406 398 eng.TeV(pow(10.0, -1.7 + 0.1*double(i)));
407 399 double aeff = rsp.aeff(0.0, 0.0, 0.0, 0.0, eng.log10TeV());
408   - //std::cout << eng << " " << eng.log10TeV() << " " << aeff << std::endl;
409 400 sum += aeff;
410 401 }
411 402 test_value(sum, ref, 0.1, "Effective area verification");
... ... @@ -2299,17 +2290,17 @@ int main(void)
2299 2290 bool success = true;
2300 2291  
2301 2292 // Create test suites and append them to the container
  2293 + TestGCTA test;
2302 2294 TestGCTAResponse rsp;
2303 2295 TestGCTAModel model;
2304 2296 TestGCTAOptimize opt;
2305 2297 TestGCTAObservation obs;
2306   - TestGCTAPointing pnt;
2307 2298 testsuites.append(rsp);
2308 2299 if (has_data) {
  2300 + testsuites.append(test);
2309 2301 testsuites.append(model);
2310 2302 testsuites.append(opt);
2311 2303 testsuites.append(obs);
2312   - testsuites.append(pnt);
2313 2304 }
2314 2305  
2315 2306 // Run the testsuites
... ...
inst/cta/test/test_CTA.hpp
  View file @63e6870
1 1 /***************************************************************************
2 2 * test_CTA.hpp - Test CTA classes *
3 3 * ----------------------------------------------------------------------- *
4   - * copyright (C) 2012-2017 by Juergen Knoedlseder *
  4 + * copyright (C) 2012-2018 by Juergen Knoedlseder *
5 5 * ----------------------------------------------------------------------- *
6 6 * *
7 7 * This program is free software: you can redistribute it and/or modify *
... ... @@ -32,6 +32,30 @@
32 32  
33 33  
34 34 /***********************************************************************//**
  35 + * @class TestGCTA
  36 + *
  37 + * @brief Test suite for miscellaneous CTA classes
  38 + *
  39 + * This class defines a unit test suite for miscellaneous CTA classes
  40 + ***************************************************************************/
  41 +class TestGCTA : public GTestSuite {
  42 +public:
  43 +
  44 + // Constructors and destructors
  45 + TestGCTA(void) : GTestSuite() {}
  46 + virtual ~TestGCTA(void) {}
  47 +
  48 + // Methods
  49 + virtual void set(void);
  50 + virtual TestGCTA* clone(void) const { return new TestGCTA(*this); }
  51 + virtual std::string classname(void) const { return "TestGCTA"; }
  52 + void test_instdir(void);
  53 + void test_pointing_load_table(void);
  54 + void test_pointing_interpolate_altaz(void);
  55 +};
  56 +
  57 +
  58 +/***********************************************************************//**
35 59 * @class TestGCTAResponse
36 60 *
37 61 * @brief Test suite for CTA response class testing
... ... @@ -46,7 +70,7 @@ public:
46 70  
47 71 // Methods
48 72 virtual void set(void);
49   - virtual TestGCTAResponse* clone(void) const;
  73 + virtual TestGCTAResponse* clone(void) const { return new TestGCTAResponse(*this); }
50 74 virtual std::string classname(void) const { return "TestGCTAResponse"; }
51 75 void test_response(void);
52 76 void test_response_aeff(void);
... ... @@ -90,7 +114,7 @@ public:
90 114  
91 115 // Methods
92 116 virtual void set(void);
93   - virtual TestGCTAModel* clone(void) const;
  117 + virtual TestGCTAModel* clone(void) const { return new TestGCTAModel(*this); }
94 118 virtual std::string classname(void) const { return "TestGCTAModel"; }
95 119 void test_model_cube_bgd(void);
96 120 void test_model_irf_bgd(void);
... ... @@ -113,7 +137,7 @@ public:
113 137  
114 138 // Methods
115 139 virtual void set(void);
116   - virtual TestGCTAObservation* clone(void) const;
  140 + virtual TestGCTAObservation* clone(void) const { return new TestGCTAObservation(*this); }
117 141 virtual std::string classname(void) const { return "TestGCTAObservation"; }
118 142 void test_event_bin(void);
119 143 void test_event_cube(void);
... ... @@ -139,7 +163,7 @@ public:
139 163  
140 164 // Methods
141 165 virtual void set(void);
142   - virtual TestGCTAOptimize* clone(void) const;
  166 + virtual TestGCTAOptimize* clone(void) const { return new TestGCTAOptimize(*this); }
143 167 virtual std::string classname(void) const { return "TestGCTAOptimize"; }
144 168 void test_unbinned_optimizer(void);
145 169 void test_binned_optimizer(void);
... ... @@ -152,27 +176,4 @@ protected:
152 176 void check_results(const GObservations& obs, const double* results);
153 177 };
154 178  
155   -
156   -/***********************************************************************//**
157   - * @class TestGCTAPointing
158   - *
159   - * @brief Test suite for CTA pointing class testing
160   - *
161   - * This class defines a unit test suite for the CTA response class.
162   - ***************************************************************************/
163   -class TestGCTAPointing : public GTestSuite {
164   -public:
165   -
166   - // Constructors and destructors
167   - TestGCTAPointing(void) : GTestSuite() {}
168   - virtual ~TestGCTAPointing(void) {}
169   -
170   - // Methods
171   - virtual void set(void);
172   - virtual TestGCTAPointing* clone(void) const;
173   - virtual std::string classname(void) const { return "TestGCTAPointing"; }
174   - void test_load_table(void);
175   - void test_interpolate_altaz(void);
176   -};
177   -
178 179 #endif /* TEST_CTA_HPP */
... ...
pyext/GMath.i
  View file @63e6870
... ... @@ -57,6 +57,7 @@ namespace gammalib {
57 57  
58 58 // Trigonometric functions
59 59 double acos(const double& arg);
  60 + double atan2(const double& y, const double& x);
60 61 double cosd(const double& angle);
61 62 double sind(const double& angle);
62 63 double tand(const double& angle);
... ...
src/numerics/GMath.cpp
  View file @63e6870
... ... @@ -92,6 +92,39 @@ double gammalib::acos(const double&amp; arg)
92 92  
93 93  
94 94 /***********************************************************************//**
  95 + * @brief Compute arc tangens in radians
  96 + *
  97 + * @param[in] y Nominator
  98 + * @param[in] x Denominator
  99 + *
  100 + * This code has been adapted from the WCSLIB function wcstrig.c::atan2d().
  101 + ***************************************************************************/
  102 +double gammalib::atan2(const double& y, const double& x)
  103 +{
  104 + // Check for rounding errors
  105 + if (y == 0.0) {
  106 + if (x >= 0.0) {
  107 + return 0.0;
  108 + }
  109 + else if (x < 0.0) {
  110 + return pi;
  111 + }
  112 + }
  113 + else if (x == 0.0) {
  114 + if (y > 0.0) {
  115 + return pihalf;
  116 + }
  117 + else if (y < 0.0) {
  118 + return -pihalf;
  119 + }
  120 + }
  121 +
  122 + // Return arc sine
  123 + return (std::atan2(y,x));
  124 +}
  125 +
  126 +
  127 +/***********************************************************************//**
95 128 * @brief Compute cosine of angle in degrees
96 129 *
97 130 * @param[in] angle Angle in degrees
... ...