Skip to content

[PWGEM] Add new event mixing logic based on photon legs #16956

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Closed
Steffimro wants to merge 5 commits into from
Closed

[PWGEM] Add new event mixing logic based on photon legs #16956

Changes from all commits
Commits
Show all changes
5 commits
Select commit Hold shift + click to select a range
674f9e7
Store leg kinematics in photon mixing pool. Replace EMPair with a loc…
Jul 7, 2026
3373a98
Bug fix: fix double filling of hSparse_DEtaDPhi_qinv_kT in MC pairing
Jul 7, 2026
3e65e30
Bug fix: e+/e- assignement in efficiency correctly implemented
Jul 7, 2026
b6a67e0
Add SE/ME same-sign leg separation QA for pair efficiency extraction
Jul 7, 2026
e3476dc
Fix o2-linter findings
Jul 7, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
129 changes: 104 additions & 25 deletions PWGEM/PhotonMeson/Tasks/photonhbt.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -302,6 +302,13 @@ struct Photonhbt {
Configurable<float> cfgMaxAsymmetry{"cfgMaxAsymmetry", -1.f, "max |p_{T, 1} - p_{T, 2}|/(p_{T, 1} + p_{T, 2}) asymmetry cut"};
} ggpaircuts;

struct : ConfigurableGroup {
std::string prefix = "pairsep_group";
Configurable<bool> cfgDoPairSepQA{"cfgDoPairSepQA", true, "fill SE/ME same-sign leg separation histograms"};
Configurable<float> cfgUSigEta{"cfgUSigEta", 0.02f, "eta scale defining u (axis units only, NOT a cut)"};
Configurable<float> cfgUSigPhi{"cfgUSigPhi", 0.05f, "phi scale defining u (axis units only, NOT a cut)"};
} pairsep;

