From cd112f8d23c56275d2ff16293a87ec0196169e89 Mon Sep 17 00:00:00 2001
From: Felix Riehn <friehn@lip.pt>
Date: Tue, 26 Jan 2021 13:01:43 +0100
Subject: [PATCH] Resolve "Hydrogen cause error in Sibyll interaction"
---
corsika/detail/modules/sibyll/Interaction.inl | 30 +++++++++++++------
corsika/modules/sibyll/Interaction.hpp | 11 ++++++-
examples/cascade_proton_example.cpp | 8 +++--
tests/modules/testSibyll.cpp | 25 ++++++++++++++++
4 files changed, 61 insertions(+), 13 deletions(-)
diff --git a/corsika/detail/modules/sibyll/Interaction.inl b/corsika/detail/modules/sibyll/Interaction.inl
index 364ba19d1..e2c5f001f 100644
--- a/corsika/detail/modules/sibyll/Interaction.inl
+++ b/corsika/detail/modules/sibyll/Interaction.inl
@@ -78,21 +78,33 @@ namespace corsika::sibyll {
const corsika::HEPEnergyType CoMenergy) const {
double sigProd, sigEla, dummy, dum1, dum3, dum4;
double dumdif[3];
- const int iBeam = corsika::sibyll::getSibyllXSCode(BeamId);
+ const int iBeam = corsika::sibyll::getSibyllXSCode(
+ BeamId); // 0 (can not interact, 1: proton-like, 2: pion-like, 3:kaon-like)
+ if (!iBeam)
+ throw std::runtime_error(
+ "Interaction: getCrossSection: interaction of beam hadron not defined in "
+ "Sibyll!");
if (!isValidCoMEnergy(CoMenergy)) {
throw std::runtime_error(
"Interaction: getCrossSection: CoM energy outside range for Sibyll!");
}
const double dEcm = CoMenergy / 1_GeV;
- if (corsika::is_nucleus(TargetId)) {
- const int iTarget = corsika::get_nucleus_A(TargetId);
- if (iTarget > maxTargetMassNumber_ || iTarget == 0)
- throw std::runtime_error(
- "Sibyll target outside range. Only nuclei with A<18 are allowed.");
- sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy, sigEla);
- } else if (TargetId == corsika::Code::Proton) {
- sib_sigma_hp_(iBeam, dEcm, dum1, sigEla, sigProd, dumdif, dum3, dum4);
+ // single nucleon target (p,n, hydrogen) or 4<=A<=18
+ if (isValidTarget(TargetId)) {
+ // single nucleon target
+ if (TargetId == corsika::Code::Proton || TargetId == Code::Hydrogen ||
+ TargetId == Code::Neutron) {
+ sib_sigma_hp_(iBeam, dEcm, dum1, sigEla, sigProd, dumdif, dum3, dum4);
+ } else {
+ // nuclear target
+ const int iTarget = corsika::get_nucleus_A(TargetId);
+ sib_sigma_hnuc_(iBeam, iTarget, dEcm, sigProd, dummy, sigEla);
+ }
} else {
+ // throw std::runtime_error(
+ // "Sibyll nuclear target outside range. Only nuclei with 4<=A<18 are
+ // allowed.");
+
// no interaction in sibyll possible, return infinite cross section? or throw?
sigProd = std::numeric_limits<double>::infinity();
sigEla = std::numeric_limits<double>::infinity();
diff --git a/corsika/modules/sibyll/Interaction.hpp b/corsika/modules/sibyll/Interaction.hpp
index 9b983eb01..605e6c184 100644
--- a/corsika/modules/sibyll/Interaction.hpp
+++ b/corsika/modules/sibyll/Interaction.hpp
@@ -25,10 +25,18 @@ namespace corsika::sibyll {
bool isValidCoMEnergy(HEPEnergyType const ecm) const {
return (minEnergyCoM_ <= ecm) && (ecm <= maxEnergyCoM_);
}
+ //! sibyll only accepts nuclei with 4<=A<=18 as targets, or protons aka Hydrogen or
+ //! neutrons (p,n == nucleon)
bool isValidTarget(Code const TargetId) const {
- return is_nucleus(TargetId) && (get_nucleus_A(TargetId) < maxTargetMassNumber_);
+ return (is_nucleus(TargetId) && (get_nucleus_A(TargetId) >= minNuclearTargetA_) &&
+ (get_nucleus_A(TargetId) < maxTargetMassNumber_)) ||
+ (TargetId == Code::Proton || TargetId == Code::Hydrogen ||
+ TargetId == Code::Neutron);
}
+ //! returns production and elastic cross section for hadrons in sibyll. Inputs are:
+ //! CorsikaId of beam particle, CorsikaId of target particle and center-of-mass
+ //! energy. Allowed targets are: nuclei or single nucleons (p,n,hydrogen).
std::tuple<CrossSectionType, CrossSectionType> getCrossSection(
Code const, Code const, HEPEnergyType const) const;
@@ -61,6 +69,7 @@ namespace corsika::sibyll {
const HEPEnergyType minEnergyCoM_ = 10. * 1e9 * electronvolt;
const HEPEnergyType maxEnergyCoM_ = 1.e6 * 1e9 * electronvolt;
const int maxTargetMassNumber_ = 18;
+ const int minNuclearTargetA_ = 4;
// data members
int count_ = 0;
diff --git a/examples/cascade_proton_example.cpp b/examples/cascade_proton_example.cpp
index 24ed9e77d..6381feff4 100644
--- a/examples/cascade_proton_example.cpp
+++ b/examples/cascade_proton_example.cpp
@@ -119,9 +119,9 @@ int main() {
RNGManager::getInstance().registerRandomStream("sibyll");
RNGManager::getInstance().registerRandomStream("pythia");
- // sibyll::Interaction sibyll(env);
+ // corsika::sibyll::Interaction sibyll;
corsika::pythia8::Interaction pythia;
- // sibyll::NuclearInteraction sibyllNuc(env, sibyll);
+ // sibyll::NuclearInteraction sibyllNuc(sibyll, env);
// sibyll::Decay decay;
corsika::pythia8::Decay decay;
ParticleCut cut(60_GeV, true, true);
@@ -135,7 +135,9 @@ int main() {
BetheBlochPDG eLoss{showerAxis};
// assemble all processes into an ordered process list
- // auto sequence = sibyll << decay << hadronicElastic << cut << trackWriter;
+ // auto sequence = make_sequence(sibyll, sibyllNuc, decay, eLoss, cut, trackWriter,
+ // stackInspect); auto sequence = make_sequence(sibyll, decay, eLoss, cut, trackWriter,
+ // stackInspect);
auto sequence = make_sequence(pythia, decay, eLoss, cut, trackWriter, stackInspect);
// define air shower object, run simulation
diff --git a/tests/modules/testSibyll.cpp b/tests/modules/testSibyll.cpp
index 6e35c6a1d..dc237b4a3 100644
--- a/tests/modules/testSibyll.cpp
+++ b/tests/modules/testSibyll.cpp
@@ -108,6 +108,31 @@ TEST_CASE("SibyllInterface", "[processes]") {
RNGManager::getInstance().registerRandomStream("sibyll");
+ SECTION("InteractionInterface - valid targets") {
+
+ Interaction model;
+ // sibyll only accepts protons or nuclei with 4<=A<=18 as targets
+ CHECK_FALSE(model.isValidTarget(Code::Electron));
+ CHECK(model.isValidTarget(Code::Hydrogen));
+ CHECK_FALSE(model.isValidTarget(Code::Deuterium));
+ CHECK(model.isValidTarget(Code::Helium));
+ CHECK_FALSE(model.isValidTarget(Code::Helium3));
+ CHECK_FALSE(model.isValidTarget(Code::Iron));
+ CHECK(model.isValidTarget(Code::Oxygen));
+
+ // hydrogen target == proton target == neutron target
+ auto const [xs_prod_pp, xs_ela_pp] =
+ model.getCrossSection(Code::Proton, Code::Proton, 100_GeV);
+ auto const [xs_prod_pn, xs_ela_pn] =
+ model.getCrossSection(Code::Proton, Code::Neutron, 100_GeV);
+ auto const [xs_prod_pHydrogen, xs_ela_pHydrogen] =
+ model.getCrossSection(Code::Proton, Code::Hydrogen, 100_GeV);
+ CHECK(xs_prod_pp == xs_prod_pHydrogen);
+ CHECK(xs_prod_pp == xs_prod_pn);
+ CHECK(xs_ela_pp == xs_ela_pHydrogen);
+ CHECK(xs_ela_pn == xs_ela_pHydrogen);
+ }
+
SECTION("InteractionInterface - low energy") {
const HEPEnergyType P0 = 60_GeV;
--
GitLab