added decay configuration routines, turned on decay for sibyll interaction

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:

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;
+      }
+      
+    };
   }
 }