diff --git a/corsika/detail/modules/proposal/ContinuousProcess.inl b/corsika/detail/modules/proposal/ContinuousProcess.inl index 16d419b0be200bb0d65b8748770818491d5e53b2..c5308315cd7bc779ea5bd7d7ad1cd0b0130a8aab 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 853302967fcc0bd460f6181f593254c9335b10bb..4874bda1f7fb7f6ca89d898cf7f1e01e8a51b791 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 b71b0fa9b1f9aff2784796588fb4ea7017a17d1f..a3eee5afa624a4b29d36498f5f79d3d579dff715 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 3bf53b145246392848f6ba09cbd293ac722c0a80..3b8976f3ac6bfe7a6d5251740782e4ca8fec89bc 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 69f27cca42b728cd52a77883cd377c6c3ccafc36..c6a389b1e2dcfcbb5308dcd53b9dd7b33589ef00 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 89037f9fead31d255b76e40180a265aaa00d3cb7..7a7eec3a14631bd2998935977d4048b31601f77e 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 d362ae14f5e3acb6ff6358ab44839229de4487f1..2afddc84ee9744909d4378eb3cf2456267655308 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 46bd6452209c2adc82f67c01c581137ff0379b1c..f012440e82deb19c3ab2b5220622c0255da9dbc5 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);