EMPhotonEventCut fEMEventCut;
struct : ConfigurableGroup {
std::string prefix = "eventcut_group";
Expand Down Expand Up @@ -588,6 +595,24 @@ struct Photonhbt {
return o;
}

struct PairSep {
float dEtaPP{0}, dPhiPP{0}, dEtaNN{0}, dPhiNN{0};
float u{999.f};
};
PairSep computePairSep(float etaP1, float phiP1, float etaN1, float phiN1,
float etaP2, float phiP2, float etaN2, float phiN2) const
{
PairSep s;
s.dEtaPP = etaP1 - etaP2;
s.dPhiPP = RecoDecay::constrainAngle(phiP1 - phiP2, -o2::constants::math::PI);
s.dEtaNN = etaN1 - etaN2;
s.dPhiNN = RecoDecay::constrainAngle(phiN1 - phiN2, -o2::constants::math::PI);
const float sE = pairsep.cfgUSigEta.value, sP = pairsep.cfgUSigPhi.value;
s.u = std::min(std::hypot(s.dEtaPP / sE, s.dPhiPP / sP),
std::hypot(s.dEtaNN / sE, s.dPhiNN / sP));
return s;
}

struct TruthGamma {
int id = -1, posId = -1, negId = -1;
float eta = 0.f, phi = 0.f, pt = 0.f;
Expand All @@ -598,6 +623,21 @@ struct Photonhbt {
float alphaTrue = 0.f;
};

struct PhotonWithLegs {
float fPt{0}, fEta{0}, fPhi{0};
float fVx{0}, fVy{0}, fVz{0};
std::array<float, 2> fLegPt{}, fLegEta{}, fLegPhi{}; // [0] = e+, [1] = e-
float pt() const { return fPt; }
float eta() const { return fEta; }
float phi() const { return fPhi; }
float vx() const { return fVx; }
float vy() const { return fVy; }
float vz() const { return fVz; }
float legEta(int i) const { return fLegEta[i]; }
float legPhi(int i) const { return fLegPhi[i]; }
float legPt(int i) const { return fLegPt[i]; }
};

std::map<std::tuple<int, int, int, int>, std::deque<std::vector<TruthGamma>>> truthGammaPool;

void addSinglePhotonQAHistogramsForStep(const std::string& path)
Expand Down Expand Up @@ -629,6 +669,19 @@ struct Photonhbt {
fRegistry.fill(HIST(base) + HIST("hSparseDeltaRDeltaZQinv"), obs.deltaR, obs.deltaZ, obs.qinv);
}

static constexpr std::string_view pairsep_types[2] = {"PairSep/same/", "PairSep/mix/"};
template <int ev_id>
inline void fillPairSep(PairSep const& s, PairQAObservables const& obs)
{
if (!pairsep.cfgDoPairSepQA.value || obs.qinv >= cfgMaxQinvForMCQA.value)
return;
fRegistry.fill(HIST(pairsep_types[ev_id]) + HIST("hDEtaDPhiKt_PP"), s.dEtaPP, s.dPhiPP, obs.kt);
fRegistry.fill(HIST(pairsep_types[ev_id]) + HIST("hDEtaDPhiKt_NN"), s.dEtaNN, s.dPhiNN, obs.kt);
fRegistry.fill(HIST(pairsep_types[ev_id]) + HIST("hU_Qinv_Kt"), s.u, obs.qinv, obs.kt);
fRegistry.fill(HIST(pairsep_types[ev_id]) + HIST("hU_Qinv_dE"), s.u, obs.qinv,
std::fabs(obs.v1.E() - obs.v2.E()));
}

template <int ev_id>
inline void fillFullRangeDeltaRCosOA(float qinv, float drOverCosOA)
{
Expand Down Expand Up @@ -779,6 +832,18 @@ struct Photonhbt {

fRegistry.addClone("Pair/same/", "Pair/mix/");

const std::array<std::string, 2> pairSepDirs = {"PairSep/same/", "PairSep/mix/"};
for (const auto& sm : pairSepDirs) {
fRegistry.add((sm + "hDEtaDPhiKt_PP").c_str(), "e^{+}e^{+};#Delta#eta;#Delta#phi (rad);k_{T} (GeV/c)",
kTHnSparseF, {{160, -0.4f, 0.4f}, {160, -0.4f, 0.4f}, axisKt}, true);
fRegistry.add((sm + "hDEtaDPhiKt_NN").c_str(), "e^{-}e^{-};#Delta#eta;#Delta#phi (rad);k_{T} (GeV/c)",
kTHnSparseF, {{160, -0.4f, 0.4f}, {160, -0.4f, 0.4f}, axisKt}, true);
fRegistry.add((sm + "hU_Qinv_Kt").c_str(), ";u;q_{inv} (GeV/c);k_{T} (GeV/c)",
kTHnSparseF, {{100, 0.f, 10.f}, {50, 0.f, 0.5f}, axisKt}, true);
fRegistry.add((sm + "hU_Qinv_dE").c_str(), ";u;q_{inv} (GeV/c);#DeltaE (GeV)",
kTHnSparseF, {{50, 0.f, 10.f}, {50, 0.f, 0.5f}, {40, 0.f, 2.f}}, true);
}

fRegistry.add("Pair/same/EtaTopology/hSparse_DEtaDPhi_kT_sameSideV0_sameSideLegs",
"both V0 same #eta-side, all legs same side;"
"#Delta#eta;#Delta#phi (rad);k_{T} (GeV/c)",
Expand Down Expand Up @@ -967,14 +1032,6 @@ struct Photonhbt {
} else {
fRegistryPairMC.fill(HIST(mcDir) + HIST("CF_1D"), hbtanalysis.cfgUseLCMS ? qabs_lcms : qinv, kt, weight);
}
float deta_pair = v1.Eta() - v2.Eta();
float dphi_pair = v1.Phi() - v2.Phi();
dphi_pair = RecoDecay::constrainAngle(dphi_pair, -o2::constants::math::PI);
if constexpr (ev_id == 0) {
fRegistry.fill(HIST("Pair/same/hSparse_DEtaDPhi_qinv_kT"), deta_pair, dphi_pair, qinv, kt, weight);
} else {
fRegistry.fill(HIST("Pair/mix/hSparse_DEtaDPhi_qinv_kT"), deta_pair, dphi_pair, qinv, kt, weight);
}
}

template <typename TG1, typename TG2>
Expand Down Expand Up @@ -1805,6 +1862,9 @@ struct Photonhbt {
continue;
idsAfterDR.insert(g1.globalIndex());
idsAfterDR.insert(g2.globalIndex());
fillPairSep<0>(computePairSep(pos1.eta(), pos1.phi(), ele1.eta(), ele1.phi(),
pos2.eta(), pos2.phi(), ele2.eta(), ele2.phi()),
obs);
if (doQA)
fillPairQAStep<0, 1>(obs, centForQA, occupancy);
if (!passRZCut(obs.deltaR, obs.deltaZ))
Expand Down Expand Up @@ -1833,15 +1893,23 @@ struct Photonhbt {
pos1.eta(), ele1.eta(),
pos2.eta(), ele2.eta()),
obs.deta, obs.dphi, obs.kt);
auto addToPool = [&](auto const& g) {
auto addToPool = [&](auto const& g, auto const& pos, auto const& ele) {
if (usedPhotonIdsPerCol.insert(g.globalIndex()).second) {
EMPair gtmp(g.pt(), g.eta(), g.phi(), 0.f);
gtmp.setConversionPointXYZ(g.vx(), g.vy(), g.vz());
emh1->AddTrackToEventPool(keyDFCollision, gtmp);
PhotonWithLegs p;
p.fPt = g.pt();
p.fEta = g.eta();
p.fPhi = g.phi();
p.fVx = g.vx();
p.fVy = g.vy();
p.fVz = g.vz();
p.fLegPt = {static_cast<float>(pos.pt()), static_cast<float>(ele.pt())};
p.fLegEta = {static_cast<float>(pos.eta()), static_cast<float>(ele.eta())};
p.fLegPhi = {static_cast<float>(pos.phi()), static_cast<float>(ele.phi())};
emh1->AddTrackToEventPool(keyDFCollision, p);
}
};
addToPool(g1);
addToPool(g2);
addToPool(g1, pos1, ele1);
addToPool(g2, pos2, ele2);
}
if (qaflags.doSinglePhotonQa) {
for (const auto& g : photons1Coll) {
Expand Down Expand Up @@ -1888,6 +1956,9 @@ struct Photonhbt {
fRegistry.fill(HIST("Pair/mix/hDeltaRCosOA"), obs.drOverCosOA);
if (obs.drOverCosOA < ggpaircuts.cfgMinDRCosOA)
continue;
fillPairSep<1>(computePairSep(g1.legEta(0), g1.legPhi(0), g1.legEta(1), g1.legPhi(1),
g2.legEta(0), g2.legPhi(0), g2.legEta(1), g2.legPhi(1)),
obs);
if (doQA)
fillPairQAStep<1, 1>(obs, centForQA, occupancy);
if (!passRZCut(obs.deltaR, obs.deltaZ))
Expand Down Expand Up @@ -2060,15 +2131,23 @@ struct Photonhbt {
}
}

auto addToPool = [&](auto const& g) {
auto addToPool = [&](auto const& g, auto const& pos, auto const& ele) {
if (usedPhotonIdsPerCol.insert(g.globalIndex()).second) {
EMPair gtmp(g.pt(), g.eta(), g.phi(), 0.f);
gtmp.setConversionPointXYZ(g.vx(), g.vy(), g.vz());
emh1->AddTrackToEventPool(keyDFCollision, gtmp);
PhotonWithLegs p;
p.fPt = g.pt();
p.fEta = g.eta();
p.fPhi = g.phi();
p.fVx = g.vx();
p.fVy = g.vy();
p.fVz = g.vz();
p.fLegPt = {static_cast<float>(pos.pt()), static_cast<float>(ele.pt())};
p.fLegEta = {static_cast<float>(pos.eta()), static_cast<float>(ele.eta())};
p.fLegPhi = {static_cast<float>(pos.phi()), static_cast<float>(ele.phi())};
emh1->AddTrackToEventPool(keyDFCollision, p);
}
};
addToPool(g1);
addToPool(g2);
addToPool(g1, pos1, ele1);
addToPool(g2, pos2, ele2);
}
if (qaflags.doSinglePhotonQa) {
for (const auto& g : photonsColl) {
Expand Down Expand Up @@ -2251,10 +2330,10 @@ struct Photonhbt {
continue;
}
const auto d = emmcParticles.iteratorAt(dId);
if (d.pdgCode() == kElectron) {
if (d.pdgCode() == kPositron) {
posId = dId;
rTrue = std::sqrt(d.vx() * d.vx() + d.vy() * d.vy());
} else if (d.pdgCode() == kPositron) {
} else if (d.pdgCode() == kElectron) {
negId = dId;
}
}
Expand Down Expand Up @@ -2512,11 +2591,11 @@ struct Photonhbt {
poolBin.pop_front();
}
} // end cfgDoTruthMix
} // end collision loop
} // end runTruthEfficiency
} // end collision loop
} // end runTruthEfficiency

using MyEMH = o2::aod::pwgem::dilepton::utils::EventMixingHandler<
std::tuple<int, int, int, int>, std::pair<int, int>, EMPair>;
std::tuple<int, int, int, int>, std::pair<int, int>, PhotonWithLegs>;

MyEMH* emh1 = nullptr;
MyEMH* emh2 = nullptr;
Expand Down
Loading