From 39dd6d182a8d2a20f5775cb785c4cf4276e61942 Mon Sep 17 00:00:00 2001 From: Alberto Caliva Date: Wed, 14 Jan 2026 14:34:56 +0100 Subject: [PATCH 01/10] [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 02/10] 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 03/10] 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()) { From 332a2ee52ebb49b0e5ad5f974e6f4c8eae308d09 Mon Sep 17 00:00:00 2001 From: Alberto Caliva Date: Thu, 15 Jan 2026 10:25:20 +0100 Subject: [PATCH 04/10] fix headers --- .../PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C | 7 +++++-- 1 file changed, 5 insertions(+), 2 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 e1d2fb705..19392d6d9 100644 --- a/MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C +++ b/MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C @@ -1,14 +1,17 @@ -#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 #include #include @@ -17,6 +20,7 @@ #include #include #include + #include #include #include @@ -24,7 +28,6 @@ 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 From f2601c55e267e592cdcd7d6013a2c4ba98b7cd75 Mon Sep 17 00:00:00 2001 From: Alberto Caliva Date: Thu, 15 Jan 2026 15:03:31 +0100 Subject: [PATCH 05/10] fix selection of physical primary particles --- .../generator_pythia8_strangeness_in_jets.C | 30 ++++++++++--------- 1 file changed, 16 insertions(+), 14 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 19392d6d9..fb92f89be 100644 --- a/MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C +++ b/MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C @@ -61,36 +61,38 @@ protected: 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; + const int absPdg = std::abs(p.id()); + + // Particle species selection + const bool isAcceptedSpecies = (absPdg == 211 || absPdg == 321 || absPdg == 2212 || absPdg == 1000010020 || absPdg == 11 || absPdg == 13); + + if (!isAcceptedSpecies) { + return false; } - // Walk up ancestry to identify charm or beauty + // Walk up ancestry int motherIdx = p.mother1(); + while (motherIdx > 0) { const auto& mother = event[motherIdx]; - int absPdg = std::abs(mother.id()); + const int absMotherPdg = std::abs(mother.id()); - // Charm or beauty hadrons - if ((absPdg / 100 == 4) || (absPdg / 100 == 5) || - (absPdg / 1000 == 4) || (absPdg / 1000 == 5)) { + // Charm or beauty hadron → accept (HF decay) + if ((absMotherPdg / 100 == 4) || (absMotherPdg / 100 == 5) || (absMotherPdg / 1000 == 4) || (absMotherPdg / 1000 == 5)) { return true; } - // Stop at beam - if (mother.mother1() <= 0) { - break; + // Weakly decaying hadron → reject (non-physical primary) + if (mother.isHadron() && mother.tau0() > 1.0) { + return false; } motherIdx = mother.mother1(); } - - return false; + return true; } bool generateEvent() override { From 30ca6bf09186a9d2dc5605807cdc57f629a2dd92 Mon Sep 17 00:00:00 2001 From: Alberto Caliva Date: Thu, 15 Jan 2026 20:48:17 +0100 Subject: [PATCH 06/10] remove redundant includes --- .../generator_pythia8_strangeness_in_jets.C | 21 ++++--------------- 1 file changed, 4 insertions(+), 17 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 fb92f89be..7692a57d3 100644 --- a/MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C +++ b/MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C @@ -4,27 +4,14 @@ #include "Pythia8/Pythia.h" -#include "TDatabasePDG.h" -#include "TMath.h" -#include "TParticlePDG.h" -#include "TRandom3.h" -#include "TSystem.h" - -#include "fairlogger/Logger.h" - -#include -#include -#include -#include #include -#include -#include -#include +#include +#include #include -#include -#include #include +#include +#include using namespace Pythia8; using namespace fastjet; From caad66015331f1015040a2794660d77e422b4ee3 Mon Sep 17 00:00:00 2001 From: Alberto Caliva Date: Fri, 16 Jan 2026 11:37:56 +0100 Subject: [PATCH 07/10] change jet definition --- .../GeneratorLF_StrangenessInJets_gap4.C | 81 +-------- .../generator_pythia8_strangeness_in_jets.C | 172 ++++++++++-------- 2 files changed, 102 insertions(+), 151 deletions(-) diff --git a/MC/config/PWGLF/ini/tests/GeneratorLF_StrangenessInJets_gap4.C b/MC/config/PWGLF/ini/tests/GeneratorLF_StrangenessInJets_gap4.C index fb7a162ee..12f069f3c 100644 --- a/MC/config/PWGLF/ini/tests/GeneratorLF_StrangenessInJets_gap4.C +++ b/MC/config/PWGLF/ini/tests/GeneratorLF_StrangenessInJets_gap4.C @@ -1,13 +1,3 @@ -#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"}; @@ -22,73 +12,16 @@ int External() { 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 - } + auto nXi = tree->Scan("MCTrack.GetPdgCode()", "TMath::Abs(MCTrack.GetPdgCode()) == 3312"); + auto nOmega = tree->Scan("MCTrack.GetPdgCode()", "TMath::Abs(MCTrack.GetPdgCode()) == 3334"); + auto nStr = nXi+nOmega; + + if (nStr == 0) { + std::cerr << "No event of interest\n"; + return 1; } - 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 index 7692a57d3..fd9d8e758 100644 --- a/MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C +++ b/MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C @@ -1,42 +1,44 @@ +#if !defined(__CLING__) || defined(__ROOTCLING__) #include "FairGenerator.h" #include "FairPrimaryGenerator.h" #include "Generators/GeneratorPythia8.h" - #include "Pythia8/Pythia.h" - -#include -#include -#include - +#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 -#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 -/// Jets built from physical primaries OR HF decay products +/// Event generator using Pythia ropes for Ξ and Ω production in jets. +/// Jets are defined as particles within a cone around the highest-pT particle +/// (approximate jet axis) and include charged final-state particles that are +/// either physical primaries or from heavy-flavor decays. +/// Jets must be fully within the acceptance and have pT > 10 GeV. +/// One event of interest is generated after 4 minimum-bias events. -class GeneratorPythia8StrangeInJet : public o2::eventgen::GeneratorPythia8 { +class GeneratorPythia8StrangenessInJet : public o2::eventgen::GeneratorPythia8 { public: /// Constructor - GeneratorPythia8StrangeInJet(double ptJetThreshold = 10.0, - double jetR = 0.4, - int gapSize = 4) + GeneratorPythia8StrangenessInJet(double ptJetMin = 10.0, double Rjet = 0.4, int gapSize = 4) : o2::eventgen::GeneratorPythia8(), - mPtJetThreshold(ptJetThreshold), - mJetR(jetR), + mptJetMin(ptJetMin), + mRjet(Rjet), mGapSize(gapSize) { fmt::printf( - ">> Pythia8 generator: Xi/Omega inside jets with ptJet > %.1f GeV, R = %.1f, gap = %d\n", - ptJetThreshold, jetR, gapSize); + ">> Pythia8 generator: Xi/Omega inside jets with pt_{jet} > %.1f GeV, R_{jet} = %.1f, gap = %d\n", + ptJetMin, Rjet, gapSize); } - - ~GeneratorPythia8StrangeInJet() = default; + /// Destructor + ~GeneratorPythia8StrangenessInJet() = default; bool Init() override { addSubGenerator(0, "Pythia8 events with Xi/Omega inside jets"); @@ -44,63 +46,81 @@ public: } protected: - /// Check if particle is physical primary OR from HF decay - bool isPhysicalPrimaryOrFromHF(const Pythia8::Particle& p, - const Pythia8::Event& event) + /// Check if particle is physical primary or from HF decay + bool isPhysicalPrimaryOrFromHF(const Pythia8::Particle& p, const Pythia8::Event& event) { + // Select only final-state particles if (!p.isFinal()) { return false; } - const int absPdg = std::abs(p.id()); - // Particle species selection - const bool isAcceptedSpecies = (absPdg == 211 || absPdg == 321 || absPdg == 2212 || absPdg == 1000010020 || absPdg == 11 || absPdg == 13); - - if (!isAcceptedSpecies) { + const int absPdg = std::abs(p.id()); + if (absPdg!=211 && absPdg!=321 && absPdg!= 2212 && absPdg!=1000010020 && absPdg!=11 && absPdg!=13) return false; - } // Walk up ancestry - int motherIdx = p.mother1(); + int motherId = p.mother1(); - while (motherIdx > 0) { - const auto& mother = event[motherIdx]; + while (motherId > 0) { + + // Get mother + const auto& mother = event[motherId]; const int absMotherPdg = std::abs(mother.id()); - // Charm or beauty hadron → accept (HF decay) + // Check if particle is from HF decay if ((absMotherPdg / 100 == 4) || (absMotherPdg / 100 == 5) || (absMotherPdg / 1000 == 4) || (absMotherPdg / 1000 == 5)) { return true; } - // Weakly decaying hadron → reject (non-physical primary) + // Reject non-physical primary hadrons if (mother.isHadron() && mother.tau0() > 1.0) { return false; } - motherIdx = mother.mother1(); + motherId = mother.mother1(); } return true; } + // Compute delta phi + double getDeltaPhi(double a1, double a2) + { + double deltaPhi{0.0}; + double phi1 = TVector2::Phi_0_2pi(a1); + double phi2 = TVector2::Phi_0_2pi(a2); + double diff = std::fabs(phi1 - phi2); + + if (diff <= M_PI) + deltaPhi = diff; + if (diff > M_PI) + deltaPhi = 2.0 * M_PI - diff; + + return deltaPhi; + } + bool generateEvent() override { - fmt::printf(">> Generating event %lu\n", mGeneratedEvents); + fmt::printf(">> Generating event %d\n", mGeneratedEvents); bool genOk = false; int localCounter{0}; + constexpr int kMaxTries{100000}; // 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) { + while (!genOk && localCounter < kMaxTries) { if (GeneratorPythia8::generateEvent()) { genOk = selectEvent(mPythia.event); } localCounter++; } - fmt::printf(">> Event accepted after %d iterations (Xi/Omega in jet)\n", - localCounter); + if (!genOk) { + fmt::printf("Failed to generate triggered event after %d tries\n",kMaxTries); + return false; + } + fmt::printf(">> Event accepted after %d iterations (Xi/Omega in jet)\n", localCounter); } notifySubGenerator(0); @@ -109,56 +129,54 @@ protected: } bool selectEvent(Pythia8::Event &event) { - const std::vector pdgXiOmega = {3312, -3312, 3334, -3334}; - const double mpi = 0.1395704; - std::vector fjParticles; + double etaJet{-999.0}, phiJet{-999.0}, ptMax{0.0}; + std::vector particleID; + bool containsXiOrOmega{false}; - for (int i = 0; i < event.size(); ++i) { + for (int i = 0; i < event.size(); i++) { const auto& p = event[i]; + if (std::abs(p.id()) == 3312 || std::abs(p.id()) == 3334) + containsXiOrOmega = true; - // --- Jet input selection --- - if (!p.isFinal()) continue; - if (!p.isCharged()) continue; + if (std::abs(p.eta()) > 0.8 || p.pT() < 0.1) continue; if (!isPhysicalPrimaryOrFromHF(p, event)) continue; - if (std::abs(p.eta()) > 0.8) continue; + particleID.emplace_back(i); - 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); - - 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 (p.pT() > ptMax) { + ptMax = p.pT(); + etaJet = p.eta(); + phiJet = p.phi(); + } } + if (ptMax == 0.0) + return false; + if (std::abs(etaJet) + Rjet > 0.8) + return false; + if (!containsXiOrOmega) + return false; - if (fjParticles.empty()) return false; - - 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()) { - int idx = c.user_index(); - int pdg = event[idx].id(); + double ptJet{0.0}; + for (int i = 0 ; i < particleID.size() ; i++) { + int id = particleID[i]; + const auto& p = event[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; - } + double deltaEta = std::abs(p.eta() - etaJet); + double deltaPhi = getDeltaPhi(p.phi(),phiJet); + double deltaR = std::sqrt(deltaEta * deltaEta + deltaPhi * deltaPhi); + if (deltaR < Rjet) { + ptJet += p.pT(); } } + if (ptJet < ptJetMin) + return false; - return false; + return true; } private: - double mPtJetThreshold{10.0}; - double mJetR{0.4}; + double mptJetMin{10.0}; + double mRjet{0.4}; int mGapSize{4}; uint64_t mGeneratedEvents{0}; }; From 0aac6c615331ffd8e17694ffb669628484b22975 Mon Sep 17 00:00:00 2001 From: Alberto Caliva Date: Fri, 16 Jan 2026 12:57:46 +0100 Subject: [PATCH 08/10] fix typo --- .../PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C | 6 +++--- 1 file changed, 3 insertions(+), 3 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 fd9d8e758..205b915a6 100644 --- a/MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C +++ b/MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C @@ -151,7 +151,7 @@ protected: } if (ptMax == 0.0) return false; - if (std::abs(etaJet) + Rjet > 0.8) + if (std::abs(etaJet) + mRjet > 0.8) return false; if (!containsXiOrOmega) return false; @@ -164,11 +164,11 @@ protected: double deltaEta = std::abs(p.eta() - etaJet); double deltaPhi = getDeltaPhi(p.phi(),phiJet); double deltaR = std::sqrt(deltaEta * deltaEta + deltaPhi * deltaPhi); - if (deltaR < Rjet) { + if (deltaR < mRjet) { ptJet += p.pT(); } } - if (ptJet < ptJetMin) + if (ptJet < mptJetMin) return false; return true; From bf6fa094a245fb8a171161aae9b4ad739088f4d6 Mon Sep 17 00:00:00 2001 From: Alberto Caliva Date: Fri, 16 Jan 2026 22:38:23 +0100 Subject: [PATCH 09/10] fix typo --- MC/config/PWGLF/ini/GeneratorLF_StrangenessInJets_gap4.ini | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/MC/config/PWGLF/ini/GeneratorLF_StrangenessInJets_gap4.ini b/MC/config/PWGLF/ini/GeneratorLF_StrangenessInJets_gap4.ini index 7143b919c..005c27ff3 100644 --- a/MC/config/PWGLF/ini/GeneratorLF_StrangenessInJets_gap4.ini +++ b/MC/config/PWGLF/ini/GeneratorLF_StrangenessInJets_gap4.ini @@ -1,6 +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) +funcName = GeneratorPythia8StrangenessInJet(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 From a5fc840a1ebc75f97e68d64886a9f5cdc5fff5b6 Mon Sep 17 00:00:00 2001 From: Alberto Caliva Date: Fri, 16 Jan 2026 22:52:28 +0100 Subject: [PATCH 10/10] add missing function --- .../PWGLF/ini/GeneratorLF_StrangenessInJets_gap4.ini | 2 +- .../pythia8/generator_pythia8_strangeness_in_jets.C | 9 +++++++++ 2 files changed, 10 insertions(+), 1 deletion(-) diff --git a/MC/config/PWGLF/ini/GeneratorLF_StrangenessInJets_gap4.ini b/MC/config/PWGLF/ini/GeneratorLF_StrangenessInJets_gap4.ini index 005c27ff3..80b0235b8 100644 --- a/MC/config/PWGLF/ini/GeneratorLF_StrangenessInJets_gap4.ini +++ b/MC/config/PWGLF/ini/GeneratorLF_StrangenessInJets_gap4.ini @@ -1,6 +1,6 @@ [GeneratorExternal] fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C -funcName = GeneratorPythia8StrangenessInJet(10.0,0.4,4) +funcName = generateStrangenessInJets(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/pythia8/generator_pythia8_strangeness_in_jets.C b/MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C index 205b915a6..8075d27cb 100644 --- a/MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C +++ b/MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C @@ -181,3 +181,12 @@ private: uint64_t mGeneratedEvents{0}; }; +///___________________________________________________________ +FairGenerator *generateStrangenessInJets(double ptJet = 10.0, double rJet = 0.4, int gap = 4) { + + auto myGenerator = new GeneratorPythia8StrangenessInJet(ptJet,rJet,gap); + auto seed = (gRandom->TRandom::GetSeed() % 900000000); + myGenerator->readString("Random:setSeed on"); + myGenerator->readString("Random:seed " + std::to_string(seed)); + return myGenerator; +}