diff --git a/Framework/Particles/ParticleProperties.h b/Framework/Particles/ParticleProperties.h index ed0e919b5acce89a7b0573489695a635a0abb104..dbac43ac444f7411d63827655fd233a1b828ebf0 100644 --- a/Framework/Particles/ParticleProperties.h +++ b/Framework/Particles/ParticleProperties.h @@ -120,14 +120,22 @@ namespace corsika::particles { } int constexpr GetNucleusA(Code const p) { + if (p == Code::Nucleus) { + throw std::runtime_error("GetNucleusA(Code::Nucleus) is impossible!"); + } return detail::nucleusA[static_cast<CodeIntType>(p)]; } int constexpr GetNucleusZ(Code const p) { + if (p == Code::Nucleus) { + throw std::runtime_error("GetNucleusZ(Code::Nucleus) is impossible!"); + } return detail::nucleusZ[static_cast<CodeIntType>(p)]; } - bool constexpr IsNucleus(Code const p) { return GetNucleusA(p) != 0; } + bool constexpr IsNucleus(Code const p) { + return (p == Code::Nucleus) || (GetNucleusA(p) != 0); + } /** * the output operator for humand-readable particle codes diff --git a/Framework/Particles/testParticles.cc b/Framework/Particles/testParticles.cc index c06ab0d11c6ba087b16fef50ef405bc4dd980583..d4021c99ce53a85ddce38dec60f830cb69daa000 100644 --- a/Framework/Particles/testParticles.cc +++ b/Framework/Particles/testParticles.cc @@ -102,6 +102,7 @@ TEST_CASE("ParticleProperties", "[Particles]") { REQUIRE(IsHadron(Code::Proton)); REQUIRE(IsHadron(Code::PiPlus)); REQUIRE(IsHadron(Code::Oxygen)); + REQUIRE(IsHadron(Code::Nucleus)); } SECTION("Particle groups: muons") { @@ -124,7 +125,6 @@ TEST_CASE("ParticleProperties", "[Particles]") { REQUIRE_FALSE(IsNeutrino(Code::Oxygen)); } - SECTION("Nuclei") { REQUIRE_FALSE(IsNucleus(Code::Gamma)); REQUIRE(IsNucleus(Code::Argon)); @@ -137,5 +137,8 @@ TEST_CASE("ParticleProperties", "[Particles]") { REQUIRE(GetNucleusA(Code::Tritium) == 3); REQUIRE(Hydrogen::GetNucleusZ() == 1); REQUIRE(Tritium::GetNucleusA() == 3); + + REQUIRE_THROWS(GetNucleusA(Code::Nucleus)); + REQUIRE_THROWS(GetNucleusZ(Code::Nucleus)); } } diff --git a/Processes/Sibyll/NuclearInteraction.cc b/Processes/Sibyll/NuclearInteraction.cc index be970455fd33701667bc375d3ffe7581c7997786..384330c636a45a3dbeed8fc11ba2193e9dd1fd31 100644 --- a/Processes/Sibyll/NuclearInteraction.cc +++ b/Processes/Sibyll/NuclearInteraction.cc @@ -227,15 +227,18 @@ namespace corsika::process::sibyll { RootCoordinateSystem::GetInstance().GetRootCoordinateSystem(); const particles::Code corsikaBeamId = vP.GetPID(); - if (!particles::IsNucleus(corsikaBeamId)) { - // no nuclear interaction - return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm); + + if (corsikaBeamId != particles::Code::Nucleus) { + // check if target-style nucleus (enum), these are not allowed as projectile + if (particles::IsNucleus(corsikaBeamId)) + throw std::runtime_error( + "NuclearInteraction: GetInteractionLength: Wrong nucleus type. Nuclear " + "projectiles should use NuclearStackExtension!"); + else { + // no nuclear interaction + return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm); + } } - // check if target-style nucleus (enum) - if (corsikaBeamId != particles::Code::Nucleus) - throw std::runtime_error( - "NuclearInteraction: GetInteractionLength: Wrong nucleus type. Nuclear " - "projectiles should use NuclearStackExtension!"); // read from cross section code table @@ -333,13 +336,6 @@ namespace corsika::process::sibyll { // const auto ProjMass = vP.GetMass(); cout << "NuclearInteraction: DoInteraction: called with:" << ProjId << endl; - if (!IsNucleus(ProjId)) { - cout << "WARNING: non nuclear projectile in NUCLIB!" << endl; - // this should not happen - // throw std::runtime_error("Non nuclear projectile in NUCLIB!"); - return process::EProcessReturn::eOk; - } - // check if target-style nucleus (enum) if (ProjId != particles::Code::Nucleus) throw std::runtime_error(