Sibyll

closes #72 (closed)

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);