Skip to content

[PWGLF] Add minimum separation variable among daughter kaons #13705

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

Merged
skundu692 merged 1 commit into from
Nov 6, 2025
Merged

[PWGLF] Add minimum separation variable among daughter kaons #13705

Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
128 changes: 38 additions & 90 deletions PWGLF/Tasks/Resonances/doublephimeson.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -28,13 +28,13 @@
#include <Math/GenVector/Boost.h>
#include <Math/Vector3D.h>
#include <Math/Vector4D.h>
#include <TLorentzVector.h>

Check failure on line 31 in PWGLF/Tasks/Resonances/doublephimeson.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector] 

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
#include <TMath.h>
#include <TVector2.h>

#include <fairlogger/Logger.h>

#include <iostream>

Check failure on line 37 in PWGLF/Tasks/Resonances/doublephimeson.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[include-iostream] 

Do not include iostream. Use O2 logging instead.
#include <iterator>
#include <string>
#include <vector>
Expand Down Expand Up @@ -106,15 +106,15 @@
const AxisSpec thnAxisCosTheta{configThnAxisCosTheta, "cos #theta"};
const AxisSpec thnAxisNumPhi{configThnAxisNumPhi, "Number of phi meson"};

histos.add("SEMassUnlike", "SEMassUnlike", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisdeltapt, thnAxisPt, thnAxisDeltaR, thnAxisInvMassDeltaPhi, thnAxisNumPhi});
histos.add("SEMassUnlike", "SEMassUnlike", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisDeltaR, thnAxisPt, thnAxisDeltaR, thnAxisInvMassDeltaPhi, thnAxisNumPhi});
// histos.add("SEMassLike", "SEMassLike", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisDeltaR, thnAxisInvMassPhi, thnAxisInvMassPhi, thnAxisNumPhi});
histos.add("MEMassUnlike", "MEMassUnlike", HistType::kTHnSparseF, {thnAxisInvMass, thnAxisPt, thnAxisDeltaR, thnAxisInvMassDeltaPhi});
}

// get kstar
TLorentzVector trackSum, PartOneCMS, PartTwoCMS, trackRelK;

Check failure on line 115 in PWGLF/Tasks/Resonances/doublephimeson.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector] 

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
float getkstar(const TLorentzVector part1,

Check failure on line 116 in PWGLF/Tasks/Resonances/doublephimeson.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector] 

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
const TLorentzVector part2)

Check failure on line 117 in PWGLF/Tasks/Resonances/doublephimeson.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector] 

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
{
// const TLorentzVector trackSum = part1 + part2;
trackSum = part1 + part2;
Expand Down Expand Up @@ -148,8 +148,8 @@
return angle;
}

float deepangle(const TLorentzVector candidate1,

Check failure on line 151 in PWGLF/Tasks/Resonances/doublephimeson.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector] 

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
const TLorentzVector candidate2)

Check failure on line 152 in PWGLF/Tasks/Resonances/doublephimeson.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector] 

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
{
double pt1, pt2, pz1, pz2, p1, p2, angle;
pt1 = candidate1.Pt();
Expand All @@ -163,10 +163,10 @@
}

// get cosTheta
TLorentzVector daughterCMS;

Check failure on line 166 in PWGLF/Tasks/Resonances/doublephimeson.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector] 

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
ROOT::Math::XYZVector threeVecDauCM, threeVecMother;
float getCosTheta(const TLorentzVector mother,

Check failure on line 168 in PWGLF/Tasks/Resonances/doublephimeson.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector] 

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
const TLorentzVector daughter)

Check failure on line 169 in PWGLF/Tasks/Resonances/doublephimeson.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[root/lorentz-vector] 

Do not use the TLorentzVector legacy class. Use std::array with RecoDecay methods or the ROOT::Math::LorentzVector template instead.
{
threeVecMother = mother.Vect();
const float beta = mother.Beta();
Expand Down Expand Up @@ -762,6 +762,7 @@
}
if (phimult < 2)
return;

// --- helpers ---
constexpr double mPhiPDG = 1.019461; // GeV/c^2

Expand All @@ -776,20 +777,35 @@
return std::sqrt(dphi * dphi + deta * deta);
};

