Merge branch '374-hydrogen-cause-error-in-sibyll-interaction' into 'master'

Resolve "Hydrogen cause error in Sibyll interaction"

Closes #374

See merge request AirShowerPhysics/corsika!312

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();

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;

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

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;