From 75bb2a016e4d534ba6bcaefa1552b6c7195fce3f Mon Sep 17 00:00:00 2001 From: jimun_lee Date: Thu, 9 Oct 2025 13:06:41 +0900 Subject: [PATCH] [PWGLF] Fixed reconstruction part of KstarInOO.cxx --- PWGLF/Tasks/Resonances/kstarInOO.cxx | 186 ++++++++++++++++++++------- 1 file changed, 141 insertions(+), 45 deletions(-) diff --git a/PWGLF/Tasks/Resonances/kstarInOO.cxx b/PWGLF/Tasks/Resonances/kstarInOO.cxx index 24ba298e3c9..bd5286986f5 100644 --- a/PWGLF/Tasks/Resonances/kstarInOO.cxx +++ b/PWGLF/Tasks/Resonances/kstarInOO.cxx @@ -179,21 +179,25 @@ struct kstarInOO { } if (cfgMcHistos) { + histos.add("hPion_PID_Purity", "hPion_PID_Purity", kTH1F, {{3, -1.5, 1.5}}); + histos.add("hKaon_PID_Purity", "hKaon_PID_Purity", kTH1F, {{3, -1.5, 1.5}}); + histos.add("hSimplePion_PID_Purity", "hSimplePion_PID_Purity", kTH1F, {{3, -1.5, 1.5}}); + histos.add("hSimpleKaon_PID_Purity", "hSimpleKaon_PID_Purity", kTH1F, {{3, -1.5, 1.5}}); + histos.add("nEvents_MC", "nEvents_MC", kTH1F, {{4, 0.0, 4.0}}); histos.add("nEvents_MC_True", "nEvents_MC_True", kTH1F, {{4, 0.0, 4.0}}); - histos.add("hMC_USS", "hMC_USS", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); - histos.add("hMC_LSS", "hMC_LSS", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); - histos.add("hMC_USS_Mix", "hMC_USS_Mix", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); - histos.add("hMC_LSS_Mix", "hMC_LSS_Mix", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); - histos.add("hMC_USS_True", "hMC_USS_True", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); histos.add("hMC_kstar_True", "hMC_kstar_True", kTHnSparseF, {cfgCentAxis, ptAxis}); - histos.add("hMC_USS_pion", "hMC_USS_pion", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); - histos.add("hMC_LSS_pion", "hMC_LSS_pion", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); - histos.add("hMC_USS_Mix_pion", "hMC_USS_Mix_pion", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); - histos.add("hMC_LSS_Mix_pion", "hMC_LSS_Mix_pion", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); - histos.add("hMC_USS_pion_True", "hMC_USS_pion_True", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); + histos.add("hMC_USS_True", "hMC_USS_True", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); + histos.add("hMC_USS_KPi", "hMC_USS_KPi", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); + histos.add("hMC_USS_PiK", "hMC_USS_PiK", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); + histos.add("hMC_LSS_KPi", "hMC_LSS_KPi", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); + histos.add("hMC_LSS_PiK", "hMC_LSS_PiK", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); + histos.add("hMC_USS_KPi_Mix", "hMC_USS_KPi_Mix", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); + histos.add("hMC_USS_PiK_Mix", "hMC_USS_PiK_Mix", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); + histos.add("hMC_USS_KPi_True", "hMC_USS_KPi_True", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); + histos.add("hMC_USS_PiK_True", "hMC_USS_PiK_True", kTHnSparseF, {cfgCentAxis, ptAxis, minvAxis}); } } // end of init @@ -226,7 +230,6 @@ struct kstarInOO { histos.fill(HIST("hPosZ_BC"), event.posZ()); histos.fill(HIST("hcentFT0C_BC"), event.centFT0C()); } - if (!event.sel8()) return false; if (std::abs(event.posZ()) > cfgEventVtxCut) @@ -264,6 +267,7 @@ struct kstarInOO { histos.fill(HIST("hTPCChi2_BC"), track.tpcChi2NCl()); histos.fill(HIST("QA_track_pT_BC"), track.pt()); } + if (cfgTrackGlobalSel && !track.isGlobalTrack()) return false; if (track.pt() < cfgTrackMinPt) @@ -278,11 +282,11 @@ struct kstarInOO { return false; if (cfgTrackGlobalWoDCATrack && !track.isGlobalTrackWoDCA()) return false; - if (track.tpcNClsFindable() < cfgTracknFindableTPCClusters) + if (cfgTracknFindableTPCClusters > 0 && track.tpcNClsFindable() < cfgTracknFindableTPCClusters) return false; if (track.tpcNClsCrossedRows() < cfgTracknTPCCrossedRows) return false; - if (track.tpcCrossedRowsOverFindableCls() > cfgTracknRowsOverFindable) + if (cfgTracknRowsOverFindable > 0 && track.tpcCrossedRowsOverFindableCls() > cfgTracknRowsOverFindable) return false; if (track.tpcChi2NCl() > cfgTracknTPCChi2) return false; @@ -358,7 +362,7 @@ struct kstarInOO { } else { tofPIDPassed = true; } - // TPC & TOF + // TPC & TOF if (tpcPIDPassed && tofPIDPassed) { if (cfgTrackCutQA && QA) { histos.fill(HIST("QA_nSigma_pion_TPC_AC"), candidate.pt(), candidate.tpcNSigmaPi()); @@ -378,10 +382,8 @@ struct kstarInOO { auto centrality = collision1.centFT0C(); for (const auto& [trk1, trk2] : combinations(o2::soa::CombinationsFullIndexPolicy(tracks1, tracks2))) { - if (std::fabs(trk1.signed1Pt()) <= 0.f || std::fabs(trk2.signed1Pt()) <= 0.f) - continue; - auto [KstarPt, Minv] = minvReconstruction(trk1, trk2, QA); + auto [KstarPt, Minv] = minvReconstruction(trk1, trk2, QA, false); double conjugate = trk1.sign() * trk2.sign(); if (cfgDataHistos) { @@ -411,28 +413,36 @@ struct kstarInOO { auto tracks2 = pionMC->sliceByCached(aod::track::collisionId, collision2.globalIndex(), cache); auto centrality = collision1.centFT0C(); + std::vector mcMemory; + std::vector PIDPurityKey_Kaon; + std::vector PIDPurityKey_Pion; + + double KstarPt_Kpi, Minv_Kpi; + for (const auto& [trk1, trk2] : combinations(o2::soa::CombinationsFullIndexPolicy(tracks1, tracks2))) { if (!trk1.has_mcParticle() || !trk2.has_mcParticle()) continue; - if (std::fabs(trk1.signed1Pt()) <= 0.f || std::fabs(trk2.signed1Pt()) <= 0.f) + + // auto [KstarPt_Kpi, Minv_Kpi] = minvReconstruction(trk1, trk2, QA, false); + // auto [KstarPt_piK, Minv_piK] = minvReconstruction(trk1, trk2, QA, true); + + std::tie(KstarPt_Kpi, Minv_Kpi) = minvReconstruction(trk1, trk2, QA, false); + std::tie(KstarPt_Kpi, Minv_Kpi) = minvReconstruction(trk1, trk2, QA, true); + + if (Minv_Kpi < 0) continue; - auto [KstarPt_Kpi, Minv_Kpi] = minvReconstruction(trk1, trk2, QA); double conjugate = trk1.sign() * trk2.sign(); if (cfgMcHistos) { - if (Minv_Kpi > 0) { - if (!IsMix) { - if (conjugate < 0) { - histos.fill(HIST("hMC_USS"), centrality, KstarPt_Kpi, Minv_Kpi); - } else if (conjugate > 0) { - histos.fill(HIST("hMC_LSS"), centrality, KstarPt_Kpi, Minv_Kpi); - } - } else { - if (conjugate < 0) { - histos.fill(HIST("hMC_USS_Mix"), centrality, KstarPt_Kpi, Minv_Kpi); - } else if (conjugate > 0) { - histos.fill(HIST("hMC_LSS_Mix"), centrality, KstarPt_Kpi, Minv_Kpi); - } + if (!IsMix) { + if (conjugate < 0) { + histos.fill(HIST("hMC_USS_KPi"), centrality, KstarPt_Kpi, Minv_Kpi); + } else if (conjugate > 0) { + histos.fill(HIST("hMC_LSS_KPi"), centrality, KstarPt_Kpi, Minv_Kpi); + } + } else { + if (conjugate < 0) { + histos.fill(HIST("hMC_USS_KPi_Mix"), centrality, KstarPt_Kpi, Minv_Kpi); } } } @@ -440,14 +450,12 @@ struct kstarInOO { // Gen MC auto particle1 = trk1.mcParticle(); auto particle2 = trk2.mcParticle(); - if (std::fabs(particle1.pdgCode()) != 321) - continue; // Not Kaon - if (std::fabs(particle2.pdgCode()) != 211) - continue; // Not Pion if (!particle1.has_mothers() || !particle2.has_mothers()) { continue; } + int mcindex1 = trk1.globalIndex(); + int mcindex2 = trk2.globalIndex(); std::vector mothers1{}; std::vector mothers1PDG{}; @@ -471,27 +479,85 @@ struct kstarInOO { if (mothers1[0] != mothers2[0]) continue; // Kaon and pion not from the same K*0 + if (std::fabs(particle1.pdgCode()) != 211 && std::fabs(particle1.pdgCode()) != 321) + continue; + if (std::fabs(particle2.pdgCode()) != 211 && std::fabs(particle2.pdgCode()) != 321) + continue; + + double track1_mass, track2_mass; + bool track1f{false}; // true means pion + + if (std::fabs(particle1.pdgCode()) == 211) { + track1f = true; + track1_mass = massPi; + } else { + track1_mass = massKa; + } + + if (std::fabs(particle2.pdgCode()) == 211) { + track2_mass = massPi; + } else { + track2_mass = massKa; + } + + if (track1_mass == track2_mass) { + return; + } + + bool exists1 = std::find(mcMemory.begin(), mcMemory.end(), mcindex1) != mcMemory.end(); + bool exists2 = std::find(mcMemory.begin(), mcMemory.end(), mcindex2) != mcMemory.end(); + if (exists1 || exists2) { + continue; + } else { + mcMemory.push_back(trk1.globalIndex()); + mcMemory.push_back(trk2.globalIndex()); + } + + TLorentzVector lDecayDaughter1, lDecayDaughter2, lResonance; + lDecayDaughter1.SetXYZM(trk1.px(), trk1.py(), trk1.pz(), track1_mass); + lDecayDaughter2.SetXYZM(trk2.px(), trk2.py(), trk2.pz(), track2_mass); + lResonance = lDecayDaughter1 + lDecayDaughter2; + if (cfgMcHistos) { - histos.fill(HIST("hMC_USS_True"), centrality, KstarPt_Kpi, Minv_Kpi); + histos.fill(HIST("hMC_USS_True"), centrality, lResonance.Pt(), lResonance.M()); + if (track1f) { + histos.fill(HIST("hMC_USS_PiK_True"), centrality, lResonance.Pt(), lResonance.M()); + } else { + histos.fill(HIST("hMC_USS_KPi_True"), centrality, lResonance.Pt(), lResonance.M()); + } } //====================== } // for } // TrackSlicingMC template - std::pair minvReconstruction(const TracksType& trk1, const TracksType& trk2, const bool QA) + std::pair minvReconstruction(const TracksType& trk1, const TracksType& trk2, const bool QA, const bool flip) { if (!trackSelection(trk1, false) || !trackSelection(trk2, false)) return {-1.0, -1.0}; - if (!trackPIDKaon(trk1, QA) || !trackPIDPion(trk2, QA)) - return {-1.0, -1.0}; - if (trk1.globalIndex() == trk2.globalIndex()) + + if (!flip) { + if (!trackPIDKaon(trk1, QA) || !trackPIDPion(trk2, QA)) { + return {-1.0, -1.0}; + } + } else { + if (!trackPIDPion(trk1, false) || !trackPIDKaon(trk2, false)) + return {-1.0, -1.0}; + } + + // if (trk1.globalIndex() == trk2.globalIndex()) + // return {-1.0, -1.0}; + if (trk1.index() >= trk2.index()) return {-1.0, -1.0}; TLorentzVector lDecayDaughter1, lDecayDaughter2, lResonance; - - lDecayDaughter1.SetXYZM(trk1.px(), trk1.py(), trk1.pz(), massKa); - lDecayDaughter2.SetXYZM(trk2.px(), trk2.py(), trk2.pz(), massPi); + if (!flip) { + lDecayDaughter1.SetXYZM(trk1.px(), trk1.py(), trk1.pz(), massKa); + lDecayDaughter2.SetXYZM(trk2.px(), trk2.py(), trk2.pz(), massPi); + } else { + lDecayDaughter1.SetXYZM(trk1.px(), trk1.py(), trk1.pz(), massPi); + lDecayDaughter2.SetXYZM(trk2.px(), trk2.py(), trk2.pz(), massKa); + } lResonance = lDecayDaughter1 + lDecayDaughter2; if (std::abs(lResonance.Eta()) > cfgTrackMaxEta) @@ -573,7 +639,7 @@ struct kstarInOO { //| //| MC STUFF (SE) //| - //======================================================= + //========================================================= int nEventsMC = 0; void processSameEventMC(EventCandidates::iterator const& collision, TrackCandidatesMC const& tracks, aod::McParticles const&) { @@ -603,6 +669,36 @@ struct kstarInOO { if (!INELgt0) return; + auto tracks1 = kaonMC->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + for (const auto& kaon : tracks1) { + if (!trackSelection(kaon, false)) + continue; + if (!trackPIDKaon(kaon, false)) + continue; + auto particle1 = kaon.mcParticle(); + if (std::fabs(particle1.pdgCode()) == 321) + histos.fill(HIST("hSimpleKaon_PID_Purity"), 1); // histogram with two bins, -1.5, 1.5 fill 1 or -1 + else if (std::fabs(particle1.pdgCode()) == 211) + histos.fill(HIST("hSimpleKaon_PID_Purity"), -1); // histogram with two bins, -1.5, 1.5 fill 1 or -1 + else + histos.fill(HIST("hSimpleKaon_PID_Purity"), 0); // histogram with two bins, -1.5, 1.5 fill 1 or -1 + } + + auto tracks2 = pionMC->sliceByCached(aod::track::collisionId, collision.globalIndex(), cache); + for (const auto& pion : tracks2) { + if (!trackSelection(pion, false)) + continue; + if (!trackPIDPion(pion, false)) + continue; + auto particle2 = pion.mcParticle(); + if (std::fabs(particle2.pdgCode()) == 211) + histos.fill(HIST("hSimplePion_PID_Purity"), 1); // histogram with two bins, -1.5, 1.5 fill 1 or -1 + else if (std::fabs(particle2.pdgCode()) == 321) + histos.fill(HIST("hSimplePion_PID_Purity"), -1); // histogram with two bins, -1.5, 1.5 fill 1 or -1 + else + histos.fill(HIST("hSimplePion_PID_Purity"), 0); // histogram with two bins, -1.5, 1.5 fill 1 or -1 + } + if (cfgMcHistos) { histos.fill(HIST("nEvents_MC"), 1.5); }