diff --git a/PWGDQ/Tasks/qaMatching.cxx b/PWGDQ/Tasks/qaMatching.cxx index ea18c0aa9d7..d756244786f 100644 --- a/PWGDQ/Tasks/qaMatching.cxx +++ b/PWGDQ/Tasks/qaMatching.cxx @@ -22,6 +22,7 @@ #include "CCDB/BasicCCDBManager.h" #include "DataFormatsParameters/GRPMagField.h" #include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisDataModel.h" #include "Framework/AnalysisTask.h" #include "Framework/runDataProcessing.h" #include "GlobalTracking/MatchGlobalFwd.h" @@ -30,6 +31,7 @@ #include #include +#include #include #include #include @@ -37,6 +39,7 @@ #include #include #include +#include #include #include @@ -44,6 +47,63 @@ using namespace o2; using namespace o2::framework; using namespace o2::aod; +namespace qamatching +{ +DECLARE_SOA_COLUMN(ReducedEventId, reducedEventId, int64_t); +DECLARE_SOA_COLUMN(P, p, float); +DECLARE_SOA_COLUMN(Pt, pt, float); +DECLARE_SOA_COLUMN(Eta, eta, float); +DECLARE_SOA_COLUMN(Phi, phi, float); +DECLARE_SOA_COLUMN(MatchLabel, matchlabel, int8_t); +DECLARE_SOA_COLUMN(TrackId, trackid, int64_t); +DECLARE_SOA_COLUMN(MatchType, matchType, int8_t); +DECLARE_SOA_COLUMN(MatchScore, matchScore, float); +DECLARE_SOA_COLUMN(MatchRanking, matchRanking, int32_t); +DECLARE_SOA_COLUMN(MftMultiplicity, mftMultiplicity, int32_t); +DECLARE_SOA_COLUMN(TrackType, tracktype, int8_t); +DECLARE_SOA_COLUMN(MftMatchAttempts, mftmatchattempts, int32_t); +DECLARE_SOA_COLUMN(X_atVtx, x_atVtx, float); +DECLARE_SOA_COLUMN(Y_atVtx, y_atVtx, float); +DECLARE_SOA_COLUMN(Z_atVtx, z_atVtx, float); +DECLARE_SOA_COLUMN(Px_atVtx, px_atVtx, float); +DECLARE_SOA_COLUMN(Py_atVtx, py_atVtx, float); +DECLARE_SOA_COLUMN(Pz_atVtx, pz_atVtx, float); +DECLARE_SOA_COLUMN(ColX, colx, float); +DECLARE_SOA_COLUMN(ColY, coly, float); +DECLARE_SOA_COLUMN(ColZ, colz, float); +} // namespace qamatching + +namespace o2::aod +{ +DECLARE_SOA_TABLE(QaMatchingEvents, "AOD", "QAMEVT", + qamatching::ReducedEventId, + qamatching::MftMultiplicity, + qamatching::ColX, + qamatching::ColY, + qamatching::ColZ); +DECLARE_SOA_TABLE(QaMatchingMCHTrack, "AOD", "QAMCHTRK", + qamatching::ReducedEventId, + qamatching::TrackId, + qamatching::TrackType, + qamatching::P, + qamatching::Pt, + qamatching::Eta, + qamatching::Phi, + qamatching::MftMatchAttempts, + qamatching::X_atVtx, + qamatching::Y_atVtx, + qamatching::Z_atVtx, + qamatching::Px_atVtx, + qamatching::Py_atVtx, + qamatching::Pz_atVtx); +DECLARE_SOA_TABLE(QaMatchingCandidates, "AOD", "QAMCAND", + qamatching::ReducedEventId, + qamatching::MatchLabel, + qamatching::TrackId, + qamatching::P, qamatching::Pt, qamatching::Eta, qamatching::Phi, + qamatching::MatchType, qamatching::MatchScore, qamatching::MatchRanking); +} // namespace o2::aod + using MyEvents = soa::Join; using MyMuons = soa::Join; using MyMuonsMC = soa::Join; @@ -153,9 +213,9 @@ struct qaMatching { /// Variables for histograms configuration Configurable fNCandidatesMax{"cfgNCandidatesMax", 5, "Number of matching candidates stored for each muon track"}; Configurable fMftTrackMultiplicityMax{"cfgMftTrackMultiplicityMax", 1000, "Maximum number of MFT tracks per collision"}; + Configurable fQaMatchingAodDebug{"cfgQaMatchingAodDebug", 0, "If >0, print AO2D filling debug (0=off, N=max collisions)"}; double mBzAtMftCenter{0}; - o2::globaltracking::MatchGlobalFwd mExtrap; using MatchingFunc_t = std::function(const o2::dataformats::GlobalFwdTrack& mchtrack, const o2::track::TrackParCovFwd& mfttrack)>; @@ -343,6 +403,10 @@ struct qaMatching { std::unordered_map matchingHistos; matrix dimuonHistos; + Produces qaMatchingEvents; + Produces qaMatchingMCHTrack; + Produces qaMatchingCandidates; + struct EfficiencyPlotter { o2::framework::HistPtr p_num; o2::framework::HistPtr p_den; @@ -1576,6 +1640,28 @@ struct qaMatching { return trueMatchIndex; } + template + int GetTrueMatchIndexTrackType(TMUON const& muonTracks, + TMUONS const& muonTracksAll, + TMFTS const& mftTracks, + const std::vector& matchCandidatesVector, + const std::vector>& matchablePairs) + { + // Same definition as GetTrueMatchIndex, but require trackType-based IsMuon. + int trueMatchIndex = 0; + for (size_t i = 0; i < matchCandidatesVector.size(); i++) { + auto const& muonTrack = muonTracks.rawIteratorAt(matchCandidatesVector[i].globalTrackId); + if (!IsMuon(muonTrack, muonTracksAll, mftTracks)) { + continue; + } + if (IsTrueGlobalMatching(muonTrack, matchablePairs)) { + trueMatchIndex = i + 1; + break; + } + } + return trueMatchIndex; + } + template bool IsMuon(const TMCH& mchTrack, const TMFT& mftTrack) @@ -1623,6 +1709,7 @@ struct qaMatching { return kMatchTypeUndefined; auto const& mchTrack = muonTrack.template matchMCHTrack_as(); + auto const& mftTrack = muonTrack.template matchMFTTrack_as(); bool isPaired = IsMatchableMCH(mchTrack.globalIndex(), matchablePairs); bool isMuon = IsMuon(muonTrack, muonTracks, mftTracks); @@ -2034,8 +2121,8 @@ struct qaMatching { // find the index of the matching candidate that corresponds to the true match // index=1 corresponds to the leading candidate // index=0 means no candidate was found that corresponds to the true match - int trueMatchIndex = GetTrueMatchIndex(muonTracks, globalTracksVector, matchablePairs); - int trueMatchIndexProd = GetTrueMatchIndex(muonTracks, matchingCandidatesProd.at(mchIndex), matchablePairs); + int trueMatchIndex = GetTrueMatchIndexTrackType(muonTracks, muonTracks, mftTracks, globalTracksVector, matchablePairs); + int trueMatchIndexProd = GetTrueMatchIndexTrackType(muonTracks, muonTracks, mftTracks, matchingCandidatesProd.at(mchIndex), matchablePairs); float mcParticleDz = -1000; if (mchTrack.has_mcParticle()) { @@ -2755,6 +2842,105 @@ struct qaMatching { FillDimuonPlotsMC(collisionInfo, collisions, muonTracks, mftTracks); } + template + void FillQaMatchingAodTablesForCollision(TCOLLISION const& collision, + TMUON const& muonTracks, + const MatchingCandidates& matchingCandidates, + int8_t matchLabel, + int64_t reducedEventId) + { + for (const auto& [mchIndex, candidates] : matchingCandidates) { + if (candidates.empty()) { + continue; + } + + const auto& mchTrack = muonTracks.rawIteratorAt(mchIndex); + if (!IsGoodGlobalMuon(mchTrack, collision)) { + continue; + } + float p = mchTrack.p(); + float pt = mchTrack.pt(); + float eta = mchTrack.eta(); + float phi = mchTrack.phi(); + + for (const auto& candidate : candidates) { + qaMatchingCandidates( + reducedEventId, + matchLabel, + mchIndex, + p, + pt, + eta, + phi, + static_cast(candidate.matchType), + static_cast(candidate.matchScore), + static_cast(candidate.matchRanking)); + } + } + } + + template + void FillQaMatchingAodEventForCollision(const CollisionInfo& collisionInfo, + TCOLLISION const& collision, + int64_t reducedEventId, + int& debugCounter) + { + int32_t mftMultiplicity = static_cast(collisionInfo.mftTracks.size()); + qaMatchingEvents( + reducedEventId, + mftMultiplicity, + static_cast(collision.posX()), + static_cast(collision.posY()), + static_cast(collision.posZ())); + + if (fQaMatchingAodDebug > 0 && debugCounter < fQaMatchingAodDebug) { + LOGF(info, "[AO2D] reducedEvent=%" PRId64 " mftMult=%d", + reducedEventId, + static_cast(mftMultiplicity)); + debugCounter += 1; + } + } + + template + void FillQaMatchingMchTracksForCollision(const CollisionInfo& collisionInfo, + TCOLLISIONS const& collisions, + TCOLLISION const& collision, + TMUON const& muonTracks, + TMFT const& mftTracks, + TBC const& bcs, + int64_t reducedEventId) + { + std::unordered_set mchIds; + for (const auto& mchIndex : collisionInfo.mchTracks) { + mchIds.insert(mchIndex); + } + for (const auto& [mchIndex, candidates] : collisionInfo.matchingCandidates) { + (void)candidates; + mchIds.insert(mchIndex); + } + + for (const auto& mchIndex : mchIds) { + auto const& mchTrack = muonTracks.rawIteratorAt(mchIndex); + int mftMchMatchAttempts = GetMftMchMatchAttempts(collisions, bcs, mchTrack, mftTracks); + auto mchTrackAtVertex = VarManager::PropagateMuon(mchTrack, collision, VarManager::kToVertex); + qaMatchingMCHTrack( + reducedEventId, + mchIndex, + static_cast(mchTrack.trackType()), + static_cast(mchTrack.p()), + static_cast(mchTrack.pt()), + static_cast(mchTrack.eta()), + static_cast(mchTrack.phi()), + static_cast(mftMchMatchAttempts), + static_cast(mchTrackAtVertex.getX()), + static_cast(mchTrackAtVertex.getY()), + static_cast(mchTrackAtVertex.getZ()), + static_cast(mchTrackAtVertex.getPx()), + static_cast(mchTrackAtVertex.getPy()), + static_cast(mchTrackAtVertex.getPz())); + } + } + void processQAMC(MyEvents const& collisions, aod::BCsWithTimestamps const& bcs, MyMuonsMC const& muonTracks, @@ -2776,6 +2962,49 @@ struct qaMatching { mftTrackCovs[mftTrackCov.matchMFTTrackId()] = mftTrackCov.globalIndex(); } + std::unordered_map reducedEventIds; + int64_t reducedEventCounter = 0; + for (auto const& [collisionIndex, collisionInfo] : fCollisionInfos) { + reducedEventIds.emplace(collisionInfo.index, reducedEventCounter); + reducedEventCounter += 1; + } + + int debugCounter = 0; + for (auto const& [collisionIndex, collisionInfo] : fCollisionInfos) { + auto it = reducedEventIds.find(collisionInfo.index); + if (it == reducedEventIds.end()) { + continue; + } + int64_t reducedEventId = it->second; + auto collision = collisions.rawIteratorAt(collisionInfo.index); + FillQaMatchingAodEventForCollision(collisionInfo, collision, reducedEventId, debugCounter); + FillQaMatchingMchTracksForCollision(collisionInfo, collisions, collision, muonTracks, mftTracks, bcs, reducedEventId); + } + + struct AodLabel { + const char* name; + int8_t id; + }; + std::array aodLabels{{{"ProdAll", 0}, {"MatchXYPhiTanl", 1}, {"MatchXYPhiTanlMom", 2}}}; + for (const auto& aodLabel : aodLabels) { + if (matchingChi2Functions.find(aodLabel.name) == matchingChi2Functions.end()) { + LOGF(warn, "[AO2D] Chi2 label not found: %s", aodLabel.name); + continue; + } + debugCounter = 0; + for (auto const& [collisionIndex, collisionInfo] : fCollisionInfos) { + auto it = reducedEventIds.find(collisionInfo.index); + if (it == reducedEventIds.end()) { + continue; + } + int64_t reducedEventId = it->second; + MatchingCandidates matchingCandidates; + RunChi2Matching(collisions, bcs, muonTracks, mftTracks, mftCovs, aodLabel.name, collisionInfo.matchablePairs, collisionInfo.matchingCandidates, matchingCandidates); + auto collision = collisions.rawIteratorAt(collisionInfo.index); + FillQaMatchingAodTablesForCollision(collision, muonTracks, matchingCandidates, aodLabel.id, reducedEventId); + } + } + for (auto const& [collisionIndex, collisionInfo] : fCollisionInfos) { ProcessCollisionMC(collisionInfo, collisions, bcs, muonTracks, mftTracks, mftCovs); }