From 110f7bc99d9d954c1b574e74ade4358724c96e43 Mon Sep 17 00:00:00 2001 From: Maurice Coquet Date: Sat, 12 Sep 2026 16:59:01 +0200 Subject: [PATCH 1/9] processPPMuonRefit in tableMaker_withAssoc --- PWGDQ/TableProducer/tableMaker_withAssoc.cxx | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/PWGDQ/TableProducer/tableMaker_withAssoc.cxx b/PWGDQ/TableProducer/tableMaker_withAssoc.cxx index a6b7fc43145..bd556f0153b 100644 --- a/PWGDQ/TableProducer/tableMaker_withAssoc.cxx +++ b/PWGDQ/TableProducer/tableMaker_withAssoc.cxx @@ -508,7 +508,7 @@ struct TableMaker { bool enableEMCalHistos = (context.mOptions.get("processPPBarrelOnlyWithEMCal") || context.mOptions.get("processPPWithEMCal")); - bool enableMuonHistos = (context.mOptions.get("processPP") || context.mOptions.get("processPPWithEMCal") || context.mOptions.get("processPPWithFilter") || context.mOptions.get("processPPWithFilterMuonOnly") || context.mOptions.get("processPPWithFilterMuonMFT") || context.mOptions.get("processPPMuonOnly") || context.mOptions.get("processPPRealignedMuonOnly") || context.mOptions.get("processPPMuonMFT") || context.mOptions.get("processPPMuonMFTWithMultsExtra") || + bool enableMuonHistos = (context.mOptions.get("processPP") || context.mOptions.get("processPPWithEMCal") || context.mOptions.get("processPPWithFilter") || context.mOptions.get("processPPWithFilterMuonOnly") || context.mOptions.get("processPPWithFilterMuonMFT") || context.mOptions.get("processPPMuonOnly") || context.mOptions.get("processPPRealignedMuonOnly") || context.mOptions.get("processPPMuonMFT") || context.mOptions.get("processPPMuonMFTWithMultsExtra") || context.mOptions.get("processPPMuonRefit") || context.mOptions.get("processPbPb") || context.mOptions.get("processPbPbMuonOnly") || context.mOptions.get("processPbPbWithFilterMuonOnly") || context.mOptions.get("processPbPbStreamMuonOnly") || context.mOptions.get("processPbPbRealignedMuonOnly") || context.mOptions.get("processPbPbMuonMFT")); if (enableBarrelHistos) { @@ -1975,6 +1975,8 @@ struct TableMaker { if constexpr (static_cast(TMFTFillMap & VarManager::ObjTypes::MFTCov)) { if (fConfigVariousOptions.fUseML.value) { skimBestMuonMatchesML(muons, mftTracks, mftCovs, collision); + } else { + skimBestMuonMatches(muons); } } else { skimBestMuonMatches(muons); @@ -2237,6 +2239,7 @@ struct TableMaker { PROCESS_SWITCH(TableMaker, processPPRealignedMuonOnly, "Build realigned muon only DQ skimmed data model typically for pp/p-Pb and UPC Pb-Pb", false); PROCESS_SWITCH(TableMaker, processPPMuonMFT, "Build muon + mft DQ skimmed data model typically for pp/p-Pb and UPC Pb-Pb", false); PROCESS_SWITCH(TableMaker, processPPMuonMFTWithMultsExtra, "Build muon + mft DQ skimmed data model typically for pp/p-Pb and UPC Pb-Pb", false); + PROCESS_SWITCH(TableMaker, processPPMuonRefit, "Build muon + mft DQ skimmed data model with MFT covariances for global muon refit, typically for pp/p-Pb", false); PROCESS_SWITCH(TableMaker, processPbPb, "Build full DQ skimmed data model typically for Pb-Pb, w/o event filtering", false); PROCESS_SWITCH(TableMaker, processPbPbBarrelOnly, "Build barrel only DQ skimmed data model typically for Pb-Pb, w/o event filtering", false); PROCESS_SWITCH(TableMaker, processPbPbBarrelOnlyNoTOF, "Build barrel only DQ skimmed data model typically for Pb-Pb, w/o event filtering, no TOF", false); From d0091aea0524b5e4a047f7f13d633b4c7653894b Mon Sep 17 00:00:00 2001 From: Maurice Coquet Date: Sat, 12 Sep 2026 17:49:29 +0200 Subject: [PATCH 2/9] Moving computation of muon dca to FillTrackCollision --- PWGDQ/Core/VarManager.h | 16 ++-------------- 1 file changed, 2 insertions(+), 14 deletions(-) diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index 4b1641d139f..d3d4e929e15 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -1921,20 +1921,6 @@ void VarManager::FillPropagateMuon(const T& muon, const C& collision, float* val values[kTgl] = propmuon.getTgl(); values[kPhi] = propmuon.getPhi(); - // Redo propagation only for muon tracks - // propagation of MFT tracks alredy done in fwdtrack-extention task - if (static_cast(muon.trackType()) > 2) { - o2::dataformats::GlobalFwdTrack propmuonAtDCA = PropagateMuon(muon, collision, kToDCA); - o2::dataformats::GlobalFwdTrack propmuonAtRabs = PropagateMuon(muon, collision, kToRabs); - float dcaX = (propmuonAtDCA.getX() - collision.posX()); - float dcaY = (propmuonAtDCA.getY() - collision.posY()); - values[kMuonDCAx] = dcaX; - values[kMuonDCAy] = dcaY; - double xAbs = propmuonAtRabs.getX(); - double yAbs = propmuonAtRabs.getY(); - values[kMuonRAtAbsorberEnd] = std::sqrt(xAbs * xAbs + yAbs * yAbs); - } - const SMatrix55& cov = propmuon.getCovariances(); values[kMuonCXX] = cov(0, 0); values[kMuonCXY] = cov(1, 0); @@ -3483,6 +3469,8 @@ void VarManager::FillTrackCollision(T const& track, C const& collision, float* v float dcaX = (propmuonAtDCA.getX() - collision.posX()); float dcaY = (propmuonAtDCA.getY() - collision.posY()); float dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY); + values[kMuonDCAx] = dcaX; + values[kMuonDCAy] = dcaY; values[kMuonPDca] = track.p() * dcaXY; } } From bd1e8e9082a3aff3c622abe20221c151ef3ce95f Mon Sep 17 00:00:00 2001 From: Maurice Coquet Date: Sat, 12 Sep 2026 17:54:09 +0200 Subject: [PATCH 3/9] Removing obsolete FillMuonPDca --- PWGDQ/Core/VarManager.h | 21 --------------------- 1 file changed, 21 deletions(-) diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index d3d4e929e15..b134a9489b5 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -1449,8 +1449,6 @@ class VarManager : public TObject template static o2::track::TrackParCovFwd PropagateFwd(const T& track, const C& cov, float z); template - static void FillMuonPDca(const T& muon, const C& collision, float* values = nullptr); - template static void FillPropagateMuon(const T& muon, const C& collision, float* values = nullptr); template static void FillBC(T const& bc, float* values = nullptr); @@ -1879,25 +1877,6 @@ o2::track::TrackParCovFwd VarManager::PropagateFwd(const T& track, const C& cov, return fwdtrack; } -template -void VarManager::FillMuonPDca(const T& muon, const C& collision, float* values) -{ - if (!values) { - values = fgValues; - } - - if constexpr ((fillMap & MuonCov) > 0 || (fillMap & ReducedMuonCov) > 0) { - - o2::dataformats::GlobalFwdTrack propmuon = PropagateMuon(muon, collision); - o2::dataformats::GlobalFwdTrack propmuonAtDCA = PropagateMuon(muon, collision, kToDCA); - - float dcaX = (propmuonAtDCA.getX() - collision.posX()); - float dcaY = (propmuonAtDCA.getY() - collision.posY()); - float dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY); - values[kMuonPDca] = muon.p() * dcaXY; - } -} - template void VarManager::FillPropagateMuon(const T& muon, const C& collision, float* values) { From 67773a6c6f35e0b5d8bff5fb27ce6f6eb0479ed0 Mon Sep 17 00:00:00 2001 From: Maurice Coquet Date: Sat, 12 Sep 2026 20:31:32 +0200 Subject: [PATCH 4/9] Reorganizing computation of dca and pdca --- PWGDQ/Core/VarManager.h | 58 +- .../TableProducer/tableMakerMC_withAssoc.cxx | 10 +- PWGDQ/TableProducer/tableMaker_withAssoc.cxx | 13 +- PWGDQ/Tasks/tests/CMakeLists.txt | 170 ++ PWGDQ/Tasks/tests/FwdTrackReAlignTables.h | 90 + .../tests/fwdtrackToCollisionAssociator.cxx | 127 ++ PWGDQ/Tasks/tests/global-muon-matcher.cxx | 1681 +++++++++++++++++ 7 files changed, 2119 insertions(+), 30 deletions(-) create mode 100644 PWGDQ/Tasks/tests/CMakeLists.txt create mode 100644 PWGDQ/Tasks/tests/FwdTrackReAlignTables.h create mode 100644 PWGDQ/Tasks/tests/fwdtrackToCollisionAssociator.cxx create mode 100644 PWGDQ/Tasks/tests/global-muon-matcher.cxx diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index b134a9489b5..9b1642287e8 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -1835,11 +1835,7 @@ o2::dataformats::GlobalFwdTrack VarManager::PropagateMuon(const T& muon, const C o2::track::TrackParCovFwd fwdtrack = o2::aod::fwdtrackutils::getTrackParCovFwd3DShift(muon, xShift, yShift, zShift, muon); o2::dataformats::GlobalFwdTrack propmuon; if (static_cast(muon.trackType()) > 2) { - o2::dataformats::GlobalFwdTrack track; - track.setParameters(fwdtrack.getParameters()); - track.setZ(fwdtrack.getZ()); - track.setCovariances(fwdtrack.getCovariances()); - auto mchTrack = mMatching.FwdtoMCH(track); + auto mchTrack = mMatching.FwdtoMCH(fwdtrack); if (endPoint == kToVertex) { o2::mch::TrackExtrap::extrapToVertex(mchTrack, collision.posX(), collision.posY(), collision.posZ(), collision.covXX(), collision.covYY()); @@ -1854,17 +1850,12 @@ o2::dataformats::GlobalFwdTrack VarManager::PropagateMuon(const T& muon, const C o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchTrack, fgzMatching); } - auto proptrack = mMatching.MCHtoFwd(mchTrack); - propmuon.setParameters(proptrack.getParameters()); - propmuon.setZ(proptrack.getZ()); - propmuon.setCovariances(proptrack.getCovariances()); + propmuon = mMatching.MCHtoFwd(mchTrack); } else if (static_cast(muon.trackType()) < 2) { std::array dcaInfOrig{999.f, 999.f, 999.f}; fwdtrack.propagateToDCAhelix(fgMagField, {collision.posX(), collision.posY(), collision.posZ()}, dcaInfOrig); - propmuon.setParameters(fwdtrack.getParameters()); - propmuon.setZ(fwdtrack.getZ()); - propmuon.setCovariances(fwdtrack.getCovariances()); + propmuon = fwdtrack; } return propmuon; } @@ -1936,6 +1927,7 @@ void VarManager::FillGlobalMuonRefit(T1 const& muontrack, T2 const& mfttrack, co double pz = propmuon.getP() * std::cos(o2::constants::math::PIHalf - std::atan(mfttrack.tgl())); double pt = std::sqrt(std::pow(px, 2) + std::pow(py, 2)); auto mftprop = o2::aod::fwdtrackutils::getTrackParCovFwd3DShift(mfttrack, xShift, yShift, zShift); + mftprop.setInvQPt(static_cast(muontrack.sign()) / pt); values[kX] = mftprop.getX(); values[kY] = mftprop.getY(); values[kZ] = mftprop.getZ(); @@ -1944,6 +1936,12 @@ void VarManager::FillGlobalMuonRefit(T1 const& muontrack, T2 const& mfttrack, co values[kPz] = pz; values[kEta] = mftprop.getEta(); values[kPhi] = mftprop.getPhi(); + + // Helix DCA of the refitted global track w.r.t. the associated collision + std::array dca{999., 999., 999.}; + mftprop.propagateToDCAhelix(fgMagField, {collision.posX(), collision.posY(), collision.posZ()}, dca); + values[kMuonDCAx] = static_cast(dca[0]); + values[kMuonDCAy] = static_cast(dca[1]); } } @@ -1971,6 +1969,12 @@ void VarManager::FillGlobalMuonRefitCov(T1 const& muontrack, T2 const& mfttrack, values[kPz] = globalRefit.getPz(); values[kEta] = globalRefit.getEta(); values[kPhi] = globalRefit.getPhi(); + + // Helix DCA of the covariance-refitted global track w.r.t. the associated collision + std::array dca{999., 999., 999.}; + globalRefit.propagateToDCAhelix(fgMagField, {collision.posX(), collision.posY(), collision.posZ()}, dca); + values[kMuonDCAx] = static_cast(dca[0]); + values[kMuonDCAy] = static_cast(dca[1]); } } } @@ -3442,15 +3446,31 @@ void VarManager::FillTrackCollision(T const& track, C const& collision, float* v } } if constexpr ((fillMap & MuonCov) > 0 || (fillMap & MuonCovRealign) > 0 || (fillMap & ReducedMuonCov) > 0) { - - o2::dataformats::GlobalFwdTrack propmuonAtDCA = PropagateMuon(track, collision, kToDCA); - - float dcaX = (propmuonAtDCA.getX() - collision.posX()); - float dcaY = (propmuonAtDCA.getY() - collision.posY()); - float dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY); + float dcaX = 999.f; + float dcaY = 999.f; + if (static_cast(track.trackType()) <= 2) { + // Global muons: helix DCA w.r.t. the associated collision + float xShift = 0.f; + float yShift = 0.f; + float zShift = 0.f; + GetFwdShiftForY(track.y(), xShift, yShift, zShift); + o2::track::TrackParCovFwd fwdtrack = o2::aod::fwdtrackutils::getTrackParCovFwd3DShift(track, xShift, yShift, zShift, track); + std::array dca{999., 999., 999.}; + fwdtrack.propagateToDCAhelix(fgMagField, {collision.posX(), collision.posY(), collision.posZ()}, dca); + dcaX = static_cast(dca[0]); + dcaY = static_cast(dca[1]); + } else { + // MCH / standalone: DCA from MCH extrapolation + o2::dataformats::GlobalFwdTrack propmuonAtDCA = PropagateMuon(track, collision, kToDCA); + dcaX = propmuonAtDCA.getX() - collision.posX(); + dcaY = propmuonAtDCA.getY() - collision.posY(); + float dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY); + values[kMuonPDca] = track.p() * dcaXY; + } + values[kMuonDCAx] = dcaX; values[kMuonDCAy] = dcaY; - values[kMuonPDca] = track.p() * dcaXY; + } } diff --git a/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx b/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx index 620de63de2f..9e005ee3b9f 100644 --- a/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx +++ b/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx @@ -1133,10 +1133,10 @@ struct TableMakerMC { VarManager::FillTrack(muon); // NOTE: If a muon is associated to multiple collisions, depending on the selections, // it may be accepted for some associations and rejected for other - if (fConfigVariousOptions.fPropMuon) { - VarManager::FillPropagateMuon(muon, collision); + if (static_cast(muon.trackType()) > 2 && fConfigVariousOptions.fPropMuon) { + VarManager::FillPropagateMuon(muon, collision); } - // recalculte pDca and global muon kinematics + // recalculate pDca / DCA and global muon kinematics if (static_cast(muon.trackType()) < 2 && fConfigVariousOptions.fRefitGlobalMuon) { auto muontrack = muon.template matchMCHTrack_as(); if (muontrack.eta() < fConfigVariousOptions.fMuonMatchEtaMin || muontrack.eta() > fConfigVariousOptions.fMuonMatchEtaMax) { @@ -1262,10 +1262,10 @@ struct TableMakerMC { } VarManager::FillTrack(muon); - if (fConfigVariousOptions.fPropMuon) { + if (static_cast(muon.trackType()) > 2 && fConfigVariousOptions.fPropMuon) { VarManager::FillPropagateMuon(muon, collision); } - // recalculte pDca and global muon kinematics + // recalculate pDca / DCA and global muon kinematics int globalClusters = muon.nClusters(); if (static_cast(muon.trackType()) < 2 && fConfigVariousOptions.fRefitGlobalMuon) { auto muontrack = muon.template matchMCHTrack_as(); diff --git a/PWGDQ/TableProducer/tableMaker_withAssoc.cxx b/PWGDQ/TableProducer/tableMaker_withAssoc.cxx index bd556f0153b..616fb49d9e1 100644 --- a/PWGDQ/TableProducer/tableMaker_withAssoc.cxx +++ b/PWGDQ/TableProducer/tableMaker_withAssoc.cxx @@ -1687,10 +1687,10 @@ struct TableMaker { // NOTE: Muons are propagated to the current associated collisions. // So if a muon is associated to multiple collisions, depending on the selections, // it may be accepted for some associations and rejected for other - if (fConfigVariousOptions.fPropMuon) { - VarManager::FillPropagateMuon(muon, collision); + if (static_cast(muon.trackType()) > 2 && fConfigVariousOptions.fPropMuon) { + VarManager::FillPropagateMuon(muon, collision); } - // recalculate pDca and global muon kinematics + // recalculate pDca / DCA and global muon kinematics if (static_cast(muon.trackType()) < 2 && fConfigVariousOptions.fRefitGlobalMuon) { auto muontrack = muon.template matchMCHTrack_as(); if (muontrack.eta() < fConfigVariousOptions.fMuonMatchEtaMin || muontrack.eta() > fConfigVariousOptions.fMuonMatchEtaMax) { @@ -1700,6 +1700,7 @@ struct TableMaker { VarManager::FillTrackCollision(muontrack, collision); // NOTE: the MFT track originally associated to the MUON track is currently used in the global muon refit // Should MUON - MFT time ambiguities be taken into account ? + // Helix DCA is filled from the refitted parameters inside FillGlobalMuonRefit(Cov) if constexpr (static_cast(TMFTFillMap & VarManager::ObjTypes::MFTCov)) { auto const& mfttrackcov = mfCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]); VarManager::FillGlobalMuonRefitCov(muontrack, mfttrack, collision, mfttrackcov); @@ -1782,10 +1783,10 @@ struct TableMaker { } VarManager::FillTrack(muon); - if (fConfigVariousOptions.fPropMuon) { - VarManager::FillPropagateMuon(muon, collision); + if (static_cast(muon.trackType()) > 2 && fConfigVariousOptions.fPropMuon) { + VarManager::FillPropagateMuon(muon, collision); } - // recalculte pDca and global muon kinematics + // recalculate pDca / DCA and global muon kinematics int globalClusters = muon.nClusters(); if (static_cast(muon.trackType()) < 2 && fConfigVariousOptions.fRefitGlobalMuon) { auto muontrack = muon.template matchMCHTrack_as(); diff --git a/PWGDQ/Tasks/tests/CMakeLists.txt b/PWGDQ/Tasks/tests/CMakeLists.txt new file mode 100644 index 00000000000..29ac8ae0ba3 --- /dev/null +++ b/PWGDQ/Tasks/tests/CMakeLists.txt @@ -0,0 +1,170 @@ +# Copyright 2019-2020 CERN and copyright holders of ALICE O2. +# See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +# All rights not expressly granted are reserved. +# +# This software is distributed under the terms of the GNU General Public +# License v3 (GPL Version 3), copied verbatim in the file "COPYING". +# +# In applying this license CERN does not waive the privileges and immunities +# granted to it by virtue of its status as an Intergovernmental Organization +# or submit itself to any jurisdiction. + +o2physics_add_dpl_workflow(table-reader + SOURCES tableReader.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore O2Physics::MLCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(table-reader-with-assoc + SOURCES tableReader_withAssoc.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2Physics::AnalysisCCDB O2Physics::PWGDQCore O2Physics::MLCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(table-reader-with-assoc-direct + SOURCES tableReader_withAssoc_direct.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2Physics::AnalysisCCDB O2Physics::PWGDQCore O2Physics::MLCore O2::ReconstructionDataFormats O2::DetectorsCommonDataFormats O2::DetectorsVertexing O2Physics::EventFilteringUtils + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(efficiency + SOURCES dqEfficiency.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(efficiency-with-assoc + SOURCES dqEfficiency_withAssoc.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(efficiency-with-assoc-direct + SOURCES dqEfficiency_withAssoc_direct.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore O2::ReconstructionDataFormats O2::DetectorsCommonDataFormats O2::DetectorsVertexing + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(energy-correlator-direct + SOURCES dqEnergyCorrelator_direct.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(filter-pp + SOURCES filterPP.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(filter-pp-with-association + SOURCES filterPPwithAssociation.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(filter-pb-pb + SOURCES filterPbPb.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore O2Physics::SGCutParHolder + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(v0-selector + SOURCES v0selector.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DCAFitter O2Physics::AnalysisCore O2Physics::PWGDQCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(dalitz-selection + SOURCES DalitzSelection.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(flow + SOURCES dqFlow.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore O2Physics::GFWCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(task-muon-mch-trk-eff + SOURCES taskMuonMchTrkEfficiency.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(task-j-psi-hf + SOURCES taskJpsiHf.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(task-muon-dca + SOURCES muonDCA.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(correlation + SOURCES dqCorrelation.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore O2Physics::GFWCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(task-mch-align-record + SOURCES mchAlignRecord.cxx + PUBLIC_LINK_LIBRARIES + O2::Framework + O2Physics::AnalysisCore + O2Physics::PWGDQCore + O2::CommonUtils + O2::MCHClustering + O2::DPLUtils + O2::CCDB + O2::DataFormatsParameters + O2::MCHBase + O2::MCHTracking + O2::DataFormatsMCH + O2::DetectorsBase + O2::MCHGeometryTransformer + O2::MathUtils + O2::MCHAlign + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(task-muon-mid-eff + SOURCES MIDefficiency.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::MIDBase + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(task-fwd-track-pid + SOURCES taskFwdTrackPid.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(quarkonia-to-hyperons + SOURCES quarkoniaToHyperons.cxx + PUBLIC_LINK_LIBRARIES O2::DetectorsBase O2::Framework O2::DCAFitter KFParticle::KFParticle O2Physics::AnalysisCore O2Physics::MLCore O2Physics::EventFilteringUtils + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(model-converter-mult-pv + SOURCES ModelConverterMultPv.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::PWGDQCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(model-converter-event-extended + SOURCES ModelConverterEventExtended.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::PWGDQCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(model-converter-mc-reduced-event + SOURCES ModelConverterReducedMCEvents.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::PWGDQCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(tag-and-probe + SOURCES TagAndProbe.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2Physics::AnalysisCCDB O2Physics::PWGDQCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(qa-matching + SOURCES qaMatching.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(mft-mch-matcher + SOURCES mftMchMatcher.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(muon-global-alignment + SOURCES muonGlobalAlignment.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore O2::MCHGeometryTransformer + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(global-muon-matcher + SOURCES global-muon-matcher.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore O2::MCHGeometryTransformer + COMPONENT_NAME Analysis) diff --git a/PWGDQ/Tasks/tests/FwdTrackReAlignTables.h b/PWGDQ/Tasks/tests/FwdTrackReAlignTables.h new file mode 100644 index 00000000000..879d3fe7486 --- /dev/null +++ b/PWGDQ/Tasks/tests/FwdTrackReAlignTables.h @@ -0,0 +1,90 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file FwdTrackReAlignTables.h +/// \brief Table definitions for re-aligned forward tracks +/// \author Chi Zhang , CEA-Saclay + +#ifndef COMMON_DATAMODEL_FWDTRACKREALIGNTABLES_H_ +#define COMMON_DATAMODEL_FWDTRACKREALIGNTABLES_H_ + +#include + +namespace o2::aod +{ +namespace fwdtrackrealign +{ +DECLARE_SOA_COLUMN(IsRemovable, isRemovable, int); //! flag to check the refit status +} + +DECLARE_SOA_TABLE_FULL(StoredFwdTracksReAlign, "FwdTracksReAlign", "AOD", "FWDTRACKREALIGN", + o2::soa::Index<>, fwdtrack::CollisionId, fwdtrack::TrackType, + fwdtrack::X, fwdtrack::Y, fwdtrack::Z, fwdtrack::Phi, fwdtrack::Tgl, + fwdtrack::Signed1Pt, fwdtrack::NClusters, fwdtrack::PDca, fwdtrack::RAtAbsorberEnd, + fwdtrackrealign::IsRemovable, + fwdtrack::Px, + fwdtrack::Py, + fwdtrack::Pz, + fwdtrack::Sign, + fwdtrack::Chi2, fwdtrack::Chi2MatchMCHMID, fwdtrack::Chi2MatchMCHMFT, + fwdtrack::MatchScoreMCHMFT, fwdtrack::MFTTrackId, fwdtrack::MCHTrackId, + fwdtrack::MCHBitMap, fwdtrack::MIDBitMap, fwdtrack::MIDBoards, + fwdtrack::TrackTime, fwdtrack::TrackTimeRes); + +// extended table with expression columns that can be used as arguments of dynamic columns +DECLARE_SOA_EXTENDED_TABLE_USER(FwdTracksReAlign, StoredFwdTracksReAlign, "FWDTRKREALIGNEXT", //! + fwdtrack::Pt, + fwdtrack::Eta, + fwdtrack::P); + +DECLARE_SOA_TABLE_FULL(StoredFwdTrksCovReAlign, "FwdCovsReAlign", "AOD", "FWDCOVREALIGN", + fwdtrack::SigmaX, fwdtrack::SigmaY, fwdtrack::SigmaPhi, fwdtrack::SigmaTgl, fwdtrack::Sigma1Pt, + fwdtrack::RhoXY, fwdtrack::RhoPhiY, fwdtrack::RhoPhiX, fwdtrack::RhoTglX, fwdtrack::RhoTglY, + fwdtrack::RhoTglPhi, fwdtrack::Rho1PtX, fwdtrack::Rho1PtY, fwdtrack::Rho1PtPhi, fwdtrack::Rho1PtTgl); + +// extended table with expression columns that can be used as arguments of dynamic columns +DECLARE_SOA_EXTENDED_TABLE_USER(FwdTrksCovReAlign, StoredFwdTrksCovReAlign, "FWDCOVREALIGNEXT", //! + fwdtrack::CXX, + fwdtrack::CXY, + fwdtrack::CYY, + fwdtrack::CPhiX, + fwdtrack::CPhiY, + fwdtrack::CPhiPhi, + fwdtrack::CTglX, + fwdtrack::CTglY, + fwdtrack::CTglPhi, + fwdtrack::CTglTgl, + fwdtrack::C1PtX, + fwdtrack::C1PtY, + fwdtrack::C1PtPhi, + fwdtrack::C1PtTgl, + fwdtrack::C1Pt21Pt2); + +using FwdTrackRealign = FwdTracksReAlign::iterator; +using FwdTrkCovRealign = FwdTrksCovReAlign::iterator; +using FullFwdTracksRealign = soa::Join; +using FullFwdTrackRealign = FullFwdTracksRealign::iterator; + +// ambiguity table for realigned muons +namespace fwdtrackrealignambiguous +{ +DECLARE_SOA_INDEX_COLUMN_FULL(FwdTrackRealign, fwdTrackRealign, int, FwdTracksReAlign, ""); //! FwdTracksReAlign index +DECLARE_SOA_SLICE_INDEX_COLUMN(BC, bc); + +} // namespace fwdtrackrealignambiguous + +DECLARE_SOA_TABLE(AmbiguousFwdTrksReAlign, "AOD", "AMBIFWDREALIGN", //! Table for FwdTracksReAlign which are not associated with a collision + o2::soa::Index<>, fwdtrackrealignambiguous::FwdTrackRealignId, fwdtrackrealignambiguous::BCIdSlice); + +using AmbiguousFwdTrkRealign = AmbiguousFwdTrksReAlign::iterator; +} // namespace o2::aod + +#endif // COMMON_DATAMODEL_FWDTRACKREALIGNTABLES_H_ diff --git a/PWGDQ/Tasks/tests/fwdtrackToCollisionAssociator.cxx b/PWGDQ/Tasks/tests/fwdtrackToCollisionAssociator.cxx new file mode 100644 index 00000000000..1de267d0e60 --- /dev/null +++ b/PWGDQ/Tasks/tests/fwdtrackToCollisionAssociator.cxx @@ -0,0 +1,127 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// \file fwdtrackToCollisionAssociator.cxx +/// \brief Associates fwd and MFT tracks to collisions considering ambiguities +/// \author Sarah Herrmann , IP2I Lyon +/// \author Maurice Coquet , CEA-Saclay/Irfu + +#include "Common/Core/CollisionAssociation.h" +#include "Common/DataModel/CollisionAssociationTables.h" +#include "Common/DataModel/FwdTrackReAlignTables.h" + +#include +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::aod; + +struct FwdTrackToCollisionAssociation { + Produces fwdassociation; + Produces fwdreverseIndices; + Produces mftassociation; + Produces mftreverseIndices; + + Configurable nSigmaForTimeCompat{"nSigmaForTimeCompat", 4.f, "number of sigmas for time compatibility"}; + Configurable timeMargin{"timeMargin", 0.f, "time margin in ns added to uncertainty because of uncalibrated TPC"}; + Configurable includeUnassigned{"includeUnassigned", false, "consider also tracks which are not assigned to any collision"}; + Configurable fillTableOfCollIdsPerTrack{"fillTableOfCollIdsPerTrack", false, "fill additional table with vector of collision ids per track"}; + Configurable bcWindowForOneSigma{"bcWindowForOneSigma", 115, "BC window to be multiplied by the number of sigmas to define maximum window to be considered"}; + + CollisionAssociation collisionAssociator; + + Preslice muonsPerCollisions = aod::fwdtrack::collisionId; + Preslice realignmuonsPerCollisions = aod::fwdtrack::collisionId; + Preslice mftsPerCollisions = aod::fwdtrack::collisionId; + + void init(InitContext const&) + { + if (doprocessFwdAssocWithTime && doprocessFwdStandardAssoc) { + LOGP(fatal, "Exactly one process function between standard and time-based association should be enabled!"); + } + if (doprocessMFTAssocWithTime && doprocessMFTStandardAssoc) { + LOGP(fatal, "Exactly one process function between standard and time-based association should be enabled!"); + } + + if (!(doprocessMFTAssocWithTime || doprocessMFTStandardAssoc || doprocessFwdAssocWithTime || doprocessFwdStandardAssoc || doprocessFwdRealignAssocWithTime || doprocessFwdRealignStandardAssoc)) { + LOGP(fatal, "At least one process function should be enabled!"); + } + + // set options in track-to-collision association + collisionAssociator.setNumSigmaForTimeCompat(nSigmaForTimeCompat); + collisionAssociator.setTimeMargin(timeMargin); + collisionAssociator.setTrackSelectionOptionForStdAssoc(track_association::TrackSelection::None); + collisionAssociator.setUsePvAssociation(track_association::PVContrReassocOpt::Disabled); + collisionAssociator.setIncludeUnassigned(includeUnassigned); + collisionAssociator.setFillTableOfCollIdsPerTrack(fillTableOfCollIdsPerTrack); + collisionAssociator.setBcWindow(bcWindowForOneSigma); + } + + void processFwdAssocWithTime(Collisions const& collisions, + FwdTracks const& muons, + AmbiguousFwdTracks const& ambiTracksFwd, + BCs const& bcs) + { + collisionAssociator.runAssocWithTime(collisions, muons, muons, ambiTracksFwd, bcs, fwdassociation, fwdreverseIndices); + } + PROCESS_SWITCH(FwdTrackToCollisionAssociation, processFwdAssocWithTime, "Use fwdtrack-to-collision association based on time", true); + + void processFwdStandardAssoc(Collisions const& collisions, + FwdTracks const& muons) + { + collisionAssociator.runStandardAssoc(collisions, muons, muonsPerCollisions, fwdassociation, fwdreverseIndices); + } + PROCESS_SWITCH(FwdTrackToCollisionAssociation, processFwdStandardAssoc, "Use standard fwdtrack-to-collision association", false); + + void processFwdRealignAssocWithTime(Collisions const& collisions, + FwdTracksReAlign const& muons, + AmbiguousFwdTrksReAlign const& ambiTracksFwd, + BCs const& bcs) + { + collisionAssociator.runAssocWithTime(collisions, muons, muons, ambiTracksFwd, bcs, fwdassociation, fwdreverseIndices); + } + PROCESS_SWITCH(FwdTrackToCollisionAssociation, processFwdRealignAssocWithTime, "Use fwdrealigntrack-to-collision association based on time", false); + + void processFwdRealignStandardAssoc(Collisions const& collisions, + FwdTracksReAlign const& muons) + { + collisionAssociator.runStandardAssoc(collisions, muons, realignmuonsPerCollisions, fwdassociation, fwdreverseIndices); + } + PROCESS_SWITCH(FwdTrackToCollisionAssociation, processFwdRealignStandardAssoc, "Use standard fwdrealigntrack-to-collision association", false); + + void processMFTAssocWithTime(Collisions const& collisions, + MFTTracks const& tracks, + AmbiguousMFTTracks const& ambiguousTracks, + BCs const& bcs) + { + collisionAssociator.runAssocWithTime(collisions, tracks, tracks, ambiguousTracks, bcs, mftassociation, mftreverseIndices); + } + PROCESS_SWITCH(FwdTrackToCollisionAssociation, processMFTAssocWithTime, "Use MFTtrack-to-collision association based on time", true); + + void processMFTStandardAssoc(Collisions const& collisions, + MFTTracks const& tracks) + { + collisionAssociator.runStandardAssoc(collisions, tracks, mftsPerCollisions, mftassociation, mftreverseIndices); + } + PROCESS_SWITCH(FwdTrackToCollisionAssociation, processMFTStandardAssoc, "Use standard mfttrack-to-collision association", false); +}; + +//________________________________________________________________________________________________________________________ +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +} diff --git a/PWGDQ/Tasks/tests/global-muon-matcher.cxx b/PWGDQ/Tasks/tests/global-muon-matcher.cxx new file mode 100644 index 00000000000..28b8ef02fcd --- /dev/null +++ b/PWGDQ/Tasks/tests/global-muon-matcher.cxx @@ -0,0 +1,1681 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +// +/// \file global-muon-matcher.cxx +/// \brief Task for analysis MFT-MCH muon matching +/// \author Andrea Ferrero +/// +#include "PWGDQ/Core/MuonMatchingMlResponse.h" +#include "PWGDQ/Core/VarManager.h" + +#include "Common/Core/fwdtrackUtilities.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/FwdTrackReAlignTables.h" +#include "Tools/ML/MlResponse.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::aod; + +namespace o2::aod::globalmuonmatching +{ +DECLARE_SOA_COLUMN(MchTrackId, mchTrackId, int64_t); +DECLARE_SOA_COLUMN(MftTrackId, mftTrackId, int64_t); +DECLARE_SOA_COLUMN(MatchChi2, matchChi2, float); +DECLARE_SOA_COLUMN(MatchScore, matchScore, float); +DECLARE_SOA_COLUMN(MatchRanking, matchRanking, int32_t); +DECLARE_SOA_COLUMN(IsTagged, isTagged, bool); +} // namespace o2::aod::globalmuonmatching + +namespace o2::aod +{ +DECLARE_SOA_TABLE(GlobalMuonMatchCandidates, "AOD", "GMCAND", + o2::soa::Index<>, + globalmuonmatching::MchTrackId, + globalmuonmatching::MftTrackId, + globalmuonmatching::MatchChi2, globalmuonmatching::MatchScore, globalmuonmatching::MatchRanking, + globalmuonmatching::IsTagged); + +namespace globalmuonmatching +{ +DECLARE_SOA_ARRAY_INDEX_COLUMN(GlobalMuonMatchCandidate, globalMuonMatchCandidate); //! Array of GlobalMuonMatchCandidates indices +} // namespace globalmuonmatching + +DECLARE_SOA_TABLE(FwdTrkMatchCands, "AOD", "FWDTRKMATCHCAND", //! Vectors of match-candidate indices stored per fwdtrack + globalmuonmatching::GlobalMuonMatchCandidateIds, o2::soa::Marker<3>); + +} // namespace o2::aod + +using MyEvents = soa::Join; +using MyMuons = soa::Join; +using MyMFTs = aod::MFTTracks; +using MyMFTCovariances = aod::MFTTracksCov; + +using SMatrix55Sym = o2::track::SMatrix55Sym; +using SMatrix55Std = o2::track::SMatrix55Std; +using SMatrix5 = o2::track::SMatrix5; + +constexpr std::array NDetElemCh = {4, 4, 4, 4, 18, 18, 26, 26, 26, 26}; +constexpr std::array SNDetElemCh = {0, 4, 8, 12, 16, 34, 52, 78, 104, 130, 156}; + +struct GlobalMuonMatching { + + static constexpr int GlobalTrackTypeMax = 2; + static constexpr int MchMidTrackType = 3; + static constexpr int NMchChambers = 10; + static constexpr int MchDetElemNumberingBase = 100; + static constexpr int NMchDetElems = 156; + static constexpr int MinRemovableTrackClusters = 10; + static constexpr int ThetaAbsBoundaryDeg = 3; + static constexpr double SlopeResolutionZ = 535.; + static constexpr float MatchingPlaneDefaultZ = -77.5; + + struct MatchingCandidate { + int64_t muonTrackId{-1}; + int64_t mftTrackId{-1}; + double matchScore{-1}; + double matchChi2{-1}; + int matchRanking{-1}; + }; + + //// Variables for selecting tagged muons + struct : ConfigurableGroup { + Configurable cfgMuonTaggingNCrossedMftPlanesLow{"cfgMuonTaggingNCrossedMftPlanesLow", 5, ""}; + Configurable cfgMuonTaggingTrackChi2MchUp{"cfgMuonTaggingTrackChi2MchUp", 5.f, ""}; + Configurable cfgMuonTaggingPMchLow{"cfgMuonTaggingPMchLow", 0.0f, ""}; + Configurable cfgMuonTaggingPtMchLow{"cfgMuonTaggingPtMchLow", 0.7f, ""}; + Configurable cfgMuonTaggingEtaMchLow{"cfgMuonTaggingEtaMchLow", -3.6f, ""}; + Configurable cfgMuonTaggingEtaMchUp{"cfgMuonTaggingEtaMchUp", -2.5f, ""}; + Configurable cfgMuonTaggingRabsLow{"cfgMuonTaggingRabsLow", 17.6f, ""}; + Configurable cfgMuonTaggingRabsUp{"cfgMuonTaggingRabsUp", 89.5f, ""}; + Configurable cfgMuonTaggingPdcaUp{"cfgMuonTaggingPdcaUp", 4.f, ""}; + Configurable cfgMuonTaggingRadiusAtMftFrontLow{"cfgMuonTaggingRadiusAtMftFrontLow", 3.f, ""}; + Configurable cfgMuonTaggingRadiusAtMftFrontUp{"cfgMuonTaggingRadiusAtMftFrontUp", 9.f, ""}; + Configurable cfgMuonTaggingRadiusAtMftBackLow{"cfgMuonTaggingRadiusAtMftBackLow", 5.f, ""}; + Configurable cfgMuonTaggingRadiusAtMftBackUp{"cfgMuonTaggingRadiusAtMftBackUp", 12.f, ""}; + } configMuonTagging; + + //// Variables for MCH realignment + struct : ConfigurableGroup { + Configurable cfgEnableMCHRealign{"cfgEnableMCHRealign", true, "Enable re-alignment of MCH clusters and tracks"}; + Configurable cfgGeoRefPath{"cfgGeoRefPath", "GLO/Config/GeometryAligned", "Path of the reference geometry file"}; + Configurable cfgGeoNewPath{"cfgGeoNewPath", "GLO/Config/GeometryAligned", "Path of the new geometry file"}; + Configurable cfgCcdbNoLaterThanRef{"cfgCcdbNoLaterThanRef", std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(), "latest acceptable timestamp of creation for the object of reference basis"}; + Configurable cfgCcdbNoLaterThanNew{"cfgCcdbNoLaterThanNew", std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(), "latest acceptable timestamp of creation for the object of new basis"}; + Configurable cfgChamberResolutionX{"cfgChamberResolutionX", 0.04, "Chamber resolution along X configuration for refit"}; // 0.4cm pp, 0.2cm PbPb + Configurable cfgChamberResolutionY{"cfgChamberResolutionY", 0.04, "Chamber resolution along Y configuration for refit"}; // 0.4cm pp, 0.2cm PbPb + Configurable cfgSigmaCutImprove{"cfgSigmaCutImprove", 6., "Sigma cut for track improvement"}; // 6 for pp, 4 for PbPb + } configMchRealign; + + //// Variables for MFT alignment corrections + struct : ConfigurableGroup { + Configurable cfgEnableMftAlignmentCorrections{"cfgEnableMftAlignmentCorrections", true, "Enable alignment corrections for the MFT tracks"}; + // slope corrections + // Configurable cfgMFTAlignmentCorrXSlopeTop{"cfgMFTAlignmentCorrXSlopeTop", (-0.0006696 - 0.0005621) / 2.f, "MFT X slope correction - top half"}; + // Configurable cfgMFTAlignmentCorrXSlopeBottom{"cfgMFTAlignmentCorrXSlopeBottom", (0.00105 + 0.001007) / 2.f, "MFT X slope correction - bottom half"}; + // Configurable cfgMFTAlignmentCorrYSlopeTop{"cfgMFTAlignmentCorrYSlopeTop", (-0.002299 - 0.002442) / 2.f, "MFT Y slope correction - top half"}; + // Configurable cfgMFTAlignmentCorrYSlopeBottom{"cfgMFTAlignmentCorrYSlopeBottom", (-0.0005339 - 0.0006921) / 2.f, "MFT Y slope correction - bottom half"}; + Configurable cfgMFTAlignmentCorrXSlopeTop{"cfgMFTAlignmentCorrXSlopeTop", 0.f, "MFT X slope correction - top half"}; + Configurable cfgMFTAlignmentCorrXSlopeBottom{"cfgMFTAlignmentCorrXSlopeBottom", 0.f, "MFT X slope correction - bottom half"}; + Configurable cfgMFTAlignmentCorrYSlopeTop{"cfgMFTAlignmentCorrYSlopeTop", 0.f, "MFT Y slope correction - top half"}; + Configurable cfgMFTAlignmentCorrYSlopeBottom{"cfgMFTAlignmentCorrYSlopeBottom", 0.f, "MFT Y slope correction - bottom half"}; + // offset corrections + Configurable cfgMFTAlignmentCorrXOffsetTop{"cfgMFTAlignmentCorrXOffsetTop", 0.f, "MFT X offset correction - top half"}; + Configurable cfgMFTAlignmentCorrXOffsetBottom{"cfgMFTAlignmentCorrXOffsetBottom", 0.f, "MFT X offset correction - bottom half"}; + Configurable cfgMFTAlignmentCorrYOffsetTop{"cfgMFTAlignmentCorrYOffsetTop", 0.f, "MFT Y offset correction - top half"}; + Configurable cfgMFTAlignmentCorrYOffsetBottom{"cfgMFTAlignmentCorrYOffsetBottom", 0.f, "MFT Y offset correction - bottom half"}; + } configMftAlignmentCorrections; + + // Variables for CCDB objects access and retrieval + struct : ConfigurableGroup { + Configurable cfgCcdbUrl{"cfgCcdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable cfgCcdbNoLaterThan{"cfgCcdbNoLaterThan", std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(), "latest acceptable timestamp of creation for the object"}; + Configurable cfgGrpPath{"cfgGrpPath", "GLO/GRP/GRP", "Path of the grp file"}; + Configurable cfgGeoPath{"cfgGeoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"}; + Configurable cfgGrpMagPath{"cfgGrpMagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; + } configCcdb; + + // Matching strategy for the *custom* matches (production baseline is always computed). + // 0 = chi2 (runChi2Matching), 1 = ML (runMlMatching) + struct : ConfigurableGroup { + Configurable cfgCustomMatchingStrategy{"cfgCustomMatchingStrategy", 0, "0=chi2, 1=ML for custom matches"}; + Configurable cfgIncludeGlobalMuonsInFwdTracks{"cfgIncludeGlobalMuonsInFwdTracks", false, "Include MFT-MCH-MID global muons in GMMCANDTRK table"}; + Configurable cfgMaxCandidatesPerMchTrack{"cfgMaxCandidatesPerMchTrack", -1, "Maximum number of match candidates stored per MCH track (-1: no limit)"}; + Configurable cfgMatchAllTracks{"cfgMatchAllTracks", false, "If true the matching is performed considering all the MFT tracks for which the covariances are available; if false the matching is performed considering only the global forward tracks stored at production"}; + } configMatching; + + double mBzAtMftCenter{0}; + + using MatchingFunc = std::function(const o2::track::TrackParCovFwd& mchtrack, const o2::track::TrackParCovFwd& mfttrack)>; + std::map mMatchingFunctionMap; ///< MFT-MCH Matching function + + // Chi2 matching interface (single configurable method) + struct : ConfigurableGroup { + Configurable cfgChi2FunctionLabel{"cfgChi2FunctionLabel", std::string{"ProdAll"}, "Text label identifying the chi2 matching method"}; + Configurable cfgChi2FunctionName{"cfgChi2FunctionName", std::string{"prod"}, "Name of the chi2 matching function"}; + Configurable cfgChi2FunctionMatchingPlaneZ{"cfgChi2FunctionMatchingPlaneZ", static_cast(o2::mft::constants::mft::LayerZCoordinate()[9]), "Z position of the matching plane"}; + } configChi2MatchingOptions; + + // ML interface (single configurable model) + struct : ConfigurableGroup { + Configurable cfgMlModelLabel{"cfgMlModelLabel", std::string{""}, "Text label identifying this ML model"}; + Configurable cfgMlModelPathCcdb{"cfgMlModelPathCcdb", "Users/m/mcoquet/MLTest", "Path of model on CCDB"}; + Configurable cfgMlModelName{"cfgMlModelName", "model.onnx", "ONNX file name (if not from CCDB full path)"}; + Configurable> cfgMlInputFeatures{"cfgMlInputFeatures", std::vector{"chi2MCHMFT"}, "Names of ML model input features"}; + Configurable cfgMlModelMatchingPlaneZ{"cfgMlModelMatchingPlaneZ", static_cast(o2::mft::constants::mft::LayerZCoordinate()[9]), "Z position of the matching plane"}; + } configMlOptions; + + std::vector binsPtMl; + std::array cutValues{}; + std::vector cutDirMl; + bool hasActiveChi2Matching{false}; + std::string activeChi2FunctionName; + double activeChi2MatchingPlaneZ{0.}; + + bool hasActiveMlMatching{false}; + o2::analysis::MlResponseMFTMuonMatch activeMlResponse; + double activeMlMatchingPlaneZ{0.}; + + int mRunNumber{0}; // needed to detect if the run changed and trigger update of magnetic field + + Service ccdbManager{}; + o2::ccdb::CcdbApi fCCDBApi; + + // 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 score + // the map key is the MCH(-MID) track global index + using MatchingCandidates = std::map>; + std::map> mMatchingCandidates; + + 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 *this; } + + private: + int nClusters{-1}; + bool removable{false}; + }; + + std::unordered_map mMchTrackPars; + std::unordered_map mMftTrackPars; + + std::unordered_map mftTrackCovs; + + Produces globalMuonMatchCandidates; + Produces fwdTrkMatchCands; + Produces gmCandidateFwdTracks; + Produces gmCandidateFwdTracksCov; + Produces gmAmbiguousFwdTracksReAlign; + + int32_t mGmmCandFwdTrackRowIndex{0}; + std::unordered_map> mAmbBcSliceByFwdTrackId; + bool mHasLastMchAmbiguousBcSlice{false}; + std::array mLastMchAmbiguousBcSlice{}; + + int32_t mMatchCandidateCounter{0}; + std::unordered_map> mMchTrackToCandidateIndices; + std::unordered_map> mMchTrackMatchingCandidates; + std::unordered_map mFwdTrackToGmmCandTrkIndex; + + mch::TrackFitter trackFitter; // Track fitter from MCH tracking library + mch::geo::TransformationCreator transformation; + std::map transformRef; // reference geometry w.r.t track data + std::map transformNew; // new geometry + double mImproveCutChi2{0.}; // Chi2 cut for track improvement. + TGeoManager* geoNew = nullptr; + TGeoManager* geoRef = nullptr; + globaltracking::MatchGlobalFwd mMatching; + + Preslice perMuon = aod::fwdtrkcl::fwdtrackId; + + template + o2::mch::TrackParam fwdToMch(const T& fwdtrack) + { + // Convert Forward Track parameters and covariances matrix to the + // MCH track format. + + // Parameter conversion + const double x2 = fwdtrack.getPhi(); + const double x3 = fwdtrack.getTanl(); + const double x4 = fwdtrack.getInvQPt(); + + const auto sinX2 = std::sin(x2); + const auto cosX2 = std::cos(x2); + + const double alpha1 = cosX2 / x3; + const double alpha3 = sinX2 / x3; + const double alpha4 = x4 / std::sqrt(x3 * x3 + sinX2 * sinX2); + + const auto kNorm = std::sqrt(x3 * x3 + sinX2 * sinX2); + const auto kNorm3 = kNorm * kNorm * kNorm; + + // Covariances matrix conversion + SMatrix55Std jacobian; + SMatrix55Sym covariances; + + covariances(0, 0) = fwdtrack.getCovariances()(0, 0); + covariances(0, 1) = fwdtrack.getCovariances()(0, 1); + covariances(0, 2) = fwdtrack.getCovariances()(0, 2); + covariances(0, 3) = fwdtrack.getCovariances()(0, 3); + covariances(0, 4) = fwdtrack.getCovariances()(0, 4); + + covariances(1, 1) = fwdtrack.getCovariances()(1, 1); + covariances(1, 2) = fwdtrack.getCovariances()(1, 2); + covariances(1, 3) = fwdtrack.getCovariances()(1, 3); + covariances(1, 4) = fwdtrack.getCovariances()(1, 4); + + covariances(2, 2) = fwdtrack.getCovariances()(2, 2); + covariances(2, 3) = fwdtrack.getCovariances()(2, 3); + covariances(2, 4) = fwdtrack.getCovariances()(2, 4); + + covariances(3, 3) = fwdtrack.getCovariances()(3, 3); + covariances(3, 4) = fwdtrack.getCovariances()(3, 4); + + covariances(4, 4) = fwdtrack.getCovariances()(4, 4); + + jacobian(0, 0) = 1; + + jacobian(1, 2) = -sinX2 / x3; + jacobian(1, 3) = -cosX2 / (x3 * x3); + + jacobian(2, 1) = 1; + + jacobian(3, 2) = cosX2 / x3; + jacobian(3, 3) = -sinX2 / (x3 * x3); + + jacobian(4, 2) = -x4 * sinX2 * cosX2 / kNorm3; + jacobian(4, 3) = -x3 * x4 / kNorm3; + jacobian(4, 4) = 1 / kNorm; + // jacobian*covariances*jacobian^T + covariances = ROOT::Math::Similarity(jacobian, covariances); + + std::array cov = {covariances(0, 0), covariances(1, 0), covariances(1, 1), covariances(2, 0), covariances(2, 1), covariances(2, 2), covariances(3, 0), covariances(3, 1), covariances(3, 2), covariances(3, 3), covariances(4, 0), covariances(4, 1), covariances(4, 2), covariances(4, 3), covariances(4, 4)}; + std::array param = {fwdtrack.getX(), alpha1, fwdtrack.getY(), alpha3, alpha4}; + + o2::mch::TrackParam convertedTrack(fwdtrack.getZ(), param.data(), cov.data()); + return {convertedTrack}; + } + + o2::track::TrackParCovFwd mchToFwd(const o2::mch::TrackParam& mchParam) + { + // Convert a MCH Track parameters and covariances matrix to the + // Forward track format. Must be called after propagation though the absorber + + o2::track::TrackParCovFwd convertedTrack; + + // Parameter conversion + const double alpha1 = mchParam.getNonBendingSlope(); + const double alpha3 = mchParam.getBendingSlope(); + const double alpha4 = mchParam.getInverseBendingMomentum(); + + const double x2 = std::atan2(-alpha3, -alpha1); + const double x3 = -1. / std::sqrt(alpha3 * alpha3 + alpha1 * alpha1); + const double x4 = alpha4 * -x3 * std::sqrt(1 + alpha3 * alpha3); + + const auto kNorm = alpha1 * alpha1 + alpha3 * alpha3; + const auto kNorm32 = kNorm * std::sqrt(kNorm); + const auto slopeLen = std::sqrt(alpha3 * alpha3 + 1); + + // Covariances matrix conversion + SMatrix55Std jacobian; + SMatrix55Sym covariances; + + covariances(0, 0) = mchParam.getCovariances()(0, 0); + covariances(0, 1) = mchParam.getCovariances()(0, 1); + covariances(0, 2) = mchParam.getCovariances()(0, 2); + covariances(0, 3) = mchParam.getCovariances()(0, 3); + covariances(0, 4) = mchParam.getCovariances()(0, 4); + + covariances(1, 1) = mchParam.getCovariances()(1, 1); + covariances(1, 2) = mchParam.getCovariances()(1, 2); + covariances(1, 3) = mchParam.getCovariances()(1, 3); + covariances(1, 4) = mchParam.getCovariances()(1, 4); + + covariances(2, 2) = mchParam.getCovariances()(2, 2); + covariances(2, 3) = mchParam.getCovariances()(2, 3); + covariances(2, 4) = mchParam.getCovariances()(2, 4); + + covariances(3, 3) = mchParam.getCovariances()(3, 3); + covariances(3, 4) = mchParam.getCovariances()(3, 4); + + covariances(4, 4) = mchParam.getCovariances()(4, 4); + + jacobian(0, 0) = 1; + + jacobian(1, 2) = 1; + + jacobian(2, 1) = -alpha3 / kNorm; + jacobian(2, 3) = alpha1 / kNorm; + + jacobian(3, 1) = alpha1 / kNorm32; + jacobian(3, 3) = alpha3 / kNorm32; + + jacobian(4, 1) = -alpha1 * alpha4 * slopeLen / kNorm32; + jacobian(4, 3) = alpha3 * alpha4 * (1 / (std::sqrt(kNorm) * slopeLen) - slopeLen / kNorm32); + jacobian(4, 4) = slopeLen / std::sqrt(kNorm); + + // jacobian*covariances*jacobian^T + covariances = ROOT::Math::Similarity(jacobian, covariances); + + // Set output + convertedTrack.setX(mchParam.getNonBendingCoor()); + convertedTrack.setY(mchParam.getBendingCoor()); + convertedTrack.setZ(mchParam.getZ()); + convertedTrack.setPhi(x2); + convertedTrack.setTanl(x3); + convertedTrack.setInvQPt(x4); + convertedTrack.setCharge(mchParam.getCharge()); + convertedTrack.setCovariances(covariances); + + return convertedTrack; + } + + int getDetElemId(int iDetElemNumber) + { + // make sure detector number is valid + if (iDetElemNumber < SNDetElemCh[0] || + iDetElemNumber >= SNDetElemCh[NMchChambers]) { + LOGF(fatal, "Invalid detector element number: %d", iDetElemNumber); + } + /// get det element number from ID + // get chamber and element number in chamber + int iCh = 0; + int iDet = 0; + for (int i = 1; i <= NMchChambers; i++) { + if (iDetElemNumber < SNDetElemCh[i]) { + iCh = i; + iDet = iDetElemNumber - SNDetElemCh[i - 1]; + break; + } + } + + // make sure detector index is valid + if (iCh <= 0 || iCh > NMchChambers || iDet >= NDetElemCh[iCh - 1]) { + LOGF(fatal, "Invalid detector element id: %d", MchDetElemNumberingBase * iCh + iDet); + } + + // add number of detectors up to this chamber + return MchDetElemNumberingBase * iCh + iDet; + } + + bool removeTrack(mch::Track& track) + { + // Refit track with re-aligned clusters + bool shouldRemoveTrack = false; + try { + trackFitter.fit(track, false); + } catch (std::exception const& e) { + shouldRemoveTrack = true; + return shouldRemoveTrack; + } + + auto itStartingParam = std::prev(track.rend()); + + while (true) { + + try { + trackFitter.fit(track, true, false, (itStartingParam == track.rbegin()) ? nullptr : &itStartingParam); + } catch (std::exception const&) { + shouldRemoveTrack = true; + break; + } + + double worstLocalChi2 = -1.0; + + track.tagRemovableClusters(0x1F, false); + + auto itWorstParam = track.end(); + + for (auto itParam = track.begin(); itParam != track.end(); ++itParam) { + if (itParam->getLocalChi2() > worstLocalChi2) { + worstLocalChi2 = itParam->getLocalChi2(); + itWorstParam = itParam; + } + } + + if (worstLocalChi2 < mImproveCutChi2) { + break; + } + + if (!itWorstParam->isRemovable()) { + shouldRemoveTrack = true; + track.removable(); + break; + } + + auto itNextParam = track.removeParamAtCluster(itWorstParam); + auto itNextToNextParam = (itNextParam == track.end()) ? itNextParam : std::next(itNextParam); + itStartingParam = track.rbegin(); + + if (track.getNClusters() < MinRemovableTrackClusters) { + shouldRemoveTrack = true; + break; + } + while (itNextToNextParam != track.end()) { + if (itNextToNextParam->getClusterPtr()->getChamberId() != itNextParam->getClusterPtr()->getChamberId()) { + itStartingParam = std::make_reverse_iterator(++itNextParam); + break; + } + ++itNextToNextParam; + } + } + + if (!shouldRemoveTrack) { + for (auto& param : track) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) + param.setParameters(param.getSmoothParameters()); + param.setCovariances(param.getSmoothCovariances()); + } + } + + return shouldRemoveTrack; + } + + template + void initCcdb(BC const& bc) + { + if (mRunNumber == bc.runNumber()) { + return; + } + + mRunNumber = bc.runNumber(); + std::map metadata; + auto soreor = o2::ccdb::BasicCCDBManager::getRunDuration(fCCDBApi, mRunNumber); + auto ts = soreor.first; + auto grpmag = fCCDBApi.retrieveFromTFileAny(configCcdb.cfgGrpMagPath, metadata, ts); + o2::base::Propagator::initFieldFromGRP(grpmag); + LOGF(info, "Set field for muons"); + VarManager::SetupMuonMagField(); + if (!o2::base::GeometryManager::isGeometryLoaded()) { + ccdbManager->get(configCcdb.cfgGeoPath); + } + mch::TrackExtrap::setField(); + mch::TrackExtrap::useExtrapV2(); + + // Load geometry information from CCDB/local + LOGF(info, "Loading reference aligned geometry from CCDB no later than %d", configMchRealign.cfgCcdbNoLaterThanRef.value); + ccdbManager->setCreatedNotAfter(configMchRealign.cfgCcdbNoLaterThanRef.value); // this timestamp has to be consistent with what has been used in reco + geoRef = ccdbManager->getForTimeStamp(configMchRealign.cfgGeoRefPath, bc.timestamp()); + ccdbManager->clearCache(configMchRealign.cfgGeoRefPath); + if (geoRef != nullptr) { + transformation = mch::geo::transformationFromTGeoManager(*geoRef); + } else { + LOGF(fatal, "Reference aligned geometry object is not available in CCDB at timestamp=%llu", bc.timestamp()); + } + for (int i = 0; i < NMchDetElems; i++) { + int iDEN = getDetElemId(i); + transformRef[iDEN] = transformation(iDEN); + } + + LOGF(info, "Loading new aligned geometry from CCDB no later than %d", configMchRealign.cfgCcdbNoLaterThanNew.value); + ccdbManager->setCreatedNotAfter(configMchRealign.cfgCcdbNoLaterThanNew.value); // make sure this timestamp can be resolved regarding the reference one + geoNew = ccdbManager->getForTimeStamp(configMchRealign.cfgGeoNewPath, bc.timestamp()); + ccdbManager->clearCache(configMchRealign.cfgGeoNewPath); + if (geoNew != nullptr) { + transformation = mch::geo::transformationFromTGeoManager(*geoNew); + } else { + LOGF(fatal, "New aligned geometry object is not available in CCDB at timestamp=%llu", bc.timestamp()); + } + for (int i = 0; i < NMchDetElems; i++) { + int iDEN = getDetElemId(i); + transformNew[iDEN] = transformation(iDEN); + } + + // Init magnetic field for MFT track extrapolation + auto* fieldB = dynamic_cast(TGeoGlobalMagField::Instance()->GetField()); + if (fieldB) { + std::array centerMft{0, 0, -61.4}; // Field at center of MFT + mBzAtMftCenter = fieldB->getBz(centerMft.data()); + // std::cout << "fieldB: " << (void*)fieldB << std::endl; + } + } + + void initMatchingFunctions() + { + using SVector2 = ROOT::Math::SVector; + using SVector4 = ROOT::Math::SVector; + using SVector5 = ROOT::Math::SVector; + + using SMatrix44 = ROOT::Math::SMatrix; + using SMatrix45 = ROOT::Math::SMatrix; + using SMatrix22 = ROOT::Math::SMatrix; + using SMatrix25 = ROOT::Math::SMatrix; + + // Define built-in matching functions + //________________________________________________________________________________ + mMatchingFunctionMap["matchALL"] = [](const o2::track::TrackParCovFwd& mchTrack, const o2::track::TrackParCovFwd& mftTrack) -> std::tuple { + // Match two tracks evaluating all parameters: X,Y, phi, tanl & q/pt + + SMatrix55Sym hK, vK; + SVector5 mK(mftTrack.getX(), mftTrack.getY(), mftTrack.getPhi(), + mftTrack.getTanl(), mftTrack.getInvQPt()), + rKKminus1; + const auto& globalMuonTrackParameters = mchTrack.getParameters(); + const auto& globalMuonTrackCovariances = mchTrack.getCovariances(); + vK(0, 0) = mftTrack.getCovariances()(0, 0); + vK(1, 1) = mftTrack.getCovariances()(1, 1); + vK(2, 2) = mftTrack.getCovariances()(2, 2); + vK(3, 3) = mftTrack.getCovariances()(3, 3); + vK(4, 4) = mftTrack.getCovariances()(4, 4); + hK(0, 0) = 1.0; + hK(1, 1) = 1.0; + hK(2, 2) = 1.0; + hK(3, 3) = 1.0; + hK(4, 4) = 1.0; + + // Covariance of residuals + SMatrix55Std invResCov = (vK + ROOT::Math::Similarity(hK, globalMuonTrackCovariances)); + invResCov.Invert(); + + // Update Parameters + rKKminus1 = mK - hK * globalMuonTrackParameters; // Residuals of prediction + + auto matchChi2Track = ROOT::Math::Similarity(rKKminus1, invResCov); + + // return chi2 and NDF + return {matchChi2Track, 5}; + }; + + //________________________________________________________________________________ + mMatchingFunctionMap["matchXYPhiTanl"] = [](const o2::track::TrackParCovFwd& mchTrack, const o2::track::TrackParCovFwd& mftTrack) -> std::tuple { + // Match two tracks evaluating positions & angles + + SMatrix45 hK; + SMatrix44 vK; + SVector4 mK(mftTrack.getX(), mftTrack.getY(), mftTrack.getPhi(), + mftTrack.getTanl()), + rKKminus1; + const auto& globalMuonTrackParameters = mchTrack.getParameters(); + const auto& globalMuonTrackCovariances = mchTrack.getCovariances(); + vK(0, 0) = mftTrack.getCovariances()(0, 0); + vK(1, 1) = mftTrack.getCovariances()(1, 1); + vK(2, 2) = mftTrack.getCovariances()(2, 2); + vK(3, 3) = mftTrack.getCovariances()(3, 3); + hK(0, 0) = 1.0; + hK(1, 1) = 1.0; + hK(2, 2) = 1.0; + hK(3, 3) = 1.0; + + // Covariance of residuals + SMatrix44 invResCov = (vK + ROOT::Math::Similarity(hK, globalMuonTrackCovariances)); + invResCov.Invert(); + + // Residuals of prediction + rKKminus1 = mK - hK * globalMuonTrackParameters; + + auto matchChi2Track = ROOT::Math::Similarity(rKKminus1, invResCov); + + // return chi2 and NDF + return {matchChi2Track, 4}; + }; + + //________________________________________________________________________________ + mMatchingFunctionMap["matchXY"] = [](const o2::track::TrackParCovFwd& mchTrack, const o2::track::TrackParCovFwd& mftTrack) -> std::tuple { + // Calculate Matching Chi2 - X and Y positions + + SMatrix25 hK; + SMatrix22 vK; + SVector2 mK(mftTrack.getX(), mftTrack.getY()), rKKminus1; + const auto& globalMuonTrackParameters = mchTrack.getParameters(); + const auto& globalMuonTrackCovariances = mchTrack.getCovariances(); + vK(0, 0) = mftTrack.getCovariances()(0, 0); + vK(1, 1) = mftTrack.getCovariances()(1, 1); + hK(0, 0) = 1.0; + hK(1, 1) = 1.0; + + // Covariance of residuals + SMatrix22 invResCov = (vK + ROOT::Math::Similarity(hK, globalMuonTrackCovariances)); + invResCov.Invert(); + + // Residuals of prediction + rKKminus1 = mK - hK * globalMuonTrackParameters; + auto matchChi2Track = ROOT::Math::Similarity(rKKminus1, invResCov); + + // return reduced chi2 + return {matchChi2Track, 2}; + }; + } + + void init(o2::framework::InitContext&) + { + // Load geometry + ccdbManager->setURL(configCcdb.cfgCcdbUrl); + ccdbManager->setCaching(true); + ccdbManager->setLocalObjectValidityChecking(); + fCCDBApi.init(configCcdb.cfgCcdbUrl); + mRunNumber = 0; + + // Configuration for track fitter + const auto& trackerParam = mch::TrackerParam::Instance(); + trackFitter.setBendingVertexDispersion(trackerParam.bendingVertexDispersion); + trackFitter.setChamberResolution(configMchRealign.cfgChamberResolutionX.value, configMchRealign.cfgChamberResolutionY.value); + trackFitter.smoothTracks(true); + trackFitter.useChamberResolution(); + mImproveCutChi2 = 2. * configMchRealign.cfgSigmaCutImprove.value * configMchRealign.cfgSigmaCutImprove.value; + + // Reset matching configuration, then populate only what we need. + hasActiveChi2Matching = false; + activeChi2FunctionName.clear(); + activeChi2MatchingPlaneZ = 0.; + + hasActiveMlMatching = false; + activeMlMatchingPlaneZ = 0.; + + if (configMatching.cfgCustomMatchingStrategy.value == 0) { + // Matching functions (custom chi2) + initMatchingFunctions(); + auto label = configChi2MatchingOptions.cfgChi2FunctionLabel.value; + auto funcName = configChi2MatchingOptions.cfgChi2FunctionName.value; + auto matchingPlaneZ = configChi2MatchingOptions.cfgChi2FunctionMatchingPlaneZ.value; + + if (!label.empty() && !funcName.empty()) { + hasActiveChi2Matching = true; + activeChi2FunctionName = funcName; + activeChi2MatchingPlaneZ = matchingPlaneZ; + } + } else { + // Matching ML models (custom ML) + // TODO : for now we use hard coded values since the current models use 1 pT bin + binsPtMl = {-1e-6, 1000.0}; + cutValues = {0.0}; + cutDirMl = {cuts_ml::CutNot}; + LabeledArray mycutsMl(cutValues.data(), 1, 1, std::vector{"pT bin 0"}, std::vector{"score"}); + + auto label = configMlOptions.cfgMlModelLabel.value; + auto modelPath = configMlOptions.cfgMlModelPathCcdb.value; + auto inputFeatures = configMlOptions.cfgMlInputFeatures.value; + auto modelName = configMlOptions.cfgMlModelName.value; + auto matchingPlaneZ = configMlOptions.cfgMlModelMatchingPlaneZ.value; + + if (!label.empty() && !modelPath.empty() && !inputFeatures.empty() && !modelName.empty()) { + activeMlResponse.configure(binsPtMl, mycutsMl, cutDirMl, 1); + activeMlResponse.setModelPathsCCDB(std::vector{modelName}, fCCDBApi, std::vector{modelPath}, configCcdb.cfgCcdbNoLaterThan.value); + activeMlResponse.cacheInputFeaturesIndices(inputFeatures); + activeMlResponse.init(); + + hasActiveMlMatching = true; + activeMlMatchingPlaneZ = matchingPlaneZ; + } + } + } + + template + bool pDcaCut(const T& mchTrack, const C& collision, double nSigmaPDCA) + { + static const double sigmaPDCA23 = 80.; + static const double sigmaPDCA310 = 54.; + static const double relPRes = 0.0004; + static const double slopeRes = 0.0005; + + constexpr double AbsorberEndZ = 505.; + constexpr double RadToDeg = 180. / o2::constants::math::PI; + double thetaAbs = std::atan(mchTrack.rAtAbsorberEnd() / AbsorberEndZ) * RadToDeg; + + // propagate muon track to vertex + auto mchTrackAtVertex = VarManager::PropagateMuon(mchTrack, collision, VarManager::kToVertex); + + // double pUncorr = mchTrack.p(); + double p = mchTrackAtVertex.getP(); + + double pDCA = mchTrack.pDca(); + double sigmaPDCA = (thetaAbs < ThetaAbsBoundaryDeg) ? sigmaPDCA23 : sigmaPDCA310; + double nrp = nSigmaPDCA * relPRes * p; + double pResEffect = sigmaPDCA / (1. - nrp / (1. + nrp)); + double slopeResEffect = SlopeResolutionZ * slopeRes * p; + double sigmaPDCAWithRes = std::sqrt(pResEffect * pResEffect + slopeResEffect * slopeResEffect); + return pDCA <= nSigmaPDCA * sigmaPDCAWithRes; + } + + template + bool isGoodMuon(const T& mchTrack, const C& collision, + double chi2Cut, + double pCut, + double pTCut, + std::array etaCut, + std::array rAbsCut, + double nSigmaPdcaCut) + { + // chi2 cut + if (mchTrack.chi2() > chi2Cut) { + return false; + } + + // momentum cut + if (mchTrack.p() < pCut) { + return false; // skip low-momentum tracks + } + + // transverse momentum cut + if (mchTrack.pt() < pTCut) { + return false; // skip low-momentum tracks + } + + // Eta cut + double eta = mchTrack.eta(); + if ((eta < etaCut[0] || eta > etaCut[1])) { + return false; + } + + // RAbs cut + double rAbs = mchTrack.rAtAbsorberEnd(); + if ((rAbs < rAbsCut[0] || rAbs > rAbsCut[1])) { + return false; + } + + // pDCA cut + return pDcaCut(mchTrack, collision, nSigmaPdcaCut); + } + + void storeFwdTrackCovariance(const SMatrix55Sym& cov) + { + const float sigX = std::sqrt(cov(0, 0)); + const float sigY = std::sqrt(cov(1, 1)); + const float sigPhi = std::sqrt(cov(2, 2)); + const float sigTgl = std::sqrt(cov(3, 3)); + const float sig1Pt = std::sqrt(cov(4, 4)); + const auto rhoXY = static_cast(128.f * cov(0, 1) / (sigX * sigY)); + const auto rhoPhiX = static_cast(128.f * cov(0, 2) / (sigPhi * sigX)); + const auto rhoPhiY = static_cast(128.f * cov(1, 2) / (sigPhi * sigY)); + const auto rhoTglX = static_cast(128.f * cov(0, 3) / (sigTgl * sigX)); + const auto rhoTglY = static_cast(128.f * cov(1, 3) / (sigTgl * sigY)); + const auto rhoTglPhi = static_cast(128.f * cov(2, 3) / (sigTgl * sigPhi)); + const auto rho1PtX = static_cast(128.f * cov(0, 4) / (sig1Pt * sigX)); + const auto rho1PtY = static_cast(128.f * cov(1, 4) / (sig1Pt * sigY)); + const auto rho1PtPhi = static_cast(128.f * cov(2, 4) / (sig1Pt * sigPhi)); + const auto rho1PtTgl = static_cast(128.f * cov(3, 4) / (sig1Pt * sigTgl)); + gmCandidateFwdTracksCov(sigX, sigY, sigPhi, sigTgl, sig1Pt, + rhoXY, rhoPhiY, rhoPhiX, rhoTglX, rhoTglY, rhoTglPhi, rho1PtX, rho1PtY, rho1PtPhi, rho1PtTgl); + } + + template + void fillBaseGmmCandFwdTrack(TMCH const& track, + TrackParExt const& trackPar, + int32_t gmmMchTrackId, + float chi2MatchMCHMFT, + float matchScoreMCHMFT) + { + const auto collisionId = track.collisionId(); + bool hasBcSlice = false; + std::array bcSlice{}; + if (collisionId < 0) { + const auto ambIt = mAmbBcSliceByFwdTrackId.find(track.globalIndex()); + if (ambIt != mAmbBcSliceByFwdTrackId.end()) { + bcSlice = ambIt->second; + hasBcSlice = true; + } + } + + gmCandidateFwdTracks( + collisionId, + track.trackType(), + trackPar.getX(), + trackPar.getY(), + trackPar.getZ(), + trackPar.getPhi(), + trackPar.getTgl(), + trackPar.getInvQPt(), + trackPar.getNClusters(), + track.pDca(), + track.rAtAbsorberEnd(), + trackPar.isRemovable(), + trackPar.getTrackChi2(), + track.chi2MatchMCHMID(), + chi2MatchMCHMFT, + matchScoreMCHMFT, + track.matchMFTTrackId(), + gmmMchTrackId, + track.mchBitMap(), + track.midBitMap(), + track.midBoards(), + track.trackTime(), + track.trackTimeRes()); + + storeFwdTrackCovariance(trackPar.getCovariances()); + if (hasBcSlice) { + gmAmbiguousFwdTracksReAlign(mGmmCandFwdTrackRowIndex, bcSlice.data()); + } + mGmmCandFwdTrackRowIndex += 1; + + mHasLastMchAmbiguousBcSlice = hasBcSlice; + if (hasBcSlice) { + mLastMchAmbiguousBcSlice = bcSlice; + } + } + + template + void fillCandidateFwdTrack(TMCH const& mchTrack, + TrackParExt const& mchPar, + int32_t gmmMchTrackId, + TMFT const& mftTrack, + TrackParExt const& mftPar, + const MatchingCandidate& candidate) + { + using o2::aod::fwdtrack::ForwardTrackTypeEnum; + using o2::aod::fwdtrackutils::propagationPoint; + + constexpr uint8_t CandidateTrackType = static_cast(ForwardTrackTypeEnum::GlobalForwardTrack); + + auto propmuonAtMft = fwdToMch(mchPar); + o2::mch::TrackExtrap::extrapToVertex(propmuonAtMft, + mftPar.getX(), + mftPar.getY(), + mftPar.getZ(), + mftPar.getSigma2X(), + mftPar.getSigma2Y()); + + const auto globalMuonRefit = o2::aod::fwdtrackutils::refitGlobalMuonCov(mchToFwd(propmuonAtMft), mftPar); + + const auto nClusters = static_cast(std::min(127, mchPar.getNClusters() + mftPar.getNClusters())); + + const float chi2 = static_cast(mchTrack.chi2()); + const int32_t collisionId = mchTrack.has_collision() ? mchTrack.collisionId() : -1; + bool hasBcSlice = false; + std::array bcSlice{}; + if (collisionId < 0) { + if (mHasLastMchAmbiguousBcSlice) { + bcSlice = mLastMchAmbiguousBcSlice; + hasBcSlice = true; + } else { + const auto ambIt = mAmbBcSliceByFwdTrackId.find(mchTrack.globalIndex()); + if (ambIt != mAmbBcSliceByFwdTrackId.end()) { + bcSlice = ambIt->second; + hasBcSlice = true; + } + } + } + + bool isRemovable = mchPar.isRemovable(); + + gmCandidateFwdTracks( + collisionId, + CandidateTrackType, + globalMuonRefit.getX(), + globalMuonRefit.getY(), + globalMuonRefit.getZ(), + globalMuonRefit.getPhi(), + globalMuonRefit.getTgl(), + globalMuonRefit.getInvQPt(), + nClusters, + mchTrack.pDca(), + mchTrack.rAtAbsorberEnd(), + isRemovable, + chi2, + mchTrack.chi2MatchMCHMID(), + static_cast(candidate.matchChi2), + static_cast(candidate.matchScore), + static_cast(mftTrack.globalIndex()), + gmmMchTrackId, + mchTrack.mchBitMap(), + mchTrack.midBitMap(), + mchTrack.midBoards(), + mchTrack.trackTime(), + mchTrack.trackTimeRes()); + + storeFwdTrackCovariance(globalMuonRefit.getCovariances()); + if (hasBcSlice) { + gmAmbiguousFwdTracksReAlign(mGmmCandFwdTrackRowIndex, bcSlice.data()); + } + mGmmCandFwdTrackRowIndex += 1; + } + + o2::track::TrackParCovFwd propagateToZMch(const o2::track::TrackParCovFwd& muon, const double z) + { + auto mchTrack = fwdToMch(muon); + + float absFront = -90.f; + float absBack = -505.f; + + if (muon.getZ() < absBack && z > absFront) { + // extrapolation through the absorber in the upstream direction + o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchTrack, z); + } else { + // all other cases + o2::mch::TrackExtrap::extrapToZCov(mchTrack, z); + } + + return mchToFwd(mchTrack); + } + + o2::track::TrackParCovFwd propagateToZMft(const o2::track::TrackParCovFwd& mftTrack, const double z) + { + o2::track::TrackParCovFwd trackExtrap{mftTrack}; + trackExtrap.propagateToZ(z, mBzAtMftCenter); + return trackExtrap; + } + + template + o2::track::TrackParCovFwd propagateToVertexMch(const TMCH& muon, + const C& collision) + { + auto mchTrack = fwdToMch(fwdtrackutils::getTrackParCovFwd(muon, muon)); + o2::mch::TrackExtrap::extrapToVertex(mchTrack, + collision.posX(), + collision.posY(), + collision.posZ(), + collision.covXX(), + collision.covYY()); + return mchToFwd(mchTrack); + } + + // tag muons based on the track quality and the track position at the front and back MFT planes + template + void getTaggedMuons(C const& collisions, + TMUON const& muonTracks, + std::vector& taggedMuons) + { + taggedMuons.clear(); + for (const auto& muonTrack : muonTracks) { + + // only consider MCH-MID matches + if (static_cast(muonTrack.trackType()) != MchMidTrackType) { + continue; + } + + // only select MCH-MID tracks associated to a collision + if (!muonTrack.has_collision()) { + continue; + } + + const auto& collision = collisions.rawIteratorAt(muonTrack.collisionId()); + + // select MCH tracks with strict quality cuts + if (!isGoodMuon(muonTrack, collision, + configMuonTagging.cfgMuonTaggingTrackChi2MchUp, + configMuonTagging.cfgMuonTaggingPMchLow, + configMuonTagging.cfgMuonTaggingPtMchLow, + {configMuonTagging.cfgMuonTaggingEtaMchLow, configMuonTagging.cfgMuonTaggingEtaMchUp}, + {configMuonTagging.cfgMuonTaggingRabsLow, configMuonTagging.cfgMuonTaggingRabsUp}, + configMuonTagging.cfgMuonTaggingPdcaUp)) { + continue; + } + + // propagate MCH track to the vertex + auto mchTrackAtVertex = propagateToVertexMch(muonTrack, collision); + + // propagate the track from the vertex to the first MFT plane + const auto& extrapToMFTfirst = propagateToZMch(mchTrackAtVertex, o2::mft::constants::mft::LayerZCoordinate()[0]); + double rFront = std::sqrt(extrapToMFTfirst.getX() * extrapToMFTfirst.getX() + extrapToMFTfirst.getY() * extrapToMFTfirst.getY()); + if (rFront < configMuonTagging.cfgMuonTaggingRadiusAtMftFrontLow.value || rFront > configMuonTagging.cfgMuonTaggingRadiusAtMftFrontUp.value) { + continue; + } + + // propagate the track from the vertex to the last MFT plane + const auto& extrapToMFTlast = propagateToZMch(mchTrackAtVertex, o2::mft::constants::mft::LayerZCoordinate()[9]); + double rBack = std::sqrt(extrapToMFTlast.getX() * extrapToMFTlast.getX() + extrapToMFTlast.getY() * extrapToMFTlast.getY()); + if (rBack < configMuonTagging.cfgMuonTaggingRadiusAtMftBackLow.value || rBack > configMuonTagging.cfgMuonTaggingRadiusAtMftBackUp.value) { + continue; + } + + int64_t muonTrackIndex = muonTrack.globalIndex(); + taggedMuons.emplace_back(muonTrackIndex); + } + } + + template + bool isMftMchTimeCompatible(EVT const& collisions, + BC const& bcs, + TMUON const& mchTrack, + TMFT const& mftTrack) + { + if (!mchTrack.has_collision() || !mftTrack.has_collision()) { + return false; + } + + const auto& collMch = collisions.rawIteratorAt(mchTrack.collisionId()); + const auto& bcMch = bcs.rawIteratorAt(collMch.bcId()); + const auto& collMft = collisions.rawIteratorAt(mftTrack.collisionId()); + const auto& bcMft = bcs.rawIteratorAt(collMft.bcId()); + + int64_t deltaBc = static_cast(bcMft.globalBC()) - static_cast(bcMch.globalBC()); + double deltaBcNS = o2::constants::lhc::LHCBunchSpacingNS * deltaBc; + double deltaTrackTime = mftTrack.trackTime() - mchTrack.trackTime() + deltaBcNS; + double trackTimeResTot = mftTrack.trackTimeRes() + mchTrack.trackTimeRes(); + + return std::fabs(deltaTrackTime) <= trackTimeResTot; + } + + template + void prepareMatchingCandidates(EVT const& collisions, + BC const& bcs, + TMUON const& muonTracks, + TMFT const& mftTracks, + MyMFTCovariances const& mftCovs) + { + mMftTrackPars.clear(); + mMchTrackPars.clear(); + mMatchingCandidates.clear(); + + LOGF(info, "Filling matching candidate tables"); + + for (const auto& muonTrack : muonTracks) { + if (static_cast(muonTrack.trackType()) <= GlobalTrackTypeMax) { + continue; + } + auto mchTrackIndex = muonTrack.globalIndex(); + + // initialize the MCH track parameters, which will be updated by the realignment if enabled + mMchTrackPars.try_emplace(mchTrackIndex, TrackParExt(fwdtrackutils::getTrackParCovFwd(muonTrack, muonTrack), muonTrack.nClusters())); + } + + for (const auto& mftTrack : mftTracks) { + auto mftTrackIndex = mftTrack.globalIndex(); + + // initialize the MFT track parameters, which will be updated by the alignment corrections if enabled + if (mftTrackCovs.contains(mftTrackIndex) && !mMftTrackPars.contains(mftTrackIndex)) { + auto const& mftTrackCov = mftCovs.rawIteratorAt(mftTrackCovs[mftTrackIndex]); + mMftTrackPars.emplace(mftTrackIndex, TrackParExt(fwdtrackutils::getTrackParCovFwd(mftTrack, mftTrackCov), mftTrack.nClusters())); + } + } + + // fill matching candidates table + if (!configMatching.cfgMatchAllTracks.value) { + // collect global MFT-MCH or MFT-MCH-MID tracks and associate them to the corresponding MCH(-MID) track + for (const auto& muonTrack : muonTracks) { + // skip MCH or MCH-MID tracks + if (static_cast(muonTrack.trackType()) > GlobalTrackTypeMax) { + continue; + } + + auto const& mchTrack = muonTrack.template matchMCHTrack_as(); + int64_t mchTrackIndex = mchTrack.globalIndex(); + auto const& mftTrack = muonTrack.template matchMFTTrack_as(); + int64_t mftTrackIndex = mftTrack.globalIndex(); + + if (!mftTrackCovs.contains(mftTrackIndex)) { + continue; + } + + mMatchingCandidates[mchTrackIndex].emplace_back(MatchingCandidate{ + .muonTrackId = muonTrack.globalIndex(), + .mftTrackId = mftTrackIndex, + .matchScore = muonTrack.matchScoreMCHMFT(), + .matchChi2 = muonTrack.chi2MatchMCHMFT()}); + } + } else { + // build matching candidates from all time-compatible MFT-MCH pairs + for (const auto& muonTrack : muonTracks) { + if (static_cast(muonTrack.trackType()) <= GlobalTrackTypeMax) { + continue; + } + auto mchTrackIndex = muonTrack.globalIndex(); + for (const auto& mftTrack : mftTracks) { + if (!isMftMchTimeCompatible(collisions, bcs, muonTrack, mftTrack)) { + continue; + } + if (!mftTrackCovs.contains(mftTrack.globalIndex())) { + continue; + } + + mMatchingCandidates[mchTrackIndex].emplace_back(MatchingCandidate{ + .mftTrackId = mftTrack.globalIndex()}); + } + } + } + + // sort the vectors of matching candidates in ascending order based on the matching chi2 value + auto compareMatchingChi2 = [](const MatchingCandidate& track1, const MatchingCandidate& track2) -> bool { + return (track1.matchChi2 < track2.matchChi2); + }; + + for (auto& [mchIndex, candidatesVector] : mMatchingCandidates) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) + std::sort(candidatesVector.begin(), candidatesVector.end(), compareMatchingChi2); + } + } + + template + o2::track::TrackParCovFwd transformMft(TMFT& mftTrack, TMFTCOV const& mftTrackCov) + { + auto track = fwdToMch(fwdtrackutils::getTrackParCovFwd(mftTrack, mftTrackCov)); + + double z = track.getZ(); + // double dZ = zMCH - z; + double x = track.getNonBendingCoor(); + double y = track.getBendingCoor(); + double xSlope = track.getNonBendingSlope(); + double ySlope = track.getBendingSlope(); + + double xSlopeCorrection = (y > 0) ? configMftAlignmentCorrections.cfgMFTAlignmentCorrXSlopeTop : configMftAlignmentCorrections.cfgMFTAlignmentCorrXSlopeBottom; + double xCorrection = xSlopeCorrection * z + + ((y > 0) ? configMftAlignmentCorrections.cfgMFTAlignmentCorrXOffsetTop : configMftAlignmentCorrections.cfgMFTAlignmentCorrXOffsetBottom); + double xNew = x + xCorrection; + double xSlopeNew = xSlope + xSlopeCorrection; + + track.setNonBendingCoor(xNew); + track.setNonBendingSlope(xSlopeNew); + + double ySlopeCorrection = (y > 0) ? configMftAlignmentCorrections.cfgMFTAlignmentCorrYSlopeTop : configMftAlignmentCorrections.cfgMFTAlignmentCorrYSlopeBottom; + double yCorrection = ySlopeCorrection * z + + ((y > 0) ? configMftAlignmentCorrections.cfgMFTAlignmentCorrYOffsetTop : configMftAlignmentCorrections.cfgMFTAlignmentCorrYOffsetBottom); + track.setBendingCoor(y + yCorrection); + track.setBendingSlope(ySlope + ySlopeCorrection); + + return mchToFwd(track); + } + + template + void runMftRealignment(TMFTs const& mftTracks, TMFTCOVs const& mftCovs) + { + for (const auto& mftTrack : mftTracks) { + auto mftTrackIndex = mftTrack.globalIndex(); + if (!mftTrackCovs.contains(mftTrackIndex)) { + continue; + } + + auto const& mftTrackCov = mftCovs.rawIteratorAt(mftTrackCovs[mftTrackIndex]); + mMftTrackPars[mftTrackIndex] = transformMft(mftTrack, mftTrackCov); + } + } + + template + void runMuonRealignment(TMuons const& muons, TMuonCls const& clusters) + { + // Loop over forward tracks + for (auto const& muon : muons) { + int mchIndex = muon.globalIndex(); + // skip global forward matches + if (muon.trackType() > GlobalTrackTypeMax) { + continue; + } + + // continue; + + auto mchTrackParIt = mMchTrackPars.find(mchIndex); + if (mchTrackParIt == mMchTrackPars.end()) { + continue; + } + + auto clustersSliced = clusters.sliceBy(perMuon, muon.globalIndex()); // Slice clusters by muon id + mch::Track convertedTrack = mch::Track(); // Temporary variable to store re-aligned clusters + + int clIndex = -1; + // Get re-aligned clusters associated to current track + for (auto const& cluster : clustersSliced) { + clIndex += 1; + + auto* clusterMCH = new mch::Cluster(); + + math_utils::Point3D local; + math_utils::Point3D master; + master.SetXYZ(cluster.x(), cluster.y(), cluster.z()); + + // Transformation from reference geometry frame to new geometry frame + transformRef[cluster.deId()].MasterToLocal(master, local); + transformNew[cluster.deId()].LocalToMaster(local, master); + + clusterMCH->x = master.x(); + clusterMCH->y = master.y(); + clusterMCH->z = master.z(); + + const 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); + // LOGF(debug, "Track %d, cluster DE%d: x:%g y:%g z:%g", muon.globalIndex(), cluster.deId(), cluster.x(), cluster.y(), cluster.z()); + // LOGF(debug, "Track %d, re-aligned cluster DE%d: x:%g y:%g z:%g", muonRealignId, cluster.deId(), clusterMCH->getX(), clusterMCH->getY(), clusterMCH->getZ()); + } + + // Refit the re-aligned track + int removable = 0; + if (convertedTrack.getNClusters() != 0) { + removable = removeTrack(convertedTrack); + } else { + LOGF(fatal, "Muon track %d has no associated clusters.", muon.globalIndex()); + } + + // Get the re-aligned track parameter: track param at the first cluster + mch::TrackParam trackParam = mch::TrackParam(convertedTrack.first()); + + // Convert MCH track to FWD track and store new parameters after realignment + mchTrackParIt->second = mchToFwd(mch::TrackParam(convertedTrack.first())); + mchTrackParIt->second.setTrackChi2(trackParam.getTrackChi2() / convertedTrack.getNDF()); + mchTrackParIt->second.setNClusters(convertedTrack.getNClusters()); + if (removable) { + mchTrackParIt->second.setRemovable(); + } + } + } + + void runChi2Matching(const std::string& funcName, + float matchingPlaneZ, + const MatchingCandidates& matchingCandidates, + MatchingCandidates& newMatchingCandidates) + { + newMatchingCandidates.clear(); + + std::string funcNameEffective = funcName; + float matchingPlaneZEffective = matchingPlaneZ; + if (funcName == "prod") { + funcNameEffective = "matchALL"; + matchingPlaneZEffective = MatchingPlaneDefaultZ; + } + + if (!mMatchingFunctionMap.contains(funcNameEffective)) { + return; + } + auto matchingFunc = mMatchingFunctionMap.at(funcNameEffective); + + for (const auto& [mchIndex, candidatesVector] : matchingCandidates) { + + // get the tracks parameters, which have been updated by the realignment if enabled + const auto mchTrackParIt = mMchTrackPars.find(mchIndex); + if (mchTrackParIt == mMchTrackPars.end()) { + continue; + } + + for (const auto& candidate : candidatesVector) { + auto mftTrackParIt = mMftTrackPars.find(candidate.mftTrackId); + if (mftTrackParIt == mMftTrackPars.end()) { + continue; + } + + auto mftTrackProp = mftTrackParIt->second.asTrackParCovFwd(); + auto mchTrackProp = mchTrackParIt->second.asTrackParCovFwd(); + + if (matchingPlaneZEffective < 0.) { + mftTrackProp = propagateToZMft(mftTrackProp, matchingPlaneZ); + mchTrackProp = propagateToZMch(mchTrackProp, matchingPlaneZ); + } + + auto matchResult = matchingFunc(mchTrackProp, mftTrackProp); + float matchChi2 = std::get<0>(matchResult); + + newMatchingCandidates[mchIndex].emplace_back(MatchingCandidate{ + .muonTrackId = candidate.muonTrackId, + .mftTrackId = candidate.mftTrackId, + .matchScore = -1, + .matchChi2 = matchChi2}); + } + } + + auto compareMatchingChi2 = [](const MatchingCandidate& track1, const MatchingCandidate& track2) -> bool { + return (track1.matchChi2 < track2.matchChi2); + }; + + for (auto& [mchIndex, globalTracksVector] : newMatchingCandidates) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) + std::sort(globalTracksVector.begin(), globalTracksVector.end(), compareMatchingChi2); + + int ranking = 1; + for (auto& candidate : globalTracksVector) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) + candidate.matchRanking = ranking; + ranking += 1; + } + } + } + + template + void runMlMatching(C const& collisions, + TMUON const& muonTracks, + TMFT const& mftTracks, + o2::analysis::MlResponseMFTMuonMatch& mlResponse, + float matchingPlaneZ, + const MatchingCandidates& matchingCandidates, + MatchingCandidates& newMatchingCandidates) + { + newMatchingCandidates.clear(); + for (const auto& [mchIndex, candidatesVector] : matchingCandidates) { + auto const& mchTrack = muonTracks.rawIteratorAt(mchIndex); + if (!mchTrack.has_collision()) { + continue; + } + + auto collision = collisions.rawIteratorAt(mchTrack.collisionId()); + + // get the tracks parameters, which have been updated by the realignment if enabled + auto mchTrackParIt = mMchTrackPars.find(mchIndex); + if (mchTrackParIt == mMchTrackPars.end()) { + continue; + } + + for (const auto& candidate : candidatesVector) { + auto const& muonTrack = (candidate.muonTrackId >= 0) ? muonTracks.rawIteratorAt(candidate.muonTrackId) : mchTrack; + auto const& mftTrack = mftTracks.rawIteratorAt(candidate.mftTrackId); + auto mftTrackParIt = mMftTrackPars.find(candidate.mftTrackId); + if (mftTrackParIt == mMftTrackPars.end()) { + continue; + } + + auto mftTrackProp = mftTrackParIt->second.asTrackParCovFwd(); + auto mchTrackProp = mchTrackParIt->second.asTrackParCovFwd(); + + if (matchingPlaneZ < 0.) { + mftTrackProp = propagateToZMft(mftTrackProp, matchingPlaneZ); + mchTrackProp = propagateToZMch(mchTrackProp, matchingPlaneZ); + } + + std::vector output; + std::vector inputML = mlResponse.getInputFeatures(muonTrack, mftTrack, mchTrack, mftTrackProp, mchTrackProp, collision); + mlResponse.isSelectedMl(inputML, 0, output); + float matchScore = output[0]; + + newMatchingCandidates[mchIndex].emplace_back(MatchingCandidate{ + .muonTrackId = candidate.muonTrackId, + .mftTrackId = candidate.mftTrackId, + .matchScore = matchScore, + .matchChi2 = -1}); + } + } + + auto compareMatchingScore = [](const MatchingCandidate& track1, const MatchingCandidate& track2) -> bool { + return (track1.matchScore > track2.matchScore); + }; + + for (auto& [mchIndex, globalTracksVector] : newMatchingCandidates) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) + std::sort(globalTracksVector.begin(), globalTracksVector.end(), compareMatchingScore); + + int ranking = 1; + for (auto& candidate : globalTracksVector) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) + candidate.matchRanking = ranking; + ranking += 1; + } + } + } + + template + void processMatchingCandidates(C const& collisions, + TMUON const& muonTracks, + TMFT const& mftTracks, + CMFT const& mftCovs, + aod::FwdTrkCls const& clusters) + { + if (configMchRealign.cfgEnableMCHRealign.value) { + runMuonRealignment(muonTracks, clusters); + } + + if (configMftAlignmentCorrections.cfgEnableMftAlignmentCorrections) { + runMftRealignment(mftTracks, mftCovs); + } + + std::vector taggedMuons; + getTaggedMuons(collisions, muonTracks, taggedMuons); + + if (configMatching.cfgCustomMatchingStrategy.value == 0) { + if (hasActiveChi2Matching) { + MatchingCandidates newMatchingCandidates; + runChi2Matching(activeChi2FunctionName, activeChi2MatchingPlaneZ, mMatchingCandidates, newMatchingCandidates); + fillMatchingCandidates(newMatchingCandidates, taggedMuons); + } + } else { + if (hasActiveMlMatching) { + MatchingCandidates newMatchingCandidates; + runMlMatching(collisions, muonTracks, mftTracks, activeMlResponse, activeMlMatchingPlaneZ, mMatchingCandidates, newMatchingCandidates); + fillMatchingCandidates(newMatchingCandidates, taggedMuons); + } + } + } + + void fillMatchingCandidates(const MatchingCandidates& matchingCandidates, + const std::vector& taggedMuons) + { + for (const auto& [mchIndex, candidates] : matchingCandidates) { + if (candidates.empty()) { + continue; + } + + bool isTagged = std::find(taggedMuons.begin(), taggedMuons.end(), mchIndex) != taggedMuons.end(); + + std::vector storedCandidates; + int nStored = 0; + for (const auto& candidate : candidates) { + if (configMatching.cfgMaxCandidatesPerMchTrack.value >= 0 && nStored >= configMatching.cfgMaxCandidatesPerMchTrack.value) { + break; + } + + int32_t candidateIndex = mMatchCandidateCounter; + globalMuonMatchCandidates( + mchIndex, + candidate.mftTrackId, + static_cast(candidate.matchChi2), + static_cast(candidate.matchScore), + static_cast(candidate.matchRanking), + isTagged); + mMatchCandidateCounter += 1; + + mMchTrackToCandidateIndices[mchIndex].push_back(candidateIndex); + storedCandidates.push_back(candidate); + nStored += 1; + } + + if (!storedCandidates.empty()) { + mMchTrackMatchingCandidates[mchIndex] = std::move(storedCandidates); + } + } + } + + int32_t countStoredCandidatesForMchTrack(int64_t mchTrackIndex) const + { + const auto candidateIterator = mMchTrackMatchingCandidates.find(mchTrackIndex); + if (candidateIterator == mMchTrackMatchingCandidates.end()) { + return 0; + } + return static_cast(candidateIterator->second.size()); + } + + template + void fillGmmCandidateFwdTracks(TMUON const& muonTracks, + TMFT const& mftTracks, + aod::AmbiguousFwdTracks const& ambFwdTracks) + { + mFwdTrackToGmmCandTrkIndex.clear(); + mGmmCandFwdTrackRowIndex = 0; + mHasLastMchAmbiguousBcSlice = false; + mAmbBcSliceByFwdTrackId.clear(); + for (const auto& ambFwdTrack : ambFwdTracks) { + const auto bcIds = ambFwdTrack.bcIds(); + mAmbBcSliceByFwdTrackId[ambFwdTrack.fwdtrackId()] = {bcIds[0], bcIds[1]}; + } + + // First pass: assign GMMCANDTRK row indices for MCH/MCH-MID base entries so that + // MCHTrackId can be remapped consistently even when global muons appear first in FwdTracks. + int32_t nextGmmCandTrkIndex = 0; + for (const auto& track : muonTracks) { + const int trackType = static_cast(track.trackType()); + if (trackType > GlobalTrackTypeMax) { + mFwdTrackToGmmCandTrkIndex[track.globalIndex()] = nextGmmCandTrkIndex; + nextGmmCandTrkIndex += 1 + countStoredCandidatesForMchTrack(track.globalIndex()); + } else if (configMatching.cfgIncludeGlobalMuonsInFwdTracks.value) { + nextGmmCandTrkIndex += 1; + } + } + + // Second pass: fill GMMCANDTRK/GMMCANDTRKCOV in FwdTracks order. + for (const auto& track : muonTracks) { + const int trackType = static_cast(track.trackType()); + + if (trackType > GlobalTrackTypeMax) { + mHasLastMchAmbiguousBcSlice = false; + const int64_t mchTrackIndex = track.globalIndex(); + const int32_t gmmMchTrackId = mFwdTrackToGmmCandTrkIndex.at(mchTrackIndex); + + const auto candidateIterator = mMchTrackMatchingCandidates.find(mchTrackIndex); + auto mchTrackParIt = mMchTrackPars.find(mchTrackIndex); + if (mchTrackParIt == mMchTrackPars.end()) { + // fill muon tracks table with original parameters + const TrackParExt trackPar{fwdtrackutils::getTrackParCovFwd(track, track)}; + fillBaseGmmCandFwdTrack(track, trackPar, gmmMchTrackId, -1.f, -1.f); + } else { + // fill muon tracks table with realignment parameters + fillBaseGmmCandFwdTrack(track, mchTrackParIt->second, gmmMchTrackId, -1.f, -1.f); + } + + if (candidateIterator != mMchTrackMatchingCandidates.end()) { + for (const auto& candidate : candidateIterator->second) { + auto mftTrackParIt = mMftTrackPars.find(candidate.mftTrackId); + if (mftTrackParIt != mMftTrackPars.end()) { + const auto& mftTrack = mftTracks.rawIteratorAt(candidate.mftTrackId); + fillCandidateFwdTrack(track, mchTrackParIt->second, gmmMchTrackId, mftTrack, mftTrackParIt->second, candidate); + } + } + } + } + + if (configMatching.cfgIncludeGlobalMuonsInFwdTracks.value && trackType <= GlobalTrackTypeMax) { + int32_t gmmMchTrackId = -1; + const auto mchIterator = mFwdTrackToGmmCandTrkIndex.find(track.matchMCHTrackId()); + if (mchIterator != mFwdTrackToGmmCandTrkIndex.end()) { + gmmMchTrackId = mchIterator->second; + } + TrackParExt parExt(fwdtrackutils::getTrackParCovFwd(track, track)); + fillBaseGmmCandFwdTrack(track, + parExt, + gmmMchTrackId, + track.chi2MatchMCHMFT(), + track.matchScoreMCHMFT()); + } + } + } + + template + void fillFwdTrkMatchCands(TMUON const& muonTracks) + { + std::vector empty{}; + for (const auto& muonTrack : muonTracks) { + if (static_cast(muonTrack.trackType()) <= GlobalTrackTypeMax) { + fwdTrkMatchCands(empty); + continue; + } + + const int64_t mchTrackIndex = muonTrack.globalIndex(); + const auto matchIterator = mMchTrackToCandidateIndices.find(mchTrackIndex); + if (matchIterator == mMchTrackToCandidateIndices.end() || matchIterator->second.empty()) { + fwdTrkMatchCands(empty); + } else { + fwdTrkMatchCands(matchIterator->second); + } + } + } + + void processData(MyEvents const& collisions, + aod::BCsWithTimestamps const& bcs, + MyMuons const& muonTracks, + MyMFTs const& mftTracks, + MyMFTCovariances const& mftCovs, + aod::FwdTrkCls const& clusters, + aod::AmbiguousFwdTracks const& ambFwdTracks) + { + auto bc = bcs.begin(); + initCcdb(bc); + + LOGF(info, "Filling MFT cov"); + mftTrackCovs.clear(); + for (const auto& mftTrackCov : mftCovs) { + mftTrackCovs[mftTrackCov.matchMFTTrackId()] = mftTrackCov.globalIndex(); + } + + mMatchCandidateCounter = 0; + mMchTrackToCandidateIndices.clear(); + mMchTrackMatchingCandidates.clear(); + mFwdTrackToGmmCandTrkIndex.clear(); + + LOGF(info, "Preparing candidates"); + prepareMatchingCandidates(collisions, bcs, muonTracks, mftTracks, mftCovs); + + LOGF(info, "Processing candidates"); + processMatchingCandidates(collisions, muonTracks, mftTracks, mftCovs, clusters); + + LOGF(info, "Filling tables"); + // fill table with track/candidates index mapping + fillFwdTrkMatchCands(muonTracks); + // fill track tables + fillGmmCandidateFwdTracks(muonTracks, mftTracks, ambFwdTracks); + } + + PROCESS_SWITCH(GlobalMuonMatching, processData, "processData", true); +}; + +// Extends the fwdtracksrealign table with expression columns +struct GlobalMuonMatchingSpawner { + Spawns realignFwdTrksCov; + Spawns realignFwdTrks; + void init(InitContext const&) {} +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc)}; +}; From 5dcad1e61a83eb833810bc7b1e74c8052e394ebe Mon Sep 17 00:00:00 2001 From: Maurice Coquet Date: Sat, 12 Sep 2026 20:45:16 +0200 Subject: [PATCH 5/9] logic fix --- PWGDQ/Core/VarManager.h | 7 +-- .../TableProducer/tableMakerMC_withAssoc.cxx | 48 ++++++++++------ PWGDQ/TableProducer/tableMaker_withAssoc.cxx | 55 ++++++++++++------- 3 files changed, 68 insertions(+), 42 deletions(-) diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index 9b1642287e8..004c7e618d4 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -3449,7 +3449,7 @@ void VarManager::FillTrackCollision(T const& track, C const& collision, float* v float dcaX = 999.f; float dcaY = 999.f; if (static_cast(track.trackType()) <= 2) { - // Global muons: helix DCA w.r.t. the associated collision + // Global / MCH-MID: helix DCA only (for globals, kMuonPDca is filled from the matched MCH in skimMuons) float xShift = 0.f; float yShift = 0.f; float zShift = 0.f; @@ -3460,17 +3460,16 @@ void VarManager::FillTrackCollision(T const& track, C const& collision, float* v dcaX = static_cast(dca[0]); dcaY = static_cast(dca[1]); } else { - // MCH / standalone: DCA from MCH extrapolation + // MCH standalone: DCA and pDCA from MCH extrapolation o2::dataformats::GlobalFwdTrack propmuonAtDCA = PropagateMuon(track, collision, kToDCA); dcaX = propmuonAtDCA.getX() - collision.posX(); dcaY = propmuonAtDCA.getY() - collision.posY(); float dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY); values[kMuonPDca] = track.p() * dcaXY; } - + values[kMuonDCAx] = dcaX; values[kMuonDCAy] = dcaY; - } } diff --git a/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx b/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx index 9e005ee3b9f..2396d797536 100644 --- a/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx +++ b/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx @@ -1134,21 +1134,28 @@ struct TableMakerMC { // NOTE: If a muon is associated to multiple collisions, depending on the selections, // it may be accepted for some associations and rejected for other if (static_cast(muon.trackType()) > 2 && fConfigVariousOptions.fPropMuon) { - VarManager::FillPropagateMuon(muon, collision); + VarManager::FillPropagateMuon(muon, collision); } // recalculate pDca / DCA and global muon kinematics - if (static_cast(muon.trackType()) < 2 && fConfigVariousOptions.fRefitGlobalMuon) { + // kMuonPDca is always taken from MCH (standalone or the MCH matched to a global) + if (static_cast(muon.trackType()) < 2) { auto muontrack = muon.template matchMCHTrack_as(); - if (muontrack.eta() < fConfigVariousOptions.fMuonMatchEtaMin || muontrack.eta() > fConfigVariousOptions.fMuonMatchEtaMax) { - continue; - } - auto mfttrack = muon.template matchMFTTrack_as(); VarManager::FillTrackCollision(muontrack, collision); - if constexpr (static_cast(TMFTFillMap & VarManager::ObjTypes::MFTCov)) { - auto const& mfttrackcov = mfCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]); - VarManager::FillGlobalMuonRefitCov(muontrack, mfttrack, collision, mfttrackcov); + if (fConfigVariousOptions.fRefitGlobalMuon) { + if (muontrack.eta() < fConfigVariousOptions.fMuonMatchEtaMin || muontrack.eta() > fConfigVariousOptions.fMuonMatchEtaMax) { + continue; + } + auto mfttrack = muon.template matchMFTTrack_as(); + // Helix DCA (kMuonDCAx/y) is filled from the refitted parameters inside FillGlobalMuonRefit(Cov) + if constexpr (static_cast(TMFTFillMap & VarManager::ObjTypes::MFTCov)) { + auto const& mfttrackcov = mfCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]); + VarManager::FillGlobalMuonRefitCov(muontrack, mfttrack, collision, mfttrackcov); + } else { + VarManager::FillGlobalMuonRefit(muontrack, mfttrack, collision); + } } else { - VarManager::FillGlobalMuonRefit(muontrack, mfttrack, collision); + // Helix DCA of the global track; leaves kMuonPDca from the matched MCH above + VarManager::FillTrackCollision(muon, collision); } } else { VarManager::FillTrackCollision(muon, collision); @@ -1266,17 +1273,24 @@ struct TableMakerMC { VarManager::FillPropagateMuon(muon, collision); } // recalculate pDca / DCA and global muon kinematics + // kMuonPDca is always taken from MCH (standalone or the MCH matched to a global) int globalClusters = muon.nClusters(); - if (static_cast(muon.trackType()) < 2 && fConfigVariousOptions.fRefitGlobalMuon) { + if (static_cast(muon.trackType()) < 2) { auto muontrack = muon.template matchMCHTrack_as(); - auto mfttrack = muon.template matchMFTTrack_as(); - globalClusters += mfttrack.nClusters(); VarManager::FillTrackCollision(muontrack, collision); - if constexpr (static_cast(TMFTFillMap & VarManager::ObjTypes::MFTCov)) { - auto const& mfttrackcov = mfCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]); - VarManager::FillGlobalMuonRefitCov(muontrack, mfttrack, collision, mfttrackcov); + if (fConfigVariousOptions.fRefitGlobalMuon) { + auto mfttrack = muon.template matchMFTTrack_as(); + globalClusters += mfttrack.nClusters(); + // Helix DCA (kMuonDCAx/y) is filled from the refitted parameters inside FillGlobalMuonRefit(Cov) + if constexpr (static_cast(TMFTFillMap & VarManager::ObjTypes::MFTCov)) { + auto const& mfttrackcov = mfCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]); + VarManager::FillGlobalMuonRefitCov(muontrack, mfttrack, collision, mfttrackcov); + } else { + VarManager::FillGlobalMuonRefit(muontrack, mfttrack, collision); + } } else { - VarManager::FillGlobalMuonRefit(muontrack, mfttrack, collision); + // Helix DCA of the global track; leaves kMuonPDca from the matched MCH above + VarManager::FillTrackCollision(muon, collision); } } else { VarManager::FillTrackCollision(muon, collision); diff --git a/PWGDQ/TableProducer/tableMaker_withAssoc.cxx b/PWGDQ/TableProducer/tableMaker_withAssoc.cxx index 616fb49d9e1..25fda45923b 100644 --- a/PWGDQ/TableProducer/tableMaker_withAssoc.cxx +++ b/PWGDQ/TableProducer/tableMaker_withAssoc.cxx @@ -1688,24 +1688,30 @@ struct TableMaker { // So if a muon is associated to multiple collisions, depending on the selections, // it may be accepted for some associations and rejected for other if (static_cast(muon.trackType()) > 2 && fConfigVariousOptions.fPropMuon) { - VarManager::FillPropagateMuon(muon, collision); + VarManager::FillPropagateMuon(muon, collision); } // recalculate pDca / DCA and global muon kinematics - if (static_cast(muon.trackType()) < 2 && fConfigVariousOptions.fRefitGlobalMuon) { + // kMuonPDca is always taken from MCH (standalone or the MCH matched to a global) + if (static_cast(muon.trackType()) < 2) { auto muontrack = muon.template matchMCHTrack_as(); - if (muontrack.eta() < fConfigVariousOptions.fMuonMatchEtaMin || muontrack.eta() > fConfigVariousOptions.fMuonMatchEtaMax) { - continue; - } - auto mfttrack = muon.template matchMFTTrack_as(); VarManager::FillTrackCollision(muontrack, collision); - // NOTE: the MFT track originally associated to the MUON track is currently used in the global muon refit - // Should MUON - MFT time ambiguities be taken into account ? - // Helix DCA is filled from the refitted parameters inside FillGlobalMuonRefit(Cov) - if constexpr (static_cast(TMFTFillMap & VarManager::ObjTypes::MFTCov)) { - auto const& mfttrackcov = mfCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]); - VarManager::FillGlobalMuonRefitCov(muontrack, mfttrack, collision, mfttrackcov); + if (fConfigVariousOptions.fRefitGlobalMuon) { + if (muontrack.eta() < fConfigVariousOptions.fMuonMatchEtaMin || muontrack.eta() > fConfigVariousOptions.fMuonMatchEtaMax) { + continue; + } + auto mfttrack = muon.template matchMFTTrack_as(); + // NOTE: the MFT track originally associated to the MUON track is currently used in the global muon refit + // Should MUON - MFT time ambiguities be taken into account ? + // Helix DCA (kMuonDCAx/y) is filled from the refitted parameters inside FillGlobalMuonRefit(Cov) + if constexpr (static_cast(TMFTFillMap & VarManager::ObjTypes::MFTCov)) { + auto const& mfttrackcov = mfCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]); + VarManager::FillGlobalMuonRefitCov(muontrack, mfttrack, collision, mfttrackcov); + } else { + VarManager::FillGlobalMuonRefit(muontrack, mfttrack, collision); + } } else { - VarManager::FillGlobalMuonRefit(muontrack, mfttrack, collision); + // Helix DCA of the global track; leaves kMuonPDca from the matched MCH above + VarManager::FillTrackCollision(muon, collision); } } else { VarManager::FillTrackCollision(muon, collision); @@ -1784,20 +1790,27 @@ struct TableMaker { VarManager::FillTrack(muon); if (static_cast(muon.trackType()) > 2 && fConfigVariousOptions.fPropMuon) { - VarManager::FillPropagateMuon(muon, collision); + VarManager::FillPropagateMuon(muon, collision); } // recalculate pDca / DCA and global muon kinematics + // kMuonPDca is always taken from MCH (standalone or the MCH matched to a global) int globalClusters = muon.nClusters(); - if (static_cast(muon.trackType()) < 2 && fConfigVariousOptions.fRefitGlobalMuon) { + if (static_cast(muon.trackType()) < 2) { auto muontrack = muon.template matchMCHTrack_as(); - auto mfttrack = muon.template matchMFTTrack_as(); - globalClusters += mfttrack.nClusters(); VarManager::FillTrackCollision(muontrack, collision); - if constexpr (static_cast(TMFTFillMap & VarManager::ObjTypes::MFTCov)) { - auto const& mfttrackcov = mfCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]); - VarManager::FillGlobalMuonRefitCov(muontrack, mfttrack, collision, mfttrackcov); + if (fConfigVariousOptions.fRefitGlobalMuon) { + auto mfttrack = muon.template matchMFTTrack_as(); + globalClusters += mfttrack.nClusters(); + // Helix DCA (kMuonDCAx/y) is filled from the refitted parameters inside FillGlobalMuonRefit(Cov) + if constexpr (static_cast(TMFTFillMap & VarManager::ObjTypes::MFTCov)) { + auto const& mfttrackcov = mfCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]); + VarManager::FillGlobalMuonRefitCov(muontrack, mfttrack, collision, mfttrackcov); + } else { + VarManager::FillGlobalMuonRefit(muontrack, mfttrack, collision); + } } else { - VarManager::FillGlobalMuonRefit(muontrack, mfttrack, collision); + // Helix DCA of the global track; leaves kMuonPDca from the matched MCH above + VarManager::FillTrackCollision(muon, collision); } } else { VarManager::FillTrackCollision(muon, collision); From 3841990c99dd4f4711ac860d340edcc9080c8abb Mon Sep 17 00:00:00 2001 From: Maurice Coquet Date: Sun, 13 Sep 2026 10:25:11 +0200 Subject: [PATCH 6/9] removing old tests --- PWGDQ/Tasks/tests/CMakeLists.txt | 170 -- PWGDQ/Tasks/tests/FwdTrackReAlignTables.h | 90 - .../tests/fwdtrackToCollisionAssociator.cxx | 127 -- PWGDQ/Tasks/tests/global-muon-matcher.cxx | 1681 ----------------- 4 files changed, 2068 deletions(-) delete mode 100644 PWGDQ/Tasks/tests/CMakeLists.txt delete mode 100644 PWGDQ/Tasks/tests/FwdTrackReAlignTables.h delete mode 100644 PWGDQ/Tasks/tests/fwdtrackToCollisionAssociator.cxx delete mode 100644 PWGDQ/Tasks/tests/global-muon-matcher.cxx diff --git a/PWGDQ/Tasks/tests/CMakeLists.txt b/PWGDQ/Tasks/tests/CMakeLists.txt deleted file mode 100644 index 29ac8ae0ba3..00000000000 --- a/PWGDQ/Tasks/tests/CMakeLists.txt +++ /dev/null @@ -1,170 +0,0 @@ -# Copyright 2019-2020 CERN and copyright holders of ALICE O2. -# See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -# All rights not expressly granted are reserved. -# -# This software is distributed under the terms of the GNU General Public -# License v3 (GPL Version 3), copied verbatim in the file "COPYING". -# -# In applying this license CERN does not waive the privileges and immunities -# granted to it by virtue of its status as an Intergovernmental Organization -# or submit itself to any jurisdiction. - -o2physics_add_dpl_workflow(table-reader - SOURCES tableReader.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore O2Physics::MLCore - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(table-reader-with-assoc - SOURCES tableReader_withAssoc.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2Physics::AnalysisCCDB O2Physics::PWGDQCore O2Physics::MLCore - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(table-reader-with-assoc-direct - SOURCES tableReader_withAssoc_direct.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2Physics::AnalysisCCDB O2Physics::PWGDQCore O2Physics::MLCore O2::ReconstructionDataFormats O2::DetectorsCommonDataFormats O2::DetectorsVertexing O2Physics::EventFilteringUtils - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(efficiency - SOURCES dqEfficiency.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(efficiency-with-assoc - SOURCES dqEfficiency_withAssoc.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(efficiency-with-assoc-direct - SOURCES dqEfficiency_withAssoc_direct.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore O2::ReconstructionDataFormats O2::DetectorsCommonDataFormats O2::DetectorsVertexing - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(energy-correlator-direct - SOURCES dqEnergyCorrelator_direct.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(filter-pp - SOURCES filterPP.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(filter-pp-with-association - SOURCES filterPPwithAssociation.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(filter-pb-pb - SOURCES filterPbPb.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore O2Physics::SGCutParHolder - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(v0-selector - SOURCES v0selector.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DCAFitter O2Physics::AnalysisCore O2Physics::PWGDQCore - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(dalitz-selection - SOURCES DalitzSelection.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(flow - SOURCES dqFlow.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore O2Physics::GFWCore - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(task-muon-mch-trk-eff - SOURCES taskMuonMchTrkEfficiency.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(task-j-psi-hf - SOURCES taskJpsiHf.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(task-muon-dca - SOURCES muonDCA.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(correlation - SOURCES dqCorrelation.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore O2Physics::GFWCore - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(task-mch-align-record - SOURCES mchAlignRecord.cxx - PUBLIC_LINK_LIBRARIES - O2::Framework - O2Physics::AnalysisCore - O2Physics::PWGDQCore - O2::CommonUtils - O2::MCHClustering - O2::DPLUtils - O2::CCDB - O2::DataFormatsParameters - O2::MCHBase - O2::MCHTracking - O2::DataFormatsMCH - O2::DetectorsBase - O2::MCHGeometryTransformer - O2::MathUtils - O2::MCHAlign - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(task-muon-mid-eff - SOURCES MIDefficiency.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::MIDBase - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(task-fwd-track-pid - SOURCES taskFwdTrackPid.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(quarkonia-to-hyperons - SOURCES quarkoniaToHyperons.cxx - PUBLIC_LINK_LIBRARIES O2::DetectorsBase O2::Framework O2::DCAFitter KFParticle::KFParticle O2Physics::AnalysisCore O2Physics::MLCore O2Physics::EventFilteringUtils - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(model-converter-mult-pv - SOURCES ModelConverterMultPv.cxx - PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::PWGDQCore - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(model-converter-event-extended - SOURCES ModelConverterEventExtended.cxx - PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::PWGDQCore - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(model-converter-mc-reduced-event - SOURCES ModelConverterReducedMCEvents.cxx - PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::PWGDQCore - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(tag-and-probe - SOURCES TagAndProbe.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2Physics::AnalysisCCDB O2Physics::PWGDQCore - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(qa-matching - SOURCES qaMatching.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(mft-mch-matcher - SOURCES mftMchMatcher.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(muon-global-alignment - SOURCES muonGlobalAlignment.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore O2::MCHGeometryTransformer - COMPONENT_NAME Analysis) - -o2physics_add_dpl_workflow(global-muon-matcher - SOURCES global-muon-matcher.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGDQCore O2::MCHGeometryTransformer - COMPONENT_NAME Analysis) diff --git a/PWGDQ/Tasks/tests/FwdTrackReAlignTables.h b/PWGDQ/Tasks/tests/FwdTrackReAlignTables.h deleted file mode 100644 index 879d3fe7486..00000000000 --- a/PWGDQ/Tasks/tests/FwdTrackReAlignTables.h +++ /dev/null @@ -1,90 +0,0 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. -// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -// All rights not expressly granted are reserved. -// -// This software is distributed under the terms of the GNU General Public -// License v3 (GPL Version 3), copied verbatim in the file "COPYING". -// -// In applying this license CERN does not waive the privileges and immunities -// granted to it by virtue of its status as an Intergovernmental Organization -// or submit itself to any jurisdiction. - -/// \file FwdTrackReAlignTables.h -/// \brief Table definitions for re-aligned forward tracks -/// \author Chi Zhang , CEA-Saclay - -#ifndef COMMON_DATAMODEL_FWDTRACKREALIGNTABLES_H_ -#define COMMON_DATAMODEL_FWDTRACKREALIGNTABLES_H_ - -#include - -namespace o2::aod -{ -namespace fwdtrackrealign -{ -DECLARE_SOA_COLUMN(IsRemovable, isRemovable, int); //! flag to check the refit status -} - -DECLARE_SOA_TABLE_FULL(StoredFwdTracksReAlign, "FwdTracksReAlign", "AOD", "FWDTRACKREALIGN", - o2::soa::Index<>, fwdtrack::CollisionId, fwdtrack::TrackType, - fwdtrack::X, fwdtrack::Y, fwdtrack::Z, fwdtrack::Phi, fwdtrack::Tgl, - fwdtrack::Signed1Pt, fwdtrack::NClusters, fwdtrack::PDca, fwdtrack::RAtAbsorberEnd, - fwdtrackrealign::IsRemovable, - fwdtrack::Px, - fwdtrack::Py, - fwdtrack::Pz, - fwdtrack::Sign, - fwdtrack::Chi2, fwdtrack::Chi2MatchMCHMID, fwdtrack::Chi2MatchMCHMFT, - fwdtrack::MatchScoreMCHMFT, fwdtrack::MFTTrackId, fwdtrack::MCHTrackId, - fwdtrack::MCHBitMap, fwdtrack::MIDBitMap, fwdtrack::MIDBoards, - fwdtrack::TrackTime, fwdtrack::TrackTimeRes); - -// extended table with expression columns that can be used as arguments of dynamic columns -DECLARE_SOA_EXTENDED_TABLE_USER(FwdTracksReAlign, StoredFwdTracksReAlign, "FWDTRKREALIGNEXT", //! - fwdtrack::Pt, - fwdtrack::Eta, - fwdtrack::P); - -DECLARE_SOA_TABLE_FULL(StoredFwdTrksCovReAlign, "FwdCovsReAlign", "AOD", "FWDCOVREALIGN", - fwdtrack::SigmaX, fwdtrack::SigmaY, fwdtrack::SigmaPhi, fwdtrack::SigmaTgl, fwdtrack::Sigma1Pt, - fwdtrack::RhoXY, fwdtrack::RhoPhiY, fwdtrack::RhoPhiX, fwdtrack::RhoTglX, fwdtrack::RhoTglY, - fwdtrack::RhoTglPhi, fwdtrack::Rho1PtX, fwdtrack::Rho1PtY, fwdtrack::Rho1PtPhi, fwdtrack::Rho1PtTgl); - -// extended table with expression columns that can be used as arguments of dynamic columns -DECLARE_SOA_EXTENDED_TABLE_USER(FwdTrksCovReAlign, StoredFwdTrksCovReAlign, "FWDCOVREALIGNEXT", //! - fwdtrack::CXX, - fwdtrack::CXY, - fwdtrack::CYY, - fwdtrack::CPhiX, - fwdtrack::CPhiY, - fwdtrack::CPhiPhi, - fwdtrack::CTglX, - fwdtrack::CTglY, - fwdtrack::CTglPhi, - fwdtrack::CTglTgl, - fwdtrack::C1PtX, - fwdtrack::C1PtY, - fwdtrack::C1PtPhi, - fwdtrack::C1PtTgl, - fwdtrack::C1Pt21Pt2); - -using FwdTrackRealign = FwdTracksReAlign::iterator; -using FwdTrkCovRealign = FwdTrksCovReAlign::iterator; -using FullFwdTracksRealign = soa::Join; -using FullFwdTrackRealign = FullFwdTracksRealign::iterator; - -// ambiguity table for realigned muons -namespace fwdtrackrealignambiguous -{ -DECLARE_SOA_INDEX_COLUMN_FULL(FwdTrackRealign, fwdTrackRealign, int, FwdTracksReAlign, ""); //! FwdTracksReAlign index -DECLARE_SOA_SLICE_INDEX_COLUMN(BC, bc); - -} // namespace fwdtrackrealignambiguous - -DECLARE_SOA_TABLE(AmbiguousFwdTrksReAlign, "AOD", "AMBIFWDREALIGN", //! Table for FwdTracksReAlign which are not associated with a collision - o2::soa::Index<>, fwdtrackrealignambiguous::FwdTrackRealignId, fwdtrackrealignambiguous::BCIdSlice); - -using AmbiguousFwdTrkRealign = AmbiguousFwdTrksReAlign::iterator; -} // namespace o2::aod - -#endif // COMMON_DATAMODEL_FWDTRACKREALIGNTABLES_H_ diff --git a/PWGDQ/Tasks/tests/fwdtrackToCollisionAssociator.cxx b/PWGDQ/Tasks/tests/fwdtrackToCollisionAssociator.cxx deleted file mode 100644 index 1de267d0e60..00000000000 --- a/PWGDQ/Tasks/tests/fwdtrackToCollisionAssociator.cxx +++ /dev/null @@ -1,127 +0,0 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. -// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -// All rights not expressly granted are reserved. -// -// This software is distributed under the terms of the GNU General Public -// License v3 (GPL Version 3), copied verbatim in the file "COPYING". -// -// In applying this license CERN does not waive the privileges and immunities -// granted to it by virtue of its status as an Intergovernmental Organization -// or submit itself to any jurisdiction. - -/// \file fwdtrackToCollisionAssociator.cxx -/// \brief Associates fwd and MFT tracks to collisions considering ambiguities -/// \author Sarah Herrmann , IP2I Lyon -/// \author Maurice Coquet , CEA-Saclay/Irfu - -#include "Common/Core/CollisionAssociation.h" -#include "Common/DataModel/CollisionAssociationTables.h" -#include "Common/DataModel/FwdTrackReAlignTables.h" - -#include -#include -#include -#include -#include -#include - -using namespace o2; -using namespace o2::framework; -using namespace o2::framework::expressions; -using namespace o2::aod; - -struct FwdTrackToCollisionAssociation { - Produces fwdassociation; - Produces fwdreverseIndices; - Produces mftassociation; - Produces mftreverseIndices; - - Configurable nSigmaForTimeCompat{"nSigmaForTimeCompat", 4.f, "number of sigmas for time compatibility"}; - Configurable timeMargin{"timeMargin", 0.f, "time margin in ns added to uncertainty because of uncalibrated TPC"}; - Configurable includeUnassigned{"includeUnassigned", false, "consider also tracks which are not assigned to any collision"}; - Configurable fillTableOfCollIdsPerTrack{"fillTableOfCollIdsPerTrack", false, "fill additional table with vector of collision ids per track"}; - Configurable bcWindowForOneSigma{"bcWindowForOneSigma", 115, "BC window to be multiplied by the number of sigmas to define maximum window to be considered"}; - - CollisionAssociation collisionAssociator; - - Preslice muonsPerCollisions = aod::fwdtrack::collisionId; - Preslice realignmuonsPerCollisions = aod::fwdtrack::collisionId; - Preslice mftsPerCollisions = aod::fwdtrack::collisionId; - - void init(InitContext const&) - { - if (doprocessFwdAssocWithTime && doprocessFwdStandardAssoc) { - LOGP(fatal, "Exactly one process function between standard and time-based association should be enabled!"); - } - if (doprocessMFTAssocWithTime && doprocessMFTStandardAssoc) { - LOGP(fatal, "Exactly one process function between standard and time-based association should be enabled!"); - } - - if (!(doprocessMFTAssocWithTime || doprocessMFTStandardAssoc || doprocessFwdAssocWithTime || doprocessFwdStandardAssoc || doprocessFwdRealignAssocWithTime || doprocessFwdRealignStandardAssoc)) { - LOGP(fatal, "At least one process function should be enabled!"); - } - - // set options in track-to-collision association - collisionAssociator.setNumSigmaForTimeCompat(nSigmaForTimeCompat); - collisionAssociator.setTimeMargin(timeMargin); - collisionAssociator.setTrackSelectionOptionForStdAssoc(track_association::TrackSelection::None); - collisionAssociator.setUsePvAssociation(track_association::PVContrReassocOpt::Disabled); - collisionAssociator.setIncludeUnassigned(includeUnassigned); - collisionAssociator.setFillTableOfCollIdsPerTrack(fillTableOfCollIdsPerTrack); - collisionAssociator.setBcWindow(bcWindowForOneSigma); - } - - void processFwdAssocWithTime(Collisions const& collisions, - FwdTracks const& muons, - AmbiguousFwdTracks const& ambiTracksFwd, - BCs const& bcs) - { - collisionAssociator.runAssocWithTime(collisions, muons, muons, ambiTracksFwd, bcs, fwdassociation, fwdreverseIndices); - } - PROCESS_SWITCH(FwdTrackToCollisionAssociation, processFwdAssocWithTime, "Use fwdtrack-to-collision association based on time", true); - - void processFwdStandardAssoc(Collisions const& collisions, - FwdTracks const& muons) - { - collisionAssociator.runStandardAssoc(collisions, muons, muonsPerCollisions, fwdassociation, fwdreverseIndices); - } - PROCESS_SWITCH(FwdTrackToCollisionAssociation, processFwdStandardAssoc, "Use standard fwdtrack-to-collision association", false); - - void processFwdRealignAssocWithTime(Collisions const& collisions, - FwdTracksReAlign const& muons, - AmbiguousFwdTrksReAlign const& ambiTracksFwd, - BCs const& bcs) - { - collisionAssociator.runAssocWithTime(collisions, muons, muons, ambiTracksFwd, bcs, fwdassociation, fwdreverseIndices); - } - PROCESS_SWITCH(FwdTrackToCollisionAssociation, processFwdRealignAssocWithTime, "Use fwdrealigntrack-to-collision association based on time", false); - - void processFwdRealignStandardAssoc(Collisions const& collisions, - FwdTracksReAlign const& muons) - { - collisionAssociator.runStandardAssoc(collisions, muons, realignmuonsPerCollisions, fwdassociation, fwdreverseIndices); - } - PROCESS_SWITCH(FwdTrackToCollisionAssociation, processFwdRealignStandardAssoc, "Use standard fwdrealigntrack-to-collision association", false); - - void processMFTAssocWithTime(Collisions const& collisions, - MFTTracks const& tracks, - AmbiguousMFTTracks const& ambiguousTracks, - BCs const& bcs) - { - collisionAssociator.runAssocWithTime(collisions, tracks, tracks, ambiguousTracks, bcs, mftassociation, mftreverseIndices); - } - PROCESS_SWITCH(FwdTrackToCollisionAssociation, processMFTAssocWithTime, "Use MFTtrack-to-collision association based on time", true); - - void processMFTStandardAssoc(Collisions const& collisions, - MFTTracks const& tracks) - { - collisionAssociator.runStandardAssoc(collisions, tracks, mftsPerCollisions, mftassociation, mftreverseIndices); - } - PROCESS_SWITCH(FwdTrackToCollisionAssociation, processMFTStandardAssoc, "Use standard mfttrack-to-collision association", false); -}; - -//________________________________________________________________________________________________________________________ -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) -{ - return WorkflowSpec{adaptAnalysisTask(cfgc)}; -} diff --git a/PWGDQ/Tasks/tests/global-muon-matcher.cxx b/PWGDQ/Tasks/tests/global-muon-matcher.cxx deleted file mode 100644 index 28b8ef02fcd..00000000000 --- a/PWGDQ/Tasks/tests/global-muon-matcher.cxx +++ /dev/null @@ -1,1681 +0,0 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. -// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. -// All rights not expressly granted are reserved. -// -// This software is distributed under the terms of the GNU General Public -// License v3 (GPL Version 3), copied verbatim in the file "COPYING". -// -// In applying this license CERN does not waive the privileges and immunities -// granted to it by virtue of its status as an Intergovernmental Organization -// or submit itself to any jurisdiction. -// -/// \file global-muon-matcher.cxx -/// \brief Task for analysis MFT-MCH muon matching -/// \author Andrea Ferrero -/// -#include "PWGDQ/Core/MuonMatchingMlResponse.h" -#include "PWGDQ/Core/VarManager.h" - -#include "Common/Core/fwdtrackUtilities.h" -#include "Common/DataModel/EventSelection.h" -#include "Common/DataModel/FwdTrackReAlignTables.h" -#include "Tools/ML/MlResponse.h" - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -#include -#include -#include -#include - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include - -using namespace o2; -using namespace o2::framework; -using namespace o2::aod; - -namespace o2::aod::globalmuonmatching -{ -DECLARE_SOA_COLUMN(MchTrackId, mchTrackId, int64_t); -DECLARE_SOA_COLUMN(MftTrackId, mftTrackId, int64_t); -DECLARE_SOA_COLUMN(MatchChi2, matchChi2, float); -DECLARE_SOA_COLUMN(MatchScore, matchScore, float); -DECLARE_SOA_COLUMN(MatchRanking, matchRanking, int32_t); -DECLARE_SOA_COLUMN(IsTagged, isTagged, bool); -} // namespace o2::aod::globalmuonmatching - -namespace o2::aod -{ -DECLARE_SOA_TABLE(GlobalMuonMatchCandidates, "AOD", "GMCAND", - o2::soa::Index<>, - globalmuonmatching::MchTrackId, - globalmuonmatching::MftTrackId, - globalmuonmatching::MatchChi2, globalmuonmatching::MatchScore, globalmuonmatching::MatchRanking, - globalmuonmatching::IsTagged); - -namespace globalmuonmatching -{ -DECLARE_SOA_ARRAY_INDEX_COLUMN(GlobalMuonMatchCandidate, globalMuonMatchCandidate); //! Array of GlobalMuonMatchCandidates indices -} // namespace globalmuonmatching - -DECLARE_SOA_TABLE(FwdTrkMatchCands, "AOD", "FWDTRKMATCHCAND", //! Vectors of match-candidate indices stored per fwdtrack - globalmuonmatching::GlobalMuonMatchCandidateIds, o2::soa::Marker<3>); - -} // namespace o2::aod - -using MyEvents = soa::Join; -using MyMuons = soa::Join; -using MyMFTs = aod::MFTTracks; -using MyMFTCovariances = aod::MFTTracksCov; - -using SMatrix55Sym = o2::track::SMatrix55Sym; -using SMatrix55Std = o2::track::SMatrix55Std; -using SMatrix5 = o2::track::SMatrix5; - -constexpr std::array NDetElemCh = {4, 4, 4, 4, 18, 18, 26, 26, 26, 26}; -constexpr std::array SNDetElemCh = {0, 4, 8, 12, 16, 34, 52, 78, 104, 130, 156}; - -struct GlobalMuonMatching { - - static constexpr int GlobalTrackTypeMax = 2; - static constexpr int MchMidTrackType = 3; - static constexpr int NMchChambers = 10; - static constexpr int MchDetElemNumberingBase = 100; - static constexpr int NMchDetElems = 156; - static constexpr int MinRemovableTrackClusters = 10; - static constexpr int ThetaAbsBoundaryDeg = 3; - static constexpr double SlopeResolutionZ = 535.; - static constexpr float MatchingPlaneDefaultZ = -77.5; - - struct MatchingCandidate { - int64_t muonTrackId{-1}; - int64_t mftTrackId{-1}; - double matchScore{-1}; - double matchChi2{-1}; - int matchRanking{-1}; - }; - - //// Variables for selecting tagged muons - struct : ConfigurableGroup { - Configurable cfgMuonTaggingNCrossedMftPlanesLow{"cfgMuonTaggingNCrossedMftPlanesLow", 5, ""}; - Configurable cfgMuonTaggingTrackChi2MchUp{"cfgMuonTaggingTrackChi2MchUp", 5.f, ""}; - Configurable cfgMuonTaggingPMchLow{"cfgMuonTaggingPMchLow", 0.0f, ""}; - Configurable cfgMuonTaggingPtMchLow{"cfgMuonTaggingPtMchLow", 0.7f, ""}; - Configurable cfgMuonTaggingEtaMchLow{"cfgMuonTaggingEtaMchLow", -3.6f, ""}; - Configurable cfgMuonTaggingEtaMchUp{"cfgMuonTaggingEtaMchUp", -2.5f, ""}; - Configurable cfgMuonTaggingRabsLow{"cfgMuonTaggingRabsLow", 17.6f, ""}; - Configurable cfgMuonTaggingRabsUp{"cfgMuonTaggingRabsUp", 89.5f, ""}; - Configurable cfgMuonTaggingPdcaUp{"cfgMuonTaggingPdcaUp", 4.f, ""}; - Configurable cfgMuonTaggingRadiusAtMftFrontLow{"cfgMuonTaggingRadiusAtMftFrontLow", 3.f, ""}; - Configurable cfgMuonTaggingRadiusAtMftFrontUp{"cfgMuonTaggingRadiusAtMftFrontUp", 9.f, ""}; - Configurable cfgMuonTaggingRadiusAtMftBackLow{"cfgMuonTaggingRadiusAtMftBackLow", 5.f, ""}; - Configurable cfgMuonTaggingRadiusAtMftBackUp{"cfgMuonTaggingRadiusAtMftBackUp", 12.f, ""}; - } configMuonTagging; - - //// Variables for MCH realignment - struct : ConfigurableGroup { - Configurable cfgEnableMCHRealign{"cfgEnableMCHRealign", true, "Enable re-alignment of MCH clusters and tracks"}; - Configurable cfgGeoRefPath{"cfgGeoRefPath", "GLO/Config/GeometryAligned", "Path of the reference geometry file"}; - Configurable cfgGeoNewPath{"cfgGeoNewPath", "GLO/Config/GeometryAligned", "Path of the new geometry file"}; - Configurable cfgCcdbNoLaterThanRef{"cfgCcdbNoLaterThanRef", std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(), "latest acceptable timestamp of creation for the object of reference basis"}; - Configurable cfgCcdbNoLaterThanNew{"cfgCcdbNoLaterThanNew", std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(), "latest acceptable timestamp of creation for the object of new basis"}; - Configurable cfgChamberResolutionX{"cfgChamberResolutionX", 0.04, "Chamber resolution along X configuration for refit"}; // 0.4cm pp, 0.2cm PbPb - Configurable cfgChamberResolutionY{"cfgChamberResolutionY", 0.04, "Chamber resolution along Y configuration for refit"}; // 0.4cm pp, 0.2cm PbPb - Configurable cfgSigmaCutImprove{"cfgSigmaCutImprove", 6., "Sigma cut for track improvement"}; // 6 for pp, 4 for PbPb - } configMchRealign; - - //// Variables for MFT alignment corrections - struct : ConfigurableGroup { - Configurable cfgEnableMftAlignmentCorrections{"cfgEnableMftAlignmentCorrections", true, "Enable alignment corrections for the MFT tracks"}; - // slope corrections - // Configurable cfgMFTAlignmentCorrXSlopeTop{"cfgMFTAlignmentCorrXSlopeTop", (-0.0006696 - 0.0005621) / 2.f, "MFT X slope correction - top half"}; - // Configurable cfgMFTAlignmentCorrXSlopeBottom{"cfgMFTAlignmentCorrXSlopeBottom", (0.00105 + 0.001007) / 2.f, "MFT X slope correction - bottom half"}; - // Configurable cfgMFTAlignmentCorrYSlopeTop{"cfgMFTAlignmentCorrYSlopeTop", (-0.002299 - 0.002442) / 2.f, "MFT Y slope correction - top half"}; - // Configurable cfgMFTAlignmentCorrYSlopeBottom{"cfgMFTAlignmentCorrYSlopeBottom", (-0.0005339 - 0.0006921) / 2.f, "MFT Y slope correction - bottom half"}; - Configurable cfgMFTAlignmentCorrXSlopeTop{"cfgMFTAlignmentCorrXSlopeTop", 0.f, "MFT X slope correction - top half"}; - Configurable cfgMFTAlignmentCorrXSlopeBottom{"cfgMFTAlignmentCorrXSlopeBottom", 0.f, "MFT X slope correction - bottom half"}; - Configurable cfgMFTAlignmentCorrYSlopeTop{"cfgMFTAlignmentCorrYSlopeTop", 0.f, "MFT Y slope correction - top half"}; - Configurable cfgMFTAlignmentCorrYSlopeBottom{"cfgMFTAlignmentCorrYSlopeBottom", 0.f, "MFT Y slope correction - bottom half"}; - // offset corrections - Configurable cfgMFTAlignmentCorrXOffsetTop{"cfgMFTAlignmentCorrXOffsetTop", 0.f, "MFT X offset correction - top half"}; - Configurable cfgMFTAlignmentCorrXOffsetBottom{"cfgMFTAlignmentCorrXOffsetBottom", 0.f, "MFT X offset correction - bottom half"}; - Configurable cfgMFTAlignmentCorrYOffsetTop{"cfgMFTAlignmentCorrYOffsetTop", 0.f, "MFT Y offset correction - top half"}; - Configurable cfgMFTAlignmentCorrYOffsetBottom{"cfgMFTAlignmentCorrYOffsetBottom", 0.f, "MFT Y offset correction - bottom half"}; - } configMftAlignmentCorrections; - - // Variables for CCDB objects access and retrieval - struct : ConfigurableGroup { - Configurable cfgCcdbUrl{"cfgCcdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; - Configurable cfgCcdbNoLaterThan{"cfgCcdbNoLaterThan", std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(), "latest acceptable timestamp of creation for the object"}; - Configurable cfgGrpPath{"cfgGrpPath", "GLO/GRP/GRP", "Path of the grp file"}; - Configurable cfgGeoPath{"cfgGeoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"}; - Configurable cfgGrpMagPath{"cfgGrpMagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; - } configCcdb; - - // Matching strategy for the *custom* matches (production baseline is always computed). - // 0 = chi2 (runChi2Matching), 1 = ML (runMlMatching) - struct : ConfigurableGroup { - Configurable cfgCustomMatchingStrategy{"cfgCustomMatchingStrategy", 0, "0=chi2, 1=ML for custom matches"}; - Configurable cfgIncludeGlobalMuonsInFwdTracks{"cfgIncludeGlobalMuonsInFwdTracks", false, "Include MFT-MCH-MID global muons in GMMCANDTRK table"}; - Configurable cfgMaxCandidatesPerMchTrack{"cfgMaxCandidatesPerMchTrack", -1, "Maximum number of match candidates stored per MCH track (-1: no limit)"}; - Configurable cfgMatchAllTracks{"cfgMatchAllTracks", false, "If true the matching is performed considering all the MFT tracks for which the covariances are available; if false the matching is performed considering only the global forward tracks stored at production"}; - } configMatching; - - double mBzAtMftCenter{0}; - - using MatchingFunc = std::function(const o2::track::TrackParCovFwd& mchtrack, const o2::track::TrackParCovFwd& mfttrack)>; - std::map mMatchingFunctionMap; ///< MFT-MCH Matching function - - // Chi2 matching interface (single configurable method) - struct : ConfigurableGroup { - Configurable cfgChi2FunctionLabel{"cfgChi2FunctionLabel", std::string{"ProdAll"}, "Text label identifying the chi2 matching method"}; - Configurable cfgChi2FunctionName{"cfgChi2FunctionName", std::string{"prod"}, "Name of the chi2 matching function"}; - Configurable cfgChi2FunctionMatchingPlaneZ{"cfgChi2FunctionMatchingPlaneZ", static_cast(o2::mft::constants::mft::LayerZCoordinate()[9]), "Z position of the matching plane"}; - } configChi2MatchingOptions; - - // ML interface (single configurable model) - struct : ConfigurableGroup { - Configurable cfgMlModelLabel{"cfgMlModelLabel", std::string{""}, "Text label identifying this ML model"}; - Configurable cfgMlModelPathCcdb{"cfgMlModelPathCcdb", "Users/m/mcoquet/MLTest", "Path of model on CCDB"}; - Configurable cfgMlModelName{"cfgMlModelName", "model.onnx", "ONNX file name (if not from CCDB full path)"}; - Configurable> cfgMlInputFeatures{"cfgMlInputFeatures", std::vector{"chi2MCHMFT"}, "Names of ML model input features"}; - Configurable cfgMlModelMatchingPlaneZ{"cfgMlModelMatchingPlaneZ", static_cast(o2::mft::constants::mft::LayerZCoordinate()[9]), "Z position of the matching plane"}; - } configMlOptions; - - std::vector binsPtMl; - std::array cutValues{}; - std::vector cutDirMl; - bool hasActiveChi2Matching{false}; - std::string activeChi2FunctionName; - double activeChi2MatchingPlaneZ{0.}; - - bool hasActiveMlMatching{false}; - o2::analysis::MlResponseMFTMuonMatch activeMlResponse; - double activeMlMatchingPlaneZ{0.}; - - int mRunNumber{0}; // needed to detect if the run changed and trigger update of magnetic field - - Service ccdbManager{}; - o2::ccdb::CcdbApi fCCDBApi; - - // 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 score - // the map key is the MCH(-MID) track global index - using MatchingCandidates = std::map>; - std::map> mMatchingCandidates; - - 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 *this; } - - private: - int nClusters{-1}; - bool removable{false}; - }; - - std::unordered_map mMchTrackPars; - std::unordered_map mMftTrackPars; - - std::unordered_map mftTrackCovs; - - Produces globalMuonMatchCandidates; - Produces fwdTrkMatchCands; - Produces gmCandidateFwdTracks; - Produces gmCandidateFwdTracksCov; - Produces gmAmbiguousFwdTracksReAlign; - - int32_t mGmmCandFwdTrackRowIndex{0}; - std::unordered_map> mAmbBcSliceByFwdTrackId; - bool mHasLastMchAmbiguousBcSlice{false}; - std::array mLastMchAmbiguousBcSlice{}; - - int32_t mMatchCandidateCounter{0}; - std::unordered_map> mMchTrackToCandidateIndices; - std::unordered_map> mMchTrackMatchingCandidates; - std::unordered_map mFwdTrackToGmmCandTrkIndex; - - mch::TrackFitter trackFitter; // Track fitter from MCH tracking library - mch::geo::TransformationCreator transformation; - std::map transformRef; // reference geometry w.r.t track data - std::map transformNew; // new geometry - double mImproveCutChi2{0.}; // Chi2 cut for track improvement. - TGeoManager* geoNew = nullptr; - TGeoManager* geoRef = nullptr; - globaltracking::MatchGlobalFwd mMatching; - - Preslice perMuon = aod::fwdtrkcl::fwdtrackId; - - template - o2::mch::TrackParam fwdToMch(const T& fwdtrack) - { - // Convert Forward Track parameters and covariances matrix to the - // MCH track format. - - // Parameter conversion - const double x2 = fwdtrack.getPhi(); - const double x3 = fwdtrack.getTanl(); - const double x4 = fwdtrack.getInvQPt(); - - const auto sinX2 = std::sin(x2); - const auto cosX2 = std::cos(x2); - - const double alpha1 = cosX2 / x3; - const double alpha3 = sinX2 / x3; - const double alpha4 = x4 / std::sqrt(x3 * x3 + sinX2 * sinX2); - - const auto kNorm = std::sqrt(x3 * x3 + sinX2 * sinX2); - const auto kNorm3 = kNorm * kNorm * kNorm; - - // Covariances matrix conversion - SMatrix55Std jacobian; - SMatrix55Sym covariances; - - covariances(0, 0) = fwdtrack.getCovariances()(0, 0); - covariances(0, 1) = fwdtrack.getCovariances()(0, 1); - covariances(0, 2) = fwdtrack.getCovariances()(0, 2); - covariances(0, 3) = fwdtrack.getCovariances()(0, 3); - covariances(0, 4) = fwdtrack.getCovariances()(0, 4); - - covariances(1, 1) = fwdtrack.getCovariances()(1, 1); - covariances(1, 2) = fwdtrack.getCovariances()(1, 2); - covariances(1, 3) = fwdtrack.getCovariances()(1, 3); - covariances(1, 4) = fwdtrack.getCovariances()(1, 4); - - covariances(2, 2) = fwdtrack.getCovariances()(2, 2); - covariances(2, 3) = fwdtrack.getCovariances()(2, 3); - covariances(2, 4) = fwdtrack.getCovariances()(2, 4); - - covariances(3, 3) = fwdtrack.getCovariances()(3, 3); - covariances(3, 4) = fwdtrack.getCovariances()(3, 4); - - covariances(4, 4) = fwdtrack.getCovariances()(4, 4); - - jacobian(0, 0) = 1; - - jacobian(1, 2) = -sinX2 / x3; - jacobian(1, 3) = -cosX2 / (x3 * x3); - - jacobian(2, 1) = 1; - - jacobian(3, 2) = cosX2 / x3; - jacobian(3, 3) = -sinX2 / (x3 * x3); - - jacobian(4, 2) = -x4 * sinX2 * cosX2 / kNorm3; - jacobian(4, 3) = -x3 * x4 / kNorm3; - jacobian(4, 4) = 1 / kNorm; - // jacobian*covariances*jacobian^T - covariances = ROOT::Math::Similarity(jacobian, covariances); - - std::array cov = {covariances(0, 0), covariances(1, 0), covariances(1, 1), covariances(2, 0), covariances(2, 1), covariances(2, 2), covariances(3, 0), covariances(3, 1), covariances(3, 2), covariances(3, 3), covariances(4, 0), covariances(4, 1), covariances(4, 2), covariances(4, 3), covariances(4, 4)}; - std::array param = {fwdtrack.getX(), alpha1, fwdtrack.getY(), alpha3, alpha4}; - - o2::mch::TrackParam convertedTrack(fwdtrack.getZ(), param.data(), cov.data()); - return {convertedTrack}; - } - - o2::track::TrackParCovFwd mchToFwd(const o2::mch::TrackParam& mchParam) - { - // Convert a MCH Track parameters and covariances matrix to the - // Forward track format. Must be called after propagation though the absorber - - o2::track::TrackParCovFwd convertedTrack; - - // Parameter conversion - const double alpha1 = mchParam.getNonBendingSlope(); - const double alpha3 = mchParam.getBendingSlope(); - const double alpha4 = mchParam.getInverseBendingMomentum(); - - const double x2 = std::atan2(-alpha3, -alpha1); - const double x3 = -1. / std::sqrt(alpha3 * alpha3 + alpha1 * alpha1); - const double x4 = alpha4 * -x3 * std::sqrt(1 + alpha3 * alpha3); - - const auto kNorm = alpha1 * alpha1 + alpha3 * alpha3; - const auto kNorm32 = kNorm * std::sqrt(kNorm); - const auto slopeLen = std::sqrt(alpha3 * alpha3 + 1); - - // Covariances matrix conversion - SMatrix55Std jacobian; - SMatrix55Sym covariances; - - covariances(0, 0) = mchParam.getCovariances()(0, 0); - covariances(0, 1) = mchParam.getCovariances()(0, 1); - covariances(0, 2) = mchParam.getCovariances()(0, 2); - covariances(0, 3) = mchParam.getCovariances()(0, 3); - covariances(0, 4) = mchParam.getCovariances()(0, 4); - - covariances(1, 1) = mchParam.getCovariances()(1, 1); - covariances(1, 2) = mchParam.getCovariances()(1, 2); - covariances(1, 3) = mchParam.getCovariances()(1, 3); - covariances(1, 4) = mchParam.getCovariances()(1, 4); - - covariances(2, 2) = mchParam.getCovariances()(2, 2); - covariances(2, 3) = mchParam.getCovariances()(2, 3); - covariances(2, 4) = mchParam.getCovariances()(2, 4); - - covariances(3, 3) = mchParam.getCovariances()(3, 3); - covariances(3, 4) = mchParam.getCovariances()(3, 4); - - covariances(4, 4) = mchParam.getCovariances()(4, 4); - - jacobian(0, 0) = 1; - - jacobian(1, 2) = 1; - - jacobian(2, 1) = -alpha3 / kNorm; - jacobian(2, 3) = alpha1 / kNorm; - - jacobian(3, 1) = alpha1 / kNorm32; - jacobian(3, 3) = alpha3 / kNorm32; - - jacobian(4, 1) = -alpha1 * alpha4 * slopeLen / kNorm32; - jacobian(4, 3) = alpha3 * alpha4 * (1 / (std::sqrt(kNorm) * slopeLen) - slopeLen / kNorm32); - jacobian(4, 4) = slopeLen / std::sqrt(kNorm); - - // jacobian*covariances*jacobian^T - covariances = ROOT::Math::Similarity(jacobian, covariances); - - // Set output - convertedTrack.setX(mchParam.getNonBendingCoor()); - convertedTrack.setY(mchParam.getBendingCoor()); - convertedTrack.setZ(mchParam.getZ()); - convertedTrack.setPhi(x2); - convertedTrack.setTanl(x3); - convertedTrack.setInvQPt(x4); - convertedTrack.setCharge(mchParam.getCharge()); - convertedTrack.setCovariances(covariances); - - return convertedTrack; - } - - int getDetElemId(int iDetElemNumber) - { - // make sure detector number is valid - if (iDetElemNumber < SNDetElemCh[0] || - iDetElemNumber >= SNDetElemCh[NMchChambers]) { - LOGF(fatal, "Invalid detector element number: %d", iDetElemNumber); - } - /// get det element number from ID - // get chamber and element number in chamber - int iCh = 0; - int iDet = 0; - for (int i = 1; i <= NMchChambers; i++) { - if (iDetElemNumber < SNDetElemCh[i]) { - iCh = i; - iDet = iDetElemNumber - SNDetElemCh[i - 1]; - break; - } - } - - // make sure detector index is valid - if (iCh <= 0 || iCh > NMchChambers || iDet >= NDetElemCh[iCh - 1]) { - LOGF(fatal, "Invalid detector element id: %d", MchDetElemNumberingBase * iCh + iDet); - } - - // add number of detectors up to this chamber - return MchDetElemNumberingBase * iCh + iDet; - } - - bool removeTrack(mch::Track& track) - { - // Refit track with re-aligned clusters - bool shouldRemoveTrack = false; - try { - trackFitter.fit(track, false); - } catch (std::exception const& e) { - shouldRemoveTrack = true; - return shouldRemoveTrack; - } - - auto itStartingParam = std::prev(track.rend()); - - while (true) { - - try { - trackFitter.fit(track, true, false, (itStartingParam == track.rbegin()) ? nullptr : &itStartingParam); - } catch (std::exception const&) { - shouldRemoveTrack = true; - break; - } - - double worstLocalChi2 = -1.0; - - track.tagRemovableClusters(0x1F, false); - - auto itWorstParam = track.end(); - - for (auto itParam = track.begin(); itParam != track.end(); ++itParam) { - if (itParam->getLocalChi2() > worstLocalChi2) { - worstLocalChi2 = itParam->getLocalChi2(); - itWorstParam = itParam; - } - } - - if (worstLocalChi2 < mImproveCutChi2) { - break; - } - - if (!itWorstParam->isRemovable()) { - shouldRemoveTrack = true; - track.removable(); - break; - } - - auto itNextParam = track.removeParamAtCluster(itWorstParam); - auto itNextToNextParam = (itNextParam == track.end()) ? itNextParam : std::next(itNextParam); - itStartingParam = track.rbegin(); - - if (track.getNClusters() < MinRemovableTrackClusters) { - shouldRemoveTrack = true; - break; - } - while (itNextToNextParam != track.end()) { - if (itNextToNextParam->getClusterPtr()->getChamberId() != itNextParam->getClusterPtr()->getChamberId()) { - itStartingParam = std::make_reverse_iterator(++itNextParam); - break; - } - ++itNextToNextParam; - } - } - - if (!shouldRemoveTrack) { - for (auto& param : track) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) - param.setParameters(param.getSmoothParameters()); - param.setCovariances(param.getSmoothCovariances()); - } - } - - return shouldRemoveTrack; - } - - template - void initCcdb(BC const& bc) - { - if (mRunNumber == bc.runNumber()) { - return; - } - - mRunNumber = bc.runNumber(); - std::map metadata; - auto soreor = o2::ccdb::BasicCCDBManager::getRunDuration(fCCDBApi, mRunNumber); - auto ts = soreor.first; - auto grpmag = fCCDBApi.retrieveFromTFileAny(configCcdb.cfgGrpMagPath, metadata, ts); - o2::base::Propagator::initFieldFromGRP(grpmag); - LOGF(info, "Set field for muons"); - VarManager::SetupMuonMagField(); - if (!o2::base::GeometryManager::isGeometryLoaded()) { - ccdbManager->get(configCcdb.cfgGeoPath); - } - mch::TrackExtrap::setField(); - mch::TrackExtrap::useExtrapV2(); - - // Load geometry information from CCDB/local - LOGF(info, "Loading reference aligned geometry from CCDB no later than %d", configMchRealign.cfgCcdbNoLaterThanRef.value); - ccdbManager->setCreatedNotAfter(configMchRealign.cfgCcdbNoLaterThanRef.value); // this timestamp has to be consistent with what has been used in reco - geoRef = ccdbManager->getForTimeStamp(configMchRealign.cfgGeoRefPath, bc.timestamp()); - ccdbManager->clearCache(configMchRealign.cfgGeoRefPath); - if (geoRef != nullptr) { - transformation = mch::geo::transformationFromTGeoManager(*geoRef); - } else { - LOGF(fatal, "Reference aligned geometry object is not available in CCDB at timestamp=%llu", bc.timestamp()); - } - for (int i = 0; i < NMchDetElems; i++) { - int iDEN = getDetElemId(i); - transformRef[iDEN] = transformation(iDEN); - } - - LOGF(info, "Loading new aligned geometry from CCDB no later than %d", configMchRealign.cfgCcdbNoLaterThanNew.value); - ccdbManager->setCreatedNotAfter(configMchRealign.cfgCcdbNoLaterThanNew.value); // make sure this timestamp can be resolved regarding the reference one - geoNew = ccdbManager->getForTimeStamp(configMchRealign.cfgGeoNewPath, bc.timestamp()); - ccdbManager->clearCache(configMchRealign.cfgGeoNewPath); - if (geoNew != nullptr) { - transformation = mch::geo::transformationFromTGeoManager(*geoNew); - } else { - LOGF(fatal, "New aligned geometry object is not available in CCDB at timestamp=%llu", bc.timestamp()); - } - for (int i = 0; i < NMchDetElems; i++) { - int iDEN = getDetElemId(i); - transformNew[iDEN] = transformation(iDEN); - } - - // Init magnetic field for MFT track extrapolation - auto* fieldB = dynamic_cast(TGeoGlobalMagField::Instance()->GetField()); - if (fieldB) { - std::array centerMft{0, 0, -61.4}; // Field at center of MFT - mBzAtMftCenter = fieldB->getBz(centerMft.data()); - // std::cout << "fieldB: " << (void*)fieldB << std::endl; - } - } - - void initMatchingFunctions() - { - using SVector2 = ROOT::Math::SVector; - using SVector4 = ROOT::Math::SVector; - using SVector5 = ROOT::Math::SVector; - - using SMatrix44 = ROOT::Math::SMatrix; - using SMatrix45 = ROOT::Math::SMatrix; - using SMatrix22 = ROOT::Math::SMatrix; - using SMatrix25 = ROOT::Math::SMatrix; - - // Define built-in matching functions - //________________________________________________________________________________ - mMatchingFunctionMap["matchALL"] = [](const o2::track::TrackParCovFwd& mchTrack, const o2::track::TrackParCovFwd& mftTrack) -> std::tuple { - // Match two tracks evaluating all parameters: X,Y, phi, tanl & q/pt - - SMatrix55Sym hK, vK; - SVector5 mK(mftTrack.getX(), mftTrack.getY(), mftTrack.getPhi(), - mftTrack.getTanl(), mftTrack.getInvQPt()), - rKKminus1; - const auto& globalMuonTrackParameters = mchTrack.getParameters(); - const auto& globalMuonTrackCovariances = mchTrack.getCovariances(); - vK(0, 0) = mftTrack.getCovariances()(0, 0); - vK(1, 1) = mftTrack.getCovariances()(1, 1); - vK(2, 2) = mftTrack.getCovariances()(2, 2); - vK(3, 3) = mftTrack.getCovariances()(3, 3); - vK(4, 4) = mftTrack.getCovariances()(4, 4); - hK(0, 0) = 1.0; - hK(1, 1) = 1.0; - hK(2, 2) = 1.0; - hK(3, 3) = 1.0; - hK(4, 4) = 1.0; - - // Covariance of residuals - SMatrix55Std invResCov = (vK + ROOT::Math::Similarity(hK, globalMuonTrackCovariances)); - invResCov.Invert(); - - // Update Parameters - rKKminus1 = mK - hK * globalMuonTrackParameters; // Residuals of prediction - - auto matchChi2Track = ROOT::Math::Similarity(rKKminus1, invResCov); - - // return chi2 and NDF - return {matchChi2Track, 5}; - }; - - //________________________________________________________________________________ - mMatchingFunctionMap["matchXYPhiTanl"] = [](const o2::track::TrackParCovFwd& mchTrack, const o2::track::TrackParCovFwd& mftTrack) -> std::tuple { - // Match two tracks evaluating positions & angles - - SMatrix45 hK; - SMatrix44 vK; - SVector4 mK(mftTrack.getX(), mftTrack.getY(), mftTrack.getPhi(), - mftTrack.getTanl()), - rKKminus1; - const auto& globalMuonTrackParameters = mchTrack.getParameters(); - const auto& globalMuonTrackCovariances = mchTrack.getCovariances(); - vK(0, 0) = mftTrack.getCovariances()(0, 0); - vK(1, 1) = mftTrack.getCovariances()(1, 1); - vK(2, 2) = mftTrack.getCovariances()(2, 2); - vK(3, 3) = mftTrack.getCovariances()(3, 3); - hK(0, 0) = 1.0; - hK(1, 1) = 1.0; - hK(2, 2) = 1.0; - hK(3, 3) = 1.0; - - // Covariance of residuals - SMatrix44 invResCov = (vK + ROOT::Math::Similarity(hK, globalMuonTrackCovariances)); - invResCov.Invert(); - - // Residuals of prediction - rKKminus1 = mK - hK * globalMuonTrackParameters; - - auto matchChi2Track = ROOT::Math::Similarity(rKKminus1, invResCov); - - // return chi2 and NDF - return {matchChi2Track, 4}; - }; - - //________________________________________________________________________________ - mMatchingFunctionMap["matchXY"] = [](const o2::track::TrackParCovFwd& mchTrack, const o2::track::TrackParCovFwd& mftTrack) -> std::tuple { - // Calculate Matching Chi2 - X and Y positions - - SMatrix25 hK; - SMatrix22 vK; - SVector2 mK(mftTrack.getX(), mftTrack.getY()), rKKminus1; - const auto& globalMuonTrackParameters = mchTrack.getParameters(); - const auto& globalMuonTrackCovariances = mchTrack.getCovariances(); - vK(0, 0) = mftTrack.getCovariances()(0, 0); - vK(1, 1) = mftTrack.getCovariances()(1, 1); - hK(0, 0) = 1.0; - hK(1, 1) = 1.0; - - // Covariance of residuals - SMatrix22 invResCov = (vK + ROOT::Math::Similarity(hK, globalMuonTrackCovariances)); - invResCov.Invert(); - - // Residuals of prediction - rKKminus1 = mK - hK * globalMuonTrackParameters; - auto matchChi2Track = ROOT::Math::Similarity(rKKminus1, invResCov); - - // return reduced chi2 - return {matchChi2Track, 2}; - }; - } - - void init(o2::framework::InitContext&) - { - // Load geometry - ccdbManager->setURL(configCcdb.cfgCcdbUrl); - ccdbManager->setCaching(true); - ccdbManager->setLocalObjectValidityChecking(); - fCCDBApi.init(configCcdb.cfgCcdbUrl); - mRunNumber = 0; - - // Configuration for track fitter - const auto& trackerParam = mch::TrackerParam::Instance(); - trackFitter.setBendingVertexDispersion(trackerParam.bendingVertexDispersion); - trackFitter.setChamberResolution(configMchRealign.cfgChamberResolutionX.value, configMchRealign.cfgChamberResolutionY.value); - trackFitter.smoothTracks(true); - trackFitter.useChamberResolution(); - mImproveCutChi2 = 2. * configMchRealign.cfgSigmaCutImprove.value * configMchRealign.cfgSigmaCutImprove.value; - - // Reset matching configuration, then populate only what we need. - hasActiveChi2Matching = false; - activeChi2FunctionName.clear(); - activeChi2MatchingPlaneZ = 0.; - - hasActiveMlMatching = false; - activeMlMatchingPlaneZ = 0.; - - if (configMatching.cfgCustomMatchingStrategy.value == 0) { - // Matching functions (custom chi2) - initMatchingFunctions(); - auto label = configChi2MatchingOptions.cfgChi2FunctionLabel.value; - auto funcName = configChi2MatchingOptions.cfgChi2FunctionName.value; - auto matchingPlaneZ = configChi2MatchingOptions.cfgChi2FunctionMatchingPlaneZ.value; - - if (!label.empty() && !funcName.empty()) { - hasActiveChi2Matching = true; - activeChi2FunctionName = funcName; - activeChi2MatchingPlaneZ = matchingPlaneZ; - } - } else { - // Matching ML models (custom ML) - // TODO : for now we use hard coded values since the current models use 1 pT bin - binsPtMl = {-1e-6, 1000.0}; - cutValues = {0.0}; - cutDirMl = {cuts_ml::CutNot}; - LabeledArray mycutsMl(cutValues.data(), 1, 1, std::vector{"pT bin 0"}, std::vector{"score"}); - - auto label = configMlOptions.cfgMlModelLabel.value; - auto modelPath = configMlOptions.cfgMlModelPathCcdb.value; - auto inputFeatures = configMlOptions.cfgMlInputFeatures.value; - auto modelName = configMlOptions.cfgMlModelName.value; - auto matchingPlaneZ = configMlOptions.cfgMlModelMatchingPlaneZ.value; - - if (!label.empty() && !modelPath.empty() && !inputFeatures.empty() && !modelName.empty()) { - activeMlResponse.configure(binsPtMl, mycutsMl, cutDirMl, 1); - activeMlResponse.setModelPathsCCDB(std::vector{modelName}, fCCDBApi, std::vector{modelPath}, configCcdb.cfgCcdbNoLaterThan.value); - activeMlResponse.cacheInputFeaturesIndices(inputFeatures); - activeMlResponse.init(); - - hasActiveMlMatching = true; - activeMlMatchingPlaneZ = matchingPlaneZ; - } - } - } - - template - bool pDcaCut(const T& mchTrack, const C& collision, double nSigmaPDCA) - { - static const double sigmaPDCA23 = 80.; - static const double sigmaPDCA310 = 54.; - static const double relPRes = 0.0004; - static const double slopeRes = 0.0005; - - constexpr double AbsorberEndZ = 505.; - constexpr double RadToDeg = 180. / o2::constants::math::PI; - double thetaAbs = std::atan(mchTrack.rAtAbsorberEnd() / AbsorberEndZ) * RadToDeg; - - // propagate muon track to vertex - auto mchTrackAtVertex = VarManager::PropagateMuon(mchTrack, collision, VarManager::kToVertex); - - // double pUncorr = mchTrack.p(); - double p = mchTrackAtVertex.getP(); - - double pDCA = mchTrack.pDca(); - double sigmaPDCA = (thetaAbs < ThetaAbsBoundaryDeg) ? sigmaPDCA23 : sigmaPDCA310; - double nrp = nSigmaPDCA * relPRes * p; - double pResEffect = sigmaPDCA / (1. - nrp / (1. + nrp)); - double slopeResEffect = SlopeResolutionZ * slopeRes * p; - double sigmaPDCAWithRes = std::sqrt(pResEffect * pResEffect + slopeResEffect * slopeResEffect); - return pDCA <= nSigmaPDCA * sigmaPDCAWithRes; - } - - template - bool isGoodMuon(const T& mchTrack, const C& collision, - double chi2Cut, - double pCut, - double pTCut, - std::array etaCut, - std::array rAbsCut, - double nSigmaPdcaCut) - { - // chi2 cut - if (mchTrack.chi2() > chi2Cut) { - return false; - } - - // momentum cut - if (mchTrack.p() < pCut) { - return false; // skip low-momentum tracks - } - - // transverse momentum cut - if (mchTrack.pt() < pTCut) { - return false; // skip low-momentum tracks - } - - // Eta cut - double eta = mchTrack.eta(); - if ((eta < etaCut[0] || eta > etaCut[1])) { - return false; - } - - // RAbs cut - double rAbs = mchTrack.rAtAbsorberEnd(); - if ((rAbs < rAbsCut[0] || rAbs > rAbsCut[1])) { - return false; - } - - // pDCA cut - return pDcaCut(mchTrack, collision, nSigmaPdcaCut); - } - - void storeFwdTrackCovariance(const SMatrix55Sym& cov) - { - const float sigX = std::sqrt(cov(0, 0)); - const float sigY = std::sqrt(cov(1, 1)); - const float sigPhi = std::sqrt(cov(2, 2)); - const float sigTgl = std::sqrt(cov(3, 3)); - const float sig1Pt = std::sqrt(cov(4, 4)); - const auto rhoXY = static_cast(128.f * cov(0, 1) / (sigX * sigY)); - const auto rhoPhiX = static_cast(128.f * cov(0, 2) / (sigPhi * sigX)); - const auto rhoPhiY = static_cast(128.f * cov(1, 2) / (sigPhi * sigY)); - const auto rhoTglX = static_cast(128.f * cov(0, 3) / (sigTgl * sigX)); - const auto rhoTglY = static_cast(128.f * cov(1, 3) / (sigTgl * sigY)); - const auto rhoTglPhi = static_cast(128.f * cov(2, 3) / (sigTgl * sigPhi)); - const auto rho1PtX = static_cast(128.f * cov(0, 4) / (sig1Pt * sigX)); - const auto rho1PtY = static_cast(128.f * cov(1, 4) / (sig1Pt * sigY)); - const auto rho1PtPhi = static_cast(128.f * cov(2, 4) / (sig1Pt * sigPhi)); - const auto rho1PtTgl = static_cast(128.f * cov(3, 4) / (sig1Pt * sigTgl)); - gmCandidateFwdTracksCov(sigX, sigY, sigPhi, sigTgl, sig1Pt, - rhoXY, rhoPhiY, rhoPhiX, rhoTglX, rhoTglY, rhoTglPhi, rho1PtX, rho1PtY, rho1PtPhi, rho1PtTgl); - } - - template - void fillBaseGmmCandFwdTrack(TMCH const& track, - TrackParExt const& trackPar, - int32_t gmmMchTrackId, - float chi2MatchMCHMFT, - float matchScoreMCHMFT) - { - const auto collisionId = track.collisionId(); - bool hasBcSlice = false; - std::array bcSlice{}; - if (collisionId < 0) { - const auto ambIt = mAmbBcSliceByFwdTrackId.find(track.globalIndex()); - if (ambIt != mAmbBcSliceByFwdTrackId.end()) { - bcSlice = ambIt->second; - hasBcSlice = true; - } - } - - gmCandidateFwdTracks( - collisionId, - track.trackType(), - trackPar.getX(), - trackPar.getY(), - trackPar.getZ(), - trackPar.getPhi(), - trackPar.getTgl(), - trackPar.getInvQPt(), - trackPar.getNClusters(), - track.pDca(), - track.rAtAbsorberEnd(), - trackPar.isRemovable(), - trackPar.getTrackChi2(), - track.chi2MatchMCHMID(), - chi2MatchMCHMFT, - matchScoreMCHMFT, - track.matchMFTTrackId(), - gmmMchTrackId, - track.mchBitMap(), - track.midBitMap(), - track.midBoards(), - track.trackTime(), - track.trackTimeRes()); - - storeFwdTrackCovariance(trackPar.getCovariances()); - if (hasBcSlice) { - gmAmbiguousFwdTracksReAlign(mGmmCandFwdTrackRowIndex, bcSlice.data()); - } - mGmmCandFwdTrackRowIndex += 1; - - mHasLastMchAmbiguousBcSlice = hasBcSlice; - if (hasBcSlice) { - mLastMchAmbiguousBcSlice = bcSlice; - } - } - - template - void fillCandidateFwdTrack(TMCH const& mchTrack, - TrackParExt const& mchPar, - int32_t gmmMchTrackId, - TMFT const& mftTrack, - TrackParExt const& mftPar, - const MatchingCandidate& candidate) - { - using o2::aod::fwdtrack::ForwardTrackTypeEnum; - using o2::aod::fwdtrackutils::propagationPoint; - - constexpr uint8_t CandidateTrackType = static_cast(ForwardTrackTypeEnum::GlobalForwardTrack); - - auto propmuonAtMft = fwdToMch(mchPar); - o2::mch::TrackExtrap::extrapToVertex(propmuonAtMft, - mftPar.getX(), - mftPar.getY(), - mftPar.getZ(), - mftPar.getSigma2X(), - mftPar.getSigma2Y()); - - const auto globalMuonRefit = o2::aod::fwdtrackutils::refitGlobalMuonCov(mchToFwd(propmuonAtMft), mftPar); - - const auto nClusters = static_cast(std::min(127, mchPar.getNClusters() + mftPar.getNClusters())); - - const float chi2 = static_cast(mchTrack.chi2()); - const int32_t collisionId = mchTrack.has_collision() ? mchTrack.collisionId() : -1; - bool hasBcSlice = false; - std::array bcSlice{}; - if (collisionId < 0) { - if (mHasLastMchAmbiguousBcSlice) { - bcSlice = mLastMchAmbiguousBcSlice; - hasBcSlice = true; - } else { - const auto ambIt = mAmbBcSliceByFwdTrackId.find(mchTrack.globalIndex()); - if (ambIt != mAmbBcSliceByFwdTrackId.end()) { - bcSlice = ambIt->second; - hasBcSlice = true; - } - } - } - - bool isRemovable = mchPar.isRemovable(); - - gmCandidateFwdTracks( - collisionId, - CandidateTrackType, - globalMuonRefit.getX(), - globalMuonRefit.getY(), - globalMuonRefit.getZ(), - globalMuonRefit.getPhi(), - globalMuonRefit.getTgl(), - globalMuonRefit.getInvQPt(), - nClusters, - mchTrack.pDca(), - mchTrack.rAtAbsorberEnd(), - isRemovable, - chi2, - mchTrack.chi2MatchMCHMID(), - static_cast(candidate.matchChi2), - static_cast(candidate.matchScore), - static_cast(mftTrack.globalIndex()), - gmmMchTrackId, - mchTrack.mchBitMap(), - mchTrack.midBitMap(), - mchTrack.midBoards(), - mchTrack.trackTime(), - mchTrack.trackTimeRes()); - - storeFwdTrackCovariance(globalMuonRefit.getCovariances()); - if (hasBcSlice) { - gmAmbiguousFwdTracksReAlign(mGmmCandFwdTrackRowIndex, bcSlice.data()); - } - mGmmCandFwdTrackRowIndex += 1; - } - - o2::track::TrackParCovFwd propagateToZMch(const o2::track::TrackParCovFwd& muon, const double z) - { - auto mchTrack = fwdToMch(muon); - - float absFront = -90.f; - float absBack = -505.f; - - if (muon.getZ() < absBack && z > absFront) { - // extrapolation through the absorber in the upstream direction - o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchTrack, z); - } else { - // all other cases - o2::mch::TrackExtrap::extrapToZCov(mchTrack, z); - } - - return mchToFwd(mchTrack); - } - - o2::track::TrackParCovFwd propagateToZMft(const o2::track::TrackParCovFwd& mftTrack, const double z) - { - o2::track::TrackParCovFwd trackExtrap{mftTrack}; - trackExtrap.propagateToZ(z, mBzAtMftCenter); - return trackExtrap; - } - - template - o2::track::TrackParCovFwd propagateToVertexMch(const TMCH& muon, - const C& collision) - { - auto mchTrack = fwdToMch(fwdtrackutils::getTrackParCovFwd(muon, muon)); - o2::mch::TrackExtrap::extrapToVertex(mchTrack, - collision.posX(), - collision.posY(), - collision.posZ(), - collision.covXX(), - collision.covYY()); - return mchToFwd(mchTrack); - } - - // tag muons based on the track quality and the track position at the front and back MFT planes - template - void getTaggedMuons(C const& collisions, - TMUON const& muonTracks, - std::vector& taggedMuons) - { - taggedMuons.clear(); - for (const auto& muonTrack : muonTracks) { - - // only consider MCH-MID matches - if (static_cast(muonTrack.trackType()) != MchMidTrackType) { - continue; - } - - // only select MCH-MID tracks associated to a collision - if (!muonTrack.has_collision()) { - continue; - } - - const auto& collision = collisions.rawIteratorAt(muonTrack.collisionId()); - - // select MCH tracks with strict quality cuts - if (!isGoodMuon(muonTrack, collision, - configMuonTagging.cfgMuonTaggingTrackChi2MchUp, - configMuonTagging.cfgMuonTaggingPMchLow, - configMuonTagging.cfgMuonTaggingPtMchLow, - {configMuonTagging.cfgMuonTaggingEtaMchLow, configMuonTagging.cfgMuonTaggingEtaMchUp}, - {configMuonTagging.cfgMuonTaggingRabsLow, configMuonTagging.cfgMuonTaggingRabsUp}, - configMuonTagging.cfgMuonTaggingPdcaUp)) { - continue; - } - - // propagate MCH track to the vertex - auto mchTrackAtVertex = propagateToVertexMch(muonTrack, collision); - - // propagate the track from the vertex to the first MFT plane - const auto& extrapToMFTfirst = propagateToZMch(mchTrackAtVertex, o2::mft::constants::mft::LayerZCoordinate()[0]); - double rFront = std::sqrt(extrapToMFTfirst.getX() * extrapToMFTfirst.getX() + extrapToMFTfirst.getY() * extrapToMFTfirst.getY()); - if (rFront < configMuonTagging.cfgMuonTaggingRadiusAtMftFrontLow.value || rFront > configMuonTagging.cfgMuonTaggingRadiusAtMftFrontUp.value) { - continue; - } - - // propagate the track from the vertex to the last MFT plane - const auto& extrapToMFTlast = propagateToZMch(mchTrackAtVertex, o2::mft::constants::mft::LayerZCoordinate()[9]); - double rBack = std::sqrt(extrapToMFTlast.getX() * extrapToMFTlast.getX() + extrapToMFTlast.getY() * extrapToMFTlast.getY()); - if (rBack < configMuonTagging.cfgMuonTaggingRadiusAtMftBackLow.value || rBack > configMuonTagging.cfgMuonTaggingRadiusAtMftBackUp.value) { - continue; - } - - int64_t muonTrackIndex = muonTrack.globalIndex(); - taggedMuons.emplace_back(muonTrackIndex); - } - } - - template - bool isMftMchTimeCompatible(EVT const& collisions, - BC const& bcs, - TMUON const& mchTrack, - TMFT const& mftTrack) - { - if (!mchTrack.has_collision() || !mftTrack.has_collision()) { - return false; - } - - const auto& collMch = collisions.rawIteratorAt(mchTrack.collisionId()); - const auto& bcMch = bcs.rawIteratorAt(collMch.bcId()); - const auto& collMft = collisions.rawIteratorAt(mftTrack.collisionId()); - const auto& bcMft = bcs.rawIteratorAt(collMft.bcId()); - - int64_t deltaBc = static_cast(bcMft.globalBC()) - static_cast(bcMch.globalBC()); - double deltaBcNS = o2::constants::lhc::LHCBunchSpacingNS * deltaBc; - double deltaTrackTime = mftTrack.trackTime() - mchTrack.trackTime() + deltaBcNS; - double trackTimeResTot = mftTrack.trackTimeRes() + mchTrack.trackTimeRes(); - - return std::fabs(deltaTrackTime) <= trackTimeResTot; - } - - template - void prepareMatchingCandidates(EVT const& collisions, - BC const& bcs, - TMUON const& muonTracks, - TMFT const& mftTracks, - MyMFTCovariances const& mftCovs) - { - mMftTrackPars.clear(); - mMchTrackPars.clear(); - mMatchingCandidates.clear(); - - LOGF(info, "Filling matching candidate tables"); - - for (const auto& muonTrack : muonTracks) { - if (static_cast(muonTrack.trackType()) <= GlobalTrackTypeMax) { - continue; - } - auto mchTrackIndex = muonTrack.globalIndex(); - - // initialize the MCH track parameters, which will be updated by the realignment if enabled - mMchTrackPars.try_emplace(mchTrackIndex, TrackParExt(fwdtrackutils::getTrackParCovFwd(muonTrack, muonTrack), muonTrack.nClusters())); - } - - for (const auto& mftTrack : mftTracks) { - auto mftTrackIndex = mftTrack.globalIndex(); - - // initialize the MFT track parameters, which will be updated by the alignment corrections if enabled - if (mftTrackCovs.contains(mftTrackIndex) && !mMftTrackPars.contains(mftTrackIndex)) { - auto const& mftTrackCov = mftCovs.rawIteratorAt(mftTrackCovs[mftTrackIndex]); - mMftTrackPars.emplace(mftTrackIndex, TrackParExt(fwdtrackutils::getTrackParCovFwd(mftTrack, mftTrackCov), mftTrack.nClusters())); - } - } - - // fill matching candidates table - if (!configMatching.cfgMatchAllTracks.value) { - // collect global MFT-MCH or MFT-MCH-MID tracks and associate them to the corresponding MCH(-MID) track - for (const auto& muonTrack : muonTracks) { - // skip MCH or MCH-MID tracks - if (static_cast(muonTrack.trackType()) > GlobalTrackTypeMax) { - continue; - } - - auto const& mchTrack = muonTrack.template matchMCHTrack_as(); - int64_t mchTrackIndex = mchTrack.globalIndex(); - auto const& mftTrack = muonTrack.template matchMFTTrack_as(); - int64_t mftTrackIndex = mftTrack.globalIndex(); - - if (!mftTrackCovs.contains(mftTrackIndex)) { - continue; - } - - mMatchingCandidates[mchTrackIndex].emplace_back(MatchingCandidate{ - .muonTrackId = muonTrack.globalIndex(), - .mftTrackId = mftTrackIndex, - .matchScore = muonTrack.matchScoreMCHMFT(), - .matchChi2 = muonTrack.chi2MatchMCHMFT()}); - } - } else { - // build matching candidates from all time-compatible MFT-MCH pairs - for (const auto& muonTrack : muonTracks) { - if (static_cast(muonTrack.trackType()) <= GlobalTrackTypeMax) { - continue; - } - auto mchTrackIndex = muonTrack.globalIndex(); - for (const auto& mftTrack : mftTracks) { - if (!isMftMchTimeCompatible(collisions, bcs, muonTrack, mftTrack)) { - continue; - } - if (!mftTrackCovs.contains(mftTrack.globalIndex())) { - continue; - } - - mMatchingCandidates[mchTrackIndex].emplace_back(MatchingCandidate{ - .mftTrackId = mftTrack.globalIndex()}); - } - } - } - - // sort the vectors of matching candidates in ascending order based on the matching chi2 value - auto compareMatchingChi2 = [](const MatchingCandidate& track1, const MatchingCandidate& track2) -> bool { - return (track1.matchChi2 < track2.matchChi2); - }; - - for (auto& [mchIndex, candidatesVector] : mMatchingCandidates) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) - std::sort(candidatesVector.begin(), candidatesVector.end(), compareMatchingChi2); - } - } - - template - o2::track::TrackParCovFwd transformMft(TMFT& mftTrack, TMFTCOV const& mftTrackCov) - { - auto track = fwdToMch(fwdtrackutils::getTrackParCovFwd(mftTrack, mftTrackCov)); - - double z = track.getZ(); - // double dZ = zMCH - z; - double x = track.getNonBendingCoor(); - double y = track.getBendingCoor(); - double xSlope = track.getNonBendingSlope(); - double ySlope = track.getBendingSlope(); - - double xSlopeCorrection = (y > 0) ? configMftAlignmentCorrections.cfgMFTAlignmentCorrXSlopeTop : configMftAlignmentCorrections.cfgMFTAlignmentCorrXSlopeBottom; - double xCorrection = xSlopeCorrection * z + - ((y > 0) ? configMftAlignmentCorrections.cfgMFTAlignmentCorrXOffsetTop : configMftAlignmentCorrections.cfgMFTAlignmentCorrXOffsetBottom); - double xNew = x + xCorrection; - double xSlopeNew = xSlope + xSlopeCorrection; - - track.setNonBendingCoor(xNew); - track.setNonBendingSlope(xSlopeNew); - - double ySlopeCorrection = (y > 0) ? configMftAlignmentCorrections.cfgMFTAlignmentCorrYSlopeTop : configMftAlignmentCorrections.cfgMFTAlignmentCorrYSlopeBottom; - double yCorrection = ySlopeCorrection * z + - ((y > 0) ? configMftAlignmentCorrections.cfgMFTAlignmentCorrYOffsetTop : configMftAlignmentCorrections.cfgMFTAlignmentCorrYOffsetBottom); - track.setBendingCoor(y + yCorrection); - track.setBendingSlope(ySlope + ySlopeCorrection); - - return mchToFwd(track); - } - - template - void runMftRealignment(TMFTs const& mftTracks, TMFTCOVs const& mftCovs) - { - for (const auto& mftTrack : mftTracks) { - auto mftTrackIndex = mftTrack.globalIndex(); - if (!mftTrackCovs.contains(mftTrackIndex)) { - continue; - } - - auto const& mftTrackCov = mftCovs.rawIteratorAt(mftTrackCovs[mftTrackIndex]); - mMftTrackPars[mftTrackIndex] = transformMft(mftTrack, mftTrackCov); - } - } - - template - void runMuonRealignment(TMuons const& muons, TMuonCls const& clusters) - { - // Loop over forward tracks - for (auto const& muon : muons) { - int mchIndex = muon.globalIndex(); - // skip global forward matches - if (muon.trackType() > GlobalTrackTypeMax) { - continue; - } - - // continue; - - auto mchTrackParIt = mMchTrackPars.find(mchIndex); - if (mchTrackParIt == mMchTrackPars.end()) { - continue; - } - - auto clustersSliced = clusters.sliceBy(perMuon, muon.globalIndex()); // Slice clusters by muon id - mch::Track convertedTrack = mch::Track(); // Temporary variable to store re-aligned clusters - - int clIndex = -1; - // Get re-aligned clusters associated to current track - for (auto const& cluster : clustersSliced) { - clIndex += 1; - - auto* clusterMCH = new mch::Cluster(); - - math_utils::Point3D local; - math_utils::Point3D master; - master.SetXYZ(cluster.x(), cluster.y(), cluster.z()); - - // Transformation from reference geometry frame to new geometry frame - transformRef[cluster.deId()].MasterToLocal(master, local); - transformNew[cluster.deId()].LocalToMaster(local, master); - - clusterMCH->x = master.x(); - clusterMCH->y = master.y(); - clusterMCH->z = master.z(); - - const 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); - // LOGF(debug, "Track %d, cluster DE%d: x:%g y:%g z:%g", muon.globalIndex(), cluster.deId(), cluster.x(), cluster.y(), cluster.z()); - // LOGF(debug, "Track %d, re-aligned cluster DE%d: x:%g y:%g z:%g", muonRealignId, cluster.deId(), clusterMCH->getX(), clusterMCH->getY(), clusterMCH->getZ()); - } - - // Refit the re-aligned track - int removable = 0; - if (convertedTrack.getNClusters() != 0) { - removable = removeTrack(convertedTrack); - } else { - LOGF(fatal, "Muon track %d has no associated clusters.", muon.globalIndex()); - } - - // Get the re-aligned track parameter: track param at the first cluster - mch::TrackParam trackParam = mch::TrackParam(convertedTrack.first()); - - // Convert MCH track to FWD track and store new parameters after realignment - mchTrackParIt->second = mchToFwd(mch::TrackParam(convertedTrack.first())); - mchTrackParIt->second.setTrackChi2(trackParam.getTrackChi2() / convertedTrack.getNDF()); - mchTrackParIt->second.setNClusters(convertedTrack.getNClusters()); - if (removable) { - mchTrackParIt->second.setRemovable(); - } - } - } - - void runChi2Matching(const std::string& funcName, - float matchingPlaneZ, - const MatchingCandidates& matchingCandidates, - MatchingCandidates& newMatchingCandidates) - { - newMatchingCandidates.clear(); - - std::string funcNameEffective = funcName; - float matchingPlaneZEffective = matchingPlaneZ; - if (funcName == "prod") { - funcNameEffective = "matchALL"; - matchingPlaneZEffective = MatchingPlaneDefaultZ; - } - - if (!mMatchingFunctionMap.contains(funcNameEffective)) { - return; - } - auto matchingFunc = mMatchingFunctionMap.at(funcNameEffective); - - for (const auto& [mchIndex, candidatesVector] : matchingCandidates) { - - // get the tracks parameters, which have been updated by the realignment if enabled - const auto mchTrackParIt = mMchTrackPars.find(mchIndex); - if (mchTrackParIt == mMchTrackPars.end()) { - continue; - } - - for (const auto& candidate : candidatesVector) { - auto mftTrackParIt = mMftTrackPars.find(candidate.mftTrackId); - if (mftTrackParIt == mMftTrackPars.end()) { - continue; - } - - auto mftTrackProp = mftTrackParIt->second.asTrackParCovFwd(); - auto mchTrackProp = mchTrackParIt->second.asTrackParCovFwd(); - - if (matchingPlaneZEffective < 0.) { - mftTrackProp = propagateToZMft(mftTrackProp, matchingPlaneZ); - mchTrackProp = propagateToZMch(mchTrackProp, matchingPlaneZ); - } - - auto matchResult = matchingFunc(mchTrackProp, mftTrackProp); - float matchChi2 = std::get<0>(matchResult); - - newMatchingCandidates[mchIndex].emplace_back(MatchingCandidate{ - .muonTrackId = candidate.muonTrackId, - .mftTrackId = candidate.mftTrackId, - .matchScore = -1, - .matchChi2 = matchChi2}); - } - } - - auto compareMatchingChi2 = [](const MatchingCandidate& track1, const MatchingCandidate& track2) -> bool { - return (track1.matchChi2 < track2.matchChi2); - }; - - for (auto& [mchIndex, globalTracksVector] : newMatchingCandidates) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) - std::sort(globalTracksVector.begin(), globalTracksVector.end(), compareMatchingChi2); - - int ranking = 1; - for (auto& candidate : globalTracksVector) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) - candidate.matchRanking = ranking; - ranking += 1; - } - } - } - - template - void runMlMatching(C const& collisions, - TMUON const& muonTracks, - TMFT const& mftTracks, - o2::analysis::MlResponseMFTMuonMatch& mlResponse, - float matchingPlaneZ, - const MatchingCandidates& matchingCandidates, - MatchingCandidates& newMatchingCandidates) - { - newMatchingCandidates.clear(); - for (const auto& [mchIndex, candidatesVector] : matchingCandidates) { - auto const& mchTrack = muonTracks.rawIteratorAt(mchIndex); - if (!mchTrack.has_collision()) { - continue; - } - - auto collision = collisions.rawIteratorAt(mchTrack.collisionId()); - - // get the tracks parameters, which have been updated by the realignment if enabled - auto mchTrackParIt = mMchTrackPars.find(mchIndex); - if (mchTrackParIt == mMchTrackPars.end()) { - continue; - } - - for (const auto& candidate : candidatesVector) { - auto const& muonTrack = (candidate.muonTrackId >= 0) ? muonTracks.rawIteratorAt(candidate.muonTrackId) : mchTrack; - auto const& mftTrack = mftTracks.rawIteratorAt(candidate.mftTrackId); - auto mftTrackParIt = mMftTrackPars.find(candidate.mftTrackId); - if (mftTrackParIt == mMftTrackPars.end()) { - continue; - } - - auto mftTrackProp = mftTrackParIt->second.asTrackParCovFwd(); - auto mchTrackProp = mchTrackParIt->second.asTrackParCovFwd(); - - if (matchingPlaneZ < 0.) { - mftTrackProp = propagateToZMft(mftTrackProp, matchingPlaneZ); - mchTrackProp = propagateToZMch(mchTrackProp, matchingPlaneZ); - } - - std::vector output; - std::vector inputML = mlResponse.getInputFeatures(muonTrack, mftTrack, mchTrack, mftTrackProp, mchTrackProp, collision); - mlResponse.isSelectedMl(inputML, 0, output); - float matchScore = output[0]; - - newMatchingCandidates[mchIndex].emplace_back(MatchingCandidate{ - .muonTrackId = candidate.muonTrackId, - .mftTrackId = candidate.mftTrackId, - .matchScore = matchScore, - .matchChi2 = -1}); - } - } - - auto compareMatchingScore = [](const MatchingCandidate& track1, const MatchingCandidate& track2) -> bool { - return (track1.matchScore > track2.matchScore); - }; - - for (auto& [mchIndex, globalTracksVector] : newMatchingCandidates) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) - std::sort(globalTracksVector.begin(), globalTracksVector.end(), compareMatchingScore); - - int ranking = 1; - for (auto& candidate : globalTracksVector) { // o2-linter: disable=const-ref-in-for-loop (object is modified in loop) - candidate.matchRanking = ranking; - ranking += 1; - } - } - } - - template - void processMatchingCandidates(C const& collisions, - TMUON const& muonTracks, - TMFT const& mftTracks, - CMFT const& mftCovs, - aod::FwdTrkCls const& clusters) - { - if (configMchRealign.cfgEnableMCHRealign.value) { - runMuonRealignment(muonTracks, clusters); - } - - if (configMftAlignmentCorrections.cfgEnableMftAlignmentCorrections) { - runMftRealignment(mftTracks, mftCovs); - } - - std::vector taggedMuons; - getTaggedMuons(collisions, muonTracks, taggedMuons); - - if (configMatching.cfgCustomMatchingStrategy.value == 0) { - if (hasActiveChi2Matching) { - MatchingCandidates newMatchingCandidates; - runChi2Matching(activeChi2FunctionName, activeChi2MatchingPlaneZ, mMatchingCandidates, newMatchingCandidates); - fillMatchingCandidates(newMatchingCandidates, taggedMuons); - } - } else { - if (hasActiveMlMatching) { - MatchingCandidates newMatchingCandidates; - runMlMatching(collisions, muonTracks, mftTracks, activeMlResponse, activeMlMatchingPlaneZ, mMatchingCandidates, newMatchingCandidates); - fillMatchingCandidates(newMatchingCandidates, taggedMuons); - } - } - } - - void fillMatchingCandidates(const MatchingCandidates& matchingCandidates, - const std::vector& taggedMuons) - { - for (const auto& [mchIndex, candidates] : matchingCandidates) { - if (candidates.empty()) { - continue; - } - - bool isTagged = std::find(taggedMuons.begin(), taggedMuons.end(), mchIndex) != taggedMuons.end(); - - std::vector storedCandidates; - int nStored = 0; - for (const auto& candidate : candidates) { - if (configMatching.cfgMaxCandidatesPerMchTrack.value >= 0 && nStored >= configMatching.cfgMaxCandidatesPerMchTrack.value) { - break; - } - - int32_t candidateIndex = mMatchCandidateCounter; - globalMuonMatchCandidates( - mchIndex, - candidate.mftTrackId, - static_cast(candidate.matchChi2), - static_cast(candidate.matchScore), - static_cast(candidate.matchRanking), - isTagged); - mMatchCandidateCounter += 1; - - mMchTrackToCandidateIndices[mchIndex].push_back(candidateIndex); - storedCandidates.push_back(candidate); - nStored += 1; - } - - if (!storedCandidates.empty()) { - mMchTrackMatchingCandidates[mchIndex] = std::move(storedCandidates); - } - } - } - - int32_t countStoredCandidatesForMchTrack(int64_t mchTrackIndex) const - { - const auto candidateIterator = mMchTrackMatchingCandidates.find(mchTrackIndex); - if (candidateIterator == mMchTrackMatchingCandidates.end()) { - return 0; - } - return static_cast(candidateIterator->second.size()); - } - - template - void fillGmmCandidateFwdTracks(TMUON const& muonTracks, - TMFT const& mftTracks, - aod::AmbiguousFwdTracks const& ambFwdTracks) - { - mFwdTrackToGmmCandTrkIndex.clear(); - mGmmCandFwdTrackRowIndex = 0; - mHasLastMchAmbiguousBcSlice = false; - mAmbBcSliceByFwdTrackId.clear(); - for (const auto& ambFwdTrack : ambFwdTracks) { - const auto bcIds = ambFwdTrack.bcIds(); - mAmbBcSliceByFwdTrackId[ambFwdTrack.fwdtrackId()] = {bcIds[0], bcIds[1]}; - } - - // First pass: assign GMMCANDTRK row indices for MCH/MCH-MID base entries so that - // MCHTrackId can be remapped consistently even when global muons appear first in FwdTracks. - int32_t nextGmmCandTrkIndex = 0; - for (const auto& track : muonTracks) { - const int trackType = static_cast(track.trackType()); - if (trackType > GlobalTrackTypeMax) { - mFwdTrackToGmmCandTrkIndex[track.globalIndex()] = nextGmmCandTrkIndex; - nextGmmCandTrkIndex += 1 + countStoredCandidatesForMchTrack(track.globalIndex()); - } else if (configMatching.cfgIncludeGlobalMuonsInFwdTracks.value) { - nextGmmCandTrkIndex += 1; - } - } - - // Second pass: fill GMMCANDTRK/GMMCANDTRKCOV in FwdTracks order. - for (const auto& track : muonTracks) { - const int trackType = static_cast(track.trackType()); - - if (trackType > GlobalTrackTypeMax) { - mHasLastMchAmbiguousBcSlice = false; - const int64_t mchTrackIndex = track.globalIndex(); - const int32_t gmmMchTrackId = mFwdTrackToGmmCandTrkIndex.at(mchTrackIndex); - - const auto candidateIterator = mMchTrackMatchingCandidates.find(mchTrackIndex); - auto mchTrackParIt = mMchTrackPars.find(mchTrackIndex); - if (mchTrackParIt == mMchTrackPars.end()) { - // fill muon tracks table with original parameters - const TrackParExt trackPar{fwdtrackutils::getTrackParCovFwd(track, track)}; - fillBaseGmmCandFwdTrack(track, trackPar, gmmMchTrackId, -1.f, -1.f); - } else { - // fill muon tracks table with realignment parameters - fillBaseGmmCandFwdTrack(track, mchTrackParIt->second, gmmMchTrackId, -1.f, -1.f); - } - - if (candidateIterator != mMchTrackMatchingCandidates.end()) { - for (const auto& candidate : candidateIterator->second) { - auto mftTrackParIt = mMftTrackPars.find(candidate.mftTrackId); - if (mftTrackParIt != mMftTrackPars.end()) { - const auto& mftTrack = mftTracks.rawIteratorAt(candidate.mftTrackId); - fillCandidateFwdTrack(track, mchTrackParIt->second, gmmMchTrackId, mftTrack, mftTrackParIt->second, candidate); - } - } - } - } - - if (configMatching.cfgIncludeGlobalMuonsInFwdTracks.value && trackType <= GlobalTrackTypeMax) { - int32_t gmmMchTrackId = -1; - const auto mchIterator = mFwdTrackToGmmCandTrkIndex.find(track.matchMCHTrackId()); - if (mchIterator != mFwdTrackToGmmCandTrkIndex.end()) { - gmmMchTrackId = mchIterator->second; - } - TrackParExt parExt(fwdtrackutils::getTrackParCovFwd(track, track)); - fillBaseGmmCandFwdTrack(track, - parExt, - gmmMchTrackId, - track.chi2MatchMCHMFT(), - track.matchScoreMCHMFT()); - } - } - } - - template - void fillFwdTrkMatchCands(TMUON const& muonTracks) - { - std::vector empty{}; - for (const auto& muonTrack : muonTracks) { - if (static_cast(muonTrack.trackType()) <= GlobalTrackTypeMax) { - fwdTrkMatchCands(empty); - continue; - } - - const int64_t mchTrackIndex = muonTrack.globalIndex(); - const auto matchIterator = mMchTrackToCandidateIndices.find(mchTrackIndex); - if (matchIterator == mMchTrackToCandidateIndices.end() || matchIterator->second.empty()) { - fwdTrkMatchCands(empty); - } else { - fwdTrkMatchCands(matchIterator->second); - } - } - } - - void processData(MyEvents const& collisions, - aod::BCsWithTimestamps const& bcs, - MyMuons const& muonTracks, - MyMFTs const& mftTracks, - MyMFTCovariances const& mftCovs, - aod::FwdTrkCls const& clusters, - aod::AmbiguousFwdTracks const& ambFwdTracks) - { - auto bc = bcs.begin(); - initCcdb(bc); - - LOGF(info, "Filling MFT cov"); - mftTrackCovs.clear(); - for (const auto& mftTrackCov : mftCovs) { - mftTrackCovs[mftTrackCov.matchMFTTrackId()] = mftTrackCov.globalIndex(); - } - - mMatchCandidateCounter = 0; - mMchTrackToCandidateIndices.clear(); - mMchTrackMatchingCandidates.clear(); - mFwdTrackToGmmCandTrkIndex.clear(); - - LOGF(info, "Preparing candidates"); - prepareMatchingCandidates(collisions, bcs, muonTracks, mftTracks, mftCovs); - - LOGF(info, "Processing candidates"); - processMatchingCandidates(collisions, muonTracks, mftTracks, mftCovs, clusters); - - LOGF(info, "Filling tables"); - // fill table with track/candidates index mapping - fillFwdTrkMatchCands(muonTracks); - // fill track tables - fillGmmCandidateFwdTracks(muonTracks, mftTracks, ambFwdTracks); - } - - PROCESS_SWITCH(GlobalMuonMatching, processData, "processData", true); -}; - -// Extends the fwdtracksrealign table with expression columns -struct GlobalMuonMatchingSpawner { - Spawns realignFwdTrksCov; - Spawns realignFwdTrks; - void init(InitContext const&) {} -}; - -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) -{ - return WorkflowSpec{ - adaptAnalysisTask(cfgc), - adaptAnalysisTask(cfgc)}; -}; From 0ec5d32535e82c847a48ebfc63d4b17ebbfee89e Mon Sep 17 00:00:00 2001 From: Maurice Coquet Date: Sun, 13 Sep 2026 10:31:57 +0200 Subject: [PATCH 7/9] removing deprecated PV shift options --- PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx | 17 ++--------------- PWGDQ/TableProducer/tableMaker_withAssoc.cxx | 9 ++------- 2 files changed, 4 insertions(+), 22 deletions(-) diff --git a/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx b/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx index 2396d797536..374fe6e666f 100644 --- a/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx +++ b/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx @@ -234,9 +234,6 @@ struct TableMakerMC { Configurable fConfigCcdbUrl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; Configurable fGeoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"}; Configurable fGrpMagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; - Configurable fZShiftPath{"zShiftPath", "Users/m/mcoquet/ZShift", "CCDB path for z shift to apply to forward tracks"}; - Configurable fUseRemoteZShift{"cfgUseRemoteZShift", false, "Enable getting Zshift from ccdb"}; - Configurable fManualZShift{"cfgManualZShift", 0.f, "Manual value for the Zshift for muons."}; Configurable fGrpMagPathRun2{"grpmagPathRun2", "GLO/GRP/GRP", "CCDB path of the GRPObject (Usage for Run 2)"}; Configurable timestampCCDB{"timestampCCDB", -1, "timestamp of the ONNX file for ML model used to query in CCDB"}; } fConfigCCDB; @@ -1138,7 +1135,7 @@ struct TableMakerMC { } // recalculate pDca / DCA and global muon kinematics // kMuonPDca is always taken from MCH (standalone or the MCH matched to a global) - if (static_cast(muon.trackType()) < 2) { + if (static_cast(muon.trackType()) <= 2) { auto muontrack = muon.template matchMCHTrack_as(); VarManager::FillTrackCollision(muontrack, collision); if (fConfigVariousOptions.fRefitGlobalMuon) { @@ -1275,7 +1272,7 @@ struct TableMakerMC { // recalculate pDca / DCA and global muon kinematics // kMuonPDca is always taken from MCH (standalone or the MCH matched to a global) int globalClusters = muon.nClusters(); - if (static_cast(muon.trackType()) < 2) { + if (static_cast(muon.trackType()) <= 2) { auto muontrack = muon.template matchMCHTrack_as(); VarManager::FillTrackCollision(muontrack, collision); if (fConfigVariousOptions.fRefitGlobalMuon) { @@ -1337,16 +1334,6 @@ struct TableMakerMC { o2::base::Propagator::initFieldFromGRP(fGrpMag); VarManager::SetMagneticField(fGrpMag->getNominalL3Field()); } - if (fConfigCCDB.fUseRemoteZShift) { - auto* fZShift = fCCDB->getForTimeStamp>(fConfigCCDB.fZShiftPath, bcs.begin().timestamp()); - if (fZShift != nullptr && !fZShift->empty()) { - VarManager::SetZShift((*fZShift)[0]); - } else { - LOG(fatal) << "Could not retrieve Z-shift value from CCDB"; - } - } else { - VarManager::SetZShift(fConfigCCDB.fManualZShift.value); - } if (fConfigVariousOptions.fPropMuon) { VarManager::SetupMuonMagField(); } diff --git a/PWGDQ/TableProducer/tableMaker_withAssoc.cxx b/PWGDQ/TableProducer/tableMaker_withAssoc.cxx index 25fda45923b..d957501b2d8 100644 --- a/PWGDQ/TableProducer/tableMaker_withAssoc.cxx +++ b/PWGDQ/TableProducer/tableMaker_withAssoc.cxx @@ -290,7 +290,6 @@ struct TableMaker { Configurable fConfigGrpMagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; Configurable fFwdShiftPath{"fwdShiftPath", "Users/m/mcoquet/ZShift", "CCDB path for the shift to apply to forward tracks: 1 (z), 3 (x,y,z), or 10 (x,y,z,slopeX,slopeY for top then bottom; slopes unused)"}; Configurable fUseRemoteFwdShift{"cfgUseRemoteFwdShift", false, "Enable getting the forward track shift from ccdb"}; - Configurable fManualZShift{"cfgManualZShift", 0.f, "Manual value for the Zshift for muons."}; Configurable fConfigGrpMagPathRun2{"grpmagPathRun2", "GLO/GRP/GRP", "CCDB path of the GRPObject (Usage for Run 2)"}; } fConfigCCDB; @@ -1692,7 +1691,7 @@ struct TableMaker { } // recalculate pDca / DCA and global muon kinematics // kMuonPDca is always taken from MCH (standalone or the MCH matched to a global) - if (static_cast(muon.trackType()) < 2) { + if (static_cast(muon.trackType()) <= 2) { auto muontrack = muon.template matchMCHTrack_as(); VarManager::FillTrackCollision(muontrack, collision); if (fConfigVariousOptions.fRefitGlobalMuon) { @@ -1700,8 +1699,6 @@ struct TableMaker { continue; } auto mfttrack = muon.template matchMFTTrack_as(); - // NOTE: the MFT track originally associated to the MUON track is currently used in the global muon refit - // Should MUON - MFT time ambiguities be taken into account ? // Helix DCA (kMuonDCAx/y) is filled from the refitted parameters inside FillGlobalMuonRefit(Cov) if constexpr (static_cast(TMFTFillMap & VarManager::ObjTypes::MFTCov)) { auto const& mfttrackcov = mfCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]); @@ -1795,7 +1792,7 @@ struct TableMaker { // recalculate pDca / DCA and global muon kinematics // kMuonPDca is always taken from MCH (standalone or the MCH matched to a global) int globalClusters = muon.nClusters(); - if (static_cast(muon.trackType()) < 2) { + if (static_cast(muon.trackType()) <= 2) { auto muontrack = muon.template matchMCHTrack_as(); VarManager::FillTrackCollision(muontrack, collision); if (fConfigVariousOptions.fRefitGlobalMuon) { @@ -1894,8 +1891,6 @@ struct TableMaker { } else { LOG(fatal) << "Unexpected number of shift values from CCDB: " << fFwdShift->size() << ", expected 1 (z), 3 (x, y, z) or 10 (top/bottom x,y,z + slopes)"; } - } else { - VarManager::SetZShift(fConfigCCDB.fManualZShift.value); } if (fConfigHistOutput.fConfigFillBcStat) { mLHCIFdata = fCCDB->getSpecific("GLO/Config/GRPLHCIF", bcs.begin().timestamp()); From fbd1553be6f0450566d2a2c428631af17f6d58d5 Mon Sep 17 00:00:00 2001 From: Maurice Coquet Date: Sun, 13 Sep 2026 11:26:03 +0200 Subject: [PATCH 8/9] fixing few errors --- PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx | 15 +++++++-------- PWGDQ/TableProducer/tableMaker_withAssoc.cxx | 2 +- 2 files changed, 8 insertions(+), 9 deletions(-) diff --git a/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx b/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx index 374fe6e666f..f369795251c 100644 --- a/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx +++ b/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx @@ -144,7 +144,7 @@ template void PrintBitMap(TMap map, int nbits) { for (int i = 0; i < nbits; i++) { - cout << ((map & (TMap(1) << i)) > 0 ? "1" : "0"); + LOG(info) << ((map & (TMap(1) << i)) > 0 ? "1" : "0"); } } */ @@ -577,17 +577,16 @@ struct TableMakerMC { /*if ((std::abs(mctrack.pdgCode())>400 && std::abs(mctrack.pdgCode())<599) || (std::abs(mctrack.pdgCode())>4000 && std::abs(mctrack.pdgCode())<5999) || (mcflags > 0)) { - cout << ">>>>>>>>>>>>>>>>>>>>>>> track idx / pdg / process / status code / HEPMC status / primary : " - << mctrack.globalIndex() << " / " << mctrack.pdgCode() << " / " - << mctrack.getProcess() << " / " << mctrack.getGenStatusCode() << " / " << mctrack.getHepMCStatusCode() << " / " << mctrack.isPhysicalPrimary() << endl; - cout << ">>>>>>>>>>>>>>>>>>>>>>> track bitmap: "; + LOG(info) << ">>>>>>>>>>>>>>>>>>>>>>> track idx / pdg / process / status code / HEPMC status / primary : " + << mctrack.globalIndex() << " / " << mctrack.pdgCode() << " / " + << mctrack.getProcess() << " / " << mctrack.getGenStatusCode() << " / " << mctrack.getHepMCStatusCode() << " / " << mctrack.isPhysicalPrimary(); + LOG(info) << ">>>>>>>>>>>>>>>>>>>>>>> track bitmap: "; PrintBitMap(mcflags, 16); - cout << endl; if (mctrack.has_mothers()) { for (const auto& m : mctrack.mothersIds()) { if (m < mcTracks.size()) { // protect against bad mother indices auto aMother = mcTracks.rawIteratorAt(m); - cout << "<<<<<< mother idx / pdg: " << m << " / " << aMother.pdgCode() << endl; + LOG(info) << "<<<<<< mother idx / pdg: " << m << " / " << aMother.pdgCode(); } } } @@ -597,7 +596,7 @@ struct TableMakerMC { if (d < mcTracks.size()) { // protect against bad daughter indices auto aDaughter = mcTracks.rawIteratorAt(d); - cout << "<<<<<< daughter idx / pdg: " << d << " / " << aDaughter.pdgCode() << endl; + LOG(info) << "<<<<<< daughter idx / pdg: " << d << " / " << aDaughter.pdgCode(); } } } diff --git a/PWGDQ/TableProducer/tableMaker_withAssoc.cxx b/PWGDQ/TableProducer/tableMaker_withAssoc.cxx index d957501b2d8..b5b6d5d76f6 100644 --- a/PWGDQ/TableProducer/tableMaker_withAssoc.cxx +++ b/PWGDQ/TableProducer/tableMaker_withAssoc.cxx @@ -1557,7 +1557,7 @@ struct TableMaker { float deltaR2 = deltaEta * deltaEta + deltaPhi * deltaPhi; auto existing = fTrackEMCalMatchMap.find(match.trackId()); if (existing == fTrackEMCalMatchMap.end() || deltaR2 < existing->second.deltaR2) { - fTrackEMCalMatchMap[match.trackId()] = EMCalMatch{static_cast(outTables.emcal.lastIndex()), deltaEta, deltaPhi, deltaR2}; + fTrackEMCalMatchMap[match.trackId()] = EMCalMatch{.clusterIdx = static_cast(outTables.emcal.lastIndex()), .deltaEta = deltaEta, .deltaPhi = deltaPhi, .deltaR2 = deltaR2}; } } } // end loop over clusters From 3a3348d754862a4573c677bb0979ac2adaa0581e Mon Sep 17 00:00:00 2001 From: Maurice Coquet Date: Sun, 13 Sep 2026 11:57:57 +0200 Subject: [PATCH 9/9] dummy commit --- PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx b/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx index f369795251c..614c1ead2eb 100644 --- a/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx +++ b/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx @@ -577,10 +577,10 @@ struct TableMakerMC { /*if ((std::abs(mctrack.pdgCode())>400 && std::abs(mctrack.pdgCode())<599) || (std::abs(mctrack.pdgCode())>4000 && std::abs(mctrack.pdgCode())<5999) || (mcflags > 0)) { - LOG(info) << ">>>>>>>>>>>>>>>>>>>>>>> track idx / pdg / process / status code / HEPMC status / primary : " + LOG(info) << ">>>>>>>>>>>>>>>>>>>>>> track idx / pdg / process / status code / HEPMC status / primary : " << mctrack.globalIndex() << " / " << mctrack.pdgCode() << " / " << mctrack.getProcess() << " / " << mctrack.getGenStatusCode() << " / " << mctrack.getHepMCStatusCode() << " / " << mctrack.isPhysicalPrimary(); - LOG(info) << ">>>>>>>>>>>>>>>>>>>>>>> track bitmap: "; + LOG(info) << ">>>>>>>>>>>>>>>>>>>>>> track bitmap: "; PrintBitMap(mcflags, 16); if (mctrack.has_mothers()) { for (const auto& m : mctrack.mothersIds()) {