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)