Skip to content

Commit 39dd6d1

Browse files
alcalivaalcaliva
committed
[PWGLF] Add generator for strangeness in jets
1 parent 6d85da1 commit 39dd6d1

File tree

3 files changed

+218
-0
lines changed

3 files changed

+218
-0
lines changed

‎MC/config/PWGLF/ini/GeneratorLF_StrangenessInJets_gap4.ini‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,6 @@
1+
[GeneratorExternal]
2+
fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C
3+
funcName = generateStrangeInJet(10.0,0.4,4)
4+
5+
[GeneratorPythia8]
6+
config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGLF/pythia8/generator/pythia8_jet_ropes_136tev.cfg

‎MC/config/PWGLF/ini/tests/GeneratorLF_StrangenessInJets_gap4.C‎

Lines changed: 94 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,94 @@
1+
#include "TFile.h"
2+
#include "TTree.h"
3+
#include "TMath.h"
4+
#include "fastjet/ClusterSequence.hh"
5+
#include <vector>
6+
#include <iostream>
7+
#include "DataFormatsMC/MCTrack.h"
8+
9+
using namespace fastjet;
10+
11+
int External() {
12+
std::string path{"o2sim_Kine.root"};
13+
14+
TFile file(path.c_str(), "READ");
15+
if (file.IsZombie()) {
16+
std::cerr << "Cannot open ROOT file " << path << "\n";
17+
return 1;
18+
}
19+
20+
auto tree = (TTree *)file.Get("o2sim");
21+
if (!tree) {
22+
std::cerr << "Cannot find tree o2sim in file " << path << "\n";
23+
return 1;
24+
}
25+
26+
std::vector<o2::MCTrack> *tracks{};
27+
tree->SetBranchAddress("MCTrack", &tracks);
28+
29+
// Jet parameters
30+
const double ptJetThreshold = 10.0;
31+
const double jetR = 0.4;
32+
const int gapSize = 4; // 4 events auto-accepted, 5th needs strange-in-jet
33+
34+
// Xi and Omega PDG codes
35+
const std::vector<int> pdgXiOmega = {3312, -3312, 3334, -3334};
36+
37+
Long64_t nEntries = tree->GetEntries();
38+
39+
for (Long64_t iEntry = 0; iEntry < nEntries; ++iEntry) {
40+
tree->GetEntry(iEntry);
41+
if (!tracks || tracks->empty()) continue;
42+
43+
bool acceptEvent = false;
44+
45+
// Gap-trigger logic
46+
if (iEntry % (gapSize + 1) < gapSize) {
47+
// Accept event automatically
48+
acceptEvent = true;
49+
std::cout << "Gap-trigger accepted event " << iEntry << " (no Xi/Omega check)\n";
50+
} else {
51+
// Require Xi/Omega inside a jet
52+
std::vector<PseudoJet> fjParticles;
53+
std::vector<int> fjIndexMap;
54+
for (size_t i = 0; i < tracks->size(); ++i) {
55+
const auto &t = tracks->at(i);
56+
if (t.GetPdgCode() == 0) continue;
57+
if (std::abs(t.GetEta()) > 0.8) continue; // acceptance cut
58+
fjParticles.emplace_back(t.GetPx(), t.GetPy(), t.GetPz(), t.GetE());
59+
fjIndexMap.push_back(i);
60+
}
61+
62+
if (!fjParticles.empty()) {
63+
JetDefinition jetDef(antikt_algorithm, jetR);
64+
ClusterSequence cs(fjParticles, jetDef);
65+
std::vector<PseudoJet> jets = sorted_by_pt(cs.inclusive_jets(ptJetThreshold));
66+
67+
for (auto &jet : jets) {
68+
auto constituents = jet.constituents();
69+
for (auto &c : constituents) {
70+
int trackIndex = fjIndexMap[c.user_index()];
71+
int pdg = tracks->at(trackIndex).GetPdgCode();
72+
if (std::find(pdgXiOmega.begin(), pdgXiOmega.end(), pdg) != pdgXiOmega.end()) {
73+
acceptEvent = true;
74+
std::cout << "Accepted event " << iEntry
75+
<< ": Jet pt = " << jet.pt()
76+
<< ", eta = " << jet.eta()
77+
<< ", phi = " << jet.phi()
78+
<< ", contains PDG " << pdg << "\n";
79+
break;
80+
}
81+
}
82+
if (acceptEvent) break;
83+
}
84+
}
85+
}
86+
87+
if (acceptEvent) {
88+
// Here you could break if you only want the first accepted event,
89+
// or continue to process all events following the gap-trigger pattern
90+
}
91+
}
92+
93+
return 0;
94+
}

‎MC/config/PWGLF/pythia8/generator_pythia8_strangeness_in_jets.C‎

