From 46eb94f0ddfd84e8ec64b178d11d7f5d2ba9cbdc Mon Sep 17 00:00:00 2001 From: aferrero2707 Date: Wed, 12 Aug 2026 10:42:27 +0200 Subject: [PATCH] [PWGDQ] add dimuon analysis in global alignment task Added comparison of dimuon properties (invariant mass, angles, DCA) for different alignment configurations. --- PWGDQ/Tasks/muonGlobalAlignment.cxx | 1259 ++++++++++++++++++--------- 1 file changed, 856 insertions(+), 403 deletions(-) diff --git a/PWGDQ/Tasks/muonGlobalAlignment.cxx b/PWGDQ/Tasks/muonGlobalAlignment.cxx index 68671d16fdf..c2577f1e05e 100644 --- a/PWGDQ/Tasks/muonGlobalAlignment.cxx +++ b/PWGDQ/Tasks/muonGlobalAlignment.cxx @@ -18,6 +18,7 @@ #include "Common/CCDB/EventSelectionParams.h" #include "Common/CCDB/RCTSelectionFlags.h" #include "Common/Core/RecoDecay.h" +#include "Common/Core/fwdtrackUtilities.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/TrackSelectionTables.h" @@ -80,6 +81,7 @@ #include #include #include +#include #include using namespace o2; @@ -110,8 +112,6 @@ using MyMFTCovariance = MyMFTCovariances::iterator; using SMatrix55 = ROOT::Math::SMatrix>; using SMatrix5 = ROOT::Math::SVector; -// static o2::globaltracking::MatchGlobalFwd sExtrap; - using o2::dataformats::GlobalFwdTrack; using o2::track::TrackParCovFwd; @@ -142,8 +142,10 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc //// Variables for selecting MCH and MFT tracks Configurable cfgTrackChi2MchUp{"cfgTrackChi2MchUp", 5.f, ""}; Configurable cfgPtMchLow{"cfgPtMchLow", 0.7f, ""}; - Configurable cfgEtaMftlow{"cfgEtaMftlow", -3.6f, ""}; - Configurable cfgEtaMftup{"cfgEtaMftup", -2.5f, ""}; + Configurable cfgEtaMchLow{"cfgEtaMchLow", -4.0f, ""}; + Configurable cfgEtaMchUp{"cfgEtaMchUp", -2.5f, ""}; + Configurable cfgEtaMftLow{"cfgEtaMftLow", -3.6f, ""}; + Configurable cfgEtaMftUp{"cfgEtaMftUp", -2.5f, ""}; Configurable cfgRabsLow{"cfgRabsLow", 17.6f, ""}; Configurable cfgRabsUp{"cfgRabsUp", 89.5f, ""}; Configurable fSigmaPdcaUp{"fSigmaPdcaUp", 6.f, ""}; @@ -151,7 +153,9 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc Configurable cfgTrackNClustMftLow{"cfgTrackNClustMftLow", 7, ""}; Configurable cfgTrackChi2MftUp{"cfgTrackChi2MftUp", 999.f, ""}; - Configurable cfgMftDcaMatchChi2Up{"cfgMftDcaMatchChi2Up", 10.f, ""}; + Configurable cfgMftDcaMatchChi2Up{"cfgMftDcaMatchChi2Up", 50.f, ""}; + Configurable cfgMftMchResidualsMatchChi2Up{"cfgMftMchResidualsMatchChi2Up", 50.f, ""}; + Configurable cfgDimuonMatchChi2Up{"cfgDimuonMatchChi2Up", 50.f, ""}; Configurable cfgMftMchResidualsPLow{"cfgMftMchResidualsPLow", 30.f, ""}; Configurable cfgMftMchResidualsPtLow{"cfgMftMchResidualsPtLow", 4.f, ""}; @@ -180,7 +184,8 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc //// Variables for re-alignment setup struct : ConfigurableGroup { - Configurable cfgEnableMCHRealign{"cfgEnableMCHRealign", true, "Enable re-alignment of MCH clusters and tracks"}; + Configurable cfgEnableMCHRefit{"cfgEnableMCHRefit", false, "Enable re-fitting of MCH tracks"}; + Configurable cfgEnableMCHRealign{"cfgEnableMCHRealign", false, "Enable re-alignment of MCH clusters and tracks"}; Configurable cfgChamberResolutionX{"cfgChamberResolutionX", 0.4, "Chamber resolution along X configuration for refit"}; // 0.4cm pp, 0.2cm PbPb Configurable cfgChamberResolutionY{"cfgChamberResolutionY", 0.4, "Chamber resolution along Y configuration for refit"}; // 0.4cm pp, 0.2cm PbPb Configurable cfgSigmaCutImprove{"cfgSigmaCutImprove", 6., "Sigma cut for track improvement"}; @@ -208,6 +213,7 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc Configurable cfgEnableMftMchResidualsAnalysis{"cfgEnableMftMchResidualsAnalysis", true, "Enable the analysis of residuals between MFT tracks and MCH clusters"}; Configurable cfgEnableMftMchResidualsExtraPlots{"cfgEnableMftMchResidualsExtraPlots", false, "Enable additional plots for the analysis of residuals between MFT tracks and MCH clusters"}; Configurable cfgEnableMftMchMatchingAnalysis{"cfgEnableMftMchMatchingAnalysis", false, "Enable the analysis of residuals between MFT and MCH tracks at reference planes"}; + Configurable cfgEnableDimuonAnalysis{"cfgEnableDimuonAnalysis", false, "Enable the analysis of di-muon pairs"}; Configurable cfgRefPlaneZMFT{"cfgRefPlaneZMFT", o2::mft::constants::mft::LayerZCoordinate()[0], "Reference plane on MFT side"}; Configurable cfgRefPlaneZMCH{"cfgRefPlaneZMCH", -526.0, "Reference plane on MCH side"}; @@ -234,6 +240,52 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc } }; + class TrackParExt : public o2::track::TrackParCovFwd + { + public: + TrackParExt() = default; + TrackParExt(const TrackParExt& t) = default; + explicit TrackParExt(o2::track::TrackParCovFwd const& t, int nc = -1, bool r = false) + : TrackParCovFwd(t), nClusters(nc), removable(r) {} + ~TrackParExt() = default; + + TrackParExt& operator=(const TrackParCovFwd& tpf) + { + o2::track::TrackParCovFwd::operator=(tpf); + return *this; + } + TrackParExt& operator=(const TrackParExt& tpe) + { + o2::track::TrackParCovFwd::operator=(tpe); + nClusters = tpe.getNClusters(); + removable = tpe.isRemovable(); + return *this; + } + + void setNClusters(int n) { nClusters = n; } + [[nodiscard]] int getNClusters() const { return nClusters; } + + void setRemovable() { removable = true; } + [[nodiscard]] bool isRemovable() const { return removable; } + + [[nodiscard]] o2::track::TrackParCovFwd asTrackParCovFwd() const + { + return {static_cast(*this)}; + } + + private: + int nClusters{-1}; + bool removable{false}; + }; + + std::unordered_map mMchTrackPars; + std::unordered_map mMftTrackPars; + std::unordered_map mMchTrackParsNew; + std::unordered_map mMftTrackParsNew; + + using MuonPair = std::pair; + using GlobalMuonPair = std::pair, std::vector>; + geo::TransformationCreator transformation; std::map transformRef; // reference geometry w.r.t track data std::map transformNew; // new geometry @@ -278,98 +330,6 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc std::map> globalMuonTracks; }; - void InitCollisions(MyEvents const& collisions, - MyBCs const& bcs, - MyMuonsWithCov const& muonTracks, - std::map& collisionInfos) - { - // fill collision information for global muon tracks (MFT-MCH-MID matches) - for (const auto& muonTrack : muonTracks) { - if (!muonTrack.has_collision()) { - continue; - } - - auto collision = collisions.rawIteratorAt(muonTrack.collisionId()); - - if (cfgRequireGoodRCT && !rctChecker(collision)) { - continue; - } - - uint64_t collisionIndex = collision.globalIndex(); - - auto bc = bcs.rawIteratorAt(collision.bcId()); - - auto& collisionInfo = collisionInfos[collisionIndex]; - collisionInfo.bc = bc.globalBC(); - collisionInfo.zVertex = collision.posZ(); - - if (static_cast(muonTrack.trackType()) > GlobalTrackTypeMax) { - // standalone MCH or MCH-MID tracks - uint64_t mchTrackIndex = muonTrack.globalIndex(); - collisionInfo.mchTracks.push_back(mchTrackIndex); - } else { - // global muon tracks (MFT-MCH or MFT-MCH-MID) - uint64_t muonTrackIndex = muonTrack.globalIndex(); - auto const& mchTrack = muonTrack.template matchMCHTrack_as(); - uint64_t mchTrackIndex = mchTrack.globalIndex(); - - // check if a vector of global muon candidates is already available for the current MCH index - // if not, initialize a new one and add the current global muon track - // bool globalMuonTrackFound = false; - auto matchingCandidateIterator = collisionInfo.globalMuonTracks.find(mchTrackIndex); - if (matchingCandidateIterator != collisionInfo.globalMuonTracks.end()) { - matchingCandidateIterator->second.push_back(muonTrackIndex); - // globalMuonTrackFound = true; - } else { - collisionInfo.globalMuonTracks[mchTrackIndex].push_back(muonTrackIndex); - } - } - } - - // sort the vectors of matching candidates in ascending order based on the matching chi2 value - auto compareChi2 = [&muonTracks](uint64_t trackIndex1, uint64_t trackIndex2) -> bool { - auto const& track1 = muonTracks.rawIteratorAt(trackIndex1); - auto const& track2 = muonTracks.rawIteratorAt(trackIndex2); - - return (track1.chi2MatchMCHMFT() < track2.chi2MatchMCHMFT()); - }; - - for (auto& [collisionIndex, collisionInfo] : collisionInfos) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) - for (auto& [mchIndex, globalTracksVector] : collisionInfo.globalMuonTracks) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) - std::sort(globalTracksVector.begin(), globalTracksVector.end(), compareChi2); - } - } - } - - void InitCollisions(MyEvents const& collisions, - MyBCs const& bcs, - MyMuonsWithCov const& muonTracks, - MyMFTs const& mftTracks, - std::map& collisionInfos) - { - InitCollisions(collisions, bcs, muonTracks, collisionInfos); - - // fill collision information for MFT standalone tracks - for (const auto& mftTrack : mftTracks) { - if (!mftTrack.has_collision()) { - continue; - } - - auto collision = collisions.rawIteratorAt(mftTrack.collisionId()); - uint64_t collisionIndex = collision.globalIndex(); - - auto bc = bcs.rawIteratorAt(collision.bcId()); - - uint64_t mftTrackIndex = mftTrack.globalIndex(); - - auto& collisionInfo = collisionInfos[collisionIndex]; - collisionInfo.bc = bc.globalBC(); - collisionInfo.zVertex = collision.posZ(); - - collisionInfo.mftTracks.push_back(mftTrackIndex); - } - } - template void initCCDB(BC const& bc) { @@ -545,8 +505,8 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc } if (cfgEnableMftMchResidualsAnalysis) { - AxisSpec dxAxis = {400, -20.0, 20.0, "#Delta x (cm)"}; - AxisSpec dyAxis = {400, -20.0, 20.0, "#Delta y (cm)"}; + AxisSpec dxAxis = {400, -10.0, 10.0, "#Delta x (cm)"}; + AxisSpec dyAxis = {400, -10.0, 10.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}}); @@ -628,6 +588,48 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc registry.add("matching/dphiAtMCH", "Tracks #Delta#phi at MCH reference plane", {HistType::kTHnSparseF, {dphiAxis, {20, -100.0, 100.0, "track x (cm)"}, {20, -100.0, 100.0, "track y (cm)"}, {4, 0, 4, "quadrant"}, {2, 0, 2, "sign"}, {20, 0, 100.0, "p (GeV/c)"}}}); } + + if (cfgEnableDimuonAnalysis) { + AxisSpec invMassAxis = {1500, 0, 15, "M_{#mu^{+}#mu^{-}} (GeV/c^{2})"}; + AxisSpec pTAxis = {30, 0, 30, "#mu^{+}#mu^{-} p_{T} (GeV/c)"}; + AxisSpec pAxis = {50, 0, 200, "#mu^{+}#mu^{-} p (GeV/c)"}; + AxisSpec muPosQuadrantAxis = {4, 0, 4, "#mu^{+} quadrant"}; + AxisSpec muNegQuadrantAxis = {4, 0, 4, "#mu^{-} quadrant"}; + AxisSpec angleAxis = {100, 0, 0.5, "#mu^{+}#mu^{-} angle (rad)"}; + AxisSpec angleDiffAxis = {500, -0.05, 0.05, "#mu^{+}#mu^{-} angle difference (rad)"}; + AxisSpec dcaxAxis = {400, -10.0, 10.0, "#mu^{+}#mu^{-} DCA_{x} (cm)"}; + AxisSpec dcayAxis = {400, -10.0, 10.0, "#mu^{+}#mu^{-} DCA_{y} (cm)"}; + AxisSpec dcaMftxAxis = {400, -0.5, 0.5, "#mu^{+}#mu^{-} DCA_{x} (cm)"}; + AxisSpec dcaMftyAxis = {400, -0.5, 0.5, "#mu^{+}#mu^{-} DCA_{y} (cm)"}; + + // di-muon invariant mass distributions + registry.add("dimuon/invariantMass_MuonKine_MuonCuts", "#mu^{+}#mu^{-} invariant mass (muon cuts)", {HistType::kTHnSparseF, {invMassAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + registry.add("dimuon/invariantMass_MuonKine_GlobalMuonCuts_GoodMatches", "#mu^{+}#mu^{-} invariant mass (global muon cuts, good matches)", {HistType::kTHnSparseF, {invMassAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + registry.add("dimuon/invariantMass_ScaledMftKine_GlobalMuonCuts_GoodMatches", "#mu^{+}#mu^{-} invariant mass (global muon cuts, rescaled MFT, good matches)", {HistType::kTHnSparseF, {invMassAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + // difference in mu+mu- opening angle between MCH and global muon tracks + registry.add("dimuon/angle_GlobalMuonCuts_GoodMatches", "#mu^{+}#mu^{-} opening angle difference (global muon cuts, good matches)", {HistType::kTHnSparseF, {angleDiffAxis, angleAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + // mu+mu- DCA + registry.add("dimuon/dcax_MuonKine_MuonCuts", "#mu^{+}#mu^{-} DCA_{x} (muon cuts)", {HistType::kTHnSparseF, {dcaxAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + registry.add("dimuon/dcay_MuonKine_MuonCuts", "#mu^{+}#mu^{-} DCA_{y} (muon cuts)", {HistType::kTHnSparseF, {dcayAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + registry.add("dimuon/dcax_MuonKine_GlobalMuonCuts_GoodMatches", "#mu^{+}#mu^{-} DCA_{x} (global muon cuts, good matches)", {HistType::kTHnSparseF, {dcaxAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + registry.add("dimuon/dcay_MuonKine_GlobalMuonCuts_GoodMatches", "#mu^{+}#mu^{-} DCA_{y} (global muon cuts, good matches)", {HistType::kTHnSparseF, {dcayAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + registry.add("dimuon/dcax_ScaledMftKine_GlobalMuonCuts_GoodMatches", "#mu^{+}#mu^{-} DCA_{x} (global muon cuts, rescaled MFT, good matches)", {HistType::kTHnSparseF, {dcaMftxAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + registry.add("dimuon/dcay_ScaledMftKine_GlobalMuonCuts_GoodMatches", "#mu^{+}#mu^{-} DCA_{y} (global muon cuts, rescaled MFT, good matches)", {HistType::kTHnSparseF, {dcaMftyAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + + // di-muon invariant mass distributions (realigned/refitted tracks) + registry.add("dimuon/realign/invariantMass_MuonKine_MuonCuts", "#mu^{+}#mu^{-} invariant mass (muon cuts)", {HistType::kTHnSparseF, {invMassAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + registry.add("dimuon/realign/invariantMass_MuonKine_GlobalMuonCuts_GoodMatches", "#mu^{+}#mu^{-} invariant mass (global muon cuts, good matches)", {HistType::kTHnSparseF, {invMassAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + registry.add("dimuon/realign/invariantMass_ScaledMftKine_GlobalMuonCuts_GoodMatches", "#mu^{+}#mu^{-} invariant mass (global muon cuts, rescaled MFT, good matches)", {HistType::kTHnSparseF, {invMassAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + // difference in mu+mu- opening angle between MCH and global muon tracks (realigned/refitted tracks) + registry.add("dimuon/realign/angle_GlobalMuonCuts_GoodMatches", "#mu^{+}#mu^{-} opening angle difference (global muon cuts, good matches)", {HistType::kTHnSparseF, {angleDiffAxis, angleAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + // mu+mu- DCA + registry.add("dimuon/realign/dcax_MuonKine_MuonCuts", "#mu^{+}#mu^{-} DCA_{x} (muon cuts)", {HistType::kTHnSparseF, {dcaxAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + registry.add("dimuon/realign/dcay_MuonKine_MuonCuts", "#mu^{+}#mu^{-} DCA_{y} (muon cuts)", {HistType::kTHnSparseF, {dcayAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + registry.add("dimuon/realign/dcax_MuonKine_GlobalMuonCuts_GoodMatches", "#mu^{+}#mu^{-} DCA_{x} (global muon cuts, good matches)", {HistType::kTHnSparseF, {dcaxAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + registry.add("dimuon/realign/dcay_MuonKine_GlobalMuonCuts_GoodMatches", "#mu^{+}#mu^{-} DCA_{y} (global muon cuts, good matches)", {HistType::kTHnSparseF, {dcayAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + registry.add("dimuon/realign/dcax_ScaledMftKine_GlobalMuonCuts_GoodMatches", "#mu^{+}#mu^{-} DCA_{x} (global muon cuts, rescaled MFT, good matches)", {HistType::kTHnSparseF, {dcaMftxAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + registry.add("dimuon/realign/dcay_ScaledMftKine_GlobalMuonCuts_GoodMatches", "#mu^{+}#mu^{-} DCA_{y} (global muon cuts, rescaled MFT, good matches)", {HistType::kTHnSparseF, {dcaMftyAxis, pAxis, pTAxis, muPosQuadrantAxis, muNegQuadrantAxis}}); + } } int GetDetElemId(int iDetElemNumber) @@ -747,8 +749,7 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc template int GetQuadrant(const T& track) { - // double phi = track.phi() * 180 / TMath::Pi(); - return GetQuadrant(track.phi()); + return GetQuadrant(static_cast(track.phi())); } template @@ -761,8 +762,6 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc // MCH track format. // Parameter conversion - // double alpha1, alpha3, alpha4, x2, x3, x4; - double x2 = fwdtrack.getPhi(); double x3 = fwdtrack.getTanl(); double x4 = fwdtrack.getInvQPt(); @@ -835,8 +834,6 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc o2::dataformats::GlobalFwdTrack convertedTrack; // Parameter conversion - // double alpha1, alpha3, alpha4, x2, x3, x4; - double alpha1 = mchParam.getNonBendingSlope(); double alpha3 = mchParam.getBendingSlope(); double alpha4 = mchParam.getInverseBendingMomentum(); @@ -1072,8 +1069,35 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc return true; } + template + bool isGoodGlobalMatching(const TMUON& muon, + double matchChi2Cut) + { + if (static_cast(muon.trackType()) > GlobalTrackTypeMax) + return false; + + // MFT-MCH match chi2 cut + if (muon.chi2MatchMCHMFT() > matchChi2Cut) + return false; + + return true; + } + + template + o2::track::TrackParCovFwd TrackToParCovFwd(const T& track) + { + double chi2 = track.chi2(); + SMatrix5 tpars(track.x(), track.y(), track.phi(), track.tgl(), track.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 trackparCov{track.z(), tpars, tcovs, chi2}; + return trackparCov; + } + template - o2::dataformats::GlobalFwdTrack FwdToTrackPar(const T& track) + o2::dataformats::GlobalFwdTrack TrackToGlobalFwd(const T& track) { double chi2 = track.chi2(); SMatrix5 tpars(track.x(), track.y(), track.phi(), track.tgl(), track.signed1Pt()); @@ -1112,26 +1136,9 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc track.setBendingSlope(ySlope + ySlopeCorrection); } - void TransformMFT(o2::dataformats::GlobalFwdTrack& track) - { - auto mchTrack = FwdtoMCH(track); - - TransformMFTPar(mchTrack); - - auto transformedTrack = MCHtoFwd(mchTrack); - track.setParameters(transformedTrack.getParameters()); - track.setZ(transformedTrack.getZ()); - track.setCovariances(transformedTrack.getCovariances()); - } - void TransformMFT(o2::track::TrackParCovFwd& fwdtrack) { - o2::dataformats::GlobalFwdTrack track; - track.setParameters(fwdtrack.getParameters()); - track.setZ(fwdtrack.getZ()); - track.setCovariances(fwdtrack.getCovariances()); - - auto mchTrack = FwdtoMCH(track); + auto mchTrack = FwdtoMCH(fwdtrack); TransformMFTPar(mchTrack); @@ -1264,6 +1271,35 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc return PropagateMCH(track, z); } + o2::dataformats::GlobalFwdTrack PropagateMCHToVertex(const o2::track::TrackParCovFwd& muon, + const double vx, const double vy, const double vz, + const double covVx, const double covVy) + { + auto mchTrack = FwdtoMCH(muon); + + o2::mch::TrackExtrap::extrapToVertex(mchTrack, vx, vy, vz, covVx, covVy); + + auto proptrack = MCHtoFwd(mchTrack); + o2::dataformats::GlobalFwdTrack propmuon; + propmuon.setParameters(proptrack.getParameters()); + propmuon.setZ(proptrack.getZ()); + propmuon.setCovariances(proptrack.getCovariances()); + + return propmuon; + } + + template + o2::dataformats::GlobalFwdTrack PropagateMCHToVertex(const o2::track::TrackParCovFwd& muon, + const C& collision) + { + return PropagateMCHToVertex(muon, + collision.posX(), + collision.posY(), + collision.posZ(), + collision.covXX(), + collision.covYY()); + } + template o2::dataformats::GlobalFwdTrack PropagateMFT(const TMFT& mftTrack, float z) { @@ -1299,27 +1335,12 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc return propmuon; } - template - o2::dataformats::GlobalFwdTrack PropagateMFTToDCA(const TMFT& mftTrack, const C& collision, float zshift) + template + o2::dataformats::GlobalFwdTrack PropagateMFTToDCA(o2::track::TrackParCovFwd mftTrack, + const C& collision, + float zshift) { // static double Bz = -10001; - double chi2 = mftTrack.chi2(); - double phiCorrDeg = 0; - double phiCorr = phiCorrDeg * o2::constants::math::Deg2Rad; - 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.cfgEnableMFTAlignmentCorrections) { - TransformMFT(fwdtrack); - } - o2::dataformats::GlobalFwdTrack propmuon; // double propVec[3] = {}; // propVec[0] = collision.posX() - mftTrack.x(); @@ -1334,40 +1355,28 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc // o2::field::MagneticField* field = static_cast(TGeoGlobalMagField::Instance()->GetField()); // Bz = field->getBz(centerZ); // } - fwdtrack.propagateToZ(collision.posZ() - zshift, mBzAtMftCenter); + mftTrack.propagateToZ(collision.posZ() - zshift, mBzAtMftCenter); - propmuon.setParameters(fwdtrack.getParameters()); - propmuon.setZ(fwdtrack.getZ()); - propmuon.setCovariances(fwdtrack.getCovariances()); + o2::dataformats::GlobalFwdTrack result; + result.setParameters(mftTrack.getParameters()); + result.setZ(mftTrack.getZ()); + result.setCovariances(mftTrack.getCovariances()); - return propmuon; + return result; } - template - o2::dataformats::GlobalFwdTrack PropagateMFTToDCA(const TMFT& mftTrack, const TMUON& mchTrack, const C& collision, float zshift) + template + o2::dataformats::GlobalFwdTrack PropagateMFTToDCA(o2::track::TrackParCovFwd mftTrack, + const o2::track::TrackParCovFwd& mchTrack, + const C& collision, + float zshift) { // static double Bz = -10001; - double chi2 = mftTrack.chi2(); - double phiCorrDeg = 0; - double phiCorr = phiCorrDeg * o2::constants::math::Deg2Rad; - 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.cfgEnableMFTAlignmentCorrections) { - TransformMFT(fwdtrack); - } // extrapolation with MCH tools - auto mchTrackAtMFT = FwdtoMCH(FwdToTrackPar(mchTrack)); - o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchTrackAtMFT, mftTrack.z()); - UpdateTrackMomentum(fwdtrack, mchTrackAtMFT); + auto mchTrackAtMFT = FwdtoMCH(mchTrack); + o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchTrackAtMFT, mftTrack.getZ()); + UpdateTrackMomentum(mftTrack, mchTrackAtMFT); // double propVec[3] = {}; // propVec[0] = collision.posX() - mftTrack.x(); @@ -1382,27 +1391,24 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc // o2::field::MagneticField* field = static_cast(TGeoGlobalMagField::Instance()->GetField()); // Bz = field->getBz(centerZ); // } - fwdtrack.propagateToZ(collision.posZ() - zshift, mBzAtMftCenter); + mftTrack.propagateToZ(collision.posZ() - zshift, mBzAtMftCenter); - o2::dataformats::GlobalFwdTrack propmuon; - propmuon.setParameters(fwdtrack.getParameters()); - propmuon.setZ(fwdtrack.getZ()); - propmuon.setCovariances(fwdtrack.getCovariances()); + o2::dataformats::GlobalFwdTrack result; + result.setParameters(mftTrack.getParameters()); + result.setZ(mftTrack.getZ()); + result.setCovariances(mftTrack.getCovariances()); - return propmuon; + return result; } - template - o2::dataformats::GlobalFwdTrack PropagateMFTtoMCH(const TMFT& mftTrack, const o2::mch::TrackParam& mchTrackPar, const double z) + o2::dataformats::GlobalFwdTrack PropagateMFTtoMCH(o2::track::TrackParCovFwd mftTrackPar, + const o2::mch::TrackParam& mchTrackPar, + const double z) { // extrapolation with MCH tools auto mchTrackAtMFT = mchTrackPar; - o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchTrackAtMFT, mftTrack.z()); + o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchTrackAtMFT, mftTrackPar.getZ()); - auto mftTrackPar = FwdToTrackPar(mftTrack); - if (configMFTAlignmentCorrections.cfgEnableMFTAlignmentCorrections) { - TransformMFT(mftTrackPar); - } auto mftTrackProp = FwdtoMCH(mftTrackPar); UpdateTrackMomentum(mftTrackProp, mchTrackAtMFT); if (z < AbsorberBackZ) { @@ -1431,76 +1437,378 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc return MCHtoFwd(mftTrackProp); } - void FillDCAPlots(MyEvents const& collisions, - MyBCs const& bcs, - MyMuonsWithCov const& muonTracks, - MyMFTs const& mftTracks, - const std::map& collisionInfos) + template + o2::dataformats::GlobalFwdTrack PropagateMFTToVertex(const o2::track::TrackParCovFwd& mftTrackPar, + const o2::track::TrackParCovFwd& mchTrackPar, + const C& collision) { - // outer loop over collisions - for (const auto& [collisionIndex, collisionInfo] : collisionInfos) { - auto const& collision = collisions.rawIteratorAt(collisionIndex); - const auto& bc = bcs.rawIteratorAt(collision.bcId()); + // extrapolation with MCH tools + auto mchTrackAtMFT = FwdtoMCH(mchTrackPar); + o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchTrackAtMFT, mftTrackPar.getZ()); - // remove TF/ROF borders and ambiguous collisions - if (!bc.selection_bit(o2::aod::evsel::kNoTimeFrameBorder) || - !bc.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) { - continue; - } + auto mftTrackProp = FwdtoMCH(mftTrackPar); - registry.get(HIST("vertex_y_vs_x"))->Fill(collision.posX(), collision.posY()); - registry.get(HIST("vertex_z"))->Fill(collision.posZ()); + // update global track momentum from the MCH track + double pRatio = mftTrackProp.p() / mchTrackAtMFT.p(); + double newInvBendMom = mftTrackProp.getInverseBendingMomentum() * pRatio; + mftTrackProp.setInverseBendingMomentum(newInvBendMom); + mftTrackProp.setCharge(mchTrackAtMFT.getCharge()); - if (cfgEnableVertexShiftAnalysis || cfgEnableMftDcaAnalysis) { - registry.get(HIST("DCA/MFT/nTracksMFT"))->Fill(collisionInfo.mftTracks.size()); - } + o2::mch::TrackExtrap::extrapToVertex(mftTrackProp, + collision.posX(), + collision.posY(), + collision.posZ(), + collision.covXX(), + collision.covYY()); - if (cfgEnableVertexShiftAnalysis || cfgEnableMftDcaAnalysis) { - // loop over MFT tracks - auto mftTrackIds = collisionInfo.mftTracks; - if (cfgMftTracksMultiplicityMax > 0 && mftTrackIds.size() > cfgMftTracksMultiplicityMax) { - auto rng = std::default_random_engine{}; - std::shuffle(std::begin(mftTrackIds), std::end(mftTrackIds), rng); - mftTrackIds.resize(cfgMftTracksMultiplicityMax); - } + return MCHtoFwd(mftTrackProp); + } - for (const auto& mftIndex : mftTrackIds) { - auto const& mftTrack = mftTracks.rawIteratorAt(mftIndex); + void getMuonPairs(const CollisionInfo& collisionInfo, + std::vector& muonPairs) + { + // outer loop over muon tracks + for (const auto& mchIndex1 : collisionInfo.mchTracks) { + // inner loop over muon tracks + for (const auto& mchIndex2 : collisionInfo.mchTracks) { + // avoid double-counting of muon pairs + if (mchIndex2 <= mchIndex1) { + continue; + } - if (mftTrack.isCA()) { - continue; - } + muonPairs.emplace_back(mchIndex1, mchIndex2); + } + } + } - bool isGoodMFT = IsGoodMFT(mftTrack, 999.f, 5); - if (!isGoodMFT) { - continue; - } + ROOT::Math::PxPyPzMVector getMuMu4Momentum(const o2::dataformats::GlobalFwdTrack& track1, const o2::dataformats::GlobalFwdTrack& track2) + { + ROOT::Math::PxPyPzMVector muon1{ + track1.getPx(), + track1.getPy(), + track1.getPz(), + o2::constants::physics::MassMuon}; + + ROOT::Math::PxPyPzMVector muon2{ + track2.getPx(), + track2.getPy(), + track2.getPz(), + o2::constants::physics::MassMuon}; + + return muon1 + muon2; + } - auto mftTrackAtDCA = PropagateMFTToDCA(mftTrack, collision, cfgVertexZshift); - double dcax = mftTrackAtDCA.getX() - collision.posX(); - double dcay = mftTrackAtDCA.getY() - collision.posY(); - double phi = mftTrack.phi() * o2::constants::math::Rad2Deg; - int mftNclusters = mftTrack.nClusters(); - double chi2NDF = static_cast(mftNclusters) * 2 - 5; + double getMuMuAngle(const o2::dataformats::GlobalFwdTrack& track1, const o2::dataformats::GlobalFwdTrack& track2) + { + ROOT::Math::XYZVector muon1{ + track1.getPx(), + track1.getPy(), + track1.getPz()}; - const int nMftLayers = 10; - std::array firedLayers{false}; - for (int layer = 0; layer < nMftLayers; layer++) { - if (((mftTrack.mftClusterSizesAndTrackFlags() >> (layer * 6)) & 0x3F) != 0) { - firedLayers[layer] = true; - } else { - firedLayers[layer] = false; - } - } + ROOT::Math::XYZVector muon2{ + track2.getPx(), + track2.getPy(), + track2.getPz()}; - if (cfgEnableMftDcaAnalysis) { - if (cfgEnableMftDcaExtraPlots) { - 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); - } - } - } + return std::acos(muon1.Unit().Dot(muon2.Unit())); + } + + double getMuMuInvariantMass(const o2::dataformats::GlobalFwdTrack& track1, const o2::dataformats::GlobalFwdTrack& track2) + { + return getMuMu4Momentum(track1, track2).M(); + } + + template + bool MchRealignTrack(const TMUON& mchTrack, const TCLUS& clusters, TrackRealigned& convertedTrack, bool applyCorrections) + { + auto mchTrackPar = FwdtoMCH(TrackToGlobalFwd(mchTrack)); + + // loop over attached clusters + int clIndex = -1; + auto clustersSliced = clusters.sliceBy(perMuon, mchTrack.globalIndex()); // Slice clusters by muon id + for (auto const& cluster : clustersSliced) { + clIndex += 1; + + int deId = cluster.deId(); + int chamber = deId / 100 - 1; + if (chamber < 0 || chamber >= NMchChambers) { + continue; + } + + math_utils::Point3D local; + math_utils::Point3D master; + + master.SetXYZ(cluster.x(), cluster.y(), cluster.z()); + + if (configRealign.cfgEnableMCHRealign) { + // Transformation from reference geometry frame to new geometry frame + transformRef[cluster.deId()].MasterToLocal(master, local); + transformNew[cluster.deId()].LocalToMaster(local, master); + } + + // shift the clusters to correct the longitudinal shift of the dipole + if (cfgDipoleZshift.value != 0) { + master.SetZ(master.z() + cfgDipoleZshift.value); + } + + if (applyCorrections) { + auto correctionsIt = mMchAlignmentCorrections.find(cluster.deId()); + if (correctionsIt != mMchAlignmentCorrections.end()) { + const auto& corrections = correctionsIt->second; + master.SetX(master.x() + corrections.x); + master.SetY(master.y() + corrections.y); + master.SetZ(master.z() + corrections.z); + } + } + + // realigned MCH cluster + auto clusterMCH = new mch::Cluster(); + clusterMCH->x = master.x(); + clusterMCH->y = master.y(); + clusterMCH->z = master.z(); + + uint32_t ClUId = mch::Cluster::buildUniqueId(static_cast(cluster.deId() / 100) - 1, cluster.deId(), clIndex); + clusterMCH->uid = ClUId; + clusterMCH->ex = cluster.isGoodX() ? 0.2 : 10.0; + clusterMCH->ey = cluster.isGoodY() ? 0.2 : 10.0; + + // Add transformed cluster into temporary variable + convertedTrack.createParamAtCluster(*clusterMCH); + } + + bool removable{false}; + // Refit the re-aligned track + if (convertedTrack.getNClusters() != 0) { + removable = RemoveTrack(convertedTrack); + } else { + LOGF(fatal, "Muon track %d has no associated clusters.", mchTrack.globalIndex()); + } + + // subtract the longitudinal shift of the dipole from the track z + if (cfgDipoleZshift.value != 0) { + auto& trackParam = *(convertedTrack.begin()); + trackParam.setZ(trackParam.getZ() - cfgDipoleZshift.value); + } + + return !removable; + } + + template + void InitCollisions(COLL const& collisions, + BC const& bcs, + TMUON const& muonTracks, + aod::FwdTrkCls const& clusters, + std::map& collisionInfos) + { + mMchTrackPars.clear(); + mMchTrackParsNew.clear(); + + // fill collision information for global muon tracks (MFT-MCH-MID matches) + for (const auto& muonTrack : muonTracks) { + if (!muonTrack.has_collision()) { + continue; + } + + auto collision = collisions.rawIteratorAt(muonTrack.collisionId()); + + if (cfgRequireGoodRCT && !rctChecker(collision)) { + continue; + } + + uint64_t collisionIndex = collision.globalIndex(); + + auto bc = bcs.rawIteratorAt(collision.bcId()); + + auto& collisionInfo = collisionInfos[collisionIndex]; + collisionInfo.bc = bc.globalBC(); + collisionInfo.zVertex = collision.posZ(); + + if (static_cast(muonTrack.trackType()) > GlobalTrackTypeMax) { + // standalone MCH or MCH-MID tracks + uint64_t mchTrackIndex = muonTrack.globalIndex(); + collisionInfo.mchTracks.push_back(mchTrackIndex); + + // initialize the original MCH track parameters + mMchTrackPars.try_emplace(mchTrackIndex, TrackParExt(fwdtrackutils::getTrackParCovFwd(muonTrack, muonTrack), muonTrack.nClusters())); + + // refit MCH track if requested + if (configRealign.cfgEnableMCHRefit || configRealign.cfgEnableMCHRealign) { + TrackRealigned convertedTrack; + bool convertedTrackOk = MchRealignTrack(muonTrack, clusters, convertedTrack, !mMchAlignmentCorrections.empty()); + + // Get the re-aligned track parameters: track param at the first cluster + mch::TrackParam trackParam = mch::TrackParam(convertedTrack.first()); + + auto mchTrackParIt = mMchTrackParsNew.try_emplace(mchTrackIndex, TrackParExt(MCHtoFwd(trackParam), convertedTrack.getNClusters())); + if (mchTrackParIt.second) { + // the insertion succeeded + mchTrackParIt.first->second.setTrackChi2(trackParam.getTrackChi2() / convertedTrack.getNDF()); + if (!convertedTrackOk) { + mchTrackParIt.first->second.setRemovable(); + } + } + } else { + // initialize the new MCH track parameters with the original ones, without refitting + mMchTrackParsNew.try_emplace(mchTrackIndex, TrackParExt(fwdtrackutils::getTrackParCovFwd(muonTrack, muonTrack), muonTrack.nClusters())); + } + } else { + // global muon tracks (MFT-MCH or MFT-MCH-MID) + uint64_t muonTrackIndex = muonTrack.globalIndex(); + auto const& mchTrack = muonTrack.template matchMCHTrack_as(); + uint64_t mchTrackIndex = mchTrack.globalIndex(); + + // check if a vector of global muon candidates is already available for the current MCH index + // if not, initialize a new one and add the current global muon track + // bool globalMuonTrackFound = false; + auto matchingCandidateIterator = collisionInfo.globalMuonTracks.find(mchTrackIndex); + if (matchingCandidateIterator != collisionInfo.globalMuonTracks.end()) { + matchingCandidateIterator->second.push_back(muonTrackIndex); + // globalMuonTrackFound = true; + } else { + collisionInfo.globalMuonTracks[mchTrackIndex].push_back(muonTrackIndex); + } + } + } + + // sort the vectors of matching candidates in ascending order based on the matching chi2 value + auto compareChi2 = [&muonTracks](uint64_t trackIndex1, uint64_t trackIndex2) -> bool { + auto const& track1 = muonTracks.rawIteratorAt(trackIndex1); + auto const& track2 = muonTracks.rawIteratorAt(trackIndex2); + + return (track1.chi2MatchMCHMFT() < track2.chi2MatchMCHMFT()); + }; + + for (auto& [collisionIndex, collisionInfo] : collisionInfos) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) + for (auto& [mchIndex, globalTracksVector] : collisionInfo.globalMuonTracks) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) + std::sort(globalTracksVector.begin(), globalTracksVector.end(), compareChi2); + } + } + } + + void InitCollisions(MyEvents const& collisions, + MyBCs const& bcs, + MyMuonsWithCov const& muonTracks, + aod::FwdTrkCls const& clusters, + MyMFTs const& mftTracks, + std::map& collisionInfos) + { + InitCollisions(collisions, bcs, muonTracks, clusters, collisionInfos); + + mMftTrackPars.clear(); + mMftTrackParsNew.clear(); + + // fill collision information for MFT standalone tracks + for (const auto& mftTrack : mftTracks) { + if (!mftTrack.has_collision()) { + continue; + } + + auto collision = collisions.rawIteratorAt(mftTrack.collisionId()); + uint64_t collisionIndex = collision.globalIndex(); + + auto bc = bcs.rawIteratorAt(collision.bcId()); + + uint64_t mftTrackIndex = mftTrack.globalIndex(); + + auto& collisionInfo = collisionInfos[collisionIndex]; + collisionInfo.bc = bc.globalBC(); + collisionInfo.zVertex = collision.posZ(); + + collisionInfo.mftTracks.push_back(mftTrackIndex); + + // initialize the original MFT track parameters + auto mftTrackFwd = TrackToParCovFwd(mftTrack); + mMftTrackPars.try_emplace(mftTrackIndex, TrackParExt(mftTrackFwd, mftTrack.nClusters())); + + // initialize the corrected MFT track parameters, if requested + if (configMFTAlignmentCorrections.cfgEnableMFTAlignmentCorrections) { + TransformMFT(mftTrackFwd); + mMftTrackParsNew.try_emplace(mftTrackIndex, TrackParExt(mftTrackFwd, mftTrack.nClusters())); + } else { + // initialize the new MFT track parameters with the original ones, without corrections + mMftTrackParsNew.try_emplace(mftTrackIndex, TrackParExt(mftTrackFwd, mftTrack.nClusters())); + } + } + } + + void FillMftPlots(MyEvents const& collisions, + MyBCs const& bcs, + MyMuonsWithCov const& muonTracks, + MyMFTs const& mftTracks, + const std::map& collisionInfos) + { + // outer loop over collisions + for (const auto& [collisionIndex, collisionInfo] : collisionInfos) { + auto const& collision = collisions.rawIteratorAt(collisionIndex); + const auto& bc = bcs.rawIteratorAt(collision.bcId()); + + // remove TF/ROF borders and ambiguous collisions + if (!bc.selection_bit(o2::aod::evsel::kNoTimeFrameBorder) || + !bc.selection_bit(o2::aod::evsel::kNoITSROFrameBorder)) { + continue; + } + + registry.get(HIST("vertex_y_vs_x"))->Fill(collision.posX(), collision.posY()); + registry.get(HIST("vertex_z"))->Fill(collision.posZ()); + + if (cfgEnableVertexShiftAnalysis || cfgEnableMftDcaAnalysis) { + registry.get(HIST("DCA/MFT/nTracksMFT"))->Fill(collisionInfo.mftTracks.size()); + } + + if (cfgEnableVertexShiftAnalysis || cfgEnableMftDcaAnalysis) { + // loop over MFT tracks + auto mftTrackIds = collisionInfo.mftTracks; + if (cfgMftTracksMultiplicityMax > 0 && mftTrackIds.size() > cfgMftTracksMultiplicityMax) { + auto rng = std::default_random_engine{}; + std::shuffle(std::begin(mftTrackIds), std::end(mftTrackIds), rng); + mftTrackIds.resize(cfgMftTracksMultiplicityMax); + } + + for (const auto& mftIndex : mftTrackIds) { + auto const& mftTrack = mftTracks.rawIteratorAt(mftIndex); + + if (mftTrack.isCA()) { + continue; + } + + bool isGoodMFT = IsGoodMFT(mftTrack, 999.f, 5); + if (!isGoodMFT) { + continue; + } + + // get the pre-stored MFT track parameters after corrections + // if MFT corrections are not enabled, the original MFT track parameters are retrieved + const auto mftTrackParIt = mMftTrackParsNew.find(mftIndex); + if (mftTrackParIt == mMftTrackParsNew.end()) { + continue; + } + const auto& mftTrackPar = mftTrackParIt->second; + + auto mftTrackAtDCA = PropagateMFTToDCA(mftTrackPar, collision, cfgVertexZshift); + double dcax = mftTrackAtDCA.getX() - collision.posX(); + double dcay = mftTrackAtDCA.getY() - collision.posY(); + double phi = mftTrack.phi() * o2::constants::math::Rad2Deg; + int mftNclusters = mftTrack.nClusters(); + double chi2NDF = static_cast(mftNclusters) * 2 - 5; + + const int nMftLayers = 10; + std::array firedLayers{false}; + for (int layer = 0; layer < nMftLayers; layer++) { + if (((mftTrack.mftClusterSizesAndTrackFlags() >> (layer * 6)) & 0x3F) != 0) { + firedLayers[layer] = true; + } else { + firedLayers[layer] = false; + } + } + + if (cfgEnableMftDcaAnalysis) { + if (cfgEnableMftDcaExtraPlots) { + 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); + } + } + } if (mftTrack.chi2() <= cfgTrackChi2MftUp) { registry.get(HIST("DCA/MFT/DCA_y_vs_x"))->Fill(dcax, dcay); @@ -1544,7 +1852,7 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc -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}; for (int zi = 0; zi < nPoints; zi++) { - auto mftTrackAtDCAshifted = PropagateMFTToDCA(mftTrack, collision, zshift[zi] / 10.f); + auto mftTrackAtDCAshifted = PropagateMFTToDCA(mftTrackPar, collision, zshift[zi] / 10.f); double dcaxShifted = mftTrackAtDCAshifted.getX() - collision.posX(); double dcayShifted = mftTrackAtDCAshifted.getY() - collision.posY(); registry.get(HIST("DCA/MFT/DCA_x_vs_phi_vs_zshift"))->Fill(zshift[zi], phi, dcaxShifted); @@ -1570,6 +1878,9 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc const auto& mchTrack = muonTrack.template matchMCHTrack_as(); const auto& mftTrack = muonTrack.template matchMFTTrack_as(); + auto mchIndex = mchTrack.globalIndex(); + auto mftIndex = mftTrack.globalIndex(); + if (muonTrack.chi2MatchMCHMFT() < cfgMftDcaMatchChi2Up.value) { continue; } @@ -1592,12 +1903,27 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc continue; } - bool isGoodMuon = IsGoodMuon(mchTrack, collision, cfgTrackChi2MchUp, 0.f, cfgPtMchLow, {cfgEtaMftlow, cfgEtaMftup}, {cfgRabsLow, cfgRabsUp}, fSigmaPdcaUp); + bool isGoodMuon = IsGoodMuon(mchTrack, collision, cfgTrackChi2MchUp, 0.f, cfgPtMchLow, {cfgEtaMftLow, cfgEtaMftUp}, {cfgRabsLow, cfgRabsUp}, fSigmaPdcaUp); if (!isGoodMuon) { continue; } - auto mftTrackAtDCA = PropagateMFTToDCA(mftTrack, mchTrack, collision, cfgVertexZshift); + // get the pre-stored MFT track parameters after corrections + // if MFT corrections are not enabled, the original MFT track parameters are retrieved + const auto mftTrackParIt = mMftTrackParsNew.find(mftIndex); + if (mftTrackParIt == mMftTrackParsNew.end()) { + continue; + } + const auto& mftTrackPar = mftTrackParIt->second; + + // get the pre-stored MCH track parameters + const auto mchTrackParIt = mMchTrackPars.find(mchIndex); + if (mchTrackParIt == mMchTrackPars.end()) { + continue; + } + const auto& mchTrackPar = mchTrackParIt->second; + + auto mftTrackAtDCA = PropagateMFTToDCA(mftTrackPar, mchTrackPar, collision, cfgVertexZshift); double dcax = mftTrackAtDCA.getX() - collision.posX(); double dcay = mftTrackAtDCA.getY() - collision.posY(); int mftNclusters = mftTrack.nClusters(); @@ -1611,84 +1937,11 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc } } - template - bool MchRealignTrack(const TMUON& mchTrack, const TCLUS& clusters, TrackRealigned& convertedTrack, bool applyCorrections) - { - // loop over attached clusters - int clIndex = -1; - auto clustersSliced = clusters.sliceBy(perMuon, mchTrack.globalIndex()); // Slice clusters by muon id - for (auto const& cluster : clustersSliced) { - clIndex += 1; - - int deId = cluster.deId(); - int chamber = deId / 100 - 1; - if (chamber < 0 || chamber >= NMchChambers) { - continue; - } - - math_utils::Point3D local; - math_utils::Point3D master; - - master.SetXYZ(cluster.x(), cluster.y(), cluster.z()); - - if (configRealign.cfgEnableMCHRealign) { - // Transformation from reference geometry frame to new geometry frame - transformRef[cluster.deId()].MasterToLocal(master, local); - transformNew[cluster.deId()].LocalToMaster(local, master); - } - - // shift the clusters to correct the longitudinal shift of the dipole - if (cfgDipoleZshift.value != 0) { - master.SetZ(master.z() + cfgDipoleZshift.value); - } - - if (applyCorrections) { - auto correctionsIt = mMchAlignmentCorrections.find(cluster.deId()); - if (correctionsIt != mMchAlignmentCorrections.end()) { - const auto& corrections = correctionsIt->second; - master.SetX(master.x() + corrections.x); - master.SetY(master.y() + corrections.y); - master.SetZ(master.z() + corrections.z); - } - } - - // realigned MCH cluster - auto clusterMCH = new mch::Cluster(); - clusterMCH->x = master.x(); - clusterMCH->y = master.y(); - clusterMCH->z = master.z(); - - uint32_t ClUId = mch::Cluster::buildUniqueId(static_cast(cluster.deId() / 100) - 1, cluster.deId(), clIndex); - clusterMCH->uid = ClUId; - clusterMCH->ex = cluster.isGoodX() ? 0.2 : 10.0; - clusterMCH->ey = cluster.isGoodY() ? 0.2 : 10.0; - - // Add transformed cluster into temporary variable - convertedTrack.createParamAtCluster(*clusterMCH); - } - - bool removable{false}; - // Refit the re-aligned track - if (convertedTrack.getNClusters() != 0) { - removable = RemoveTrack(convertedTrack); - } else { - LOGF(fatal, "Muon track %d has no associated clusters.", mchTrack.globalIndex()); - } - - // subtract the longitudinal shift of the dipole from the track z - if (cfgDipoleZshift.value != 0) { - auto& trackParam = *(convertedTrack.begin()); - trackParam.setZ(trackParam.getZ() - cfgDipoleZshift.value); - } - - return !removable; - } - - void FillResidualsPlots(MyEvents const& collisions, - MyBCs const& bcs, - MyMuonsWithCov const& muonTracks, - aod::FwdTrkCls const& clusters, - const std::map& collisionInfos) + void FillMchPlots(MyEvents const& collisions, + MyBCs const& bcs, + MyMuonsWithCov const& muonTracks, + aod::FwdTrkCls const& clusters, + const std::map& collisionInfos) { if (!cfgEnableMftMchResidualsAnalysis && !cfgEnableMftMchMatchingAnalysis) { return; @@ -1714,7 +1967,10 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc int quadrant = GetQuadrant(mftTrack); int posNeg = (mchTrack.sign() >= 0) ? 0 : 1; - bool isGoodMuon = IsGoodMuon(mchTrack, collision, cfgTrackChi2MchUp, cfgMftMchResidualsPLow, cfgMftMchResidualsPtLow, {cfgEtaMftlow, cfgEtaMftup}, {cfgRabsLow, cfgRabsUp}, fSigmaPdcaUp); + auto mchIndex = mchTrack.globalIndex(); + auto mftIndex = mftTrack.globalIndex(); + + bool isGoodMuon = IsGoodMuon(mchTrack, collision, cfgTrackChi2MchUp, cfgMftMchResidualsPLow, cfgMftMchResidualsPtLow, {cfgEtaMftLow, cfgEtaMftUp}, {cfgRabsLow, cfgRabsUp}, fSigmaPdcaUp); if (!isGoodMuon) { continue; } @@ -1724,26 +1980,39 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc continue; } - double matchChi2 = muonTrack.chi2MatchMCHMFT(); - if (matchChi2 > cfgMftDcaMatchChi2Up.value) { + // get the pre-stored MFT track parameters + const auto mftTrackParIt = mMftTrackPars.find(mftIndex); + if (mftTrackParIt == mMftTrackPars.end()) { + continue; + } + const auto& mftTrackPar = mftTrackParIt->second; + + // get the pre-stored MFT track parameters after corrections + // if MFT corrections are not enabled, the original MFT track parameters are retrieved + const auto mftTrackParNewIt = mMftTrackParsNew.find(mftIndex); + if (mftTrackParNewIt == mMftTrackParsNew.end()) { continue; } + const auto& mftTrackParNew = mftTrackParNewIt->second; - // refit MCH track if enabled - TrackRealigned convertedTrack; - bool convertedTrackOk = false; - if (configRealign.cfgEnableMCHRealign) { - convertedTrackOk = MchRealignTrack(mchTrack, clusters, convertedTrack, false); + // get the pre-stored MCH track parameters + const auto mchTrackParIt = mMchTrackPars.find(mchIndex); + if (mchTrackParIt == mMchTrackPars.end()) { + continue; } + const auto& mchTrackPar = mchTrackParIt->second; - // apply alignment corrections if available - TrackRealigned convertedTrackWithCorr; - bool convertedTrackWithCorrOk = false; - if (!mMchAlignmentCorrections.empty()) { - convertedTrackWithCorrOk = MchRealignTrack(mchTrack, clusters, convertedTrackWithCorr, true); + // get the pre-stored MCH track parameters after refit + const auto mchTrackParNewIt = mMchTrackParsNew.find(mchIndex); + if (mchTrackParNewIt == mMchTrackParsNew.end()) { + continue; } + const auto& mchTrackParNew = mchTrackParNewIt->second; + + double matchChi2 = muonTrack.chi2MatchMCHMFT(); - if (cfgEnableMftMchResidualsAnalysis) { + // Residuals analysis between MFT tracks and MCH clusters + if (cfgEnableMftMchResidualsAnalysis && (matchChi2 <= cfgMftMchResidualsMatchChi2Up.value)) { // loop over attached clusters auto clustersSliced = clusters.sliceBy(perMuon, mchTrack.globalIndex()); // Slice clusters by muon id for (auto const& cluster : clustersSliced) { @@ -1755,18 +2024,17 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc int deIndex = getDEindex(deId); math_utils::Point3D local; - math_utils::Point3D master; - math_utils::Point3D masterWithCorr; + math_utils::Point3D master; // original cluster position + math_utils::Point3D masterRealign; // cluster position after realignment master.SetXYZ(cluster.x(), cluster.y(), cluster.z()); - masterWithCorr.SetXYZ(cluster.x(), cluster.y(), cluster.z()); + masterRealign.SetXYZ(cluster.x(), cluster.y(), cluster.z()); // apply realignment to MCH cluster if (configRealign.cfgEnableMCHRealign) { // Transformation from reference geometry frame to new geometry frame transformRef[cluster.deId()].MasterToLocal(master, local); - transformNew[cluster.deId()].LocalToMaster(local, master); - transformNew[cluster.deId()].LocalToMaster(local, masterWithCorr); + transformNew[cluster.deId()].LocalToMaster(local, masterRealign); } // apply alignment corrections to MCH cluster (if available) @@ -1774,78 +2042,72 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc auto correctionsIt = mMchAlignmentCorrections.find(cluster.deId()); if (correctionsIt != mMchAlignmentCorrections.end()) { const auto& corrections = correctionsIt->second; - masterWithCorr.SetX(masterWithCorr.x() + corrections.x); - masterWithCorr.SetY(masterWithCorr.y() + corrections.y); - masterWithCorr.SetZ(masterWithCorr.z() + corrections.z); + masterRealign.SetX(masterRealign.x() + corrections.x); + masterRealign.SetY(masterRealign.y() + corrections.y); + masterRealign.SetZ(masterRealign.z() + corrections.z); } } - // MFT-MCH residuals (MCH cluster is realigned if enabled) - // if the realignment is enabled and successful, the MFT track is extrpolated - // by taking the momentum from the MCH track refitted with the new alignment - if (!configRealign.cfgEnableMCHRealign || convertedTrackOk) { - auto mftTrackAtCluster = configRealign.cfgEnableMCHRealign ? PropagateMFTtoMCH(mftTrack, mch::TrackParam(convertedTrack.first()), master.z()) : PropagateMFTtoMCH(mftTrack, FwdtoMCH(FwdToTrackPar(mchTrack)), master.z()); - auto mftTrackParamAtCluster = FwdtoMCH(mftTrackAtCluster); + // MFT-MCH residuals from original alignment + const auto mftTrackAtCluster = PropagateMFTtoMCH(mftTrackPar, FwdtoMCH(mchTrackPar), master.z()); + const auto mftTrackParamAtCluster = FwdtoMCH(mftTrackAtCluster); - std::array xPos{master.x(), mftTrackAtCluster.getX()}; - std::array yPos{master.y(), mftTrackAtCluster.getY()}; + const std::array xPos{master.x(), mftTrackAtCluster.getX()}; + const std::array yPos{master.y(), mftTrackAtCluster.getY()}; - 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_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, 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()); - } + 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 // if the alignment corrections are available and the refitting is successful, the MFT track is extrpolated // by taking the momentum from the MCH track refitted with the alignment corrections and the new // alignment (if realignment is enabled) - if (convertedTrackWithCorrOk) { - auto mftTrackAtClusterWithCorr = PropagateMFTtoMCH(mftTrack, mch::TrackParam(convertedTrackWithCorr.first()), masterWithCorr.z()); - auto mftTrackParamAtClusterWithCorr = FwdtoMCH(mftTrackAtClusterWithCorr); + if (!mchTrackParNew.isRemovable()) { + const auto mftTrackAtClusterWithCorr = PropagateMFTtoMCH(mftTrackParNew, FwdtoMCH(mchTrackParNew), masterRealign.z()); + const auto mftTrackParamAtClusterWithCorr = FwdtoMCH(mftTrackAtClusterWithCorr); - std::array xPos{masterWithCorr.x(), mftTrackAtClusterWithCorr.getX()}; - std::array yPos{masterWithCorr.y(), mftTrackAtClusterWithCorr.getY()}; + const std::array xPosWithCorr{masterRealign.x(), mftTrackAtClusterWithCorr.getX()}; + const std::array yPosWithCorr{masterRealign.y(), mftTrackAtClusterWithCorr.getY()}; - 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_chamber_corr"))->Fill(chamber + 1, quadrant, posNeg, xPosWithCorr[0] - xPosWithCorr[1]); + registry.get(HIST("residuals/dy_vs_chamber_corr"))->Fill(chamber + 1, quadrant, posNeg, yPosWithCorr[0] - yPosWithCorr[1]); - 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()); + registry.get(HIST("residuals/dx_vs_de_corr"))->Fill(xPosWithCorr[0] - xPosWithCorr[1], deIndex, quadrant, posNeg, mchTrack.p(), mftTrackParamAtClusterWithCorr.getNonBendingSlope()); + registry.get(HIST("residuals/dy_vs_de_corr"))->Fill(yPosWithCorr[0] - yPosWithCorr[1], deIndex, quadrant, posNeg, mchTrack.p(), mftTrackParamAtClusterWithCorr.getBendingSlope()); } } - if (!configRealign.cfgEnableMCHRealign || convertedTrackOk) { - auto mchTrackAtDCA = configRealign.cfgEnableMCHRealign ? PropagateMCHRealigned(convertedTrack, collision.posZ()) : PropagateMCH(mchTrack, collision.posZ()); - auto dcax = mchTrackAtDCA.getX() - collision.posX(); - 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(), quadrant, posNeg, dcax); - registry.get(HIST("DCA/MCH/DCA_y_vs_sign_vs_quadrant_vs_mom"))->Fill(mchTrack.p(), quadrant, posNeg, dcay); - - if (cfgEnableMftMchResidualsExtraPlots) { - 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.cfgEnableMCHRealign ? 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()); - } + const auto mchTrackAtDCA = PropagateMCHParam(FwdtoMCH(mchTrackPar), collision.posZ()); + const auto dcax = mchTrackAtDCA.getX() - collision.posX(); + const 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(mchTrackPar.getP(), quadrant, posNeg, dcax); + registry.get(HIST("DCA/MCH/DCA_y_vs_sign_vs_quadrant_vs_mom"))->Fill(mchTrackPar.getP(), quadrant, posNeg, dcay); + + if (cfgEnableMftMchResidualsExtraPlots) { + 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); + const auto mchTrackAtMFT = PropagateMCHParam(FwdtoMCH(mchTrackPar), mftTrack.z()); + const double deltaPhi = mchTrackAtMFT.getPhi() - mftTrack.phi(); + registry.get(HIST("residuals/dphi_at_mft"))->Fill(deltaPhi, mftTrack.x(), mftTrack.y(), posNeg, mchTrackAtMFT.getP()); } - if (convertedTrackWithCorrOk) { - auto mchTrackAtDCA = PropagateMCHRealigned(convertedTrackWithCorr, collision.posZ()); - auto dcax = mchTrackAtDCA.getX() - collision.posX(); - auto dcay = mchTrackAtDCA.getY() - collision.posY(); + if (!mchTrackParNew.isRemovable()) { + const auto mchTrackAtDCAWithCorr = PropagateMCHParam(FwdtoMCH(mchTrackParNew), collision.posZ()); + const auto dcaxWithCorr = mchTrackAtDCAWithCorr.getX() - collision.posX(); + const auto dcayWithCorr = mchTrackAtDCAWithCorr.getY() - collision.posY(); - 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); + registry.get(HIST("DCA/MCH/DCA_x_vs_sign_vs_quadrant_vs_mom_corr"))->Fill(mchTrackParNew.getP(), quadrant, posNeg, dcaxWithCorr); + registry.get(HIST("DCA/MCH/DCA_y_vs_sign_vs_quadrant_vs_mom_corr"))->Fill(mchTrackParNew.getP(), quadrant, posNeg, dcayWithCorr); } } // MFT-MCH track residuals analysis - if (cfgEnableMftMchMatchingAnalysis && convertedTrackWithCorrOk) { + if (cfgEnableMftMchMatchingAnalysis && !mchTrackParNew.isRemovable() && (matchChi2 <= cfgMftMchResidualsMatchChi2Up.value)) { static constexpr int nRefPlanes = 2; const std::array refPlaneZ{cfgRefPlaneZMFT, cfgRefPlaneZMCH}; @@ -1856,26 +2118,215 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc std::array, 2> dphiPlots{registry.get(HIST("matching/dphiAtMFT")), registry.get(HIST("matching/dphiAtMCH"))}; for (int iRefPlane = 0; iRefPlane < nRefPlanes; iRefPlane++) { - const auto mftTrackAtRefPlane = configRealign.cfgEnableMCHRealign ? PropagateMFTtoMCH(mftTrack, mch::TrackParam(convertedTrackWithCorr.first()), refPlaneZ[iRefPlane]) : PropagateMFTtoMCH(mftTrack, FwdtoMCH(FwdToTrackPar(mchTrack)), refPlaneZ[iRefPlane]); - const auto mchTrackAtRefPlane = configRealign.cfgEnableMCHRealign ? PropagateMCHRealigned(convertedTrackWithCorr, refPlaneZ[iRefPlane]) : PropagateMCH(mchTrack, refPlaneZ[iRefPlane]); + const auto mftTrackAtRefPlane = PropagateMFTtoMCH(mftTrackParNew, FwdtoMCH(mchTrackParNew), refPlaneZ[iRefPlane]); + const auto mchTrackAtRefPlane = PropagateMCHParam(FwdtoMCH(mchTrackParNew), refPlaneZ[iRefPlane]); const auto& refTrackAtRefPlane = (iRefPlane == 0) ? mftTrackAtRefPlane : mchTrackAtRefPlane; auto dx = mchTrackAtRefPlane.getX() - mftTrackAtRefPlane.getX(); - dxPlots[iRefPlane]->Fill(dx, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrant, posNeg, mchTrack.p()); + dxPlots[iRefPlane]->Fill(dx, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrant, posNeg, mchTrackParNew.getP()); auto dy = mchTrackAtRefPlane.getY() - mftTrackAtRefPlane.getY(); - dyPlots[iRefPlane]->Fill(dy, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrant, posNeg, mchTrack.p()); + dyPlots[iRefPlane]->Fill(dy, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrant, posNeg, mchTrackParNew.getP()); - auto mftParamAtRefPlane = FwdtoMCH(mftTrackAtRefPlane); - auto mchParamAtRefPlane = FwdtoMCH(mchTrackAtRefPlane); + const auto mftParamAtRefPlane = FwdtoMCH(mftTrackAtRefPlane); + const auto mchParamAtRefPlane = FwdtoMCH(mchTrackAtRefPlane); auto dsx = mchParamAtRefPlane.getNonBendingSlope() - mftParamAtRefPlane.getNonBendingSlope(); - dsxPlots[iRefPlane]->Fill(dsx, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrant, posNeg, mchTrack.p()); + dsxPlots[iRefPlane]->Fill(dsx, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrant, posNeg, mchTrackParNew.getP()); auto dsy = mchParamAtRefPlane.getBendingSlope() - mftParamAtRefPlane.getBendingSlope(); - dsyPlots[iRefPlane]->Fill(dsy, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrant, posNeg, mchTrack.p()); + dsyPlots[iRefPlane]->Fill(dsy, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrant, posNeg, mchTrackParNew.getP()); auto dphi = RecoDecay::constrainAngle(mchTrackAtRefPlane.getPhi() - mftTrackAtRefPlane.getPhi(), -o2::constants::math::PI); - dphiPlots[iRefPlane]->Fill(dphi, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrant, posNeg, mchTrack.p()); + dphiPlots[iRefPlane]->Fill(dphi, refTrackAtRefPlane.getX(), refTrackAtRefPlane.getY(), quadrant, posNeg, mchTrackParNew.getP()); + } + } + } + } + } + +#define FILL_DIMUON_PLOT(trackPar1, trackPar2, trackPar1AtVertex, trackPar2AtVertex, histName) \ + { \ + auto mumu4mom = getMuMu4Momentum(trackPar1AtVertex, trackPar2AtVertex); \ + double p = mumu4mom.P(); \ + double pT = mumu4mom.Pt(); \ + double mass = mumu4mom.M(); \ + int quadrant1 = GetQuadrant(static_cast(std::atan2(trackPar1.getY(), trackPar1.getX()))); \ + int quadrant2 = GetQuadrant(static_cast(std::atan2(trackPar2.getY(), trackPar2.getX()))); \ + registry.get(HIST(histName))->Fill(mass, p, pT, quadrant1, quadrant2); \ + } + +#define FILL_DIMUON_DCA_PLOTS(trackPar1, trackPar2, trackPar1AtVertex, trackPar2AtVertex, trackPar1AtDca, trackPar2AtDca, histNameX, histNameY) \ + { \ + auto mumu4mom = getMuMu4Momentum(trackPar1AtVertex, trackPar2AtVertex); \ + double p = mumu4mom.P(); \ + double pT = mumu4mom.Pt(); \ + double dcax = trackPar1AtDca.getX() - trackPar2AtDca.getX(); \ + double dcay = trackPar1AtDca.getY() - trackPar2AtDca.getY(); \ + int quadrant1 = GetQuadrant(static_cast(std::atan2(trackPar1.getY(), trackPar1.getX()))); \ + int quadrant2 = GetQuadrant(static_cast(std::atan2(trackPar2.getY(), trackPar2.getX()))); \ + registry.get(HIST(histNameX))->Fill(dcax, p, pT, quadrant1, quadrant2); \ + registry.get(HIST(histNameY))->Fill(dcay, p, pT, quadrant1, quadrant2); \ + } + +#define FILL_DIMUON_ANGLE_PLOT(trackPar1, trackPar2, trackPar1AtVertex, trackPar2AtVertex, mchAngle, fwdAngle, histName) \ + { \ + auto mumu4mom = getMuMu4Momentum(trackPar1AtVertex, trackPar2AtVertex); \ + double p = mumu4mom.P(); \ + double pT = mumu4mom.Pt(); \ + double dAngle = mchAngle - fwdAngle; \ + int quadrant1 = GetQuadrant(static_cast(std::atan2(trackPar1.getY(), trackPar1.getX()))); \ + int quadrant2 = GetQuadrant(static_cast(std::atan2(trackPar2.getY(), trackPar2.getX()))); \ + registry.get(HIST(histName))->Fill(dAngle, fwdAngle, p, pT, quadrant1, quadrant2); \ + } + + void FillDimuonPlots(MyEvents const& collisions, + MyMuonsWithCov const& muonTracks, + const std::map& collisionInfos) + { + if (!cfgEnableDimuonAnalysis) { + return; + } + + for (const auto& [collisionIndex, collisionInfo] : collisionInfos) { + auto const& collision = collisions.rawIteratorAt(collisionIndex); + + std::vector muonPairs; + getMuonPairs(collisionInfo, muonPairs); + + for (const auto& [mchIndex1, mchIndex2] : muonPairs) { + + auto const& muonTrack1 = muonTracks.rawIteratorAt(mchIndex1); + auto const& muonTrack2 = muonTracks.rawIteratorAt(mchIndex2); + int sign1 = muonTrack1.sign(); + int sign2 = muonTrack2.sign(); + + // only consider opposite-sign pairs + if ((sign1 * sign2) >= 0) + continue; + + bool isGoodMuon1 = IsGoodMuon(muonTrack1, collision, cfgTrackChi2MchUp, 0.f, cfgPtMchLow, {cfgEtaMchLow, cfgEtaMchUp}, {cfgRabsLow, cfgRabsUp}, fSigmaPdcaUp); + bool isGoodMuon2 = IsGoodMuon(muonTrack2, collision, cfgTrackChi2MchUp, 0.f, cfgPtMchLow, {cfgEtaMchLow, cfgEtaMchUp}, {cfgRabsLow, cfgRabsUp}, fSigmaPdcaUp); + bool goodMuonTracks = (isGoodMuon1 && isGoodMuon2); + + if (!goodMuonTracks) { + continue; + } + + // get the pre-stored MCH track parameters + const auto mchTrackParIt1 = mMchTrackPars.find(mchIndex1); + if (mchTrackParIt1 == mMchTrackPars.end()) { + continue; + } + const auto mchTrackParIt2 = mMchTrackPars.find(mchIndex2); + if (mchTrackParIt2 == mMchTrackPars.end()) { + continue; + } + const auto& mchTrackPar1 = mchTrackParIt1->second; + auto mchTrackPar1AtDca = PropagateMCHParam(FwdtoMCH(mchTrackPar1), collision.posZ()); + auto mchTrackPar1AtVertex = PropagateMCHToVertex(mchTrackPar1, collision); + const auto& mchTrackPar2 = mchTrackParIt2->second; + auto mchTrackPar2AtDca = PropagateMCHParam(FwdtoMCH(mchTrackPar2), collision.posZ()); + auto mchTrackPar2AtVertex = PropagateMCHToVertex(mchTrackPar2, collision); + + // get the pre-stored MCH track parameters after refit + const auto mchTrackParNewIt1 = mMchTrackParsNew.find(mchIndex1); + if (mchTrackParNewIt1 == mMchTrackParsNew.end()) { + continue; + } + const auto mchTrackParNewIt2 = mMchTrackParsNew.find(mchIndex2); + if (mchTrackParNewIt2 == mMchTrackParsNew.end()) { + continue; + } + const auto& mchTrackParNew1 = mchTrackParNewIt1->second; + auto mchTrackParNew1AtDca = PropagateMCHParam(FwdtoMCH(mchTrackParNew1), collision.posZ()); + auto mchTrackParNew1AtVertex = PropagateMCHToVertex(mchTrackParNew1, collision); + const auto& mchTrackParNew2 = mchTrackParNewIt2->second; + auto mchTrackParNew2AtDca = PropagateMCHParam(FwdtoMCH(mchTrackParNew2), collision.posZ()); + auto mchTrackParNew2AtVertex = PropagateMCHToVertex(mchTrackParNew2, collision); + + FILL_DIMUON_PLOT(mchTrackPar1, mchTrackPar2, mchTrackPar1AtVertex, mchTrackPar2AtVertex, "dimuon/invariantMass_MuonKine_MuonCuts"); + FILL_DIMUON_PLOT(mchTrackParNew1, mchTrackParNew2, mchTrackParNew1AtVertex, mchTrackParNew2AtVertex, "dimuon/realign/invariantMass_MuonKine_MuonCuts"); + + FILL_DIMUON_DCA_PLOTS(mchTrackPar1, mchTrackPar2, + mchTrackPar1AtVertex, mchTrackPar2AtVertex, + mchTrackPar1AtDca, mchTrackPar2AtDca, + "dimuon/dcax_MuonKine_MuonCuts", "dimuon/dcay_MuonKine_MuonCuts"); + FILL_DIMUON_DCA_PLOTS(mchTrackParNew1, mchTrackParNew2, + mchTrackParNew1AtVertex, mchTrackParNew2AtVertex, + mchTrackParNew1AtDca, mchTrackParNew2AtDca, + "dimuon/realign/dcax_MuonKine_MuonCuts", "dimuon/realign/dcay_MuonKine_MuonCuts"); + + double mchAngle = getMuMuAngle(mchTrackPar1AtVertex, mchTrackPar2AtVertex); + double mchAngleNew = getMuMuAngle(mchTrackParNew1AtVertex, mchTrackParNew2AtVertex); + + try { + const auto& candidates1 = collisionInfo.globalMuonTracks.at(mchIndex1); + const auto& candidates2 = collisionInfo.globalMuonTracks.at(mchIndex2); + + auto fwdIndex1 = candidates1[0]; + auto fwdIndex2 = candidates2[0]; + + auto const& fwdTrack1 = muonTracks.rawIteratorAt(fwdIndex1); + auto const& fwdTrack2 = muonTracks.rawIteratorAt(fwdIndex2); + + bool isGoodGlobalMuon1 = IsGoodMuon(muonTrack1, collision, cfgTrackChi2MchUp, 0.f, cfgPtMchLow, {cfgEtaMftLow, cfgEtaMftUp}, {cfgRabsLow, cfgRabsUp}, fSigmaPdcaUp); + bool isGoodGlobalMuon2 = IsGoodMuon(muonTrack2, collision, cfgTrackChi2MchUp, 0.f, cfgPtMchLow, {cfgEtaMftLow, cfgEtaMftUp}, {cfgRabsLow, cfgRabsUp}, fSigmaPdcaUp); + bool goodGlobalMuonTracks = (isGoodGlobalMuon1 && isGoodGlobalMuon2); + + bool isGoodMatch1 = isGoodGlobalMatching(fwdTrack1, cfgDimuonMatchChi2Up.value); + bool isGoodMatch2 = isGoodGlobalMatching(fwdTrack2, cfgDimuonMatchChi2Up.value); + bool goodGlobalMuonMatches = (isGoodMatch1 && isGoodMatch2); + + if (!goodGlobalMuonTracks || !goodGlobalMuonMatches) { + continue; } + + FILL_DIMUON_PLOT(mchTrackPar1, mchTrackPar2, mchTrackPar1AtVertex, mchTrackPar2AtVertex, "dimuon/invariantMass_MuonKine_GlobalMuonCuts_GoodMatches"); + FILL_DIMUON_PLOT(mchTrackParNew1, mchTrackParNew2, mchTrackParNew1AtVertex, mchTrackParNew2AtVertex, "dimuon/realign/invariantMass_MuonKine_GlobalMuonCuts_GoodMatches"); + + FILL_DIMUON_DCA_PLOTS(mchTrackPar1, mchTrackPar2, + mchTrackPar1AtVertex, mchTrackPar2AtVertex, + mchTrackPar1AtDca, mchTrackPar2AtDca, + "dimuon/dcax_MuonKine_GlobalMuonCuts_GoodMatches", "dimuon/dcay_MuonKine_GlobalMuonCuts_GoodMatches"); + FILL_DIMUON_DCA_PLOTS(mchTrackParNew1, mchTrackParNew2, + mchTrackParNew1AtVertex, mchTrackParNew2AtVertex, + mchTrackParNew1AtDca, mchTrackParNew2AtDca, + "dimuon/realign/dcax_MuonKine_GlobalMuonCuts_GoodMatches", "dimuon/realign/dcay_MuonKine_GlobalMuonCuts_GoodMatches"); + + auto mftIndex1 = fwdTrack1.matchMFTTrackId(); + auto mftIndex2 = fwdTrack2.matchMFTTrackId(); + + const auto mftTrackPar1 = mMftTrackPars.at(mftIndex1); + auto fwdTrackPar1AtDca = PropagateMFTToDCA(mftTrackPar1, mchTrackPar1, collision, cfgVertexZshift); + auto fwdTrackPar1AtVertex = PropagateMFTToVertex(mftTrackPar1, mchTrackPar1, collision); + const auto mftTrackPar2 = mMftTrackPars.at(mftIndex2); + auto fwdTrackPar2AtDca = PropagateMFTToDCA(mftTrackPar2, mchTrackPar2, collision, cfgVertexZshift); + auto fwdTrackPar2AtVertex = PropagateMFTToVertex(mftTrackPar2, mchTrackPar2, collision); + const auto mftTrackParNew1 = mMftTrackParsNew.at(mftIndex1); + auto fwdTrackParNew1AtDca = PropagateMFTToDCA(mftTrackParNew1, mchTrackParNew1, collision, cfgVertexZshift); + auto fwdTrackParNew1AtVertex = PropagateMFTToVertex(mftTrackParNew1, mchTrackParNew1, collision); + const auto mftTrackParNew2 = mMftTrackParsNew.at(mftIndex2); + auto fwdTrackParNew2AtDca = PropagateMFTToDCA(mftTrackParNew2, mchTrackParNew2, collision, cfgVertexZshift); + auto fwdTrackParNew2AtVertex = PropagateMFTToVertex(mftTrackParNew2, mchTrackParNew2, collision); + + FILL_DIMUON_PLOT(mchTrackPar1, mchTrackPar2, fwdTrackPar1AtVertex, fwdTrackPar2AtVertex, "dimuon/invariantMass_ScaledMftKine_GlobalMuonCuts_GoodMatches"); + FILL_DIMUON_PLOT(mchTrackParNew1, mchTrackParNew2, fwdTrackParNew1AtVertex, fwdTrackParNew2AtVertex, "dimuon/realign/invariantMass_ScaledMftKine_GlobalMuonCuts_GoodMatches"); + + FILL_DIMUON_DCA_PLOTS(mchTrackPar1, mchTrackPar2, + fwdTrackPar1AtVertex, fwdTrackPar2AtVertex, + fwdTrackPar1AtDca, fwdTrackPar2AtDca, + "dimuon/dcax_ScaledMftKine_GlobalMuonCuts_GoodMatches", "dimuon/dcay_ScaledMftKine_GlobalMuonCuts_GoodMatches"); + FILL_DIMUON_DCA_PLOTS(mchTrackParNew1, mchTrackParNew2, + fwdTrackParNew1AtVertex, fwdTrackParNew2AtVertex, + fwdTrackParNew1AtDca, fwdTrackParNew2AtDca, + "dimuon/realign/dcax_ScaledMftKine_GlobalMuonCuts_GoodMatches", "dimuon/realign/dcay_ScaledMftKine_GlobalMuonCuts_GoodMatches"); + + double fwdAngle = getMuMuAngle(fwdTrackPar1AtVertex, fwdTrackPar2AtVertex); + double fwdAngleNew = getMuMuAngle(fwdTrackParNew1AtVertex, fwdTrackParNew2AtVertex); + + FILL_DIMUON_ANGLE_PLOT(mchTrackPar1, mchTrackPar2, fwdTrackPar1AtVertex, fwdTrackPar2AtVertex, mchAngle, fwdAngle, "dimuon/angle_GlobalMuonCuts_GoodMatches") + FILL_DIMUON_ANGLE_PLOT(mchTrackParNew1, mchTrackParNew2, fwdTrackParNew1AtVertex, fwdTrackParNew2AtVertex, mchAngleNew, fwdAngleNew, "dimuon/realign/angle_GlobalMuonCuts_GoodMatches") + } catch (const std::exception&) { + continue; } } } @@ -1897,14 +2348,16 @@ struct muonGlobalAlignment { // o2-linter: disable=name/workflow-file,name/struc } std::map collisionInfos; - InitCollisions(collisions, bcs, muonTracks, mftTracks, collisionInfos); + InitCollisions(collisions, bcs, muonTracks, clusters, mftTracks, collisionInfos); + + FillMftPlots(collisions, bcs, muonTracks, mftTracks, collisionInfos); - FillDCAPlots(collisions, bcs, muonTracks, mftTracks, collisionInfos); + FillMchPlots(collisions, bcs, muonTracks, clusters, collisionInfos); - FillResidualsPlots(collisions, bcs, muonTracks, clusters, collisionInfos); + FillDimuonPlots(collisions, muonTracks, collisionInfos); } - PROCESS_SWITCH(muonGlobalAlignment, processQA, "process qa", true); + PROCESS_SWITCH(muonGlobalAlignment, processQA, "processQA", true); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)