added interaction length unit to cascade example

Documentation/Examples/cascade_example.cc

@@ -20,6 +20,8 @@
 #include <corsika/random/RNGManager.h>
 #include <corsika/cascade/sibyll2.3c.h>
 
+//#include <corsika/units/PhysicalConstants.h>
+#include <corsika/units/PhysicalUnits.h>
 using namespace corsika;
 using namespace corsika::process;
 using namespace corsika::units;
@@ -39,21 +41,40 @@ public:
   double MinStepLength(Particle& p) const {
     // beam particles for sibyll : 1, 2, 3 for p, pi, k
     int kBeam   = 1;
+
+    /* 
+       the target should be defined by the Environment,
+       ideally as full particle object so that the four momenta 
+       and the boosts can be defined..
+     */
     // target nuclei: A < 18
-    int kTarget = 0.4*16 + 0.6*14;
+    // FOR NOW: assume target is oxygen
+    int kTarget = 16;
     double beamEnergy =  p.GetEnergy() / 1_GeV; 
     std::cout << "ProcessSplit: " << "MinStep: en: " << beamEnergy << " pid:" << kBeam << std::endl;
     double prodCrossSection,dummy;
+    
     sib_sigma_hnuc_(kBeam, kTarget, beamEnergy, prodCrossSection, dummy );
+    
     std::cout << "ProcessSplit: " << "MinStep: sibyll return: " << prodCrossSection << std::endl;
     CrossSectionType sig = prodCrossSection  / 1000. * barn;
-    std::cout << "ProcessSplit: " << "MinStep: CrossSection= " << sig << std::endl;
+    std::cout << "ProcessSplit: " << "MinStep: CrossSection (mb): " << sig / 1_mbarn << std::endl;
 
+    const MassType nucleon_mass = 0.93827_GeV / corsika::units::si::constants::cSquared;
+    std::cout << "ProcessSplit: " << "nucleon mass " << nucleon_mass << std::endl;
     // calculate interaction length in medium
-
+    double int_length = kTarget * ( nucleon_mass / 1_g ) / ( sig / 1_cmeter / 1_cmeter );
     // pick random step lenth
-    
-    return prodCrossSection;
+    std::cout << "ProcessSplit: " << "interaction length (g/cm2): " << int_length << std::endl;
+    // add exponential sampling
+    int a = 0;
+    const double next_step = -int_length * log(s_rndm_(a));
+    /*
+      what are the units of the output? slant depth or 3space length?
+      
+    */
+    std::cout << "ProcessSplit: " << "next step (g/cm2): " << next_step << std::endl;
+    return next_step;
   }
 
   template <typename Particle, typename Trajectory, typename Stack>
@@ -65,12 +86,21 @@ public:
   template <typename Particle, typename Stack>
   void DoDiscrete(Particle& p, Stack& s) const {
     // get energy of particle from stack
-    // stack is in GeV in lab. frame
-    // convert to GeV in cm. frame (assuming proton at rest as target)
+    /*
+      stack is in GeV in lab. frame
+      convert to GeV in cm. frame 
+      (assuming proton at rest as target AND 
+      assuming no pT, i.e. shower frame-z is aligned with hadron-int-frame-z)
+    */
     EnergyType E   = p.GetEnergy();
     EnergyType Ecm = sqrt( 2. * E * 0.93827_GeV );
     // FOR NOW: set beam to proton
     int kBeam   = 13; //p.GetPID();
+    /* 
+       the target should be defined by the Environment,
+       ideally as full particle object so that the four momenta 
+       and the boosts can be defined..
+     */
     // FOR NOW: set target to proton
     int kTarget = 1; //p.GetPID();
     std::cout << "ProcessSplit: " << " DoDiscrete: E(GeV):" << E / 1_GeV << " Ecm(GeV): " << Ecm / 1_GeV << std::endl;
@@ -83,6 +113,8 @@ public:
       double sqs = Ecm / 1_GeV;
       // running sibyll
       sibyll_( kBeam, kTarget, sqs);
+      // running decays
+      //decsib_();
       // print final state
       int print_unit = 6;
       sib_list_( print_unit );
@@ -96,7 +128,7 @@ public:
       
       // add particles from sibyll to stack
       for(int i=0; i<s_plist_.np; ++i){
-	//transform to lab. frame
+	//transform to lab. frame, primitve
 	const double en_lab = gambet * s_plist_.p[2][i] + gamma * s_plist_.p[3][i];
 	// add to corsika stack
 	s.NewParticle().SetEnergy( en_lab * 1_GeV );
@@ -108,6 +140,8 @@ public:
   {
     fCount = 0;
 
+    // define reference frame? --> defines boosts between corsika stack and model stack.
+    
     // initialize random numbers for sibyll
     // FOR NOW USE SIBYLL INTERNAL !!!
     rnd_ini_();
@@ -127,6 +161,8 @@ public:
     
     //initialize Sibyll
     sibyll_ini_();
+
+    // set particles stable / unstable
   }
   
   int GetCount() { return fCount; }