From 5a80b6d4c8ff4c4e940a78d47f14306ac128ebc3 Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Thu, 1 Nov 2018 20:42:32 +0100 Subject: [PATCH] Fix precision issue with GCTABackground3D Monte-Carlo simulations (#2708) --- ChangeLog | 3 ++- NEWS | 8 +++++++- doc/source/admin/release_history/1.6.rst | 2 ++ inst/cta/src/GCTABackground3D.cpp | 62 ++++++++++++++++++++++++++++++++++++++++++++------------------ 4 files changed, 55 insertions(+), 20 deletions(-) diff --git a/ChangeLog b/ChangeLog index 1bb3dcc..760a2b3 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,8 +1,9 @@ -2018-10-31 +2018-11-01 * Version 1.6.0 released ======================== + Fix precision issue with GCTABackground3D Monte-Carlo simulations (#2708) Revise GCTAEdisp class interface (#2704) GCTAEdisp returns now energy dispersion values per MeV (#2704) Protect GCTAInstDir from access to invalid information (#2694) diff --git a/NEWS b/NEWS index 960f196..9dac06d 100644 --- a/NEWS +++ b/NEWS @@ -1,6 +1,6 @@ New Features and Important Changes in GammaLib 1.6.0 -31 October 2018 +1 November 2018 1. Introduction @@ -306,6 +306,12 @@ None 20. CTA interface ----------------- +GCTABackground3D::init_mc_cache now pre-computes the Monte Carlo cache at energy +nodes of the background cube, and in addition, inserts 8 logarithmically spaced +nodes between each pair of nodes. This assures that the energy variations of the +background rate are precisely captured and reproduced in the Monte Carlo +simulations (#2708). + The GCTAEdisp::operator() now returns energy dispersion values per MeV. The units are now consistent with the expectation in GResponse, and this changes assurance the coherence between the various classes. In particular, no re- diff --git a/doc/source/admin/release_history/1.6.rst b/doc/source/admin/release_history/1.6.rst index 77e785d..3d342a6 100644 --- a/doc/source/admin/release_history/1.6.rst +++ b/doc/source/admin/release_history/1.6.rst @@ -19,6 +19,8 @@ In particular, this release provides: Bug fixes --------- +* [`2708 `_] - + Fix precision issue with ``GCTABackground3D`` Monte-Carlo simulations * [`2683 `_] - Correct handling of stacked On/Off observations without background IRF. Fix bug that potentially wrote the wrong exposure time in a PHA file. diff --git a/inst/cta/src/GCTABackground3D.cpp b/inst/cta/src/GCTABackground3D.cpp index 5a6a938..55cc961 100644 --- a/inst/cta/src/GCTABackground3D.cpp +++ b/inst/cta/src/GCTABackground3D.cpp @@ -912,6 +912,9 @@ int GCTABackground3D::index(const int& idetx, ***************************************************************************/ void GCTABackground3D::init_mc_cache(void) const { + // Set number of subbins per energy bin of response cube + const int nsubbins = 10; + // Initialise cache m_mc_max.clear(); m_mc_spectrum.clear(); @@ -936,23 +939,50 @@ void GCTABackground3D::init_mc_cache(void) const } m_mc_one_minus_costheta = 1.0 - std::cos(theta); - // Set number of energy bins. Make sure that we compute at least 100 energy - // bins - int nbins = m_num_energy * 10; - if (nbins < 100) { - nbins = 100; - } - - // Compute DETX and DETY binsize in radians and energy binsize in log10(TeV) + // Compute DETX and DETY binsize in radians double detx_bin = (m_detx_max - m_detx_min) / m_num_detx; double dety_bin = (m_dety_max - m_dety_min) / m_num_dety; - double logE_bin = (m_logE_max - m_logE_min) / (nbins-1); + + // Set energy bins for Monte Carlo cache spectrum. Make sure that the + // actual energy nodes are part of the energy bins to avoid that sharp + // structures in the background rate variations are missed. + GEnergies energies; + for (int i = 0; i < m_num_energy; ++i) { + + // If this is the first node then add energies below the first node + // down to the lower energy limit of the background response + if (i == 0) { + GEnergy emin(m_background.axis_lo(m_inx_energy, i), + m_background.axis_lo_unit(m_inx_energy)); + GEnergies enodes(nsubbins/2, emin, m_energy[i]); + enodes.remove(enodes.size()-1); // avoids duplication of energy + energies.extend(enodes); + } + + // If this is not the last node then append nodes between current + // node and next node + if (i < m_num_energy-1) { + GEnergies enodes(nsubbins, m_energy[i], m_energy[i+1]); + enodes.remove(enodes.size()-1); // avoids duplication of energy + energies.extend(enodes); + } + + // ... otherwise this is the last node. Add energies above the last + // node up to the upper energy limit of the background response + else { + GEnergy emax(m_background.axis_hi(m_inx_energy, i), + m_background.axis_hi_unit(m_inx_energy)); + GEnergies enodes(nsubbins/2, m_energy[i], emax); + energies.extend(enodes); + } + + } // endfor: looped over energy bins of background response // Loop over all energy bins - for (int i = 0; i < nbins; ++i) { + for (int i = 0; i < energies.size(); ++i) { // Get logE of energy bin - double logE = m_logE_min + i * logE_bin; + double logE = energies[i].log10TeV(); // Initialise cache with cumulative pixel fluxes and compute total flux // in response table for normalization. Negative pixels are excluded @@ -1007,18 +1037,14 @@ void GCTABackground3D::init_mc_cache(void) const } // endfor: loop over DETX - // Set energy value (unit independent) - GEnergy energy; - energy.log10(logE, m_background.axis_lo_unit(m_inx_energy)); - - // Only append node if flux is positive + // Append node if flux is positive if (total_flux > 0.0) { - m_mc_spectrum.append(energy, total_flux); + m_mc_spectrum.append(energies[i], total_flux); } // Dump spectrum for debugging #if defined(G_DEBUG_MC_INIT) - std::cout << "Energy=" << energy; + std::cout << "Energy=" << energies[i]; std::cout << " Rate_total=" << total_flux; std::cout << " events/s/MeV" << std::endl; #endif -- libgit2 0.21.4