Skip to content

PWGEM/PhotonMeson: update for HF->ee MC #5722

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
dsekihat merged 1 commit into from
Apr 19, 2024
Merged

PWGEM/PhotonMeson: update for HF->ee MC #5722

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
166 changes: 134 additions & 32 deletions PWGEM/Dilepton/Utils/MCUtilities.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -30,14 +30,69 @@ enum class EM_HFeeType : int {
kBCe_Be_SameB = 3, // ULS
kBCe_Be_DiffB = 4, // LS
};

template <typename TMCParticle, typename TMCParticles>
int IsFromBeauty(TMCParticle const& p, TMCParticles const& mcparticles)
{
if (!p.has_mothers()) {
return -999;
}

int motherid = p.mothersIds()[0]; // first mother index
while (motherid > -1) {
if (motherid < mcparticles.size()) { // protect against bad mother indices. why is this needed?
auto mp = mcparticles.iteratorAt(motherid);
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')) {
return motherid;
}
if (mp.has_mothers()) {
motherid = mp.mothersIds()[0];
} else {
return -999;
}
} else {
LOGF(info, "Mother label(%d) exceeds the McParticles size(%d)", motherid, mcparticles.size());
}
}

return -999;
}

template <typename TMCParticle, typename TMCParticles>
int IsFromCharm(TMCParticle const& p, TMCParticles const& mcparticles)
{
if (!p.has_mothers()) {
return -999;
}

int motherid = p.mothersIds()[0]; // first mother index
while (motherid > -1) {
if (motherid < mcparticles.size()) { // protect against bad mother indices. why is this needed?
auto mp = mcparticles.iteratorAt(motherid);
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')) {
return motherid;
}
if (mp.has_mothers()) {
motherid = mp.mothersIds()[0];
} else {
return -999;
}
} else {
LOGF(info, "Mother label(%d) exceeds the McParticles size(%d)", motherid, mcparticles.size());
}
}

return -999;
}

template <typename TMCParticle1, typename TMCParticle2, typename TMCParticles>
int IsHF(TMCParticle1 const& p1, TMCParticle2 const& p2, TMCParticles const& mcparticles)
{
if (!p1.has_mothers() || !p2.has_mothers()) {
return static_cast<int>(EM_HFeeType::kUndef);
}

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

Expand Down Expand Up @@ -81,44 +136,91 @@ int IsHF(TMCParticle1 const& p1, TMCParticle2 const& p2, TMCParticles const& mcp
}
}

if (std::to_string(mothers_pdg1[0]).find("5") != std::string::npos && std::to_string(mothers_pdg2[0]).find("5") != std::string::npos) {
return static_cast<int>(EM_HFeeType::kBe_Be); // bb->ee, decay type = 2
// this is easy. first mother is b hadron for both leg.
}

