Skip to content

[PWGLF] Appended the injected event in MB event correctly #2481

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
njacazio merged 3 commits into from
Sep 28, 2026
Merged

[PWGLF] Appended the injected event in MB event correctly #2481

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
8912261
Appended the injected event in MB correctly
sawankumawat Sep 28, 2026
3e9d3ee
Made phi max pT configurable
sawankumawat Sep 28, 2026
5ca7311
Custom phi injected resonance
sawankumawat Sep 28, 2026
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
2 changes: 1 addition & 1 deletion MC/config/PWGLF/ini/GeneratorLF_phiphiInvMass.ini
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,6 +1,6 @@
[GeneratorExternal]
fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/pythia8/GeneratorLF_phiphi.C
funcName=generatePhiResonanceGun(999999, 0.0, 50.0, -1.0, 1.0, "${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/pythia8/generator/pythia8_inel_136tev.cfg", 3)
funcName=generatePhiResonanceGun(999999, 888888, 0.0, 50.0, 100.0, -1.0, 1.0, "${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/pythia8/generator/pythia8_inel_136tev.cfg", 3)

[GeneratorPythia8]
config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/pythia8/generator/pythia8_inel_136tev.cfg
Expand Down
60 changes: 31 additions & 29 deletions MC/config/PWGLF/ini/tests/GeneratorLF_phiphiInvMass.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -76,14 +76,34 @@ int External()
for (int k = phiTrack.getFirstDaughterTrackId(); k <= phiTrack.getLastDaughterTrackId(); ++k)
{
auto grandDauPdg = tracks->at(k).GetPdgCode();
if (grandDauPdg == 321) ++nKPlusFromResonancePhi;
if (grandDauPdg == -321) ++nKMinusFromResonancePhi;
if (grandDauPdg == 321)
++nKPlusFromResonancePhi;
if (grandDauPdg == -321)
++nKMinusFromResonancePhi;
}
}
}
}
}

else if (pdg == 888888)
{
++nDirectInjectedPhi;
hasInjection = true;

if (track.getFirstDaughterTrackId() >= 0)
{
for (int j = track.getFirstDaughterTrackId(); j <= track.getLastDaughterTrackId(); ++j)
{
auto dauPdg = tracks->at(j).GetPdgCode();
if (dauPdg == 321)
++nKPlusFromDirectPhi;
if (dauPdg == -321)
++nKMinusFromDirectPhi;
}
}
}

// 2. Process Phi (333) Mesons
else if (pdg == 333)
{
Expand All @@ -96,36 +116,18 @@ int External()
continue;
}

bool isDirectInjected = (motherId < 0);
// Minimum Bias Phi
++nMBPhi;

if (isDirectInjected)
if (track.getFirstDaughterTrackId() >= 0)
{
++nDirectInjectedPhi;
hasInjection = true;

if (track.getFirstDaughterTrackId() >= 0)
for (int j = track.getFirstDaughterTrackId(); j <= track.getLastDaughterTrackId(); ++j)
{
for (int j = track.getFirstDaughterTrackId(); j <= track.getLastDaughterTrackId(); ++j)
{
auto dauPdg = tracks->at(j).GetPdgCode();
if (dauPdg == 321) ++nKPlusFromDirectPhi;
if (dauPdg == -321) ++nKMinusFromDirectPhi;
}
}
}
else
{
// Minimum Bias Phi
++nMBPhi;

if (track.getFirstDaughterTrackId() >= 0)
{
for (int j = track.getFirstDaughterTrackId(); j <= track.getLastDaughterTrackId(); ++j)
{
auto dauPdg = tracks->at(j).GetPdgCode();
if (dauPdg == 321) ++nKPlusFromMBPhi;
if (dauPdg == -321) ++nKMinusFromMBPhi;
}
auto dauPdg = tracks->at(j).GetPdgCode();
if (dauPdg == 321)
++nKPlusFromMBPhi;
if (dauPdg == -321)
++nKMinusFromMBPhi;
}
}
}
Expand Down
155 changes: 73 additions & 82 deletions MC/config/PWGLF/pythia8/GeneratorLF_phiphi.C
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -6,7 +6,6 @@
#include "Pythia8/Pythia.h"
#include "TRandom3.h"
#include "TMath.h"
// #include "TF1.h"
#include "TParticle.h"
#include "TSystem.h"
#if __has_include("SimulationDataFormat/MCGenStatus.h")
Expand All @@ -19,40 +18,51 @@
#endif
#include <cmath>
#include <string>
#include <vector>
#endif

