From 4bc50090adcbc60b9d1de9978e2e2f5ee905f9c1 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Sat, 1 Jun 2019 11:21:29 +0100
Subject: [PATCH] added exception for particle decaying into itself
---
Processes/Sibyll/Decay.cc | 12 ++++++++++--
Processes/Sibyll/testSibyll.cc | 2 +-
2 files changed, 11 insertions(+), 3 deletions(-)
diff --git a/Processes/Sibyll/Decay.cc b/Processes/Sibyll/Decay.cc
index 09ac820cf..d39aff32b 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 51d438b19..f4f8c77d8 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();
--
GitLab