added test case to document the interface inconsistency

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