diff --git a/ChangeLog b/ChangeLog index 119f6b5..f1ac4a9 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,11 @@ +2020-12-02 + + * Bug fix version 1.7.3 released + ================================ + + Fix inaccuracy in CTA background model evaluation (#3472) + + 2020-10-15 * Bug fix version 1.7.2 released diff --git a/NEWS b/NEWS index 4ccde08..27ab2e9 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ -New Features and Important Changes in GammaLib 1.7.2 +New Features and Important Changes in GammaLib 1.7.3 -15 October 2020 +2 December 2020 1. Introduction @@ -63,11 +63,7 @@ None 8. Model module --------------- -The copy constructors of the GModelSpatialDiffuse, GModelSpatialElliptical -and GModelSpatialRadial classes did not properly copy the base class members -in their copy constructors. This prevented specifically the fitting of -diffuse maps in the stacked analysis. The code was modified to fix the -issue. +None 9. Numerics module @@ -127,7 +123,12 @@ None 20. CTA interface ----------------- -None +The GCTAModelIrfBackground::eval() method uses now the GCTABackground3D::rate_ebin() +integrator method instead of the GCTABackground3D::operator() operator so that the +background model evaluation is no longer sensitive to the number of energy bins. +The GCTACubeBackground::fill() method was simplified, removing the power law +integration which was not accurate enough, and using directly the +GCTAModelIrfBackground::eval() method which is now exact (#3472). 21. COMPTEL interface diff --git a/README.md b/README.md index ef39f61..f4b3550 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,6 @@ GammaLib information ==================== -* Version: 1.7.2 (15 October 2020) +* Version: 1.7.3 (2 December 2020) [![Build Status](https://cta-jenkins.irap.omp.eu/buildStatus/icon?job=gammalib-integrate-os)](https://cta-jenkins.irap.omp.eu/job/gammalib-integrate-os/) diff --git a/doc/source/admin/release_history/1.7.3.rst b/doc/source/admin/release_history/1.7.3.rst new file mode 100644 index 0000000..30a9b2b --- /dev/null +++ b/doc/source/admin/release_history/1.7.3.rst @@ -0,0 +1,14 @@ +.. _1.7.3: + +GammaLib 1.7.3 release (2 December 2020) +======================================== + +Overview +-------- + +GammaLib 1.7.3 is a bug fix release for GammaLib release 1.7. + +The following bugs were fixed: + +* [`3472 `_] - + Fix inaccuracy in CTA background model evaluation diff --git a/doc/source/admin/release_history/index.rst b/doc/source/admin/release_history/index.rst index 56c767e..fcab8cd 100644 --- a/doc/source/admin/release_history/index.rst +++ b/doc/source/admin/release_history/index.rst @@ -6,6 +6,7 @@ Release History .. toctree:: :maxdepth: 1 + GammaLib 1.7.3 <1.7.3> GammaLib 1.7.2 <1.7.2> GammaLib 1.7.1 <1.7.1> GammaLib 1.7 <1.7> diff --git a/inst/cta/include/GCTAEventBin.hpp b/inst/cta/include/GCTAEventBin.hpp index b334fa2..b1cbeb4 100644 --- a/inst/cta/include/GCTAEventBin.hpp +++ b/inst/cta/include/GCTAEventBin.hpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAEventBin.hpp - CTA event bin class * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2017 by Juergen Knoedlseder * + * copyright (C) 2010-2020 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -80,6 +80,8 @@ public: const int& ipix(void) const; const int& ieng(void) const; const double& solidangle(void) const; + GEnergy emin(void) const; + GEnergy emax(void) const; const GEnergy& ewidth(void) const; const double& ontime(void) const; const double& weight(void) const; @@ -185,4 +187,18 @@ void GCTAEventBin::ieng(const int& ieng) return; } + +/***********************************************************************//** + * @brief Return maximum energy of event bin + * + * @return Maximum energy of event bin. + * + * Returns the maximum energy of event bin. + ***************************************************************************/ +inline +GEnergy GCTAEventBin::emax(void) const +{ + return (emin() + ewidth()); +} + #endif /* GCTAEVENTBIN_HPP */ diff --git a/inst/cta/pyext/GCTAEventBin.i b/inst/cta/pyext/GCTAEventBin.i index b577e19..2b6e217 100644 --- a/inst/cta/pyext/GCTAEventBin.i +++ b/inst/cta/pyext/GCTAEventBin.i @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAEventBin.i - CTA event bin class * * ----------------------------------------------------------------------- * - * copyright (C) 2011-2018 by Juergen Knoedlseder * + * copyright (C) 2011-2020 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -61,6 +61,8 @@ public: const int& ipix(void) const; const int& ieng(void) const; const double& solidangle(void) const; + GEnergy emin(void) const; + GEnergy emax(void) const; const GEnergy& ewidth(void) const; const double& ontime(void) const; const double& weight(void) const; diff --git a/inst/cta/src/GCTACubeBackground.cpp b/inst/cta/src/GCTACubeBackground.cpp index 96938bc..7bee4ff 100644 --- a/inst/cta/src/GCTACubeBackground.cpp +++ b/inst/cta/src/GCTACubeBackground.cpp @@ -466,44 +466,23 @@ void GCTACubeBackground::fill(const GObservations& obs, GLog* log) // DETX/DETY coordinates GCTAInstDir dir = cta->pointing().instdir(bin->dir().dir()); - // Initialise values for integration over energy bin - double f1 = 0.0; - // Loop over all energy bins of the cube for (int iebin = 0, ibin = ipix; iebin < nebins; ++iebin, ibin += npix) { - // If this is the first bin then compute f1 at the lower bin edge - if (iebin == 0) { - GCTAEventBin* bin = eventcube[ibin]; - bin->dir(dir); // Set instrument direction - bin->energy(m_ebounds.emin(iebin)); - f1 = obs.models().eval(*bin, *cta); - } - - // Compute f2 at the upper bin edge - GCTAEventBin* bin = eventcube[ibin]; - bin->dir(dir); // Set instrument direction - bin->energy(m_ebounds.emax(iebin)); - double f2 = obs.models().eval(*bin, *cta); - - // Add background model value to cube if the energy bin is - // fully container in the observation interval and if both f1 - // and f2 are positive. Add a small margin here to the energy - // bin boundaries to take provision for rounding errors. + // Compute background model value if the energy bin is fully + // contained in the observation interval. Add a small margin + // here to the energy bin boundaries to take provision for + // rounding errors. if (obs_ebounds.contains(m_ebounds.emin(iebin) + energy_margin, - m_ebounds.emax(iebin) - energy_margin) && - (f1 > 0.0 && f2 > 0.0)) { + m_ebounds.emax(iebin) - energy_margin)) { - // Compute background model value by integrating over the - // bin and deviding the result by the bin width - double x1 = m_ebounds.emin(iebin).MeV(); - double x2 = m_ebounds.emax(iebin).MeV(); - double model = gammalib::plaw_integral(x1, f1, x2, f2) / (x2 - x1); + // Set energy bin + GCTAEventBin* bin = eventcube[ibin]; + bin->dir(dir); // Set instrument direction + bin->energy(m_ebounds.elogmean(iebin)); - // Multiply by ontime to get the correct weighting for each - // observation. We divide by the total ontime later to get - // the background model in units of counts/(MeV s sr). - model *= ontime; + // Compute model + double model = obs.models().eval(*bin, *cta) * ontime; // Add existing number of counts model += bin->counts(); @@ -513,10 +492,7 @@ void GCTACubeBackground::fill(const GObservations& obs, GLog* log) } // endif: added background model value to cube - // Set f1 to f2 for next energy layer - f1 = f2; - - } // endfor: looped over energy layers + } // endfor: looped over energy bins } // endfor: looped over spatial pixels diff --git a/inst/cta/src/GCTAEventBin.cpp b/inst/cta/src/GCTAEventBin.cpp index 7b6298d..81c3f78 100644 --- a/inst/cta/src/GCTAEventBin.cpp +++ b/inst/cta/src/GCTAEventBin.cpp @@ -1,7 +1,7 @@ /*************************************************************************** * GCTAEventBin.cpp - CTA event bin class * * ----------------------------------------------------------------------- * - * copyright (C) 2010-2016 by Juergen Knoedlseder * + * copyright (C) 2010-2020 by Juergen Knoedlseder * * ----------------------------------------------------------------------- * * * * This program is free software: you can redistribute it and/or modify * @@ -342,6 +342,40 @@ const double& GCTAEventBin::solidangle(void) const /***********************************************************************//** + * @brief Return minimum energy of event bin + * + * @return Minimum energy of event bin + * + * Returns minimum energy of event bin, computed using + * + * \f[ + * E_{\rm min} = \frac{-\Delta E + sqrt{\Delta E^2 + 4 E^2}}{2} + * \f] + * + * where + * \f$\Delta E\f$ is the energy bin width, returned by ewidth() and + * \f$E\f$ is the energy, returned by energy(). + ***************************************************************************/ +GEnergy GCTAEventBin::emin(void) const +{ + // Get energy and energy width in MeV + double energy = this->energy().MeV(); + double ewidth = this->ewidth().MeV(); + + // Compute minimum energy in MeV + double arg = ewidth * ewidth + 4.0 * energy * energy; + double e_min = 0.5 * (-ewidth + std::sqrt(arg)); + + // Set minimum energy + GEnergy emin; + emin.MeV(e_min); + + // Return minimum energy + return emin; +} + + +/***********************************************************************//** * @brief Return energy width of event bin * * @return Energy width of event bin diff --git a/inst/cta/src/GCTAModelIrfBackground.cpp b/inst/cta/src/GCTAModelIrfBackground.cpp index 2a0e10d..250bca6 100644 --- a/inst/cta/src/GCTAModelIrfBackground.cpp +++ b/inst/cta/src/GCTAModelIrfBackground.cpp @@ -40,6 +40,7 @@ #include "GCTAObservation.hpp" #include "GCTAResponseIrf.hpp" #include "GCTABackground.hpp" +#include "GCTAEventBin.hpp" #include "GCTASupport.hpp" /* __ Globals ____________________________________________________________ */ @@ -58,9 +59,10 @@ const GModelRegistry g_cta_inst_background_registry(&g_cta_inst_backgrou /* __ Coding definitions _________________________________________________ */ #define G_USE_NPRED_CACHE +#define G_USE_RATE_EBIN /* __ Debug definitions __________________________________________________ */ -//#define G_DUMP_MC +#define G_DUMP_MC //#define G_DEBUG_NPRED /* __ Constants __________________________________________________________ */ @@ -349,10 +351,27 @@ double GCTAModelIrfBackground::eval(const GEvent& event, // Set DETX and DETY in instrument direction GCTAInstDir inst_dir = pnt.instdir(dir.dir()); - - // Evaluate function + + // Evaluate spatial component + #if defined(G_USE_RATE_EBIN) + double spat = 0.0; + if (event.is_atom()) { + double logE = event.energy().log10TeV(); + spat = bgd(logE, inst_dir.detx(), inst_dir.dety()); + } + else { + const GCTAEventBin* bin = static_cast(&event); + GEnergy emin = bin->emin(); + GEnergy emax = emin + bin->ewidth(); + double norm = 1.0 / bin->ewidth().MeV(); + spat = bgd.rate_ebin(inst_dir, emin, emax) * norm; + } + #else double logE = event.energy().log10TeV(); double spat = bgd(logE, inst_dir.detx(), inst_dir.dety()); + #endif + + // Evaluate function double spec = (spectral() != NULL) ? spectral()->eval(event.energy(), event.time(), gradients) : 1.0; diff --git a/inst/cta/test/test_CTA.cpp b/inst/cta/test/test_CTA.cpp index 32fc81d..f11978d 100644 --- a/inst/cta/test/test_CTA.cpp +++ b/inst/cta/test/test_CTA.cpp @@ -2229,6 +2229,8 @@ void TestGCTAObservation::test_event_bin(void) test_value(bin.ieng(), 7, "Check ieng() for filled bin"); test_value(bin.solidangle(), 3.14, 1.0e-10, "Check solidangle() for filled bin"); + test_value(bin.emin().TeV(), 1.996503, "Check emin() for filled bin"); + test_value(bin.emax().TeV(), 2.003503, "Check emax() for filled bin"); test_value(bin.ewidth().GeV(), 7.0, 1.0e-10, "Check ewidth() for filled bin"); test_value(bin.ontime(), 101.0, 1.0e-10, "Check ontime() for filled bin");