Skip to content

Commit fed51ce

Browse files
jackal1-66jackal1-66
committed
Implemented Performance Generator
1 parent 1c4641e commit fed51ce

File tree

6 files changed

+410
-0
lines changed

6 files changed

+410
-0
lines changed

‎MC/config/common/external/generator/perfConf.json‎

Lines changed: 23 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,23 @@
1+
{
2+
"generators": [
3+
{
4+
"cocktail": [
5+
{
6+
"name": "pythia8pp",
7+
"config": ""
8+
},
9+
{
10+
"name": "external",
11+
"config": {
12+
"fileName": "${O2DPG_MC_CONFIG_ROOT}/MC/config/common/external/generator/performanceGenerator.C",
13+
"funcName": "Generator_Performance()",
14+
"iniFile": ""
15+
}
16+
}
17+
]
18+
}
19+
],
20+
"fractions": [
21+
1
22+
]
23+
}

‎MC/config/common/external/generator/performanceGenerator.C‎

Lines changed: 289 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,289 @@
1+
// External generator requested in https://its.cern.ch/jira/browse/O2-6235
2+
// for multidimensional performance studies
3+
namespace o2
4+
{
5+
namespace eventgen
6+
{
7+
8+
class GenPerf : public Generator
9+
{
10+
public:
11+
GenPerf(float fraction = 0.03f, unsigned short int nsig = 100, unsigned short int tag = 1)
12+
{
13+
if (fraction == -1) {
14+
LOG(info) << nsig << " Signal particles will be generated in each event";
15+
mNSig = nsig;
16+
mFraction = -1.f;
17+
} else if (fraction >= 0) {
18+
LOG(info) << "Fraction based signal generation is enabled";
19+
LOG(info) << fraction << "*nUE tracks per event will be generated";
20+
mFraction = fraction;
21+
} else {
22+
LOG(fatal) << "Wrong fraction selected. Accepted values are:";
23+
LOG(fatal) << "\t -1 => fixed number of tracks per event";
24+
LOG(fatal) << ">=0 => fraction based signal generation over the number of UE tracks per event";
25+
exit(1);
26+
}
27+
initGenMap();
28+
if (genMap.find(tag) == genMap.end()) {
29+
LOG(fatal) << "Wrong tag selected. Accepted values are:";
30+
for (const auto& [key, _] : genMap) {
31+
LOG(fatal) << "\t" << key;
32+
}
33+
exit(1);
34+
} else {
35+
mTag = tag;
36+
LOG(info) << "Generator with tag " << mTag << " is selected";
37+
}
38+
Generator::setTimeUnit(1.0);
39+
Generator::setPositionUnit(1.0);
40+
}
41+
42+
Bool_t generateEvent() override
43+
{
44+
return kTRUE;
45+
}
46+
47+
Bool_t importParticles() override
48+
{
49+
mNUE = 0;
50+
if ( mFraction != -1) {
51+
// This line assumes that the current generator is run in a cocktail with another generator
52+
// which is run before the current one in a sequential way
53+
mNUE = genList.front()->getParticles().size();
54+
LOG(debug) << "Number of tracks from UE is " << mNUE;
55+
}
56+
unsigned short nSig = (mFraction == -1) ? mNSig : std::lround(mFraction * mNUE);
57+
LOG(debug) << "Generating additional " << nSig << " particles";
58+
for (int k = 0; k < nSig; k++){
59+
mParticles.push_back(genMap[mTag]());
60+
}
61+
return kTRUE;
62+
}
63+
64+
private:
65+
float mFraction = 0.03f; // Fraction based generation
66+
unsigned short int mNSig = 0; // Number of particles to generate
67+
unsigned int mNUE = 0; // Number of tracks in the Underlying event
68+
unsigned short int mTag = 1; // Tag to select the generation function
69+
const std::vector<std::shared_ptr<o2::eventgen::Generator>> genList = GeneratorHybrid::getGenerators();
70+
std::map<unsigned short int, std::function<TParticle()>> genMap;
71+
UInt_t mGenID = 42;
72+
73+
// This is performance test generator with uniform weighting for PDG
74+
TParticle generateParticle0()
75+
{
76+
// 1. Get the singleton instances
77+
TDatabasePDG *pdg = TDatabasePDG::Instance();
78+
// 2. Define the list of PDG codes
79+
const int ncodes = 13;
80+
const int pdgCodes[ncodes] = {
81+
310, // K0_s
82+
421, // D0
83+
3122, // Lambda
84+
-3122, // Anti-Lambda
85+
443, // J/psi
86+
13, // mu-
87+
22, // gamma
88+
23, // Z0
89+
1, 2, 3, 4, 5 // Quarks: d, u, s, c, b (t-quark is 6, often excluded for kinematics)
90+
};
91+
// 3. Randomly select and validate a PDG code
92+
// TMath::Nint(gRandom->Rndm() * ncodes) selects an index from 0 to ncodes-1 safely.
93+
int index = TMath::Nint(gRandom->Rndm() * (ncodes - 1));
94+
int pdgCode = pdgCodes[index];
95+
// Check if the particle exists and switch to antiparticle if needed
96+
if (pdg->GetParticle(pdgCode) == nullptr)
97+
{
98+
if (pdg->GetParticle(-pdgCode) != nullptr)
99+
{
100+
pdgCode *= -1; // Use the negative code (antiparticle)
101+
}
102+
else
103+
{
104+
LOG(error) << "Error: PDG code " << pdgCode << " not found in TDatabasePDG. Using Muon (13).";
105+
pdgCode = 13;
106+
}
107+
}
108+
// 4. Generate Kinematics (p_T, phi, eta)
109+
float pt = 1 / (gRandom->Rndm()); // flat 1/pt distribution
110+
float phi = gRandom->Rndm() * 2.0f * TMath::Pi();
111+
float eta = 3.0f * (gRandom->Rndm() - 0.5f); // eta from -1.5 to 1.5
112+
// Initial position (origin)
113+
float xyz[3] = {0.0f, 0.0f, 0.0f};
114+
// if cosmic, you might want to randomize the vertex position
115+
if (pdgCode == 13 || pdgCode == -13)
116+
{
117+
xyz[0] = (gRandom->Rndm() - 0.5) * 300.0f; // x from -100 to 100 cm
118+
xyz[1] = (gRandom->Rndm() - 0.5) * 300.0f; // y from -100 to 100 cm
119+
xyz[2] = 400;
120+
pt = 1 / (gRandom->Rndm() + 0.01);
121+
eta = gRandom->Gaus() * 0.2;
122+
}
123+
//
124+
// Convert spherical coordinates (pt, phi, eta) to Cartesian (px, py, pz)
125+
float pz = pt * TMath::SinH(eta);
126+
float px = pt * TMath::Cos(phi);
127+
float py = pt * TMath::Sin(phi);
128+
// 5. Calculate Energy (E) from Mass (M)
129+
TParticlePDG *particleInfo = pdg->GetParticle(pdgCode);
130+
double mass = particleInfo ? particleInfo->Mass() : 0.1056; // Default to muon mass if lookup fails
131+
double energy = TMath::Sqrt(px * px + py * py + pz * pz + mass * mass);
132+
133+
// 6. Create and return the TParticle object by value
134+
// TParticle(pdgCode, trackIndex, Mother, Daughter1, Daughter2, Px, Py, Pz, E, Vx, Vy, Vz, Time)
135+
int status = -1; // Status code, -1 for undefined
136+
// Set your custom performance generator ID (e.g., ID 42)
137+
TParticle generatedParticle(pdgCode, status, -1, -1, -1, -1, px, py, pz, energy, xyz[0], xyz[1], xyz[2], 0.0);
138+
generatedParticle.SetStatusCode(o2::mcgenstatus::MCGenStatusEncoding(generatedParticle.GetStatusCode(), 0).fullEncoding);
139+
generatedParticle.SetUniqueID(mGenID);
140+
generatedParticle.SetBit(ParticleStatus::kToBeDone, //
141+
o2::mcgenstatus::getHepMCStatusCode(generatedParticle.GetStatusCode()) == 1);
142+
return generatedParticle;
143+
}
144+
145+
// Particle configuration for ALICE O2 performance testing
146+
struct ParticleSpec
147+
{
148+
int pdgCode;
149+
float fraction; // Relative probability for probe statistics
150+
float pTScale; // Scales pt
151+
};
152+
153+
// Optimized for rare probes (J/psi, D0, jets) with flat distributions
154+
const std::vector<ParticleSpec> g_particle_specs = {
155+
// PDG | Fraction | pTScale
156+
{22, 1.0f, 1.0f}, // Photon: High yield for PID/calo
157+
{13, 1.f, 1.0f}, // Muon: Cosmic override applied
158+
{-13, 1.f, 1.0f}, // Anti-muon
159+
{23, 0.1f, 10.0f}, // Z0: Rare,
160+
{310, 1.f, 1.0f}, // K0_s: Common hadron
161+
{421, 0.2f, 1.5f}, // D0
162+
{443, 0.1f, 5.0f}, // J/psi: Boosted for candle
163+
{3122, 0.5f, 1.0f}, // Lambda
164+
{-3122, 0.5f, 1.0f}, // Anti-Lambda
165+
{211, 1.0f, 1.0f}, // Pi+
166+
{-211, 1.0f, 1.0f}, // Pi-:
167+
//
168+
{21, 0.1f, 3.0f}, // Gluon: Jet proxy (status=11)
169+
{1, 0.1f, 3.0f}, // d quark: Jet proxy
170+
{-1, 0.1f, 3.0f}, // anti-d
171+
{2, 0.1f, 3.0f}, // u quark: Jet proxy
172+
{-2, 0.1f, 3.0f}, // anti-u
173+
{3, 0.1f, 5.0f}, // s quark: Strange
174+
{-3, 0.1f, 5.0f}, // anti-s
175+
{4, 0.1f, 5.0f}, // c quark: Heavy flavor
176+
{-4, 0.1f, 5.0f}, // anti-c
177+
{5, 0.1f, 8.0f}, // b quark: Very hard
178+
{-5, 0.1f, 8.0f} // anti-b
179+
};
180+
181+
// pT bounds: Max pT ~5 TeV (ALICE Pb-Pb energy)
182+
const float kMaxInvPt = 1.0f; // Min pT = 1 GeV
183+
const float kBaseMinInvPt = 2e-4f; // Max pT = 5000 GeV (unscaled)
184+
185+
// Check if particle is a parton (quark/gluon, status=11)
186+
bool isParton(int& pdgCode)
187+
{
188+
int absCode = TMath::Abs(pdgCode);
189+
return (absCode >= 1 && absCode <= 5) || absCode == 21;
190+
}
191+
192+
// Generator for flat distributions in pT, eta for calibration
193+
TParticle generateParticle1()
194+
{
195+
TDatabasePDG *pdg = TDatabasePDG::Instance();
196+
// 1. Weighted Random Selection
197+
static float totalWeight = 0.0f;
198+
if (totalWeight == 0.0f)
199+
{
200+
totalWeight = std::accumulate(g_particle_specs.begin(), g_particle_specs.end(), 0.0f,
201+
[](float sum, const ParticleSpec &spec)
202+
{ return sum + spec.fraction; });
203+
}
204+
float randVal = gRandom->Rndm() * totalWeight;
205+
float cumulativeWeight = 0.0f;
206+
const ParticleSpec *selectedSpec = nullptr;
207+
for (const auto &spec : g_particle_specs)
208+
{
209+
cumulativeWeight += spec.fraction;
210+
if (randVal <= cumulativeWeight)
211+
{
212+
selectedSpec = &spec;
213+
break;
214+
}
215+
}
216+
if (!selectedSpec)
217+
selectedSpec = &g_particle_specs.back();
218+
int pdgCode = selectedSpec->pdgCode;
219+
float pTScale = selectedSpec->pTScale;
220+
// 2. PDG Validation
221+
if (!pdg->GetParticle(pdgCode))
222+
{
223+
if (pdg->GetParticle(-pdgCode))
224+
pdgCode *= -1;
225+
else
226+
{
227+
LOG(error) << "Error: PDG " << pdgCode << " not found. Using muon (13).\n";
228+
pdgCode = 13;
229+
pTScale = 1.0f;
230+
}
231+
}
232+
// 3. Status: 11 for partons (jets), 1 for final-state
233+
int status = isParton(pdgCode) ? 11 : 1;
234+
// 4. Kinematics (flat 1/pT, max ~5000 GeV / pTScale)
235+
float min_inv_pt = kBaseMinInvPt / pTScale; // E.g., max pT=40,000 GeV for b quarks
236+
float inv_pt = (gRandom->Rndm() / pTScale) * (kMaxInvPt - min_inv_pt) + min_inv_pt;
237+
float pt = 1.0f / inv_pt;
238+
float phi = gRandom->Rndm() * 2.0f * TMath::Pi();
239+
float eta = gRandom->Rndm() * 3.0f - 1.5f; // ALICE TPC: -1.5 to 1.5
240+
// Vertex: Delta (embedding handles smearing)
241+
float xyz[3] = {0.0f, 0.0f, 0.0f};
242+
// 5. Cosmic Muon Override
243+
if (TMath::Abs(pdgCode) == 13)
244+
{
245+
xyz[0] = (gRandom->Rndm() - 0.5f) * 300.0f;
246+
xyz[1] = (gRandom->Rndm() - 0.5f) * 300.0f;
247+
xyz[2] = 400.0f;
248+
inv_pt = (gRandom->Rndm() + 0.01f) / pTScale; // Apply pTScale
249+
pt = 1.0f / inv_pt;
250+
eta = TMath::Max(-4.0, TMath::Min(4.0, gRandom->Gaus(0.0, 0.2)));
251+
status = 1;
252+
}
253+
// 6. Momentum and Energy
254+
float pz = pt * TMath::SinH(eta);
255+
float px = pt * TMath::Cos(phi);
256+
float py = pt * TMath::Sin(phi);
257+
TParticlePDG *particleInfo = pdg->GetParticle(pdgCode);
258+
double mass = particleInfo ? particleInfo->Mass() : 0.1056;
259+
double energy = TMath::Sqrt(px * px + py * py + pz * pz + mass * mass);
260+
// 7. TParticle Creation (quarks/gluons need fragmentation in O2)
261+
TParticle generatedParticle(pdgCode, status, -1, -1, -1, -1, px, py, pz, energy, xyz[0], xyz[1], xyz[2], 0.0);
262+
generatedParticle.SetStatusCode(o2::mcgenstatus::MCGenStatusEncoding(generatedParticle.GetStatusCode(), 0).fullEncoding);
263+
generatedParticle.SetUniqueID(mGenID);
264+
generatedParticle.SetBit(ParticleStatus::kToBeDone, //
265+
o2::mcgenstatus::getHepMCStatusCode(generatedParticle.GetStatusCode()) == 1);
266+
return generatedParticle;
267+
}
268+
269+
void initGenMap()
270+
{
271+
genMap[0] = [this]()
272+
{ return generateParticle0(); };
273+
genMap[1] = [this]()
274+
{ return generateParticle1(); };
275+
}
276+
};
277+
278+
} // namespace eventgen
279+
} // namespace o2
280+
281+
// Performance test generator
282+
// fraction == -1 enables the fixed number of signal particles per event (nsig)
283+
// tag selects the generator type to be used
284+
FairGenerator *
285+
Generator_Performance(const float fraction = 0.03f, const unsigned short int nsig = 100, unsigned short int tag = 1)
286+
{
287+
auto generator = new o2::eventgen::GenPerf(fraction, nsig, tag);
288+
return generator;
289+
}

