diff --git a/inst/cta/include/GCTAEdisp2D.hpp b/inst/cta/include/GCTAEdisp2D.hpp index 6ff3f80..6032932 100644 --- a/inst/cta/include/GCTAEdisp2D.hpp +++ b/inst/cta/include/GCTAEdisp2D.hpp @@ -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 > m_cumul; + + mutable bool m_cdf_computed; + mutable double m_theta; + mutable std::vector > > m_cumul; }; diff --git a/inst/cta/include/GCTAEdispRmf.hpp b/inst/cta/include/GCTAEdispRmf.hpp index 73dcc74..1e4cbca 100644 --- a/inst/cta/include/GCTAEdispRmf.hpp +++ b/inst/cta/include/GCTAEdispRmf.hpp @@ -129,7 +129,6 @@ private: mutable std::vector m_mc_measured_cdf; mutable bool m_cdf_computed; mutable double m_theta; - mutable std::vector > > m_cumul; }; diff --git a/inst/cta/src/GCTAEdisp2D.cpp b/inst/cta/src/GCTAEdisp2D.cpp index 261aa4c..c1cb89d 100644 --- a/inst/cta/src/GCTAEdisp2D.cpp +++ b/inst/cta/src/GCTAEdisp2D.cpp @@ -183,7 +183,7 @@ double GCTAEdisp2D::operator()(const double& 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& 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& 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 >()); + + //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 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 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 diff --git a/inst/cta/src/GCTAEdispRmf.cpp b/inst/cta/src/GCTAEdispRmf.cpp index 7891a9a..73532ba 100644 --- a/inst/cta/src/GCTAEdispRmf.cpp +++ b/inst/cta/src/GCTAEdispRmf.cpp @@ -373,7 +373,7 @@ GEnergy GCTAEdispRmf::mc(GRan& 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& 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; } diff --git a/inst/cta/test/example_sim_edisp.py b/inst/cta/test/example_sim_edisp.py index 73cfcf6..3b6d385 100755 --- a/inst/cta/test/example_sim_edisp.py +++ b/inst/cta/test/example_sim_edisp.py @@ -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__ == '__main__': 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 diff --git a/inst/cta/test/test_CTA.cpp b/inst/cta/test/test_CTA.cpp index 69613e4..a5ea1fb 100644 --- a/inst/cta/test/test_CTA.cpp +++ b/inst/cta/test/test_CTA.cpp @@ -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; }