diff --git a/Processes/UrQMD/UrQMD.cc b/Processes/UrQMD/UrQMD.cc
index 293d0c478273c6ff1907fea27050a7b539981f97..2f48668b725623bfca612e5e0ad6d7bbeb9c68a1 100644
--- a/Processes/UrQMD/UrQMD.cc
+++ b/Processes/UrQMD/UrQMD.cc
@@ -50,50 +50,46 @@ corsika::process::EProcessReturn UrQMD::DoInteraction(SetupProjectile& projectil
auto const& projectilePosition = projectile.GetPosition();
auto const projectileTime = projectile.GetTime();
+ // sample target particle
+ auto const& mediumComposition =
+ projectile.GetNode()->GetModelProperties().GetNuclearComposition();
+ auto const componentCrossSections = std::invoke([&]() {
+ auto const& components = mediumComposition.GetComponents();
+ std::vector<CrossSectionType> crossSections;
+ crossSections.reserve(components.size());
+
+ for (auto const c : components) {
+ crossSections.push_back(GetCrossSection(projectileCode, c, projectileEnergyLab));
+ }
+
+ return crossSections;
+ });
+
+ auto const targetCode = mediumComposition.SampleTarget(componentCrossSections, fRNG);
+ auto const targetA = particles::GetNucleusA(targetCode);
+ auto const targetZ = particles::GetNucleusZ(targetCode);
+
inputs_.nevents = 1;
sys_.eos = 0; // could be configurable in principle
inputs_.outsteps = 1;
sys_.nsteps = 1;
- // todo: sample target
-
- int const Atarget = 14;
- int tableIndex = 0; // 0: nitrogen, 1: oxygen, 2: argon target
-
- /*
- // corsika 7
- bmin = 0.d0
- CTOption(5) = 1
- if ( iflbmax.eq.1 ) then
- bdist = BIM(LIT)
- else
- bdist=nucrad(Ap)+nucrad(At)+2*CTParam(30)
- endif
-
- // conex
- CTOption(5)=1
- if ( prspflg.eq.0 ) then
- bdist = BIM(LT)
- else
- bdist = xsbmax
- endif
- */
-
// initialization regarding projectile
if (particles::Code::Nucleus == projectileCode) {
// is this everything?
inputs_.prspflg = 0;
- rsys_.bdist = cxs_u2_.bim[tableIndex];
sys_.Ap = projectile.GetNuclearA();
sys_.Zp = projectile.GetNuclearZ();
+ rsys_.bdist = nucrad_(targetA) + nucrad_(sys_.Ap) + 2 * options_.CTParam[30 - 1];
+
int const id = 1;
cascinit_(sys_.Zp, sys_.Ap, id);
} else {
inputs_.prspflg = 1;
sys_.Ap = 1; // even for non-baryons this has to be set, see vanilla UrQMD.f
- rsys_.bdist = nucrad_(Atarget) + nucrad_(1) + 2 * options_.CTParam[30 - 1];
+ rsys_.bdist = nucrad_(targetA) + nucrad_(1) + 2 * options_.CTParam[30 - 1];
auto const [ityp, iso3] = ConvertToUrQMD(projectileCode);
// todo: conversion of K_long/short into strong eigenstates;
@@ -104,25 +100,9 @@ corsika::process::EProcessReturn UrQMD::DoInteraction(SetupProjectile& projectil
rsys_.ebeam = (projectileEnergyLab - projectile.GetMass()) * (1 / 1_GeV);
// initilazation regarding target
- auto const& mediumComposition =
- projectile.GetNode()->GetModelProperties().GetNuclearComposition();
- auto const componentCrossSections = std::invoke([&]() {
- auto const& components = mediumComposition.GetComponents();
- std::vector<CrossSectionType> crossSections;
- crossSections.reserve(components.size());
-
- for (auto const c : components) {
- crossSections.push_back(GetCrossSection(projectileCode, c, projectileEnergyLab));
- }
-
- return crossSections;
- });
-
- auto const targetCode = mediumComposition.SampleTarget(componentCrossSections, fRNG);
-
if (particles::IsNucleus(targetCode)) {
- sys_.Zt = particles::GetNucleusZ(targetCode);
- sys_.At = particles::GetNucleusA(targetCode);
+ sys_.Zt = targetZ;
+ sys_.At = targetA;
inputs_.trspflg = 0; // nucleus as target
int const id = 2;
cascinit_(sys_.Zt, sys_.At, id);
diff --git a/Processes/UrQMD/testUrQMD.cc b/Processes/UrQMD/testUrQMD.cc
index 2dad439de480f1ccb01408693addb3504fdef3db..c715f11808437b22a5e3fbbf346d14df9d912dbb 100644
--- a/Processes/UrQMD/testUrQMD.cc
+++ b/Processes/UrQMD/testUrQMD.cc
@@ -18,6 +18,8 @@
#include <corsika/units/PhysicalConstants.h>
#include <corsika/units/PhysicalUnits.h>
+#include <corsika/utl/CorsikaFenv.h>
+
#include <corsika/particles/ParticleProperties.h>
#include <corsika/setup/SetupStack.h>
#include <corsika/setup/SetupTrajectory.h>
@@ -34,7 +36,23 @@ using namespace corsika;
using namespace corsika::process::UrQMD;
using namespace corsika::units::si;
+template <typename TStack>
+auto sumCharge(TStack& stack) {
+ int totalCharge = 0;
+ int count = 0;
+ for (auto& p : stack) {
+ count++;
+ totalCharge += particles::GetChargeNumber(p.GetPID());
+ std::cout << p.GetPID() << " " << particles::GetChargeNumber(p.GetPID()) << std::endl;
+ }
+
+ std::cout << count << " particles on stack" << std::endl;
+
+ return totalCharge;
+}
+
TEST_CASE("UrQMD") {
+ feenableexcept(FE_INVALID);
corsika::random::RNGManager::GetInstance().RegisterRandomStream("UrQMD");
UrQMD urqmd;
@@ -63,21 +81,48 @@ TEST_CASE("UrQMD") {
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();
-
- [[maybe_unused]] const process::EProcessReturn ret = urqmd.DoInteraction(projectile);
+ const HEPEnergyType P0 = 1000_GeV;
+ auto pLab = corsika::stack::MomentumVector(cs, {P0, 0_GeV, 0_GeV});
+
+ SECTION("nucleon projectile") {
+ setup::Stack stack;
+
+ unsigned short constexpr A = 16, Z = 8;
+ auto constexpr mN = corsika::units::constants::nucleonMass;
+ HEPMomentumType E0 = sqrt(A * A * mN * mN + P0 * P0);
+ 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});
+
+ particle.SetNode(nodePtr);
+ corsika::stack::SecondaryView view(particle);
+ auto projectile = view.GetProjectile();
+
+ [[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(projectile);
+
+ REQUIRE(sumCharge(stack) == Z + particles::GetChargeNumber(particles::Code::Oxygen));
+ }
+
+ SECTION("\"special\" projectile") {
+
+ setup::Stack stack;
+
+ auto constexpr code = particles::Code::Neutron;
+ auto constexpr mass = particles::GetMass(code);
+ HEPMomentumType E0 = sqrt(mass * mass + pLab.squaredNorm());
+ auto particle = stack.AddParticle(
+ std::tuple<particles::Code, units::si::HEPEnergyType,
+ corsika::stack::MomentumVector, geometry::Point, units::si::TimeType>{
+ code, E0, pLab, origin, 0_ns});
+ particle.SetNode(nodePtr);
+ corsika::stack::SecondaryView view(particle);
+ auto projectile = view.GetProjectile();
+
+ [[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(projectile);
+
+ REQUIRE(sumCharge(stack) == particles::GetChargeNumber(particles::Code::PiPlus) +
+ particles::GetChargeNumber(particles::Code::Oxygen));
+ }
}