Added inv_ln10, bisection searches for MC, computation of cumulated probability for each true energy

Florent Forest
1 parent cdd1386a
Showing 7 changed files with 171 additions and 107 deletions
include/GMath.hpp
inst/cta/include/GCTAEdisp2D.hpp
inst/cta/include/GCTAEdispRmf.hpp
inst/cta/src/GCTAEdisp2D.cpp
inst/cta/src/GCTAEdispPerfTable.cpp
inst/cta/src/GCTAEdispRmf.cpp
inst/cta/test/example_sim_edisp.py
@@ -51,6 +51,7 @@ namespace gammalib {
     const double sqrt_onehalf = std::sqrt(1.0/2.0);
     const double sqrt_pihalf  = std::sqrt(pihalf);
     const double sqrt_two     = std::sqrt(2.0);
+    const double inv_ln10     =  1.0/ln10;
 }
  
 /* __ Prototypes _________________________________________________________ */
@@ -103,10 +103,18 @@ private:
     void update(const double& logEobs,
                 const double& logEsrc,
                 const double& theta) const;
+    void update_cumul(const double& logEsrc,
+                      const double& theta = 0.0,
+                      const double& phi = 0.0,
+                      const double& zenith = 0.0,
+                      const double& azimuth = 0.0) const;
  
     // Members
     std::string       m_filename;  //!< Name of Edisp response file
     GCTAResponseTable m_edisp;     //!< Edisp response table
+    mutable double                                  m_logEsrc;
+    mutable double                                  m_theta;
+    mutable std::vector<std::pair<double, double> > m_cumul;
 };
  
  
@@ -100,11 +100,10 @@ private:
     void set_cache(void) const;
     void set_mc_cache(void) const;
     void update(const double& arg1, const double& arg2) const;
-    void update_cumul(const double& logEsrc,
-                      const double& theta = 0.0,
-                      const double& phi = 0.0,
-                      const double& zenith = 0.0,
-                      const double& azimuth = 0.0) const;
+    void compute_cumul(const double& theta = 0.0,
+                       const double& phi = 0.0,
+                       const double& zenith = 0.0,
+                       const double& azimuth = 0.0) const;
  
     // Members
     std::string   m_filename;  //!< Name of response file
@@ -128,9 +127,11 @@ private:
     // Monte Carlo cache
     mutable std::vector<int>                        m_mc_measured_start;
     mutable std::vector<GVector>                    m_mc_measured_cdf;
-    mutable double                                  m_logEsrc;
+    mutable bool                                    m_cdf_computed;
     mutable double                                  m_theta;
-    mutable std::vector<std::pair<double, double> > m_cumul;
+
+    mutable std::vector<std::vector<std::pair<double, double> > > m_cumul;
+
 };
  
  