if (std::to_string(mothers_pdg1[0]).find("4") != std::string::npos && std::to_string(mothers_pdg2[0]).find("4") != std::string::npos) {
// mother is c hadron. next, check c is prompt or non-prompt.
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');
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');
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;
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;
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;
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;

bool is_c_from_b1 = false;
for (unsigned int i1 = 1; i1 < mothers_pdg1.size(); i1++) {
if (std::to_string(mothers_pdg1[i1]).find("5") != std::string::npos) {
is_c_from_b1 = true;
break;
}
}
bool is_c_from_b2 = false;
for (unsigned int i2 = 1; i2 < mothers_pdg2.size(); i2++) {
if (std::to_string(mothers_pdg2[i2]).find("5") != std::string::npos) {
is_c_from_b2 = true;
break;
}
}

if (!is_c_from_b1 && !is_c_from_b2) {
return static_cast<int>(EM_HFeeType::kCe_Ce); // prompt cc->ee, decay type = 0
} else if (is_c_from_b1 && is_c_from_b2) {
return static_cast<int>(EM_HFeeType::kBCe_BCe); // b->c->e and b->c->e, decay type = 1
} else {
if (is_prompt_c1 && is_prompt_c2 && p1.pdgCode() * p2.pdgCode() < 0) {
mothers_id1.clear();
mothers_pdg1.clear();
mothers_id2.clear();
mothers_pdg2.clear();
mothers_id1.shrink_to_fit();
mothers_pdg1.shrink_to_fit();
mothers_id2.shrink_to_fit();
mothers_pdg2.shrink_to_fit();
return static_cast<int>(EM_HFeeType::kCe_Ce); // cc->ee, decay type = 0
} else if (is_direct_from_b1 && is_direct_from_b2 && p1.pdgCode() * p2.pdgCode() < 0) {
mothers_id1.clear();
mothers_pdg1.clear();
mothers_id2.clear();
mothers_pdg2.clear();
mothers_id1.shrink_to_fit();
mothers_pdg1.shrink_to_fit();
mothers_id2.shrink_to_fit();
mothers_pdg2.shrink_to_fit();
return static_cast<int>(EM_HFeeType::kBe_Be); // bb->ee, decay type = 2
} else if (is_c_from_b1 && is_c_from_b2 && p1.pdgCode() * p2.pdgCode() < 0) {
mothers_id1.clear();
mothers_pdg1.clear();
mothers_id2.clear();
mothers_pdg2.clear();
mothers_id1.shrink_to_fit();
mothers_pdg1.shrink_to_fit();
mothers_id2.shrink_to_fit();
mothers_pdg2.shrink_to_fit();
return static_cast<int>(EM_HFeeType::kBCe_BCe); // b->c->e and b->c->e, decay type = 1
} else if ((is_direct_from_b1 && is_c_from_b2) || (is_direct_from_b2 && is_c_from_b1)) {
if (p1.pdgCode() * p2.pdgCode() < 0) { // ULS
for (auto& mid1 : mothers_id1) {
for (auto& mid2 : mothers_id2) {
if (mid1 == mid2) {
auto common_mp = mcparticles.iteratorAt(mid1);
int mp_pdg = common_mp.pdgCode();
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')) {
mothers_id1.clear();
mothers_pdg1.clear();
mothers_id2.clear();
mothers_pdg2.clear();
mothers_id1.shrink_to_fit();
mothers_pdg1.shrink_to_fit();
mothers_id2.shrink_to_fit();
mothers_pdg2.shrink_to_fit();
return static_cast<int>(EM_HFeeType::kBCe_Be_SameB); // b->c->e and b->e, decay type = 3. this should happen only in ULS.
}
}
} // end of motherid2
} // end of motherid1
} else { // LS
bool is_same_mother_found = false;
for (auto& mid1 : mothers_id1) {
for (auto& mid2 : mothers_id2) {
if (mid1 == mid2) {
return static_cast<int>(EM_HFeeType::kBCe_Be_SameB); // b->c->e and c->e, decay type = 3. this should happen only in ULS.
auto common_mp = mcparticles.iteratorAt(mid1);
int mp_pdg = common_mp.pdgCode();
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')) {
is_same_mother_found = true;
}
}
} // end of mother id 2
} // end of mother id 1
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
} // end of motherid2
} // end of motherid1
if (!is_same_mother_found) {
mothers_id1.clear();
mothers_pdg1.clear();
mothers_id2.clear();
mothers_pdg2.clear();
mothers_id1.shrink_to_fit();
mothers_pdg1.shrink_to_fit();
mothers_id2.shrink_to_fit();
mothers_pdg2.shrink_to_fit();
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
}
}
}

mothers_id1.clear();
mothers_pdg1.clear();
mothers_id2.clear();
Expand Down
53 changes: 37 additions & 16 deletions PWGEM/PhotonMeson/Core/CutsLibrary.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -41,35 +41,56 @@ EMEventCut* o2::aod::pwgem::photon::eventcuts::GetCut(const char* cutName)
EMEventCut* cut = new EMEventCut(cutName, cutName);
std::string nameStr = cutName;

if (!nameStr.compare("minbias")) {
cut->SetRequireFT0AND(true);
cut->SetZvtxRange(-10.f, +10.f);
cut->SetRequireNoTFB(false);
cut->SetRequireNoITSROFB(false);
return cut;
}

if (!nameStr.compare("nocut")) {
cut->SetRequireSel8(false);
cut->SetRequireFT0AND(false);
cut->SetZvtxRange(-1e+10, +1e+10);
cut->SetRequireNoTFB(false);
cut->SetRequireNoITSROFB(false);
cut->SetRequireNoSameBunchPileup(false);
cut->SetRequireVertexITSTPC(false);
cut->SetRequireIsGoodZvtxFT0vsPV(false);
return cut;
}

