From eae25d7a9c2b4366546e89008152bcb456c35065 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Fri, 15 Apr 2022 11:27:49 +0100
Subject: [PATCH] fix interface
---
.../modules/proposal/HadronicPhotonModel.inl | 2 +-
.../modules/proposal/InteractionModel.inl | 56 +++++++++----------
tests/modules/testProposal.cpp | 7 ++-
3 files changed, 34 insertions(+), 31 deletions(-)
diff --git a/corsika/detail/modules/proposal/HadronicPhotonModel.inl b/corsika/detail/modules/proposal/HadronicPhotonModel.inl
index 1d47b2f8f..5aa5aded3 100644
--- a/corsika/detail/modules/proposal/HadronicPhotonModel.inl
+++ b/corsika/detail/modules/proposal/HadronicPhotonModel.inl
@@ -47,7 +47,7 @@ namespace corsika::proposal {
// call inner hadronic event generator
CORSIKA_LOG_TRACE("calling HadronicInteraction...");
- CORSIKA_LOG_INFO("{} + {} interactions. Ekinlab = {}", hadPhotonCode, targetId,
+ CORSIKA_LOG_INFO("{} + {} interactions. 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 44db596d8..d1f63c685 100644
--- a/corsika/detail/modules/proposal/InteractionModel.inl
+++ b/corsika/detail/modules/proposal/InteractionModel.inl
@@ -77,7 +77,8 @@ namespace corsika::proposal {
CORSIKA_LOG_WARN(
"PROPOSAL: No particle interaction possible. "
"Put initial particle back on stack.");
- view.addSecondary(std::make_tuple(projectileId, projectile.getEnergy(),
+ view.addSecondary(std::make_tuple(projectileId,
+ projectile.getEnergy() - get_mass(projectileId),
projectile.getDirection()));
return ProcessReturn::Ok;
}
@@ -104,34 +105,31 @@ namespace corsika::proposal {
if (type != PROPOSAL::InteractionType::Ioniz)
target = PROPOSAL::Component::GetComponentForHash(target_hash);
- auto const photonEnergy = projectile.getEnergy() * v;
-
- if (type == PROPOSAL::InteractionType::Photonuclear) {
- auto const photonDirection =
- projectileP4.getSpaceLikeComponents()
- .normalized(); // photon collinear with projectile
- FourMomentum const photonP4(photonEnergy, photonEnergy * photonDirection);
- Code const targetId =
- get_nucleus_code(target.GetAtomicNum(), target.GetNucCharge());
- CORSIKA_LOG_INFO(
- "photo-hadronic interaction ({} + {})! Energy = "
- "{} GeV, v = {}, Photon energy (v*E) = {} GeV",
- projectileId, targetId, projectile.getEnergy() / 1_GeV, v,
- photonEnergy / 1_GeV);
- this->doHadronicPhotonInteraction(view, labCS, photonP4, targetId);
- if (projectileId != Code::Photon)
- // add lepton, apply energy loss to kinetic energy
- view.addSecondary(std::make_tuple(
- projectileId, (1 - v) * (projectile.getEnergy() - get_mass(projectileId)),
- photonDirection));
- } else {
- auto sec =
- std::get<eSECONDARIES>(c->second)->CalculateSecondaries(loss, target, rnd);
- for (auto& s : sec) {
- auto E = s.energy * 1_MeV;
- auto vecProposal = s.direction;
- auto dir = DirectionVector(
- labCS, {vecProposal.GetX(), vecProposal.GetY(), vecProposal.GetZ()});
+ auto sec =
+ std::get<eSECONDARIES>(c->second)->CalculateSecondaries(loss, target, rnd);
+ for (auto& s : sec) {
+ auto E = s.energy * 1_MeV;
+ auto vecProposal = s.direction;
+ auto dir = DirectionVector(
+ labCS, {vecProposal.GetX(), vecProposal.GetY(), vecProposal.GetZ()});
+
+ if (static_cast<PROPOSAL::ParticleType>(s.type) ==
+ PROPOSAL::ParticleType::Hadron) {
+ FourMomentum const photonP4(E, E * dir);
+ // for PROPOSAL media target.GetAtomicNum() returns the atomic number in units
+ // of atomic mass, so Nitrogen is 14.0067. When using media in CORSIKA the same
+ // field is filled with the number of nucleons (ie. 14 for Nitrogen) To be sure
+ // we explicitly cast to int
+ auto const A = int(target.GetAtomicNum());
+ 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);
+ this->doHadronicPhotonInteraction(view, labCS, photonP4, targetId);
+ } else {
auto sec_code = convert_from_PDG(static_cast<PDGCode>(s.type));
view.addSecondary(std::make_tuple(sec_code, E - get_mass(sec_code), dir));
}
diff --git a/tests/modules/testProposal.cpp b/tests/modules/testProposal.cpp
index 650cfaab7..d4b3c62dd 100644
--- a/tests/modules/testProposal.cpp
+++ b/tests/modules/testProposal.cpp
@@ -65,9 +65,14 @@ TEST_CASE("ProposalInterface", "modules") {
auto& stack = *stackPtr;
auto particle = stack.first();
FourMomentum P4(
- 100_GeV,
+ 100_MeV,
{cs, {sqrt(static_pow<2>(100_MeV) - static_pow<2>(Proton::mass)), 0_eV, 0_eV}});
CHECK(emModel.getCrossSection(particle, Code::Proton, P4) == 0_mb);
+
+ FourMomentum eleP4(
+ 100_MeV,
+ {cs, {sqrt(static_pow<2>(100_MeV) - static_pow<2>(Electron::mass)), 0_eV, 0_eV}});
+ CHECK(emModel.getCrossSection(particle, Code::Electron, eleP4) > 0_mb);
}
SECTION("InteractionInterface - LE hadronic photon interaction") {
--
GitLab