infinite interaction length for non-hadrons

Processes/UrQMD/UrQMD.cc

@@ -7,6 +7,7 @@
 #include <corsika/setup/SetupTrajectory.h>
 #include <corsika/units/PhysicalUnits.h>
 
+#include <algorithm>
 #include <functional>
 #include <iostream>
 #include <random>
@@ -23,10 +24,8 @@ using SetupTrack = corsika::setup::Trajectory;
 
 template <typename TParticle> // need template here, as this is called both with
                               // SetupParticle as well as SetupProjectile
-                              CrossSectionType UrQMD::GetCrossSection(
-                                  TParticle const& vProjectile,
-                                  corsika::particles::Code vTargetCode)
-                                  const {
+CrossSectionType UrQMD::GetCrossSection(TParticle const& vProjectile,
+                                        corsika::particles::Code vTargetCode) const {
   using namespace units::si;
 
   // TODO: energy cuts, return 0 for non-hadrons
@@ -39,9 +38,10 @@ template <typename TParticle> // need template here, as this is called both with
       !IsNucleus(vTargetCode)) { // both particles are "special"
     auto const mProj = particles::GetMass(projectileCode);
     auto const mTar = particles::GetMass(vTargetCode);
-    double sqrtS = sqrt(units::si::detail::static_pow<2>(mProj) + units::si::detail::static_pow<2>(mTar) +
-                        2 * projectileEnergyLab * mTar) *
-                   (1 / 1_GeV);
+    double sqrtS =
+        sqrt(units::si::detail::static_pow<2>(mProj) +
+             units::si::detail::static_pow<2>(mTar) + 2 * projectileEnergyLab * mTar) *
+        (1 / 1_GeV);
 
     // we must set some UrQMD globals first...
     auto const [ityp, iso3] = ConvertToUrQMD(projectileCode);
@@ -66,7 +66,29 @@ template <typename TParticle> // need template here, as this is called both with
   }
 }
 
+bool UrQMD::CanInteract(particles::Code vCode) const {
+  // According to the manual, UrQMD can use all mesons, baryons and nucleons
+  // which are modeled also as input particles. I think it is safer to accept
+  // only the usual long-lived species as input.
+  // TODO: Charmed mesons should be added to the list, too
+
+  static particles::Code const validProjectileCodes[] = {
+      particles::Code::Nucleus, particles::Code::Proton,      particles::Code::AntiProton,
+      particles::Code::Neutron, particles::Code::AntiNeutron, particles::Code::PiPlus,
+      particles::Code::PiMinus, particles::Code::KPlus,       particles::Code::KMinus,
+      particles::Code::K0,      particles::Code::K0Bar};
+
+  return std::find(std::cbegin(validProjectileCodes), std::cend(validProjectileCodes),
+                   vCode) != std::cend(validProjectileCodes);
+}
+
 GrammageType UrQMD::GetInteractionLength(SetupParticle& vParticle, SetupTrack&) const {
+  if (!CanInteract(vParticle.GetPID())) {
+    // we could do the canInteract check in GetCrossSection, too but if
+    // we do it here we have the advantage of avoiding the loop
+    return std::numeric_limits<double>::infinity() * 1_g / (1_cm * 1_cm);
+  }
+
   auto const& mediumComposition =
       vParticle.GetNode()->GetModelProperties().GetNuclearComposition();
   using namespace std::placeholders;

Processes/UrQMD/UrQMD.h

@@ -19,12 +19,14 @@ namespace corsika::process::UrQMD {
         corsika::setup::Stack::StackIterator&, corsika::setup::Trajectory&) const;
 
     template <typename TParticle>
-    corsika::units::si::CrossSectionType GetCrossSection(
-        TParticle const&, corsika::particles::Code) const;
+    corsika::units::si::CrossSectionType GetCrossSection(TParticle const&,
+                                                         corsika::particles::Code) const;
 
     corsika::process::EProcessReturn DoInteraction(
         corsika::setup::StackView::StackIterator&);
 
+    bool CanInteract(particles::Code) const;
+
   private:
     corsika::random::RNG& fRNG =
         corsika::random::RNGManager::GetInstance().GetRandomStream("UrQMD");

Processes/UrQMD/testUrQMD.cc

@@ -42,10 +42,8 @@ using namespace corsika::units::si;
 template <typename TStackView>
 auto sumCharge(TStackView const& view) {
   int totalCharge = 0;
-  
-  for (auto const& p : view) {
-    totalCharge += particles::GetChargeNumber(p.GetPID());
-  }
+
+  for (auto const& p : view) { totalCharge += particles::GetChargeNumber(p.GetPID()); }
 
   return totalCharge;
 }
@@ -165,7 +163,7 @@ TEST_CASE("UrQMD") {
     auto [stackPtr, secViewPtr] = setupStack(A, Z, 400_GeV, nodePtr, *csPtr);
 
     // must be assigned to variable, cannot be used as rvalue?!
-    auto projectile =    secViewPtr     ->GetProjectile();
+    auto projectile = secViewPtr->GetProjectile();
     [[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(projectile);
 
     REQUIRE(sumCharge(*secViewPtr) ==
@@ -175,14 +173,14 @@ TEST_CASE("UrQMD") {
   SECTION("\"special\" projectile") {
     auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Oxygen);
     auto [stackPtr, secViewPtr] =
-        setupStack(particles::Code::PiPlus, 400_GeV, nodePtr, *csPtr);
+        setupStack(particles::Code::K0, 400_GeV, nodePtr, *csPtr);
 
     // must be assigned to variable, cannot be used as rvalue?!
-    auto projectile = secViewPtr->GetProjectile(); 
+    auto projectile = secViewPtr->GetProjectile();
     [[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(projectile);
 
     REQUIRE(sumCharge(*secViewPtr) ==
-            particles::GetChargeNumber(particles::Code::PiPlus) +
+            particles::GetChargeNumber(particles::Code::K0) +
                 particles::GetChargeNumber(particles::Code::Oxygen));
   }
 }