Skip to content

Commit ed4aea6

Browse files
dsekihatdsekihat
authored
PWGEM/PhotonMeson: update for HF->ee MC (#5722)
1 parent cc95e49 commit ed4aea6

7 files changed

Lines changed: 675 additions & 90 deletions

File tree

‎PWGEM/Dilepton/Utils/MCUtilities.h‎

Lines changed: 134 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -30,14 +30,69 @@ enum class EM_HFeeType : int {
3030
kBCe_Be_SameB = 3, // ULS
3131
kBCe_Be_DiffB = 4, // LS
3232
};
33+
34+
template <typename TMCParticle, typename TMCParticles>
35+
int IsFromBeauty(TMCParticle const& p, TMCParticles const& mcparticles)
36+
{
37+
if (!p.has_mothers()) {
38+
return -999;
39+
}
40+
41+
int motherid = p.mothersIds()[0]; // first mother index
42+
while (motherid > -1) {
43+
if (motherid < mcparticles.size()) { // protect against bad mother indices. why is this needed?
44+
auto mp = mcparticles.iteratorAt(motherid);
45+
if (abs(mp.pdgCode()) < 1e+9 && (std::to_string(abs(mp.pdgCode()))[std::to_string(abs(mp.pdgCode())).length() - 3] == '5' || std::to_string(abs(mp.pdgCode()))[std::to_string(abs(mp.pdgCode())).length() - 4] == '5')) {
46+
return motherid;
47+
}
48+
if (mp.has_mothers()) {
49+
motherid = mp.mothersIds()[0];
50+
} else {
51+
return -999;
52+
}
53+
} else {
54+
LOGF(info, "Mother label(%d) exceeds the McParticles size(%d)", motherid, mcparticles.size());
55+
}
56+
}
57+
58+
return -999;
59+
}
60+
61+
template <typename TMCParticle, typename TMCParticles>
62+
int IsFromCharm(TMCParticle const& p, TMCParticles const& mcparticles)
63+
{
64+
if (!p.has_mothers()) {
65+
return -999;
66+
}
67+
68+
int motherid = p.mothersIds()[0]; // first mother index
69+
while (motherid > -1) {
70+
if (motherid < mcparticles.size()) { // protect against bad mother indices. why is this needed?
71+
auto mp = mcparticles.iteratorAt(motherid);
72+
if (abs(mp.pdgCode()) < 1e+9 && (std::to_string(abs(mp.pdgCode()))[std::to_string(abs(mp.pdgCode())).length() - 3] == '4' || std::to_string(abs(mp.pdgCode()))[std::to_string(abs(mp.pdgCode())).length() - 4] == '4')) {
73+
return motherid;
74+
}
75+
if (mp.has_mothers()) {
76+
motherid = mp.mothersIds()[0];
77+
} else {
78+
return -999;
79+
}
80+
} else {
81+
LOGF(info, "Mother label(%d) exceeds the McParticles size(%d)", motherid, mcparticles.size());
82+
}
83+
}
84+
85+
return -999;
86+
}
87+
3388
template <typename TMCParticle1, typename TMCParticle2, typename TMCParticles>
3489
int IsHF(TMCParticle1 const& p1, TMCParticle2 const& p2, TMCParticles const& mcparticles)
3590
{
3691
if (!p1.has_mothers() || !p2.has_mothers()) {
3792
return static_cast<int>(EM_HFeeType::kUndef);
3893
}
3994

40-
if (p1.mothersIds()[0] == p2.mothersIds()[0]) { // same mother
95+
if (p1.mothersIds()[0] == p2.mothersIds()[0]) { // reject same mother. e.g. jspi 443
4196
return static_cast<int>(EM_HFeeType::kUndef); // this never happens in correlated HF->ee decays
4297
}
4398

@@ -81,44 +136,91 @@ int IsHF(TMCParticle1 const& p1, TMCParticle2 const& p2, TMCParticles const& mcp
81136
}
82137
}
83138

84-
if (std::to_string(mothers_pdg1[0]).find("5") != std::string::npos && std::to_string(mothers_pdg2[0]).find("5") != std::string::npos) {
85-
return static_cast<int>(EM_HFeeType::kBe_Be); // bb->ee, decay type = 2
86-
// this is easy. first mother is b hadron for both leg.
87-
}
88-
89-
if (std::to_string(mothers_pdg1[0]).find("4") != std::string::npos && std::to_string(mothers_pdg2[0]).find("4") != std::string::npos) {
90-
// mother is c hadron. next, check c is prompt or non-prompt.
139+
bool is_direct_from_b1 = abs(mothers_pdg1[0]) < 1e+9 && (std::to_string(mothers_pdg1[0])[std::to_string(mothers_pdg1[0]).length() - 3] == '5' || std::to_string(mothers_pdg1[0])[std::to_string(mothers_pdg1[0]).length() - 4] == '5');
140+
bool is_direct_from_b2 = abs(mothers_pdg2[0]) < 1e+9 && (std::to_string(mothers_pdg2[0])[std::to_string(mothers_pdg2[0]).length() - 3] == '5' || std::to_string(mothers_pdg2[0])[std::to_string(mothers_pdg2[0]).length() - 4] == '5');
141+
bool is_prompt_c1 = abs(mothers_pdg1[0]) < 1e+9 && (std::to_string(mothers_pdg1[0])[std::to_string(mothers_pdg1[0]).length() - 3] == '4' || std::to_string(mothers_pdg1[0])[std::to_string(mothers_pdg1[0]).length() - 4] == '4') && IsFromBeauty(p1, mcparticles) < 0;
142+
bool is_prompt_c2 = abs(mothers_pdg2[0]) < 1e+9 && (std::to_string(mothers_pdg2[0])[std::to_string(mothers_pdg2[0]).length() - 3] == '4' || std::to_string(mothers_pdg2[0])[std::to_string(mothers_pdg2[0]).length() - 4] == '4') && IsFromBeauty(p2, mcparticles) < 0;
143+
bool is_c_from_b1 = abs(mothers_pdg1[0]) < 1e+9 && (std::to_string(mothers_pdg1[0])[std::to_string(mothers_pdg1[0]).length() - 3] == '4' || std::to_string(mothers_pdg1[0])[std::to_string(mothers_pdg1[0]).length() - 4] == '4') && IsFromBeauty(p1, mcparticles) > 0;
144+
bool is_c_from_b2 = abs(mothers_pdg2[0]) < 1e+9 && (std::to_string(mothers_pdg2[0])[std::to_string(mothers_pdg2[0]).length() - 3] == '4' || std::to_string(mothers_pdg2[0])[std::to_string(mothers_pdg2[0]).length() - 4] == '4') && IsFromBeauty(p2, mcparticles) > 0;
91145

92-
bool is_c_from_b1 = false;
93-
for (unsigned int i1 = 1; i1 < mothers_pdg1.size(); i1++) {
94-
if (std::to_string(mothers_pdg1[i1]).find("5") != std::string::npos) {
95-
is_c_from_b1 = true;
96-
break;
97-
}
98-
}
99-
bool is_c_from_b2 = false;
100-
for (unsigned int i2 = 1; i2 < mothers_pdg2.size(); i2++) {
101-
if (std::to_string(mothers_pdg2[i2]).find("5") != std::string::npos) {
102-
is_c_from_b2 = true;
103-
break;
104-
}
105-
}
106-
107-
if (!is_c_from_b1 && !is_c_from_b2) {
108-
return static_cast<int>(EM_HFeeType::kCe_Ce); // prompt cc->ee, decay type = 0
109-
} else if (is_c_from_b1 && is_c_from_b2) {
110-
return static_cast<int>(EM_HFeeType::kBCe_BCe); // b->c->e and b->c->e, decay type = 1
111-
} else {
146+
if (is_prompt_c1 && is_prompt_c2 && p1.pdgCode() * p2.pdgCode() < 0) {
147+
mothers_id1.clear();
148+
mothers_pdg1.clear();
149+
mothers_id2.clear();
150+
mothers_pdg2.clear();
151+
mothers_id1.shrink_to_fit();
152+
mothers_pdg1.shrink_to_fit();
153+
mothers_id2.shrink_to_fit();
154+
mothers_pdg2.shrink_to_fit();
155+
return static_cast<int>(EM_HFeeType::kCe_Ce); // cc->ee, decay type = 0
156+
} else if (is_direct_from_b1 && is_direct_from_b2 && p1.pdgCode() * p2.pdgCode() < 0) {
157+
mothers_id1.clear();
158+
mothers_pdg1.clear();
159+
mothers_id2.clear();
160+
mothers_pdg2.clear();
161+
mothers_id1.shrink_to_fit();
162+
mothers_pdg1.shrink_to_fit();
163+
mothers_id2.shrink_to_fit();
164+
mothers_pdg2.shrink_to_fit();
165+
return static_cast<int>(EM_HFeeType::kBe_Be); // bb->ee, decay type = 2
166+
} else if (is_c_from_b1 && is_c_from_b2 && p1.pdgCode() * p2.pdgCode() < 0) {
167+
mothers_id1.clear();
168+
mothers_pdg1.clear();
169+
mothers_id2.clear();
170+
mothers_pdg2.clear();
171+
mothers_id1.shrink_to_fit();
172+
mothers_pdg1.shrink_to_fit();
173+
mothers_id2.shrink_to_fit();
174+
mothers_pdg2.shrink_to_fit();
175+
return static_cast<int>(EM_HFeeType::kBCe_BCe); // b->c->e and b->c->e, decay type = 1
176+
} else if ((is_direct_from_b1 && is_c_from_b2) || (is_direct_from_b2 && is_c_from_b1)) {
177+
if (p1.pdgCode() * p2.pdgCode() < 0) { // ULS
178+
for (auto& mid1 : mothers_id1) {
179+
for (auto& mid2 : mothers_id2) {
180+
if (mid1 == mid2) {
181+
auto common_mp = mcparticles.iteratorAt(mid1);
182+
int mp_pdg = common_mp.pdgCode();
183+
if (abs(mp_pdg) < 1e+9 && (std::to_string(abs(mp_pdg))[std::to_string(abs(mp_pdg)).length() - 3] == '5' || std::to_string(abs(mp_pdg))[std::to_string(abs(mp_pdg)).length() - 4] == '5')) {
184+
mothers_id1.clear();
185+
mothers_pdg1.clear();
186+
mothers_id2.clear();
187+
mothers_pdg2.clear();
188+
mothers_id1.shrink_to_fit();
189+
mothers_pdg1.shrink_to_fit();
190+
mothers_id2.shrink_to_fit();
191+
mothers_pdg2.shrink_to_fit();
192+
return static_cast<int>(EM_HFeeType::kBCe_Be_SameB); // b->c->e and b->e, decay type = 3. this should happen only in ULS.
193+
}
194+
}
195+
} // end of motherid2
196+
} // end of motherid1
197+
} else { // LS
198+
bool is_same_mother_found = false;
112199
for (auto& mid1 : mothers_id1) {
113200
for (auto& mid2 : mothers_id2) {
114201
if (mid1 == mid2) {
115-
return static_cast<int>(EM_HFeeType::kBCe_Be_SameB); // b->c->e and c->e, decay type = 3. this should happen only in ULS.
202+
auto common_mp = mcparticles.iteratorAt(mid1);
203+
int mp_pdg = common_mp.pdgCode();
204+
if (abs(mp_pdg) < 1e+9 && (std::to_string(abs(mp_pdg))[std::to_string(abs(mp_pdg)).length() - 3] == '5' || std::to_string(abs(mp_pdg))[std::to_string(abs(mp_pdg)).length() - 4] == '5')) {
205+
is_same_mother_found = true;
206+
}
116207
}
117-
} // end of mother id 2
118-
} // end of mother id 1
119-
return static_cast<int>(EM_HFeeType::kBCe_Be_DiffB); // b->c->e and c->e, decay type = 4. this should happen only in LS. But, this may happen, when ele/pos is reconstructed as pos/ele wrongly. and create LS pair
208+
} // end of motherid2
209+
} // end of motherid1
210+
if (!is_same_mother_found) {
211+
mothers_id1.clear();
212+
mothers_pdg1.clear();
213+
mothers_id2.clear();
214+
mothers_pdg2.clear();
215+
mothers_id1.shrink_to_fit();
216+
mothers_pdg1.shrink_to_fit();
217+
mothers_id2.shrink_to_fit();
218+
mothers_pdg2.shrink_to_fit();
219+
return static_cast<int>(EM_HFeeType::kBCe_Be_DiffB); // b->c->e and b->e, decay type = 4. this should happen only in LS. But, this may happen, when ele/pos is reconstructed as pos/ele wrongly. and create LS pair
220+
}
120221
}
121222
}
223+
122224
mothers_id1.clear();
123225
mothers_pdg1.clear();
124226
mothers_id2.clear();

‎PWGEM/PhotonMeson/Core/CutsLibrary.cxx‎

Lines changed: 37 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -41,35 +41,56 @@ EMEventCut* o2::aod::pwgem::photon::eventcuts::GetCut(const char* cutName)
4141
EMEventCut* cut = new EMEventCut(cutName, cutName);
4242
std::string nameStr = cutName;
4343

44-
if (!nameStr.compare("minbias")) {
45-
cut->SetRequireFT0AND(true);
46-
cut->SetZvtxRange(-10.f, +10.f);
47-
cut->SetRequireNoTFB(false);
48-
cut->SetRequireNoITSROFB(false);
49-
return cut;
50-
}
51-
5244
if (!nameStr.compare("nocut")) {
45+
cut->SetRequireSel8(false);
5346
cut->SetRequireFT0AND(false);
5447
cut->SetZvtxRange(-1e+10, +1e+10);
5548
cut->SetRequireNoTFB(false);
5649
cut->SetRequireNoITSROFB(false);
50+
cut->SetRequireNoSameBunchPileup(false);
51+
cut->SetRequireVertexITSTPC(false);
52+
cut->SetRequireIsGoodZvtxFT0vsPV(false);
5753
return cut;
5854
}
5955

60-
if (!nameStr.compare("minbias_notfb")) {
56+
if (!nameStr.compare("ft0and")) {
57+
cut->SetRequireSel8(false);
6158
cut->SetRequireFT0AND(true);
6259
cut->SetZvtxRange(-10.f, +10.f);
63-
cut->SetRequireNoTFB(true);
60+
cut->SetRequireNoTFB(false);
6461
cut->SetRequireNoITSROFB(false);
62+
cut->SetRequireNoSameBunchPileup(false);
63+
cut->SetRequireVertexITSTPC(false);
64+
cut->SetRequireIsGoodZvtxFT0vsPV(false);
6565
return cut;
6666
}
6767

68-
if (!nameStr.compare("minbias_notfb_noitsrofb")) {
68+
if (!nameStr.compare("minbias")) {
69+
cut->SetRequireSel8(true);
6970
cut->SetRequireFT0AND(true);
7071
cut->SetZvtxRange(-10.f, +10.f);
71-
cut->SetRequireNoTFB(true);
72-
cut->SetRequireNoITSROFB(true);
72+
cut->SetRequireNoTFB(false); // included in sel8
73+
cut->SetRequireNoITSROFB(false); // included in sel8
74+
cut->SetRequireNoSameBunchPileup(false);
75+
cut->SetRequireVertexITSTPC(false);
76+
cut->SetRequireIsGoodZvtxFT0vsPV(false);
77+
78+
if (nameStr.find("notfb") != std::string::npos) {
79+
cut->SetRequireNoTFB(true);
80+
}
81+
if (nameStr.find("noitsrofb") != std::string::npos) {
82+
cut->SetRequireNoITSROFB(true);
83+
}
84+
if (nameStr.find("nosbp") != std::string::npos) {
85+
cut->SetRequireNoSameBunchPileup(true);
86+
}
87+
if (nameStr.find("vtxitstpc") != std::string::npos) {
88+
cut->SetRequireVertexITSTPC(true);
89+
}
90+
if (nameStr.find("goodvtx") != std::string::npos) {
91+
cut->SetRequireIsGoodZvtxFT0vsPV(true);
92+
}
93+
7394
return cut;
7495
}
7596

@@ -424,7 +445,7 @@ DalitzEECut* o2::aod::pwgem::photon::dalitzeecuts::GetCut(const char* cutName)
424445

425446
if (!nameStr.compare("nocut")) {
426447
// apply kinetic cuts
427-
cut->SetTrackPtRange(0.05, 1e+10f);
448+
cut->SetTrackPtRange(0.1, 1e+10f);
428449
cut->SetTrackEtaRange(-0.9, +0.9);
429450

430451
// for pair
@@ -434,11 +455,11 @@ DalitzEECut* o2::aod::pwgem::photon::dalitzeecuts::GetCut(const char* cutName)
434455
cut->ApplyPrefilter(false);
435456

436457
// for track cuts
437-
cut->SetMinNCrossedRowsTPC(100);
458+
cut->SetMinNCrossedRowsTPC(40);
438459
cut->SetMinNCrossedRowsOverFindableClustersTPC(0.8);
439460
cut->SetChi2PerClusterTPC(0.0, 4.0);
440461
cut->SetChi2PerClusterITS(0.0, 5.0);
441-
cut->SetNClustersITS(5, 7);
462+
cut->SetNClustersITS(4, 7);
442463
cut->SetMaxDcaXY(1.0);
443464
cut->SetMaxDcaZ(1.0);
444465
return cut;

‎PWGEM/PhotonMeson/Core/EMEventCut.cxx‎

Lines changed: 25 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,13 @@
1818

1919
ClassImp(EMEventCut);
2020

21-
const char* EMEventCut::mCutNames[static_cast<int>(EMEventCut::EMEventCuts::kNCuts)] = {"RequireFT0AND", "Zvtx", "RequireNoTFB", "RequireNoITSROFB"};
21+
const char* EMEventCut::mCutNames[static_cast<int>(EMEventCut::EMEventCuts::kNCuts)] = {"Sel8", "FT0AND", "Zvtx", "eNoTFB", "RequireNoITSROFB", "NoSameBunchPileup", "GoodVertexITSTPC", "GoodZvtxFT0vsPV"};
22+
23+
void EMEventCut::SetRequireSel8(bool flag)
24+
{
25+
mRequireSel8 = flag;
26+
LOG(info) << "EM Event Cut, require sel8: " << mRequireSel8;
27+
}
2228

2329
void EMEventCut::SetRequireFT0AND(bool flag)
2430
{
@@ -45,6 +51,24 @@ void EMEventCut::SetRequireNoITSROFB(bool flag)
4551
LOG(info) << "EM Event Cut, require No ITS ROF border: " << mRequireNoITSROFB;
4652
}
4753

54+
void EMEventCut::SetRequireNoSameBunchPileup(bool flag)
55+
{
56+
mRequireNoSameBunchPileup = flag;
57+
LOG(info) << "EM Event Cut, require No same bunch pileup: " << mRequireNoSameBunchPileup;
58+
}
59+
60+
void EMEventCut::SetRequireVertexITSTPC(bool flag)
61+
{
62+
mRequireVertexITSTPC = flag;
63+
LOG(info) << "EM Event Cut, require vertex reconstructed by ITS-TPC matched track: " << mRequireVertexITSTPC;
64+
}
65+
66+
void EMEventCut::SetRequireIsGoodZvtxFT0vsPV(bool flag)
67+
{
68+
mRequireGoodZvtxFT0vsPV = flag;
69+
LOG(info) << "EM Event Cut, require good Zvtx between FT0 vs. PV: " << mRequireGoodZvtxFT0vsPV;
70+
}
71+
4872
void EMEventCut::print() const
4973
{
5074
LOG(info) << "EM Event Cut:";

0 commit comments

Comments
 (0)