Correct MC simulations for normalized phase curves (#2172)

Jürgen Knödlseder
1 parent 45ed9e48
Showing 3 changed files with 14 additions and 8 deletions
ChangeLog
include/GModelTemporalPhaseCurve.hpp
src/model/GModelTemporalPhaseCurve.cpp
-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)
@@ -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<double> m_values;   //!< Function values at nodes
-    GFilename           m_filename; //!< Name of phase file function
+    mutable GNodeArray  m_nodes;     //!< Phase values of nodes
+    std::vector<double> m_values;    //!< Function values at nodes
+    double              m_scale;     //!< Scale factor due to normalization
+    GFilename           m_filename;  //!< Name of phase file function
 };
@@ -316,7 +316,7 @@ GTimes GModelTemporalPhaseCurve::mc(const double&amp; rate, const GTime&amp;  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&amp; rate, const GTime&amp;  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&amp; 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