From d36ec1ca9161e2385a7c6f227a31bc2b4143f319 Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@tu-dortmund.de>
Date: Sat, 27 Apr 2019 22:19:23 -0300
Subject: [PATCH] infinite interaction length for non-hadrons
---
Processes/UrQMD/UrQMD.cc | 36 +++++++++++++++++++++++++++++-------
Processes/UrQMD/UrQMD.h | 6 ++++--
Processes/UrQMD/testUrQMD.cc | 14 ++++++--------
3 files changed, 39 insertions(+), 17 deletions(-)
diff --git a/Processes/UrQMD/UrQMD.cc b/Processes/UrQMD/UrQMD.cc
index 9a88c8588..1f7e44215 100644
--- a/Processes/UrQMD/UrQMD.cc
+++ b/Processes/UrQMD/UrQMD.cc
@@ -7,6 +7,7 @@
#include <corsika/setup/SetupTrajectory.h>
#include <corsika/units/PhysicalUnits.h>
+#include <algorithm>
#include <functional>
#include <iostream>
#include <random>
@@ -23,10 +24,8 @@ using SetupTrack = corsika::setup::Trajectory;
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 {
+CrossSectionType UrQMD::GetCrossSection(TParticle const& vProjectile,
+ corsika::particles::Code vTargetCode) const {
using namespace units::si;
// TODO: energy cuts, return 0 for non-hadrons
@@ -39,9 +38,10 @@ template <typename TParticle> // need template here, as this is called both with
!IsNucleus(vTargetCode)) { // both particles are "special"
auto const mProj = particles::GetMass(projectileCode);
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 * projectileEnergyLab * mTar) *
+ (1 / 1_GeV);
// we must set some UrQMD globals first...
auto const [ityp, iso3] = ConvertToUrQMD(projectileCode);
@@ -66,7 +66,29 @@ template <typename TParticle> // need template here, as this is called both with
}
}
+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
+ // only the usual long-lived species as input.
+ // TODO: Charmed mesons should be added to the list, too
+
+ static particles::Code const validProjectileCodes[] = {
+ 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};
+
+ return std::find(std::cbegin(validProjectileCodes), std::cend(validProjectileCodes),
+ vCode) != std::cend(validProjectileCodes);
+}
+
GrammageType UrQMD::GetInteractionLength(SetupParticle& vParticle, SetupTrack&) const {
+ 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);
+ }
+
auto const& mediumComposition =
vParticle.GetNode()->GetModelProperties().GetNuclearComposition();
using namespace std::placeholders;
diff --git a/Processes/UrQMD/UrQMD.h b/Processes/UrQMD/UrQMD.h
index 798d86d7a..0b34e2e14 100644
--- a/Processes/UrQMD/UrQMD.h
+++ b/Processes/UrQMD/UrQMD.h
@@ -19,12 +19,14 @@ namespace corsika::process::UrQMD {
corsika::setup::Stack::StackIterator&, corsika::setup::Trajectory&) const;
template <typename TParticle>
- corsika::units::si::CrossSectionType GetCrossSection(
- TParticle const&, corsika::particles::Code) const;
+ corsika::units::si::CrossSectionType GetCrossSection(TParticle const&,
+ corsika::particles::Code) const;
corsika::process::EProcessReturn DoInteraction(
corsika::setup::StackView::StackIterator&);
+ bool CanInteract(particles::Code) const;
+
private:
corsika::random::RNG& fRNG =
corsika::random::RNGManager::GetInstance().GetRandomStream("UrQMD");
diff --git a/Processes/UrQMD/testUrQMD.cc b/Processes/UrQMD/testUrQMD.cc
index f65838d32..0a9af8749 100644
--- a/Processes/UrQMD/testUrQMD.cc
+++ b/Processes/UrQMD/testUrQMD.cc
@@ -42,10 +42,8 @@ using namespace corsika::units::si;
template <typename TStackView>
auto sumCharge(TStackView const& view) {
int totalCharge = 0;
-
- for (auto const& p : view) {
- totalCharge += particles::GetChargeNumber(p.GetPID());
- }
+
+ for (auto const& p : view) { totalCharge += particles::GetChargeNumber(p.GetPID()); }
return totalCharge;
}
@@ -165,7 +163,7 @@ TEST_CASE("UrQMD") {
auto [stackPtr, secViewPtr] = setupStack(A, Z, 400_GeV, nodePtr, *csPtr);
// must be assigned to variable, cannot be used as rvalue?!
- auto projectile = secViewPtr ->GetProjectile();
+ auto projectile = secViewPtr->GetProjectile();
[[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(projectile);
REQUIRE(sumCharge(*secViewPtr) ==
@@ -175,14 +173,14 @@ TEST_CASE("UrQMD") {
SECTION("\"special\" projectile") {
auto [env, csPtr, nodePtr] = setupEnvironment(particles::Code::Oxygen);
auto [stackPtr, secViewPtr] =
- setupStack(particles::Code::PiPlus, 400_GeV, nodePtr, *csPtr);
+ setupStack(particles::Code::K0, 400_GeV, nodePtr, *csPtr);
// must be assigned to variable, cannot be used as rvalue?!
- auto projectile = secViewPtr->GetProjectile();
+ auto projectile = secViewPtr->GetProjectile();
[[maybe_unused]] process::EProcessReturn const ret = urqmd.DoInteraction(projectile);
REQUIRE(sumCharge(*secViewPtr) ==
- particles::GetChargeNumber(particles::Code::PiPlus) +
+ particles::GetChargeNumber(particles::Code::K0) +
particles::GetChargeNumber(particles::Code::Oxygen));
}
}
--
GitLab