Edisp2D MC method is now working !

Florent Forest
1 parent 4c2771c1
Showing 6 changed files with 164 additions and 57 deletions
inst/cta/include/GCTAEdisp2D.hpp
inst/cta/include/GCTAEdispRmf.hpp
inst/cta/src/GCTAEdisp2D.cpp
inst/cta/src/GCTAEdispRmf.cpp
inst/cta/test/example_sim_edisp.py
inst/cta/test/test_CTA.cpp
@@ -103,8 +103,7 @@ private:
     void update(const double& logEobs,
                 const double& logEsrc,
                 const double& theta) const;
-    void update_cumul(const double& logEsrc,
-                      const double& theta = 0.0,
+    void compute_cumul(const double& theta = 0.0,
                       const double& phi = 0.0,
                       const double& zenith = 0.0,
                       const double& azimuth = 0.0) const;
@@ -112,9 +111,10 @@ private:
     // 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;
+
+    mutable bool      m_cdf_computed;
+    mutable double    m_theta;
+    mutable std::vector<std::vector<std::pair<double, double> > > m_cumul;
 };
  
  
@@ -129,7 +129,6 @@ private:
     mutable std::vector<GVector>                    m_mc_measured_cdf;
     mutable bool                                    m_cdf_computed;
     mutable double                                  m_theta;
-
     mutable std::vector<std::vector<std::pair<double, double> > > m_cumul;
  
 };
