diff --git a/Processes/Sibyll/CMakeLists.txt b/Processes/Sibyll/CMakeLists.txt
index bf89db756b5aaca2be901596c1015b5d724245e7..215685c6c38255f72ae7dd21924d214b1d66d5a1 100644
--- a/Processes/Sibyll/CMakeLists.txt
+++ b/Processes/Sibyll/CMakeLists.txt
@@ -27,6 +27,7 @@ set (
SibStack.h
Decay.h
Interaction.h
+ NuclearInteraction.h
${PROJECT_BINARY_DIR}/Processes/Sibyll/Generated.inc
)
diff --git a/Processes/Sibyll/NuclearInteraction.h b/Processes/Sibyll/NuclearInteraction.h
new file mode 100644
index 0000000000000000000000000000000000000000..5f0d29781e1b5260f61cb65ef3a9e25f209f279b
--- /dev/null
+++ b/Processes/Sibyll/NuclearInteraction.h
@@ -0,0 +1,220 @@
+
+/**
+ * (c) Copyright 2018 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * See file AUTHORS for a list of contributors.
+ *
+ * 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.
+ */
+
+#ifndef _corsika_process_sibyll_nuclearinteraction_h_
+#define _corsika_process_sibyll_nuclearinteraction_h_
+
+#include <corsika/process/InteractionProcess.h>
+
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/NuclearComposition.h>
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/process/sibyll/ParticleConversion.h>
+#include <corsika/process/sibyll/SibStack.h>
+#include <corsika/process/sibyll/sibyll2.3c.h>
+#include <corsika/random/RNGManager.h>
+#include <corsika/units/PhysicalUnits.h>
+#include <corsika/utl/COMBoost.h>
+
+namespace corsika::process::sibyll {
+
+ class NuclearInteraction : public corsika::process::InteractionProcess<NuclearInteraction> {
+
+ int fCount = 0;
+ int fNucCount = 0;
+
+ public:
+ NuclearInteraction(corsika::environment::Environment const& env)
+ : fEnvironment(env) {}
+ ~NuclearInteraction() {
+ std::cout << "Sibyll::NuclearInteraction n=" << fCount << " Nnuc=" << fNucCount
+ << std::endl;
+ }
+
+ void Init() {
+
+ using corsika::random::RNGManager;
+ using std::cout;
+ using std::endl;
+
+ // initialize hadronic interaction module
+ sibyll_ini_();
+
+ // initialize nuclib
+ nuc_nuc_ini_();
+ }
+
+ std::tuple<corsika::units::si::CrossSectionType, int> GetCrossSection(
+ const corsika::particles::Code BeamId, const corsika::particles::Code TargetId,
+ const corsika::units::si::HEPEnergyType CoMenergy) {
+ using namespace corsika::units::si;
+ double sigProd, dummy, dum1, dum2, dum3, dum4;
+ double dumdif[3];
+
+ if(!corsika::particles::IsNucleus(BeamId)){
+ sigProd = std::numeric_limits<double>::infinity();
+ return std::make_tuple(sigProd * 1_mbarn, 1);
+ }
+
+ // TODO: use nuclib to calc. nuclear cross sections
+ // FOR NOW: use proton cross section for nuclei
+ auto const BeamIdToUse = corsika::particles::Proton::GetCode();
+ std::cout << "WARNING: replacing beam nucleus with proton!" << std::endl;
+
+ const int iBeam = process::sibyll::GetSibyllXSCode(BeamIdToUse);
+ const double dEcm = CoMenergy / 1_GeV;
+ if (corsika::particles::IsNucleus(TargetId)) {
+ const int iTarget = corsika::particles::GetNucleusA(TargetId);
+ if (iTarget > 18 || 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);
+ return std::make_tuple(sigProd * 1_mbarn, iTarget);
+ } else if (TargetId == corsika::particles::Proton::GetCode()) {
+ sib_sigma_hp_(iBeam, dEcm, dum1, dum2, sigProd, dumdif, dum3, dum4);
+ return std::make_tuple(sigProd * 1_mbarn, 1);
+ } else {
+ throw std::runtime_error("GetCrossSection: no interaction in sibyll possible");
+ // no interaction in sibyll possible, return infinite cross section? or throw?
+ sigProd = std::numeric_limits<double>::infinity();
+ return std::make_tuple(sigProd * 1_mbarn, 0);
+ }
+ }
+
+ template <typename Particle, typename Track>
+ corsika::units::si::GrammageType GetInteractionLength(Particle const& p, Track&) {
+
+ using namespace corsika::units;
+ using namespace corsika::units::si;
+ using namespace corsika::geometry;
+ using std::cout;
+ using std::endl;
+
+ // coordinate system, get global frame of reference
+ CoordinateSystem& rootCS =
+ RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
+
+ const particles::Code corsikaBeamId = p.GetPID();
+
+ if(!corsika::particles::IsNucleus(corsikaBeamId))
+ // no nuclear interaction
+ return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
+
+ // beam particles for sibyll : 1, 2, 3 for p, pi, k
+ // read from cross section code table
+
+ const HEPMassType nucleon_mass = 0.5 * (corsika::particles::Proton::GetMass() +
+ corsika::particles::Neutron::GetMass());
+
+ // FOR NOW: assume target is at rest
+ MomentumVector pTarget(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+
+ // total momentum and energy
+ HEPEnergyType Elab = p.GetEnergy() + nucleon_mass;
+ MomentumVector pTotLab(rootCS, {0.0_GeV, 0.0_GeV, 0.0_GeV});
+ pTotLab += p.GetMomentum();
+ pTotLab += pTarget;
+ auto const pTotLabNorm = pTotLab.norm();
+ // calculate cm. energy
+ const HEPEnergyType ECoM = sqrt(
+ (Elab + pTotLabNorm) * (Elab - pTotLabNorm)); // binomial for numerical accuracy
+
+ std::cout << "NuclearInteraction: LambdaInt: \n"
+ << " input energy: " << p.GetEnergy() / 1_GeV << endl
+ << " beam pid:" << p.GetPID() << endl;
+
+ if ( Elab >= 8.5_GeV && ECoM >= 10_GeV) {
+
+ // get target from environment
+ /*
+ the target should be defined by the Environment,
+ ideally as full particle object so that the four momenta
+ and the boosts can be defined..
+ */
+ const auto currentNode =
+ fEnvironment.GetUniverse()->GetContainingNode(p.GetPosition());
+ const auto mediumComposition =
+ currentNode->GetModelProperties().GetNuclearComposition();
+ // determine average interaction length
+ // weighted sum
+ int i = -1;
+ double avgTargetMassNumber = 0.;
+ si::CrossSectionType weightedProdCrossSection = 0_mbarn;
+ // get weights of components from environment/medium
+ const auto w = mediumComposition.GetFractions();
+ // loop over components in medium
+ for (auto const targetId : mediumComposition.GetComponents()) {
+ i++;
+ cout << "NuclearInteraction: get interaction length for target: " << targetId << endl;
+
+ auto const [productionCrossSection, numberOfNucleons] =
+ GetCrossSection(corsikaBeamId, targetId, ECoM);
+
+ std::cout << "NuclearInteraction: "
+ << " IntLength: sibyll return (mb): "
+ << productionCrossSection / 1_mbarn << std::endl;
+ weightedProdCrossSection += w[i] * productionCrossSection;
+ avgTargetMassNumber += w[i] * numberOfNucleons;
+ }
+ cout << "NuclearInteraction: "
+ << "IntLength: weighted CrossSection (mb): "
+ << weightedProdCrossSection / 1_mbarn << endl;
+
+ // calculate interaction length in medium
+ GrammageType const int_length =
+ avgTargetMassNumber * corsika::units::constants::u / weightedProdCrossSection;
+ std::cout << "NuclearInteraction: "
+ << "interaction length (g/cm2): "
+ << int_length / (0.001_kg) * 1_cm * 1_cm << std::endl;
+
+ return int_length;
+ } else {
+ return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
+ }
+ }
+
+ template <typename Particle, typename Stack>
+ corsika::process::EProcessReturn DoInteraction(Particle& p, Stack& s) {
+
+ using namespace corsika::units;
+ using namespace corsika::utl;
+ using namespace corsika::units::si;
+ using namespace corsika::geometry;
+ using std::cout;
+ using std::endl;
+
+ const auto corsikaBeamId = p.GetPID();
+ const auto kNucleus = IsNucleus(corsikaBeamId);
+ if(!kNucleus)
+ return process::EProcessReturn::eOk;
+
+
+ // add proton instead
+ auto pnew = s.NewParticle();
+ pnew.SetPID( corsika::particles::Proton::GetCode() );
+ pnew.SetMomentum( p.GetMomentum() );
+ pnew.SetEnergy( p.GetEnergy() );
+
+ // delete current particle
+ p.Delete();
+
+ return process::EProcessReturn::eOk;
+ }
+
+ private:
+ corsika::environment::Environment const& fEnvironment;
+ corsika::random::RNG& fRNG =
+ corsika::random::RNGManager::GetInstance().GetRandomStream("s_rndm");
+ };
+
+} // namespace corsika::process::sibyll
+
+#endif
diff --git a/Processes/Sibyll/sibyll2.3c.h b/Processes/Sibyll/sibyll2.3c.h
index 104b55eb0bf2e6fc12277f6d62bf1adcf567acb4..cd13ebb67bc9642310f6cf95f8194d6fc6b551f5 100644
--- a/Processes/Sibyll/sibyll2.3c.h
+++ b/Processes/Sibyll/sibyll2.3c.h
@@ -79,6 +79,9 @@ void sibyll_(const int&, const int&, const double&);
// subroutine to initiate sibyll
void sibyll_ini_();
+// subroutine to initiate nuclib
+void nuc_nuc_ini_();
+
// subroutine to SET DECAYS
void dec_ini_();
diff --git a/Processes/Sibyll/testSibyll.cc b/Processes/Sibyll/testSibyll.cc
index 9b078fd087c16d6238837ce26783468640777054..16449b87c96f57f9eba63c92f4aa30973b597978 100644
--- a/Processes/Sibyll/testSibyll.cc
+++ b/Processes/Sibyll/testSibyll.cc
@@ -11,6 +11,7 @@
#include <corsika/process/sibyll/Decay.h>
#include <corsika/process/sibyll/Interaction.h>
+#include <corsika/process/sibyll/NuclearInteraction.h>
#include <corsika/process/sibyll/ParticleConversion.h>
#include <corsika/random/RNGManager.h>
@@ -127,6 +128,20 @@ TEST_CASE("SibyllInterface", "[processes]") {
model.GetInteractionLength(particle, track);
}
+ SECTION("NuclearInteractionInterface") {
+
+ setup::Stack stack;
+ auto particle = stack.NewParticle();
+
+ Interaction model(env);
+
+ model.Init();
+ [[maybe_unused]] const process::EProcessReturn ret =
+ model.DoInteraction(particle, stack);
+ [[maybe_unused]] const GrammageType length =
+ model.GetInteractionLength(particle, track);
+ }
+
SECTION("DecayInterface") {
setup::Stack stack;