added momentum to sibyll stack, implemented boost between sibyll stack and corsika stack

Documentation/Examples/cascade_example.cc

@@ -23,6 +23,7 @@
 #include <corsika/process/sibyll/ParticleConversion.h>
 
 #include <corsika/units/PhysicalUnits.h>
+
 using namespace corsika;
 using namespace corsika::process;
 using namespace corsika::units;
@@ -42,9 +43,13 @@ public:
   template <typename Particle>
   double MinStepLength(Particle& p) const {
 
+    const Code corsikaBeamId = p.GetPID();
+    
     // beam particles for sibyll : 1, 2, 3 for p, pi, k
     // read from cross section code table
-    int kBeam   =  1;
+    int kBeam   =  process::sibyll::GetSibyllXSCode( corsikaBeamId );   
+
+    bool kInteraction = process::sibyll::CanInteract( corsikaBeamId );
     
     /* 
        the target should be defined by the Environment,
@@ -53,30 +58,44 @@ public:
      */
     // target nuclei: A < 18
     // FOR NOW: assume target is oxygen
-    int kTarget = 1;
-    double beamEnergy =  p.GetEnergy() / 1_GeV; 
-    std::cout << "ProcessSplit: " << "MinStep: en: " << beamEnergy << " pid:" << kBeam << std::endl;
-    double prodCrossSection,dummy,dum1,dum2,dum3,dum4;
-    double dumdif[3];
-
-    if(kTarget==1)
-      sib_sigma_hp_(kBeam, beamEnergy, dum1, dum2, prodCrossSection, dumdif,dum3, dum4 );
-    else
-      sib_sigma_hnuc_(kBeam, kTarget, beamEnergy, prodCrossSection, dummy );
+    int kTarget = 16;
+    double beamEnergy =  p.GetEnergy() / 1_GeV;
+#warning boost to cm. still missing, sibyll cross section input is cm. energy!
+    
+    std::cout << "ProcessSplit: " << "MinStep: input en: " << beamEnergy
+      	      << " beam can interact:" << kBeam
+	      << " beam XS code:" << kBeam
+      	      << " beam pid:" << p.GetPID()
+	      << " target mass number:" << kTarget << std::endl;
+
+    double next_step;
+    if(kInteraction){
+      
+      double prodCrossSection,dummy,dum1,dum2,dum3,dum4;
+      double dumdif[3];
+      
+      if(kTarget==1)
+	sib_sigma_hp_(kBeam, beamEnergy, dum1, dum2, prodCrossSection, dumdif,dum3, dum4 );
+      else
+	sib_sigma_hnuc_(kBeam, kTarget, beamEnergy, prodCrossSection, dummy );
+    
+      std::cout << "ProcessSplit: " << "MinStep: sibyll return: " << prodCrossSection << std::endl;
+      CrossSectionType sig = prodCrossSection * 1_mbarn;
+      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
+      std::cout << "ProcessSplit: " << "interaction length (g/cm2): " << int_length << std::endl;
+      // add exponential sampling
+      int a = 0;
+      next_step = -int_length * log(s_rndm_(a));
+    }else
+#warning define infinite interaction length? then we can skip the test in DoDiscrete()
+      next_step = 1.e8;
     
-    std::cout << "ProcessSplit: " << "MinStep: sibyll return: " << prodCrossSection << std::endl;
-    CrossSectionType sig = prodCrossSection * 1_mbarn;
-    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
-    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?
       
@@ -95,28 +114,75 @@ public:
   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
-      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 );
+      cout << "defining coordinates" << endl;
+      // coordinate system, get global frame of reference
+      CoordinateSystem rootCS = CoordinateSystem::CreateRootCS();
+
+      QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m};
+      Point pOrig(rootCS, coordinates);
+      
+      /* 
+	 the target should be defined by the Environment,
+	 ideally as full particle object so that the four momenta 
+	 and the boosts can be defined..
+
+	 here we need: GetTargetMassNumber() or GetTargetPID()??
+	               GetTargetMomentum() (zero in EAS)
+      */
+      // FOR NOW: set target to proton
+      int kTarget = 1; //p.GetPID();
 
