Skip to content

Commit 0509601

Browse files
jackal1-66jackal1-66
authored
Implemented first external generator version of AMPT (#2479)
1 parent e956bae commit 0509601

4 files changed

Lines changed: 471 additions & 0 deletions

File tree

‎MC/config/common/ampt/generator/PbPb_536TeV_AMPT.ampt‎

Lines changed: 151 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,151 @@
1+
5360 ! EFRM (sqrt(S_NN) in GeV if FRAME is CMS)
2+
CMS ! FRAME
3+
A ! PROJ
4+
A ! TARG
5+
208 ! IAP (projectile A number)
6+
82 ! IZP (projectile Z number)
7+
208 ! IAT (target A number)
8+
82 ! IZT (target Z number)
9+
1 ! NEVNT (total number of events) - overwritten by generator_AMPT.C
10+
0. ! BMIN (mininum impact parameter in fm)
11+
20. ! BMAX (maximum impact parameter in fm, also see below)
12+
4 ! ISOFT (D=4): select Default AMPT or String Melting(see below)
13+
1000 ! NTMAX: number of timesteps (D=150), see below
14+
0.2 ! DT: timestep in fm (hadron cascade time= DT*NTMAX) (D=0.2)
15+
0.30 ! PARJ(41): parameter a in Lund symmetric splitting function
16+
0.15 ! PARJ(42): parameter b in Lund symmetric splitting function
17+
1 ! (D=1,yes;0,no) flag for popcorn mechanism(netbaryon stopping)
18+
1.0 ! PARJ(5) to control BMBbar vs BBbar in popcorn (D=1.0)
19+
1 ! shadowing flag (Default=1,yes; 0,no)
20+
0 ! quenching flag (D=0,no; 1,yes)
21+
2.0 ! quenching parameter -dE/dx (GeV/fm) in case quenching flag=1
22+
2.0 ! p0 cutoff in HIJING for minijet productions (D=2.0)
23+
3.2264d0 ! parton screening mass in fm^(-1) (D=2.265d0), see below
24+
0 ! IZPC: (D=0 forward-angle parton scatterings; 100,isotropic)
25+
0.33d0 ! alpha in parton cascade (D=0.33d0), see parton screening mass
26+
1d6 ! dpcoal in GeV
27+
1d6 ! drcoal in fm
28+
0 ! ihjsed: take HIJING seed from below (D=0)or at runtime(11) - overwritten by generator_AMPT.C
29+
13150909 ! random seed for HIJING - overwritten by generator_AMPT.C
30+
8 ! random seed for parton cascade - overwritten by generator_AMPT.C
31+
0 ! flag for K0s weak decays (D=0,no; 1,yes)
32+
1 ! flag for phi decays at end of hadron cascade (D=1,yes; 0,no)
33+
0 ! flag for pi0 decays at end of hadron cascade (D=0,no; 1,yes)
34+
0 ! optional OSCAR output (D=0,no; 1,yes; 2&3,more parton info)
35+
0 ! flag for perturbative deuteron calculation (D=0,no; 1or2,yes)
36+
1 ! integer factor for perturbative deuterons(>=1 & <=10000)
37+
1 ! choice of cross section assumptions for deuteron reactions
38+
-7. ! Pt in GeV: generate events with >=1 minijet above this value
39+
1000 ! maxmiss (D=1000): maximum # of tries to repeat a HIJING event
40+
3 ! flag on initial and final state radiation (D=3,both yes; 0,no)
41+
1 ! flag on Kt kick (D=1,yes; 0,no)
42+
0 ! flag to turn on quark pair embedding (D=0,no; 1,yes)
43+
7., 0. ! Initial Px and Py values (GeV) of the embedded quark (u or d)
44+
0., 0. ! Initial x & y values (fm) of the embedded back-to-back q/qbar
45+
1, 5., 0. ! nsembd(D=0), psembd (in GeV),tmaxembd (in radian).
46+
0 ! Flag to enable users to modify shadowing (D=0,no; 1,yes)
47+
1.d0 ! Factor used to modify nuclear shadowing
48+
1 ! Flag for random orientation of reaction plane (D=0,no; 1,yes)
49+
50+
%%%%%%%%%% O2DPG notes:
51+
Pb-Pb at sqrt(s_NN) = 5.36 TeV, minimum bias (b = 0-20 fm), String Melting
52+
with the LHC settings of arXiv:1403.6321 (a=0.30, b=0.15/GeV^2, 1.5 mb parton
53+
cross section: alpha=0.33 and screening mass 3.2264/fm).
54+
Values are read by AMPT line by line: do NOT add or remove lines above.
55+
NEVNT (line 9), ihjsed (line 28) and the two seeds (lines 29-30) are
56+
overwritten by MC/config/common/external/generator/generator_AMPT.C.
57+
%%%%%%%%%% Further explanations:
58+
BMAX: the upper limit HIPR1(34)+HIPR1(35)=19.87fm (dAu), 25.60fm(AuAu).
59+
ISOFT: 1 Default,
60+
4 String Melting.
61+
PARJ(41) & (42): for string melting AMPT, 0.55 & 0.15/GeV^2 are recommended
62+
for top RHIC energies and 0.30 & 0.15/GeV^2 are recommended for
63+
LHC energies (see arXiv:1403.6321 for details).
64+
NTMAX: number of time-steps for hadron cascade.
65+
Use a large value (e.g. 1000) for LHC studies or HBT studies at RHIC.
66+
Using NTMAX=2 or 3 effectively turns off hadronic cascade.
67+
parton screening mass (in 1/fm): its square is inversely proportional to
68+
the parton cross section. Use D=2.265d0 for 3mb cross section
69+
when alpha in parton cascade is set to 0.33;
70+
(note: 3.2264d0 for 3mb cross section when alpha is set to 0.47).
71+
Using 1d4 effectively turns off parton cascade.
72+
ihjsed: if =11, take HIJING random seed at runtime so that
73+
every run may be automatically different (see file 'exec').
74+
iksdcy: flag for K0s weak decays for comparison with data.
75+
iphidcy: flag for phi meson decays at the end of hadron cascade for comparison
76+
with data; default is yes; use 0 to turn off these decays.
77+
Note: phi meson decay during hadron cascade is always enabled.
78+
ipi0dcy: flag for pi0 electromagnetic decays at the end of hadron cascade for
79+
comparison with data; set to 1 to turn on pi0 decays.
80+
ioscar: 0 Dafault,
81+
1 Write output in the OSCAR format,
82+
2 Write out the complete parton information
83+
(ana/parton-initial-afterPropagation.dat)
84+
right after string melting (before parton cascade),
85+
3 Write out several more files on parton information (see readme).
86+
idpert: flag for perturbative deuteron and antideuteron calculations
87+
with results in ana/ampt_pert.dat:
88+
0 No perturbative calculations,
89+
1 Trigger a production of NPERTD perturbative deuterons
90+
in each NN collision,
91+
2 Trigger a production of NPERTD perturbative deuterons only in
92+
an NN collision where a conventional deuteron is produced.
93+
Note: conventional deuteron calculations are always performed
94+
with results in ana/ampt.dat.
95+
NPERTD: number of perturbative deuterons produced in each triggered collision;
96+
setting it to 0 turns off perturbative deuteron productions.
97+
idxsec: choose a cross section model for deuteron inelastic/elastic collisions:
98+
1: same |matrix element|**2/s (after averaging over initial spins
99+
and isospins) for B+B -> deuteron+meson at the same sqrt(s);
100+
2: same |matrix element|**2/s for B+B -> deuteron+meson
101+
at the same sqrt(s)-threshold;
102+
3: same |matrix element|**2/s for deuteron+meson -> B+B
103+
at the same sqrt(s);
104+
4: same |matrix element|**2/s for deuteron+meson -> B+B
105+
at the same sqrt(s)-threshold;
106+
1 or 3 also chooses the same cross section for deuteron+meson or baryon
107+
elastic collision at the same sqrt(s);
108+
2 or 4 also chooses the same cross section for deuteron+meson or baryon
109+
elastic collision at the same sqrt(s)-threshold.
110+
%%%%%%%%%% For jet studies:
111+
pttrig: generate events with at least 1 initial minijet parton above this Pt
112+
value, otherwise repeat HIJING event until reaching maxmiss tries;
113+
use a negative value to disable this requirement and get normal events.
114+
maxmiss: maximum number of tries for the repetition of a HIJING event to obtain
115+
a minijet above the Pt value of pttrig; increase maxmiss if some events
116+
fail to generate at least 1 initial minijet parton above pttrig.
117+
it is safer to set a large value for high pttrig and/or large b value
118+
and/or smaller colliding nuclei.
119+
IHPR2(2): flag to turn off initial and final state radiation:
120+
0 both radiation off, 1 only final off, 2 only initial off, 3 both on.
121+
IHPR2(5): flag to turn off Pt kick due to soft interactions: 0 off, 1 on.
122+
Setting both IHPR2(2) and IHPR2(5) to zero makes it more likely to
123+
have two high-Pt minijet partons that are close to back-to-back.
124+
%%%%%%%%%% To embed a back-to-back light q/qbar jet pair
125+
%%%%%%%%%% and a given number of soft pions along each jet into each event:
126+
iembed: flag to turn on quark pair embedding:
127+
1: on with fixed position(xembd,pembd) and Pt(pxqembd,pyqembd);
128+
2: on with fixed position(xembd,pembd) and random azimuthal angle
129+
with Pt-magnitude given by sqrt(pxqembd^2+pyqembd^2);
130+
3: on with random position and fixed Pt(pxqembd,pyqembd);
131+
4: on with random position and random random azimuthal angle
132+
with Pt-magnitude given by sqrt(pxqembd^2+pyqembd^2);
133+
for iembed=3 or 4: need a position file "embed-jet-xy.txt";
134+
Other integers: off.
135+
pxqembd, pyqembd: sqrt(pxqembd^2+pyqembd^2) > 70MeV/c is required;
136+
the embedded quark and antiquark have pz=0.
137+
xembd, yembd: the embedded quark and antiquark jets have z=0 initially. Note:
138+
the x-axis is defined as the direction along the impact parameter.
139+
nsembd: number of soft pions to be embedded with each high-Pt parton
140+
in the embedded jet pair.
141+
psembd: Momentum of each embedded soft pion in GeV.
142+
tmaxembd: maximum angle(rad) of embedded soft pions relative to high-Pt parton.
143+
%%%%%%%%%% User modification of nuclear shadowing:
144+
ishadow: set to 1 to enable users to adjust nuclear shadowing
145+
provided the shadowing flag IHPR2(6) is turned on; default value is 0.
146+
dshadow: valid when ishadow=1; this parameter modifies the HIJING shadowing
147+
parameterization Ra(x,r)==1+fa(x,r) via Ra(x,r)==1+fa(x,r)*dshadow,
148+
so the value of 0.d0 turns off shadowing
149+
and the value of 1.d0 uses the default HIJING shadowing;
150+
currently limited to 0.d0<=dshadow<=1.d0 to make sure Ra(x,r)>0.
151+
iphirp: set to 1 to turn on random orientation of reaction plane (D=0)

‎MC/config/common/external/generator/generator_AMPT.C‎

Lines changed: 245 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,245 @@
1+
#include <cmath>
2+
#include <filesystem>
3+
#include <fstream>
4+
#include <map>
5+
#include <sstream>
6+
#include <string>
7+
#include <vector>
8+
#include <fairlogger/Logger.h>
9+
#include "TRandom.h"
10+
#include "TParticle.h"
11+
#include "TString.h"
12+
#include "Generators/Generator.h"
13+
#include "Generators/GeneratorFileOrCmd.h"
14+
#include "CommonUtils/FileSystemUtils.h"
15+
#include "SimulationDataFormat/MCEventHeader.h"
16+
#include "SimulationDataFormat/MCUtils.h"
17+
18+
/// AMPT adapted external generator
19+
///
20+
/// @author Marco Giacalone (marco.giacalone@cern.ch)
21+
/// @date 09/26
22+
// The AMPT output (ana/ampt.dat) is converted directly into TParticles, without any HepMC step, as done in the past.
23+
// Two modes are available:
24+
// - "run" : the path is an AMPT input card (input.ampt). AMPT is started in the background in a dedicated
25+
// working directory, where ana/ampt.dat is a named pipe, so events are streamed to O2 as soon as
26+
// they are produced. NEVNT and the random seeds of the card are replaced by the macro.
27+
// - "file" : the path is an existing ampt.dat file, whose events are read sequentially.
28+
// All the particles are transported, including the spectator nucleons (needed by the ZDC).
29+
// This could be changed in case it's actually not needed. To verify
30+
31+
// o2-sim -g external --noGeant -n 2 --configFile ${O2DPG_MC_CONFIG_ROOT}/MC/config/common/ini/GeneratorAMPTPbPb536TeV.ini
32+
// or, reading an existing AMPT output
33+
// o2-sim -g external --noGeant -n 2 --configKeyValues "GeneratorExternal.fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/common/external/generator/generator_AMPT.C;GeneratorExternal.funcName=generateAMPT(\"/path/to/ampt.dat\",\"file\")"
34+
35+
namespace o2
36+
{
37+
namespace eventgen
38+
{
39+
40+
// Inheriting from GeneratorFileOrCmd as well to use the pipe mechanism
41+
class GeneratorAMPT : public Generator, public GeneratorFileOrCmd
42+
{
43+
public:
44+
GeneratorAMPT(const std::string& path, bool runAMPT, unsigned int nEvents) : mPath(path)
45+
{
46+
setNEvents(nEvents);
47+
if (runAMPT) {
48+
// AMPT does NSEED = 2 * NSEED + 1 internally, hence the seed must stay below 2^30
49+
setSeed(gRandom->Integer(536870911) + 1);
50+
mWorkDir = Form("ampt_%lu", mSeed);
51+
// AMPT always reads a number from stdin (used as seed only when ihjsed = 11). If AMPT fails,
52+
// the pipe is opened and closed by the shell, so the reader gets an end of file instead of hanging.
53+
// The timeout avoids a shell blocked forever on the pipe when the reader is gone (this was seen if o2-sim crashed mid-run for some reason)
54+
setCmd("cd " + mWorkDir + " && { echo 0 | \"$AMPT_ROOT/bin/ampt\" > ampt.log 2>&1 || timeout 60 sh -c ': > ana/ampt.dat'; }");
55+
}
56+
}
57+
~GeneratorAMPT() override
58+
{
59+
stop();
60+
if (not mCmd.empty()) {
61+
removeTemp();
62+
}
63+
}
64+
65+
Bool_t Init() override
66+
{
67+
if (mCmd.empty()) {
68+
mInput.open(mPath);
69+
} else {
70+
if (!getenv("AMPT_ROOT")) {
71+
LOG(fatal) << "AMPT_ROOT is not set, load the AMPT package";
72+
return false;
73+
}
74+
if (!prepareWorkDir() || !makeFifo() || !executeCmdLine(mCmd)) {
75+
return false;
76+
}
77+
LOG(info) << "AMPT started in " << mWorkDir << " with seed " << mSeed << " for " << mNEvents << " events";
78+
mInput.open(mTemporary); // blocks until AMPT opens the pipe for writing
79+
}
80+
if (!mInput.is_open()) {
81+
LOG(fatal) << "Cannot open AMPT output " << (mCmd.empty() ? mPath : mTemporary);
82+
return false;
83+
}
84+
return Generator::Init();
85+
}
86+
87+
// Reads the next event from ampt.dat: one header line with 11 fields followed by
88+
// exactly nParticles lines with 9 fields (pdg px py pz m x y z t)
89+
Bool_t generateEvent() override
90+
{
91+
std::string line;
92+
if (!nextLine(line)) {
93+
LOG(fatal) << "AMPT output ended after " << mEventCounter << " events (" << mNEvents << " requested)"
94+
<< (mCmd.empty() ? "" : ", see " + mWorkDir + "/ampt.log");
95+
return false;
96+
}
97+
std::istringstream header(line);
98+
header >> mHeader.event >> mHeader.run >> mHeader.nParticles >> mHeader.b >> mHeader.nPartProj >> mHeader.nPartTarg >> mHeader.nElP >> mHeader.nInP >> mHeader.nElT >> mHeader.nInT >> mHeader.phiRP;
99+
if (header.fail() || mHeader.nParticles < 0) {
100+
LOG(fatal) << "Malformed AMPT event header: " << line;
101+
return false;
102+
}
103+
mTracks.clear();
104+
mTracks.reserve(mHeader.nParticles);
105+
for (int i = 0; i < mHeader.nParticles; ++i) {
106+
AMPTTrack t;
107+
if (!nextLine(line)) {
108+
LOG(fatal) << "AMPT event " << mHeader.event << " is truncated: " << i << " out of " << mHeader.nParticles << " particles";
109+
return false;
110+
}
111+
std::istringstream particle(line);
112+
particle >> t.pdg >> t.px >> t.py >> t.pz >> t.m;
113+
if (particle.fail()) {
114+
LOG(fatal) << "Malformed AMPT particle line: " << line;
115+
return false;
116+
}
117+
mTracks.push_back(t);
118+
}
119+
mEventCounter++;
120+
LOG(info) << "AMPT event " << mEventCounter << "/" << mNEvents << ": " << mHeader.nParticles << " particles, b = " << mHeader.b << " fm";
121+
return true;
122+
}
123+
124+
Bool_t importParticles() override
125+
{
126+
mParticles.clear();
127+
for (const auto& t : mTracks) {
128+
// AMPT uses PDG codes, apart from (anti)deuterons
129+
int pdg = std::abs(t.pdg) == 42 ? (t.pdg > 0 ? 1 : -1) * 1000010020 : t.pdg;
130+
double e = std::sqrt(t.px * t.px + t.py * t.py + t.pz * t.pz + t.m * t.m);
131+
// Freeze-out coordinates are in fm, hence negligible: particles are produced at the interaction vertex
132+
TParticle particle(pdg, 1, -1, -1, -1, -1, t.px, t.py, t.pz, e, 0., 0., 0., 0.);
133+
o2::mcutils::MCGenHelper::encodeParticleStatusAndTracking(particle, true);
134+
mParticles.push_back(particle);
135+
}
136+
return true;
137+
}
138+
139+
void updateHeader(o2::dataformats::MCEventHeader* eventHeader) override
140+
{
141+
using Key = o2::dataformats::MCInfoKeys;
142+
eventHeader->putInfo<std::string>(Key::generator, "ampt");
143+
eventHeader->SetB(mHeader.b);
144+
eventHeader->putInfo<float>(Key::impactParameter, mHeader.b);
145+
eventHeader->putInfo<int>(Key::nPart, mHeader.nPartProj + mHeader.nPartTarg);
146+
eventHeader->putInfo<int>(Key::nPartProjectile, mHeader.nPartProj);
147+
eventHeader->putInfo<int>(Key::nPartTarget, mHeader.nPartTarg);
148+
eventHeader->putInfo<double>(Key::planeAngle, mHeader.phiRP);
149+
// Participant nucleons from elastic and inelastic collisions in projectile and target
150+
eventHeader->putInfo<int>("ampt_NELP", mHeader.nElP);
151+
eventHeader->putInfo<int>("ampt_NINP", mHeader.nInP);
152+
eventHeader->putInfo<int>("ampt_NELT", mHeader.nElT);
153+
eventHeader->putInfo<int>("ampt_NINTHJ", mHeader.nInT);
154+
}
155+
156+
void stop() override
157+
{
158+
mInput.close();
159+
if (not mCmd.empty()) {
160+
// AMPT exits by itself once all the events are written, otherwise it is terminated
161+
terminateCmd(sStopGraceMillis);
162+
}
163+
}
164+
165+
private:
166+
struct AMPTHeader {
167+
int event = 0, run = 0, nParticles = 0;
168+
double b = 0.;
169+
int nPartProj = 0, nPartTarg = 0, nElP = 0, nInP = 0, nElT = 0, nInT = 0;
170+
double phiRP = 0.;
171+
};
172+
struct AMPTTrack {
173+
int pdg = 0;
174+
double px = 0., py = 0., pz = 0., m = 0.;
175+
};
176+
177+
// Creates the working directory with the edited card. ana/ampt.dat becomes the named pipe, while
178+
// the outputs growing with the number of events and not needed here points to /dev/null
179+
bool prepareWorkDir()
180+
{
181+
std::filesystem::create_directories(mWorkDir + "/ana");
182+
// AMPT reads the card sequentially, one value per line: lines are replaced by their position
183+
const std::map<int, std::string> replacements = {
184+
{9, std::to_string(mNEvents) + "\t\t! NEVNT (set by generator_AMPT.C)"},
185+
{28, "0\t\t! ihjsed (set by generator_AMPT.C)"},
186+
{29, std::to_string(mSeed) + "\t\t! random seed for HIJING (set by generator_AMPT.C)"},
187+
{30, std::to_string(gRandom->Integer(536870911) + 1) + "\t\t! random seed for parton cascade (set by generator_AMPT.C)"}};
188+
std::ifstream src(mPath);
189+
std::ofstream dst(mWorkDir + "/input.ampt");
190+
std::string line;
191+
for (int lineNumber = 1; std::getline(src, line); ++lineNumber) {
192+
auto replacement = replacements.find(lineNumber);
193+
dst << (replacement != replacements.end() ? replacement->second : line) << "\n";
194+
}
195+
for (const auto& file : {"zpc.dat", "npart-xy.dat"}) {
196+
std::filesystem::remove(mWorkDir + "/ana/" + file);
197+
std::filesystem::create_symlink("/dev/null", mWorkDir + "/ana/" + file);
198+
}
199+
// Named pipe used by makeFifo and removed by removeTemp
200+
mTemporary = mWorkDir + "/ana/ampt.dat";
201+
mFileNames = {mTemporary};
202+
return true;
203+
}
204+
205+
bool nextLine(std::string& line)
206+
{
207+
while (std::getline(mInput, line)) {
208+
if (line.find_first_not_of(" \t\r") != std::string::npos) {
209+
return true;
210+
}
211+
}
212+
return false;
213+
}
214+
215+
std::string mPath;
216+
std::string mWorkDir;
217+
std::ifstream mInput;
218+
unsigned int mEventCounter = 0;
219+
AMPTHeader mHeader;
220+
std::vector<AMPTTrack> mTracks;
221+
};
222+
223+
} // namespace eventgen
224+
} // namespace o2
225+
226+
// path : AMPT input card (mode "run") or existing ampt.dat file (mode "file")
227+
// maxEvents : number of events for AMPT when the total number of events is not known (e.g. hyperloop)
228+
FairGenerator* generateAMPT(std::string path, std::string mode = "run", int maxEvents = 2147483647)
229+
{
230+
path = o2::utils::expandShellVarsInFileName(path);
231+
if (mode != "run" && mode != "file") {
232+
LOG(fatal) << "Unknown AMPT generator mode " << mode << ", use \"run\" or \"file\"";
233+
return nullptr;
234+
}
235+
if (!std::filesystem::exists(path)) {
236+
LOG(fatal) << "AMPT " << (mode == "run" ? "input card " : "output file ") << path << " does not exist";
237+
return nullptr;
238+
}
239+
unsigned int nEvents = o2::eventgen::Generator::getTotalNEvents();
240+
if (nEvents == 0) {
241+
nEvents = maxEvents;
242+
}
243+
LOG(info) << "AMPT generator in \"" << mode << "\" mode using " << path;
244+
return new o2::eventgen::GeneratorAMPT(path, mode == "run", nEvents);
245+
}

0 commit comments

Comments
 (0)