diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index fdc7069499e80599c92ef25e80b63b62d175be60..5c631f5a02b0ba7fd8ac6810e990adf867653f36 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -43,6 +43,35 @@ class ProcessSplit : public corsika::process::BaseProcess<ProcessSplit> {
public:
ProcessSplit() {}
+ void setTrackedParticlesStable() const {
+ /*
+ Sibyll is hadronic generator
+ only hadrons decay
+ */
+ // set particles unstable
+ corsika::process::sibyll::setHadronsUnstable();
+ // make tracked particles stable
+ std::cout << "ProcessSplit: setting tracked hadrons stable.." << std::endl;
+ setup::Stack ds;
+ ds.NewParticle().SetPID(Code::PiPlus);
+ ds.NewParticle().SetPID(Code::PiMinus);
+ ds.NewParticle().SetPID(Code::KPlus);
+ ds.NewParticle().SetPID(Code::KMinus);
+ ds.NewParticle().SetPID(Code::K0Long);
+ ds.NewParticle().SetPID(Code::K0Short);
+
+ for (auto& p : ds) {
+ int s_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
+ // set particle stable by setting table value negative
+ // cout << "ProcessSplit: doDiscrete: setting " << p.GetPID() << "(" << s_id << ")"
+ // << " stable in Sibyll .." << endl;
+ s_csydec_.idb[s_id - 1] = (-1) * abs( s_csydec_.idb[s_id - 1] );
+ p.Delete();
+ }
+
+ }
+
+
template <typename Particle>
double MinStepLength(Particle& p, setup::Trajectory&) const {
@@ -211,7 +240,8 @@ public:
// running sibyll, filling stack
sibyll_( kBeam, kTarget, sqs);
// running decays
- //decsib_();
+ setTrackedParticlesStable();
+ decsib_();
// print final state
int print_unit = 6;
sib_list_( print_unit );
@@ -228,6 +258,9 @@ public:
super_stupid::MomentumVector Ptot_final(rootCS, {0.0_newton_second, 0.0_newton_second, 0.0_newton_second});
for (auto &psib: ss){
++i;
+ // skip particles that have decayed in Sibyll
+ if( abs(s_plist_.llist[ i ]) > 100 ) continue;
+
//transform energy to lab. frame, primitve
// compute beta_vec * p_vec
// arbitrary Lorentz transformation based on sibyll routines
@@ -300,29 +333,10 @@ public:
// initialize Sibyll
sibyll_ini_();
- // set particles stable / unstable
- // use stack to loop over particles
- setup::Stack ds;
- ds.NewParticle().SetPID(Code::Proton);
- ds.NewParticle().SetPID(Code::Neutron);
- ds.NewParticle().SetPID(Code::PiPlus);
- ds.NewParticle().SetPID(Code::PiMinus);
- ds.NewParticle().SetPID(Code::Pi0);
- ds.NewParticle().SetPID(Code::KPlus);
- ds.NewParticle().SetPID(Code::KMinus);
- ds.NewParticle().SetPID(Code::K0Long);
- ds.NewParticle().SetPID(Code::K0Short);
-
- for (auto& p : ds) {
- int s_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
- // set particle stable by setting table value negative
- cout << "ProcessSplit: Init: setting " << p.GetPID() << "(" << s_id << ")"
- << " stable in Sibyll .." << endl;
- s_csydec_.idb[s_id] = -s_csydec_.idb[s_id - 1];
- p.Delete();
- }
+ setTrackedParticlesStable();
}
+
int GetCount() { return fCount; }
private:
diff --git a/Processes/Sibyll/ProcessDecay.h b/Processes/Sibyll/ProcessDecay.h
index 4989ab90a19f3b7f07080dc2647dad1bfa5cf21f..ecae303e79169ce0961793c187f90253f135f219 100644
--- a/Processes/Sibyll/ProcessDecay.h
+++ b/Processes/Sibyll/ProcessDecay.h
@@ -12,62 +12,107 @@ namespace corsika::process {
namespace sibyll {
-class ProcessDecay : public corsika::process::BaseProcess<ProcessDecay> {
-public:
- ProcessDecay() {}
- void Init() {}
-
- template <typename Particle>
- double MinStepLength(Particle& p, setup::Trajectory&) const {
- EnergyType E = p.GetEnergy();
- MassType m = corsika::particles::GetMass(p.GetPID());
- // env.GetDensity();
- const MassDensityType density = 1.25e-3 * kilogram / ( 1_cm * 1_cm * 1_cm );
+
+ void setHadronsUnstable(){
+ // name? also makes EM particles stable
+
+ // loop over all particles in sibyll
+ // should be changed to loop over human readable list
+ // i.e. corsika::particles::ListOfParticles()
+ std::cout << "Sibyll: setting hadrons unstable.." << std::endl;
+ // make ALL particles unstable, then set EM stable
+ for(auto &p: corsika2sibyll){
+ //std::cout << (int)p << std::endl;
+ const int sibCode = (int)p;
+ // skip unknown and antiparticles
+ if( sibCode< 1) continue;
+ //std::cout << "Sibyll: Decay: setting " << ConvertFromSibyll( static_cast<SibyllCode> ( sibCode ) ) << " unstable" << std::endl;
+ s_csydec_.idb[ sibCode - 1 ] = abs( s_csydec_.idb[ sibCode - 1 ] );
+ //std::cout << "decay table value: " << s_csydec_.idb[ sibCode - 1 ] << std::endl;
+ }
+ // set Leptons and Proton and Neutron stable
+ // use stack to loop over particles
+ setup::Stack ds;
+ ds.NewParticle().SetPID(corsika::particles::Code::Proton);
+ ds.NewParticle().SetPID(corsika::particles::Code::Neutron);
+ ds.NewParticle().SetPID(corsika::particles::Code::Electron);
+ ds.NewParticle().SetPID(corsika::particles::Code::Positron);
+ ds.NewParticle().SetPID(corsika::particles::Code::NuE);
+ ds.NewParticle().SetPID(corsika::particles::Code::NuEBar);
+ ds.NewParticle().SetPID(corsika::particles::Code::MuMinus);
+ ds.NewParticle().SetPID(corsika::particles::Code::MuPlus);
+ ds.NewParticle().SetPID(corsika::particles::Code::NuMu);
+ ds.NewParticle().SetPID(corsika::particles::Code::NuMuBar);
+
+ for (auto& p : ds) {
+ int s_id = process::sibyll::ConvertToSibyllRaw(p.GetPID());
+ // set particle stable by setting table value negative
+ // cout << "Sibyll: setting " << p.GetPID() << "(" << s_id << ")"
+ // << " stable in Sibyll .." << endl;
+ s_csydec_.idb[ s_id - 1 ] = (-1) * abs(s_csydec_.idb[ s_id - 1 ]);
+ p.Delete();
+ }
+
+ }
+
- const double gamma = E / m / constants::cSquared;
- // lifetimes not implemented yet
- TimeType t0;
- switch( p.GetPID() ){
- case corsika::particles::Code::PiPlus :
- t0 = 2.6e-8 * 1_s;
- break;
-
- case corsika::particles::Code::PiMinus :
- t0 = 2.6e-8 * 1_s;
- break;
+ class ProcessDecay : public corsika::process::BaseProcess<ProcessDecay> {
+ public:
+ ProcessDecay() {}
+ void Init() {
+ //setHadronsUnstable();
+ }
- case corsika::particles::Code::KPlus :
- t0 = 1.e-5 * 1_s;
- break;
+ void setAllStable(){
+ // name? also makes EM particles stable
+
+ // loop over all particles in sibyll
+ // should be changed to loop over human readable list
+ // i.e. corsika::particles::ListOfParticles()
+ for(auto &p: corsika2sibyll){
+ //std::cout << (int)p << std::endl;
+ const int sibCode = (int)p;
+ // skip unknown and antiparticles
+ if( sibCode< 1) continue;
+ std::cout << "Sibyll: Decay: setting " << ConvertFromSibyll( static_cast<SibyllCode> ( sibCode ) ) << " stable" << std::endl;
+ s_csydec_.idb[ sibCode - 1 ] = -1 * abs( s_csydec_.idb[ sibCode - 1 ] );
+ std::cout << "decay table value: " << s_csydec_.idb[ sibCode - 1 ] << std::endl;
+ }
+ }
- case corsika::particles::Code::KMinus :
- t0 = 1.e-5 * 1_s;
- break;
+ friend void setHadronsUnstable();
- default:
- t0 = 1.e8 * 1_s;
- break;
- }
- cout << "ProcessDecay: MinStep: t0: " << t0 << endl;
- cout << "ProcessDecay: MinStep: gamma: " << gamma << endl;
- cout << "ProcessDecay: MinStep: density: " << density << endl;
- // return as column density
- const double x0 = density * t0 * gamma * constants::c / kilogram * 1_cm * 1_cm;
- cout << "ProcessDecay: MinStep: x0: " << x0 << endl;
- return x0;
- }
-
- template <typename Particle, typename Stack>
- void DoDiscrete(Particle& p, Stack& s) const {
+ template <typename Particle>
+ double MinStepLength(Particle& p, setup::Trajectory&) const {
+ EnergyType E = p.GetEnergy();
+ MassType m = corsika::particles::GetMass(p.GetPID());
+ // env.GetDensity();
+
+ const MassDensityType density = 1.25e-3 * kilogram / ( 1_cm * 1_cm * 1_cm );
- }
-
- template <typename Particle, typename Stack>
- EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
- return EProcessReturn::eOk;
- }
+ const double gamma = E / m / constants::cSquared;
-};
+ TimeType t0 = GetLifetime( p.GetPID() );
+ cout << "ProcessDecay: MinStep: t0: " << t0 << endl;
+ cout << "ProcessDecay: MinStep: gamma: " << gamma << endl;
+ cout << "ProcessDecay: MinStep: density: " << density << endl;
+ // return as column density
+ const double x0 = density * t0 * gamma * constants::c / kilogram * 1_cm * 1_cm;
+ cout << "ProcessDecay: MinStep: x0: " << x0 << endl;
+ return x0;
+ }
+
+ template <typename Particle, typename Stack>
+ void DoDiscrete(Particle& p, Stack& s) const {
+
+ }
+
+ template <typename Particle, typename Stack>
+ EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
+ return EProcessReturn::eOk;
+ }
+
+ };
}
}