Skip to content

Commit 2bd4888

Browse files
DanielSamitzDanielSamitz
authored
PWGEM: hf cocktail task updates (#6235)
* check partonic history for pairs * fixed searchMoterhs, more info in hfTable * clang-format * fix partonic check for beauty * fix partonic checks * MegaLInter * clang-format
1 parent f516797 commit 2bd4888

2 files changed

Lines changed: 249 additions & 55 deletions

File tree

‎PWGEM/Dilepton/Tasks/lmeeHFCocktail.cxx‎

Lines changed: 150 additions & 55 deletions
Original file line numberDiff line numberDiff line change
@@ -43,9 +43,21 @@ namespace o2::aod
4343
namespace hftable
4444
{
4545
DECLARE_SOA_COLUMN(IsHF, isHF, int);
46-
}
46+
DECLARE_SOA_COLUMN(BHadronId, bHadronId, int);
47+
DECLARE_SOA_COLUMN(CHadronId, cHadronId, int);
48+
DECLARE_SOA_COLUMN(BQuarkId, bQuarkId, int);
49+
DECLARE_SOA_COLUMN(CQuarkId, cQuarkId, int);
50+
DECLARE_SOA_COLUMN(BQuarkOriginId, bQuarkOriginId, int);
51+
DECLARE_SOA_COLUMN(CQuarkOriginId, cQuarkOriginId, int);
52+
} // namespace hftable
4753
DECLARE_SOA_TABLE(HfTable, "AOD", "HFTABLE",
48-
hftable::IsHF);
54+
hftable::IsHF,
55+
hftable::BHadronId,
56+
hftable::CHadronId,
57+
hftable::BQuarkId,
58+
hftable::CQuarkId,
59+
hftable::BQuarkOriginId,
60+
hftable::CQuarkOriginId);
4961
} // namespace o2::aod
5062

5163
const char* stageNames[3] = {"gen", "eff", "eff_and_acc"};
@@ -96,27 +108,78 @@ void doPair(T& p1, T& p2, std::vector<std::shared_ptr<TH1>> hMee, std::vector<st
96108
}
97109
}
98110

