improved test of Pythia decay

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};

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") {