if (!nameStr.compare("minbias_notfb")) {
if (!nameStr.compare("ft0and")) {
cut->SetRequireSel8(false);
cut->SetRequireFT0AND(true);
cut->SetZvtxRange(-10.f, +10.f);
cut->SetRequireNoTFB(true);
cut->SetRequireNoTFB(false);
cut->SetRequireNoITSROFB(false);
cut->SetRequireNoSameBunchPileup(false);
cut->SetRequireVertexITSTPC(false);
cut->SetRequireIsGoodZvtxFT0vsPV(false);
return cut;
}

if (!nameStr.compare("minbias_notfb_noitsrofb")) {
if (!nameStr.compare("minbias")) {
cut->SetRequireSel8(true);
cut->SetRequireFT0AND(true);
cut->SetZvtxRange(-10.f, +10.f);
cut->SetRequireNoTFB(true);
cut->SetRequireNoITSROFB(true);
cut->SetRequireNoTFB(false); // included in sel8
cut->SetRequireNoITSROFB(false); // included in sel8
cut->SetRequireNoSameBunchPileup(false);
cut->SetRequireVertexITSTPC(false);
cut->SetRequireIsGoodZvtxFT0vsPV(false);

if (nameStr.find("notfb") != std::string::npos) {
cut->SetRequireNoTFB(true);
}
if (nameStr.find("noitsrofb") != std::string::npos) {
cut->SetRequireNoITSROFB(true);
}
if (nameStr.find("nosbp") != std::string::npos) {
cut->SetRequireNoSameBunchPileup(true);
}
if (nameStr.find("vtxitstpc") != std::string::npos) {
cut->SetRequireVertexITSTPC(true);
}
if (nameStr.find("goodvtx") != std::string::npos) {
cut->SetRequireIsGoodZvtxFT0vsPV(true);
}

return cut;
}

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

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

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

// for track cuts
cut->SetMinNCrossedRowsTPC(100);
cut->SetMinNCrossedRowsTPC(40);
cut->SetMinNCrossedRowsOverFindableClustersTPC(0.8);
cut->SetChi2PerClusterTPC(0.0, 4.0);
cut->SetChi2PerClusterITS(0.0, 5.0);
cut->SetNClustersITS(5, 7);
cut->SetNClustersITS(4, 7);
cut->SetMaxDcaXY(1.0);
cut->SetMaxDcaZ(1.0);
return cut;
Expand Down
26 changes: 25 additions & 1 deletion PWGEM/PhotonMeson/Core/EMEventCut.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,7 +18,13 @@

ClassImp(EMEventCut);

const char* EMEventCut::mCutNames[static_cast<int>(EMEventCut::EMEventCuts::kNCuts)] = {"RequireFT0AND", "Zvtx", "RequireNoTFB", "RequireNoITSROFB"};
const char* EMEventCut::mCutNames[static_cast<int>(EMEventCut::EMEventCuts::kNCuts)] = {"Sel8", "FT0AND", "Zvtx", "eNoTFB", "RequireNoITSROFB", "NoSameBunchPileup", "GoodVertexITSTPC", "GoodZvtxFT0vsPV"};

void EMEventCut::SetRequireSel8(bool flag)
{
mRequireSel8 = flag;
LOG(info) << "EM Event Cut, require sel8: " << mRequireSel8;
}

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

void EMEventCut::SetRequireNoSameBunchPileup(bool flag)
{
mRequireNoSameBunchPileup = flag;
LOG(info) << "EM Event Cut, require No same bunch pileup: " << mRequireNoSameBunchPileup;
}

void EMEventCut::SetRequireVertexITSTPC(bool flag)
{
mRequireVertexITSTPC = flag;
LOG(info) << "EM Event Cut, require vertex reconstructed by ITS-TPC matched track: " << mRequireVertexITSTPC;
}

void EMEventCut::SetRequireIsGoodZvtxFT0vsPV(bool flag)
{
mRequireGoodZvtxFT0vsPV = flag;
LOG(info) << "EM Event Cut, require good Zvtx between FT0 vs. PV: " << mRequireGoodZvtxFT0vsPV;
}

void EMEventCut::print() const
{
LOG(info) << "EM Event Cut:";
Expand Down
Loading