diff --git a/Processes/UrQMD/UrQMD.cc b/Processes/UrQMD/UrQMD.cc index 6f462878111a12b619e37718d18b769b8cd634b0..c91dca3e53c817535e05a4b97524052307ca9143 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 402b16ca8e3f2b88b20ca6e4acca80066a040f04..d74fed7b8ca28168b4d5c6e00d7a9253a5fd3b1c 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 4fe9c044e447f577559444be28e0798ab3cd678e..259e3fd4739344c2fe5f67a8fa14bf4cd7ac17dc 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