Skip to content

Commit ef14ae5

Browse files
aferrero2707aferrero2707
committed
[PWGDQ] add histograms of MFT-MCH match features in QA task
1 parent 5165955 commit ef14ae5

1 file changed

Lines changed: 135 additions & 14 deletions

File tree

‎PWGDQ/Tasks/qaMatching.cxx‎

Lines changed: 135 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -184,7 +184,7 @@ static float chi2ToScore(float chi2, int ndf, float chi2max)
184184
return static_cast<float>(result);
185185
}
186186

187-
static void SetMatchTypeAxisLabels(TAxis* axis)
187+
static void setMatchTypeAxisLabels(TAxis* axis)
188188
{
189189
axis->SetBinLabel(1, "true (leading)");
190190
axis->SetBinLabel(2, "wrong (leading)");
@@ -230,12 +230,15 @@ struct QaMatching {
230230
static constexpr int ExtrapolationMethodVertex = 3;
231231
static constexpr int ExtrapolationMethodMftDca = 4;
232232
static constexpr int DecayRankingDirect = 2;
233+
static constexpr float MatchingPlaneDefaultZ = -77.5;
233234

234235
struct MatchingCandidate {
235236
int64_t collisionId{-1};
236237
int64_t globalTrackId{-1};
237238
int64_t muonTrackId{-1};
238239
int64_t mftTrackId{-1};
240+
o2::track::TrackParCovFwd mftTrackProp;
241+
o2::track::TrackParCovFwd mchTrackProp;
239242
double matchScore{-1};
240243
double matchChi2{-1};
241244
int matchRanking{-1};
@@ -649,6 +652,42 @@ struct QaMatching {
649652
}
650653
};
651654

655+
struct MatchFeaturesHistos {
656+
o2::framework::HistPtr hDeltaP;
657+
o2::framework::HistPtr hDeltaPt;
658+
o2::framework::HistPtr hDeltaX;
659+
o2::framework::HistPtr hDeltaY;
660+
o2::framework::HistPtr hDeltaPhi;
661+
o2::framework::HistPtr hDeltaTanl;
662+
o2::framework::HistPtr hDeltaEta;
663+
o2::framework::HistPtr hRabs;
664+
665+
MatchFeaturesHistos(std::string path, HistogramRegistry* registry, int numCandidates)
666+
{
667+
AxisSpec indexAxis = {numCandidates + 1, 0, static_cast<double>(numCandidates + 1), "ranking index"};
668+
AxisSpec scoreAxis = {100, 0, 1, "match score"};
669+
int matchTypeMax = static_cast<int>(kMatchTypeUndefined) + 1;
670+
AxisSpec matchTypeAxis = {matchTypeMax, 0, static_cast<double>(matchTypeMax), "match type"};
671+
AxisSpec dxAxis = {100, -10, 10, "#Deltax (cm)"};
672+
AxisSpec dyAxis = {100, -10, 10, "#Deltay (cm)"};
673+
AxisSpec dpAxis = {100, -10, 10, "#Deltap (GeV/c)"};
674+
AxisSpec dptAxis = {100, -1, 1, "#Deltap_{T} (GeV/c)"};
675+
AxisSpec dphiAxis = {100, -1, 1, "#Delta#phi (rad)"};
676+
AxisSpec dtanlAxis = {100, -10, 10, "#Deltatanl"};
677+
AxisSpec detaAxis = {100, -1, 1, "#Delta#eta"};
678+
AxisSpec rabsAxis = {100, 0, 100, "R_{abs}"};
679+
680+
hDeltaP = registry->add((path + "/deltaP").c_str(), "MFT-MCH #Deltap", {HistType::kTHnSparseF, {dpAxis, scoreAxis, indexAxis, matchTypeAxis}});
681+
hDeltaPt = registry->add((path + "/deltaPt").c_str(), "MFT-MCH #Deltap_{T}", {HistType::kTHnSparseF, {dptAxis, scoreAxis, indexAxis, matchTypeAxis}});
682+
hDeltaX = registry->add((path + "/deltaX").c_str(), "MFT-MCH #Deltax", {HistType::kTHnSparseF, {dxAxis, scoreAxis, indexAxis, matchTypeAxis}});
683+
hDeltaY = registry->add((path + "/deltaY").c_str(), "MFT-MCH #Deltay", {HistType::kTHnSparseF, {dyAxis, scoreAxis, indexAxis, matchTypeAxis}});
684+
hDeltaPhi = registry->add((path + "/deltaPhi").c_str(), "MFT-MCH #Delta#phi", {HistType::kTHnSparseF, {dphiAxis, scoreAxis, indexAxis, matchTypeAxis}});
685+
hDeltaTanl = registry->add((path + "/deltaTanl").c_str(), "MFT-MCH #DeltaTanl", {HistType::kTHnSparseF, {dtanlAxis, scoreAxis, indexAxis, matchTypeAxis}});
686+
hDeltaEta = registry->add((path + "/deltaEta").c_str(), "MFT-MCH #Delta#eta", {HistType::kTHnSparseF, {detaAxis, scoreAxis, indexAxis, matchTypeAxis}});
687+
hRabs = registry->add((path + "/Rabs").c_str(), "MFT-MCH R_{abs}", {HistType::kTHnSparseF, {rabsAxis, scoreAxis, indexAxis, matchTypeAxis}});
688+
}
689+
};
690+
652691
struct MatchRankingHistos {
653692
o2::framework::HistPtr hist;
654693
o2::framework::HistPtr histVsP;
@@ -695,6 +734,8 @@ struct QaMatching {
695734
std::unique_ptr<MatchRankingHistos> fMatchRankingPaired;
696735
std::unique_ptr<MatchRankingHistos> fMatchRankingPairedGoodMCH;
697736

737+
std::unique_ptr<MatchFeaturesHistos> fMatchFeaturesGoodMCH;
738+
698739
//-
699740
o2::framework::HistPtr fMissedMatches;
700741
o2::framework::HistPtr fMissedMatchesGoodMCH;
@@ -757,6 +798,8 @@ struct QaMatching {
757798
std::string histName;
758799
std::string histTitle;
759800

801+
fMatchFeaturesGoodMCH = std::make_unique<MatchFeaturesHistos>(path + "matchFeaturesGoodMCH", registry, numCandidates);
802+
760803
if (isMc) {
761804
fMatchRanking = std::make_unique<MatchRankingHistos>(path + "matchRanking", "True match ranking", registry, mftMultMax, numCandidates);
762805
fMatchRankingGoodMCH = std::make_unique<MatchRankingHistos>(path + "matchRankingGoodMCH", "True match ranking (good MCH tracks)", registry, mftMultMax, numCandidates);
@@ -808,41 +851,41 @@ struct QaMatching {
808851
histName = path + "matchType";
809852
histTitle = "Match type";
810853
fMatchType = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {matchTypeAxis}});
811-
SetMatchTypeAxisLabels(std::get<std::shared_ptr<TH1>>(fMatchType)->GetXaxis());
854+
setMatchTypeAxisLabels(std::get<std::shared_ptr<TH1>>(fMatchType)->GetXaxis());
812855
histName = path + "matchTypeVsP";
813856
histTitle = "Match type vs. p";
814857
fMatchTypeVsP = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {pAxis, matchTypeAxis}});
815-
SetMatchTypeAxisLabels(std::get<std::shared_ptr<TH2>>(fMatchTypeVsP)->GetYaxis());
858+
setMatchTypeAxisLabels(std::get<std::shared_ptr<TH2>>(fMatchTypeVsP)->GetYaxis());
816859
histName = path + "matchTypeVsPt";
817860
histTitle = "Match type vs. p_{T}";
818861
fMatchTypeVsPt = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {ptAxis, matchTypeAxis}});
819-
SetMatchTypeAxisLabels(std::get<std::shared_ptr<TH2>>(fMatchTypeVsPt)->GetYaxis());
862+
setMatchTypeAxisLabels(std::get<std::shared_ptr<TH2>>(fMatchTypeVsPt)->GetYaxis());
820863

821864
histName = path + "matchChi2VsType";
822865
histTitle = "Match #chi^{2} vs. match type";
823866
fMatchChi2VsType = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {matchTypeAxis, chi2Axis}});
824-
SetMatchTypeAxisLabels(std::get<std::shared_ptr<TH2>>(fMatchChi2VsType)->GetXaxis());
867+
setMatchTypeAxisLabels(std::get<std::shared_ptr<TH2>>(fMatchChi2VsType)->GetXaxis());
825868
histName = path + "matchChi2VsTypeVsP";
826869
histTitle = "Match #chi^{2} vs. match type vs. p";
827870
fMatchChi2VsTypeVsP = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH3F, {pAxis, matchTypeAxis, chi2Axis}});
828-
SetMatchTypeAxisLabels(std::get<std::shared_ptr<TH3>>(fMatchChi2VsTypeVsP)->GetYaxis());
871+
setMatchTypeAxisLabels(std::get<std::shared_ptr<TH3>>(fMatchChi2VsTypeVsP)->GetYaxis());
829872
histName = path + "matchChi2VsTypeVsPt";
830873
histTitle = "Match #chi^{2} vs. match type vs. p_{T}";
831874
fMatchChi2VsTypeVsPt = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH3F, {ptAxis, matchTypeAxis, chi2Axis}});
832-
SetMatchTypeAxisLabels(std::get<std::shared_ptr<TH3>>(fMatchChi2VsTypeVsPt)->GetYaxis());
875+
setMatchTypeAxisLabels(std::get<std::shared_ptr<TH3>>(fMatchChi2VsTypeVsPt)->GetYaxis());
833876
//-
834877
histName = path + "matchScoreVsType";
835878
histTitle = "Match score vs. match type";
836879
fMatchScoreVsType = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {matchTypeAxis, scoreAxis}});
837-
SetMatchTypeAxisLabels(std::get<std::shared_ptr<TH2>>(fMatchScoreVsType)->GetXaxis());
880+
setMatchTypeAxisLabels(std::get<std::shared_ptr<TH2>>(fMatchScoreVsType)->GetXaxis());
838881
histName = path + "matchScoreVsTypeVsP";
839882
histTitle = "Match score vs. match type vs. p";
840883
fMatchScoreVsTypeVsP = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH3F, {pAxis, matchTypeAxis, scoreAxis}});
841-
SetMatchTypeAxisLabels(std::get<std::shared_ptr<TH3>>(fMatchScoreVsTypeVsP)->GetYaxis());
884+
setMatchTypeAxisLabels(std::get<std::shared_ptr<TH3>>(fMatchScoreVsTypeVsP)->GetYaxis());
842885
histName = path + "matchScoreVsTypeVsPt";
843886
histTitle = "Match score vs. match type vs. p_{T}";
844887
fMatchScoreVsTypeVsPt = registry->add(histName.c_str(), histTitle.c_str(), {HistType::kTH3F, {ptAxis, matchTypeAxis, scoreAxis}});
845-
SetMatchTypeAxisLabels(std::get<std::shared_ptr<TH3>>(fMatchScoreVsTypeVsPt)->GetYaxis());
888+
setMatchTypeAxisLabels(std::get<std::shared_ptr<TH3>>(fMatchScoreVsTypeVsPt)->GetYaxis());
846889
}
847890

