diff git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 328f952036e6a3a119d41628677463001beee75e..62aea277473939304d00f631184f3d76414017bf 100644
 a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ 140,38 +140,90 @@ TEST_CASE("SibyllInterface", "[processes]") {
/*
Interactions between hadrons (h) and nuclei (A) in Sibyll are treated in the
 hadronnucleon centerofmass frame (hnCoM). In addition the incoming hadron (h) and
 nucleon (n) are assumed massless, such that the energy and momentum in the hnCoM are
+ hadronnucleon centerofmass frame (hnCoM). The incoming hadron (h) and
+ nucleon (N) are assumed massless, such that the energy and momentum in the hnCoM are
: E_i_cm = 0.5 * SQS and P_i_cm = + 0.5 * SQS where i is either the projectile
hadron or the target nucleon and SQS is the hadronnucleon centerofmass energy.
+ The true energies and momenta, accounting for the hadron masses, are: E_i = ( S +
+ m_i**2  m_j**2 ) / (2 * SQS) and Pcm = +
+ sqrt( (S(m_j+m_i)**2) * (s(m_jm_i)**2) ) / (2*SQS) where m_i is the projectiles
+ mass and m_j is the target particles mass. In terms of lab. frame variables Pcm =
+ m_j * Plab_i / SQS, where Plab_i is the momentum of the projectile (i) in the lab.
+ and m_j is the mass of the target, i.e. the particle at rest (usually a nucleon).
+
Any hadronnucleus event can contain several nucleon interactions. In case of Nw
 (number of wounded nucleons) in the hadronnucleus interaction the total energy and
 momentum in the hadronnucleon centerofmass frame are: momentum: p_projectile +
 p_nucleon_1 + p_nucleon_2 + .... p_nucleon_Nw = (Nw1) Pcm with centerofmass
 momentum Pcm = p_projectile =  p_nucleon_i energy: E_projectile + E_nucleon_1 + ...
 E_nucleon_Nw = (Nw+1) / 2 * SQS with SQS the hadronnucleon centerofmass energy.
+ (number of wounded nucleons) nucleons interacting in the hadronnucleus interaction,
+ the total energy and momentum in the hadron(i)nucleon(N) centerofmass frame are:
+ momentum: p_projectile + p_nucleon_1 + p_nucleon_2 + .... p_nucleon_Nw = (Nw1) *
+ Pcm with centerofmass momentum Pcm = p_projectile =  p_nucleon_i. For the energy:
+ E_projectile + E_nucleon_1 + ... E_nucleon_Nw = E_projectile + Nw * E_nucleon.
+
+ Using the above definitions of centerofmass energies and momenta this leads to the
+ total energy: E_tot = SQS/2 * (1+Nw) + (m_N**2m_i**2)/(2*SQS) * (Nw1) and P_tot
+ = m_N * Plab_i / SQS * (Nw1).
+
+ A Lorentztransformation of these quantities to the lab. frame recovers Plab_i for
+ the total momentum, so momentum is exactly conserved, and Elab_i + Nw * m_N for the
+ total energy. Not surprisingly the total energy differs from the total energy before
+ the collision by the mass of the additional nucleons (Nw1)*m_N. In relative terms
+ the additional energy is entirely negligible and as it is not kinetic energy there
+ is zero influence on the shower development.
+
+ Due to the ommission of the hadron masses in Sibyll, the total energy and momentum
+ in the centerofmass system after the collision are just: E_tot = SQS/2 * (1+Nw)
+ and P_tot = SQS/2 * (1Nw). After the Lorentztransformation the total momentum in
+ the lab. thus differs from the initial value by (1Nw)/2 * ( m_N + m_i**2 / (2 *
+ Plab_i) ) and momentum is NOT conserved. Note however that the second term quickly
+ vanishes as the lab. momentum of the projectile increases. The first term is fixed
+ as it depends only on the number of additional nucleons, in relative terms it is
+ always small at high energies.
+
+ For this reason the numerical precision in these tests is limited to 5% to still
+ pass at low energies and no absolute check is implemented, e.g.
+
+ CHECK(pSum.GetComponents(cs).GetX() / P0 == Approx(1).margin(0.05));
+ CHECK((pSum  plab).norm()/1_GeV == Approx(0).margin(plab.norm() * 0.05/1_GeV));
+
+ /FR'2020
+
+ See also:
 In case of a single interaction (Nw=1), by definition, the total momentum is zero in
 the hadronnucleon frame. After the boost to the lab. frame the lab. momentum of the
 projectile before the interaction is recoverred.
+ Issue 272 / MR 204
+ https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika//merge_requests/204
 In case of multiple interactions (Nw>1), the momentum is not zero and the total
 momentum in the lab. frame after the boost is different from the original projectile
 (momentum violation).
+ */
 The level of violation of momentum conservation is further enhanced due to the
 approximation of massless hadrons. At low energies (~10GeV), where the masses can
 not be neglected the violation is at the level of percent.
+ CHECK(pSum.GetComponents(cs).GetX() / P0 == Approx(1).margin(0.05));
+ CHECK(pSum.GetComponents(cs).GetY() / 1_GeV == Approx(0).margin(1e4));
+ CHECK(pSum.GetComponents(cs).GetZ() / 1_GeV == Approx(0).margin(1e4));
 For this reason the numerical precision in these tests is limited to 5% at low
 energies. see also:
+ CHECK((pSum  plab).norm() / 1_GeV == Approx(0).margin(plab.norm() * 0.05 / 1_GeV));
+ CHECK(pSum.norm() / P0 == Approx(1).margin(0.05));
+ [[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
+ }
 Issue 272 / MR 204
 https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika//merge_requests/204
+ SECTION("InteractionInterface  high energy") {
 */
+ setup::Stack stack;
+ const HEPEnergyType E0 = 60_EeV;
+ HEPMomentumType P0 =
+ sqrt(E0 * E0  particles::Proton::GetMass() * particles::Proton::GetMass());
+ auto plab = corsika::stack::MomentumVector(cs, {P0, 0_eV, 0_eV});
+ geometry::Point pos(cs, 0_m, 0_m, 0_m);
+ auto particle = stack.AddParticle(
+ std::tuple{
+ particles::Code::Proton, E0, plab, pos, 0_ns});
+ particle.SetNode(nodePtr);
+ corsika::stack::SecondaryView view(particle);
+ auto projectile = view.GetProjectile();
+
+ Interaction model;
+ model.Init();
+ [[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
+ auto const pSum = sumMomentum(view, cs);
CHECK(pSum.GetComponents(cs).GetX() / P0 == Approx(1).margin(0.05));
CHECK(pSum.GetComponents(cs).GetY() / 1_GeV == Approx(0).margin(1e4));
CHECK(pSum.GetComponents(cs).GetZ() / 1_GeV == Approx(0).margin(1e4));
@@ 180,7 +232,7 @@ TEST_CASE("SibyllInterface", "[processes]") {
CHECK(pSum.norm() / P0 == Approx(1).margin(0.05));
[[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
}

+
SECTION("NuclearInteractionInterface") {
setup::Stack stack;