use CoordinateSystem::RotateToZ in COMBoost and add rest frame boost

Framework/Utilities/COMBoost.cc

@@ -23,36 +23,11 @@ using namespace corsika::geometry;
 
 COMBoost::COMBoost(FourVector<HEPEnergyType, Vector<hepmomentum_d>> const& Pprojectile,
                    const HEPMassType massTarget)
-    : fCS(Pprojectile.GetSpaceLikeComponents().GetCoordinateSystem()) {
+    : originalCS_{Pprojectile.GetSpaceLikeComponents().GetCoordinateSystem()}
+    , rotatedCS_{originalCS_.RotateToZ(Pprojectile.GetSpaceLikeComponents())} {
   auto const pProjectile = Pprojectile.GetSpaceLikeComponents();
   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 sign = sgn(a(2));
-  auto const c = 1 / (1 + sign * a(2));
-
-  Eigen::Matrix3d A, B;
-
-  if (sign > 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; // .
-
-  } 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;
 
   auto const eProjectile = Pprojectile.GetTimeLikeComponent();
   auto const massProjectileSquared = eProjectile * eProjectile - pProjNormSquared;
@@ -63,15 +38,27 @@ COMBoost::COMBoost(FourVector<HEPEnergyType, Vector<hepmomentum_d>> const& Pproj
   auto const sinhEta = -pProjNorm / sqrtS;
   auto const coshEta = sqrt(1 + pProjNormSquared / s);
 
+  setBoost(coshEta, sinhEta);
+
   std::cout << "COMBoost (1-beta)=" << 1 - sinhEta / coshEta << " gamma=" << coshEta
             << std::endl;
-  std::cout << "  det = " << fRotation.determinant() - 1 << std::endl;
+  std::cout << "  det = " << boost_.determinant() - 1 << std::endl;
+}
 
-  fBoost << coshEta, sinhEta, sinhEta, coshEta;
+COMBoost::COMBoost(geometry::Vector<units::si::hepmomentum_d> const& momentum,
+                   units::si::HEPEnergyType mass)
+    : originalCS_{momentum.GetCoordinateSystem()}
+    , rotatedCS_{originalCS_.RotateToZ(momentum)} {
+  auto const squaredNorm = momentum.squaredNorm();
+  auto const norm = sqrt(squaredNorm);
+  auto const sinhEta = -norm / mass;
+  auto const coshEta = sqrt(1 + squaredNorm / (mass * mass));
+  setBoost(coshEta, sinhEta);
+}
 
-  fInverseBoost << coshEta, -sinhEta, -sinhEta, coshEta;
+void COMBoost::setBoost(double coshEta, double sinhEta) {
+  boost_ << coshEta, sinhEta, sinhEta, coshEta;
+  inverseBoost_ << coshEta, -sinhEta, -sinhEta, coshEta;
 }
 
-/*
-  Here we instantiate all physically meaningful versions of COMBoost
- */
+CoordinateSystem const& COMBoost::GetRotatedCS() const { return rotatedCS_; }

Framework/Utilities/COMBoost.h

@@ -25,72 +25,77 @@ namespace corsika::utl {
    */
 
   class COMBoost {
-    Eigen::Matrix3d fRotation;
-    Eigen::Matrix2d fBoost, fInverseBoost;
-    corsika::geometry::CoordinateSystem const& fCS;
+    Eigen::Matrix2d boost_, inverseBoost_;
+    corsika::geometry::CoordinateSystem const &originalCS_, rotatedCS_;
+
+    void setBoost(double coshEta, double sinhEta);
 
   public:
-    //! construct a COMBoost given four-vector of prjectile and mass of target
+    //! construct a COMBoost given four-vector of projectile and mass of target
     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; }
+    //! construct a COMBoost to boost into the rest frame given a 3-momentum and mass
+    COMBoost(geometry::Vector<units::si::hepmomentum_d> const& momentum,
+             units::si::HEPEnergyType mass);
 
     //! transforms a 4-momentum from lab frame to the center-of-mass frame
     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;
+      auto pComponents = p.GetSpaceLikeComponents().GetComponents(rotatedCS_);
+      Eigen::Vector3d eVecRotated = pComponents.eVector;
       Eigen::Vector2d lab;
 
       lab << (p.GetTimeLikeComponent() * (1 / 1_GeV)),
           (eVecRotated(2) * (1 / 1_GeV).magnitude());
 
-      auto const boostedZ = fBoost * lab;
+      auto const boostedZ = boost_ * 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));
+      return FourVector(
+          E_CoM, corsika::geometry::Vector<hepmomentum_d>(rotatedCS_, eVecRotated));
     }
 
     //! transforms a 4-momentum from the center-of-mass frame back to lab frame
     template <typename FourVector>
     FourVector fromCoM(const FourVector& p) const {
       using namespace corsika::units::si;
+      auto pCM = p.GetSpaceLikeComponents().GetComponents(rotatedCS_);
+      auto const Ecm = p.GetTimeLikeComponent();
+
       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
+      com << (Ecm * (1 / 1_GeV)), (pCM.eVector(2) * (1 / 1_GeV).magnitude());
+
+      std::cout << "COMBoost::fromCoM Ecm=" << Ecm / 1_GeV << " GeV, "
+                << " pcm = " << pCM / 1_GeV << " (norm = " << pCM.norm() / 1_GeV
                 << " GeV), invariant mass = " << p.GetNorm() / 1_GeV << " GeV"
                 << std::endl;
 
-      auto const boostedZ = fInverseBoost * com;
+      auto const boostedZ = inverseBoost_ * 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;
+      pCM.eVector(2) = boostedZ(1) * (1_GeV).magnitude();
 
-      FourVector f(E_lab, corsika::geometry::Vector(fCS, pLab));
+      geometry::Vector<typename decltype(pCM)::dimension> pLab{rotatedCS_, pCM};
+      pLab.rebase(originalCS_);
+
+      FourVector f(E_lab, pLab);
 
       std::cout << "COMBoost::fromCoM --> Elab=" << E_lab / 1_GeV << "GeV, "
-                << " pcm = " << pLab / 1_GeV << " (norm =" << pLab.norm() / 1_GeV
+                << " plab = " << pLab.GetComponents() << " (norm =" << pLab.norm() / 1_GeV
                 << " GeV), invariant mass = " << f.GetNorm() / 1_GeV << " GeV"
                 << std::endl;
 
       return f;
     }
