From 2febbb3e94bcdf3b7b4b8d5c7e4b7b3e2ef2d18f Mon Sep 17 00:00:00 2001
From: ralfulrich <ralf.ulrich@kit.edu>
Date: Mon, 25 Oct 2021 20:49:45 +0200
Subject: [PATCH] added urqmd
---
corsika/detail/framework/core/Cascade.inl | 8 +-
.../framework/process/ProcessSequence.inl | 27 +-
.../process/SwitchProcessSequence.inl | 22 +-
.../modules/urqmd/ParticleConversion.inl | 104 +++++++
corsika/detail/modules/urqmd/UrQMD.inl | 288 ++++++------------
corsika/framework/process/ProcessSequence.hpp | 4 +-
.../process/SwitchProcessSequence.hpp | 2 +-
corsika/modules/urqmd/ParticleConversion.hpp | 38 +++
corsika/modules/urqmd/UrQMD.hpp | 38 +--
tests/framework/testCascade.cpp | 8 +-
tests/modules/CMakeLists.txt | 2 +-
tests/modules/testSibyll.cpp | 3 +-
tests/modules/testUrQMD.cpp | 147 +++------
13 files changed, 331 insertions(+), 360 deletions(-)
create mode 100644 corsika/detail/modules/urqmd/ParticleConversion.inl
create mode 100644 corsika/modules/urqmd/ParticleConversion.hpp
diff --git a/corsika/detail/framework/core/Cascade.inl b/corsika/detail/framework/core/Cascade.inl
index 0e00520bb..7aabfed10 100644
--- a/corsika/detail/framework/core/Cascade.inl
+++ b/corsika/detail/framework/core/Cascade.inl
@@ -79,7 +79,6 @@ namespace corsika {
currentLogicalNode->getModelProperties().getNuclearComposition();
// determine sqrtS per nucleon pair, sqrtS_NN
- Code const projectileId = particle.getPID();
HEPEnergyType const Elab = particle.getEnergy();
FourMomentum const projectileP4{Elab, particle.getMomentum()};
CrossSectionType const sigma =
@@ -88,8 +87,7 @@ namespace corsika {
get_mass(targetId),
MomentumVector(particle.getMomentum().getCoordinateSystem(),
{0_GeV, 0_GeV, 0_GeV}));
- return sequence_.getCrossSection(projectileId, targetId, projectileP4,
- targetP4);
+ return sequence_.getCrossSection(particle, targetId, targetP4);
});
interaction(secondaries, projectileP4, composition, sigma);
sequence_.doSecondaries(secondaries);
@@ -113,7 +111,6 @@ namespace corsika {
currentLogicalNode->getModelProperties().getNuclearComposition();
// determine sqrtS per nucleon pair, sqrtS_NN
- Code const projectileId = particle.getPID();
HEPEnergyType const Elab = particle.getEnergy();
FourMomentum const projectileP4{Elab, particle.getMomentum()};
@@ -125,8 +122,7 @@ namespace corsika {
get_mass(targetId),
MomentumVector(particle.getMomentum().getCoordinateSystem(),
{0_GeV, 0_GeV, 0_GeV}));
- return sequence_.getCrossSection(projectileId, targetId, projectileP4,
- targetP4);
+ return sequence_.getCrossSection(particle, targetId, targetP4);
});
// calculate interaction length in medium
diff --git a/corsika/detail/framework/process/ProcessSequence.inl b/corsika/detail/framework/process/ProcessSequence.inl
index fb977f252..024c82809 100644
--- a/corsika/detail/framework/process/ProcessSequence.inl
+++ b/corsika/detail/framework/process/ProcessSequence.inl
@@ -310,30 +310,36 @@ namespace corsika {
template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
int IndexProcess2>
+ template <typename TParticle>
inline CrossSectionType
ProcessSequence<TProcess1, TProcess2, IndexStart, IndexProcess1,
- IndexProcess2>::getCrossSection(Code const projectileId,
+ IndexProcess2>::getCrossSection(TParticle const& projectile,
Code const targetId,
- FourMomentum const& projectileP4,
FourMomentum const& targetP4) const {
CrossSectionType tot = CrossSectionType::zero();
if constexpr (is_process_v<process1_type>) { // to protect from further compiler
// errors if process1_type is invalid
- if constexpr (is_interaction_process_v<process1_type> ||
- process1_type::is_process_sequence) {
- tot += A_.getCrossSection(projectileId, targetId, projectileP4, targetP4);
+ if constexpr (is_interaction_process_v<process1_type>) {
+ tot += A_.getCrossSection(projectile.getPID(), targetId,
+ {projectile.getEnergy(), projectile.getMomentum()},
+ targetP4);
+ } else if constexpr (process1_type::is_process_sequence) {
+ tot += A_.getCrossSection(projectile, targetId, targetP4);
}
}
if constexpr (is_process_v<process2_type>) { // to protect from further compiler
// errors if process2_type is invalid
- if constexpr (is_interaction_process_v<process2_type> ||
- process2_type::is_process_sequence) {
- tot += B_.getCrossSection(projectileId, targetId, projectileP4, targetP4);
+ if constexpr (is_interaction_process_v<process2_type>) {
+ tot += B_.getCrossSection(projectile.getPID(), targetId,
+ {projectile.getEnergy(), projectile.getMomentum()},
+ targetP4);
+ } else if constexpr (process2_type::is_process_sequence) {
+ tot += B_.getCrossSection(projectile, targetId, targetP4);
}
+ return tot;
}
- return tot;
}
template <typename TProcess1, typename TProcess2, int IndexStart, int IndexProcess1,
@@ -436,6 +442,7 @@ namespace corsika {
Code const projectileId = projectile.getPID();
// get cross section vector for all material components
+ // for selected process A
static_assert(
has_method_getCrossSection_v<TProcess1, // process object
CrossSectionType, // return type
@@ -500,7 +507,7 @@ namespace corsika {
auto const& projectile = view.parent();
Code const projectileId = projectile.getPID();
- // get cross section vector for all material components
+ // get cross section vector for all material components, for selected process B
static_assert(has_method_getCrossSection_v<TProcess2, // process object
CrossSectionType, // return type
Code, Code, FourMomentum const&,
diff --git a/corsika/detail/framework/process/SwitchProcessSequence.inl b/corsika/detail/framework/process/SwitchProcessSequence.inl
index b7cf27db0..02e197ec0 100644
--- a/corsika/detail/framework/process/SwitchProcessSequence.inl
+++ b/corsika/detail/framework/process/SwitchProcessSequence.inl
@@ -213,26 +213,24 @@ namespace corsika {
CrossSectionType SwitchProcessSequence<
TCondition, TSequence, USequence, IndexStart, IndexProcess1,
IndexProcess2>::getCrossSection(TParticle const& projectile, Code const targetId,
- HEPEnergyType const sqrtSnn) const {
+ FourMomentum const& targetP4) const {
- if (select_(projectile.parent())) {
+ if (select_(projectile)) {
if constexpr (is_interaction_process_v<process1_type>) {
- return A_.getCrossSection(projectile.getPID(), targetId, sqrtSnn,
- projectile.getNuclearA(),
- is_nucleus(targetId) ? get_nucleus_A(targetId) : 0);
+ return A_.getCrossSection(projectile.getPID(), targetId,
+ {projectile.getEnergy(), projectile.getMomentum()},
+ targetP4);
} else if (process1_type::is_process_sequence) {
- return A_.getCrossSection(projectile, targetId, sqrtSnn,
- is_nucleus(targetId) ? get_nucleus_A(targetId) : 0);
+ return A_.getCrossSection(projectile, targetId, targetP4);
}
} else {
if constexpr (is_interaction_process_v<process2_type>) {
- return B_.getCrossSection(projectile.getPID(), targetId, sqrtSnn,
- projectile.getNuclearA(),
- is_nucleus(targetId) ? get_nucleus_A(targetId) : 0);
+ return B_.getCrossSection(projectile.getPID(), targetId,
+ {projectile.getEnergy(), projectile.getMomentum()},
+ targetP4);
} else if (process2_type::is_process_sequence) {
- return B_.getCrossSection(projectile, targetId, sqrtSnn,
- is_nucleus(targetId) ? get_nucleus_A(targetId) : 0);
+ return B_.getCrossSection(projectile, targetId, targetP4);
}
}
return CrossSectionType::zero(); // default value
diff --git a/corsika/detail/modules/urqmd/ParticleConversion.inl b/corsika/detail/modules/urqmd/ParticleConversion.inl
new file mode 100644
index 000000000..de533d048
--- /dev/null
+++ b/corsika/detail/modules/urqmd/ParticleConversion.inl
@@ -0,0 +1,104 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/modules/urqmd/ParticleConversion.hpp>
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+#include <corsika/framework/core/PhysicalUnits.hpp>
+
+#include <urqmd.hpp>
+
+namespace corsika::urqmd {
+
+ inline bool canInteract(Code const vCode) {
+ // 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 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);
+ }
+
+ inline std::pair<int, int> convertToUrQMD(Code const code) {
+ static const std::map<int, std::pair<int, int>> mapPDGToUrQMD{
+ // data mostly from github.com/afedynitch/ParticleDataTool
+ {22, {100, 0}}, // photon
+ {111, {101, 0}}, // pi0
+ {211, {101, 2}}, // pi+
+ {-211, {101, -2}}, // pi-
+ {321, {106, 1}}, // K+
+ {-321, {-106, -1}}, // K-
+ {311, {106, -1}}, // K0
+ {-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-
+ {113, {104, 0}}, // rho0
+ {323, {108, 2}}, // K*+
+ {-323, {108, -2}}, // K*-
+ {313, {108, 0}}, // K*0
+ {-313, {-108, 0}}, // K*0-bar
+ {223, {103, 0}}, // omega
+ {333, {109, 0}}, // phi
+ {3222, {40, 2}}, // Sigma+
+ {3212, {40, 0}}, // Sigma0
+ {3112, {40, -2}}, // Sigma-
+ {3322, {49, 0}}, // Xi0
+ {3312, {49, -1}}, // Xi-
+ {3122, {27, 0}}, // Lambda0
+ {2224, {17, 4}}, // Delta++
+ {2214, {17, 2}}, // Delta+
+ {2114, {17, 0}}, // Delta0
+ {1114, {17, -2}}, // Delta-
+ {3224, {41, 2}}, // Sigma*+
+ {3214, {41, 0}}, // Sigma*0
+ {3114, {41, -2}}, // Sigma*-
+ {3324, {50, 0}}, // Xi*0
+ {3314, {50, -1}}, // Xi*-
+ {3334, {55, 0}}, // Omega-
+ {411, {133, 2}}, // D+
+ {-411, {133, -2}}, // D-
+ {421, {133, 0}}, // D0
+ {-421, {-133, 0}}, // D0-bar
+ {441, {107, 0}}, // etaC
+ {431, {138, 1}}, // Ds+
+ {-431, {138, -1}}, // Ds-
+ {433, {139, 1}}, // Ds*+
+ {-433, {139, -1}}, // Ds*-
+ {413, {134, 1}}, // D*+
+ {-413, {134, -1}}, // D*-
+ {10421, {134, 0}}, // D*0
+ {-10421, {-134, 0}}, // D*0-bar
+ {443, {135, 0}}, // jpsi
+ };
+
+ return mapPDGToUrQMD.at(static_cast<int>(get_PDG(code)));
+ }
+
+ inline Code convertFromUrQMD(int vItyp, int vIso3) {
+ int const pdgInt =
+ ::urqmd::pdgid_(vItyp, vIso3); // use the conversion function provided by UrQMD
+ if (pdgInt == 0) { // ::urqmd::pdgid_ returns 0 on error
+ throw std::runtime_error("UrQMD pdgid() returned 0");
+ }
+ auto const pdg = static_cast<PDGCode>(pdgInt);
+ return convert_from_PDG(pdg);
+ }
+
+} // namespace corsika::urqmd
diff --git a/corsika/detail/modules/urqmd/UrQMD.inl b/corsika/detail/modules/urqmd/UrQMD.inl
index 80da719ed..d5c2bff04 100644
--- a/corsika/detail/modules/urqmd/UrQMD.inl
+++ b/corsika/detail/modules/urqmd/UrQMD.inl
@@ -9,11 +9,13 @@
#pragma once
#include <corsika/modules/urqmd/UrQMD.hpp>
+#include <corsika/modules/urqmd/ParticleConversion.hpp>
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
#include <corsika/framework/geometry/QuantityVector.hpp>
#include <corsika/framework/geometry/Vector.hpp>
+#include <corsika/framework/utility/COMBoost.hpp>
#include <boost/filesystem.hpp>
#include <boost/multi_array.hpp>
@@ -35,12 +37,22 @@ namespace corsika::urqmd {
::urqmd::iniurqmdc8_();
}
- inline CrossSectionType UrQMD::getTabulatedCrossSection(Code projectileCode,
- Code targetCode,
- HEPEnergyType labEnergy) const {
+ inline void UrQMD::isValid(Code const projectileId, Code const targetId) const {
+
+ if (!is_hadron(projectileId) || !corsika::urqmd::canInteract(projectileId)) {
+ throw std::runtime_error("UrQMD projectile is not a compatible hadron.");
+ }
+ if (!is_nucleus(targetId)) {
+ throw std::runtime_error("UrQMD target is not a nucleus.");
+ }
+ }
+
+ inline CrossSectionType UrQMD::getTabulatedCrossSection(
+ Code const projectileId, Code const targetId, HEPEnergyType const labEnergy) const {
+
// translated to C++ from CORSIKA 7 subroutine cxtot_u
- auto const kinEnergy = labEnergy - get_mass(projectileCode);
+ auto const kinEnergy = labEnergy - get_mass(projectileId);
assert(kinEnergy >= HEPEnergyType::zero());
@@ -54,7 +66,7 @@ namespace corsika::urqmd {
w[2 - 1] = w[2 - 1] - 2 * w[3 - 1];
int projectileIndex;
- switch (projectileCode) {
+ switch (projectileId) {
case Code::Proton:
projectileIndex = 0;
break;
@@ -88,14 +100,14 @@ namespace corsika::urqmd {
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);
+ CORSIKA_LOG_WARN("UrQMD cross-section not tabulated for {}", projectileId);
return CrossSectionType::zero();
// LCOV_EXCL_STOP
}
}
int targetIndex;
- switch (targetCode) {
+ switch (targetId) {
case Code::Nitrogen:
targetIndex = 0;
break;
@@ -107,7 +119,7 @@ namespace corsika::urqmd {
break;
default:
std::stringstream ss;
- ss << "UrQMD cross-section not tabluated for target " << targetCode;
+ ss << "UrQMD cross-section not tabluated for target " << targetId;
throw std::runtime_error(ss.str().data());
}
@@ -119,53 +131,17 @@ namespace corsika::urqmd {
CORSIKA_LOG_TRACE(
"UrQMD::GetTabulatedCrossSection proj={}, targ={}, E={} GeV, sigma={}",
- get_name(projectileCode), get_name(targetCode), labEnergy / 1_GeV, result);
+ get_name(projectileId), get_name(targetId), labEnergy / 1_GeV, result);
return result;
}
- inline CrossSectionType UrQMD::getCrossSection(Code vProjectileCode, Code vTargetCode,
- HEPEnergyType vLabEnergy,
- int vAProjectile = 1) {
-
- // the following is a translation of ptsigtot() into C++
- if (!is_nucleus(vProjectileCode) &&
- !is_nucleus(vTargetCode)) { // both particles are "special"
- auto const mProj = get_mass(vProjectileCode);
- auto const mTar = get_mass(vTargetCode);
- double sqrtS =
- sqrt(static_pow<2>(mProj) + static_pow<2>(mTar) + 2 * vLabEnergy * mTar) *
- (1 / 1_GeV);
-
- // we must set some UrQMD globals first...
- auto const [ityp, iso3] = convertToUrQMD(vProjectileCode);
- ::urqmd::inputs_.spityp[0] = ityp;
- ::urqmd::inputs_.spiso3[0] = iso3;
-
- auto const [itypTar, iso3Tar] = convertToUrQMD(vTargetCode);
- ::urqmd::inputs_.spityp[1] = itypTar;
- ::urqmd::inputs_.spiso3[1] = iso3Tar;
-
- int one = 1;
- int two = 2;
- return ::urqmd::sigtot_(one, two, sqrtS) * 1_mb;
- } else {
- int const Ap = vAProjectile;
- int const At = is_nucleus(vTargetCode) ? get_nucleus_A(vTargetCode) : 1;
-
- double const maxImpact = ::urqmd::nucrad_(Ap) + ::urqmd::nucrad_(At) +
- 2 * ::urqmd::options_.CTParam[30 - 1];
- return 10_mb * M_PI * static_pow<2>(maxImpact);
- // is a constant cross-section really reasonable?
- }
- }
-
- template <typename TParticle>
- inline CrossSectionType UrQMD::getCrossSection(TParticle const& projectile,
- Code targetCode) const {
- auto const projectileCode = projectile.getPID();
+ inline CrossSectionType UrQMD::getCrossSection(Code const projectileId,
+ Code const targetId,
+ FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) const {
- if (is_nucleus(projectileCode)) {
+ if (is_nucleus(projectileId)) {
/*
* unfortunately unavoidable at the moment until we have tools to get the actual
* inealstic cross-section from UrQMD
@@ -173,72 +149,61 @@ namespace corsika::urqmd {
return CrossSectionType::zero();
}
- return getTabulatedCrossSection(projectileCode, targetCode, projectile.getEnergy());
- }
+ isValid(projectileId, targetId); // throws
- inline bool UrQMD::canInteract(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
+ // define projectile, in lab frame
+ auto const sqrtS2 = (projectileP4 + targetP4).getNormSqr();
+ HEPEnergyType const Elab = (sqrtS2 - static_pow<2>(get_mass(projectileId)) -
+ static_pow<2>(get_mass(targetId))) /
+ (2 * get_mass(targetId));
- 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};
+ bool const tabulated = true;
+ if (tabulated) { return getTabulatedCrossSection(projectileId, targetId, Elab); }
- return std::find(std::cbegin(validProjectileCodes), std::cend(validProjectileCodes),
- vCode) != std::cend(validProjectileCodes);
- }
+ // the following is a translation of ptsigtot() into C++
+ if (!is_nucleus(projectileId) &&
+ !is_nucleus(targetId)) { // both particles are "special"
- template <typename TParticle>
- inline GrammageType UrQMD::getInteractionLength(TParticle const& vParticle) const {
+ double sqrtS = sqrt(sqrtS2) / 1_GeV;
- 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);
- }
+ // we must set some UrQMD globals first...
+ auto const [ityp, iso3] = corsika::urqmd::convertToUrQMD(projectileId);
+ ::urqmd::inputs_.spityp[0] = ityp;
+ ::urqmd::inputs_.spiso3[0] = iso3;
- auto const& mediumComposition =
- vParticle.getNode()->getModelProperties().getNuclearComposition();
- using namespace std::placeholders;
+ auto const [itypTar, iso3Tar] = corsika::urqmd::convertToUrQMD(targetId);
+ ::urqmd::inputs_.spityp[1] = itypTar;
+ ::urqmd::inputs_.spiso3[1] = iso3Tar;
- CrossSectionType const weightedProdCrossSection = mediumComposition.getWeightedSum(
- std::bind(&UrQMD::getCrossSection<decltype(vParticle)>, this, vParticle, _1));
+ int one = 1;
+ int two = 2;
+ return ::urqmd::sigtot_(one, two, sqrtS) * 1_mb;
+ }
- return mediumComposition.getAverageMassNumber() * constants::u /
- weightedProdCrossSection;
+ // at least one of them is a nucleus
+ int const Ap = is_nucleus(projectileId) ? get_nucleus_A(projectileId) : 1;
+ int const At = is_nucleus(targetId) ? get_nucleus_A(targetId) : 1;
+ double const maxImpact = ::urqmd::nucrad_(Ap) + ::urqmd::nucrad_(At) +
+ 2 * ::urqmd::options_.CTParam[30 - 1];
+ return 10_mb * M_PI * static_pow<2>(maxImpact);
+ // is a constant cross-section really reasonable?
}
template <typename TView>
- inline void UrQMD::doInteraction(TView& view) {
-
- auto projectile = view.getProjectile();
-
- Code projectileCode = projectile.getPID();
- auto const projectileEnergyLab = projectile.getEnergy();
- auto const& projectileMomentumLab = projectile.getMomentum();
- 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());
+ inline void UrQMD::doInteraction(TView& view, Code const projectileId,
+ Code const targetId, FourMomentum const& projectileP4,
+ FourMomentum const& targetP4) {
- for (auto const c : components) {
- crossSections.push_back(getCrossSection(projectile, c));
- }
+ // define projectile, in lab frame
+ auto const sqrtS2 = (projectileP4 + targetP4).getNormSqr();
+ HEPEnergyType const Elab = (sqrtS2 - static_pow<2>(get_mass(projectileId)) -
+ static_pow<2>(get_mass(targetId))) /
+ (2 * get_mass(targetId));
- return crossSections;
- });
+ isValid(projectileId, targetId); // throws
- auto const targetCode = mediumComposition.sampleTarget(componentCrossSections, RNG_);
- auto const targetA = get_nucleus_A(targetCode);
- auto const targetZ = get_nucleus_Z(targetCode);
+ size_t const targetA = get_nucleus_A(targetId);
+ size_t const targetZ = get_nucleus_Z(targetId);
::urqmd::inputs_.nevents = 1;
::urqmd::sys_.eos = 0; // could be configurable in principle
@@ -246,14 +211,13 @@ namespace corsika::urqmd {
::urqmd::sys_.nsteps = 1;
// initialization regarding projectile
- if (is_nucleus(projectileCode)) {
+ if (is_nucleus(projectileId)) {
// is this everything?
::urqmd::inputs_.prspflg = 0;
- ::urqmd::sys_.Ap = get_nucleus_A(projectileCode);
- ::urqmd::sys_.Zp = get_nucleus_Z(projectileCode);
- ::urqmd::rsys_.ebeam =
- (projectileEnergyLab - projectile.getMass()) * (1 / 1_GeV) / ::urqmd::sys_.Ap;
+ ::urqmd::sys_.Ap = get_nucleus_A(projectileId);
+ ::urqmd::sys_.Zp = get_nucleus_Z(projectileId);
+ ::urqmd::rsys_.ebeam = (Elab - get_mass(projectileId)) / 1_GeV / ::urqmd::sys_.Ap;
::urqmd::rsys_.bdist = ::urqmd::nucrad_(targetA) +
::urqmd::nucrad_(::urqmd::sys_.Ap) +
@@ -267,20 +231,19 @@ namespace corsika::urqmd {
1; // even for non-baryons this has to be set, see vanilla UrQMD.f
::urqmd::rsys_.bdist = ::urqmd::nucrad_(targetA) + ::urqmd::nucrad_(1) +
2 * ::urqmd::options_.CTParam[30 - 1];
- ::urqmd::rsys_.ebeam = (projectileEnergyLab - projectile.getMass()) * (1 / 1_GeV);
+ ::urqmd::rsys_.ebeam = (Elab - get_mass(projectileId)) / 1_GeV;
- if (projectileCode == Code::K0Long || projectileCode == Code::K0Short) {
- projectileCode = booleanDist_(RNG_) ? Code::K0 : Code::K0Bar;
- }
-
- auto const [ityp, iso3] = convertToUrQMD(projectileCode);
+ auto const [ityp, iso3] = corsika::urqmd::convertToUrQMD(
+ (projectileId == Code::K0Long || projectileId == Code::K0Short)
+ ? (booleanDist_(RNG_) ? Code::K0 : Code::K0Bar)
+ : projectileId);
// todo: conversion of K_long/short into strong eigenstates;
::urqmd::inputs_.spityp[0] = ityp;
::urqmd::inputs_.spiso3[0] = iso3;
}
- // initilazation regarding target
- if (is_nucleus(targetCode)) {
+ // initialization regarding target
+ if (is_nucleus(targetId)) {
::urqmd::sys_.Zt = targetZ;
::urqmd::sys_.At = targetA;
::urqmd::inputs_.trspflg = 0; // nucleus as target
@@ -288,7 +251,7 @@ namespace corsika::urqmd {
::urqmd::cascinit_(::urqmd::sys_.Zt, ::urqmd::sys_.At, id);
} else {
::urqmd::inputs_.trspflg = 1; // special particle as target
- auto const [ityp, iso3] = convertToUrQMD(targetCode);
+ auto const [ityp, iso3] = corsika::urqmd::convertToUrQMD(targetId);
::urqmd::inputs_.spityp[1] = ityp;
::urqmd::inputs_.spiso3[1] = iso3;
}
@@ -297,103 +260,36 @@ namespace corsika::urqmd {
iflb_; // flag for retrying interaction in case of empty event, 0 means retry
::urqmd::urqmd_(iflb);
+ auto projectile = view.getProjectile();
+ auto const& projectilePosition = projectile.getPosition();
+ auto const projectileTime = projectile.getTime();
+
// now retrieve secondaries from UrQMD
- auto const& originalCS = projectileMomentumLab.getCoordinateSystem();
- CoordinateSystemPtr const& zAxisFrame =
- make_rotationToZ(originalCS, projectileMomentumLab);
+ COMBoost const boost(projectileP4, targetP4);
+ auto const& originalCS = boost.getOriginalCS();
+ auto const& csPrime = boost.getRotatedCS();
for (int i = 0; i < ::urqmd::sys_.npart; ++i) {
- auto code = convertFromUrQMD(::urqmd::isys_.ityp[i], ::urqmd::isys_.iso3[i]);
+ auto code = corsika::urqmd::convertFromUrQMD(::urqmd::isys_.ityp[i],
+ ::urqmd::isys_.iso3[i]);
if (code == Code::K0 || code == Code::K0Bar) {
code = booleanDist_(RNG_) ? Code::K0Short : Code::K0Long;
}
// "coor_.p0[i] * 1_GeV" is likely off-shell as UrQMD doesn't preserve masses well
- auto momentum =
- Vector(zAxisFrame,
- QuantityVector<dimensionless_d>{
- ::urqmd::coor_.px[i], ::urqmd::coor_.py[i], ::urqmd::coor_.pz[i]} *
- 1_GeV);
+ MomentumVector momentum{csPrime,
+ {::urqmd::coor_.px[i] * 1_GeV, ::urqmd::coor_.py[i] * 1_GeV,
+ ::urqmd::coor_.pz[i] * 1_GeV}};
momentum.rebase(originalCS); // transform back into standard lab frame
CORSIKA_LOG_DEBUG(" {} {} {} ", i, code, momentum.getComponents());
- projectile.addSecondary(
+ view.addSecondary(
std::make_tuple(code, momentum, projectilePosition, projectileTime));
}
CORSIKA_LOG_DEBUG("UrQMD generated {} secondaries!", ::urqmd::sys_.npart);
}
- inline Code convertFromUrQMD(int vItyp, int vIso3) {
- int const pdgInt =
- ::urqmd::pdgid_(vItyp, vIso3); // use the conversion function provided by UrQMD
- if (pdgInt == 0) { // ::urqmd::pdgid_ returns 0 on error
- throw std::runtime_error("UrQMD pdgid() returned 0");
- }
- auto const pdg = static_cast<PDGCode>(pdgInt);
- return convert_from_PDG(pdg);
- }
-
- inline std::pair<int, int> convertToUrQMD(Code code) {
- static const std::map<int, std::pair<int, int>> mapPDGToUrQMD{
- // data mostly from github.com/afedynitch/ParticleDataTool
- {22, {100, 0}}, // photon
- {111, {101, 0}}, // pi0
- {211, {101, 2}}, // pi+
- {-211, {101, -2}}, // pi-
- {321, {106, 1}}, // K+
- {-321, {-106, -1}}, // K-
- {311, {106, -1}}, // K0
- {-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-
- {113, {104, 0}}, // rho0
- {323, {108, 2}}, // K*+
- {-323, {108, -2}}, // K*-
- {313, {108, 0}}, // K*0
- {-313, {-108, 0}}, // K*0-bar
- {223, {103, 0}}, // omega
- {333, {109, 0}}, // phi
- {3222, {40, 2}}, // Sigma+
- {3212, {40, 0}}, // Sigma0
- {3112, {40, -2}}, // Sigma-
- {3322, {49, 0}}, // Xi0
- {3312, {49, -1}}, // Xi-
- {3122, {27, 0}}, // Lambda0
- {2224, {17, 4}}, // Delta++
- {2214, {17, 2}}, // Delta+
- {2114, {17, 0}}, // Delta0
- {1114, {17, -2}}, // Delta-
- {3224, {41, 2}}, // Sigma*+
- {3214, {41, 0}}, // Sigma*0
- {3114, {41, -2}}, // Sigma*-
- {3324, {50, 0}}, // Xi*0
- {3314, {50, -1}}, // Xi*-
- {3334, {55, 0}}, // Omega-
- {411, {133, 2}}, // D+
- {-411, {133, -2}}, // D-
- {421, {133, 0}}, // D0
- {-421, {-133, 0}}, // D0-bar
- {441, {107, 0}}, // etaC
- {431, {138, 1}}, // Ds+
- {-431, {138, -1}}, // Ds-
- {433, {139, 1}}, // Ds*+
- {-433, {139, -1}}, // Ds*-
- {413, {134, 1}}, // D*+
- {-413, {134, -1}}, // D*-
- {10421, {134, 0}}, // D*0
- {-10421, {-134, 0}}, // D*0-bar
- {443, {135, 0}}, // jpsi
- };
-
- 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);
diff --git a/corsika/framework/process/ProcessSequence.hpp b/corsika/framework/process/ProcessSequence.hpp
index 4ecf3d0cd..06111acf8 100644
--- a/corsika/framework/process/ProcessSequence.hpp
+++ b/corsika/framework/process/ProcessSequence.hpp
@@ -243,8 +243,8 @@ namespace corsika {
template <typename TParticle, typename TTrack>
ContinuousProcessStepLength getMaxStepLength(TParticle&& particle, TTrack&& vTrack);
- CrossSectionType getCrossSection(Code const projectileId, Code const targetId,
- FourMomentum const& projectilP4,
+ template <typename TParticle>
+ CrossSectionType getCrossSection(TParticle const& projectile, Code const targetId,
FourMomentum const& targetP4) const;
template <typename TParticle>
diff --git a/corsika/framework/process/SwitchProcessSequence.hpp b/corsika/framework/process/SwitchProcessSequence.hpp
index 74f07fff1..2b6190369 100644
--- a/corsika/framework/process/SwitchProcessSequence.hpp
+++ b/corsika/framework/process/SwitchProcessSequence.hpp
@@ -147,7 +147,7 @@ namespace corsika {
template <typename TParticle>
CrossSectionType getCrossSection(TParticle const& projectile, Code const targetId,
- HEPEnergyType const sqrtSnn) const;
+ FourMomentum const& targetP4) const;
template <typename TSecondaryView, typename TRNG>
ProcessReturn selectInteraction(TSecondaryView& view,
diff --git a/corsika/modules/urqmd/ParticleConversion.hpp b/corsika/modules/urqmd/ParticleConversion.hpp
new file mode 100644
index 000000000..85ae95d0d
--- /dev/null
+++ b/corsika/modules/urqmd/ParticleConversion.hpp
@@ -0,0 +1,38 @@
+/*
+ * (c) Copyright 2020 CORSIKA Project, corsika-project@lists.kit.edu
+ *
+ * This software is distributed under the terms of the GNU General Public
+ * Licence version 3 (GPL Version 3). See file LICENSE for a full version of
+ * the license.
+ */
+
+#pragma once
+
+#include <corsika/framework/core/ParticleProperties.hpp>
+
+#include <array>
+#include <utility>
+#include <string>
+
+namespace corsika::urqmd {
+
+ /**
+ * Checks if particle with Code can interaction in UrQMD.
+ */
+ bool canInteract(Code const);
+
+ /**
+ * convert CORSIKA code to UrQMD code tuple.
+ *
+ * In the current implementation a detour via the PDG code is made.
+ */
+ std::pair<int, int> convertToUrQMD(Code const);
+
+ /**
+ * convert UrQMD code to CORSIKA code.
+ */
+ Code convertFromUrQMD(int vItyp, int vIso3);
+
+} // namespace corsika::urqmd
+
+#include <corsika/detail/modules/urqmd/ParticleConversion.inl>
diff --git a/corsika/modules/urqmd/UrQMD.hpp b/corsika/modules/urqmd/UrQMD.hpp
index 2f3c27804..86d238e0f 100644
--- a/corsika/modules/urqmd/UrQMD.hpp
+++ b/corsika/modules/urqmd/UrQMD.hpp
@@ -10,9 +10,9 @@
#include <corsika/framework/core/ParticleProperties.hpp>
#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/framework/geometry/FourVector.hpp>
#include <boost/filesystem/path.hpp>
#include <boost/multi_array.hpp>
@@ -23,32 +23,30 @@
namespace corsika::urqmd {
- class UrQMD : public InteractionProcess<UrQMD> {
+ class UrQMD {
public:
/**
- * @param path Location of UrQMD XS data file
+ * The UrQMD interaction model.
+ *
+ * @param path Location of UrQMD XS data file.
* @param retryFlag Internal UrQMD flag for retrying interaction in case of empty
- * event, 0 means retry
+ * event, 0 means retry.
*/
UrQMD(boost::filesystem::path const path = corsika_data("UrQMD/UrQMD-1.3.1-xs.dat"),
int const retryFlag = 0);
- template <typename TParticle>
- GrammageType getInteractionLength(TParticle const&) const;
+ void isValid(Code const projectileId, Code const targetId) const;
- CrossSectionType getTabulatedCrossSection(Code, Code, HEPEnergyType) const;
+ CrossSectionType getTabulatedCrossSection(Code const, Code const,
+ HEPEnergyType const) const;
- template <typename TParticle>
- CrossSectionType getCrossSection(TParticle const&, Code) const;
+ CrossSectionType getCrossSection(Code const projectileId, Code const targetId,
+ FourMomentum const& projP4,
+ FourMomentum const& targP4) const;
template <typename TView>
- void doInteraction(TView&);
-
- bool canInteract(Code) const;
-
- void blob(int) {}
-
- static CrossSectionType getCrossSection(Code, Code, HEPEnergyType, int);
+ void doInteraction(TView&, Code const projectile, Code const targetId,
+ FourMomentum const& projP4, FourMomentum const& targP4);
private:
void readXSFile(boost::filesystem::path);
@@ -60,14 +58,6 @@ namespace corsika::urqmd {
boost::multi_array<CrossSectionType, 3> xs_interp_support_table_;
};
- /**
- * convert CORSIKA code to UrQMD code tuple
- *
- * In the current implementation a detour via the PDG code is made.
- */
- std::pair<int, int> convertToUrQMD(Code);
- Code convertFromUrQMD(int vItyp, int vIso3);
-
} // namespace corsika::urqmd
#include <corsika/detail/modules/urqmd/UrQMD.inl>
diff --git a/tests/framework/testCascade.cpp b/tests/framework/testCascade.cpp
index f2ac504e3..18732cb1f 100644
--- a/tests/framework/testCascade.cpp
+++ b/tests/framework/testCascade.cpp
@@ -113,8 +113,8 @@ private:
class ProcessCut : public SecondariesProcess<ProcessCut> {
public:
- ProcessCut(HEPEnergyType e)
- : fEcrit(e) {}
+ ProcessCut(HEPEnergyType const e)
+ : Ecrit_(e) {}
template <typename TStack>
void doSecondaries(TStack& vS) {
@@ -122,7 +122,7 @@ public:
auto p = vS.begin();
while (p != vS.end()) {
HEPEnergyType E = p.getEnergy();
- if (E < fEcrit) {
+ if (E < Ecrit_) {
p.erase();
count_++;
}
@@ -138,7 +138,7 @@ public:
private:
int count_ = 0;
int calls_ = 0;
- HEPEnergyType fEcrit;
+ HEPEnergyType Ecrit_;
};
TEST_CASE("Cascade", "[Cascade]") {
diff --git a/tests/modules/CMakeLists.txt b/tests/modules/CMakeLists.txt
index 5d6ec770f..1f36ee5f4 100644
--- a/tests/modules/CMakeLists.txt
+++ b/tests/modules/CMakeLists.txt
@@ -6,7 +6,7 @@ set (test_modules_sources
testObservationPlane.cpp
testQGSJetII.cpp
#testPythia8.cpp
- #testUrQMD.cpp
+ testUrQMD.cpp
#testCONEX.cpp
## testOnShellCheck.cpp
testParticleCut.cpp
diff --git a/tests/modules/testSibyll.cpp b/tests/modules/testSibyll.cpp
index 718a43601..94fb90792 100644
--- a/tests/modules/testSibyll.cpp
+++ b/tests/modules/testSibyll.cpp
@@ -29,7 +29,6 @@
#include <corsika/modules/sibyll/Random.hpp>
// TODODODODODODODODOODOD THIS MUST BE REMOVED
-#include <corsika/modules/urqmd/Random.hpp>
#include <corsika/modules/epos/Random.hpp>
// TODODODODODODODODOODOD THIS MUST BE REMOVED
@@ -279,7 +278,7 @@ TEST_CASE("SibyllInterface", "modules") {
// CHECK(view.getSize() == 11);
CHECK(cx / 1_mb == Approx(1100).margin(100)); // this is not physics validation
// CHECK(view.getSize() == 20); // also sibyll not stable wrt. to compiler changes
- CHECK(view.getSize() == Approx(200).margin(90)); // this is not physics validation
+ CHECK(view.getSize() == Approx(90).margin(10)); // this is not physics validation
}
}
diff --git a/tests/modules/testUrQMD.cpp b/tests/modules/testUrQMD.cpp
index 89bde00d2..0cabc8b72 100644
--- a/tests/modules/testUrQMD.cpp
+++ b/tests/modules/testUrQMD.cpp
@@ -11,8 +11,8 @@
#include <corsika/framework/core/ParticleProperties.hpp>
#include <corsika/framework/core/PhysicalConstants.hpp>
#include <corsika/framework/core/PhysicalUnits.hpp>
-#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/geometry/RootCoordinateSystem.hpp>
+#include <corsika/framework/geometry/Point.hpp>
#include <corsika/framework/geometry/Vector.hpp>
#include <corsika/framework/random/RNGManager.hpp>
#include <corsika/framework/utility/CorsikaFenv.hpp>
@@ -71,73 +71,46 @@ TEST_CASE("UrQMD") {
RNGManager<>::getInstance().registerRandomStream("urqmd");
UrQMD urqmd;
- SECTION("interaction length") {
- auto [env, csPtr, nodePtr] = setup::testing::setup_environment(Code::Nitrogen);
- auto const& cs = *csPtr;
- { [[maybe_unused]] auto const& env_dummy = env; }
-
- 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, 100_GeV, nodePtr, cs);
- CHECK(stack->getEntries() == 1);
- CHECK(view->getEntries() == 0);
-
- // simple check whether the cross-section is non-vanishing
- // 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, 100_GeV, nodePtr, cs);
- [[maybe_unused]] setup::StackView& viewRef = *(view.get());
- 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, 100_GeV, nodePtr, cs);
- [[maybe_unused]] setup::StackView& viewRef = *(view.get());
- CHECK_THROWS(urqmd.getInteractionLength(stack->getNextParticle()));
+ auto const rootCS = get_root_CoordinateSystem();
+
+ SECTION("valid") {
+ // this is how it is currently done
+ CHECK_THROWS(urqmd.isValid(Code::K0, Code::Proton));
+ CHECK_THROWS(urqmd.isValid(Code::DPlus, Code::Proton));
+ CHECK_THROWS(urqmd.isValid(Code::Electron, Code::Proton));
+ CHECK_THROWS(urqmd.isValid(Code::Proton, Code::Electron));
+ CHECK_THROWS(urqmd.isValid(Code::Oxygen, Code::Oxygen));
+ CHECK_THROWS(urqmd.isValid(Code::PiPlus, Code::Omega));
+ CHECK_THROWS(
+ urqmd.isValid(Code::PiPlus, Code::Proton)); // Proton is not a valid target....
+
+ CHECK_NOTHROW(urqmd.isValid(Code::Proton, Code::Oxygen));
+ CHECK_NOTHROW(urqmd.isValid(Code::PiPlus, Code::Argon));
}
- SECTION("nucleus projectile") {
- 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
+ SECTION("cross sections") {
+ FourMomentum const targetP4{Nitrogen::mass, {rootCS, {0_eV, 0_eV, 0_eV}}};
- unsigned short constexpr A = 14, Z = 7;
- auto [stackPtr, secViewPtr] = setup::testing::setup_stack(
- get_nucleus_code(A, Z), 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr,
- *csPtr);
- [[maybe_unused]] setup::StackView& viewRef = *(secViewPtr.get());
- CHECK(stackPtr->getEntries() == 1);
- CHECK(secViewPtr->getEntries() == 0);
+ HEPMomentumType const P0 = 100_GeV;
+ Code const validProjectileCodes[] = {
+ Code::PiPlus, Code::PiMinus, Code::Proton, Code::AntiProton, Code::AntiNeutron,
+ Code::Neutron, Code::KPlus, Code::KMinus, Code::K0Long};
+ // Code::K0, Code::K0Bar are not valid projectiles (no mass eigenstates)
+ CrossSectionType const checkCX[] = {219_mb, 222_mb, 303_mb, 324_mb, 324_mb,
+ 303_mb, 189_mb, 198_mb, 172_mb};
- // 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) == Z + get_charge_number(Code::Oxygen));
-
- auto const secMomSum =
- sumMomentum(*secViewPtr, projectileMomentum.getCoordinateSystem());
- CHECK((secMomSum - projectileMomentum).getNorm() / projectileMomentum.getNorm() ==
- Approx(0).margin(1e-2));
+ int i = 0;
+ for (auto code : validProjectileCodes) {
+ FourMomentum const projectileP4{
+ sqrt(static_pow<2>(get_mass(code)) + static_pow<2>(P0)),
+ {rootCS, {0_GeV, 0_GeV, P0}}};
+ auto const cx = urqmd.getCrossSection(code, Code::Nitrogen, projectileP4, targetP4);
+ CORSIKA_LOG_INFO("UrQMD cross seciton for {} is {} mb", code, cx / 1_mb);
+ CHECK(cx / 1_mb == Approx(checkCX[i++] / 1_mb).margin(1));
+ }
}
- SECTION("\"special\" projectile") {
+ SECTION("pion+ projectile") {
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
@@ -151,7 +124,11 @@ TEST_CASE("UrQMD") {
auto projectile = secViewPtr->getProjectile();
auto const projectileMomentum = projectile.getMomentum();
- urqmd.doInteraction(*secViewPtr);
+ FourMomentum const projectileP4{
+ sqrt(static_pow<2>(PiPlus::mass) + static_pow<2>(40_GeV)),
+ {rootCS, {40_GeV, 0_GeV, 0_GeV}}};
+ FourMomentum const targetP4{Oxygen::mass, {rootCS, {0_GeV, 0_GeV, 0_GeV}}};
+ urqmd.doInteraction(*secViewPtr, Code::PiPlus, Code::Oxygen, projectileP4, targetP4);
CHECK(sumCharge(*secViewPtr) ==
get_charge_number(Code::PiPlus) + get_charge_number(Code::Oxygen));
@@ -162,44 +139,6 @@ 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, 40_GeV, (setup::Environment::BaseNodeType* const)nodePtr, *csPtr);
- [[maybe_unused]] auto particle = stackPtr->first();
- 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, 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
@@ -214,7 +153,11 @@ TEST_CASE("UrQMD") {
auto projectile = secViewPtr->getProjectile();
auto const projectileMomentum = projectile.getMomentum();
- urqmd.doInteraction(*secViewPtr);
+ FourMomentum const projectileP4{
+ sqrt(static_pow<2>(K0Long::mass) + static_pow<2>(40_GeV)),
+ {rootCS, {40_GeV, 0_GeV, 0_GeV}}};
+ FourMomentum const targetP4{Oxygen::mass, {rootCS, {0_GeV, 0_GeV, 0_GeV}}};
+ urqmd.doInteraction(*secViewPtr, Code::K0Long, Code::Oxygen, projectileP4, targetP4);
CHECK(sumCharge(*secViewPtr) ==
get_charge_number(Code::K0Long) + get_charge_number(Code::Oxygen));
--
GitLab