From 5f7e1d1dcd5de700e610b15e761d515696573f23 Mon Sep 17 00:00:00 2001 From: aferrero2707 Date: Tue, 28 Apr 2026 16:27:29 +0200 Subject: [PATCH] [PWGDQ] added track slope dependence to residals plots The X and Y residual plots now have an additional dimension filled with the track slope in the same coordinate. The residual vs. slope dependencies are needed to identify and correct possible z misalignments of the MCH detection elements. In addition, a cut has been added to the maximum MFT-MCH matching chi2 in order to select only good matches for the global alignment analysis. Finally, the quadrants are now computed from the MFT track position instead of the MCH one, since we noticed some bias in the MCH DCA in MC simulations. --- PWGDQ/Tasks/muonGlobalAlignment.cxx | 67 ++++++++++++++++------------- 1 file changed, 37 insertions(+), 30 deletions(-) diff --git a/PWGDQ/Tasks/muonGlobalAlignment.cxx b/PWGDQ/Tasks/muonGlobalAlignment.cxx index 00bd2d964d7..5f9dcbea5a5 100644 --- a/PWGDQ/Tasks/muonGlobalAlignment.cxx +++ b/PWGDQ/Tasks/muonGlobalAlignment.cxx @@ -528,8 +528,8 @@ struct muonGlobalAlignment { } 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}}); @@ -543,20 +543,20 @@ struct muonGlobalAlignment { 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(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(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(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(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(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(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(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(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(getNumDE()), "DE"}}}); @@ -1656,7 +1656,8 @@ struct muonGlobalAlignment { auto const& muonTrack = muonTracks.rawIteratorAt(globalTracksVector[0]); const auto& mchTrack = muonTrack.template matchMCHTrack_as(); const auto& mftTrack = muonTrack.template matchMFTTrack_as(); - 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); @@ -1667,6 +1668,10 @@ struct muonGlobalAlignment { if (!isGoodMFT) continue; + double matchChi2 = muonTrack.chi2MatchMCHMFT() / 5.f; + if (matchChi2 > 10.f) + continue; + // refit MCH track if enabled TrackRealigned convertedTrack; bool convertedTrackOk = false; @@ -1722,15 +1727,16 @@ struct muonGlobalAlignment { // 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 xPos{master.x(), mftTrackAtCluster.getX()}; std::array yPos{master.y(), mftTrackAtCluster.getY()}; - registry.get(HIST("residuals/dx_vs_chamber"))->Fill(chamber + 1, quadrantMch, posNeg, xPos[0] - xPos[1]); - registry.get(HIST("residuals/dy_vs_chamber"))->Fill(chamber + 1, quadrantMch, posNeg, yPos[0] - yPos[1]); + registry.get(HIST("residuals/dx_vs_chamber"))->Fill(chamber + 1, quadrant, posNeg, xPos[0] - xPos[1]); + registry.get(HIST("residuals/dy_vs_chamber"))->Fill(chamber + 1, quadrant, posNeg, yPos[0] - yPos[1]); - registry.get(HIST("residuals/dx_vs_de"))->Fill(xPos[0] - xPos[1], deIndex, quadrantMch, posNeg, mchTrack.p()); - registry.get(HIST("residuals/dy_vs_de"))->Fill(yPos[0] - yPos[1], deIndex, quadrantMch, posNeg, mchTrack.p()); + registry.get(HIST("residuals/dx_vs_de"))->Fill(xPos[0] - xPos[1], deIndex, quadrant, posNeg, mchTrack.p(), mftTrackParamAtCluster.getNonBendingSlope()); + registry.get(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 @@ -1739,15 +1745,16 @@ struct muonGlobalAlignment { // alignment (if realignment is enabled) if (convertedTrackWithCorrOk) { auto mftTrackAtClusterWithCorr = PropagateMFTtoMCH(mftTrack, mch::TrackParam(convertedTrackWithCorr.first()), masterWithCorr.z()); + auto mftTrackParamAtClusterWithCorr = FwdtoMCH(mftTrackAtClusterWithCorr); std::array xPos{masterWithCorr.x(), mftTrackAtClusterWithCorr.getX()}; std::array yPos{masterWithCorr.y(), mftTrackAtClusterWithCorr.getY()}; - registry.get(HIST("residuals/dx_vs_chamber_corr"))->Fill(chamber + 1, quadrantMch, posNeg, xPos[0] - xPos[1]); - registry.get(HIST("residuals/dy_vs_chamber_corr"))->Fill(chamber + 1, quadrantMch, posNeg, yPos[0] - yPos[1]); + registry.get(HIST("residuals/dx_vs_chamber_corr"))->Fill(chamber + 1, quadrant, posNeg, xPos[0] - xPos[1]); + registry.get(HIST("residuals/dy_vs_chamber_corr"))->Fill(chamber + 1, quadrant, posNeg, yPos[0] - yPos[1]); - registry.get(HIST("residuals/dx_vs_de_corr"))->Fill(xPos[0] - xPos[1], deIndex, quadrantMch, posNeg, mchTrack.p()); - registry.get(HIST("residuals/dy_vs_de_corr"))->Fill(yPos[0] - yPos[1], deIndex, quadrantMch, posNeg, mchTrack.p()); + registry.get(HIST("residuals/dx_vs_de_corr"))->Fill(xPos[0] - xPos[1], deIndex, quadrant, posNeg, mchTrack.p(), mftTrackParamAtClusterWithCorr.getNonBendingSlope()); + registry.get(HIST("residuals/dy_vs_de_corr"))->Fill(yPos[0] - yPos[1], deIndex, quadrant, posNeg, mchTrack.p(), mftTrackParamAtClusterWithCorr.getBendingSlope()); } } @@ -1757,12 +1764,12 @@ struct muonGlobalAlignment { auto dcay = mchTrackAtDCA.getY() - collision.posY(); registry.get(HIST("DCA/MCH/DCA_y_vs_x"))->Fill(dcax, dcay); - registry.get(HIST("DCA/MCH/DCA_x_vs_sign_vs_quadrant_vs_mom"))->Fill(mchTrack.p(), quadrantMch, posNeg, dcax); - registry.get(HIST("DCA/MCH/DCA_y_vs_sign_vs_quadrant_vs_mom"))->Fill(mchTrack.p(), quadrantMch, posNeg, dcay); + registry.get(HIST("DCA/MCH/DCA_x_vs_sign_vs_quadrant_vs_mom"))->Fill(mchTrack.p(), quadrant, posNeg, dcax); + registry.get(HIST("DCA/MCH/DCA_y_vs_sign_vs_quadrant_vs_mom"))->Fill(mchTrack.p(), quadrant, posNeg, dcay); if (fEnableMftMchResidualsExtraPlots) { - registry.get(HIST("DCA/MCH/DCA_x_vs_sign_vs_quadrant_vs_vz"))->Fill(collision.posZ(), quadrantMch, posNeg, dcax); - registry.get(HIST("DCA/MCH/DCA_y_vs_sign_vs_quadrant_vs_vz"))->Fill(collision.posZ(), quadrantMch, posNeg, dcay); + registry.get(HIST("DCA/MCH/DCA_x_vs_sign_vs_quadrant_vs_vz"))->Fill(collision.posZ(), quadrant, posNeg, dcax); + registry.get(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(HIST("residuals/dphi_at_mft"))->Fill(deltaPhi, mftTrack.x(), mftTrack.y(), posNeg, mchTrackAtMFT.getP()); @@ -1774,8 +1781,8 @@ struct muonGlobalAlignment { auto dcax = mchTrackAtDCA.getX() - collision.posX(); auto dcay = mchTrackAtDCA.getY() - collision.posY(); - registry.get(HIST("DCA/MCH/DCA_x_vs_sign_vs_quadrant_vs_mom_corr"))->Fill(mchTrack.p(), quadrantMch, posNeg, dcax); - registry.get(HIST("DCA/MCH/DCA_y_vs_sign_vs_quadrant_vs_mom_corr"))->Fill(mchTrack.p(), quadrantMch, posNeg, dcay); + registry.get(HIST("DCA/MCH/DCA_x_vs_sign_vs_quadrant_vs_mom_corr"))->Fill(mchTrack.p(), quadrant, posNeg, dcax); + registry.get(HIST("DCA/MCH/DCA_y_vs_sign_vs_quadrant_vs_mom_corr"))->Fill(mchTrack.p(), quadrant, posNeg, dcay); } } @@ -1795,17 +1802,17 @@ struct muonGlobalAlignment { 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()) { @@ -1813,7 +1820,7 @@ struct muonGlobalAlignment { } 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()); } } }