From 8337dd0ad2c8830e199b27ca6cdb062eda2943a0 Mon Sep 17 00:00:00 2001 From: Marco Giacalone Date: Mon, 28 Sep 2026 17:50:01 +0200 Subject: [PATCH] Implemented first external generator version of AMPT --- .../ampt/generator/PbPb_536TeV_AMPT.ampt | 151 +++++++++++ .../external/generator/generator_AMPT.C | 245 ++++++++++++++++++ .../common/ini/GeneratorAMPTPbPb536TeV.ini | 6 + .../ini/tests/GeneratorAMPTPbPb536TeV.C | 69 +++++ 4 files changed, 471 insertions(+) create mode 100644 MC/config/common/ampt/generator/PbPb_536TeV_AMPT.ampt create mode 100644 MC/config/common/external/generator/generator_AMPT.C create mode 100644 MC/config/common/ini/GeneratorAMPTPbPb536TeV.ini create mode 100644 MC/config/common/ini/tests/GeneratorAMPTPbPb536TeV.C diff --git a/MC/config/common/ampt/generator/PbPb_536TeV_AMPT.ampt b/MC/config/common/ampt/generator/PbPb_536TeV_AMPT.ampt new file mode 100644 index 000000000..31eb92c85 --- /dev/null +++ b/MC/config/common/ampt/generator/PbPb_536TeV_AMPT.ampt @@ -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) diff --git a/MC/config/common/external/generator/generator_AMPT.C b/MC/config/common/external/generator/generator_AMPT.C new file mode 100644 index 000000000..2de19ffeb --- /dev/null +++ b/MC/config/common/external/generator/generator_AMPT.C @@ -0,0 +1,245 @@ +#include +#include +#include +#include +#include +#include +#include +#include +#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(Key::generator, "ampt"); + eventHeader->SetB(mHeader.b); + eventHeader->putInfo(Key::impactParameter, mHeader.b); + eventHeader->putInfo(Key::nPart, mHeader.nPartProj + mHeader.nPartTarg); + eventHeader->putInfo(Key::nPartProjectile, mHeader.nPartProj); + eventHeader->putInfo(Key::nPartTarget, mHeader.nPartTarg); + eventHeader->putInfo(Key::planeAngle, mHeader.phiRP); + // Participant nucleons from elastic and inelastic collisions in projectile and target + eventHeader->putInfo("ampt_NELP", mHeader.nElP); + eventHeader->putInfo("ampt_NINP", mHeader.nInP); + eventHeader->putInfo("ampt_NELT", mHeader.nElT); + eventHeader->putInfo("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 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 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); +} diff --git a/MC/config/common/ini/GeneratorAMPTPbPb536TeV.ini b/MC/config/common/ini/GeneratorAMPTPbPb536TeV.ini new file mode 100644 index 000000000..2900caf2b --- /dev/null +++ b/MC/config/common/ini/GeneratorAMPTPbPb536TeV.ini @@ -0,0 +1,6 @@ +#NEV_TEST> 2 +### AMPT Pb-Pb 5.36 TeV (String Melting), converted directly to MCTracks without HepMC +### Use generateAMPT("/path/to/ampt.dat", "file") to read an existing AMPT output instead +[GeneratorExternal] +fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/common/external/generator/generator_AMPT.C +funcName=generateAMPT("${O2DPG_MC_CONFIG_ROOT}/MC/config/common/ampt/generator/PbPb_536TeV_AMPT.ampt") diff --git a/MC/config/common/ini/tests/GeneratorAMPTPbPb536TeV.C b/MC/config/common/ini/tests/GeneratorAMPTPbPb536TeV.C new file mode 100644 index 000000000..1d68673f4 --- /dev/null +++ b/MC/config/common/ini/tests/GeneratorAMPTPbPb536TeV.C @@ -0,0 +1,69 @@ +int External() +{ + std::string path{"o2sim_Kine.root"}; + // Check that file exists, can be opened and has the correct tree + TFile file(path.c_str(), "READ"); + if (file.IsZombie()) { + std::cerr << "Cannot open ROOT file " << path << "\n"; + return 1; + } + auto tree = (TTree*)file.Get("o2sim"); + if (!tree) { + std::cerr << "Cannot find tree o2sim in file " << path << "\n"; + return 1; + } + std::vector* tracks{}; + tree->SetBranchAddress("MCTrack", &tracks); + o2::dataformats::MCEventHeader* eventHeader = nullptr; + tree->SetBranchAddress("MCEventHeader.", &eventHeader); + + // Check if there are 2 events, as simulated in the o2dpg-test + auto nEvents = tree->GetEntries(); + if (nEvents != 2) { + std::cerr << "Expected 2 events, got " << nEvents << "\n"; + return 1; + } + // Beam momentum per nucleon for Pb-Pb at sqrt(s_NN) = 5.36 TeV + const double beamMomentum = std::sqrt(2680. * 2680. - 0.938 * 0.938); + for (Long64_t i = 0; i < nEvents; ++i) { + tree->GetEntry(i); + if (tracks->empty()) { + std::cerr << "Empty entry found at event " << i << "\n"; + return 1; + } + // Check the heavy-ion information of the header + bool isValid = false; + auto nPart = eventHeader->getInfo(o2::dataformats::MCInfoKeys::nPart, isValid); + if (!isValid || nPart <= 0) { + std::cerr << "Event " << i << " has no valid number of participants\n"; + return 1; + } + if (eventHeader->GetB() < 0. || eventHeader->GetB() > 20.) { + std::cerr << "Event " << i << " has impact parameter " << eventHeader->GetB() << " outside the [0, 20] fm range\n"; + return 1; + } + // Spectator nucleons (pT = 0, beam momentum) must be kept for both beams. + // AMPT writes momenta with 4 decimals, hence the tolerance + int nSpectatorsA = 0, nSpectatorsC = 0; + for (const auto& track : *tracks) { + if (track.GetPdgCode() != 2212 && track.GetPdgCode() != 2112) { + continue; + } + if (track.GetPt() == 0. && std::abs(std::abs(track.GetStartVertexMomentumZ()) - beamMomentum) < 1e-2) { + (track.GetStartVertexMomentumZ() > 0 ? nSpectatorsA : nSpectatorsC)++; + } + } + // The total number of nucleons in each nucleus must be conserved + if (nSpectatorsA + nSpectatorsC + nPart > 2 * 208) { + std::cerr << "Event " << i << " has more nucleons than the colliding nuclei\n"; + return 1; + } + if (nSpectatorsA == 0 || nSpectatorsC == 0) { + std::cerr << "Event " << i << " has no spectator nucleons (A side: " << nSpectatorsA << ", C side: " << nSpectatorsC << ")\n"; + return 1; + } + std::cout << "Event " << i << ": " << tracks->size() << " tracks, b = " << eventHeader->GetB() << " fm, Npart = " << nPart + << ", spectators A/C = " << nSpectatorsA << "/" << nSpectatorsC << "\n"; + } + return 0; +}