From 9e9d3cfc8632f9b812532790adbe3539e38477ea Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@kit.edu>
Date: Thu, 25 Apr 2019 19:11:59 -0300
Subject: [PATCH] implemented UrQMD::GetCrossSection
---
Processes/UrQMD/UrQMD.cc | 84 ++++++++++++++++++++++++++----------
Processes/UrQMD/UrQMD.h | 4 +-
Processes/UrQMD/init.f | 3 +-
Processes/UrQMD/testUrQMD.cc | 5 ++-
4 files changed, 69 insertions(+), 27 deletions(-)
diff --git a/Processes/UrQMD/UrQMD.cc b/Processes/UrQMD/UrQMD.cc
index 2f48668b7..452714fd8 100644
--- a/Processes/UrQMD/UrQMD.cc
+++ b/Processes/UrQMD/UrQMD.cc
@@ -21,45 +21,83 @@ using SetupParticle = corsika::setup::Stack::StackIterator;
using SetupProjectile = corsika::setup::StackView::StackIterator;
using SetupTrack = corsika::setup::Trajectory;
-CrossSectionType UrQMD::GetCrossSection(
- corsika::particles::Code projectileID, corsika::particles::Code targetID,
- corsika::units::si::HEPEnergyType projectileEnergy) const {
+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,
+ corsika::units::si::HEPEnergyType vProjectileEnergyLab)
+ const {
using namespace units::si;
- return 10_mb; // TODO: implement
+
+ // TODO: energy cuts, return 0 for non-hadrons
+
+ auto const projectileCode = vProjectile.GetPID();
+
+ // the following is a translation of ptsigtot() into C++
+ if (projectileCode != particles::Code::Nucleus &&
+ !IsNucleus(vTargetCode)) { // both particles are "special"
+ auto const mProj = particles::GetMass(projectileCode);
+ auto const mTar = particles::GetMass(vTargetCode);
+ double sqrtS = sqrt(detail::static_pow<2>(mProj) + detail::static_pow<2>(mTar) +
+ 2 * vProjectileEnergyLab * mTar) *
+ (1 / 1_GeV);
+
+ // we must set some UrQMD globals first...
+ auto const [ityp, iso3] = ConvertToUrQMD(projectileCode);
+ inputs_.spityp[0] = ityp;
+ inputs_.spiso3[0] = iso3;
+
+ auto const [itypTar, iso3Tar] = ConvertToUrQMD(vTargetCode);
+ inputs_.spityp[1] = itypTar;
+ inputs_.spiso3[1] = iso3Tar;
+
+ int one = 1;
+ int two = 2;
+ return sigtot_(one, two, sqrtS) * 1_mb;
+ } else {
+ int const Ap =
+ (projectileCode == particles::Code::Nucleus) ? vProjectile.GetNuclearA() : 1;
+ int const At = IsNucleus(vTargetCode) ? particles::GetNucleusA(vTargetCode) : 1;
+
+ double const maxImpact = nucrad_(Ap) + nucrad_(At) + 2 * options_.CTParam[30 - 1];
+ return 10_mb * M_PI * detail::static_pow<2>(maxImpact);
+ // is a constant cross-section really reasonable?
+ }
}
-GrammageType UrQMD::GetInteractionLength(SetupParticle& particle, SetupTrack&) const {
+GrammageType UrQMD::GetInteractionLength(SetupParticle& vParticle, SetupTrack&) const {
auto const& mediumComposition =
- particle.GetNode()->GetModelProperties().GetNuclearComposition();
+ vParticle.GetNode()->GetModelProperties().GetNuclearComposition();
using namespace std::placeholders;
- CrossSectionType const weightedProdCrossSection =
- mediumComposition.WeightedSum(std::bind(
- &UrQMD::GetCrossSection, this, particle.GetPID(), _1, particle.GetEnergy()));
+ CrossSectionType const weightedProdCrossSection = mediumComposition.WeightedSum(
+ std::bind(&UrQMD::GetCrossSection<decltype(vParticle)>, this, vParticle, _1,
+ vParticle.GetEnergy()));
return mediumComposition.GetAverageMassNumber() * units::constants::u /
weightedProdCrossSection;
}
-corsika::process::EProcessReturn UrQMD::DoInteraction(SetupProjectile& projectile) {
+corsika::process::EProcessReturn UrQMD::DoInteraction(SetupProjectile& vProjectile) {
using namespace units::si;
- auto const projectileCode = projectile.GetPID();
- auto const projectileEnergyLab = projectile.GetEnergy();
- auto const& projectileMomentumLab = projectile.GetMomentum();
- auto const& projectilePosition = projectile.GetPosition();
- auto const projectileTime = projectile.GetTime();
+ auto const projectileCode = vProjectile.GetPID();
+ auto const projectileEnergyLab = vProjectile.GetEnergy();
+ auto const& projectileMomentumLab = vProjectile.GetMomentum();
+ auto const& projectilePosition = vProjectile.GetPosition();
+ auto const projectileTime = vProjectile.GetTime();
// sample target particle
auto const& mediumComposition =
- projectile.GetNode()->GetModelProperties().GetNuclearComposition();
+ vProjectile.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));
+ crossSections.push_back(GetCrossSection(vProjectile, c, projectileEnergyLab));
}
return crossSections;
@@ -79,8 +117,8 @@ corsika::process::EProcessReturn UrQMD::DoInteraction(SetupProjectile& projectil
// is this everything?
inputs_.prspflg = 0;
- sys_.Ap = projectile.GetNuclearA();
- sys_.Zp = projectile.GetNuclearZ();
+ sys_.Ap = vProjectile.GetNuclearA();
+ sys_.Zp = vProjectile.GetNuclearZ();
rsys_.bdist = nucrad_(targetA) + nucrad_(sys_.Ap) + 2 * options_.CTParam[30 - 1];
@@ -97,7 +135,7 @@ corsika::process::EProcessReturn UrQMD::DoInteraction(SetupProjectile& projectil
inputs_.spiso3[0] = iso3;
}
- rsys_.ebeam = (projectileEnergyLab - projectile.GetMass()) * (1 / 1_GeV);
+ rsys_.ebeam = (projectileEnergyLab - vProjectile.GetMass()) * (1 / 1_GeV);
// initilazation regarding target
if (particles::IsNucleus(targetCode)) {
@@ -113,7 +151,7 @@ corsika::process::EProcessReturn UrQMD::DoInteraction(SetupProjectile& projectil
inputs_.spiso3[1] = iso3;
}
- int iflb = 0; // flag for retrying interaction in case of elastic scattering
+ int iflb = 0; // flag for retrying interaction in case of empty event, 0 means retry
urqmd_(iflb);
// now retrieve secondaries from UrQMD
@@ -132,7 +170,7 @@ corsika::process::EProcessReturn UrQMD::DoInteraction(SetupProjectile& projectil
momentum.rebase(originalCS); // transform back into standard lab frame
std::cout << i << " " << code << " " << momentum.GetComponents() << std::endl;
- projectile.AddSecondary(
+ vProjectile.AddSecondary(
std::tuple<particles::Code, HEPEnergyType, stack::MomentumVector, geometry::Point,
TimeType>{code, energy, momentum, projectilePosition, projectileTime});
}
@@ -140,7 +178,7 @@ corsika::process::EProcessReturn UrQMD::DoInteraction(SetupProjectile& projectil
std::cout << "UrQMD generated " << sys_.npart << " secondaries!" << std::endl;
if (sys_.npart > 0) // delete only in case of inelastic collision, otherwise keep
- projectile.Delete();
+ vProjectile.Delete();
return process::EProcessReturn::eOk;
}
diff --git a/Processes/UrQMD/UrQMD.h b/Processes/UrQMD/UrQMD.h
index ca3b4b398..c32082a83 100644
--- a/Processes/UrQMD/UrQMD.h
+++ b/Processes/UrQMD/UrQMD.h
@@ -18,8 +18,9 @@ namespace corsika::process::UrQMD {
corsika::units::si::GrammageType GetInteractionLength(
corsika::setup::Stack::StackIterator&, corsika::setup::Trajectory&) const;
+ template <typename TParticle>
corsika::units::si::CrossSectionType GetCrossSection(
- corsika::particles::Code, corsika::particles::Code,
+ TParticle const&, corsika::particles::Code,
corsika::units::si::HEPEnergyType) const;
corsika::process::EProcessReturn DoInteraction(
@@ -55,6 +56,7 @@ namespace corsika::process::UrQMD {
double nucrad_(int const&);
void urqmd_(int&);
int pdgid_(int&, int&);
+ double sigtot_(int&, int&, double&);
// defined in coms.f
extern struct {
diff --git a/Processes/UrQMD/init.f b/Processes/UrQMD/init.f
index d934f7f74..76afc9dfa 100644
--- a/Processes/UrQMD/init.f
+++ b/Processes/UrQMD/init.f
@@ -338,11 +338,12 @@ c determine impact parameter
if(CTOption(5).eq.0) then
bimp=bdist
elseif(CTOption(5).eq.1) then
+C M.R. we don't truncate bdist here, logic happens in CORSIKA
c hjd1
c if(bdist.gt.(nucrad(Ap)+nucrad(At)+2*CTParam(30)))
c & bdist=nucrad(Ap)+nucrad(At)+2*CTParam(30)
c hjd1
-c ! M.R. 2019-04-24: updated sampling procedure from UrQMD 3.4
+C M.R. 2019-04-24: updated sampling procedure from UrQMD 3.4
bimp=sqrt(bmin**2 + ranf(0) * (bdist**2 - bmin**2))
cdh if (bimp<bmin) goto 215
diff --git a/Processes/UrQMD/testUrQMD.cc b/Processes/UrQMD/testUrQMD.cc
index 6cddf3e92..fec1f3ea0 100644
--- a/Processes/UrQMD/testUrQMD.cc
+++ b/Processes/UrQMD/testUrQMD.cc
@@ -41,9 +41,10 @@ auto sumCharge(TStack& stack) {
int totalCharge = 0;
int count = 0;
for (auto& p : stack) {
- count++;
+ std::cout << count++ << " ";
totalCharge += particles::GetChargeNumber(p.GetPID());
- std::cout << p.GetPID() << " " << particles::GetChargeNumber(p.GetPID()) << std::endl;
+ std::cout << p.GetPID() << " " << particles::GetChargeNumber(p.GetPID()) << ' '
+ << p.GetMomentum().GetComponents() << std::endl;
}
std::cout << count << " particles on stack" << std::endl;
--
GitLab