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); + } + } + } } }