diff --git a/MC/config/PWGLF/ini/GeneratorLF_doublephi_triggerMasspTcut.ini b/MC/config/PWGLF/ini/GeneratorLF_doublephi_triggerMasspTcut.ini new file mode 100644 index 000000000..a51496b44 --- /dev/null +++ b/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 diff --git a/MC/config/PWGLF/ini/tests/GeneratorLF_doublephi_triggerMasspTcut.C b/MC/config/PWGLF/ini/tests/GeneratorLF_doublephi_triggerMasspTcut.C new file mode 100644 index 000000000..58767f843 --- /dev/null +++ b/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 *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(); +} diff --git a/MC/config/PWGLF/pythia8/generator_pythia8_twophi_trigger.C b/MC/config/PWGLF/pythia8/generator_pythia8_twophi_trigger.C index 8de46d445..af6366971 100644 --- a/MC/config/PWGLF/pythia8/generator_pythia8_twophi_trigger.C +++ b/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++; } diff --git a/MC/config/PWGLF/pythia8/generator_pythia8_twophi_triggerMassCut.C b/MC/config/PWGLF/pythia8/generator_pythia8_twophi_triggerMassCut.C new file mode 100644 index 000000000..4cb177bc3 --- /dev/null +++ b/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 +#include +#include +#include "TLorentzVector.h" +#include +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 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; +} \ No newline at end of file