From 2647a2a9b0acf518911729439f80f06ea23908ab Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Mon, 29 Apr 2019 13:14:31 -0300
Subject: [PATCH] test momentum conservation
---
Processes/UrQMD/testUrQMD.cc | 28 +++++++++++++++++++++++++---
1 file changed, 25 insertions(+), 3 deletions(-)
diff --git a/Processes/UrQMD/testUrQMD.cc b/Processes/UrQMD/testUrQMD.cc
index 0a9af8749..3e347a827 100644
--- a/Processes/UrQMD/testUrQMD.cc
+++ b/Processes/UrQMD/testUrQMD.cc
@@ -48,6 +48,15 @@ auto sumCharge(TStackView const& view) {
return totalCharge;
}
+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;
+}
+
auto setupEnvironment(particles::Code vTargetCode) {
// setup environment, geometry
auto env = std::make_unique<environment::Environment<environment::IMediumModel>>();
@@ -159,28 +168,41 @@ TEST_CASE("UrQMD") {
SECTION("nucleon projectile") {
auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Oxygen);
- unsigned short constexpr A = 4, Z = 2;
+ unsigned short constexpr A = 14, Z = 7;
auto [stackPtr, secViewPtr] = setupStack(A, Z, 400_GeV, nodePtr, *csPtr);
// must be assigned to variable, cannot be used as rvalue?!
auto projectile = secViewPtr->GetProjectile();
+ auto const projectileMomentum = projectile.GetMomentum();
[[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(projectile);
REQUIRE(sumCharge(*secViewPtr) ==
Z + particles::GetChargeNumber(particles::Code::Oxygen));
+
+ auto const secMomSum =
+ sumMomentum(*secViewPtr, projectileMomentum.GetCoordinateSystem());
+ REQUIRE((secMomSum - projectileMomentum).norm() / projectileMomentum.norm() ==
+ Approx(0).margin(1e-2));
}
SECTION("\"special\" projectile") {
auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Oxygen);
auto [stackPtr, secViewPtr] =
- setupStack(particles::Code::K0, 400_GeV, nodePtr, *csPtr);
+ setupStack(particles::Code::PiPlus, 400_GeV, nodePtr, *csPtr);
// must be assigned to variable, cannot be used as rvalue?!
auto projectile = secViewPtr->GetProjectile();
+ auto const projectileMomentum = projectile.GetMomentum();
+
[[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(projectile);
REQUIRE(sumCharge(*secViewPtr) ==
- particles::GetChargeNumber(particles::Code::K0) +
+ particles::GetChargeNumber(particles::Code::PiPlus) +
particles::GetChargeNumber(particles::Code::Oxygen));
+
+ auto const secMomSum =
+ sumMomentum(*secViewPtr, projectileMomentum.GetCoordinateSystem());
+ REQUIRE((secMomSum - projectileMomentum).norm() / projectileMomentum.norm() ==
+ Approx(0).margin(1e-2));
}
}
--
GitLab