new calculation of beta*gamma with better accuracy

Framework/Utilities/COMBoost.cc

@@ -15,6 +15,8 @@
 #include <corsika/utl/COMBoost.h>
 #include <corsika/utl/sgn.h>
 
+#include <cmath>
+
 using namespace corsika::utl;
 using namespace corsika::units::si;
 using namespace corsika::geometry;
@@ -23,16 +25,17 @@ COMBoost::COMBoost(FourVector<HEPEnergyType, Vector<hepmomentum_d>> const& Pproj
                    const HEPMassType massTarget)
     : fCS(Pprojectile.GetSpaceLikeComponents().GetCoordinateSystem()) {
   auto const pProjectile = Pprojectile.GetSpaceLikeComponents();
-  auto const pProjNorm = pProjectile.norm();
+  auto const pProjNormSquared = pProjectile.squaredNorm();
+  auto const pProjNorm = sqrt(pProjNormSquared);
   auto const a = (pProjectile / pProjNorm).GetComponents().eVector;
   auto const a1 = a(0), a2 = a(1);
 
-  auto const s = sgn(a(2));
-  auto const c = 1 / (1 + s * a(2));
+  auto const sign = sgn(a(2));
+  auto const c = 1 / (1 + sign * a(2));
 
   Eigen::Matrix3d A, B;
 
-  if (s > 0) {
+  if (sign > 0) {
     A << 1, 0, -a1,                     // comment to prevent clang-format
         0, 1, -a2,                      // .
         a1, a2, 1;                      // .
@@ -51,15 +54,17 @@ COMBoost::COMBoost(FourVector<HEPEnergyType, Vector<hepmomentum_d>> const& Pproj
 
   fRotation = A + B;
 
-  // calculate boost
-  double const beta = pProjNorm / (Pprojectile.GetTimeLikeComponent() + massTarget);
+  auto const eProjectile = Pprojectile.GetTimeLikeComponent();
+  auto const massProjectileSquared = eProjectile * eProjectile - pProjNormSquared;
+  auto const s =
+      massTarget * massTarget + massProjectileSquared + 2 * eProjectile * massTarget;
 
-  /* Accurracy matters here, beta = 1 - epsilon for ultra-relativistic boosts */
-  double const coshEta = 1 / std::sqrt((1 + beta) * (1 - beta));
-  //~ double const coshEta = 1 / std::sqrt((1-beta*beta));
-  double const sinhEta = -beta * coshEta;
+  auto const sqrtS = sqrt(s);
+  auto const sinhEta = -pProjNorm / sqrtS;
+  auto const coshEta = sqrt(1 + pProjNormSquared / s);
 
-  std::cout << "COMBoost (1-beta)=" << 1 - beta << " gamma=" << coshEta << std::endl;
+  std::cout << "COMBoost (1-beta)=" << 1 - sinhEta / coshEta << " gamma=" << coshEta
+            << std::endl;
   std::cout << "  det = " << fRotation.determinant() - 1 << std::endl;
 
   fBoost << coshEta, sinhEta, sinhEta, coshEta;