Merge branch '201-cascade_example-with-nucleus-primary-doesn-t-work' into 'master'

Resolve "cascade_example with nucleus primary doesn't work"

Closes #201

See merge request AirShowerPhysics/corsika!146

Framework/Particles/ParticleProperties.h

@@ -120,14 +120,22 @@ namespace corsika::particles {
   }
 
   int constexpr GetNucleusA(Code const p) {
+    if (p == Code::Nucleus) {
+      throw std::runtime_error("GetNucleusA(Code::Nucleus) is impossible!");
+    }
     return detail::nucleusA[static_cast<CodeIntType>(p)];
   }
 
   int constexpr GetNucleusZ(Code const p) {
+    if (p == Code::Nucleus) {
+      throw std::runtime_error("GetNucleusZ(Code::Nucleus) is impossible!");
+    }
     return detail::nucleusZ[static_cast<CodeIntType>(p)];
   }
 
-  bool constexpr IsNucleus(Code const p) { return GetNucleusA(p) != 0; }
+  bool constexpr IsNucleus(Code const p) {
+    return (p == Code::Nucleus) || (GetNucleusA(p) != 0);
+  }
 
   /**
    * the output operator for humand-readable particle codes

Framework/Particles/testParticles.cc

@@ -102,6 +102,7 @@ TEST_CASE("ParticleProperties", "[Particles]") {
     REQUIRE(IsHadron(Code::Proton));
     REQUIRE(IsHadron(Code::PiPlus));
     REQUIRE(IsHadron(Code::Oxygen));
+    REQUIRE(IsHadron(Code::Nucleus));
   }
 
   SECTION("Particle groups: muons") {
@@ -124,7 +125,6 @@ TEST_CASE("ParticleProperties", "[Particles]") {
     REQUIRE_FALSE(IsNeutrino(Code::Oxygen));
   }
 
-
   SECTION("Nuclei") {
     REQUIRE_FALSE(IsNucleus(Code::Gamma));
     REQUIRE(IsNucleus(Code::Argon));
@@ -137,5 +137,8 @@ TEST_CASE("ParticleProperties", "[Particles]") {
     REQUIRE(GetNucleusA(Code::Tritium) == 3);
     REQUIRE(Hydrogen::GetNucleusZ() == 1);
     REQUIRE(Tritium::GetNucleusA() == 3);
+
+    REQUIRE_THROWS(GetNucleusA(Code::Nucleus));
+    REQUIRE_THROWS(GetNucleusZ(Code::Nucleus));
   }
 }

Processes/Sibyll/NuclearInteraction.cc

@@ -227,15 +227,18 @@ namespace corsika::process::sibyll {
         RootCoordinateSystem::GetInstance().GetRootCoordinateSystem();
 
     const particles::Code corsikaBeamId = vP.GetPID();
-    if (!particles::IsNucleus(corsikaBeamId)) {
-      // no nuclear interaction
-      return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
+
+    if (corsikaBeamId != particles::Code::Nucleus) {
+      // check if target-style nucleus (enum), these are not allowed as projectile
+      if (particles::IsNucleus(corsikaBeamId))
+        throw std::runtime_error(
+            "NuclearInteraction: GetInteractionLength: Wrong nucleus type. Nuclear "
+            "projectiles should use NuclearStackExtension!");
+      else {
+        // no nuclear interaction
+        return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
+      }
     }
-    // check if target-style nucleus (enum)
-    if (corsikaBeamId != particles::Code::Nucleus)
-      throw std::runtime_error(
-          "NuclearInteraction: GetInteractionLength: Wrong nucleus type. Nuclear "
-          "projectiles should use NuclearStackExtension!");
 
     // read from cross section code table
 
@@ -333,13 +336,6 @@ namespace corsika::process::sibyll {
     //      const auto ProjMass = vP.GetMass();
     cout << "NuclearInteraction: DoInteraction: called with:" << ProjId << endl;
 
-    if (!IsNucleus(ProjId)) {
-      cout << "WARNING: non nuclear projectile in NUCLIB!" << endl;
-      // this should not happen
-      // throw std::runtime_error("Non nuclear projectile in NUCLIB!");
-      return process::EProcessReturn::eOk;
-    }
-
     // check if target-style nucleus (enum)
     if (ProjId != particles::Code::Nucleus)
       throw std::runtime_error(