fixes in particleConversion for sibyll, added decay table to interface, added pid to sibstack

Documentation/Examples/cascade_example.cc

@@ -42,19 +42,11 @@ public:
 
   template <typename Particle>
   double MinStepLength(Particle& p) const {
+
     // beam particles for sibyll : 1, 2, 3 for p, pi, k
-    corsika::particles::Code c_id = p.GetPID();
-    corsika::process::sibyll::Code s_id = process::sibyll::ConvertToSibyll( p.GetPID() );
-    corsika::particles::Code p_id_2 = process::sibyll::ConvertFromSibyll( corsika::process::sibyll::Code::Proton );
-    cout << p_id_2 << endl;
-    std::cout << "MinStepLength: particle input " << "corsika id: " << c_id << std::endl;
-    auto test = static_cast<corsika::process::sibyll::SibyllCodeIntType>(s_id);
-    std::cout << "MinStepLength: particle input " << "sibyll id: |" << (int)test << "|" <<std::endl;
-											  // std::cout << "MinStepLength: particle input " << "sibyll id: " << process::sibyll::ConvertToSibyll( p.GetPID() ) << std::endl;
-    cout << p.GetPID() << " --> " << process::sibyll::ConvertToSibyllRaw( p.GetPID() ) << endl;
+    // read from cross section code table
+    int kBeam   =  1;
     
-    int kBeam   = 1;
-
     /* 
        the target should be defined by the Environment,
        ideally as full particle object so that the four momenta 
@@ -98,6 +90,9 @@ public:
 
   template <typename Particle, typename Stack>
   void DoDiscrete(Particle& p, Stack& s) const {
+    cout << "DoDiscrete: " << p.GetPID() << " interaction? " << process::sibyll::CanInteract( p.GetPID() )  << endl; 
+    if( process::sibyll::CanInteract( p.GetPID() ) ){
+    
     // get energy of particle from stack
     /*
       stack is in GeV in lab. frame
@@ -107,8 +102,9 @@ public:
     */
     EnergyType E   = p.GetEnergy();
     EnergyType Ecm = sqrt( 2. * E * 0.93827_GeV );
-    // FOR NOW: set beam to proton
-    int kBeam   = 13; //p.GetPID();
+
+    int kBeam   = process::sibyll::ConvertToSibyllRaw( p.GetPID() );
+    
     /* 
        the target should be defined by the Environment,
        ideally as full particle object so that the four momenta 
@@ -156,10 +152,15 @@ public:
 	//transform to lab. frame, primitve
 	const double en_lab = gambet * s_plist_.p[2][i] + gamma * p.GetEnergy();	
 	// add to corsika stack
-	s.NewParticle().SetEnergy( en_lab * 1_GeV );
+	auto pnew = s.NewParticle();
+	pnew.SetEnergy( en_lab * 1_GeV );
+	pnew.SetPID( process::sibyll::ConvertFromSibyll( p.GetPID() ) );
+	
       }
      
     }
+    }else
+      p.Delete();
   }
   
   
@@ -183,14 +184,31 @@ public:
     // test random number generator
     std::cout << "ProcessSplit: " << " test sequence of random numbers."  << std::endl;
     int a = 0;
-    for(int i=0; i<5; ++i)
+    for(int i=0; i<8; ++i)
       std::cout << i << " " << s_rndm_(a) << std::endl;
     
     //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::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();
+    }
   }
   
   int GetCount() { return fCount; }

Framework/Cascade/SibStack.h

@@ -6,6 +6,7 @@
 
 #include <corsika/stack/Stack.h>
 #include <corsika/cascade/sibyll2.3c.h>
+#include <corsika/process/sibyll/ParticleConversion.h>
 
 using namespace std;
 using namespace corsika::stack;
@@ -17,8 +18,8 @@ class SibStackData {
   
   void Clear() { s_plist_.np = 0; }
   
-  int GetSize() const { return s_plist_.np-1;  }
-  int GetCapacity() const { return s_plist_.np-1; }
+  int GetSize() const { return s_plist_.np;  }
+  int GetCapacity() const { return 8000; }
 
   
   void SetId(const int i, const int v) { s_plist_.llist[i]=v; }
@@ -45,6 +46,9 @@ class ParticleInterface : public ParticleBase<StackIteratorInterface> {
  public:
   void SetEnergy(const double v) { GetStackData().SetEnergy(GetIndex(), v); }
   double 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())); }
+  
 };
 
 

Framework/Cascade/sibyll2.3c.h

@@ -23,6 +23,14 @@ extern"C"{
   }s_plist_;
 
 
+  extern struct {
+    double  cbr[223+16+12+8];
+    int    kdec[1338+6*(16+12+8)];
+    int    lbarp[99];
+    int     idb[99];
+  }s_csydec_;
+
+  
   // additional particle stack for the mother particles of unstable particles
   // stable particles have entry zero
   extern struct { 

Processes/Sibyll/code_generator.py

@@ -28,7 +28,7 @@ def read_sibyll_codes(filename, pythia_db):
             line = line.strip()
             if line[0] == '#':
                 continue            
-            identifier, sib_code, canInteractFlag = line.split()
+            identifier, sib_code, canInteractFlag, xsctnId = line.split()
             try:
                 pythia_db[identifier]["sibyll_code"] = int(sib_code)
                 pythia_db[identifier]["sibyll_canInteract"] = int(canInteractFlag)
@@ -148,3 +148,4 @@ if __name__ == "__main__":
         print(generate_corsika2sibyll(pythia_db), file=f)
         print(generate_known_particle(pythia_db), file=f)
         print(generate_sibyll2corsika(pythia_db), file=f)
+        print(generate_interacting_particle(pythia_db), file=f)

Processes/Sibyll/sibyll_codes.dat

@@ -12,7 +12,8 @@ NuTau 92 0 0
 NuTauBar 93 0 0
 Gamma 1 0 0
 Pi0 6 1 2
-Rho0 27 1 2
+# rho0 could interact but sibyll has no cross section/interaction length. was used for gamma had int
+Rho0 27 0 0
 K0Long 11 1 3
 PiPlus 7 1 2
 PiMinus 8 1 2

Processes/Sibyll/testSibyll.cc

@@ -23,20 +23,20 @@ TEST_CASE("Sibyll", "[processes]") {
 
   SECTION("Sibyll -> Corsika") {
     REQUIRE(corsika::particles::Electron::GetCode() ==
-            process::sibyll::ConvertFromSibyll(process::sibyll::Code::Electron));
+            process::sibyll::ConvertFromSibyll(process::sibyll::SibyllCode::Electron));
   }
 
   SECTION("Corsika -> Sibyll") {
     REQUIRE(process::sibyll::ConvertToSibyll(corsika::particles::Electron::GetCode()) ==
-            process::sibyll::Code::Electron);
+            process::sibyll::SibyllCode::Electron);
     REQUIRE(process::sibyll::ConvertToSibyllRaw(corsika::particles::Proton::GetCode()) ==
             13 );
   }
 
-  SECTION("handledBySibyll") {
-    REQUIRE(process::sibyll::HandledBySibyll(corsika::particles::Electron::GetCode()));
+  SECTION("KnownBySibyll") {
+    REQUIRE(process::sibyll::KnownBySibyll(corsika::particles::Electron::GetCode()));
 
     REQUIRE_FALSE(
-        process::sibyll::HandledBySibyll(corsika::particles::XiPrimeC0::GetCode()));
+        process::sibyll::KnownBySibyll(corsika::particles::XiPrimeC0::GetCode()));
   }
 }