diff --git a/Documentation/Examples/cascade_example.cc b/Documentation/Examples/cascade_example.cc
index d274cac0de63eab87c8a0b5c9c6700ddfe495b08..367b7700c91cf74d29121e06013f2bc795c0c322 100644
--- a/Documentation/Examples/cascade_example.cc
+++ b/Documentation/Examples/cascade_example.cc
@@ -23,6 +23,8 @@
#include <corsika/cascade/sibyll2.3c.h>
#include <corsika/process/sibyll/ParticleConversion.h>
+#include <corsika/process/sibyll/ProcessDecay.h>
+
#include <corsika/units/PhysicalUnits.h>
using namespace corsika;
@@ -42,6 +44,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 {
@@ -210,7 +241,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 );
@@ -227,6 +259,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
@@ -299,29 +334,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:
@@ -334,6 +350,7 @@ double s_rndm_(int&) {
return rmng() / (double)rmng.max();
}
+
int main() {
// coordinate system, get global frame of reference
@@ -346,7 +363,8 @@ int main() {
stack_inspector::StackInspector<setup::Stack> p0(true);
ProcessSplit p1;
- const auto sequence = p0 + p1;
+ corsika::process::sibyll::ProcessDecay p2;
+ const auto sequence = p0 + p1 + p2;
setup::Stack stack;
corsika::cascade::Cascade EAS(tracking, sequence, stack);
diff --git a/Framework/Cascade/SibStack.h b/Framework/Cascade/SibStack.h
index 57a3f9c16b099c6d00e71b405223f6c3f7eb41c5..2275cdf491ae34433f8afba51041d867d301b8e2 100644
--- a/Framework/Cascade/SibStack.h
+++ b/Framework/Cascade/SibStack.h
@@ -6,6 +6,7 @@
#include <corsika/cascade/sibyll2.3c.h>
#include <corsika/process/sibyll/ParticleConversion.h>
+#include <corsika/units/PhysicalUnits.h>
#include <corsika/stack/Stack.h>
using namespace std;
@@ -64,11 +65,12 @@ class ParticleInterface : public ParticleBase<StackIteratorInterface> {
using ParticleBase<StackIteratorInterface>::GetIndex;
public:
- void SetEnergy(const double v) { GetStackData().SetEnergy(GetIndex(), v); }
+ void SetEnergy(const EnergyType v) { GetStackData().SetEnergy(GetIndex(), v ); }
EnergyType GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); }
void SetPID(const int v) { GetStackData().SetId(GetIndex(), v); }
corsika::process::sibyll::SibyllCode GetPID() const { return static_cast<corsika::process::sibyll::SibyllCode> (GetStackData().GetId(GetIndex())); }
super_stupid::MomentumVector GetMomentum() const { return GetStackData().GetMomentum(GetIndex()); }
+ void SetMomentum(const super_stupid::MomentumVector& v) { GetStackData().SetMomentum(GetIndex(), v); }
};
diff --git a/Framework/Units/PhysicalUnits.h b/Framework/Units/PhysicalUnits.h
index e97003eee0a0e649021e20c420776afd6e586e59..dfddec6277470476ee5fada52c546352e9061fde 100644
--- a/Framework/Units/PhysicalUnits.h
+++ b/Framework/Units/PhysicalUnits.h
@@ -36,9 +36,10 @@ namespace corsika::units::si {
phys::units::quantity<phys::units::electric_charge_d, double>;
using EnergyType = phys::units::quantity<phys::units::energy_d, double>;
using MassType = phys::units::quantity<phys::units::mass_d, double>;
- using MomentumType = phys::units::quantity<momentum_d, double>;
+ using MassDensityType = phys::units::quantity<phys::units::mass_density_d, double>;
using CrossSectionType = phys::units::quantity<sigma_d, double>;
-
+ using MomentumType = phys::units::quantity<momentum_d, double>;
+
} // end namespace corsika::units::si
/**
diff --git a/Processes/Sibyll/CMakeLists.txt b/Processes/Sibyll/CMakeLists.txt
index f2e4a9fee1aabbbb5196a597db36de06e6385594..a0985a1ad2304e7e5bd4f1024cf43b5cefaf8bc4 100644
--- a/Processes/Sibyll/CMakeLists.txt
+++ b/Processes/Sibyll/CMakeLists.txt
@@ -20,6 +20,7 @@ set (
set (
MODEL_HEADERS
ParticleConversion.h
+ ProcessDecay.h
${PROJECT_BINARY_DIR}/Processes/Sibyll/Generated.inc
)
diff --git a/Processes/Sibyll/ProcessDecay.h b/Processes/Sibyll/ProcessDecay.h
new file mode 100644
index 0000000000000000000000000000000000000000..449a281b43ee2ff83c889441f5b567e53989adac
--- /dev/null
+++ b/Processes/Sibyll/ProcessDecay.h
@@ -0,0 +1,151 @@
+#ifndef _include_ProcessDecay_h_
+#define _include_ProcessDecay_h_
+
+#include <corsika/process/ContinuousProcess.h>
+
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/process/sibyll/ParticleConversion.h>
+#include <corsika/cascade/SibStack.h>
+
+//using namespace corsika::particles;
+
+namespace corsika::process {
+
+ namespace sibyll {
+
+
+ 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();
+ }
+
+ }
+
+
+ class ProcessDecay : public corsika::process::BaseProcess<ProcessDecay> {
+ public:
+ ProcessDecay() {}
+ void Init() {
+ //setHadronsUnstable();
+ }
+
+ 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;
+ }
+ }
+
+ friend void setHadronsUnstable();
+
+ 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 );
+
+ 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 {
+ SibStack ss;
+ ss.Clear();
+ // copy particle to sibyll stack
+ auto pin = ss.NewParticle();
+ pin.SetPID( process::sibyll::ConvertToSibyllRaw( p.GetPID() ) );
+ pin.SetEnergy( p.GetEnergy() );
+ pin.SetMomentum( p.GetMomentum() );
+ // set all particles/hadrons unstable
+ setHadronsUnstable();
+ // call sibyll decay
+ std::cout << "calling Sibyll decay routine.." << std::endl;
+ decsib_();
+ // print output
+ //int print_unit = 6;
+ //sib_list_( print_unit );
+ // copy particles from sibyll stack to corsika
+ int i = -1;
+ for (auto &psib: ss){
+ ++i;
+ // FOR NOW: skip particles that have decayed in Sibyll, move to iterator?
+ if( abs(s_plist_.llist[ i ]) > 100 ) continue;
+ // add to corsika stack
+ //cout << "decay product: " << process::sibyll::ConvertFromSibyll( psib.GetPID() ) << endl;
+ auto pnew = s.NewParticle();
+ pnew.SetEnergy( psib.GetEnergy() );
+ pnew.SetPID( process::sibyll::ConvertFromSibyll( psib.GetPID() ) );
+ pnew.SetMomentum( psib.GetMomentum() );
+ }
+ // empty sibyll stack
+ ss.Clear();
+ // remove original particle from stack
+ p.Delete();
+ }
+
+ template <typename Particle, typename Stack>
+ EProcessReturn DoContinuous(Particle&, setup::Trajectory&, Stack&) const {
+ return EProcessReturn::eOk;
+ }
+
+ };
+ }
+}
+
+#endif