diff --git a/Documentation/Examples/boundary_example.cc b/Documentation/Examples/boundary_example.cc
index c8532185a0484c1df13a353fc15b4362b6793840..ff4885b378251e0e93e3137c5357dc0f1858d211 100644
--- a/Documentation/Examples/boundary_example.cc
+++ b/Documentation/Examples/boundary_example.cc
@@ -121,9 +121,7 @@ int main() {
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
process::sibyll::Interaction sibyll;
- process::sibyll::Decay decay{{particles::Code::PiPlus, particles::Code::PiMinus,
- particles::Code::KPlus, particles::Code::KMinus,
- particles::Code::K0Long, particles::Code::K0Short}};
+ process::sibyll::Decay decay;
process::particle_cut::ParticleCut cut(20_GeV);
@@ -180,4 +178,8 @@ int main() {
cut.GetCutEnergy() + cut.GetInvEnergy() + cut.GetEmEnergy();
cout << "total energy (GeV): " << Efinal / 1_GeV << endl
<< "relative difference (%): " << (Efinal / E0 - 1.) * 100 << endl;
+
+ // basic check for unit-tests
+ assert(cut.GetNumberEmParticles() == 29785);
+ assert(cut.GetNumberInvParticles() == 26697);
}
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index dc08e7167e9a58309f4ad5296b0173d856a21c45..c2cc3c4f49f0ce60b089a9e3b893151e6b7464e8 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -130,15 +130,11 @@ int main() {
tracking_line::TrackingLine tracking;
stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0);
- const std::vector<particles::Code> trackedHadrons = {
- particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
- particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
-
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
random::RNGManager::GetInstance().RegisterRandomStream("pythia");
process::sibyll::Interaction sibyll;
process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
- process::sibyll::Decay decay(trackedHadrons);
+ process::sibyll::Decay decay;
process::particle_cut::ParticleCut cut(20_GeV);
process::track_writer::TrackWriter trackWriter("tracks.dat");
@@ -162,4 +158,8 @@ int main() {
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
cout << "total dEdX energy (GeV): " << eLoss.GetTotal() / 1_GeV << endl
<< "relative difference (%): " << eLoss.GetTotal() / E0 * 100 << endl;
+
+ // basic check for unit-tests
+ assert(cut.GetNumberEmParticles() == 127);
+ assert(cut.GetNumberInvParticles() == 116);
}
diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index e2c27ab30f68623703f90b59fb51616482e757e6..7745c16289848a8fe8d58c122ae2f607b1bace7c 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -125,14 +125,10 @@ int main() {
tracking_line::TrackingLine tracking;
stack_inspector::StackInspector<setup::Stack> stackInspect(1, true, E0);
- const std::vector<particles::Code> trackedHadrons = {
- particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
- particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
-
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
process::sibyll::Interaction sibyll;
process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
- process::sibyll::Decay decay(trackedHadrons);
+ process::sibyll::Decay decay;
process::particle_cut::ParticleCut cut(20_GeV);
process::track_writer::TrackWriter trackWriter("tracks.dat");
@@ -155,4 +151,8 @@ int main() {
<< "relative difference (%): " << (Efinal / E0 - 1) * 100 << endl;
cout << "total dEdX energy (GeV): " << eLoss.GetTotal() / 1_GeV << endl
<< "relative difference (%): " << eLoss.GetTotal() / E0 * 100 << endl;
+
+ // basic check for unit-tests
+ assert(cut.GetNumberEmParticles() == 526);
+ assert(cut.GetNumberInvParticles() == 645);
}
diff --git a/Framework/Cascade/Cascade.h b/Framework/Cascade/Cascade.h
index 873e47817c20ade54e96c517df4662857a1247a0..99c28b6c97819a7d6c2a92bb21590aa9244fcae6 100644
--- a/Framework/Cascade/Cascade.h
+++ b/Framework/Cascade/Cascade.h
@@ -263,6 +263,10 @@ namespace corsika::cascade {
InverseTimeType inv_decay_count = 0 / second;
fProcessSequence.SelectDecay(vParticle, projectile, sample_process,
inv_decay_count);
+ // make sure particle actually did decay if it should have done so
+ if (secondaries.GetSize() == 1 &&
+ projectile.GetPID() == secondaries.GetNextParticle().GetPID())
+ throw std::runtime_error("Cascade::Step: Particle decays into itself!");
}
fProcessSequence.DoSecondaries(secondaries);
@@ -273,6 +277,7 @@ namespace corsika::cascade {
// be "protected" and not accessible to physics
} else { // step-length limitation within volume
+
std::cout << "step-length limitation" << std::endl;
fProcessSequence.DoSecondaries(secondaries);
}
diff --git a/Processes/ParticleCut/ParticleCut.h b/Processes/ParticleCut/ParticleCut.h
index 6d5f7a685cccbce7f2e079a3c291f90dbf20039b..739a1919ebc20a8c080759baedf4eae5cb942e75 100644
--- a/Processes/ParticleCut/ParticleCut.h
+++ b/Processes/ParticleCut/ParticleCut.h
@@ -24,9 +24,9 @@ namespace corsika::process {
units::si::HEPEnergyType fEnergy = 0 * units::si::electronvolt;
units::si::HEPEnergyType fEmEnergy = 0 * units::si::electronvolt;
- int fEmCount = 0;
+ unsigned int fEmCount = 0;
units::si::HEPEnergyType fInvEnergy = 0 * units::si::electronvolt;
- int fInvCount = 0;
+ unsigned int fInvCount = 0;
public:
ParticleCut(const units::si::HEPEnergyType vCut)
@@ -46,6 +46,8 @@ namespace corsika::process {
units::si::HEPEnergyType GetInvEnergy() const { return fInvEnergy; }
units::si::HEPEnergyType GetCutEnergy() const { return fEnergy; }
units::si::HEPEnergyType GetEmEnergy() const { return fEmEnergy; }
+ unsigned int GetNumberEmParticles() const { return fEmCount; }
+ unsigned int GetNumberInvParticles() const { return fInvCount; }
};
} // namespace particle_cut
} // namespace corsika::process
diff --git a/Processes/Sibyll/Decay.cc b/Processes/Sibyll/Decay.cc
index b2ac4a9efebeea382157dc3af3329d8ad878b4c4..ea50d4477715e7d765abb000087fb56ef6f3e722 100644
--- a/Processes/Sibyll/Decay.cc
+++ b/Processes/Sibyll/Decay.cc
@@ -23,82 +23,69 @@ 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;
namespace corsika::process::sibyll {
- Decay::Decay(vector<particles::Code> pParticles)
- : fTrackedParticles(pParticles) {}
+ Decay::Decay() {}
Decay::~Decay() { cout << "Sibyll::Decay n=" << fCount << endl; }
void Decay::Init() {
- SetHadronsUnstable();
- SetParticleListStable(fTrackedParticles);
+ // switch off decays to avoid internal decay chains
+ SetAllStable();
+ }
+
+ void Decay::SetStable(const vector<particles::Code> vParticleList) {
+ for (auto p : vParticleList) Decay::SetStable(p);
+ }
+
+ void Decay::SetUnstable(const vector<particles::Code> vParticleList) {
+ for (auto p : vParticleList) Decay::SetUnstable(p);
+ }
+
+ bool Decay::IsStable(const particles::Code vCode) {
+ return abs(process::sibyll::ConvertToSibyllRaw(vCode)) <= 0 ? true : false;
+ }
+
+ bool Decay::IsUnstable(const particles::Code vCode) {
+ return abs(process::sibyll::ConvertToSibyllRaw(vCode)) > 0 ? true : false;
}
- void Decay::SetParticleListStable(const vector<particles::Code> particleList) {
- /*
- Sibyll is hadronic generator
- only hadrons decay
- */
- // set particles unstable
- SetHadronsUnstable();
- cout << "Interaction: setting tracked hadrons stable.." << endl;
- for (auto p : particleList) Decay::SetStable(p);
+ void Decay::SetDecay(const particles::Code vCode, const bool vMakeUnstable) {
+ vMakeUnstable ? SetUnstable(vCode) : SetStable(vCode);
}
- void Decay::SetUnstable(const particles::Code pCode) {
- int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
+ void Decay::SetUnstable(const particles::Code vCode) {
+ cout << "Sibyll::Interaction: setting " << vCode << " unstable.." << endl;
+ const int s_id = abs(process::sibyll::ConvertToSibyllRaw(vCode));
s_csydec_.idb[s_id - 1] = abs(s_csydec_.idb[s_id - 1]);
}
- void Decay::SetStable(const particles::Code pCode) {
- int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
+ void Decay::SetStable(const particles::Code vCode) {
+ cout << "Sibyll::Interaction: setting " << vCode << " stable.." << endl;
+ const int s_id = abs(process::sibyll::ConvertToSibyllRaw(vCode));
s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
}
void Decay::SetAllStable() {
- // name? also makes EM particles stable
-
- cout << "Decay: setting all particles stable.." << endl;
-
- // loop over all particles in sibyll
- // should be changed to loop over human readable list
- // i.e. particles::ListOfParticles()
- for (auto& p : corsika2sibyll) {
- // cout << (int)p << endl;
- const int sibCode = static_cast<int>(p);
- // skip unknown and antiparticles
- if (sibCode < 1) continue;
- s_csydec_.idb[sibCode - 1] = -1 * abs(s_csydec_.idb[sibCode - 1]);
- }
+ for (int i = 0; i < 99; ++i) s_csydec_.idb[i] = -1 * abs(s_csydec_.idb[i]);
}
- void Decay::SetHadronsUnstable() {
-
- // name? also makes EM particles stable
+ void Decay::SetAllUnstable() {
+ for (int i = 0; i < 99; ++i) s_csydec_.idb[i] = abs(s_csydec_.idb[i]);
+ }
- // loop over all particles in sibyll
- // should be changed to loop over human readable list
- // i.e. particles::ListOfParticles()
- cout << "Sibyll: setting hadrons unstable.." << endl;
- // make ALL particles unstable, then set EM stable
- for (int sibCode : corsika2sibyll) {
- if (sibCode < 1) continue;
- s_csydec_.idb[sibCode - 1] = abs(s_csydec_.idb[sibCode - 1]);
- }
- // set Leptons and Proton and Neutron stable
- // use stack to loop over particles
- constexpr particles::Code particleList[] = {
- particles::Code::Proton, particles::Code::Neutron, particles::Code::Electron,
- particles::Code::Positron, particles::Code::NuE, particles::Code::NuEBar,
- particles::Code::MuMinus, particles::Code::MuPlus, particles::Code::NuMu,
- particles::Code::NuMuBar};
-
- for (auto p : particleList) {
- const int sibid = process::sibyll::ConvertToSibyllRaw(p);
- s_csydec_.idb[sibid - 1] = (-1) * abs(s_csydec_.idb[sibid - 1]);
- }
+ void Decay::PrintDecayConfig(const particles::Code vCode) {
+ cout << "Decay: Sibyll decay configuration:" << endl;
+ const int sibCode = process::sibyll::ConvertToSibyllRaw(vCode);
+ const int absSibCode = abs(sibCode);
+ cout << vCode << " is ";
+ if (s_csydec_.idb[absSibCode - 1] <= 0)
+ cout << "stable" << endl;
+ else
+ cout << "unstable" << endl;
}
template <>
@@ -148,12 +135,16 @@ namespace corsika::process::sibyll {
// remember position
Point const decayPoint = vP.GetPosition();
TimeType const t0 = vP.GetTime();
- // set all particles/hadrons unstable
- // setHadronsUnstable();
+ // remember if particles is unstable
+ // auto const priorIsUnstable = IsUnstable(pCode);
+ // switch on decay for this particle
SetUnstable(pCode);
+ PrintDecayConfig(pCode);
+
// call sibyll decay
cout << "Decay: calling Sibyll decay routine.." << endl;
decsib_();
+
// reset to stable
SetStable(pCode);
// print output
diff --git a/Processes/Sibyll/Decay.h b/Processes/Sibyll/Decay.h
index 0416bf6e24c191ef68f6714c29844d8f13b9cab7..05bc0186ee6b2b4231d540cf8a1ccb209f65994d 100644
--- a/Processes/Sibyll/Decay.h
+++ b/Processes/Sibyll/Decay.h
@@ -13,6 +13,7 @@
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/DecayProcess.h>
+#include <corsika/process/SecondariesProcess.h>
#include <vector>
@@ -21,27 +22,42 @@ namespace corsika::process {
namespace sibyll {
class Decay : public corsika::process::DecayProcess<Decay> {
- std::vector<particles::Code> const fTrackedParticles;
int fCount = 0;
public:
- Decay(std::vector<particles::Code>);
+ Decay();
~Decay();
void Init();
- void SetParticleListStable(const std::vector<particles::Code>);
+ void SetStable(const std::vector<particles::Code>);
+ void SetUnstable(const std::vector<particles::Code>);
void SetUnstable(const corsika::particles::Code);
void SetStable(const corsika::particles::Code);
+ void SetAllUnstable();
void SetAllStable();
+ bool IsStable(const corsika::particles::Code);
+ bool IsUnstable(const corsika::particles::Code);
+ void SetDecay(const particles::Code, const bool);
+
+ void PrintDecayConfig(const corsika::particles::Code);
void SetHadronsUnstable();
template <typename TParticle>
corsika::units::si::TimeType GetLifetime(TParticle const&) const;
+ /**
+ In this function SIBYLL is called to produce to decay the input particle.
+ */
+
template <typename TProjectile>
void DoDecay(TProjectile&);
+
+ template <typename TParticleView>
+ EProcessReturn DoSecondaries(TParticleView&);
};
+
} // namespace sibyll
+
} // namespace corsika::process
#endif
diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc
index d55506688b9e6c8d607ee431171088f3c8237de7..def7ca1cc996c2629613d8bef2f76efe234ee6f2 100644
--- a/Processes/Sibyll/Interaction.cc
+++ b/Processes/Sibyll/Interaction.cc
@@ -47,46 +47,38 @@ namespace corsika::process::sibyll {
// initialize Sibyll
if (!fInitialized) {
sibyll_ini_();
-
- // any decays in sibyll? if yes need to define which particles
- if (fInternalDecays) {
- // define which particles are passed to corsika, i.e. which particles make it into
- // history even very shortlived particles like charm or pi0 are of interest here
- const std::vector<particles::Code> hadronsWeWantTrackedByCorsika = {
- particles::Code::PiPlus, particles::Code::PiMinus,
- particles::Code::Pi0, particles::Code::KMinus,
- particles::Code::KPlus, particles::Code::K0Long,
- particles::Code::K0Short, particles::Code::SigmaPlus,
- particles::Code::SigmaMinus, particles::Code::Lambda0,
- particles::Code::Xi0, particles::Code::XiMinus,
- particles::Code::OmegaMinus, particles::Code::DPlus,
- particles::Code::DMinus, particles::Code::D0,
- particles::Code::D0Bar};
-
- Interaction::SetParticleListStable(hadronsWeWantTrackedByCorsika);
- }
-
fInitialized = true;
}
}
- void Interaction::SetParticleListStable(
- std::vector<particles::Code> const& particleList) {
- for (auto p : particleList) Interaction::SetStable(p);
+ void Interaction::SetStable(std::vector<particles::Code> const& vParticleList) {
+ for (auto p : vParticleList) Interaction::SetStable(p);
+ }
+
+ void Interaction::SetUnstable(std::vector<particles::Code> const& vParticleList) {
+ for (auto p : vParticleList) Interaction::SetUnstable(p);
}
- void Interaction::SetUnstable(const particles::Code pCode) {
- cout << "Sibyll::Interaction: setting " << pCode << " unstable.." << endl;
- int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
+ void Interaction::SetUnstable(const particles::Code vCode) {
+ cout << "Sibyll::Interaction: setting " << vCode << " unstable.." << endl;
+ const int s_id = abs(process::sibyll::ConvertToSibyllRaw(vCode));
s_csydec_.idb[s_id - 1] = abs(s_csydec_.idb[s_id - 1]);
}
- void Interaction::SetStable(const particles::Code pCode) {
- cout << "Sibyll::Interaction: setting " << pCode << " stable.." << endl;
- int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
+ void Interaction::SetStable(const particles::Code vCode) {
+ cout << "Sibyll::Interaction: setting " << vCode << " stable.." << endl;
+ const int s_id = abs(process::sibyll::ConvertToSibyllRaw(vCode));
s_csydec_.idb[s_id - 1] = (-1) * abs(s_csydec_.idb[s_id - 1]);
}
+ void Interaction::SetAllStable() {
+ for (int i = 0; i < 99; ++i) s_csydec_.idb[i] = -1 * abs(s_csydec_.idb[i]);
+ }
+
+ void Interaction::SetAllUnstable() {
+ for (int i = 0; i < 99; ++i) s_csydec_.idb[i] = abs(s_csydec_.idb[i]);
+ }
+
tuple<units::si::CrossSectionType, units::si::CrossSectionType>
Interaction::GetCrossSection(const particles::Code BeamId,
const particles::Code TargetId,
@@ -325,8 +317,16 @@ namespace corsika::process::sibyll {
const double sqs = Ecm / 1_GeV;
// running sibyll, filling stack
sibyll_(kBeam, targetSibCode, sqs);
- // running decays
- if (fInternalDecays) decsib_();
+ if (fInternalDecays) {
+ // particles that decay internally will never appear on the corsika stack
+ // switch on all decays except for the particles we want to take part in the
+ // tracking
+ SetAllUnstable();
+ SetStable(fTrackedParticles);
+ decsib_();
+ // reset
+ SetAllStable();
+ }
// print final state
int print_unit = 6;
sib_list_(print_unit);
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index 0eb65029c8efeb5a909f517ea0dc385e206ae44a..d442de3cdb1256e9ebe332d108e52541122ffde7 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -31,9 +31,13 @@ namespace corsika::process::sibyll {
void Init();
- void SetParticleListStable(std::vector<particles::Code> const&);
+ void SetStable(std::vector<particles::Code> const&);
+ void SetUnstable(std::vector<particles::Code> const&);
+
void SetUnstable(const corsika::particles::Code);
void SetStable(const corsika::particles::Code);
+ void SetAllUnstable();
+ void SetAllStable();
bool WasInitialized() { return fInitialized; }
bool IsValidCoMEnergy(corsika::units::si::HEPEnergyType ecm) const {
@@ -65,7 +69,18 @@ namespace corsika::process::sibyll {
private:
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
-
+ // FOR NOW keep trackedParticles private, could be configurable
+ std::vector<particles::Code> const fTrackedParticles = {
+ particles::Code::PiPlus, particles::Code::PiMinus,
+ particles::Code::Pi0, particles::Code::KMinus,
+ particles::Code::KPlus, particles::Code::K0Long,
+ particles::Code::K0Short, particles::Code::SigmaPlus,
+ particles::Code::SigmaMinus, particles::Code::Lambda0,
+ particles::Code::Xi0, particles::Code::XiMinus,
+ particles::Code::OmegaMinus, particles::Code::DPlus,
+ particles::Code::DMinus, particles::Code::D0,
+ particles::Code::MuMinus, particles::Code::MuPlus,
+ particles::Code::D0Bar};
const bool fInternalDecays = true;
const corsika::units::si::HEPEnergyType fMinEnergyCoM =
10. * 1e9 * corsika::units::si::electronvolt;
diff --git a/Processes/Sibyll/sibyll2.3c.f b/Processes/Sibyll/sibyll2.3c.f
index e1f07491be0e017d84ac1f809ea00301615f4356..457b0e53bc18f0b74b4a42d3c359efd391359518 100644
--- a/Processes/Sibyll/sibyll2.3c.f
+++ b/Processes/Sibyll/sibyll2.3c.f
@@ -478,7 +478,7 @@ C-----------------------------------------------------------------------
CALL BLOCK_INI
CALL NUC_GEOM_INI
CALL SIG_AIR_INI
- CALL DEC_INI
+c CALL DEC_INI
c... charm frag. normalisation
CALL ZNORMAL
diff --git a/Processes/Sibyll/sibyll_codes.dat b/Processes/Sibyll/sibyll_codes.dat
index 08209d0f9f27249932c850199dc7cb5c12c92ef9..45aa507f964304a941b2ccf5c257f52e11f46c26 100644
--- a/Processes/Sibyll/sibyll_codes.dat
+++ b/Processes/Sibyll/sibyll_codes.dat
@@ -17,6 +17,8 @@ Pi0 6 1 2
# rho0 could interact but sibyll has no cross section/interaction length. was used for gamma had int
Rho0 27 0 0
K0Long 11 1 3
+K0 21 0 3
+K0Bar 22 0 3
PiPlus 7 1 2
PiMinus 8 1 2
RhoPlus 25 0 0
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index a736614e35ac17c2b2291d40ea33705aa9a8e4bb..b0099568795e074ce6bc72e2635eaa1a35c84344 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -69,6 +69,7 @@ TEST_CASE("Sibyll", "[processes]") {
#include <corsika/environment/Environment.h>
#include <corsika/environment/HomogeneousMedium.h>
#include <corsika/environment/NuclearComposition.h>
+#include <corsika/process/sibyll/sibyll2.3c.h>
using namespace corsika::units::si;
using namespace corsika::units;
@@ -157,18 +158,34 @@ 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();
- const std::vector<particles::Code> particleList = {
- particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
- particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
-
- Decay model(particleList);
+ Decay model;
model.Init();
/*[[maybe_unused]] const process::EProcessReturn ret =*/model.DoDecay(projectile);
+ // run checks
[[maybe_unused]] const TimeType time = model.GetLifetime(particle);
}
+
+ SECTION("DecayConfiguration") {
+
+ Decay model;
+
+ const std::vector<particles::Code> particleTestList = {
+ particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
+ particles::Code::Lambda0Bar, particles::Code::NuE, particles::Code::D0Bar};
+
+ for (auto& pCode : particleTestList) {
+ model.SetUnstable(pCode);
+ // get state of sibyll internal config
+ REQUIRE(0 <= s_csydec_.idb[abs(process::sibyll::ConvertToSibyllRaw(pCode)) - 1]);
+
+ model.SetStable(pCode);
+ // get state of sibyll internal config
+ REQUIRE(0 >= s_csydec_.idb[abs(process::sibyll::ConvertToSibyllRaw(pCode)) - 1]);
+ }
+ }
}