111+
struct MyConfigs : ConfigurableGroup {
112+
Configurable<float> fConfigPtMin{"cfgPtMin", 0.2, "min. pT"};
113+
Configurable<float> fConfigEtaMax{"cfgEtaMax", 0.8, "max. |eta|"};
114+
ConfigurableAxis fConfigPtBins{"cfgPtBins", {200, 0.f, 10.f}, "pT binning"};
115+
ConfigurableAxis fConfigEtaBins{"cfgEtaBins", {200, -10.f, 10.f}, "eta binning"};
116+
ConfigurableAxis fConfigMeeBins{"cfgMeeBins", {800, 0.f, 8.f}, "Mee binning"};
117+
ConfigurableAxis fConfigPteeBins{"cfgPteeBins", {400, 0.f, 10.f}, "pTee binning"};
118+
Configurable<bool> fConfigCheckPartonic{"cfgCheckPartonic", true, "check entire partonic history for pairs"};
119+
};
120+
99121
struct lmeehfcocktailprefilter {
100122

101123
Produces<o2::aod::HfTable> hfTable;
102124
void process(aod::McParticles const& mcParticles)
103125
{
104126
for (auto const& p : mcParticles) {
105127

106-
if (abs(p.pdgCode()) != 11 || o2::mcgenstatus::getHepMCStatusCode(p.statusCode()) != 1) {
107-
hfTable(EFromHFType::kNoE);
128+
if (abs(p.pdgCode()) != 11 || o2::mcgenstatus::getHepMCStatusCode(p.statusCode()) != 1 || !p.has_mothers()) {
129+
hfTable(EFromHFType::kNoE, -1, -1, -1, -1, -1, -1);
108130
continue;
109131
}
110132

111-
// no HF = 0; c->e = 1; b->e = 2; b->c->e = 3;
112-
int isHF = 0;
113-
if (IsFromCharm(p, mcParticles) > -1) {
114-
isHF = isHF + 1;
133+
int mother_pdg = mcParticles.iteratorAt(p.mothersIds()[0]).pdgCode();
134+
bool direct_charm_mother = abs(mother_pdg) < 1e+9 && (std::to_string(mother_pdg)[std::to_string(mother_pdg).length() - 3] == '4' || std::to_string(mother_pdg)[std::to_string(mother_pdg).length() - 4] == '4');
135+
if (abs(mother_pdg) == 443) {
136+
direct_charm_mother = false; // we don't want JPsi here
115137
}
116-
if (IsFromBeauty(p, mcParticles) > -1) {
117-
isHF = isHF + 2;
138+
int cHadronId = -1;
139+
if (direct_charm_mother) {
140+
cHadronId = p.mothersIds()[0];
141+
}
142+
bool direct_beauty_mother = abs(mother_pdg) < 1e+9 && (std::to_string(mother_pdg)[std::to_string(mother_pdg).length() - 3] == '5' || std::to_string(mother_pdg)[std::to_string(mother_pdg).length() - 4] == '5');
143+
int bHadronId = IsFromBeauty(p, mcParticles);
144+
145+
int bQuarkId = -1;
146+
int cQuarkId = -1;
147+
int cQuarkOriginId = -1;
148+
int bQuarkOriginId = -1;
149+
150+
int isHF = EFromHFType::kNoHFE;
151+
if (direct_charm_mother && bHadronId < 0) { // c->e
152+
isHF = EFromHFType::kCE;
153+
auto cHadron = mcParticles.iteratorAt(cHadronId);
154+
cQuarkId = searchMothers(cHadron, mcParticles, 4, true);
155+
if (cQuarkId > -1) {
156+
auto cQuark = mcParticles.iteratorAt(cQuarkId);
157+
cQuarkOriginId = searchMothers(cQuark, mcParticles, 4, false);
158+
}
159+
} else if (direct_charm_mother && bHadronId > -1) { // b->c->e
160+
isHF = EFromHFType::kBCE;
161+
auto bHadron = mcParticles.iteratorAt(bHadronId);
162+
bQuarkId = searchMothers(bHadron, mcParticles, 5, true);
163+
if (bQuarkId > -1) {
164+
auto bQuark = mcParticles.iteratorAt(bQuarkId);
165+
bQuarkOriginId = searchMothers(bQuark, mcParticles, 5, false);
166+
}
167+
auto cHadron = mcParticles.iteratorAt(cHadronId);
168+
cQuarkId = searchMothers(cHadron, mcParticles, 4, true);
169+
if (cQuarkId > -1) {
170+
auto cQuark = mcParticles.iteratorAt(cQuarkId);
171+
cQuarkOriginId = searchMothers(cQuark, mcParticles, 4, false);
172+
}
173+
} else if (direct_beauty_mother) { // b->e
174+
isHF = EFromHFType::kBE;
175+
auto bHadron = mcParticles.iteratorAt(bHadronId);
176+
bQuarkId = searchMothers(bHadron, mcParticles, 5, true);
177+
if (bQuarkId > -1) {
178+
auto bQuark = mcParticles.iteratorAt(bQuarkId);
179+
bQuarkOriginId = searchMothers(bQuark, mcParticles, 5, false);
180+
}
118181
}
119-
hfTable(isHF);
182+
hfTable(isHF, bHadronId, cHadronId, bQuarkId, cQuarkId, bQuarkOriginId, cQuarkOriginId);
120183
}
121184
}
122185
};
@@ -130,12 +193,7 @@ struct lmeehfcocktailbeauty {
130193
std::vector<std::shared_ptr<TH1>> hLSpp_Mee, hLSmm_Mee;
131194
std::vector<std::shared_ptr<TH2>> hLSpp_MeePtee, hLSmm_MeePtee;
132195

133-
Configurable<float> fConfigPtMin{"cfgPtMin", 0.2, "min. pT"};
134-
Configurable<float> fConfigEtaMax{"cfgEtaMax", 0.8, "max. |eta|"};
135-
ConfigurableAxis fConfigPtBins{"cfgPtBins", {200, 0.f, 10.f}, "pT binning"};
136-
ConfigurableAxis fConfigEtaBins{"cfgEtaBins", {200, -10.f, 10.f}, "eta binning"};
137-
ConfigurableAxis fConfigMeeBins{"cfgMeeBins", {800, 0.f, 8.f}, "Mee binning"};
138-
ConfigurableAxis fConfigPteeBins{"cfgPteeBins", {400, 0.f, 10.f}, "pTee binning"};
196+
MyConfigs myConfigs;
139197

140198
Filter hfFilter = o2::aod::hftable::isHF == static_cast<int>(EFromHFType::kBE) || o2::aod::hftable::isHF == static_cast<int>(EFromHFType::kBCE);
141199
using MyFilteredMcParticlesSmeared = soa::Filtered<soa::Join<aod::McParticles, aod::SmearedTracks, aod::HfTable>>;
@@ -154,10 +212,10 @@ struct lmeehfcocktailbeauty {
154212
const char* typeNamesSingle[2] = {"be", "bce"};
155213
const char* typeTitlesSingle[2] = {"b->e", "b->c->e"};
156214

157-
AxisSpec eta_axis = {fConfigEtaBins, "#eta"};
158-
AxisSpec pt_axis = {fConfigPtBins, "#it{p}_{T} (GeV/c)"};
159-
AxisSpec mass_axis = {fConfigMeeBins, "m_{ee} (GeV/c^{2})"};
160-
AxisSpec ptee_axis = {fConfigPteeBins, "#it{p}_{T,ee} (GeV/c)"};
215+
AxisSpec eta_axis = {myConfigs.fConfigEtaBins, "#eta"};
216+
AxisSpec pt_axis = {myConfigs.fConfigPtBins, "#it{p}_{T} (GeV/c)"};
217+
AxisSpec mass_axis = {myConfigs.fConfigMeeBins, "m_{ee} (GeV/c^{2})"};
218+
AxisSpec ptee_axis = {myConfigs.fConfigPteeBins, "#it{p}_{T,ee} (GeV/c)"};
161219

162220
// single histograms
163221
for (int i = 0; i < 2; i++) {
@@ -196,46 +254,65 @@ struct lmeehfcocktailbeauty {
196254

197255
void processBeauty(aod::McCollisions const& collisions, MyFilteredMcParticlesSmeared const& mcParticles, aod::McParticles const& mcParticlesAll)
198256
{
257+
for (auto const& p : mcParticles) {
258+
if (myConfigs.fConfigCheckPartonic && p.bQuarkOriginId() < 0) {
259+
continue;
260+
}
261+
int from_quark = p.isHF() - 2;
262+
doSingle(p, hEta[from_quark], hPt[from_quark], hPtEta[from_quark], myConfigs.fConfigPtMin, myConfigs.fConfigEtaMax);
263+
}
264+
199265
for (auto const& collision : collisions) {
200266

201267
registry.fill(HIST("NEvents"), 0.5);
202268

203-
for (auto const& p : mcParticles) {
204-
int from_quark = p.isHF() - 2;
205-
doSingle(p, hEta[from_quark], hPt[from_quark], hPtEta[from_quark], fConfigPtMin, fConfigEtaMax);
206-
}
207-
208269
auto const electronsGrouped = Electrons->sliceBy(perCollision, collision.globalIndex());
209270
auto const positronsGrouped = Positrons->sliceBy(perCollision, collision.globalIndex());
210271
// ULS spectrum
211272
for (auto const& [particle1, particle2] : combinations(o2::soa::CombinationsFullIndexPolicy(electronsGrouped, positronsGrouped))) {
273+
if (myConfigs.fConfigCheckPartonic) {
274+
if (particle1.bQuarkOriginId() < 0 || particle2.bQuarkOriginId() < 0 || particle1.bQuarkOriginId() != particle2.bQuarkOriginId())
275+
continue;
276+
}
212277
int type = IsHF(particle1, particle2, mcParticlesAll);
213278
if (type == static_cast<int>(EM_HFeeType::kUndef))
214279
continue;
215280
if ((type < static_cast<int>(EM_HFeeType::kBe_Be)) || (type > static_cast<int>(EM_HFeeType::kBCe_Be_SameB))) {
216281
LOG(error) << "Something is wrong here. There should only be pairs of type kBe_Be = 1, kBCe_BCe = 2 and kBCe_Be_SameB = 3 left at this point.";
217282
}
218-
doPair(particle1, particle2, hULS_Mee[type - 1], hULS_MeePtee[type - 1], fConfigPtMin, fConfigEtaMax);
219-
doPair(particle1, particle2, hULS_Mee[3], hULS_MeePtee[3], fConfigPtMin, fConfigEtaMax); // fill the 'allB' histograms that holds the sum of the others
283+
doPair(particle1, particle2, hULS_Mee[type - 1], hULS_MeePtee[type - 1], myConfigs.fConfigPtMin, myConfigs.fConfigEtaMax);
284+
doPair(particle1, particle2, hULS_Mee[3], hULS_MeePtee[3], myConfigs.fConfigPtMin, myConfigs.fConfigEtaMax); // fill the 'allB' histograms that holds the sum of the others
220285
}
221286
// LS spectrum
222287
for (auto const& [particle1, particle2] : combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(electronsGrouped, electronsGrouped))) {
288+
if (myConfigs.fConfigCheckPartonic) {
289+
if (particle1.bQuarkOriginId() < 0 || particle2.bQuarkOriginId() < 0 || particle1.bQuarkOriginId() != particle2.bQuarkOriginId())
290+
continue;
291+
}
223292
int type = IsHF(particle1, particle2, mcParticlesAll);
224293
if (type == static_cast<int>(EM_HFeeType::kUndef))
225294
continue;
226-
if (type != static_cast<int>(EM_HFeeType::kBCe_Be_DiffB)) {
227-
LOG(error) << "Something is wrong here. There should only be pairs of type kBCe_Be_DiffB = 4 left at this point.";
295+
if (myConfigs.fConfigCheckPartonic) {
296+
if (type != static_cast<int>(EM_HFeeType::kBCe_Be_DiffB)) {
297+
LOG(error) << "Something is wrong here. There should only be pairs of type kBCe_Be_DiffB = 4 left at this point.";
298+
}
228299
}
229-
doPair(particle1, particle2, hLSmm_Mee, hLSmm_MeePtee, fConfigPtMin, fConfigEtaMax);
300+
doPair(particle1, particle2, hLSmm_Mee, hLSmm_MeePtee, myConfigs.fConfigPtMin, myConfigs.fConfigEtaMax);
230301
}
231302
for (auto const& [particle1, particle2] : combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(positronsGrouped, positronsGrouped))) {
303+
if (myConfigs.fConfigCheckPartonic) {
304+
if (particle1.bQuarkOriginId() < 0 || particle2.bQuarkOriginId() < 0 || particle1.bQuarkOriginId() != particle2.bQuarkOriginId())
305+
continue;
306+
}
232307
int type = IsHF(particle1, particle2, mcParticlesAll);
233308
if (type == static_cast<int>(EM_HFeeType::kUndef))
234309
continue;
235-
if (type != static_cast<int>(EM_HFeeType::kBCe_Be_DiffB)) {
236-
LOG(error) << "Something is wrong here. There should only be pairs of type kBCe_Be_DiffB = 4 left at this point.";
310+
if (myConfigs.fConfigCheckPartonic) {
311+
if (type != static_cast<int>(EM_HFeeType::kBCe_Be_DiffB)) {
312+
LOG(error) << "Something is wrong here. There should only be pairs of type kBCe_Be_DiffB = 4 left at this point.";
313+
}
237314
}
238-
doPair(particle1, particle2, hLSpp_Mee, hLSpp_MeePtee, fConfigPtMin, fConfigEtaMax);
315+
doPair(particle1, particle2, hLSpp_Mee, hLSpp_MeePtee, myConfigs.fConfigPtMin, myConfigs.fConfigEtaMax);
239316
}
240317
}
241318
}
@@ -253,15 +330,10 @@ struct lmeehfcocktailcharm {
253330

254331
HistogramRegistry registry{"registry", {}};
255332

256-
std::vector<std::shared_ptr<TH1>> hEta, hPt, hULS_Mee;
257-
std::vector<std::shared_ptr<TH2>> hPtEta, hULS_MeePtee;
333+
std::vector<std::shared_ptr<TH1>> hEta, hPt, hULS_Mee, hLS_Mee;
334+
std::vector<std::shared_ptr<TH2>> hPtEta, hULS_MeePtee, hLS_MeePtee;
258335

259-
Configurable<float> fConfigPtMin{"cfgPtMin", 0.2f, "min. pT"};
260-
Configurable<float> fConfigEtaMax{"cfgEtaMax", 0.8f, "max. |eta|"};
261-
ConfigurableAxis fConfigPtBins{"cfgPtBins", {200, 0.f, 10.f}, "pT binning"};
262-
ConfigurableAxis fConfigEtaBins{"cfgEtaBins", {200, -10.f, 10.f}, "eta binning"};
263-
ConfigurableAxis fConfigMeeBins{"cfgMeeBins", {800, 0.f, 8.f}, "Mee binning"};
264-
ConfigurableAxis fConfigPteeBins{"cfgPteeBins", {400, 0.f, 10.f}, "pTee binning"};
336+
MyConfigs myConfigs;
265337

266338
Filter hfFilter = o2::aod::hftable::isHF == static_cast<int>(EFromHFType::kCE);
267339
using MyFilteredMcParticlesSmeared = soa::Filtered<soa::Join<aod::McParticles, aod::SmearedTracks, aod::HfTable>>;
@@ -279,10 +351,10 @@ struct lmeehfcocktailcharm {
279351
const char* typeNamesSingle = "ce";
280352
const char* typeTitlesSingle = "c->e";
281353

282-
AxisSpec eta_axis = {fConfigEtaBins, "#eta"};
283-
AxisSpec pt_axis = {fConfigPtBins, "#it{p}_{T} (GeV/c)"};
284-
AxisSpec mass_axis = {fConfigMeeBins, "m_{ee} (GeV/c^{2})"};
285-
AxisSpec ptee_axis = {fConfigPteeBins, "#it{p}_{T,ee} (GeV/c)"};
354+
AxisSpec eta_axis = {myConfigs.fConfigEtaBins, "#eta"};
355+
AxisSpec pt_axis = {myConfigs.fConfigPtBins, "#it{p}_{T} (GeV/c)"};
356+
AxisSpec mass_axis = {myConfigs.fConfigMeeBins, "m_{ee} (GeV/c^{2})"};
357+
AxisSpec ptee_axis = {myConfigs.fConfigPteeBins, "#it{p}_{T,ee} (GeV/c)"};
286358

287359
// single histograms
288360
for (int j = 0; j < 3; j++) {
@@ -297,29 +369,52 @@ struct lmeehfcocktailcharm {
297369
hULS_Mee.push_back(registry.add<TH1>(Form("ULS_Mee_%s_%s", typeNamesPairULS, stageNames[j]), Form("ULS Mee %s %s", typeNamesPairULS, stageNames[j]), HistType::kTH1F, {mass_axis}, true));
298370
hULS_MeePtee.push_back(registry.add<TH2>(Form("ULS_MeePtee_%s_%s", typeNamesPairULS, stageNames[j]), Form("ULS MeePtee %s %s", typeNamesPairULS, stageNames[j]), HistType::kTH2F, {mass_axis, ptee_axis}, true));
299371
}
372+
for (int j = 0; j < 3; j++) {
373+
hLS_Mee.push_back(registry.add<TH1>(Form("LS_Mee_%s_%s", typeNamesPairULS, stageNames[j]), Form("LS Mee %s %s", typeNamesPairULS, stageNames[j]), HistType::kTH1F, {mass_axis}, true));
374+
hLS_MeePtee.push_back(registry.add<TH2>(Form("LS_MeePtee_%s_%s", typeNamesPairULS, stageNames[j]), Form("LS MeePtee %s %s", typeNamesPairULS, stageNames[j]), HistType::kTH2F, {mass_axis, ptee_axis}, true));
375+
}
300376
}
301377

302378
void processCharm(aod::McCollisions const& collisions, MyFilteredMcParticlesSmeared const& mcParticles, aod::McParticles const& mcParticlesAll)
303379
{
380+
for (auto const& p : mcParticles) {
381+
if (myConfigs.fConfigCheckPartonic && p.cQuarkOriginId() < 0) {
382+
continue;
383+
}
384+
doSingle(p, hEta, hPt, hPtEta, myConfigs.fConfigPtMin, myConfigs.fConfigEtaMax);
385+
}
386+
304387
for (auto const& collision : collisions) {
305388

306389
registry.fill(HIST("NEvents"), 0.5);
307390

308-
for (auto const& p : mcParticles) {
309-
doSingle(p, hEta, hPt, hPtEta, fConfigPtMin, fConfigEtaMax);
310-
}
311-
312391
auto const electronsGrouped = Electrons->sliceBy(perCollision, collision.globalIndex());
313392
auto const positronsGrouped = Positrons->sliceBy(perCollision, collision.globalIndex());
314393
// ULS spectrum
315394
for (auto const& [particle1, particle2] : combinations(o2::soa::CombinationsFullIndexPolicy(electronsGrouped, positronsGrouped))) {
316-
int type = IsHF(particle1, particle2, mcParticlesAll);
317-
if (type == static_cast<int>(EM_HFeeType::kUndef))
318-
continue;
319-
if (type != static_cast<int>(EM_HFeeType::kCe_Ce)) {
395+
if (myConfigs.fConfigCheckPartonic) {
396+
if (particle1.cQuarkOriginId() < 0 || particle2.cQuarkOriginId() < 0 || particle1.cQuarkOriginId() != particle2.cQuarkOriginId())
397+
continue;
398+
}
399+
if (IsHF(particle1, particle2, mcParticlesAll) != static_cast<int>(EM_HFeeType::kCe_Ce)) {
320400
LOG(error) << "Something is wrong here. There should only be pairs of type kCe_Ce = 0 left at this point.";
321401
}
322-
doPair(particle1, particle2, hULS_Mee, hULS_MeePtee, fConfigPtMin, fConfigEtaMax);
402+
doPair(particle1, particle2, hULS_Mee, hULS_MeePtee, myConfigs.fConfigPtMin, myConfigs.fConfigEtaMax);
403+
}
404+
// LS
405+
for (auto const& [particle1, particle2] : combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(electronsGrouped, electronsGrouped))) {
406+
if (myConfigs.fConfigCheckPartonic) {
407+
if (particle1.cQuarkOriginId() < 0 || particle2.cQuarkOriginId() < 0 || particle1.cQuarkOriginId() != particle2.cQuarkOriginId())
408+
continue;
409+
}
410+
doPair(particle1, particle2, hLS_Mee, hLS_MeePtee, myConfigs.fConfigPtMin, myConfigs.fConfigEtaMax);
411+
}
412+
for (auto const& [particle1, particle2] : combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(positronsGrouped, positronsGrouped))) {
413+
if (myConfigs.fConfigCheckPartonic) {
414+
if (particle1.cQuarkOriginId() < 0 || particle2.cQuarkOriginId() < 0 || particle1.cQuarkOriginId() != particle2.cQuarkOriginId())
415+
continue;
416+
}
417+
doPair(particle1, particle2, hLS_Mee, hLS_MeePtee, myConfigs.fConfigPtMin, myConfigs.fConfigEtaMax);
323418
}
324419
}
325420
}

0 commit comments

Comments
 (0)