diff --git a/corsika/detail/modules/sibyll/Interaction.inl b/corsika/detail/modules/sibyll/Interaction.inl
index bbcacf4e5e7e90fca3fa5a21a0aa639848ccdc3c..39d240df8c513ec007302ff146bef00bd56cfff7 100644
--- a/corsika/detail/modules/sibyll/Interaction.inl
+++ b/corsika/detail/modules/sibyll/Interaction.inl
@@ -44,34 +44,6 @@ namespace corsika::sibyll {
CORSIKA_LOG_DEBUG("Sibyll::Interaction n={}, Nnuc={}", count_, nucCount_);
}
- inline void Interaction::setStable(std::vector<corsika::Code> const& vParticleList) {
- for (auto p : vParticleList) Interaction::setStable(p);
- }
-
- inline void Interaction::setUnstable(std::vector<corsika::Code> const& vParticleList) {
- for (auto p : vParticleList) Interaction::setUnstable(p);
- }
-
- inline void Interaction::setUnstable(const corsika::Code vCode) {
- CORSIKA_LOG_DEBUG("Sibyll::Interaction: setting {} unstable..", vCode);
- const int s_id = abs(corsika::sibyll::convertToSibyllRaw(vCode));
- s_csydec_.idb[s_id - 1] = abs(s_csydec_.idb[s_id - 1]);
- }
-
- inline void Interaction::setStable(const corsika::Code vCode) {
- CORSIKA_LOG_DEBUG("Sibyll::Interaction: setting {} stable..", vCode);
- const int s_id = abs(corsika::sibyll::convertToSibyllRaw(vCode));
- s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
- }
-
- inline void Interaction::setAllUnstable() {
- for (int i = 0; i < 99; ++i) s_csydec_.idb[i] = abs(s_csydec_.idb[i]);
- }
-
- inline void Interaction::setAllStable() {
- for (int i = 0; i < 99; ++i) s_csydec_.idb[i] = -1 * abs(s_csydec_.idb[i]);
- }
-
inline std::tuple<corsika::CrossSectionType, corsika::CrossSectionType>
Interaction::getCrossSection(const corsika::Code BeamId, const corsika::Code TargetId,
const corsika::HEPEnergyType CoMenergy) const {
diff --git a/corsika/modules/sibyll/Interaction.hpp b/corsika/modules/sibyll/Interaction.hpp
index 605e6c1843bd8be4ef481e8f4698ea1bd3610c62..81f0aa68571101dec0a48bb0aae72055882980f2 100644
--- a/corsika/modules/sibyll/Interaction.hpp
+++ b/corsika/modules/sibyll/Interaction.hpp
@@ -52,14 +52,6 @@ namespace corsika::sibyll {
void doInteraction(TSecondaries&);
private:
- void setStable(std::vector<Code> const&);
- void setUnstable(std::vector<Code> const&);
-
- void setUnstable(Code const);
- void setStable(Code const);
- void setAllUnstable();
- void setAllStable();
-
int getMaxTargetMassNumber() const { return maxTargetMassNumber_; }
HEPEnergyType getMinEnergyCoM() const { return minEnergyCoM_; }
HEPEnergyType getMaxEnergyCoM() const { return maxEnergyCoM_; }
diff --git a/tests/modules/testSibyll.cpp b/tests/modules/testSibyll.cpp
index 85b512740c1471e9a83d40947330ddbb7694dcc0..f6a517e1a8822b36a30b743cfc15ce1b9883305f 100644
--- a/tests/modules/testSibyll.cpp
+++ b/tests/modules/testSibyll.cpp
@@ -68,8 +68,11 @@ TEST_CASE("Sibyll", "[processes]") {
}
SECTION("sibyll mass") {
-
CHECK_FALSE(corsika::sibyll::getSibyllMass(Code::Electron) == 0_GeV);
+ // Nucleus not a particle
+ CHECK_THROWS(corsika::sibyll::getSibyllMass(Code::Nucleus));
+ // Higgs not a particle in Sibyll
+ CHECK_THROWS(corsika::sibyll::getSibyllMass(Code::H0));
}
}
@@ -129,6 +132,16 @@ TEST_CASE("SibyllInterface", "[processes]") {
CHECK(xs_prod_pp == xs_prod_pn);
CHECK(xs_ela_pp == xs_ela_pHydrogen);
CHECK(xs_ela_pn == xs_ela_pHydrogen);
+
+ // out of range
+ // beam particle
+ CHECK_THROWS(
+ std::get<0>(model.getCrossSection(Code::Electron, Code::Hydrogen, 100_GeV)));
+ // target particle
+ CHECK(std::get<0>(model.getCrossSection(Code::Proton, Code::Electron, 100_GeV)) ==
+ std::numeric_limits<double>::infinity() * 1_mb);
+ // energy out of range
+ CHECK_THROWS(std::get<0>(model.getCrossSection(Code::Proton, Code::Hydrogen, 5_GeV)));
}
SECTION("InteractionInterface - low energy") {
@@ -142,7 +155,9 @@ TEST_CASE("SibyllInterface", "[processes]") {
auto particle = stack->first();
- Interaction model;
+ // also print particles after sibyll was called
+ Interaction model(true);
+
model.doInteraction(view);
auto const pSum = sumMomentum(view, cs);
@@ -209,13 +224,59 @@ TEST_CASE("SibyllInterface", "[processes]") {
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);
+ [[maybe_unused]] GrammageType const length = model.getInteractionLength(particle);
CHECK(length / 1_g * 1_cm * 1_cm == Approx(88.7).margin(0.1));
// CHECK(view.getEntries() == 9); //! \todo: this was 20 before refactory-2020: check
// "also sibyll not stable wrt. to compiler
// changes"
}
+ SECTION("InteractionInterface - energy too low") {
+
+ const HEPEnergyType P0 = 5_GeV;
+ auto [stack, viewPtr] = setup::testing::setup_stack(
+ Code::Proton, 0, 0, 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;
+
+ auto particle = stack->first();
+
+ Interaction model;
+ CHECK_THROWS(model.doInteraction(view));
+
+ [[maybe_unused]] GrammageType const length = model.getInteractionLength(particle);
+ CHECK(model.getInteractionLength(particle) / 1_g * 1_cm * 1_cm ==
+ std::numeric_limits<double>::infinity());
+ }
+
+ SECTION("InteractionInterface - energy too high") {
+
+ const HEPEnergyType P0 = 1000_EeV;
+ auto [stack, viewPtr] = setup::testing::setup_stack(
+ Code::Proton, 0, 0, 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;
+
+ Interaction model;
+ CHECK_THROWS(model.doInteraction(view));
+ }
+
+ SECTION("InteractionInterface - target nucleus out of range") {
+ auto [env1, csPtr1, nodePtr1] = setup::testing::setup_environment(Code::Argon);
+ auto const& cs1 = *csPtr1;
+ const HEPEnergyType P0 = 150_GeV;
+ auto [stack, viewPtr] = setup::testing::setup_stack(
+ Code::Electron, 0, 0, P0, (setup::Environment::BaseNodeType* const)nodePtr1, cs1);
+ MomentumVector plab = MomentumVector(
+ cs1, {P0, 0_eV, 0_eV}); // this is secret knowledge about setupStack
+ setup::StackView& view = *viewPtr;
+
+ Interaction model;
+ CHECK_THROWS(model.doInteraction(view));
+ }
+
SECTION("NuclearInteractionInterface") {
auto [stack, viewPtr] =
@@ -248,16 +309,32 @@ TEST_CASE("SibyllInterface", "[processes]") {
Decay model;
model.printDecayConfig();
[[maybe_unused]] const TimeType time = model.getLifetime(particle);
-
+ auto const gamma = particle.getEnergy() / particle.getMass();
+ CHECK(time == get_lifetime(Code::Lambda0) * gamma);
model.doDecay(view);
// run checks
// lambda decays into proton and pi- or neutron and pi+
CHECK(stack.getEntries() == 3);
}
+ SECTION("DecayInterface - decay not handled") {
+ // sibyll does not know the higgs for example
+ auto [stackPtr, viewPtr] = setup::testing::setup_stack(
+ Code::H0, 0, 0, 10_GeV, (setup::Environment::BaseNodeType* const)nodePtr, cs);
+ setup::StackView& view = *viewPtr;
+ auto& stack = *stackPtr;
+ auto particle = stack.first();
+
+ Decay model;
+
+ CHECK(model.getLifetime(particle) == std::numeric_limits<double>::infinity() * 1_s);
+ CHECK_THROWS(model.doDecay(view));
+ }
+
SECTION("DecayConfiguration") {
Decay model({Code::PiPlus, Code::PiMinus});
+ model.printDecayConfig();
CHECK(model.isDecayHandled(Code::PiPlus));
CHECK(model.isDecayHandled(Code::PiMinus));
CHECK_FALSE(model.isDecayHandled(Code::KPlus));
@@ -269,10 +346,13 @@ TEST_CASE("SibyllInterface", "[processes]") {
model.setHandleDecay(particleTestList);
for (auto& pCode : particleTestList) CHECK(model.isDecayHandled(pCode));
- // individually
+ // set decay individually
model.setHandleDecay(Code::KMinus);
+ // fail
+ CHECK_THROWS(model.setHandleDecay(Code::H0));
// possible decays
+ CHECK_FALSE(model.canHandleDecay(Code::H0));
CHECK_FALSE(model.canHandleDecay(Code::Proton));
CHECK_FALSE(model.canHandleDecay(Code::Electron));
CHECK(model.canHandleDecay(Code::PiPlus));