diff --git a/corsika/detail/modules/proposal/ContinuousProcess.inl b/corsika/detail/modules/proposal/ContinuousProcess.inl
index b190619b39b88f0b81ceafbb7fff4a73992d4862..36f277270a6e5a41f738589ac182b0ba69b8103f 100644
--- a/corsika/detail/modules/proposal/ContinuousProcess.inl
+++ b/corsika/detail/modules/proposal/ContinuousProcess.inl
@@ -28,11 +28,9 @@ namespace corsika::proposal {
if (code == Code::Photon) return; // no continuous builders needed for photons
// 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
+ // take some minutes if you have to build the tables and cannot read the tables
// from disk
- auto const emCut = get_energy_production_threshold(
- code); //! energy resolutions globally defined for individual particles
- auto c = p_cross->second(media.at(comp.getHash()), emCut);
+ auto c = p_cross->second(media.at(comp.getHash()), proposal_energycutsettings[code]);
// choose multiple scattering model
static constexpr auto ms_type = PROPOSAL::MultipleScatteringType::MoliereInterpol;
diff --git a/corsika/detail/modules/proposal/InteractionModel.inl b/corsika/detail/modules/proposal/InteractionModel.inl
index b5fb6cc8f6e8a2d62580dcf8ab9da10fe413cb49..05c44ba2fb725ddcf109e1cdef71793a81d7aa74 100644
--- a/corsika/detail/modules/proposal/InteractionModel.inl
+++ b/corsika/detail/modules/proposal/InteractionModel.inl
@@ -37,12 +37,9 @@ namespace corsika::proposal {
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
+ // take some minutes if you have to build the tables and cannot read the tables
// from disk
- auto const emCut = get_energy_production_threshold(
- code); //! energy resolutions globally defined for individual particles
-
- auto c = p_cross->second(media.at(comp.getHash()), emCut);
+ auto c = p_cross->second(media.at(comp.getHash()), proposal_energycutsettings[code]);
// Look which interactions take place and build the corresponding
// interaction and secondary builder. The interaction integral will
diff --git a/corsika/detail/modules/proposal/ProposalProcessBase.inl b/corsika/detail/modules/proposal/ProposalProcessBase.inl
index c7c38c10423ab6fe270c0851f7b77ca21d8ab73c..0ea18a4bf9cb6ed42a3a59bda43175110008506e 100644
--- a/corsika/detail/modules/proposal/ProposalProcessBase.inl
+++ b/corsika/detail/modules/proposal/ProposalProcessBase.inl
@@ -27,6 +27,27 @@ namespace corsika::proposal {
return false;
}
+ inline HEPEnergyType ProposalProcessBase::getOptimizedEmCut(Code code) const {
+ // get energy above which energy losses need to be considered
+ auto const production_threshold = get_energy_production_threshold(code);
+
+ HEPEnergyType lowest_table_value = 0_GeV;
+
+ // find tables for EnergyCuts closest (but still smaller than) production_threshold
+ for (auto table_energy : energycut_table_values) {
+ if (table_energy <= production_threshold && table_energy > lowest_table_value) {
+ lowest_table_value = table_energy;
+ }
+ }
+
+ if (lowest_table_value == 0_GeV) {
+ // no appropriate table available
+ return production_threshold;
+ }
+
+ return lowest_table_value;
+ };
+
template <typename TEnvironment>
inline ProposalProcessBase::ProposalProcessBase(TEnvironment const& _env) {
_env.getUniverse()->walk([&](auto& vtn) {
@@ -54,6 +75,34 @@ namespace corsika::proposal {
//! path, otherwise interpolation tables would only stored in main memory if
//! no explicit intrpolation def is specified.
PROPOSAL::InterpolationSettings::TABLES_PATH = corsika_data("PROPOSAL").c_str();
+
+ //! Initialize EnergyCutSettings
+ for (auto particle_code : tracked) {
+ if (particle_code == Code::Photon) {
+ // no EnergyCut for photon, only-stochastic propagation
+ continue;
+ }
+ // use optimized emcut for which tables should be available
+ auto optimized_ecut = getOptimizedEmCut(particle_code);
+ CORSIKA_LOG_INFO("PROPOSAL: Use tables with EnergyCut of {} for particle {}.",
+ optimized_ecut, particle_code);
+ proposal_energycutsettings[particle_code] = optimized_ecut;
+ }
+
+ //! Initialize PROPOSAL tables for all media and all particles
+ for (auto medium : media) {
+ for (auto particle_code : tracked) {
+ buildTables(medium.second, particle_code,
+ proposal_energycutsettings[particle_code]);
+ }
+ }
+ }
+
+ void ProposalProcessBase::buildTables(PROPOSAL::Medium medium, Code particle_code,
+ HEPEnergyType emcut) {
+ CORSIKA_LOG_DEBUG("PROPOSAL: Initialize tables for particle {} in {}", particle_code,
+ medium.GetName());
+ cross[particle_code](medium, emcut);
}
inline size_t ProposalProcessBase::hash::operator()(const calc_key_t& p) const
diff --git a/corsika/modules/proposal/ProposalProcessBase.hpp b/corsika/modules/proposal/ProposalProcessBase.hpp
index 0f13f54b4c2b8303b75afca262c67d761d7165c2..db56b02594db7870b11a97994fb4ea02e349eb68 100644
--- a/corsika/modules/proposal/ProposalProcessBase.hpp
+++ b/corsika/modules/proposal/ProposalProcessBase.hpp
@@ -32,6 +32,14 @@ namespace corsika::proposal {
//!
static constexpr double v_cut = 0.01;
+ //!
+ //! List of EnergyCut values for which CORSIKA will, by default, create and provide
+ //! PROPOSAL tables.
+ //!
+ static constexpr std::array<HEPEnergyType, 10> energycut_table_values{
+ 1000_MeV, 100_MeV, 20_MeV, 10_MeV, 3_MeV,
+ 1_MeV, 0.4_MeV, 0.25_MeV, 0.15_MeV, 0.05_MeV};
+
//!
//! Internal map from particle codes to particle properties required for
//! crosssections, decay and scattering algorithms. In the future the
@@ -223,6 +231,10 @@ namespace corsika::proposal {
media; //!< maps nuclear composition from univers to media to produce
//!< crosssections, which requires further ionization constants.
+ //!< save emcut for tracked particles
+ std::unordered_map<Code, corsika::units::si::HEPEnergyType>
+ proposal_energycutsettings;
+
//!
//! Store cut and nuclear composition of the whole universe in media which are
//! required for creating crosssections by proposal.
@@ -235,6 +247,12 @@ namespace corsika::proposal {
//!
bool canInteract(Code pcode) const;
+ //!
+ //! Finds the optimal EnergyCut for which PROPOSAL tables should (by default) be
+ //! available.
+ //!
+ HEPEnergyType getOptimizedEmCut(Code code) const;
+
using calc_key_t = std::pair<std::size_t, Code>;
//!
@@ -250,6 +268,11 @@ namespace corsika::proposal {
//!
virtual void buildCalculator(Code, NuclearComposition const&) = 0;
+ //!
+ //! Initialize PROPOSAL tables for given medium, code, and energy cut
+ //!
+ void buildTables(PROPOSAL::Medium, Code, HEPEnergyType);
+
//!
//! Searches the particle dependet calculator dependent of actuall medium composition
//! and particle type. If no calculator is found, the corresponding new calculator is
diff --git a/examples/corsika.cpp b/examples/corsika.cpp
index 7f40521ceff1cb5c3d822a8bc5f7ff5ff2b3f077..df55471074e45dc97306e3471dc00074b97182b1 100644
--- a/examples/corsika.cpp
+++ b/examples/corsika.cpp
@@ -398,12 +398,14 @@ int main(int argc, char** argv) {
HEPEnergyType const mucut = 1_GeV * app["--mucut"]->as<double>();
ParticleCut<SubWriter<decltype(dEdX)>> cut(emcut, emcut, hadcut, mucut, true, dEdX);
- // tell proposal that we are interested in all energy losses above the emcut
- set_energy_production_threshold(Code::Electron, emcut);
- set_energy_production_threshold(Code::Positron, emcut);
- set_energy_production_threshold(Code::Photon, emcut);
- set_energy_production_threshold(Code::MuMinus, mucut);
- set_energy_production_threshold(Code::MuPlus, mucut);
+ // tell proposal that we are interested in all energy losses above the particle cut
+ set_energy_production_threshold(Code::Electron, std::min({emcut, hadcut, mucut}));
+ set_energy_production_threshold(Code::Positron, std::min({emcut, hadcut, mucut}));
+ set_energy_production_threshold(Code::Photon, std::min({emcut, hadcut, mucut}));
+ set_energy_production_threshold(Code::MuMinus, std::min({emcut, hadcut, mucut}));
+ set_energy_production_threshold(Code::MuPlus, std::min({emcut, hadcut, mucut}));
+ set_energy_production_threshold(Code::TauMinus, std::min({emcut, hadcut, mucut}));
+ set_energy_production_threshold(Code::TauPlus, std::min({emcut, hadcut, mucut}));
// energy threshold for high energy hadronic model. Affects LE/HE switch for
// hadron interactions and the hadronic photon model in proposal