Skip to content

Commit 05fde30

Browse files
fgrosafgrosa
authored
[PWGHF] Add possibility to replace particles in the stack for HF generator (#1985)
* Add possibility to replace particles in the stack for HF generator * Minor change * Add Xic0 excited states * Fix test * Add protection * Fix file name * Minor fix * Minor fix * Make test more robust against small amount of test events
1 parent 7279641 commit 05fde30

File tree

4 files changed

+466
-26
lines changed

4 files changed

+466
-26
lines changed

‎MC/config/PWGHF/external/generator/generator_pythia8_gaptriggered_hf.C‎

Lines changed: 136 additions & 26 deletions
Original file line numberDiff line numberDiff line change
@@ -2,6 +2,7 @@
22
#include "Generators/GeneratorPythia8.h"
33
#include "Pythia8/Pythia.h"
44
#include "TRandom.h"
5+
#include "TDatabasePDG.h"
56
#include <fairlogger/Logger.h>
67

78
#include <string>
@@ -16,21 +17,39 @@ public:
1617
GeneratorPythia8GapTriggeredHF() = default;
1718

1819
/// constructor
19-
GeneratorPythia8GapTriggeredHF(int inputTriggerRatio = 5, std::vector<int> quarkPdgList = {}, std::vector<int> hadronPdgList = {})
20+
GeneratorPythia8GapTriggeredHF(int inputTriggerRatio = 5, std::vector<int> quarkPdgList = {}, std::vector<int> hadronPdgList = {}, std::vector<std::array<int, 2>> partPdgToReplaceList = {}, std::vector<float> freqReplaceList = {})
2021
{
2122

2223
mGeneratedEvents = 0;
2324
mInverseTriggerRatio = inputTriggerRatio;
2425
mQuarkRapidityMin = -1.5;
25-
mQuarkRapidityMax = -1.5;
26+
mQuarkRapidityMax = 1.5;
2627
mHadRapidityMin = -1.5;
27-
mHadRapidityMax = -1.5;
28+
mHadRapidityMax = 1.5;
2829
mQuarkPdg = 0;
2930
mHadronPdg = 0;
3031
mQuarkPdgList = quarkPdgList;
3132
mHadronPdgList = hadronPdgList;
33+
mPartPdgToReplaceList = partPdgToReplaceList;
34+
mFreqReplaceList = freqReplaceList;
35+
// Ds1*(2700), Ds1*(2860), Ds3*(2860), Xic(3055)+, Xic(3080)+, Xic(3055)0, Xic(3080)0
36+
mCustomPartPdgs = {30433, 40433, 437, 4315, 4316, 4325, 4326};
37+
mCustomPartMasses[30433] = 2.714f;
38+
mCustomPartMasses[40433] = 2.859f;
39+
mCustomPartMasses[437] = 2.860f;
40+
mCustomPartMasses[4315] = 3.0590f;
41+
mCustomPartMasses[4316] = 3.0799f;
42+
mCustomPartMasses[4325] = 3.0559f;
43+
mCustomPartMasses[4326] = 3.0772f;
44+
mCustomPartWidths[30433] = 0.122f;
45+
mCustomPartWidths[40433] = 0.160f;
46+
mCustomPartWidths[437] = 0.053f;
47+
mCustomPartWidths[4315] = 0.0064f;
48+
mCustomPartWidths[4316] = 0.0056f;
49+
mCustomPartWidths[4325] = 0.0078f;
50+
mCustomPartWidths[4326] = 0.0036f;
3251
Print();
33-
}
52+
}
3453

3554
/// Destructor
3655
~GeneratorPythia8GapTriggeredHF() = default;
@@ -54,6 +73,11 @@ public:
5473
{
5574
LOG(info)<<Form("* %d ", pdg);
5675
}
76+
LOG(info)<<Form("* Replacements: ");
77+
for (auto iRepl{0u}; iRepl<mPartPdgToReplaceList.size(); ++iRepl)
78+
{
79+
LOGP(info, "* {} -> {} (freq. {})", mPartPdgToReplaceList[iRepl].at(0), mPartPdgToReplaceList[iRepl].at(1), mFreqReplaceList[iRepl]);
80+
}
5781
LOG(info)<<"***********************************************************************";
5882
}
5983