// Double_t FuncLavy(Double_t *x, Double_t *par)
// {

// Double_t p = (par[0] - 1) * (par[0] - 2) * par[1] * x[0] / (((pow((1 + (((sqrt((par[2] * par[2]) + (x[0] * x[0]))) - par[2]) / (par[0] * par[3]))), par[0]) * (par[0] * par[3] * ((par[0] * par[3]) + (par[2] * (par[0] - 2)))))));
// return (p);
// }

class GeneratorPhiResonance : public o2::eventgen::GeneratorPythia8
{
public:
GeneratorPhiResonance(int resoPDG = 999999,
float ptMin = 0.0, float ptMax = 50.0,
int customPhiPDG = 888888,
float ptMin = 0.0, float ptMax = 50.0, float ptMaxPhi = 100.0,
float yMin = -1.0, float yMax = 1.0,
std::string pythiaCfgMb = "${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/pythia8/generator/pythia8_inel_136tev.cfg",
int signalInterval = 3)
: GeneratorPythia8(), mResoPDG(resoPDG), mPtMin(ptMin), mPtMaxPhiPhi(ptMax), mYMin(yMin), mYMax(yMax), mSignalInterval(signalInterval)
: GeneratorPythia8(), mResoPDG(resoPDG), mCustomPhiPDG(customPhiPDG), mPtMin(ptMin), mPtMaxPhiPhi(ptMax), mPtMaxPhi(ptMaxPhi), mYMin(yMin), mYMax(yMax), mSignalInterval(signalInterval)
{
// 1. Initialize Gun Pythia object & define custom resonance
// # id::all = name antiName spinType chargeType colType m0 mWidth mMin mMax tau0
// 1. Define Custom Directly Injected Phi (PDG: 888888) with mass, width, and decay to kaons
std::string createCustomPhi = std::to_string(mCustomPhiPDG) + ":new = custom_phi custom_phi 3 0 0 1.019461 0.004249 0.980 1.100 0.0";
std::string customPhiMayDecay = std::to_string(mCustomPhiPDG) + ":mayDecay = on";
std::string addPhiDecayKPlusKMinus = std::to_string(mCustomPhiPDG) + ":addChannel = 1 0.492 0 321 -321";

// 2. Define Custom Signal Resonance (PDG: 999999) decay into standard Phis (333 333)
std::string createReso = std::to_string(mResoPDG) + ":new = f2_Custom void 5 0 0 2.714 0.012 2.05 3.50 0.0";
mPythiaGun.readString(createReso);
mPythiaGun.readString(std::to_string(mResoPDG) + ":mayDecay = on");
// id:addChannel = onMode bRatio meMode product1 product2, (onMode = 1: allow decay, bRatio = branching ratio, meMode = matrix element mode where 0 is isotropic decay)
std::string addDecay = std::to_string(mResoPDG) + ":addChannel = 1 1.0 0 333 333";
mPythiaGun.readString(addDecay);
std::string resoMayDecay = std::to_string(mResoPDG) + ":mayDecay = on";
std::string addResoDecay = std::to_string(mResoPDG) + ":addChannel = 1 1.0 0 333 333";

// Helper lambda to load custom particle definitions across ALL Pythia engines
auto applyCustomParticles = [&](Pythia8::Pythia &pythiaInst)
{
pythiaInst.readString(createCustomPhi);
pythiaInst.readString(customPhiMayDecay);
pythiaInst.readString(addPhiDecayKPlusKMinus);
pythiaInst.readString(createReso);
pythiaInst.readString(resoMayDecay);
pythiaInst.readString(addResoDecay);
};

// 1: Apply particle definitions to mPythia, mPythiaGun, and pythiaObjectMinimumBias
applyCustomParticles(mPythia);
applyCustomParticles(mPythiaGun);

mPythiaGun.readString("ProcessLevel:all off");
mPythiaGun.readString("Random:setSeed = on");
mPythiaGun.readString("Random:seed = " + std::to_string(1 + gRandom->Integer(900000000)));
mPythiaGun.init();

// 2. Initialize Minimum Bias Pythia Engine
// 3. Initialize Minimum Bias Pythia Engine
if (pythiaCfgMb.empty())
{
auto &param = o2::eventgen::GeneratorPythia8Param::Instance();
Expand All @@ -67,25 +77,8 @@ public:
pythiaObjectMinimumBias.readString("Random:setSeed = on");
pythiaObjectMinimumBias.readString("Random:seed = " + std::to_string(1 + gRandom->Integer(900000000)));

// Add custom particle definition to MB instance so particle table matches
pythiaObjectMinimumBias.readString(createReso);
pythiaObjectMinimumBias.readString(std::to_string(mResoPDG) + ":mayDecay = on");
pythiaObjectMinimumBias.readString(addDecay);

applyCustomParticles(pythiaObjectMinimumBias);
pythiaObjectMinimumBias.init();

// // Thermal pT distribution for phi-phi resonance
// mThermal = new TF1("mThermal", "x*sqrt(x*x+[0]*[0])*exp(-sqrt(x*x+[0]*[0])/[1])", mPtMin, mPtMaxPhiPhi);

// // Lévy-Tsallis pT distribution for direct phi
// mLevyTsallis = new TF1("mLevyTsallis", FuncLavy, mPtMin, 100.0, 4);

// mLevyTsallis->SetParameters(
// 7.60279, // n
// 0.0374237, // dN/dy
// 1.01946, // mass
// 0.338379 // T
// );
}

Bool_t generateEvent() override
Expand All @@ -99,69 +92,77 @@ public:
mbOK = pythiaObjectMinimumBias.next();
}

// Copy MB event into mPythia
mPythia.event = pythiaObjectMinimumBias.event;
// 2: Copy MB event using exact pointer binding from generator_pythia8_LF_rapidity_width.C
copyMinimumBiasEventForInjection();

// 2. Clear Gun event container
mPythiaGun.event.reset();

// 3. Inject Signal Gun Particles into mPythiaGun
if (mEventCounter % mSignalInterval == 0)
{
// Theramal distribution
// Resonant signal -> Decays into 333 333 (Standard Phis)
injectParticle(mResoPDG, 1, true);
}
else
{
// From published
injectParticle(333, 2, false);
// Directly injected uncorrelated Phi -> Uses Custom PDG 888888
injectParticle(mCustomPhiPDG, 2, false);
}

// 4. Force Decay of injected particles using Pythia's Decayer
for (int i = 1; i < mPythiaGun.event.size(); ++i)
{
if (mPythiaGun.event[i].status() > 0)
{ // Active injected particles
mPythiaGun.particleData.mayDecay(mPythiaGun.event[i].id(), true);
mPythiaGun.moreDecays();
}
}
mPythiaGun.moreDecays();
mPythiaGun.next();

// 5. Merge mPythiaGun event into mPythia.event
// 3: Index Mapping during Event Merging
int offset = mPythia.event.size();
std::vector<int> indexMap(mPythiaGun.event.size(), 0);

for (int i = 1; i < mPythiaGun.event.size(); ++i)
{ // Skip system particle 0
{
indexMap[i] = mPythia.event.size();
Pythia8::Particle p = mPythiaGun.event[i];
mPythia.event.append(p);
}

// Adjust history indices accurately
int mother1 = (p.mother1() > 0) ? p.mother1() + offset - 1 : p.mother1();
int mother2 = (p.mother2() > 0) ? p.mother2() + offset - 1 : p.mother2();
int daughter1 = (p.daughter1() > 0) ? p.daughter1() + offset - 1 : p.daughter1();
int daughter2 = (p.daughter2() > 0) ? p.daughter2() + offset - 1 : p.daughter2();
// Re-link mother and daughter index relationships accurately
for (int i = 1; i < mPythiaGun.event.size(); ++i)
{
int newIdx = indexMap[i];
Pythia8::Particle &p = mPythia.event[newIdx];

p.mothers(mother1, mother2);
p.daughters(daughter1, daughter2);
int m1 = p.mother1();
int m2 = p.mother2();
int d1 = p.daughter1();
int d2 = p.daughter2();

mPythia.event.append(p);
p.mothers((m1 > 0 && m1 < (int)indexMap.size()) ? indexMap[m1] : 0,
(m2 > 0 && m2 < (int)indexMap.size()) ? indexMap[m2] : 0);

p.daughters((d1 > 0 && d1 < (int)indexMap.size()) ? indexMap[d1] : 0,
(d2 > 0 && d2 < (int)indexMap.size()) ? indexMap[d2] : 0);
}

// 6. CRITICAL: Restore Pythia particleData pointers for O2 exporter
// 4: Restore Pythia particleData pointers for O2 exporter
mPythia.event.restorePtrs();

// 7. Invoke base generator hooks to sync O2 event record
// return GeneratorPythia8::generateEvent();
return true; // Skip base generator processing to avoid overwriting injected particles
return true;
}

private:
void copyMinimumBiasEventForInjection()
{
mPythia.event = pythiaObjectMinimumBias.event;
mPythia.event.init("Minimum-bias event with injected particles", &mPythia.particleData);
mPythia.event.restorePtrs();
}

void injectParticle(int pdg, int nParticles, bool thermalPt)
{
const double phiMass = 1.019461;

for (int i = 0; i < nParticles; ++i)
{
// const double pt = gRandom->Uniform(mPtMin, mPtMaxPhiPhi);
const double y = gRandom->Uniform(mYMin, mYMax);
const double phi = gRandom->Uniform(0, TMath::TwoPi());

Expand All @@ -175,26 +176,17 @@ private:
}
else
{
mass = mPythiaGun.particleData.mSel(pdg);
mass = mPythiaGun.particleData.mSel(333); // Use standard phi mass for directly injected custom phi
}

double pt;

if (thermalPt)
{
// const double T = 0.160;

// mThermal->SetParameter(0, mass);
// mThermal->SetParameter(1, T);

// pt = mThermal->GetRandom();

pt = gRandom->Uniform(mPtMin, mPtMaxPhiPhi); // Falling back to flat pT due to low statistics in high pT
pt = gRandom->Uniform(mPtMin, mPtMaxPhiPhi);
}
else
{
// pt = mLevyTsallis->GetRandom();
pt = gRandom->Uniform(mPtMin, 100.0); // Falling back to flat pT due to low statistics in high pT
pt = gRandom->Uniform(mPtMin, mPtMaxPhi);
}

const double px = pt * std::cos(phi);
Expand All @@ -215,24 +207,23 @@ private:
particle.yProd(0.);
particle.zProd(0.);

mPythiaGun.particleData.mayDecay(pdg, true);
mPythiaGun.event.append(particle);
}
}