@@ -77,6 +77,7 @@ GCTAEdisp2D::GCTAEdisp2D(void) : GCTAEdisp()
  ***************************************************************************/
 GCTAEdisp2D::GCTAEdisp2D(const std::string& filename) : GCTAEdisp()
 {
+
     // Initialise class members
     init_members();
  
@@ -182,8 +183,6 @@ double GCTAEdisp2D::operator()(const double&amp; logEobs,
     double EobsOverEsrc = std::exp((logEobs-logEsrc) * gammalib::ln10);
  
     // Compute edisp
-    //edisp = m_edisp(0, logEsrc, std::exp((logEobs-logEsrc)*std::log(10.0)),
-    //                theta);
     edisp = m_edisp(0, logEsrc, EobsOverEsrc, theta);
  
     // Return
@@ -237,6 +236,7 @@ GCTAEdisp2D* GCTAEdisp2D::clone(void) const
  ***************************************************************************/
 void GCTAEdisp2D::load(const std::string& filename)
 {
+
     // Open FITS file
     GFits fits(filename);
  
@@ -269,6 +269,7 @@ void GCTAEdisp2D::load(const std::string&amp; filename)
  ***************************************************************************/
 void GCTAEdisp2D::read(const GFits& fits)
 {
+
     // Clear response table
     m_edisp.clear();
  
@@ -310,10 +311,10 @@ void GCTAEdisp2D::read(const GFits&amp; fits)
             if (sum != 0.0) {
                 for (int i_migra = 0, i = offset; i_migra < migra_size;
                      ++i_migra, i += etrue_size) {
-                    m_edisp(0, i) /= (sum / gammalib::ln10);
+                    m_edisp(0, i) *= (gammalib::ln10 / sum);
                 }
             }
-            
+
         } // endfor: looped over theta
     } // endfor: looped over Etrue
  
@@ -391,66 +392,32 @@ GEnergy GCTAEdisp2D::mc(GRan&amp;         ran,
                         const double& zenith,
                         const double& azimuth) const
 {
-    // Vector representing cumulative probability distribution
-    std::vector<std::pair<double, double> > cumul;
-
-    // Initialize cumulative probability
-    double sum = 0.0;
-
-    // Loop over EobsOverEsrc
-    // TODO: We should precompute the cumulative distributions for all
-    //       tabluated Esrc values and put them in a Monte Carlo cache.
-    //       It's probably sufficient to select the nearest vector and
-    //       use it to draw energies. Otherwise this method will be
-    //       very slow as it computes the individual vectors for each
-    //       call.
-    for (int i = 0; i < m_edisp.axis(1); ++i) {
-
-        // Stop if we have reached the maximum
-        if (sum == 1.0) {
-            break;
-        }
-
-        // Compute delta(Eobs/Etrue) and Eobs/Etrue values
-        double delta        = m_edisp.axis_hi(1, i) -
-                              m_edisp.axis_lo(1, i);
-        double EobsOverEsrc = 0.5 * (m_edisp.axis_hi(1, i) +
-                                     m_edisp.axis_lo(1, i));
-
-        // Compute cumulative probability
-        double add = m_edisp(0, logEsrc, EobsOverEsrc, theta) *
-                     delta / EobsOverEsrc / gammalib::ln10;
-
-        // Limit cumulative probability to 1
-        sum += add;
-        if (sum > 1.0) {
-            sum = 1.0;
-        }
  
-        // Create pair containing EobsOverEsrc and cumulated probability
-        std::pair<double, double> pair(EobsOverEsrc, sum);
-
-        // Add to vector
-        cumul.push_back(pair);
-
-    } // endfor: looped over EobsOverEsrc
+    // Update cumul
+    update_cumul(logEsrc, theta, phi, zenith, azimuth);
  
     // Draw random number between 0 and 1 from uniform distribution
     double p = ran.uniform();
  
-    // Find right index
-    // TODO: This can lead to index = m_edisp.axis(1) on axit, which
-    //       leads to an array indexing problem in the interpolation.
-    //       Also, a bisection search would be more efficient.
-    int index = 0;
-    while(index < m_edisp.axis(1) - 2 && cumul[index+1].second < p) {
-        index++;
-    } // index found
+    // Find right index with bisection search
+    int low = 0;
+    int high = m_cumul.size();
+    while ((high-low) > 1) {
+        int mid = (low+high) / 2;
+        if (p < m_cumul[mid].second) {
+            high = mid;
+        }
+        else if (m_cumul[mid].second <= p) {
+            low = mid;
+        }
+    }
+    // Index found
+    int index = low;
  
     // Interpolate EobsOverEsrc value
-    double result =   (cumul[index+1].second - p)*cumul[index].first
-                    + (p - cumul[index].second)*cumul[index+1].first;
-    result       /=   (cumul[index+1].second - cumul[index].second);
+    double result =   (m_cumul[index+1].second - p)*m_cumul[index].first
+                    + (p - m_cumul[index].second)*m_cumul[index+1].first;
+    result       /=   (m_cumul[index+1].second - m_cumul[index].second);
  
     // Final result
     double logEnergy = std::log10(result) + logEsrc;
@@ -520,7 +487,7 @@ GEbounds GCTAEdisp2D::ebounds_obs(const double&amp; logEsrc,
         }
  
     } // endfor: looped over EobsOverEsrc
-    
+
     // If energy dispersion has never become negligible until end of loop,
     // reset logEobsMax
     if (!maxFound) {
@@ -566,7 +533,7 @@ GEbounds GCTAEdisp2D::ebounds_src(const double&amp; logEobs,
     double logEsrcMax =  30.0;
     bool   minFound   = false;
     bool   maxFound   = false;
-    
+
     // Find boundaries
     for (int i = 0; i < m_edisp.axis(0); ++i) {
  
@@ -580,7 +547,7 @@ GEbounds GCTAEdisp2D::ebounds_src(const double&amp; logEobs,
             minFound   = true;
             logEsrcMin = logEsrc;
         }
-        
+
         // Find last non-negligible matrix term
         else if (minFound && !maxFound &&
                  m_edisp(0, logEsrc, EobsOverEsrc, theta) < eps) {
@@ -593,15 +560,15 @@ GEbounds GCTAEdisp2D::ebounds_src(const double&amp; logEobs,
                  m_edisp(0, logEsrc, EobsOverEsrc, theta) >= eps) {
             maxFound = false;
         }
-        
+
     } // endfor: looped over true energy
-    
+
     // If energy dispersion has never become negligible until end of loop,
     // reset logEobsMax
     if (!maxFound) {
         logEsrcMax = 30.0;
     }
-    
+
     // Compute energy boundaries
     GEnergy emin;
     GEnergy emax;
@@ -671,6 +638,11 @@ void GCTAEdisp2D::init_members(void)
     m_filename.clear();
     m_edisp.clear();
  
+    // Initialise Monte Carlo cache
+    m_logEsrc     = -30.0;
+    m_theta       = 0.0;
+    m_cumul.clear();
+
     // Return
     return;
 }
@@ -700,3 +672,55 @@ void GCTAEdisp2D::free_members(void)
     // Return
     return;
 }
+
+/***********************************************************************//**
+ * @brief Update cumulative probability
+ *
+ * TODO : vector of vectors for each Esrc and select nearest vector in mc
+ *
+ ***************************************************************************/
+void GCTAEdisp2D::update_cumul(const double& logEsrc,
+                               const double& theta,
+                               const double& phi,
+                               const double& zenith,
+                               const double& azimuth) const
+{
+
+    if (logEsrc != m_logEsrc || theta != m_theta) {
+
+        m_cumul.clear();
+
+        m_logEsrc = logEsrc;
+        m_theta = theta;
+
+        // Initialize cumulative probability
+        double sum = 0.0;
+
+        for (int i = 0; i < m_edisp.axis(1); ++i) {
+
+            // Compute delta(Eobs/Etrue) and Eobs/Etrue values
+            double delta        = m_edisp.axis_hi(1, i) -
+                                  m_edisp.axis_lo(1, i);
+            double EobsOverEsrc = 0.5 * (m_edisp.axis_hi(1, i) +
+                                         m_edisp.axis_lo(1, i));
+
+            // Compute cumulative probability
+            double add = m_edisp(0, logEsrc, EobsOverEsrc, theta) *
+                         delta / EobsOverEsrc * gammalib::inv_ln10;
+
+            // Add to sum (and keep sum lower or equal to 1)
+            sum = sum+add >= 1.0 ? 1.0 : sum+add;
+
+            // Create pair containing EobsOverEsrc and cumulated probability
+            std::pair<double, double> pair(EobsOverEsrc, sum);
+            std::cout << "Eobs/Esrc = " << EobsOverEsrc << ", sum = " << sum << std::endl;
+
+            // Add to vector
+            m_cumul.push_back(pair);
+
+        } // endfor: looped over EobsOverEsrc
+    }
+
+    // Return
+    return;
+}
@@ -190,7 +190,7 @@ double GCTAEdispPerfTable::operator()(const double&amp; logEobs,
     // Compute energy dispersion value
     double delta = logEobs - logEsrc;
     double edisp = m_par_scale * std::exp(m_par_width * delta * delta);
-    
+
     // Return energy dispersion
     return edisp;
 }
@@ -249,8 +249,6 @@ GCTAEdispPerfTable* GCTAEdispPerfTable::clone(void) const
  ***************************************************************************/
 void GCTAEdispPerfTable::load(const std::string& filename)
 {
-    // Set conversion factor from RMS(ln(Eest/Etrue)) to RMS(log10(Eest/Etrue))
-    const double conv = 0.434294481903;
  
     // Clear arrays
     m_logE.clear();
@@ -292,7 +290,7 @@ void GCTAEdispPerfTable::load(const std::string&amp; filename)
  
         // Push elements in node array and vector
         m_logE.append(gammalib::todouble(elements[0]));
-        m_sigma.push_back(gammalib::todouble(elements[4]) * conv);
+        m_sigma.push_back(gammalib::todouble(elements[4]) * gammalib::inv_ln10);
  
     } // endwhile: looped over lines
  
@@ -40,8 +40,8 @@
 /* __ Macros _____________________________________________________________ */
  
 /* __ Coding definitions _________________________________________________ */
-#define G_OLD_MC_CODE
-#define G_MC_REJECTION
+//#define G_OLD_MC_CODE
+//#define G_MC_REJECTION
  
 /* __ Debug definitions __________________________________________________ */
  
@@ -333,7 +333,7 @@ GEnergy GCTAEdispRmf::mc(GRan&amp;         ran,
                 ftest  = ran.uniform() * fmax;
             }
             energy.log10TeV(e);
-            
+
             #else
             // Determine measured energy index from Monte-Carlo cache
             int imeasured = ran.cdf(m_mc_measured_cdf[itrue]) + offset;
@@ -348,7 +348,7 @@ GEnergy GCTAEdispRmf::mc(GRan&amp;         ran,
             // Set interval
             energy.log10TeV(e);
             #endif
-        
+
         } // endif: offset was valid
  
     } // endif: there was redistribution information
@@ -365,26 +365,46 @@ GEnergy GCTAEdispRmf::mc(GRan&amp;         ran,
                          const double& azimuth) const
 {
  
-    // Update cumul
-    update_cumul(logEsrc, theta, phi, zenith, azimuth);
+    // Compute cumulative probability
+    compute_cumul(theta, phi, zenith, azimuth);
+
+    // Find right logEsrc index with bisection search
+    int low = 0;
+    int high = m_rmf.ntrue();
+    while ((high-low) > 1) {
+        int mid = (low+high) / 2;
+        if (logEsrc < m_rmf.etrue().elogmean(mid).log10TeV()) {
+            high = mid;
+        }
+        else {
+            low = mid;
+        }
+    }
+    // Index found
+    int isrc = low;
  
     // Draw random number between 0 and 1 from uniform distribution
     double p = ran.uniform();
  
-    // Find right index
-    // TODO: This can lead to index = m_edisp.axis(1) on exit, which
-    //       leads to an array indexing problem in the interpolation.
-    //       Also, a bisection search would be more efficient.
-    int index = 0;
-    while(index < m_cumul.size() - 2 && m_cumul[index+1].second < p) {
-        index++;
-    } // index found
+    // Find right index with bisection search
+    low = 0;
+    high = m_cumul[isrc].size();
+    while ((high-low) > 1) {
+        int mid = (low+high) / 2;
+        if (p < m_cumul[isrc][mid].second) {
+            high = mid;
+        }
+        else if (m_cumul[isrc][mid].second <= p) {
+            low = mid;
+        }
+    }
+    // Index found
+    int index = low;
  
     // Interpolate Eobs value
-    double Eobs =   (m_cumul[index+1].second - p)*m_cumul[index].first
-                  + (p - m_cumul[index].second)*m_cumul[index+1].first;
-    Eobs       /=   (m_cumul[index+1].second - m_cumul[index].second);
-
+    double Eobs =   (m_cumul[isrc][index+1].second - p)*m_cumul[isrc][index].first
+                  + (p - m_cumul[isrc][index].second)*m_cumul[isrc][index+1].first;
+    Eobs       /=   (m_cumul[isrc][index+1].second - m_cumul[isrc][index].second);
  
     // Set energy
     GEnergy energy;
@@ -524,8 +544,8 @@ void GCTAEdispRmf::init_members(void)
     // Initialise Monte Carlo cache
     m_mc_measured_start.clear();
     m_mc_measured_cdf.clear();
-    m_logEsrc     = -30.0;
-    m_theta       = 0.0;
+    m_cdf_computed   = false;
+    m_theta          = 0.0;
     m_cumul.clear();
  
     // Return
@@ -562,9 +582,9 @@ void GCTAEdispRmf::copy_members(const GCTAEdispRmf&amp; edisp)
     // Copy Monte Carlo cache
     m_mc_measured_start = edisp.m_mc_measured_start;
     m_mc_measured_cdf   = edisp.m_mc_measured_cdf;
-    m_logEsrc           = edisp.m_logEsrc;
     m_theta             = edisp.m_theta;
     m_cumul             = edisp.m_cumul;
+    m_cdf_computed      = edisp.m_cdf_computed;
  
     // Return
     return;
@@ -770,19 +790,29 @@ void GCTAEdispRmf::update(const double&amp; etrue, const double&amp; emeasured) const
  
  
 /***********************************************************************//**
- * @brief Update cumulative probability
+ * @brief Compute cumulative probability
  *
  ***************************************************************************/
-void GCTAEdispRmf::update_cumul(const double& logEsrc,
-                                const double& theta,
-                                const double& phi,
-                                const double& zenith,
-                                const double& azimuth) const
+void GCTAEdispRmf::compute_cumul(const double& theta,
+                                 const double& phi,
+                                 const double& zenith,
+                                 const double& azimuth) const
 {
-    if (logEsrc != m_logEsrc || theta != m_theta) {
+    // TODO ? compute also for each theta ??
+
+    if (!m_cdf_computed) {
  
-        m_logEsrc = logEsrc;
-        m_theta = theta;
+        m_cumul.clear();
+
+        //m_theta = theta;
+        m_cdf_computed = true;
+
+    // Loop over Esrc
+    for (int isrc = 0; isrc < m_rmf.ntrue(); ++isrc) {
+
+        m_cumul.push_back(std::vector<std::pair<double, double> >());
+
+        double logEsrc = m_rmf.etrue().elogmean(isrc).log10TeV();
  
         // Initialize cumulative probability
         double sum = 0.0;
@@ -802,21 +832,22 @@ void GCTAEdispRmf::update_cumul(const double&amp; logEsrc,
                 double Eobs = Eobsmin+i*deltaEobs/n;
  
                 // Compute cumulative probability
-                double add = GCTAEdispRmf::operator()(logEsrc, std::log10(Eobs), theta) * deltaEobs / n / Eobs / std::log(10.0);
+                double add = GCTAEdispRmf::operator()(logEsrc, std::log10(Eobs), theta) * deltaEobs / (n * Eobs * gammalib::ln10);
  
-                //sum = sum+add/n >= 1.0 ? 1.0 : sum+add/n;
+                // Add to sum (and keep sum lower or equal to 1)
                 sum = sum+add >= 1.0 ? 1.0 : sum+add;
  
                 // Create pair containing Eobs and cumulative probability
                 std::pair<double, double> pair(Eobs, sum);
  
                 // Add to vector
-                m_cumul.push_back(pair);
+                m_cumul[isrc].push_back(pair);
             }
  
         }
     }
  
+    }
     // Return
     return;
 }
@@ -140,6 +140,7 @@ if __name__ == &#39;__main__&#39;:
     # Load edisp
     edisp = gammalib.GCTAEdispRmf("./caldb/dc1/rmf.fits")
     #edisp = gammalib.GCTAEdispPerfTable("./caldb/cta_dummy_irf.dat")
+    # NE MARCHE PAS /!\
     #edisp = gammalib.GCTAEdisp2D("./caldb/edisp_matrix.fits")
  
     # Simulate Edisp