// anti-merging: same-sign kaons must be separated
const auto passDaughterDR = [&](const ROOT::Math::PtEtaPhiMVector& kplusA,
const ROOT::Math::PtEtaPhiMVector& kplusB,
const ROOT::Math::PtEtaPhiMVector& kminusA,
const ROOT::Math::PtEtaPhiMVector& kminusB,
double thr) {
// minimum ΔR among all kaons in the candidate (4 kaons → 6 combinations)
const auto minKaonDeltaR = [&](const ROOT::Math::PtEtaPhiMVector& kplusA,
const ROOT::Math::PtEtaPhiMVector& kplusB,
const ROOT::Math::PtEtaPhiMVector& kminusA,
const ROOT::Math::PtEtaPhiMVector& kminusB) {
// same-sign first (keep your QA histos)
const double dRkplus = deltaR(kplusA.Phi(), kplusA.Eta(), kplusB.Phi(), kplusB.Eta());
const double dRkminus = deltaR(kminusA.Phi(), kminusA.Eta(), kminusB.Phi(), kminusB.Eta());
histos.fill(HIST("hDeltaRkaonplus"), dRkplus);
histos.fill(HIST("hDeltaRkaonminus"), dRkminus);
return (dRkplus > thr) && (dRkminus > thr);

// all other combinations
const double dR_k1p_k1m = deltaR(kplusA.Phi(), kplusA.Eta(), kminusA.Phi(), kminusA.Eta());
const double dR_k1p_k2m = deltaR(kplusA.Phi(), kplusA.Eta(), kminusB.Phi(), kminusB.Eta());
const double dR_k2p_k1m = deltaR(kplusB.Phi(), kplusB.Eta(), kminusA.Phi(), kminusA.Eta());
const double dR_k2p_k2m = deltaR(kplusB.Phi(), kplusB.Eta(), kminusB.Phi(), kminusB.Eta());

double minDR = dRkplus;
minDR = std::min(minDR, dRkminus);
minDR = std::min(minDR, dR_k1p_k1m);
minDR = std::min(minDR, dR_k1p_k2m);
minDR = std::min(minDR, dR_k2p_k1m);
minDR = std::min(minDR, dR_k2p_k2m);
return minDR;
};

// --- collect candidates once ---
std::vector<ROOT::Math::PtEtaPhiMVector> pairV, phi1V, phi2V, kplus1V, kplus2V, kminus1V, kminus2V;
std::vector<double> minDRV; // store minimum ΔR for each pair

for (auto const& t1 : phitracks) {
const double kplus1pt = std::hypot(t1.phid1Px(), t1.phid1Py());
Expand Down Expand Up @@ -858,13 +874,12 @@
if (pair.M() < minExoticMass || pair.M() > maxExoticMass)
continue;

// daughter anti-merging
// daughter ΔR QA and minΔR (NO CUT anymore)
ROOT::Math::PtEtaPhiMVector k1pV(k1p.Pt(), k1p.Eta(), k1p.Phi(), 0.493);
ROOT::Math::PtEtaPhiMVector k1mV(k1m.Pt(), k1m.Eta(), k1m.Phi(), 0.493);
ROOT::Math::PtEtaPhiMVector k2pV(k2p.Pt(), k2p.Eta(), k2p.Phi(), 0.493);
ROOT::Math::PtEtaPhiMVector k2mV(k2m.Pt(), k2m.Eta(), k2m.Phi(), 0.493);
if (!passDaughterDR(k1pV, k2pV, k1mV, k2mV, daughterDeltaR))
continue;
const double minDR = minKaonDeltaR(k1pV, k2pV, k1mV, k2mV);

// store for one-pass fill
pairV.emplace_back(pair.Pt(), pair.Eta(), pair.Phi(), pair.M());
Expand All @@ -874,6 +889,7 @@
kminus1V.emplace_back(k1m.Pt(), k1m.Eta(), k1m.Phi(), 0.493);
kplus2V.emplace_back(k2p.Pt(), k2p.Eta(), k2p.Phi(), 0.493);
kminus2V.emplace_back(k2m.Pt(), k2m.Eta(), k2m.Phi(), 0.493);
minDRV.emplace_back(minDR); // per-candidate minimum ΔR of kaons
}
}

Expand All @@ -886,91 +902,23 @@
p1.SetPtEtaPhiM(phi1V[i].Pt(), phi1V[i].Eta(), phi1V[i].Phi(), phi1V[i].M());
p2.SetPtEtaPhiM(phi2V[i].Pt(), phi2V[i].Eta(), phi2V[i].Phi(), phi2V[i].M());
pair.SetPtEtaPhiM(pairV[i].Pt(), pairV[i].Eta(), pairV[i].Phi(), pairV[i].M());

const double dM = deltaMPhi(p1.M(), p2.M());
const double M = pair.M();
// const double pT = pair.Pt();
const double dR = deltaR(p1.Phi(), p1.Eta(), p2.Phi(), p2.Eta());
histos.fill(HIST("SEMassUnlike"), M, std::abs(p1.Pt() - p2.Pt()) / pair.Pt(), pair.Pt(), dR, dM, pairV.size());
const double minDR = minDRV[i];

// NOTE: second axis is now minΔR(all kaons), ΔpT/pT has been removed
histos.fill(HIST("SEMassUnlike"),
M,
minDR,
pair.Pt(),
dR,
dM,
pairV.size());
}
}
PROCESS_SWITCH(doublephimeson, processopti3, "Process Optimized same event", false);

