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

* Appended the injected event in MB correctly

* Made phi max pT configurable

* Custom phi injected resonance

Author: sawankumawat

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

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;
                     }
                 }
             }

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 <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();
@@ -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,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());
 
@@ -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);
 }