Skip to content

Commit 410c38a

Browse files
skundu692skundu692
committed
Add new process function to remove cross pairing background
1 parent f296787 commit 410c38a

1 file changed

Lines changed: 209 additions & 1 deletion

File tree

‎PWGLF/Tasks/Resonances/doublephimeson.cxx‎

Lines changed: 209 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -806,7 +806,7 @@ struct doublephimeson {
806806
const double dR_k1p_k2m = deltaR(kplusA.Phi(), kplusA.Eta(), kminusB.Phi(), kminusB.Eta());
807807
const double dR_k2p_k1m = deltaR(kplusB.Phi(), kplusB.Eta(), kminusA.Phi(), kminusA.Eta());
808808
const double dR_k2p_k2m = deltaR(kplusB.Phi(), kplusB.Eta(), kminusB.Phi(), kminusB.Eta());
809-
809+
810810
double minDR = dRkplus;
811811
minDR = std::min(minDR, dRkminus);
812812
minDR = std::min(minDR, dR_k1p_k1m);
@@ -926,6 +926,214 @@ struct doublephimeson {
926926
}
927927
PROCESS_SWITCH(doublephimeson, processopti3, "Process Optimized same event", false);
928928

929+
void processopti4(aod::RedPhiEvents::iterator const& collision, aod::PhiTracks const& phitracks)
930+
{
931+
if (additionalEvsel && (collision.numPos() < 2 || collision.numNeg() < 2))
932+
return;
933+
934+
// --- φ multiplicity with PID ---
935+
int phimult = 0;
936+
for (auto const& t : phitracks) {
937+
const double kpluspt = std::hypot(t.phid1Px(), t.phid1Py());
938+
const double kminuspt = std::hypot(t.phid2Px(), t.phid2Py());
939+
940+
// PID QA before
941+
histos.fill(HIST("hnsigmaTPCTOFKaonBefore"), t.phid1TPC(), t.phid1TOF(), kpluspt);
942+
histos.fill(HIST("hnsigmaTPCKaonPlusBefore"), t.phid1TPC(), kpluspt);
943+
histos.fill(HIST("hnsigmaTPCKaonMinusBefore"), t.phid2TPC(), kminuspt);
944+
945+
if (t.phiMass() < minPhiMass1 || t.phiMass() > maxPhiMass1)
946+
continue;
947+
if (kpluspt > maxKaonPt || kminuspt > maxKaonPt)
948+
continue;
949+
if (!selectionPID(t.phid1TPC(), t.phid1TOF(), t.phid1TOFHit(), strategyPID1, kpluspt))
950+
continue;
951+
if (!selectionPID(t.phid2TPC(), t.phid2TOF(), t.phid2TOFHit(), strategyPID2, kminuspt))
952+
continue;
953+
954+
// PID QA after
955+
histos.fill(HIST("hnsigmaTPCTOFKaon"), t.phid1TPC(), t.phid1TOF(), kpluspt);
956+
histos.fill(HIST("hnsigmaTPCKaonPlus"), t.phid1TPC(), kpluspt);
957+
histos.fill(HIST("hnsigmaTPCKaonMinus"), t.phid2TPC(), kminuspt);
958+
959+
++phimult;
960+
}
961+
if (phimult < 2)
962+
return;
963+
964+
// --- helpers ---
965+
constexpr double mPhiPDG = 1.019461; // GeV/c^2
966+
constexpr double mKPDG = 0.493677; // GeV/c^2
967+
968+
const auto deltaMPhi = [=](double m1, double m2) {
969+
const double d1 = m1 - mPhiPDG, d2 = m2 - mPhiPDG;
970+
return std::sqrt(d1 * d1 + d2 * d2);
971+
};
972+
973+
const auto deltaR = [](double phi1, double eta1, double phi2, double eta2) {
974+
const double dphi = TVector2::Phi_mpi_pi(phi1 - phi2);
975+
const double deta = eta1 - eta2;
976+
return std::sqrt(dphi * dphi + deta * deta);
977+
};
978+
979+
// minimum ΔR among all kaons in the candidate (4 kaons → 6 combinations)
980+
const auto minKaonDeltaR = [&](const ROOT::Math::PtEtaPhiMVector& kplusA,
981+
const ROOT::Math::PtEtaPhiMVector& kplusB,
982+
const ROOT::Math::PtEtaPhiMVector& kminusA,
983+
const ROOT::Math::PtEtaPhiMVector& kminusB) {
984+
// same-sign first (keep your QA histos)
985+
const double dRkplus = deltaR(kplusA.Phi(), kplusA.Eta(), kplusB.Phi(), kplusB.Eta());
986+
const double dRkminus = deltaR(kminusA.Phi(), kminusA.Eta(), kminusB.Phi(), kminusB.Eta());
987+
histos.fill(HIST("hDeltaRkaonplus"), dRkplus);
988+
histos.fill(HIST("hDeltaRkaonminus"), dRkminus);
989+
990+
// all other combinations
991+
const double dR_k1p_k1m = deltaR(kplusA.Phi(), kplusA.Eta(), kminusA.Phi(), kminusA.Eta());
992+
const double dR_k1p_k2m = deltaR(kplusA.Phi(), kplusA.Eta(), kminusB.Phi(), kminusB.Eta());
993+
const double dR_k2p_k1m = deltaR(kplusB.Phi(), kplusB.Eta(), kminusA.Phi(), kminusA.Eta());
994+
const double dR_k2p_k2m = deltaR(kplusB.Phi(), kplusB.Eta(), kminusB.Phi(), kminusB.Eta());
995+
996+
double minDR = dRkplus;
997+
minDR = std::min(minDR, dRkminus);
998+
minDR = std::min(minDR, dR_k1p_k1m);
999+
minDR = std::min(minDR, dR_k1p_k2m);
1000+
minDR = std::min(minDR, dR_k2p_k1m);
1001+
minDR = std::min(minDR, dR_k2p_k2m);
1002+
return minDR;
1003+
};
1004+
1005+
// --- collect candidates once ---
1006+
std::vector<ROOT::Math::PtEtaPhiMVector> pairV, phi1V, phi2V;
1007+
std::vector<double> minDRV; // store minimum ΔR for each pair
1008+
1009+
// optional: swapped-cross-mass veto window (minimal new knobs)
1010+
const double crossPhiLow = 1.01; // or a dedicated config
1011+
const double crossPhiHigh = 1.03; // or a dedicated config
1012+
1013+
for (auto const& t1 : phitracks) {
1014+
const double kplus1pt = std::hypot(t1.phid1Px(), t1.phid1Py());
1015+
const double kminus1pt = std::hypot(t1.phid2Px(), t1.phid2Py());
1016+
1017+
if (kplus1pt > maxKaonPt || kminus1pt > maxKaonPt)
1018+
continue;
1019+
if (!selectionPID(t1.phid1TPC(), t1.phid1TOF(), t1.phid1TOFHit(), strategyPID1, kplus1pt))
1020+
continue;
1021+
if (!selectionPID(t1.phid2TPC(), t1.phid2TOF(), t1.phid2TOFHit(), strategyPID2, kminus1pt))
1022+
continue;
1023+
1024+
TLorentzVector phi1, k1p, k1m;
1025+
phi1.SetXYZM(t1.phiPx(), t1.phiPy(), t1.phiPz(), t1.phiMass());
1026+
k1p.SetXYZM(t1.phid1Px(), t1.phid1Py(), t1.phid1Pz(), mKPDG);
1027+
k1m.SetXYZM(t1.phid2Px(), t1.phid2Py(), t1.phid2Pz(), mKPDG);
1028+
1029+
// φ1 mass window + φ1 pT
1030+
if (t1.phiMass() < minPhiMass1 || t1.phiMass() > maxPhiMass1)
1031+
continue;
1032+
if (phi1.Pt() < minPhiPt || phi1.Pt() > maxPhiPt)
1033+
continue;
1034+
1035+
const auto id1 = t1.index();
1036+
1037+
for (auto const& t2 : phitracks) {
1038+
const auto id2 = t2.index();
1039+
if (id2 <= id1)
1040+
continue;
1041+
1042+
const double kplus2pt = std::hypot(t2.phid1Px(), t2.phid1Py());
1043+
const double kminus2pt = std::hypot(t2.phid2Px(), t2.phid2Py());
1044+
if (kplus2pt > maxKaonPt || kminus2pt > maxKaonPt)
1045+
continue;
1046+
if (!selectionPID(t2.phid1TPC(), t2.phid1TOF(), t2.phid1TOFHit(), strategyPID1, kplus2pt))
1047+
continue;
1048+
if (!selectionPID(t2.phid2TPC(), t2.phid2TOF(), t2.phid2TOFHit(), strategyPID2, kminus2pt))
1049+
continue;
1050+
1051+
// FIX + robust: block ANY shared daughter (4-way)
1052+
if (t1.phid1Index() == t2.phid1Index() || t1.phid1Index() == t2.phid2Index() ||
1053+
t1.phid2Index() == t2.phid1Index() || t1.phid2Index() == t2.phid2Index())
1054+
continue;
1055+
1056+
TLorentzVector phi2, k2p, k2m;
1057+
phi2.SetXYZM(t2.phiPx(), t2.phiPy(), t2.phiPz(), t2.phiMass());
1058+
k2p.SetXYZM(t2.phid1Px(), t2.phid1Py(), t2.phid1Pz(), mKPDG);
1059+
k2m.SetXYZM(t2.phid2Px(), t2.phid2Py(), t2.phid2Pz(), mKPDG);
1060+
1061+
// φ2 mass window + FIX: apply pT cut to phi2 (not phi1)
1062+
if (t2.phiMass() < minPhiMass2 || t2.phiMass() > maxPhiMass2)
1063+
continue;
1064+
if (phi2.Pt() < minPhiPt || phi2.Pt() > maxPhiPt)
1065+
continue;
1066+
1067+
// NEW: cross (swapped) K+K- mass veto
1068+
// veto if either cross-pair lands in φ window (tune to be looser/tighter)
1069+
const double mCross12 = (k1p + k2m).M(); // K+1 + K-2
1070+
const double mCross21 = (k2p + k1m).M(); // K+2 + K-1
1071+
if ((mCross12 > crossPhiLow && mCross12 < crossPhiHigh) ||
1072+
(mCross21 > crossPhiLow && mCross21 < crossPhiHigh)) {
1073+
continue;
1074+
}
1075+
1076+
// Δm cut (configurable)
1077+
const double dM = deltaMPhi(phi1.M(), phi2.M());
1078+
if (dM > maxDeltaMPhi)
1079+
continue;
1080+
1081+
TLorentzVector pair = phi1 + phi2;
1082+
if (pair.M() < minExoticMass || pair.M() > maxExoticMass)
1083+
continue;
1084+
1085+
histos.fill(HIST("hPhiMass"), phi1.M(), phi2.M(), pair.Pt());
1086+
1087+
// daughter ΔR QA and minΔR (NO CUT)
1088+
ROOT::Math::PtEtaPhiMVector k1pV(k1p.Pt(), k1p.Eta(), k1p.Phi(), mKPDG);
1089+
ROOT::Math::PtEtaPhiMVector k1mV(k1m.Pt(), k1m.Eta(), k1m.Phi(), mKPDG);
1090+
ROOT::Math::PtEtaPhiMVector k2pV(k2p.Pt(), k2p.Eta(), k2p.Phi(), mKPDG);
1091+
ROOT::Math::PtEtaPhiMVector k2mV(k2m.Pt(), k2m.Eta(), k2m.Phi(), mKPDG);
1092+
const double minDR = minKaonDeltaR(k1pV, k2pV, k1mV, k2mV);
1093+
1094+
// store for one-pass fill
1095+
pairV.emplace_back(pair.Pt(), pair.Eta(), pair.Phi(), pair.M());
1096+
phi1V.emplace_back(phi1.Pt(), phi1.Eta(), phi1.Phi(), phi1.M());
1097+
phi2V.emplace_back(phi2.Pt(), phi2.Eta(), phi2.Phi(), phi2.M());
1098+
minDRV.emplace_back(minDR);
1099+
}
1100+
}
1101+
1102+
if (pairV.empty())
1103+
return;
1104+
1105+
// --- fill the single THnSparse ---
1106+
for (size_t i = 0; i < pairV.size(); ++i) {
1107+
TLorentzVector p1, p2, pair;
1108+
p1.SetPtEtaPhiM(phi1V[i].Pt(), phi1V[i].Eta(), phi1V[i].Phi(), phi1V[i].M());
1109+
p2.SetPtEtaPhiM(phi2V[i].Pt(), phi2V[i].Eta(), phi2V[i].Phi(), phi2V[i].M());
1110+
pair.SetPtEtaPhiM(pairV[i].Pt(), pairV[i].Eta(), pairV[i].Phi(), pairV[i].M());
1111+
1112+
const double dM = deltaMPhi(p1.M(), p2.M());
1113+
const double M = pair.M();
1114+
const double dR = deltaR(p1.Phi(), p1.Eta(), p2.Phi(), p2.Eta());
1115+
const double minDR = minDRV[i];
1116+
1117+
// (optional but recommended) protect ptcorr from blow-ups
1118+
const double denom = (pair.Pt() - p1.Pt());
1119+
if (std::abs(denom) < 1e-9)
1120+
continue;
1121+
const double ptcorr = p1.Pt() / denom;
1122+
1123+
histos.fill(HIST("hPtCorrelation"), pair.Pt(), ptcorr);
1124+
1125+
// NOTE: second axis is minΔR(all kaons), ΔpT/pT has been removed
1126+
histos.fill(HIST("SEMassUnlike"),
1127+
M,
1128+
minDR,
1129+
pair.Pt(),
1130+
dR,
1131+
dM,
1132+
ptcorr);
1133+
}
1134+
}
1135+
PROCESS_SWITCH(doublephimeson, processopti4, "Process Optimized same event", true);
1136+
9291137
SliceCache cache;
9301138
using BinningTypeVertexContributor = ColumnBinningPolicy<aod::collision::PosZ, aod::collision::NumContrib>;
9311139

0 commit comments

Comments
 (0)