Skip to content

Commit 7127f73

Browse files
mhemmer-cernmhemmer-cern
authored
[PWGEM] PM: Add TooCloseV0 cut to AreSelectedRunning function call: (#15442)
1 parent 12546fe commit 7127f73

File tree

3 files changed

+246
-26
lines changed

3 files changed

+246
-26
lines changed

‎PWGEM/PhotonMeson/Core/V0PhotonCut.cxx‎

Lines changed: 22 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -115,6 +115,28 @@ void V0PhotonCut::RejectITSib(bool flag)
115115
mRejectITSib = flag;
116116
LOG(info) << "V0 Photon Cut, reject photon on ITSib: " << mRejectITSib;
117117
}
118+
119+
void V0PhotonCut::setTooCloseType(V0PhotonCut::TooCloseCuts type)
120+
{
121+
mTooCloseType = type;
122+
LOG(info) << "V0 Photon Cut, TooCloseV0 cut type: " << static_cast<uint>(mTooCloseType);
123+
}
124+
void V0PhotonCut::setMinV0DistSquared(float value)
125+
{
126+
mMinV0DistSquared = value;
127+
LOG(info) << "V0 Photon Cut, min V0 distance squared: " << mMinV0DistSquared;
128+
}
129+
void V0PhotonCut::setDeltaR(float value)
130+
{
131+
mDeltaR = value;
132+
LOG(info) << "V0 Photon Cut, delta R for too close V0: " << mDeltaR;
133+
}
134+
void V0PhotonCut::setMinOpeningAngle(float value)
135+
{
136+
mMinOpeningAngle = value;
137+
LOG(info) << "V0 Photon Cut, min opening angle for too close V0: " << mMinOpeningAngle;
138+
}
139+
118140
void V0PhotonCut::SetTPCNsigmaElRange(float min, float max)
119141
{
120142
mMinTPCNsigmaEl = min;

‎PWGEM/PhotonMeson/Core/V0PhotonCut.h‎

Lines changed: 186 additions & 24 deletions
Original file line numberDiff line numberDiff line change
@@ -27,6 +27,9 @@
2727
#include <Framework/Array2D.h>
2828
#include <Framework/HistogramRegistry.h>
2929

30+
#include <Math/Vector3D.h> // IWYU pragma: keep
31+
#include <Math/Vector3Dfwd.h>
32+
#include <Math/VectorUtil.h>
3033
#include <TMath.h>
3134
#include <TNamed.h>
3235

@@ -183,6 +186,7 @@ class V0PhotonCut : public TNamed
183186
kRZLine,
184187
kOnWwireIB,
185188
kOnWwireOB,
189+
kIsTooClose,
186190
// leg cut
187191
kTrackPtRange,
188192
kTrackEtaRange,
@@ -206,6 +210,12 @@ class V0PhotonCut : public TNamed
206210
kNCuts
207211
};
208212

