From f7a786e2248917ffbf2978a6ecb37d11218efc5a Mon Sep 17 00:00:00 2001 From: Daiki Sekihata Date: Wed, 5 Nov 2025 16:31:45 +0100 Subject: [PATCH] PWGEM/Dilepton: update HFll in MC with status code --- PWGEM/Dilepton/Core/SingleTrackQCMC.h | 36 ++---- PWGEM/Dilepton/Utils/MCUtilities.h | 172 +++++++++++++++++++------- 2 files changed, 139 insertions(+), 69 deletions(-) diff --git a/PWGEM/Dilepton/Core/SingleTrackQCMC.h b/PWGEM/Dilepton/Core/SingleTrackQCMC.h index 5f5fce49270..9e2cdbf8847 100644 --- a/PWGEM/Dilepton/Core/SingleTrackQCMC.h +++ b/PWGEM/Dilepton/Core/SingleTrackQCMC.h @@ -317,10 +317,6 @@ struct SingleTrackQCMC { fRegistry.add("Track/PID/positive/hTPCNsigmaPr", "TPC n sigma pr;p_{in} (GeV/c);n #sigma_{p}^{TPC}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); fRegistry.add("Track/PID/positive/hTOFbeta", "TOF #beta;p_{pv} (GeV/c);#beta", kTH2F, {{1000, 0, 10}, {240, 0, 1.2}}, false); fRegistry.add("Track/PID/positive/hTOFNsigmaEl", "TOF n sigma el;p_{pv} (GeV/c);n #sigma_{e}^{TOF}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); - // fRegistry.add("Track/PID/positive/hTOFNsigmaMu", "TOF n sigma mu;p_{pv} (GeV/c);n #sigma_{#mu}^{TOF}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); - // fRegistry.add("Track/PID/positive/hTOFNsigmaPi", "TOF n sigma pi;p_{pv} (GeV/c);n #sigma_{#pi}^{TOF}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); - // fRegistry.add("Track/PID/positive/hTOFNsigmaKa", "TOF n sigma ka;p_{pv} (GeV/c);n #sigma_{K}^{TOF}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); - // fRegistry.add("Track/PID/positive/hTOFNsigmaPr", "TOF n sigma pr;p_{pv} (GeV/c);n #sigma_{p}^{TOF}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); fRegistry.add("Track/PID/positive/hMeanClusterSizeITS", "mean cluster size ITS;p_{pv} (GeV/c); on ITS #times cos(#lambda)", kTH2F, {{1000, 0.f, 10.f}, {150, 0, 15}}, false); fRegistry.add("Track/PID/positive/hMeanClusterSizeITSib", "mean cluster size ITS inner barrel;p_{pv} (GeV/c); on ITS #times cos(#lambda)", kTH2F, {{1000, 0.f, 10.f}, {150, 0, 15}}, false); fRegistry.add("Track/PID/positive/hMeanClusterSizeITSob", "mean cluster size ITS outer barrel;p_{pv} (GeV/c); on ITS #times cos(#lambda)", kTH2F, {{1000, 0.f, 10.f}, {150, 0, 15}}, false); @@ -693,10 +689,6 @@ struct SingleTrackQCMC { fRegistry.fill(HIST("Track/PID/positive/hTPCNsigmaKa"), track.tpcInnerParam(), track.tpcNSigmaKa()); fRegistry.fill(HIST("Track/PID/positive/hTPCNsigmaPr"), track.tpcInnerParam(), track.tpcNSigmaPr()); fRegistry.fill(HIST("Track/PID/positive/hTOFNsigmaEl"), track.p(), track.tofNSigmaEl()); - // fRegistry.fill(HIST("Track/PID/positive/hTOFNsigmaMu"), track.p(), track.tofNSigmaMu()); - // fRegistry.fill(HIST("Track/PID/positive/hTOFNsigmaPi"), track.p(), track.tofNSigmaPi()); - // fRegistry.fill(HIST("Track/PID/positive/hTOFNsigmaKa"), track.p(), track.tofNSigmaKa()); - // fRegistry.fill(HIST("Track/PID/positive/hTOFNsigmaPr"), track.p(), track.tofNSigmaPr()); } } else { fRegistry.fill(HIST("Track/") + HIST(lepton_source_types[lepton_source_id]) + HIST("negative/hs"), track.pt(), track.eta(), track.phi(), dca3D, dcaXY, dcaZ, -mctrack.pdgCode() / pdg_lepton, weight); @@ -742,10 +734,6 @@ struct SingleTrackQCMC { fRegistry.fill(HIST("Track/PID/negative/hTPCNsigmaKa"), track.tpcInnerParam(), track.tpcNSigmaKa()); fRegistry.fill(HIST("Track/PID/negative/hTPCNsigmaPr"), track.tpcInnerParam(), track.tpcNSigmaPr()); fRegistry.fill(HIST("Track/PID/negative/hTOFNsigmaEl"), track.p(), track.tofNSigmaEl()); - // fRegistry.fill(HIST("Track/PID/negative/hTOFNsigmaMu"), track.p(), track.tofNSigmaMu()); - // fRegistry.fill(HIST("Track/PID/negative/hTOFNsigmaPi"), track.p(), track.tofNSigmaPi()); - // fRegistry.fill(HIST("Track/PID/negative/hTOFNsigmaKa"), track.p(), track.tofNSigmaKa()); - // fRegistry.fill(HIST("Track/PID/negative/hTOFNsigmaPr"), track.p(), track.tofNSigmaPr()); } } } @@ -908,14 +896,14 @@ struct SingleTrackQCMC { } else { fillTrackInfo<5, TMCParticles>(track); } - } else if (IsFromBeauty(mctrack, mcparticles) > 0) { // b is found in full decay chain. - if (IsFromCharm(mctrack, mcparticles) > 0) { // c is found in full decay chain. - fillTrackInfo<9, TMCParticles>(track); + } else if (isWeakDecayFromBeautyHadron(mctrack, mcparticles)) { // hb->l is found in full decay chain. + fillTrackInfo<8, TMCParticles>(track); + } else if (isWeakDecayFromCharmHadron(mctrack, mcparticles)) { // hc->l is found in full decay chain. + if (IsFromBeauty(mcmother, mcparticles) > 0) { + fillTrackInfo<9, TMCParticles>(track); // hb->hc->l is fond. } else { - fillTrackInfo<8, TMCParticles>(track); + fillTrackInfo<7, TMCParticles>(track); // prompt hc->l is found. } - } else if (IsFromCharm(mctrack, mcparticles) > 0) { // c is found in full decay chain. Not from b. - fillTrackInfo<7, TMCParticles>(track); } } else { if (pdg_mother == 22) { // photon conversion @@ -1013,14 +1001,14 @@ struct SingleTrackQCMC { } else { fRegistry.fill(HIST("Generated/PromptPsi2S/hs"), pt, eta, phi, -lepton.pdgCode() / pdg_lepton); } - } else if (IsFromBeauty(lepton, mcparticles) > 0) { // b is found in full decay chain. - if (IsFromCharm(lepton, mcparticles) > 0) { // c is found in full decay chain. - fRegistry.fill(HIST("Generated/b2c2l/hs"), pt, eta, phi, -lepton.pdgCode() / pdg_lepton); + } else if (isWeakDecayFromBeautyHadron(lepton, mcparticles)) { // hb->l is found + fRegistry.fill(HIST("Generated/b2l/hs"), pt, eta, phi, -lepton.pdgCode() / pdg_lepton); + } else if (isWeakDecayFromCharmHadron(lepton, mcparticles)) { // hc->l is found in full decay chain. + if (IsFromBeauty(mcmother, mcparticles) > 0) { + fRegistry.fill(HIST("Generated/b2c2l/hs"), pt, eta, phi, -lepton.pdgCode() / pdg_lepton); // hb->hc->l is found in full decay chain. } else { - fRegistry.fill(HIST("Generated/b2l/hs"), pt, eta, phi, -lepton.pdgCode() / pdg_lepton); + fRegistry.fill(HIST("Generated/c2l/hs"), pt, eta, phi, -lepton.pdgCode() / pdg_lepton); // prompt hc->l is found in full decay chain. } - } else if (IsFromCharm(lepton, mcparticles) > 0) { // c is found in full decay chain. Not from b. - fRegistry.fill(HIST("Generated/c2l/hs"), pt, eta, phi, -lepton.pdgCode() / pdg_lepton); } } // end of mc lepton loop per collision diff --git a/PWGEM/Dilepton/Utils/MCUtilities.h b/PWGEM/Dilepton/Utils/MCUtilities.h index fdbe349ff62..c5928c0560b 100644 --- a/PWGEM/Dilepton/Utils/MCUtilities.h +++ b/PWGEM/Dilepton/Utils/MCUtilities.h @@ -71,6 +71,127 @@ int hasFakeMatchMFTMCH(TTrack const& track) } } //_______________________________________________________________________ +template +bool isCharmonia(T const& track) +{ + if (std::abs(track.pdgCode()) < 100) { + return false; + } + + std::string pdgStr = std::to_string(std::abs(track.pdgCode())); + int n = pdgStr.length(); + int pdg3 = std::stoi(pdgStr.substr(n - 3, 3)); + + if (pdg3 == 441 || pdg3 == 443 || pdg3 == 445 || pdg3 == 447) { + return true; + } else { + return false; + } +} +//_______________________________________________________________________ +template +bool isCharmMeson(T const& track) +{ + if (isCharmonia(track)) { + return false; + } + + if (400 < std::abs(track.pdgCode()) && std::abs(track.pdgCode()) < 500) { + return true; + } else { + return false; + } +} +//_______________________________________________________________________ +template +bool isCharmBaryon(T const& track) +{ + if (4000 < std::abs(track.pdgCode()) && std::abs(track.pdgCode()) < 5000) { + return true; + } else { + return false; + } +} +//_______________________________________________________________________ +template +bool isBottomonia(T const& track) +{ + if (std::abs(track.pdgCode()) < 100) { + return false; + } + + std::string pdgStr = std::to_string(std::abs(track.pdgCode())); + int n = pdgStr.length(); + int pdg3 = std::stoi(pdgStr.substr(n - 3, 3)); + + if (pdg3 == 551 || pdg3 == 553 || pdg3 == 555 || pdg3 == 557) { + return true; + } else { + return false; + } +} +//_______________________________________________________________________ +template +bool isBeautyMeson(T const& track) +{ + if (isBottomonia(track)) { + return false; + } + + if (500 < std::abs(track.pdgCode()) && std::abs(track.pdgCode()) < 600) { + return true; + } else { + return false; + } +} +//_______________________________________________________________________ +template +bool isBeautyBaryon(T const& track) +{ + if (5000 < std::abs(track.pdgCode()) && std::abs(track.pdgCode()) < 6000) { + return true; + } else { + return false; + } +} +//_______________________________________________________________________ +template +bool isWeakDecayFromCharmHadron(T const& mcParticle, U const& mcParticles) +{ + // require that the direct mother is charm hadron via semileptonic. e.g. hc->e, not hc->X->pi0->eegamma + if (!mcParticle.has_mothers()) { + return false; + } + if (mcParticle.getProcess() != 4) { // weak decay + return false; + } + auto mp = mcParticles.iteratorAt(mcParticle.mothersIds()[0]); + if (isCharmMeson(mp) || isCharmBaryon(mp)) { + return true; + } else { + return false; + } +} +//_______________________________________________________________________ +template +bool isWeakDecayFromBeautyHadron(T const& mcParticle, U const& mcParticles) +{ + // require that the direct mother is beauty hadron via semileptonice decay. e.g. hb->e, not hb->X->pi0->eegamma + if (!mcParticle.has_mothers()) { + return false; + } + if (mcParticle.getProcess() != 4) { // weak decay + return false; + } + auto mp = mcParticles.iteratorAt(mcParticle.mothersIds()[0]); + if (isBeautyMeson(mp) || isBeautyBaryon(mp)) { + return true; + } else { + return false; + } +} +//_______________________________________________________________________ +//_______________________________________________________________________ template int FindCommonMotherFrom2ProngsWithoutPDG(TMCParticle1 const& p1, TMCParticle2 const& p2) { @@ -436,12 +557,12 @@ int IsHF(TMCParticle1 const& p1, TMCParticle2 const& p2, TMCParticles const& mcp bool isFOFound1 = findFlavorOscillationB(p1, mcparticles); bool isFOFound2 = findFlavorOscillationB(p2, mcparticles); - bool is_direct_from_b1 = IsFromBeauty(p1, mcparticles) > 0 && IsFromCharm(p1, mcparticles) < 0; - bool is_direct_from_b2 = IsFromBeauty(p2, mcparticles) > 0 && IsFromCharm(p2, mcparticles) < 0; - bool is_prompt_c1 = IsFromBeauty(p1, mcparticles) < 0 && IsFromCharm(p1, mcparticles) > 0; - bool is_prompt_c2 = IsFromBeauty(p2, mcparticles) < 0 && IsFromCharm(p2, mcparticles) > 0; - bool is_c_from_b1 = IsFromBeauty(p1, mcparticles) > 0 && IsFromCharm(p1, mcparticles) > 0; - bool is_c_from_b2 = IsFromBeauty(p2, mcparticles) > 0 && IsFromCharm(p2, mcparticles) > 0; + bool is_direct_from_b1 = isWeakDecayFromBeautyHadron(p1, mcparticles); + bool is_direct_from_b2 = isWeakDecayFromBeautyHadron(p2, mcparticles); + bool is_prompt_c1 = isWeakDecayFromCharmHadron(p1, mcparticles) && IsFromBeauty(p1, mcparticles) < 0; + bool is_prompt_c2 = isWeakDecayFromCharmHadron(p2, mcparticles) && IsFromBeauty(p2, mcparticles) < 0; + bool is_c_from_b1 = isWeakDecayFromCharmHadron(p1, mcparticles) && IsFromBeauty(p1, mcparticles) > 0; + bool is_c_from_b2 = isWeakDecayFromCharmHadron(p2, mcparticles) && IsFromBeauty(p2, mcparticles) > 0; if (is_prompt_c1 && is_prompt_c2 && mpfh1.pdgCode() * mpfh2.pdgCode() < 0) { // charmed mesons never oscillate. only ULS mothers_id1.clear(); @@ -640,45 +761,6 @@ bool checkFromSameQuarkPair(T& p1, T& p2, U& mcParticles, int pdg) return id1 == id2 && id1 > -1 && id2 > -1; } //_______________________________________________________________________ -template -bool isCharmMeson(T const& track) -{ - if (400 < std::abs(track.pdgCode()) && std::abs(track.pdgCode()) < 500) { - return true; - } else { - return false; - } -} -//_______________________________________________________________________ -template -bool isCharmBaryon(T const& track) -{ - if (4000 < std::abs(track.pdgCode()) && std::abs(track.pdgCode()) < 5000) { - return true; - } else { - return false; - } -} -//_______________________________________________________________________ -template -bool isBeautyMeson(T const& track) -{ - if (500 < std::abs(track.pdgCode()) && std::abs(track.pdgCode()) < 600) { - return true; - } else { - return false; - } -} -//_______________________________________________________________________ -template -bool isBeautyBaryon(T const& track) -{ - if (5000 < std::abs(track.pdgCode()) && std::abs(track.pdgCode()) < 6000) { - return true; - } else { - return false; - } -} //_______________________________________________________________________ //_______________________________________________________________________ } // namespace o2::aod::pwgem::dilepton::utils::mcutil