From 71872b3d24e9e623b79d30615d07836671b9733c Mon Sep 17 00:00:00 2001 From: Domenik Ehlert Date: Wed, 6 Mar 2024 18:08:46 +0100 Subject: [PATCH 1/2] Implement optional logarithmic sampling of the photon energies --- src/module/PhotoPionProduction.cpp | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/src/module/PhotoPionProduction.cpp b/src/module/PhotoPionProduction.cpp index b8f6a7db2..a3f49c981 100644 --- a/src/module/PhotoPionProduction.cpp +++ b/src/module/PhotoPionProduction.cpp @@ -422,7 +422,11 @@ double PhotoPionProduction::sampleEps(bool onProton, double E, double z) const { Random &random = Random::instance(); for (int i = 0; i < 1000000; i++) { - double eps = epsMin + random.rand() * (epsMax - epsMin); + double eps; + if (sampleLog){ + eps = epsMin * pow(epsMax / epsMin, random.rand()); + } else + eps = epsMin + random.rand() * (epsMax - epsMin); double pEps = probEps(eps, onProton, Ein, z); if (random.rand() * pEpsMax < pEps) return eps * eV; From 12518ec93a6e00da3e9e248095b7aa4fbd83b2c8 Mon Sep 17 00:00:00 2001 From: Domenik Ehlert Date: Fri, 15 Mar 2024 09:21:12 +0100 Subject: [PATCH 2/2] Photon sampling with power law proposal function --- include/crpropa/module/PhotoPionProduction.h | 4 ++++ src/module/PhotoPionProduction.cpp | 25 ++++++++++++++++---- 2 files changed, 25 insertions(+), 4 deletions(-) diff --git a/include/crpropa/module/PhotoPionProduction.h b/include/crpropa/module/PhotoPionProduction.h index f736f86b7..3f2ff6772 100644 --- a/include/crpropa/module/PhotoPionProduction.h +++ b/include/crpropa/module/PhotoPionProduction.h @@ -43,6 +43,10 @@ class PhotoPionProduction: public Module { // - output: (s-p^2) * sigma_(nucleon/gamma) [GeV^2 * mubarn] double functs(double s, bool onProton) const; + double proposal_pdf(double eps, double epsMin, double epsMax) const; + + double proposal_inv_cdf(double prob, double epsMin, double epsMax) const; + // called by: sampleEps, gaussInt // - input: photon energy eps [eV], Ein [GeV] // - output: probability to encounter photon of energy eps diff --git a/src/module/PhotoPionProduction.cpp b/src/module/PhotoPionProduction.cpp index a3f49c981..b57388593 100644 --- a/src/module/PhotoPionProduction.cpp +++ b/src/module/PhotoPionProduction.cpp @@ -413,6 +413,16 @@ SophiaEventOutput PhotoPionProduction::sophiaEvent(bool onProton, double Ein, do return output; } +double PhotoPionProduction::proposal_pdf(double eps, double epsMin, double epsMax) const { + double prob = 1 / (eps * log(epsMax / epsMin)); + return prob; +} + +double PhotoPionProduction::proposal_inv_cdf(double prob, double epsMin, double epsMax) const { + double eps = epsMin * pow(epsMax / epsMin, prob); + return eps; +} + double PhotoPionProduction::sampleEps(bool onProton, double E, double z) const { // sample eps between epsMin ... epsMax double Ein = E / GeV; @@ -423,13 +433,20 @@ double PhotoPionProduction::sampleEps(bool onProton, double E, double z) const { Random &random = Random::instance(); for (int i = 0; i < 1000000; i++) { double eps; + double pEps; if (sampleLog){ - eps = epsMin * pow(epsMax / epsMin, random.rand()); + // convert uniformly-sampled probability to eps via inverse CDF + eps = proposal_inv_cdf(random.rand(), epsMin, epsMax); + // accept / reject based on ratio of pdfs + pEps = probEps(eps, onProton, Ein, z); + double pEps_proposal = proposal_pdf(eps, epsMin, epsMax); + if (pEps / pEps_proposal > random.rand() * pEpsMax) + return eps * eV; } else eps = epsMin + random.rand() * (epsMax - epsMin); - double pEps = probEps(eps, onProton, Ein, z); - if (random.rand() * pEpsMax < pEps) - return eps * eV; + pEps = probEps(eps, onProton, Ein, z); + if (random.rand() * pEpsMax < pEps) + return eps * eV; } throw std::runtime_error("error: no photon found in sampleEps, please make sure that photon field provides photons for the interaction by adapting the energy range of the tabulated photon field."); }