From f02bf3f306d209454f202b4579c5073388bb0823 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Thu, 8 Nov 2018 14:40:50 +0100 Subject: [PATCH] Revise GCTAEdispCube::operator() interface (#2721) --- ChangeLog | 3 ++- NEWS | 8 ++++++-- doc/source/admin/release_history/1.6.rst | 2 ++ inst/cta/include/GCTACubeEdisp.hpp | 10 +++++----- inst/cta/pyext/GCTACubeEdisp.i | 6 +++--- inst/cta/src/GCTACubeEdisp.cpp | 59 ++++++++++++++++++++++++++++++++++++----------------------- inst/cta/src/GCTAResponseCube.cpp | 65 ++++++++++++++++++++--------------------------------------------- inst/cta/test/test_CTA.cpp | 80 +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++--------- inst/cta/test/test_CTA.hpp | 3 +++ 9 files changed, 148 insertions(+), 88 deletions(-) diff --git a/ChangeLog b/ChangeLog index eb2e5ff..5ab9cbd 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -2018-11-06 +2018-11-08 * Version 1.6.0 released ======================== + Revise GCTAEdispCube::operator() interface (#2721) Rename GCTAResponseIrf::lo_save_thres() and GCTAResponseIrf::hi_save_thres() to GCTAResponseIrf::lo_safe_thres() and GCTAResponseIrf::hi_safe_thres() (#2716) diff --git a/NEWS b/NEWS index 0ed8274..0763d6d 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in GammaLib 1.6.0 -6 November 2018 +8 November 2018 1. Introduction @@ -151,12 +151,12 @@ The arguments for the following methods have been changed: - GCTAModelRadial::mc(GRan&) - GCTAEdisp::operator(GEnergy&, GEnergy&, double&, double&, double&, double&) - GCTAEdisp::mc(GRan&, GEnergy&, double&, double&, double&, double&) +- GCTACubeEdisp::operator()(GEnergy&, GEnergy&, GSkyDir&) - GCTAResponseIrf::edisp(GEnergy&, GEnergy&, double&, double&, double&, double&) - GLATAeff(GFits&, std::string&) - GLATEdisp(GFits&, std::string&) - GLATPsf(GFits&, std::string&) - The return value of the following methods has been changed: - none @@ -314,6 +314,10 @@ None 20. CTA interface ----------------- +The arguments of the GCTACubeEdisp::operator() were changed from, passing +now reconstructed and true energy objects, followed by a true sky direction. +This is more consistent with the GCTAEdisp interface (#2721). + Rename GCTAResponseIrf::lo_save_thres() and GCTAResponseIrf::hi_save_thres() to GCTAResponseIrf::lo_safe_thres() and GCTAResponseIrf::hi_safe_thres() (#2716). diff --git a/doc/source/admin/release_history/1.6.rst b/doc/source/admin/release_history/1.6.rst index 1d776e8..bf6e3f2 100644 --- a/doc/source/admin/release_history/1.6.rst +++ b/doc/source/admin/release_history/1.6.rst @@ -47,6 +47,8 @@ Bug fixes Improvements ------------ +* [`2721 `_] - + Revise ``GCTAEdispCube::operator()`` interface * [`2715 `_] - Add generic CTA response cache * [`2704 `_] - diff --git a/inst/cta/include/GCTACubeEdisp.hpp b/inst/cta/include/GCTACubeEdisp.hpp index f41fa43..221523e 100644 --- a/inst/cta/include/GCTACubeEdisp.hpp +++ b/inst/cta/include/GCTACubeEdisp.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTACubeEdisp.hpp - CTA cube analysis energy dispersion class * * ----------------------------------------------------------------------- * - * copyright (C) 2016-2017 by Michael Mayer * + * copyright (C) 2016-2018 by Michael Mayer * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -84,9 +84,9 @@ public: // Operators GCTACubeEdisp& operator=(const GCTACubeEdisp& cube); - double operator()(const GSkyDir& dir, - const double& migra, - const GEnergy& energy) const; + double operator()(const GEnergy& ereco, + const GEnergy& etrue, + const GSkyDir& dir) const; // Methods void clear(void); @@ -116,7 +116,7 @@ protected: void clear_cube(void); void fill_cube(const GCTAObservation& obs, GSkyMap* exposure = NULL, GLog* log = NULL); - void update(const double& migra, const double& logEsrc) const; + void update(const GEnergy& ereco, const GEnergy& etrue) const; void set_eng_axis(void); void compute_ebounds(void) const; diff --git a/inst/cta/pyext/GCTACubeEdisp.i b/inst/cta/pyext/GCTACubeEdisp.i index e524b52..f13d9b8 100644 --- a/inst/cta/pyext/GCTACubeEdisp.i +++ b/inst/cta/pyext/GCTACubeEdisp.i @@ -64,9 +64,9 @@ public: virtual ~GCTACubeEdisp(void); // Interpolation operator - double operator()(const GSkyDir& dir, - const double& migra, - const GEnergy& energy) const; + double operator()(const GEnergy& ereco, + const GEnergy& etrue, + const GSkyDir& dir) const; // Methods void clear(void); diff --git a/inst/cta/src/GCTACubeEdisp.cpp b/inst/cta/src/GCTACubeEdisp.cpp index 75026d7..18fff38 100644 --- a/inst/cta/src/GCTACubeEdisp.cpp +++ b/inst/cta/src/GCTACubeEdisp.cpp @@ -323,23 +323,23 @@ GCTACubeEdisp& GCTACubeEdisp::operator=(const GCTACubeEdisp& cube) /***********************************************************************//** - * @brief Return energy dispersion (in units of MeV^-1) + * @brief Return energy dispersion in units of MeV\f$^{-1}\f$ * + * @param[in] ereco Reconstructed event energy. + * @param[in] etrue True photon energy. * @param[in] dir Coordinate of the true photon position. - * @param[in] migra Energy migration (reconstructed over true energy). - * @param[in] srcEng True photon energy. - * @return Energy dispersion (in units of MeV^-1) + * @return Energy dispersion (MeV\f$^{-1}\f$) * * Returns the energy dispersion for a given energy migration (reconstructed * over true energy), true photon energy, and sky direction in units of - * MeV^-1. + * MeV\f$^{-1}\f$. ***************************************************************************/ -double GCTACubeEdisp::operator()(const GSkyDir& dir, - const double& migra, - const GEnergy& srcEng) const +double GCTACubeEdisp::operator()(const GEnergy& ereco, + const GEnergy& etrue, + const GSkyDir& dir) const { // Update indices and weighting factors for interpolation - update(migra, srcEng.log10TeV()); + update(ereco, etrue); // Perform bi-linear interpolation double edisp = m_wgt1 * m_cube(dir, m_inx1) + @@ -837,12 +837,21 @@ void GCTACubeEdisp::clear_cube(void) /***********************************************************************//** * @brief Fill energy dispersion cube from observation container * - * @param[in] obs Observation. + * @param[in] obs CTA observation. * @param[in] exposure Pointer towards exposure map. * @param[in] log Pointer towards logger. * * @exception GException::invalid_value * No RoI or response found in CTA observation. + * + * Fills the energy dispersion cube from the observation container. + * + * If the @p exposure argument is not NULL, the energy dispersion is weighted + * by the exposure, and the exposure weighting is added to the exposure map. + * + * If the @p log argument is not NULL, the method puts information about + * exclusion and inclusion of a CTA observation into the energy dispersion + * cube. ***************************************************************************/ void GCTACubeEdisp::fill_cube(const GCTAObservation& obs, GSkyMap* exposure, @@ -965,12 +974,14 @@ void GCTACubeEdisp::fill_cube(const GCTAObservation& obs, // Get log10 of true energy in TeV double logEsrc = m_energies[iebin].log10TeV(); - // Compute exposure weight - double weight = rsp->aeff(theta, 0.0, 0.0, 0.0, logEsrc) * - obs.livetime(); - - // If available, accumulate weights + // Compute exposure weight. If no exposure map is + // specified, the weight is one. Otherwise, the weight + // is the effective area for this offset angle and + // source energy times the livetime of the observation. + double weight = 1.0; if (exposure != NULL) { + weight = rsp->aeff(theta, 0.0, 0.0, 0.0, logEsrc) * + obs.livetime(); (*exposure)(pixel, iebin) += weight; } @@ -985,9 +996,8 @@ void GCTACubeEdisp::fill_cube(const GCTAObservation& obs, continue; } - // Compute log10 of reconstructed energy in TeV - double logEobs = std::log10(migra) + logEsrc; - GEnergy Eobs(std::pow(10.0, logEobs), "TeV"); + // Compute reconstructed energy + GEnergy Eobs = migra * m_energies[iebin]; // Set map index int imap = offset(imigra, iebin); @@ -1016,19 +1026,22 @@ void GCTACubeEdisp::fill_cube(const GCTAObservation& obs, /***********************************************************************//** * @brief Update energy dispersion cube indices and weights cache * - * @param[in] migra Reconstructed over true energy. - * @param[in] logEsrc Log10 of true photon energy in TeV. + * @param[in] ereco Reconstructed event energy. + * @param[in] etrue True photon energy. * * Updates the energy dispersion cube indices, stored in the members m_inx1, * m_inx2, m_inx3, and m_inx4, and weights cache, stored in m_wgt1, m_wgt2, * m_wgt3, and m_wgt4. ***************************************************************************/ -void GCTACubeEdisp::update(const double& migra, const double& logEsrc) const +void GCTACubeEdisp::update(const GEnergy& ereco, const GEnergy& etrue) const { // Set node array for energy interpolation - m_elogmeans.set_value(logEsrc); + m_elogmeans.set_value(etrue.log10TeV()); + + // Compute migration value + double migra = ereco/etrue; - // Set node array for delta interpolation + // Set node array for migra interpolation m_migras.set_value(migra); // Set indices for bi-linear interpolation diff --git a/inst/cta/src/GCTAResponseCube.cpp b/inst/cta/src/GCTAResponseCube.cpp index 365c14a..4187849 100644 --- a/inst/cta/src/GCTAResponseCube.cpp +++ b/inst/cta/src/GCTAResponseCube.cpp @@ -354,14 +354,8 @@ double GCTAResponseCube::irf(const GEvent& event, // Multiply-in energy dispersion if (use_edisp() && irf > 0.0) { - - // Determine fraction of reconstructed and true energy - double migra = obsEng.TeV() / srcEng.TeV(); - - // Multiply-in energy dispersion - irf *= edisp()(srcDir, migra, srcEng); - - } // endif: energy dispersion was available and psf was non-zero + irf *= edisp()(obsEng, srcEng, srcDir); + } // Divide by livetime irf /= livetime; @@ -1135,9 +1129,6 @@ double GCTAResponseCube::irf_ptsrc(const GEvent& event, // direction in radians double delta = bin->dir().dir().dist(srcDir); - // Determine fraction of reconstructed and true energy - double migra = bin->energy().TeV() / srcEng.TeV(); - // Get maximum angular separation for PSF (in radians) double delta_max = psf().delta_max(); @@ -1162,11 +1153,8 @@ double GCTAResponseCube::irf_ptsrc(const GEvent& event, // Multiply-in energy dispersion if (use_edisp() && irf > 0.0) { - - // Multiply-in energy dispersion - irf *= edisp()(srcDir, migra, srcEng); - - } // endif: energy dispersion was available and psf was non-zero + irf *= edisp()(bin->energy(), srcEng, srcDir); + } // Apply deadtime correction irf *= exposure().deadc(); @@ -1236,9 +1224,6 @@ double GCTAResponseCube::irf_radial(const GEvent& event, // (radians) double rho_obs = model->dir().dist(obsDir); - // Determine fraction of reconstructed and true energy - double migra = bin->energy().TeV() / srcEng.TeV(); - // Get livetime (in seconds) double livetime = exposure().livetime(); @@ -1264,17 +1249,14 @@ double GCTAResponseCube::irf_radial(const GEvent& event, irf *= psf_radial(model, rho_obs, obsDir, srcEng, obsTime); // Multiply-in energy dispersion + // + // The current code assumes that the energy dispersion at the + // observed and true event direction does not vary + // significantly. In other words, the code assumes that the + // energy dispersion is constant over the size of the PSF. if (use_edisp() && irf > 0.0) { - - // Multiply-in energy dispersion - // - // The current code assumes that the energy dispersion at the - // observed and true event direction does not vary - // significantly. In other words, the code assumes that the - // energy dispersion is constant over the size of the PSF. - irf *= edisp()(obsDir, migra, srcEng); - - } // endif: energy dispersion was available and psf was non-zero + irf *= edisp()(bin->energy(), srcEng, obsDir); + } // Apply deadtime correction irf *= exposure().deadc(); @@ -1344,9 +1326,6 @@ double GCTAResponseCube::irf_elliptical(const GEvent& event, double rho_obs = model->dir().dist(obsDir); double posangle_obs = model->dir().posang(obsDir); - // Determine fraction of reconstructed and true energy - double migra = bin->energy().TeV() / srcEng.TeV(); - // Get livetime (in seconds) double livetime = exposure().livetime(); @@ -1372,23 +1351,20 @@ double GCTAResponseCube::irf_elliptical(const GEvent& event, irf *= psf_elliptical(model, rho_obs, posangle_obs, obsDir, srcEng, obsTime); // Multiply-in energy dispersion + // + // The current code assumes that the energy dispersion at the + // observed and true event direction does not vary + // significantly. In other words, the code assumes that the + // energy dispersion is constant over the size of the PSF. if (use_edisp() && irf > 0.0) { - - // Multiply-in energy dispersion - // - // The current code assumes that the energy dispersion at the - // observed and true event direction does not vary - // significantly. In other words, the code assumes that the - // energy dispersion is constant over the size of the PSF. - irf *= edisp()(obsDir, migra, srcEng); - - } // endif: energy dispersion was available and psf was non-zero + irf *= edisp()(bin->energy(), srcEng, obsDir); + } // Apply deadtime correction irf *= exposure().deadc(); } // endif: exposure was positive - + } // endif: we were sufficiently close and livetime >0 // Compile option: Check for NaN/Inf @@ -1471,8 +1447,7 @@ double GCTAResponseCube::irf_diffuse(const GEvent& event, if (use_edisp()) { irf = cache->irf(bin->ipix(), source.energy()); if (irf > 0.0) { - irf *= edisp()(bin->dir().dir(), bin->energy()/source.energy(), - source.energy()); + irf *= edisp()(bin->energy(), source.energy(), bin->dir().dir()); } } diff --git a/inst/cta/test/test_CTA.cpp b/inst/cta/test/test_CTA.cpp index f26ee9e..28329e0 100644 --- a/inst/cta/test/test_CTA.cpp +++ b/inst/cta/test/test_CTA.cpp @@ -1364,8 +1364,8 @@ void TestGCTAResponse::test_response_bkgcube(void) ***************************************************************************/ void TestGCTAResponse::test_response_edispcube(void) { - // Test PSF cube constructors - test_try("CTA Edisp cube void constructor"); + // Test energy dispersion cube constructors + test_try("CTA energy dispersion cube void constructor"); try { GCTACubeEdisp cube; test_try_success(); @@ -1373,7 +1373,7 @@ void TestGCTAResponse::test_response_edispcube(void) catch (std::exception &e) { test_try_failure(e); } - test_try("CTA Edisp cube map constructor"); + test_try("CTA energy dispersion cube map constructor"); try { GEnergies energies(20, GEnergy(0.1, "TeV"), GEnergy(100.0, "TeV")); GCTACubeEdisp cube("CAR", "CEL", 83.63, 22.01, 0.4, 0.4, 10, 10, @@ -1384,24 +1384,41 @@ void TestGCTAResponse::test_response_edispcube(void) test_try_failure(e); } - // Test set method + // Load IRF for energy dispersion computation GCTAObservation obs_cta; obs_cta.load(cta_events); obs_cta.response(cta_irf, GCaldb(cta_caldb)); - GEnergies energies(20, GEnergy(0.1, "TeV"), GEnergy(100.0, "TeV")); + + // Initialise energy dispersion cube + GEnergies energies(100, GEnergy(0.1, "TeV"), GEnergy(100.0, "TeV")); GCTACubeEdisp cube("CAR", "CEL", 83.63, 22.01, 0.4, 0.4, 10, 10, - energies, 2.0, 20); + energies, 3.0, 100); + + // Test setting of observation cube.set(obs_cta); - cube.save("test_cta_edispcube_one.fits", true); - // Test fill method + // Test energy dispersion integration + test_edispcube_integration(cube, 0.1, 10.0); + + // Test filling of single observation GObservations obs; obs_cta.id("000001"); obs.append(obs_cta); + cube.fill(obs); + + // Test energy dispersion integration + test_edispcube_integration(cube, 0.1, 10.0); + + // Test fill of two observations obs_cta.id("000002"); obs.append(obs_cta); cube.fill(obs); - cube.save("test_cta_edispcube_two.fits", true); + + // Test energy dispersion integration + test_edispcube_integration(cube, 0.1, 10.0); + + // Test saving + cube.save("test_cta_edispcube.fits", true); // Return return; @@ -1688,6 +1705,51 @@ void TestGCTAResponse::test_edisp_integration(const GCTAEdisp& edisp, /***********************************************************************//** + * @brief Utility function for energy dispersion cube integration tests + * + * @param[in] edisp Energy dispersion cube instance. + * @param[in] e_src_min Minimum true energy (TeV). + * @param[in] e_src_max Maximum true energy (TeV). + * + * Tests whether the energy dispersion cube integrated over observed + * energies gives unity. + ***************************************************************************/ +void TestGCTAResponse::test_edispcube_integration(const GCTACubeEdisp& edisp, + const double& e_src_min, + const double& e_src_max) +{ + // Set sky direction + GSkyDir dir; + dir.radec_deg(83.63, 22.01); + + // Test energy dispersion cube integration + for (double e_src = e_src_min; e_src <= e_src_max; e_src *= 2.0) { + + // Compute log10 of true energy + GEnergy etrue(e_src, "TeV"); + + // Perform numerical integration by summing + const int nE = 10000; + double emin = 1000.0; //!< 1 GeV + double emax = 3.0 * etrue.MeV(); //!< Migra_max = 3.0 + double dE = (emax-emin)/nE; + double sum = 0.0; + double E_obs = emin; + for (int i = 0; i < nE; ++i, E_obs += dE) { + GEnergy ereco(E_obs,"MeV"); + double dp = edisp(ereco, etrue, dir); + sum += dp * dE; + } + test_value(sum, 1.0, 0.001, "Response cube energy dispersion integration " + " for "+etrue.print()); + } + + // Return + return; +} + + +/***********************************************************************//** * @brief Test CTA background ***************************************************************************/ void TestGCTAModel::test_model_bgd(void) diff --git a/inst/cta/test/test_CTA.hpp b/inst/cta/test/test_CTA.hpp index 46caff6..0204417 100644 --- a/inst/cta/test/test_CTA.hpp +++ b/inst/cta/test/test_CTA.hpp @@ -99,6 +99,9 @@ public: void test_edisp_integration(const GCTAEdisp& edisp, const double& e_src_min = 0.1, const double& e_src_max = 10.0); + void test_edispcube_integration(const GCTACubeEdisp& edisp, + const double& e_src_min = 0.1, + const double& e_src_max = 10.0); }; -- libgit2 0.21.4