From b33b0f65c74a9749628ce3d692abcae29ed28c1e Mon Sep 17 00:00:00 2001 From: Jürgen Knödlseder Date: Fri, 28 Jul 2017 10:11:08 +0200 Subject: [PATCH] Correct MC simulations for normalized phase curves (#2172) --- ChangeLog | 4 ++-- include/GModelTemporalPhaseCurve.hpp | 7 ++++--- src/model/GModelTemporalPhaseCurve.cpp | 11 ++++++++--- 3 files changed, 14 insertions(+), 8 deletions(-) diff --git a/ChangeLog b/ChangeLog index 1ba180d..9d7203c 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,11 +1,11 @@ -2017-07-27 +2017-07-28 * Version 1.4.0 released ========================= Add GCTAEventBin::ipix() and ieng() setter methods (#2171) Write response information into simulated observation definition file (#2170) - Normalise phase curves by default (#2154) + Normalise phase curves by default (#2154, #2172) Support Python list slicing for container classes (#2128) Add GSkyRegions::overlaps() and GSkyMap::overlaps() classes (#1980) diff --git a/include/GModelTemporalPhaseCurve.hpp b/include/GModelTemporalPhaseCurve.hpp index 7f28230..e43ee6c 100644 --- a/include/GModelTemporalPhaseCurve.hpp +++ b/include/GModelTemporalPhaseCurve.hpp @@ -144,9 +144,10 @@ protected: bool m_has_normalize; //!< XML has normalize attribute // Protected members - mutable GNodeArray m_nodes; //!< Phase values of nodes - std::vector m_values; //!< Function values at nodes - GFilename m_filename; //!< Name of phase file function + mutable GNodeArray m_nodes; //!< Phase values of nodes + std::vector m_values; //!< Function values at nodes + double m_scale; //!< Scale factor due to normalization + GFilename m_filename; //!< Name of phase file function }; diff --git a/src/model/GModelTemporalPhaseCurve.cpp b/src/model/GModelTemporalPhaseCurve.cpp index 77990aa..eb05452 100644 --- a/src/model/GModelTemporalPhaseCurve.cpp +++ b/src/model/GModelTemporalPhaseCurve.cpp @@ -316,7 +316,7 @@ GTimes GModelTemporalPhaseCurve::mc(const double& rate, const GTime& tmin, GTimes times; // Compute event rate (in events per seconds) - double lambda = rate * norm(); + double lambda = rate * norm() * m_scale; // Initialise start and stop times in seconds double time = tmin.secs(); @@ -341,7 +341,7 @@ GTimes GModelTemporalPhaseCurve::mc(const double& rate, const GTime& tmin, double value = eval(srcTime); // Get uniform random number - double uniform = ran.uniform(); + double uniform = ran.uniform() * m_scale; // If the random number is not larger than the phase curve value // then accept the time @@ -686,6 +686,7 @@ void GModelTemporalPhaseCurve::init_members(void) m_nodes.clear(); m_values.clear(); m_filename.clear(); + m_scale = 1.0; m_normalize = true; m_has_normalize = false; @@ -711,6 +712,7 @@ void GModelTemporalPhaseCurve::copy_members(const GModelTemporalPhaseCurve& mode m_nodes = model.m_nodes; m_values = model.m_values; m_filename = model.m_filename; + m_scale = model.m_scale; m_normalize = model.m_normalize; m_has_normalize = model.m_has_normalize; @@ -909,9 +911,12 @@ void GModelTemporalPhaseCurve::normalize_nodes(void) throw GException::invalid_value(G_NORMALISE_NODES, msg); } + // Set scale factor + m_scale = 1.0 / sum; + // Normalise the node values for (int i = 0; i < m_values.size(); ++i) { - m_values[i] /= sum; + m_values[i] *= m_scale; } // Return -- libgit2 0.21.4