Merge branch '403-urqmd-cross-section-regression' into 'master'

Resolve "UrQMD cross-section regression"

Closes #403

See merge request !335

corsika/detail/framework/utility/CorsikaData.inl

@@ -7,14 +7,15 @@
  */
 #pragma once
 
+#include <boost/filesystem/path.hpp>
+
 #include <cstdlib>
 #include <stdexcept>
 #include <string>
 
-inline std::string corsika::corsika_data(std::string const& key) {
+inline boost::filesystem::path corsika::corsika_data(boost::filesystem::path const& key) {
   if (auto const* p = std::getenv("CORSIKA_DATA"); p != nullptr) {
-    auto const path = std::string(p) + "/" + key;
-    return path;
+    return boost::filesystem::path(p) / key;
   } else {
     throw std::runtime_error("CORSIKA_DATA not set");
   }

corsika/detail/modules/qgsjetII/Interaction.inl

@@ -239,10 +239,14 @@ namespace corsika::qgsjetII {
       if (is_nucleus(targetCode)) { // nucleus
         targetMassNumber = get_nucleus_A(targetCode);
         if (targetMassNumber > maxMassNumber_)
-          throw std::runtime_error("QgsjetII target mass outside range.");
+          throw std::runtime_error(
+              "QgsjetII target mass outside range."); // LCOV_EXCL_LINE there is no
+                                                      // allowed path here
       } else {
-        if (targetCode != Proton::code)
-          throw std::runtime_error("QgsjetII Taget not possible.");
+        if (targetCode != Proton::code) // LCOV_EXCL_LINE there is no allowed path here
+          throw std::runtime_error(
+              "QgsjetII Taget not possible."); // LCOV_EXCL_LINE there is no allowed path
+                                               // here
       }
       CORSIKA_LOG_DEBUG("Interaction: target qgsjetII code/A: {}", targetMassNumber);
 
@@ -252,7 +256,9 @@ namespace corsika::qgsjetII {
       if (qgsjet_hadron_type == QgsjetIIHadronType::NucleusType) {
         projectileMassNumber = projectile.getNuclearA();
         if (projectileMassNumber > maxMassNumber_)
-          throw std::runtime_error("QgsjetII projectile mass outside range.");
+          throw std::runtime_error(
+              "QgsjetII projectile mass outside range."); // LCOV_EXCL_LINE there is no
+                                                          // allowed path here
         std::array<QgsjetIIHadronType, 2> constexpr nucleons = {
             QgsjetIIHadronType::ProtonType, QgsjetIIHadronType::NeutronType};
         std::uniform_int_distribution select(0, 1);

corsika/detail/modules/urqmd/UrQMD.inl

@@ -15,15 +15,112 @@
 #include <corsika/framework/geometry/QuantityVector.hpp>
 #include <corsika/framework/geometry/Vector.hpp>
 
+#include <boost/filesystem.hpp>
+#include <boost/multi_array.hpp>
+
 #include <algorithm>
+#include <cassert>
 #include <functional>
 #include <iostream>
+#include <fstream>
+#include <sstream>
 
 #include <urqmd.hpp>
 
 namespace corsika::urqmd {
 
-  inline UrQMD::UrQMD() { ::urqmd::iniurqmdc8_(); }
+  inline UrQMD::UrQMD(boost::filesystem::path const& xs_file) {
+    readXSFile(xs_file);
+    ::urqmd::iniurqmdc8_();
+  }
+
+  inline CrossSectionType UrQMD::getTabulatedCrossSection(Code projectileCode,
+                                                          Code targetCode,
+                                                          HEPEnergyType labEnergy) const {
+    // translated to C++ from CORSIKA 7 subroutine cxtot_u
+
+    auto const kinEnergy = labEnergy - get_mass(projectileCode);
+
+    assert(kinEnergy >= HEPEnergyType::zero());
+
+    double const logKinEnergy = std::log10(kinEnergy * (1 / 1_GeV));
+    double const ye = std::max(10 * logKinEnergy + 10.5, 1.);
+    int const je = std::min(int(ye), int(xs_interp_support_table_.shape()[2] - 2));
+    std::array<double, 3> w;
+    w[2 - 1] = ye - je;
+    w[3 - 1] = w[2 - 1] * (w[2 - 1] - 1.) * .5;
+    w[1 - 1] = 1 - w[2 - 1] + w[3 - 1];
+    w[2 - 1] = w[2 - 1] - 2 * w[3 - 1];
+
+    int projectileIndex;
+    switch (projectileCode) {
+      case Code::Proton:
+        projectileIndex = 0;
+        break;
+      case Code::AntiProton:
+        projectileIndex = 1;
+        break;
+      case Code::Neutron:
+        projectileIndex = 2;
+        break;
+      case Code::AntiNeutron:
+        projectileIndex = 3;
+        break;
+      case Code::PiPlus:
+        projectileIndex = 4;
+        break;
+      case Code::PiMinus:
+        projectileIndex = 5;
+        break;
+      case Code::KPlus:
+        projectileIndex = 6;
+        break;
+      case Code::KMinus:
+        projectileIndex = 7;
+        break;
+      case Code::K0Short:
+      case Code::K0Long:
+      /* since K0Short and K0Long are treated the same, we can also add K0 and K0Bar
+       * to the list. This is a deviation from CORSIKA 7. */
+      case Code::K0:
+      case Code::K0Bar:
+        projectileIndex = 8;
+        break;
+      default: { // LCOV_EXCL_START since this can never happen due to canInteract
+        CORSIKA_LOG_WARN("UrQMD cross-section not tabulated for {}", projectileCode);
+        return CrossSectionType::zero(); // LCOV_EXCL_STOP
+      }
+    }
+
+    int targetIndex;
+    switch (targetCode) {
+      case Code::Nitrogen:
+        targetIndex = 0;
+        break;
+      case Code::Oxygen:
+        targetIndex = 1;
+        break;
+      case Code::Argon:
+        targetIndex = 2;
+        break;
+      default:
+        std::stringstream ss;
+        ss << "UrQMD cross-section not tabluated for target " << targetCode;
+        throw std::runtime_error(ss.str().data());
+    }
+
+    auto result = CrossSectionType::zero();
+    for (int i = 0; i < 3; ++i) {
+      result +=
+          xs_interp_support_table_[projectileIndex][targetIndex][je + i - 1 - 1] * w[i];
+    }
+
+    CORSIKA_LOG_TRACE(
+        "UrQMD::GetTabulatedCrossSection proj={}, targ={}, E={} GeV, sigma={}",
+        get_name(projectileCode), get_name(targetCode), labEnergy / 1_GeV, result);
+
+    return result;
+  }
 
   inline CrossSectionType UrQMD::getCrossSection(Code vProjectileCode, Code vTargetCode,
                                                  HEPEnergyType vLabEnergy,
@@ -63,20 +160,19 @@ namespace corsika::urqmd {
 
   template <typename TParticle> // need template here, as this is called both with
                                 // SetupParticle as well as SetupProjectile
-  inline CrossSectionType UrQMD::getCrossSection(TParticle const& vProjectile,
-                                                 Code vTargetCode) const {
-    // TODO: return 0 for non-hadrons?
-
-    auto const projectileCode = vProjectile.getPID();
-    auto const projectileEnergyLab = vProjectile.getEnergy();
-
-    if (projectileCode == Code::K0Long) {
-      return 0.5 * (getCrossSection(Code::K0, vTargetCode, projectileEnergyLab) +
-                    getCrossSection(Code::K0Bar, vTargetCode, projectileEnergyLab));
+  inline CrossSectionType UrQMD::getCrossSection(TParticle const& projectile,
+                                                 Code targetCode) const {
+    auto const projectileCode = projectile.getPID();
+
+    if (projectileCode == Code::Nucleus) {
+      /*
+       * unfortunately unavoidable at the moment until we have tools to get the actual
+       * inealstic cross-section from UrQMD
+       */
+      return CrossSectionType::zero();
     }
 
-    int const Ap = (projectileCode == Code::Nucleus) ? vProjectile.getNuclearA() : 1;
-    return getCrossSection(projectileCode, vTargetCode, projectileEnergyLab, Ap);
+    return getTabulatedCrossSection(projectileCode, targetCode, projectile.getEnergy());
   }
 
   inline bool UrQMD::canInteract(Code vCode) const {
@@ -85,10 +181,9 @@ namespace corsika::urqmd {
     // only the usual long-lived species as input.
     // TODO: Charmed mesons should be added to the list, too
 
-    static Code const validProjectileCodes[] = {
-        Code::Nucleus,     Code::Proton, Code::AntiProton, Code::Neutron,
-        Code::AntiNeutron, Code::PiPlus, Code::PiMinus,    Code::KPlus,
-        Code::KMinus,      Code::K0,     Code::K0Bar,      Code::K0Long};
+    static std::array constexpr validProjectileCodes{
+        Code::Proton,  Code::AntiProton, Code::Neutron, Code::AntiNeutron, Code::PiPlus,
+        Code::PiMinus, Code::KPlus,      Code::KMinus,  Code::K0Short,     Code::K0Long};
 
     return std::find(std::cbegin(validProjectileCodes), std::cend(validProjectileCodes),
                      vCode) != std::cend(validProjectileCodes);
@@ -301,4 +396,41 @@ namespace corsika::urqmd {
     return mapPDGToUrQMD.at(static_cast<int>(get_PDG(code)));
   }
 
+  inline void UrQMD::readXSFile(boost::filesystem::path const& filename) {
+    boost::filesystem::ifstream file(filename, std::ios::in);
+
+    if (!file.is_open()) {
+      throw std::runtime_error(
+          filename.native() +
+          " could not be opened."); // LCOV_EXCL_LINE since this is pointless to test
+    }
+
+    std::string line;
+
+    std::getline(file, line);
+    std::stringstream ss(line);
+
+    char dummy;
+    int nTargets, nProjectiles, nSupports;
+    ss >> dummy >> nTargets >> nProjectiles >> nSupports;
+
+    decltype(xs_interp_support_table_)::extent_gen extents;
+    xs_interp_support_table_.resize(extents[nProjectiles][nTargets][nSupports]);
+
+    for (int i = 0; i < nTargets; ++i) {
+      for (int j = 0; j < nProjectiles; ++j) {
+        for (int k = 0; k < nSupports; ++k) {
+          std::getline(file, line);
+          std::stringstream s(line);
+          double energy, sigma;
+          s >> energy >> sigma;
+          xs_interp_support_table_[j][i][k] = sigma * 1_mb;
+        }
+
+        std::getline(file, line);
+        std::getline(file, line);
+      }
+    }
+  }
+
 } // namespace corsika::urqmd

corsika/framework/utility/COMBoost.hpp

@@ -19,14 +19,36 @@
 namespace corsika {
 
   /**
+     @defgroup Utilities
+
+     Collection of classes and methods to perform recurring tasks.
+   **/
+
+  /**
+     @class
+     @ingroup Utilities
+
      This utility class handles Lorentz boost between different
      referenence frames, using FourVector.
+
+     The class is initialized with projectile and optionally target
+     energy/momentum data. During initialization, a rotation matrix is
+     calculated to represent the projectile movement (and thus the
+     boost) along the z-axis. Also the inverse of this rotation is
+     calculated. The Lorentz boost matrix and its inverse are
+     determined as 2x2 matrices considering the energy and
+     pz-momentum.
+
+     Different constructors are offered with different specialization
+     for the cases of collisions (projectile-target) or just decays
+     (projectile only).
    */
 
   class COMBoost {
 
   public:
-    //! construct a COMBoost given four-vector of projectile and mass of target
+    //! construct a COMBoost given four-vector of projectile and mass of target (target at
+    //! rest)
     COMBoost(FourVector<HEPEnergyType, MomentumVector> const& Pprojectile,
              HEPEnergyType const massTarget);
 
@@ -41,9 +63,11 @@ namespace corsika {
     template <typename FourVector>
     FourVector fromCoM(FourVector const& p) const;
 
+    //! returns the rotated coordinate system
     CoordinateSystemPtr getRotatedCS() const;
 
   protected:
+    //! internal method
     void setBoost(double coshEta, double sinhEta);
 
   private:

corsika/framework/utility/CorsikaData.hpp

@@ -8,13 +8,18 @@
 
 #pragma once
 
-#include <string>
+#include <boost/filesystem/path.hpp>
 
 namespace corsika {
   /**
+   * @file CorsikaData.hpp
+   * @ingroup Utilities
+   * @{
    * returns the full path of the file \p filename within the CORSIKA_DATA directory
    */
-  std::string corsika_data(std::string const& filename);
+  boost::filesystem::path corsika_data(boost::filesystem::path const& filename);
+
+  //! @}
 } // namespace corsika
 
 #include <corsika/detail/framework/utility/CorsikaData.inl>

corsika/framework/utility/QuarticSolver.hpp

@@ -11,6 +11,9 @@
 #include <complex>
 
 /**
+ * @file QuarticSolver.hpp
+ * @ingroup Utilities
+ * @{
  * \todo convert to class
  */
 
@@ -56,6 +59,8 @@ namespace corsika::quartic_solver {
   // afterwards.
   DComplex* solve_quartic(double a, double b, double c, double d);
 
+  //! @}
+
 } // namespace corsika::quartic_solver
 
 #include <corsika/detail/framework/utility/QuarticSolver.inl>

corsika/framework/utility/SaveBoostHistogram.hpp

@@ -13,6 +13,8 @@
 namespace corsika {
 
   /**
+   * @ingroup Utilities
+   *
    * This functions saves a boost::histogram into a numpy file. Only rather basic axis
    * types are supported: regular, variable, integer, category<int>. Only "ordinary" bin
    * counts (i.e. a double or int) are supported, nothing fancy like profiles.

corsika/modules/urqmd/UrQMD.hpp

@@ -12,21 +12,28 @@
 #include <corsika/framework/core/PhysicalUnits.hpp>
 #include <corsika/framework/process/InteractionProcess.hpp>
 #include <corsika/framework/random/RNGManager.hpp>
+#include <corsika/framework/utility/CorsikaData.hpp>
 
 #include <corsika/setup/SetupStack.hpp>
 
+#include <boost/filesystem/path.hpp>
+#include <boost/multi_array.hpp>
+
 #include <array>
 #include <utility>
+#include <string>
 
 namespace corsika::urqmd {
 
   class UrQMD : public InteractionProcess<UrQMD> {
   public:
-    UrQMD();
+    UrQMD(boost::filesystem::path const& path = corsika_data("UrQMD/UrQMD-1.3.1-xs.dat"));
 
     template <typename TParticle>
     GrammageType getInteractionLength(TParticle const&) const;
 
+    CrossSectionType getTabulatedCrossSection(Code, Code, HEPEnergyType) const;
+
     template <typename TParticle>
     CrossSectionType getCrossSection(TParticle const&, Code) const;
 
@@ -39,11 +46,12 @@ namespace corsika::urqmd {
 
   private:
     static CrossSectionType getCrossSection(Code, Code, HEPEnergyType, int);
+    void readXSFile(boost::filesystem::path const&);
 
     // data members
     default_prng_type& RNG_ = RNGManager::getInstance().getRandomStream("urqmd");
-
     std::uniform_int_distribution<int> booleanDist_{0, 1};
+    boost::multi_array<CrossSectionType, 3> xs_interp_support_table_;
   };
 
   /**

documentation/index.rst

@@ -13,6 +13,7 @@ Welcome to the CORSIKA 8 air shower simulation framework.
    units
    environment
    stack
+   utilities
    api
    
 

documentation/utilities.rst

+Utilities
+==========
+
+.. doxygengroup:: Utilities
+   :project: CORSIKA8
+   :members:
+
+
+      

modules/urqmd/urqmdInterface.F

@@ -432,16 +432,3 @@ c~       xsegymin=0.25d0
 
       end
 
-c
-c M. Reininghaus, 2020-04-08
-c
-      integer function ReadSigmaLn(ia, ib, ic)
-      implicit none
-
-      include 'comres.f'
-
-      integer :: ia, ib, ic
-
-      ReadSigmaLn = SigmaLn(ia, ib, ic)
-
-      end function ReadSigmaLn

tests/modules/testQGSJetII.cpp

@@ -175,7 +175,7 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
   SECTION("InteractionInterface Nuclei") {
 
     auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
-        Code::Nucleus, 20, 10, 10100_GeV,
+        Code::Nucleus, 60, 30, 20100_GeV,
         (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
     setup::StackView& view = *(secViewPtr.get());
     auto particle = stackPtr->first();
@@ -183,10 +183,15 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
     auto const projectileMomentum = projectile.getMomentum();
 
     corsika::qgsjetII::Interaction model;
-    model.doInteraction(view);
+    model.doInteraction(view); // this also should produce some fragments
+    CHECK(view.getSize() == Approx(188).margin(2)); // this is not physics validation
+    int countFragments = 0;
+    for (auto const& sec : view) { countFragments += (sec.getPID() == Code::Nucleus); }
+    CHECK(countFragments == Approx(2).margin(1)); // this is not physics validation
     [[maybe_unused]] const GrammageType length = model.getInteractionLength(particle);
 
-    CHECK(length / (1_g / square(1_cm)) == Approx(20.13).margin(0.1));
+    CHECK(length / (1_g / square(1_cm)) ==
+          Approx(12).margin(2)); // this is not physics validation
   }
 
   SECTION("Heavy nuclei") {
@@ -210,17 +215,50 @@ TEST_CASE("QgsjetIIInterface", "[processes]") {
   }
 
   SECTION("Allowed Particles") {
-
-    auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
-        Code::Electron, 0, 0, 1100_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
-        *csPtr);
-    auto particle = stackPtr->first();
-    auto projectile = secViewPtr->getProjectile();
-    auto const projectileMomentum = projectile.getMomentum();
-
-    corsika::qgsjetII::Interaction model;
-
-    GrammageType const length = model.getInteractionLength(particle);
-    CHECK(length / (1_g / square(1_cm)) == std::numeric_limits<double>::infinity());
+    { // electron
+      auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+          Code::Electron, 0, 0, 100_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+          *csPtr);
+      auto particle = stackPtr->first();
+      corsika::qgsjetII::Interaction model;
+      GrammageType const length = model.getInteractionLength(particle);
+      CHECK(length / (1_g / square(1_cm)) == std::numeric_limits<double>::infinity());
+    }
+    { // pi0 is internally converted into pi+/pi-
+      auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+          Code::Pi0, 0, 0, 1000_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+          *csPtr);
+      setup::StackView& view = *(secViewPtr.get());
+      corsika::qgsjetII::Interaction model;
+      model.doInteraction(view);
+      CHECK(view.getSize() == Approx(18).margin(2)); // this is not physics validation
+    }
+    { // rho0 is internally converted into pi-/pi+
+      auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+          Code::Rho0, 0, 0, 1000_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+          *csPtr);
+      setup::StackView& view = *(secViewPtr.get());
+      corsika::qgsjetII::Interaction model;
+      model.doInteraction(view);
+      CHECK(view.getSize() == Approx(7).margin(2)); // this is not physics validation
+    }
+    { // Lambda is internally converted into neutron
+      auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+          Code::Lambda0, 0, 0, 100_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+          *csPtr);
+      setup::StackView& view = *(secViewPtr.get());
+      corsika::qgsjetII::Interaction model;
+      model.doInteraction(view);
+      CHECK(view.getSize() == Approx(25).margin(3)); // this is not physics validation
+    }
+    { // AntiLambda is internally converted into anti neutron
+      auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+          Code::Lambda0Bar, 0, 0, 100_GeV,
+          (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
+      setup::StackView& view = *(secViewPtr.get());
+      corsika::qgsjetII::Interaction model;
+      model.doInteraction(view);
+      CHECK(view.getSize() == Approx(25).margin(3)); // this is not physics validation
+    }
   }
 }

tests/modules/testUrQMD.cpp

@@ -77,24 +77,41 @@ TEST_CASE("UrQMD") {
     auto const& cs = *csPtr;
     { [[maybe_unused]] auto const& env_dummy = env; }
 
-    Code validProjectileCodes[] = {Code::PiPlus,  Code::PiMinus, Code::Proton,
-                                   Code::Neutron, Code::KPlus,   Code::KMinus,
-                                   Code::K0,      Code::K0Bar,   Code::K0Long};
+    Code validProjectileCodes[] = {Code::PiPlus,     Code::PiMinus,     Code::Proton,
+                                   Code::AntiProton, Code::AntiNeutron, Code::Neutron,
+                                   Code::KPlus,      Code::KMinus,      Code::K0,
+                                   Code::K0Bar,      Code::K0Long};
 
     for (auto code : validProjectileCodes) {
-      auto [stack, view] = setup::testing::setup_stack(
-          code, 0, 0, 100_GeV, (setup::Environment::BaseNodeType* const)nodePtr, cs);
+      auto [stack, view] = setup::testing::setup_stack(code, 0, 0, 100_GeV, nodePtr, cs);
       CHECK(stack->getEntries() == 1);
       CHECK(view->getEntries() == 0);
 
       // simple check whether the cross-section is non-vanishing
-      CHECK(urqmd.getCrossSection(view->getProjectile(), Code::Proton) / 1_mb > 0);
-      CHECK(urqmd.getCrossSection(view->getProjectile(), Code::Nitrogen) / 1_mb > 0);
-      CHECK(urqmd.getCrossSection(view->getProjectile(), Code::Oxygen) / 1_mb > 0);
-      CHECK(urqmd.getCrossSection(view->getProjectile(), Code::Argon) / 1_mb > 0);
+      // only nuclei with available tabluated data so far
+      CHECK(urqmd.getInteractionLength(stack->getNextParticle()) > 1_g / square(1_cm));
     }
   }
 
+  SECTION("targets options") {
+    auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Argon);
+    auto const& cs = *csPtr;
+    { [[maybe_unused]] auto const& env_dummy = env; }
+    auto [stack, view] =
+        setup::testing::setup_stack(Code::Proton, 0, 0, 100_GeV, nodePtr, cs);
+    CHECK(urqmd.getInteractionLength(stack->getNextParticle()) / 1_g * square(1_cm) ==
+          Approx(105).margin(5));
+  }
+
+  SECTION("invalid targets options") {
+    auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Omega);
+    auto const& cs = *csPtr;
+    { [[maybe_unused]] auto const& env_dummy = env; }
+    auto [stack, view] =
+        setup::testing::setup_stack(Code::Neutron, 0, 0, 100_GeV, nodePtr, cs);
+    CHECK_THROWS(urqmd.getInteractionLength(stack->getNextParticle()));
+  }
+
   SECTION("nucleus projectile") {
     auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
     [[maybe_unused]] auto const& env_dummy = env;      // against warnings
@@ -146,6 +163,45 @@ TEST_CASE("UrQMD") {
           Approx(0).margin(1e-2));
   }
 
+  SECTION("\"special\" projectile and target") {
+    {
+      auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Proton);
+      [[maybe_unused]] auto const& env_dummy = env;      // against warnings
+      [[maybe_unused]] auto const& node_dummy = nodePtr; // against warnings
+
+      auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+          Code::PiPlus, 0, 0, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+          *csPtr);
+      CHECK_THROWS(urqmd.doInteraction(*secViewPtr)); // Code::Proton not a valid target
+    }
+
+    {
+      auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
+      [[maybe_unused]] auto const& env_dummy = env;      // against warnings
+      [[maybe_unused]] auto const& node_dummy = nodePtr; // against warnings
+
+      auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+          Code::PiPlus, 0, 0, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+          *csPtr);
+      CHECK(stackPtr->getEntries() == 1);
+      CHECK(secViewPtr->getEntries() == 0);
+
+      // must be assigned to variable, cannot be used as rvalue?!
+      auto projectile = secViewPtr->getProjectile();
+      auto const projectileMomentum = projectile.getMomentum();
+
+      urqmd.doInteraction(*secViewPtr);
+
+      CHECK(sumCharge(*secViewPtr) ==
+            get_charge_number(Code::PiPlus) + get_charge_number(Code::Oxygen));
+
+      auto const secMomSum =
+          sumMomentum(*secViewPtr, projectileMomentum.getCoordinateSystem());
+      CHECK((secMomSum - projectileMomentum).getNorm() / projectileMomentum.getNorm() ==
+            Approx(0).margin(1e-2));
+    }
+  }
+
   SECTION("K0Long projectile") {
     auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Oxygen);
     [[maybe_unused]] auto const& env_dummy = env;      // against warnings
@@ -171,4 +227,22 @@ TEST_CASE("UrQMD") {
     CHECK((secMomSum - projectileMomentum).getNorm() / projectileMomentum.getNorm() ==
           Approx(0).margin(1e-2));
   }
+
+  SECTION("K0Short projectile") {
+    auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Argon);
+    [[maybe_unused]] auto const& env_dummy = env;      // against warnings
+    [[maybe_unused]] auto const& node_dummy = nodePtr; // against warnings
+
+    auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
+        Code::K0Short, 0, 0, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
+        *csPtr);
+    CHECK(stackPtr->getEntries() == 1);
+    CHECK(secViewPtr->getEntries() == 0);
+
+    // must be assigned to variable, cannot be used as rvalue?!
+    auto projectile = secViewPtr->getProjectile();
+    auto const projectileMomentum = projectile.getMomentum();
+
+    CHECK_THROWS(urqmd.doInteraction(*secViewPtr));
+  }
 }