diff --git a/Processes/Sibyll/Interaction.cc b/Processes/Sibyll/Interaction.cc
index 8300462c6de387e9f0ed62d8adca49db96c3fcd3..d738dc81caa233416e74b1ca56d5e6e3d5662473 100644
--- a/Processes/Sibyll/Interaction.cc
+++ b/Processes/Sibyll/Interaction.cc
@@ -312,26 +312,26 @@ namespace corsika::process::sibyll {
auto const pCoM = Vector<hepmomentum_d>(csPrime, tmp);
HEPEnergyType const eCoM = psib.GetEnergy();
auto const Plab = boost.fromCoM(FourVector(eCoM, pCoM));
-
- auto const pid = process::sibyll::ConvertFromSibyll(psib.GetPID());
+ auto const p3lab = Plab.GetSpaceLikeComponents();
+ assert(p3lab.GetCoordinateSystem() == originalCS); // just to be sure!
// add to corsika stack
auto pnew = vP.AddSecondary(
tuple<particles::Code, units::si::HEPEnergyType, stack::MomentumVector,
geometry::Point, units::si::TimeType>{
- pid, Plab.GetTimeLikeComponent(), Plab.GetSpaceLikeComponents(), pOrig,
- tOrig});
+ process::sibyll::ConvertFromSibyll(psib.GetPID()),
+ Plab.GetTimeLikeComponent(), p3lab, pOrig, tOrig});
Plab_final += pnew.GetMomentum();
Elab_final += pnew.GetEnergy();
Ecm_final += psib.GetEnergy();
}
cout << "conservation (all GeV):" << endl
- << "Ecm_initial=" << Ecm / 1_GeV << " Ecm_final=" << Ecm_final / 1_GeV
- << endl
- << "Elab_initial=" << eProjectileLab / 1_GeV
- << " Elab_final=" << Elab_final / 1_GeV
- << " diff (%)=" << (Elab_final / eProjectileLab / get_nwounded() - 1) * 100
+ << "Ecm_initial(per nucleon)=" << Ecm / 1_GeV << " Ecm_final(per nucleon)="
+ << Ecm_final * 2. / (get_nwounded() + 1) / 1_GeV << endl
+ << "Elab_initial=" << Etot / 1_GeV << " Elab_final=" << Elab_final / 1_GeV
+ << " diff (%)=" << (Elab_final / Etot / get_nwounded() - 1) * 100
+ << " E in nucleons=" << constants::nucleonMass * get_nwounded() / 1_GeV
<< endl
<< "Plab_initial=" << (pProjectileLab / 1_GeV).GetComponents()
<< ", Plab_final=" << (Plab_final / 1_GeV).GetComponents() << endl;
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 308e8ea1835d8bea50f57c79d4074024d595ba3f..dee7ee31b722aaabdde09e864c0bbfe80d0b36b3 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -117,7 +117,7 @@ TEST_CASE("SibyllInterface", "[processes]") {
random::RNGManager::GetInstance().RegisterRandomStream("sibyll");
- SECTION("InteractionInterface") {
+ SECTION("InteractionInterface - low energy") {
setup::Stack stack;
const HEPEnergyType E0 = 60_GeV;
@@ -136,8 +136,99 @@ TEST_CASE("SibyllInterface", "[processes]") {
Interaction model;
[[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
- [[maybe_unused]] auto const pSum = sumMomentum(view, cs);
+ auto const pSum = sumMomentum(view, cs);
+
+ /*
+ Interactions between hadrons (h) and nuclei (A) in Sibyll are treated in the
+ hadron-nucleon center-of-mass 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 hadron-nucleon center-of-mass 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_j-m_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 hadron-nucleus event can contain several nucleon interactions. In case of Nw
+ (number of wounded nucleons) nucleons interacting in the hadron-nucleus interaction,
+ the total energy and momentum in the hadron(i)-nucleon(N) center-of-mass frame are:
+ momentum: p_projectile + p_nucleon_1 + p_nucleon_2 + .... p_nucleon_Nw = -(Nw-1) *
+ Pcm with center-of-mass 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 center-of-mass energies and momenta this leads to the
+ total energy: E_tot = SQS/2 * (1+Nw) + (m_N**2-m_i**2)/(2*SQS) * (Nw-1) and P_tot
+ = -m_N * Plab_i / SQS * (Nw-1).
+
+ 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 (Nw-1)*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 center-of-mass system after the collision are just: E_tot = SQS/2 * (1+Nw)
+ and P_tot = SQS/2 * (1-Nw). After the Lorentztransformation the total momentum in
+ the lab. thus differs from the initial value by (1-Nw)/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:
+
+ Issue 272 / MR 204
+ https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/-/merge_requests/204
+
+ */
+
+ CHECK(pSum.GetComponents(cs).GetX() / P0 == Approx(1).margin(0.05));
+ CHECK(pSum.GetComponents(cs).GetY() / 1_GeV == Approx(0).margin(1e-4));
+ CHECK(pSum.GetComponents(cs).GetZ() / 1_GeV == Approx(0).margin(1e-4));
+
+ 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);
+ }
+
+ 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, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ particles::Code::Proton, E0, plab, pos, 0_ns});
+ particle.SetNode(nodePtr);
+ corsika::stack::SecondaryView view(particle);
+ auto projectile = view.GetProjectile();
+
+ Interaction model;
+
+ [[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.001));
+ CHECK(pSum.GetComponents(cs).GetY() / 1_GeV == Approx(0).margin(1e-4));
+ CHECK(pSum.GetComponents(cs).GetZ() / 1_GeV == Approx(0).margin(1e-4));
+ CHECK((pSum - plab).norm() / 1_GeV == Approx(0).margin(plab.norm() * 0.001 / 1_GeV));
+ CHECK(pSum.norm() / P0 == Approx(1).margin(0.05));
[[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
}