From 007fbbb3ce051124da2e946f445a5b80c24cfff6 Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@tu-dortmund.de>
Date: Sat, 29 Dec 2018 00:20:17 +0100
Subject: [PATCH] rotation in COMBoost for (0,0,-z) case
---
Framework/Utilities/COMBoost.cc | 18 ++++---
Framework/Utilities/testCOMBoost.cc | 79 +++++++++++++++--------------
2 files changed, 54 insertions(+), 43 deletions(-)
diff --git a/Framework/Utilities/COMBoost.cc b/Framework/Utilities/COMBoost.cc
index 4e8a9688..42d5995d 100644
--- a/Framework/Utilities/COMBoost.cc
+++ b/Framework/Utilities/COMBoost.cc
@@ -18,17 +18,23 @@ COMBoost::COMBoost(EnergyType eProjectile, COMBoost::MomentumVector const& pProj
: fRotation(Eigen::Matrix3d::Identity())
, fCS(pProjectile.GetCoordinateSystem()) {
// calculate matrix for rotating pProjectile to z-axis first
- // TODO: handle the case when pProjectile ~ (0, 0, -1)
auto const pProjNorm = pProjectile.norm();
auto const a = (pProjectile / pProjNorm).GetComponents().eVector;
- Eigen::Vector3d const b{0, 0, 1};
- auto const v = a.cross(b);
+ if (a(0) == 0 && a(1) == 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:
+ 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 vHat;
- vHat << 0, -v(2), v(1), v(2), 0, -v(0), -v(1), v(0), 0;
+ Eigen::Matrix3d vHat;
+ vHat << 0, -v(2), v(1), v(2), 0, -v(0), -v(1), v(0), 0;
- fRotation += vHat + vHat * vHat / (1 + a.dot(b));
+ fRotation += vHat + vHat * vHat / (1 + a.dot(b));
+ }
// calculate boost
double const x = pProjNorm * units::constants::c /
diff --git a/Framework/Utilities/testCOMBoost.cc b/Framework/Utilities/testCOMBoost.cc
index faac8c6d..3af25856 100644
--- a/Framework/Utilities/testCOMBoost.cc
+++ b/Framework/Utilities/testCOMBoost.cc
@@ -41,41 +41,46 @@ TEST_CASE("boosts") {
// define projectile kinematics in lab frame
MassType const projectileMass = 1._GeV / cSquared;
- Vector<momentum_d> pProjectileLab{rootCS, {0_GeV / c, 1_PeV / c, 0_GeV / c}};
- EnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
-
- // define target kinematics in lab frame
- MassType const targetMass = 1_GeV / cSquared;
- Vector<momentum_d> pTargetLab{rootCS, {0_Ns, 0_Ns, 0_Ns}};
- EnergyType const eTargetLab = energy(targetMass, pTargetLab);
-
- // define boost to com frame
- COMBoost boost(eProjectileLab, pProjectileLab, targetMass);
-
- // boost projecticle
- auto const [eProjectileCoM, pProjectileCoM] =
- boost.toCoM(eProjectileLab, pProjectileLab);
-
- // boost target
- auto const [eTargetCoM, pTargetCoM] = boost.toCoM(eTargetLab, pTargetLab);
-
- // sum of momenta in CoM, should be 0
- auto const sumPCoM = pProjectileCoM + pTargetCoM;
- CHECK(sumPCoM[2] / (1_GeV / c) == Approx(0).margin(absMargin));
-
- // mandelstam-s should be invariant under transformation
- CHECK(s(eProjectileLab + eTargetLab,
- pProjectileLab.GetComponents() + pTargetLab.GetComponents()) /
- (1_GeV / c) / (1_GeV / c) ==
- Approx(s(eProjectileCoM + eTargetCoM, pProjectileCoM + pTargetCoM) / (1_GeV / c) /
- (1_GeV / c)));
-
- // boost back...
- auto const [eProjectileBack, pProjectileBack] =
- boost.fromCoM(eProjectileCoM, pProjectileCoM);
-
- // ...should yield original values before the boosts
- CHECK(eProjectileBack / eProjectileLab == Approx(1));
- CHECK((pProjectileBack - pProjectileLab).norm() / pProjectileLab.norm() ==
- Approx(0).margin(absMargin));
+ std::vector<Vector<momentum_d>> labProjectiles{
+ {rootCS, {0_GeV / c, 1_PeV / c, 0_GeV / c}}, // standard case
+ {rootCS, {0_GeV / c, 0_GeV / c, -1_GeV / c}}}; // "special" case
+
+ for (auto const& pProjectileLab : labProjectiles) {
+ EnergyType const eProjectileLab = energy(projectileMass, pProjectileLab);
+
+ // define target kinematics in lab frame
+ MassType const targetMass = 1_GeV / cSquared;
+ Vector<momentum_d> pTargetLab{rootCS, {0_Ns, 0_Ns, 0_Ns}};
+ EnergyType const eTargetLab = energy(targetMass, pTargetLab);
+
+ // define boost to com frame
+ COMBoost boost(eProjectileLab, pProjectileLab, targetMass);
+
+ // boost projecticle
+ auto const [eProjectileCoM, pProjectileCoM] =
+ boost.toCoM(eProjectileLab, pProjectileLab);
+
+ // boost target
+ auto const [eTargetCoM, pTargetCoM] = boost.toCoM(eTargetLab, pTargetLab);
+
+ // sum of momenta in CoM, should be 0
+ auto const sumPCoM = pProjectileCoM + pTargetCoM;
+ CHECK(sumPCoM[2] / (1_GeV / c) == Approx(0).margin(absMargin));
+
+ // mandelstam-s should be invariant under transformation
+ CHECK(s(eProjectileLab + eTargetLab,
+ pProjectileLab.GetComponents() + pTargetLab.GetComponents()) /
+ (1_GeV / c) / (1_GeV / c) ==
+ Approx(s(eProjectileCoM + eTargetCoM, pProjectileCoM + pTargetCoM) /
+ (1_GeV / c) / (1_GeV / c)));
+
+ // boost back...
+ auto const [eProjectileBack, pProjectileBack] =
+ boost.fromCoM(eProjectileCoM, pProjectileCoM);
+
+ // ...should yield original values before the boosts
+ CHECK(eProjectileBack / eProjectileLab == Approx(1));
+ CHECK((pProjectileBack - pProjectileLab).norm() / pProjectileLab.norm() ==
+ Approx(0).margin(absMargin));
+ }
}
--
GitLab