Lines changed: 118 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,118 @@
1+
#if !defined(__CLING__) || defined(__ROOTCLING__)
2+
#include "FairGenerator.h"
3+
#include "FairPrimaryGenerator.h"
4+
#include "Generators/GeneratorPythia8.h"
5+
#include "Pythia8/Pythia.h"
6+
#include "TDatabasePDG.h"
7+
#include "TMath.h"
8+
#include "TParticlePDG.h"
9+
#include "TRandom3.h"
10+
#include "TSystem.h"
11+
#include "fairlogger/Logger.h"
12+
#include "fastjet/ClusterSequence.hh"
13+
#include <cmath>
14+
#include <fstream>
15+
#include <string>
16+
#include <vector>
17+
18+
using namespace Pythia8;
19+
using namespace fastjet;
20+
#endif
21+
22+
/// Pythia8 event generator for pp collisions
23+
/// Selection of events with Xi or Omega inside jets with pt > 10 GeV
24+
25+
class GeneratorPythia8StrangeInJet : public o2::eventgen::GeneratorPythia8 {
26+
public:
27+
/// Constructor
28+
GeneratorPythia8StrangeInJet(double ptJetThreshold = 10.0, double jetR = 0.4, int gapSize = 4)
29+
: o2::eventgen::GeneratorPythia8(),
30+
mPtJetThreshold(ptJetThreshold),
31+
mJetR(jetR),
32+
mGapSize(gapSize)
33+
{
34+
fmt::printf(
35+
">> Pythia8 generator: Xi/Omega inside jets with ptJet > %.1f GeV, R = %.1f, gap = %d\n",
36+
ptJetThreshold, jetR, gapSize);
37+
}
38+
39+
~GeneratorPythia8StrangeInJet() = default;
40+
41+
bool Init() override {
42+
addSubGenerator(0, "Pythia8 events with Xi/Omega inside jets");
43+
return o2::eventgen::GeneratorPythia8::Init();
44+
}
45+
46+
protected:
47+
bool generateEvent() override {
48+
fmt::printf(">> Generating event %lu\n", mGeneratedEvents);
49+
50+
bool genOk = false;
51+
int localCounter{0};
52+
53+
// Accept mGapSize events unconditionally, then one triggered event
54+
if (mGeneratedEvents % (mGapSize + 1) < mGapSize) {
55+
genOk = GeneratorPythia8::generateEvent();
56+
fmt::printf(">> Gap-trigger accepted event (no strangeness check)\n");
57+
} else {
58+
while (!genOk) {
59+
if (GeneratorPythia8::generateEvent()) {
60+
genOk = selectEvent(mPythia.event);
61+
}
62+
localCounter++;
63+
}
64+
fmt::printf(">> Event accepted after %d iterations (Xi/Omega in jet)\n",
65+
localCounter);
66+
}
67+
68+
notifySubGenerator(0);
69+
mGeneratedEvents++;
70+
return true;
71+
}
72+
73+
bool selectEvent(Pythia8::Event &event) {
74+
const std::vector<int> pdgXiOmega = {3312, -3312, 3334, -3334};
75+
76+
std::vector<PseudoJet> fjParticles;
77+
78+
for (int i = 0; i < event.size(); ++i) {
79+
const auto &p = event[i];
80+
81+
// --- Jet input selection ---
82+
if (!p.isFinal()) continue;
83+
if (!p.isPrimary()) continue;
84+
if (!p.isCharged()) continue;
85+
if (std::abs(p.eta()) > 0.8) continue;
86+
87+
PseudoJet pj(p.px(), p.py(), p.pz(), p.e());
88+
pj.set_user_index(i); // map back to Pythia index
89+
fjParticles.push_back(pj);
90+
}
91+
92+
if (fjParticles.empty()) return false;
93+
94+
JetDefinition jetDef(antikt_algorithm, mJetR);
95+
ClusterSequence cs(fjParticles, jetDef);
96+
auto jets = sorted_by_pt(cs.inclusive_jets(mPtJetThreshold));
97+
98+
for (const auto &jet : jets) {
99+
for (const auto &c : jet.constituents()) {
100+
int idx = c.user_index();
101+
int pdg = event[idx].id();
102+
if (std::find(pdgXiOmega.begin(), pdgXiOmega.end(), pdg) != pdgXiOmega.end()) {
103+
fmt::printf(
104+
">> Accepted jet: pt = %.2f, eta = %.2f, phi = %.2f, contains PDG %d\n",
105+
jet.pt(), jet.eta(), jet.phi(), pdg);
106+
return true;
107+
}
108+
}
109+
}
110+
return false;
111+
}
112+
113+
private:
114+
double mPtJetThreshold{10.0};
115+
double mJetR{0.4};
116+
int mGapSize{4};
117+
uint64_t mGeneratedEvents{0};
118+
};

0 commit comments

Comments
 (0)