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..641fc6782 --- /dev/null +++ b/MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C @@ -0,0 +1,166 @@ +#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 +/// 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) + : 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: + /// 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); + + 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}; + const double mpi = 0.1395704; + + 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.isCharged()) continue; + if (!isPhysicalPrimaryOrFromHF(p, event)) continue; + if (std::abs(p.eta()) > 0.8) continue; + + 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); + } + + 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}; +}; +