Skip to content

Implemented first external generator version of AMPT #2479

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
jackal1-66 merged 1 commit into from
Sep 30, 2026
Merged

Implemented first external generator version of AMPT #2479

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
151 changes: 151 additions & 0 deletions MC/config/common/ampt/generator/PbPb_536TeV_AMPT.ampt
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,151 @@
5360 ! EFRM (sqrt(S_NN) in GeV if FRAME is CMS)
CMS ! FRAME
A ! PROJ
A ! TARG
208 ! IAP (projectile A number)
82 ! IZP (projectile Z number)
208 ! IAT (target A number)
82 ! IZT (target Z number)
1 ! NEVNT (total number of events) - overwritten by generator_AMPT.C
0. ! BMIN (mininum impact parameter in fm)
20. ! BMAX (maximum impact parameter in fm, also see below)
4 ! ISOFT (D=4): select Default AMPT or String Melting(see below)
1000 ! NTMAX: number of timesteps (D=150), see below
0.2 ! DT: timestep in fm (hadron cascade time= DT*NTMAX) (D=0.2)
0.30 ! PARJ(41): parameter a in Lund symmetric splitting function
0.15 ! PARJ(42): parameter b in Lund symmetric splitting function
1 ! (D=1,yes;0,no) flag for popcorn mechanism(netbaryon stopping)
1.0 ! PARJ(5) to control BMBbar vs BBbar in popcorn (D=1.0)
1 ! shadowing flag (Default=1,yes; 0,no)
0 ! quenching flag (D=0,no; 1,yes)
2.0 ! quenching parameter -dE/dx (GeV/fm) in case quenching flag=1
2.0 ! p0 cutoff in HIJING for minijet productions (D=2.0)
3.2264d0 ! parton screening mass in fm^(-1) (D=2.265d0), see below
0 ! IZPC: (D=0 forward-angle parton scatterings; 100,isotropic)
0.33d0 ! alpha in parton cascade (D=0.33d0), see parton screening mass
1d6 ! dpcoal in GeV
1d6 ! drcoal in fm
0 ! ihjsed: take HIJING seed from below (D=0)or at runtime(11) - overwritten by generator_AMPT.C
13150909 ! random seed for HIJING - overwritten by generator_AMPT.C
8 ! random seed for parton cascade - overwritten by generator_AMPT.C
0 ! flag for K0s weak decays (D=0,no; 1,yes)
1 ! flag for phi decays at end of hadron cascade (D=1,yes; 0,no)
0 ! flag for pi0 decays at end of hadron cascade (D=0,no; 1,yes)
0 ! optional OSCAR output (D=0,no; 1,yes; 2&3,more parton info)
0 ! flag for perturbative deuteron calculation (D=0,no; 1or2,yes)
1 ! integer factor for perturbative deuterons(>=1 & <=10000)
1 ! choice of cross section assumptions for deuteron reactions
-7. ! Pt in GeV: generate events with >=1 minijet above this value
1000 ! maxmiss (D=1000): maximum # of tries to repeat a HIJING event
3 ! flag on initial and final state radiation (D=3,both yes; 0,no)
1 ! flag on Kt kick (D=1,yes; 0,no)
0 ! flag to turn on quark pair embedding (D=0,no; 1,yes)
7., 0. ! Initial Px and Py values (GeV) of the embedded quark (u or d)
0., 0. ! Initial x & y values (fm) of the embedded back-to-back q/qbar
1, 5., 0. ! nsembd(D=0), psembd (in GeV),tmaxembd (in radian).
0 ! Flag to enable users to modify shadowing (D=0,no; 1,yes)
1.d0 ! Factor used to modify nuclear shadowing
1 ! Flag for random orientation of reaction plane (D=0,no; 1,yes)

