diff --git a/corsika/detail/modules/proposal/InteractionModel.inl b/corsika/detail/modules/proposal/InteractionModel.inl
index 0c92b39b9684fa244d1bcb5fb1ae611f9f073194..2838c1676035cdb672313b48dea43f7bbfcaac1b 100644
--- a/corsika/detail/modules/proposal/InteractionModel.inl
+++ b/corsika/detail/modules/proposal/InteractionModel.inl
@@ -22,9 +22,9 @@ namespace corsika::proposal {
template <typename THadronicModel>
template <typename TEnvironment>
inline InteractionModel<THadronicModel>::InteractionModel(TEnvironment const& _env,
- THadronicModel& hadint)
+ THadronicModel& _hadint)
: ProposalProcessBase(_env)
- , hadronicInteraction_(hadint) {}
+ , HadronicPhotonModel<THadronicModel>(_hadint) {}
template <typename THadronicModel>
inline void InteractionModel<THadronicModel>::buildCalculator(
@@ -107,14 +107,11 @@ namespace corsika::proposal {
auto const photonEnergy = projectile.getEnergy() * v;
- if (type == PROPOSAL::InteractionType::Photonuclear)
+ if (type == PROPOSAL::InteractionType::Photonuclear) {
CORSIKA_LOG_INFO(
"HE photo-hadronic interaction! Energy = "
"{} GeV, v = {}, v * E = {}",
projectile.getEnergy() / 1_GeV, v, photonEnergy);
-
- if (type == PROPOSAL::InteractionType::Photonuclear &&
- photonEnergy > heHadronicModelThresholdLab_) {
auto const photonDirection =
projectileP4.getSpaceLikeComponents()
.normalized(); // photon collinear with projectile
@@ -143,45 +140,54 @@ namespace corsika::proposal {
return ProcessReturn::Ok;
}
+ template <typename THadronicModel>
+ inline HadronicPhotonModel<THadronicModel>::HadronicPhotonModel(THadronicModel& _hadint)
+ : heHadronicInteraction_(_hadint){};
+
template <typename THadronicModel>
template <typename TStackView>
- inline ProcessReturn InteractionModel<THadronicModel>::doHadronicPhotonInteraction(
+ inline ProcessReturn HadronicPhotonModel<THadronicModel>::doHadronicPhotonInteraction(
TStackView& view, CoordinateSystemPtr const& labCS, FourMomentum const& photonP4,
Code const& targetId) {
- CORSIKA_LOG_INFO(
- "HE photo-hadronic interaction! calling hadronic interaction model..");
-
- // copy from sibyll::NuclearInteractionModel
- // temporarily add to stack, will be removed after interaction in DoInteraction
- typename TStackView::inner_stack_value_type photonStack;
- Point const pDummy(labCS, {0_m, 0_m, 0_m});
- TimeType const tDummy = 0_ns;
- Code const hadPhotonCode = Code::Rho0; // stand in for hadronic-photon
- // target at rest
- FourMomentum const targetP4(get_mass(targetId),
- MomentumVector(labCS, {0_GeV, 0_GeV, 0_GeV}));
- auto hadronicPhoton = photonStack.addParticle(
- std::make_tuple(hadPhotonCode, photonP4.getTimeLikeComponent(),
- photonP4.getSpaceLikeComponents().normalized(), pDummy, tDummy));
- hadronicPhoton.setNode(view.getProjectile().getNode());
- // create inelastic interaction of the hadronic photon
- // create new StackView for the photon
- TStackView photon_secondaries(hadronicPhoton);
-
- // call inner hadronic event generator
- CORSIKA_LOG_TRACE("calling HadronicInteraction...");
- CORSIKA_LOG_INFO("{} + {} interactions. Ekinlab = {}", hadPhotonCode, targetId,
- photonP4.getTimeLikeComponent() / 1_GeV);
- hadronicInteraction_.doInteraction(photon_secondaries, hadPhotonCode, targetId,
- photonP4, targetP4);
- for (const auto& pSec : photon_secondaries) {
- auto const p3lab = pSec.getMomentum();
- Code const pid = pSec.getPID();
- HEPEnergyType const secEkin =
- calculate_kinetic_energy(p3lab.getNorm(), get_mass(pid));
- view.addSecondary(std::make_tuple(pid, secEkin, p3lab.normalized()));
+ if (photonP4.getTimeLikeComponent() > heHadronicModelThresholdLab_) {
+ CORSIKA_LOG_INFO(
+ "HE photo-hadronic interaction! calling hadronic interaction model..");
+
+ // copy from sibyll::NuclearInteractionModel
+ // temporarily add to stack, will be removed after interaction in DoInteraction
+ typename TStackView::inner_stack_value_type photonStack;
+ Point const pDummy(labCS, {0_m, 0_m, 0_m});
+ TimeType const tDummy = 0_ns;
+ Code const hadPhotonCode = Code::Rho0; // stand in for hadronic-photon
+ // target at rest
+ FourMomentum const targetP4(get_mass(targetId),
+ MomentumVector(labCS, {0_GeV, 0_GeV, 0_GeV}));
+ auto hadronicPhoton = photonStack.addParticle(std::make_tuple(
+ hadPhotonCode, photonP4.getTimeLikeComponent(),
+ photonP4.getSpaceLikeComponents().normalized(), pDummy, tDummy));
+ hadronicPhoton.setNode(view.getProjectile().getNode());
+ // create inelastic interaction of the hadronic photon
+ // create new StackView for the photon
+ TStackView photon_secondaries(hadronicPhoton);
+
+ // call inner hadronic event generator
+ CORSIKA_LOG_TRACE("calling HadronicInteraction...");
+ CORSIKA_LOG_INFO("{} + {} interactions. Ekinlab = {}", hadPhotonCode, targetId,
+ photonP4.getTimeLikeComponent() / 1_GeV);
+ heHadronicInteraction_.doInteraction(photon_secondaries, hadPhotonCode, targetId,
+ photonP4, targetP4);
+ for (const auto& pSec : photon_secondaries) {
+ auto const p3lab = pSec.getMomentum();
+ Code const pid = pSec.getPID();
+ HEPEnergyType const secEkin =
+ calculate_kinetic_energy(p3lab.getNorm(), get_mass(pid));
+ view.addSecondary(std::make_tuple(pid, secEkin, p3lab.normalized()));
+ }
+ CORSIKA_LOG_INFO("number of particles produced: {}", view.getEntries());
+ } else {
+ CORSIKA_LOG_INFO(
+ "LE photo-hadronic interaction! production of secondaries not implemented..");
}
- CORSIKA_LOG_INFO("number of particles produced: {}", view.getEntries());
return ProcessReturn::Ok;
}
diff --git a/corsika/modules/proposal/InteractionModel.hpp b/corsika/modules/proposal/InteractionModel.hpp
index e5663108d194b824ada0f08f2009622484976169..9c3381db2a527d4fbbd010599e65a7025214168d 100644
--- a/corsika/modules/proposal/InteractionModel.hpp
+++ b/corsika/modules/proposal/InteractionModel.hpp
@@ -33,7 +33,26 @@ namespace corsika::proposal {
//!
template <class THadronicModel>
- class InteractionModel : public ProposalProcessBase {
+ class HadronicPhotonModel {
+ public:
+ HadronicPhotonModel(THadronicModel&);
+ //!
+ //! Calculate produce the hadronic secondaries in a hadronic photon interaction and
+ //! store them on the particle stack.
+ //!
+ template <typename TSecondaryView>
+ ProcessReturn doHadronicPhotonInteraction(TSecondaryView&, CoordinateSystemPtr const&,
+ FourMomentum const&, Code const&);
+
+ private:
+ THadronicModel& heHadronicInteraction_;
+ static HEPEnergyType constexpr heHadronicModelThresholdLab_ =
+ 80. * 1e9 * electronvolt;
+ };
+
+ template <class THadronicModel>
+ class InteractionModel : public ProposalProcessBase,
+ public HadronicPhotonModel<THadronicModel> {
enum { eSECONDARIES, eINTERACTION };
using calculator_t = std::tuple<std::unique_ptr<PROPOSAL::SecondariesCalculator>,
@@ -64,25 +83,12 @@ namespace corsika::proposal {
ProcessReturn doInteraction(TSecondaryView&, Code const projectileId,
FourMomentum const& projectileP4);
- //!
- //! Calculate produce the hadronic secondaries in a hadronic photon interaction and
- //! store them on the particle stack.
- //!
- template <typename TSecondaryView>
- ProcessReturn doHadronicPhotonInteraction(TSecondaryView&, CoordinateSystemPtr const&,
- FourMomentum const&, Code const&);
-
//!
//! Calculates and returns the cross section.
//!
template <typename TParticle>
CrossSectionType getCrossSection(TParticle const& p, Code const projectileId,
FourMomentum const& projectileP4);
-
- private:
- THadronicModel& hadronicInteraction_;
- static HEPEnergyType constexpr heHadronicModelThresholdLab_ =
- 80. * 1e9 * electronvolt;
};
} // namespace corsika::proposal
diff --git a/tests/modules/testProposal.cpp b/tests/modules/testProposal.cpp
index 4575eeb08f32fbf0e4063fec81606834d5412c2a..8718f1c52966ef52791bfa66b6e24b50351a0711 100644
--- a/tests/modules/testProposal.cpp
+++ b/tests/modules/testProposal.cpp
@@ -70,7 +70,20 @@ TEST_CASE("ProposalInterface", "modules") {
CHECK(emModel.getCrossSection(particle, Code::Proton, P4) == 0_mb);
}
- SECTION("InteractionInterface - hadronic photon interaction") {
+ SECTION("InteractionInterface - LE hadronic photon interaction") {
+ auto& stack = *stackPtr;
+ // auto particle = stack.first();
+ FourMomentum P4(10_GeV, {cs, {10_GeV, 0_eV, 0_eV}});
+ // finish successfully
+ CHECK(emModel.doHadronicPhotonInteraction(view, cs, P4, Code::Oxygen) ==
+ ProcessReturn::Ok);
+ // no LE interactions
+ CHECK(stack.getEntries() == 1);
+ CORSIKA_LOG_INFO("Number of particles produced in hadronic photon interaction: {}",
+ stack.getEntries()-1);
+ }
+
+ SECTION("InteractionInterface - HE hadronic photon interaction") {
auto& stack = *stackPtr;
// auto particle = stack.first();
FourMomentum P4(100_TeV, {cs, {100_TeV, 0_eV, 0_eV}});
@@ -79,6 +92,6 @@ TEST_CASE("ProposalInterface", "modules") {
ProcessReturn::Ok);
CHECK(stack.getEntries() > 1);
CORSIKA_LOG_INFO("Number of particles produced in hadronic photon interaction: {}",
- stack.getEntries());
+ stack.getEntries()-1);
}
}
\ No newline at end of file