[PWGLF] Added phi-phi trigger code with additional selections

MC/config/PWGLF/ini/GeneratorLF_doublephi_triggerMasspTcut.ini

@@ -0,0 +1,10 @@
+[GeneratorExternal]
+fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/pythia8/generator_pythia8_twophi_triggerMassCut.C
+funcName=generateDoublePhi(0, 0.0, 100.0, 0.8)
+
+[GeneratorPythia8]
+config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/pythia8/generator/pythia8_inel_136tev.cfg
+
+[DecayerPythia8]
+config[0]=${O2DPG_MC_CONFIG_ROOT}/MC/config/common/pythia8/decayer/base.cfg
+config[1]=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/pythia8/generator/resonances.cfg

MC/config/PWGLF/ini/tests/GeneratorLF_doublephi_triggerMasspTcut.C

@@ -0,0 +1,76 @@
+int External()
+{
+    const 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<o2::MCTrack> *tracks{};
+    tree->SetBranchAddress("MCTrack", &tracks);
+
+    // Counters
+    int nMBPhi = 0;
+    int nKPlusFromMBPhi = 0;
+    int nKMinusFromMBPhi = 0;
+    int numberOfEventsProcessed = 0;
+
+    for (Long64_t i = 0; i < tree->GetEntries(); ++i)
+    {
+        tree->GetEntry(i);
+        ++numberOfEventsProcessed;
+
+        for (size_t idx = 0; idx < tracks->size(); ++idx)
+        {
+            const auto &track = tracks->at(idx);
+            const auto pdg = track.GetPdgCode();
+
+            if (pdg == 333)
+            {
+                ++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;
+                        }
+                    }
+                }
+            }
+        }
+    }
+
+    // --------------------------- Output ---------------------------
+    std::cout << "=================================================\n";
+    std::cout << "Total Events: " << tree->GetEntries() << "\n\n";
+    std::cout << "Total events processed: " << numberOfEventsProcessed << "\n";
+    std::cout << "Total Minimum Bias Phi (333): " << nMBPhi << "\n";
+    std::cout << "  -> Decayed to K+: " << nKPlusFromMBPhi << "\n";
+    std::cout << "  -> Decayed to K-: " << nKMinusFromMBPhi << "\n";
+    std::cout << "=================================================\n";
+
+    return 0;
+}
+
+void GeneratorLF_doublephi_trigger()
+{
+    External();
+}

MC/config/PWGLF/pythia8/generator_pythia8_twophi_trigger.C

@@ -127,15 +127,15 @@ protected:
             if (isPhiFromHFDecay(p, event))
                 continue;
 
-            // // Avoid double-counting copy/history entries:
-            // // Ensure this is the physical produced phi (e.g. check if its daughter is a copy of itself)
-            // int d1 = p.daughter1();
-            // int d2 = p.daughter2();
-            // if (d1 > 0 && d1 == d2 && std::abs(event[d1].id()) == 333)
-            // {
-            //     // p decayed into another copy of phi, so skip this intermediate entry
-            //     continue;
-            // }
+            // Avoid double-counting copy/history entries:
+            // Ensure this is the physical produced phi (e.g. check if its daughter is a copy of itself)
+            int d1 = p.daughter1();
+            int d2 = p.daughter2();
+            if (d1 > 0 && d1 == d2 && std::abs(event[d1].id()) == 333)
+            {
+                // p decayed into another copy of phi, so skip this intermediate entry
+                continue;
+            }
 
             nPhi++;
         }

MC/config/PWGLF/pythia8/generator_pythia8_twophi_triggerMassCut.C

