From 407c2b3100ecf18d5c23b7017cb827dba207ae78 Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@tu-dortmund.de>
Date: Sun, 24 May 2020 18:41:22 +0200
Subject: [PATCH] improved test of Pythia decay
---
Framework/Geometry/testGeometry.cc | 2 --
Processes/Pythia/testPythia.cc | 19 ++++++++++++++++---
2 files changed, 16 insertions(+), 5 deletions(-)
diff --git a/Framework/Geometry/testGeometry.cc b/Framework/Geometry/testGeometry.cc
index be4dccbc7..6e054d336 100644
--- a/Framework/Geometry/testGeometry.cc
+++ b/Framework/Geometry/testGeometry.cc
@@ -127,7 +127,6 @@ TEST_CASE("transformations between CoordinateSystems") {
SECTION("RotateToZ positive") {
Vector const v{rootCS, 0_m, 1_m, 1_m};
auto const csPrime = rootCS.RotateToZ(v);
- auto const transform = csPrime.GetTransform().matrix();
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};
@@ -155,7 +154,6 @@ TEST_CASE("transformations between CoordinateSystems") {
SECTION("RotateToZ negative") {
Vector const v{rootCS, 0_m, 0_m, -1_m};
auto const csPrime = rootCS.RotateToZ(v);
- auto const transform = csPrime.GetTransform().matrix();
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};
diff --git a/Processes/Pythia/testPythia.cc b/Processes/Pythia/testPythia.cc
index 5dd0cc8bc..daef1f9ef 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") {
--
GitLab