int mEventCounter = 0;
int mResoPDG;
int mCustomPhiPDG;
int mSignalInterval;
float mPtMin, mPtMaxPhiPhi, mYMin, mYMax;
float mPtMin, mPtMaxPhiPhi, mPtMaxPhi, mYMin, mYMax;

Pythia8::Pythia mPythiaGun;
Pythia8::Pythia pythiaObjectMinimumBias;

// TF1 *mThermal;
// TF1 *mLevyTsallis;
};

/// Entry point for o2-sim
FairGenerator *generatePhiResonanceGun(int resoPDG = 999999, float ptMin = 0.0, float ptMax = 50.0, float yMin = -1.0, float yMax = 1.0, std::string pythiaCfgMb = "", int signalInterval = 3)
FairGenerator *generatePhiResonanceGun(int resoPDG = 999999, int customPhiPDG = 888888, float ptMin = 0.0, float ptMax = 50.0, float ptMaxPhi = 100.0, float yMin = -1.0, float yMax = 1.0, std::string pythiaCfgMb = "", int signalInterval = 3)
{
return new GeneratorPhiResonance(resoPDG, ptMin, ptMax, yMin, yMax, pythiaCfgMb, signalInterval);
return new GeneratorPhiResonance(resoPDG, customPhiPDG, ptMin, ptMax, ptMaxPhi, yMin, yMax, pythiaCfgMb, signalInterval);
}
Loading