diff --git a/PWGDQ/Tasks/global-muon-matcher.cxx b/PWGDQ/Tasks/global-muon-matcher.cxx index c16f2917b6e..390466dddc3 100644 --- a/PWGDQ/Tasks/global-muon-matcher.cxx +++ b/PWGDQ/Tasks/global-muon-matcher.cxx @@ -9,13 +9,14 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. // -/// \file global-muon-matcher.cxx +/// \file global-muon-matcher.cxx // o2-linter: disable=name/file-cpp,name/workflow-file (legacy workflow executable name) /// \brief Task for analysis MFT-MCH muon matching /// \author Andrea Ferrero /// #include "PWGDQ/Core/MuonMatchingMlResponse.h" #include "PWGDQ/Core/VarManager.h" +#include "Common/Core/RecoDecay.h" #include "Common/Core/fwdtrackUtilities.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/FwdTrackReAlignTables.h" @@ -59,6 +60,7 @@ #include #include #include +#include #include #include #include @@ -74,37 +76,19 @@ using namespace o2; using namespace o2::framework; using namespace o2::aod; -namespace o2::aod::globalmuonmatching -{ -DECLARE_SOA_COLUMN(MchTrackId, mchTrackId, int64_t); -DECLARE_SOA_COLUMN(MftTrackId, mftTrackId, int64_t); -DECLARE_SOA_COLUMN(MatchChi2, matchChi2, float); -DECLARE_SOA_COLUMN(MatchScore, matchScore, float); -DECLARE_SOA_COLUMN(MatchRanking, matchRanking, int32_t); -DECLARE_SOA_COLUMN(IsTagged, isTagged, bool); -} // namespace o2::aod::globalmuonmatching - namespace o2::aod { -DECLARE_SOA_TABLE(GlobalMuonMatchCandidates, "AOD", "GMCAND", - o2::soa::Index<>, - globalmuonmatching::MchTrackId, - globalmuonmatching::MftTrackId, - globalmuonmatching::MatchChi2, globalmuonmatching::MatchScore, globalmuonmatching::MatchRanking, - globalmuonmatching::IsTagged); - namespace globalmuonmatching { -DECLARE_SOA_ARRAY_INDEX_COLUMN(GlobalMuonMatchCandidate, globalMuonMatchCandidate); //! Array of GlobalMuonMatchCandidates indices -DECLARE_SOA_INDEX_COLUMN_FULL(FwdTrackRealign, fwdTrackRealign, int, FwdTracksReAlign, ""); //! Index of ambiguous FwdTracksReAlign entry -DECLARE_SOA_SLICE_INDEX_COLUMN_FULL(Bc, bc, int32_t, BCs, ""); //! BC index slice compatible with the track time window +DECLARE_SOA_COLUMN(IsTagged, isTagged, bool); //! Whether the MCH(-MID) track passes tagging cuts +DECLARE_SOA_COLUMN(MatchRanking, matchRanking, int32_t); //! Match candidate ranking (-1 for base MCH entries) +DECLARE_SOA_COLUMN(MixedGroupIndex, mixedGroupIndex, int32_t); //! Mixed-event group index (-1 for same-event candidates) } // namespace globalmuonmatching -DECLARE_SOA_TABLE(FwdTrkMatchCands, "AOD", "FWDTRKMATCHCAND", //! Vectors of match-candidate indices stored per fwdtrack - globalmuonmatching::GlobalMuonMatchCandidateIds, o2::soa::Marker<3>); - -DECLARE_SOA_TABLE(AmbiguousFwdTracksReAlign, "AOD", "AMBIGFWDREALIGN", //! FwdTracksReAlign entries without a unique collision association - o2::soa::Index<>, globalmuonmatching::FwdTrackRealignId, globalmuonmatching::BcIdSlice); +DECLARE_SOA_TABLE(GmmCandFwdTrkExtras, "AOD", "GMMCANDEXTRA", //! Extra info joinable to FwdTracksReAlign + globalmuonmatching::IsTagged, + globalmuonmatching::MatchRanking, + globalmuonmatching::MixedGroupIndex); } // namespace o2::aod using MyEvents = soa::Join; @@ -119,6 +103,12 @@ using SMatrix5 = o2::track::SMatrix5; constexpr std::array NDetElemCh = {4, 4, 4, 4, 18, 18, 26, 26, 26, 26}; constexpr std::array SNDetElemCh = {0, 4, 8, 12, 16, 34, 52, 78, 104, 130, 156}; +// compute minimum difference between azimuthal angles +static float getDeltaPhi(float phi1, float phi2) +{ + return RecoDecay::constrainAngle(phi1 - phi2, -o2::constants::math::PI); +} + struct GlobalMuonMatching { static constexpr int GlobalTrackTypeMax = 2; @@ -137,6 +127,16 @@ struct GlobalMuonMatching { double matchScore{-1}; double matchChi2{-1}; int matchRanking{-1}; + int32_t mixedGroupIndex{-1}; + }; + + struct MchTrackInfo { + int nMatchAttempts{-1}; + bool isTagged{false}; + // vector of MFT-MCH matching candidates + std::vector matchingCandidates; + // vector of vectors of MFT-MCH matching candidates from mixed events + std::vector> mixedMatchingCandidates; }; //// Variables for selecting tagged muons @@ -205,6 +205,15 @@ struct GlobalMuonMatching { Configurable cfgMatchAllTracks{"cfgMatchAllTracks", false, "If true the matching is performed considering all the MFT tracks for which the covariances are available; if false the matching is performed considering only the global forward tracks stored at production"}; } configMatching; + struct : ConfigurableGroup { + Configurable cfgMixingDepth{"cfgMixingDepth", -1, "Maximum number of mixed candidate groups per MCH track (-1: no limit)"}; + Configurable cfgMinDeltaBc{"cfgMinDeltaBc", 3564, "Minimum DeltaBc between mixed collisions"}; + Configurable cfgMaxDeltaPhi{"cfgMaxDeltaPhi", static_cast(o2::constants::math::PI / 10), "Maximum DelptaPhi between mixed MCH tracks (rad)"}; + Configurable cfgMaxDeltaR{"cfgMaxDeltaR", 10.f, "Maximum DeltaR between mixed MCH tracks"}; + Configurable cfgMaxDeltaAttempts{"cfgMaxDeltaAttempts", 0.1f, "Maximum relative difference in match attempts"}; + Configurable cfgMaxDeltaZ{"cfgMaxDeltaZ", 1.f, "Maximum deltaZ between mixed collisions"}; + } configEventMixing; + double mBzAtMftCenter{0}; using MatchingFunc = std::function(const o2::track::TrackParCovFwd& mchtrack, const o2::track::TrackParCovFwd& mfttrack)>; @@ -246,7 +255,6 @@ struct GlobalMuonMatching { // to be sorted in descending order with respect to the matching score // the map key is the MCH(-MID) track global index using MatchingCandidates = std::map>; - std::map> mMatchingCandidates; class TrackParExt : public o2::track::TrackParCovFwd { @@ -276,7 +284,10 @@ struct GlobalMuonMatching { void setRemovable() { removable = true; } [[nodiscard]] bool isRemovable() const { return removable; } - [[nodiscard]] o2::track::TrackParCovFwd asTrackParCovFwd() const { return *this; } + [[nodiscard]] o2::track::TrackParCovFwd asTrackParCovFwd() const + { + return {static_cast(*this)}; + } private: int nClusters{-1}; @@ -288,20 +299,18 @@ struct GlobalMuonMatching { std::unordered_map mftTrackCovs; - Produces globalMuonMatchCandidates; - Produces fwdTrkMatchCands; Produces gmCandidateFwdTracks; Produces gmCandidateFwdTracksCov; - Produces gmAmbiguousFwdTracksReAlign; + Produces gmCandidateFwdTrackExtras; + Produces gmAmbiguousFwdTracksReAlign; int32_t mGmmCandFwdTrackRowIndex{0}; std::unordered_map> mAmbBcSliceByFwdTrackId; bool mHasLastMchAmbiguousBcSlice{false}; std::array mLastMchAmbiguousBcSlice{}; - int32_t mMatchCandidateCounter{0}; - std::unordered_map> mMchTrackToCandidateIndices; - std::unordered_map> mMchTrackMatchingCandidates; + std::unordered_map mMchTrackInfos; + std::unordered_map> mStoredMatchingCandidates; std::unordered_map mFwdTrackToGmmCandTrkIndex; mch::TrackFitter trackFitter; // Track fitter from MCH tracking library @@ -873,12 +882,19 @@ struct GlobalMuonMatching { rhoXY, rhoPhiY, rhoPhiX, rhoTglX, rhoTglY, rhoTglPhi, rho1PtX, rho1PtY, rho1PtPhi, rho1PtTgl); } + bool isMchTrackTagged(int64_t mchTrackIndex) const + { + const auto it = mMchTrackInfos.find(mchTrackIndex); + return it != mMchTrackInfos.end() && it->second.isTagged; + } + template void fillBaseGmmCandFwdTrack(TMCH const& track, TrackParExt const& trackPar, int32_t gmmMchTrackId, float chi2MatchMCHMFT, - float matchScoreMCHMFT) + float matchScoreMCHMFT, + bool isTagged) { const auto collisionId = track.collisionId(); bool hasBcSlice = false; @@ -917,6 +933,7 @@ struct GlobalMuonMatching { track.trackTimeRes()); storeFwdTrackCovariance(trackPar.getCovariances()); + gmCandidateFwdTrackExtras(isTagged, -1, -1); if (hasBcSlice) { gmAmbiguousFwdTracksReAlign(mGmmCandFwdTrackRowIndex, bcSlice.data()); } @@ -939,7 +956,7 @@ struct GlobalMuonMatching { using o2::aod::fwdtrack::ForwardTrackTypeEnum; using o2::aod::fwdtrackutils::propagationPoint; - constexpr uint8_t CandidateTrackType = static_cast(ForwardTrackTypeEnum::GlobalForwardTrack); + constexpr auto CandidateTrackType = static_cast(ForwardTrackTypeEnum::GlobalForwardTrack); auto propmuonAtMft = fwdToMch(mchPar); o2::mch::TrackExtrap::extrapToVertex(propmuonAtMft, @@ -953,7 +970,7 @@ struct GlobalMuonMatching { const auto nClusters = static_cast(std::min(127, mchPar.getNClusters() + mftPar.getNClusters())); - const float chi2 = static_cast(mchTrack.chi2()); + const auto chi2 = static_cast(mchTrack.chi2()); const int32_t collisionId = mchTrack.collisionId(); bool hasBcSlice = false; std::array bcSlice{}; @@ -998,6 +1015,7 @@ struct GlobalMuonMatching { mchTrack.trackTimeRes()); storeFwdTrackCovariance(globalMuonRefit.getCovariances()); + gmCandidateFwdTrackExtras(isMchTrackTagged(mchTrack.globalIndex()), candidate.matchRanking, candidate.mixedGroupIndex); if (hasBcSlice) { gmAmbiguousFwdTracksReAlign(mGmmCandFwdTrackRowIndex, bcSlice.data()); } @@ -1120,6 +1138,41 @@ struct GlobalMuonMatching { return std::fabs(deltaTrackTime) <= trackTimeResTot; } + template + int getMftMchMatchAttempts(EVT const& collisions, + BC const& bcs, + TMUON const& mchTrack, + TMFTS const& mftTracks) + { + if (!mchTrack.has_collision()) { + return 0; + } + const auto& collMch = collisions.rawIteratorAt(mchTrack.collisionId()); + const auto& bcMch = bcs.rawIteratorAt(collMch.bcId()); + + int attempts{0}; + for (const auto& mftTrack : mftTracks) { + if (!mftTrack.has_collision()) { + continue; + } + + const auto& collMft = collisions.rawIteratorAt(mftTrack.collisionId()); + const auto& bcMft = bcs.rawIteratorAt(collMft.bcId()); + + int64_t deltaBc = static_cast(bcMft.globalBC()) - static_cast(bcMch.globalBC()); + double deltaBcNS = o2::constants::lhc::LHCBunchSpacingNS * deltaBc; + double deltaTrackTime = mftTrack.trackTime() - mchTrack.trackTime() + deltaBcNS; + double trackTimeResTot = mftTrack.trackTimeRes() + mchTrack.trackTimeRes(); + + if (std::fabs(deltaTrackTime) > trackTimeResTot) { + continue; + } + attempts += 1; + } + + return attempts; + } + template void prepareMatchingCandidates(EVT const& collisions, BC const& bcs, @@ -1129,7 +1182,6 @@ struct GlobalMuonMatching { { mMftTrackPars.clear(); mMchTrackPars.clear(); - mMatchingCandidates.clear(); LOGF(info, "Filling matching candidate tables"); @@ -1141,6 +1193,11 @@ struct GlobalMuonMatching { // initialize the MCH track parameters, which will be updated by the realignment if enabled mMchTrackPars.try_emplace(mchTrackIndex, TrackParExt(fwdtrackutils::getTrackParCovFwd(muonTrack, muonTrack), muonTrack.nClusters())); + + mMchTrackInfos.try_emplace(mchTrackIndex, MchTrackInfo{ + .nMatchAttempts = 0, + .matchingCandidates = std::vector(), + .mixedMatchingCandidates = std::vector>()}); } for (const auto& mftTrack : mftTracks) { @@ -1171,12 +1228,19 @@ struct GlobalMuonMatching { continue; } - mMatchingCandidates[mchTrackIndex].emplace_back(MatchingCandidate{ + auto& mchTrackInfo = mMchTrackInfos[mchTrackIndex]; + mchTrackInfo.matchingCandidates.emplace_back(MatchingCandidate{ .muonTrackId = muonTrack.globalIndex(), .mftTrackId = mftTrackIndex, .matchScore = muonTrack.matchScoreMCHMFT(), .matchChi2 = muonTrack.chi2MatchMCHMFT()}); } + + // set the number of match attempts for this track + for (auto& mchTrackInfo : mMchTrackInfos) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) + const auto& mchTrack = muonTracks.rawIteratorAt(mchTrackInfo.first); + mchTrackInfo.second.nMatchAttempts = getMftMchMatchAttempts(collisions, bcs, mchTrack, mftTracks); + } } else { // build matching candidates from all time-compatible MFT-MCH pairs for (const auto& muonTrack : muonTracks) { @@ -1192,19 +1256,111 @@ struct GlobalMuonMatching { continue; } - mMatchingCandidates[mchTrackIndex].emplace_back(MatchingCandidate{ + auto& mchTrackInfo = mMchTrackInfos[mchTrackIndex]; + mchTrackInfo.matchingCandidates.emplace_back(MatchingCandidate{ .mftTrackId = mftTrack.globalIndex()}); } } + + // set the number of match attempts for this track + for (auto& mchTrackInfo : mMchTrackInfos) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) + mchTrackInfo.second.nMatchAttempts = mchTrackInfo.second.matchingCandidates.size(); + } } - // sort the vectors of matching candidates in ascending order based on the matching chi2 value - auto compareMatchingChi2 = [](const MatchingCandidate& track1, const MatchingCandidate& track2) -> bool { - return (track1.matchChi2 < track2.matchChi2); - }; + std::vector taggedMuons; + getTaggedMuons(collisions, muonTracks, taggedMuons); + for (const auto& mchTrackIndex : taggedMuons) { + auto it = mMchTrackInfos.find(mchTrackIndex); + if (it != mMchTrackInfos.end()) { + it->second.isTagged = true; + } + } + } + + template + void prepareEventMixingMatchingCandidates(EVT const& collisions, + BC const& bcs, + TMUON const& muonTracks) + { + LOGF(info, "Filling mixed matching candidate tables"); + + const int mixingDepth = configEventMixing.cfgMixingDepth.value; + const int64_t minDeltaBc = configEventMixing.cfgMinDeltaBc.value; + const float maxDeltaPhi = configEventMixing.cfgMaxDeltaPhi.value; + const float maxDeltaR = configEventMixing.cfgMaxDeltaR.value; + const float maxDeltaAttemptsRel = configEventMixing.cfgMaxDeltaAttempts.value; + const float maxDeltaZ = configEventMixing.cfgMaxDeltaZ.value; + + for (auto& [mchIndex1, mchTrackInfo1] : mMchTrackInfos) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) + const auto& mchTrack1 = muonTracks.rawIteratorAt(mchIndex1); + + if (!mchTrack1.has_collision()) { + continue; + } + + const auto& collision1 = collisions.rawIteratorAt(mchTrack1.collisionId()); + const auto& bc1 = bcs.rawIteratorAt(collision1.bcId()); + + auto phi1 = std::atan2(mchTrack1.y(), mchTrack1.x()); + auto r1 = std::hypot(mchTrack1.x(), mchTrack1.y()); + + auto nAttempts1 = mchTrackInfo1.nMatchAttempts; + + auto vz1 = collision1.posZ(); + + for (const auto& [mchIndex2, mchTrackInfo2] : mMchTrackInfos) { + if (mixingDepth >= 0 && static_cast(mchTrackInfo1.mixedMatchingCandidates.size()) >= mixingDepth) { + break; + } + + const auto& mchTrack2 = muonTracks.rawIteratorAt(mchIndex2); + + if (!mchTrack2.has_collision()) { + continue; + } + + const auto& collision2 = collisions.rawIteratorAt(mchTrack2.collisionId()); + const auto& bc2 = bcs.rawIteratorAt(collision2.bcId()); + + auto deltaBc = std::abs(static_cast(bc2.globalBC()) - static_cast(bc1.globalBC())); + if (deltaBc < minDeltaBc) { + continue; + } + + auto phi2 = std::atan2(mchTrack2.y(), mchTrack2.x()); + auto deltaPhi = std::fabs(getDeltaPhi(phi1, phi2)); + if (deltaPhi > maxDeltaPhi) { + continue; + } + + auto r2 = std::hypot(mchTrack2.x(), mchTrack2.y()); + auto deltaR = r2 - r1; + if (deltaR > maxDeltaR) { + continue; + } + + auto nAttempts2 = mchTrackInfo2.nMatchAttempts; + + float deltaAttempts = nAttempts2 - nAttempts1; + float deltaAttemptsRel = (nAttempts1 > 0) ? deltaAttempts / nAttempts1 : 0; + if (deltaAttemptsRel > maxDeltaAttemptsRel) { + continue; + } - for (auto& [mchIndex, candidatesVector] : mMatchingCandidates) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) - std::sort(candidatesVector.begin(), candidatesVector.end(), compareMatchingChi2); + auto vz2 = collision2.posZ(); + auto deltaZ = std::fabs(vz2 - vz1); + if (deltaZ > maxDeltaZ) { + continue; + } + + // add the candidates of MCH track #2 to the list of mixed candidates of track #1 + mchTrackInfo1.mixedMatchingCandidates.push_back(mchTrackInfo2.matchingCandidates); + // update the muon track index of the mixed candidates to the index of track #1 + for (auto& candidate : mchTrackInfo1.mixedMatchingCandidates.back()) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) + candidate.muonTrackId = mchIndex1; + } + } } } @@ -1326,8 +1482,8 @@ struct GlobalMuonMatching { void runChi2Matching(const std::string& funcName, float matchingPlaneZ, - const MatchingCandidates& matchingCandidates, - MatchingCandidates& newMatchingCandidates) + MatchingCandidates& newMatchingCandidates, + bool useMixedMatchingCandidates) { newMatchingCandidates.clear(); @@ -1343,50 +1499,74 @@ struct GlobalMuonMatching { } auto matchingFunc = mMatchingFunctionMap.at(funcNameEffective); - for (const auto& [mchIndex, candidatesVector] : matchingCandidates) { - + for (const auto& [mchIndex, mchTrackInfo] : mMchTrackInfos) { // get the tracks parameters, which have been updated by the realignment if enabled const auto mchTrackParIt = mMchTrackPars.find(mchIndex); if (mchTrackParIt == mMchTrackPars.end()) { continue; } - for (const auto& candidate : candidatesVector) { - auto mftTrackParIt = mMftTrackPars.find(candidate.mftTrackId); - if (mftTrackParIt == mMftTrackPars.end()) { - continue; - } + auto processGroup = [&, mchIndex](const std::vector& candidatesGroup, int32_t mixedGroupIndex) { + std::vector groupResults; + groupResults.reserve(candidatesGroup.size()); + + for (const auto& candidate : candidatesGroup) { + auto mftTrackParIt = mMftTrackPars.find(candidate.mftTrackId); + if (mftTrackParIt == mMftTrackPars.end()) { + continue; + } + + auto mftTrackProp = mftTrackParIt->second.asTrackParCovFwd(); + auto mchTrackProp = mchTrackParIt->second.asTrackParCovFwd(); + + if (matchingPlaneZEffective < 0.) { + mftTrackProp = propagateToZMft(mftTrackProp, matchingPlaneZ); + mchTrackProp = propagateToZMch(mchTrackProp, matchingPlaneZ); + } - auto mftTrackProp = mftTrackParIt->second.asTrackParCovFwd(); - auto mchTrackProp = mchTrackParIt->second.asTrackParCovFwd(); + auto matchResult = matchingFunc(mchTrackProp, mftTrackProp); + float matchChi2 = std::get<0>(matchResult); - if (matchingPlaneZEffective < 0.) { - mftTrackProp = propagateToZMft(mftTrackProp, matchingPlaneZ); - mchTrackProp = propagateToZMch(mchTrackProp, matchingPlaneZ); + groupResults.emplace_back(MatchingCandidate{ + .muonTrackId = candidate.muonTrackId, + .mftTrackId = candidate.mftTrackId, + .matchScore = -1, + .matchChi2 = matchChi2}); } - auto matchResult = matchingFunc(mchTrackProp, mftTrackProp); - float matchChi2 = std::get<0>(matchResult); + auto compareMatchingChi2 = [](const MatchingCandidate& track1, const MatchingCandidate& track2) -> bool { + return (track1.matchChi2 < track2.matchChi2); + }; + std::sort(groupResults.begin(), groupResults.end(), compareMatchingChi2); - newMatchingCandidates[mchIndex].emplace_back(MatchingCandidate{ - .muonTrackId = candidate.muonTrackId, - .mftTrackId = candidate.mftTrackId, - .matchScore = -1, - .matchChi2 = matchChi2}); - } - } + int ranking = 1; + for (auto& matchedCandidate : groupResults) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) + matchedCandidate.matchRanking = ranking; + ranking += 1; + } - auto compareMatchingChi2 = [](const MatchingCandidate& track1, const MatchingCandidate& track2) -> bool { - return (track1.matchChi2 < track2.matchChi2); - }; + const int maxCandidates = configMatching.cfgMaxCandidatesPerMchTrack.value; - for (auto& [mchIndex, globalTracksVector] : newMatchingCandidates) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) - std::sort(globalTracksVector.begin(), globalTracksVector.end(), compareMatchingChi2); + auto& storedCandidates = newMatchingCandidates[mchIndex]; + size_t nStoredThisGroup = 0; + for (auto& result : groupResults) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) + if (maxCandidates >= 0 && nStoredThisGroup >= static_cast(maxCandidates)) { + break; + } + result.mixedGroupIndex = useMixedMatchingCandidates ? mixedGroupIndex : -1; + storedCandidates.push_back(result); + ++nStoredThisGroup; + } + }; - int ranking = 1; - for (auto& candidate : globalTracksVector) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) - candidate.matchRanking = ranking; - ranking += 1; + if (useMixedMatchingCandidates) { + int32_t groupIdx = 0; + for (const auto& candidatesGroup : mchTrackInfo.mixedMatchingCandidates) { + processGroup(candidatesGroup, groupIdx); + groupIdx += 1; + } + } else { + processGroup(mchTrackInfo.matchingCandidates, -1); } } } @@ -1397,11 +1577,11 @@ struct GlobalMuonMatching { TMFT const& mftTracks, o2::analysis::MlResponseMFTMuonMatch& mlResponse, float matchingPlaneZ, - const MatchingCandidates& matchingCandidates, - MatchingCandidates& newMatchingCandidates) + MatchingCandidates& newMatchingCandidates, + bool useMixedMatchingCandidates) { newMatchingCandidates.clear(); - for (const auto& [mchIndex, candidatesVector] : matchingCandidates) { + for (const auto& [mchIndex, mchTrackInfo] : mMchTrackInfos) { auto const& mchTrack = muonTracks.rawIteratorAt(mchIndex); if (!mchTrack.has_collision()) { continue; @@ -1415,46 +1595,71 @@ struct GlobalMuonMatching { continue; } - for (const auto& candidate : candidatesVector) { - auto const& muonTrack = (candidate.muonTrackId >= 0) ? muonTracks.rawIteratorAt(candidate.muonTrackId) : mchTrack; - auto const& mftTrack = mftTracks.rawIteratorAt(candidate.mftTrackId); - auto mftTrackParIt = mMftTrackPars.find(candidate.mftTrackId); - if (mftTrackParIt == mMftTrackPars.end()) { - continue; - } + auto processGroup = [&, mchIndex](const std::vector& candidatesGroup, int32_t mixedGroupIndex) { + std::vector groupResults; + groupResults.reserve(candidatesGroup.size()); - auto mftTrackProp = mftTrackParIt->second.asTrackParCovFwd(); - auto mchTrackProp = mchTrackParIt->second.asTrackParCovFwd(); + for (const auto& candidate : candidatesGroup) { + auto const& muonTrack = (candidate.muonTrackId >= 0) ? muonTracks.rawIteratorAt(candidate.muonTrackId) : mchTrack; + auto const& mftTrack = mftTracks.rawIteratorAt(candidate.mftTrackId); + auto mftTrackParIt = mMftTrackPars.find(candidate.mftTrackId); + if (mftTrackParIt == mMftTrackPars.end()) { + continue; + } + + auto mftTrackProp = mftTrackParIt->second.asTrackParCovFwd(); + auto mchTrackProp = mchTrackParIt->second.asTrackParCovFwd(); + + if (matchingPlaneZ < 0.) { + mftTrackProp = propagateToZMft(mftTrackProp, matchingPlaneZ); + mchTrackProp = propagateToZMch(mchTrackProp, matchingPlaneZ); + } - if (matchingPlaneZ < 0.) { - mftTrackProp = propagateToZMft(mftTrackProp, matchingPlaneZ); - mchTrackProp = propagateToZMch(mchTrackProp, matchingPlaneZ); + std::vector output; + std::vector inputML = mlResponse.getInputFeatures(muonTrack, mftTrack, mchTrack, mftTrackProp, mchTrackProp, collision); + mlResponse.isSelectedMl(inputML, 0, output); + float matchScore = output[0]; + + groupResults.emplace_back(MatchingCandidate{ + .muonTrackId = candidate.muonTrackId, + .mftTrackId = candidate.mftTrackId, + .matchScore = matchScore, + .matchChi2 = -1}); } - std::vector output; - std::vector inputML = mlResponse.getInputFeatures(muonTrack, mftTrack, mchTrack, mftTrackProp, mchTrackProp, collision); - mlResponse.isSelectedMl(inputML, 0, output); - float matchScore = output[0]; + auto compareMatchingScore = [](const MatchingCandidate& track1, const MatchingCandidate& track2) -> bool { + return (track1.matchScore > track2.matchScore); + }; + std::sort(groupResults.begin(), groupResults.end(), compareMatchingScore); - newMatchingCandidates[mchIndex].emplace_back(MatchingCandidate{ - .muonTrackId = candidate.muonTrackId, - .mftTrackId = candidate.mftTrackId, - .matchScore = matchScore, - .matchChi2 = -1}); - } - } + int ranking = 1; + for (auto& matchedCandidate : groupResults) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) + matchedCandidate.matchRanking = ranking; + ranking += 1; + } - auto compareMatchingScore = [](const MatchingCandidate& track1, const MatchingCandidate& track2) -> bool { - return (track1.matchScore > track2.matchScore); - }; + const int maxCandidates = configMatching.cfgMaxCandidatesPerMchTrack.value; - for (auto& [mchIndex, globalTracksVector] : newMatchingCandidates) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) - std::sort(globalTracksVector.begin(), globalTracksVector.end(), compareMatchingScore); + auto& storedCandidates = newMatchingCandidates[mchIndex]; + size_t nStoredThisGroup = 0; + for (auto& result : groupResults) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) + if (maxCandidates >= 0 && nStoredThisGroup >= static_cast(maxCandidates)) { + break; + } + result.mixedGroupIndex = useMixedMatchingCandidates ? mixedGroupIndex : -1; + storedCandidates.push_back(result); + ++nStoredThisGroup; + } + }; - int ranking = 1; - for (auto& candidate : globalTracksVector) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) - candidate.matchRanking = ranking; - ranking += 1; + if (useMixedMatchingCandidates) { + int32_t groupIdx = 0; + for (const auto& candidatesGroup : mchTrackInfo.mixedMatchingCandidates) { + processGroup(candidatesGroup, groupIdx); + groupIdx += 1; + } + } else { + processGroup(mchTrackInfo.matchingCandidates, -1); } } } @@ -1464,7 +1669,8 @@ struct GlobalMuonMatching { TMUON const& muonTracks, TMFT const& mftTracks, CMFT const& mftCovs, - aod::FwdTrkCls const& clusters) + aod::FwdTrkCls const& clusters, + bool useMixedMatchingCandidates) { if (configMchRealign.cfgEnableMCHRealign.value) { runMuonRealignment(muonTracks, clusters); @@ -1474,66 +1680,30 @@ struct GlobalMuonMatching { runMftRealignment(mftTracks, mftCovs); } - std::vector taggedMuons; - getTaggedMuons(collisions, muonTracks, taggedMuons); - + MatchingCandidates newMatchingCandidates; if (configMatching.cfgCustomMatchingStrategy.value == 0) { if (hasActiveChi2Matching) { - MatchingCandidates newMatchingCandidates; - runChi2Matching(activeChi2FunctionName, activeChi2MatchingPlaneZ, mMatchingCandidates, newMatchingCandidates); - fillMatchingCandidates(newMatchingCandidates, taggedMuons); + runChi2Matching(activeChi2FunctionName, activeChi2MatchingPlaneZ, newMatchingCandidates, useMixedMatchingCandidates); } } else { if (hasActiveMlMatching) { - MatchingCandidates newMatchingCandidates; - runMlMatching(collisions, muonTracks, mftTracks, activeMlResponse, activeMlMatchingPlaneZ, mMatchingCandidates, newMatchingCandidates); - fillMatchingCandidates(newMatchingCandidates, taggedMuons); + runMlMatching(collisions, muonTracks, mftTracks, activeMlResponse, activeMlMatchingPlaneZ, newMatchingCandidates, useMixedMatchingCandidates); } } - } - void fillMatchingCandidates(const MatchingCandidates& matchingCandidates, - const std::vector& taggedMuons) - { - for (const auto& [mchIndex, candidates] : matchingCandidates) { + for (const auto& [mchIndex, candidates] : newMatchingCandidates) { if (candidates.empty()) { continue; } - bool isTagged = std::find(taggedMuons.begin(), taggedMuons.end(), mchIndex) != taggedMuons.end(); - - std::vector storedCandidates; - int nStored = 0; - for (const auto& candidate : candidates) { - if (configMatching.cfgMaxCandidatesPerMchTrack.value >= 0 && nStored >= configMatching.cfgMaxCandidatesPerMchTrack.value) { - break; - } - - int32_t candidateIndex = mMatchCandidateCounter; - globalMuonMatchCandidates( - mchIndex, - candidate.mftTrackId, - static_cast(candidate.matchChi2), - static_cast(candidate.matchScore), - static_cast(candidate.matchRanking), - isTagged); - mMatchCandidateCounter += 1; - - mMchTrackToCandidateIndices[mchIndex].push_back(candidateIndex); - storedCandidates.push_back(candidate); - nStored += 1; - } - - if (!storedCandidates.empty()) { - mMchTrackMatchingCandidates[mchIndex] = std::move(storedCandidates); - } + mStoredMatchingCandidates[mchIndex] = candidates; } } int32_t countStoredCandidatesForMchTrack(int64_t mchTrackIndex) const { - const auto candidateIterator = mMchTrackMatchingCandidates.find(mchTrackIndex); - if (candidateIterator == mMchTrackMatchingCandidates.end()) { + const auto candidateIterator = mStoredMatchingCandidates.find(mchTrackIndex); + if (candidateIterator == mStoredMatchingCandidates.end()) { return 0; } return static_cast(candidateIterator->second.size()); @@ -1575,18 +1745,12 @@ struct GlobalMuonMatching { const int64_t mchTrackIndex = track.globalIndex(); const int32_t gmmMchTrackId = mFwdTrackToGmmCandTrkIndex.at(mchTrackIndex); - const auto candidateIterator = mMchTrackMatchingCandidates.find(mchTrackIndex); + const auto candidateIterator = mStoredMatchingCandidates.find(mchTrackIndex); auto mchTrackParIt = mMchTrackPars.find(mchTrackIndex); - if (mchTrackParIt == mMchTrackPars.end()) { - // fill muon tracks table with original parameters - const TrackParExt trackPar{fwdtrackutils::getTrackParCovFwd(track, track)}; - fillBaseGmmCandFwdTrack(track, trackPar, gmmMchTrackId, -1.f, -1.f); - } else { - // fill muon tracks table with realignment parameters - fillBaseGmmCandFwdTrack(track, mchTrackParIt->second, gmmMchTrackId, -1.f, -1.f); - } + fillBaseGmmCandFwdTrack(track, mchTrackParIt->second, gmmMchTrackId, -1.f, -1.f, + isMchTrackTagged(mchTrackIndex)); - if (candidateIterator != mMchTrackMatchingCandidates.end()) { + if (candidateIterator != mStoredMatchingCandidates.end()) { for (const auto& candidate : candidateIterator->second) { auto mftTrackParIt = mMftTrackPars.find(candidate.mftTrackId); if (mftTrackParIt != mMftTrackPars.end()) { @@ -1608,38 +1772,20 @@ struct GlobalMuonMatching { parExt, gmmMchTrackId, track.chi2MatchMCHMFT(), - track.matchScoreMCHMFT()); + track.matchScoreMCHMFT(), + isMchTrackTagged(track.matchMCHTrackId())); } } } - template - void fillFwdTrkMatchCands(TMUON const& muonTracks) - { - std::vector empty{}; - for (const auto& muonTrack : muonTracks) { - if (static_cast(muonTrack.trackType()) <= GlobalTrackTypeMax) { - fwdTrkMatchCands(empty); - continue; - } - - const int64_t mchTrackIndex = muonTrack.globalIndex(); - const auto matchIterator = mMchTrackToCandidateIndices.find(mchTrackIndex); - if (matchIterator == mMchTrackToCandidateIndices.end() || matchIterator->second.empty()) { - fwdTrkMatchCands(empty); - } else { - fwdTrkMatchCands(matchIterator->second); - } - } - } - - void processData(MyEvents const& collisions, - aod::BCsWithTimestamps const& bcs, - MyMuons const& muonTracks, - MyMFTs const& mftTracks, - MyMFTCovariances const& mftCovs, - aod::FwdTrkCls const& clusters, - aod::AmbiguousFwdTracks const& ambFwdTracks) + void runCandidateProcessing(MyEvents const& collisions, + aod::BCsWithTimestamps const& bcs, + MyMuons const& muonTracks, + MyMFTs const& mftTracks, + MyMFTCovariances const& mftCovs, + aod::FwdTrkCls const& clusters, + aod::AmbiguousFwdTracks const& ambFwdTracks, + bool useMixedMatchingCandidates) { auto bc = bcs.begin(); initCcdb(bc); @@ -1650,25 +1796,47 @@ struct GlobalMuonMatching { mftTrackCovs[mftTrackCov.matchMFTTrackId()] = mftTrackCov.globalIndex(); } - mMatchCandidateCounter = 0; - mMchTrackToCandidateIndices.clear(); - mMchTrackMatchingCandidates.clear(); + mStoredMatchingCandidates.clear(); mFwdTrackToGmmCandTrkIndex.clear(); + mMchTrackInfos.clear(); LOGF(info, "Preparing candidates"); prepareMatchingCandidates(collisions, bcs, muonTracks, mftTracks, mftCovs); + if (useMixedMatchingCandidates) { + prepareEventMixingMatchingCandidates(collisions, bcs, muonTracks); + } LOGF(info, "Processing candidates"); - processMatchingCandidates(collisions, muonTracks, mftTracks, mftCovs, clusters); + processMatchingCandidates(collisions, muonTracks, mftTracks, mftCovs, clusters, useMixedMatchingCandidates); LOGF(info, "Filling tables"); - // fill table with track/candidates index mapping - fillFwdTrkMatchCands(muonTracks); - // fill track tables fillGmmCandidateFwdTracks(muonTracks, mftTracks, ambFwdTracks); } + void processData(MyEvents const& collisions, + aod::BCsWithTimestamps const& bcs, + MyMuons const& muonTracks, + MyMFTs const& mftTracks, + MyMFTCovariances const& mftCovs, + aod::FwdTrkCls const& clusters, + aod::AmbiguousFwdTracks const& ambFwdTracks) + { + runCandidateProcessing(collisions, bcs, muonTracks, mftTracks, mftCovs, clusters, ambFwdTracks, false); + } + + void processMixedData(MyEvents const& collisions, + aod::BCsWithTimestamps const& bcs, + MyMuons const& muonTracks, + MyMFTs const& mftTracks, + MyMFTCovariances const& mftCovs, + aod::FwdTrkCls const& clusters, + aod::AmbiguousFwdTracks const& ambFwdTracks) + { + runCandidateProcessing(collisions, bcs, muonTracks, mftTracks, mftCovs, clusters, ambFwdTracks, true); + } + PROCESS_SWITCH(GlobalMuonMatching, processData, "processData", true); + PROCESS_SWITCH(GlobalMuonMatching, processMixedData, "process event-mixed matching candidates", false); }; // Extends the fwdtracksrealign table with expression columns