diff --git a/Framework/Geometry/CoordinateSystem.h b/Framework/Geometry/CoordinateSystem.h
index 26da1a701139a5ed1b887d516582888807ae9d21..614c76da9d2f32925fa9d8c18828196fbbe8f19f 100644
--- a/Framework/Geometry/CoordinateSystem.h
+++ b/Framework/Geometry/CoordinateSystem.h
@@ -72,19 +72,19 @@ namespace corsika::geometry {
Eigen::Matrix3d A, B;
if (s > 0) {
- A << 1, 0, -a1, // comment to prevent clang-format
- 0, 1, -a2, // .
- a1, a2, 1; // .
+ 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, // .
- 0, -1, -a2, // .
- a1, a2, -1; // .
- B << -a1 * a1 * c, -a1 * a2 * c, 0, // .
- +a1 * a2 * c, +a2 * a2 * c, 0, // .
+ 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; // .
}
diff --git a/Framework/Geometry/testGeometry.cc b/Framework/Geometry/testGeometry.cc
index 440e771e73446de8729d6c57a9b8aee69ca80414..6e054d336b1a2a66e6e615b506079eae4f4aebbf 100644
--- a/Framework/Geometry/testGeometry.cc
+++ b/Framework/Geometry/testGeometry.cc
@@ -123,6 +123,61 @@ TEST_CASE("transformations between CoordinateSystems") {
auto comp3 = v1.GetComponents(combined);
REQUIRE((comp1 - comp3).norm().magnitude() == Approx(0).margin(absMargin));
}
+
+ SECTION("RotateToZ positive") {
+ Vector const v{rootCS, 0_m, 1_m, 1_m};
+ auto const csPrime = rootCS.RotateToZ(v);
+ Vector const zPrime{csPrime, 0_m, 0_m, 5_m};
+ Vector const xPrime{csPrime, 5_m, 0_m, 0_m};
+ Vector const yPrime{csPrime, 0_m, 5_m, 0_m};
+
+ CHECK(zPrime.dot(v).magnitude() > 0);
+ CHECK(zPrime.GetComponents(rootCS)[1].magnitude() ==
+ Approx(zPrime.GetComponents(rootCS)[2].magnitude()));
+ CHECK(zPrime.GetComponents(rootCS)[0].magnitude() == Approx(0));
+
+ CHECK(xPrime.GetComponents(rootCS).eVector.dot(
+ yPrime.GetComponents(rootCS).eVector) == Approx(0));
+ CHECK(zPrime.GetComponents(rootCS).eVector.dot(
+ xPrime.GetComponents(rootCS).eVector) == Approx(0));
+ CHECK(yPrime.GetComponents(rootCS).eVector.dot(
+ zPrime.GetComponents(rootCS).eVector) == Approx(0));
+
+ CHECK(yPrime.GetComponents(rootCS).eVector.dot(
+ yPrime.GetComponents(rootCS).eVector) == Approx((5_m * 5_m).magnitude()));
+ CHECK(xPrime.GetComponents(rootCS).eVector.dot(
+ xPrime.GetComponents(rootCS).eVector) == Approx((5_m * 5_m).magnitude()));
+ CHECK(zPrime.GetComponents(rootCS).eVector.dot(
+ zPrime.GetComponents(rootCS).eVector) == Approx((5_m * 5_m).magnitude()));
+ }
+
+ SECTION("RotateToZ negative") {
+ Vector const v{rootCS, 0_m, 0_m, -1_m};
+ auto const csPrime = rootCS.RotateToZ(v);
+ Vector const zPrime{csPrime, 0_m, 0_m, 5_m};
+ Vector const xPrime{csPrime, 5_m, 0_m, 0_m};
+ Vector const yPrime{csPrime, 0_m, 5_m, 0_m};
+
+ CHECK(zPrime.dot(v).magnitude() > 0);
+ CHECK(xPrime.GetComponents(rootCS).eVector.dot(v.GetComponents().eVector) ==
+ Approx(0));
+ CHECK(yPrime.GetComponents(rootCS).eVector.dot(v.GetComponents().eVector) ==
+ Approx(0));
+
+ CHECK(xPrime.GetComponents(rootCS).eVector.dot(
+ yPrime.GetComponents(rootCS).eVector) == Approx(0));
+ CHECK(zPrime.GetComponents(rootCS).eVector.dot(
+ xPrime.GetComponents(rootCS).eVector) == Approx(0));
+ CHECK(yPrime.GetComponents(rootCS).eVector.dot(
+ zPrime.GetComponents(rootCS).eVector) == Approx(0));
+
+ CHECK(yPrime.GetComponents(rootCS).eVector.dot(
+ yPrime.GetComponents(rootCS).eVector) == Approx((5_m * 5_m).magnitude()));
+ CHECK(xPrime.GetComponents(rootCS).eVector.dot(
+ xPrime.GetComponents(rootCS).eVector) == Approx((5_m * 5_m).magnitude()));
+ CHECK(zPrime.GetComponents(rootCS).eVector.dot(
+ zPrime.GetComponents(rootCS).eVector) == Approx((5_m * 5_m).magnitude()));
+ }
}
TEST_CASE("Sphere") {
diff --git a/Framework/Utilities/COMBoost.cc b/Framework/Utilities/COMBoost.cc
index 1d65d6439a2731a70988d440b76ed4e68a77262f..abdb7a48b88f733aba4c03dc4d8521a741aa27b5 100644
--- a/Framework/Utilities/COMBoost.cc
+++ b/Framework/Utilities/COMBoost.cc
@@ -15,58 +15,50 @@
#include <corsika/utl/COMBoost.h>
#include <corsika/utl/sgn.h>
+#include <cmath>
+
using namespace corsika::utl;
using namespace corsika::units::si;
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 pProjNorm = pProjectile.norm();
- auto const a = (pProjectile / pProjNorm).GetComponents().eVector;
- auto const a1 = a(0), a2 = a(1);
-
- auto const s = sgn(a(2));
- auto const c = 1 / (1 + s * a(2));
-
- Eigen::Matrix3d A, B;
-
- 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; // .
-
- } 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; // .
- }
+ auto const pProjNormSquared = pProjectile.squaredNorm();
+ auto const pProjNorm = sqrt(pProjNormSquared);
- fRotation = A + B;
+ auto const eProjectile = Pprojectile.GetTimeLikeComponent();
+ auto const massProjectileSquared = eProjectile * eProjectile - pProjNormSquared;
+ auto const s =
+ massTarget * massTarget + massProjectileSquared + 2 * eProjectile * massTarget;
- // calculate boost
- double const beta = pProjNorm / (Pprojectile.GetTimeLikeComponent() + massTarget);
+ auto const sqrtS = sqrt(s);
+ auto const sinhEta = -pProjNorm / sqrtS;
+ auto const coshEta = sqrt(1 + pProjNormSquared / s);
- /* Accurracy matters here, beta = 1 - epsilon for ultra-relativistic boosts */
- double const coshEta = 1 / std::sqrt((1 + beta) * (1 - beta));
- //~ double const coshEta = 1 / std::sqrt((1-beta*beta));
- double const sinhEta = -beta * coshEta;
+ setBoost(coshEta, sinhEta);
- std::cout << "COMBoost (1-beta)=" << 1 - beta << " gamma=" << coshEta << std::endl;
- std::cout << " det = " << fRotation.determinant() - 1 << std::endl;
+ std::cout << "COMBoost (1-beta)=" << 1 - sinhEta / coshEta << " gamma=" << coshEta
+ << 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()));
+ }
}
diff --git a/Processes/Pythia/Decay.cc b/Processes/Pythia/Decay.cc
index 7586f0b42c63e88f532ba6bba17f27348a09ae9c..27890d28bfe1ce42576fc81d08550ea80daf70de 100644
--- a/Processes/Pythia/Decay.cc
+++ b/Processes/Pythia/Decay.cc
@@ -12,8 +12,10 @@
#include <corsika/process/pythia/Decay.h>
#include <corsika/process/pythia/Random.h>
+#include <corsika/geometry/FourVector.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/utl/COMBoost.h>
using std::cout;
using std::endl;
@@ -103,16 +105,14 @@ namespace corsika::process::pythia {
void Decay::SetHandleDecay(const vector<particles::Code> vParticleList) {
handleAllDecays_ = false;
- for (auto p : vParticleList) Decay::SetHandleDecay(p);
+ for (auto p : vParticleList) SetHandleDecay(p);
}
bool Decay::IsDecayHandled(const corsika::particles::Code vParticleCode) {
- if (handleAllDecays_ && Decay::CanHandleDecay(vParticleCode))
+ if (handleAllDecays_ && CanHandleDecay(vParticleCode))
return true;
else
- return Decay::handledDecays_.find(vParticleCode) != Decay::handledDecays_.end()
- ? true
- : false;
+ return handledDecays_.find(vParticleCode) != Decay::handledDecays_.end();
}
bool Decay::IsStable(const particles::Code vCode) {
@@ -192,12 +192,17 @@ namespace corsika::process::pythia {
using namespace units;
using namespace units::si;
- auto const decayPoint = vP.GetPosition();
+ auto const& decayPoint = vP.GetPosition();
auto const t0 = vP.GetTime();
- // coordinate system, get global frame of reference
- geometry::CoordinateSystem& rootCS =
- geometry::RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+ auto const& labMomentum = vP.GetMomentum();
+ geometry::CoordinateSystem const& labCS = labMomentum.GetCoordinateSystem();
+
+ // define target kinematics in lab frame
+ // define boost to and from CoM frame
+ // CoM frame definition in Pythia projectile: +z
+ utl::COMBoost const boost(labMomentum, vP.GetMass());
+ auto const& rotatedCS = boost.GetRotatedCS();
fCount++;
@@ -213,12 +218,11 @@ namespace corsika::process::pythia {
// input particle PDG
auto const pdgCode = static_cast<int>(particles::GetPDG(particleId));
- auto const pcomp = vP.GetMomentum().GetComponents();
- double px = pcomp[0] / 1_GeV;
- double py = pcomp[1] / 1_GeV;
- double pz = pcomp[2] / 1_GeV;
- double en = vP.GetEnergy() / 1_GeV;
- double m = particles::GetMass(particleId) / 1_GeV;
+ double constexpr px = 0;
+ double constexpr py = 0;
+ double constexpr pz = 0;
+ double const en = vP.GetMass() / 1_GeV;
+ double const m = en;
// add particle to pythia stack
event.append(pdgCode, 1, 0, 0, px, py, pz, en, m);
@@ -236,17 +240,22 @@ namespace corsika::process::pythia {
if (event[i].isFinal()) {
auto const pyId =
particles::ConvertFromPDG(static_cast<particles::PDGCode>(event[i].id()));
- HEPEnergyType pyEn = event[i].e() * 1_GeV;
- MomentumVector pyP(rootCS, {event[i].px() * 1_GeV, event[i].py() * 1_GeV,
- event[i].pz() * 1_GeV});
+ HEPEnergyType const Erest = event[i].e() * 1_GeV;
+ MomentumVector const pRest(
+ rotatedCS,
+ {event[i].px() * 1_GeV, event[i].py() * 1_GeV, event[i].pz() * 1_GeV});
+ geometry::FourVector const fourMomRest{Erest, pRest};
+ auto const fourMomLab = boost.fromCoM(fourMomRest);
- cout << "particle: id=" << pyId << " momentum=" << pyP.GetComponents() / 1_GeV
- << " energy=" << pyEn << endl;
+ cout << "particle: id=" << pyId << " momentum="
+ << fourMomLab.GetSpaceLikeComponents().GetComponents(labCS) / 1_GeV
+ << " energy=" << fourMomLab.GetTimeLikeComponent() << endl;
vP.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType,
corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
- pyId, pyEn, pyP, decayPoint, t0});
+ pyId, fourMomLab.GetTimeLikeComponent(),
+ fourMomLab.GetSpaceLikeComponents(), decayPoint, t0});
}
// set particle stable
diff --git a/Processes/Pythia/testPythia.cc b/Processes/Pythia/testPythia.cc
index 5dd0cc8bc6c60e23697d722c8d983ce0f44efc89..daef1f9ef1648a65a184a0b075296bac3b7e16df 100644
--- a/Processes/Pythia/testPythia.cc
+++ b/Processes/Pythia/testPythia.cc
@@ -19,6 +19,7 @@
#include <corsika/geometry/Point.h>
#include <corsika/units/PhysicalUnits.h>
+#include <corsika/utl/CorsikaFenv.h>
#include <catch2/catch.hpp>
TEST_CASE("Pythia", "[processes]") {
@@ -88,6 +89,15 @@ TEST_CASE("Pythia", "[processes]") {
using namespace corsika;
using namespace corsika::units::si;
+template <typename TStackView>
+auto sumMomentum(TStackView const& view, geometry::CoordinateSystem const& vCS) {
+ geometry::Vector<hepenergy_d> sum{vCS, 0_eV, 0_eV, 0_eV};
+
+ for (auto const& p : view) { sum += p.GetMomentum(); }
+
+ return sum;
+}
+
TEST_CASE("pythia process") {
// setup environment, geometry
@@ -110,12 +120,12 @@ TEST_CASE("pythia process") {
auto const* nodePtr = theMedium.get(); // save the medium for later use before moving it
SECTION("pythia decay") {
-
+ feenableexcept(FE_INVALID);
setup::Stack stack;
const HEPEnergyType E0 = 10_GeV;
HEPMomentumType P0 =
sqrt(E0 * E0 - particles::PiPlus::GetMass() * particles::PiPlus::GetMass());
- auto plab = corsika::stack::MomentumVector(cs, {0_GeV, 0_GeV, -P0});
+ auto plab = corsika::stack::MomentumVector(cs, {0_GeV, 0_GeV, P0});
geometry::Point pos(cs, 0_m, 0_m, 0_m);
auto particle = stack.AddParticle(
std::tuple<particles::Code, units::si::HEPEnergyType,
@@ -131,7 +141,10 @@ TEST_CASE("pythia process") {
model.Init();
[[maybe_unused]] const TimeType time = model.GetLifetime(particle);
model.DoDecay(projectile);
- REQUIRE(stack.GetSize() == 3);
+ CHECK(stack.GetSize() == 3);
+ auto const pSum = sumMomentum(view, cs);
+ CHECK((pSum - plab).norm() / 1_GeV == Approx(0).margin(1e-4));
+ CHECK((pSum.norm() - plab.norm()) / 1_GeV == Approx(0).margin(1e-4));
}
SECTION("pythia decay config") {