Appended the injected event in MB correctly

Author: sawankumawat

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,15 +18,9 @@
 #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:
@@ -38,21 +31,30 @@ public:
                           int signalInterval = 3)
         : GeneratorPythia8(), mResoPDG(resoPDG), mPtMin(ptMin), mPtMaxPhiPhi(ptMax), 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 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(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 +69,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 +84,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
+            // Directly injected uncorrelated Phi
             injectParticle(333, 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);
+        }
+
+        // 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];
 
-            // 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();
+            int m1 = p.mother1();
+            int m2 = p.mother2();
+            int d1 = p.daughter1();
+            int d2 = p.daughter2();
 
-            p.mothers(mother1, mother2);
-            p.daughters(daughter1, daughter2);
+            p.mothers((m1 > 0 && m1 < (int)indexMap.size()) ? indexMap[m1] : 0,
+                      (m2 > 0 && m2 < (int)indexMap.size()) ? indexMap[m2] : 0);
 
-            mPythia.event.append(p);
+            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());
 
@@ -179,22 +172,13 @@ private:
             }
 
             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, 100.0);
             }
 
             const double px = pt * std::cos(phi);
@@ -215,6 +199,7 @@ private:
             particle.yProd(0.);
             particle.zProd(0.);
 
+            mPythiaGun.particleData.mayDecay(pdg, true);
             mPythiaGun.event.append(particle);
         }
     }
@@ -226,9 +211,6 @@ private:
 
     Pythia8::Pythia mPythiaGun;
     Pythia8::Pythia pythiaObjectMinimumBias;
-
-    // TF1 *mThermal;
-    // TF1 *mLevyTsallis;
 };
 
 /// Entry point for o2-sim