very basic test with NN interactions

Processes/UrQMD/CMakeLists.txt

@@ -63,6 +63,7 @@ target_link_libraries (
   CORSIKAunits
   CORSIKAgeometry
   CORSIKArandom
+  CORSIKAsetup
   )
 
 target_include_directories (
@@ -90,6 +91,10 @@ add_executable (testUrQMD
 target_link_libraries (
   testUrQMD
   ProcessUrQMD
+  CORSIKAsetup
+  CORSIKArandom
+  CORSIKAgeometry
+  CORSIKAunits
   CORSIKAthirdparty # for catch2
   )
 CORSIKA_ADD_TEST(testUrQMD)

Processes/UrQMD/UrQMD.cc

@@ -17,7 +17,8 @@ using namespace corsika::units::si;
 UrQMD::UrQMD() { iniurqmd_(); }
 
 using SetupStack = corsika::setup::Stack;
-using SetupParticle = SetupStack::StackIterator;
+using SetupParticle = corsika::setup::Stack::StackIterator;
+using SetupProjectile = corsika::setup::StackView::StackIterator;
 using SetupTrack = corsika::setup::Trajectory;
 
 CrossSectionType UrQMD::GetCrossSection(
@@ -27,7 +28,6 @@ CrossSectionType UrQMD::GetCrossSection(
   return 10_mb; // TODO: implement
 }
 
-template <>
 GrammageType UrQMD::GetInteractionLength(SetupParticle& particle, SetupTrack&) const {
   auto const& mediumComposition =
       particle.GetNode()->GetModelProperties().GetNuclearComposition();
@@ -41,9 +41,7 @@ GrammageType UrQMD::GetInteractionLength(SetupParticle& particle, SetupTrack&) c
          weightedProdCrossSection;
 }
 
-template <>
-corsika::process::EProcessReturn UrQMD::DoInteraction(SetupParticle& projectile,
-                                                      SetupStack&) {
+corsika::process::EProcessReturn UrQMD::DoInteraction(SetupProjectile& projectile) {
   using namespace units::si;
 
   auto const projectileCode = projectile.GetPID();
@@ -104,7 +102,7 @@ corsika::process::EProcessReturn UrQMD::DoInteraction(SetupParticle& projectile,
     inputs_.spiso3[0] = iso3;
   }
 
-  rsys_.ebeam = (projectileEnergyLab - particles::GetMass(projectileCode)) * (1 / 1_GeV);
+  rsys_.ebeam = (projectileEnergyLab - projectile.GetMass()) * (1 / 1_GeV);
 
   // initilazation regarding target
   auto const& mediumComposition =
@@ -159,6 +157,8 @@ corsika::process::EProcessReturn UrQMD::DoInteraction(SetupParticle& projectile,
                    TimeType>{code, energy, momentum, projectilePosition, projectileTime});
   }
 
+  std::cout << "UrQMD generated " << sys_.npart << " secondaries!" << std::endl;
+
   if (sys_.npart > 0) // delete only in case of inelastic collision, otherwise keep
     projectile.Delete();
 }

Processes/UrQMD/UrQMD.h

@@ -4,6 +4,8 @@
 #include <corsika/particles/ParticleProperties.h>
 #include <corsika/process/InteractionProcess.h>
 #include <corsika/random/RNGManager.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
 #include <corsika/units/PhysicalUnits.h>
 
 #include <array>
@@ -13,16 +15,15 @@ namespace corsika::process::UrQMD {
   class UrQMD : public corsika::process::InteractionProcess<UrQMD> {
   public:
     UrQMD();
-
-    template <typename Particle, typename Track>
-    corsika::units::si::GrammageType GetInteractionLength(Particle&, Track&) const;
+    corsika::units::si::GrammageType GetInteractionLength(
+        corsika::setup::Stack::StackIterator&, corsika::setup::Trajectory&) const;
 
     corsika::units::si::CrossSectionType GetCrossSection(
         corsika::particles::Code, corsika::particles::Code,
         corsika::units::si::HEPEnergyType) const;
 
-    template <typename Particle, typename Stack>
-    corsika::process::EProcessReturn DoInteraction(Particle&, Stack&);
+    corsika::process::EProcessReturn DoInteraction(
+        corsika::setup::StackView::StackIterator&);
 
   private:
     corsika::random::RNG& fRNG =

Processes/UrQMD/testUrQMD.cc

@@ -9,11 +9,75 @@
  */
 
 #include <corsika/process/urqmd/UrQMD.h>
+#include <corsika/random/RNGManager.h>
+
+#include <corsika/geometry/Point.h>
+#include <corsika/geometry/RootCoordinateSystem.h>
+#include <corsika/geometry/Vector.h>
+
+#include <corsika/units/PhysicalConstants.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+
+#include <corsika/environment/Environment.h>
+#include <corsika/environment/HomogeneousMedium.h>
+#include <corsika/environment/NuclearComposition.h>
 
 #define CATCH_CONFIG_MAIN // This tells Catch to provide a main() - only do this in one
                           // cpp file
 #include <catch2/catch.hpp>
 
+using namespace corsika;
 using namespace corsika::process::UrQMD;
+using namespace corsika::units::si;
+
+TEST_CASE("UrQMD") {
+  corsika::random::RNGManager::GetInstance().RegisterRandomStream("UrQMD");
+  UrQMD urqmd;
+
+  // setup environment, geometry
+  environment::Environment<environment::IMediumModel> env;
+  auto& universe = *(env.GetUniverse());
+  const geometry::CoordinateSystem& cs = env.GetCoordinateSystem();
+
+  auto theMedium =
+      environment::Environment<environment::IMediumModel>::CreateNode<geometry::Sphere>(
+          geometry::Point{cs, 0_m, 0_m, 0_m},
+          1_km * std::numeric_limits<double>::infinity());
+
+  using MyHomogeneousModel = environment::HomogeneousMedium<environment::IMediumModel>;
+  theMedium->SetModelProperties<MyHomogeneousModel>(
+      1_kg / (1_m * 1_m * 1_m),
+      environment::NuclearComposition(
+          std::vector<particles::Code>{particles::Code::Oxygen}, std::vector<float>{1.}));
+
+  auto const* nodePtr = theMedium.get();
+  universe.AddChild(std::move(theMedium));
+
+  geometry::Point const origin(cs, {0_m, 0_m, 0_m});
+  geometry::Vector<units::si::SpeedType::dimension_type> v(cs, 0_m / second, 0_m / second,
+                                                           1_m / second);
+  geometry::Line line(origin, v);
+  geometry::Trajectory<geometry::Line> track(line, 10_s);
+
+  auto constexpr mN = corsika::units::constants::nucleonMass;
+
+  setup::Stack stack;
+  const HEPEnergyType P0 = 50_GeV;
+  unsigned short constexpr A = 56, Z = 26;
+  HEPMomentumType E0 = sqrt(A * A * mN * mN + P0 * P0);
+  auto plab = corsika::stack::MomentumVector(cs, {P0, 0_GeV, 0_GeV});
+  auto particle =
+      stack.AddParticle(std::tuple<particles::Code, units::si::HEPEnergyType,
+                                   corsika::stack::MomentumVector, geometry::Point,
+                                   units::si::TimeType, unsigned short, unsigned short>{
+          particles::Code::Nucleus, E0, plab, origin, 0_ns, A, Z}); // iron
+  particle.SetNode(nodePtr);
+  corsika::stack::SecondaryView view(particle);
+  auto projectile = view.GetProjectile();
 
-TEST_CASE("UrQMD") { UrQMD proc; }
+  [[maybe_unused]] const process::EProcessReturn ret = urqmd.DoInteraction(projectile);
+}