From d7ef0ae65a85e5e5dc9927336e35d7c980272109 Mon Sep 17 00:00:00 2001 From: aferrero2707 Date: Mon, 25 May 2026 13:59:16 +0200 Subject: [PATCH] [PWGDQ] loop over MFT covariances for MFT-MCH matching The matching is performed by looping directly over the available MFT covariances, instead of the stored candidates. It allows to process AO2Ds in which all the MFT covariances are saved, regardless of the number of stored candidates. --- PWGDQ/Tasks/qaMatching.cxx | 896 +++++++++++++++++++++++-------------- 1 file changed, 555 insertions(+), 341 deletions(-) diff --git a/PWGDQ/Tasks/qaMatching.cxx b/PWGDQ/Tasks/qaMatching.cxx index 6f2056d9c24..88eb74a2b88 100644 --- a/PWGDQ/Tasks/qaMatching.cxx +++ b/PWGDQ/Tasks/qaMatching.cxx @@ -95,6 +95,7 @@ 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(MatchChi2, matchChi2, float); DECLARE_SOA_COLUMN(MatchScore, matchScore, float); DECLARE_SOA_COLUMN(MatchRanking, matchRanking, int32_t); DECLARE_SOA_COLUMN(MftMultiplicity, mftMultiplicity, int32_t); @@ -152,7 +153,7 @@ DECLARE_SOA_TABLE(QaMatchingCandidates, "AOD", "QAMCAND", qamatching::MatchLabel, qamatching::TrackId, qamatching::P, qamatching::Pt, qamatching::Eta, qamatching::Phi, - qamatching::MatchType, qamatching::MatchScore, qamatching::MatchRanking, + qamatching::MatchType, qamatching::MatchChi2, qamatching::MatchScore, qamatching::MatchRanking, qamatching::XAtVtx, qamatching::YAtVtx, qamatching::ZAtVtx, @@ -184,6 +185,19 @@ static float chi2ToScore(float chi2, int ndf, float chi2max) return static_cast(result); } +static void setMatchTypeAxisLabels(TAxis* axis) +{ + axis->SetBinLabel(1, "true (leading)"); + axis->SetBinLabel(2, "wrong (leading)"); + axis->SetBinLabel(3, "decay (leading)"); + axis->SetBinLabel(4, "fake (leading)"); + axis->SetBinLabel(5, "true (non leading)"); + axis->SetBinLabel(6, "wrong (non leading)"); + axis->SetBinLabel(7, "decay (non leading)"); + axis->SetBinLabel(8, "fake (non leading)"); + axis->SetBinLabel(9, "unknown"); +} + struct QaMatching { template @@ -217,12 +231,16 @@ struct QaMatching { static constexpr int ExtrapolationMethodVertex = 3; static constexpr int ExtrapolationMethodMftDca = 4; static constexpr int DecayRankingDirect = 2; + static constexpr float MatchingPlaneDefaultZ = -77.5; struct MatchingCandidate { int64_t collisionId{-1}; int64_t globalTrackId{-1}; int64_t muonTrackId{-1}; int64_t mftTrackId{-1}; + int trackType{-1}; + o2::track::TrackParCovFwd mftTrackProp; + o2::track::TrackParCovFwd mchTrackProp; double matchScore{-1}; double matchChi2{-1}; int matchRanking{-1}; @@ -233,6 +251,8 @@ struct QaMatching { MuonMatchType matchType{kMatchTypeUndefined}; }; + Configurable cfgIsMc{"cfgIsMc", true, "Wheter the processed data is from MC simulations"}; + //// Variables for selecting muon tracks Configurable cfgPMchLow{"cfgPMchLow", 0.0f, ""}; Configurable cfgPtMchLow{"cfgPtMchLow", 0.7f, ""}; @@ -634,6 +654,42 @@ struct QaMatching { } }; + struct MatchFeaturesHistos { + o2::framework::HistPtr hDeltaP; + o2::framework::HistPtr hDeltaPt; + o2::framework::HistPtr hDeltaX; + o2::framework::HistPtr hDeltaY; + o2::framework::HistPtr hDeltaPhi; + o2::framework::HistPtr hDeltaTanl; + o2::framework::HistPtr hDeltaEta; + o2::framework::HistPtr hRabs; + + MatchFeaturesHistos(std::string path, HistogramRegistry* registry, int numCandidates) + { + AxisSpec indexAxis = {numCandidates + 1, 0, static_cast(numCandidates + 1), "ranking index"}; + AxisSpec scoreAxis = {100, 0, 1, "match score"}; + int matchTypeMax = static_cast(kMatchTypeUndefined) + 1; + AxisSpec matchTypeAxis = {matchTypeMax, 0, static_cast(matchTypeMax), "match type"}; + AxisSpec dxAxis = {100, -10, 10, "#Deltax (cm)"}; + AxisSpec dyAxis = {100, -10, 10, "#Deltay (cm)"}; + AxisSpec dpAxis = {100, -10, 10, "#Deltap (GeV/c)"}; + AxisSpec dptAxis = {100, -1, 1, "#Deltap_{T} (GeV/c)"}; + AxisSpec dphiAxis = {100, -1, 1, "#Delta#phi (rad)"}; + AxisSpec dtanlAxis = {100, -10, 10, "#Deltatanl"}; + AxisSpec detaAxis = {100, -1, 1, "#Delta#eta"}; + AxisSpec rabsAxis = {100, 0, 100, "R_{abs}"}; + + hDeltaP = registry->add((path + "/deltaP").c_str(), "MFT-MCH #Deltap", {HistType::kTHnSparseF, {dpAxis, scoreAxis, indexAxis, matchTypeAxis}}); + hDeltaPt = registry->add((path + "/deltaPt").c_str(), "MFT-MCH #Deltap_{T}", {HistType::kTHnSparseF, {dptAxis, scoreAxis, indexAxis, matchTypeAxis}}); + hDeltaX = registry->add((path + "/deltaX").c_str(), "MFT-MCH #Deltax", {HistType::kTHnSparseF, {dxAxis, scoreAxis, indexAxis, matchTypeAxis}}); + hDeltaY = registry->add((path + "/deltaY").c_str(), "MFT-MCH #Deltay", {HistType::kTHnSparseF, {dyAxis, scoreAxis, indexAxis, matchTypeAxis}}); + hDeltaPhi = registry->add((path + "/deltaPhi").c_str(), "MFT-MCH #Delta#phi", {HistType::kTHnSparseF, {dphiAxis, scoreAxis, indexAxis, matchTypeAxis}}); + hDeltaTanl = registry->add((path + "/deltaTanl").c_str(), "MFT-MCH #DeltaTanl", {HistType::kTHnSparseF, {dtanlAxis, scoreAxis, indexAxis, matchTypeAxis}}); + hDeltaEta = registry->add((path + "/deltaEta").c_str(), "MFT-MCH #Delta#eta", {HistType::kTHnSparseF, {detaAxis, scoreAxis, indexAxis, matchTypeAxis}}); + hRabs = registry->add((path + "/Rabs").c_str(), "MFT-MCH R_{abs}", {HistType::kTHnSparseF, {rabsAxis, scoreAxis, indexAxis, matchTypeAxis}}); + } + }; + struct MatchRankingHistos { o2::framework::HistPtr hist; o2::framework::HistPtr histVsP; @@ -680,6 +736,8 @@ struct QaMatching { std::unique_ptr fMatchRankingPaired; std::unique_ptr fMatchRankingPairedGoodMCH; + std::unique_ptr fMatchFeaturesGoodMCH; + //- o2::framework::HistPtr fMissedMatches; o2::framework::HistPtr fMissedMatchesGoodMCH; @@ -717,10 +775,10 @@ struct QaMatching { o2::framework::HistPtr fTrueMatchChi2VsProd; //- - EfficiencyPlotter fMatchingPurityPlotter; - EfficiencyPlotter fPairingEfficiencyPlotter; - EfficiencyPlotter fMatchingEfficiencyPlotter; - EfficiencyPlotter fFakeMatchingEfficiencyPlotter; + std::unique_ptr fMatchingPurityPlotter; + std::unique_ptr fPairingEfficiencyPlotter; + std::unique_ptr fMatchingEfficiencyPlotter; + std::unique_ptr fFakeMatchingEfficiencyPlotter; HistogramRegistry* registry; @@ -728,191 +786,129 @@ struct QaMatching { HistogramRegistry* reg, bool createPdgMomHistograms, int mftMultMax, - int numCandidates) - : fMatchingPurityPlotter(path + "matching-purity/", "Matching purity", *reg, createPdgMomHistograms), - fPairingEfficiencyPlotter(path + "pairing-efficiency/", "Pairing efficiency", *reg, createPdgMomHistograms), - fMatchingEfficiencyPlotter(path + "matching-efficiency/", "Matching efficiency", *reg, createPdgMomHistograms), - fFakeMatchingEfficiencyPlotter(path + "fake-matching-efficiency/", "Fake matching efficiency", *reg, createPdgMomHistograms) + int numCandidates, + bool isMc) { registry = reg; AxisSpec pAxis = {100, 0, 100, "p (GeV/c)"}; AxisSpec ptAxis = {100, 0, 10, "p_{T} (GeV/c)"}; AxisSpec dzAxis = {100, 0, 50, "#Deltaz (cm)"}; AxisSpec indexAxis = {6, 0, 6, "ranking index"}; - - std::string histName = path + "matchRanking"; - std::string histTitle = "True match ranking"; - - fMatchRanking = std::make_unique(path + "matchRanking", "True match ranking", registry, mftMultMax, numCandidates); - fMatchRankingGoodMCH = std::make_unique(path + "matchRankingGoodMCH", "True match ranking (good MCH tracks)", registry, mftMultMax, numCandidates); - fMatchRankingPaired = std::make_unique(path + "matchRankingPaired", "True match ranking (paired MCH tracks)", registry, mftMultMax, numCandidates); - fMatchRankingPairedGoodMCH = std::make_unique(path + "matchRankingPairedGoodMCH", "True match ranking (good paired MCH tracks)", registry, mftMultMax, numCandidates); - - //- - AxisSpec missedMatchAxis = {5, 0, 5, ""}; - histName = path + "missedMatches"; - histTitle = "Missed matches"; - fMissedMatches = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {missedMatchAxis}}); - std::get>(fMissedMatches)->GetXaxis()->SetBinLabel(1, "not paired"); - std::get>(fMissedMatches)->GetXaxis()->SetBinLabel(2, "fake MCH"); - std::get>(fMissedMatches)->GetXaxis()->SetBinLabel(3, "not stored"); - histName = path + "missedMatchesGoodMCH"; - histTitle = "Missed matches - good MCH tracks"; - fMissedMatchesGoodMCH = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {missedMatchAxis}}); - std::get>(fMissedMatchesGoodMCH)->GetXaxis()->SetBinLabel(1, "not paired"); - std::get>(fMissedMatchesGoodMCH)->GetXaxis()->SetBinLabel(2, "fake MCH"); - std::get>(fMissedMatchesGoodMCH)->GetXaxis()->SetBinLabel(3, "not stored"); - - AxisSpec decayRankingAxis = {5, 0, 5, "decay ranking"}; - histName = path + "decayRankingGoodMatches"; - histTitle = "Decay ranking - good matches"; - fDecayRankingGoodMatches = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {decayRankingAxis}}); - histName = path + "decayRankingNonLeadingMatches"; - histTitle = "Decay ranking - non-leading matches"; - fDecayRankingNonLeadingMatches = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {decayRankingAxis}}); - histName = path + "decayRankingMissedMatches"; - histTitle = "Decay ranking - missed matches"; - fDecayRankingMissedMatches = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {decayRankingAxis}}); - - AxisSpec scoreGapAxis = {100, 0, 1, "match score difference"}; - histName = path + "scoreGapLeadingTrueMatches"; - histTitle = "Score gap between leading and subleading matches - good matches"; - fScoreGapLeadingTrueMatches = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {scoreGapAxis}}); - histName = path + "scoreGapNonLeadingTrueMatches"; - histTitle = "Score gap between leading and subleading matches - non-leading matches"; - fScoreGapNonLeadingTrueMatches = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {scoreGapAxis}}); - - //- AxisSpec chi2Axis = {100, 0, 100, "matching #chi^{2}/NDF"}; AxisSpec scoreAxis = {100, 0, 1, "matching score"}; - int matchTypeMax = static_cast(kMatchTypeUndefined) + 1; - AxisSpec matchTypeAxis = {matchTypeMax, 0, static_cast(matchTypeMax), ""}; - histName = path + "matchType"; - histTitle = "Match type"; - fMatchType = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {matchTypeAxis}}); - std::get>(fMatchType)->GetXaxis()->SetBinLabel(1, "true (leading)"); - std::get>(fMatchType)->GetXaxis()->SetBinLabel(2, "wrong (leading)"); - std::get>(fMatchType)->GetXaxis()->SetBinLabel(3, "decay (leading)"); - std::get>(fMatchType)->GetXaxis()->SetBinLabel(4, "fake (leading)"); - std::get>(fMatchType)->GetXaxis()->SetBinLabel(5, "true (non leading)"); - std::get>(fMatchType)->GetXaxis()->SetBinLabel(6, "wrong (non leading)"); - std::get>(fMatchType)->GetXaxis()->SetBinLabel(7, "decay (non leading)"); - std::get>(fMatchType)->GetXaxis()->SetBinLabel(8, "fake (non leading)"); - std::get>(fMatchType)->GetXaxis()->SetBinLabel(9, "undefined"); - histName = path + "matchTypeVsP"; - histTitle = "Match type vs. p"; - fMatchTypeVsP = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {pAxis, matchTypeAxis}}); - std::get>(fMatchTypeVsP)->GetYaxis()->SetBinLabel(1, "true (leading)"); - std::get>(fMatchTypeVsP)->GetYaxis()->SetBinLabel(2, "wrong (leading)"); - std::get>(fMatchTypeVsP)->GetYaxis()->SetBinLabel(3, "decay (leading)"); - std::get>(fMatchTypeVsP)->GetYaxis()->SetBinLabel(4, "fake (leading)"); - std::get>(fMatchTypeVsP)->GetYaxis()->SetBinLabel(5, "true (non leading)"); - std::get>(fMatchTypeVsP)->GetYaxis()->SetBinLabel(6, "wrong (non leading)"); - std::get>(fMatchTypeVsP)->GetYaxis()->SetBinLabel(7, "decay (non leading)"); - std::get>(fMatchTypeVsP)->GetYaxis()->SetBinLabel(8, "fake (non leading)"); - std::get>(fMatchTypeVsP)->GetYaxis()->SetBinLabel(9, "undefined"); - histName = path + "matchTypeVsPt"; - histTitle = "Match type vs. p_{T}"; - fMatchTypeVsPt = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {ptAxis, matchTypeAxis}}); - std::get>(fMatchTypeVsPt)->GetYaxis()->SetBinLabel(1, "true (leading)"); - std::get>(fMatchTypeVsPt)->GetYaxis()->SetBinLabel(2, "wrong (leading)"); - std::get>(fMatchTypeVsPt)->GetYaxis()->SetBinLabel(3, "decay (leading)"); - std::get>(fMatchTypeVsPt)->GetYaxis()->SetBinLabel(4, "fake (leading)"); - std::get>(fMatchTypeVsPt)->GetYaxis()->SetBinLabel(5, "true (non leading)"); - std::get>(fMatchTypeVsPt)->GetYaxis()->SetBinLabel(6, "wrong (non leading)"); - std::get>(fMatchTypeVsPt)->GetYaxis()->SetBinLabel(7, "decay (non leading)"); - std::get>(fMatchTypeVsPt)->GetYaxis()->SetBinLabel(8, "fake (non leading)"); - std::get>(fMatchTypeVsPt)->GetYaxis()->SetBinLabel(9, "undefined"); - - histName = path + "matchChi2VsType"; - histTitle = "Match #chi^{2} vs. match type"; - fMatchChi2VsType = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {matchTypeAxis, chi2Axis}}); - std::get>(fMatchChi2VsType)->GetXaxis()->SetBinLabel(1, "true (leading)"); - std::get>(fMatchChi2VsType)->GetXaxis()->SetBinLabel(2, "wrong (leading)"); - std::get>(fMatchChi2VsType)->GetXaxis()->SetBinLabel(3, "decay (leading)"); - std::get>(fMatchChi2VsType)->GetXaxis()->SetBinLabel(4, "fake (leading)"); - std::get>(fMatchChi2VsType)->GetXaxis()->SetBinLabel(5, "true (non leading)"); - std::get>(fMatchChi2VsType)->GetXaxis()->SetBinLabel(6, "wrong (non leading)"); - std::get>(fMatchChi2VsType)->GetXaxis()->SetBinLabel(7, "decay (non leading)"); - std::get>(fMatchChi2VsType)->GetXaxis()->SetBinLabel(8, "fake (non leading)"); - std::get>(fMatchChi2VsType)->GetXaxis()->SetBinLabel(9, "undefined"); - histName = path + "matchChi2VsTypeVsP"; - histTitle = "Match #chi^{2} vs. match type vs. p"; - fMatchChi2VsTypeVsP = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH3F, {pAxis, matchTypeAxis, chi2Axis}}); - std::get>(fMatchChi2VsTypeVsP)->GetYaxis()->SetBinLabel(1, "true (leading)"); - std::get>(fMatchChi2VsTypeVsP)->GetYaxis()->SetBinLabel(2, "wrong (leading)"); - std::get>(fMatchChi2VsTypeVsP)->GetYaxis()->SetBinLabel(3, "decay (leading)"); - std::get>(fMatchChi2VsTypeVsP)->GetYaxis()->SetBinLabel(4, "fake (leading)"); - std::get>(fMatchChi2VsTypeVsP)->GetYaxis()->SetBinLabel(5, "true (non leading)"); - std::get>(fMatchChi2VsTypeVsP)->GetYaxis()->SetBinLabel(6, "wrong (non leading)"); - std::get>(fMatchChi2VsTypeVsP)->GetYaxis()->SetBinLabel(7, "decay (non leading)"); - std::get>(fMatchChi2VsTypeVsP)->GetYaxis()->SetBinLabel(8, "fake (non leading)"); - std::get>(fMatchChi2VsTypeVsP)->GetYaxis()->SetBinLabel(9, "undefined"); - histName = path + "matchChi2VsTypeVsPt"; - histTitle = "Match #chi^{2} vs. match type vs. p_{T}"; - fMatchChi2VsTypeVsPt = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH3F, {ptAxis, matchTypeAxis, chi2Axis}}); - std::get>(fMatchChi2VsTypeVsPt)->GetYaxis()->SetBinLabel(1, "true (leading)"); - std::get>(fMatchChi2VsTypeVsPt)->GetYaxis()->SetBinLabel(2, "wrong (leading)"); - std::get>(fMatchChi2VsTypeVsPt)->GetYaxis()->SetBinLabel(3, "decay (leading)"); - std::get>(fMatchChi2VsTypeVsPt)->GetYaxis()->SetBinLabel(4, "fake (leading)"); - std::get>(fMatchChi2VsTypeVsPt)->GetYaxis()->SetBinLabel(5, "true (non leading)"); - std::get>(fMatchChi2VsTypeVsPt)->GetYaxis()->SetBinLabel(6, "wrong (non leading)"); - std::get>(fMatchChi2VsTypeVsPt)->GetYaxis()->SetBinLabel(7, "decay (non leading)"); - std::get>(fMatchChi2VsTypeVsPt)->GetYaxis()->SetBinLabel(8, "fake (non leading)"); - std::get>(fMatchChi2VsTypeVsPt)->GetYaxis()->SetBinLabel(9, "undefined"); - //- - histName = path + "matchScoreVsType"; - histTitle = "Match score vs. match type"; - fMatchScoreVsType = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {matchTypeAxis, scoreAxis}}); - std::get>(fMatchScoreVsType)->GetXaxis()->SetBinLabel(1, "true (leading)"); - std::get>(fMatchScoreVsType)->GetXaxis()->SetBinLabel(2, "wrong (leading)"); - std::get>(fMatchScoreVsType)->GetXaxis()->SetBinLabel(3, "decay (leading)"); - std::get>(fMatchScoreVsType)->GetXaxis()->SetBinLabel(4, "fake (leading)"); - std::get>(fMatchScoreVsType)->GetXaxis()->SetBinLabel(5, "true (non leading)"); - std::get>(fMatchScoreVsType)->GetXaxis()->SetBinLabel(6, "wrong (non leading)"); - std::get>(fMatchScoreVsType)->GetXaxis()->SetBinLabel(7, "decay (non leading)"); - std::get>(fMatchScoreVsType)->GetXaxis()->SetBinLabel(8, "fake (non leading)"); - std::get>(fMatchScoreVsType)->GetXaxis()->SetBinLabel(9, "undefined"); - histName = path + "matchScoreVsTypeVsP"; - histTitle = "Match score vs. match type vs. p"; - fMatchScoreVsTypeVsP = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH3F, {pAxis, matchTypeAxis, scoreAxis}}); - std::get>(fMatchScoreVsTypeVsP)->GetYaxis()->SetBinLabel(1, "true (leading)"); - std::get>(fMatchScoreVsTypeVsP)->GetYaxis()->SetBinLabel(2, "wrong (leading)"); - std::get>(fMatchScoreVsTypeVsP)->GetYaxis()->SetBinLabel(3, "decay (leading)"); - std::get>(fMatchScoreVsTypeVsP)->GetYaxis()->SetBinLabel(4, "fake (leading)"); - std::get>(fMatchScoreVsTypeVsP)->GetYaxis()->SetBinLabel(5, "true (non leading)"); - std::get>(fMatchScoreVsTypeVsP)->GetYaxis()->SetBinLabel(6, "wrong (non leading)"); - std::get>(fMatchScoreVsTypeVsP)->GetYaxis()->SetBinLabel(7, "decay (non leading)"); - std::get>(fMatchScoreVsTypeVsP)->GetYaxis()->SetBinLabel(8, "fake (non leading)"); - std::get>(fMatchScoreVsTypeVsP)->GetYaxis()->SetBinLabel(9, "undefined"); - histName = path + "matchScoreVsTypeVsPt"; - histTitle = "Match score vs. match type vs. p_{T}"; - fMatchScoreVsTypeVsPt = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH3F, {ptAxis, matchTypeAxis, scoreAxis}}); - std::get>(fMatchScoreVsTypeVsPt)->GetYaxis()->SetBinLabel(1, "true (leading)"); - std::get>(fMatchScoreVsTypeVsPt)->GetYaxis()->SetBinLabel(2, "wrong (leading)"); - std::get>(fMatchScoreVsTypeVsPt)->GetYaxis()->SetBinLabel(3, "decay (leading)"); - std::get>(fMatchScoreVsTypeVsPt)->GetYaxis()->SetBinLabel(4, "fake (leading)"); - std::get>(fMatchScoreVsTypeVsPt)->GetYaxis()->SetBinLabel(5, "true (non leading)"); - std::get>(fMatchScoreVsTypeVsPt)->GetYaxis()->SetBinLabel(6, "wrong (non leading)"); - std::get>(fMatchScoreVsTypeVsPt)->GetYaxis()->SetBinLabel(7, "decay (non leading)"); - std::get>(fMatchScoreVsTypeVsPt)->GetYaxis()->SetBinLabel(8, "fake (non leading)"); - std::get>(fMatchScoreVsTypeVsPt)->GetYaxis()->SetBinLabel(9, "undefined"); + + std::string histName; + std::string histTitle; + + fMatchFeaturesGoodMCH = std::make_unique(path + "matchFeaturesGoodMCH", registry, numCandidates); + + if (isMc) { + fMatchRanking = std::make_unique(path + "matchRanking", "True match ranking", registry, mftMultMax, numCandidates); + fMatchRankingGoodMCH = std::make_unique(path + "matchRankingGoodMCH", "True match ranking (good MCH tracks)", registry, mftMultMax, numCandidates); + fMatchRankingPaired = std::make_unique(path + "matchRankingPaired", "True match ranking (paired MCH tracks)", registry, mftMultMax, numCandidates); + fMatchRankingPairedGoodMCH = std::make_unique(path + "matchRankingPairedGoodMCH", "True match ranking (good paired MCH tracks)", registry, mftMultMax, numCandidates); + + fMatchingPurityPlotter = std::make_unique(path + "matching-purity/", "Matching purity", *reg, createPdgMomHistograms); + fPairingEfficiencyPlotter = std::make_unique(path + "pairing-efficiency/", "Pairing efficiency", *reg, createPdgMomHistograms); + fMatchingEfficiencyPlotter = std::make_unique(path + "matching-efficiency/", "Matching efficiency", *reg, createPdgMomHistograms); + fFakeMatchingEfficiencyPlotter = std::make_unique(path + "fake-matching-efficiency/", "Fake matching efficiency", *reg, createPdgMomHistograms); + + //- + AxisSpec missedMatchAxis = {5, 0, 5, ""}; + histName = path + "missedMatches"; + histTitle = "Missed matches"; + fMissedMatches = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {missedMatchAxis}}); + std::get>(fMissedMatches)->GetXaxis()->SetBinLabel(1, "not paired"); + std::get>(fMissedMatches)->GetXaxis()->SetBinLabel(2, "fake MCH"); + std::get>(fMissedMatches)->GetXaxis()->SetBinLabel(3, "not stored"); + histName = path + "missedMatchesGoodMCH"; + histTitle = "Missed matches - good MCH tracks"; + fMissedMatchesGoodMCH = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {missedMatchAxis}}); + std::get>(fMissedMatchesGoodMCH)->GetXaxis()->SetBinLabel(1, "not paired"); + std::get>(fMissedMatchesGoodMCH)->GetXaxis()->SetBinLabel(2, "fake MCH"); + std::get>(fMissedMatchesGoodMCH)->GetXaxis()->SetBinLabel(3, "not stored"); + + AxisSpec decayRankingAxis = {5, 0, 5, "decay ranking"}; + histName = path + "decayRankingGoodMatches"; + histTitle = "Decay ranking - good matches"; + fDecayRankingGoodMatches = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {decayRankingAxis}}); + histName = path + "decayRankingNonLeadingMatches"; + histTitle = "Decay ranking - non-leading matches"; + fDecayRankingNonLeadingMatches = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {decayRankingAxis}}); + histName = path + "decayRankingMissedMatches"; + histTitle = "Decay ranking - missed matches"; + fDecayRankingMissedMatches = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {decayRankingAxis}}); + + AxisSpec scoreGapAxis = {100, 0, 1, "match score difference"}; + histName = path + "scoreGapLeadingTrueMatches"; + histTitle = "Score gap between leading and subleading matches - good matches"; + fScoreGapLeadingTrueMatches = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {scoreGapAxis}}); + histName = path + "scoreGapNonLeadingTrueMatches"; + histTitle = "Score gap between leading and subleading matches - non-leading matches"; + fScoreGapNonLeadingTrueMatches = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {scoreGapAxis}}); + + //- + int matchTypeMax = static_cast(kMatchTypeUndefined) + 1; + AxisSpec matchTypeAxis = {matchTypeMax, 0, static_cast(matchTypeMax), ""}; + histName = path + "matchType"; + histTitle = "Match type"; + fMatchType = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {matchTypeAxis}}); + setMatchTypeAxisLabels(std::get>(fMatchType)->GetXaxis()); + histName = path + "matchTypeVsP"; + histTitle = "Match type vs. p"; + fMatchTypeVsP = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {pAxis, matchTypeAxis}}); + setMatchTypeAxisLabels(std::get>(fMatchTypeVsP)->GetYaxis()); + histName = path + "matchTypeVsPt"; + histTitle = "Match type vs. p_{T}"; + fMatchTypeVsPt = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {ptAxis, matchTypeAxis}}); + setMatchTypeAxisLabels(std::get>(fMatchTypeVsPt)->GetYaxis()); + + histName = path + "matchChi2VsType"; + histTitle = "Match #chi^{2} vs. match type"; + fMatchChi2VsType = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {matchTypeAxis, chi2Axis}}); + setMatchTypeAxisLabels(std::get>(fMatchChi2VsType)->GetXaxis()); + histName = path + "matchChi2VsTypeVsP"; + histTitle = "Match #chi^{2} vs. match type vs. p"; + fMatchChi2VsTypeVsP = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH3F, {pAxis, matchTypeAxis, chi2Axis}}); + setMatchTypeAxisLabels(std::get>(fMatchChi2VsTypeVsP)->GetYaxis()); + histName = path + "matchChi2VsTypeVsPt"; + histTitle = "Match #chi^{2} vs. match type vs. p_{T}"; + fMatchChi2VsTypeVsPt = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH3F, {ptAxis, matchTypeAxis, chi2Axis}}); + setMatchTypeAxisLabels(std::get>(fMatchChi2VsTypeVsPt)->GetYaxis()); + //- + histName = path + "matchScoreVsType"; + histTitle = "Match score vs. match type"; + fMatchScoreVsType = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {matchTypeAxis, scoreAxis}}); + setMatchTypeAxisLabels(std::get>(fMatchScoreVsType)->GetXaxis()); + histName = path + "matchScoreVsTypeVsP"; + histTitle = "Match score vs. match type vs. p"; + fMatchScoreVsTypeVsP = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH3F, {pAxis, matchTypeAxis, scoreAxis}}); + setMatchTypeAxisLabels(std::get>(fMatchScoreVsTypeVsP)->GetYaxis()); + histName = path + "matchScoreVsTypeVsPt"; + histTitle = "Match score vs. match type vs. p_{T}"; + fMatchScoreVsTypeVsPt = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH3F, {ptAxis, matchTypeAxis, scoreAxis}}); + setMatchTypeAxisLabels(std::get>(fMatchScoreVsTypeVsPt)->GetYaxis()); + } AxisSpec prodScoreAxis = {100, 0, 1, "matching score (prod)"}; histName = path + "matchScoreVsProd"; histTitle = "Match score vs. production"; fMatchScoreVsProd = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {prodScoreAxis, scoreAxis}}); - histName = path + "trueMatchScoreVsProd"; - histTitle = "Match score vs. production - true match"; - fTrueMatchScoreVsProd = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {prodScoreAxis, scoreAxis}}); + if (isMc) { + histName = path + "trueMatchScoreVsProd"; + histTitle = "Match score vs. production - true match"; + fTrueMatchScoreVsProd = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {prodScoreAxis, scoreAxis}}); + } AxisSpec prodChi2Axis = {100, 0, 100, "matching #chi^{2}/NDF (prod)"}; histName = path + "matchChi2VsProd"; histTitle = "Match #chi^{2} vs. production"; fMatchChi2VsProd = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {prodChi2Axis, chi2Axis}}); - histName = path + "trueMatchChi2VsProd"; - histTitle = "Match #chi^{2} vs. production - true match"; - fTrueMatchChi2VsProd = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {{100, 0, 10, "matching #chi^{2} (prod)"}, {100, 0, 10, "matching #chi^{2}"}}}); + if (isMc) { + histName = path + "trueMatchChi2VsProd"; + histTitle = "Match #chi^{2} vs. production - true match"; + fTrueMatchChi2VsProd = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {{100, 0, 10, "matching #chi^{2} (prod)"}, {100, 0, 10, "matching #chi^{2}"}}}); + } } }; std::unique_ptr fChi2MatchingPlotter; @@ -948,7 +944,7 @@ struct QaMatching { } } - void createMatchingHistosMc() + void createMatchingHistos() { AxisSpec chi2Axis = {1000, 0, 1000, "chi^{2}"}; AxisSpec chi2AxisSmall = {200, 0, 100, "chi^{2}"}; @@ -956,7 +952,7 @@ struct QaMatching { AxisSpec pTAxis = {100, 0, 10, "p_{T} (GeV/c)"}; AxisSpec etaAxis = {100, -4, -2, "#eta"}; AxisSpec phiAxis = {90, -180, 180, "#phi (degrees)"}; - std::string histPath = "matching/MC/"; + std::string histPath = cfgIsMc.value ? "matching/MC/" : "matching/"; AxisSpec trackPositionXAtMftAxis = {100, -15, 15, "MFT x (cm)"}; AxisSpec trackPositionYAtMftAxis = {100, -15, 15, "MFT y (cm)"}; @@ -966,18 +962,18 @@ struct QaMatching { registry.add((histPath + "selectedMCHTracksAtMFTTrue").c_str(), "Selected MCH tracks position at MFT end - true", {HistType::kTH2F, {trackPositionXAtMftAxis, trackPositionYAtMftAxis}}); registry.add((histPath + "selectedMCHTracksAtMFTFake").c_str(), "Selected MCH tracks position at MFT end - fake", {HistType::kTH2F, {trackPositionXAtMftAxis, trackPositionYAtMftAxis}}); - fChi2MatchingPlotter = std::make_unique(histPath + "Prod/", ®istryMatching, configQas.cfgCreatePdgMomHistograms, cfgMftTrackMultiplicityMax, cfgNCandidatesMax); + fChi2MatchingPlotter = std::make_unique(histPath + "Prod/", ®istryMatching, configQas.cfgCreatePdgMomHistograms, cfgMftTrackMultiplicityMax, cfgNCandidatesMax, cfgIsMc.value); int registryIndex = 0; for (const auto& [label, func] : matchingChi2Functions) { - fMatchingPlotters[label] = std::make_unique(histPath + label + "/", registryMatchingVec[registryIndex], configQas.cfgCreatePdgMomHistograms, cfgMftTrackMultiplicityMax, cfgNCandidatesMax); + fMatchingPlotters[label] = std::make_unique(histPath + label + "/", registryMatchingVec[registryIndex], configQas.cfgCreatePdgMomHistograms, cfgMftTrackMultiplicityMax, cfgNCandidatesMax, cfgIsMc.value); registryIndex += 1; } for (const auto& [label, response] : matchingMlResponses) { - fMatchingPlotters[label] = std::make_unique(histPath + label + "/", (registryMatchingVec[registryIndex]), configQas.cfgCreatePdgMomHistograms, cfgMftTrackMultiplicityMax, cfgNCandidatesMax); + fMatchingPlotters[label] = std::make_unique(histPath + label + "/", (registryMatchingVec[registryIndex]), configQas.cfgCreatePdgMomHistograms, cfgMftTrackMultiplicityMax, cfgNCandidatesMax, cfgIsMc.value); registryIndex += 1; } - fTaggedMuonsMatchingPlotter = std::make_unique(histPath + "Tagged/", ®istryMatching, configQas.cfgCreatePdgMomHistograms, cfgMftTrackMultiplicityMax, cfgNCandidatesMax); + fTaggedMuonsMatchingPlotter = std::make_unique(histPath + "Tagged/", ®istryMatching, configQas.cfgCreatePdgMomHistograms, cfgMftTrackMultiplicityMax, cfgNCandidatesMax, cfgIsMc.value); } void createDimuonHistos() @@ -1203,7 +1199,7 @@ struct QaMatching { std::get>(matchMethodsHist)->GetXaxis()->SetBinLabel(iLabel + 1, matchMethodLabels[iLabel].c_str()); } - createMatchingHistosMc(); + createMatchingHistos(); createDimuonHistos(); } @@ -1318,35 +1314,31 @@ struct QaMatching { return isGoodMft(mftTrack, cfgTrackChi2MftUp, cfgTrackNClustMftLow); } - template - bool isGoodGlobalMatching(const TMUON& muonTrack, - double matchingScore, + bool isGoodGlobalMatching(const MatchingCandidate& candidate, double matchingScoreCut) { - if (static_cast(muonTrack.trackType()) > GlobalTrackTypeMax) + if (static_cast(candidate.trackType) > GlobalTrackTypeMax) return false; // MFT-MCH matching score cut - if (matchingScore < matchingScoreCut) + if (candidate.matchScore < matchingScoreCut) return false; return true; } - template - bool isGoodGlobalMatching(const TMUON& muonTrack, double matchingScore) + bool isGoodGlobalMatching(const MatchingCandidate& candidate) { - return isGoodGlobalMatching(muonTrack, matchingScore, cfgMatchingChi2ScoreMftMchLow); + return isGoodGlobalMatching(candidate, cfgMatchingChi2ScoreMftMchLow); } - template - bool isTrueGlobalMatching(const TMUON& muonTrack, const std::vector>& matchablePairs) + bool isTrueGlobalMatching(const MatchingCandidate& candidate, const std::vector>& matchablePairs) { - if (static_cast(muonTrack.trackType()) > GlobalTrackTypeMax) + if (candidate.trackType > GlobalTrackTypeMax) return false; - int64_t mchTrackId = static_cast(muonTrack.matchMCHTrackId()); - int64_t mftTrackId = static_cast(muonTrack.matchMFTTrackId()); + int64_t mchTrackId = candidate.muonTrackId; + int64_t mftTrackId = candidate.mftTrackId; std::pair trackIndexes = std::make_pair(mchTrackId, mftTrackId); @@ -1599,30 +1591,30 @@ struct QaMatching { } template - o2::dataformats::GlobalFwdTrack propagateToVertexMft(const TMFT& muon, + o2::dataformats::GlobalFwdTrack propagateToVertexMft(const TMFT& mftTrack, const TMCH& mchTrack, const C& collision) { // extrapolation with MCH tools auto mchTrackAtMFT = mExtrap.FwdtoMCH(fwdToTrackPar(mchTrack)); - o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchTrackAtMFT, muon.z()); + o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchTrackAtMFT, mftTrack.z()); - auto muonTrackProp = mExtrap.FwdtoMCH(fwdToTrackPar(muon)); + auto mftTrackProp = mExtrap.FwdtoMCH(fwdToTrackPar(mftTrack)); // update global track momentum from the MCH track - double pRatio = muonTrackProp.p() / mchTrackAtMFT.p(); - double newInvBendMom = muonTrackProp.getInverseBendingMomentum() * pRatio; - muonTrackProp.setInverseBendingMomentum(newInvBendMom); - muonTrackProp.setCharge(mchTrackAtMFT.getCharge()); + double pRatio = mftTrackProp.p() / mchTrackAtMFT.p(); + double newInvBendMom = mftTrackProp.getInverseBendingMomentum() * pRatio; + mftTrackProp.setInverseBendingMomentum(newInvBendMom); + mftTrackProp.setCharge(mchTrackAtMFT.getCharge()); - o2::mch::TrackExtrap::extrapToVertex(muonTrackProp, + o2::mch::TrackExtrap::extrapToVertex(mftTrackProp, collision.posX(), collision.posY(), collision.posZ(), collision.covXX(), collision.covYY()); - return mExtrap.MCHtoFwd(muonTrackProp); + return mExtrap.MCHtoFwd(mftTrackProp); } template @@ -1725,9 +1717,7 @@ struct QaMatching { } } - template - int getTrueMatchIndex(TMUON const& muonTracks, - const std::vector& matchCandidatesVector, + int getTrueMatchIndex(const std::vector& matchCandidatesVector, const std::vector>& matchablePairs) { // find the index of the matching candidate that corresponds to the true match @@ -1735,9 +1725,7 @@ struct QaMatching { // index=0 means no candidate was found that corresponds to the true match int trueMatchIndex = 0; for (size_t i = 0; i < matchCandidatesVector.size(); i++) { - auto const& muonTrack = muonTracks.rawIteratorAt(matchCandidatesVector[i].globalTrackId); - - if (isTrueGlobalMatching(muonTrack, matchablePairs)) { + if (isTrueGlobalMatching(matchCandidatesVector[i], matchablePairs)) { trueMatchIndex = i + 1; break; } @@ -1781,20 +1769,20 @@ struct QaMatching { return isMuon(mchTrack, mftTrack); } - template - MuonMatchType getMatchType(const TMUON& muonTrack, - TMUONS const& /*muonTracks*/, + template + MuonMatchType getMatchType(const MatchingCandidate& candidate, + TMUONS const& muonTracks, TMFTS const& mftTracks, const std::vector>& matchablePairs, int ranking) { - if (static_cast(muonTrack.trackType()) > GlobalTrackTypeMax) + if (candidate.trackType > GlobalTrackTypeMax) return kMatchTypeUndefined; - auto const& mchTrack = muonTrack.template matchMCHTrack_as(); + auto const& mchTrack = muonTracks.rawIteratorAt(candidate.muonTrackId); - bool isPairable = isMatchableMch(mchTrack.globalIndex(), matchablePairs); - bool isTrueMatch = isTrueGlobalMatching(muonTrack, matchablePairs); + bool isPairable = isMatchableMch(candidate.muonTrackId, matchablePairs); + bool isTrueMatch = isTrueGlobalMatching(candidate, matchablePairs); int decayRanking = getDecayRanking(mchTrack, mftTracks); MuonMatchType result{kMatchTypeUndefined}; @@ -1973,11 +1961,12 @@ struct QaMatching { return attempts; } - template + template void fillCollisions(EVT const& collisions, BC const& bcs, TMUON const& muonTracks, TMFT const& mftTracks, + MyMFTCovariances const& mftCovs, CollisionInfos& collisionInfos) { collisionInfos.clear(); @@ -2009,8 +1998,10 @@ struct QaMatching { collisionInfo.bc = bc.globalBC(); collisionInfo.zVertex = collision.posZ(); - if (collisionInfo.matchablePairs.empty()) { - fillMatchablePairs(collisionInfo, muonTracks, mftTracks); + if constexpr (isMC) { + if (collisionInfo.matchablePairs.empty()) { + fillMatchablePairs(collisionInfo, muonTracks, mftTracks); + } } if (static_cast(muonTrack.trackType()) > GlobalTrackTypeMax) { @@ -2027,6 +2018,18 @@ struct QaMatching { auto const& mftTrack = muonTrack.template matchMFTTrack_as(); int64_t mftTrackIndex = mftTrack.globalIndex(); + // get MFT track covariances + if (mftTrackCovs.count(mftTrack.globalIndex()) < 1) { + continue; + } + auto const& mftTrackCov = mftCovs.rawIteratorAt(mftTrackCovs[mftTrack.globalIndex()]); + + // propagate MCH and MFT tracks to matching plane + auto mchTrackProp = fwdToTrackPar(mchTrack, mchTrack); + mchTrackProp = propagateToMatchingPlaneMch(mchTrack, mftTrack, mftTrackCov, collision, MatchingPlaneDefaultZ, 0); + auto mftTrackProp = fwdToTrackPar(mftTrack, mftTrackCov); + mftTrackProp = propagateToMatchingPlaneMft(mchTrack, mftTrack, mftTrackCov, collision, MatchingPlaneDefaultZ, 0); + // check if a vector of global muon candidates is already available for the current MCH index // if not, initialize a new one and add the current global muon track // bool globalMuonTrackFound = false; @@ -2037,6 +2040,9 @@ struct QaMatching { muonTrackIndex, mchTrackIndex, mftTrackIndex, + static_cast(muonTrack.trackType()), + mftTrackProp, + mchTrackProp, matchScore, matchChi2, -1, @@ -2051,6 +2057,9 @@ struct QaMatching { muonTrackIndex, mchTrackIndex, mftTrackIndex, + static_cast(muonTrack.trackType()), + mftTrackProp, + mchTrackProp, matchScore, matchChi2, -1, @@ -2100,11 +2109,13 @@ struct QaMatching { auto mftMchMatchAttempts = getMftMchMatchAttempts(collisions, bcs, mchTrack, mftTracks); int ranking = 1; for (auto& candidate : globalTracksVector) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) - const auto& muonTrack = muonTracks.rawIteratorAt(candidate.globalTrackId); - candidate.matchRanking = ranking; candidate.matchRankingProd = ranking; - candidate.matchType = getMatchType(muonTrack, muonTracks, mftTracks, collisionInfo.matchablePairs, ranking); + if constexpr (isMC) { + candidate.matchType = getMatchType(candidate, muonTracks, mftTracks, collisionInfo.matchablePairs, ranking); + } else { + candidate.matchType = kMatchTypeUndefined; + } candidate.mftMchMatchAttempts = mftMchMatchAttempts; ranking += 1; } @@ -2112,6 +2123,65 @@ struct QaMatching { } } + template + void fillMatchingPlots(C const& collision, + TMUON const& muonTracks, + const MatchingCandidates& matchingCandidates, + MatchingPlotter* plotter) + { + // ==================================== + // Matching properties + for (const auto& [mchIndex, globalTracksVector] : matchingCandidates) { + if (globalTracksVector.size() < 1) + continue; + + // get the standalone MCH track + auto const& mchTrack = muonTracks.rawIteratorAt(mchIndex); + + // skip global muon tracks that do not pass the MCH and MFT quality cuts + if (!isGoodGlobalMuon(mchTrack, collision)) + continue; + + for (const auto& candidate : globalTracksVector) { + double dp = candidate.mchTrackProp.getP() - candidate.mftTrackProp.getP(); + double dpt = candidate.mchTrackProp.getPt() - candidate.mftTrackProp.getPt(); + double dx = candidate.mchTrackProp.getX() - candidate.mftTrackProp.getX(); + double dy = candidate.mchTrackProp.getY() - candidate.mftTrackProp.getY(); + double dphi = candidate.mchTrackProp.getPhi() - candidate.mftTrackProp.getPhi(); + double dtanl = candidate.mchTrackProp.getTanl() - candidate.mftTrackProp.getTanl(); + double deta = candidate.mchTrackProp.getEta() - candidate.mftTrackProp.getEta(); + int matchType = static_cast(candidate.matchType); + std::get>(plotter->fMatchFeaturesGoodMCH->hDeltaP)->Fill(dp, candidate.matchScore, candidate.matchRanking, matchType); + std::get>(plotter->fMatchFeaturesGoodMCH->hDeltaPt)->Fill(dpt, candidate.matchScore, candidate.matchRanking, matchType); + std::get>(plotter->fMatchFeaturesGoodMCH->hDeltaX)->Fill(dx, candidate.matchScore, candidate.matchRanking, matchType); + std::get>(plotter->fMatchFeaturesGoodMCH->hDeltaY)->Fill(dy, candidate.matchScore, candidate.matchRanking, matchType); + std::get>(plotter->fMatchFeaturesGoodMCH->hDeltaPhi)->Fill(dphi, candidate.matchScore, candidate.matchRanking, matchType); + std::get>(plotter->fMatchFeaturesGoodMCH->hDeltaTanl)->Fill(dtanl, candidate.matchScore, candidate.matchRanking, matchType); + std::get>(plotter->fMatchFeaturesGoodMCH->hDeltaEta)->Fill(deta, candidate.matchScore, candidate.matchRanking, matchType); + std::get>(plotter->fMatchFeaturesGoodMCH->hRabs)->Fill(mchTrack.rAtAbsorberEnd(), candidate.matchScore, candidate.matchRanking, matchType); + } + } + + // ==================================== + // Matching chi2 and score vs. production values + for (const auto& [mchIndex, globalTracksVector] : matchingCandidates) { + if (globalTracksVector.size() < 1) + continue; + + // loop over candidates + for (const auto& candidate : globalTracksVector) { + float matchChi2 = candidate.matchChi2; + float matchScore = candidate.matchScore; + + float matchChi2Prod = candidate.matchChi2Prod; + float matchScoreProd = candidate.matchScoreProd; + + std::get>(plotter->fMatchScoreVsProd)->Fill(matchScoreProd, matchScore); + std::get>(plotter->fMatchChi2VsProd)->Fill(matchChi2Prod, matchChi2); + } + } + } + template void fillMatchingPlotsMc(C const& collision, const CollisionInfo& collisionInfo, @@ -2155,8 +2225,9 @@ 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 = getTrueMatchIndex(globalTracksVector, matchablePairs); + const auto prodCandidatesIt = matchingCandidatesProd.find(mchIndex); + int trueMatchIndexProd = (prodCandidatesIt != matchingCandidatesProd.end()) ? getTrueMatchIndex(prodCandidatesIt->second, matchablePairs) : 0; float mcParticleDz = -1000; if (mchTrack.has_mcParticle()) { @@ -2216,6 +2287,27 @@ struct QaMatching { std::get>(plotter->fMatchRankingPairedGoodMCH->histVsDeltaChi2)->Fill(dchi2, trueMatchIndex); } + if (isGoodMCH) { + for (const auto& candidate : globalTracksVector) { + double dp = candidate.mchTrackProp.getP() - candidate.mftTrackProp.getP(); + double dpt = candidate.mchTrackProp.getPt() - candidate.mftTrackProp.getPt(); + double dx = candidate.mchTrackProp.getX() - candidate.mftTrackProp.getX(); + double dy = candidate.mchTrackProp.getY() - candidate.mftTrackProp.getY(); + double dphi = candidate.mchTrackProp.getPhi() - candidate.mftTrackProp.getPhi(); + double dtanl = candidate.mchTrackProp.getTanl() - candidate.mftTrackProp.getTanl(); + double deta = candidate.mchTrackProp.getEta() - candidate.mftTrackProp.getEta(); + int matchType = static_cast(candidate.matchType); + std::get>(plotter->fMatchFeaturesGoodMCH->hDeltaP)->Fill(dp, candidate.matchScore, candidate.matchRanking, matchType); + std::get>(plotter->fMatchFeaturesGoodMCH->hDeltaPt)->Fill(dpt, candidate.matchScore, candidate.matchRanking, matchType); + std::get>(plotter->fMatchFeaturesGoodMCH->hDeltaX)->Fill(dx, candidate.matchScore, candidate.matchRanking, matchType); + std::get>(plotter->fMatchFeaturesGoodMCH->hDeltaY)->Fill(dy, candidate.matchScore, candidate.matchRanking, matchType); + std::get>(plotter->fMatchFeaturesGoodMCH->hDeltaPhi)->Fill(dphi, candidate.matchScore, candidate.matchRanking, matchType); + std::get>(plotter->fMatchFeaturesGoodMCH->hDeltaTanl)->Fill(dtanl, candidate.matchScore, candidate.matchRanking, matchType); + std::get>(plotter->fMatchFeaturesGoodMCH->hDeltaEta)->Fill(deta, candidate.matchScore, candidate.matchRanking, matchType); + std::get>(plotter->fMatchFeaturesGoodMCH->hRabs)->Fill(mchTrack.rAtAbsorberEnd(), candidate.matchScore, candidate.matchRanking, matchType); + } + } + if (trueMatchIndex == 0) { // missed matches if (!isPairedMCH) { @@ -2304,18 +2396,16 @@ struct QaMatching { if (globalTracksVector.size() < 1) continue; - int trueMatchIndex = getTrueMatchIndex(muonTracks, globalTracksVector, matchablePairs); + int trueMatchIndex = getTrueMatchIndex(globalTracksVector, matchablePairs); // loop over candidates int candidateIndex = 1; for (const auto& candidate : globalTracksVector) { - auto const& muonTrack = muonTracks.rawIteratorAt(candidate.globalTrackId); - - float matchScore = candidate.matchScore; float matchChi2 = candidate.matchChi2; + float matchScore = candidate.matchScore; - float matchChi2Prod = muonTrack.chi2MatchMCHMFT() / 5.f; - float matchScoreProd = chi2ToScore(muonTrack.chi2MatchMCHMFT(), 5, 50.f); + float matchChi2Prod = candidate.matchChi2Prod; + float matchScoreProd = candidate.matchScoreProd; std::get>(plotter->fMatchScoreVsProd)->Fill(matchScoreProd, matchScore); std::get>(plotter->fMatchChi2VsProd)->Fill(matchChi2Prod, matchChi2); @@ -2335,37 +2425,30 @@ struct QaMatching { if (globalTracksVector.size() < 1) continue; - // get the leading matching candidate - auto const& muonTrack = muonTracks.rawIteratorAt(globalTracksVector[0].globalTrackId); - double matchingScore = globalTracksVector[0].matchScore; - // get the standalone MCH and MFT tracks auto const& mchTrack = muonTracks.rawIteratorAt(mchIndex); - auto const& mftTrack = muonTrack.template matchMFTTrack_as(); // skip global muon tracks that do not pass the MCH and MFT quality cuts if (!isGoodGlobalMuon(mchTrack, collision)) continue; - if (!isGoodMft(mftTrack)) - continue; // skip candidates that do not pass the matching quality cuts - if (!isGoodGlobalMatching(muonTrack, matchingScore, matchingScoreCut)) + if (!isGoodGlobalMatching(globalTracksVector[0], matchingScoreCut)) continue; // check if the matching candidate is a true one - bool isTrueMatch = isTrueGlobalMatching(muonTrack, matchablePairs); + bool isTrueMatch = isTrueGlobalMatching(globalTracksVector[0], matchablePairs); // ---- MC ancestry ---- - auto motherParticles = getMotherParticles(muonTrack); + auto motherParticles = getMotherParticles(mchTrack); int motherPDG = 0; if (motherParticles.size() > 1) { motherPDG = motherParticles[1].first; } // fill matching purity plots - plotter->fMatchingPurityPlotter.fill(mchTrack, isTrueMatch); + plotter->fMatchingPurityPlotter->fill(mchTrack, isTrueMatch); if (configQas.cfgCreatePdgMomHistograms) { - plotter->fMatchingPurityPlotter.fill(mchTrack, motherPDG, isTrueMatch); + plotter->fMatchingPurityPlotter->fill(mchTrack, motherPDG, isTrueMatch); } } @@ -2389,15 +2472,10 @@ struct QaMatching { if (matchingCandidates.count(matchableMchIndex) > 0) { const auto& globalTracksVector = matchingCandidates.at(static_cast(matchableMchIndex)); if (!globalTracksVector.empty()) { - // get the leading matching candidate - auto const& muonTrack = muonTracks.rawIteratorAt(globalTracksVector[0].globalTrackId); - double matchingScore = globalTracksVector[0].matchScore; + // get the standalone MFT track index + auto mftIndex = globalTracksVector[0].mftTrackId; - // get the standalone MFT track - auto const& mftTrack = muonTrack.template matchMFTTrack_as(); - auto mftIndex = mftTrack.globalIndex(); - - goodMatchFound = isGoodGlobalMatching(muonTrack, matchingScore, matchingScoreCut); + goodMatchFound = isGoodGlobalMatching(globalTracksVector[0], matchingScoreCut); isTrueMatch = (mftIndex == matchableMftIndex); } } @@ -2421,17 +2499,17 @@ struct QaMatching { } // fill matching efficiency plots - plotter->fPairingEfficiencyPlotter.fill(mchTrack, goodMatchFound); + plotter->fPairingEfficiencyPlotter->fill(mchTrack, goodMatchFound); if (configQas.cfgCreatePdgMomHistograms) { - plotter->fPairingEfficiencyPlotter.fill(mchTrack, motherPDG, goodMatchFound); + plotter->fPairingEfficiencyPlotter->fill(mchTrack, motherPDG, goodMatchFound); } - plotter->fMatchingEfficiencyPlotter.fill(mchTrack, (goodMatchFound && isTrueMatch)); + plotter->fMatchingEfficiencyPlotter->fill(mchTrack, (goodMatchFound && isTrueMatch)); if (configQas.cfgCreatePdgMomHistograms) { - plotter->fMatchingEfficiencyPlotter.fill(mchTrack, motherPDG, (goodMatchFound && isTrueMatch)); + plotter->fMatchingEfficiencyPlotter->fill(mchTrack, motherPDG, (goodMatchFound && isTrueMatch)); } - plotter->fFakeMatchingEfficiencyPlotter.fill(mchTrack, (goodMatchFound && !isTrueMatch)); + plotter->fFakeMatchingEfficiencyPlotter->fill(mchTrack, (goodMatchFound && !isTrueMatch)); if (configQas.cfgCreatePdgMomHistograms) { - plotter->fFakeMatchingEfficiencyPlotter.fill(mchTrack, motherPDG, (goodMatchFound && !isTrueMatch)); + plotter->fFakeMatchingEfficiencyPlotter->fill(mchTrack, motherPDG, (goodMatchFound && !isTrueMatch)); } } } @@ -2440,7 +2518,7 @@ struct QaMatching { void fillDimuonPlotsMc(const CollisionInfo& collisionInfo, C const& collisions, TMUON const& muonTracks, - TMFT const& /*mftTracks*/) + TMFT const& mftTracks) { std::vector muonPairs; std::vector globalMuonPairs; @@ -2476,12 +2554,10 @@ struct QaMatching { auto const& collision = collisions.rawIteratorAt(muon1.first); - auto const& muonTrack1 = muonTracks.rawIteratorAt(candidates1[0].globalTrackId); - auto const& muonTrack2 = muonTracks.rawIteratorAt(candidates2[0].globalTrackId); - auto matchScore1 = candidates1[0].matchScore; - auto matchScore2 = candidates2[0].matchScore; - auto const& mchTrack1 = muonTrack1.template matchMCHTrack_as(); - auto const& mchTrack2 = muonTrack2.template matchMCHTrack_as(); + auto const& mchTrack1 = muonTracks.rawIteratorAt(candidates1[0].muonTrackId); + auto const& mchTrack2 = muonTracks.rawIteratorAt(candidates2[0].muonTrackId); + auto const& mftTrack1 = mftTracks.rawIteratorAt(candidates1[0].mftTrackId); + auto const& mftTrack2 = mftTracks.rawIteratorAt(candidates2[0].mftTrackId); int sign1 = mchTrack1.sign(); int sign2 = mchTrack2.sign(); @@ -2503,12 +2579,12 @@ struct QaMatching { continue; } - bool goodGlobalMuonMatches = (isGoodGlobalMatching(muonTrack1, matchScore1) && isGoodGlobalMatching(muonTrack2, matchScore2)); + bool goodGlobalMuonMatches = (isGoodGlobalMatching(candidates1[0]) && isGoodGlobalMatching(candidates2[0])); double massMCH = getMuMuInvariantMass(propagateToVertexMch(mchTrack1, collision), propagateToVertexMch(mchTrack2, collision)); - double mass = getMuMuInvariantMass(propagateToVertexMch(muonTrack1, collision), - propagateToVertexMch(muonTrack2, collision)); + double mass = getMuMuInvariantMass(propagateToVertexMft(mftTrack1, mchTrack1, collision), + propagateToVertexMft(mftTrack2, mchTrack2, collision)); registryDimuon.get(HIST("dimuon/invariantMass_MuonKine_GlobalMuonCuts"))->Fill(massMCH); registryDimuon.get(HIST("dimuon/invariantMass_ScaledMftKine_GlobalMuonCuts"))->Fill(mass); registryDimuon.get(HIST("dimuon/MC/invariantMass_MuonKine_GlobalMuonCuts_vs_match_type"))->Fill(massMCH, matchType); @@ -2523,7 +2599,7 @@ struct QaMatching { } } - template + template void runChi2Matching(C const& collisions, BC const& bcs, TMUON const& muonTracks, @@ -2547,40 +2623,76 @@ struct QaMatching { return; auto matchingFunc = mMatchingFunctionMap.at(funcName); - for (const auto& [mchIndex, globalTracksVector] : matchingCandidates) { - auto const& mchTrack = muonTracks.rawIteratorAt(mchIndex); + for (const auto& muonTrack : muonTracks) { + if (!muonTrack.has_collision()) { + continue; + } - for (const auto& candidate : globalTracksVector) { - auto const& muonTrack = muonTracks.rawIteratorAt(candidate.globalTrackId); - if (!muonTrack.has_collision()) + // skip global tracks + if (static_cast(muonTrack.trackType()) <= GlobalTrackTypeMax) { + continue; + } + + int64_t mchIndex = muonTrack.globalIndex(); + auto collision = collisions.rawIteratorAt(muonTrack.collisionId()); + + const auto& collMch = collisions.rawIteratorAt(muonTrack.collisionId()); + const auto& bcMch = bcs.rawIteratorAt(collMch.bcId()); + + // loop over MFT track/covariance mapping + for (const auto& [mftIndex, mftCovIndex] : mftTrackCovs) { + auto const& mftTrack = mftTracks.rawIteratorAt(mftIndex); + auto const& mftTrackCov = mftCovs.rawIteratorAt(mftCovIndex); + + if (!mftTrack.has_collision()) { continue; + } + + const auto& collMft = collisions.rawIteratorAt(mftTrack.collisionId()); + const auto& bcMft = bcs.rawIteratorAt(collMft.bcId()); - auto collision = collisions.rawIteratorAt(muonTrack.collisionId()); + // check time compatibility between MCH and MFT tracks + int64_t deltaBc = static_cast(bcMft.globalBC()) - static_cast(bcMch.globalBC()); + double deltaBcNS = o2::constants::lhc::LHCBunchSpacingNS * deltaBc; + double deltaTrackTime = mftTrack.trackTime() - muonTrack.trackTime() + deltaBcNS; + double trackTimeResTot = mftTrack.trackTimeRes() + muonTrack.trackTimeRes(); - // get MCH and MFT standalone tracks - // auto mchTrack = muonTrack.template matchMCHTrack_as(); - auto const& mftTrack = muonTrack.template matchMFTTrack_as(); - if (mftTrackCovs.count(mftTrack.globalIndex()) < 1) { - // std::cout << std::format("Covariance matrix for MFT track #{} not found", mftTrack.globalIndex()) << std::endl; + if (std::fabs(deltaTrackTime) > trackTimeResTot) { continue; } - auto const& mftTrackCov = mftCovs.rawIteratorAt(mftTrackCovs[mftTrack.globalIndex()]); // get tracks parameters in O2 format auto mftTrackProp = fwdToTrackPar(mftTrack, mftTrackCov); - auto mchTrackProp = fwdToTrackPar(mchTrack, mchTrack); + auto mchTrackProp = fwdToTrackPar(muonTrack, muonTrack); if (matchingPlaneZ < 0.) { - mftTrackProp = propagateToMatchingPlaneMft(mchTrack, mftTrack, mftTrackCov, collision, matchingPlaneZ, extrapMethod); - mchTrackProp = propagateToMatchingPlaneMch(mchTrack, mftTrack, mftTrackCov, collision, matchingPlaneZ, extrapMethod); + mftTrackProp = propagateToMatchingPlaneMft(muonTrack, mftTrack, mftTrackCov, collision, matchingPlaneZ, extrapMethod); + mchTrackProp = propagateToMatchingPlaneMch(muonTrack, mftTrack, mftTrackCov, collision, matchingPlaneZ, extrapMethod); + } + + float matchChi2Prod = -1; + float matchScoreProd = -1; + int globalTrackType = -1; + int64_t globalTrackIndex = -1; + // search candidate from production + auto candidatesIt = matchingCandidates.find(mchIndex); + if (candidatesIt != matchingCandidates.end()) { + for (const auto& candidate : candidatesIt->second) { + if (candidate.mftTrackId != mftIndex) { + continue; + } + matchChi2Prod = candidate.matchChi2; + matchScoreProd = candidate.matchScore; + globalTrackType = candidate.trackType; + globalTrackIndex = candidate.globalTrackId; + break; + } } // run the chi2 matching function auto matchResult = matchingFunc(mchTrackProp, mftTrackProp); float matchChi2 = std::get<0>(matchResult) / std::get<1>(matchResult); float matchScore = chi2ToScore(std::get<0>(matchResult), std::get<1>(matchResult), 10.f * std::get<1>(matchResult)); - float matchChi2Prod = muonTrack.chi2MatchMCHMFT() / 5.f; - float matchScoreProd = chi2ToScore(muonTrack.chi2MatchMCHMFT(), 5, 50.f); // check if a vector of global muon candidates is already available for the current MCH index // if not, initialize a new one and add the current global muon track @@ -2588,9 +2700,12 @@ struct QaMatching { if (matchingCandidateIterator != newMatchingCandidates.end()) { matchingCandidateIterator->second.emplace_back(MatchingCandidate{ muonTrack.collisionId(), - candidate.globalTrackId, + globalTrackIndex, mchIndex, mftTrack.globalIndex(), + globalTrackType, + mftTrackProp, + mchTrackProp, matchScore, matchChi2, -1, @@ -2601,9 +2716,12 @@ struct QaMatching { } else { newMatchingCandidates[mchIndex].emplace_back(MatchingCandidate{ muonTrack.collisionId(), - candidate.globalTrackId, + globalTrackIndex, mchIndex, mftTrack.globalIndex(), + globalTrackType, + mftTrackProp, + mchTrackProp, matchScore, matchChi2, -1, @@ -2627,17 +2745,19 @@ struct QaMatching { auto mftMchMatchAttempts = getMftMchMatchAttempts(collisions, bcs, mchTrack, mftTracks); int ranking = 1; for (auto& candidate : globalTracksVector) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) - const auto& muonTrack = muonTracks.rawIteratorAt(candidate.globalTrackId); - candidate.matchRanking = ranking; - candidate.matchType = getMatchType(muonTrack, muonTracks, mftTracks, matchablePairs, ranking); + if constexpr (isMC) { + candidate.matchType = getMatchType(candidate, muonTracks, mftTracks, matchablePairs, ranking); + } else { + candidate.matchType = kMatchTypeUndefined; + } candidate.mftMchMatchAttempts = mftMchMatchAttempts; ranking += 1; } } } - template + template void runChi2Matching(C const& collisions, BC const& bcs, TMUON const& muonTracks, @@ -2668,10 +2788,10 @@ struct QaMatching { auto matchingPlaneZ = matchingPlanesZ[label]; auto extrapMethod = matchingExtrapMethod[label]; - runChi2Matching(collisions, bcs, muonTracks, mftTracks, mftCovs, funcName, matchingPlaneZ, extrapMethod, matchablePairs, matchingCandidates, newMatchingCandidates); + runChi2Matching(collisions, bcs, muonTracks, mftTracks, mftCovs, funcName, matchingPlaneZ, extrapMethod, matchablePairs, matchingCandidates, newMatchingCandidates); } - template + template void runMlMatching(C const& collisions, BC const& bcs, TMUON const& muonTracks, @@ -2688,28 +2808,47 @@ struct QaMatching { return; auto& mlResponse = mlIter->second; - for (const auto& [mchIndex, globalTracksVector] : matchingCandidates) { - auto const& mchTrack = muonTracks.rawIteratorAt(mchIndex); - for (const auto& candidate : globalTracksVector) { - auto const& muonTrack = muonTracks.rawIteratorAt(candidate.globalTrackId); - if (!muonTrack.has_collision()) + for (const auto& muonTrack : muonTracks) { + if (!muonTrack.has_collision()) { + continue; + } + + // skip global tracks + if (static_cast(muonTrack.trackType()) <= GlobalTrackTypeMax) { + continue; + } + + int64_t mchIndex = muonTrack.globalIndex(); + auto collision = collisions.rawIteratorAt(muonTrack.collisionId()); + + const auto& collMch = collisions.rawIteratorAt(muonTrack.collisionId()); + const auto& bcMch = bcs.rawIteratorAt(collMch.bcId()); + + // loop over MFT track/covariance mapping + for (const auto& [mftIndex, mftCovIndex] : mftTrackCovs) { + auto const& mftTrack = mftTracks.rawIteratorAt(mftIndex); + auto const& mftTrackCov = mftCovs.rawIteratorAt(mftCovIndex); + + if (!mftTrack.has_collision()) { continue; + } + + const auto& collMft = collisions.rawIteratorAt(mftTrack.collisionId()); + const auto& bcMft = bcs.rawIteratorAt(collMft.bcId()); - auto collision = collisions.rawIteratorAt(muonTrack.collisionId()); + // check time compatibility between MCH and MFT tracks + int64_t deltaBc = static_cast(bcMft.globalBC()) - static_cast(bcMch.globalBC()); + double deltaBcNS = o2::constants::lhc::LHCBunchSpacingNS * deltaBc; + double deltaTrackTime = mftTrack.trackTime() - muonTrack.trackTime() + deltaBcNS; + double trackTimeResTot = mftTrack.trackTimeRes() + muonTrack.trackTimeRes(); - // get MFT standalone track - auto const& mftTrack = muonTrack.template matchMFTTrack_as(); - if (mftTrackCovs.count(mftTrack.globalIndex()) < 1) { - // std::cout << std::format("Covariance matrix for MFT track #{} not found", mftTrack.globalIndex()) << std::endl; + if (std::fabs(deltaTrackTime) > trackTimeResTot) { continue; } - // std::cout << fmt::format("Getting covariance matrix for MFT track #{} -> {}", mftTrack.globalIndex(), mftTrackCovs[mftTrack.globalIndex()]) << std::endl; - auto const& mftTrackCov = mftCovs.rawIteratorAt(mftTrackCovs[mftTrack.globalIndex()]); - // std::cout << fmt::format("Covariance matrix for MFT track #{} retrieved", mftTrack.globalIndex()) << std::endl; // get tracks parameters in O2 format auto mftTrackProp = fwdToTrackPar(mftTrack, mftTrackCov); - auto mchTrackProp = fwdToTrackPar(mchTrack, mchTrack); + auto mchTrackProp = fwdToTrackPar(muonTrack, muonTrack); // extrapolate to the matching plane auto matchingPlaneZ = matchingPlanesZ[label]; @@ -2718,14 +2857,36 @@ struct QaMatching { mchTrackProp = propagateToZMch(mchTrackProp, matchingPlaneZ); } + float matchChi2Prod = -1; + float matchScoreProd = -1; + int globalTrackType = -1; + int64_t globalTrackIndex = -1; + // search candidate from production + auto candidatesIt = matchingCandidates.find(mchIndex); + if (candidatesIt != matchingCandidates.end()) { + for (const auto& candidate : candidatesIt->second) { + if (candidate.mftTrackId != mftIndex) { + continue; + } + matchChi2Prod = candidate.matchChi2; + matchScoreProd = candidate.matchScore; + globalTrackType = candidate.trackType; + globalTrackIndex = candidate.globalTrackId; + break; + } + } + // run the ML model std::vector output; - std::vector inputML = mlResponse.getInputFeatures(muonTrack, mftTrack, mchTrack, mftTrackProp, mchTrackProp, collision); + std::vector inputML; + if (globalTrackIndex < 0) { + inputML = mlResponse.getInputFeatures(muonTrack, mftTrack, muonTrack, mftTrackProp, mchTrackProp, collision); + } else { + const auto& globalMuonTrack = muonTracks.rawIteratorAt(globalTrackIndex); + inputML = mlResponse.getInputFeatures(globalMuonTrack, mftTrack, muonTrack, mftTrackProp, mchTrackProp, collision); + } mlResponse.isSelectedMl(inputML, 0, output); float matchScore = output[0]; - float matchChi2Prod = muonTrack.chi2MatchMCHMFT() / MatchingDegreesOfFreedom; - float matchScoreProd = chi2ToScore(muonTrack.chi2MatchMCHMFT(), MatchingDegreesOfFreedom, MatchingScoreChi2Max); - // std::cout << std::format("Matching score: {}, Chi2: {}", matchingScore, muonTrack.chi2MatchMCHMFT()) << std::endl; // check if a vector of global muon candidates is already available for the current MCH index // if not, initialize a new one and add the current global muon track @@ -2733,9 +2894,12 @@ struct QaMatching { if (matchingCandidateIterator != newMatchingCandidates.end()) { matchingCandidateIterator->second.emplace_back(MatchingCandidate{ muonTrack.collisionId(), - candidate.globalTrackId, + globalTrackIndex, mchIndex, mftTrack.globalIndex(), + globalTrackType, + mftTrackProp, + mchTrackProp, matchScore, -1, -1, @@ -2746,9 +2910,12 @@ struct QaMatching { } else { newMatchingCandidates[mchIndex].emplace_back(MatchingCandidate{ muonTrack.collisionId(), - candidate.globalTrackId, + globalTrackIndex, mchIndex, mftTrack.globalIndex(), + globalTrackType, + mftTrackProp, + mchTrackProp, matchScore, -1, -1, @@ -2772,23 +2939,25 @@ struct QaMatching { auto mftMchMatchAttempts = getMftMchMatchAttempts(collisions, bcs, mchTrack, mftTracks); int ranking = 1; for (auto& candidate : globalTracksVector) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) - const auto& muonTrack = muonTracks.rawIteratorAt(candidate.globalTrackId); - candidate.matchRanking = ranking; - candidate.matchType = getMatchType(muonTrack, muonTracks, mftTracks, matchablePairs, ranking); + if constexpr (isMC) { + candidate.matchType = getMatchType(candidate, muonTracks, mftTracks, matchablePairs, ranking); + } else { + candidate.matchType = kMatchTypeUndefined; + } candidate.mftMchMatchAttempts = mftMchMatchAttempts; ranking += 1; } } } - template - void processCollisionMc(const CollisionInfo& collisionInfo, - C const& collisions, - BC const& bcs, - TMUON const& muonTracks, - TMFT const& mftTracks, - CMFT const& mftCovs) + template + void processCollision(const CollisionInfo& collisionInfo, + C const& collisions, + BC const& bcs, + TMUON const& muonTracks, + TMFT const& mftTracks, + CMFT const& mftCovs) { auto collision = collisions.rawIteratorAt(collisionInfo.index); @@ -2806,8 +2975,12 @@ struct QaMatching { //------------------------------- // Chi2-based matching from production - fillQaMatchingAodTablesForCollision(collision, muonTracks, collisionInfo.matchingCandidates, matchingMethodCounter, collisionInfo.reducedEventId); - fillMatchingPlotsMc(collision, collisionInfo, muonTracks, mftTracks, collisionInfo.matchingCandidates, collisionInfo.matchingCandidates, collisionInfo.matchablePairs, cfgMatchingChi2ScoreMftMchLow, fChi2MatchingPlotter.get(), false); + fillQaMatchingAodTablesForCollision(collision, muonTracks, mftTracks, collisionInfo.matchingCandidates, matchingMethodCounter, collisionInfo.reducedEventId); + if constexpr (isMC) { + fillMatchingPlotsMc(collision, collisionInfo, muonTracks, mftTracks, collisionInfo.matchingCandidates, collisionInfo.matchingCandidates, collisionInfo.matchablePairs, cfgMatchingChi2ScoreMftMchLow, fChi2MatchingPlotter.get(), false); + } else { + fillMatchingPlots(collision, muonTracks, collisionInfo.matchingCandidates, fChi2MatchingPlotter.get()); + } //------------------------------- // Tagged muons @@ -2840,35 +3013,46 @@ struct QaMatching { taggedMatchingCandidates[mchIndex] = globalTracksVector; } } - matchingMethodCounter += 1; - fillQaMatchingAodTablesForCollision(collision, muonTracks, collisionInfo.matchingCandidates, matchingMethodCounter, collisionInfo.reducedEventId); - fillMatchingPlotsMc(collision, collisionInfo, muonTracks, mftTracks, taggedMatchingCandidates, collisionInfo.matchingCandidates, collisionInfo.matchablePairs, cfgMatchingChi2ScoreMftMchLow, fTaggedMuonsMatchingPlotter.get()); + if constexpr (isMC) { + fillMatchingPlotsMc(collision, collisionInfo, muonTracks, mftTracks, taggedMatchingCandidates, collisionInfo.matchingCandidates, collisionInfo.matchablePairs, cfgMatchingChi2ScoreMftMchLow, fTaggedMuonsMatchingPlotter.get()); + } else { + fillMatchingPlots(collision, muonTracks, taggedMatchingCandidates, fTaggedMuonsMatchingPlotter.get()); + } //------------------------------- // Custom chi2-based matching methods for (const auto& [label, func] : matchingChi2Functions) { MatchingCandidates matchingCandidates; - runChi2Matching(collisions, bcs, muonTracks, mftTracks, mftCovs, label, collisionInfo.matchablePairs, collisionInfo.matchingCandidates, matchingCandidates); + runChi2Matching(collisions, bcs, muonTracks, mftTracks, mftCovs, label, collisionInfo.matchablePairs, collisionInfo.matchingCandidates, matchingCandidates); auto* plotter = fMatchingPlotters.at(label).get(); double matchingScoreCut = matchingScoreCuts.at(label); matchingMethodCounter += 1; - fillQaMatchingAodTablesForCollision(collision, muonTracks, matchingCandidates, matchingMethodCounter, collisionInfo.reducedEventId); - fillMatchingPlotsMc(collision, collisionInfo, muonTracks, mftTracks, matchingCandidates, collisionInfo.matchingCandidates, collisionInfo.matchablePairs, matchingScoreCut, plotter, false); + fillQaMatchingAodTablesForCollision(collision, muonTracks, mftTracks, matchingCandidates, matchingMethodCounter, collisionInfo.reducedEventId); + if constexpr (isMC) { + fillMatchingPlotsMc(collision, collisionInfo, muonTracks, mftTracks, matchingCandidates, collisionInfo.matchingCandidates, collisionInfo.matchablePairs, matchingScoreCut, plotter, false); + } else { + fillMatchingPlots(collision, muonTracks, matchingCandidates, plotter); + } } - // ML-based matching analysis + //------------------------------- + // Custom ML-based matching methods for (const auto& [label, mlResponse] : matchingMlResponses) { MatchingCandidates matchingCandidates; - runMlMatching(collisions, bcs, muonTracks, mftTracks, mftCovs, label, collisionInfo.matchablePairs, collisionInfo.matchingCandidates, matchingCandidates); + runMlMatching(collisions, bcs, muonTracks, mftTracks, mftCovs, label, collisionInfo.matchablePairs, collisionInfo.matchingCandidates, matchingCandidates); auto* plotter = fMatchingPlotters.at(label).get(); double matchingScoreCut = matchingScoreCuts.at(label); matchingMethodCounter += 1; - fillQaMatchingAodTablesForCollision(collision, muonTracks, matchingCandidates, matchingMethodCounter, collisionInfo.reducedEventId); - fillMatchingPlotsMc(collision, collisionInfo, muonTracks, mftTracks, matchingCandidates, collisionInfo.matchingCandidates, collisionInfo.matchablePairs, matchingScoreCut, plotter); + fillQaMatchingAodTablesForCollision(collision, muonTracks, mftTracks, matchingCandidates, matchingMethodCounter, collisionInfo.reducedEventId); + if constexpr (isMC) { + fillMatchingPlotsMc(collision, collisionInfo, muonTracks, mftTracks, matchingCandidates, collisionInfo.matchingCandidates, collisionInfo.matchablePairs, matchingScoreCut, plotter); + } else { + fillMatchingPlots(collision, muonTracks, matchingCandidates, plotter); + } } //------------------------------- @@ -2876,9 +3060,10 @@ struct QaMatching { fillDimuonPlotsMc(collisionInfo, collisions, muonTracks, mftTracks); } - template + template void fillQaMatchingAodTablesForCollision(TCOLLISION const& collision, TMUON const& muonTracks, + TMFT const& mftTracks, const MatchingCandidates& matchingCandidates, int8_t matchLabel, int32_t reducedEventId) @@ -2894,17 +3079,19 @@ struct QaMatching { } for (const auto& candidate : candidates) { - const auto& candidateTrack = muonTracks.rawIteratorAt(candidate.globalTrackId); - auto candidateTrackAtVertex = VarManager::PropagateMuon(candidateTrack, collision, VarManager::kToVertex); + const auto& mftTrack = mftTracks.rawIteratorAt(candidate.mftTrackId); + // propagate global forward track to vertex using momentum rescaling method + auto candidateTrackAtVertex = propagateToVertexMft(mftTrack, mchTrack, collision); qaMatchingCandidates( reducedEventId, matchLabel, mchIndex, - static_cast(candidateTrack.p()), - static_cast(candidateTrack.pt()), - static_cast(candidateTrack.eta()), - static_cast(candidateTrack.phi()), + static_cast(mchTrack.p()), + static_cast(mchTrack.pt()), + static_cast(mchTrack.eta()), + static_cast(mchTrack.phi()), static_cast(candidate.matchType), + static_cast(candidate.matchChi2), static_cast(candidate.matchScore), static_cast(candidate.matchRanking), static_cast(candidateTrackAtVertex.getX()), @@ -3005,14 +3192,41 @@ struct QaMatching { mftTrackCovs[mftTrackCov.matchMFTTrackId()] = mftTrackCov.globalIndex(); } - fillCollisions(collisions, bcs, muonTracks, mftTracks, fCollisionInfos); + fillCollisions(collisions, bcs, muonTracks, mftTracks, mftCovs, fCollisionInfos); for (auto const& [collisionIndex, collisionInfo] : fCollisionInfos) { - processCollisionMc(collisionInfo, collisions, bcs, muonTracks, mftTracks, mftCovs); + processCollision(collisionInfo, collisions, bcs, muonTracks, mftTracks, mftCovs); } } PROCESS_SWITCH(QaMatching, processQAMC, "processQAMC", true); + + void processQA(MyEvents const& collisions, + aod::BCsWithTimestamps const& bcs, + MyMuons const& muonTracks, + MyMFTs const& mftTracks, + MyMFTCovariances const& mftCovs) + { + auto bc = bcs.begin(); + initCcdb(bc); + + for (const auto& muon : muonTracks) { + registry.get(HIST("nTracksPerType"))->Fill(static_cast(muon.trackType())); + } + + mftTrackCovs.clear(); + for (const auto& mftTrackCov : mftCovs) { + mftTrackCovs[mftTrackCov.matchMFTTrackId()] = mftTrackCov.globalIndex(); + } + + fillCollisions(collisions, bcs, muonTracks, mftTracks, mftCovs, fCollisionInfos); + + for (auto const& [collisionIndex, collisionInfo] : fCollisionInfos) { + processCollision(collisionInfo, collisions, bcs, muonTracks, mftTracks, mftCovs); + } + } + + PROCESS_SWITCH(QaMatching, processQA, "processQA", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)