Skip to content

[PWGDQ] added track slope dependence to residals plots #15999

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
mcoquet642 merged 1 commit into from
Apr 28, 2026
+37 −30
Merged

[PWGDQ] added track slope dependence to residals plots #15999

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
67 changes: 37 additions & 30 deletions PWGDQ/Tasks/muonGlobalAlignment.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -272,7 +272,7 @@
std::map<uint64_t, CollisionInfo>& collisionInfos)
{
// fill collision information for global muon tracks (MFT-MCH-MID matches)
for (auto muonTrack : muonTracks) {

Check failure on line 275 in PWGDQ/Tasks/muonGlobalAlignment.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.
if (!muonTrack.has_collision())
continue;

Expand Down Expand Up @@ -320,8 +320,8 @@
return (track1.chi2MatchMCHMFT() < track2.chi2MatchMCHMFT());
};

for (auto& [collisionIndex, collisionInfo] : collisionInfos) {

Check failure on line 323 in PWGDQ/Tasks/muonGlobalAlignment.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.
for (auto& [mchIndex, globalTracksVector] : collisionInfo.globalMuonTracks) {

Check failure on line 324 in PWGDQ/Tasks/muonGlobalAlignment.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.
std::sort(globalTracksVector.begin(), globalTracksVector.end(), compareChi2);
}
}
Expand All @@ -336,7 +336,7 @@
InitCollisions(collisions, bcs, muonTracks, collisionInfos);

// fill collision information for MFT standalone tracks
for (auto mftTrack : mftTracks) {

Check failure on line 339 in PWGDQ/Tasks/muonGlobalAlignment.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[const-ref-in-for-loop] 

Use constant references for non-modified iterators in range-based for loops.
if (!mftTrack.has_collision())
continue;

Expand Down Expand Up @@ -528,8 +528,8 @@
}

if (fEnableMftMchResidualsAnalysis) {
AxisSpec dxAxis = {600, -30.0, 30.0, "#Delta x (cm)"};
AxisSpec dyAxis = {600, -30.0, 30.0, "#Delta y (cm)"};
AxisSpec dxAxis = {400, -20.0, 20.0, "#Delta x (cm)"};
AxisSpec dyAxis = {400, -20.0, 20.0, "#Delta y (cm)"};

registry.add("DCA/MCH/DCA_y_vs_x", std::format("DCA y vs. x").c_str(), {HistType::kTH2F, {dcaxMCHAxis, dcayMCHAxis}});
registry.add("DCA/MCH/DCA_x_vs_sign_vs_quadrant_vs_mom", std::format("DCA(x) vs. p, quadrant, chargeSign").c_str(), {HistType::kTHnSparseF, {{20, 0, 100.0, "p (GeV/c)"}, {4, 0, 4, "quadrant"}, {2, 0, 2, "sign"}, dcaxMCHAxis}});
Expand All @@ -543,20 +543,20 @@
registry.add("residuals/dy_vs_chamber", "Cluster y residual vs. chamber, quadrant, chargeSign",
{HistType::kTHnSparseF, {{10, 1, 11, "chamber"}, {4, 0, 4, "quadrant"}, {2, 0, 2, "sign"}, dyAxis}});

registry.add("residuals/dx_vs_de", "Cluster x residual vs. DE, quadrant, chargeSign, momentum",
{HistType::kTHnSparseF, {dxAxis, {getNumDE(), 0, static_cast<double>(getNumDE()), "DE"}, {4, 0, 4, "quadrant"}, {2, 0, 2, "sign"}, {20, 0, 100.0, "p (GeV/c)"}}});
registry.add("residuals/dy_vs_de", "Cluster y residual vs. DE, quadrant, chargeSign, momentum",
{HistType::kTHnSparseF, {dyAxis, {getNumDE(), 0, static_cast<double>(getNumDE()), "DE"}, {4, 0, 4, "quadrant"}, {2, 0, 2, "sign"}, {20, 0, 100.0, "p (GeV/c)"}}});
registry.add("residuals/dx_vs_de", "Cluster x residual vs. DE, quadrant, chargeSign, momentum, xslope",
{HistType::kTHnSparseF, {dxAxis, {getNumDE(), 0, static_cast<double>(getNumDE()), "DE"}, {4, 0, 4, "quadrant"}, {2, 0, 2, "sign"}, {20, 0, 100.0, "p (GeV/c)"}, sxAxis}});
registry.add("residuals/dy_vs_de", "Cluster y residual vs. DE, quadrant, chargeSign, momentum, yslope",
{HistType::kTHnSparseF, {dyAxis, {getNumDE(), 0, static_cast<double>(getNumDE()), "DE"}, {4, 0, 4, "quadrant"}, {2, 0, 2, "sign"}, {20, 0, 100.0, "p (GeV/c)"}, syAxis}});

registry.add("residuals/dx_vs_chamber_corr", "Cluster x residual vs. chamber, quadrant, chargeSign (with corrections)",
{HistType::kTHnSparseF, {{10, 1, 11, "chamber"}, {4, 0, 4, "quadrant"}, {2, 0, 2, "sign"}, dxAxis}});
registry.add("residuals/dy_vs_chamber_corr", "Cluster y residual vs. chamber, quadrant, chargeSign (with corrections)",
{HistType::kTHnSparseF, {{10, 1, 11, "chamber"}, {4, 0, 4, "quadrant"}, {2, 0, 2, "sign"}, dyAxis}});

registry.add("residuals/dx_vs_de_corr", "Cluster x residual vs. DE, quadrant, chargeSign, momentum (with corrections)",
{HistType::kTHnSparseF, {dxAxis, {getNumDE(), 0, static_cast<double>(getNumDE()), "DE"}, {4, 0, 4, "quadrant"}, {2, 0, 2, "sign"}, {20, 0, 100.0, "p (GeV/c)"}}});
registry.add("residuals/dy_vs_de_corr", "Cluster y residual vs. DE, quadrant, chargeSign, momentum (with corrections)",
{HistType::kTHnSparseF, {dyAxis, {getNumDE(), 0, static_cast<double>(getNumDE()), "DE"}, {4, 0, 4, "quadrant"}, {2, 0, 2, "sign"}, {20, 0, 100.0, "p (GeV/c)"}}});
registry.add("residuals/dx_vs_de_corr", "Cluster x residual vs. DE, quadrant, chargeSign, momentum, xslope (with corrections)",
{HistType::kTHnSparseF, {dxAxis, {getNumDE(), 0, static_cast<double>(getNumDE()), "DE"}, {4, 0, 4, "quadrant"}, {2, 0, 2, "sign"}, {20, 0, 100.0, "p (GeV/c)"}, sxAxis}});
registry.add("residuals/dy_vs_de_corr", "Cluster y residual vs. DE, quadrant, chargeSign, momentum, yslope (with corrections)",
{HistType::kTHnSparseF, {dyAxis, {getNumDE(), 0, static_cast<double>(getNumDE()), "DE"}, {4, 0, 4, "quadrant"}, {2, 0, 2, "sign"}, {20, 0, 100.0, "p (GeV/c)"}, syAxis}});

registry.add("residuals/de_alignment_corrections_x", "DE alignment corrections - X coordinate",
{HistType::kTH1F, {{getNumDE(), 0, static_cast<double>(getNumDE()), "DE"}}});
Expand Down Expand Up @@ -1094,8 +1094,8 @@
track.setBendingCoor(y + yCorrection);
track.setBendingSlope(ySlope + ySlopeCorrection);
/*
std::cout << std::format("[TOTO] MFT position: pos={:0.3f},{:0.3f}", x, y) << std::endl;

Check failure on line 1097 in PWGDQ/Tasks/muonGlobalAlignment.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[logging] 

Use O2 logging (LOG, LOGF, LOGP).
std::cout << std::format("[TOTO] MFT corrections: pos={:0.3f},{:0.3f} slope={:0.12f},{:0.12f} angle={:0.12f},{:0.12f}",

Check failure on line 1098 in PWGDQ/Tasks/muonGlobalAlignment.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[logging] 

Use O2 logging (LOG, LOGF, LOGP).
xCorrection, yCorrection, xSlopeCorrection, ySlopeCorrection,
std::atan2(xSlopeCorrection, 1), std::atan2(ySlopeCorrection, 1)) << std::endl;
*/
Expand Down Expand Up @@ -1133,8 +1133,8 @@
template <typename T>
T UpdateTrackMomentum(const T& track, const double p, int sign)
{
double px = p * sin(M_PI / 2 - atan(track.tgl())) * cos(track.phi());

Check failure on line 1136 in PWGDQ/Tasks/muonGlobalAlignment.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[std-prefix] 

Use std:: prefix for names from the std namespace.
double py = p * sin(M_PI / 2 - atan(track.tgl())) * sin(track.phi());

Check failure on line 1137 in PWGDQ/Tasks/muonGlobalAlignment.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[std-prefix] 

Use std:: prefix for names from the std namespace.
double pt = std::sqrt(std::pow(px, 2) + std::pow(py, 2));

SMatrix5 tpars = {track.x(), track.y(), track.phi(), track.tgl(), sign / pt};
Expand All @@ -1154,8 +1154,8 @@
template <typename T>
T UpdateTrackMomentum(const T& track, const o2::mch::TrackParam& track4mom)
{
double px = track4mom.p() * sin(M_PI / 2 - atan(track.tgl())) * cos(track.phi());

Check failure on line 1157 in PWGDQ/Tasks/muonGlobalAlignment.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[std-prefix] 

Use std:: prefix for names from the std namespace.
double py = track4mom.p() * sin(M_PI / 2 - atan(track.tgl())) * sin(track.phi());

Check failure on line 1158 in PWGDQ/Tasks/muonGlobalAlignment.cxx

View workflow job for this annotation

GitHub Actions / O2 linter

[std-prefix] 

Use std:: prefix for names from the std namespace.
double pt = std::sqrt(std::pow(px, 2) + std::pow(py, 2));
double sign = track4mom.getCharge();

Expand Down Expand Up @@ -1656,7 +1656,8 @@
auto const& muonTrack = muonTracks.rawIteratorAt(globalTracksVector[0]);
const auto& mchTrack = muonTrack.template matchMCHTrack_as<MyMuonsWithCov>();
const auto& mftTrack = muonTrack.template matchMFTTrack_as<MyMFTs>();
int quadrantMch = GetQuadrant(mchTrack);
// int quadrant = GetQuadrant(mchTrack);
int quadrant = GetQuadrant(mftTrack);
int posNeg = (mchTrack.sign() >= 0) ? 0 : 1;

bool isGoodMuon = IsGoodMuon(mchTrack, collision, fTrackChi2MchUp, fMftMchResidualsPLow, fMftMchResidualsPtLow, {fEtaMftLow, fEtaMftUp}, {fRabsLow, fRabsUp}, fSigmaPdcaUp);
Expand All @@ -1667,6 +1668,10 @@
if (!isGoodMFT)
continue;

double matchChi2 = muonTrack.chi2MatchMCHMFT() / 5.f;
if (matchChi2 > 10.f)
continue;

// refit MCH track if enabled
TrackRealigned convertedTrack;
bool convertedTrackOk = false;
Expand Down Expand Up @@ -1722,15 +1727,16 @@
// by taking the momentum from the MCH track refitted with the new alignment
if (!configRealign.fEnableMCHRealign || convertedTrackOk) {
auto mftTrackAtCluster = configRealign.fEnableMCHRealign ? PropagateMFTtoMCH(mftTrack, mch::TrackParam(convertedTrack.first()), master.z()) : PropagateMFTtoMCH(mftTrack, FwdtoMCH(FwdToTrackPar(mchTrack)), master.z());
auto mftTrackParamAtCluster = FwdtoMCH(mftTrackAtCluster);

std::array<double, 2> xPos{master.x(), mftTrackAtCluster.getX()};
std::array<double, 2> yPos{master.y(), mftTrackAtCluster.getY()};

registry.get<THnSparse>(HIST("residuals/dx_vs_chamber"))->Fill(chamber + 1, quadrantMch, posNeg, xPos[0] - xPos[1]);
registry.get<THnSparse>(HIST("residuals/dy_vs_chamber"))->Fill(chamber + 1, quadrantMch, posNeg, yPos[0] - yPos[1]);
registry.get<THnSparse>(HIST("residuals/dx_vs_chamber"))->Fill(chamber + 1, quadrant, posNeg, xPos[0] - xPos[1]);
registry.get<THnSparse>(HIST("residuals/dy_vs_chamber"))->Fill(chamber + 1, quadrant, posNeg, yPos[0] - yPos[1]);

registry.get<THnSparse>(HIST("residuals/dx_vs_de"))->Fill(xPos[0] - xPos[1], deIndex, quadrantMch, posNeg, mchTrack.p());
registry.get<THnSparse>(HIST("residuals/dy_vs_de"))->Fill(yPos[0] - yPos[1], deIndex, quadrantMch, posNeg, mchTrack.p());
registry.get<THnSparse>(HIST("residuals/dx_vs_de"))->Fill(xPos[0] - xPos[1], deIndex, quadrant, posNeg, mchTrack.p(), mftTrackParamAtCluster.getNonBendingSlope());
registry.get<THnSparse>(HIST("residuals/dy_vs_de"))->Fill(yPos[0] - yPos[1], deIndex, quadrant, posNeg, mchTrack.p(), mftTrackParamAtCluster.getBendingSlope());
}

// MFT-MCH residuals with realigned and/or corrected MCH clusters
Expand All @@ -1739,15 +1745,16 @@
// alignment (if realignment is enabled)
if (convertedTrackWithCorrOk) {
auto mftTrackAtClusterWithCorr = PropagateMFTtoMCH(mftTrack, mch::TrackParam(convertedTrackWithCorr.first()), masterWithCorr.z());
auto mftTrackParamAtClusterWithCorr = FwdtoMCH(mftTrackAtClusterWithCorr);

std::array<double, 2> xPos{masterWithCorr.x(), mftTrackAtClusterWithCorr.getX()};
std::array<double, 2> yPos{masterWithCorr.y(), mftTrackAtClusterWithCorr.getY()};

registry.get<THnSparse>(HIST("residuals/dx_vs_chamber_corr"))->Fill(chamber + 1, quadrantMch, posNeg, xPos[0] - xPos[1]);
registry.get<THnSparse>(HIST("residuals/dy_vs_chamber_corr"))->Fill(chamber + 1, quadrantMch, posNeg, yPos[0] - yPos[1]);
registry.get<THnSparse>(HIST("residuals/dx_vs_chamber_corr"))->Fill(chamber + 1, quadrant, posNeg, xPos[0] - xPos[1]);
registry.get<THnSparse>(HIST("residuals/dy_vs_chamber_corr"))->Fill(chamber + 1, quadrant, posNeg, yPos[0] - yPos[1]);

registry.get<THnSparse>(HIST("residuals/dx_vs_de_corr"))->Fill(xPos[0] - xPos[1], deIndex, quadrantMch, posNeg, mchTrack.p());
registry.get<THnSparse>(HIST("residuals/dy_vs_de_corr"))->Fill(yPos[0] - yPos[1], deIndex, quadrantMch, posNeg, mchTrack.p());
registry.get<THnSparse>(HIST("residuals/dx_vs_de_corr"))->Fill(xPos[0] - xPos[1], deIndex, quadrant, posNeg, mchTrack.p(), mftTrackParamAtClusterWithCorr.getNonBendingSlope());
registry.get<THnSparse>(HIST("residuals/dy_vs_de_corr"))->Fill(yPos[0] - yPos[1], deIndex, quadrant, posNeg, mchTrack.p(), mftTrackParamAtClusterWithCorr.getBendingSlope());
}
}

Expand All @@ -1757,12 +1764,12 @@
auto dcay = mchTrackAtDCA.getY() - collision.posY();

registry.get<TH2>(HIST("DCA/MCH/DCA_y_vs_x"))->Fill(dcax, dcay);
registry.get<THnSparse>(HIST("DCA/MCH/DCA_x_vs_sign_vs_quadrant_vs_mom"))->Fill(mchTrack.p(), quadrantMch, posNeg, dcax);
registry.get<THnSparse>(HIST("DCA/MCH/DCA_y_vs_sign_vs_quadrant_vs_mom"))->Fill(mchTrack.p(), quadrantMch, posNeg, dcay);
registry.get<THnSparse>(HIST("DCA/MCH/DCA_x_vs_sign_vs_quadrant_vs_mom"))->Fill(mchTrack.p(), quadrant, posNeg, dcax);
registry.get<THnSparse>(HIST("DCA/MCH/DCA_y_vs_sign_vs_quadrant_vs_mom"))->Fill(mchTrack.p(), quadrant, posNeg, dcay);

if (fEnableMftMchResidualsExtraPlots) {
registry.get<THnSparse>(HIST("DCA/MCH/DCA_x_vs_sign_vs_quadrant_vs_vz"))->Fill(collision.posZ(), quadrantMch, posNeg, dcax);
registry.get<THnSparse>(HIST("DCA/MCH/DCA_y_vs_sign_vs_quadrant_vs_vz"))->Fill(collision.posZ(), quadrantMch, posNeg, dcay);
registry.get<THnSparse>(HIST("DCA/MCH/DCA_x_vs_sign_vs_quadrant_vs_vz"))->Fill(collision.posZ(), quadrant, posNeg, dcax);
registry.get<THnSparse>(HIST("DCA/MCH/DCA_y_vs_sign_vs_quadrant_vs_vz"))->Fill(collision.posZ(), quadrant, posNeg, dcay);
auto mchTrackAtMFT = configRealign.fEnableMCHRealign ? PropagateMCHRealigned(convertedTrack, mftTrack.z()) : PropagateMCH(mchTrack, mftTrack.z());
double deltaPhi = mchTrackAtMFT.getPhi() - mftTrack.phi();
registry.get<THnSparse>(HIST("residuals/dphi_at_mft"))->Fill(deltaPhi, mftTrack.x(), mftTrack.y(), posNeg, mchTrackAtMFT.getP());
Expand All @@ -1774,8 +1781,8 @@
auto dcax = mchTrackAtDCA.getX() - collision.posX();
auto dcay = mchTrackAtDCA.getY() - collision.posY();

registry.get<THnSparse>(HIST("DCA/MCH/DCA_x_vs_sign_vs_quadrant_vs_mom_corr"))->Fill(mchTrack.p(), quadrantMch, posNeg, dcax);
registry.get<THnSparse>(HIST("DCA/MCH/DCA_y_vs_sign_vs_quadrant_vs_mom_corr"))->Fill(mchTrack.p(), quadrantMch, posNeg, dcay);
registry.get<THnSparse>(HIST("DCA/MCH/DCA_x_vs_sign_vs_quadrant_vs_mom_corr"))->Fill(mchTrack.p(), quadrant, posNeg, dcax);
registry.get<THnSparse>(HIST("DCA/MCH/DCA_y_vs_sign_vs_quadrant_vs_mom_corr"))->Fill(mchTrack.p(), quadrant, posNeg, dcay);
}
}

Expand All @@ -1795,25 +1802,25 @@
const auto& refTrackAtRefPlane = (iRefPlane == 0) ? mftTrackAtRefPlane : mchTrackAtRefPlane;

auto dx = mchTrackAtRefPlane.getX() - mftTrackAtRefPlane.getX();
dxPlots[iRefPlane]->Fill(dx, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrantMch, posNeg, mchTrack.p());
dxPlots[iRefPlane]->Fill(dx, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrant, posNeg, mchTrack.p());
auto dy = mchTrackAtRefPlane.getY() - mftTrackAtRefPlane.getY();
dyPlots[iRefPlane]->Fill(dy, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrantMch, posNeg, mchTrack.p());
dyPlots[iRefPlane]->Fill(dy, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrant, posNeg, mchTrack.p());

auto mftParamAtRefPlane = FwdtoMCH(mftTrackAtRefPlane);
auto mchParamAtRefPlane = FwdtoMCH(mchTrackAtRefPlane);

auto dsx = mchParamAtRefPlane.getNonBendingSlope() - mftParamAtRefPlane.getNonBendingSlope();
dsxPlots[iRefPlane]->Fill(dsx, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrantMch, posNeg, mchTrack.p());
dsxPlots[iRefPlane]->Fill(dsx, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrant, posNeg, mchTrack.p());
auto dsy = mchParamAtRefPlane.getBendingSlope() - mftParamAtRefPlane.getBendingSlope();
dsyPlots[iRefPlane]->Fill(dsy, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrantMch, posNeg, mchTrack.p());
dsyPlots[iRefPlane]->Fill(dsy, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrant, posNeg, mchTrack.p());

auto dphi = mchTrackAtRefPlane.getPhi() - mftTrackAtRefPlane.getPhi();
if (dphi < -TMath::Pi()) {
dphi += TMath::Pi() * 2.0;
} else if (dphi > TMath::Pi()) {
dphi -= TMath::Pi() * 2.0;
}
dphiPlots[iRefPlane]->Fill(dphi, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrantMch, posNeg, mchTrack.p());
dphiPlots[iRefPlane]->Fill(dphi, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrant, posNeg, mchTrack.p());
}
}
}
Expand Down
Loading