diff --git a/corsika/detail/modules/epos/InteractionModel.inl b/corsika/detail/modules/epos/InteractionModel.inl
index 97fe7ef4b077063650d3d0cae0798c22e0bdf735..039d9635de189d5c6a96309e6a96b4bd07365dd7 100644
--- a/corsika/detail/modules/epos/InteractionModel.inl
+++ b/corsika/detail/modules/epos/InteractionModel.inl
@@ -28,13 +28,12 @@ namespace corsika::epos {
: data_path_(dataPath)
, epos_listing_(epos_printout_on) {
// initialize Eposlhc
- static bool initialized = false;
- if (!initialized) {
+ if (!isInitialized_) {
+ isInitialized_ = true;
if (dataPath == "") {
data_path_ = (std::string(corsika_data("EPOS").c_str()) + "/").c_str();
}
initialize();
- initialized = true;
}
setParticlesStable();
}
@@ -136,9 +135,9 @@ namespace corsika::epos {
strcpy(::epos::fname_.fncs, CS.data);
::epos::nfname_.nfncs = CS.length;
- // dummy event (prepare commons)
- initializeEventLab(Code::Iron, Iron::nucleus_A, Iron::nucleus_Z, Code::Argon,
- Argon::nucleus_A, Argon::nucleus_Z, 100_GeV);
+ // initialiazes maximum energy and mass
+ initializeEventLab(Code::Lead, Lead::nucleus_A, Lead::nucleus_Z, Code::Lead,
+ Lead::nucleus_A, Lead::nucleus_Z, 1_TeV);
}
inline void InteractionModel::initializeEventCoM(Code const idBeam, int const iBeamA,
@@ -157,11 +156,8 @@ namespace corsika::epos {
::epos::enrgy_.ecms = Ecm / 1_GeV; // -> c.m.s. frame
- CORSIKA_LOGGER_TRACE(logger_,
- "inside EPOS: "
- "Ecm={}, "
- "Elab={}",
- ::epos::enrgy_.ecms, ::epos::enrgy_.elab);
+ CORSIKA_LOGGER_TRACE(logger_, "inside EPOS: Ecm={}, Elab={}", ::epos::enrgy_.ecms,
+ ::epos::enrgy_.elab);
configureParticles(idBeam, iBeamA, iBeamZ, idTarget, iTargetA, iTargetZ);
::epos::ainit_();
@@ -184,11 +180,7 @@ namespace corsika::epos {
// hadron-nucleon momentum
::epos::hadr1_.pnll = float(Plab / 1_GeV); // -> lab frame
- CORSIKA_LOGGER_TRACE(logger_,
- "inside EPOS: "
- "Ecm={}, "
- "Elab={}, "
- "Pnll={}",
+ CORSIKA_LOGGER_TRACE(logger_, "inside EPOS: Ecm={}, Elab={}, Pnll={}",
::epos::enrgy_.ecms, ::epos::enrgy_.elab, ::epos::hadr1_.pnll);
configureParticles(idBeam, iBeamA, iBeamZ, idTarget, iTargetA, iTargetZ);
@@ -449,13 +441,12 @@ namespace corsika::epos {
// create event
int iarg = 1;
::epos::aepos_(iarg);
-
::epos::afinal_();
- if (epos_listing_) {
+ if (epos_listing_) { // LCOV_EXCL_START
char nam[9] = "EPOSLHC&";
::epos::alistf_(nam, 9);
- }
+ } // LCOV_EXCL_STOP
// NSTORE-part
diff --git a/corsika/detail/modules/sibyll/InteractionModel.inl b/corsika/detail/modules/sibyll/InteractionModel.inl
index 9ecb434bc62cb5a0bbccc6a13ce889032b4bdcda..99f078673c2d80df3c6251bd6d22eed25a077c4f 100644
--- a/corsika/detail/modules/sibyll/InteractionModel.inl
+++ b/corsika/detail/modules/sibyll/InteractionModel.inl
@@ -151,8 +151,9 @@ namespace corsika::sibyll {
HEPEnergyType Elab_final = 0_GeV, Ecm_final = 0_GeV;
for (auto& psib : ss) {
// abort on particles that have decayed in Sibyll. Should not happen!
- if (psib.hasDecayed())
+ if (psib.hasDecayed()) { // LCOV_EXCL_START
throw std::runtime_error("found particle that decayed in SIBYLL!");
+ } // LCOV_EXCL_STOP
// transform 4-momentum to lab. frame
// note that the momentum needs to be rotated back
diff --git a/corsika/detail/modules/sibyll/NuclearInteractionModel.inl b/corsika/detail/modules/sibyll/NuclearInteractionModel.inl
index c2bfecd92a79c33c4a108460ba0fd4d3d710cf2b..4a564c7e96bcc64e059dfec58b43eb3f5b8ed7ac 100644
--- a/corsika/detail/modules/sibyll/NuclearInteractionModel.inl
+++ b/corsika/detail/modules/sibyll/NuclearInteractionModel.inl
@@ -57,10 +57,10 @@ namespace corsika::sibyll {
inline void
NuclearInteractionModel<TEnvironment, TNucleonModel>::printCrossSectionTable(
Code const pCode) const {
- if (!hadronicInteraction_.isValid(Code::Proton, pCode, 100_GeV)) {
+ if (!hadronicInteraction_.isValid(Code::Proton, pCode, 100_GeV)) { // LCOV_EXCL_START
CORSIKA_LOG_ERROR("Invalid target type {} for hadron interaction model.", pCode);
return;
- }
+ } // LCOV_EXCL_STOP
int const k = targetComponentsIndex_.at(pCode);
Code const pNuclei[] = {Code::Helium, Code::Lithium7, Code::Oxygen,
@@ -200,6 +200,9 @@ namespace corsika::sibyll {
// model is only designed for projectile nuclei. Collisions are broken down into
// "nucleon-target" collisions.
+ if (!is_nucleus(projectileId)) {
+ throw std::runtime_error("Can only handle nuclear projectiles.");
+ }
size_t const projectileA = get_nucleus_A(projectileId);
// this is center-of-mass for projectile_nucleon - target
diff --git a/corsika/modules/epos/InteractionModel.hpp b/corsika/modules/epos/InteractionModel.hpp
index 5a2c4a2901844988cb6d696c0d390022ae794b7d..0f08a919b7cc640d20351235dfc2164b155dbe50 100644
--- a/corsika/modules/epos/InteractionModel.hpp
+++ b/corsika/modules/epos/InteractionModel.hpp
@@ -83,6 +83,8 @@ namespace corsika::epos {
void setParticlesStable() const;
private:
+ inline static bool isInitialized_ = false;
+
std::string data_path_;
unsigned int count_ = 0;
bool epos_listing_;
diff --git a/tests/modules/testEpos.cpp b/tests/modules/testEpos.cpp
index 4700eab1091cf545a887b5b19ed06ae7227d00b6..d0f13d94fa941ae6931b9874b43d208240d9b5ad 100644
--- a/tests/modules/testEpos.cpp
+++ b/tests/modules/testEpos.cpp
@@ -29,7 +29,7 @@
using namespace corsika;
using namespace corsika::epos;
-TEST_CASE("Epos", "module,process") {
+TEST_CASE("EposBasics", "module,process") {
logging::set_level(logging::level::trace);
@@ -71,6 +71,7 @@ TEST_CASE("Epos", "module,process") {
SECTION("epos mass") {
CHECK_FALSE(corsika::epos::getEposMass(Code::Electron) / 1_GeV == Approx(0));
+ CHECK_THROWS(corsika::epos::getEposMass(Code::Unknown));
}
/*
@@ -88,6 +89,7 @@ TEST_CASE("Epos", "module,process") {
CHECK(p == convert_from_PDG(getEposPDGId(p)));
}
}
+ CHECK_THROWS(getEposPDGId(Code::Oxygen));
}
}
@@ -118,16 +120,17 @@ auto sqs2elab(HEPEnergyType const sqs, HEPEnergyType const ma, HEPEnergyType con
return (sqs * sqs - ma * ma - mb * mb) / 2. / mb;
}
-TEST_CASE("EposInterface", "modules") {
+TEST_CASE("Epos", "modules") {
logging::set_level(logging::level::trace);
+ RNGManager<>::getInstance().registerRandomStream("epos");
+ InteractionModel model;
+
auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
auto const& cs = *csPtr;
[[maybe_unused]] auto const& env_dummy = env;
- RNGManager<>::getInstance().registerRandomStream("epos");
-
SECTION("InteractionInterface - random number") {
auto const rndm = ::epos::rangen_();
CHECK(rndm > 0);
@@ -136,8 +139,6 @@ TEST_CASE("EposInterface", "modules") {
SECTION("InteractionInterface - isValid") {
- InteractionModel model;
-
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));
@@ -147,8 +148,6 @@ TEST_CASE("EposInterface", "modules") {
SECTION("InteractionInterface - getCrossSectionInelEla") {
- InteractionModel model;
-
// hydrogen target == proton target == neutron target
auto const [xs_prod_pp, xs_ela_pp] = model.getCrossSectionInelEla(
Code::Proton, Code::Proton,
@@ -172,12 +171,19 @@ TEST_CASE("EposInterface", "modules") {
CHECK(xs_prod_pp == xs_prod_pn);
CHECK(xs_ela_pp == xs_ela_pHydrogen);
CHECK(xs_ela_pn == xs_ela_pHydrogen);
+
+ // invalid system
+ auto const [xs_prod_0, xs_ela_0] = model.getCrossSectionInelEla(
+ Code::Electron, Code::Electron,
+ {sqrt(static_pow<2>(100_GeV) + static_pow<2>(Electron::mass)),
+ {cs, 100_GeV, 0_GeV, 0_GeV}},
+ {Electron::mass, {cs, 0_GeV, 0_GeV, 0_GeV}});
+ CHECK(xs_prod_0 / 1_mb == Approx(0));
+ CHECK(xs_ela_0 / 1_mb == Approx(0));
}
SECTION("InteractionModelInterface - hadron cross sections") {
- InteractionModel model;
-
// p-p at 7TeV around 70mb according to LHC
auto const xs_prod = model.getCrossSection(
Code::Proton, Code::Proton,
@@ -211,8 +217,6 @@ TEST_CASE("EposInterface", "modules") {
SECTION("InteractionInterface - nuclear cross sections") {
- InteractionModel model;
-
auto const xs_prod = model.getCrossSection(
Code::Proton, Code::Oxygen,
{100_GeV,
@@ -228,35 +232,67 @@ TEST_CASE("EposInterface", "modules") {
{Oxygen::mass, {cs, 0_GeV, 0_GeV, 0_GeV}});
CHECK(xs_prod2 / 1_mb == Approx(1076.7).margin(3.1));
}
-}
-
-TEST_CASE("EposDoInteractionNucleus", "module") {
-
- logging::set_level(logging::level::trace);
-
- auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
- auto const& cs = *csPtr;
- [[maybe_unused]] auto const& env_dummy = env;
-
- RNGManager<>::getInstance().registerRandomStream("epos");
-
- InteractionModel model;
-
- SECTION("InteractionInterface - nuclear projectile") {
+ /*
+ 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")
+ {
HEPEnergyType const P0 = 10_TeV;
- Code const pid = get_nucleus_code(40, 20);
+ Code const pid = Code::Proton;
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
+ MomentumVector(cs, {P0, 0_eV, 0_eV}); // this is secret knowledge about
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))), {cs, P0, 0_GeV, 0_GeV}},
- {Oxygen::mass, {cs, 0_GeV, 0_GeV, 0_GeV}});
+ 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);
@@ -274,4 +310,4 @@ TEST_CASE("EposDoInteractionNucleus", "module") {
// CHECK(length / 1_g * 1_cm * 1_cm ==
// Approx(30).margin(20)); // this is no physics validation
}
-}
+}
\ No newline at end of file
diff --git a/tests/modules/testSibyll.cpp b/tests/modules/testSibyll.cpp
index 7cd05aa6ef07d0fb936f833cf6cdaf073dc15489..042dc72a873d657b3e7da682fdead5d15df67e86 100644
--- a/tests/modules/testSibyll.cpp
+++ b/tests/modules/testSibyll.cpp
@@ -100,27 +100,9 @@ auto sumMomentum(TStackView const& view, CoordinateSystemPtr const& vCS) {
return sum;
}
-/*
-calculate COM boost object assuming fixed target collision with projectile
-*/
-COMBoost getCOMboost(HEPEnergyType const& eProjectileLab,
- MomentumVector const& pProjectileLab,
- CoordinateSystemPtr const& cs) {
- // define target
- // for Sibyll is always a single nucleon
- // FOR NOW: target is always at rest
- auto const pTargetLab = MomentumVector(cs, 0_GeV, 0_GeV, 0_GeV);
- FourVector const P4projLab(eProjectileLab, pProjectileLab);
- // define target kinematics in lab frame
- // define boost to and from CoM frame
- // CoM frame definition in Sibyll projectile: +z
- COMBoost const boost(P4projLab, constants::nucleonMass);
- return boost;
-}
-
TEST_CASE("SibyllInterface", "modules") {
- logging::set_level(logging::level::info);
+ logging::set_level(logging::level::trace);
// the environment and stack should eventually disappear from here
auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
@@ -168,6 +150,11 @@ TEST_CASE("SibyllInterface", "modules") {
CHECK(xs_ela_pp == xs_ela_pHydrogen);
CHECK(xs_ela_pn == xs_ela_pHydrogen);
+ // invalids
+ auto const xs_prod_0 = model.getCrossSection(Code::Electron, Code::Proton, aP4, bP4);
+ CHECK(xs_prod_0 / 1_mb == Approx(0));
+ CHECK_THROWS(model.doInteraction(view, Code::Electron, Code::Proton, aP4, bP4));
+
CHECK_THROWS(convertFromSibyll(corsika::sibyll::SibyllCode::Unknown));
}
@@ -274,8 +261,18 @@ TEST_CASE("SibyllInterface", "modules") {
MomentumVector(cs, {0_eV, 0_eV, 0_eV}));
model.doInteraction(view, pid, Code::Oxygen, P4, targetP4);
CrossSectionType const cx = model.getCrossSection(pid, Code::Oxygen, P4, targetP4);
- CHECK(cx / 1_mb == Approx(1250).margin(100)); // this is not physics validation
- CHECK(view.getSize() == Approx(40).margin(30)); // this is not physics validation
+ CHECK(cx / 1_mb == Approx(1250).margin(100)); // this is not physics validation
+ CHECK(view.getSize() == Approx(150).margin(140)); // this is not physics validation
+
+ // invalid to underlying model
+ FourMomentum P4mu(
+ 100_GeV,
+ {cs, {sqrt(static_pow<2>(100_GeV) - static_pow<2>(MuPlus::mass)), 0_eV, 0_eV}});
+ CrossSectionType const cx0 =
+ model.getCrossSection(Code::MuPlus, Code::Oxygen, P4mu, targetP4);
+ CHECK(cx0 / 1_mb == Approx(0));
+
+ CHECK_THROWS(model.doInteraction(view, Code::MuPlus, Code::Oxygen, P4mu, targetP4));
}
}
diff --git a/tests/modules/testUrQMD.cpp b/tests/modules/testUrQMD.cpp
index 3b3da197bf501a9fdf73b11fcf851c39fe953182..f75d59b75b44a41d25376d013c9d9e0d813efa32 100644
--- a/tests/modules/testUrQMD.cpp
+++ b/tests/modules/testUrQMD.cpp
@@ -89,6 +89,7 @@ TEST_CASE("UrQMD") {
}
SECTION("cross sections") {
+
FourMomentum const targetP4{Nitrogen::mass, {rootCS, {0_eV, 0_eV, 0_eV}}};
HEPMomentumType const P0 = 100_GeV;
@@ -108,6 +109,9 @@ TEST_CASE("UrQMD") {
CORSIKA_LOG_INFO("UrQMD cross seciton for {} is {} mb", code, cx / 1_mb);
CHECK(cx / 1_mb == Approx(checkCX[i++] / 1_mb).margin(1));
}
+
+ // invalid
+ CHECK_THROWS(urqmd.getTabulatedCrossSection(Code::Proton, Code::Proton, 100_GeV));
}
SECTION("pion+ projectile") {