diff --git a/Framework/Utilities/COMBoost.cc b/Framework/Utilities/COMBoost.cc
index 5d8b00857b66afb73c7cbbe922b65e87c57e767a..a12b533bbc0ef756737920b750a5c00a9556a58f 100644
--- a/Framework/Utilities/COMBoost.cc
+++ b/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
diff --git a/Framework/Utilities/COMBoost.h b/Framework/Utilities/COMBoost.h
index eebc8358ad5f5891587c3ee9bda036e7c6b98784..57867769dc420c77630ebe3d446f0834999dd988 100644
--- a/Framework/Utilities/COMBoost.h
+++ b/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
diff --git a/Framework/Utilities/testCOMBoost.cc b/Framework/Utilities/testCOMBoost.cc
index 501db86435b535f14b47e6f072cbf885878e8293..159f97a7e598fe269d566fa11f53e36bab3cb24a 100644
--- a/Framework/Utilities/testCOMBoost.cc
+++ b/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}};