213+
enum class TooCloseCuts : uint8_t {
214+
kNoCut = 0,
215+
kDistance3D = 1,
216+
kRadAndAngle = 2
217+
};
218+
209219
/// \brief add histograms to registry
210220
/// \param fRegistry pointer to histogram registry
211221
void addQAHistograms(o2::framework::HistogramRegistry* fRegistry = nullptr) const
@@ -268,26 +278,27 @@ class V0PhotonCut : public TNamed
268278
hPhotonQualityCuts->GetXaxis()->SetBinLabel(12, "RZ_{line}");
269279
hPhotonQualityCuts->GetXaxis()->SetBinLabel(13, "Wire_{IB}");
270280
hPhotonQualityCuts->GetXaxis()->SetBinLabel(14, "Wire_{OB}");
271-
hPhotonQualityCuts->GetXaxis()->SetBinLabel(15, "#it{p}_{T,leg}");
272-
hPhotonQualityCuts->GetXaxis()->SetBinLabel(16, "#it{#eta}_{leg}");
273-
hPhotonQualityCuts->GetXaxis()->SetBinLabel(17, "#it{N}_{cl,TPC}");
274-
hPhotonQualityCuts->GetXaxis()->SetBinLabel(18, "#it{N}_{cr,TPC}");
275-
hPhotonQualityCuts->GetXaxis()->SetBinLabel(19, "#it{N}_{cr,TPC}/#it{N}_{cl,TPC}");
276-
hPhotonQualityCuts->GetXaxis()->SetBinLabel(20, "FracSharedCl");
277-
hPhotonQualityCuts->GetXaxis()->SetBinLabel(21, "#chi^{2}_{TPC}/NDF");
278-
hPhotonQualityCuts->GetXaxis()->SetBinLabel(22, "#it{N#sigma}_{e,TPC}");
279-
hPhotonQualityCuts->GetXaxis()->SetBinLabel(23, "#it{N#sigma}_{#pi,TPC}");
280-
hPhotonQualityCuts->GetXaxis()->SetBinLabel(24, "DCA_{xy}");
281-
hPhotonQualityCuts->GetXaxis()->SetBinLabel(25, "DCA_{z}");
282-
hPhotonQualityCuts->GetXaxis()->SetBinLabel(26, "#it{N}_{cl,ITS}");
283-
hPhotonQualityCuts->GetXaxis()->SetBinLabel(27, "#chi^{2}_{ITS}/NDF");
284-
hPhotonQualityCuts->GetXaxis()->SetBinLabel(28, "size_{ITS}");
285-
hPhotonQualityCuts->GetXaxis()->SetBinLabel(29, "ITSTPC");
286-
hPhotonQualityCuts->GetXaxis()->SetBinLabel(30, "ITSOnly");
287-
hPhotonQualityCuts->GetXaxis()->SetBinLabel(31, "TPCOnly");
288-
hPhotonQualityCuts->GetXaxis()->SetBinLabel(32, "TPCTRD");
289-
hPhotonQualityCuts->GetXaxis()->SetBinLabel(33, "TPCTOF");
290-
hPhotonQualityCuts->GetXaxis()->SetBinLabel(34, "Out");
281+
hPhotonQualityCuts->GetXaxis()->SetBinLabel(15, "IsTooClose");
282+
hPhotonQualityCuts->GetXaxis()->SetBinLabel(16, "#it{p}_{T,leg}");
283+
hPhotonQualityCuts->GetXaxis()->SetBinLabel(17, "#it{#eta}_{leg}");
284+
hPhotonQualityCuts->GetXaxis()->SetBinLabel(18, "#it{N}_{cl,TPC}");
285+
hPhotonQualityCuts->GetXaxis()->SetBinLabel(19, "#it{N}_{cr,TPC}");
286+
hPhotonQualityCuts->GetXaxis()->SetBinLabel(20, "#it{N}_{cr,TPC}/#it{N}_{cl,TPC}");
287+
hPhotonQualityCuts->GetXaxis()->SetBinLabel(21, "FracSharedCl");
288+
hPhotonQualityCuts->GetXaxis()->SetBinLabel(22, "#chi^{2}_{TPC}/NDF");
289+
hPhotonQualityCuts->GetXaxis()->SetBinLabel(23, "#it{N#sigma}_{e,TPC}");
290+
hPhotonQualityCuts->GetXaxis()->SetBinLabel(24, "#it{N#sigma}_{#pi,TPC}");
291+
hPhotonQualityCuts->GetXaxis()->SetBinLabel(25, "DCA_{xy}");
292+
hPhotonQualityCuts->GetXaxis()->SetBinLabel(26, "DCA_{z}");
293+
hPhotonQualityCuts->GetXaxis()->SetBinLabel(27, "#it{N}_{cl,ITS}");
294+
hPhotonQualityCuts->GetXaxis()->SetBinLabel(28, "#chi^{2}_{ITS}/NDF");
295+
hPhotonQualityCuts->GetXaxis()->SetBinLabel(29, "size_{ITS}");
296+
hPhotonQualityCuts->GetXaxis()->SetBinLabel(30, "ITSTPC");
297+
hPhotonQualityCuts->GetXaxis()->SetBinLabel(31, "ITSOnly");
298+
hPhotonQualityCuts->GetXaxis()->SetBinLabel(32, "TPCOnly");
299+
hPhotonQualityCuts->GetXaxis()->SetBinLabel(33, "TPCTRD");
300+
hPhotonQualityCuts->GetXaxis()->SetBinLabel(34, "TPCTOF");
301+
hPhotonQualityCuts->GetXaxis()->SetBinLabel(35, "Out");
291302
}
292303
}
293304

@@ -351,6 +362,119 @@ class V0PhotonCut : public TNamed
351362
fRegistry->fill(HIST("QA/V0Photon/after/Neg/hTPCHits"), ele.tpcNClsFound(), ele.tpcNClsCrossedRows());
352363
}
353364

