diff --git a/Documentation/Examples/boundary_example.cc b/Documentation/Examples/boundary_example.cc
index 169cb1c6a30a46a16494ac7e37b5bae2e875c256..c046ba18af3dbac225c46d520ad289b7fddc32f8 100644
--- a/Documentation/Examples/boundary_example.cc
+++ b/Documentation/Examples/boundary_example.cc
@@ -123,9 +123,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);
diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index 7690d664811f4192d224044ab39ebea7b92d6136..c336ee8a54cdda75e35d905283ff437e3b73f44c 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -100,12 +100,12 @@ int main() {
const std::vector<particles::Code> trackedHadrons = {
particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
- particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
+ particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short, particles::Code::MuMinus,particles::Code::MuPlus};
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");
diff --git a/Documentation/Examples/vertical_EAS.cc b/Documentation/Examples/vertical_EAS.cc
index 12b1cac0e59c082da7f374e49c777bb3b96ecee2..69771bee47663940c3d255bcc74f3bd3c013a11f 100644
--- a/Documentation/Examples/vertical_EAS.cc
+++ b/Documentation/Examples/vertical_EAS.cc
@@ -169,13 +169,9 @@ int main() {
// setup processes, decays and interactions
tracking_line::TrackingLine tracking;
- const std::vector<particles::Code> trackedHadrons = {
- particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::KPlus,
- particles::Code::KMinus, particles::Code::K0Long, particles::Code::K0Short};
-
process::sibyll::Interaction sibyll;
process::sibyll::NuclearInteraction sibyllNuc(sibyll, env);
- process::sibyll::Decay decay(trackedHadrons);
+ process::sibyll::Decay decay;
//~ process::pythia::Decay decay(trackedHadrons);
process::particle_cut::ParticleCut cut(2_GeV);
diff --git a/Processes/Sibyll/Decay.cc b/Processes/Sibyll/Decay.cc
index 9cb42f0c5c2f20db356218450684b41d5d6f1bd5..9e93748328754fd0e028fc9e70aab95a1d9b8911 100644
--- a/Processes/Sibyll/Decay.cc
+++ b/Processes/Sibyll/Decay.cc
@@ -29,77 +29,62 @@ 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::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::SetStable(const vector<particles::Code> vParticleList) {
+ for (auto p : vParticleList) Decay::SetStable(p);
}
- void Decay::SetUnstable(const particles::Code pCode) {
- int s_id = process::sibyll::ConvertToSibyllRaw(pCode);
+ 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::SetDecay(const particles::Code vCode, const bool vMakeUnstable) {
+ vMakeUnstable ? SetUnstable(vCode) : SetStable(vCode);
+ }
+
+ 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 <>
@@ -149,12 +134,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 03efcb97138f2525b676ea66e0bbe121fc4045dc..dd2c48ed1d1e149a9d4fbe320e6ed56fe59a9b85 100644
--- a/Processes/Sibyll/Decay.h
+++ b/Processes/Sibyll/Decay.h
@@ -22,23 +22,33 @@ 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&);
};
diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc
index c089ec083c9c7d9ec1d91136cb0c7796dc958c57..6fc86272b1e9acee7a9e08bf06ade2a940bb6e5b 100644
--- a/Processes/Sibyll/Interaction.cc
+++ b/Processes/Sibyll/Interaction.cc
@@ -48,46 +48,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,
@@ -326,8 +318,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 643bd0d51bd4f387aafa3870f2dddfb0ed1fe9d5..b9d06b9d4479cfa2be00aad315daef37fa058517 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -32,9 +32,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 {
@@ -66,7 +70,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/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 0c42c90abb86a758c45f5ad4a79ebe87d3f1d3fa..7c1b88031625c95096ef72f5d7103f97dfa25441 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -72,6 +72,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;
@@ -164,14 +165,29 @@ TEST_CASE("SibyllInterface", "[processes]") {
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);
[[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]);
+ }
+ }
}