@@ -62,6 +86,21 @@ public:
6286
addSubGenerator(0, "Minimum bias");
6387
addSubGenerator(4, "Charm injected");
6488
addSubGenerator(5, "Beauty injected");
89+
90+
std::vector<int> pdgToReplace{};
91+
for (auto iRepl{0u}; iRepl<mPartPdgToReplaceList.size(); ++iRepl)
92+
{
93+
for (auto iPdgRep{0u}; iPdgRep<pdgToReplace.size(); ++iPdgRep) {
94+
if (mPartPdgToReplaceList[iRepl].at(0) == pdgToReplace[iPdgRep]) {
95+
mFreqReplaceList[iRepl] += mFreqReplaceList[iPdgRep];
96+
}
97+
}
98+
if (mFreqReplaceList[iRepl] > 1.f) {
99+
LOGP(fatal, "Replacing more than 100% of a particle!");
100+
}
101+
pdgToReplace.push_back(mPartPdgToReplaceList[iRepl].at(0));
102+
}
103+
65104
return o2::eventgen::GeneratorPythia8::Init();
66105
}
67106

@@ -115,7 +154,7 @@ protected:
115154
{
116155
if (GeneratorPythia8::generateEvent())
117156
{
118-
genOk = selectEvent(mPythia.event);
157+
genOk = selectEvent();
119158
}
120159
}
121160
notifySubGenerator(mQuarkPdg);
@@ -136,34 +175,35 @@ protected:
136175
return true;
137176
}
138177

