diff --git a/corsika/detail/modules/sibyll/InteractionModel.inl b/corsika/detail/modules/sibyll/InteractionModel.inl
index 1af79a9d2ba646901a0d81f16289abe2a1a7c6cf..1541e8cb9c9cab72b71c9b0b84b90f79246639da 100644
--- a/corsika/detail/modules/sibyll/InteractionModel.inl
+++ b/corsika/detail/modules/sibyll/InteractionModel.inl
@@ -51,14 +51,14 @@ namespace corsika::sibyll {
targA = get_nucleus_A(targetId);
}
- if ((is_nucleus(targetId) &&
- (targA < minNuclearTargetA_ || targA >= maxTargetMassNumber_)) &&
- targetId != Code::Proton && targetId != Code::Hydrogen &&
- targetId != Code::Neutron) {
+ if (is_nucleus(targetId)) {
+ if (targA != 1 && (targA < minNuclearTargetA_ || targA >= maxTargetMassNumber_)) {
+ throw std::runtime_error("Target outside of allowed range for SIBYLL");
+ }
+ } else if (targetId != Code::Proton && targetId != Code::Neutron) {
throw std::runtime_error("Target cannot be handled by SIBYLL");
}
-
- if (is_nucleus(projectileId) && !corsika::sibyll::canInteract(projectileId)) {
+ if (is_nucleus(projectileId) || !corsika::sibyll::canInteract(projectileId)) {
throw std::runtime_error("Projectile cannot be handled by SIBYLL");
}
}
@@ -73,7 +73,8 @@ namespace corsika::sibyll {
double dummy, dum1, dum3, dum4, dumdif[3]; // dummies needed for fortran call
int const iBeam = corsika::sibyll::getSibyllXSCode(
- projectileId); // 0 (can not interact, 1: proton-like, 2: pion-like, 3:kaon-like)
+ projectileId); // 0 (can not interact, 1: proton-like, 2: pion-like,
+ // 3:kaon-like)
double const dEcm = sqrtSnn / 1_GeV;
// single nucleon target (p,n, hydrogen) or 4<=A<=18
diff --git a/corsika/detail/modules/sibyll/NuclearInteractionModel.inl b/corsika/detail/modules/sibyll/NuclearInteractionModel.inl
index 8811e0eeafc791e63d8e98f08462280ed56e4d0c..820839ff065697e42512eb43f54e1c07a98bb62d 100644
--- a/corsika/detail/modules/sibyll/NuclearInteractionModel.inl
+++ b/corsika/detail/modules/sibyll/NuclearInteractionModel.inl
@@ -49,7 +49,7 @@ namespace corsika::sibyll {
hadronicInteraction_.isValid(Code::Proton, targetId, sqrtSnn, 1, targetA); // throws
// projectile limits:
- if (is_nucleus(projectileId)) {
+ if (!is_nucleus(projectileId)) {
throw std::runtime_error("can only handle nuclear projectile");
}
if (projectileA >= getMaxNucleusAProjectile() || projectileA < 2) {
diff --git a/corsika/modules/sibyll/InteractionModel.hpp b/corsika/modules/sibyll/InteractionModel.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..4cba07fa5ee3ae4a920dca118a64d187c890853e
--- /dev/null
+++ b/corsika/modules/sibyll/InteractionModel.hpp
@@ -0,0 +1,125 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+#include <tuple>
+
+namespace corsika::sibyll {
+
+ /**
+ * @brief sibyll::InteractionModel provides the SIBYLL proton-nucleus interaction model.
+ *
+ * This is a TModel argument for InteractionProcess<TModel>.
+ */
+
+ class InteractionModel {
+
+ public:
+ InteractionModel();
+ ~InteractionModel();
+
+ /**
+ * @brief Set the Verbose flag.
+ *
+ * If flag is true, SIBYLL will printout additional secondary particle information
+ * lists, etc.
+ *
+ * @param flag to switch.
+ */
+ void setVerbose(bool const flag);
+
+ /**
+ * @brief evaluated validity of collision system.
+ *
+ * sibyll only accepts nuclei with 4<=A<=18 as targets, or protons aka Hydrogen or
+ * neutrons (p,n == nucleon).
+ */
+ void constexpr isValid(Code const projectileId, Code const targetId,
+ HEPEnergyType const sqrtSnn, unsigned int const projectileA,
+ unsigned int const targetA) const;
+
+ /**
+ * @brief returns inelastic AND elastic cross sections.
+ *
+ * These cross sections must correspond to the process described in doInteraction
+ * AND elastic scattering (sigma_tot = sigma_inel + sigma_el). Allowed targets are:
+ * nuclei or single nucleons (p,n,hydrogen). This "InelEla" method is used since
+ * Sibyll must be useful inside the NuclearInteraction model, which requires that.
+ *
+ * @param projectile is the Code of the projectile
+ * @param target is the Code of the target
+ * @param sqrtSnn is the center-of-mass energy (per nucleon pair)
+ * @param Aprojectil is the mass number of the projectils, if it is a nucleus
+ * @param Atarget is the mass number of the target, if it is a nucleus
+ *
+ * @return a tuple of: inelastic cross section, elastic cross section
+ */
+ std::tuple<CrossSectionType, CrossSectionType> getCrossSectionInelEla(
+ Code const projectile, Code const target, HEPEnergyType const sqrtSnn,
+ unsigned int const Aprojectile = 1, unsigned int const Atarget = 1) const;
+
+ /**
+ * @brief returns inelastic (production) cross section.
+ *
+ * This cross section must correspond to the process described in doInteraction.
+ * Allowed targets are: nuclei or single nucleons (p,n,hydrogen).
+ *
+ * @param projectile is the Code of the projectile
+ * @param target is the Code of the target
+ * @param sqrtSnn is the center-of-mass energy (per nucleon pair)
+ * @param Aprojectil is the mass number of the projectils, if it is a nucleus
+ * @param Atarget is the mass number of the target, if it is a nucleus
+ *
+ * @return inelastic cross section
+ * elastic cross section
+ */
+ CrossSectionType getCrossSection(Code const projectile, Code const target,
+ HEPEnergyType const sqrtSnn,
+ unsigned int const Aprojectile = 1,
+ unsigned int const Atarget = 1) const {
+ return std::get<0>(
+ getCrossSectionInelEla(projectile, target, sqrtSnn, Aprojectile, Atarget));
+ }
+
+ /**
+ * In this function SIBYLL is called to produce one event. The
+ * event is copied (and boosted) into the shower lab frame.
+ */
+
+ template <typename TSecondaries>
+ void doInteraction(TSecondaries&, COMBoost const& boost, Code const projectile,
+ Code const target, HEPEnergyType const sqrtSnn,
+ unsigned int const Aprojectile = 1,
+ unsigned int const Atarget = 1);
+
+ private:
+ HEPEnergyType constexpr getMinEnergyCoM() const { return minEnergyCoM_; }
+ HEPEnergyType constexpr getMaxEnergyCoM() const { return maxEnergyCoM_; }
+
+ // hard model limits
+ static HEPEnergyType constexpr minEnergyCoM_ = 10. * 1e9 * electronvolt;
+ static HEPEnergyType constexpr maxEnergyCoM_ = 1.e6 * 1e9 * electronvolt;
+ static unsigned int constexpr maxTargetMassNumber_ = 18;
+ static unsigned int constexpr minNuclearTargetA_ = 4;
+
+ default_prng_type& RNG_ = RNGManager<>::getInstance().getRandomStream("sibyll");
+
+ // data members
+ int count_ = 0;
+ int nucCount_ = 0;
+ bool sibyll_listing_;
+ };
+
+} // namespace corsika::sibyll
+
+#include <corsika/detail/modules/sibyll/InteractionModel.inl>
diff --git a/corsika/modules/sibyll/NuclearInteractionModel.hpp b/corsika/modules/sibyll/NuclearInteractionModel.hpp
new file mode 100644
index 0000000000000000000000000000000000000000..61b2f807e68190fc1f948174fb3087175359f5e1
--- /dev/null
+++ b/corsika/modules/sibyll/NuclearInteractionModel.hpp
@@ -0,0 +1,78 @@
+/*
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
+
+namespace corsika::sibyll {
+
+ /**
+ * The sibyll::NuclearInteractionModel provides the SIBYLL semi superposition model.
+ *
+ * can transform a proton-nucleus interaction model into a nucleus-nucleus interaction
+ * model.
+ *
+ * @tparam TNucleonModel
+ */
+ template <class TEnvironment, class TNucleonModel>
+ class NuclearInteractionModel {
+
+ public:
+ NuclearInteractionModel(TNucleonModel&, TEnvironment const&);
+ ~NuclearInteractionModel();
+
+ void constexpr isValid(Code const projectileId, Code const targetId,
+ HEPEnergyType const sqrtSnn, unsigned int const projectileA,
+ unsigned int const targetA) const;
+
+ void initializeNuclearCrossSections();
+ void printCrossSectionTable(Code) const;
+ CrossSectionType readCrossSectionTable(int const, Code const,
+ HEPEnergyType const) const;
+ HEPEnergyType getMinEnergyPerNucleonCoM() const { return gMinEnergyPerNucleonCoM_; }
+ HEPEnergyType getMaxEnergyPerNucleonCoM() const { return gMaxEnergyPerNucleonCoM_; }
+ unsigned int constexpr getMaxNucleusAProjectile() const {
+ return gMaxNucleusAProjectile_;
+ }
+ unsigned int constexpr getMaxNFragments() const { return gMaxNFragments_; }
+ unsigned int constexpr getNEnergyBins() const { return gNEnBins_; }
+
+ CrossSectionType getCrossSection(Code const, Code const, HEPEnergyType const,
+ unsigned int const projectileA = 1,
+ unsigned int const targetA = 1) const;
+
+ template <typename TSecondaryView>
+ void doInteraction(TSecondaryView&, COMBoost const& boost, Code const, Code const,
+ HEPEnergyType const, unsigned int const projectileA = 1,
+ unsigned int const targetA = 1);
+
+ private:
+ int count_ = 0;
+ int nucCount_ = 0;
+
+ TEnvironment const& environment_;
+ TNucleonModel& hadronicInteraction_;
+ std::map<Code, int> targetComponentsIndex_;
+ default_prng_type& RNG_ = RNGManager<>::getInstance().getRandomStream("sibyll");
+ static unsigned int constexpr gNSample_ =
+ 500; // number of samples in MC estimation of cross section
+ static unsigned int constexpr gMaxNucleusAProjectile_ = 56;
+ static unsigned int constexpr gNEnBins_ = 6;
+ static unsigned int constexpr gMaxNFragments_ = 60;
+ // energy limits defined by table used for cross section in signuc.f
+ // 10**1 GeV to 10**6 GeV
+ static HEPEnergyType constexpr gMinEnergyPerNucleonCoM_ = 10. * 1e9 * electronvolt;
+ static HEPEnergyType constexpr gMaxEnergyPerNucleonCoM_ = 1.e6 * 1e9 * electronvolt;
+ };
+
+} // namespace corsika::sibyll
+
+#include <corsika/detail/modules/sibyll/NuclearInteractionModel.inl>
diff --git a/tests/modules/testSibyll.cpp b/tests/modules/testSibyll.cpp
index 8bc4879f93c320d63774fe5c7579e932eb705fc5..477f1899b1970b8f478933640427cd5ba9791b4a 100644
--- a/tests/modules/testSibyll.cpp
+++ b/tests/modules/testSibyll.cpp
@@ -146,8 +146,8 @@ TEST_CASE("SibyllInterface", "modules") {
CHECK_THROWS(model.isValid(Code::Proton, Code::Iron, 100_GeV, 1, 56));
CHECK_NOTHROW(model.isValid(Code::Proton, Code::Oxygen, 100_GeV, 1, 16));
// beam particles
- CHECK_NOTHROW(model.isValid(Code::Electron, Code::Oxygen, 100_GeV, 1, 1));
- CHECK_NOTHROW(model.isValid(Code::Nucleus, Code::Oxygen, 100_GeV, 1, 20));
+ CHECK_THROWS(model.isValid(Code::Electron, Code::Oxygen, 100_GeV, 1, 1));
+ CHECK_THROWS(model.isValid(Code::Nucleus, Code::Oxygen, 100_GeV, 1, 20));
// energy too low
CHECK_THROWS(model.isValid(Code::Proton, Code::Proton, 9_GeV, 1, 1));
CHECK_NOTHROW(model.isValid(Code::Proton, Code::Proton, 11_GeV, 1, 1));