diff --git a/PWGEM/PhotonMeson/Core/V0PhotonCut.cxx b/PWGEM/PhotonMeson/Core/V0PhotonCut.cxx index 5abc97c5d1c..0511e72d0c7 100644 --- a/PWGEM/PhotonMeson/Core/V0PhotonCut.cxx +++ b/PWGEM/PhotonMeson/Core/V0PhotonCut.cxx @@ -115,6 +115,28 @@ void V0PhotonCut::RejectITSib(bool flag) mRejectITSib = flag; LOG(info) << "V0 Photon Cut, reject photon on ITSib: " << mRejectITSib; } + +void V0PhotonCut::setTooCloseType(V0PhotonCut::TooCloseCuts type) +{ + mTooCloseType = type; + LOG(info) << "V0 Photon Cut, TooCloseV0 cut type: " << static_cast(mTooCloseType); +} +void V0PhotonCut::setMinV0DistSquared(float value) +{ + mMinV0DistSquared = value; + LOG(info) << "V0 Photon Cut, min V0 distance squared: " << mMinV0DistSquared; +} +void V0PhotonCut::setDeltaR(float value) +{ + mDeltaR = value; + LOG(info) << "V0 Photon Cut, delta R for too close V0: " << mDeltaR; +} +void V0PhotonCut::setMinOpeningAngle(float value) +{ + mMinOpeningAngle = value; + LOG(info) << "V0 Photon Cut, min opening angle for too close V0: " << mMinOpeningAngle; +} + void V0PhotonCut::SetTPCNsigmaElRange(float min, float max) { mMinTPCNsigmaEl = min; diff --git a/PWGEM/PhotonMeson/Core/V0PhotonCut.h b/PWGEM/PhotonMeson/Core/V0PhotonCut.h index f73cbf6d16a..664981df21e 100644 --- a/PWGEM/PhotonMeson/Core/V0PhotonCut.h +++ b/PWGEM/PhotonMeson/Core/V0PhotonCut.h @@ -27,6 +27,9 @@ #include #include +#include // IWYU pragma: keep +#include +#include #include #include @@ -183,6 +186,7 @@ class V0PhotonCut : public TNamed kRZLine, kOnWwireIB, kOnWwireOB, + kIsTooClose, // leg cut kTrackPtRange, kTrackEtaRange, @@ -206,6 +210,12 @@ class V0PhotonCut : public TNamed kNCuts }; + enum class TooCloseCuts : uint8_t { + kNoCut = 0, + kDistance3D = 1, + kRadAndAngle = 2 + }; + /// \brief add histograms to registry /// \param fRegistry pointer to histogram registry void addQAHistograms(o2::framework::HistogramRegistry* fRegistry = nullptr) const @@ -268,26 +278,27 @@ class V0PhotonCut : public TNamed hPhotonQualityCuts->GetXaxis()->SetBinLabel(12, "RZ_{line}"); hPhotonQualityCuts->GetXaxis()->SetBinLabel(13, "Wire_{IB}"); hPhotonQualityCuts->GetXaxis()->SetBinLabel(14, "Wire_{OB}"); - hPhotonQualityCuts->GetXaxis()->SetBinLabel(15, "#it{p}_{T,leg}"); - hPhotonQualityCuts->GetXaxis()->SetBinLabel(16, "#it{#eta}_{leg}"); - hPhotonQualityCuts->GetXaxis()->SetBinLabel(17, "#it{N}_{cl,TPC}"); - hPhotonQualityCuts->GetXaxis()->SetBinLabel(18, "#it{N}_{cr,TPC}"); - hPhotonQualityCuts->GetXaxis()->SetBinLabel(19, "#it{N}_{cr,TPC}/#it{N}_{cl,TPC}"); - hPhotonQualityCuts->GetXaxis()->SetBinLabel(20, "FracSharedCl"); - hPhotonQualityCuts->GetXaxis()->SetBinLabel(21, "#chi^{2}_{TPC}/NDF"); - hPhotonQualityCuts->GetXaxis()->SetBinLabel(22, "#it{N#sigma}_{e,TPC}"); - hPhotonQualityCuts->GetXaxis()->SetBinLabel(23, "#it{N#sigma}_{#pi,TPC}"); - hPhotonQualityCuts->GetXaxis()->SetBinLabel(24, "DCA_{xy}"); - hPhotonQualityCuts->GetXaxis()->SetBinLabel(25, "DCA_{z}"); - hPhotonQualityCuts->GetXaxis()->SetBinLabel(26, "#it{N}_{cl,ITS}"); - hPhotonQualityCuts->GetXaxis()->SetBinLabel(27, "#chi^{2}_{ITS}/NDF"); - hPhotonQualityCuts->GetXaxis()->SetBinLabel(28, "size_{ITS}"); - hPhotonQualityCuts->GetXaxis()->SetBinLabel(29, "ITSTPC"); - hPhotonQualityCuts->GetXaxis()->SetBinLabel(30, "ITSOnly"); - hPhotonQualityCuts->GetXaxis()->SetBinLabel(31, "TPCOnly"); - hPhotonQualityCuts->GetXaxis()->SetBinLabel(32, "TPCTRD"); - hPhotonQualityCuts->GetXaxis()->SetBinLabel(33, "TPCTOF"); - hPhotonQualityCuts->GetXaxis()->SetBinLabel(34, "Out"); + hPhotonQualityCuts->GetXaxis()->SetBinLabel(15, "IsTooClose"); + hPhotonQualityCuts->GetXaxis()->SetBinLabel(16, "#it{p}_{T,leg}"); + hPhotonQualityCuts->GetXaxis()->SetBinLabel(17, "#it{#eta}_{leg}"); + hPhotonQualityCuts->GetXaxis()->SetBinLabel(18, "#it{N}_{cl,TPC}"); + hPhotonQualityCuts->GetXaxis()->SetBinLabel(19, "#it{N}_{cr,TPC}"); + hPhotonQualityCuts->GetXaxis()->SetBinLabel(20, "#it{N}_{cr,TPC}/#it{N}_{cl,TPC}"); + hPhotonQualityCuts->GetXaxis()->SetBinLabel(21, "FracSharedCl"); + hPhotonQualityCuts->GetXaxis()->SetBinLabel(22, "#chi^{2}_{TPC}/NDF"); + hPhotonQualityCuts->GetXaxis()->SetBinLabel(23, "#it{N#sigma}_{e,TPC}"); + hPhotonQualityCuts->GetXaxis()->SetBinLabel(24, "#it{N#sigma}_{#pi,TPC}"); + hPhotonQualityCuts->GetXaxis()->SetBinLabel(25, "DCA_{xy}"); + hPhotonQualityCuts->GetXaxis()->SetBinLabel(26, "DCA_{z}"); + hPhotonQualityCuts->GetXaxis()->SetBinLabel(27, "#it{N}_{cl,ITS}"); + hPhotonQualityCuts->GetXaxis()->SetBinLabel(28, "#chi^{2}_{ITS}/NDF"); + hPhotonQualityCuts->GetXaxis()->SetBinLabel(29, "size_{ITS}"); + hPhotonQualityCuts->GetXaxis()->SetBinLabel(30, "ITSTPC"); + hPhotonQualityCuts->GetXaxis()->SetBinLabel(31, "ITSOnly"); + hPhotonQualityCuts->GetXaxis()->SetBinLabel(32, "TPCOnly"); + hPhotonQualityCuts->GetXaxis()->SetBinLabel(33, "TPCTRD"); + hPhotonQualityCuts->GetXaxis()->SetBinLabel(34, "TPCTOF"); + hPhotonQualityCuts->GetXaxis()->SetBinLabel(35, "Out"); } } @@ -351,6 +362,119 @@ class V0PhotonCut : public TNamed fRegistry->fill(HIST("QA/V0Photon/after/Neg/hTPCHits"), ele.tpcNClsFound(), ele.tpcNClsCrossedRows()); } + /// \brief creates a mask for the V0s if they are too close to another V0 and have higher chi^2 + /// \param v0s V0 table + template + void createCloseV0CutMask(TV0 const& v0s) const + { + const bool useDistance3D = (mTooCloseType == TooCloseCuts::kDistance3D); + const float windowWidth = useDistance3D ? std::sqrt(mMinV0DistSquared) : mDeltaR; + const float cosMinAngle = std::cos(mMinOpeningAngle); + + int tableSize = v0s.size(); + std::vector rejectMask(tableSize, 0); + + if (mTooCloseType == TooCloseCuts::kNoCut) { + mRejectMask = rejectMask; + return; + } + + auto currentV0Iter = v0s.begin(); + auto otherV0Iter = v0s.begin(); + + int groupStart = 0; + while (groupStart < tableSize) { + + // --- find the end of this collision's group --- + int currentCollisionId = v0s.iteratorAt(groupStart).collisionId(); + int groupEnd = groupStart + 1; + while (groupEnd < tableSize && v0s.iteratorAt(groupEnd).collisionId() == currentCollisionId) { + groupEnd++; + } + + int groupSize = groupEnd - groupStart; + + std::vector> indexedRadii(groupSize); + for (int k = 0; k < groupSize; k++) { + currentV0Iter.setCursor(groupStart + k); + indexedRadii[k] = {groupStart + k, currentV0Iter.v0radius()}; + } + std::sort(indexedRadii.begin(), indexedRadii.end(), [](const auto& a, const auto& b) { + return a.second < b.second; + }); + // extract sorted indices and pre-sorted radii + std::vector sortedIndices(groupSize); + std::vector sortedRadii(groupSize); + for (int k = 0; k < groupSize; k++) { + sortedIndices[k] = indexedRadii[k].first; + sortedRadii[k] = indexedRadii[k].second; + } + + // --- sliding window within this group --- + int windowStart = 0; // reset per group + for (int i = 0; i < groupSize; i++) { + + float currentRadius = sortedRadii[i]; + while (windowStart < groupSize && sortedRadii[windowStart] < currentRadius - windowWidth) { + windowStart++; + } + + currentV0Iter.setCursor(sortedIndices[i]); + + float vx1 = currentV0Iter.vx(); + float vy1 = currentV0Iter.vy(); + float vz1 = currentV0Iter.vz(); + + float px1 = currentV0Iter.px(); + float py1 = currentV0Iter.py(); + float pz1 = currentV0Iter.pz(); + float chi2I = currentV0Iter.chiSquareNDF(); + + for (int j = windowStart; j < groupSize; j++) { + if (j == i) { + continue; + } + if (sortedRadii[j] > currentRadius + windowWidth) { + break; + } + + otherV0Iter.setCursor(sortedIndices[j]); + + bool tooClose = false; + if (useDistance3D) { + float dx = vx1 - otherV0Iter.vx(); + float dy = vy1 - otherV0Iter.vy(); + float dz = vz1 - otherV0Iter.vz(); + float distSquared = dx * dx + dy * dy + dz * dz; + tooClose = distSquared < mMinV0DistSquared; + } else { + float px2 = otherV0Iter.px(), py2 = otherV0Iter.py(), pz2 = otherV0Iter.pz(); + float dot = px1 * px2 + py1 * py2 + pz1 * pz2; + float mag1 = px1 * px1 + py1 * py1 + pz1 * pz1; + float mag2 = px2 * px2 + py2 * py2 + pz2 * pz2; + float denom = std::sqrt(mag1 * mag2); + if (denom > 0) { + float cosAngle = dot / denom; + cosAngle = std::clamp(cosAngle, -1.0f, 1.0f); + tooClose = cosAngle > cosMinAngle; + } + } + + if (tooClose) { + if (chi2I > otherV0Iter.chiSquareNDF()) { + rejectMask[sortedIndices[i]] = 1; + + } else { + rejectMask[sortedIndices[j]] = 1; + } + } + } + } + groupStart = groupEnd; + } + mRejectMask = rejectMask; + } + /// \brief check if given v0 photon survives all cuts /// \param flags EMBitFlags where results will be stored /// \param v0s v0 photon table to check @@ -360,6 +484,8 @@ class V0PhotonCut : public TNamed if (v0s.size() <= 0) { return; } + + createCloseV0CutMask(v0s); // auto legIter = legs.begin(); // auto legEnd = legs.end(); size_t iV0 = 0; @@ -489,6 +615,13 @@ class V0PhotonCut : public TNamed } } + if (!IsSelectedV0(v0, V0PhotonCuts::kIsTooClose)) { + if (doQA) { + fRegistry->fill(HIST("QA/V0Photon/hPhotonQualityCuts"), static_cast(V0PhotonCuts::kIsTooClose) + 1, v0Pt); + } + return false; + } + for (const auto& track : {pos, ele}) { if (!IsSelectedTrack(track, V0PhotonCuts::kTrackPtRange)) { if (doQA) { @@ -707,11 +840,24 @@ class V0PhotonCut : public TNamed case V0PhotonCuts::kAP: return std::pow(v0.alpha() / mMaxAlpha, 2) + std::pow(v0.qtarm() / mMaxQt, 2) < 1.0; - case V0PhotonCuts::kPsiPair: - return true; + // TODO: implement fully + case V0PhotonCuts::kPsiPair: { + if constexpr (requires { v0.psipair(); }) { + // return (std::fabs(v0.psipair() < 0.18f * std::exp( -0.55f * v0.chiSquareNDF()))); + return true; + } else { + return true; + } + } - case V0PhotonCuts::kPhiV: - return true; + // TODO: implement fully + case V0PhotonCuts::kPhiV: { + if constexpr (requires { v0.phiv(); }) { + return true; + } else { + return true; + } + } case V0PhotonCuts::kRxy: { if (v0.v0radius() < mMinRxy || mMaxRxy < v0.v0radius()) { @@ -775,6 +921,12 @@ class V0PhotonCut : public TNamed float dxy = std::sqrt(std::pow(x - x_exp, 2) + std::pow(y - y_exp, 2)); return !(dxy > margin_xy); } + case V0PhotonCuts::kIsTooClose: { + if (mRejectMask.size() == 0) { + return true; + } + return (mRejectMask[v0.globalIndex()] == 0); + } default: return false; } @@ -950,6 +1102,11 @@ class V0PhotonCut : public TNamed void SetOnWwireOB(bool flag = false); void RejectITSib(bool flag = false); + void setTooCloseType(V0PhotonCut::TooCloseCuts type); + void setMinV0DistSquared(float value); + void setDeltaR(float value); + void setMinOpeningAngle(float value); + void SetTrackPtRange(float minPt = 0.f, float maxPt = 1e10f); void SetTrackEtaRange(float minEta = -1e10f, float maxEta = 1e10f); void SetMinNClustersTPC(int minNClustersTPC); @@ -1017,6 +1174,11 @@ class V0PhotonCut : public TNamed bool mIsOnWwireIB{false}; bool mIsOnWwireOB{false}; bool mRejectITSib{false}; + TooCloseCuts mTooCloseType{V0PhotonCut::TooCloseCuts::kRadAndAngle}; // for TooCloseV0Cut: either squared distance between conversion points OR opening angle and deltaR + float mMinV0DistSquared{1.}; // for TooCloseV0Cut: cut value when using squared distance between conversion points + float mDeltaR{6.}; // for TooCloseV0Cut: V0PhotonCut::TooCloseCuts::kRadAndAngle when deltaR < this -> compare chi2 + float mMinOpeningAngle{0.02}; // for TooCloseV0Cut: V0PhotonCut::TooCloseCuts::kRadAndAngle when opening angle < this -> compare chi2 + mutable std::vector mRejectMask{}; // ML cuts bool mApplyMlCuts{false}; diff --git a/PWGEM/PhotonMeson/Tasks/calibTaskEmc.cxx b/PWGEM/PhotonMeson/Tasks/calibTaskEmc.cxx index 1d9ca1dbfeb..8d0b0ee7b3a 100644 --- a/PWGEM/PhotonMeson/Tasks/calibTaskEmc.cxx +++ b/PWGEM/PhotonMeson/Tasks/calibTaskEmc.cxx @@ -183,6 +183,11 @@ struct CalibTaskEmc { o2::framework::Configurable rejectV0onITSib{"rejectV0onITSib", true, "flag to reject V0s on ITSib"}; o2::framework::Configurable applyPrefilter{"applyPrefilter", false, "flag to apply prefilter to V0"}; + o2::framework::Configurable mTooCloseType{"mTooCloseType", 2, "type of cut for too close (0 = no, 1 = squared distance, 2 = opening angle + dR) "}; + o2::framework::Configurable mMinV0DistSquared{"mMinV0DistSquared", 4.f, "min squared distance for mTooCloseType == 1"}; + o2::framework::Configurable mDeltaR{"mDeltaR", 6.f, "deltaR for mTooCloseType == 2"}; + o2::framework::Configurable mMinOpeningAngle{"mMinOpeningAngle", 0.02, "min opening angle for mTooCloseType == 2"}; + o2::framework::Configurable minNClusterTPC{"minNClusterTPC", 0, "min NCluster TPC"}; o2::framework::Configurable minNCrossedRowsTPC{"minNCrossedRowsTPC", 40, "min ncrossed rows in TPC"}; o2::framework::Configurable minNCrossedRowsOverFindableClustersTPC{"minNCrossedRowsOverFindableClustersTPC", 0.8f, "min fraction of crossed rows over findable clusters in TPC"}; @@ -324,6 +329,10 @@ struct CalibTaskEmc { fV0PhotonCut.SetRxyRange(pcmcuts.minRV0, pcmcuts.maxRV0); fV0PhotonCut.SetAPRange(pcmcuts.maxAlphaAP, pcmcuts.maxQtAP); fV0PhotonCut.RejectITSib(pcmcuts.rejectV0onITSib); + fV0PhotonCut.setTooCloseType(static_cast(pcmcuts.mTooCloseType.value)); + fV0PhotonCut.setMinV0DistSquared(pcmcuts.mMinV0DistSquared); + fV0PhotonCut.setDeltaR(pcmcuts.mDeltaR); + fV0PhotonCut.setMinOpeningAngle(pcmcuts.mMinOpeningAngle); // for track fV0PhotonCut.SetMinNClustersTPC(pcmcuts.minNClusterTPC); @@ -364,20 +373,25 @@ struct CalibTaskEmc { const AxisSpec thAxisEnergy{1000, 0., 100., "#it{E}_{clus} (GeV)"}; const AxisSpec thAxisEta{320, -0.8, 0.8, "#eta"}; const AxisSpec thAxisPhi{500, 0, 2 * 3.14159, "phi"}; + const AxisSpec thAxisOpeningAngle{180, 0, o2::constants::math::PI, "opening angle (rad)"}; const AxisSpec thnAxisMixingVtx{mixingConfig.cfgVtxBins, "#it{z} (cm)"}; const AxisSpec thnAxisMixingCent{mixingConfig.cfgCentBins, "Centrality (%)"}; if (doprocessEMCal || doprocessEMCalPCMC) { - registry.add("hSparsePi0", "m_{inv} vs p_T vs cent for same event", HistType::kTH3D, {thnAxisInvMass, thAxisEnergyCalib, thnAxisCent}); + registry.add("hSparsePi0", "m_{inv} vs E vs cent for same event", HistType::kTH3D, {thnAxisInvMass, thAxisEnergyCalib, thnAxisCent}); + registry.add("hOpeningAngleSE", "opening vs E vs cent for same event", HistType::kTH3D, {thAxisOpeningAngle, thAxisEnergyCalib, thnAxisCent}); } else if (doprocessPCM) { registry.add("hSparsePi0", "m_{inv} vs p_T vs cent for same event", HistType::kTH3D, {thnAxisInvMass, thnAxisPtCalib, thnAxisCent}); + registry.add("hOpeningAngleSE", "opening vs p_T vs cent for same event", HistType::kTH3D, {thAxisOpeningAngle, thnAxisPtCalib, thnAxisCent}); } if (doprocessEMCalMixed || doprocessEMCalPCMMixed) { registry.add("hSparseBkgMix", "m_{inv} vs p_T vs cent for mixed event", HistType::kTH3D, {thnAxisInvMass, thAxisEnergyCalib, thnAxisCent}); + registry.add("hOpeningAngleME", "opening vs E vs cent for same event", HistType::kTH3D, {thAxisOpeningAngle, thAxisEnergyCalib, thnAxisCent}); } else if (doprocessPCMMixed) { registry.add("hSparseBkgMix", "m_{inv} vs p_T vs cent for mixed event", HistType::kTH3D, {thnAxisInvMass, thnAxisPtCalib, thnAxisCent}); + registry.add("hOpeningAngleME", "opening vs p_T vs cent for same event", HistType::kTH3D, {thAxisOpeningAngle, thnAxisPtCalib, thnAxisCent}); } if (doprocessEMCalMixed || doprocessEMCalPCMMixed || doprocessPCMMixed) { @@ -428,7 +442,6 @@ struct CalibTaskEmc { /// \param mass is the invariant mass of the candidate /// \param pt is the transverse momentum of the candidate /// \param cent is the centrality of the collision - /// \param sp is the scalar product template void fillThn(const float mass, const float pt, const float cent) { @@ -436,6 +449,17 @@ struct CalibTaskEmc { registry.fill(HIST(HistTypes[histType]), mass, pt, cent); } + /// Fill THnSparse + /// \param openingAngle opening angle between the two photons + /// \param pt is the transverse momentum of the candidate + /// \param cent is the centrality of the collision + template + void fillOpeningAngleHisto(const float openingAngle, const float pt, const float cent) + { + static constexpr std::string_view HistTypes[3] = {"hOpeningAngleSE", "hOpeningAngleRot", "hOpeningAngleME"}; + registry.fill(HIST(HistTypes[histType]), openingAngle, pt, cent); + } + /// Get the centrality /// \param collision is the collision with the centrality information template @@ -636,6 +660,8 @@ struct CalibTaskEmc { continue; } registry.fill(HIST("hMesonCuts"), 6); + float cent = getCentrality(collision); + fillOpeningAngleHisto<0>(openingAngle, vMeson.Pt(), cent); runFlowAnalysis<0>(collision, vMeson, g1.e()); } } @@ -747,6 +773,8 @@ struct CalibTaskEmc { continue; } registry.fill(HIST("hMesonCutsMixed"), 6); + float cent = getCentrality(c1); + fillOpeningAngleHisto<2>(openingAngle, vMeson.Pt(), cent); runFlowAnalysis<2>(c1, vMeson, g1.e()); } } @@ -823,6 +851,8 @@ struct CalibTaskEmc { registry.fill(HIST("hMesonCuts"), 5); continue; } + float cent = getCentrality(collision); + fillOpeningAngleHisto<0>(openingAngle, vMeson.Pt(), cent); runFlowAnalysis<0>(collision, vMeson, g1.corrE()); } } @@ -913,6 +943,8 @@ struct CalibTaskEmc { continue; } registry.fill(HIST("hMesonCutsMixed"), 6); + float cent = getCentrality(c1); + fillOpeningAngleHisto<2>(openingAngle, vMeson.Pt(), cent); runFlowAnalysis<2>(c1, vMeson, g1.corrE()); } } @@ -973,6 +1005,8 @@ struct CalibTaskEmc { registry.fill(HIST("mesonQA/hAlphaPt"), asymmetry, photon1Pt); } registry.fill(HIST("hMesonCuts"), 6); + float cent = getCentrality(collision); + fillOpeningAngleHisto<0>(openingAngle, vMeson.Pt(), cent); runFlowAnalysis<0>(collision, vMeson, photon1Pt); } } // end of loop over collisions @@ -1047,6 +1081,8 @@ struct CalibTaskEmc { registry.fill(HIST("mesonQA/hAlphaPtMixed"), asymmetry, photon1Pt); } registry.fill(HIST("hMesonCutsMixed"), 6); + float cent = getCentrality(c1); + fillOpeningAngleHisto<2>(openingAngle, vMeson.Pt(), cent); runFlowAnalysis<2>(c1, vMeson, photon1Pt); } }