-    int kBeam   = process::sibyll::ConvertToSibyllRaw( p.GetPID() );
+      // proton mass in units of energy
+      const EnergyType proton_mass_en = 0.93827_GeV ; //0.93827_GeV / si::constants::cSquared ;
+
+      cout << "defining target momentum.." << endl;
+      // FOR NOW: target is always at rest
+      const EnergyType Etarget = 0. * 1_GeV + proton_mass_en;      
+      const auto pTarget = super_stupid::MomentumVector(rootCS, 0. * 1_GeV / si::constants::c,  0. * 1_GeV / si::constants::c, 0. * 1_GeV / si::constants::c);
+      cout << "target momentum (GeV/c): " << pTarget.GetComponents() / 1_GeV * si::constants::c << endl;
+      //      const auto pBeam = super_stupid::MomentumVector(rootCS, 0. * 1_GeV / si::constants::c,  0. * 1_GeV / si::constants::c, 0. * 1_GeV / si::constants::c);
+      //      cout << "beam momentum: " << pBeam.GetComponents() << endl;
+      cout << "beam momentum (GeV/c): " << p.GetMomentum().GetComponents() / 1_GeV * si::constants::c << endl;
+      
+      // 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 AND 
+	assuming no pT, i.e. shower frame-z is aligned with hadron-int-frame-z)
+      */
+      // cout << "defining total energy" << endl;
+      // total energy: E_beam + E_target
+      // in lab. frame: E_beam + m_target*c**2
+      EnergyType E   = p.GetEnergy();
+      EnergyType Etot   = E + Etarget;
+      // cout << "tot. energy: " << Etot / 1_GeV << endl;
+      // cout << "defining total momentum" << endl;
+      // total momentum
+      super_stupid::MomentumVector Ptot = p.GetMomentum(); // + pTarget;
+      // cout << "tot. momentum: " << Ptot.GetComponents() / 1_GeV * si::constants::c << endl;
+      // cout << "inv. mass.." << endl;
+      // invariant mass, i.e. cm. energy
+      EnergyType Ecm = sqrt( Etot * Etot - Ptot.squaredNorm() * si::constants::cSquared ); //sqrt( 2. * E * 0.93827_GeV );
+      // cout << "inv. mass: " << Ecm / 1_GeV << endl; 
+      // cout << "boost parameters.." << endl;
+       /*
+	get transformation between Stack-frame and SibStack-frame
+	for EAS Stack-frame is lab. frame, could be different for CRMC-mode
+	the transformation should be derived from the input momenta
+      */      
+      //      const double gamma  = ( E + proton_mass * si::constants::cSquared ) / Ecm ;
+      // const double gambet =  sqrt( E * E - proton_mass * proton_mass ) / Ecm;
+      const double gamma  = Etot / Ecm ;
+      const auto gambet = Ptot / (Ecm / si::constants::c );      
+
+      std::cout << "ProcessSplit: " << " DoDiscrete: gamma:" << gamma << endl;
+      std::cout << "ProcessSplit: " << " DoDiscrete: gambet:" << gambet.GetComponents() << endl;
+      
+      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 
-       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;
+      std::cout << "ProcessSplit: " << " DoDiscrete: E(GeV):" << E / 1_GeV << " Ecm(GeV): " << Ecm / 1_GeV << std::endl;
     if (E < 8.5_GeV || Ecm < 10_GeV ) {
       std::cout << "ProcessSplit: " << " DoDiscrete: dropping particle.." << std::endl;
       p.Delete();
@@ -138,27 +204,53 @@ public:
       // add particles from sibyll to stack
       // link to sibyll stack
       SibStack ss;
-      /*
-	get transformation between Stack-frame and SibStack-frame
-	for EAS Stack-frame is lab. frame, could be different for CRMC-mode
-	the transformation should be derived from the input momenta
-	in general transformation is rotation + boost
-      */
-      const EnergyType proton_mass = 0.93827_GeV;
-      const double gamma  = ( E + proton_mass ) / ( Ecm );
-      const double gambet =  sqrt( E * E - proton_mass * proton_mass ) / Ecm;
 
       // SibStack does not know about momentum yet so we need counter to access momentum array in Sibyll
       int i = -1;
-      for (auto &p: ss){
+      super_stupid::MomentumVector Ptot_final(rootCS, {0.0_newton_second, 0.0_newton_second, 0.0_newton_second});
+      for (auto &psib: ss){
 	++i;
-	//transform to lab. frame, primitve
-	const double en_lab = gambet * s_plist_.p[2][i] + gamma * p.GetEnergy();	
+	//transform energy to lab. frame, primitve
+	// compute beta_vec * p_vec
+	// arbitrary Lorentz transformation based on sibyll routines
+	const auto gammaBetaComponents = gambet.GetComponents();
+	const auto pSibyllComponents = psib.GetMomentum().GetComponents();	
+	EnergyType en_lab = 0. * 1_GeV;	
+	MomentumType p_lab_components[3];
+	en_lab =  psib.GetEnergy() * gamma;
+	EnergyType pnorm = 0. * 1_GeV;
+	for(int j=0; j<3; ++j)
+	  pnorm += ( pSibyllComponents[j] * gammaBetaComponents[j] * si::constants::c ) / ( gamma + 1.);
+	pnorm += psib.GetEnergy();
+	
+	for(int j=0; j<3; ++j){
+	  p_lab_components[j] = pSibyllComponents[j] - (-1) * pnorm * gammaBetaComponents[j] /  si::constants::c;
+	  // cout << "p:" << j << " pSib (GeV/c): " << pSibyllComponents[j] / 1_GeV * si::constants::c
+	  //      << " gb: " << gammaBetaComponents[j] << endl;
+	  en_lab -= (-1) * pSibyllComponents[j] * gammaBetaComponents[j] * si::constants::c;
+	}
+	//	const EnergyType en_lab = psib.GetEnergy()*gamma + gambet * psib.GetMomentum() * si::constants::c );
+	// cout << " en cm  (GeV): " << psib.GetEnergy() / 1_GeV << endl
+	//      << " en lab (GeV): " << en_lab / 1_GeV << endl;
+	// cout << " pz cm  (GeV/c): " << psib.GetMomentum().GetComponents()[2] / 1_GeV * si::constants::c << endl
+	//      << " pz lab (GeV/c): " << p_lab_components[2] / 1_GeV * si::constants::c << endl;
+	
 	// add to corsika stack
 	auto pnew = s.NewParticle();
-	pnew.SetEnergy( en_lab * 1_GeV );
-	pnew.SetPID( process::sibyll::ConvertFromSibyll( p.GetPID() ) );
+	pnew.SetEnergy( en_lab );
+	pnew.SetPID( process::sibyll::ConvertFromSibyll( psib.GetPID() ) );
+	//cout << "momentum sib (cm): " << psib.GetMomentum().GetComponents() / 1_GeV * si::constants::c << endl;
+	
+	corsika::geometry::QuantityVector<momentum_d> p_lab_c{ p_lab_components[0],
+							       p_lab_components[1],
+							       p_lab_components[2]};
+	pnew.SetMomentum( super_stupid::MomentumVector( rootCS, p_lab_c) );
+	//cout << "momentum sib (lab): " << pnew.GetMomentum().GetComponents() / 1_GeV * si::constants::c << endl;
+	//cout << "s_cm  (GeV2): " << (psib.GetEnergy() * psib.GetEnergy() - psib.GetMomentum().squaredNorm() * si::constants::cSquared ) / 1_GeV / 1_GeV << endl;
+	//cout << "s_lab (GeV2): " << (pnew.GetEnergy() * pnew.GetEnergy() - pnew.GetMomentum().squaredNorm() * si::constants::cSquared ) / 1_GeV / 1_GeV << endl;
+	Ptot_final += pnew.GetMomentum();
       }
+      //cout << "tot. momentum final (GeV/c): " << Ptot_final.GetComponents() / 1_GeV * si::constants::c << endl;
     }
     }else
       p.Delete();
