diff --git a/corsika/detail/modules/proposal/HadronicPhotonModel.inl b/corsika/detail/modules/proposal/HadronicPhotonModel.inl
index 5aa5aded31ffe21ca19146055e829161581c0721..c6e30a4d34b905a94bc1f2f00d02cc4d52df197a 100644
--- a/corsika/detail/modules/proposal/HadronicPhotonModel.inl
+++ b/corsika/detail/modules/proposal/HadronicPhotonModel.inl
@@ -16,15 +16,32 @@
namespace corsika::proposal {
template <typename THadronicModel>
- inline HadronicPhotonModel<THadronicModel>::HadronicPhotonModel(THadronicModel& _hadint)
- : heHadronicInteraction_(_hadint){};
+ inline HadronicPhotonModel<THadronicModel>::HadronicPhotonModel(
+ THadronicModel& _hadint, HEPEnergyType const& _heenthresholdNN)
+ : heHadronicInteraction_(_hadint)
+ , heHadronicModelThresholdLabNN_(_heenthresholdNN) {
+ // check validity of threshold assuming photon-nucleon
+ // sqrtS per target nucleon
+ HEPEnergyType const sqrtS =
+ calculate_com_energy(_heenthresholdNN, Rho0::mass, Proton::mass);
+ if (!heHadronicInteraction_.isValid(Code::Rho0, Code::Proton, sqrtS)) {
+ CORSIKA_LOG_ERROR(
+ "Invalid energy threshold for hadron interaction model. theshold_lab= {} GeV, "
+ "theshold_com={} GeV",
+ _heenthresholdNN / 1_GeV, sqrtS / 1_GeV);
+ throw std::runtime_error("Configuration error!");
+ }
+ CORSIKA_LOG_INFO(
+ "Threshold for HE hadronic interactions in proposal set to Elab={} GeV",
+ _heenthresholdNN / 1_GeV);
+ };
template <typename THadronicModel>
template <typename TStackView>
inline ProcessReturn HadronicPhotonModel<THadronicModel>::doHadronicPhotonInteraction(
TStackView& view, CoordinateSystemPtr const& labCS, FourMomentum const& photonP4,
Code const& targetId) {
- if (photonP4.getTimeLikeComponent() > heHadronicModelThresholdLab_) {
+ if (photonP4.getTimeLikeComponent() > heHadronicModelThresholdLabNN_) {
CORSIKA_LOG_INFO(
"HE photo-hadronic interaction! calling hadronic interaction model..");
@@ -46,8 +63,7 @@ namespace corsika::proposal {
TStackView photon_secondaries(hadronicPhoton);
// call inner hadronic event generator
- CORSIKA_LOG_TRACE("calling HadronicInteraction...");
- CORSIKA_LOG_INFO("{} + {} interactions. Ekinlab = {} GeV", hadPhotonCode, targetId,
+ CORSIKA_LOG_INFO("{} + {} interaction. Ekinlab = {} GeV", hadPhotonCode, targetId,
photonP4.getTimeLikeComponent() / 1_GeV);
heHadronicInteraction_.doInteraction(photon_secondaries, hadPhotonCode, targetId,
photonP4, targetP4);
diff --git a/corsika/detail/modules/proposal/InteractionModel.inl b/corsika/detail/modules/proposal/InteractionModel.inl
index d1f63c685eb827d78d83321bbe49e3aa014c6be9..867b0385836d47f68f844dd449974d67e9469444 100644
--- a/corsika/detail/modules/proposal/InteractionModel.inl
+++ b/corsika/detail/modules/proposal/InteractionModel.inl
@@ -20,10 +20,11 @@ namespace corsika::proposal {
template <typename THadronicModel>
template <typename TEnvironment>
- inline InteractionModel<THadronicModel>::InteractionModel(TEnvironment const& _env,
- THadronicModel& _hadint)
+ inline InteractionModel<THadronicModel>::InteractionModel(
+ TEnvironment const& _env, THadronicModel& _hadint,
+ HEPEnergyType const& _enthreshold)
: ProposalProcessBase(_env)
- , HadronicPhotonModel<THadronicModel>(_hadint) {}
+ , HadronicPhotonModel<THadronicModel>(_hadint, _enthreshold) {}
template <typename THadronicModel>
inline void InteractionModel<THadronicModel>::buildCalculator(
@@ -124,10 +125,8 @@ namespace corsika::proposal {
auto const Z = int(target.GetNucCharge());
Code const targetId = get_nucleus_code(A, Z);
CORSIKA_LOG_INFO(
- "photo-hadronic interaction! projectile={} target={} energy={} GeV",
- projectileId,
- (is_nucleus(targetId) ? get_nucleus_name(targetId) : get_name(targetId)),
- E / 1_GeV);
+ "photo-hadronic interaction of projectile={} with target={}! Energy={} GeV",
+ projectileId, targetId, E / 1_GeV);
this->doHadronicPhotonInteraction(view, labCS, photonP4, targetId);
} else {
auto sec_code = convert_from_PDG(static_cast<PDGCode>(s.type));
diff --git a/corsika/modules/PROPOSAL.hpp b/corsika/modules/PROPOSAL.hpp
index a019ac63d0ad3714999394c604a4cf3381345fcb..13e18c3de7ed14f8dc6a7296f95bf4077ee59e9a 100644
--- a/corsika/modules/PROPOSAL.hpp
+++ b/corsika/modules/PROPOSAL.hpp
@@ -18,7 +18,7 @@ namespace corsika::proposal {
public InteractionProcess<Interaction<THadronicModel>> {
public:
template <typename TEnvironment>
- Interaction(TEnvironment const& env, THadronicModel& model)
- : InteractionModel<THadronicModel>(env, model) {}
+ Interaction(TEnvironment const& env, THadronicModel& model, HEPEnergyType const& thr)
+ : InteractionModel<THadronicModel>(env, model, thr) {}
};
} // namespace corsika::proposal
diff --git a/corsika/modules/proposal/HadronicPhotonModel.hpp b/corsika/modules/proposal/HadronicPhotonModel.hpp
index a3d9b606436491868cd6979c7822f594df86416b..befc2bd1b4b8b37df33df1defa0d66901f9e7c51 100644
--- a/corsika/modules/proposal/HadronicPhotonModel.hpp
+++ b/corsika/modules/proposal/HadronicPhotonModel.hpp
@@ -24,7 +24,7 @@ namespace corsika::proposal {
template <class THadronicModel>
class HadronicPhotonModel {
public:
- HadronicPhotonModel(THadronicModel&);
+ HadronicPhotonModel(THadronicModel&, HEPEnergyType const&);
//!
//! Calculate produce the hadronic secondaries in a hadronic photon interaction and
//! store them on the particle stack.
@@ -35,8 +35,8 @@ namespace corsika::proposal {
private:
THadronicModel& heHadronicInteraction_;
- static HEPEnergyType constexpr heHadronicModelThresholdLab_ =
- 80. * 1e9 * electronvolt; //!< energy threshold between LE and HE model
+ //! threshold for high energy hadronic interaction model. Lab. energy per nucleon
+ HEPEnergyType heHadronicModelThresholdLabNN_;
};
} // namespace corsika::proposal
diff --git a/corsika/modules/proposal/InteractionModel.hpp b/corsika/modules/proposal/InteractionModel.hpp
index ec5739225130ad1ada78f2740fc89aac12a51f36..135f685d1c14cf9d655567ecd639c163c594a4f6 100644
--- a/corsika/modules/proposal/InteractionModel.hpp
+++ b/corsika/modules/proposal/InteractionModel.hpp
@@ -55,7 +55,7 @@ namespace corsika::proposal {
//! compositions and stochastic description limited by the particle cut.
//!
template <typename TEnvironment>
- InteractionModel(TEnvironment const& env, THadronicModel&);
+ InteractionModel(TEnvironment const& env, THadronicModel&, HEPEnergyType const&);
//!
//! Calculate the rates for the different targets and interactions. Sample a
diff --git a/examples/corsika.cpp b/examples/corsika.cpp
index 68d428ab083a1cb28bfd85ea00e702d24f2aa936..9d4e05656025716094f2dd53d93fba482931a70a 100644
--- a/examples/corsika.cpp
+++ b/examples/corsika.cpp
@@ -321,7 +321,10 @@ int main(int argc, char** argv) {
HEPEnergyType const hadcut = 50_GeV;
ParticleCut<SubWriter<decltype(dEdX)>> cut(emcut, emcut, hadcut, hadcut, true, dEdX);
- corsika::proposal::Interaction emCascade(env, sibyll);
+ // energy threshold for high energy hadronic model. Affects LE/HE switch for hadron
+ // interactions and the hadronic photon model in proposal
+ HEPEnergyType heHadronModelThreshold = 63.1_GeV;
+ corsika::proposal::Interaction emCascade(env, sibyll, heHadronModelThreshold);
// NOT available for PROPOSAL due to interface trouble:
// InteractionCounter emCascadeCounted(emCascade);
// corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuous(env);
@@ -342,7 +345,8 @@ int main(int argc, char** argv) {
: cutE_(cutE) {}
bool operator()(const Particle& p) const { return (p.getKineticEnergy() < cutE_); }
};
- auto hadronSequence = make_select(EnergySwitch(63.1_GeV), urqmdCounted, heModelCounted);
+ auto hadronSequence =
+ make_select(EnergySwitch(heHadronModelThreshold), urqmdCounted, heModelCounted);
auto decaySequence = make_sequence(decayPythia, decaySibyll);
TrackWriter trackWriter{tracks};
diff --git a/examples/em_shower.cpp b/examples/em_shower.cpp
index f4e2c37a4917a9873b99257e50be738202a27369..aec2e1a3e355ad90dfdee2c8228cf2ed759ae3b4 100644
--- a/examples/em_shower.cpp
+++ b/examples/em_shower.cpp
@@ -164,7 +164,8 @@ int main(int argc, char** argv) {
ParticleCut<SubWriter<decltype(dEdX)>> cut(60_GeV, 60_GeV, 100_PeV, 100_PeV, true,
dEdX);
corsika::sibyll::Interaction sibyll;
- corsika::proposal::Interaction emCascade(env, sibyll);
+ HEPEnergyType heThresholdNN = 80_GeV;
+ corsika::proposal::Interaction emCascade(env, sibyll, heThresholdNN);
corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>> emContinuous(env, dEdX);
// BetheBlochPDG<SubWriter<decltype(dEdX)>> emContinuous{dEdX};
diff --git a/examples/mars.cpp b/examples/mars.cpp
index 23acc590eee759ea643cca0356588d6224449ecc..df9e531fa906a687801e49d529824442f6b54314 100644
--- a/examples/mars.cpp
+++ b/examples/mars.cpp
@@ -353,7 +353,11 @@ int main(int argc, char** argv) {
// decaySibyll.printDecayConfig();
- corsika::proposal::Interaction emCascade(env, sibyll);
+ // energy threshold for high energy hadronic model. Affects LE/HE switch for hadron
+ // interactions and the hadronic photon model in proposal
+ HEPEnergyType heHadronModelThreshold = 63.1_GeV;
+
+ corsika::proposal::Interaction emCascade(env, sibyll, heHadronModelThreshold);
// NOT possible right now, due to interface difference for PROPOSAL:
// InteractionCounter emCascadeCounted(emCascade);
// corsika::proposal::ContinuousProcess<SubWriter<decltype(dEdX)>>
@@ -375,7 +379,8 @@ int main(int argc, char** argv) {
: cutE_(cutE) {}
bool operator()(Particle const& p) const { return (p.getKineticEnergy() < cutE_); }
};
- auto hadronSequence = make_select(EnergySwitch(63.1_GeV), urqmdCounted, heModelCounted);
+ auto hadronSequence =
+ make_select(EnergySwitch(heHadronModelThreshold), urqmdCounted, heModelCounted);
auto decaySequence = make_sequence(decayPythia, decaySibyll);
// track writer
diff --git a/tests/modules/testProposal.cpp b/tests/modules/testProposal.cpp
index d4b3c62dd8d70bc033d2ee134460a11148e8f93f..7b6f148120891113bcd840688e2fe3e6ce4071dc 100644
--- a/tests/modules/testProposal.cpp
+++ b/tests/modules/testProposal.cpp
@@ -41,6 +41,9 @@ public:
MomentumVector(csPrime, {0_GeV, 0_GeV, 0_GeV}).normalized()));
}
}
+ bool constexpr isValid(Code const, Code const, HEPEnergyType const) const {
+ return true;
+ };
};
TEST_CASE("ProposalInterface", "modules") {
@@ -58,8 +61,8 @@ TEST_CASE("ProposalInterface", "modules") {
RNGManager<>::getInstance().registerRandomStream("proposal");
DummyHadronicModel hadModel;
-
- corsika::proposal::InteractionModel emModel(*env, hadModel);
+ HEPEnergyType heThresholdLab = 80_GeV;
+ corsika::proposal::InteractionModel emModel(*env, hadModel, heThresholdLab);
SECTION("InteractionInterface - cross section") {
auto& stack = *stackPtr;