diff --git a/PWGEM/Dilepton/Utils/MCUtilities.h b/PWGEM/Dilepton/Utils/MCUtilities.h index d1b54e7b068..a4cdb572dce 100644 --- a/PWGEM/Dilepton/Utils/MCUtilities.h +++ b/PWGEM/Dilepton/Utils/MCUtilities.h @@ -30,6 +30,61 @@ enum class EM_HFeeType : int { kBCe_Be_SameB = 3, // ULS kBCe_Be_DiffB = 4, // LS }; + +template +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 +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 int IsHF(TMCParticle1 const& p1, TMCParticle2 const& p2, TMCParticles const& mcparticles) { @@ -37,7 +92,7 @@ int IsHF(TMCParticle1 const& p1, TMCParticle2 const& p2, TMCParticles const& mcp return static_cast(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(EM_HFeeType::kUndef); // this never happens in correlated HF->ee decays } @@ -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(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(EM_HFeeType::kCe_Ce); // prompt cc->ee, decay type = 0 - } else if (is_c_from_b1 && is_c_from_b2) { - return static_cast(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(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(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(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(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(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(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(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(); diff --git a/PWGEM/PhotonMeson/Core/CutsLibrary.cxx b/PWGEM/PhotonMeson/Core/CutsLibrary.cxx index 6583fe589ba..be3a6297fd4 100644 --- a/PWGEM/PhotonMeson/Core/CutsLibrary.cxx +++ b/PWGEM/PhotonMeson/Core/CutsLibrary.cxx @@ -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; } @@ -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 @@ -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; diff --git a/PWGEM/PhotonMeson/Core/EMEventCut.cxx b/PWGEM/PhotonMeson/Core/EMEventCut.cxx index b605d2c7e8b..244dc762dc5 100644 --- a/PWGEM/PhotonMeson/Core/EMEventCut.cxx +++ b/PWGEM/PhotonMeson/Core/EMEventCut.cxx @@ -18,7 +18,13 @@ ClassImp(EMEventCut); -const char* EMEventCut::mCutNames[static_cast(EMEventCut::EMEventCuts::kNCuts)] = {"RequireFT0AND", "Zvtx", "RequireNoTFB", "RequireNoITSROFB"}; +const char* EMEventCut::mCutNames[static_cast(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) { @@ -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:"; diff --git a/PWGEM/PhotonMeson/Core/EMEventCut.h b/PWGEM/PhotonMeson/Core/EMEventCut.h index 8fc86e0bb93..cb225ad1c24 100644 --- a/PWGEM/PhotonMeson/Core/EMEventCut.h +++ b/PWGEM/PhotonMeson/Core/EMEventCut.h @@ -29,10 +29,14 @@ class EMEventCut : public TNamed EMEventCut(const char* name, const char* title) : TNamed(name, title) {} enum class EMEventCuts : int { - kFT0AND = 0, // i.e. sel8 + kSel8 = 0, + kFT0AND, kZvtx, kNoTFB, // no time frame border kNoITSROFB, // no ITS read out frame border + kNoSameBunchPileup, + kIsVertexITSTPC, + kIsGoodZvtxFT0vsPV, kNCuts }; @@ -41,6 +45,9 @@ class EMEventCut : public TNamed template bool IsSelected(T const& collision) const { + if (mRequireSel8 && !IsSelected(collision, EMEventCuts::kSel8)) { + return false; + } if (mRequireFT0AND && !IsSelected(collision, EMEventCuts::kFT0AND)) { return false; } @@ -53,6 +60,15 @@ class EMEventCut : public TNamed if (mRequireNoITSROFB && !IsSelected(collision, EMEventCuts::kNoITSROFB)) { return false; } + if (mRequireNoSameBunchPileup && !IsSelected(collision, EMEventCuts::kNoSameBunchPileup)) { + return false; + } + if (mRequireVertexITSTPC && !IsSelected(collision, EMEventCuts::kIsVertexITSTPC)) { + return false; + } + if (mRequireGoodZvtxFT0vsPV && !IsSelected(collision, EMEventCuts::kIsGoodZvtxFT0vsPV)) { + return false; + } return true; } @@ -60,9 +76,11 @@ class EMEventCut : public TNamed bool IsSelected(T const& collision, const EMEventCuts& cut) const { switch (cut) { - case EMEventCuts::kFT0AND: + case EMEventCuts::kSel8: return collision.sel8(); - // return collision.selection_bit(o2::aod::evsel::kIsTriggerTVX); // alternative way. + + case EMEventCuts::kFT0AND: + return collision.selection_bit(o2::aod::evsel::kIsTriggerTVX); case EMEventCuts::kZvtx: return mMinZvtx < collision.posZ() && collision.posZ() < mMaxZvtx; @@ -73,25 +91,42 @@ class EMEventCut : public TNamed case EMEventCuts::kNoITSROFB: return collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder); + case EMEventCuts::kNoSameBunchPileup: + return collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup); + + case EMEventCuts::kIsVertexITSTPC: + return collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC); + + case EMEventCuts::kIsGoodZvtxFT0vsPV: + return collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV); + default: return false; } } // Setters + void SetRequireSel8(bool flag); void SetRequireFT0AND(bool flag); void SetZvtxRange(float min, float max); void SetRequireNoTFB(bool flag); void SetRequireNoITSROFB(bool flag); + void SetRequireNoSameBunchPileup(bool flag); + void SetRequireVertexITSTPC(bool flag); + void SetRequireIsGoodZvtxFT0vsPV(bool flag); /// @brief Print the track selection void print() const; private: + bool mRequireSel8{true}; bool mRequireFT0AND{true}; float mMinZvtx{-10.f}, mMaxZvtx{+10.f}; bool mRequireNoTFB{true}; bool mRequireNoITSROFB{true}; + bool mRequireNoSameBunchPileup{false}; + bool mRequireVertexITSTPC{false}; + bool mRequireGoodZvtxFT0vsPV{false}; ClassDef(EMEventCut, 1); }; diff --git a/PWGEM/PhotonMeson/Core/HistogramsLibrary.cxx b/PWGEM/PhotonMeson/Core/HistogramsLibrary.cxx index e413a235f3d..2f8f40e4088 100644 --- a/PWGEM/PhotonMeson/Core/HistogramsLibrary.cxx +++ b/PWGEM/PhotonMeson/Core/HistogramsLibrary.cxx @@ -299,7 +299,18 @@ void o2::aod::pwgem::photon::histogram::DefineHistograms(THashList* list, const list->Add(new TH2F("hMvsOPA_Pi0", "m_{ee} vs. opening angle;opening angle (rad.);m_{ee} (GeV/c^{2})", 500, 0, 0.5, 100, 0.0f, 0.1f)); // ee from pi0 dalitz decay list->Add(new TH2F("hMvsOPA_Eta", "m_{ee} vs. opening angle;opening angle (rad.);m_{ee} (GeV/c^{2})", 500, 0, 0.5, 100, 0.0f, 0.1f)); // ee from eta dalitz decay list->Add(new TH2F("hMvsOPA_Photon", "m_{ee} vs. opening angle;opening angle (rad.);m_{ee} (GeV/c^{2})", 500, 0, 0.5, 100, 0.0f, 0.1f)); // ee from photon conversion - } // end of mc + + static constexpr std::string_view parnames_LMEE[] = {"Pi0", "Eta", "EtaPrime", "Rho", "Omega", "Phi", "PromptJpsi", "NonPromptJpsi", "Ce_Ce", "Be_Be", "BCe_BCe", "BCe_Be_SameB", "BCe_Be_DiffB"}; + const int npar_lmee = sizeof(parnames_LMEE) / sizeof(parnames_LMEE[0]); + for (int i = 0; i < npar_lmee; i++) { + THnSparseF* hs_dilepton_mc_rec = new THnSparseF(Form("hs_dilepton_mc_rec_%s", parnames_LMEE[i].data()), Form("hs_dilepton_mc_rec from %s;m_{ee} (GeV/c^{2});p_{T,ee} (GeV/c);DCA_{ee}^{3D} (#sigma);#varphi_{V} (rad.);", parnames_LMEE[i].data()), ndim, nbins, xmin, xmax); + hs_dilepton_mc_rec->SetBinEdges(0, mee); + hs_dilepton_mc_rec->SetBinEdges(1, pt); + hs_dilepton_mc_rec->SetBinEdges(2, dca); + hs_dilepton_mc_rec->Sumw2(); + list->Add(hs_dilepton_mc_rec); + } + } // end of mc } else if (TString(histClass).Contains("MuMu")) { const int ndim = 4; // m, pt, dca, phiv const int nbins[ndim] = {90, 20, ndca - 1, 1}; @@ -597,6 +608,15 @@ void o2::aod::pwgem::photon::histogram::DefineHistograms(THashList* list, const TH2F* hMvsPt = new TH2F("hMvsPt", "m_{ee} vs. p_{T,ee};m_{ee} (GeV/c^{2});p_{T,ee} (GeV/c)", 400, 0, 4.0f, 1000, 0, 10.f); hMvsPt->Sumw2(); list->Add(hMvsPt); + + static constexpr std::string_view parnames_LMEE[] = {"Pi0", "Eta", "EtaPrime", "Rho", "Omega", "Phi", "PromptJpsi", "NonPromptJpsi", "Ce_Ce", "Be_Be", "BCe_BCe", "BCe_Be_SameB", "BCe_Be_DiffB"}; + const int npar_lmee = sizeof(parnames_LMEE) / sizeof(parnames_LMEE[0]); + for (int i = 0; i < npar_lmee; i++) { + TH2F* hMvsPt = new TH2F(Form("hMvsPt_%s", parnames_LMEE[i].data()), Form("m_{ee} vs. p_{T,ee} from MC %s;m_{ee} (GeV/c^{2});p_{T,ee} (GeV/c)", parnames_LMEE[i].data()), 400, 0, 4.0f, 1000, 0, 10.f); + hMvsPt->Sumw2(); + list->Add(hMvsPt); + } + } else if (TString(subGroup) == "dimuon") { TH2F* hMvsPt = new TH2F("hMvsPt", "m_{#mu#mu} vs. p_{T,#mu#mu};m_{#mu#mu} (GeV/c^{2});p_{T,#mu#mu} (GeV/c)", 90, 0.2, 1.1f, 200, 0, 2.f); hMvsPt->Sumw2(); diff --git a/PWGEM/PhotonMeson/TableProducer/associateMCinfo.cxx b/PWGEM/PhotonMeson/TableProducer/associateMCinfo.cxx index 707b29b16ec..8d4e4fea83a 100644 --- a/PWGEM/PhotonMeson/TableProducer/associateMCinfo.cxx +++ b/PWGEM/PhotonMeson/TableProducer/associateMCinfo.cxx @@ -216,6 +216,49 @@ struct AssociateMCInfo { fEventIdx[mctrack.globalIndex()] = fEventLabels.find(mcCollision.globalIndex())->second; fCounters[0]++; } + + bool is_hf_l = false; + if ((mctrack.isPhysicalPrimary() || mctrack.producedByGenerator()) && (abs(pdg) == 11 || abs(pdg) == 13) && mctrack.has_mothers()) { + auto mp = mctrack.template mothers_first_as(); // mother particle of electron + int pdg_mother = abs(mp.pdgCode()); + if (std::to_string(pdg_mother)[std::to_string(pdg_mother).length() - 3] == '4' || std::to_string(pdg_mother)[std::to_string(pdg_mother).length() - 3] == '5' || std::to_string(pdg_mother)[std::to_string(pdg_mother).length() - 4] == '4' || std::to_string(pdg_mother)[std::to_string(pdg_mother).length() - 4] == '5') { + is_hf_l = true; + } + } + + if (is_hf_l) { + // Next, store mother-chain for only HF->l, because HF->ll analysis requires correlation between HF hadrons or quarks. PLEASE DON'T do this for other partices. + int motherid = -999; // first mother index + if (mctrack.has_mothers()) { + motherid = mctrack.mothersIds()[0]; // first mother index + } + while (motherid > -1) { + if (motherid < mcTracks.size()) { // protect against bad mother indices. why is this needed? + auto mp = mcTracks.iteratorAt(motherid); + + if (abs(mp.pdgCode()) < 100) { // don't store quark/gluon informaiton, because data size explodes. + break; + } + + // if the MC truth particle corresponding to this reconstructed track which is not already written, add it to the skimmed MC stack + if (!(fNewLabels.find(mp.globalIndex()) != fNewLabels.end())) { + fNewLabels[mp.globalIndex()] = fCounters[0]; + fNewLabelsReversed[fCounters[0]] = mp.globalIndex(); + // fMCFlags[mp.globalIndex()] = mcflags; + fEventIdx[mp.globalIndex()] = fEventLabels.find(mcCollision.globalIndex())->second; + fCounters[0]++; + } + + if (mp.has_mothers()) { + motherid = mp.mothersIds()[0]; // first mother index + } else { + motherid = -999; + } + } else { + motherid = -999; + } + } // end of mother chain loop + } } // end of mc track loop if constexpr (static_cast(system & kPCM)) { diff --git a/PWGEM/PhotonMeson/Tasks/dalitzEEQCMC.cxx b/PWGEM/PhotonMeson/Tasks/dalitzEEQCMC.cxx index 043b21f9843..88a54ab702f 100644 --- a/PWGEM/PhotonMeson/Tasks/dalitzEEQCMC.cxx +++ b/PWGEM/PhotonMeson/Tasks/dalitzEEQCMC.cxx @@ -29,6 +29,7 @@ #include "PWGEM/PhotonMeson/Core/CutsLibrary.h" #include "PWGEM/PhotonMeson/Core/HistogramsLibrary.h" #include "PWGEM/PhotonMeson/Utils/MCUtilities.h" +#include "PWGEM/Dilepton/Utils/MCUtilities.h" using namespace o2; using namespace o2::aod; @@ -37,6 +38,7 @@ using namespace o2::framework::expressions; using namespace o2::soa; using namespace o2::aod::pwgem::mcutil; using namespace o2::aod::pwgem::photon; +using namespace o2::aod::pwgem::dilepton::mcutil; using std::array; using MyCollisions = soa::Join; @@ -148,20 +150,20 @@ struct DalitzEEQCMC { int FindLF(TTrack const& posmc, TTrack const& elemc, TMCTracks const& mcparticles) { int arr[] = { - FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 111, mcparticles), FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 221, mcparticles), FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 331, mcparticles), FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 113, mcparticles), FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 223, mcparticles), FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 333, mcparticles)}; + FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 111, mcparticles), FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 221, mcparticles), FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 331, mcparticles), FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 113, mcparticles), FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 223, mcparticles), FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 333, mcparticles), FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 443, mcparticles)}; int size = sizeof(arr) / sizeof(*arr); int max = *std::max_element(arr, arr + size); return max; } Partition uls_pairs = o2::aod::dalitzee::sign == 0; + Partition lspp_pairs = o2::aod::dalitzee::sign == +1; + Partition lsmm_pairs = o2::aod::dalitzee::sign == -1; SliceCache cache; Preslice perCollision = aod::dalitzee::emeventId; Partition grouped_collisions = (cfgCentMin < o2::aod::cent::centFT0M && o2::aod::cent::centFT0M < cfgCentMax) || (cfgCentMin < o2::aod::cent::centFT0A && o2::aod::cent::centFT0A < cfgCentMax) || (cfgCentMin < o2::aod::cent::centFT0C && o2::aod::cent::centFT0C < cfgCentMax); // this goes to same event. - std::vector used_trackIds; - void processQCMC(MyCollisions const&, MyDalitzEEs const&, MyMCTracks const&, aod::EMMCParticles const& mcparticles, aod::EMMCEvents const&) { THashList* list_ev_before = static_cast(fMainList->FindObject("Event")->FindObject(event_types[0].data())); @@ -187,11 +189,15 @@ struct DalitzEEQCMC { reinterpret_cast(list_ev_after->FindObject("hCollisionCounter"))->Fill("accepted", 1.f); auto uls_pairs_per_coll = uls_pairs->sliceByCached(o2::aod::dalitzee::emeventId, collision.globalIndex(), cache); + auto lspp_pairs_per_coll = lspp_pairs->sliceByCached(o2::aod::dalitzee::emeventId, collision.globalIndex(), cache); + auto lsmm_pairs_per_coll = lsmm_pairs->sliceByCached(o2::aod::dalitzee::emeventId, collision.globalIndex(), cache); + + std::vector used_trackIds; + used_trackIds.reserve(uls_pairs_per_coll.size() * 2); for (const auto& cut : fDalitzEECuts) { THashList* list_dalitzee_cut = static_cast(list_dalitzee->FindObject(cut.GetName())); THashList* list_track_cut = static_cast(list_track->FindObject(cut.GetName())); - used_trackIds.reserve(uls_pairs_per_coll.size() * 2); int nuls = 0; for (auto& uls_pair : uls_pairs_per_coll) { @@ -204,17 +210,29 @@ struct DalitzEEQCMC { auto posmc = pos.template emmcparticle_as(); auto elemc = ele.template emmcparticle_as(); + + if (posmc.pdgCode() != -11 || elemc.pdgCode() != 11) { + continue; + } + int mother_id = FindLF(posmc, elemc, mcparticles); + int hfee_type = IsHF(posmc, elemc, mcparticles); int photonid = FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 22, mcparticles); - if (mother_id < 0 && photonid < 0) { + if (mother_id < 0 && hfee_type < 0 && photonid < 0) { continue; } - if (mother_id > 0) { + + dca_pos_3d = pos.dca3DinSigma(); + dca_ele_3d = ele.dca3DinSigma(); + dca_ee_3d = std::sqrt((dca_pos_3d * dca_pos_3d + dca_ele_3d * dca_ele_3d) / 2.); + values[0] = uls_pair.mass(); + values[1] = uls_pair.pt(); + values[2] = dca_ee_3d; + values[3] = uls_pair.phiv(); + + if (mother_id > -1) { auto mcmother = mcparticles.iteratorAt(mother_id); if (mcmother.isPhysicalPrimary() || mcmother.producedByGenerator()) { - dca_pos_3d = pos.dca3DinSigma(); - dca_ele_3d = ele.dca3DinSigma(); - dca_ee_3d = std::sqrt((dca_pos_3d * dca_pos_3d + dca_ele_3d * dca_ele_3d) / 2.); if (cfgDoDCAstudy) { values_single[0] = uls_pair.mass(); @@ -224,10 +242,6 @@ struct DalitzEEQCMC { reinterpret_cast(list_dalitzee_cut->FindObject("hs_dilepton_uls_dca_same"))->Fill(values_single); } - values[0] = uls_pair.mass(); - values[1] = uls_pair.pt(); - values[2] = dca_ee_3d; - values[3] = uls_pair.phiv(); reinterpret_cast(list_dalitzee_cut->FindObject("hs_dilepton_uls_same"))->Fill(values); if (mcmother.pdgCode() == 111) { @@ -238,6 +252,37 @@ struct DalitzEEQCMC { reinterpret_cast(list_dalitzee_cut->FindObject("hMvsOPA_Eta"))->Fill(uls_pair.opangle(), uls_pair.mass()); } + switch (abs(mcmother.pdgCode())) { + case 111: + reinterpret_cast(list_dalitzee_cut->FindObject("hs_dilepton_mc_rec_Pi0"))->Fill(values); + break; + case 221: + reinterpret_cast(list_dalitzee_cut->FindObject("hs_dilepton_mc_rec_Eta"))->Fill(values); + break; + case 331: + reinterpret_cast(list_dalitzee_cut->FindObject("hs_dilepton_mc_rec_EtaPrime"))->Fill(values); + break; + case 113: + reinterpret_cast(list_dalitzee_cut->FindObject("hs_dilepton_mc_rec_Rho"))->Fill(values); + break; + case 223: + reinterpret_cast(list_dalitzee_cut->FindObject("hs_dilepton_mc_rec_Omega"))->Fill(values); + break; + case 333: + reinterpret_cast(list_dalitzee_cut->FindObject("hs_dilepton_mc_rec_Phi"))->Fill(values); + break; + case 443: { + if (IsFromBeauty(mcmother, mcparticles) > 0) { + reinterpret_cast(list_dalitzee_cut->FindObject("hs_dilepton_mc_rec_NonPromptJpsi"))->Fill(values); + } else { + reinterpret_cast(list_dalitzee_cut->FindObject("hs_dilepton_mc_rec_PromptJpsi"))->Fill(values); + } + break; + } + default: + break; + } + nuls++; for (auto& track : {pos, ele}) { if (std::find(used_trackIds.begin(), used_trackIds.end(), track.globalIndex()) == used_trackIds.end()) { @@ -250,7 +295,30 @@ struct DalitzEEQCMC { } } } // end of LF - } else if (photonid > 0) { + } else if (hfee_type > -1) { + if ((posmc.isPhysicalPrimary() || posmc.producedByGenerator()) && (elemc.isPhysicalPrimary() || elemc.producedByGenerator())) { + switch (hfee_type) { + case static_cast(EM_HFeeType::kCe_Ce): + reinterpret_cast(list_dalitzee_cut->FindObject("hs_dilepton_mc_rec_Ce_Ce"))->Fill(values); + break; + case static_cast(EM_HFeeType::kBe_Be): + reinterpret_cast(list_dalitzee_cut->FindObject("hs_dilepton_mc_rec_Be_Be"))->Fill(values); + break; + case static_cast(EM_HFeeType::kBCe_BCe): + reinterpret_cast(list_dalitzee_cut->FindObject("hs_dilepton_mc_rec_BCe_BCe"))->Fill(values); + break; + case static_cast(EM_HFeeType::kBCe_Be_SameB): // ULS + reinterpret_cast(list_dalitzee_cut->FindObject("hs_dilepton_mc_rec_BCe_Be_SameB"))->Fill(values); + break; + case static_cast(EM_HFeeType::kBCe_Be_DiffB): // LS + LOGF(info, "You should not see kBCe_Be_DiffB in ULS. Good luck."); + break; + default: + break; + } + nuls++; + } + } else if (photonid > -1) { auto mcphoton = mcparticles.iteratorAt(photonid); if ((mcphoton.isPhysicalPrimary() || mcphoton.producedByGenerator()) && IsEleFromPC(elemc, mcparticles) && IsEleFromPC(posmc, mcparticles)) { reinterpret_cast(list_dalitzee_cut->FindObject("hMvsPhiV_Photon"))->Fill(uls_pair.phiv(), uls_pair.mass()); @@ -260,6 +328,134 @@ struct DalitzEEQCMC { } // end of uls pair loop reinterpret_cast(list_dalitzee_cut->FindObject("hNpair_uls"))->Fill(nuls); + int nlspp = 0; + for (auto& lspp_pair : lspp_pairs_per_coll) { + auto pos = lspp_pair.template posTrack_as(); + auto ele = lspp_pair.template negTrack_as(); + + if (!cut.IsSelected(lspp_pair)) { + continue; + } + + auto posmc = pos.template emmcparticle_as(); + auto elemc = ele.template emmcparticle_as(); + + if (posmc.pdgCode() != -11 || elemc.pdgCode() != -11) { + continue; + } + + int hfee_type = IsHF(posmc, elemc, mcparticles); + if (hfee_type < 0) { + continue; + } + + dca_pos_3d = pos.dca3DinSigma(); + dca_ele_3d = ele.dca3DinSigma(); + dca_ee_3d = std::sqrt((dca_pos_3d * dca_pos_3d + dca_ele_3d * dca_ele_3d) / 2.); + values[0] = lspp_pair.mass(); + values[1] = lspp_pair.pt(); + values[2] = dca_ee_3d; + values[3] = lspp_pair.phiv(); + + if ((posmc.isPhysicalPrimary() || posmc.producedByGenerator()) && (elemc.isPhysicalPrimary() || elemc.producedByGenerator())) { + switch (hfee_type) { + case static_cast(EM_HFeeType::kCe_Ce): + LOGF(info, "You should not see kCe_Ce in LS++. Good luck."); + break; + case static_cast(EM_HFeeType::kBe_Be): + LOGF(info, "You should not see kBe_Be in LS++. Good luck."); + break; + case static_cast(EM_HFeeType::kBCe_BCe): + LOGF(info, "You should not see kBCe_BCe in LS++. Good luck."); + break; + case static_cast(EM_HFeeType::kBCe_Be_SameB): // ULS + LOGF(info, "You should not see kBCe_Be_SameB in LS++. Good luck."); + break; + case static_cast(EM_HFeeType::kBCe_Be_DiffB): // LS + reinterpret_cast(list_dalitzee_cut->FindObject("hs_dilepton_mc_rec_BCe_Be_DiffB"))->Fill(values); + break; + default: + break; + } + } + nlspp++; + for (auto& track : {pos, ele}) { + if (std::find(used_trackIds.begin(), used_trackIds.end(), track.globalIndex()) == used_trackIds.end()) { + o2::aod::pwgem::photon::histogram::FillHistClass(list_track_cut, "", track); + used_trackIds.emplace_back(track.globalIndex()); + auto mctrack = track.template emmcparticle_as(); + reinterpret_cast(fMainList->FindObject("Track")->FindObject(cut.GetName())->FindObject("hPtGen_DeltaPtOverPtGen"))->Fill(mctrack.pt(), (track.pt() - mctrack.pt()) / mctrack.pt()); + reinterpret_cast(fMainList->FindObject("Track")->FindObject(cut.GetName())->FindObject("hPtGen_DeltaEta"))->Fill(mctrack.pt(), track.eta() - mctrack.eta()); + reinterpret_cast(fMainList->FindObject("Track")->FindObject(cut.GetName())->FindObject("hPtGen_DeltaPhi"))->Fill(mctrack.pt(), track.phi() - mctrack.phi()); + } + } + } // end of lspp pair loop + reinterpret_cast(list_dalitzee_cut->FindObject("hNpair_lspp"))->Fill(nlspp); + + int nlsmm = 0; + for (auto& lsmm_pair : lsmm_pairs_per_coll) { + auto pos = lsmm_pair.template posTrack_as(); + auto ele = lsmm_pair.template negTrack_as(); + + if (!cut.IsSelected(lsmm_pair)) { + continue; + } + + auto posmc = pos.template emmcparticle_as(); + auto elemc = ele.template emmcparticle_as(); + + if (posmc.pdgCode() != 11 || elemc.pdgCode() != 11) { + continue; + } + + int hfee_type = IsHF(posmc, elemc, mcparticles); + if (hfee_type < 0) { + continue; + } + + dca_pos_3d = pos.dca3DinSigma(); + dca_ele_3d = ele.dca3DinSigma(); + dca_ee_3d = std::sqrt((dca_pos_3d * dca_pos_3d + dca_ele_3d * dca_ele_3d) / 2.); + values[0] = lsmm_pair.mass(); + values[1] = lsmm_pair.pt(); + values[2] = dca_ee_3d; + values[3] = lsmm_pair.phiv(); + + if ((posmc.isPhysicalPrimary() || posmc.producedByGenerator()) && (elemc.isPhysicalPrimary() || elemc.producedByGenerator())) { + switch (hfee_type) { + case static_cast(EM_HFeeType::kCe_Ce): + LOGF(info, "You should not see kCe_Ce in LS--. Good luck."); + break; + case static_cast(EM_HFeeType::kBe_Be): + LOGF(info, "You should not see kBe_Be in LS--. Good luck."); + break; + case static_cast(EM_HFeeType::kBCe_BCe): + LOGF(info, "You should not see kBCe_BCe in LS--. Good luck."); + break; + case static_cast(EM_HFeeType::kBCe_Be_SameB): // ULS + LOGF(info, "You should not see kBCe_Be_SameB in LS--. Good luck."); + break; + case static_cast(EM_HFeeType::kBCe_Be_DiffB): // LS + reinterpret_cast(list_dalitzee_cut->FindObject("hs_dilepton_mc_rec_BCe_Be_DiffB"))->Fill(values); + break; + default: + break; + } + } + nlsmm++; + for (auto& track : {pos, ele}) { + if (std::find(used_trackIds.begin(), used_trackIds.end(), track.globalIndex()) == used_trackIds.end()) { + o2::aod::pwgem::photon::histogram::FillHistClass(list_track_cut, "", track); + used_trackIds.emplace_back(track.globalIndex()); + auto mctrack = track.template emmcparticle_as(); + reinterpret_cast(fMainList->FindObject("Track")->FindObject(cut.GetName())->FindObject("hPtGen_DeltaPtOverPtGen"))->Fill(mctrack.pt(), (track.pt() - mctrack.pt()) / mctrack.pt()); + reinterpret_cast(fMainList->FindObject("Track")->FindObject(cut.GetName())->FindObject("hPtGen_DeltaEta"))->Fill(mctrack.pt(), track.eta() - mctrack.eta()); + reinterpret_cast(fMainList->FindObject("Track")->FindObject(cut.GetName())->FindObject("hPtGen_DeltaPhi"))->Fill(mctrack.pt(), track.phi() - mctrack.phi()); + } + } + } // end of lsmm pair loop + reinterpret_cast(list_dalitzee_cut->FindObject("hNpair_lsmm"))->Fill(nlsmm); + used_trackIds.clear(); used_trackIds.shrink_to_fit(); } // end of cut loop @@ -288,19 +484,6 @@ struct DalitzEEQCMC { reinterpret_cast(fMainList->FindObject("Generated")->FindObject("hCollisionCounter"))->Fill(1.0); reinterpret_cast(fMainList->FindObject("Generated")->FindObject("hZvtx_before"))->Fill(mccollision.posZ()); - if (!collision.sel8()) { - continue; - } - reinterpret_cast(fMainList->FindObject("Generated")->FindObject("hCollisionCounter"))->Fill(2.0); - - if (collision.numContrib() < 0.5) { - continue; - } - reinterpret_cast(fMainList->FindObject("Generated")->FindObject("hCollisionCounter"))->Fill(3.0); - - if (abs(collision.posZ()) > 10.0) { - continue; - } if (!fEMEventCut.IsSelected(collision)) { continue; @@ -320,19 +503,176 @@ struct DalitzEEQCMC { continue; } + if (!t1.isPhysicalPrimary() && !t1.producedByGenerator()) { + continue; + } + if (!t2.isPhysicalPrimary() && !t2.producedByGenerator()) { + continue; + } + int mother_id = FindLF(t1, t2, mcparticles); - if (mother_id < 0) { + int hfee_type = IsHF(t1, t2, mcparticles); + if (mother_id < 0 && hfee_type < 0) { + continue; + } + ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), o2::constants::physics::MassElectron); + ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), o2::constants::physics::MassElectron); + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + + if (mother_id > -1) { + auto mcmother = mcparticles.iteratorAt(mother_id); + if (mcmother.isPhysicalPrimary() || mcmother.producedByGenerator()) { + reinterpret_cast(fMainList->FindObject("Generated")->FindObject("hMvsPt"))->Fill(v12.M(), v12.Pt()); + + switch (abs(mcmother.pdgCode())) { + case 111: + reinterpret_cast(fMainList->FindObject("Generated")->FindObject("hMvsPt_Pi0"))->Fill(v12.M(), v12.Pt()); + break; + case 221: + reinterpret_cast(fMainList->FindObject("Generated")->FindObject("hMvsPt_Eta"))->Fill(v12.M(), v12.Pt()); + break; + case 331: + reinterpret_cast(fMainList->FindObject("Generated")->FindObject("hMvsPt_EtaPrime"))->Fill(v12.M(), v12.Pt()); + break; + case 113: + reinterpret_cast(fMainList->FindObject("Generated")->FindObject("hMvsPt_Rho"))->Fill(v12.M(), v12.Pt()); + break; + case 223: + reinterpret_cast(fMainList->FindObject("Generated")->FindObject("hMvsPt_Omega"))->Fill(v12.M(), v12.Pt()); + break; + case 333: + reinterpret_cast(fMainList->FindObject("Generated")->FindObject("hMvsPt_Phi"))->Fill(v12.M(), v12.Pt()); + break; + case 443: { + if (IsFromBeauty(mcmother, mcparticles) > 0) { + reinterpret_cast(fMainList->FindObject("Generated")->FindObject("hMvsPt_NonPromptJpsi"))->Fill(v12.M(), v12.Pt()); + } else { + reinterpret_cast(fMainList->FindObject("Generated")->FindObject("hMvsPt_PromptJpsi"))->Fill(v12.M(), v12.Pt()); + } + break; + } + default: + break; + } + } + } else if (hfee_type > -1) { + switch (hfee_type) { + case static_cast(EM_HFeeType::kCe_Ce): + reinterpret_cast(fMainList->FindObject("Generated")->FindObject("hMvsPt_Ce_Ce"))->Fill(v12.M(), v12.Pt()); + break; + case static_cast(EM_HFeeType::kBe_Be): + reinterpret_cast(fMainList->FindObject("Generated")->FindObject("hMvsPt_Be_Be"))->Fill(v12.M(), v12.Pt()); + break; + case static_cast(EM_HFeeType::kBCe_BCe): + reinterpret_cast(fMainList->FindObject("Generated")->FindObject("hMvsPt_BCe_BCe"))->Fill(v12.M(), v12.Pt()); + break; + case static_cast(EM_HFeeType::kBCe_Be_SameB): // ULS + reinterpret_cast(fMainList->FindObject("Generated")->FindObject("hMvsPt_BCe_Be_SameB"))->Fill(v12.M(), v12.Pt()); + break; + case static_cast(EM_HFeeType::kBCe_Be_DiffB): // LS + LOGF(info, "You should not see kBCe_Be_DiffB in ULS. Good luck."); + break; + default: + break; + } + } + } // end of true ULS pair loop + + for (auto& [t1, t2] : combinations(CombinationsStrictlyUpperIndexPolicy(posTracks_per_coll, posTracks_per_coll))) { + // LOGF(info, "pdg1 = %d, pdg2 = %d", t1.pdgCode(), t2.pdgCode()); + + if (t1.pt() < min_mcPt || max_mcPt < t1.pt() || abs(t1.eta()) > max_mcEta) { + continue; + } + if (t2.pt() < min_mcPt || max_mcPt < t2.pt() || abs(t2.eta()) > max_mcEta) { + continue; + } + + if (!t1.isPhysicalPrimary() && !t1.producedByGenerator()) { continue; } - auto mcmother = mcparticles.iteratorAt(mother_id); - if (mcmother.isPhysicalPrimary() || mcmother.producedByGenerator()) { - ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), o2::constants::physics::MassElectron); - ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), o2::constants::physics::MassElectron); - ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; - reinterpret_cast(fMainList->FindObject("Generated")->FindObject("hMvsPt"))->Fill(v12.M(), v12.Pt()); - } // end of LF - } // end of true ULS pair loop - } // end of collision loop + if (!t2.isPhysicalPrimary() && !t2.producedByGenerator()) { + continue; + } + + int hfee_type = IsHF(t1, t2, mcparticles); + if (hfee_type < 0) { + continue; + } + ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), o2::constants::physics::MassElectron); + ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), o2::constants::physics::MassElectron); + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + if (hfee_type > -1) { + switch (hfee_type) { + case static_cast(EM_HFeeType::kCe_Ce): + LOGF(info, "You should not see kCe_Ce in LS++. Good luck."); + break; + case static_cast(EM_HFeeType::kBe_Be): + LOGF(info, "You should not see kBe_Be in LS++. Good luck."); + break; + case static_cast(EM_HFeeType::kBCe_BCe): + LOGF(info, "You should not see kBCe_BCe in LS++. Good luck."); + break; + case static_cast(EM_HFeeType::kBCe_Be_SameB): // ULS + LOGF(info, "You should not see kBCe_Be_SameB in LS++. Good luck."); + break; + case static_cast(EM_HFeeType::kBCe_Be_DiffB): // LS + reinterpret_cast(fMainList->FindObject("Generated")->FindObject("hMvsPt_BCe_Be_DiffB"))->Fill(v12.M(), v12.Pt()); + break; + default: + break; + } + } + } // end of true LS++ pair loop + + for (auto& [t1, t2] : combinations(CombinationsStrictlyUpperIndexPolicy(negTracks_per_coll, negTracks_per_coll))) { + // LOGF(info, "pdg1 = %d, pdg2 = %d", t1.pdgCode(), t2.pdgCode()); + + if (t1.pt() < min_mcPt || max_mcPt < t1.pt() || abs(t1.eta()) > max_mcEta) { + continue; + } + if (t2.pt() < min_mcPt || max_mcPt < t2.pt() || abs(t2.eta()) > max_mcEta) { + continue; + } + + if (!t1.isPhysicalPrimary() && !t1.producedByGenerator()) { + continue; + } + if (!t2.isPhysicalPrimary() && !t2.producedByGenerator()) { + continue; + } + + int hfee_type = IsHF(t1, t2, mcparticles); + if (hfee_type < 0) { + continue; + } + ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), o2::constants::physics::MassElectron); + ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), o2::constants::physics::MassElectron); + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + if (hfee_type > -1) { + switch (hfee_type) { + case static_cast(EM_HFeeType::kCe_Ce): + LOGF(info, "You should not see kCe_Ce in LS--. Good luck."); + break; + case static_cast(EM_HFeeType::kBe_Be): + LOGF(info, "You should not see kBe_Be in LS--. Good luck."); + break; + case static_cast(EM_HFeeType::kBCe_BCe): + LOGF(info, "You should not see kBCe_BCe in LS--. Good luck."); + break; + case static_cast(EM_HFeeType::kBCe_Be_SameB): // ULS + LOGF(info, "You should not see kBCe_Be_SameB in LS--. Good luck."); + break; + case static_cast(EM_HFeeType::kBCe_Be_DiffB): // LS + reinterpret_cast(fMainList->FindObject("Generated")->FindObject("hMvsPt_BCe_Be_DiffB"))->Fill(v12.M(), v12.Pt()); + break; + default: + break; + } + } + } // end of true LS++ pair loop + + } // end of collision loop } PROCESS_SWITCH(DalitzEEQCMC, processGen, "run genrated info", true);