diff --git a/Processes/Sibyll/Decay.cc b/Processes/Sibyll/Decay.cc index 09ac820cf43c487f0c415e70784542495d814d87..d39aff32b8296b9fc5664e7fc08bffb9902846ef 100644 --- a/Processes/Sibyll/Decay.cc +++ b/Processes/Sibyll/Decay.cc @@ -23,6 +23,7 @@ using std::vector; using namespace corsika; using namespace corsika::setup; + using SetupProjectile = corsika::setup::StackView::ParticleType; using SetupParticle = corsika::setup::Stack::ParticleType; @@ -149,19 +150,26 @@ namespace corsika::process::sibyll { int print_unit = 6; sib_list_(print_unit); + int nSecondaries = 0; + particles::Code pCodeSecondary; // copy particles from sibyll stack to corsika for (auto& psib : ss) { // FOR NOW: skip particles that have decayed in Sibyll, move to iterator? if (psib.HasDecayed()) continue; + nSecondaries++; + pCodeSecondary = process::sibyll::ConvertFromSibyll(psib.GetPID()); // add to corsika stack vP.AddSecondary( tuple<particles::Code, units::si::HEPEnergyType, corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{ - process::sibyll::ConvertFromSibyll(psib.GetPID()), psib.GetEnergy(), - psib.GetMomentum(), decayPoint, t0}); + pCodeSecondary, psib.GetEnergy(), psib.GetMomentum(), decayPoint, t0}); } // empty sibyll stack ss.Clear(); + + // check if particle decays into itself + if (1 == nSecondaries && pCode == pCodeSecondary) + throw std::runtime_error("Sibyll::Decay: Particle decays into itself!"); } } // namespace corsika::process::sibyll diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc index 51d438b193c5847b7ac3fc00c964d413e54407c1..f4f8c77d8a80c990dad2a95b6550f2719a52a090 100644 --- a/Processes/Sibyll/testSibyll.cc +++ b/Processes/Sibyll/testSibyll.cc @@ -158,7 +158,7 @@ TEST_CASE("SibyllInterface", "[processes]") { auto particle = stack.AddParticle( std::tuple<particles::Code, units::si::HEPEnergyType, corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{ - particles::Code::Proton, E0, plab, pos, 0_ns}); + particles::Code::Lambda0, E0, plab, pos, 0_ns}); corsika::stack::SecondaryView view(particle); auto projectile = view.GetProjectile();