139-
bool selectEvent(const Pythia8::Event &event)
178+
bool selectEvent()
140179
{
141180

142181
bool isGoodAtPartonLevel{mQuarkPdgList.size() == 0};
143182
bool isGoodAtHadronLevel{mHadronPdgList.size() == 0};
183+
bool anyPartToReplace{mPartPdgToReplaceList.size() > 0};
144184

145-
for (auto iPart{0}; iPart < event.size(); ++iPart)
185+
for (auto iPart{0}; iPart < mPythia.event.size(); ++iPart)
146186
{
147187

148188
// search for Q-Qbar mother with at least one Q in rapidity window
149189
if (!isGoodAtPartonLevel)
150190
{
151-
auto daughterList = event[iPart].daughterList();
191+
auto daughterList = mPythia.event[iPart].daughterList();
152192
bool hasQ = false, hasQbar = false, atSelectedY = false;
153193
for (auto iDau : daughterList)
154194
{
155-
if (event[iDau].id() == mQuarkPdg)
195+
if (mPythia.event[iDau].id() == mQuarkPdg)
156196
{
157197
hasQ = true;
158-
if ((event[iDau].y() > mQuarkRapidityMin) && (event[iDau].y() < mQuarkRapidityMax))
198+
if ((mPythia.event[iDau].y() > mQuarkRapidityMin) && (mPythia.event[iDau].y() < mQuarkRapidityMax))
159199
{
160200
atSelectedY = true;
161201
}
162202
}
163-
if (event[iDau].id() == -mQuarkPdg)
203+
if (mPythia.event[iDau].id() == -mQuarkPdg)
164204
{
165205
hasQbar = true;
166-
if ((event[iDau].y() > mQuarkRapidityMin) && (event[iDau].y() < mQuarkRapidityMax))
206+
if ((mPythia.event[iDau].y() > mQuarkRapidityMin) && (mPythia.event[iDau].y() < mQuarkRapidityMax))
167207
{
168208
atSelectedY = true;
169209
}
@@ -178,26 +218,91 @@ protected:
178218
// search for hadron in rapidity window
179219
if (!isGoodAtHadronLevel)
180220
{
181-
int id = std::abs(event[iPart].id());
182-
float rap = event[iPart].y();
221+
int id = std::abs(mPythia.event[iPart].id());
222+
float rap = mPythia.event[iPart].y();
183223
if (id == mHadronPdg && rap > mHadRapidityMin && rap < mHadRapidityMax)
184224
{
185225
isGoodAtHadronLevel = true;
186226
}
187227
}
188228

189-
// we send the trigger
190-
if (isGoodAtPartonLevel && isGoodAtHadronLevel)
229+
// if requested, we replace the particle
230+
const double pseudoRndm = mPythia.event[iPart].pT() * 1000. - (int64_t)(mPythia.event[iPart].pT() * 1000);
231+
for (auto iPartToReplace{0u}; iPartToReplace<mPartPdgToReplaceList.size(); ++iPartToReplace) {
232+
if (std::abs(mPythia.event[iPart].id()) == mPartPdgToReplaceList[iPartToReplace][0] && pseudoRndm < mFreqReplaceList[iPartToReplace]) {
233+
LOGP(debug, "REPLACING PARTICLE {} WITH {}, BEING RNDM {}", mPartPdgToReplaceList[iPartToReplace][0], mPartPdgToReplaceList[iPartToReplace][1], pseudoRndm);
234+
replaceParticle(iPart, mPartPdgToReplaceList[iPartToReplace][1]);
235+
break;
236+
}
237+
}
238+
239+
// we send the trigger immediately (if there are no particles to replace, that can be different from the trigger ones)
240+
if (isGoodAtPartonLevel && isGoodAtHadronLevel && !anyPartToReplace)
191241
{
192-
LOG(debug)<<"EVENT SELECTED: Found particle "<<event[iPart].id()<<" at rapidity "<<event[iPart].y()<<"\n";
242+
LOG(debug)<<"EVENT SELECTED: Found particle "<<mPythia.event[iPart].id()<<" at rapidity "<<mPythia.event[iPart].y()<<"\n";
193243
return true;
194244
}
195245
}
196246

247+
// we send the trigger
248+
if (isGoodAtPartonLevel && isGoodAtHadronLevel) {
249+
return true;
250+
}
251+
197252
return false;
198253
};
199254

200-
private:
255+
bool replaceParticle(int iPartToReplace, int pdgCodeNew) {
256+
257+
std::array<int, 2> mothers = {mPythia.event[iPartToReplace].mother1(), mPythia.event[iPartToReplace].mother2()};
258+
259+
std::array<int, 25> pdgDiquarks = {1103, 2101, 2103, 2203, 3101, 3103, 3201, 3203, 3303, 4101, 4103, 4201, 4203, 4301, 4303, 4403, 5101, 5103, 5201, 5203, 5301, 5303, 5401, 5403, 5503};
260+
261+
for (auto const& mother: mothers) {
262+
auto pdgMother = std::abs(mPythia.event[mother].id());
263+
if (pdgMother > 100 && std::find(pdgDiquarks.begin(), pdgDiquarks.end(), pdgMother) == pdgDiquarks.end()) {
264+
return false;
265+
}
266+
}
267+
268+
int charge = mPythia.event[iPartToReplace].id() / std::abs(mPythia.event[iPartToReplace].id());
269+
float px = mPythia.event[iPartToReplace].px();
270+
float py = mPythia.event[iPartToReplace].py();
271+
float pz = mPythia.event[iPartToReplace].pz();
272+
float mass = 0.f;
273+
274+
if (std::find(mCustomPartPdgs.begin(), mCustomPartPdgs.end(), pdgCodeNew) == mCustomPartPdgs.end()) { // not a custom particle
275+
float width = TDatabasePDG::Instance()->GetParticle(pdgCodeNew)->Width();
276+
float massRest = TDatabasePDG::Instance()->GetParticle(pdgCodeNew)->Mass();
277+
if (width > 0.f) {
278+
mass = gRandom->BreitWigner(massRest, width);
279+
} else {
280+
mass = massRest;
281+
}
282+
} else {
283+
if (mCustomPartWidths[pdgCodeNew] > 0.f) {
284+
mass = gRandom->BreitWigner(mCustomPartMasses[pdgCodeNew], mCustomPartWidths[pdgCodeNew]);
285+
} else {
286+
mass = mCustomPartMasses[pdgCodeNew];
287+
}
288+
}
289+
float energy = std::sqrt(px*px + py*py + pz*pz + mass*mass);
290+
291+
// we remove particle to replace and its daughters
292+
mPythia.event[iPartToReplace].undoDecay();
293+
mPythia.event.remove(iPartToReplace, iPartToReplace, true); // we remove the original particle
294+
295+
int status = std::abs(mPythia.event[iPartToReplace].status());
296+
if (status < 81 || status > 89) {
297+
status = 81;
298+
}
299+
mPythia.event.append(charge * pdgCodeNew, status, mothers[0], mothers[1], 0, 0, 0, 0, px, py, pz, energy, mass);
300+
mPythia.moreDecays(); // we need to decay the new particle
301+
302+
return true;
303+
}
304+
305+
private:
201306
// Interface to override import particles
202307
Pythia8::Event mOutputEvent;
203308

@@ -209,6 +314,11 @@ private:
209314
float mHadRapidityMin;
210315
float mHadRapidityMax;
211316
unsigned int mUsedSeed;
317+
std::vector<std::array<int, 2>> mPartPdgToReplaceList;
318+
std::vector<float> mFreqReplaceList;
319+
std::vector<int> mCustomPartPdgs;
320+
std::map<int, float> mCustomPartMasses;
321+
std::map<int, float> mCustomPartWidths;
212322

213323
// Control gap-triggering
214324
unsigned long long mGeneratedEvents;
@@ -223,9 +333,9 @@ private:
223333

224334
// Predefined generators:
225335
// Charm-enriched
226-
FairGenerator *GeneratorPythia8GapTriggeredCharm(int inputTriggerRatio, float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector<int> hadronPdgList = {})
336+
FairGenerator *GeneratorPythia8GapTriggeredCharm(int inputTriggerRatio, float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector<int> hadronPdgList = {}, std::vector<std::array<int, 2>> partPdgToReplaceList = {}, std::vector<float> freqReplaceList = {})
227337
{
228-
auto myGen = new GeneratorPythia8GapTriggeredHF(inputTriggerRatio, std::vector<int>{4}, hadronPdgList);
338+
auto myGen = new GeneratorPythia8GapTriggeredHF(inputTriggerRatio, std::vector<int>{4}, hadronPdgList, partPdgToReplaceList, freqReplaceList);
229339
auto seed = (gRandom->TRandom::GetSeed() % 900000000);
230340
myGen->setUsedSeed(seed);
231341
myGen->readString("Random:setSeed on");
@@ -239,9 +349,9 @@ FairGenerator *GeneratorPythia8GapTriggeredCharm(int inputTriggerRatio, float yQ
239349
}
240350

241351
// Beauty-enriched
242-
FairGenerator *GeneratorPythia8GapTriggeredBeauty(int inputTriggerRatio, float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector<int> hadronPdgList = {})
352+
FairGenerator *GeneratorPythia8GapTriggeredBeauty(int inputTriggerRatio, float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector<int> hadronPdgList = {}, std::vector<std::array<int, 2>> partPdgToReplaceList = {}, std::vector<float> freqReplaceList = {})
243353
{
244-
auto myGen = new GeneratorPythia8GapTriggeredHF(inputTriggerRatio, std::vector<int>{5}, hadronPdgList);
354+
auto myGen = new GeneratorPythia8GapTriggeredHF(inputTriggerRatio, std::vector<int>{5}, hadronPdgList, partPdgToReplaceList, freqReplaceList);
245355
auto seed = (gRandom->TRandom::GetSeed() % 900000000);
246356
myGen->setUsedSeed(seed);
247357
myGen->readString("Random:setSeed on");
@@ -255,9 +365,9 @@ FairGenerator *GeneratorPythia8GapTriggeredBeauty(int inputTriggerRatio, float y
255365
}
256366

257367
// Charm and beauty enriched (with same ratio)
258-
FairGenerator *GeneratorPythia8GapTriggeredCharmAndBeauty(int inputTriggerRatio, float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector<int> hadronPdgList = {})
368+
FairGenerator *GeneratorPythia8GapTriggeredCharmAndBeauty(int inputTriggerRatio, float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector<int> hadronPdgList = {}, std::vector<std::array<int, 2>> partPdgToReplaceList = {}, std::vector<float> freqReplaceList = {})
259369
{
260-
auto myGen = new GeneratorPythia8GapTriggeredHF(inputTriggerRatio, std::vector<int>{4, 5}, hadronPdgList);
370+
auto myGen = new GeneratorPythia8GapTriggeredHF(inputTriggerRatio, std::vector<int>{4, 5}, hadronPdgList, partPdgToReplaceList, freqReplaceList);
261371
auto seed = (gRandom->TRandom::GetSeed() % 900000000);
262372
myGen->setUsedSeed(seed);
263373
myGen->readString("Random:setSeed on");
@@ -270,13 +380,13 @@ FairGenerator *GeneratorPythia8GapTriggeredCharmAndBeauty(int inputTriggerRatio,
270380
return myGen;
271381
}
272382

273-
FairGenerator *GeneratorPythia8GapHF(int inputTriggerRatio, float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector<int> quarkPdgList = {}, std::vector<int> hadronPdgList = {})
383+
FairGenerator *GeneratorPythia8GapHF(int inputTriggerRatio, float yQuarkMin = -1.5, float yQuarkMax = 1.5, float yHadronMin = -1.5, float yHadronMax = 1.5, std::vector<int> quarkPdgList = {}, std::vector<int> hadronPdgList = {}, std::vector<std::array<int, 2>> partPdgToReplaceList = {}, std::vector<float> freqReplaceList = {})
274384
{
275385
if (hadronPdgList.size() == 0 && quarkPdgList.size() == 0)
276386
{
277387
LOG(fatal) << "GeneratorPythia8GapHF: At least one quark or hadron PDG code must be specified";
278388
}
279-
auto myGen = new GeneratorPythia8GapTriggeredHF(inputTriggerRatio, quarkPdgList, hadronPdgList);
389+
auto myGen = new GeneratorPythia8GapTriggeredHF(inputTriggerRatio, quarkPdgList, hadronPdgList, partPdgToReplaceList, freqReplaceList);
280390
auto seed = (gRandom->TRandom::GetSeed() % 900000000);
281391
myGen->setUsedSeed(seed);
282392
myGen->readString("Random:setSeed on");

‎MC/config/PWGHF/ini/GeneratorHF_D2H_ccbar_and_bbbar_gap5_DResoTrigger.ini‎

Lines changed: 8 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,8 @@
1+
### The external generator derives from GeneratorPythia8.
2+
[GeneratorExternal]
3+
fileName=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/external/generator/generator_pythia8_gaptriggered_hf.C
4+
funcName=GeneratorPythia8GapTriggeredCharmAndBeauty(5, -1.5, 1.5, -1.5, 1.5, {10433, 435}, {{10433, 30433}, {10433, 437}, {435, 4325}, {435, 4326}, {425, 4315}, {425, 4316}}, {0.1, 0.1, 0.1, 0.1, 0.5, 0.5})
5+
6+
[GeneratorPythia8]
7+
config=${O2DPG_MC_CONFIG_ROOT}/MC/config/PWGHF/pythia8/generator/pythia8_charmhadronic_with_decays_DResoTrigger.cfg
8+
includePartonEvent=true

0 commit comments

Comments
 (0)