Merge branch '107-implement-elasticmodel' of...

Merge branch '107-implement-elasticmodel' of gitlab.ikp.kit.edu:AirShowerPhysics/corsika into 107-implement-elasticmodel

Documentation/Examples/cascade_example.cc

@@ -251,7 +251,6 @@ int main() {
   // setup particle stack, and add primary particle
   setup::Stack stack;
   stack.Clear();
-
   const HEPEnergyType E0 = 100_TeV; // 1_PeV crashes with bad COMboost in second interaction (crash later)
   double theta = 0.;
   double phi = 0.;

Framework/Utilities/COMBoost.cc

@@ -82,8 +82,8 @@ FourVector COMBoost<FourVector>::fromCoM(const FourVector& p) const {
   com << (p.GetTimeLikeComponent() * (1 / 1_GeV)),
       (p.GetSpaceLikeComponents().GetComponents().eVector(2) * (1 / 1_GeV).magnitude());
 
-  std::cout << "COMBoost::fromCoM Ecm=" << p.GetTimeLikeComponent() / 1_GeV << "GeV, "
-            << " pcm=" << p.GetSpaceLikeComponents().GetComponents() / 1_GeV << "GeV"
+  std::cout << "COMBoost::fromCoM Ecm=" << p.GetTimeLikeComponent() / 1_GeV << " GeV, "
+            << " pcm=" << p.GetSpaceLikeComponents().GetComponents().squaredNorm() / 1_GeV << " GeV"
             << std::endl;
 
   auto const boostedZ = fInverseBoost * com;
@@ -94,7 +94,7 @@ FourVector COMBoost<FourVector>::fromCoM(const FourVector& p) const {
   pLab.eVector = fRotation.transpose() * pLab.eVector;
 
   std::cout << "COMBoost::fromCoM --> Elab=" << E_lab / 1_GeV << "GeV, "
-            << " pcm=" << pLab / 1_GeV << "GeV" << std::endl;
+            << " pcm=" << pLab.squaredNorm() / 1_GeV << "GeV" << std::endl;
 
   return FourVector(E_lab, corsika::geometry::Vector(fCS, pLab));
 }

Processes/HadronicElasticModel/HadronicElasticModel.h

@@ -53,29 +53,35 @@ namespace corsika::process::HadronicElasticModel {
     auto B(decltype(units::si::detail::static_pow<2>(units::si::electronvolt)) s) const {
       using namespace corsika::units::constants;
       auto constexpr b_p = 2.3;
-      return (2 * b_p + 2 * b_p + 4 * pow(s * invGeVsq, gfEpsilon) - 4.2) * invGeVsq;
+      auto const result =
+          (2 * b_p + 2 * b_p + 4 * pow(s * invGeVsq, gfEpsilon) - 4.2) * invGeVsq;
+      std::cout << "B(" << s << ") = " << result / invGeVsq << " GeV¯²" << std::endl;
+      return result;
     }
 
     corsika::units::si::CrossSectionType CrossSection(SquaredHEPEnergyType s) const {
       using namespace corsika::units::si;
       using namespace corsika::units::constants;
       // assuming every target behaves like a proton, fX and fY are universal
-      CrossSectionType const sigmaTot =
+      CrossSectionType const sigmaTotal =
           fX * pow(s * invGeVsq, gfEpsilon) + fY * pow(s * invGeVsq, -gfEta);
 
       // according to Schuler & Sjöstrand, PRD 49, 2257 (1994)
       // (we ignore rho because rho^2 is just ~2 %)
       auto const sigmaElastic =
-          units::si::detail::static_pow<2>(sigmaTot) /
+          units::si::detail::static_pow<2>(sigmaTotal) /
           (16 * M_PI * ConvertHEPToSI<CrossSectionType::dimension_type>(B(s)));
+
+      std::cout << "HEM sigmaTot = " << sigmaTotal / 1_mbarn << " mb" << std::endl;
+      std::cout << "HEM sigmaElastic = " << sigmaElastic / 1_mbarn << " mb" << std::endl;
       return sigmaElastic;
     }
 
   public:
     HadronicElasticInteraction(corsika::environment::Environment const&,
-                               // x & y values taken from Pythia8 for pp collisions
-                               units::si::CrossSectionType x = 0.217 * units::si::barn,
-                               units::si::CrossSectionType y = 0.5608 * units::si::barn);
+                               // x & y values taken from DL for pp collisions
+                               units::si::CrossSectionType x = 0.0217 * units::si::barn,
+                               units::si::CrossSectionType y = 0.05608 * units::si::barn);
     void Init();
 
     template <typename Particle, typename Track>