848891
AxisSpec prodScoreAxis = {100, 0, 1, "matching score (prod)"};
@@ -909,8 +952,6 @@ struct QaMatching {
909952
AxisSpec phiAxis = {90, -180, 180, "#phi (degrees)"};
910953
std::string histPath = cfgIsMc.value ? "matching/MC/" : "matching/";
911954

912-
std::cout << std::format("TOTO cfgIsMc: {}", cfgIsMc.value) << std::endl;
913-
914955
AxisSpec trackPositionXAtMftAxis = {100, -15, 15, "MFT x (cm)"};
915956
AxisSpec trackPositionYAtMftAxis = {100, -15, 15, "MFT y (cm)"};
916957
registry.add((histPath + "pairedMCHTracksAtMFT").c_str(), "Paired MCH tracks position at MFT end", {HistType::kTH2F, {trackPositionXAtMftAxis, trackPositionYAtMftAxis}});
@@ -1931,6 +1972,7 @@ struct QaMatching {
19311972
BC const& bcs,
19321973
TMUON const& muonTracks,
19331974
TMFT const& mftTracks,
1975+
MyMFTCovariances const& mftCovs,
19341976
CollisionInfos& collisionInfos)
19351977
{
19361978
collisionInfos.clear();
@@ -1982,6 +2024,18 @@ struct QaMatching {
19822024
auto const& mftTrack = muonTrack.template matchMFTTrack_as<TMFT>();
19832025
int64_t mftTrackIndex = mftTrack.globalIndex();
19842026

2027+
// get MFT track covariances
2028+
if (mftTrackCovs.count(mftTrack.globalIndex()) < 1) {
2029+
continue;
2030+
}
2031+
auto const& mftTrackCov = mftCovs.rawIteratorAt(mftTrackCovs[mftTrack.globalIndex()]);
2032+
2033+
// propagate MCH and MFT tracks to matching plane
2034+
auto mchTrackProp = fwdToTrackPar(mchTrack, mchTrack);
2035+
mchTrackProp = propagateToMatchingPlaneMch(mchTrack, mftTrack, mftTrackCov, collision, MatchingPlaneDefaultZ, 0);
2036+
auto mftTrackProp = fwdToTrackPar(mftTrack, mftTrackCov);
2037+
mftTrackProp = propagateToMatchingPlaneMft(mchTrack, mftTrack, mftTrackCov, collision, MatchingPlaneDefaultZ, 0);
2038+
19852039
// check if a vector of global muon candidates is already available for the current MCH index
19862040
// if not, initialize a new one and add the current global muon track
19872041
// bool globalMuonTrackFound = false;
@@ -1992,6 +2046,8 @@ struct QaMatching {
19922046
muonTrackIndex,
19932047
mchTrackIndex,
19942048
mftTrackIndex,
2049+
mftTrackProp,
2050+
mchTrackProp,
19952051
matchScore,
19962052
matchChi2,
19972053
-1,
@@ -2006,6 +2062,8 @@ struct QaMatching {
20062062
muonTrackIndex,
20072063
mchTrackIndex,
20082064
mftTrackIndex,
2065+
mftTrackProp,
2066+
mchTrackProp,
20092067
matchScore,
20102068
matchChi2,
20112069
-1,
@@ -2077,6 +2135,39 @@ struct QaMatching {
20772135
const MatchingCandidates& matchingCandidates,
20782136
MatchingPlotter* plotter)
20792137
{
2138+
// ====================================
2139+
// Matching properties
2140+
for (const auto& [mchIndex, globalTracksVector] : matchingCandidates) {
2141+
if (globalTracksVector.size() < 1)
2142+
continue;
2143+
2144+
// get the standalone MCH track
2145+
auto const& mchTrack = muonTracks.rawIteratorAt(mchIndex);
2146+
2147+
// skip global muon tracks that do not pass the MCH and MFT quality cuts
2148+
if (!isGoodGlobalMuon(mchTrack, collision))
2149+
continue;
2150+
2151+
for (const auto& candidate : globalTracksVector) {
2152+
double dp = candidate.mchTrackProp.getP() - candidate.mftTrackProp.getP();
2153+
double dpt = candidate.mchTrackProp.getPt() - candidate.mftTrackProp.getPt();
2154+
double dx = candidate.mchTrackProp.getX() - candidate.mftTrackProp.getX();
2155+
double dy = candidate.mchTrackProp.getY() - candidate.mftTrackProp.getY();
2156+
double dphi = candidate.mchTrackProp.getPhi() - candidate.mftTrackProp.getPhi();
2157+
double dtanl = candidate.mchTrackProp.getTanl() - candidate.mftTrackProp.getTanl();
2158+
double deta = candidate.mchTrackProp.getEta() - candidate.mftTrackProp.getEta();
2159+
int matchType = static_cast<int>(candidate.matchType);
2160+
std::get<std::shared_ptr<THnSparse>>(plotter->fMatchFeaturesGoodMCH->hDeltaP)->Fill(dp, candidate.matchScore, candidate.matchRanking, matchType);
2161+
std::get<std::shared_ptr<THnSparse>>(plotter->fMatchFeaturesGoodMCH->hDeltaPt)->Fill(dpt, candidate.matchScore, candidate.matchRanking, matchType);
2162+
std::get<std::shared_ptr<THnSparse>>(plotter->fMatchFeaturesGoodMCH->hDeltaX)->Fill(dx, candidate.matchScore, candidate.matchRanking, matchType);
2163+
std::get<std::shared_ptr<THnSparse>>(plotter->fMatchFeaturesGoodMCH->hDeltaY)->Fill(dy, candidate.matchScore, candidate.matchRanking, matchType);
2164+
std::get<std::shared_ptr<THnSparse>>(plotter->fMatchFeaturesGoodMCH->hDeltaPhi)->Fill(dphi, candidate.matchScore, candidate.matchRanking, matchType);
2165+
std::get<std::shared_ptr<THnSparse>>(plotter->fMatchFeaturesGoodMCH->hDeltaTanl)->Fill(dtanl, candidate.matchScore, candidate.matchRanking, matchType);
2166+
std::get<std::shared_ptr<THnSparse>>(plotter->fMatchFeaturesGoodMCH->hDeltaEta)->Fill(deta, candidate.matchScore, candidate.matchRanking, matchType);
2167+
std::get<std::shared_ptr<THnSparse>>(plotter->fMatchFeaturesGoodMCH->hRabs)->Fill(mchTrack.rAtAbsorberEnd(), candidate.matchScore, candidate.matchRanking, matchType);
2168+
}
2169+
}
2170+
20802171
// ====================================
20812172
// Matching chi2 and score vs. production values
20822173
for (const auto& [mchIndex, globalTracksVector] : matchingCandidates) {
@@ -2206,6 +2297,28 @@ struct QaMatching {
22062297
std::get<std::shared_ptr<TH2>>(plotter->fMatchRankingPairedGoodMCH->histVsDeltaChi2)->Fill(dchi2, trueMatchIndex);
22072298
}
22082299

2300+
if (isGoodMCH) {
2301+
for (const auto& candidate : globalTracksVector) {
2302+
double dp = candidate.mchTrackProp.getP() - candidate.mftTrackProp.getP();
2303+
double dpt = candidate.mchTrackProp.getPt() - candidate.mftTrackProp.getPt();
2304+
double dx = candidate.mchTrackProp.getX() - candidate.mftTrackProp.getX();
2305+
double dy = candidate.mchTrackProp.getY() - candidate.mftTrackProp.getY();
2306+
double dphi = candidate.mchTrackProp.getPhi() - candidate.mftTrackProp.getPhi();
2307+
double dtanl = candidate.mchTrackProp.getTanl() - candidate.mftTrackProp.getTanl();
2308+
double deta = candidate.mchTrackProp.getEta() - candidate.mftTrackProp.getEta();
2309+
int matchType = static_cast<int>(candidate.matchType);
2310+
std::get<std::shared_ptr<THnSparse>>(plotter->fMatchFeaturesGoodMCH->hDeltaP)->Fill(dp, candidate.matchScore, candidate.matchRanking, matchType);
2311+
std::get<std::shared_ptr<THnSparse>>(plotter->fMatchFeaturesGoodMCH->hDeltaPt)->Fill(dpt, candidate.matchScore, candidate.matchRanking, matchType);
2312+
std::get<std::shared_ptr<THnSparse>>(plotter->fMatchFeaturesGoodMCH->hDeltaX)->Fill(dx, candidate.matchScore, candidate.matchRanking, matchType);
2313+
std::get<std::shared_ptr<THnSparse>>(plotter->fMatchFeaturesGoodMCH->hDeltaY)->Fill(dy, candidate.matchScore, candidate.matchRanking, matchType);
2314+
std::get<std::shared_ptr<THnSparse>>(plotter->fMatchFeaturesGoodMCH->hDeltaPhi)->Fill(dphi, candidate.matchScore, candidate.matchRanking, matchType);
2315+
std::get<std::shared_ptr<THnSparse>>(plotter->fMatchFeaturesGoodMCH->hDeltaTanl)->Fill(dtanl, candidate.matchScore, candidate.matchRanking, matchType);
2316+
std::get<std::shared_ptr<THnSparse>>(plotter->fMatchFeaturesGoodMCH->hDeltaEta)->Fill(deta, candidate.matchScore, candidate.matchRanking, matchType);
2317+
std::get<std::shared_ptr<THnSparse>>(plotter->fMatchFeaturesGoodMCH->hRabs)->Fill(deta, candidate.matchScore, candidate.matchRanking, matchType);
2318+
std::get<std::shared_ptr<THnSparse>>(plotter->fMatchFeaturesGoodMCH->hRabs)->Fill(mchTrack.rAtAbsorberEnd(), candidate.matchScore, candidate.matchRanking, matchType);
2319+
}
2320+
}
2321+
22092322
if (trueMatchIndex == 0) {
22102323
// missed matches
22112324
if (!isPairedMCH) {
@@ -2577,6 +2690,8 @@ struct QaMatching {
25772690
candidate.globalTrackId,
25782691
mchIndex,
25792692
mftTrack.globalIndex(),
2693+
mftTrackProp,
2694+
mchTrackProp,
25802695
matchScore,
25812696
matchChi2,
25822697
-1,
@@ -2590,6 +2705,8 @@ struct QaMatching {
25902705
candidate.globalTrackId,
25912706
mchIndex,
25922707
mftTrack.globalIndex(),
2708+
mftTrackProp,
2709+
mchTrackProp,
25932710
matchScore,
25942711
matchChi2,
25952712
-1,
@@ -2722,6 +2839,8 @@ struct QaMatching {
27222839
candidate.globalTrackId,
27232840
mchIndex,
27242841
mftTrack.globalIndex(),
2842+
mftTrackProp,
2843+
mchTrackProp,
27252844
matchScore,
27262845
-1,
27272846
-1,
@@ -2735,6 +2854,8 @@ struct QaMatching {
27352854
candidate.globalTrackId,
27362855
mchIndex,
27372856
mftTrack.globalIndex(),
2857+
mftTrackProp,
2858+
mchTrackProp,
27382859
matchScore,
27392860
-1,
27402861
-1,
@@ -3010,7 +3131,7 @@ struct QaMatching {
30103131
mftTrackCovs[mftTrackCov.matchMFTTrackId()] = mftTrackCov.globalIndex();
30113132
}
30123133

3013-
fillCollisions<true>(collisions, bcs, muonTracks, mftTracks, fCollisionInfos);
3134+
fillCollisions<true>(collisions, bcs, muonTracks, mftTracks, mftCovs, fCollisionInfos);
30143135

30153136
for (auto const& [collisionIndex, collisionInfo] : fCollisionInfos) {
30163137
processCollision<true>(collisionInfo, collisions, bcs, muonTracks, mftTracks, mftCovs);
@@ -3037,7 +3158,7 @@ struct QaMatching {
30373158
mftTrackCovs[mftTrackCov.matchMFTTrackId()] = mftTrackCov.globalIndex();
30383159
}
30393160

3040-
fillCollisions<false>(collisions, bcs, muonTracks, mftTracks, fCollisionInfos);
3161+
fillCollisions<false>(collisions, bcs, muonTracks, mftTracks, mftCovs, fCollisionInfos);
30413162

30423163
for (auto const& [collisionIndex, collisionInfo] : fCollisionInfos) {
30433164
processCollision<false>(collisionInfo, collisions, bcs, muonTracks, mftTracks, mftCovs);

0 commit comments

Comments
 (0)