fixed inaccuracy in COMBoost

Framework/Utilities/COMBoost.cc

@@ -10,41 +10,53 @@
  */
 
 #include <corsika/geometry/CoordinateSystem.h>
+#include <corsika/geometry/FourVector.h>
 #include <corsika/geometry/Vector.h>
 #include <corsika/units/PhysicalUnits.h>
 #include <corsika/utl/COMBoost.h>
 
 using namespace corsika::utl;
 using namespace corsika::units::si;
+using namespace corsika::geometry;
+
+//! sign function without branches
+template <typename T>
+static int sgn(T val) {
+  return (T(0) < val) - (val < T(0));
+}
 
-template <typename FourVector>
-COMBoost<FourVector>::COMBoost(const FourVector& Pprojectile,
-                               const HEPMassType massTarget)
-    : fRotation(Eigen::Matrix3d::Identity())
-    , fCS(Pprojectile.GetSpaceLikeComponents().GetCoordinateSystem()) {
-  // calculate matrix for rotating pProjectile to z-axis first
+COMBoost::COMBoost(FourVector<HEPEnergyType, Vector<hepmomentum_d>> const& Pprojectile,
+                   const HEPMassType massTarget)
+    : fCS(Pprojectile.GetSpaceLikeComponents().GetCoordinateSystem()) {
   auto const pProjectile = Pprojectile.GetSpaceLikeComponents();
   auto const pProjNorm = pProjectile.norm();
   auto const a = (pProjectile / pProjNorm).GetComponents().eVector;
+  auto const a1 = a(0), a2 = a(1);
 
-  if (a(0) == 0 && a(1) == 0) {
-    if (a(2) < 0) {
-      // if pProjectile ~ (0, 0, -1), the standard formula for the rotation matrix breaks
-      // down but we can easily define a suitable rotation manually. We just need some
-      // SO(3) matrix that reverses the z-axis and I like this one:
+  auto const s = sgn(a(2));
+  auto const c = 1 / (1 + s * a(2));
 
-      fRotation << 1, 0, 0, 0, -1, 0, 0, 0, -1;
-    }
-  } else {
-    Eigen::Vector3d const b{0, 0, 1};
-    auto const v = a.cross(b);
+  Eigen::Matrix3d A, B;
 
-    Eigen::Matrix3d vHat;
-    vHat << 0, -v(2), v(1), v(2), 0, -v(0), -v(1), v(0), 0;
+  if (s > 0) {
+    A << 1, 0, -a1,                     // comment to prevent clang-format
+        0, 1, -a2,                      // .
+        a1, a2, 1;                      // .
+    B << -a1 * a1 * c, -a1 * a2 * c, 0, // .
+        -a1 * a2 * c, -a2 * a2 * c, 0,  // .
+        0, 0, -(a1 * a1 + a2 * a2) * c; // .
 
-    fRotation += vHat + vHat * vHat / (1 + a.dot(b));
+  } else {
+    A << 1, 0, a1,                      // comment to prevent clang-format
+        0, -1, -a2,                     // .
+        a1, a2, -1;                     // .
+    B << -a1 * a1 * c, -a1 * a2 * c, 0, // .
+        +a1 * a2 * c, +a2 * a2 * c, 0,  // .
+        0, 0, (a1 * a1 + a2 * a2) * c;  // .
   }
 
+  fRotation = A + B;
+
   // calculate boost
   double const beta = pProjNorm / (Pprojectile.GetTimeLikeComponent() + massTarget);
 
@@ -54,68 +66,13 @@ COMBoost<FourVector>::COMBoost(const FourVector& Pprojectile,
   double const sinhEta = -beta * coshEta;
 
   std::cout << "COMBoost (1-beta)=" << 1 - beta << " gamma=" << coshEta << std::endl;
+  std::cout << "  det = " << fRotation.determinant() - 1 << std::endl;
 
   fBoost << coshEta, sinhEta, sinhEta, coshEta;
 
   fInverseBoost << coshEta, -sinhEta, -sinhEta, coshEta;
 }
 
-template <typename FourVector>
-FourVector COMBoost<FourVector>::toCoM(const FourVector& p) const {
-  auto pComponents = p.GetSpaceLikeComponents().GetComponents(fCS);
-  Eigen::Vector3d eVecRotated = fRotation * pComponents.eVector;
-  Eigen::Vector2d lab;
-
-  lab << (p.GetTimeLikeComponent() * (1 / 1_GeV)),
-      (eVecRotated(2) * (1 / 1_GeV).magnitude());
-
-  auto const boostedZ = fBoost * lab;
-  auto const E_CoM = boostedZ(0) * 1_GeV;
-
-  eVecRotated(2) = boostedZ(1) * (1_GeV).magnitude();
-
-  return FourVector(E_CoM, corsika::geometry::Vector<hepmomentum_d>(fCS, eVecRotated));
-}
-
-template <typename FourVector>
-FourVector COMBoost<FourVector>::fromCoM(const FourVector& p) const {
-  Eigen::Vector2d com;
-  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().squaredNorm() / 1_GeV
-            << " GeV" << std::endl;
-
-  auto const boostedZ = fInverseBoost * com;
-  auto const E_lab = boostedZ(0) * 1_GeV;
-
-  auto pLab = p.GetSpaceLikeComponents().GetComponents();
-  pLab.eVector(2) = boostedZ(1) * (1_GeV).magnitude();
-  pLab.eVector = fRotation.transpose() * pLab.eVector;
-
-  std::cout << "COMBoost::fromCoM --> Elab=" << E_lab / 1_GeV << "GeV, "
-            << " pcm=" << pLab.squaredNorm() / 1_GeV << "GeV" << std::endl;
-
-  return FourVector(E_lab, corsika::geometry::Vector(fCS, pLab));
-}
-
 /*
   Here we instantiate all physically meaningful versions of COMBoost
  */
-
-#include <corsika/geometry/FourVector.h>
-
-namespace corsika::utl {
-
-  using corsika::geometry::FourVector;
-  using corsika::geometry::Vector;
-  using namespace corsika::units;
-
-  // for normal HEP energy/momentum units in [GeV]
-  template class COMBoost<FourVector<HEPEnergyType, Vector<hepmomentum_d>>>;
-  // template class COMBoost<FourVector<HEPEnergyType&, Vector<hepmomentum_d>&>>;
-
-  // template class COMBoost<FourVector<TimeType, Vector<length_d>>>;
-  // template class COMBoost<FourVector<TimeType&, Vector<length_d>&>>;
-} // namespace corsika::utl

Framework/Utilities/COMBoost.h

@@ -13,6 +13,7 @@
 #define _include_corsika_utilties_comboost_h_
 
 #include <corsika/geometry/CoordinateSystem.h>
+#include <corsika/geometry/FourVector.h>
 #include <corsika/units/PhysicalUnits.h>
 
 #include <Eigen/Dense>
@@ -24,7 +25,6 @@ namespace corsika::utl {
      referenence frames, using FourVectors.
    */
 
-  template <typename FourVector>
   class COMBoost {
     Eigen::Matrix3d fRotation;
     Eigen::Matrix2d fBoost, fInverseBoost;
@@ -32,14 +32,66 @@ namespace corsika::utl {
 
   public:
     //! construct a COMBoost given four-vector of prjectile and mass of target
-    COMBoost(const FourVector& Pprojectile,
-             const corsika::units::si::HEPEnergyType massTarget);
+    COMBoost(
+        const corsika::geometry::FourVector<
+            corsika::units::si::HEPEnergyType,
+            corsika::geometry::Vector<corsika::units::si::hepmomentum_d>>& Pprojectile,
+        const corsika::units::si::HEPEnergyType massTarget);
+
+    auto const& GetRotationMatrix() const { return fRotation; }
 
     //! transforms a 4-momentum from lab frame to the center-of-mass frame
-    FourVector toCoM(const FourVector& p) const;
+    template <typename FourVector>
+    FourVector toCoM(const FourVector& p) const {
+      using namespace corsika::units::si;
+      auto pComponents = p.GetSpaceLikeComponents().GetComponents(fCS);
+      Eigen::Vector3d eVecRotated = fRotation * pComponents.eVector;
+      Eigen::Vector2d lab;
+
+      lab << (p.GetTimeLikeComponent() * (1 / 1_GeV)),
+          (eVecRotated(2) * (1 / 1_GeV).magnitude());
+
+      auto const boostedZ = fBoost * lab;
+      auto const E_CoM = boostedZ(0) * 1_GeV;
+
+      eVecRotated(2) = boostedZ(1) * (1_GeV).magnitude();
+
+      return FourVector(E_CoM,
+                        corsika::geometry::Vector<hepmomentum_d>(fCS, eVecRotated));
+    }
 
     //! transforms a 4-momentum from the center-of-mass frame back to lab frame
-    FourVector fromCoM(const FourVector& pCoM) const;
+    template <typename FourVector>
+    FourVector fromCoM(const FourVector& p) const {
+      using namespace corsika::units::si;
+      Eigen::Vector2d com;
+      com << (p.GetTimeLikeComponent() * (1 / 1_GeV)),
+          (p.GetSpaceLikeComponents().GetComponents().eVector(2) *
+           (1 / 1_GeV).magnitude());
+
+      auto const plab = p.GetSpaceLikeComponents().GetComponents();
+      std::cout << "COMBoost::fromCoM Ecm=" << p.GetTimeLikeComponent() / 1_GeV
+                << " GeV, "
+                << " pcm = " << plab / 1_GeV << " (norm = " << plab.norm() / 1_GeV
+                << " GeV), invariant mass = " << p.GetNorm() / 1_GeV << " GeV"
+                << std::endl;
+
+      auto const boostedZ = fInverseBoost * com;
+      auto const E_lab = boostedZ(0) * 1_GeV;
+
+      auto pLab = p.GetSpaceLikeComponents().GetComponents();
+      pLab.eVector(2) = boostedZ(1) * (1_GeV).magnitude();
+      pLab.eVector = fRotation.transpose() * pLab.eVector;
+
+      FourVector f(E_lab, corsika::geometry::Vector(fCS, pLab));
+
+      std::cout << "COMBoost::fromCoM --> Elab=" << E_lab / 1_GeV << "GeV, "
+                << " pcm = " << pLab / 1_GeV << " (norm =" << pLab.norm() / 1_GeV
+                << " GeV), invariant mass = " << f.GetNorm() / 1_GeV << " GeV"
+                << std::endl;
+
+      return f;
+    }
   };
 } // namespace corsika::utl
 

Framework/Utilities/testCOMBoost.cc

@@ -18,6 +18,9 @@
 #include <corsika/geometry/Vector.h>
 #include <corsika/units/PhysicalUnits.h>
 #include <corsika/utl/COMBoost.h>
+
+#include <Eigen/Dense>
+
 #include <iostream>
 
 using namespace corsika::geometry;
@@ -28,21 +31,103 @@ using corsika::units::constants::cSquared;
 
 double constexpr absMargin = 1e-6;
 
-TEST_CASE("boosts") {
-  CoordinateSystem& rootCS =
+CoordinateSystem const& rootCS =
       RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
-
-  // helper function for energy-momentum
+      
+ // helper function for energy-momentum
   // relativistic energy
-  auto energy = [](HEPMassType m, Vector<hepmomentum_d> const& p) {
+  auto const energy = [](HEPMassType m, Vector<hepmomentum_d> const& p) {
     return sqrt(m * m + p.squaredNorm());
   };
 
   // helper function for mandelstam-s
-  auto s = [](HEPEnergyType E, QuantityVector<hepmomentum_d> const& p) {
+  auto const s = [](HEPEnergyType E, QuantityVector<hepmomentum_d> const& p) {
     return E * E - p.squaredNorm();
   };
 
+
+TEST_CASE("rotation") {
+    // define projectile kinematics in lab frame
+    HEPMassType const projectileMass = 1_GeV;
+    HEPMassType const targetMass = 1.0e300_eV;
+    Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, 0_PeV, 1_GeV}};
+    HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
+    const FourVector PprojLab(eProjectileLab, pProjectileLab);
+    
+    Eigen::Vector3d e1, e2, e3;
+        e1 << 1, 0, 0;
+        e2 << 0, 1, 0;
+        e3 << 0, 0, 1;
+
+    // define boost to com frame
+    SECTION("pos. z-axis") {
+        COMBoost boost({eProjectileLab, {rootCS, {0_GeV, 0_GeV, 1_GeV}}}, targetMass);
+        auto const& rot = boost.GetRotationMatrix();
+        
+        CHECK((rot * e3 - e3).norm() == Approx(0).margin(absMargin));
+        CHECK((rot * e1).norm() == Approx(1));
+        CHECK((rot * e2).norm() == Approx(1));
+        CHECK((rot * e3).norm() == Approx(1));
+        CHECK(rot.determinant() == Approx(1));
+    }
+    
+    SECTION("y-axis in upper half") {
+        COMBoost boost({eProjectileLab, {rootCS, {0_GeV, 1_GeV, 1_meV}}}, targetMass);
+        auto const& rot = boost.GetRotationMatrix();
+        
+        CHECK((rot * e2 - e3).norm() == Approx(0).margin(absMargin));
+        CHECK((rot * e1).norm() == Approx(1));
+        CHECK((rot * e2).norm() == Approx(1));
+        CHECK((rot * e3).norm() == Approx(1));
+        CHECK(rot.determinant() == Approx(1));
+    }
+    
+    SECTION("x-axis in upper half") {
+        COMBoost boost({eProjectileLab, {rootCS, {1_GeV, 0_GeV, 1_meV}}}, targetMass);
+        auto const& rot = boost.GetRotationMatrix();
+                
+        CHECK((rot * e1 - e3).norm() == Approx(0).margin(absMargin));
+        CHECK((rot * e1).norm() == Approx(1));
+        CHECK((rot * e2).norm() == Approx(1));
+        CHECK((rot * e3).norm() == Approx(1));
+        CHECK(rot.determinant() == Approx(1));
+    }
+    
+    SECTION("neg. z-axis") {
+        COMBoost boost({eProjectileLab, {rootCS, {0_GeV, 0_GeV, -1_GeV}}}, targetMass);
+        auto const& rot = boost.GetRotationMatrix();
+        
+        CHECK((rot * (-e3) - e3).norm() == Approx(0).margin(absMargin));
+        CHECK((rot * e1).norm() == Approx(1));
+        CHECK((rot * e2).norm() == Approx(1));
+        CHECK((rot * e3).norm() == Approx(1));
+        CHECK(rot.determinant() == Approx(1));
+    }
+    
+    SECTION("x-axis lower half") {
+        COMBoost boost({eProjectileLab, {rootCS, {1_GeV, 0_GeV, -1_meV}}}, targetMass);
+        auto const& rot = boost.GetRotationMatrix();
+        
+        CHECK((rot * e1 - e3).norm() == Approx(0).margin(absMargin));
+        CHECK((rot * e1).norm() == Approx(1));
+        CHECK((rot * e2).norm() == Approx(1));
+        CHECK((rot * e3).norm() == Approx(1));
+        CHECK(rot.determinant() == Approx(1));
+    }
+    
+    SECTION("y-axis lower half") {
+        COMBoost boost({eProjectileLab, {rootCS, {0_GeV, 1_GeV, -1_meV}}}, targetMass);
+        auto const& rot = boost.GetRotationMatrix();
+        
+        CHECK((rot * e2 - e3).norm() == Approx(0).margin(absMargin));
+        CHECK((rot * e1).norm() == Approx(1));
+        CHECK((rot * e2).norm() == Approx(1));
+        CHECK((rot * e3).norm() == Approx(1));
+        CHECK(rot.determinant() == Approx(1));
+    }
+}
+
+TEST_CASE("boosts") {
   // define target kinematics in lab frame
   HEPMassType const targetMass = 1_GeV;
   Vector<hepmomentum_d> pTargetLab{rootCS, {0_eV, 0_eV, 0_eV}};