@@ -183,7 +183,7 @@ double GCTAEdisp2D::operator()(const double&amp; logEobs,
     double EobsOverEsrc = std::exp((logEobs-logEsrc) * gammalib::ln10);
  
     // Compute edisp
-    edisp = m_edisp(0, logEsrc, EobsOverEsrc, theta);
+    edisp = m_edisp(0, logEsrc, EobsOverEsrc, theta) * gammalib::ln10;
  
     // Return
     return edisp;
@@ -321,6 +321,9 @@ void GCTAEdisp2D::read(const GFits&amp; fits)
     // Set true energy axis to logarithmic scale
     m_edisp.axis_log10(0);
  
+    // Set emigra axis to log scale
+    m_edisp.axis_log10(1);
+
     // Set offset angle axis to radians
     m_edisp.axis_radians(2);
  
@@ -393,38 +396,70 @@ GEnergy GCTAEdisp2D::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_edisp.axis(0);
+    double Esrc = std::exp(logEsrc*gammalib::ln10);
+    while ((high-low) > 1) {
+        int mid = (low+high) / 2;
+        if (Esrc < m_edisp.axis_lo(0, mid)) {
+            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 with bisection search
-    int low = 0;
-    int high = m_cumul.size();
+    low = 0;
+    high = m_cumul[isrc].size();
     while ((high-low) > 1) {
         int mid = (low+high) / 2;
-        if (p < m_cumul[mid].second) {
+        if (p < m_cumul[isrc][mid].second) {
             high = mid;
         }
-        else if (m_cumul[mid].second <= p) {
+        else if (m_cumul[isrc][mid].second <= p) {
             low = mid;
         }
     }
+
     // Index found
-    int index = low;
+    int index = 0;
+    if (low >= m_cumul[isrc].size() - 1) {
+        index = m_cumul[isrc].size() - 2;
+    }
+    else {
+        index = low;
+    }
+
+    //std::cout << "p = " << p << std::endl;
+    //std::cout << "index = " << index << std::endl;
  
     // Interpolate EobsOverEsrc value
-    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);
+    double EobsOverEsrc = 0.0;
+
+    if (m_cumul[isrc][index+1].second == m_cumul[isrc][index].second) {
+        EobsOverEsrc = m_cumul[isrc][index].first;
+    }
+    else {
+        EobsOverEsrc  =   (m_cumul[isrc][index+1].second - p)*m_cumul[isrc][index].first
+                        + (p - m_cumul[isrc][index].second)*m_cumul[isrc][index+1].first;
+        EobsOverEsrc /=   (m_cumul[isrc][index+1].second - m_cumul[isrc][index].second);
+    }
  
-    // Final result
-    double logEnergy = std::log10(result) + logEsrc;
  
     // Set energy
     GEnergy energy;
-    energy.log10TeV(logEnergy);
+    energy.TeV(EobsOverEsrc * Esrc);
+    //std::cout << "eobs = " << EobsOverEsrc * Esrc << std::endl;
  
     // Return energy
     return energy;
@@ -639,8 +674,8 @@ void GCTAEdisp2D::init_members(void)
     m_edisp.clear();
  
     // Initialise Monte Carlo cache
-    m_logEsrc     = -30.0;
-    m_theta       = 0.0;
+    m_theta        = 0.0;
+    m_cdf_computed = false;
     m_cumul.clear();
  
     // Return
@@ -674,51 +709,100 @@ void GCTAEdisp2D::free_members(void)
 }
  
 /***********************************************************************//**
- * @brief Update cumulative probability
- *
- * TODO : vector of vectors for each Esrc and select nearest vector in mc
+ * @brief Compute cumulative probability
  *
  ***************************************************************************/
-void GCTAEdisp2D::update_cumul(const double& logEsrc,
-                               const double& theta,
-                               const double& phi,
-                               const double& zenith,
-                               const double& azimuth) const
+void GCTAEdisp2D::compute_cumul(const double& theta,
+                                const double& phi,
+                                const double& zenith,
+                                const double& azimuth) const
 {
+    // TODO ? compute also for each theta ??
  
-    if (logEsrc != m_logEsrc || theta != m_theta) {
+    if (!m_cdf_computed) {
  
         m_cumul.clear();
  
-        m_logEsrc = logEsrc;
-        m_theta = theta;
+        //m_theta = theta;
+        m_cdf_computed = true;
+
+    // Loop over Esrc
+    for (int isrc = 0; isrc < m_edisp.axis(0); ++isrc) {
+
+        m_cumul.push_back(std::vector<std::pair<double, double> >());
+
+        //double Esrc = 0.5 * (m_edisp.axis_lo(0, isrc) + m_edisp.axis_hi(0, isrc));
+        double Esrc = m_edisp.axis_lo(0, isrc);
  
         // Initialize cumulative probability
         double sum = 0.0;
  
-        for (int i = 0; i < m_edisp.axis(1); ++i) {
+        // Loop over Eobs from 0.1 TeV to 100.0 TeV
+        //int nsteps = (int) (m_edisp.axis_hi(1, m_edisp.axis(1)-1) * Esrc / deltaEobs);
+        int nsteps = 10000;
+        double deltaEobs = 100.0/nsteps;
+        double Eobs = 0.1;
  
-            // 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));
+        for (int iobs = 0; iobs < nsteps; ++iobs) {
+
+            double EobsOverEsrc = Eobs/Esrc;
  
             // Compute cumulative probability
-            double add = m_edisp(0, logEsrc, EobsOverEsrc, theta) *
-                         delta / EobsOverEsrc * gammalib::inv_ln10;
+            double add = m_edisp(0, std::log10(Esrc), EobsOverEsrc, theta) *
+                         deltaEobs / Eobs;
  
             // 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);
+            m_cumul[isrc].push_back(pair);
+
+            // Compute new Eobs value
+            Eobs += deltaEobs;
+
+        } // endfor: looped over Eobs
+
+        /* // Loop over EobsOverEtrue
+
+        // Divide bin into n sub-bins
+        const int n = 10;
+
+        for (int imigra = 0; imigra < m_edisp.axis(1); ++imigra) {
+
+            // Compute delta(Eobs/Etrue) and Eobs/Etrue values
+            double delta        = m_edisp.axis_hi(1, imigra) -
+                                  m_edisp.axis_lo(1, imigra);
+            double EobsOverEsrcmin = m_edisp.axis_lo(1, imigra);
+
+            for (int i = 0; i < n; ++i) {
+
+                double EobsOverEsrc = EobsOverEsrcmin + i*delta/n;
+
+                // Compute cumulative probability
+                double add = m_edisp(0, std::log10(Esrc), EobsOverEsrc, theta) *
+                             delta / (n * EobsOverEsrc * gammalib::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);
+
+                //if(Esrc == 10.0 && EobsOverEsrc < 2.0)
+                //std::cout << "Eobs/Esrc = " << EobsOverEsrc << ", add = " << add << ", sum = " << sum << std::endl;
+
+                // Add to vector
+                m_cumul[isrc].push_back(pair);
+
+            }
  
         } // endfor: looped over EobsOverEsrc
+        */
+    }
+
     }
  
     // Return
@@ -373,7 +373,7 @@ GEnergy GCTAEdispRmf::mc(GRan&amp;         ran,
     int high = m_rmf.ntrue();
     while ((high-low) > 1) {
         int mid = (low+high) / 2;
-        if (logEsrc < m_rmf.etrue().elogmean(mid).log10TeV()) {
+        if (logEsrc < m_rmf.etrue().emin(mid).log10TeV()) {
             high = mid;
         }
         else {
@@ -398,19 +398,33 @@ GEnergy GCTAEdispRmf::mc(GRan&amp;         ran,
             low = mid;
         }
     }
+
     // Index found
-    int index = low;
+    int index = 0;
+    if (low >= m_cumul[isrc].size() - 1) {
+        index = m_cumul[isrc].size() - 2;
+    }
+    else {
+        index = low;
+    }
  
     // Interpolate Eobs value
-    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);
+    double Eobs = 0.0;
+
+    if (m_cumul[isrc][index+1].second == m_cumul[isrc][index].second) {
+        Eobs = m_cumul[isrc][index].first;
+    }
+    else {
+        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;
     energy.TeV(Eobs);
  
-
     // Return energy
     return energy;
 }
@@ -31,7 +31,7 @@ import time
 # ============== #
 # Simulate Edisp #
 # ============== #
-def sim_edisp(edisp, etrue, eobs_max=40.0, ebins=500, nmc=100000):
+def sim_edisp(edisp, etrue, eobs_max=40.0, ebins=1000, nmc=100000):
     """
     Simulate Edisp and show results using matplotlib (if available).
     """
@@ -64,8 +64,8 @@ def sim_edisp(edisp, etrue, eobs_max=40.0, ebins=500, nmc=100000):
         # Simulate observed energies
         tstart = time.clock()
         for i in range(nmc):
-            #eobs     = edisp.mc(ran, logEtrue).TeV()
-            eobs     = edisp.mc(ran, logEtrue, float(i)).TeV()
+            eobs     = edisp.mc(ran, logEtrue).TeV()
+            #eobs     = edisp.mc(ran, logEtrue, float(i)).TeV()
             index    = int(eobs/dEobs)
             if (index < ebins):
                 counts[index] += 1.0
@@ -89,7 +89,7 @@ def sim_edisp(edisp, etrue, eobs_max=40.0, ebins=500, nmc=100000):
             sum += value
             exp_edisp.append(value)
             cumul.append(sum/nmc)
-        print(sum)
+        print(sum/nmc)
  
         # Create figure
         fig, ax1 = plt.subplots()
@@ -138,11 +138,9 @@ if __name__ == &#39;__main__&#39;:
     print("******************************")
  
     # Load edisp
-    edisp = gammalib.GCTAEdispRmf("./caldb/dc1/rmf.fits")
+    #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")
+    edisp = gammalib.GCTAEdisp2D("./caldb/edisp_matrix.fits")
  
     # Simulate Edisp
     sim_edisp(edisp, 20.0) # Energy in TeV
-    #sim_edisp(edisp, 1.0) # Energy in TeV
@@ -611,6 +611,12 @@ void TestGCTAResponse::test_response_edisp_RMF(void)
     // Test normalisation
     test_edisp_integration(edisp, 0.5, 10.0);
  
+    std::cout << "Eobs = 1 TeV : " << edisp(0.0, 0.0) << std::endl;
+    std::cout << "Eobs = 10 TeV : " << edisp(1.0, 0.0) << std::endl;
+
+    std::cout << "Eobs = 10 TeV : " << edisp(1.0, 1.0) << std::endl;
+    std::cout << "Eobs = 100 TeV : " << edisp(2.0, 1.0) << std::endl;
+
     // Return
     return;
 }
@@ -656,6 +662,12 @@ void TestGCTAResponse::test_response_edisp_2D(void)
         test_try_failure(e);
     }
  
+    std::cout << "Eobs = 1 TeV : " << edisp(0.0, 0.0) << std::endl;
+    std::cout << "Eobs = 10 TeV : " << edisp(1.0, 0.0) << std::endl;
+
+    std::cout << "Eobs = 10 TeV : " << edisp(1.0, 1.0) << std::endl;
+    std::cout << "Eobs = 100 TeV : " << edisp(2.0, 1.0) << std::endl;
+
     // Return
     return;
 }