From 03c28a611cbd1c2bde895480c9eba616bdaad294 Mon Sep 17 00:00:00 2001 From: Florent Forest Date: Thu, 26 Feb 2015 10:14:40 +0100 Subject: [PATCH] Corrected GCTAEdisp2D methods Tested Edisp2D and EdispPerfTable methods --- inst/cta/include/GCTAEdisp2D.hpp | 9 --------- inst/cta/src/GCTAEdisp2D.cpp | 181 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------------------------------------------------------------------------------------- inst/cta/test/test_CTA.cpp | 143 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++---------------------------- 3 files changed, 208 insertions(+), 125 deletions(-) diff --git a/inst/cta/include/GCTAEdisp2D.hpp b/inst/cta/include/GCTAEdisp2D.hpp index b2209ed..034e04c 100644 --- a/inst/cta/include/GCTAEdisp2D.hpp +++ b/inst/cta/include/GCTAEdisp2D.hpp @@ -100,15 +100,6 @@ private: // Members std::string m_filename; //!< Name of Edisp response file GCTAResponseTable m_edisp; //!< Edisp response table - - // Precomputation cache - mutable double m_par_logEobs; - mutable double m_par_logEsrc; - mutable double m_par_theta; - mutable double m_norm; - mutable double m_EobsOnEsrc; //!< linear ratio Eobs/Esrc - mutable double m_eres; - mutable double m_ebias; }; /***********************************************************************//** diff --git a/inst/cta/src/GCTAEdisp2D.cpp b/inst/cta/src/GCTAEdisp2D.cpp index 7001e4a..1fc5489 100644 --- a/inst/cta/src/GCTAEdisp2D.cpp +++ b/inst/cta/src/GCTAEdisp2D.cpp @@ -178,7 +178,8 @@ double GCTAEdisp2D::operator()(const double& logEobs, double edisp = 0.0; // Compute edisp - edisp = m_edisp(0, logEsrc, Eobs / Esrc, theta); + edisp = m_edisp(0, logEsrc, std::exp((logEobs-logEsrc)*std::log(10.0)), theta); + // Return return edisp; @@ -335,19 +336,20 @@ void GCTAEdisp2D::read(const GFits& fits) // Normalize vectors of migration matrix for all etrue and theta for (int i_etrue = 0; i_etrue < etrue_size; ++i_etrue) { - double EobsOnEtrue = 0.5 * (m_edisp.axis_hi(1, i_etrue) + m_edisp.axis_lo(1, i_etrue)); for (int i_theta = 0; i_theta < theta_size; ++i_theta) { // Compute sum double sum = 0.0; for (int i_migra = 0; i_migra < migra_size; ++i_migra) { - // Compute delta(Eobs/Etrue) + // Compute delta(Eobs/Etrue) and Eobs/Etrue double delta = m_edisp.axis_hi(1, i_migra) - m_edisp.axis_lo(1, i_migra); + double EobsOnEtrue = 0.5 * (m_edisp.axis_hi(1, i_migra) + m_edisp.axis_lo(1, i_migra)); + // Add dispersion to sum sum += m_edisp(0, i_etrue + (i_migra + i_theta*migra_size)*etrue_size) * delta / EobsOnEtrue; } if (sum > 0.0) { for (int i_migra = 0; i_migra < migra_size; ++i_migra) { - // Normalize by sum - m_edisp(0, i_etrue + (i_migra + i_theta*migra_size)*etrue_size) /= sum; + // Normalize by sum/log(10.0) + m_edisp(0, i_etrue + (i_migra + i_theta*migra_size)*etrue_size) /= sum / std::log(10.0); } } } @@ -355,9 +357,6 @@ void GCTAEdisp2D::read(const GFits& fits) // Set true energy axis to logarithmic scale m_edisp.axis_log10(0); - - // Set migration energy ratio axis to logarithmic scale - m_edisp.axis_log10(1); // Set offset angle axis to radians m_edisp.axis_radians(2); @@ -389,37 +388,39 @@ GEnergy GCTAEdisp2D::mc(GRan& ran, // Vector representing cumulative probability distribution std::vector > cumul; - const double deltaLogE = 0.01; - double logEobs = 0.01; + // Initialize cumulative probability double sum = 0.0; - // Loop through Eobs - for(int i = 0; i < 3.0 / deltaLogE; ++i) { - // Increment Eobs - logEobs += deltaLogE; + + // Loop through EobsOnEtrue + for (int i = 0; i < m_edisp.axis(1); ++i) { + + // Compute delta(Eobs/Etrue) and Eobs/Etrue values + double delta = m_edisp.axis_hi(1, i) - m_edisp.axis_lo(1, i); + double EobsOnEtrue = 0.5 * (m_edisp.axis_hi(1, i) + m_edisp.axis_lo(1, i)); // Compute cumulated probability - sum += GCTAEdisp2D::operator()(logEobs, logEsrc, theta) * deltaLogE; + sum += m_edisp(0, logEsrc, EobsOnEtrue, theta) * delta; - // Create pair containing Eobs/Esrc and cumulated probability - std::pair pair(std::exp((logEobs - logEsrc) * std::log(10.0)), sum); + // Create pair containing EobsOnEtrue and cumulated probability + std::pair pair(EobsOnEtrue, sum); // Add to vector cumul.push_back(pair); } - // Draw random number between 0 and 1 + // Draw random number between 0 and 1 from uniform distribution double p = ran.uniform(); // Find right index int index = 0; - while(index < 3.0 / deltaLogE - 1.0 && cumul[index+1].second < p) { + while(index < m_edisp.axis(1) - 1 && cumul[index+1].second < p) { index++; } // index found - // Compute result Eobs/Esrc value + // Compute result EobsOnEtrue value double result = (cumul[index+1].second - p)*cumul[index].first + (p - cumul[index].second)*cumul[index+1].first; - result *= (cumul[index+1].second - cumul[index].second); + result /= (cumul[index+1].second - cumul[index].second); // Final result @@ -438,15 +439,14 @@ GEnergy GCTAEdisp2D::mc(GRan& ran, * @brief Return observed energy interval that contains the energy dispersion. * * @param[in] logEsrc Log10 of the true photon energy (TeV). - * @param[in] theta Offset angle in camera system (rad). Not used. + * @param[in] theta Offset angle in camera system (rad). * @param[in] phi Azimuth angle in camera system (rad). Not used. * @param[in] zenith Zenith angle in Earth system (rad). Not used. * @param[in] azimuth Azimuth angle in Earth system (rad). Not used. * * Returns the band of observed energies outside of which the energy - * dispersion becomes negligible for a given true energy @p logEsrc. This - * band is set to \f$\pm 5 \times \sigma\f$, where \f$\sigma\f$ is the - * Gaussian width of the energy dispersion. + * dispersion becomes negligible for a given true energy @p logEsrc and offset + * angle @p theta. An energy in considered negligible if inferior to eps. ***************************************************************************/ GEbounds GCTAEdisp2D::ebounds_obs(const double& logEsrc, const double& theta, @@ -454,17 +454,43 @@ GEbounds GCTAEdisp2D::ebounds_obs(const double& logEsrc, const double& zenith, const double& azimuth) const { - // Set energy band constant - const double number_of_sigmas = 5.0; + // Set epsilon + const double eps = 1.0e-06; - // Update the parameter cache - update(m_par_logEobs, logEsrc, theta); + double logEobsMin = -10.0; + double logEobsMax = 30.0; + bool minFound = false; + bool maxFound = false; + + for (int i = 0; i < m_edisp.axis(1); ++i) { + // Compute value EobsOnEtrue + double EobsOnEtrue = 0.5 * (m_edisp.axis_hi(1, i) + m_edisp.axis_lo(1, i)); + + // Find first non-negligible matrix term + if (!minFound && m_edisp(0, logEsrc, EobsOnEtrue, theta) >= eps) { + minFound = true; + logEobsMin = std::log10(EobsOnEtrue) + logEsrc; + } + // Find last non-negligible matrix term + else if (minFound && !maxFound && m_edisp(0, logEsrc, EobsOnEtrue, theta) < eps) { + maxFound = true; + logEobsMax = std::log10(EobsOnEtrue) + logEsrc; + } + else if (minFound && maxFound && m_edisp(0, logEsrc, EobsOnEtrue, theta) >= eps) { + maxFound = false; + } + } + + // If energy dispersion has never become negligible until end of loop, reset logEobsMax + if (!maxFound) { + logEobsMax = 30.0; + } // Compute energy boundaries GEnergy emin; GEnergy emax; - emin.log10TeV(logEsrc - number_of_sigmas * m_eres); - emax.log10TeV(logEsrc + number_of_sigmas * m_eres); + emin.log10TeV(logEobsMin); + emax.log10TeV(logEobsMax); // Return energy boundaries return (GEbounds(emin, emax)); @@ -475,15 +501,14 @@ GEbounds GCTAEdisp2D::ebounds_obs(const double& logEsrc, * @brief Return true energy interval that contains the energy dispersion. * * @param[in] logEobs Log10 of the observed event energy (TeV). - * @param[in] theta Offset angle in camera system (rad). Not used. + * @param[in] theta Offset angle in camera system (rad). * @param[in] phi Azimuth angle in camera system (rad). Not used. * @param[in] zenith Zenith angle in Earth system (rad). Not used. * @param[in] azimuth Azimuth angle in Earth system (rad). Not used. * * Returns the band of true photon energies outside of which the energy - * dispersion becomes negligible for a given observed energy @p logEobs. This - * band is set to \f$\pm 5 \times \sigma\f$, where \f$\sigma\f$ is the - * Gaussian width of the energy dispersion. + * dispersion becomes negligible for a given observed energy @p logEobs and + * offset angle @p theta. An energy in considered negligible if inferior to eps. ***************************************************************************/ GEbounds GCTAEdisp2D::ebounds_src(const double& logEobs, const double& theta, @@ -491,17 +516,43 @@ GEbounds GCTAEdisp2D::ebounds_src(const double& logEobs, const double& zenith, const double& azimuth) const { - // Set energy band constant - const double number_of_sigmas = 5.0; - - // Update the parameter cache - update(logEobs, m_par_logEsrc, theta); + // Set epsilon + const double eps = 1.0e-06; + double logEsrcMin = -10.0; + double logEsrcMax = 30.0; + bool minFound = false; + bool maxFound = false; + + for (int i = 0; i < m_edisp.axis(0); ++i) { + // Compute value of corresponding logEsrc + double logEsrc = 0.5 * (m_edisp.axis_hi(0, i) + m_edisp.axis_lo(0, i)); + double EobsOnEtrue = std::exp((logEobs-logEsrc)*std::log(10.0)); + // Find first non-negligible matrix term + if (!minFound && m_edisp(0, logEsrc, EobsOnEtrue, theta) >= eps) { + minFound = true; + logEsrcMin = logEsrc; + } + // Find last non-negligible matrix term + else if (minFound && !maxFound && m_edisp(0, logEsrc, EobsOnEtrue, theta) < eps) { + maxFound = true; + logEsrcMax = logEsrc; + } + else if (minFound && maxFound && m_edisp(0, logEsrc, EobsOnEtrue, theta) >= eps) { + maxFound = false; + } + } + + // If energy dispersion has never become negligible until end of loop, reset logEobsMax + if (!maxFound) { + logEsrcMax = 30.0; + } + // Compute energy boundaries GEnergy emin; GEnergy emax; - emin.log10TeV(logEobs - number_of_sigmas * m_eres); - emax.log10TeV(logEobs + number_of_sigmas * m_eres); + emin.log10TeV(logEsrcMin); + emax.log10TeV(logEsrcMax); // Return energy boundaries return (GEbounds(emin, emax)); @@ -564,13 +615,6 @@ void GCTAEdisp2D::init_members(void) // Initialise members m_filename.clear(); m_edisp.clear(); - m_par_logEobs = -1.0e30; - m_par_logEsrc = -1.0e30; - m_par_theta = -1.0; - m_norm = 1.0; - m_EobsOnEsrc = 0.0; - m_eres = 0.0; - m_ebias = 0.0; // Return return; @@ -601,42 +645,3 @@ void GCTAEdisp2D::free_members(void) // Return return; } - -/***********************************************************************//** - * @brief Update Energy dispersion parameter cache - * - * @param[in] logEobs Log10 of the observed photon energy (TeV). - * @param[in] logEsrc Log10 of the true photon energy (TeV). - * @param[in] theta Offset angle in camera system (rad). - * - * This method updates the Energy dispersion parameter cache. - ***************************************************************************/ -void GCTAEdisp2D::update(const double& logEobs, const double& logEsrc, const double& theta) const -{ - /* - // Only compute Edisp parameters if arguments have changed - if (logEobs != m_par_logEobs || logEsrc != m_par_logEsrc || theta != m_par_theta) { - - // Save parameters - m_par_logEobs = logEobs; - m_par_logEsrc = logEsrc; - m_par_theta = theta; - - // Compute linear Eobs/Esrc ratio - m_EobsOnEsrc = std::exp( (logEobs - logEsrc) * std::log(10)); - - - - - // Get ERES and EBIAS values from response table - m_eres = m_edisp(0, logEsrc, theta); - m_ebias = m_edisp(1, logEsrc, theta); - - // Compute global normalization parameter - m_norm = (m_eres > 0.0) ? gammalib::inv_sqrt2pi / m_eres : 0.0; - - } - */ - // Return - return; -} diff --git a/inst/cta/test/test_CTA.cpp b/inst/cta/test/test_CTA.cpp index b9bef7c..8fa0849 100644 --- a/inst/cta/test/test_CTA.cpp +++ b/inst/cta/test/test_CTA.cpp @@ -37,7 +37,9 @@ #include "GNodeArray.hpp" #include "test_CTA.hpp" #include "GCTAEdisp2D.hpp" +#include "GCTAEdispPerfTable.hpp" #include "GCTAResponseTable.hpp" +#include "GMath.hpp" /* __ Namespaces _________________________________________________________ */ @@ -590,45 +592,60 @@ void TestGCTAResponse::test_response_edisp2D(void) catch (std::exception &e) { test_try_failure(e); } - + // Test normalization + test_try("GCTAEdisp2D normalization"); try { - GCTAEdisp2D edisp2D_2(cta_irf_matrix); - - GCTAResponseTable m_edisp = edisp2D_2.table(); - int etrue_size = m_edisp.axis(0); - int migra_size = m_edisp.axis(1); - int theta_size = m_edisp.axis(2); - - for (int i_etrue = 0; i_etrue < etrue_size; ++i_etrue) { - double EobsOnEtrue = 0.5 * (m_edisp.axis_hi(1, i_etrue) + m_edisp.axis_lo(1, i_etrue)); - for (int i_theta = 0; i_theta < theta_size; ++i_theta) { - // Compute sum - double sum = 0.0; - for (int i_migra = 0; i_migra < migra_size; ++i_migra) { - // Compute delta(Eobs/Etrue) - double delta = m_edisp.axis_hi(1, i_migra) - m_edisp.axis_lo(1, i_migra); - sum += m_edisp(0, i_etrue + (i_migra + i_theta*migra_size)*etrue_size) * delta / EobsOnEtrue; + GCTAEdisp2D edisp2D_2(cta_irf_matrix); + //std::cout << "TEST NORM 2D" << std::endl; + + const double deltaLogE = 0.01; + for(double logEsrc = -3.0; logEsrc < 4.0; logEsrc += deltaLogE) { + double sum = 0.0; + for(double logEobs = -3.0; logEobs < 7.0; logEobs += deltaLogE) { + sum += edisp2D_2(logEobs, logEsrc) * deltaLogE; } - std::cout << sum << std::endl; - } + //std::cout << sum << std::endl; + } + + test_try_success(); + } + catch (std::exception &e) { + test_try_failure(e); } + // Test ebounds_obs and ebounds_src methods + test_try("GCTAEdisp2D ebounds_src and ebounds_obs methods"); + try { + GCTAEdisp2D edisp2D_3(cta_irf_matrix); + //std::cout << "TEST EBOUNDS 2D" << std::endl; + // ebounds_obs (for a given logEsrc) + //std::cout << edisp2D_3.ebounds_obs(-1.0).print(EXPLICIT) << std::endl; + // ebounds_src (for a given logEobs) + //std::cout << edisp2D_3.ebounds_src(2.0).print(EXPLICIT) << std::endl; + + test_try_success(); + } + catch (std::exception &e) { + test_try_failure(e); + } - const double deltaLogE = 0.01; - for(double logEsrc = -1.0; logEsrc < 2.0; logEsrc += deltaLogE) { - double sum = 0.0; - for(double logEobs = -1.0; logEobs < 2.0; logEobs += deltaLogE) { - sum += edisp2D_2(logEobs, logEsrc) * deltaLogE; + // Test mc method + test_try("GCTAEdisp2D mc method"); + try { + GCTAEdisp2D edisp2D_4(cta_irf_matrix); + GRan ran; + //std::cout << "TEST MC 2D" << std::endl; + for (int i = 0; i < 20; ++i) { + //std::cout << edisp2D_4.mc(ran, -1.0).GeV() << std::endl; } - std::cout << sum << std::endl; - } - + + test_try_success(); } catch (std::exception &e) { test_try_failure(e); } - + // Return return; } @@ -638,7 +655,77 @@ void TestGCTAResponse::test_response_edisp2D(void) ***************************************************************************/ void TestGCTAResponse::test_response_edispPerfTable(void) { + /// /// TO DO + + // Load response + GCTAResponseIrf rsp; + rsp.caldb(GCaldb(cta_caldb)); + // Test Energy Dispersion + test_energy_integration(rsp); + + + // Test GCTAEdispPerfTable constructor + test_try("GCTAEdispPerfTable constructor"); + try { + GCTAEdispPerfTable test(cta_caldb+"/"+cta_irf+".dat"); + test_try_success(); + } + catch (std::exception &e) { + test_try_failure(e); + } + + GCTAEdispPerfTable test2(cta_caldb+"/"+cta_irf+".dat"); + + // Test ebounds_obs and ebounds_src methods + test_try("GCTAEdispPerfTable ebounds_obs and ebounds_src methods"); + try { + //std::cout << "TEST EBOUNDS PERFTABLE" << std::endl; + // ebounds_obs (for a given logEsrc) + //std::cout << test2.ebounds_obs(-1.0).print(EXPLICIT) << std::endl; + // ebounds_src (for a given logEobs) + //std::cout << test2.ebounds_src(2.0).print(EXPLICIT) << std::endl; + test_try_success(); + } + catch (std::exception &e) { + test_try_failure(e); + } + + // Test mc method + test_try("GCTAEdispPerfTable mc method"); + try { + GRan ran; + //std::cout << "TEST MC PERFTABLE" << std::endl; + for (int i = 0; i < 20; ++i) { + //std::cout << test2.mc(ran, -1.0).GeV() << std::endl; + } + test_try_success(); + } + catch (std::exception &e) { + test_try_failure(e); + } + + // Test normalization + test_try("GCTAEdispPerfTable normalization"); + try { + //std::cout << "TEST NORM PERFTABLE" << std::endl; + const double deltaLogE = 0.01; + for(double logEsrc = -3.0; logEsrc < 4.0; logEsrc += deltaLogE) { + double sum = 0.0; + for(double logEobs = -3.0; logEobs < 7.0; logEobs += deltaLogE) { + sum += test2(logEobs, logEsrc) * deltaLogE; + } + //std::cout << sum << std::endl; + } + + test_try_success(); + } + catch (std::exception &e) { + test_try_failure(e); + } + + + // Return return; } -- libgit2 0.21.4