Skip to content

Commit 765080f

Browse files
jackal1-66jackal1-66
committed
Implemented two modes for Z0
1 parent 6248725 commit 765080f

File tree

1 file changed

+99
-65
lines changed

1 file changed

+99
-65
lines changed

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

Lines changed: 99 additions & 65 deletions
Original file line numberDiff line numberDiff line change
@@ -8,7 +8,8 @@ namespace o2
88
class GenPerf : public Generator
99
{
1010
public:
11-
GenPerf(float fraction = 0.03f, unsigned short int nsig = 100, unsigned short int tag = 1)
11+
GenPerf(float fraction = 0.03f, unsigned short int nsig = 100, unsigned short int tag = 1, bool setDecayZ0 = false)
12+
: mDecayZ0(setDecayZ0)
1213
{
1314
if (fraction == -1) {
1415
LOG(info) << nsig << " Signal particles will be generated in each event";
@@ -35,8 +36,15 @@ namespace o2
3536
mTag = tag;
3637
LOG(info) << "Generator with tag " << mTag << " is selected";
3738
}
38-
mDecayer = std::make_unique<DecayerPythia8>();
39-
mDecayer->Init();
39+
if (mDecayZ0){
40+
LOG(info) << "Z0 decays will be handled with Pythia8";
41+
mPythia = std::make_unique<Pythia8::Pythia>();
42+
// Configure Pythia8 with minimal beam setup for standalone decays
43+
mPythia->readString("WeakSingleBoson:ffbar2gmZ = on");
44+
mPythia->init(); // Initialize
45+
} else {
46+
LOG(info) << "Z0 decays will be created but not transported (incompatible with Geant4 physics list)";
47+
}
4048
Generator::setTimeUnit(1.0);
4149
Generator::setPositionUnit(1.0);
4250
}
@@ -64,7 +72,20 @@ namespace o2
6472
unsigned short nSig = (mFraction == -1) ? mNSig : std::lround(mFraction * mNUE);
6573
LOG(debug) << "Generating additional " << nSig << " particles";
6674
for (int k = 0; k < nSig; k++){
67-
mParticles.push_back(genMap[mTag]());
75+
auto part = genMap[mTag]();
76+
if(part.GetPdgCode() == 23) {
77+
if(!mDecayZ0) {
78+
mParticles.push_back(part);
79+
} else {
80+
auto daughters = decayZ0(part);
81+
for (auto &dau : daughters)
82+
{
83+
mParticles.push_back(dau);
84+
}
85+
}
86+
} else {
87+
mParticles.push_back(part);
88+
}
6889
}
6990
return kTRUE;
7091
}
@@ -74,7 +95,8 @@ namespace o2
7495
unsigned short int mNSig = 0; // Number of particles to generate
7596
unsigned int mNUE = 0; // Number of tracks in the Underlying event
7697
unsigned short int mTag = 1; // Tag to select the generation function
77-
std::unique_ptr<DecayerPythia8> mDecayer; // Pythia8 decayer for particles not present in the physics list of Geant4 (like Z0)
98+
bool mDecayZ0 = false; // Whether to decay Z0 with Pythia8
99+
std::unique_ptr<Pythia8::Pythia> mPythia; // Pythia8 instance for particle decays not present in the physics list of Geant4 (like Z0)
78100
const std::vector<std::shared_ptr<o2::eventgen::Generator>>* mGenList = nullptr; // Cached generators list
79101
std::map<unsigned short int, std::function<TParticle()>> genMap;
80102
UInt_t mGenID = 42;
@@ -148,36 +170,6 @@ namespace o2
148170
generatedParticle.SetUniqueID(mGenID);
149171
if (pdgCode == 23) {
150172
generatedParticle.SetBit(ParticleStatus::kToBeDone, false); // Force Z0 to be decayed by the transport
151-
LOG(debug) << "Processing Z0 with DecayerPythia8";
152-
TLorentzVector *pDec = new TLorentzVector(px, py, pz, energy);
153-
mDecayer->Decay(pdgCode, pDec);
154-
TClonesArray *daughters = new TClonesArray("TParticle");
155-
mDecayer->ImportParticles(daughters);
156-
unsigned short int nDaughters = daughters->GetEntriesFast();
157-
if (daughters && nDaughters > 0) {
158-
for (int i = 0; i < daughters->GetEntriesFast(); i++) {
159-
TParticle* daughter = (TParticle*)daughters->At(i);
160-
daughter->SetUniqueID(mGenID);
161-
if (i > 0)
162-
{
163-
daughter->SetBit(ParticleStatus::kToBeDone, //
164-
o2::mcgenstatus::getHepMCStatusCode(generatedParticle.GetStatusCode()) == 1);
165-
}
166-
else
167-
{
168-
// First daughter is the mother (Z0)
169-
daughter->SetBit(ParticleStatus::kToBeDone, false);
170-
}
171-
LOG(debug) << "Daughter " << i << ": PDG=" << daughter->GetPdgCode() << ", E=" << daughter->Energy() << ", p=(" << daughter->Px() << "," << daughter->Py() << "," << daughter->Pz() << ")";
172-
mParticles.push_back(*daughter);
173-
}
174-
LOG(debug) << "Z0 decayed into " << daughters->GetEntriesFast() << " particles";
175-
daughters->Clear("C");
176-
delete daughters;
177-
} else {
178-
LOG(warn) << "DecayerPythia8 failed to decay Z0 or no daughters found";
179-
}
180-
delete pDec;
181173
} else {
182174
generatedParticle.SetBit(ParticleStatus::kToBeDone, //
183175
o2::mcgenstatus::getHepMCStatusCode(generatedParticle.GetStatusCode()) == 1);
@@ -305,34 +297,8 @@ namespace o2
305297
generatedParticle.SetStatusCode(o2::mcgenstatus::MCGenStatusEncoding(generatedParticle.GetStatusCode(), 0).fullEncoding);
306298
generatedParticle.SetUniqueID(mGenID);
307299
if (pdgCode == 23) {
308-
generatedParticle.SetBit(ParticleStatus::kToBeDone, false); // Force Z0 to be decayed by the transport
309-
LOG(debug) << "Processing Z0 with DecayerPythia8";
310-
TLorentzVector *pDec = new TLorentzVector(px, py, pz, energy);
311-
mDecayer->Decay(pdgCode, pDec);
312-
TClonesArray *daughters = new TClonesArray("TParticle");
313-
mDecayer->ImportParticles(daughters);
314-
unsigned short int nDaughters = daughters->GetEntriesFast();
315-
if (daughters && nDaughters > 0) {
316-
for (int i = 0; i < daughters->GetEntriesFast(); i++) {
317-
TParticle* daughter = (TParticle*)daughters->At(i);
318-
daughter->SetUniqueID(mGenID);
319-
if (i > 0) {
320-
daughter->SetBit(ParticleStatus::kToBeDone, //
321-
o2::mcgenstatus::getHepMCStatusCode(generatedParticle.GetStatusCode()) == 1);
322-
} else {
323-
// First daughter is the mother (Z0)
324-
daughter->SetBit(ParticleStatus::kToBeDone, false);
325-
}
326-
LOG(debug) << "Daughter " << i << ": PDG=" << daughter->GetPdgCode() << ", E=" << daughter->Energy() << ", p=(" << daughter->Px() << "," << daughter->Py() << "," << daughter->Pz() << ")";
327-
mParticles.push_back(*daughter);
328-
}
329-
LOG(debug) << "Z0 decayed into " << daughters->GetEntriesFast() << " particles";
330-
daughters->Clear("C");
331-
delete daughters;
332-
} else {
333-
LOG(warn) << "DecayerPythia8 failed to decay Z0 or no daughters found";
334-
}
335-
delete pDec;
300+
generatedParticle.SetBit(ParticleStatus::kToBeDone, false);
301+
// Z0 will follow another decay procedure
336302
} else {
337303
generatedParticle.SetBit(ParticleStatus::kToBeDone, //
338304
o2::mcgenstatus::getHepMCStatusCode(generatedParticle.GetStatusCode()) == 1);
@@ -347,6 +313,74 @@ namespace o2
347313
genMap[1] = [this]()
348314
{ return generateParticle1(); };
349315
}
316+
317+
std::vector<TParticle> decayZ0(TParticle &z0)
318+
{
319+
std::vector<TParticle> subparts;
320+
// Start a Pythia8 event with Z0
321+
mPythia->next();
322+
// Find the Z0 and its final products
323+
auto &event = mPythia->event;
324+
for (int j = 0; j < event.size(); ++j)
325+
{
326+
const Pythia8::Particle &p = event[j];
327+
if (p.id() == 23) // PDG code for Z0
328+
{
329+
// Push Z0 itself
330+
subparts.push_back(TParticle(p.id(), p.status(),
331+
-1, -1, -1, -1,
332+
p.px(), p.py(),
333+
p.pz(), p.e(),
334+
z0.Vx(), z0.Vy(), z0.Vz(), 0.0));
335+
subparts.back().SetStatusCode(o2::mcgenstatus::MCGenStatusEncoding(p.status(), 0).fullEncoding);
336+
subparts.back().SetUniqueID(mGenID);
337+
subparts.back().SetBit(ParticleStatus::kToBeDone, false);
338+
// Navigate through intermediate Z0s to find final decay products
339+
int iZ0 = j;
340+
while (event[iZ0].daughter1() != 0 &&
341+
event[event[iZ0].daughter1()].id() == 23)
342+
{
343+
iZ0 = event[iZ0].daughter1();
344+
}
345+
// Recursively collect all final-state descendants
346+
std::function<void(int)> collectFinalState = [&](int idx)
347+
{
348+
const Pythia8::Particle &particle = event[idx];
349+
350+
if (particle.isFinal())
351+
{
352+
subparts.push_back(TParticle(particle.id(), particle.status(),
353+
-1, -1, -1, -1,
354+
particle.px(), particle.py(),
355+
particle.pz(), particle.e(),
356+
p.xProd() + z0.Vx(), p.yProd() + z0.Vy(), p.zProd() + z0.Vz(), 0.0));
357+
subparts.back().SetStatusCode(o2::mcgenstatus::MCGenStatusEncoding(particle.status(), 0).fullEncoding);
358+
subparts.back().SetUniqueID(mGenID+1);
359+
subparts.back().SetBit(ParticleStatus::kToBeDone,
360+
o2::mcgenstatus::getHepMCStatusCode(subparts.back().GetStatusCode()) == 1);
361+
}
362+
else
363+
{
364+
// Not final-state, recurse through daughters
365+
int d1 = particle.daughter1();
366+
int d2 = particle.daughter2();
367+
if (d1 > 0)
368+
{
369+
for (int k = d1; k <= d2; ++k)
370+
{
371+
collectFinalState(k);
372+
}
373+
}
374+
}
375+
};
376+
377+
// Start collecting from the final Z0
378+
collectFinalState(iZ0);
379+
break; // Found and processed the Z0
380+
}
381+
}
382+
return subparts;
383+
}
350384
};
351385

352386
} // namespace eventgen
@@ -356,8 +390,8 @@ namespace o2
356390
// fraction == -1 enables the fixed number of signal particles per event (nsig)
357391
// tag selects the generator type to be used
358392
FairGenerator *
359-
Generator_Performance(const float fraction = 0.03f, const unsigned short int nsig = 100, unsigned short int tag = 1)
393+
Generator_Performance(const float fraction = 0.03f, const unsigned short int nsig = 100, unsigned short int tag = 1, bool setDecayZ0 = false)
360394
{
361-
auto generator = new o2::eventgen::GenPerf(fraction, nsig, tag);
395+
auto generator = new o2::eventgen::GenPerf(fraction, nsig, tag, setDecayZ0);
362396
return generator;
363397
}

0 commit comments

Comments
 (0)