SliceCache cache;
using BinningTypeVertexContributor = ColumnBinningPolicy<aod::collision::PosZ, aod::collision::NumContrib>;
void processMixedEvent(aod::RedPhiEvents& collisions, aod::PhiTracks& phitracks)
{
auto tracksTuple = std::make_tuple(phitracks);
BinningTypeVertexContributor binningOnPositions{{CfgVtxBins, CfgMultBins}, true};
SameKindPair<aod::RedPhiEvents, aod::PhiTracks, BinningTypeVertexContributor> pair{binningOnPositions, nEvtMixing, -1, collisions, tracksTuple, &cache};

for (auto& [collision1, tracks1, collision2, tracks2] : pair) {
if (collision1.index() == collision2.index()) {
continue;
}
for (auto& [phitrackd1, phitrackd2] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(tracks1, tracks2))) {
if (phitrackd1.phiMass() < minPhiMass1 || phitrackd1.phiMass() > maxPhiMass1) {
continue;
}
if (phitrackd2.phiMass() < minPhiMass2 || phitrackd2.phiMass() > maxPhiMass2) {
continue;
}
auto kaonplusd1pt = TMath::Sqrt(phitrackd1.phid1Px() * phitrackd1.phid1Px() + phitrackd1.phid1Py() * phitrackd1.phid1Py());
auto kaonminusd1pt = TMath::Sqrt(phitrackd1.phid2Px() * phitrackd1.phid2Px() + phitrackd1.phid2Py() * phitrackd1.phid2Py());
auto kaonplusd2pt = TMath::Sqrt(phitrackd2.phid1Px() * phitrackd2.phid1Px() + phitrackd2.phid1Py() * phitrackd2.phid1Py());
auto kaonminusd2pt = TMath::Sqrt(phitrackd2.phid2Px() * phitrackd2.phid2Px() + phitrackd2.phid2Py() * phitrackd2.phid2Py());
if (kaonplusd1pt > maxKaonPt) {
continue;
}
if (kaonminusd1pt > maxKaonPt) {
continue;
}
if (kaonplusd2pt > maxKaonPt) {
continue;
}
if (kaonminusd2pt > maxKaonPt) {
continue;
}
if (!selectionPID(phitrackd1.phid1TPC(), phitrackd1.phid1TOF(), phitrackd1.phid1TOFHit(), strategyPID1, kaonplusd1pt)) {
continue;
}
if (!selectionPID(phitrackd1.phid2TPC(), phitrackd1.phid2TOF(), phitrackd1.phid2TOFHit(), strategyPID1, kaonminusd1pt)) {
continue;
}
Phid1.SetXYZM(phitrackd1.phiPx(), phitrackd1.phiPy(), phitrackd1.phiPz(), phitrackd1.phiMass());
if (!selectionPID(phitrackd2.phid1TPC(), phitrackd2.phid1TOF(), phitrackd2.phid1TOFHit(), strategyPID2, kaonplusd2pt)) {
continue;
}
if (!selectionPID(phitrackd2.phid2TPC(), phitrackd2.phid2TOF(), phitrackd2.phid2TOFHit(), strategyPID2, kaonminusd2pt)) {
continue;
}
Phid2.SetXYZM(phitrackd2.phiPx(), phitrackd2.phiPy(), phitrackd2.phiPz(), phitrackd2.phiMass());
exotic = Phid1 + Phid2;

Phi1kaonplus.SetXYZM(phitrackd1.phid1Px(), phitrackd1.phid1Py(), phitrackd1.phid1Pz(), 0.493);
Phi1kaonminus.SetXYZM(phitrackd1.phid2Px(), phitrackd1.phid2Py(), phitrackd1.phid2Pz(), 0.493);
Phi2kaonplus.SetXYZM(phitrackd2.phid1Px(), phitrackd2.phid1Py(), phitrackd2.phid1Pz(), 0.493);
Phi2kaonminus.SetXYZM(phitrackd2.phid2Px(), phitrackd2.phid2Py(), phitrackd2.phid2Pz(), 0.493);
auto deltaRd1 = TMath::Sqrt(TMath::Power(Phi1kaonplus.Phi() - Phi2kaonplus.Phi(), 2.0) + TMath::Power(Phi1kaonplus.Eta() - Phi2kaonplus.Eta(), 2.0));
auto deltaRd2 = TMath::Sqrt(TMath::Power(Phi1kaonminus.Phi() - Phi2kaonminus.Phi(), 2.0) + TMath::Power(Phi1kaonminus.Eta() - Phi2kaonminus.Eta(), 2.0));
auto deltaR = TMath::Sqrt(TMath::Power(Phid1.Phi() - Phid2.Phi(), 2.0) + TMath::Power(Phid1.Eta() - Phid2.Eta(), 2.0));
// auto costheta = (Phid1.Px() * Phid2.Px() + Phid1.Py() * Phid2.Py() + Phid1.Pz() * Phid2.Pz()) / (Phid1.P() * Phid2.P());
auto deltam = TMath::Sqrt(TMath::Power(Phid1.M() - 1.0192, 2.0) + TMath::Power(Phid2.M() - 1.0192, 2.0));
if (deltaRd1 < daughterDeltaR) {
continue;
}
if (deltaRd2 < daughterDeltaR) {
continue;
}
if (!isDeep) {
histos.fill(HIST("MEMassUnlike"), exotic.M(), exotic.Pt(), deltaR, deltam);
}
if (isDeep) {
histos.fill(HIST("MEMassUnlike"), exotic.M(), exotic.Pt(), deltaR, deltam);
}
}
}
}
PROCESS_SWITCH(doublephimeson, processMixedEvent, "Process EventMixing for combinatorial background", false);
};

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