From 2aed5bcf292161afb88b7a682d797e2e39bc2c40 Mon Sep 17 00:00:00 2001
From: Felix Riehn <felix@matilda>
Date: Mon, 27 Jul 2020 18:05:11 +0100
Subject: [PATCH] added expected precision to sibyll
---
Processes/Sibyll/Interaction.h | 26 ++++++++++++++++++++++++++
Processes/Sibyll/testSibyll.cc | 9 ++++++---
2 files changed, 32 insertions(+), 3 deletions(-)
diff --git a/Processes/Sibyll/Interaction.h b/Processes/Sibyll/Interaction.h
index fcb3b94e5..fd867204b 100644
--- a/Processes/Sibyll/Interaction.h
+++ b/Processes/Sibyll/Interaction.h
@@ -10,6 +10,7 @@
#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/InteractionProcess.h>
+#include <corsika/process/sibyll/sibyll2.3d.h>
#include <corsika/random/RNGManager.h>
#include <corsika/units/PhysicalUnits.h>
#include <tuple>
@@ -35,6 +36,18 @@ namespace corsika::process::sibyll {
int GetMaxTargetMassNumber() const { return maxTargetMassNumber_; }
corsika::units::si::HEPEnergyType GetMinEnergyCoM() const { return minEnergyCoM_; }
corsika::units::si::HEPEnergyType GetMaxEnergyCoM() const { return maxEnergyCoM_; }
+ double get_relative_precision_momentum() const {
+ if (get_nwounded() == 1)
+ return precision_momentum_single_;
+ else
+ return precision_momentum_ * get_nwounded();
+ }
+ double get_relative_precision_energy() const {
+ if (get_nwounded() == 1)
+ return precision_energy_single_;
+ else
+ return precision_energy_ * get_nwounded();
+ }
bool IsValidTarget(corsika::particles::Code TargetId) const {
return (corsika::particles::GetNucleusA(TargetId) < maxTargetMassNumber_) &&
corsika::particles::IsNucleus(TargetId);
@@ -64,6 +77,19 @@ namespace corsika::process::sibyll {
const corsika::units::si::HEPEnergyType maxEnergyCoM_ =
1.e6 * 1e9 * corsika::units::si::electronvolt;
const int maxTargetMassNumber_ = 18;
+ /*
+ for interactions with a single target nucleon energy and momentum conservation
+ are fullfilled
+ for more than one target nucleon, conservation is only approximately true in the
+ lab frame in addition there seems to be a bug in sibyll that leads to violation of
+ momentum conservation already in the nuc-nuc frame. see Issue #272
+
+ https://gitlab.ikp.kit.edu/AirShowerPhysics/corsika/-/issues/272
+ */
+ const double precision_energy_single_ = 1.e-8;
+ const double precision_momentum_single_ = 1.e-8;
+ const double precision_energy_ = 2.e-2;
+ const double precision_momentum_ = 5.e-2;
};
} // namespace corsika::process::sibyll
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index e4ec476d7..bfc3456a1 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -142,12 +142,15 @@ TEST_CASE("SibyllInterface", "[processes]") {
std::cout << pSum.GetComponents(cs) << std::endl;
std::cout << plab.GetComponents(cs) << std::endl;
- CHECK(pSum.GetComponents(cs).GetX() / P0 == Approx(1).margin(1e-4));
+ CHECK(pSum.GetComponents(cs).GetX() / P0 ==
+ Approx(1).margin(model.get_relative_precision_momentum()));
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(1e-4));
- CHECK(pSum.norm() / P0 == Approx(1).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()));
[[maybe_unused]] const GrammageType length = model.GetInteractionLength(particle);
}
--
GitLab