diff --git a/PWGDQ/Tasks/qaMatching.cxx b/PWGDQ/Tasks/qaMatching.cxx index b5c0bb2ef65..9560ec54751 100644 --- a/PWGDQ/Tasks/qaMatching.cxx +++ b/PWGDQ/Tasks/qaMatching.cxx @@ -84,6 +84,77 @@ using namespace o2; using namespace o2::framework; using namespace o2::aod; +namespace qamatching +{ +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(XAtVtx, xAtVtx, float); +DECLARE_SOA_COLUMN(YAtVtx, yAtVtx, float); +DECLARE_SOA_COLUMN(ZAtVtx, zAtVtx, float); +DECLARE_SOA_COLUMN(PxAtVtx, pxAtVtx, float); +DECLARE_SOA_COLUMN(PyAtVtx, pyAtVtx, float); +DECLARE_SOA_COLUMN(PzAtVtx, pzAtVtx, 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", + o2::soa::Index<>, + qamatching::MftMultiplicity, + qamatching::ColX, + qamatching::ColY, + qamatching::ColZ); +} // namespace o2::aod + +namespace qamatching +{ +DECLARE_SOA_INDEX_COLUMN_FULL(ReducedEvent, reducedEvent, int32_t, o2::aod::QaMatchingEvents, ""); +} // namespace qamatching + +namespace o2::aod +{ +DECLARE_SOA_TABLE(QaMatchingMCHTrack, "AOD", "QAMCHTRK", + qamatching::ReducedEventId, + qamatching::TrackId, + qamatching::TrackType, + qamatching::P, + qamatching::Pt, + qamatching::Eta, + qamatching::Phi, + qamatching::MftMatchAttempts, + qamatching::XAtVtx, + qamatching::YAtVtx, + qamatching::ZAtVtx, + qamatching::PxAtVtx, + qamatching::PyAtVtx, + qamatching::PzAtVtx); +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, + qamatching::XAtVtx, + qamatching::YAtVtx, + qamatching::ZAtVtx, + qamatching::PxAtVtx, + qamatching::PyAtVtx, + qamatching::PzAtVtx); +} // namespace o2::aod + using MyEvents = soa::Join; using MyMuons = soa::Join; using MyMuonsMC = soa::Join; @@ -139,6 +210,7 @@ struct QaMatching { static constexpr int ExtrapolationMethodMftFirstPoint = 2; static constexpr int ExtrapolationMethodVertex = 3; static constexpr int ExtrapolationMethodMftDca = 4; + static constexpr int DecayRankingDirect = 2; struct MatchingCandidate { int64_t collisionId{-1}; @@ -210,6 +282,8 @@ struct QaMatching { Configurable cfgNCandidatesMax{"cfgNCandidatesMax", 5, "Number of matching candidates stored for each muon track"}; Configurable cfgMftTrackMultiplicityMax{"cfgMftTrackMultiplicityMax", 1000, "Maximum number of MFT tracks per collision"}; + Configurable cfgQaMatchingAodDebug{"cfgQaMatchingAodDebug", 0, "If >0, print AO2D filling debug (0=off, N=max collisions)"}; + double mBzAtMftCenter{0}; o2::globaltracking::MatchGlobalFwd mExtrap; @@ -399,6 +473,10 @@ struct QaMatching { std::unordered_map matchingHistos; Matrix dimuonHistos; + Produces qaMatchingEvents; + Produces qaMatchingMCHTrack; + Produces qaMatchingCandidates; + struct EfficiencyPlotter { o2::framework::HistPtr pNum; o2::framework::HistPtr pDen; @@ -1696,7 +1774,7 @@ struct QaMatching { } else { result = (ranking == 1) ? kMatchTypeWrongLeading : kMatchTypeWrongNonLeading; } - } else if (decayRanking == 1) { + } else if (decayRanking == DecayRankingDirect) { result = (ranking == 1) ? kMatchTypeDecayLeading : kMatchTypeDecayNonLeading; } else { result = (ranking == 1) ? kMatchTypeFakeLeading : kMatchTypeFakeNonLeading; @@ -2329,6 +2407,17 @@ struct QaMatching { } } + if (cfgQaMatchingAodDebug > 0 && goodMatchFound && isTrueMatch) { + LOGF(info, + "[good&true] mchId=%lld trackType=%d p=%.3f pt=%.3f eta=%.3f phi=%.3f", + static_cast(mchTrack.globalIndex()), + static_cast(mchTrack.trackType()), + mchTrack.p(), + mchTrack.pt(), + mchTrack.eta(), + mchTrack.phi()); + } + // ---- MC ancestry ---- auto motherParticles = getMotherParticles(mchTrack); int motherPDG = 0; @@ -2790,6 +2879,110 @@ struct QaMatching { fillDimuonPlotsMc(collisionInfo, collisions, muonTracks, mftTracks); } + template + void fillQaMatchingAodTablesForCollision(TCOLLISION const& collision, + TMUON const& muonTracks, + const MatchingCandidates& matchingCandidates, + int8_t matchLabel, + int32_t reducedEventId) + { + for (const auto& [mchIndex, candidates] : matchingCandidates) { + if (candidates.empty()) { + continue; + } + + const auto& mchTrack = muonTracks.rawIteratorAt(mchIndex); + if (!isGoodGlobalMuon(mchTrack, collision)) { + continue; + } + + for (const auto& candidate : candidates) { + const auto& candidateTrack = muonTracks.rawIteratorAt(candidate.globalTrackId); + auto candidateTrackAtVertex = VarManager::PropagateMuon(candidateTrack, collision, VarManager::kToVertex); + qaMatchingCandidates( + reducedEventId, + matchLabel, + mchIndex, + static_cast(candidateTrack.p()), + static_cast(candidateTrack.pt()), + static_cast(candidateTrack.eta()), + static_cast(candidateTrack.phi()), + static_cast(candidate.matchType), + static_cast(candidate.matchScore), + static_cast(candidate.matchRanking), + static_cast(candidateTrackAtVertex.getX()), + static_cast(candidateTrackAtVertex.getY()), + static_cast(candidateTrackAtVertex.getZ()), + static_cast(candidateTrackAtVertex.getPx()), + static_cast(candidateTrackAtVertex.getPy()), + static_cast(candidateTrackAtVertex.getPz())); + } + } + } + + template + void fillQaMatchingAodEventForCollision(const CollisionInfo& collisionInfo, + TCOLLISION const& collision, + int32_t reducedEventId, + int& debugCounter) + { + int32_t mftMultiplicity = static_cast(collisionInfo.mftTracks.size()); + qaMatchingEvents( + mftMultiplicity, + static_cast(collision.posX()), + static_cast(collision.posY()), + static_cast(collision.posZ())); + + if (cfgQaMatchingAodDebug > 0 && debugCounter < cfgQaMatchingAodDebug) { + LOGF(info, "[AO2D] reducedEvent=%", reducedEventId); + debugCounter += 1; + } + } + + template + void fillQaMatchingMchTracksForCollision(const CollisionInfo& collisionInfo, + TCOLLISIONS const& collisions, + TCOLLISION const& collision, + TMUON const& muonTracks, + TMFT const& mftTracks, + TBC const& bcs, + int32_t reducedEventId) + { + std::vector mchIds; + for (const auto& mchIndex : collisionInfo.mchTracks) { + if (std::find(mchIds.begin(), mchIds.end(), mchIndex) == mchIds.end()) { + mchIds.emplace_back(mchIndex); + } + } + for (const auto& [mchIndex, candidates] : collisionInfo.matchingCandidates) { + (void)candidates; + if (std::find(mchIds.begin(), mchIds.end(), mchIndex) == mchIds.end()) { + mchIds.emplace_back(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, @@ -2811,6 +3004,49 @@ struct QaMatching { mftTrackCovs[mftTrackCov.matchMFTTrackId()] = mftTrackCov.globalIndex(); } + std::unordered_map reducedEventIds; + int32_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; + } + int32_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; + } + int32_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); }