diff --git a/Processes/UrQMD/CMakeLists.txt b/Processes/UrQMD/CMakeLists.txt
index 3ba71c209f3d3e746458c2a80e4c1623f4684981..4c8d8e2cbe6eacf4a1049e6a0d49529d56a89e89 100644
--- a/Processes/UrQMD/CMakeLists.txt
+++ b/Processes/UrQMD/CMakeLists.txt
@@ -58,6 +58,7 @@ set_target_properties (
# target dependencies on other libraries (also the header onlys)
target_link_libraries (
ProcessUrQMD
+ CORSIKAprocesssequence
CORSIKAparticles
CORSIKAunits
CORSIKAgeometry
diff --git a/Processes/UrQMD/UrQMD.cc b/Processes/UrQMD/UrQMD.cc
index 81f034de544423b58e51c3b56bcc60a86eed65f2..72d773cd95397e964f5e7b85ea76770baadc6231 100644
--- a/Processes/UrQMD/UrQMD.cc
+++ b/Processes/UrQMD/UrQMD.cc
@@ -1,13 +1,289 @@
+#include <corsika/geometry/QuantityVector.h>
+#include <corsika/geometry/Vector.h>
+#include <corsika/particles/ParticleProperties.h>
#include <corsika/process/urqmd/UrQMD.h>
#include <corsika/random/RNGManager.h>
+#include <corsika/setup/SetupStack.h>
+#include <corsika/setup/SetupTrajectory.h>
+#include <corsika/units/PhysicalUnits.h>
+
+#include <functional>
+#include <iostream>
#include <random>
-corsika::process::UrQMD::UrQMD::UrQMD() { iniurqmd_(); }
+using namespace corsika::process::UrQMD;
+using namespace corsika::units::si;
+
+UrQMD::UrQMD() { iniurqmd_(); }
+
+using SetupStack = corsika::setup::Stack;
+using SetupParticle = SetupStack::StackIterator;
+using SetupTrack = corsika::setup::Trajectory;
+
+CrossSectionType UrQMD::GetCrossSection(
+ corsika::particles::Code projectileID, corsika::particles::Code targetID,
+ corsika::units::si::HEPEnergyType projectileEnergy) const {
+ using namespace units::si;
+ return 10_mb; // TODO: implement
+}
+
+template <>
+GrammageType UrQMD::GetInteractionLength(SetupParticle& particle, SetupTrack&) const {
+ auto const& mediumComposition =
+ particle.GetNode()->GetModelProperties().GetNuclearComposition();
+ using namespace std::placeholders;
+
+ CrossSectionType const weightedProdCrossSection =
+ mediumComposition.WeightedSum(std::bind(
+ &UrQMD::GetCrossSection, this, particle.GetPID(), _1, particle.GetEnergy()));
+
+ return mediumComposition.GetAverageMassNumber() * units::constants::u /
+ weightedProdCrossSection;
+}
+
+template <>
+corsika::process::EProcessReturn UrQMD::DoInteraction(SetupParticle& projectile,
+ SetupStack&) {
+ 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();
+
+ inputs_.nevents = 1;
+ sys_.eos = 0; // could be configurable in principle
+ inputs_.outsteps = 1;
+ sys_.nsteps = 1;
+
+ // todo: sample target
+
+ int const Atarget = 14;
+ int tableIndex = 0; // 0: nitrogen, 1: oxygen, 2: argon target
+
+ /*
+ // corsika 7
+ bmin = 0.d0
+ CTOption(5) = 1
+ if ( iflbmax.eq.1 ) then
+ bdist = BIM(LIT)
+ else
+ bdist=nucrad(Ap)+nucrad(At)+2*CTParam(30)
+ endif
+
+ // conex
+ CTOption(5)=1
+ if ( prspflg.eq.0 ) then
+ bdist = BIM(LT)
+ else
+ bdist = xsbmax
+ endif
+ */
+
+ // initialization regarding projectile
+ if (particles::Code::Nucleus == projectileCode) {
+ // is this everything?
+ inputs_.prspflg = 0;
+ rsys_.bdist = cxs_u2_.bim[tableIndex];
+
+ sys_.Ap = projectile.GetNuclearA();
+ sys_.Zp = projectile.GetNuclearZ();
+
+ int const id = 1;
+ cascinit_(sys_.Zp, sys_.Ap, id);
+ } else {
+ inputs_.prspflg = 1;
+ int const Ap =
+ 1; // what value to use here for non-baryons??? see CONEX UrQMD interface
+ rsys_.bdist = nucrad_(Ap) + nucrad_(Atarget) + 2 * options_.CTParam[30 - 1];
+
+ auto const [ityp, iso3] = ConvertToUrQMD(projectileCode);
+ // todo: conversion of K_long/short into strong eigenstates;
+ inputs_.spityp[0] = ityp;
+ inputs_.spiso3[0] = iso3;
+ }
-double ranf_(int*) {
+ rsys_.ebeam = (projectileEnergyLab - particles::GetMass(projectileCode)) * (1 / 1_GeV);
+
+ // initilazation regarding target
+ 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());
+
+ for (auto const c : components) {
+ crossSections.push_back(GetCrossSection(projectileCode, c, projectileEnergyLab));
+ }
+
+ return crossSections;
+ });
+
+ auto const targetCode = mediumComposition.SampleTarget(componentCrossSections, fRNG);
+
+ if (particles::IsNucleus(targetCode)) {
+ sys_.Zt = particles::GetNucleusZ(targetCode);
+ sys_.At = particles::GetNucleusA(targetCode);
+ inputs_.trspflg = 0; // nucleus as target
+ int const id = 2;
+ cascinit_(sys_.Zt, sys_.At, id);
+ } else {
+ inputs_.trspflg = 1; // special particle as target
+ auto const [ityp, iso3] = ConvertToUrQMD(targetCode);
+ inputs_.spityp[1] = ityp;
+ inputs_.spiso3[1] = iso3;
+ }
+
+ int iflb = 0; // flag for retrying interaction in case of elastic scattering
+ urqmd_(iflb);
+
+ // now retrieve secondaries from UrQMD
+ auto const& originalCS = projectileMomentumLab.GetCoordinateSystem();
+ geometry::CoordinateSystem const zAxisFrame =
+ originalCS.RotateToZ(projectileMomentumLab);
+
+ for (int i = 0; i < sys_.npart; ++i) {
+ auto const code = ConvertFromUrQMD(isys_.ityp[i], isys_.iso3[i]);
+ auto const energy = coor_.p0[i] * 1_GeV;
+ auto momentum = geometry::Vector(
+ zAxisFrame,
+ geometry::QuantityVector<dimensionless_d>{coor_.px[i], coor_.py[i], coor_.pz[i]} *
+ 1_GeV);
+
+ momentum.rebase(originalCS); // transform back into standard lab frame
+
+ projectile.AddSecondary(
+ std::tuple<particles::Code, HEPEnergyType, stack::MomentumVector, geometry::Point,
+ TimeType>{code, energy, momentum, projectilePosition, projectileTime});
+ }
+
+ if (sys_.npart > 0) // delete only in case of inelastic collision, otherwise keep
+ projectile.Delete();
+}
+
+/**
+ * the random number generator function of UrQMD
+ */
+double corsika::process::UrQMD::ranf_(int&) {
static corsika::random::RNG& rng =
- corsika::random::RNGManager::GetInstance().GetRandomStream("ranf");
+ corsika::random::RNGManager::GetInstance().GetRandomStream("UrQMD");
+ static std::uniform_real_distribution<double> dist;
- std::uniform_real_distribution<double> dist;
return dist(rng);
}
+
+corsika::particles::Code corsika::process::UrQMD::ConvertFromUrQMD(int vItyp, int vIso3) {
+ static const std::map<std::pair<int, int>, int> mapUrQMDToPDG{
+ // data from github.com/afedynitch/ParticleDataTool
+ {{100, 0}, 22}, // gamma
+ {{101, 0}, 111}, // pi0
+ {{101, 2}, 211}, // pi+
+ {{101, -2}, -211}, // pi-
+ {{106, 1}, 321}, // K+
+ {{106, -1}, -321}, // K-
+ {{1, 1}, 2212}, // p
+ {{1, -1}, 2112}, // n
+ {{102, 0}, 221}, // eta
+ {{104, 2}, 213}, // rho+
+ {{104, -2}, -213}, // rho-
+ {{104, 0}, 113}, // rho0
+ {{108, 2}, 323}, // K*+
+ {{108, -2}, -323}, // K*-
+ {{108, 0}, 313}, // K*0
+ {{-108, 0}, -313}, // K*0-bar
+ {{103, 0}, 223}, // omega
+ {{109, 0}, 333}, // phi
+ {{40, 2}, 3222}, // Sigma+
+ {{40, 0}, 3212}, // Sigma0
+ {{40, -2}, 3112}, // Sigma-
+ {{49, 0}, 3322}, // Xi0
+ {{49, -1}, 3312}, // Xi-
+ {{27, 0}, 3122}, // Lambda0
+ {{17, 4}, 2224}, // Delta++
+ {{17, 2}, 2214}, // Delta+
+ {{17, 0}, 2114}, // Delta0
+ {{17, -2}, 1114}, // Delta-
+ {{41, 2}, 3224}, // Sigma*+
+ {{41, 0}, 3214}, // Sigma*0
+ {{41, -2}, 3114}, // Sigma*-
+ {{50, 0}, 3324}, // Xi*0
+ {{50, -1}, 3314}, // Xi*-
+ {{55, 0}, 3334}, // Omega-
+ {{133, 2}, 411}, // D+
+ {{133, -2}, -411}, // D-
+ {{133, 0}, 421}, // D0
+ {{-133, 0}, -421}, // D0-bar
+ {{107, 0}, 441}, // etaC
+ {{138, 1}, 431}, // Ds+
+ {{138, -1}, -431}, // Ds-
+ {{139, 1}, 433}, // Ds*+
+ {{139, -1}, -433}, // Ds*-
+ {{134, 1}, 413}, // D*+
+ {{134, -1}, -413}, // D*-
+ {{134, 0}, 10421}, // D*0
+ {{-134, 0}, -10421}, // D*0-bar
+ {{135, 0}, 443}, // jpsi
+ };
+
+ auto const pdg = static_cast<particles::PDGCode>(mapUrQMDToPDG.at({vItyp, vIso3}));
+ return particles::ConvertFromPDG(pdg);
+}
+
+std::pair<int, int> corsika::process::UrQMD::ConvertToUrQMD(
+ corsika::particles::Code code) {
+ static const std::map<int, std::pair<int, int>> mapPDGToUrQMD{
+ // data from github.com/afedynitch/ParticleDataTool
+ {22, {100, 0}}, // gamma
+ {111, {101, 0}}, // pi0
+ {211, {101, 2}}, // pi+
+ {-211, {101, -2}}, // pi-
+ {321, {106, 1}}, // K+
+ {-321, {106, -1}}, // K-
+ {2212, {1, 1}}, // p
+ {2112, {1, -1}}, // n
+ {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>(GetPDG(code)));
+}
diff --git a/Processes/UrQMD/UrQMD.h b/Processes/UrQMD/UrQMD.h
index aad8b45e3cf5eb715e3c32e64a5dd6455314ac01..9693d6eccac2606872caabb8a28dd25f1dba63a9 100644
--- a/Processes/UrQMD/UrQMD.h
+++ b/Processes/UrQMD/UrQMD.h
@@ -1,118 +1,122 @@
#ifndef _Processes_UrQMD_UrQMD_h
#define _Processes_UrQMD_UrQMD_h
+#include <corsika/particles/ParticleProperties.h>
+#include <corsika/process/InteractionProcess.h>
+#include <corsika/random/RNGManager.h>
+#include <corsika/units/PhysicalUnits.h>
+
#include <array>
+#include <utility>
namespace corsika::process::UrQMD {
- class UrQMD {
+ class UrQMD : public corsika::process::InteractionProcess<UrQMD> {
public:
UrQMD();
+
+ template <typename Particle, typename Track>
+ corsika::units::si::GrammageType GetInteractionLength(Particle&, Track&) const;
+
+ corsika::units::si::CrossSectionType GetCrossSection(
+ corsika::particles::Code, corsika::particles::Code,
+ corsika::units::si::HEPEnergyType) const;
+
+ template <typename Particle, typename Stack>
+ corsika::process::EProcessReturn DoInteraction(Particle&, Stack&);
+
+ private:
+ corsika::random::RNG& fRNG =
+ corsika::random::RNGManager::GetInstance().GetRandomStream("UrQMD");
};
namespace constants {
// from coms.f
int constexpr nmax = 500;
- int constexpr nspl = 500;
- int constexpr smax = 500;
- // from comres.f
- int constexpr minnuc = 1;
- int constexpr minmes = 100;
- int constexpr maxmes = 132;
- int constexpr numnuc = 16;
- int constexpr numdel = 10;
- int constexpr maxnuc = minnuc + numnuc - 1;
- int constexpr mindel = minnuc + maxnuc;
- int constexpr maxdel = mindel + numdel - 1;
- int constexpr minres = minnuc + 1;
- int constexpr maxres = maxdel;
- int constexpr numlam = 13;
- int constexpr numsig = 9;
- int constexpr numcas = 6;
- int constexpr numome = 1;
- int constexpr minlam = mindel + numdel;
- int constexpr maxlam = minlam + numlam - 1;
- int constexpr minsig = minlam + numlam;
- int constexpr maxsig = minsig + numsig - 1;
- int constexpr mincas = minsig + numsig;
- int constexpr maxcas = mincas + numcas - 1;
- int constexpr minome = mincas + numcas;
- int constexpr maxome = minome + numome - 1;
- int constexpr minbar = minnuc;
- int constexpr maxbar = maxome;
- int constexpr offmeson = minmes;
- int constexpr maxmeson = maxmes;
- int constexpr maxbra = 11;
- int constexpr maxbrm = 25;
- int constexpr maxbrs1 = 10;
- int constexpr maxbrs2 = 3;
- int constexpr nsigs = 10;
- int constexpr itblsz = 100;
- int constexpr maxreac = 13;
- int constexpr maxpsig = 12;
-
- // from comwid.f
- int constexpr widnsp = 120;
- double constexpr mintab = 0.10;
- double constexpr maxtab1 = 5.0;
- double constexpr maxtab2 = 50.0;
- int constexpr tabver = 9;
// from options.f
int constexpr numcto = 400;
int constexpr numctp = 400;
- int constexpr maxstables = 20;
-
- // from colltab.f
- int constexpr ncollmax = 100;
// from inputs.f
int constexpr aamax = 300;
- // from newpart.f
- int constexpr mprt = 200;
- int constexpr oprt = 2;
-
- // from boxinc.f
- int constexpr bptmax = 20;
-
- // from comnorm.f
- int constexpr n = 400;
-
- // from comstr.f
- int constexpr njspin = 8;
-
- // iso
- int constexpr jmax = 7;
} // namespace constants
- using nmaxIntArray = std::array<int, constants::nmax>;
-} // namespace corsika::process::UrQMD
-
-extern "C" {
-void iniurqmd_();
-double ranf_(int*);
-
-struct {
- int npart, nbar, nmes, ctag, nsteps, uid_cnt, ranseed, event, Ap, At, Zp, Zt, eos,
- dectag, NHardRes, NSoftRes, NDecRes, NElColl, NBlColl;
-} sys_;
+ template <typename T>
+ using nmaxArray = std::array<T, constants::nmax>;
+ using nmaxIntArray = nmaxArray<int>;
+ using nmaxDoubleArray = nmaxArray<double>;
+
+ extern "C" {
+ void iniurqmd_();
+ double ranf_(int&);
+ void cascinit_(int const&, int const&, int const&);
+ double nucrad_(int const&);
+ void urqmd_(int&);
+
+ // defined in coms.f
+ extern struct {
+ int npart, nbar, nmes, ctag, nsteps, uid_cnt, ranseed, event;
+ int Ap; // projectile mass number (in case of nucleus)
+ int At; // target mass number (in case of nucleus)
+ int Zp; // projectile charge number (in case of nucleus)
+ int Zt; // target charge number (in case of nucleus)
+ int eos, dectag, NHardRes, NSoftRes, NDecRes, NElColl, NBlColl;
+ } sys_;
+
+ extern struct {
+ double time, acttime, bdist, bimp, bmin;
+ double ebeam; // lab-frame energy of projectile
+ double ecm;
+ } rsys_;
+
+ // defined in coms.f
+ extern struct {
+ nmaxIntArray spin, ncoll, charge, ityp, lstcoll, iso3, origin, strid, uid;
+ } isys_;
+
+ // defined in coor.f
+ extern struct {
+ nmaxDoubleArray r0, rx, ry, rz, p0, px, py, pz, fmass, rww, dectime;
+ } coor_;
+
+ // defined in inputs.f
+ extern struct {
+ int nevents;
+ std::array<int, 2> spityp; // particle codes of: [0]: projectile, [1]: target
+ int prspflg; // projectile special flag
+ int trspflg; // target special flag, set to 1 unless target is nucleus > H
+ std::array<int, 2> spiso3; // particle codes of: [0]: projectile, [1]: target
+ int outsteps, bflag, srtflag, efuncflag, nsrt, npb, firstev;
+ } inputs_;
+
+ // defined in inputs.f
+ extern struct {
+ double srtmin, srtmax, pbeam, betann, betatar, betapro, pbmin, pbmax;
+ } input2_;
+
+ // defined in options.f
+ extern struct {
+ std::array<double, constants::numcto> CTOption;
+ std::array<double, constants::numctp> CTParam;
+ } options_;
+
+ extern struct {
+ int fixedseed, bf13, bf14, bf15, bf16, bf17, bf18, bf19, bf20;
+ } loptions_;
+
+ // defined in urqmdInterface.F
+ extern struct { std::array<double, 3> xs, bim; } cxs_u2_;
+ }
+
+ /**
+ * convert CORSIKA code to UrQMD code tuple
+ *
+ * In the current implementation a detour via the PDG code is made.
+ */
+ std::pair<int, int> ConvertToUrQMD(particles::Code);
+ particles::Code ConvertFromUrQMD(int vItyp, int vIso3);
-struct {
- double time, acttime, bdist, bimp, bmin, ebeam, ecm;
-} rsys_;
-
-struct {
- int firstseed;
-} comseed_;
-
-struct {
- corsika::process::UrQMD::nmaxIntArray lsct;
- int logSky, logYuk, logCb, logPau;
-} logic_;
-
-struct {
- corsika::process::UrQMD::nmaxIntArray spin, ncoll, charge, ityp, lstcoll, iso3, origin, strid, uid;
-} isys_;
-}
+} // namespace corsika::process::UrQMD
#endif
diff --git a/Processes/UrQMD/urqmdInterface.F b/Processes/UrQMD/urqmdInterface.F
index f22e5d94d7d83f490dabf3f6433a98e918fe30db..8eaf1f740aab02052962e8d151c5a84d91bfd784 100644
--- a/Processes/UrQMD/urqmdInterface.F
+++ b/Processes/UrQMD/urqmdInterface.F
@@ -46,7 +46,8 @@ c local
integer iamaxu,idmaxu,iemaxu
common /cxs_u1/ sig_u1(mxie,mxid,mxia),iamaxu,idmaxu,iemaxu
double precision xs(3),bim(3)
- common /cxs_u2/ xs
+c M.R.: bim added to cxs_u2
+ common /cxs_u2/ xs,bim
integer iudebug
data bim/6.d0,6.d0,7.d0/
integer init
@@ -335,7 +336,7 @@ C DEFAULT SETTINGS FOR CTOption
CTOStrng(26)=' use z -> 1-z for diquark-pairs'
CTOption(27)=1 ! hjd1
CTOStrng(27)=' reference frame (1=target, 2=projectile, else=cms)'
- CTOption(28)=0
+ CTOption(28)=1 ! M.R. 2019-04-15
CTOStrng(28)=' propagate spectators also '
CTOption(29)=2
CTOStrng(29)=' no transverse momentum in clustr '