‎MC/config/common/ini/GeneratorPerformance.ini‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,4 @@
1+
# Test performance generator for multidimensional studies using fraction based signal particles per event compared to underlying event
2+
# generated with PYTHIA8 pp at 14 TeV
3+
[GeneratorHybrid]
4+
configFile = ${O2DPG_MC_CONFIG_ROOT}/MC/config/common/external/generator/perfConf.json

‎MC/config/common/ini/GeneratorPerformanceFix.ini‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,6 @@
1+
# Test performance generator for multidimensional studies using fix number of signal particles per event
2+
# Parameters are in order: fraction of signal particles, fixed number of signal particles per event, tag to select the generator type
3+
# fraction = -1 enables the fixed number of signal particles per event (nsig)
4+
[GeneratorExternal]
5+
fileName = ${O2DPG_MC_CONFIG_ROOT}/MC/config/common/external/generator/performanceGenerator.C
6+
funcName = Generator_Performance(-1, 100, 1)

‎MC/config/common/ini/tests/GeneratorPerformance.C‎

Lines changed: 45 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,45 @@
1+
int External() {
2+
std::string path{"o2sim_Kine.root"};
3+
TFile file(path.c_str(), "READ");
4+
if (file.IsZombie()) {
5+
std::cerr << "Cannot open ROOT file " << path << "\n";
6+
return 1;
7+
}
8+
auto tree = (TTree *)file.Get("o2sim");
9+
if (!tree) {
10+
std::cerr << "Cannot find tree 'o2sim' in file " << path << "\n";
11+
return 1;
12+
}
13+
// Get the MCTrack branch
14+
std::vector<o2::MCTrack> *tracks{};
15+
tree->SetBranchAddress("MCTrack", &tracks);
16+
// Check if processes with ID 42 are available
17+
// And are 0.03 of the UE tracks
18+
const int processID = 42; // Performance test particle custom process ID
19+
int nEvents = tree->GetEntries();
20+
short int count_perf = 0;
21+
22+
for (int i = 0; i < nEvents; i++) {
23+
tree->GetEntry(i);
24+
int nTracks = tracks->size();
25+
count_perf = 0;
26+
for (auto &track : *tracks)
27+
{
28+
const auto& process = track.getProcess();
29+
if (process == 42)
30+
{
31+
count_perf++;
32+
}
33+
}
34+
int UEtracks = nTracks - count_perf;
35+
unsigned short int expSig = std::lround(0.03 * UEtracks);
36+
if (count_perf != expSig)
37+
{
38+
std::cerr << "Event " << i << ": Expected " << expSig << " performance test particles, found " << count_perf << "\n";
39+
return 1;
40+
}
41+
}
42+
43+
file.Close();
44+
return 0;
45+
}

0 commit comments

Comments
 (0)