365+
/// \brief creates a mask for the V0s if they are too close to another V0 and have higher chi^2
366+
/// \param v0s V0 table
367+
template <o2::soa::is_table TV0>
368+
void createCloseV0CutMask(TV0 const& v0s) const
369+
{
370+
const bool useDistance3D = (mTooCloseType == TooCloseCuts::kDistance3D);
371+
const float windowWidth = useDistance3D ? std::sqrt(mMinV0DistSquared) : mDeltaR;
372+
const float cosMinAngle = std::cos(mMinOpeningAngle);
373+
374+
int tableSize = v0s.size();
375+
std::vector<uint8_t> rejectMask(tableSize, 0);
376+
377+
if (mTooCloseType == TooCloseCuts::kNoCut) {
378+
mRejectMask = rejectMask;
379+
return;
380+
}
381+
382+
auto currentV0Iter = v0s.begin();
383+
auto otherV0Iter = v0s.begin();
384+
385+
int groupStart = 0;
386+
while (groupStart < tableSize) {
387+
388+
// --- find the end of this collision's group ---
389+
int currentCollisionId = v0s.iteratorAt(groupStart).collisionId();
390+
int groupEnd = groupStart + 1;
391+
while (groupEnd < tableSize && v0s.iteratorAt(groupEnd).collisionId() == currentCollisionId) {
392+
groupEnd++;
393+
}
394+
395+
int groupSize = groupEnd - groupStart;
396+
397+
std::vector<std::pair<int, float>> indexedRadii(groupSize);
398+
for (int k = 0; k < groupSize; k++) {
399+
currentV0Iter.setCursor(groupStart + k);
400+
indexedRadii[k] = {groupStart + k, currentV0Iter.v0radius()};
401+
}
402+
std::sort(indexedRadii.begin(), indexedRadii.end(), [](const auto& a, const auto& b) {
403+
return a.second < b.second;
404+
});
405+
// extract sorted indices and pre-sorted radii
406+
std::vector<int> sortedIndices(groupSize);
407+
std::vector<float> sortedRadii(groupSize);
408+
for (int k = 0; k < groupSize; k++) {
409+
sortedIndices[k] = indexedRadii[k].first;
410+
sortedRadii[k] = indexedRadii[k].second;
411+
}
412+
413+
// --- sliding window within this group ---
414+
int windowStart = 0; // reset per group
415+
for (int i = 0; i < groupSize; i++) {
416+
417+
float currentRadius = sortedRadii[i];
418+
while (windowStart < groupSize && sortedRadii[windowStart] < currentRadius - windowWidth) {
419+
windowStart++;
420+
}
421+
422+
currentV0Iter.setCursor(sortedIndices[i]);
423+
424+
float vx1 = currentV0Iter.vx();
425+
float vy1 = currentV0Iter.vy();
426+
float vz1 = currentV0Iter.vz();
427+
428+
float px1 = currentV0Iter.px();
429+
float py1 = currentV0Iter.py();
430+
float pz1 = currentV0Iter.pz();
431+
float chi2I = currentV0Iter.chiSquareNDF();
432+
433+
for (int j = windowStart; j < groupSize; j++) {
434+
if (j == i) {
435+
continue;
436+
}
437+
if (sortedRadii[j] > currentRadius + windowWidth) {
438+
break;
439+
}
440+
441+
otherV0Iter.setCursor(sortedIndices[j]);
442+
443+
bool tooClose = false;
444+
if (useDistance3D) {
445+
float dx = vx1 - otherV0Iter.vx();
446+
float dy = vy1 - otherV0Iter.vy();
447+
float dz = vz1 - otherV0Iter.vz();
448+
float distSquared = dx * dx + dy * dy + dz * dz;
449+
tooClose = distSquared < mMinV0DistSquared;
450+
} else {
451+
float px2 = otherV0Iter.px(), py2 = otherV0Iter.py(), pz2 = otherV0Iter.pz();
452+
float dot = px1 * px2 + py1 * py2 + pz1 * pz2;
453+
float mag1 = px1 * px1 + py1 * py1 + pz1 * pz1;
454+
float mag2 = px2 * px2 + py2 * py2 + pz2 * pz2;
455+
float denom = std::sqrt(mag1 * mag2);
456+
if (denom > 0) {
457+
float cosAngle = dot / denom;
458+
cosAngle = std::clamp(cosAngle, -1.0f, 1.0f);
459+
tooClose = cosAngle > cosMinAngle;
460+
}
461+
}
462+
463+
if (tooClose) {
464+
if (chi2I > otherV0Iter.chiSquareNDF()) {
465+
rejectMask[sortedIndices[i]] = 1;
466+
467+
} else {
468+
rejectMask[sortedIndices[j]] = 1;
469+
}
470+
}
471+
}
472+
}
473+
groupStart = groupEnd;
474+
}
475+
mRejectMask = rejectMask;
476+
}
477+
354478
/// \brief check if given v0 photon survives all cuts
355479
/// \param flags EMBitFlags where results will be stored
356480
/// \param v0s v0 photon table to check
@@ -360,6 +484,8 @@ class V0PhotonCut : public TNamed
360484
if (v0s.size() <= 0) {
361485
return;
362486
}
487+
488+
createCloseV0CutMask(v0s);
363489
// auto legIter = legs.begin();
364490
// auto legEnd = legs.end();
365491
size_t iV0 = 0;
@@ -489,6 +615,13 @@ class V0PhotonCut : public TNamed
489615
}
490616
}
491617

