diff --git a/MC/config/PWGLF/ini/GeneratorLF_phiphiInvMass.ini b/MC/config/PWGLF/ini/GeneratorLF_phiphiInvMass.ini index 097c63c3e..d9ad4041c 100644 --- a/MC/config/PWGLF/ini/GeneratorLF_phiphiInvMass.ini +++ b/MC/config/PWGLF/ini/GeneratorLF_phiphiInvMass.ini @@ -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 diff --git a/MC/config/PWGLF/ini/tests/GeneratorLF_phiphiInvMass.C b/MC/config/PWGLF/ini/tests/GeneratorLF_phiphiInvMass.C index 1ae3c2e12..4f08d32f8 100644 --- a/MC/config/PWGLF/ini/tests/GeneratorLF_phiphiInvMass.C +++ b/MC/config/PWGLF/ini/tests/GeneratorLF_phiphiInvMass.C @@ -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) { @@ -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; } } } diff --git a/MC/config/PWGLF/pythia8/GeneratorLF_phiphi.C b/MC/config/PWGLF/pythia8/GeneratorLF_phiphi.C index e8f7dad5e..8537ebd3a 100644 --- a/MC/config/PWGLF/pythia8/GeneratorLF_phiphi.C +++ b/MC/config/PWGLF/pythia8/GeneratorLF_phiphi.C @@ -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") @@ -19,40 +18,51 @@ #endif #include #include +#include #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 ¶m = o2::eventgen::GeneratorPythia8Param::Instance(); @@ -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 @@ -99,8 +92,8 @@ 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(); @@ -108,60 +101,68 @@ public: // 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 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()); @@ -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); @@ -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); } \ No newline at end of file