%%%%%%%%%% O2DPG notes:
Pb-Pb at sqrt(s_NN) = 5.36 TeV, minimum bias (b = 0-20 fm), String Melting
with the LHC settings of arXiv:1403.6321 (a=0.30, b=0.15/GeV^2, 1.5 mb parton
cross section: alpha=0.33 and screening mass 3.2264/fm).
Values are read by AMPT line by line: do NOT add or remove lines above.
NEVNT (line 9), ihjsed (line 28) and the two seeds (lines 29-30) are
overwritten by MC/config/common/external/generator/generator_AMPT.C.
%%%%%%%%%% Further explanations:
BMAX: the upper limit HIPR1(34)+HIPR1(35)=19.87fm (dAu), 25.60fm(AuAu).
ISOFT: 1 Default,
4 String Melting.
PARJ(41) & (42): for string melting AMPT, 0.55 & 0.15/GeV^2 are recommended
for top RHIC energies and 0.30 & 0.15/GeV^2 are recommended for
LHC energies (see arXiv:1403.6321 for details).
NTMAX: number of time-steps for hadron cascade.
Use a large value (e.g. 1000) for LHC studies or HBT studies at RHIC.
Using NTMAX=2 or 3 effectively turns off hadronic cascade.
parton screening mass (in 1/fm): its square is inversely proportional to
the parton cross section. Use D=2.265d0 for 3mb cross section
when alpha in parton cascade is set to 0.33;
(note: 3.2264d0 for 3mb cross section when alpha is set to 0.47).
Using 1d4 effectively turns off parton cascade.
ihjsed: if =11, take HIJING random seed at runtime so that
every run may be automatically different (see file 'exec').
iksdcy: flag for K0s weak decays for comparison with data.
iphidcy: flag for phi meson decays at the end of hadron cascade for comparison
with data; default is yes; use 0 to turn off these decays.
Note: phi meson decay during hadron cascade is always enabled.
ipi0dcy: flag for pi0 electromagnetic decays at the end of hadron cascade for
comparison with data; set to 1 to turn on pi0 decays.
ioscar: 0 Dafault,
1 Write output in the OSCAR format,
2 Write out the complete parton information
(ana/parton-initial-afterPropagation.dat)
right after string melting (before parton cascade),
3 Write out several more files on parton information (see readme).
idpert: flag for perturbative deuteron and antideuteron calculations
with results in ana/ampt_pert.dat:
0 No perturbative calculations,
1 Trigger a production of NPERTD perturbative deuterons
in each NN collision,
2 Trigger a production of NPERTD perturbative deuterons only in
an NN collision where a conventional deuteron is produced.
Note: conventional deuteron calculations are always performed
with results in ana/ampt.dat.
NPERTD: number of perturbative deuterons produced in each triggered collision;
setting it to 0 turns off perturbative deuteron productions.
idxsec: choose a cross section model for deuteron inelastic/elastic collisions:
1: same |matrix element|**2/s (after averaging over initial spins
and isospins) for B+B -> deuteron+meson at the same sqrt(s);
2: same |matrix element|**2/s for B+B -> deuteron+meson
at the same sqrt(s)-threshold;
3: same |matrix element|**2/s for deuteron+meson -> B+B
at the same sqrt(s);
4: same |matrix element|**2/s for deuteron+meson -> B+B
at the same sqrt(s)-threshold;
1 or 3 also chooses the same cross section for deuteron+meson or baryon
elastic collision at the same sqrt(s);
2 or 4 also chooses the same cross section for deuteron+meson or baryon
elastic collision at the same sqrt(s)-threshold.
%%%%%%%%%% For jet studies:
pttrig: generate events with at least 1 initial minijet parton above this Pt
value, otherwise repeat HIJING event until reaching maxmiss tries;
use a negative value to disable this requirement and get normal events.
maxmiss: maximum number of tries for the repetition of a HIJING event to obtain
a minijet above the Pt value of pttrig; increase maxmiss if some events
fail to generate at least 1 initial minijet parton above pttrig.
it is safer to set a large value for high pttrig and/or large b value
and/or smaller colliding nuclei.
IHPR2(2): flag to turn off initial and final state radiation:
0 both radiation off, 1 only final off, 2 only initial off, 3 both on.
IHPR2(5): flag to turn off Pt kick due to soft interactions: 0 off, 1 on.
Setting both IHPR2(2) and IHPR2(5) to zero makes it more likely to
have two high-Pt minijet partons that are close to back-to-back.
%%%%%%%%%% To embed a back-to-back light q/qbar jet pair
%%%%%%%%%% and a given number of soft pions along each jet into each event:
iembed: flag to turn on quark pair embedding:
1: on with fixed position(xembd,pembd) and Pt(pxqembd,pyqembd);
2: on with fixed position(xembd,pembd) and random azimuthal angle
with Pt-magnitude given by sqrt(pxqembd^2+pyqembd^2);
3: on with random position and fixed Pt(pxqembd,pyqembd);
4: on with random position and random random azimuthal angle
with Pt-magnitude given by sqrt(pxqembd^2+pyqembd^2);
for iembed=3 or 4: need a position file "embed-jet-xy.txt";
Other integers: off.
pxqembd, pyqembd: sqrt(pxqembd^2+pyqembd^2) > 70MeV/c is required;
the embedded quark and antiquark have pz=0.
xembd, yembd: the embedded quark and antiquark jets have z=0 initially. Note:
the x-axis is defined as the direction along the impact parameter.
nsembd: number of soft pions to be embedded with each high-Pt parton
in the embedded jet pair.
psembd: Momentum of each embedded soft pion in GeV.
tmaxembd: maximum angle(rad) of embedded soft pions relative to high-Pt parton.
%%%%%%%%%% User modification of nuclear shadowing:
ishadow: set to 1 to enable users to adjust nuclear shadowing
provided the shadowing flag IHPR2(6) is turned on; default value is 0.
dshadow: valid when ishadow=1; this parameter modifies the HIJING shadowing
parameterization Ra(x,r)==1+fa(x,r) via Ra(x,r)==1+fa(x,r)*dshadow,
so the value of 0.d0 turns off shadowing
and the value of 1.d0 uses the default HIJING shadowing;
currently limited to 0.d0<=dshadow<=1.d0 to make sure Ra(x,r)>0.
iphirp: set to 1 to turn on random orientation of reaction plane (D=0)
245 changes: 245 additions & 0 deletions MC/config/common/external/generator/generator_AMPT.C
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,245 @@
#include <cmath>
#include <filesystem>
#include <fstream>
#include <map>
#include <sstream>
#include <string>
#include <vector>
#include <fairlogger/Logger.h>
#include "TRandom.h"
#include "TParticle.h"
#include "TString.h"
#include "Generators/Generator.h"
#include "Generators/GeneratorFileOrCmd.h"
#include "CommonUtils/FileSystemUtils.h"
#include "SimulationDataFormat/MCEventHeader.h"
#include "SimulationDataFormat/MCUtils.h"

