From 63e6870b722829f42acd1144528e533c5afc3476 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Fri, 9 Feb 2018 10:29:59 +0100 Subject: [PATCH] Add GCTAInstDir theta() and phi() methods (#2308) --- NEWS | 15 +++++++++++---- include/GMath.hpp | 1 + inst/cta/include/GCTAInstDir.hpp | 82 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------- inst/cta/pyext/GCTAInstDir.i | 4 +++- inst/cta/src/GCTAModelAeffBackground.cpp | 12 ++++-------- inst/cta/test/test_CTA.cpp | 151 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------------------------------------------------------------------- inst/cta/test/test_CTA.hpp | 57 +++++++++++++++++++++++++++++---------------------------- pyext/GMath.i | 1 + src/numerics/GMath.cpp | 33 +++++++++++++++++++++++++++++++++ 9 files changed, 215 insertions(+), 141 deletions(-) diff --git a/NEWS b/NEWS index f9ba2e7..cbdf5a7 100644 --- a/NEWS +++ b/NEWS @@ -21,7 +21,10 @@ The following classes have been renamed: - none The following methods have been added: -- gammalib::plaw_integral +- gammalib::atan2(double&, double&) +- gammalib::plaw_integral(double&, double&, double&, double&) +- GCTAInstDir::theta() +- GCTAInstDir::phi() The following methods have been removed: - none @@ -68,8 +71,11 @@ None 9. Numerics module ------------------ -The gammalib::plaw_integral() function has been added to compute the area under -a power law between two boundaries. +The gammalib::atan2() function was added to compute the arc tangens of a ratio +in radians (#2308). + +The gammalib::plaw_integral() function was added to compute the area under a +power law between two boundaries (#2308). 10. Observation module @@ -124,7 +130,8 @@ None 20. CTA interface ----------------- -None +A GCTAInstDir::theta() and GCTAInstDir::phi() methods were added that return the +distance to the camera centre in radians and the azimuth angle, respectively (#2308). 21. COMPTEL interface diff --git a/include/GMath.hpp b/include/GMath.hpp index 4308dd1..2d88585 100644 --- a/include/GMath.hpp +++ b/include/GMath.hpp @@ -61,6 +61,7 @@ namespace gammalib { // Trigonometric functions double acos(const double& arg); + double atan2(const double& y, const double& x); double cosd(const double& angle); double sind(const double& angle); double tand(const double& angle); diff --git a/inst/cta/include/GCTAInstDir.hpp b/inst/cta/include/GCTAInstDir.hpp index 0f5e435..03f7967 100644 --- a/inst/cta/include/GCTAInstDir.hpp +++ b/inst/cta/include/GCTAInstDir.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAInstDir.hpp - CTA instrument direction class * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2014 by Juergen Knoedlseder * + * copyright (C) 2010-2018 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -38,8 +38,9 @@ * * @brief CTA instrument direction class. * - * The CTA instrument direction is an encapsulation of a sky direction - * as CTA is an imaging device. + * The CTA instrument direction defines the measured (or reconstructed) + * direction of an event. The instrument direction is given in sky and + * detector coordinates. ***************************************************************************/ class GCTAInstDir : public GInstDir { @@ -67,6 +68,8 @@ public: void dety(const double &y); const double& detx(void) const; const double& dety(void) const; + double theta(void) const; + double phi(void) const; protected: // Protected methods @@ -94,6 +97,21 @@ std::string GCTAInstDir::classname(void) const /***********************************************************************//** + * @brief Set sky direction + * + * @param[in] dir Sky direction. + * + * Set the sky direction. + ***************************************************************************/ +inline +void GCTAInstDir::dir(const GSkyDir& dir) +{ + m_dir = dir; + return; +} + + +/***********************************************************************//** * @brief Return reference to sky direction * * @return Reference to sky direction. @@ -122,16 +140,31 @@ const GSkyDir& GCTAInstDir::dir(void) const /***********************************************************************//** - * @brief Set sky direction + * @brief Set DETX coordinate (in radians) * - * @param[in] dir Sky direction. + * @param[in] x DETX coordinate (in radians). * - * Set the sky direction. + * Set DETX coordinate. ***************************************************************************/ inline -void GCTAInstDir::dir(const GSkyDir& dir) +void GCTAInstDir::detx(const double &x) { - m_dir = dir; + m_detx = x; + return; +} + + +/***********************************************************************//** + * @brief Set DETY coordinate (in radians) + * + * @param[in] y DETY coordinate (in radians). + * + * Set DETY coordinate. + ***************************************************************************/ +inline +void GCTAInstDir::dety(const double &y) +{ + m_dety = y; return; } @@ -165,32 +198,41 @@ const double& GCTAInstDir::dety(void) const /***********************************************************************//** - * @brief Set DETX coordinate (in radians) + * @brief Return offset angle (in radians) * - * @param[in] x DETX coordinate (in radians). + * @return Offset angle (in radians). * - * Set DETX coordinate. + * Returns the offset angle from the camera centre (in radians). The offset + * angle \f$\theta\f$ is computed using + * + * \f[\theta = \sqrt{{\rm DETX}^2 + {\rm DETY}^2}\f] + * + * where \f${\rm DETX}\f$ and \f${\rm DETY}\f$ are the detector coordinates. ***************************************************************************/ inline -void GCTAInstDir::detx(const double &x) +double GCTAInstDir::theta(void) const { - m_detx = x; - return; + return (std::sqrt(m_detx*m_detx + m_dety*m_dety)); } /***********************************************************************//** - * @brief Set DETY coordinate (in radians) + * @brief Return azimuth angle (in radians) * - * @param[in] y DETY coordinate (in radians). + * @return Azimuth angle (in radians). * - * Set DETY coordinate. + * Returns the azimuth angle from the camera centre (in radians). The + * azimuth angle \f$\phi\f$ is computed using + * + * \f[\phi = \arctan \left( \frac{{\rm DETY}}{{\rm DETX}} \right)\f] + * + * where \f${\rm DETX}\f$ and \f${\rm DETY}\f$ are the detector coordinates. + * The \c atan2(DETY,DETX) function is used for the computation. ***************************************************************************/ inline -void GCTAInstDir::dety(const double &y) +double GCTAInstDir::phi(void) const { - m_dety = y; - return; + return (gammalib::atan2(m_dety, m_detx)); } #endif /* GCTAINSTDIR_HPP */ diff --git a/inst/cta/pyext/GCTAInstDir.i b/inst/cta/pyext/GCTAInstDir.i index ac2c2cc..2b45639 100644 --- a/inst/cta/pyext/GCTAInstDir.i +++ b/inst/cta/pyext/GCTAInstDir.i @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAInstDir.i - CTA instrument direction class * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2015 by Juergen Knoedlseder * + * copyright (C) 2010-2018 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -55,6 +55,8 @@ public: void dety(const double &y); const double& detx(void) const; const double& dety(void) const; + double theta(void) const; + double phi(void) const; }; diff --git a/inst/cta/src/GCTAModelAeffBackground.cpp b/inst/cta/src/GCTAModelAeffBackground.cpp index ba139b7..0a01d7d 100644 --- a/inst/cta/src/GCTAModelAeffBackground.cpp +++ b/inst/cta/src/GCTAModelAeffBackground.cpp @@ -318,18 +318,14 @@ double GCTAModelAeffBackground::eval(const GEvent& event, throw GException::invalid_argument(G_EVAL, msg); } - // Set DETX and DETY in instrument direction + // Set instrument direction GCTAInstDir inst_dir = cta->pointing().instdir(dir->dir()); - // Set theta and phi from instrument coordinates - double theta = std::sqrt(inst_dir.detx() * inst_dir.detx() + - inst_dir.dety() * inst_dir.dety()); - double phi = gammalib::atan2d(inst_dir.dety(), inst_dir.detx()) * - gammalib::deg2rad; - // Evaluate function double logE = event.energy().log10TeV(); - double spat = (*aeff)(logE, theta, phi, + double spat = (*aeff)(logE, + inst_dir.theta(), + inst_dir.phi(), cta->pointing().zenith(), cta->pointing().azimuth(), false); double spec = (spectral() != NULL) diff --git a/inst/cta/test/test_CTA.cpp b/inst/cta/test/test_CTA.cpp index 628fc79..a7fca15 100644 --- a/inst/cta/test/test_CTA.cpp +++ b/inst/cta/test/test_CTA.cpp @@ -86,6 +86,29 @@ const std::string cta_onoff_onreg = datadir+"/onoff_region_on.reg"; /***********************************************************************//** + * @brief Set miscellaneous CTA classes test methods + ***************************************************************************/ +void TestGCTA::set(void) +{ + // Set test name + name("Miscellaneous CTA classes"); + + // Append GCTAInstDir tests to test suite + append(static_cast(&TestGCTA::test_instdir), + "Test GCTAInstDir class"); + + // Append GCTAPointing tests to test suite + append(static_cast(&TestGCTA::test_pointing_load_table), + "Test loading of table into GCTAPointing instance"); + append(static_cast(&TestGCTA::test_pointing_interpolate_altaz), + "Test alt/az interpolation given a time"); + + // Return + return; +} + + +/***********************************************************************//** * @brief Set CTA response test methods ***************************************************************************/ void TestGCTAResponse::set(void) @@ -133,18 +156,6 @@ void TestGCTAResponse::set(void) /***********************************************************************//** - * @brief Clone test suite - * - * @return Pointer to deep copy of test suite. - ***************************************************************************/ -TestGCTAResponse* TestGCTAResponse::clone(void) const -{ - // Clone test suite - return new TestGCTAResponse(*this); -} - - -/***********************************************************************//** * @brief Set CTA model test methods ***************************************************************************/ void TestGCTAModel::set(void) @@ -166,18 +177,6 @@ void TestGCTAModel::set(void) /***********************************************************************//** - * @brief Clone test suite - * - * @return Pointer to deep copy of test suite. - ***************************************************************************/ -TestGCTAModel* TestGCTAModel::clone(void) const -{ - // Clone test suite - return new TestGCTAModel(*this); -} - - -/***********************************************************************//** * @brief Set CTA observation test methods ***************************************************************************/ void TestGCTAObservation::set(void) @@ -205,18 +204,6 @@ void TestGCTAObservation::set(void) /***********************************************************************//** - * @brief Clone test suite - * - * @return Pointer to deep copy of test suite. - ***************************************************************************/ -TestGCTAObservation* TestGCTAObservation::clone(void) const -{ - // Clone test suite - return new TestGCTAObservation(*this); -} - - -/***********************************************************************//** * @brief Set CTA optimizer test methods ***************************************************************************/ void TestGCTAOptimize::set(void) @@ -242,33 +229,46 @@ void TestGCTAOptimize::set(void) /***********************************************************************//** - * @brief Clone test suite - * - * @return Pointer to deep copy of test suite. + * @brief Test GCTAInstDir class ***************************************************************************/ -TestGCTAOptimize* TestGCTAOptimize::clone(void) const +void TestGCTA::test_instdir(void) { - // Clone test suite - return new TestGCTAOptimize(*this); -} - - - - -/***********************************************************************//** - * @brief Set CTA pointing test methods - ***************************************************************************/ -void TestGCTAPointing::set(void) -{ - // Set test name - name("GCTAPointing"); - - // Append tests to test suite - append(static_cast(&TestGCTAPointing::test_load_table), - "Test load pointing from table"); - - append(static_cast(&TestGCTAPointing::test_interpolate_altaz), - "Test alt/az interpolation given a time"); + // Check content of empty instance + GCTAInstDir instdir1; + test_value(instdir1.dir().ra(), 0.0, "Right Ascension of empty instance"); + test_value(instdir1.dir().dec(), 0.0, "Declination of empty instance"); + test_value(instdir1.detx(), 0.0, "DETX of empty instance"); + test_value(instdir1.dety(), 0.0, "DETY of empty instance"); + test_value(instdir1.theta(), 0.0, "Theta of empty instance"); + test_value(instdir1.phi(), 0.0, "Phi of empty instance"); + + // Set content + double ra = 83.6331 * gammalib::deg2rad; + double dec = 22.0145 * gammalib::deg2rad; + double detx = 0.5 * gammalib::deg2rad; + double dety = 1.5 * gammalib::deg2rad; + GSkyDir dir; + dir.radec(ra, dec); + instdir1.dir(dir); + instdir1.detx(detx); + instdir1.dety(dety); + + // Check content of filled instance + test_value(instdir1.dir().ra(), ra, "Right Ascension of filled instance"); + test_value(instdir1.dir().dec(), dec, "Declination of filled instance"); + test_value(instdir1.detx(), detx, "DETX of filled instance"); + test_value(instdir1.dety(), dety, "DETY of filled instance"); + test_value(instdir1.theta(), 0.02759607852, "Theta of filled instance"); + test_value(instdir1.phi(), 1.24904577, "Phi of filled instance"); + + // Check copying of instance + GCTAInstDir instdir2(instdir1); + test_value(instdir2.dir().ra(), ra, "Right Ascension of copied instance"); + test_value(instdir2.dir().dec(), dec, "Declination of copied instance"); + test_value(instdir2.detx(), detx, "DETX of copied instance"); + test_value(instdir2.dety(), dety, "DETY of copied instance"); + test_value(instdir2.theta(), 0.02759607852, "Theta of copied instance"); + test_value(instdir2.phi(), 1.24904577, "Phi of copied instance"); // Return return; @@ -276,23 +276,14 @@ void TestGCTAPointing::set(void) /***********************************************************************//** - * @brief Clone test suite - * - * @return Pointer to deep copy of test suite. - ***************************************************************************/ -TestGCTAPointing* TestGCTAPointing::clone(void) const -{ - // Clone test suite - return new TestGCTAPointing(*this); -} - - -/***********************************************************************//** * @brief Test ability to load a CTA pointing table ***************************************************************************/ -void TestGCTAPointing::test_load_table(void) +void TestGCTA::test_pointing_load_table(void) { + // Allocate classes GCTAPointing pnt; + + // Load pointing table pnt.load(cta_point_table); // Return @@ -303,12 +294,13 @@ void TestGCTAPointing::test_load_table(void) /***********************************************************************//** * @brief Test interpolation of alt/az pointing dir as a function of time ***************************************************************************/ -void TestGCTAPointing::test_interpolate_altaz(void) +void TestGCTA::test_pointing_interpolate_altaz(void) { - + // Allocate classes GCTAObservation run; + GCTAPointing pnt; - GCTAPointing pnt; + // Load pointing table pnt.load(cta_point_table); // Test an out-of bounds time @@ -405,7 +397,6 @@ void TestGCTAResponse::test_response_aeff(void) for (int i = 0; i < 30; ++i) { eng.TeV(pow(10.0, -1.7 + 0.1*double(i))); double aeff = rsp.aeff(0.0, 0.0, 0.0, 0.0, eng.log10TeV()); - //std::cout << eng << " " << eng.log10TeV() << " " << aeff << std::endl; sum += aeff; } test_value(sum, ref, 0.1, "Effective area verification"); @@ -2299,17 +2290,17 @@ int main(void) bool success = true; // Create test suites and append them to the container + TestGCTA test; TestGCTAResponse rsp; TestGCTAModel model; TestGCTAOptimize opt; TestGCTAObservation obs; - TestGCTAPointing pnt; testsuites.append(rsp); if (has_data) { + testsuites.append(test); testsuites.append(model); testsuites.append(opt); testsuites.append(obs); - testsuites.append(pnt); } // Run the testsuites diff --git a/inst/cta/test/test_CTA.hpp b/inst/cta/test/test_CTA.hpp index fe6c850..f21aa0f 100644 --- a/inst/cta/test/test_CTA.hpp +++ b/inst/cta/test/test_CTA.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * test_CTA.hpp - Test CTA classes * * ----------------------------------------------------------------------- * - * copyright (C) 2012-2017 by Juergen Knoedlseder * + * copyright (C) 2012-2018 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -32,6 +32,30 @@ /***********************************************************************//** + * @class TestGCTA + * + * @brief Test suite for miscellaneous CTA classes + * + * This class defines a unit test suite for miscellaneous CTA classes + ***************************************************************************/ +class TestGCTA : public GTestSuite { +public: + + // Constructors and destructors + TestGCTA(void) : GTestSuite() {} + virtual ~TestGCTA(void) {} + + // Methods + virtual void set(void); + virtual TestGCTA* clone(void) const { return new TestGCTA(*this); } + virtual std::string classname(void) const { return "TestGCTA"; } + void test_instdir(void); + void test_pointing_load_table(void); + void test_pointing_interpolate_altaz(void); +}; + + +/***********************************************************************//** * @class TestGCTAResponse * * @brief Test suite for CTA response class testing @@ -46,7 +70,7 @@ public: // Methods virtual void set(void); - virtual TestGCTAResponse* clone(void) const; + virtual TestGCTAResponse* clone(void) const { return new TestGCTAResponse(*this); } virtual std::string classname(void) const { return "TestGCTAResponse"; } void test_response(void); void test_response_aeff(void); @@ -90,7 +114,7 @@ public: // Methods virtual void set(void); - virtual TestGCTAModel* clone(void) const; + virtual TestGCTAModel* clone(void) const { return new TestGCTAModel(*this); } virtual std::string classname(void) const { return "TestGCTAModel"; } void test_model_cube_bgd(void); void test_model_irf_bgd(void); @@ -113,7 +137,7 @@ public: // Methods virtual void set(void); - virtual TestGCTAObservation* clone(void) const; + virtual TestGCTAObservation* clone(void) const { return new TestGCTAObservation(*this); } virtual std::string classname(void) const { return "TestGCTAObservation"; } void test_event_bin(void); void test_event_cube(void); @@ -139,7 +163,7 @@ public: // Methods virtual void set(void); - virtual TestGCTAOptimize* clone(void) const; + virtual TestGCTAOptimize* clone(void) const { return new TestGCTAOptimize(*this); } virtual std::string classname(void) const { return "TestGCTAOptimize"; } void test_unbinned_optimizer(void); void test_binned_optimizer(void); @@ -152,27 +176,4 @@ protected: void check_results(const GObservations& obs, const double* results); }; - -/***********************************************************************//** - * @class TestGCTAPointing - * - * @brief Test suite for CTA pointing class testing - * - * This class defines a unit test suite for the CTA response class. - ***************************************************************************/ -class TestGCTAPointing : public GTestSuite { -public: - - // Constructors and destructors - TestGCTAPointing(void) : GTestSuite() {} - virtual ~TestGCTAPointing(void) {} - - // Methods - virtual void set(void); - virtual TestGCTAPointing* clone(void) const; - virtual std::string classname(void) const { return "TestGCTAPointing"; } - void test_load_table(void); - void test_interpolate_altaz(void); -}; - #endif /* TEST_CTA_HPP */ diff --git a/pyext/GMath.i b/pyext/GMath.i index 05aae7a..b43dcb8 100644 --- a/pyext/GMath.i +++ b/pyext/GMath.i @@ -57,6 +57,7 @@ namespace gammalib { // Trigonometric functions double acos(const double& arg); + double atan2(const double& y, const double& x); double cosd(const double& angle); double sind(const double& angle); double tand(const double& angle); diff --git a/src/numerics/GMath.cpp b/src/numerics/GMath.cpp index 6d0b252..e4d6b99 100644 --- a/src/numerics/GMath.cpp +++ b/src/numerics/GMath.cpp @@ -92,6 +92,39 @@ double gammalib::acos(const double& arg) /***********************************************************************//** + * @brief Compute arc tangens in radians + * + * @param[in] y Nominator + * @param[in] x Denominator + * + * This code has been adapted from the WCSLIB function wcstrig.c::atan2d(). + ***************************************************************************/ +double gammalib::atan2(const double& y, const double& x) +{ + // Check for rounding errors + if (y == 0.0) { + if (x >= 0.0) { + return 0.0; + } + else if (x < 0.0) { + return pi; + } + } + else if (x == 0.0) { + if (y > 0.0) { + return pihalf; + } + else if (y < 0.0) { + return -pihalf; + } + } + + // Return arc sine + return (std::atan2(y,x)); +} + + +/***********************************************************************//** * @brief Compute cosine of angle in degrees * * @param[in] angle Angle in degrees -- libgit2 0.21.4