From 77404e91b26e65f2668572db3802db095d5fbe2b Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Mon, 18 Jan 2021 18:49:17 +0000
Subject: [PATCH] removed explicit threshold from proposal
---
.../modules/proposal/ContinuousProcess.inl | 25 +++++++++++--------
.../detail/modules/proposal/Interaction.inl | 9 ++++---
.../modules/proposal/ProposalProcessBase.inl | 7 +++---
.../modules/proposal/ContinuousProcess.hpp | 2 +-
corsika/modules/proposal/Interaction.hpp | 2 +-
.../modules/proposal/ProposalProcessBase.hpp | 10 +++++---
examples/em_shower.cpp | 4 +--
examples/vertical_EAS.cpp | 4 +--
8 files changed, 36 insertions(+), 27 deletions(-)
diff --git a/corsika/detail/modules/proposal/ContinuousProcess.inl b/corsika/detail/modules/proposal/ContinuousProcess.inl
index 16d419b0b..c5308315c 100644
--- a/corsika/detail/modules/proposal/ContinuousProcess.inl
+++ b/corsika/detail/modules/proposal/ContinuousProcess.inl
@@ -28,11 +28,13 @@ namespace corsika::proposal {
auto p_cross = cross.find(code);
if (p_cross == cross.end())
throw std::runtime_error("PROPOSAL could not find corresponding builder");
-
+
// interpolate the crosssection for given media and energy cut. These may
// take some minutes if you have to build the tables and cannot read the
// from disk
- auto c = p_cross->second(media.at(comp.getHash()), emCut_);
+ auto const emCut = get_energy_threshold(
+ code); //! energy thresholds globally defined for individual particles
+ auto c = p_cross->second(media.at(comp.getHash()), emCut);
// Build displacement integral and scattering object and interpolate them too and
// saved in the calc map by a key build out of a hash of composed of the component and
@@ -45,9 +47,9 @@ namespace corsika::proposal {
}
template <>
- ContinuousProcess::ContinuousProcess(setup::Environment const& _env,
- HEPEnergyType _emCut)
- : ProposalProcessBase(_env, _emCut) {}
+ ContinuousProcess::ContinuousProcess(setup::Environment const& _env
+ )
+ : ProposalProcessBase(_env) {}
template <>
void ContinuousProcess::scatter(setup::Stack::particle_type& vP,
@@ -115,21 +117,24 @@ namespace corsika::proposal {
template <>
LengthType ContinuousProcess::getMaxStepLength(setup::Stack::particle_type const& vP,
setup::Trajectory const& vT) {
-
- if (!canInteract(vP.getPID())) return meter * std::numeric_limits<double>::infinity();
+ auto const code = vP.getPID();
+ if (!canInteract(code)) return meter * std::numeric_limits<double>::infinity();
// Limit the step size of a conitnuous loss. The maximal continuous loss seems to be a
// hyper parameter which must be adjusted.
//
- // in any case: never go below 0.99*emCut_ This needs to be
- // slightly smaller than emCut_ since, either this Step is limited
+ auto const emCut = get_energy_threshold(
+ code); //! energy thresholds globally defined for individual particles
+
+ // in any case: never go below 0.99*emCut This needs to be
+ // slightly smaller than emCut since, either this Step is limited
// by energy_lim, then the particle is stopped in a very short
// range (before doing anythin else) and is then removed
// instantly. The exact position where it reaches emCut is not
// important, the important fact is that its E_kin is zero
// afterwards.
//
- auto energy_lim = std::max(0.9 * vP.getEnergy(), 0.99 * emCut_);
+ auto energy_lim = std::max(0.9 * vP.getEnergy(), 0.99 * emCut);
// solving the track integral for giving energy lim
auto c = getCalculator(vP, calc);
diff --git a/corsika/detail/modules/proposal/Interaction.inl b/corsika/detail/modules/proposal/Interaction.inl
index 853302967..4874bda1f 100644
--- a/corsika/detail/modules/proposal/Interaction.inl
+++ b/corsika/detail/modules/proposal/Interaction.inl
@@ -24,8 +24,8 @@
namespace corsika::proposal {
template <>
- Interaction::Interaction(setup::Environment const& _env, HEPEnergyType _emCut)
- : ProposalProcessBase(_env, _emCut) {}
+ Interaction::Interaction(setup::Environment const& _env)
+ : ProposalProcessBase(_env) {}
void Interaction::buildCalculator(Code code, NuclearComposition const& comp) {
// search crosssection builder for given particle
@@ -36,7 +36,10 @@ namespace corsika::proposal {
// interpolate the crosssection for given media and energy cut. These may
// take some minutes if you have to build the tables and cannot read the
// from disk
- auto c = p_cross->second(media.at(comp.getHash()), emCut_);
+ auto const emCut = get_energy_threshold(
+ code); //! energy thresholds globally defined for individual particles
+
+ auto c = p_cross->second(media.at(comp.getHash()), emCut);
// Look which interactions take place and build the corresponding
// interaction and secondarie builder. The interaction integral will
diff --git a/corsika/detail/modules/proposal/ProposalProcessBase.inl b/corsika/detail/modules/proposal/ProposalProcessBase.inl
index b71b0fa9b..a3eee5afa 100644
--- a/corsika/detail/modules/proposal/ProposalProcessBase.inl
+++ b/corsika/detail/modules/proposal/ProposalProcessBase.inl
@@ -30,10 +30,9 @@ namespace corsika::proposal {
return false;
}
- ProposalProcessBase::ProposalProcessBase(setup::Environment const& _env,
- HEPEnergyType _emCut)
- : emCut_(_emCut)
- , RNG_(RNGManager::getInstance().getRandomStream("proposal")) {
+ ProposalProcessBase::ProposalProcessBase(setup::Environment const& _env
+ )
+ : RNG_(RNGManager::getInstance().getRandomStream("proposal")) {
_env.getUniverse()->walk([&](auto& vtn) {
if (vtn.hasModelProperties()) {
const auto& prop = vtn.getModelProperties();
diff --git a/corsika/modules/proposal/ContinuousProcess.hpp b/corsika/modules/proposal/ContinuousProcess.hpp
index 3bf53b145..3b8976f3a 100644
--- a/corsika/modules/proposal/ContinuousProcess.hpp
+++ b/corsika/modules/proposal/ContinuousProcess.hpp
@@ -52,7 +52,7 @@ namespace corsika::proposal {
//! compositions and stochastic description limited by the particle cut.
//!
template <typename TEnvironment>
- ContinuousProcess(TEnvironment const&, HEPEnergyType _emCut);
+ ContinuousProcess(TEnvironment const&);
//!
//! Multiple Scattering of the lepton. Stochastic deflection is not yet taken into
diff --git a/corsika/modules/proposal/Interaction.hpp b/corsika/modules/proposal/Interaction.hpp
index 69f27cca4..c6a389b1e 100644
--- a/corsika/modules/proposal/Interaction.hpp
+++ b/corsika/modules/proposal/Interaction.hpp
@@ -48,7 +48,7 @@ namespace corsika::proposal {
//! compositions and stochastic description limited by the particle cut.
//!
template <typename TEnvironment>
- Interaction(TEnvironment const& env, HEPEnergyType emCut);
+ Interaction(TEnvironment const& env);
//!
//! Calculate the rates for the different targets and interactions. Sample a
diff --git a/corsika/modules/proposal/ProposalProcessBase.hpp b/corsika/modules/proposal/ProposalProcessBase.hpp
index 89037f9fe..7a7eec3a1 100644
--- a/corsika/modules/proposal/ProposalProcessBase.hpp
+++ b/corsika/modules/proposal/ProposalProcessBase.hpp
@@ -44,7 +44,11 @@ namespace corsika::proposal {
//!
template <typename T>
static auto cross_builder =
- [](PROPOSAL::Medium& m, corsika::units::si::HEPEnergyType emCut) {
+ [](PROPOSAL::Medium& m,
+ corsika::units::si::HEPEnergyType
+ emCut) { //!< Stochastic losses smaller than the given cut
+ //!< will be handeled continuously.
+
using namespace corsika::units::si;
auto p_cut =
std::make_shared<const PROPOSAL::EnergyCutSettings>(emCut / 1_MeV, 1, true);
@@ -73,8 +77,6 @@ namespace corsika::proposal {
//!
class ProposalProcessBase {
protected:
- HEPEnergyType emCut_; //!< Stochastic losses smaller than the given cut
- //!< will be handeled continuously.
RNGManager::prng_type RNG_; //!< random number generator used by proposal
std::unordered_map<std::size_t, PROPOSAL::Medium>
@@ -85,7 +87,7 @@ namespace corsika::proposal {
//! Store cut and nuclear composition of the whole universe in media which are
//! required for creating crosssections by proposal.
//!
- ProposalProcessBase(corsika::setup::Environment const& _env, HEPEnergyType _emCut);
+ ProposalProcessBase(corsika::setup::Environment const& _env);
//!
//! Checks if a particle can be processed by proposal
diff --git a/examples/em_shower.cpp b/examples/em_shower.cpp
index d362ae14f..2afddc84e 100644
--- a/examples/em_shower.cpp
+++ b/examples/em_shower.cpp
@@ -143,8 +143,8 @@ int main(int argc, char** argv) {
// PROPOSAL processs proposal{...};
ParticleCut cut(10_GeV, 10_GeV, 100_PeV, 100_PeV, true);
- corsika::proposal::Interaction proposal(env, cut.getElectronECut());
- corsika::proposal::ContinuousProcess em_continuous(env, cut.getElectronECut());
+ corsika::proposal::Interaction proposal(env);
+ corsika::proposal::ContinuousProcess em_continuous(env);
InteractionCounter proposalCounted(proposal);
TrackWriter trackWriter("tracks.dat");
diff --git a/examples/vertical_EAS.cpp b/examples/vertical_EAS.cpp
index 46bd64522..f012440e8 100644
--- a/examples/vertical_EAS.cpp
+++ b/examples/vertical_EAS.cpp
@@ -222,8 +222,8 @@ int main(int argc, char** argv) {
decaySibyll.printDecayConfig();
ParticleCut cut{60_GeV, 60_GeV, 60_GeV, 60_GeV, true};
- corsika::proposal::Interaction proposal(env, cut.getElectronECut());
- corsika::proposal::ContinuousProcess em_continuous(env, cut.getElectronECut());
+ corsika::proposal::Interaction proposal(env);
+ corsika::proposal::ContinuousProcess em_continuous(env);
InteractionCounter proposalCounted(proposal);
OnShellCheck reset_particle_mass(1.e-3, 1.e-1, false);
--
GitLab