Merge branch 'sibyll' of gitlab.ikp.kit.edu:AirShowerPhysics/corsika into sibyll

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,35 +86,52 @@ 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
+    /*
+      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();
-    double Ecm = sqrt( 2. * E * 0.93827_GeV ) / 1_GeV ;
+    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 << std::endl;
-    if (E < 8.5_GeV || Ecm < 10. ) {
+    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();
       fCount++;
     } else {
-      //p.SetEnergy(E / 2);
-      //s.NewParticle().SetEnergy(E / 2);
-      
+      // Sibyll does not know about units..
+      double sqs = Ecm / 1_GeV;
       // running sibyll
-      sibyll_( kBeam, kTarget, Ecm);
+      sibyll_( kBeam, kTarget, sqs);
+      // running decays
+      //decsib_();
+      // print final state
       int print_unit = 6;
       sib_list_( print_unit );
       
       // delete current particle
       p.Delete();
+
+      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;
       
       // add particles from sibyll to stack
       for(int i=0; i<s_plist_.np; ++i){
-	
-	s.NewParticle().SetEnergy( s_plist_.p[3][i] * 1_GeV );
+	//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 );
       }
     }
   }
@@ -102,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_();
@@ -121,6 +161,8 @@ public:
     
     //initialize Sibyll
     sibyll_ini_();
+
+    // set particles stable / unstable
   }
   
   int GetCount() { return fCount; }

Framework/Units/PhysicalConstants.h

@@ -54,7 +54,7 @@ namespace corsika::units::si::constants {
   // unified atomic mass unit
   constexpr quantity<mass_d> u{Rep(1.6605402e-27L) * kilogram};
 
-  // millibarn
+  // barn
   constexpr quantity<area_d> barn{Rep(1.e-28L) * meter * meter};
   
   // etc.

Framework/Units/PhysicalUnits.h

@@ -24,6 +24,31 @@ namespace phys {
   }   // namespace units
 } // namespace phys
 
+namespace phys {
+  namespace units {
+    namespace literals {
+      QUANTITY_DEFINE_SCALING_LITERALS(barn, area_d,
+                                       magnitude(corsika::units::si::constants::barn))
+
+      // phys::units::quantity<energy_d/mass_d> Joule2Kg = c2; // 1_Joule / 1_kg;
+
+    } // namespace literals
+  }   // namespace units
+}
+
+namespace phys {
+  namespace units {
+    namespace literals {
+      QUANTITY_DEFINE_SCALING_LITERALS(meter, length_d,
+                                       magnitude(corsika::units::si::constants::meter))
+
+      // phys::units::quantity<energy_d/mass_d> Joule2Kg = c2; // 1_Joule / 1_kg;
+
+    } // namespace literals
+  }   // namespace units
+}
+
+
 namespace corsika::units::si {
   using namespace phys::units;
   using namespace phys::units::literals;

Framework/Units/testUnits.cc

@@ -54,7 +54,8 @@ TEST_CASE("PhysicalUnits", "[Units]") {
     REQUIRE(1_mol / 1_amol == Approx(1e18));
     REQUIRE(1_K / 1_zK == Approx(1e21));
     REQUIRE(1_K / 1_yK == Approx(1e24));
-
+    //    REQUIRE(1_barn / 1_mbarn == Approx(1e3));
+    
     REQUIRE(1_A / 1_hA == Approx(1e-2));
     REQUIRE(1_m / 1_km == Approx(1e-3));
     REQUIRE(1_m / 1_Mm == Approx(1e-6));

Processes/StackInspector/StackInspector.cc

@@ -38,14 +38,15 @@ process::EProcessReturn StackInspector<Stack, Trajectory>::DoContinuous(Particle
   // std::cout << "generation  " << countStep << std::endl;
   [[maybe_unused]] int i = 0;
   EnergyType Etot = 0_GeV;
+
   for (auto& iterP : s) {
     EnergyType E = iterP.GetEnergy();
     Etot += E;
-    // std::cout << " particle data: " << i++ << ", id=" << iterP << " | " << std::endl;
+    std::cout << " particle data: " << i++ << ", id=" << iterP.GetPID()<< ", en=" << iterP.GetEnergy() / 1_GeV << " | " << std::endl;
   }
   countStep++;
   // cout << "#=" << countStep << " " << s.GetSize() << " " << Etot/1_GeV << endl;
-  cout << countStep << " " << s.GetSize() << " " << Etot / 1_GeV << " " << endl;
+  cout << "call: " << countStep << " stack size: " << s.GetSize() << " energy: " << Etot / 1_GeV << endl;
   return EProcessReturn::eOk;
 }