From 4bd3d7d42bea706610c46aaf84bb7fe04344272a Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Tue, 21 Jun 2016 23:28:54 +0200 Subject: [PATCH] GCTAAeff2D no longer extrapolates beyond response table boundaries (#1784) --- ChangeLog | 1 + NEWS | 3 +++ inst/cta/include/GCTAAeff2D.hpp | 22 +++++++++------------- inst/cta/src/GCTAAeff2D.cpp | 155 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++-------------------------------- inst/cta/test/test_CTA.py | 26 ++++++++++++++++++++++---- 5 files changed, 158 insertions(+), 49 deletions(-) diff --git a/ChangeLog b/ChangeLog index 2880c72..a1aad15 100644 --- a/ChangeLog +++ b/ChangeLog @@ -3,6 +3,7 @@ * Version 1.1.0 released ========================= + GCTAAeff2D no longer extrapolates beyond response table (#1784) Add GWcsSIN projection class (#1798) CTA cube response no longer selects on RoI and energies (#1789) Add weights to CTA event cube and event bin (#1789) diff --git a/NEWS b/NEWS index 8060c05..3d17e25 100644 --- a/NEWS +++ b/NEWS @@ -416,6 +416,9 @@ The cube response classes GCTACubeExposure, GCTACubePsf, GCTACubeEdisp and GCTACubeBackground do no longer select on RoI and energy boundaries as this is now handled by the weight array (#1789). +The GCTAAeff2D class no longer extrapolates the effective area beyond the boundaries +of the response table, but sets the effective area values to zero (#1784). + 21. COMPTEL interface --------------------- diff --git a/inst/cta/include/GCTAAeff2D.hpp b/inst/cta/include/GCTAAeff2D.hpp index aa92a9e..5127074 100644 --- a/inst/cta/include/GCTAAeff2D.hpp +++ b/inst/cta/include/GCTAAeff2D.hpp @@ -89,6 +89,8 @@ private: void init_members(void); void copy_members(const GCTAAeff2D& aeff); void free_members(void); + void set_indices(void); + void set_boundaries(void); // Members GFilename m_filename; //!< Name of Aeff response file @@ -97,6 +99,10 @@ private: int m_inx_theta; //!< Theta index int m_inx_aeff; //!< Effective area (true energy) int m_inx_aeff_reco; //!< Effective area (reconstructed energy) + double m_logE_min; //!< Minimum logE (log10(E/TeV)) + double m_logE_max; //!< Maximum logE (log10(E/TeV)) + double m_theta_min; //!< Minimum theta (radians) + double m_theta_max; //!< Maximum theta (radians) }; @@ -120,7 +126,6 @@ std::string GCTAAeff2D::classname(void) const inline GFilename GCTAAeff2D::filename(void) const { - // Return filename return m_filename; } @@ -129,6 +134,9 @@ GFilename GCTAAeff2D::filename(void) const * @brief Return response table * * @return Response table. + * + * Returns the response table of the effective area. The effective area + * values are given in units of cm2. ***************************************************************************/ inline const GCTAResponseTable& GCTAAeff2D::table(void) const @@ -136,16 +144,4 @@ const GCTAResponseTable& GCTAAeff2D::table(void) const return m_aeff; } - -/***********************************************************************//** - * @brief Assign response table - * - * @param[in] table Response table. - ***************************************************************************/ -inline -void GCTAAeff2D::table(const GCTAResponseTable& table) -{ - m_aeff = table; -} - #endif /* GCTAAEFF2D_HPP */ diff --git a/inst/cta/src/GCTAAeff2D.cpp b/inst/cta/src/GCTAAeff2D.cpp index 25e9bf9..1eee94a 100644 --- a/inst/cta/src/GCTAAeff2D.cpp +++ b/inst/cta/src/GCTAAeff2D.cpp @@ -29,6 +29,7 @@ #include #endif #include "GTools.hpp" +#include "GMath.hpp" #include "GException.hpp" #include "GFilename.hpp" #include "GFits.hpp" @@ -170,13 +171,14 @@ GCTAAeff2D& GCTAAeff2D::operator=(const GCTAAeff2D& aeff) * @param[in] phi Azimuth angle in camera system (rad). * @param[in] zenith Zenith angle in Earth system (rad). * @param[in] azimuth Azimuth angle in Earth system (rad). - * @param[in] etrue Use true energy (default: true). + * @param[in] etrue Use true energy? * @return Effective area in cm2. * * Returns the effective area in units of cm2 for a given energy and * offset angle. The effective area is bi-linearly interpolated in the * log10(energy) - offset angle plane. The method assures that the effective - * area value never becomes negative. + * area value never becomes negative. Outside the logE - theta range covered + * by the response table the effective area will be set to zero. * * The method supports true and reconstructed energies for logE. To access * the effective area as function of true energy, specify etrue=true @@ -190,17 +192,26 @@ double GCTAAeff2D::operator()(const double& logE, const double& azimuth, const bool& etrue) const { - // Set parameter index - int index = (etrue) ? m_inx_aeff : m_inx_aeff_reco; + // Initialise effective area + double aeff(0.0); - // Get effective area value in cm2 - double aeff = (m_inx_energy == 0) ? m_aeff(index, logE, theta) : - m_aeff(index, theta, logE); + // Continue only if logE and theta are in validity range + if ((logE >= m_logE_min) && (logE <= m_logE_max) && + (theta >= m_theta_min) && (theta <= m_theta_max)) { - // Make sure that effective area is not negative - if (aeff < 0.0) { - aeff = 0.0; - } + // Set parameter index + int index = (etrue) ? m_inx_aeff : m_inx_aeff_reco; + + // Get effective area value in cm2 + aeff = (m_inx_energy == 0) ? m_aeff(index, logE, theta) : + m_aeff(index, theta, logE); + + // Make sure that effective area is not negative + if (aeff < 0.0) { + aeff = 0.0; + } + + } // endif: logE and theta in validity range // Return effective area value return aeff; @@ -281,15 +292,6 @@ void GCTAAeff2D::read(const GFitsTable& table) // Read effective area table m_aeff.read(table); - // Get mandatory indices (throw exception if not found) - m_inx_energy = m_aeff.axis("ENERG"); - m_inx_theta = m_aeff.axis("THETA"); - m_inx_aeff = m_aeff.table("EFFAREA"); - - // Get optional index (use "EFFAREA" if "EFFAREA_RECO" does not exist) - m_inx_aeff_reco = (m_aeff.has_table("EFFAREA_RECO")) ? - m_aeff.table("EFFAREA_RECO") : m_inx_aeff; - // Throw an exception if the table is not two-dimensional if (m_aeff.axes() != 2) { std::string msg = "Expected two-dimensional effective area response " @@ -299,6 +301,9 @@ void GCTAAeff2D::read(const GFitsTable& table) throw GException::invalid_value(G_READ, msg); } + // Set table indices + set_indices(); + // Set energy axis to logarithmic scale m_aeff.axis_log10(m_inx_energy); @@ -310,6 +315,9 @@ void GCTAAeff2D::read(const GFitsTable& table) m_aeff.scale(i, 1.0e4); } + // Set table boundaries + set_boundaries(); + // Return return; } @@ -347,10 +355,8 @@ void GCTAAeff2D::write(GFitsBinTable& table) const * * @param[in] filename FITS file name. * - * Loads the effective area from a FITS file. - * - * If no extension name is provided, the effective area will be loaded from - * the "EFFECTIVE AREA" extension. + * Loads the effective area from a FITS file. If no extension name is given + * the effective area will be loaded from the "EFFECTIVE AREA" extension. ***************************************************************************/ void GCTAAeff2D::load(const GFilename& filename) { @@ -375,15 +381,15 @@ void GCTAAeff2D::load(const GFilename& filename) /***********************************************************************//** - * @brief Save effectiva area into FITS file + * @brief Save effective area into FITS file * * @param[in] filename FITS file name. - * @param[in] clobber Overwrite existing file? (default: false) + * @param[in] clobber Overwrite existing file? * * Saves effectiva area into a FITS file. If a file with the given * @p filename does not yet exist it will be created, otherwise the method * opens the existing file. The method will create a (or replace an existing) - * effectiva area extension. The extension name can be specified as part + * effective area extension. The extension name can be specified as part * of the @p filename, or if no extension name is given, is assumed to be * "EFFECTIVE AREA". * @@ -430,7 +436,7 @@ void GCTAAeff2D::save(const GFilename& filename, const bool& clobber) const * @param[in] logE Log10 of the true photon energy (TeV). * @param[in] zenith Zenith angle in Earth system (rad). * @param[in] azimuth Azimuth angle in Earth system (rad). - * @param[in] etrue Use true energy (default: true). + * @param[in] etrue Use true energy? * @return Maximum effective area (cm2). * * Returns the maximum effective area for a given energy, zenith and azimuth @@ -453,9 +459,9 @@ double GCTAAeff2D::max(const double& logE, // Loop over theta values for (int i = 0; i < n_theta; ++i) { - // Compute lower and upper theta bin values - double theta_lo = m_aeff.axis_lo(m_inx_theta, i); - double theta_hi = m_aeff.axis_hi(m_inx_theta, i); + // Compute lower and upper theta bin values (radians) + double theta_lo = m_aeff.axis_lo(m_inx_theta, i) * gammalib::deg2rad; + double theta_hi = m_aeff.axis_hi(m_inx_theta, i) * gammalib::deg2rad; // Get effective area value in cm2 double aeff_lo = (m_inx_energy == 0) ? m_aeff(index, logE, theta_lo) : @@ -480,9 +486,39 @@ double GCTAAeff2D::max(const double& logE, /***********************************************************************//** + * @brief Assign response table + * + * @param[in] table Response table. + * + * Assigns the response table for an effective area. The effective area + * values are given in units of cm2. + ***************************************************************************/ +void GCTAAeff2D::table(const GCTAResponseTable& table) +{ + // Assign response table + m_aeff = table; + + // Set indices + set_indices(); + + // Set energy axis to logarithmic scale + m_aeff.axis_log10(m_inx_energy); + + // Set offset angle axis to radians + m_aeff.axis_radians(m_inx_theta); + + // Set table boundaries + set_boundaries(); + + // Return + return; +} + + +/***********************************************************************//** * @brief Print effective area information * - * @param[in] chatter Chattiness (defaults to NORMAL). + * @param[in] chatter Chattiness. * @return String containing effective area information. ***************************************************************************/ std::string GCTAAeff2D::print(const GChatter& chatter) const @@ -542,6 +578,10 @@ void GCTAAeff2D::init_members(void) m_inx_theta = 1; m_inx_aeff = 0; m_inx_aeff_reco = 1; + m_logE_min = 0.0; + m_logE_max = 0.0; + m_theta_min = 0.0; + m_theta_max = 0.0; // Return return; @@ -562,6 +602,10 @@ void GCTAAeff2D::copy_members(const GCTAAeff2D& aeff) m_inx_theta = aeff.m_inx_theta; m_inx_aeff = aeff.m_inx_aeff; m_inx_aeff_reco = aeff.m_inx_aeff_reco; + m_logE_min = aeff.m_logE_min; + m_logE_max = aeff.m_logE_max; + m_theta_min = aeff.m_theta_min; + m_theta_max = aeff.m_theta_max; // Return return; @@ -576,3 +620,50 @@ void GCTAAeff2D::free_members(void) // Return return; } + + +/***********************************************************************//** + * @brief Set table indices + * + * Sets the data members m_inx_energy, m_inx_theta, m_inx_aeff and optionally + * m_inx_aeff_reco. + ***************************************************************************/ +void GCTAAeff2D::set_indices(void) +{ + // Get mandatory indices (throw exception if not found) + m_inx_energy = m_aeff.axis("ENERG"); + m_inx_theta = m_aeff.axis("THETA"); + m_inx_aeff = m_aeff.table("EFFAREA"); + + // Get optional index (use "EFFAREA" if "EFFAREA_RECO" does not exist) + m_inx_aeff_reco = (m_aeff.has_table("EFFAREA_RECO")) ? + m_aeff.table("EFFAREA_RECO") : m_inx_aeff; + + // Return + return; +} + + +/***********************************************************************//** + * @brief Set effective area boundaries + * + * Sets the data members m_logE_min, m_logE_max, m_theta_min and m_theta_max + * that define the validity range of the effective area. + ***************************************************************************/ +void GCTAAeff2D::set_boundaries(void) +{ + // Compute minimum and maximum logE boundaries + m_logE_min = std::log10(m_aeff.axis_lo(m_inx_aeff, 0)); + m_logE_max = std::log10(m_aeff.axis_hi(m_inx_aeff, + m_aeff.axis_bins(m_inx_aeff)-1)); + + // Compute minimum and maximum theta boundaries + m_theta_min = m_aeff.axis_lo(m_inx_theta, 0) * + gammalib::deg2rad; + m_theta_max = m_aeff.axis_hi(m_inx_theta, m_aeff.axis_bins(m_inx_theta)-1) * + gammalib::deg2rad; + + // Return + return; +} + diff --git a/inst/cta/test/test_CTA.py b/inst/cta/test/test_CTA.py index 442c2b4..2cfc969 100755 --- a/inst/cta/test/test_CTA.py +++ b/inst/cta/test/test_CTA.py @@ -82,10 +82,28 @@ class Test(gammalib.GPythonTestSuite): aeff = gammalib.GCTAAeff2D(filename) # Test Aeff values - self.test_value(aeff(0.0, 0.0), 5535774176.75, 0.1) - self.test_value(aeff(1.0, 0.0), 20732069462.7, 0.1) - self.test_value(aeff(0.0, 0.01745), 5682897797.76, 0.1) - self.test_value(aeff(1.0, 0.01745), 18446656815.1, 0.1) + self.test_value(aeff(0.0, 0.0), 5535774176.75, 0.1, + 'Test reference effective area value') + self.test_value(aeff(1.0, 0.0), 20732069462.7, 0.1, + 'Test reference effective area value') + self.test_value(aeff(0.0, 0.01745), 5682897797.76, 0.1, + 'Test reference effective area value') + self.test_value(aeff(1.0, 0.01745), 18446656815.1, 0.1, + 'Test reference effective area value') + + # Test that Aeff values outside boundaries are zero + self.test_value(aeff(-1.80001, 0.0), 0.0, 1.0e-6, + 'Test that effective area is zero for energy below' + ' minimum energy') + self.test_value(aeff(+2.20001, 0.0), 0.0, 1.0e-6, + 'Test that effective area is zero for energy above' + ' maximum energy') + self.test_value(aeff(0.0, -0.00001), 0.0, 1.0e-6, + 'Test that effective area is zero for offset angle' + ' below minimum offset angle') + self.test_value(aeff(0.0, 0.13963), 0.0, 1.0e-6, + 'Test that effective area is zero for offset angle' + ' above maximum offset angle') # Test GCTAAeffPerfTable file constructor filename = self.caldb + '/cta_dummy_irf.dat' -- libgit2 0.21.4