From 24ff1d8d2ae52fa8f8d315b6f4ac447ee020f4b8 Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Mon, 1 Nov 2021 15:25:04 +0100
Subject: [PATCH] catch special case in COMBoost
---
corsika/detail/framework/utility/COMBoost.inl | 27 +++++----
corsika/detail/modules/epos/Interaction.inl | 10 ++--
tests/framework/testCOMBoost.cpp | 55 +++++++++++++++++--
3 files changed, 71 insertions(+), 21 deletions(-)
diff --git a/corsika/detail/framework/utility/COMBoost.inl b/corsika/detail/framework/utility/COMBoost.inl
index 7551b6e0a..4905b940b 100644
--- a/corsika/detail/framework/utility/COMBoost.inl
+++ b/corsika/detail/framework/utility/COMBoost.inl
@@ -22,9 +22,7 @@ namespace corsika {
inline COMBoost::COMBoost(FourMomentum const& P4projectile,
HEPMassType const massTarget)
: originalCS_{P4projectile.getSpaceLikeComponents().getCoordinateSystem()}
- , rotatedCS_{
- make_rotationToZ(P4projectile.getSpaceLikeComponents().getCoordinateSystem(),
- P4projectile.getSpaceLikeComponents())} {
+ , rotatedCS_{make_rotationToZ(originalCS_, P4projectile.getSpaceLikeComponents())} {
auto const pProjectile = P4projectile.getSpaceLikeComponents();
auto const pProjNormSquared = pProjectile.getSquaredNorm();
auto const pProjNorm = sqrt(pProjNormSquared);
@@ -45,15 +43,20 @@ namespace corsika {
inline COMBoost::COMBoost(FourMomentum const& P4projectile,
FourMomentum const& P4target)
- : originalCS_{P4projectile.getSpaceLikeComponents().getCoordinateSystem()}
- , rotatedCS_{make_rotationToZ(
- P4projectile.getSpaceLikeComponents().getCoordinateSystem(),
- P4projectile.getSpaceLikeComponents() + P4target.getSpaceLikeComponents())} {
+ : originalCS_{P4projectile.getSpaceLikeComponents().getCoordinateSystem()} {
+
// this is the center-of-momentum CM frame
auto const pCM =
P4projectile.getSpaceLikeComponents() + P4target.getSpaceLikeComponents();
auto const pCM2 = pCM.getSquaredNorm();
auto const pCMnorm = sqrt(pCM2);
+ if (pCMnorm == 0_eV) {
+ // CM is at reset
+ rotatedCS_ = originalCS_;
+ } else {
+ rotatedCS_ = make_rotationToZ(originalCS_, P4projectile.getSpaceLikeComponents() +
+ P4target.getSpaceLikeComponents());
+ }
auto const s = (P4projectile + P4target).getNormSqr();
auto const sqrtS = sqrt(s);
@@ -91,6 +94,9 @@ namespace corsika {
eVecRotated(2) = boostedZ(1) * (1_GeV).magnitude();
+ CORSIKA_LOG_TRACE("E0={}, p={}, E0'={}, p'={}", p4.getTimeLikeComponent() / 1_GeV,
+ eVecRotated(2) * (1 / 1_GeV).magnitude(), E_CoM / 1_GeV, boostedZ);
+
return FourVector(E_CoM, MomentumVector(rotatedCS_, eVecRotated));
}
@@ -102,10 +108,9 @@ namespace corsika {
Eigen::Vector2d com;
com << (Ecm * (1 / 1_GeV)), (pCM.getEigenVector()(2) * (1 / 1_GeV).magnitude());
- CORSIKA_LOG_TRACE(
- "COMBoost::fromCoM Ecm={} GeV"
- " pcm={} GeV (norm = {} GeV), invariant mass={} GeV",
- Ecm / 1_GeV, pCM / 1_GeV, pCM.getNorm() / 1_GeV, p4.getNorm() / 1_GeV);
+ CORSIKA_LOG_TRACE("Ecm={} GeV, pcm={} GeV (norm = {} GeV), invariant mass={} GeV",
+ Ecm / 1_GeV, pCM / 1_GeV, pCM.getNorm() / 1_GeV,
+ p4.getNorm() / 1_GeV);
auto const boostedZ = inverseBoost_ * com;
auto const E_lab = boostedZ(0) * 1_GeV;
diff --git a/corsika/detail/modules/epos/Interaction.inl b/corsika/detail/modules/epos/Interaction.inl
index 33ece71dc..7346b1846 100644
--- a/corsika/detail/modules/epos/Interaction.inl
+++ b/corsika/detail/modules/epos/Interaction.inl
@@ -406,14 +406,14 @@ namespace corsika::epos {
auto const sqrtS2 = (projectileP4 + targetP4).getNormSqr();
auto const sqrtS = sqrt(sqrtS2);
if (!isValid(projectileId, targetId, sqrtS)) {
- throw std::runtime_error("invalid projectiel/target/energy combination.");
+ throw std::runtime_error("invalid projectile/target/energy combination.");
}
HEPEnergyType const Elab = (sqrtS2 - static_pow<2>(get_mass(projectileId)) -
static_pow<2>(get_mass(targetId))) /
(2 * get_mass(targetId));
// system of initial-state
- COMBoost boost(projectileP4, targetP4);
+ COMBoost const boost(projectileP4, targetP4);
auto const& originalCS = boost.getOriginalCS();
auto const& csPrime =
@@ -424,7 +424,7 @@ namespace corsika::epos {
MomentumVector pLab(csPrime, {0_eV, 0_eV, pLabMag});
// internal EPOS lab system
- COMBoost boostInternal({Elab, pLab}, get_mass(targetId));
+ COMBoost const boostInternal({Elab, pLab}, get_mass(targetId));
CORSIKA_LOGGER_DEBUG(logger_, "doInteraction: {} interaction, Elab={} ", projectileId,
Elab);
@@ -467,8 +467,8 @@ namespace corsika::epos {
HEPEnergyType Elab_final = 0_GeV;
// position and time of interaction, not used in QgsjetII
- auto const projectile = view.getProjectile();
- Point const pOrig = projectile.getPosition();
+ auto const& projectile = view.getProjectile();
+ Point const& pOrig = projectile.getPosition();
TimeType const tOrig = projectile.getTime();
// secondaries
diff --git a/tests/framework/testCOMBoost.cpp b/tests/framework/testCOMBoost.cpp
index d092a4b61..f9542e394 100644
--- a/tests/framework/testCOMBoost.cpp
+++ b/tests/framework/testCOMBoost.cpp
@@ -23,17 +23,19 @@ CoordinateSystemPtr rootCS = get_root_CoordinateSystem();
/**
* \todo such helper functions should be moved to the FourVector class:
- **/
+ */
// helper function for energy-momentum
// relativistic energy
-auto const energy = [](HEPMassType m, MomentumVector const& p) {
+auto const energy = [](HEPMassType const m, MomentumVector const& p) {
return sqrt(m * m + p.getSquaredNorm());
};
-auto const momentum = [](HEPEnergyType E, HEPMassType m) { return sqrt(E * E - m * m); };
+auto const momentum = [](HEPEnergyType const E, HEPMassType const m) {
+ return sqrt(E * E - m * m);
+};
// helper function for mandelstam-s
-auto const s = [](HEPEnergyType E, QuantityVector<hepmomentum_d> const& p) {
+auto const s = [](HEPEnergyType const E, QuantityVector<hepmomentum_d> const& p) {
return E * E - p.getSquaredNorm();
};
@@ -170,7 +172,7 @@ TEST_CASE("rotation") {
TEST_CASE("boosts") {
- logging::set_level(logging::level::info);
+ logging::set_level(logging::level::trace);
// define target kinematics in lab frame
HEPMassType const targetMass = 1_GeV;
@@ -322,6 +324,49 @@ TEST_CASE("boosts") {
PprojCoM.getSpaceLikeComponents() + PtargCoM.getSpaceLikeComponents();
CHECK(sumPCoM.getNorm() / P0 == Approx(0).margin(absMargin)); // MAKE RELATIVE CHECK
}
+
+ SECTION("CoM system") {
+
+ MomentumVector pCM{rootCS, 0_GeV, 0_GeV, 5_GeV};
+
+ COMBoost boostCMS({energy(1_GeV, pCM), pCM}, {energy(1_GeV, pCM), -pCM});
+
+ auto test1 = boostCMS.fromCoM(FourMomentum{
+ 0_GeV, MomentumVector(boostCMS.getOriginalCS(), {0_GeV, 0_GeV, 0_GeV})});
+ CHECK(test1.getNorm() == 0_GeV);
+ auto test2 = boostCMS.fromCoM(FourMomentum{
+ 0_GeV, MomentumVector(boostCMS.getRotatedCS(), {0_GeV, 0_GeV, 0_GeV})});
+ CHECK(test2.getNorm() == 0_GeV);
+
+ auto test3 = boostCMS.toCoM(FourMomentum{
+ 0_GeV, MomentumVector(boostCMS.getOriginalCS(), {0_GeV, 0_GeV, 0_GeV})});
+ CHECK(test3.getNorm() == 0_GeV);
+ auto test4 = boostCMS.toCoM(FourMomentum{
+ 0_GeV, MomentumVector(boostCMS.getRotatedCS(), {0_GeV, 0_GeV, 0_GeV})});
+ CHECK(test4.getNorm() == 0_GeV);
+
+ HEPEnergyType const sqrtS =
+ (FourMomentum{energy(1_GeV, pCM), pCM} + FourMomentum{energy(1_GeV, pCM), -pCM})
+ .getNorm();
+ HEPEnergyType const eLab =
+ (static_pow<2>(sqrtS) - 2 * static_pow<2>(1_GeV)) / (2 * 1_GeV);
+ COMBoost boostLab({eLab, MomentumVector{rootCS, momentum(eLab, 1_GeV), 0_eV, 0_eV}},
+ {1_GeV, MomentumVector{rootCS, 0_eV, 0_eV, 0_eV}});
+
+ FourMomentum p4lab_trans(
+ 10_GeV,
+ MomentumVector(boostLab.getOriginalCS(), {0_eV, momentum(10_GeV, 1_GeV), 0_eV}));
+ FourMomentum p4lab_long(
+ 10_GeV,
+ MomentumVector(boostLab.getOriginalCS(), {momentum(10_GeV, 1_GeV), 0_GeV, 0_eV}));
+ // boost of transverse momentum
+ CHECK(boostLab.toCoM(p4lab_trans).getNorm() / 1_GeV == Approx(1));
+ CHECK(boostLab.toCoM(p4lab_trans).getTimeLikeComponent() / 1_GeV == Approx(50.99));
+ // boost of longitudinal momentum
+ CHECK(boostLab.toCoM(p4lab_long).getNorm() / 1_GeV == Approx(1));
+ CHECK(boostLab.toCoM(p4lab_long).getTimeLikeComponent() / 1_GeV ==
+ Approx(1.24).margin(0.1));
+ }
}
TEST_CASE("rest frame") {
--
GitLab