From f0c363ef736de4f124ea6af505e72767da053e27 Mon Sep 17 00:00:00 2001 From: Marco Giacalone Date: Mon, 28 Sep 2026 22:57:12 +0200 Subject: [PATCH] First implementation of HGPythia --- .../external/generator/GeneratorHGPythia8.C | 460 ++++++++++++++++++ .../ini/GeneratorHGPythia8PbPb536TeV.ini | 8 + .../ini/tests/GeneratorHGPythia8PbPb536TeV.C | 56 +++ 3 files changed, 524 insertions(+) create mode 100644 MC/config/common/external/generator/GeneratorHGPythia8.C create mode 100644 MC/config/common/ini/GeneratorHGPythia8PbPb536TeV.ini create mode 100644 MC/config/common/ini/tests/GeneratorHGPythia8PbPb536TeV.C diff --git a/MC/config/common/external/generator/GeneratorHGPythia8.C b/MC/config/common/external/generator/GeneratorHGPythia8.C new file mode 100644 index 000000000..0b9fbd915 --- /dev/null +++ b/MC/config/common/external/generator/GeneratorHGPythia8.C @@ -0,0 +1,460 @@ +// HG-PYTHIA: heavy-ion events built from a Glauber model with eikonal NN +// interactions and minijet (hard scattering) counting, where each NN +// collision is represented by a PYTHIA 8 pp event with the matching number +// of multiparton interactions. The heavy-ion event is the sum of the pp +// events. +// +// Model: C. Loizides, A. Morsch, Phys. Lett. B773 (2017) 408 (arXiv:1705.08856) +// Reference implementation: https://github.com/abaty/HGPythia (A. Baty) +// +// The PYTHIA pp configuration (beams at sqrt(s_NN), SoftQCD:inelastic, +// decays...) is given through the usual GeneratorPythia8 parameters, e.g. +// GeneratorPythia8.config=${O2DPG_MC_CONFIG_ROOT}/MC/config/common/pythia8/generator/pythia8_inel_536.cfg +// The collision energy used for the Glauber model is taken from PYTHIA. +// +// Differences with respect to the reference implementation: +// - the unused J/psi and bookkeeping parts are not ported +// +/// @author Marco Giacalone (marco.giacalone@cern.ch) +/// @date 09-2026 + +#include "Generators/GeneratorPythia8.h" +#include "SimConfig/SimConfig.h" +#include "SimulationDataFormat/MCEventHeader.h" +#include "SimulationDataFormat/MCGenProperties.h" +#include "SimulationDataFormat/ParticleStatus.h" +#include "Pythia8/Pythia.h" +#include "TF1.h" +#include "TF2.h" +#include "TMath.h" +#include "TParticle.h" +#include "TRandom3.h" +#include +#include +#include +#include +#include + +namespace o2 +{ +namespace eventgen +{ + +class GeneratorHGPythia8 : public GeneratorPythia8 +{ + public: + /// Number of bins for the hard scatterings per NN collision (and MPIs per PYTHIA event), as in the + /// reference implementation: collisions with 20 or more hard scatterings are counted in the last bin (19) + static constexpr int kMaxMPI = 20; + /// Maximum number of consecutive PYTHIA events not matching any requested number of MPIs + static constexpr long kMaxRejected = 10000000; + + GeneratorHGPythia8(int A = 208, int B = 208) : GeneratorPythia8(), mA(A), mB(B) + { + // the PYTHIA interface only describes the last NN collision: do not expose it as "pythia8" + mInterfaceName = "hgpythia8"; + // defaults, overridden by the PYTHIA configuration file (GeneratorPythia8.config) read in Init + readString("Beams:idA 2212"); + readString("Beams:idB 2212"); + readString("SoftQCD:inelastic on"); + readString("ParticleDecays:limitTau0 on"); + readString("ParticleDecays:tau0Max 10."); + } + ~GeneratorHGPythia8() override = default; + + /// Impact parameter range in fm. A non-positive bMax selects o2-sim --bMax, + /// if given, otherwise 20 fm (10 fm for p-A, 5 fm for pp) + void setImpactParameterRange(double bMin, double bMax) + { + mBMin = bMin; + mBMax = bMax; + } + /// Hard (minijet) NN cross-section in mb, <0 takes it from the built-in energy table + void setSigmaHard(double sigma) { mSigmaHardIn = sigma; } + /// Soft NN cross-section in mb + void setSigmaSoft(double sigma) { mSigmaSoft = sigma; } + /// Impact-parameter dependent shadowing of the hard cross-section + void setShadowing(bool val) { mShadowing = val; } + /// Include elastic NN scattering in the eikonal + void setElastic(bool val) { mElastic = val; } + + /// Glauber information of the current event + double getImpactParameter() const { return mImpactParameter; } + int getNcoll() const { return mNcoll; } + int getNcollHard() const { return mNcollHard; } + int getNhard() const { return mNhard; } + int getNpartProjectile() const { return mNpartProj; } + int getNpartTarget() const { return mNpartTarg; } + int getNpartBlackDisc() const { return mNpartBlackDisc; } + int getNcollBlackDisc() const { return mNcollBlackDisc; } + double getEccentricity() const { return mEccentricity; } + + Bool_t Init() override + { + if (mA < 1 || mA > 208 || mB < 1 || mB > 208) { + LOG(fatal) << "GeneratorHGPythia8: mass numbers must be between 1 and 208, got A = " << mA << ", B = " << mB; + return false; + } + // PYTHIA configuration, seeding and initialisation + if (!GeneratorPythia8::Init()) { + return false; + } + if (mPythia.settings.mode("Beams:idA") != 2212 || mPythia.settings.mode("Beams:idB") != 2212) { + LOG(warn) << "GeneratorHGPythia8: PYTHIA is expected to generate pp collisions"; + } + // Glauber random numbers seeded from the PYTHIA seed (0 means time dependent in both cases) + mRandom = std::make_unique(mPythia.settings.mode("Random:seed")); + + const double energy = mPythia.info.eCM(); + mSigmaHard = mSigmaHardIn >= 0 ? mSigmaHardIn : sigmaHardFromTable(energy); + if (mSigmaHard < 0) { + LOG(fatal) << "GeneratorHGPythia8: no hard cross-section available for sqrt(s_NN) = " << energy + << " GeV, please set it with setSigmaHard()"; + return false; + } + if (mBMax <= 0) { + const auto simBMax = o2::conf::SimConfig::Instance().getBMax(); + mBMax = simBMax > 0 ? simBMax : ((mA == 1 && mB == 1) ? 5. : ((mA == 1 || mB == 1) ? 10. : 20.)); + } + mDensity[0].reset(makeNucleonDensity(mA, "HGPythia8A")); + mDensity[1].reset(makeNucleonDensity(mB, "HGPythia8B")); + for (int j = 0; j < 2; ++j) { + mX[j].resize(j == 0 ? mA : mB); + mY[j].resize(j == 0 ? mA : mB); + mWounded[j].resize(j == 0 ? mA : mB); + mWoundedBlackDisc[j].resize(j == 0 ? mA : mB); + } + mHIEvent.init("HG-PYTHIA event", &mPythia.particleData); + + LOG(info) << "GeneratorHGPythia8: A = " << mA << ", B = " << mB << ", sqrt(s_NN) = " << energy + << " GeV, sigma_hard = " << mSigmaHard << " mb, sigma_soft = " << mSigmaSoft + << " mb, b in [" << mBMin << ", " << mBMax << "] fm, shadowing " << mShadowing + << ", elastic " << mElastic; + return true; + } + + Bool_t generateEvent() override + { + std::array nMPI{}; + sampleGlauber(nMPI); + + // Generate PYTHIA events until each NN collision has a partner with the + // matching number of MPIs, and sum them + mHIEvent.reset(); + int nNeeded = mNcoll; + long nRejected = 0; + const bool prune = !mGenConfig.includePartonEvent; + auto select = [this](const Pythia8::Particle& p) { + const int st = p.statusHepMC(); + return (st == 1 || st == 2 || st == 4) && mUserFilterFcn(p); + }; + while (nNeeded > 0) { + if (!mPythia.next()) { + continue; + } + const int mpi = mPythia.info.nMPI(); + if (mpi >= kMaxMPI || nMPI[mpi] == 0) { + if (++nRejected > kMaxRejected) { + LOG(fatal) << "GeneratorHGPythia8: " << kMaxRejected << " consecutive PYTHIA events do not match the " + << "requested number of MPIs, check the PYTHIA configuration (SoftQCD:inelastic needed)"; + return false; + } + continue; + } + nRejected = 0; + nMPI[mpi]--; + nNeeded--; + if (prune) { + pruneEvent(mPythia.event, select); + } + mHIEvent += mPythia.event; + } + return true; + } + + Bool_t importParticles() override + { + // same conversion as GeneratorPythia8::importParticles, without pruning again the summed event + for (int i = 1; i < mHIEvent.size(); ++i) { + const auto& particle = mHIEvent[i]; + auto st = o2::mcgenstatus::MCGenStatusEncoding(particle.statusHepMC(), particle.status()).fullEncoding; + mParticles.push_back(TParticle(particle.id(), st, + particle.mother1() - 1, particle.mother2() - 1, + particle.daughter1() - 1, particle.daughter2() - 1, + particle.px(), particle.py(), particle.pz(), particle.e(), + particle.xProd(), particle.yProd(), particle.zProd(), particle.tProd())); + mParticles.back().SetBit(ParticleStatus::kToBeDone, particle.statusHepMC() == 1); + } + return true; + } + + void updateHeader(o2::dataformats::MCEventHeader* eventHeader) override + { + using Key = o2::dataformats::MCInfoKeys; + eventHeader->putInfo(Key::generator, "hgpythia8"); + eventHeader->putInfo(Key::generatorVersion, PYTHIA_VERSION_INTEGER); + eventHeader->SetB(mImpactParameter); + eventHeader->putInfo(Key::impactParameter, mImpactParameter); + eventHeader->putInfo(Key::planeAngle, 0.); // impact parameter along x + eventHeader->putInfo(Key::nColl, mNcoll); + eventHeader->putInfo(Key::nCollHard, mNcollHard); + eventHeader->putInfo(Key::nPart, mNpartProj + mNpartTarg); + eventHeader->putInfo(Key::nPartProjectile, mNpartProj); + eventHeader->putInfo(Key::nPartTarget, mNpartTarg); + eventHeader->putInfo(Key::sigmaInelNN, mSigmaSoft + mSigmaHard); + eventHeader->putInfo("nHard", mNhard); + eventHeader->putInfo("Npart_blackdisc", mNpartBlackDisc); + eventHeader->putInfo("Ncoll_blackdisc", mNcollBlackDisc); + eventHeader->putInfo("eccentricity", mEccentricity); + } + + private: + /// Hard cross-section (pT > 2 GeV) in mb as a function of sqrt(s_NN), from the reference implementation + double sigmaHardFromTable(double energy) const + { + if (mShadowing) { + if (energy < 201) return 7.71968269; + if (energy < 2800) return 39.5451393; + if (energy < 5100) return 54.2068253; + if (energy < 8100) return 69.7713470; + return -1; + } + if (energy < 20) return 0.161440969; + if (energy < 40) return 1.07414019; + if (energy < 64) return 2.32993174; + if (energy < 201) return 11.6641903; + if (energy < 2800) return 85.2298813; + if (energy < 5100) return 124.296341; + if (energy < 5500) return 130.82; + if (energy < 6400) return 144.17; + if (energy < 8100) return 166.184998; + return -1; + } + + /// Radial (radial and polar for deformed nuclei) nucleon density; a single nucleon (A = 1) is + /// sampled from the exponential proton profile with R = 0.234 fm (rms charge radius), as in TGlauberMC + static TF1* makeNucleonDensity(int A, const std::string& suffix) + { + auto name = [&suffix](const char* n) { return std::string(n) + suffix; }; + if (A == 208) { + return new TF1(name("wsPb").c_str(), "7.208e-4*4.*TMath::Pi()*x^2/(1+exp((x-6.62)/0.546))", 0., 20.); + } + if (A == 197) { + return new TF1(name("wsAu").c_str(), "8.596e-04*4.*TMath::Pi()*x^2/(1+exp((x-6.38)/0.535))", 0., 20.); + } + if (A == 129) { + auto f = new TF2(name("wsXe2a").c_str(), "x*x*TMath::Sin(y)/(1+exp((x-[0]*(1+[2]*0.315*(3*pow(cos(y),2)-1.0)+[3]*0.105*(35*pow(cos(y),4)-30*pow(cos(y),2)+3)))/[1]))", 0, 15, 0.0, TMath::Pi()); + f->SetNpx(120); + f->SetNpy(120); + f->SetParameters(5.36, 0.59, 0.18, 0); + return f; + } + if (A == 40) { + return new TF1(name("wsAr").c_str(), "1.*TMath::Pi()*x^2/(1+exp((x-3.53)/0.542))", 0., 15.); + } + if (A == 20) { + auto f = new TF1(name("wsNe").c_str(), "x*x*(1+[2]*(x/[0])**2)/(1+exp((x-[0])/[1]))", 0, 10.); + f->SetParameters(2.791, 0.698, -0.168); + return f; + } + if (A == 16) { + auto f = new TF1(name("wsO").c_str(), "x*x*(1+[2]*(x/[0])**2)/(1+exp((x-[0])/[1]))", 0, 10.); + f->SetParameters(2.608, 0.513, -0.051); + return f; + } + if (A == 6) { + return new TF1(name("wsC").c_str(), "7.208e-4*4.*TMath::Pi()*x^2*(1.-0.149*(x/2.46)**2)/(1+exp((x-2.46)/0.522))", 0., 10.); + } + if (A == 3) { + return new TF1(name("wsHe").c_str(), "7.208e-4*4.*TMath::Pi()*x^2*(1.+0.517*(x/0.964)**2)/(1+exp((x-0.964)/0.322))", 0., 10.); + } + if (A == 1) { + return new TF1(name("prot").c_str(), "x*x*exp(-x/0.234)", 0., 5.); + } + LOG(fatal) << "GeneratorHGPythia8: nucleus with A = " << A << " not supported (1, 3, 6, 16, 20, 40, 129, 197, 208)"; + return nullptr; + } + + /// Matter distribution in the proton (eikonal) + static double eikonal(double x) + { + constexpr double p0 = 3.9, p1 = 96.; + return p0 * p0 / p1 * TMath::Power(p0 * x, 3) * TMath::BesselK(3, p0 * x); + } + + /// Sample the nucleon positions in the transverse plane for nucleus j, centred at x = dx + void sampleNucleus(int j, double dx) + { + auto density = mDensity[j].get(); + auto f2 = dynamic_cast(density); + for (size_t k = 0; k < mX[j].size(); ++k) { + double x = 0., y = 0.; + if (f2) { + double r, theta; + f2->GetRandom2(r, theta, mRandom.get()); + const double phi = 2. * TMath::Pi() * mRandom->Rndm(); + x = r * TMath::Sin(phi) * TMath::Sin(theta); + y = r * TMath::Cos(phi) * TMath::Sin(theta); + } else if (density) { + const double r = density->GetRandom(mRandom.get()); + const double phi = 2. * TMath::Pi() * mRandom->Rndm(); + const double costh = 2. * mRandom->Rndm() - 1.; + const double sinth = costh * costh < 1. ? TMath::Sqrt(1. - costh * costh) : 0.; + x = r * sinth * TMath::Cos(phi); + y = r * sinth * TMath::Sin(phi); + } + mX[j][k] = x + dx; + mY[j][k] = y; + } + } + + /// Sample a Glauber configuration with at least one inelastic NN collision and count, for + /// each NN collision, the number of hard scatterings (nMPI[0]: collisions without any) + void sampleGlauber(std::array& nMPI) + { + constexpr double dmax = 1.43; // black-disc NN distance (fm) + const double b02 = 0.5 * mSigmaSoft * 0.1 / TMath::Pi(); + do { + nMPI.fill(0); + const double b = TMath::Sqrt(mBMin * mBMin + mRandom->Rndm() * (mBMax * mBMax - mBMin * mBMin)); + sampleNucleus(0, b / 2.); + sampleNucleus(1, -b / 2.); + for (int j = 0; j < 2; ++j) { + std::fill(mWounded[j].begin(), mWounded[j].end(), 0); + std::fill(mWoundedBlackDisc[j].begin(), mWoundedBlackDisc[j].end(), 0); + } + // shadowing of the hard cross-section depends on b only + const double rrb = TMath::Min(1., b * b / 35.2 / 1.44); + const double aphx = 0.1 * 4. / 3. * 4.92 * TMath::Sqrt(1. - rrb); + const double sigHS = mShadowing ? mSigmaHard - aphx * 103.65 : mSigmaHard; + const double gstot0 = mElastic ? 2. * (1. - TMath::Exp(-(mSigmaSoft + sigHS) / mSigmaSoft * eikonal(0.001))) : 1.; + + mImpactParameter = b; + mNcoll = mNcollHard = mNhard = mNcollBlackDisc = 0; + for (int i = 0; i < mA; ++i) { + for (int j = 0; j < mB; ++j) { + const double dx = mX[0][i] - mX[1][j]; + const double dy = mY[0][i] - mY[1][j]; + double r2 = dx * dx + dy * dy; + if (r2 < dmax * dmax) { + mWoundedBlackDisc[0][i] = 1; + mWoundedBlackDisc[1][j] = 1; + mNcollBlackDisc++; + } + if (r2 > 25.) { + continue; + } + // interaction probability + r2 /= b02; + r2 /= gstot0; + const double chi = eikonal(TMath::Sqrt(r2)); + const double gs = 1. - TMath::Exp(-2. * (mSigmaSoft + sigHS) / mSigmaSoft * chi); + const double gstot = 2. * (1. - TMath::Sqrt(1. - gs)); + const double rantot = mRandom->Rndm() * gstot0; + if (rantot > gstot && mElastic) { + continue; + } + if (rantot > gs) { + continue; + } + mWounded[0][i] = 1; + mWounded[1][j] = 1; + mNcoll++; + // minijets + const double tt = 2. * chi * sigHS / mSigmaSoft; + const double ts = 2. * chi; + if (rantot < TMath::Exp(-tt) * (1. - TMath::Exp(-ts))) { + nMPI[0]++; + continue; + } + double xr = -TMath::Log(TMath::Exp(-tt) + mRandom->Rndm() * (1. - TMath::Exp(-tt))); + int njet = 0; + while (true) { + njet++; + xr -= TMath::Log(mRandom->Rndm()); + if (xr > tt) { + break; + } + } + njet = TMath::Min(njet, kMaxMPI - 1); + nMPI[njet]++; + mNhard += njet; + mNcollHard++; + } + } + } while (mNcoll < 1); + + // participants and participant eccentricity + mNpartProj = mNpartTarg = mNpartBlackDisc = 0; + double mx = 0., my = 0., mx2 = 0., my2 = 0., mxy = 0.; + for (int j = 0; j < 2; ++j) { + for (size_t k = 0; k < mX[j].size(); ++k) { + mNpartBlackDisc += mWoundedBlackDisc[j][k]; + if (!mWounded[j][k]) { + continue; + } + (j == 0 ? mNpartProj : mNpartTarg)++; + mx += mX[j][k]; + my += mY[j][k]; + mx2 += mX[j][k] * mX[j][k]; + my2 += mY[j][k] * mY[j][k]; + mxy += mX[j][k] * mY[j][k]; + } + } + const double iw = mNpartProj + mNpartTarg; + mx2 -= mx * mx / iw; + my2 -= my * my / iw; + mxy -= mx * my / iw; + mEccentricity = (mx2 + my2) > 0 ? TMath::Sqrt((my2 - mx2) * (my2 - mx2) + 4. * mxy * mxy) / (mx2 + my2) : 0.; + } + + // configuration + int mA = 208; + int mB = 208; + double mBMin = 0.; + double mBMax = -1.; + double mSigmaHardIn = -1.; + double mSigmaHard = -1.; + double mSigmaSoft = 57.; + bool mShadowing = false; + bool mElastic = false; + + // Glauber state + std::unique_ptr mRandom; + std::array, 2> mDensity; + std::array, 2> mX; + std::array, 2> mY; + std::array, 2> mWounded; + std::array, 2> mWoundedBlackDisc; + double mImpactParameter = 0.; + int mNcoll = 0; + int mNcollHard = 0; + int mNhard = 0; + int mNpartProj = 0; + int mNpartTarg = 0; + int mNpartBlackDisc = 0; + int mNcollBlackDisc = 0; + double mEccentricity = 0.; + + // sum of the PYTHIA events of the current heavy-ion event + Pythia8::Event mHIEvent; +}; + +} // namespace eventgen +} // namespace o2 + +/// HG-PYTHIA generator for nuclei with mass numbers A (along +z, PYTHIA beam A) and B. +/// bMax <= 0: o2-sim --bMax if given, otherwise 20 fm (10 fm for p-A, 5 fm for pp). +/// sigmaHard < 0: hard cross-section from the built-in sqrt(s_NN) table (up to 8.1 TeV). +FairGenerator* generateHGPythia8(int A = 208, int B = 208, double bMin = 0., double bMax = -1., + double sigmaHard = -1., double sigmaSoft = 57., + bool shadowing = false, bool elastic = false) +{ + auto gen = new o2::eventgen::GeneratorHGPythia8(A, B); + gen->setImpactParameterRange(bMin, bMax); + gen->setSigmaHard(sigmaHard); + gen->setSigmaSoft(sigmaSoft); + gen->setShadowing(shadowing); + gen->setElastic(elastic); + return gen; +} diff --git a/MC/config/common/ini/GeneratorHGPythia8PbPb536TeV.ini b/MC/config/common/ini/GeneratorHGPythia8PbPb536TeV.ini new file mode 100644 index 000000000..ea1e57c76 --- /dev/null +++ b/MC/config/common/ini/GeneratorHGPythia8PbPb536TeV.ini @@ -0,0 +1,8 @@ +#NEV_TEST> 5 +# HG-PYTHIA (Glauber + PYTHIA 8 pp collisions) Pb-Pb at sqrt(s_NN) = 5.36 TeV, minimum bias +[GeneratorExternal] +fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/common/external/generator/GeneratorHGPythia8.C +funcName=generateHGPythia8(208, 208) + +[GeneratorPythia8] +config=${O2DPG_MC_CONFIG_ROOT}/MC/config/common/pythia8/generator/pythia8_inel_536.cfg diff --git a/MC/config/common/ini/tests/GeneratorHGPythia8PbPb536TeV.C b/MC/config/common/ini/tests/GeneratorHGPythia8PbPb536TeV.C new file mode 100644 index 000000000..a54dcabba --- /dev/null +++ b/MC/config/common/ini/tests/GeneratorHGPythia8PbPb536TeV.C @@ -0,0 +1,56 @@ +int External() +{ + const int A = 208; // mass number of both nuclei + std::string path{"o2sim_Kine.root"}; + 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* header = nullptr; + tree->SetBranchAddress("MCEventHeader.", &header); + + using Key = o2::dataformats::MCInfoKeys; + const auto nEvents = tree->GetEntries(); + if (nEvents == 0) { + std::cerr << "No events found\n"; + return 1; + } + for (Long64_t i = 0; i < nEvents; ++i) { + tree->GetEntry(i); + bool valid = false; + const int nColl = header->getInfo(Key::nColl, valid); + if (!valid || nColl < 1 || nColl > A * A) { + std::cerr << "Missing or invalid Ncoll in event " << i << "\n"; + return 1; + } + const int nPart = header->getInfo(Key::nPart, valid); + if (!valid || nPart < 2 || nPart > 2 * A) { + std::cerr << "Missing or invalid Npart in event " << i << "\n"; + return 1; + } + // each NN collision is a PYTHIA pp event: two beam protons with half of sqrt(s_NN) each + int nBeams = 0; + for (const auto& track : *tracks) { + if (o2::mcgenstatus::getHepMCStatusCode(track.getStatusCode()) == 4) { + if (track.GetPdgCode() != 2212 || std::abs(track.GetEnergy() - 2680.) > 1e-3) { + std::cerr << "Unexpected beam particle " << track.GetPdgCode() << " with energy " << track.GetEnergy() << " in event " << i << "\n"; + return 1; + } + nBeams++; + } + } + if (nBeams != 2 * nColl) { + std::cerr << "Found " << nBeams << " beam protons for " << nColl << " NN collisions in event " << i << "\n"; + return 1; + } + } + return 0; +}