@@ -226,6 +318,12 @@ double s_rndm_(int &)
 
 int main(){
 
+  // coordinate system, get global frame of reference
+  CoordinateSystem rootCS = CoordinateSystem::CreateRootCS();
+  
+  QuantityVector<length_d> const coordinates{0_m, 0_m, 0_m};
+  Point pOrig(rootCS, coordinates);    
+  
   stack_inspector::StackInspector<setup::Stack, setup::Trajectory> p0(true);
   ProcessSplit p1;
   const auto sequence = p0 + p1;
@@ -236,8 +334,14 @@ int main(){
 
   stack.Clear();
   auto particle = stack.NewParticle();
-  EnergyType E0 = 100_GeV;
+  EnergyType E0 = 500_GeV;
+  MomentumType P0 = sqrt( E0*E0 - 0.93827_GeV * 0.93827_GeV ) / si::constants::c;
+  auto plab = super_stupid::MomentumVector(rootCS,
+					   0. * 1_GeV / si::constants::c,
+					   0. * 1_GeV / si::constants::c,
+					   P0);
   particle.SetEnergy(E0);
+  particle.SetMomentum(plab);
   particle.SetPID( Code::Proton );
   EAS.Init();
   EAS.Run();

Framework/Cascade/SibStack.h

@@ -10,6 +10,8 @@
 
 using namespace std;
 using namespace corsika::stack;
+using namespace corsika::units;
+using namespace corsika::geometry;
 
 class SibStackData {
   
@@ -19,14 +21,30 @@ class SibStackData {
   void Clear() { s_plist_.np = 0; }
   
   int GetSize() const { return s_plist_.np;  }
+#warning check actual capacity of sibyll stack  
   int GetCapacity() const { return 8000; }
 
   
   void SetId(const int i, const int v) { s_plist_.llist[i]=v; }
-  void SetEnergy(const int i, const double v) { s_plist_.p[3][i]=v;  }
-
+  void SetEnergy(const int i, const EnergyType v) { s_plist_.p[3][i] = v / 1_GeV;  }
+  void SetMomentum(const int i, const super_stupid::MomentumVector& v)
+  {
+    auto tmp = v.GetComponents();
+    for(int idx=0; idx<3; ++idx)      
+      s_plist_.p[idx][i] = tmp[idx] / 1_GeV * si::constants::c;    
+  }
+  
   int GetId(const int i) const { return s_plist_.llist[i]; }
-  double GetEnergy(const int i) const { return s_plist_.p[3][i]; }
+  
+  EnergyType GetEnergy(const int i) const { return s_plist_.p[3][i] * 1_GeV; }
+  
+  super_stupid::MomentumVector GetMomentum(const int i) const
+  {
+    CoordinateSystem rootCS = CoordinateSystem::CreateRootCS();
+    corsika::geometry::QuantityVector<momentum_d> components{ s_plist_.p[0][i] * 1_GeV / si::constants::c , s_plist_.p[1][i] * 1_GeV / si::constants::c, s_plist_.p[2][i] * 1_GeV / si::constants::c};
+    super_stupid::MomentumVector v1(rootCS,components);
+    return v1;
+  }
   
   void Copy(const int i1, const int i2) {
     s_plist_.llist[i1] = s_plist_.llist[i2];
@@ -45,9 +63,10 @@ class ParticleInterface : public ParticleBase<StackIteratorInterface> {
   using ParticleBase<StackIteratorInterface>::GetIndex;
  public:
   void SetEnergy(const double v) { GetStackData().SetEnergy(GetIndex(), v); }
-  double GetEnergy() const { return GetStackData().GetEnergy(GetIndex()); }
+  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()); }
   
 };