diff --git a/Framework/Utilities/COMBoost.cc b/Framework/Utilities/COMBoost.cc
index eedabc6f99b5f5df2cda9fa2354ca8e474f054d2..abdb7a48b88f733aba4c03dc4d8521a741aa27b5 100644
--- a/Framework/Utilities/COMBoost.cc
+++ b/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_; }
diff --git a/Framework/Utilities/COMBoost.h b/Framework/Utilities/COMBoost.h
index 48532531e5caa981f1d01daf258af5d737fa1216..e25d74e828c9140815dd98ec426871d8055a6415 100644
--- a/Framework/Utilities/COMBoost.h
+++ b/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
diff --git a/Framework/Utilities/testCOMBoost.cc b/Framework/Utilities/testCOMBoost.cc
index 84226cf86dbd4cda644168f21224784407646183..b665c5fc004e938a370d2befd3f050a76a57ca4d 100644
--- a/Framework/Utilities/testCOMBoost.cc
+++ b/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()));
+ }
}