From d748d69a6956429b594ab5698bb01bbe352f8245 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Fri, 31 Jul 2020 11:03:47 +0100
Subject: [PATCH] comment on limited precision
---
Processes/Sibyll/testSibyll.cc | 49 ++++++++++++++++++++++++++--------
1 file changed, 38 insertions(+), 11 deletions(-)
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 94aa6394c..23c9eb043 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -119,7 +119,7 @@ TEST_CASE("SibyllInterface", "[processes]") {
random::RNGManager::GetInstance().RegisterRandomStream("s_rndm");
- SECTION("InteractionInterface") {
+ SECTION("InteractionInterface - low energy") {
setup::Stack stack;
const HEPEnergyType E0 = 60_GeV;
@@ -141,22 +141,49 @@ TEST_CASE("SibyllInterface", "[processes]") {
[[maybe_unused]] const process::EProcessReturn ret = model.DoInteraction(projectile);
auto const pSum = sumMomentum(view, cs);
- // for debugging!
- std::cout << pSum.GetComponents(cs) << std::endl;
- std::cout << plab.GetComponents(cs) << std::endl;
+ /*
+ Interactions between hadrons (h) and nuclei (A) in Sibyll are treated in the
+ hadron-nucleon center-of-mass frame (hnCoM). In addition 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.
- CHECK(pSum.GetComponents(cs).GetX() / P0 ==
- Approx(1).margin(model.get_relative_precision_momentum()));
+ Any hadron-nucleus event can contain several nucleon interactions. In case of Nw
+ (number of wounded nucleons) in the hadron-nucleus interaction the total energy and
+ momentum in the hadron-nucleon 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 energy: E_projectile + E_nucleon_1 + ...
+ E_nucleon_Nw = (Nw+1) / 2 * SQS with SQS the hadron-nucleon center-of-mass energy.
+
+ In case of a single interaction (Nw=1), by definition, the total momentum is zero in
+ the hadron-nucleon frame. After the boost to the lab. frame the lab. momentum of the
+ projectile before the interaction is recoverred.
+
+ 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.
+
+ For this reason the numerical precision in these tests is limited to 5% at low
+ energies. 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() * model.get_relative_precision_momentum() / 1_GeV));
- CHECK(pSum.norm() / P0 == Approx(1).margin(model.get_relative_precision_momentum()));
+ 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("NuclearInteractionInterface") {
setup::Stack stack;
--
GitLab