+
+    geometry::CoordinateSystem const& GetRotatedCS() const;
   };
 } // namespace corsika::utl
 

Framework/Utilities/testCOMBoost.cc

@@ -42,87 +42,6 @@ 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;
@@ -266,4 +185,31 @@ TEST_CASE("boosts") {
         PprojCoM.GetSpaceLikeComponents() + PtargCoM.GetSpaceLikeComponents();
     CHECK(sumPCoM.norm() / P0 == Approx(0).margin(absMargin)); // MAKE RELATIVE CHECK
   }
+
+  SECTION("rest frame") {
+    HEPMassType const projectileMass = 1_GeV;
+    HEPMomentumType const P0 = 1_TeV;
+    Vector<hepmomentum_d> pProjectileLab{rootCS, {0_GeV, P0, 0_GeV}};
+    HEPEnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
+    const FourVector PprojLab(eProjectileLab, pProjectileLab);
+
+    COMBoost boostRest(pProjectileLab, projectileMass);
+    auto const& csPrime = boostRest.GetRotatedCS();
+    FourVector const rest4Mom = boostRest.toCoM(PprojLab);
+
+    CHECK(rest4Mom.GetTimeLikeComponent() / 1_GeV == Approx(projectileMass / 1_GeV));
+    CHECK(rest4Mom.GetSpaceLikeComponents().norm() / 1_GeV ==
+          Approx(0).margin(absMargin));
+
+    FourVector const a{0_eV, Vector{csPrime, 0_eV, 5_GeV, 0_eV}};
+    FourVector const b{0_eV, Vector{rootCS, 3_GeV, 0_eV, 0_eV}};
+    auto const aLab = boostRest.fromCoM(a);
+    auto const bLab = boostRest.fromCoM(b);
+
+    CHECK(aLab.GetNorm() / a.GetNorm() == Approx(1));
+    CHECK(aLab.GetSpaceLikeComponents().GetComponents(csPrime)[1].magnitude() ==
+          Approx((5_GeV).magnitude()));
+    CHECK(bLab.GetSpaceLikeComponents().GetComponents(rootCS)[0].magnitude() ==
+          Approx((3_GeV).magnitude()));
+  }
 }