From 39dd6d182a8d2a20f5775cb785c4cf4276e61942 Mon Sep 17 00:00:00 2001 From: Alberto Caliva Date: Wed, 14 Jan 2026 14:34:56 +0100 Subject: [PATCH 1/3] [PWGLF] Add generator for strangeness in jets --- .../GeneratorLF_StrangenessInJets_gap4.ini | 6 + .../GeneratorLF_StrangenessInJets_gap4.C | 94 ++++++++++++++ .../generator_pythia8_strangeness_in_jets.C | 118 ++++++++++++++++++ 3 files changed, 218 insertions(+) create mode 100644 MC/config/PWGLF/ini/GeneratorLF_StrangenessInJets_gap4.ini create mode 100644 MC/config/PWGLF/ini/tests/GeneratorLF_StrangenessInJets_gap4.C create mode 100644 MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C diff --git a/MC/config/PWGLF/ini/GeneratorLF_StrangenessInJets_gap4.ini b/MC/config/PWGLF/ini/GeneratorLF_StrangenessInJets_gap4.ini new file mode 100644 index 000000000..7143b919c --- /dev/null +++ b/MC/config/PWGLF/ini/GeneratorLF_StrangenessInJets_gap4.ini @@ -0,0 +1,6 @@ +[GeneratorExternal] +fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C +funcName = generateStrangeInJet(10.0,0.4,4) + +[GeneratorPythia8] +config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/pythia8/generator/pythia8_jet_ropes_136tev.cfg \ No newline at end of file diff --git a/MC/config/PWGLF/ini/tests/GeneratorLF_StrangenessInJets_gap4.C b/MC/config/PWGLF/ini/tests/GeneratorLF_StrangenessInJets_gap4.C new file mode 100644 index 000000000..fb7a162ee --- /dev/null +++ b/MC/config/PWGLF/ini/tests/GeneratorLF_StrangenessInJets_gap4.C @@ -0,0 +1,94 @@ +#include "TFile.h" +#include "TTree.h" +#include "TMath.h" +#include "fastjet/ClusterSequence.hh" +#include +#include +#include "DataFormatsMC/MCTrack.h" + +using namespace fastjet; + +int External() { + 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); + + // Jet parameters + const double ptJetThreshold = 10.0; + const double jetR = 0.4; + const int gapSize = 4; // 4 events auto-accepted, 5th needs strange-in-jet + + // Xi and Omega PDG codes + const std::vector pdgXiOmega = {3312, -3312, 3334, -3334}; + + Long64_t nEntries = tree->GetEntries(); + + for (Long64_t iEntry = 0; iEntry < nEntries; ++iEntry) { + tree->GetEntry(iEntry); + if (!tracks || tracks->empty()) continue; + + bool acceptEvent = false; + + // Gap-trigger logic + if (iEntry % (gapSize + 1) < gapSize) { + // Accept event automatically + acceptEvent = true; + std::cout << "Gap-trigger accepted event " << iEntry << " (no Xi/Omega check)\n"; + } else { + // Require Xi/Omega inside a jet + std::vector fjParticles; + std::vector fjIndexMap; + for (size_t i = 0; i < tracks->size(); ++i) { + const auto &t = tracks->at(i); + if (t.GetPdgCode() == 0) continue; + if (std::abs(t.GetEta()) > 0.8) continue; // acceptance cut + fjParticles.emplace_back(t.GetPx(), t.GetPy(), t.GetPz(), t.GetE()); + fjIndexMap.push_back(i); + } + + if (!fjParticles.empty()) { + JetDefinition jetDef(antikt_algorithm, jetR); + ClusterSequence cs(fjParticles, jetDef); + std::vector jets = sorted_by_pt(cs.inclusive_jets(ptJetThreshold)); + + for (auto &jet : jets) { + auto constituents = jet.constituents(); + for (auto &c : constituents) { + int trackIndex = fjIndexMap[c.user_index()]; + int pdg = tracks->at(trackIndex).GetPdgCode(); + if (std::find(pdgXiOmega.begin(), pdgXiOmega.end(), pdg) != pdgXiOmega.end()) { + acceptEvent = true; + std::cout << "Accepted event " << iEntry + << ": Jet pt = " << jet.pt() + << ", eta = " << jet.eta() + << ", phi = " << jet.phi() + << ", contains PDG " << pdg << "\n"; + break; + } + } + if (acceptEvent) break; + } + } + } + + if (acceptEvent) { + // Here you could break if you only want the first accepted event, + // or continue to process all events following the gap-trigger pattern + } + } + + return 0; +} diff --git a/MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C b/MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C new file mode 100644 index 000000000..61679ea6e --- /dev/null +++ b/MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C @@ -0,0 +1,118 @@ +#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 "fairlogger/Logger.h" +#include "fastjet/ClusterSequence.hh" +#include +#include +#include +#include + +using namespace Pythia8; +using namespace fastjet; +#endif + +/// Pythia8 event generator for pp collisions +/// Selection of events with Xi or Omega inside jets with pt > 10 GeV + +class GeneratorPythia8StrangeInJet : public o2::eventgen::GeneratorPythia8 { +public: + /// Constructor + GeneratorPythia8StrangeInJet(double ptJetThreshold = 10.0, double jetR = 0.4, int gapSize = 4) + : o2::eventgen::GeneratorPythia8(), + mPtJetThreshold(ptJetThreshold), + mJetR(jetR), + mGapSize(gapSize) + { + fmt::printf( + ">> Pythia8 generator: Xi/Omega inside jets with ptJet > %.1f GeV, R = %.1f, gap = %d\n", + ptJetThreshold, jetR, gapSize); + } + + ~GeneratorPythia8StrangeInJet() = default; + + bool Init() override { + addSubGenerator(0, "Pythia8 events with Xi/Omega inside jets"); + return o2::eventgen::GeneratorPythia8::Init(); + } + +protected: + bool generateEvent() override { + fmt::printf(">> Generating event %lu\n", mGeneratedEvents); + + bool genOk = false; + int localCounter{0}; + + // Accept mGapSize events unconditionally, then one triggered event + if (mGeneratedEvents % (mGapSize + 1) < mGapSize) { + genOk = GeneratorPythia8::generateEvent(); + fmt::printf(">> Gap-trigger accepted event (no strangeness check)\n"); + } else { + while (!genOk) { + if (GeneratorPythia8::generateEvent()) { + genOk = selectEvent(mPythia.event); + } + localCounter++; + } + fmt::printf(">> Event accepted after %d iterations (Xi/Omega in jet)\n", + localCounter); + } + + notifySubGenerator(0); + mGeneratedEvents++; + return true; + } + + bool selectEvent(Pythia8::Event &event) { + const std::vector pdgXiOmega = {3312, -3312, 3334, -3334}; + + std::vector fjParticles; + + for (int i = 0; i < event.size(); ++i) { + const auto &p = event[i]; + + // --- Jet input selection --- + if (!p.isFinal()) continue; + if (!p.isPrimary()) continue; + if (!p.isCharged()) continue; + if (std::abs(p.eta()) > 0.8) continue; + + PseudoJet pj(p.px(), p.py(), p.pz(), p.e()); + pj.set_user_index(i); // map back to Pythia index + fjParticles.push_back(pj); + } + + if (fjParticles.empty()) return false; + + JetDefinition jetDef(antikt_algorithm, mJetR); + ClusterSequence cs(fjParticles, jetDef); + auto jets = sorted_by_pt(cs.inclusive_jets(mPtJetThreshold)); + + for (const auto &jet : jets) { + for (const auto &c : jet.constituents()) { + int idx = c.user_index(); + int pdg = event[idx].id(); + if (std::find(pdgXiOmega.begin(), pdgXiOmega.end(), pdg) != pdgXiOmega.end()) { + fmt::printf( + ">> Accepted jet: pt = %.2f, eta = %.2f, phi = %.2f, contains PDG %d\n", + jet.pt(), jet.eta(), jet.phi(), pdg); + return true; + } + } + } + return false; + } + +private: + double mPtJetThreshold{10.0}; + double mJetR{0.4}; + int mGapSize{4}; + uint64_t mGeneratedEvents{0}; +}; From 207d32224b7db11f87d7c9884b9204e63e712f75 Mon Sep 17 00:00:00 2001 From: Alberto Caliva Date: Wed, 14 Jan 2026 16:47:29 +0100 Subject: [PATCH 2/3] replace primary with physical primary or from HF decays --- .../generator_pythia8_strangeness_in_jets.C | 66 ++++++++++++++++--- 1 file changed, 57 insertions(+), 9 deletions(-) diff --git a/MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C b/MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C index 61679ea6e..641fc6782 100644 --- a/MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C +++ b/MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C @@ -21,15 +21,18 @@ using namespace fastjet; /// Pythia8 event generator for pp collisions /// Selection of events with Xi or Omega inside jets with pt > 10 GeV +/// Jets built from physical primaries OR HF decay products class GeneratorPythia8StrangeInJet : public o2::eventgen::GeneratorPythia8 { public: /// Constructor - GeneratorPythia8StrangeInJet(double ptJetThreshold = 10.0, double jetR = 0.4, int gapSize = 4) + GeneratorPythia8StrangeInJet(double ptJetThreshold = 10.0, + double jetR = 0.4, + int gapSize = 4) : o2::eventgen::GeneratorPythia8(), - mPtJetThreshold(ptJetThreshold), - mJetR(jetR), - mGapSize(gapSize) + mPtJetThreshold(ptJetThreshold), + mJetR(jetR), + mGapSize(gapSize) { fmt::printf( ">> Pythia8 generator: Xi/Omega inside jets with ptJet > %.1f GeV, R = %.1f, gap = %d\n", @@ -44,6 +47,42 @@ public: } protected: + /// Check if particle is physical primary OR from HF decay + bool isPhysicalPrimaryOrFromHF(const Pythia8::Particle& p, + const Pythia8::Event& event) + { + // Must be final + if (!p.isFinal()) { + return false; + } + + // Physical primary: no real mother (or beam) + if (p.mother1() <= 0) { + return true; + } + + // Walk up ancestry to identify charm or beauty + int motherIdx = p.mother1(); + while (motherIdx > 0) { + const auto& mother = event[motherIdx]; + int absPdg = std::abs(mother.id()); + + // Charm or beauty hadrons + if ((absPdg / 100 == 4) || (absPdg / 100 == 5) || + (absPdg / 1000 == 4) || (absPdg / 1000 == 5)) { + return true; + } + + // Stop at beam + if (mother.mother1() <= 0) { + break; + } + motherIdx = mother.mother1(); + } + + return false; + } + bool generateEvent() override { fmt::printf(">> Generating event %lu\n", mGeneratedEvents); @@ -72,19 +111,25 @@ protected: bool selectEvent(Pythia8::Event &event) { const std::vector pdgXiOmega = {3312, -3312, 3334, -3334}; + const double mpi = 0.1395704; std::vector fjParticles; for (int i = 0; i < event.size(); ++i) { - const auto &p = event[i]; + const auto& p = event[i]; // --- Jet input selection --- if (!p.isFinal()) continue; - if (!p.isPrimary()) continue; if (!p.isCharged()) continue; + if (!isPhysicalPrimaryOrFromHF(p, event)) continue; if (std::abs(p.eta()) > 0.8) continue; - PseudoJet pj(p.px(), p.py(), p.pz(), p.e()); + double pt = std::sqrt(p.px() * p.px() + p.py() * p.py()); + if (pt < 0.1) continue; + + double energy = std::sqrt(p.p() * p.p() + mpi * mpi); + + PseudoJet pj(p.px(), p.py(), p.pz(), energy); pj.set_user_index(i); // map back to Pythia index fjParticles.push_back(pj); } @@ -95,10 +140,11 @@ protected: ClusterSequence cs(fjParticles, jetDef); auto jets = sorted_by_pt(cs.inclusive_jets(mPtJetThreshold)); - for (const auto &jet : jets) { - for (const auto &c : jet.constituents()) { + for (const auto& jet : jets) { + for (const auto& c : jet.constituents()) { int idx = c.user_index(); int pdg = event[idx].id(); + if (std::find(pdgXiOmega.begin(), pdgXiOmega.end(), pdg) != pdgXiOmega.end()) { fmt::printf( ">> Accepted jet: pt = %.2f, eta = %.2f, phi = %.2f, contains PDG %d\n", @@ -107,6 +153,7 @@ protected: } } } + return false; } @@ -116,3 +163,4 @@ private: int mGapSize{4}; uint64_t mGeneratedEvents{0}; }; + From e9a736fe4e5e59ad20e95c775e7114e78967c4a4 Mon Sep 17 00:00:00 2001 From: Alberto Caliva Date: Thu, 15 Jan 2026 07:41:39 +0100 Subject: [PATCH 3/3] added missing namespace qualifier --- .../generator_pythia8_strangeness_in_jets.C | 19 +++++++++++++------ 1 file changed, 13 insertions(+), 6 deletions(-) diff --git a/MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C b/MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C index 641fc6782..e1d2fb705 100644 --- a/MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C +++ b/MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C @@ -9,7 +9,14 @@ #include "TRandom3.h" #include "TSystem.h" #include "fairlogger/Logger.h" -#include "fastjet/ClusterSequence.hh" +#include +#include +#include +#include +#include +#include +#include +#include #include #include #include @@ -113,7 +120,7 @@ protected: const std::vector pdgXiOmega = {3312, -3312, 3334, -3334}; const double mpi = 0.1395704; - std::vector fjParticles; + std::vector fjParticles; for (int i = 0; i < event.size(); ++i) { const auto& p = event[i]; @@ -129,16 +136,16 @@ protected: double energy = std::sqrt(p.p() * p.p() + mpi * mpi); - PseudoJet pj(p.px(), p.py(), p.pz(), energy); + fastjet::PseudoJet pj(p.px(), p.py(), p.pz(), energy); pj.set_user_index(i); // map back to Pythia index fjParticles.push_back(pj); } if (fjParticles.empty()) return false; - JetDefinition jetDef(antikt_algorithm, mJetR); - ClusterSequence cs(fjParticles, jetDef); - auto jets = sorted_by_pt(cs.inclusive_jets(mPtJetThreshold)); + fastjet::JetDefinition jetDef(fastjet::antikt_algorithm, mJetR); + fastjet::ClusterSequence cs(fjParticles, jetDef); + auto jets = fastjet::sorted_by_pt(cs.inclusive_jets(mPtJetThreshold)); for (const auto& jet : jets) { for (const auto& c : jet.constituents()) {