Skip to content

Commit d53549f

Browse files
skundu692skundu692
authored
[PWGLF] Add mixed event process function (#13841)
1 parent 1e603cd commit d53549f

File tree

1 file changed

+189
-1
lines changed

1 file changed

+189
-1
lines changed

‎PWGLF/Tasks/Resonances/doublephimeson.cxx‎

Lines changed: 189 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -112,7 +112,7 @@ struct doublephimeson {
112112

113113
histos.add("SEMassUnlike", "SEMassUnlike", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisDeltaR, thnAxisPt, thnAxisDeltaR, thnAxisInvMassDeltaPhi, thnAxisPtCorr});
114114
// histos.add("SEMassLike", "SEMassLike", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisDeltaR, thnAxisInvMassPhi, thnAxisInvMassPhi, thnAxisNumPhi});
115-
histos.add("MEMassUnlike", "MEMassUnlike", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisDeltaR, thnAxisInvMassDeltaPhi});
115+
histos.add("MEMassUnlike", "MEMassUnlike", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisDeltaR, thnAxisPt, thnAxisDeltaR, thnAxisInvMassDeltaPhi, thnAxisPtCorr});
116116
}
117117

118118
// get kstar
@@ -925,6 +925,194 @@ struct doublephimeson {
925925
}
926926
}
927927
PROCESS_SWITCH(doublephimeson, processopti3, "Process Optimized same event", false);
928+
929+
SliceCache cache;
930+
using BinningTypeVertexContributor = ColumnBinningPolicy<aod::collision::PosZ, aod::collision::NumContrib>;
931+
932+
void processMixedEvent(aod::RedPhiEvents& collisions, aod::PhiTracks& phitracks)
933+
{
934+
auto tracksTuple = std::make_tuple(phitracks);
935+
BinningTypeVertexContributor binningOnPositions{{CfgVtxBins, CfgMultBins}, true};
936+
SameKindPair<aod::RedPhiEvents, aod::PhiTracks, BinningTypeVertexContributor> pair{
937+
binningOnPositions, nEvtMixing, -1, collisions, tracksTuple, &cache};
938+
939+
// --- helpers (same as in processopti3) ---
940+
constexpr double mPhiPDG = 1.019461; // GeV/c^2
941+
942+
const auto deltaMPhi = [=](double m1, double m2) {
943+
const double d1 = m1 - mPhiPDG;
944+
const double d2 = m2 - mPhiPDG;
945+
return std::sqrt(d1 * d1 + d2 * d2);
946+
};
947+
948+
const auto deltaR = [](double phi1, double eta1, double phi2, double eta2) {
949+
const double dphi = TVector2::Phi_mpi_pi(phi1 - phi2);
950+
const double deta = eta1 - eta2;
951+
return std::sqrt(dphi * dphi + deta * deta);
952+
};
953+
954+
const auto minKaonDeltaR =
955+
[&](const ROOT::Math::PtEtaPhiMVector& kplusA,
956+
const ROOT::Math::PtEtaPhiMVector& kplusB,
957+
const ROOT::Math::PtEtaPhiMVector& kminusA,
958+
const ROOT::Math::PtEtaPhiMVector& kminusB) {
959+
// same-sign first (keep QA as in SE)
960+
const double dRkplus =
961+
deltaR(kplusA.Phi(), kplusA.Eta(), kplusB.Phi(), kplusB.Eta());
962+
const double dRkminus =
963+
deltaR(kminusA.Phi(), kminusA.Eta(), kminusB.Phi(), kminusB.Eta());
964+
histos.fill(HIST("hDeltaRkaonplus"), dRkplus);
965+
histos.fill(HIST("hDeltaRkaonminus"), dRkminus);
966+
967+
// all other combinations
968+
const double dR_k1p_k1m =
969+
deltaR(kplusA.Phi(), kplusA.Eta(), kminusA.Phi(), kminusA.Eta());
970+
const double dR_k1p_k2m =
971+
deltaR(kplusA.Phi(), kplusA.Eta(), kminusB.Phi(), kminusB.Eta());
972+
const double dR_k2p_k1m =
973+
deltaR(kplusB.Phi(), kplusB.Eta(), kminusA.Phi(), kminusA.Eta());
974+
const double dR_k2p_k2m =
975+
deltaR(kplusB.Phi(), kplusB.Eta(), kminusB.Phi(), kminusB.Eta());
976+
977+
double minDR = dRkplus;
978+
minDR = std::min(minDR, dRkminus);
979+
minDR = std::min(minDR, dR_k1p_k1m);
980+
minDR = std::min(minDR, dR_k1p_k2m);
981+
minDR = std::min(minDR, dR_k2p_k1m);
982+
minDR = std::min(minDR, dR_k2p_k2m);
983+
return minDR;
984+
};
985+
986+
struct PhiCand {
987+
ROOT::Math::PtEtaPhiMVector phi;
988+
ROOT::Math::PtEtaPhiMVector kplus;
989+
ROOT::Math::PtEtaPhiMVector kminus;
990+
};
991+
992+
for (auto& [collision1, tracks1, collision2, tracks2] : pair) {
993+
// safety: should never happen but keep it
994+
if (collision1.index() == collision2.index()) {
995+
continue;
996+
}
997+
998+
// optional event-level selection (same idea as in SE)
999+
if (additionalEvsel) {
1000+
if (collision1.numPos() < 2 || collision1.numNeg() < 2) {
1001+
continue;
1002+
}
1003+
if (collision2.numPos() < 2 || collision2.numNeg() < 2) {
1004+
continue;
1005+
}
1006+
}
1007+
1008+
std::vector<PhiCand> cands1, cands2;
1009+
1010+
// --- build φ candidates for event 1 (φ1) ---
1011+
for (auto const& t1 : tracks1) {
1012+
const double kplus1pt = std::hypot(t1.phid1Px(), t1.phid1Py());
1013+
const double kminus1pt = std::hypot(t1.phid2Px(), t1.phid2Py());
1014+
1015+
if (kplus1pt > maxKaonPt || kminus1pt > maxKaonPt)
1016+
continue;
1017+
if (!selectionPID(t1.phid1TPC(), t1.phid1TOF(), t1.phid1TOFHit(), strategyPID1, kplus1pt))
1018+
continue;
1019+
if (!selectionPID(t1.phid2TPC(), t1.phid2TOF(), t1.phid2TOFHit(), strategyPID2, kminus1pt))
1020+
continue;
1021+
1022+
TLorentzVector phi1, k1p, k1m;
1023+
phi1.SetXYZM(t1.phiPx(), t1.phiPy(), t1.phiPz(), t1.phiMass());
1024+
k1p.SetXYZM(t1.phid1Px(), t1.phid1Py(), t1.phid1Pz(), 0.493);
1025+
k1m.SetXYZM(t1.phid2Px(), t1.phid2Py(), t1.phid2Pz(), 0.493);
1026+
1027+
if (t1.phiMass() < minPhiMass1 || t1.phiMass() > maxPhiMass1)
1028+
continue;
1029+
if (phi1.Pt() < minPhiPt || phi1.Pt() > maxPhiPt)
1030+
continue;
1031+
1032+
histos.fill(HIST("hPhiMass"), phi1.M(), phi1.Pt());
1033+
1034+
PhiCand cand;
1035+
cand.phi = ROOT::Math::PtEtaPhiMVector(phi1.Pt(), phi1.Eta(), phi1.Phi(), phi1.M());
1036+
cand.kplus = ROOT::Math::PtEtaPhiMVector(k1p.Pt(), k1p.Eta(), k1p.Phi(), 0.493);
1037+
cand.kminus = ROOT::Math::PtEtaPhiMVector(k1m.Pt(), k1m.Eta(), k1m.Phi(), 0.493);
1038+
1039+
cands1.emplace_back(std::move(cand));
1040+
}
1041+
1042+
// --- build φ candidates for event 2 (φ2) ---
1043+
for (auto const& t2 : tracks2) {
1044+
const double kplus2pt = std::hypot(t2.phid1Px(), t2.phid1Py());
1045+
const double kminus2pt = std::hypot(t2.phid2Px(), t2.phid2Py());
1046+
1047+
if (kplus2pt > maxKaonPt || kminus2pt > maxKaonPt)
1048+
continue;
1049+
if (!selectionPID(t2.phid1TPC(), t2.phid1TOF(), t2.phid1TOFHit(), strategyPID1, kplus2pt))
1050+
continue;
1051+
if (!selectionPID(t2.phid2TPC(), t2.phid2TOF(), t2.phid2TOFHit(), strategyPID2, kminus2pt))
1052+
continue;
1053+
1054+
TLorentzVector phi2, k2p, k2m;
1055+
phi2.SetXYZM(t2.phiPx(), t2.phiPy(), t2.phiPz(), t2.phiMass());
1056+
k2p.SetXYZM(t2.phid1Px(), t2.phid1Py(), t2.phid1Pz(), 0.493);
1057+
k2m.SetXYZM(t2.phid2Px(), t2.phid2Py(), t2.phid2Pz(), 0.493);
1058+
1059+
if (t2.phiMass() < minPhiMass2 || t2.phiMass() > maxPhiMass2)
1060+
continue;
1061+
if (phi2.Pt() < minPhiPt || phi2.Pt() > maxPhiPt)
1062+
continue;
1063+
1064+
histos.fill(HIST("hPhiMass"), phi2.M(), phi2.Pt());
1065+
1066+
PhiCand cand;
1067+
cand.phi = ROOT::Math::PtEtaPhiMVector(phi2.Pt(), phi2.Eta(), phi2.Phi(), phi2.M());
1068+
cand.kplus = ROOT::Math::PtEtaPhiMVector(k2p.Pt(), k2p.Eta(), k2p.Phi(), 0.493);
1069+
cand.kminus = ROOT::Math::PtEtaPhiMVector(k2m.Pt(), k2m.Eta(), k2m.Phi(), 0.493);
1070+
1071+
cands2.emplace_back(std::move(cand));
1072+
}
1073+
1074+
if (cands1.empty() || cands2.empty()) {
1075+
continue;
1076+
}
1077+
1078+
// --- build mixed-event pairs and fill MEMassUnlike ---
1079+
for (auto const& c1 : cands1) {
1080+
TLorentzVector phi1;
1081+
phi1.SetPtEtaPhiM(c1.phi.Pt(), c1.phi.Eta(), c1.phi.Phi(), c1.phi.M());
1082+
1083+
for (auto const& c2 : cands2) {
1084+
TLorentzVector phi2;
1085+
phi2.SetPtEtaPhiM(c2.phi.Pt(), c2.phi.Eta(), c2.phi.Phi(), c2.phi.M());
1086+
1087+
const double dM = deltaMPhi(phi1.M(), phi2.M());
1088+
if (dM > maxDeltaMPhi)
1089+
continue;
1090+
1091+
TLorentzVector pair = phi1 + phi2;
1092+
if (pair.M() < minExoticMass || pair.M() > maxExoticMass)
1093+
continue;
1094+
1095+
const double minDR = minKaonDeltaR(c1.kplus, c2.kplus, c1.kminus, c2.kminus);
1096+
const double dR = deltaR(phi1.Phi(), phi1.Eta(), phi2.Phi(), phi2.Eta());
1097+
1098+
// same definition as SE
1099+
const double ptcorr = (pair.Pt() - phi1.Pt() != 0.)
1100+
? phi1.Pt() / (pair.Pt() - phi1.Pt())
1101+
: 0.;
1102+
1103+
histos.fill(HIST("MEMassUnlike"),
1104+
pair.M(), // M(phi-phi)
1105+
minDR, // min ΔR among all kaon pairs
1106+
pair.Pt(), // pT(phi-phi)
1107+
dR, // ΔR(phi1, phi2)
1108+
dM, // Δm(phi)
1109+
ptcorr); // pT correlation
1110+
}
1111+
}
1112+
}
1113+
}
1114+
PROCESS_SWITCH(doublephimeson, processMixedEvent,
1115+
"Process EventMixing for combinatorial background", false);
9281116
};
9291117

9301118
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask<doublephimeson>(cfgc)}; }

0 commit comments

Comments
 (0)