@@ -0,0 +1,188 @@
+#if !defined(__CLING__) || defined(__ROOTCLING__)
+#include "FairGenerator.h"
+#include "FairPrimaryGenerator.h"
+#include "Generators/GeneratorPythia8.h"
+#include "Pythia8/Pythia.h"
+#include "TDatabasePDG.h"
+#include "TMath.h"
+#include "TParticlePDG.h"
+#include "TRandom3.h"
+#include "TSystem.h"
+#include "TVector2.h"
+#include "fairlogger/Logger.h"
+#include <cmath>
+#include <fstream>
+#include <string>
+#include "TLorentzVector.h"
+#include <vector>
+using namespace Pythia8;
+#endif
+
+/// Event generator using Pythia ropes
+/// Triggers events containing at least two generated phi(1020) mesons.
+
+class GeneratorPythia8DoublePhi : public o2::eventgen::GeneratorPythia8
+{
+public:
+    /// Constructor
+    GeneratorPythia8DoublePhi(int gapSize = 0, double minPt = 0.0, double maxPt = 100.0, double maxEta = 0.8)
+        : o2::eventgen::GeneratorPythia8(),
+          mGapSize(gapSize),
+          mMinPt(minPt),
+          mMaxPt(maxPt),
+          mMaxEta(maxEta)
+    {
+        fmt::printf(">> Pythia8 generator: two phi(1020) mesons, gap = %d, minPtPhi = %f, maxPtPhi = %f, |etaPhi| < %f\n", gapSize, minPt, maxPt, maxEta);
+    }
+    /// Destructor
+    ~GeneratorPythia8DoublePhi() = default;
+
+    bool Init() override
+    {
+        addSubGenerator(0, "Pythia8 events with two phi(1020) mesons");
+        return o2::eventgen::GeneratorPythia8::Init();
+    }
+
+protected:
+    bool isPhiFromHFDecay(const Pythia8::Particle &p, const Pythia8::Event &event)
+    {
+
+        // Walk up ancestry
+        int motherId = p.mother1();
+
+        while (motherId > 0)
+        {
+            // Get mother
+            const auto &mother = event[motherId];
+            const int absMotherPdg = std::abs(mother.id());
+
+            // Check if particle is from HF decay
+            if (((absMotherPdg / 100) % 10 == 4) ||
+                ((absMotherPdg / 100) % 10 == 5) ||
+                ((absMotherPdg / 1000) % 10 == 4) ||
+                ((absMotherPdg / 1000) % 10 == 5))
+            {
+                return true;
+            }
+
+            motherId = mother.mother1();
+        }
+        return false;
+    }
+
+    bool generateEvent() override
+    {
+        // fmt::printf(">> Generating event %d\n", mGeneratedEvents);
+
+        bool genOk = false;
+        int localCounter{0};
+        constexpr int kMaxTries{300000};
+
+        // If mGapSize <= 0, filter ALL events to contain two phis.
+        // Otherwise, generate mGapSize gap events before 1 triggered event.
+        if (mGapSize > 0 && (mGeneratedEvents % (mGapSize + 1) < mGapSize))
+        {
+            genOk = GeneratorPythia8::generateEvent();
+            // fmt::printf(">> Gap-event (no phi check)\n");
+        }
+        else
+        {
+            while (!genOk && localCounter < kMaxTries)
+            {
+                if (GeneratorPythia8::generateEvent())
+                {
+                    genOk = selectEvent(mPythia.event);
+                }
+                localCounter++;
+            }
+            if (!genOk)
+            {
+                fmt::printf("Failed to generate triggered event after %d tries\n", kMaxTries);
+                return false;
+            }
+            fmt::printf(">> Triggered event: event accepted after %d iterations (double phi(1020))\n", localCounter);
+        }
+
+        notifySubGenerator(0);
+        mGeneratedEvents++;
+        return true;
+    }
+
+    bool selectEvent(Pythia8::Event &event)
+    {
+        std::vector<TLorentzVector> phiCandidates;
+
+        for (int i = 0; i < event.size(); i++)
+        {
+            const auto &p = event[i];
+
+            if (std::abs(p.id()) != 333)
+                continue;
+
+            if (p.pT() < mMinPt || p.pT() > mMaxPt)
+                continue;
+
+            if (std::abs(p.eta()) > mMaxEta)
+                continue;
+
+            if (isPhiFromHFDecay(p, event))
+                continue;
+
+            // Avoid double-counting copy/history entries:
+            // Ensure this is the physical produced phi (e.g. check if its daughter is a copy of itself)
+            int d1 = p.daughter1();
+            int d2 = p.daughter2();
+            if (d1 > 0 && d1 == d2 && std::abs(event[d1].id()) == 333)
+            {
+                // p decayed into another copy of phi, so skip this intermediate entry
+                continue;
+            }
+
+            TLorentzVector phi;
+            phi.SetPtEtaPhiM(p.pT(), p.eta(), p.phi(), p.m());
+
+            phiCandidates.push_back(phi);
+        }
+        if (phiCandidates.size() < 2)
+            return false;
+
+        // Check all possible phi-phi pairs
+        for (size_t i = 0; i < phiCandidates.size(); i++)
+        {
+            for (size_t j = i + 1; j < phiCandidates.size(); j++)
+            {
+                TLorentzVector phiPhi = phiCandidates[i] + phiCandidates[j];
+
+                double mass = phiPhi.M();
+                double pt = phiPhi.Pt();
+
+                if (mass > 2.4 && pt > 4.0)
+                    return true;
+            }
+        }
+
+        return false;
+    }
+
+private:
+    int mGapSize{0};
+    double mMinPt{0.0};
+    double mMaxPt{100.0};
+    double mMaxEta{0.8};
+    uint64_t mGeneratedEvents{0};
+};
+
+///___________________________________________________________
+FairGenerator *generateDoublePhi(int gap = 0, double minPt = 0.0, double maxPt = 100.0, double maxEta = 0.8)
+{
+    auto myGenerator = new GeneratorPythia8DoublePhi(gap, minPt, maxPt, maxEta);
+
+    myGenerator->readString("333:onMode = off");
+    myGenerator->readString("333:onIfMatch = 321 -321");
+
+    auto seed = (gRandom->TRandom::GetSeed() % 900000000);
+    myGenerator->readString("Random:setSeed on");
+    myGenerator->readString("Random:seed " + std::to_string(seed));
+
+    return myGenerator;
+}