From 9e7af6af26c02500cb1e6d7a4f0cc0674f69864b Mon Sep 17 00:00:00 2001 From: aferrero2707 Date: Mon, 16 Feb 2026 08:58:54 +0100 Subject: [PATCH] [PWGDQ] added DCA analysis with global forward tracks The DCA for global forward tracks is computed by extrapolating the matched MFT track using the momentum measured by the matched MCH track. The commit also adds some extra optional plots for the MFT standalone DCA analysis. --- PWGDQ/Tasks/muonGlobalAlignment.cxx | 285 ++++++++++++++++++++++------ 1 file changed, 228 insertions(+), 57 deletions(-) diff --git a/PWGDQ/Tasks/muonGlobalAlignment.cxx b/PWGDQ/Tasks/muonGlobalAlignment.cxx index a3fa52cfc24..e91112a629a 100644 --- a/PWGDQ/Tasks/muonGlobalAlignment.cxx +++ b/PWGDQ/Tasks/muonGlobalAlignment.cxx @@ -132,6 +132,8 @@ struct muonGlobalAlignment { Configurable fEnableVertexShiftAnalysis{"cfgEnableVertexShiftAnalysis", true, "Enable the analysis of vertex shift"}; Configurable fEnableMftDcaAnalysis{"cfgEnableMftDcaAnalysis", true, "Enable the analysis of DCA-based MFT alignment"}; + Configurable fEnableMftDcaExtraPlots{"cfgEnableMftDcaExtraPlots", false, "Enable additional plots for the analysis of DCA-based MFT alignment"}; + Configurable fEnableGlobalFwdDcaAnalysis{"cfgEnableGlobalFwdDcaAnalysis", true, "Enable the analysis of DCA-based MFT alignment using global forward tracks"}; Configurable fEnableMftMchResidualsAnalysis{"cfgEnableMftMchResidualsAnalysis", true, "Enable the analysis of residuals between MFT tracks and MCH clusters"}; int mRunNumber{0}; // needed to detect if the run changed and trigger update of magnetic field @@ -171,6 +173,8 @@ struct muonGlobalAlignment { double mBzAtMftCenter{0}; HistogramRegistry registry{"registry", {}}; + std::array mMftTrackEffNum; + std::array mMftTrackEffDen; // vector of all MFT-MCH(-MID) matching candidates associated to the same MCH(-MID) track, // to be sorted in descending order with respect to the matching quality @@ -362,17 +366,22 @@ struct muonGlobalAlignment { AxisSpec dcaxMCHAxis = {400, -10.0, 10.0, "DCA_{x} (cm)"}; AxisSpec dcayMCHAxis = {400, -10.0, 10.0, "DCA_{y} (cm)"}; AxisSpec dcazAxis = {20, -10.0, 10.0, "v_{z} (cm)"}; - AxisSpec txAxis = {30, -mftLadderWidth * 15.f / 2.f, mftLadderWidth * 15.f / 2.f, "track_{x} (cm)"}; - AxisSpec tyAxis = {20, -10.f, 10.f, "track_{y} (cm)"}; + AxisSpec txAxis = {30 * 4, -mftLadderWidth * 15.f / 2.f, mftLadderWidth * 15.f / 2.f, "track_{x} (cm)"}; + AxisSpec tyAxis = {24 * 4, -12.f, 12.f, "track_{y} (cm)"}; + AxisSpec txFineAxis = {1500, -15.f, 15.f, "track_{x} (cm)"}; + AxisSpec tyFineAxis = {1500, -15.f, 15.f, "track_{y} (cm)"}; AxisSpec vxAxis = {400, -0.5, 0.5, "vtx_{x} (cm)"}; AxisSpec vyAxis = {400, -0.5, 0.5, "vtx_{y} (cm)"}; AxisSpec vzAxis = {1000, -10.0, 10.0, "vtx_{z} (cm)"}; AxisSpec phiAxis = {36, -180.0, 180.0, "#phi (degrees)"}; + AxisSpec momAxis = {500, 0, 100.0, "p (GeV/c)"}; AxisSpec nMftClustersAxis = {6, 5, 11, "# of clusters"}; AxisSpec mftTrackTypeAxis = {2, 0, 2, "track type"}; + AxisSpec mftLayerAxis = {10, 0, 10, "layer"}; AxisSpec trackChargeSignAxis = {2, 0, 0, "sign"}; AxisSpec layersPatternAxis = {1024, 0, 1024, "layers pattern"}; AxisSpec zshiftAxis = {21, -5.25, 5.25, "z shift (mm)"}; + AxisSpec chi2Axis = {500, 0, 500, "chi2"}; registry.add("vertex_y_vs_x", std::format("Vertex y vs. x").c_str(), {HistType::kTH2F, {vxAxis, vyAxis}}); registry.add("vertex_z", std::format("Vertex z").c_str(), {HistType::kTH1F, {vzAxis}}); @@ -388,12 +397,36 @@ struct muonGlobalAlignment { } if (fEnableMftDcaAnalysis) { - registry.add("DCA/MFT/DCA_x", "DCA(x) vs. vz, tx, ty, nclus, trackType", - HistType::kTHnSparseF, {dcaxMFTAxis, dcazAxis, txAxis, tyAxis, nMftClustersAxis, mftTrackTypeAxis}); - registry.add("DCA/MFT/DCA_y", "DCA(y) vs. vz, tx, ty, nclus, trackType", - HistType::kTHnSparseF, {dcayMFTAxis, dcazAxis, txAxis, tyAxis, nMftClustersAxis, mftTrackTypeAxis}); - registry.add("DCA/MFT/layers", "Layers pattern vs. tx, ty, nclus, trackType", - HistType::kTHnSparseF, {layersPatternAxis, txAxis, tyAxis, nMftClustersAxis, mftTrackTypeAxis}); + registry.add("DCA/MFT/DCA_x", "DCA(x) vs. vz, tx, ty, nclus", + HistType::kTHnSparseF, {dcaxMFTAxis, dcazAxis, txAxis, tyAxis, nMftClustersAxis}); + registry.add("DCA/MFT/DCA_y", "DCA(y) vs. vz, tx, ty, nclus", + HistType::kTHnSparseF, {dcayMFTAxis, dcazAxis, txAxis, tyAxis, nMftClustersAxis}); + + if (fEnableMftDcaExtraPlots) { + registry.add("DCA/MFT/layers", "Layers vs. tx, ty, nclus", + HistType::kTHnSparseF, {mftLayerAxis, txAxis, tyAxis, nMftClustersAxis}); + registry.add("DCA/MFT/trackChi2", "Track #chi^{2} vs. tx, ty, nclus, layer", + HistType::kTHnSparseF, {chi2Axis, txAxis, tyAxis, nMftClustersAxis}); + registry.add("DCA/MFT/trackMomentum", "Track momentum vs. tx, ty, nclus, layer", + HistType::kTHnSparseF, {momAxis, txAxis, tyAxis, nMftClustersAxis}); + + const int nMftLayers = 10; + for (int i = 0; i < nMftLayers; i++) { + mMftTrackEffNum[i] = registry.add((std::string("DCA/MFT/mftTrackEffNum_") + std::to_string(i)).c_str(), + (std::string("Track efficiency num. - layer ") + std::to_string(i)).c_str(), + HistType::kTH2F, {txFineAxis, tyFineAxis}); + mMftTrackEffDen[i] = registry.add((std::string("DCA/MFT/mftTrackEffDen_") + std::to_string(i)).c_str(), + (std::string("Track efficiency den. - layer ") + std::to_string(i)).c_str(), + HistType::kTH2F, {txFineAxis, tyFineAxis}); + } + } + } + + if (fEnableGlobalFwdDcaAnalysis) { + registry.add("DCA/GlobalFwd/DCA_x", "DCA(x) vs. vz, tx, ty, nclus", + HistType::kTHnSparseF, {dcaxMFTAxis, dcazAxis, txAxis, tyAxis, nMftClustersAxis}); + registry.add("DCA/GlobalFwd/DCA_y", "DCA(y) vs. vz, tx, ty, nclus", + HistType::kTHnSparseF, {dcayMFTAxis, dcazAxis, txAxis, tyAxis, nMftClustersAxis}); } if (fEnableMftMchResidualsAnalysis) { @@ -934,6 +967,68 @@ struct muonGlobalAlignment { fwdtrack.setCovariances(transformedTrack.getCovariances()); } + template + T UpdateTrackMomentum(const T& track, const double p, int sign) + { + double px = p * sin(M_PI / 2 - atan(track.tgl())) * cos(track.phi()); + double py = p * sin(M_PI / 2 - atan(track.tgl())) * sin(track.phi()); + double pt = std::sqrt(std::pow(px, 2) + std::pow(py, 2)); + + SMatrix5 tpars = {track.x(), track.y(), track.phi(), track.tgl(), sign / pt}; + std::vector v1{0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0}; + SMatrix55 tcovs(v1.begin(), v1.end()); + + T newTrack; + newTrack.setParameters(tpars); + newTrack.setZ(track.z()); + newTrack.setCovariances(tcovs); + + return newTrack; + } + + template + T UpdateTrackMomentum(const T& track, const o2::mch::TrackParam& track4mom) + { + double px = track4mom.p() * sin(M_PI / 2 - atan(track.tgl())) * cos(track.phi()); + double py = track4mom.p() * sin(M_PI / 2 - atan(track.tgl())) * sin(track.phi()); + double pt = std::sqrt(std::pow(px, 2) + std::pow(py, 2)); + double sign = track4mom.getCharge(); + + SMatrix5 tpars = {track.x(), track.y(), track.phi(), track.tgl(), sign / pt}; + std::vector v1{0, 0, 0, 0, 0, + 0, 0, 0, 0, 0, + 0, 0, 0, 0, 0}; + SMatrix55 tcovs(v1.begin(), v1.end()); + + T newTrack; + newTrack.setParameters(tpars); + newTrack.setZ(track.z()); + newTrack.setCovariances(tcovs); + + return track; + } + + void UpdateTrackMomentum(o2::mch::TrackParam& track, const o2::mch::TrackParam& track4mom) + { + double pRatio = track.p() / track4mom.p(); + double newInvBendMom = track.getInverseBendingMomentum() * pRatio; + track.setInverseBendingMomentum(newInvBendMom); + track.setCharge(track4mom.getCharge()); + } + + void UpdateTrackMomentum(o2::track::TrackParCovFwd& track, const o2::mch::TrackParam& track4mom) + { + auto newTrackMCH = FwdtoMCH(track); + UpdateTrackMomentum(newTrackMCH, track4mom); + auto newTrack = MCHtoFwd(newTrackMCH); + + track.setParameters(newTrack.getParameters()); + track.setZ(newTrack.getZ()); + track.setCovariances(newTrack.getCovariances()); + } + o2::dataformats::GlobalFwdTrack PropagateMCHParam(mch::TrackParam mchTrack, const double z) { float absFront = -90.f; @@ -1041,7 +1136,7 @@ struct muonGlobalAlignment { } template - o2::dataformats::GlobalFwdTrack PropagateMFTToDCA(const TMFT& mftTrack, const C& collision, float zshift = 0) + o2::dataformats::GlobalFwdTrack PropagateMFTToDCA(const TMFT& mftTrack, const C& collision, float zshift) { static double Bz = -10001; double chi2 = mftTrack.chi2(); @@ -1084,55 +1179,53 @@ struct muonGlobalAlignment { return propmuon; } - template - T UpdateTrackMomentum(const T& track, const double p, int sign) + template + o2::dataformats::GlobalFwdTrack PropagateMFTToDCA(const TMFT& mftTrack, const TMUON& mchTrack, const C& collision, float zshift) { - double px = p * sin(M_PI / 2 - atan(track.tgl())) * cos(track.phi()); - double py = p * sin(M_PI / 2 - atan(track.tgl())) * sin(track.phi()); - double pt = std::sqrt(std::pow(px, 2) + std::pow(py, 2)); - - SMatrix5 tpars = {track.x(), track.y(), track.phi(), track.tgl(), sign / pt}; + static double Bz = -10001; + double chi2 = mftTrack.chi2(); + double phiCorrDeg = 0; + double phiCorr = phiCorrDeg * TMath::Pi() / 180.f; + double tR = std::hypot(mftTrack.x(), mftTrack.y()); + double tphi = std::atan2(mftTrack.y(), mftTrack.x()); + double tx = std::cos(tphi + phiCorr) * tR; + double ty = std::sin(tphi + phiCorr) * tR; + SMatrix5 tpars = {tx, ty, mftTrack.phi() + phiCorr, mftTrack.tgl(), mftTrack.signed1Pt()}; std::vector v1{0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}; SMatrix55 tcovs(v1.begin(), v1.end()); + o2::track::TrackParCovFwd fwdtrack{mftTrack.z(), tpars, tcovs, chi2}; + if (configMFTAlignmentCorrections.fEnableMFTAlignmentCorrections) { + TransformMFT(fwdtrack); + } - T newTrack; - newTrack.setParameters(tpars); - newTrack.setZ(track.z()); - newTrack.setCovariances(tcovs); - - return newTrack; - } - - template - T UpdateTrackMomentum(const T& track, const o2::mch::TrackParam& track4mom) - { - double px = track4mom.p() * sin(M_PI / 2 - atan(track.tgl())) * cos(track.phi()); - double py = track4mom.p() * sin(M_PI / 2 - atan(track.tgl())) * sin(track.phi()); - double pt = std::sqrt(std::pow(px, 2) + std::pow(py, 2)); - double sign = track4mom.getCharge(); + // extrapolation with MCH tools + auto mchTrackAtMFT = FwdtoMCH(FwdToTrackPar(mchTrack)); + o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchTrackAtMFT, mftTrack.z()); + UpdateTrackMomentum(fwdtrack, mchTrackAtMFT); - SMatrix5 tpars = {track.x(), track.y(), track.phi(), track.tgl(), sign / pt}; - std::vector v1{0, 0, 0, 0, 0, - 0, 0, 0, 0, 0, - 0, 0, 0, 0, 0}; - SMatrix55 tcovs(v1.begin(), v1.end()); + // double propVec[3] = {}; + // propVec[0] = collision.posX() - mftTrack.x(); + // propVec[1] = collision.posY() - mftTrack.y(); + // propVec[2] = collision.posZ() - mftTrack.z(); - T newTrack; - newTrack.setParameters(tpars); - newTrack.setZ(track.z()); - newTrack.setCovariances(tcovs); + // double centerZ[3] = {mftTrack.x() + propVec[0] / 2., + // mftTrack.y() + propVec[1] / 2., + // mftTrack.z() + propVec[2] / 2.}; + if (Bz < -10000) { + double centerZ[3] = {0, 0, -45.f / 2.f}; + o2::field::MagneticField* field = static_cast(TGeoGlobalMagField::Instance()->GetField()); + Bz = field->getBz(centerZ); + } + fwdtrack.propagateToZ(collision.posZ() - zshift, Bz); - return track; - } + o2::dataformats::GlobalFwdTrack propmuon; + propmuon.setParameters(fwdtrack.getParameters()); + propmuon.setZ(fwdtrack.getZ()); + propmuon.setCovariances(fwdtrack.getCovariances()); - void UpdateTrackMomentum(o2::mch::TrackParam& track, const o2::mch::TrackParam& track4mom) - { - double pRatio = track.p() / track4mom.p(); - double newInvBendMom = track.getInverseBendingMomentum() * pRatio; - track.setInverseBendingMomentum(newInvBendMom); - track.setCharge(track4mom.getCharge()); + return propmuon; } template @@ -1170,7 +1263,7 @@ struct muonGlobalAlignment { void FillDCAPlots(MyEvents const& collisions, MyBCs const& bcs, - MyMuonsWithCov const& /*muonTracks*/, + MyMuonsWithCov const& muonTracks, MyMFTs const& mftTracks, const std::map& collisionInfos) { @@ -1203,6 +1296,10 @@ struct muonGlobalAlignment { for (auto mftIndex : mftTrackIds) { auto const& mftTrack = mftTracks.rawIteratorAt(mftIndex); + if (mftTrack.isCA()) { + continue; + } + bool isGoodMFT = IsGoodMFT(mftTrack, 999.f, 5); if (!isGoodMFT) continue; @@ -1212,26 +1309,55 @@ struct muonGlobalAlignment { double dcay = mftTrackAtDCA.getY() - collision.posY(); double phi = mftTrack.phi() * 180 / TMath::Pi(); int mftNclusters = mftTrack.nClusters(); - int mftTrackType = mftTrack.isCA() ? 1 : 0; + double chi2NDF = static_cast(mftNclusters) * 2 - 5; const int nMftLayers = 10; - int layerPattern = 0; + std::array firedLayers; for (int layer = 0; layer < nMftLayers; layer++) { if ((mftTrack.mftClusterSizesAndTrackFlags() >> (layer * 6)) & 0x3F) { - layerPattern += (1 << layer); + firedLayers[layer] = true; + } else { + firedLayers[layer] = false; } } if (fEnableMftDcaAnalysis) { - registry.get(HIST("DCA/MFT/DCA_y_vs_x"))->Fill(dcax, dcay); - registry.get(HIST("DCA/MFT/DCA_x"))->Fill(dcax, collision.posZ(), mftTrack.x(), mftTrack.y(), mftNclusters, mftTrackType); - registry.get(HIST("DCA/MFT/DCA_y"))->Fill(dcay, collision.posZ(), mftTrack.x(), mftTrack.y(), mftNclusters, mftTrackType); + const int nMftLayers = 10; + if (fEnableMftDcaExtraPlots) { + for (int i = 0; i < nMftLayers; i++) { + if (firedLayers[i]) { + registry.get(HIST("DCA/MFT/trackChi2"))->Fill(mftTrack.chi2() / chi2NDF, mftTrack.x(), mftTrack.y(), mftNclusters, i); + } + } + } - registry.get(HIST("DCA/MFT/layers"))->Fill(layerPattern, mftTrack.x(), mftTrack.y(), mftNclusters, mftTrackType); + if (mftTrack.chi2() <= fTrackChi2MftUp) { + registry.get(HIST("DCA/MFT/DCA_y_vs_x"))->Fill(dcax, dcay); + registry.get(HIST("DCA/MFT/DCA_x"))->Fill(dcax, collision.posZ(), mftTrack.x(), mftTrack.y(), mftNclusters); + registry.get(HIST("DCA/MFT/DCA_y"))->Fill(dcay, collision.posZ(), mftTrack.x(), mftTrack.y(), mftNclusters); + + if (fEnableMftDcaExtraPlots) { + if (mftNclusters >= 6) { + for (int i = 0; i < nMftLayers; i++) { + auto mftTrackAtLayer = PropagateMFT(mftTrack, o2::mft::constants::mft::LayerZCoordinate()[i]); + std::get>(mMftTrackEffDen[i])->Fill(mftTrackAtLayer.getX(), mftTrackAtLayer.getY()); + if (firedLayers[i]) { + std::get>(mMftTrackEffNum[i])->Fill(mftTrackAtLayer.getX(), mftTrackAtLayer.getY()); + registry.get(HIST("DCA/MFT/layers"))->Fill(i, mftTrack.x(), mftTrack.y(), mftNclusters); + } + } + } + for (int i = 0; i < nMftLayers; i++) { + if (firedLayers[i]) { + registry.get(HIST("DCA/MFT/trackMomentum"))->Fill(mftTrack.p(), mftTrack.x(), mftTrack.y(), mftNclusters, i); + } + } + } + } } if (fEnableVertexShiftAnalysis) { - if (std::fabs(collision.posZ()) < 1.f) { + if (mftTrack.chi2() <= fTrackChi2MftUp && std::fabs(collision.posZ()) < 1.f) { float zshift[21] = {// in millimeters -5.0, -4.5, -4.0, -3.5, -3.0, -2.5, -2.0, -1.5, -1.0, -0.5, 0.0, 0.5, 1.0, 1.5, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0}; @@ -1246,6 +1372,51 @@ struct muonGlobalAlignment { } } } + + if (fEnableGlobalFwdDcaAnalysis) { + // loop over global muon tracks + for (auto& [muonIndex, globalTracksVector] : collisionInfo.globalMuonTracks) { + auto const& muonTrack = muonTracks.rawIteratorAt(globalTracksVector[0]); + const auto& mchTrack = muonTrack.template matchMCHTrack_as(); + const auto& mftTrack = muonTrack.template matchMFTTrack_as(); + + if (muonTrack.chi2MatchMCHMFT() < 50) { + continue; + } + + if (globalTracksVector.size() > 1) { + auto const& muonTrack2 = muonTracks.rawIteratorAt(globalTracksVector[1]); + double dchi2 = muonTrack2.chi2MatchMCHMFT() - muonTrack.chi2MatchMCHMFT(); + + if (dchi2 < 50) { + continue; + } + } + + if (mftTrack.isCA()) { + continue; + } + + bool isGoodMFT = IsGoodMFT(mftTrack, 999.f, 5); + if (!isGoodMFT) { + continue; + } + + bool isGoodMuon = IsGoodMuon(mchTrack, collision, fTrackChi2MchUp, 0.f, fPtMchLow, {fEtaMftLow, fEtaMftUp}, {fRabsLow, fRabsUp}, fSigmaPdcaUp); + if (!isGoodMuon) + continue; + + auto mftTrackAtDCA = PropagateMFTToDCA(mftTrack, mchTrack, collision, fVertexZshift); + double dcax = mftTrackAtDCA.getX() - collision.posX(); + double dcay = mftTrackAtDCA.getY() - collision.posY(); + int mftNclusters = mftTrack.nClusters(); + + if (mftTrack.chi2() <= fTrackChi2MftUp) { + registry.get(HIST("DCA/GlobalFwd/DCA_x"))->Fill(dcax, collision.posZ(), mftTrack.x(), mftTrack.y(), mftNclusters); + registry.get(HIST("DCA/GlobalFwd/DCA_y"))->Fill(dcay, collision.posZ(), mftTrack.x(), mftTrack.y(), mftNclusters); + } + } + } } }