From d00feece03d51fa047d114ef644e63266f040417 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Thu, 23 Mar 2023 12:02:06 +0000
Subject: [PATCH] added test case to document the interface inconsistency
---
tests/modules/testSophia.cpp | 36 ++++++++++++++++++++++++++++++++++++
1 file changed, 36 insertions(+)
diff --git a/tests/modules/testSophia.cpp b/tests/modules/testSophia.cpp
index 8043ddbcb..521b315eb 100644
--- a/tests/modules/testSophia.cpp
+++ b/tests/modules/testSophia.cpp
@@ -28,6 +28,7 @@
link them togehter, it will fail.
*/
#include <corsika/modules/sophia/Random.hpp>
+#include "corsika/framework/core/PhysicalConstants.hpp"
using namespace corsika;
using namespace corsika::sophia;
@@ -121,6 +122,41 @@ TEST_CASE("SophiaInterface", "modules") {
CHECK_THROWS(model.doInteraction(view, Code::Electron, Code::Proton, aP4, bP4));
}
+ SECTION("InteractionInterface - sqrtSNN inconsistency") {
+ corsika::sophia::InteractionModel model;
+ // the pion-production threshold is given by the CoM energy of the neutron-pi+ system.
+ // sth = m_n**2 + 2 m_n * m_pi+ + m_pi+**2 = 1.1646 GeV**2
+ HEPEnergyType const eGamma =
+ 0.15145_GeV; // energy of photon at pion-production threshold w proton target
+ FourMomentum const photonP4(eGamma, {cs, eGamma, 0_GeV, 0_GeV});
+
+ // nucleon case (this was implemented in the HadronicPhotonModel)
+ // it passes isValid and then rejects event generation when the proton is selected as
+ // target nucleon
+ FourMomentum const nucleonP4(constants::nucleonMass, {cs, 0_GeV, 0_GeV, 0_GeV});
+ auto const sqrtSNN_nuc = (photonP4 + nucleonP4).getNorm();
+ CHECK(model.isValid(Code::Photon, Code::Neutron, sqrtSNN_nuc));
+ view.clear();
+ // this will not throw but SOPHIA will abort inernally. catch by checking stack size
+ model.doInteraction(view, Code::Photon, Code::Proton, photonP4, nucleonP4);
+ CHECK_FALSE(view.getSize() > 0);
+
+ // proton case (will throw in interface)
+ FourMomentum const protonP4(Proton::mass, {cs, 0_GeV, 0_GeV, 0_GeV});
+ auto const sqrtSNN_p = (photonP4 + protonP4).getNorm();
+ CHECK_FALSE(model.isValid(Code::Photon, Code::Proton, sqrtSNN_p));
+ view.clear();
+ CHECK_THROWS(
+ model.doInteraction(view, Code::Photon, Code::Proton, photonP4, protonP4));
+
+ // neutron case
+ FourMomentum const neutronP4(Neutron::mass, {cs, 0_GeV, 0_GeV, 0_GeV});
+ auto const sqrtSNN_n = (photonP4 + neutronP4).getNorm();
+ CHECK(model.isValid(Code::Photon, Code::Neutron, sqrtSNN_n));
+ view.clear();
+ model.doInteraction(view, Code::Photon, Code::Neutron, photonP4, neutronP4);
+ }
+
SECTION("InteractionInterface - interaction") {
const HEPEnergyType P0 = 1.2_GeV;
MomentumVector const plab = MomentumVector(cs, {P0, 0_eV, 0_eV});
--
GitLab