/// AMPT adapted external generator
///
/// @author Marco Giacalone (marco.giacalone@cern.ch)
/// @date 09/26
// The AMPT output (ana/ampt.dat) is converted directly into TParticles, without any HepMC step, as done in the past.
// Two modes are available:
// - "run" : the path is an AMPT input card (input.ampt). AMPT is started in the background in a dedicated
// working directory, where ana/ampt.dat is a named pipe, so events are streamed to O2 as soon as
// they are produced. NEVNT and the random seeds of the card are replaced by the macro.
// - "file" : the path is an existing ampt.dat file, whose events are read sequentially.
// All the particles are transported, including the spectator nucleons (needed by the ZDC).
// This could be changed in case it's actually not needed. To verify

// o2-sim -g external --noGeant -n 2 --configFile ${O2DPG_MC_CONFIG_ROOT}/MC/config/common/ini/GeneratorAMPTPbPb536TeV.ini
// or, reading an existing AMPT output
// 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\")"

namespace o2
{
namespace eventgen
{

// Inheriting from GeneratorFileOrCmd as well to use the pipe mechanism
class GeneratorAMPT : public Generator, public GeneratorFileOrCmd
{
public:
GeneratorAMPT(const std::string& path, bool runAMPT, unsigned int nEvents) : mPath(path)
{
setNEvents(nEvents);
if (runAMPT) {
// AMPT does NSEED = 2 * NSEED + 1 internally, hence the seed must stay below 2^30
setSeed(gRandom->Integer(536870911) + 1);
mWorkDir = Form("ampt_%lu", mSeed);
// AMPT always reads a number from stdin (used as seed only when ihjsed = 11). If AMPT fails,
// the pipe is opened and closed by the shell, so the reader gets an end of file instead of hanging.
// 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)
setCmd("cd " + mWorkDir + " && { echo 0 | \"$AMPT_ROOT/bin/ampt\" > ampt.log 2>&1 || timeout 60 sh -c ': > ana/ampt.dat'; }");
}
}
~GeneratorAMPT() override
{
stop();
if (not mCmd.empty()) {
removeTemp();
}
}

Bool_t Init() override
{
if (mCmd.empty()) {
mInput.open(mPath);
} else {
if (!getenv("AMPT_ROOT")) {
LOG(fatal) << "AMPT_ROOT is not set, load the AMPT package";
return false;
}
if (!prepareWorkDir() || !makeFifo() || !executeCmdLine(mCmd)) {
return false;
}
LOG(info) << "AMPT started in " << mWorkDir << " with seed " << mSeed << " for " << mNEvents << " events";
mInput.open(mTemporary); // blocks until AMPT opens the pipe for writing
}
if (!mInput.is_open()) {
LOG(fatal) << "Cannot open AMPT output " << (mCmd.empty() ? mPath : mTemporary);
return false;
}
return Generator::Init();
}

// Reads the next event from ampt.dat: one header line with 11 fields followed by
// exactly nParticles lines with 9 fields (pdg px py pz m x y z t)
Bool_t generateEvent() override
{
std::string line;
if (!nextLine(line)) {
LOG(fatal) << "AMPT output ended after " << mEventCounter << " events (" << mNEvents << " requested)"
<< (mCmd.empty() ? "" : ", see " + mWorkDir + "/ampt.log");
return false;
}
std::istringstream header(line);
header >> mHeader.event >> mHeader.run >> mHeader.nParticles >> mHeader.b >> mHeader.nPartProj >> mHeader.nPartTarg >> mHeader.nElP >> mHeader.nInP >> mHeader.nElT >> mHeader.nInT >> mHeader.phiRP;
if (header.fail() || mHeader.nParticles < 0) {
LOG(fatal) << "Malformed AMPT event header: " << line;
return false;
}
mTracks.clear();
mTracks.reserve(mHeader.nParticles);
for (int i = 0; i < mHeader.nParticles; ++i) {
AMPTTrack t;
if (!nextLine(line)) {
LOG(fatal) << "AMPT event " << mHeader.event << " is truncated: " << i << " out of " << mHeader.nParticles << " particles";
return false;
}
std::istringstream particle(line);
particle >> t.pdg >> t.px >> t.py >> t.pz >> t.m;
if (particle.fail()) {
LOG(fatal) << "Malformed AMPT particle line: " << line;
return false;
}
mTracks.push_back(t);
}
mEventCounter++;
LOG(info) << "AMPT event " << mEventCounter << "/" << mNEvents << ": " << mHeader.nParticles << " particles, b = " << mHeader.b << " fm";
return true;
}

Bool_t importParticles() override
{
mParticles.clear();
for (const auto& t : mTracks) {
// AMPT uses PDG codes, apart from (anti)deuterons
int pdg = std::abs(t.pdg) == 42 ? (t.pdg > 0 ? 1 : -1) * 1000010020 : t.pdg;
double e = std::sqrt(t.px * t.px + t.py * t.py + t.pz * t.pz + t.m * t.m);
// Freeze-out coordinates are in fm, hence negligible: particles are produced at the interaction vertex
TParticle particle(pdg, 1, -1, -1, -1, -1, t.px, t.py, t.pz, e, 0., 0., 0., 0.);
o2::mcutils::MCGenHelper::encodeParticleStatusAndTracking(particle, true);
mParticles.push_back(particle);
}
return true;
}

void updateHeader(o2::dataformats::MCEventHeader* eventHeader) override
{
using Key = o2::dataformats::MCInfoKeys;
eventHeader->putInfo<std::string>(Key::generator, "ampt");
eventHeader->SetB(mHeader.b);
eventHeader->putInfo<float>(Key::impactParameter, mHeader.b);
eventHeader->putInfo<int>(Key::nPart, mHeader.nPartProj + mHeader.nPartTarg);
eventHeader->putInfo<int>(Key::nPartProjectile, mHeader.nPartProj);
eventHeader->putInfo<int>(Key::nPartTarget, mHeader.nPartTarg);
eventHeader->putInfo<double>(Key::planeAngle, mHeader.phiRP);
// Participant nucleons from elastic and inelastic collisions in projectile and target
eventHeader->putInfo<int>("ampt_NELP", mHeader.nElP);
eventHeader->putInfo<int>("ampt_NINP", mHeader.nInP);
eventHeader->putInfo<int>("ampt_NELT", mHeader.nElT);
eventHeader->putInfo<int>("ampt_NINTHJ", mHeader.nInT);
}

void stop() override
{
mInput.close();
if (not mCmd.empty()) {
// AMPT exits by itself once all the events are written, otherwise it is terminated
terminateCmd(sStopGraceMillis);
}
}

private:
struct AMPTHeader {
int event = 0, run = 0, nParticles = 0;
double b = 0.;
int nPartProj = 0, nPartTarg = 0, nElP = 0, nInP = 0, nElT = 0, nInT = 0;
double phiRP = 0.;
};
struct AMPTTrack {
int pdg = 0;
double px = 0., py = 0., pz = 0., m = 0.;
};

// Creates the working directory with the edited card. ana/ampt.dat becomes the named pipe, while
// the outputs growing with the number of events and not needed here points to /dev/null
bool prepareWorkDir()
{
std::filesystem::create_directories(mWorkDir + "/ana");
// AMPT reads the card sequentially, one value per line: lines are replaced by their position
const std::map<int, std::string> replacements = {
{9, std::to_string(mNEvents) + "\t\t! NEVNT (set by generator_AMPT.C)"},
{28, "0\t\t! ihjsed (set by generator_AMPT.C)"},
{29, std::to_string(mSeed) + "\t\t! random seed for HIJING (set by generator_AMPT.C)"},
{30, std::to_string(gRandom->Integer(536870911) + 1) + "\t\t! random seed for parton cascade (set by generator_AMPT.C)"}};
std::ifstream src(mPath);
std::ofstream dst(mWorkDir + "/input.ampt");
std::string line;
for (int lineNumber = 1; std::getline(src, line); ++lineNumber) {
auto replacement = replacements.find(lineNumber);
dst << (replacement != replacements.end() ? replacement->second : line) << "\n";
}
for (const auto& file : {"zpc.dat", "npart-xy.dat"}) {
std::filesystem::remove(mWorkDir + "/ana/" + file);
std::filesystem::create_symlink("/dev/null", mWorkDir + "/ana/" + file);
}
// Named pipe used by makeFifo and removed by removeTemp
mTemporary = mWorkDir + "/ana/ampt.dat";
mFileNames = {mTemporary};
return true;
}

bool nextLine(std::string& line)
{
while (std::getline(mInput, line)) {
if (line.find_first_not_of(" \t\r") != std::string::npos) {
return true;
}
}
return false;
}

std::string mPath;
std::string mWorkDir;
std::ifstream mInput;
unsigned int mEventCounter = 0;
AMPTHeader mHeader;
std::vector<AMPTTrack> mTracks;
};

} // namespace eventgen
} // namespace o2

// path : AMPT input card (mode "run") or existing ampt.dat file (mode "file")
// maxEvents : number of events for AMPT when the total number of events is not known (e.g. hyperloop)
FairGenerator* generateAMPT(std::string path, std::string mode = "run", int maxEvents = 2147483647)
{
path = o2::utils::expandShellVarsInFileName(path);
if (mode != "run" && mode != "file") {
LOG(fatal) << "Unknown AMPT generator mode " << mode << ", use \"run\" or \"file\"";
return nullptr;
}
if (!std::filesystem::exists(path)) {
LOG(fatal) << "AMPT " << (mode == "run" ? "input card " : "output file ") << path << " does not exist";
return nullptr;
}
unsigned int nEvents = o2::eventgen::Generator::getTotalNEvents();
if (nEvents == 0) {
nEvents = maxEvents;
}
LOG(info) << "AMPT generator in \"" << mode << "\" mode using " << path;
return new o2::eventgen::GeneratorAMPT(path, mode == "run", nEvents);
}
Loading
Loading