diff --git a/Processes/Sibyll/Decay.cc b/Processes/Sibyll/Decay.cc index d39aff32b8296b9fc5664e7fc08bffb9902846ef..f7f0c1f360279bb4e968da747dd59d53cbe79b9f 100644 --- a/Processes/Sibyll/Decay.cc +++ b/Processes/Sibyll/Decay.cc @@ -24,6 +24,7 @@ using std::vector; using namespace corsika; using namespace corsika::setup; +using SetupView = corsika::setup::StackView; using SetupProjectile = corsika::setup::StackView::ParticleType; using SetupParticle = corsika::setup::Stack::ParticleType; @@ -150,26 +151,36 @@ 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>{ - pCodeSecondary, psib.GetEnergy(), psib.GetMomentum(), decayPoint, t0}); + process::sibyll::ConvertFromSibyll(psib.GetPID()), psib.GetEnergy(), + psib.GetMomentum(), decayPoint, t0}); } // empty sibyll stack ss.Clear(); + } - // check if particle decays into itself - if (1 == nSecondaries && pCode == pCodeSecondary) + template <> + EProcessReturn Decay::DoSecondaries(SetupView& vS) { // corsika::setup::StackView&vS){} + auto pCode = vS.GetProjectile().GetPID(); + if (vS.GetSize() == 1 && pCode == vS.GetNextParticle().GetPID()) throw std::runtime_error("Sibyll::Decay: Particle decays into itself!"); - } + /* + here we could also post-process the decay products and let short-lived resonances + decay, etc + + See Issue 196 + + https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/issues/196 + + */ + return EProcessReturn::eOk; + } } // namespace corsika::process::sibyll diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h index 8c630de8df90b2aaa304ebbd592c96cc5f338994..745a83072814f60fb017703ab2a8c286eaa11b7e 100644 --- a/Processes/Sibyll/Decay.h +++ b/Processes/Sibyll/Decay.h @@ -50,6 +50,9 @@ namespace corsika::process { template <typename TProjectile> void DoDecay(TProjectile&); + + template <typename TParticleView> + EProcessReturn DoSecondaries(TParticleView&); }; } // namespace sibyll } // namespace corsika::process diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc index f4f8c77d8a80c990dad2a95b6550f2719a52a090..3a0f365d587a95680f8464a6dc3464f25bbafc40 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::Lambda0, E0, plab, pos, 0_ns}); + particles::Code::Proton, E0, plab, pos, 0_ns}); corsika::stack::SecondaryView view(particle); auto projectile = view.GetProjectile(); @@ -166,6 +166,8 @@ TEST_CASE("SibyllInterface", "[processes]") { model.Init(); /*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoDecay(projectile); + // run checks + model.DoSecondaries(view); [[maybe_unused]] const TimeType time = model.GetLifetime(particle); }