From f98c436d4838e5e794073850b2be967047d221fc Mon Sep 17 00:00:00 2001
From: Maximilian Reininghaus <maximilian.reininghaus@tu-dortmund.de>
Date: Wed, 8 Apr 2020 21:54:31 +0200
Subject: [PATCH] resemble ptsigtot() better
---
Processes/UrQMD/UrQMD.cc | 39 +++++++++++++++++++++++++++++---
Processes/UrQMD/UrQMD.h | 13 ++++++++---
Processes/UrQMD/urqmdInterface.F | 13 +++++++++++
3 files changed, 59 insertions(+), 6 deletions(-)
diff --git a/Processes/UrQMD/UrQMD.cc b/Processes/UrQMD/UrQMD.cc
index 6f4628781..c91dca3e5 100644
--- a/Processes/UrQMD/UrQMD.cc
+++ b/Processes/UrQMD/UrQMD.cc
@@ -30,8 +30,9 @@ using SetupProjectile = corsika::setup::StackView::StackIterator;
CrossSectionType UrQMD::GetCrossSection(particles::Code vProjectileCode,
corsika::particles::Code vTargetCode,
- HEPEnergyType vLabEnergy, int vAProjectile = 1) {
- // the following is a translation of ptsigtot() into C++
+ HEPEnergyType vLabEnergy,
+ int vAProjectile) const {
+ // the following is a (incomplete!) translation of ptsigtot() into C++
if (vProjectileCode != particles::Code::Nucleus &&
!IsNucleus(vTargetCode)) { // both particles are "special"
auto const mProj = particles::GetMass(vProjectileCode);
@@ -51,7 +52,39 @@ CrossSectionType UrQMD::GetCrossSection(particles::Code vProjectileCode,
int one = 1;
int two = 2;
- return sigtot_(one, two, sqrtS) * 1_mb;
+ int three = 3;
+
+ double const totalXS = sigtot_(one, two, sqrtS);
+
+ // subtract elastic cross-section as in ptsigtot()
+ int itypmn, itypmx, iso3mn, iso3mx;
+ if (ityp < itypTar) {
+ itypmn = ityp;
+ itypmx = itypTar;
+
+ iso3mn = iso3;
+ iso3mx = iso3Tar;
+ } else {
+ itypmx = ityp;
+ itypmn = itypTar;
+
+ iso3mx = iso3;
+ iso3mn = iso3Tar;
+ }
+
+ int isigline = collclass_(itypmx, iso3mx, itypmn, iso3mn);
+ int iline = readsigmaln_(three, one, isigline);
+ double sigEl;
+ double massProj = mProj / 1_GeV;
+ double massTar = mTar / 1_GeV;
+
+ crossx_(iline, sqrtS, ityp, iso3, massProj, itypTar, iso3Tar, massTar, sigEl);
+
+ if (totalXS > sigEl) {
+ return (totalXS - sigEl) * 1_mb;
+ } else {
+ return sigEl * 0_mb;
+ }
} else {
int const Ap = vAProjectile;
int const At = IsNucleus(vTargetCode) ? particles::GetNucleusA(vTargetCode) : 1;
diff --git a/Processes/UrQMD/UrQMD.h b/Processes/UrQMD/UrQMD.h
index 402b16ca8..d74fed7b8 100644
--- a/Processes/UrQMD/UrQMD.h
+++ b/Processes/UrQMD/UrQMD.h
@@ -33,14 +33,16 @@ namespace corsika::process::UrQMD {
corsika::units::si::CrossSectionType GetCrossSection(TParticle const&,
corsika::particles::Code) const;
+ corsika::units::si::CrossSectionType GetCrossSection(
+ particles::Code, particles::Code, corsika::units::si::HEPEnergyType,
+ int Ap = 1) const;
+
corsika::process::EProcessReturn DoInteraction(
corsika::setup::StackView::StackIterator&);
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");
@@ -66,13 +68,18 @@ namespace corsika::process::UrQMD {
using nmaxDoubleArray = nmaxArray<double>;
extern "C" {
+ // FORTRAN functions defined in UrQMD
void iniurqmd_();
double ranf_(int&);
void cascinit_(int const&, int const&, int const&);
double nucrad_(int const&);
void urqmd_(int&);
- int pdgid_(int&, int&);
+ int pdgid_(int const&, int const&);
double sigtot_(int&, int&, double&);
+ int collclass_(int&, int&, int&, int&);
+ double crossx_(int const&, double const&, int const&, int const&, double const&,
+ int const&, int const&, double const&, double&);
+ int readsigmaln_(int const&, int const&, int const&);
// defined in coms.f
extern struct {
diff --git a/Processes/UrQMD/urqmdInterface.F b/Processes/UrQMD/urqmdInterface.F
index 4fe9c044e..259e3fd47 100644
--- a/Processes/UrQMD/urqmdInterface.F
+++ b/Processes/UrQMD/urqmdInterface.F
@@ -432,3 +432,16 @@ c~ xsegymin=0.25d0
end
+c
+c M. Reininghaus, 2020-04-08
+c
+ integer function ReadSigmaLn(ia, ib, ic)
+ implicit none
+
+ include 'comres.f'
+
+ integer :: ia, ib, ic
+
+ ReadSigmaLn = SigmaLn(ia, ib, ic)
+
+ end function ReadSigmaLn
--
GitLab