618+
if (!IsSelectedV0(v0, V0PhotonCuts::kIsTooClose)) {
619+
if (doQA) {
620+
fRegistry->fill(HIST("QA/V0Photon/hPhotonQualityCuts"), static_cast<int>(V0PhotonCuts::kIsTooClose) + 1, v0Pt);
621+
}
622+
return false;
623+
}
624+
492625
for (const auto& track : {pos, ele}) {
493626
if (!IsSelectedTrack(track, V0PhotonCuts::kTrackPtRange)) {
494627
if (doQA) {
@@ -707,11 +840,24 @@ class V0PhotonCut : public TNamed
707840
case V0PhotonCuts::kAP:
708841
return std::pow(v0.alpha() / mMaxAlpha, 2) + std::pow(v0.qtarm() / mMaxQt, 2) < 1.0;
709842

710-
case V0PhotonCuts::kPsiPair:
711-
return true;
843+
// TODO: implement fully
844+
case V0PhotonCuts::kPsiPair: {
845+
if constexpr (requires { v0.psipair(); }) {
846+
// return (std::fabs(v0.psipair() < 0.18f * std::exp( -0.55f * v0.chiSquareNDF())));
847+
return true;
848+
} else {
849+
return true;
850+
}
851+
}
712852

713-
case V0PhotonCuts::kPhiV:
714-
return true;
853+
// TODO: implement fully
854+
case V0PhotonCuts::kPhiV: {
855+
if constexpr (requires { v0.phiv(); }) {
856+
return true;
857+
} else {
858+
return true;
859+
}
860+
}
715861

716862
case V0PhotonCuts::kRxy: {
717863
if (v0.v0radius() < mMinRxy || mMaxRxy < v0.v0radius()) {
@@ -775,6 +921,12 @@ class V0PhotonCut : public TNamed
775921
float dxy = std::sqrt(std::pow(x - x_exp, 2) + std::pow(y - y_exp, 2));
776922
return !(dxy > margin_xy);
777923
}
924+
case V0PhotonCuts::kIsTooClose: {
925+
if (mRejectMask.size() == 0) {
926+
return true;
927+
}
928+
return (mRejectMask[v0.globalIndex()] == 0);
929+
}
778930
default:
779931
return false;
780932
}
@@ -950,6 +1102,11 @@ class V0PhotonCut : public TNamed
9501102
void SetOnWwireOB(bool flag = false);
9511103
void RejectITSib(bool flag = false);
9521104

1105+
void setTooCloseType(V0PhotonCut::TooCloseCuts type);
1106+
void setMinV0DistSquared(float value);
1107+
void setDeltaR(float value);
1108+
void setMinOpeningAngle(float value);
1109+
9531110
void SetTrackPtRange(float minPt = 0.f, float maxPt = 1e10f);
9541111
void SetTrackEtaRange(float minEta = -1e10f, float maxEta = 1e10f);
9551112
void SetMinNClustersTPC(int minNClustersTPC);
@@ -1017,6 +1174,11 @@ class V0PhotonCut : public TNamed
10171174
bool mIsOnWwireIB{false};
10181175
bool mIsOnWwireOB{false};
10191176
bool mRejectITSib{false};
1177+
TooCloseCuts mTooCloseType{V0PhotonCut::TooCloseCuts::kRadAndAngle}; // for TooCloseV0Cut: either squared distance between conversion points OR opening angle and deltaR
1178+
float mMinV0DistSquared{1.}; // for TooCloseV0Cut: cut value when using squared distance between conversion points
1179+
float mDeltaR{6.}; // for TooCloseV0Cut: V0PhotonCut::TooCloseCuts::kRadAndAngle when deltaR < this -> compare chi2
1180+
float mMinOpeningAngle{0.02}; // for TooCloseV0Cut: V0PhotonCut::TooCloseCuts::kRadAndAngle when opening angle < this -> compare chi2
1181+
mutable std::vector<uint8_t> mRejectMask{};
10201182

10211183
// ML cuts
10221184
bool mApplyMlCuts{false};

0 commit comments

Comments
 (0)