diff --git a/Processes/UrQMD/UrQMD.cc b/Processes/UrQMD/UrQMD.cc
index b60edba1b5fdc636723c70cbfcee76e92281af4f..fba5dc1833be9f813d14b6dfd9cb3c872192594c 100644
--- a/Processes/UrQMD/UrQMD.cc
+++ b/Processes/UrQMD/UrQMD.cc
@@ -28,29 +28,20 @@ using SetupStack = corsika::setup::Stack;
using SetupParticle = corsika::setup::Stack::StackIterator;
using SetupProjectile = corsika::setup::StackView::StackIterator;
-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 {
- using namespace units::si;
-
- // TODO: energy cuts, return 0 for non-hadrons
-
- auto const projectileCode = vProjectile.GetPID();
- auto const projectileEnergyLab = vProjectile.GetEnergy();
-
+CrossSectionType UrQMD::GetCrossSection(particles::Code vProjectileCode,
+ corsika::particles::Code vTargetCode,
+ HEPEnergyType vLabEnergy, int vAProjectile = 1) {
// the following is a translation of ptsigtot() into C++
- if (projectileCode != particles::Code::Nucleus &&
+ if (vProjectileCode != particles::Code::Nucleus &&
!IsNucleus(vTargetCode)) { // both particles are "special"
- auto const mProj = particles::GetMass(projectileCode);
+ auto const mProj = particles::GetMass(vProjectileCode);
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 * vLabEnergy * mTar) *
+ (1 / 1_GeV);
// we must set some UrQMD globals first...
- auto const [ityp, iso3] = ConvertToUrQMD(projectileCode);
+ auto const [ityp, iso3] = ConvertToUrQMD(vProjectileCode);
inputs_.spityp[0] = ityp;
inputs_.spiso3[0] = iso3;
@@ -62,8 +53,7 @@ CrossSectionType UrQMD::GetCrossSection(TParticle const& vProjectile,
int two = 2;
return sigtot_(one, two, sqrtS) * 1_mb;
} else {
- int const Ap =
- (projectileCode == particles::Code::Nucleus) ? vProjectile.GetNuclearA() : 1;
+ int const Ap = vAProjectile;
int const At = IsNucleus(vTargetCode) ? particles::GetNucleusA(vTargetCode) : 1;
double const maxImpact = nucrad_(Ap) + nucrad_(At) + 2 * options_.CTParam[30 - 1];
@@ -72,6 +62,26 @@ CrossSectionType UrQMD::GetCrossSection(TParticle const& vProjectile,
}
}
+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 {
+ // TODO: return 0 for non-hadrons?
+
+ auto const projectileCode = vProjectile.GetPID();
+ auto const projectileEnergyLab = vProjectile.GetEnergy();
+
+ if (projectileCode == particles::Code::K0Long) {
+ return 0.5 *
+ (GetCrossSection(particles::Code::K0, vTargetCode, projectileEnergyLab) +
+ GetCrossSection(particles::Code::K0Bar, vTargetCode, projectileEnergyLab));
+ }
+
+ int const Ap =
+ (projectileCode == particles::Code::Nucleus) ? vProjectile.GetNuclearA() : 1;
+ return GetCrossSection(projectileCode, vTargetCode, projectileEnergyLab, Ap);
+}
+
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
@@ -82,7 +92,7 @@ bool UrQMD::CanInteract(particles::Code vCode) const {
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};
+ particles::Code::K0, particles::Code::K0Bar, particles::Code::K0Long};
return std::find(std::cbegin(validProjectileCodes), std::cend(validProjectileCodes),
vCode) != std::cend(validProjectileCodes);
@@ -109,7 +119,7 @@ GrammageType UrQMD::GetInteractionLength(SetupParticle& vParticle) const {
corsika::process::EProcessReturn UrQMD::DoInteraction(SetupProjectile& vProjectile) {
using namespace units::si;
- auto const projectileCode = vProjectile.GetPID();
+ auto projectileCode = vProjectile.GetPID();
auto const projectileEnergyLab = vProjectile.GetEnergy();
auto const& projectileMomentumLab = vProjectile.GetMomentum();
auto const& projectilePosition = vProjectile.GetPosition();
@@ -159,6 +169,12 @@ corsika::process::EProcessReturn UrQMD::DoInteraction(SetupProjectile& vProjecti
rsys_.bdist = nucrad_(targetA) + nucrad_(1) + 2 * options_.CTParam[30 - 1];
rsys_.ebeam = (projectileEnergyLab - vProjectile.GetMass()) * (1 / 1_GeV);
+ if (projectileCode == particles::Code::K0Long) {
+ projectileCode = fBooleanDist(fRNG) ? particles::Code::K0 : particles::Code::K0Bar;
+ } else if (projectileCode == particles::Code::K0Short) {
+ throw std::runtime_error("K0Short should not interact");
+ }
+
auto const [ityp, iso3] = ConvertToUrQMD(projectileCode);
// todo: conversion of K_long/short into strong eigenstates;
inputs_.spityp[0] = ityp;
@@ -188,7 +204,11 @@ corsika::process::EProcessReturn UrQMD::DoInteraction(SetupProjectile& vProjecti
originalCS.RotateToZ(projectileMomentumLab);
for (int i = 0; i < sys_.npart; ++i) {
- auto const code = ConvertFromUrQMD(isys_.ityp[i], isys_.iso3[i]);
+ auto code = ConvertFromUrQMD(isys_.ityp[i], isys_.iso3[i]);
+ if (code == particles::Code::K0 || code == particles::Code::K0Bar) {
+ code = fBooleanDist(fRNG) ? particles::Code::K0Short : particles::Code::K0Long;
+ }
+
// "coor_.p0[i] * 1_GeV" is likely off-shell as UrQMD doesn't preserve masses well
auto momentum = geometry::Vector(
zAxisFrame,
@@ -245,6 +265,8 @@ std::pair<int, int> corsika::process::UrQMD::ConvertToUrQMD(
{-311, {-106, 1}}, // K0bar
{2212, {1, 1}}, // p
{2112, {1, -1}}, // n
+ {-2212, {-1, -1}}, // pbar
+ {-2112, {-1, 1}}, // nbar
{221, {102, 0}}, // eta
{213, {104, 2}}, // rho+
{-213, {104, -2}}, // rho-
diff --git a/Processes/UrQMD/UrQMD.h b/Processes/UrQMD/UrQMD.h
index 3be0168ac3aaf327a79a952869b0d60fec12e0a0..bdbca146ea015ce66174b7c2b2b99f619636c466 100644
--- a/Processes/UrQMD/UrQMD.h
+++ b/Processes/UrQMD/UrQMD.h
@@ -18,6 +18,7 @@
#include <corsika/units/PhysicalUnits.h>
#include <array>
+#include <random>
#include <utility>
namespace corsika::process::UrQMD {
@@ -38,8 +39,12 @@ namespace corsika::process::UrQMD {
bool CanInteract(particles::Code) const;
private:
+ static corsika::units::si::CrossSectionType GetCrossSection(
+ particles::Code, particles::Code, corsika::units::si::HEPEnergyType, int);
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("UrQMD");
+
+ std::uniform_int_distribution<int> fBooleanDist{0, 1};
};
namespace constants {
diff --git a/Processes/UrQMD/testUrQMD.cc b/Processes/UrQMD/testUrQMD.cc
index 3e347a827eaf5564836d5f82cabb55dbce471ffe..81b4215475cf2d927ed6a7e86cb20462798c5c3b 100644
--- a/Processes/UrQMD/testUrQMD.cc
+++ b/Processes/UrQMD/testUrQMD.cc
@@ -144,7 +144,7 @@ TEST_CASE("UrQMD") {
particles::Code validProjectileCodes[] = {
particles::Code::PiPlus, particles::Code::PiMinus, particles::Code::Proton,
particles::Code::Neutron, particles::Code::KPlus, particles::Code::KMinus,
- particles::Code::K0, particles::Code::K0Bar};
+ particles::Code::K0, particles::Code::K0Bar, particles::Code::K0Long};
for (auto code : validProjectileCodes) {
auto [stack, view] = setupStack(code, 100_GeV, nodePtr, cs);
@@ -205,4 +205,25 @@ TEST_CASE("UrQMD") {
REQUIRE((secMomSum - projectileMomentum).norm() / projectileMomentum.norm() ==
Approx(0).margin(1e-2));
}
+
+ SECTION("K0Long projectile") {
+ auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Oxygen);
+ auto [stackPtr, secViewPtr] =
+ setupStack(particles::Code::K0Long, 400_GeV, nodePtr, *csPtr);
+
+ // must be assigned to variable, cannot be used as rvalue?!
+ auto projectile = secViewPtr->GetProjectile();
+ auto const projectileMomentum = projectile.GetMomentum();
+
+ [[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(projectile);
+
+ REQUIRE(sumCharge(*secViewPtr) ==
+ particles::GetChargeNumber(particles::Code::K0Long) +
+ particles::GetChargeNumber(particles::Code::Oxygen));
+
+ auto const secMomSum =
+ sumMomentum(*secViewPtr, projectileMomentum.GetCoordinateSystem());
+ REQUIRE((secMomSum - projectileMomentum).norm() / projectileMomentum.norm() ==
+ Approx(0).margin(1e-2));
+ }
}