diff --git a/corsika/detail/modules/epos/InteractionModel.inl b/corsika/detail/modules/epos/InteractionModel.inl
index dbaeba716a2aba70d7d3bf5a55dd617d2585609f..5d7996f9ae406c180cdf7dd149b181676b2e954e 100644
--- a/corsika/detail/modules/epos/InteractionModel.inl
+++ b/corsika/detail/modules/epos/InteractionModel.inl
@@ -10,9 +10,8 @@
#include <corsika/modules/epos/InteractionModel.hpp>
#include <corsika/modules/epos/EposStack.hpp>
-
#include <corsika/framework/geometry/Point.hpp>
-
+#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/utility/COMBoost.hpp>
#include <corsika/framework/utility/CorsikaData.hpp>
@@ -35,8 +34,8 @@ namespace corsika::epos {
data_path_ = (std::string(corsika_data("EPOS").c_str()) + "/").c_str();
}
initialize();
+ setParticlesStable();
}
- setParticlesStable();
}
inline void InteractionModel::setParticlesStable() const {
@@ -46,10 +45,30 @@ namespace corsika::epos {
if (!is_hadron(p)) continue;
int const eid = convertToEposRaw(p);
if (eid != 0) {
- ::epos::nodcy_.nrnody = ::epos::nodcy_.nrnody + 1;
- ::epos::nodcy_.nody[::epos::nodcy_.nrnody - 1] = eid;
+ // LCOV_EXCL_START
+ // this is only a safeguard against messing up the epos internals by initializing
+ // more than once.
+ unsigned int const n_particles_stable_epos =
+ ::epos::nodcy_.nrnody; // avoid waring -Wsign-compare
+ if (n_particles_stable_epos < ::epos::mxnody) {
+ CORSIKA_LOGGER_TRACE(logger_, "setting {} with EposId={} stable inside EPOS.",
+ p, eid);
+ ::epos::nodcy_.nrnody = ::epos::nodcy_.nrnody + 1;
+ ::epos::nodcy_.nody[::epos::nodcy_.nrnody - 1] = eid;
+ } else {
+ CORSIKA_LOGGER_ERROR(logger_, "List of stable particles too long for Epos!");
+ throw std::runtime_error("Epos initialization error!");
+ }
+ // LCOV_EXCL_STOP
+ } else {
+ CORSIKA_LOG_TRACE(
+ "particle conversion Corsika-->Epos not known for {}. Using {}. Setting "
+ "unstable in Epos!",
+ p, eid);
}
}
+ CORSIKA_LOGGER_DEBUG(logger_, "set {} particles stable inside Epos",
+ ::epos::nodcy_.nrnody);
}
inline bool InteractionModel::isValid(Code const projectileId, Code const targetId,
@@ -59,7 +78,7 @@ namespace corsika::epos {
if (!is_nucleus(targetId) && targetId != Code::Neutron && targetId != Code::Proton) {
return false;
}
- if (is_nucleus(targetId) && (get_nucleus_A(targetId) >= maxTargetMassNumber_)) {
+ if (is_nucleus(targetId) && (get_nucleus_A(targetId) >= get_nucleus_A(maxNucleus_))) {
return false;
}
if ((minEnergyCoM_ > sqrtS) || (sqrtS > maxEnergyCoM_)) { return false; }
@@ -136,9 +155,10 @@ namespace corsika::epos {
strcpy(::epos::fname_.fncs, CS.data);
::epos::nfname_.nfncs = CS.length;
- // initialiazes maximum energy and mass
- initializeEventCoM(Code::Lead, Lead::nucleus_A, Lead::nucleus_Z, Code::Lead,
- Lead::nucleus_A, Lead::nucleus_Z, 1_PeV);
+ // initializes maximum energy and mass
+ initializeEventCoM(
+ maxNucleus_, get_nucleus_A(maxNucleus_), get_nucleus_Z(maxNucleus_), maxNucleus_,
+ get_nucleus_A(maxNucleus_), get_nucleus_Z(maxNucleus_), maxEnergyCoM_);
}
inline void InteractionModel::initializeEventCoM(Code const idBeam, int const iBeamA,
@@ -284,7 +304,7 @@ namespace corsika::epos {
"beamA={}, "
"beamZ={} "
"targetId={}, "
- "ELab={:4.3f} GeV,",
+ "ELab={:12.2f} GeV,",
BeamId, BeamA, BeamZ, TargetId, EnergyLab / 1_GeV);
// read cross section from epos internal tables
@@ -300,7 +320,7 @@ namespace corsika::epos {
int const iBeam = epos::getEposXSCode(
BeamId); // 0 (can not interact, 1: pion-like, 2: proton-like, 3:kaon-like)
CORSIKA_LOGGER_TRACE(logger_,
- "projectile cross section type={} "
+ "readCrossSectionTableLab: projectile cross section type={} "
"(0: cannot interact, 1:pion, 2:baryon, 3:kaon)",
iBeam);
@@ -314,16 +334,25 @@ namespace corsika::epos {
int iMode = 3; // 0: air, >0 not air
- CORSIKA_LOGGER_DEBUG(logger_,
- "inside Epos "
+ CORSIKA_LOGGER_TRACE(logger_,
+ "readCrossSectionTableLab: inside Epos "
"beamId={}, beamXS={}",
::epos::hadr2_.idproj, ::epos::had10_.iclpro);
+ CORSIKA_LOGGER_TRACE(logger_,
+ "readCrossSectionTableLab: calling Epos cross section with"
+ "Ekin = {}, Abeam = {}, Atarget = {}, iMode = {}",
+ Ekin, Abeam, Atarget, iMode);
+
// cross section from table, FAST
float sigProdEpos = ::epos::eposcrse_(Ekin, Abeam, Atarget, iMode);
// sig-el from analytic calculation, no fast
float sigElaEpos = ::epos::eposelacrse_(Ekin, Abeam, Atarget, iMode);
+ CORSIKA_LOGGER_TRACE(logger_,
+ "readCrossSectionTableLab: result: sigProd = {}, sigEla = {}",
+ sigProdEpos, sigElaEpos);
+
return std::make_tuple(sigProdEpos * 1_mb, sigElaEpos * 1_mb);
}
diff --git a/corsika/modules/epos/InteractionModel.hpp b/corsika/modules/epos/InteractionModel.hpp
index b6701bfb5423669018975ea6edc47ff33a6b192b..e9ab49ae5f5c4723c3489e4bd984f798c9466749 100644
--- a/corsika/modules/epos/InteractionModel.hpp
+++ b/corsika/modules/epos/InteractionModel.hpp
@@ -89,16 +89,18 @@ namespace corsika::epos {
void doInteraction(TSecondaries&, Code const projectileId, Code const targetId,
FourMomentum const& projectileP4, FourMomentum const& targetP4);
- void initialize() const;
void initializeEventCoM(Code const, int const, int const, Code const, int const,
int const, HEPEnergyType const) const;
void initializeEventLab(Code const, int const, int const, Code const, int const,
int const, HEPEnergyType const) const;
void configureParticles(Code const, int const, int const, Code const, int const,
int const) const;
- void setParticlesStable() const;
private:
+ // initialize and setParticlesStable are private since they can only be called once at
+ // the beginning and are already called in the constructor!
+ void initialize() const;
+ void setParticlesStable() const;
inline static bool isInitialized_ = false;
std::string data_path_;
@@ -109,8 +111,7 @@ namespace corsika::epos {
std::shared_ptr<spdlog::logger> logger_ = get_logger("corsika_epos_Interaction");
HEPEnergyType const minEnergyCoM_ = 6 * 1e9 * electronvolt;
HEPEnergyType const maxEnergyCoM_ = 2.e6 * 1e9 * electronvolt;
- static unsigned int constexpr maxTargetMassNumber_ = 20;
- static unsigned int constexpr minNuclearTargetA_ = 4;
+ static Code constexpr maxNucleus_ = Code::Lead;
};
} // namespace corsika::epos
diff --git a/modules/epos/epos-sem-lhc.f b/modules/epos/epos-sem-lhc.f
index c02b67f582929499c1b89a80ac5f138a069e5985..218b4223090945b2b641f57cec73836ae52e3fa4 100644
--- a/modules/epos/epos-sem-lhc.f
+++ b/modules/epos/epos-sem-lhc.f
@@ -2430,6 +2430,7 @@ c p=sqrt(max(0.,egy**2-ampro**2))
maproj=mapr
matarg=matg
engy=egy
+ call Class("pscrse-1")
if(matg.eq.1.and.mapr.eq.1)then
if(iqq.eq.1)then !sig ela
call psfz(1,gz2,0.)
@@ -2479,6 +2480,7 @@ c p=sqrt(max(0.,egy**2-ampro**2))
maproj=maprojsave
matarg=matargsave
engy=engysave
+ call Class("pscrse-2")
else
ye=log10(max(1.,egy/1.5))+1.
je=min(5,int(ye))
@@ -2597,7 +2599,7 @@ c------------------------------------------------------------------------------
else
eposelacrse=pscrse(ek,mapro,matar,1)
endif
-
+
return
end
diff --git a/src/modules/epos/epos_codes.dat b/src/modules/epos/epos_codes.dat
index 2a88f095014e6c2288da567add66c012bc2b0d90..a5ce4c6808969a36d2b6476a2b1931b5cbc93eb7 100644
--- a/src/modules/epos/epos_codes.dat
+++ b/src/modules/epos/epos_codes.dat
@@ -19,20 +19,30 @@ PiPlus 120 1 Pion
PiMinus -120 1 Pion
KPlus 130 1 Kaon
KMinus -130 1 Kaon
+K0 230 0 CannotInteract
+K0Bar -230 0 CannotInteract
+K0Long -20 1 Kaon
+K0Short 20 1 Kaon
Neutron 1220 1 Baryon
AntiNeutron -1220 1 Baryon
Proton 1120 1 Baryon
AntiProton -1120 1 Baryon
-Eta 220 0 CannotInteract
-Lambda0 2130 0 CannotInteract
-Lambda0Bar -2130 0 CannotInteract
+Lambda0 2130 1 Baryon
+Lambda0Bar -2130 1 Baryon
+
+SigmaPlus 1130 1 Baryon
+SigmaMinusBar -1130 1 Baryon
+Sigma0 1230 1 Baryon
+Sigma0Bar -1230 1 Baryon
+SigmaMinus 2230 1 Baryon
+SigmaPlusBar -2230 1 Baryon
-SigmaPlus 1130 0 CannotInteract
-SigmaMinusBar -1130 0 CannotInteract
-Sigma0 1230 0 CannotInteract
-Sigma0Bar -1230 0 CannotInteract
-SigmaMinus 2230 0 CannotInteract
-SigmaPlusBar -2230 0 CannotInteract
+SigmaStarPlus 1131 0 CannotInteract
+SigmaStarMinusBar -1131 0 CannotInteract
+SigmaStar0 1231 0 CannotInteract
+SigmaStar0Bar -1231 0 CannotInteract
+SigmaStarMinus 2231 0 CannotInteract
+SigmaStarPlusBar -2231 0 CannotInteract
DeltaPlus 1121 0 CannotInteract
DeltaMinusBar -1121 0 CannotInteract
@@ -41,14 +51,24 @@ Delta0Bar -1221 0 CannotInteract
DeltaMinus 2221 0 CannotInteract
DeltaPlusBar -2221 0 CannotInteract
+DeltaPlusPlus 1111 0 CannotInteract
+DeltaMinusMinusBar -1111 0 CannotInteract
+
Xi0 1330 0 CannotInteract
XiMinus 2330 0 CannotInteract
Xi0Bar -1330 0 CannotInteract
XiPlusBar -2330 0 CannotInteract
+XiStar0 1331 0 CannotInteract
+XiStarMinus 2331 0 CannotInteract
+XiStar0Bar -1331 0 CannotInteract
+XiStarPlusBar -2331 0 CannotInteract
+
OmegaMinus 3331 0 CannotInteract
OmegaPlusBar -3331 0 CannotInteract
+
+Eta 220 0 CannotInteract
EtaPrime 330 0 CannotInteract
Phi 331 0 CannotInteract
Omega 221 0 CannotInteract
@@ -56,9 +76,6 @@ Rho0 111 1 Pion
RhoPlus 121 0 CannotInteract
RhoMinus -121 0 CannotInteract
-DeltaPlusPlus 1111 0 CannotInteract
-DeltaMinusMinusBar -1111 0 CannotInteract
-
KStarPlus 131 0 CannotInteract
KStarMinus -131 0 CannotInteract
KStar0 231 0 CannotInteract
diff --git a/tests/modules/testEpos.cpp b/tests/modules/testEpos.cpp
index 2c9c079c17a718b4c91975fa45f83cab19fec867..3b219df6a9edeec534f25d9cbaa3a9fa83214bec 100644
--- a/tests/modules/testEpos.cpp
+++ b/tests/modules/testEpos.cpp
@@ -66,8 +66,8 @@ TEST_CASE("EposBasics", "module,process") {
SECTION("cross-section type") {
CHECK(corsika::epos::getEposXSCode(Code::Electron) == 0);
- CHECK(corsika::epos::getEposXSCode(Code::K0Long) == 0);
- CHECK(corsika::epos::getEposXSCode(Code::SigmaPlus) == 0);
+ CHECK(corsika::epos::getEposXSCode(Code::K0Long) == 3);
+ CHECK(corsika::epos::getEposXSCode(Code::SigmaPlus) == 2);
CHECK(corsika::epos::getEposXSCode(Code::KMinus) == 3);
CHECK(corsika::epos::getEposXSCode(Code::PiMinus) == 1);
CHECK(corsika::epos::getEposXSCode(Code::Proton) == 2);
@@ -148,7 +148,8 @@ TEST_CASE("Epos", "modules") {
CHECK_FALSE(model.isValid(Code::Proton, Code::Electron, 100_GeV));
CHECK(model.isValid(Code::Proton, Code::Hydrogen, 100_GeV));
CHECK(model.isValid(Code::Proton, Code::Helium, 100_GeV));
- CHECK_FALSE(model.isValid(Code::Proton, Code::Iron, 100_GeV));
+ CHECK_FALSE(model.isValid(Code::Proton, Code::Iron, 10_EeV));
+ CHECK_FALSE(model.isValid(Code::Proton, get_nucleus_code(240, 120), 10_EeV));
CHECK(model.isValid(Code::Proton, Code::Oxygen, 100_GeV));
}
@@ -239,51 +240,53 @@ TEST_CASE("Epos", "modules") {
CHECK(xs_prod2 / 1_mb == Approx(1076.7).margin(3.1));
}
- /*
- SECTION("InteractionInterface - invalid") {
- Code const pid = Code::Electron;
- HEPEnergyType const P0 = 10_TeV;
- auto [stack, viewPtr] = setup::testing::setup_stack(
- pid, P0, (setup::Environment::BaseNodeType* const)nodePtr, cs);
- setup::StackView& view = *viewPtr;
- CHECK_THROWS(model.doInteraction(
- view, pid, Code::Oxygen,
- {sqrt(static_pow<2>(P0) + static_pow<2>(get_mass(pid))), {cs, P0, 0_GeV,
- 0_GeV}}, {Oxygen::mass, {cs, 0_GeV, 0_GeV, 0_GeV}}));
- }
- */
- /*
- SECTION("InteractionInterface - nuclear projectile") {
-
- HEPEnergyType const P0 = 10_TeV;
- Code const pid = get_nucleus_code(40, 20);
- auto [stack, viewPtr] = setup::testing::setup_stack(
- pid, P0, (setup::Environment::BaseNodeType* const)nodePtr, cs);
- MomentumVector plab =
- MomentumVector(cs, {P0, 0_eV, 0_eV}); // this is secret knowledge about
- setupStack setup::StackView& view = *viewPtr;
-
- // @todo This is very obscure since it fails for -O2, but for both clang and gcc ???
- model.doInteraction(view, pid, Code::Oxygen,
- {sqrt(static_pow<2>(P0) + static_pow<2>(get_mass(pid))), plab},
- {Oxygen::mass, {cs, 0_GeV, 0_GeV, 0_GeV}});
-
- auto const pSum = sumMomentum(view, cs);
-
- CHECK(pSum.getComponents(cs).getX() / P0 == Approx(1).margin(0.05));
- CHECK(pSum.getComponents(cs).getY() / 1_GeV ==
- Approx(0).margin(0.5)); // this is not physics validation
- CHECK(pSum.getComponents(cs).getZ() / 1_GeV ==
- Approx(0).margin(0.5)); // this is not physics validation
-
- CHECK((pSum - plab).getNorm() / 1_GeV ==
- Approx(0).margin(plab.getNorm() * 0.05 / 1_GeV));
- CHECK(pSum.getNorm() / P0 == Approx(1).margin(0.05));
- // [[maybe_unused]] const GrammageType length =
- // model.getInteractionLength(particle);
- // CHECK(length / 1_g * 1_cm * 1_cm ==
- // Approx(30).margin(20)); // this is no physics validation
- }*/
+ SECTION("InteractionInterface - invalid") {
+ Code const pid = Code::Electron;
+ HEPEnergyType const P0 = 10_TeV;
+ auto [stack, viewPtr] = setup::testing::setup_stack(
+ pid, P0, (DummyEnvironment::BaseNodeType* const)nodePtr, cs);
+ test::StackView& view = *viewPtr;
+ CHECK_THROWS(model.doInteraction(
+ view, pid, Code::Oxygen,
+ {sqrt(static_pow<2>(P0) + static_pow<2>(get_mass(pid))), {cs, P0, 0_GeV, 0_GeV}},
+ {Oxygen::mass, {cs, 0_GeV, 0_GeV, 0_GeV}}));
+ }
+
+ SECTION("InteractionInterface - nuclear projectile") {
+
+ HEPEnergyType const P0 = 10_TeV;
+ Code const pid = get_nucleus_code(40, 20);
+ auto [stack, viewPtr] = setup::testing::setup_stack(
+ pid, P0, (DummyEnvironment::BaseNodeType* const)nodePtr, cs);
+ MomentumVector plab =
+ MomentumVector(cs, {P0, 0_eV, 0_eV}); // this is secret knowledge about
+ test::StackView& view = *viewPtr;
+
+ // @todo This is very obscure since it fails for -O2, but for both clang and gcc ???
+ model.doInteraction(view, pid, Code::Oxygen,
+ {sqrt(static_pow<2>(P0) + static_pow<2>(get_mass(pid))), plab},
+ {Oxygen::mass, {cs, 0_GeV, 0_GeV, 0_GeV}});
+
+ // simply check if stack is not empty after the event. Energy and momentum
+ // conservation will be tested elsewhere
+ CHECK(view.getSize() > 0);
+
+ // auto const pSum = sumMomentum(view, cs);
+
+ // CHECK(pSum.getComponents(cs).getX() / P0 == Approx(1).margin(0.05));
+ // CHECK(pSum.getComponents(cs).getY() / 1_GeV ==
+ // Approx(0).margin(0.5)); // this is not physics validation
+ // CHECK(pSum.getComponents(cs).getZ() / 1_GeV ==
+ // Approx(0).margin(0.5)); // this is not physics validation
+
+ // CHECK((pSum - plab).getNorm() / 1_GeV ==
+ // Approx(0).margin(plab.getNorm() * 0.05 / 1_GeV));
+ // CHECK(pSum.getNorm() / P0 == Approx(1).margin(0.05));
+ // [[maybe_unused]] const GrammageType length =
+ // model.getInteractionLength(particle);
+ // CHECK(length / 1_g * 1_cm * 1_cm ==
+ // Approx(30).margin(20)); // this is no physics validation
+ }
// SECTION("InteractionInterface")
{
@@ -316,4 +319,4 @@ TEST_CASE("Epos", "modules") {
// CHECK(length / 1_g * 1_cm * 1_cm ==
// Approx(30).margin(20)); // this is no physics validation
}
-}
\ No newline at end of file
+}