From a5a2ba73e6e67ee470373d5336bbc341f7dea74c Mon Sep 17 00:00:00 2001 From: Maurice Coquet Date: Fri, 2 Jan 2026 10:24:26 +0100 Subject: [PATCH 1/2] Loading z shift of fwdtracks from ccdb --- PWGDQ/Core/VarManager.cxx | 1 + PWGDQ/Core/VarManager.h | 57 +- .../TableProducer/tableMakerMC_withAssoc.cxx | 5 + PWGDQ/TableProducer/tableMaker_withAssoc.cxx | 5 + PWGDQ/Tasks/mftMchMatcher.cxx | 536 ++++++++++++++++++ 5 files changed, 589 insertions(+), 15 deletions(-) create mode 100644 PWGDQ/Tasks/mftMchMatcher.cxx diff --git a/PWGDQ/Core/VarManager.cxx b/PWGDQ/Core/VarManager.cxx index 3d5475649a9..34ae57c8475 100644 --- a/PWGDQ/Core/VarManager.cxx +++ b/PWGDQ/Core/VarManager.cxx @@ -30,6 +30,7 @@ bool VarManager::fgUsedVars[VarManager::kNVars] = {false}; bool VarManager::fgUsedKF = false; float VarManager::fgMagField = 0.5; float VarManager::fgzMatching = -77.5; +float VarManager::fgzShiftFwd = 0.0; float VarManager::fgValues[VarManager::kNVars] = {0.0f}; float VarManager::fgTPCInterSectorBoundary = 1.0; // cm int VarManager::fgITSROFbias = 0; diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index da047fbaed4..d79a19d3b36 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -1051,6 +1051,12 @@ class VarManager : public TObject return fgzMatching; } + // Set z shift for forward tracks + static void SetZShift(float z) + { + fgzShiftFwd = z; + } + // Setup the 2 prong KFParticle static void SetupTwoProngKFParticle(float magField) { @@ -1332,6 +1338,7 @@ class VarManager : public TObject static float fgMagField; static float fgzMatching; + static float fgzShiftFwd; static float fgCenterOfMassEnergy; // collision energy static float fgMassofCollidingParticle; // mass of the colliding particle static float fgTPCInterSectorBoundary; // TPC inter-sector border size at the TPC outer radius, in cm @@ -1473,7 +1480,7 @@ KFPVertex VarManager::createKFPVertexFromCollision(const T& collision) template o2::dataformats::GlobalFwdTrack VarManager::PropagateMuon(const T& muon, const C& collision, const int endPoint) { - o2::track::TrackParCovFwd fwdtrack = FwdToTrackPar(muon, muon); + o2::track::TrackParCovFwd fwdtrack = o2::aod::fwdtrackutils::getTrackParCovFwdShift(muon, fgzShiftFwd, muon); o2::dataformats::GlobalFwdTrack propmuon; if (static_cast(muon.trackType()) > 2) { o2::dataformats::GlobalFwdTrack track; @@ -1501,12 +1508,8 @@ o2::dataformats::GlobalFwdTrack VarManager::PropagateMuon(const T& muon, const C propmuon.setCovariances(proptrack.getCovariances()); } else if (static_cast(muon.trackType()) < 2) { - double centerMFT[3] = {0, 0, -61.4}; - o2::field::MagneticField* field = static_cast(TGeoGlobalMagField::Instance()->GetField()); - auto Bz = field->getBz(centerMFT); // Get field at centre of MFT - auto geoMan = o2::base::GeometryManager::meanMaterialBudget(muon.x(), muon.y(), muon.z(), collision.posX(), collision.posY(), collision.posZ()); - auto x2x0 = static_cast(geoMan.meanX2X0); - fwdtrack.propagateToVtxhelixWithMCS(collision.posZ(), {collision.posX(), collision.posY()}, {collision.covXX(), collision.covYY()}, Bz, x2x0); + 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()); @@ -1609,14 +1612,15 @@ void VarManager::FillGlobalMuonRefit(T1 const& muontrack, T2 const& mfttrack, co double py = propmuon.getP() * sin(M_PI / 2 - atan(mfttrack.tgl())) * sin(mfttrack.phi()); double pz = propmuon.getP() * cos(M_PI / 2 - atan(mfttrack.tgl())); double pt = std::sqrt(std::pow(px, 2) + std::pow(py, 2)); - values[kX] = mfttrack.x(); - values[kY] = mfttrack.y(); - values[kZ] = mfttrack.z(); - values[kTgl] = mfttrack.tgl(); + auto mftprop = o2::aod::fwdtrackutils::getTrackParCovFwdShift(mfttrack, fgzShiftFwd); + values[kX] = mftprop.getX(); + values[kY] = mftprop.getY(); + values[kZ] = mftprop.getZ(); + values[kTgl] = mftprop.getTgl(); values[kPt] = pt; values[kPz] = pz; - values[kEta] = mfttrack.eta(); - values[kPhi] = mfttrack.phi(); + values[kEta] = mftprop.getEta(); + values[kPhi] = mftprop.getPhi(); } } @@ -1629,7 +1633,7 @@ void VarManager::FillGlobalMuonRefitCov(T1 const& muontrack, T2 const& mfttrack, if constexpr ((MuonfillMap & MuonCov) > 0) { if constexpr ((MFTfillMap & MFTCov) > 0) { o2::dataformats::GlobalFwdTrack propmuon = PropagateMuon(muontrack, collision); - o2::track::TrackParCovFwd mft = FwdToTrackPar(mfttrack, mftcov); + auto mft = o2::aod::fwdtrackutils::getTrackParCovFwdShift(mfttrack, fgzShiftFwd, mftcov); o2::dataformats::GlobalFwdTrack globalRefit = o2::aod::fwdtrackutils::refitGlobalMuonCov(propmuon, mft); values[kX] = globalRefit.getX(); @@ -2694,7 +2698,7 @@ void VarManager::FillTrack(T const& track, float* values) values[kMuonTimeRes] = track.trackTimeRes(); } // Quantities based on the muon covariance table - if constexpr ((fillMap & ReducedMuonCov) > 0 || (fillMap & MuonCov) > 0 || (fillMap & MuonCovRealign) > 0) { + if constexpr ((fillMap & ReducedMuonCov) > 0) { values[kX] = track.x(); values[kY] = track.y(); values[kZ] = track.z(); @@ -2715,6 +2719,29 @@ void VarManager::FillTrack(T const& track, float* values) values[kMuonC1Pt2Tgl] = track.c1PtTgl(); values[kMuonC1Pt21Pt2] = track.c1Pt21Pt2(); } + if constexpr ((fillMap & MuonCov) > 0 || (fillMap & MuonCovRealign) > 0) { + auto muonTrack = o2::aod::fwdtrackutils::getTrackParCovFwdShift(track, fgzShiftFwd, track); + auto muonCov = muonTrack.getCovariances(); + values[kX] = muonTrack.getX(); + values[kY] = muonTrack.getY(); + values[kZ] = muonTrack.getZ(); + values[kTgl] = muonTrack.getTgl(); + values[kMuonCXX] = muonCov(0, 0); + values[kMuonCXY] = muonCov(0, 1); + values[kMuonCYY] = muonCov(1, 1); + values[kMuonCPhiX] = muonCov(2, 0); + values[kMuonCPhiY] = muonCov(2, 1); + values[kMuonCPhiPhi] = muonCov(2, 2); + values[kMuonCTglX] = muonCov(3, 0); + values[kMuonCTglY] = muonCov(3, 1); + values[kMuonCTglPhi] = muonCov(3, 2); + values[kMuonCTglTgl] = muonCov(3, 3); + values[kMuonC1Pt2X] = muonCov(4, 0); + values[kMuonC1Pt2Y] = muonCov(4, 1); + values[kMuonC1Pt2Phi] = muonCov(4, 2); + values[kMuonC1Pt2Tgl] = muonCov(4, 3); + values[kMuonC1Pt21Pt2] = muonCov(4, 4); + } // Quantities based on the pair table(s) if constexpr ((fillMap & Pair) > 0) { diff --git a/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx b/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx index a0f39c9e34b..2b99db8386f 100644 --- a/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx +++ b/PWGDQ/TableProducer/tableMakerMC_withAssoc.cxx @@ -207,6 +207,7 @@ 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 fo"}; 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; @@ -1179,10 +1180,14 @@ struct TableMakerMC { } } else { fGrpMag = fCCDB->getForTimeStamp(fConfigCCDB.fGrpMagPath, bcs.begin().timestamp()); + auto* fZShift = fCCDB->getForTimeStamp>(fConfigCCDB.fZShiftPath, bcs.begin().timestamp()); if (fGrpMag != nullptr) { o2::base::Propagator::initFieldFromGRP(fGrpMag); VarManager::SetMagneticField(fGrpMag->getNominalL3Field()); } + if (fZShift != nullptr && !fZShift->empty()) { + VarManager::SetZShift((*fZShift)[0]); + } if (fConfigVariousOptions.fPropMuon) { VarManager::SetupMuonMagField(); } diff --git a/PWGDQ/TableProducer/tableMaker_withAssoc.cxx b/PWGDQ/TableProducer/tableMaker_withAssoc.cxx index 1637cbb1d87..019c85b821b 100644 --- a/PWGDQ/TableProducer/tableMaker_withAssoc.cxx +++ b/PWGDQ/TableProducer/tableMaker_withAssoc.cxx @@ -241,6 +241,7 @@ struct TableMaker { Configurable fConfigNoLaterThan{"ccdb-no-later-than", std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(), "latest acceptable timestamp of creation for the object"}; Configurable fConfigGeoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"}; Configurable fConfigGrpMagPath{"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 fo"}; Configurable fConfigGrpMagPathRun2{"grpmagPathRun2", "GLO/GRP/GRP", "CCDB path of the GRPObject (Usage for Run 2)"}; } fConfigCCDB; @@ -1480,10 +1481,14 @@ struct TableMaker { } } else { fGrpMag = fCCDB->getForTimeStamp(fConfigCCDB.fConfigGrpMagPath, bcs.begin().timestamp()); + auto* fZShift = fCCDB->getForTimeStamp>(fConfigCCDB.fZShiftPath, bcs.begin().timestamp()); if (fGrpMag != nullptr) { o2::base::Propagator::initFieldFromGRP(fGrpMag); VarManager::SetMagneticField(fGrpMag->getNominalL3Field()); } + if (fZShift != nullptr && !fZShift->empty()) { + VarManager::SetZShift((*fZShift)[0]); + } if (fConfigVariousOptions.fPropMuon) { VarManager::SetupMuonMagField(); } diff --git a/PWGDQ/Tasks/mftMchMatcher.cxx b/PWGDQ/Tasks/mftMchMatcher.cxx new file mode 100644 index 00000000000..efa3b297821 --- /dev/null +++ b/PWGDQ/Tasks/mftMchMatcher.cxx @@ -0,0 +1,536 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copdyright 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 mftMchMatcher.cxx +/// \brief MFT-MCH matching tool for data preparation + +#include "Common/DataModel/EventSelection.h" + +#include "CCDB/BasicCCDBManager.h" +#include "DataFormatsParameters/GRPMagField.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisTask.h" +#include "Framework/runDataProcessing.h" +#include "GlobalTracking/MatchGlobalFwd.h" +#include "MFTTracking/Constants.h" +#include "DetectorsBase/Propagator.h" +#include "Field/MagneticField.h" +#include + +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::aod; + +using MyEvents = soa::Join; +using MyMuonsWithCov = soa::Join; +using MyMuonsMC = soa::Join; +using MyMFTs = aod::MFTTracks; +using MyMFTCovariances = aod::MFTTracksCov; +using MyMFTsMC = soa::Join; + +using MyMuon = MyMuonsWithCov::iterator; +using MyMuonMC = MyMuonsMC::iterator; +using MyMFT = MyMFTs::iterator; +using MyMFTCovariance = MyMFTCovariances::iterator; + + +using namespace std; + +using SMatrix55 = ROOT::Math::SMatrix>; +using SMatrix55Std = ROOT::Math::SMatrix; +using SMatrix5 = ROOT::Math::SVector; + +namespace o2::aod +{ +namespace fwdmatchcandidates +{ +DECLARE_SOA_COLUMN(XMCH, xMCH, float); +DECLARE_SOA_COLUMN(YMCH, yMCH, float); +DECLARE_SOA_COLUMN(PhiMCH, phiMCH, float); +DECLARE_SOA_COLUMN(TanlMCH, tanlMCH, float); +DECLARE_SOA_COLUMN(InvQPtMCH, invQPtMCH, float); +DECLARE_SOA_COLUMN(TimeMCH, timeMCH, float); +DECLARE_SOA_COLUMN(TimeResMCH, timeResMCH, float); +DECLARE_SOA_COLUMN(Chi2MCH, chi2MCH, float); +DECLARE_SOA_COLUMN(Rabs, rabs, float); +DECLARE_SOA_COLUMN(PDCA, pDCA, float); +DECLARE_SOA_COLUMN(McMaskMCH, mcMaskMCH, int); + +DECLARE_SOA_COLUMN(SigmaXMCH, sigmaXMCH, float); +DECLARE_SOA_COLUMN(SigmaYMCH, sigmaYMCH, float); +DECLARE_SOA_COLUMN(SigmaPhiMCH, sigmaPhiMCH, float); +DECLARE_SOA_COLUMN(SigmaTglMCH, sigmaTglMCH, float); +DECLARE_SOA_COLUMN(Sigma1PtMCH, sigma1PtMCH, float); +DECLARE_SOA_COLUMN(RhoXYMCH, rhoXYMCH, float); +DECLARE_SOA_COLUMN(RhoPhiXMCH, rhoPhiXMCH, float); +DECLARE_SOA_COLUMN(RhoPhiYMCH, rhoPhiYMCH, float); +DECLARE_SOA_COLUMN(RhoTglXMCH, rhoTglXMCH, float); +DECLARE_SOA_COLUMN(RhoTglYMCH, rhoTglYMCH, float); +DECLARE_SOA_COLUMN(RhoTglPhiMCH, rhoTglPhiMCH, float); +DECLARE_SOA_COLUMN(Rho1PtXMCH, rho1PtXMCH, float); +DECLARE_SOA_COLUMN(Rho1PtYMCH, rho1PtYMCH, float); +DECLARE_SOA_COLUMN(Rho1PtPhiMCH, rho1PtPhiMCH, float); +DECLARE_SOA_COLUMN(Rho1PtTglMCH, rho1PtTglMCH, float); + +DECLARE_SOA_COLUMN(XMFT, xMFT, float); +DECLARE_SOA_COLUMN(YMFT, yMFT, float); +DECLARE_SOA_COLUMN(PhiMFT, phiMFT, float); +DECLARE_SOA_COLUMN(TanlMFT, tanlMFT, float); +DECLARE_SOA_COLUMN(InvQPtMFT, invQPtMFT, float); +DECLARE_SOA_COLUMN(TimeMFT, timeMFT, float); +DECLARE_SOA_COLUMN(TimeResMFT, timeResMFT, float); +DECLARE_SOA_COLUMN(Chi2MFT, chi2MFT, float); +DECLARE_SOA_COLUMN(McMaskMFT, mcMaskMFT, int); + +DECLARE_SOA_COLUMN(SigmaXMFT, sigmaXMFT, float); +DECLARE_SOA_COLUMN(SigmaYMFT, sigmaYMFT, float); +DECLARE_SOA_COLUMN(SigmaPhiMFT, sigmaPhiMFT, float); +DECLARE_SOA_COLUMN(SigmaTglMFT, sigmaTglMFT, float); +DECLARE_SOA_COLUMN(Sigma1PtMFT, sigma1PtMFT, float); +DECLARE_SOA_COLUMN(RhoXYMFT, rhoXYMFT, float); +DECLARE_SOA_COLUMN(RhoPhiYMFT, rhoPhiYMFT, float); +DECLARE_SOA_COLUMN(RhoPhiXMFT, rhoPhiXMFT, float); +DECLARE_SOA_COLUMN(RhoTglXMFT, rhoTglXMFT, float); +DECLARE_SOA_COLUMN(RhoTglYMFT, rhoTglYMFT, float); +DECLARE_SOA_COLUMN(RhoTglPhiMFT, rhoTglPhiMFT, float); +DECLARE_SOA_COLUMN(Rho1PtXMFT, rho1PtXMFT, float); +DECLARE_SOA_COLUMN(Rho1PtYMFT, rho1PtYMFT, float); +DECLARE_SOA_COLUMN(Rho1PtPhiMFT, rho1PtPhiMFT, float); +DECLARE_SOA_COLUMN(Rho1PtTglMFT, rho1PtTglMFT, float); + +DECLARE_SOA_COLUMN(Chi2Glob, chi2Glob, float); +DECLARE_SOA_COLUMN(Chi2Match, chi2Match, float); +DECLARE_SOA_COLUMN(IsAmbig, isAmbig, bool); +DECLARE_SOA_COLUMN(MFTMult, mftMult, int); +DECLARE_SOA_COLUMN(DCAX, dcaX, float); +DECLARE_SOA_COLUMN(DCAY, dcaY, float); +DECLARE_SOA_COLUMN(McMaskGlob, mcMaskGlob, int); +} + +DECLARE_SOA_TABLE(FwdMatchMLCandidates, "AOD", "FWDMLCAND", + fwdmatchcandidates::XMCH, + fwdmatchcandidates::YMCH, + fwdmatchcandidates::PhiMCH, + fwdmatchcandidates::TanlMCH, + fwdmatchcandidates::InvQPtMCH, + fwdmatchcandidates::TimeMCH, + fwdmatchcandidates::TimeResMCH, + fwdmatchcandidates::Chi2MCH, + fwdmatchcandidates::PDCA, + fwdmatchcandidates::Rabs, + fwdmatchcandidates::CXXMCH, + fwdmatchcandidates::CYYMCH, + fwdmatchcandidates::CPhiPhiMCH, + fwdmatchcandidates::CTglTglMCH, + fwdmatchcandidates::C1Pt1PtMCH, + fwdmatchcandidates::CXYMCH, + fwdmatchcandidates::CPhiYMCH, + fwdmatchcandidates::CPhiXMCH, + fwdmatchcandidates::CTglXMCH, + fwdmatchcandidates::CTglYMCH, + fwdmatchcandidates::CTglPhiMCH, + fwdmatchcandidates::C1PtXMCH, + fwdmatchcandidates::C1PtYMCH, + fwdmatchcandidates::C1PtPhiMCH, + fwdmatchcandidates::C1PtTglMCH, + fwdmatchcandidates::XMFT, + fwdmatchcandidates::YMFT, + fwdmatchcandidates::PhiMFT, + fwdmatchcandidates::TanlMFT, + fwdmatchcandidates::InvQPtMFT, + fwdmatchcandidates::TimeMFT, + fwdmatchcandidates::TimeResMFT, + fwdmatchcandidates::Chi2MFT, + fwdmatchcandidates::CXXMFT, + fwdmatchcandidates::CYYMFT, + fwdmatchcandidates::CPhiPhiMFT, + fwdmatchcandidates::CTglTglMFT, + fwdmatchcandidates::C1Pt1PtMFT, + fwdmatchcandidates::CXYMFT, + fwdmatchcandidates::CPhiYMFT, + fwdmatchcandidates::CPhiXMFT, + fwdmatchcandidates::CTglXMFT, + fwdmatchcandidates::CTglYMFT, + fwdmatchcandidates::CTglPhiMFT, + fwdmatchcandidates::C1PtXMFT, + fwdmatchcandidates::C1PtYMFT, + fwdmatchcandidates::C1PtPhiMFT, + fwdmatchcandidates::C1PtTglMFT, + fwdmatchcandidates::Chi2Glob, + fwdmatchcandidates::Chi2Match, + fwdmatchcandidates::DCAX, + fwdmatchcandidates::DCAY, + fwdmatchcandidates::IsAmbig, + fwdmatchcandidates::MFTMult, + fwdmatchcandidates::McMaskMCH, + fwdmatchcandidates::McMaskMFT, + fwdmatchcandidates::McMaskGlob); +} + +struct mftMchMatcher { + Produces fwdMatchMLCandidates; + //// Variables for selecting muon tracks + Configurable fPMchLow{"cfgPMchLow", 0.0f, ""}; + Configurable fPtMchLow{"cfgPtMchLow", 0.7f, ""}; + Configurable fEtaMchLow{"cfgEtaMchLow", -4.0f, ""}; + Configurable fEtaMchUp{"cfgEtaMchUp", -2.5f, ""}; + Configurable fRabsLow{"cfgRabsLow", 17.6f, ""}; + Configurable fRabsUp{"cfgRabsUp", 89.5f, ""}; + Configurable fSigmaPdcaUp{"cfgPdcaUp", 6.f, ""}; + Configurable fTrackChi2MchUp{"cfgTrackChi2MchUp", 5.f, ""}; + Configurable fMatchingChi2MchMidUp{"cfgMatchingChi2MchMidUp", 999.f, ""}; + + //// Variables for selecting mft tracks + Configurable fEtaMftLow{"cfgEtaMftlow", -3.6f, ""}; + Configurable fEtaMftUp{"cfgEtaMftup", -2.5f, ""}; + Configurable fTrackChi2MFTUp{"cfgTrackChi2MchUp", 10.f, ""}; + Configurable fPtMFTLow{"cfgPtMchLow", 0.1f, ""}; + + //// Variables for matching configuration + Configurable fMatchingPlaneZ{"cfgMatchingPlaneZ", -77.5f, ""}; + Configurable fMaxCandidates{"cfgMaxCandidates", 0, ""}; + + //// Variables for ccdb + Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable grpPath{"grpPath", "GLO/GRP/GRP", "Path of the grp file"}; + Configurable grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; + Configurable geoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"}; + + double mBzAtMftCenter{0}; + o2::globaltracking::MatchGlobalFwd mExtrap; + + int mRunNumber{0}; // needed to detect if the run changed and trigger update of magnetic field + Service ccdbManager; + o2::ccdb::CcdbApi fCCDBApi; + + std::unordered_map mftCovIndexes; + + HistogramRegistry registry{"registry", {}}; + + 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; + + double thetaAbs = TMath::ATan(mchTrack.rAtAbsorberEnd() / 505.) * TMath::RadToDeg(); + + // propagate muon track to vertex + auto mchTrackAtVertex = FwdtoMCH(FwdToTrackPar(mchTrack, mchTrack)); + o2::mch::TrackExtrap::extrapToVertex(mchTrackAtVertex, collision.posX(), collision.posY(), collision.posZ(), collision.covXX(), collision.covYY()); + + // double pUncorr = mchTrack.p(); + double p = mchTrackAtVertex.p(); + + double pDCA = mchTrack.pDca(); + double sigmaPDCA = (thetaAbs < 3) ? sigmaPDCA23 : sigmaPDCA310; + double nrp = nSigmaPDCA * relPRes * p; + double pResEffect = sigmaPDCA / (1. - nrp / (1. + nrp)); + double slopeResEffect = 535. * slopeRes * p; + double sigmaPDCAWithRes = TMath::Sqrt(pResEffect * pResEffect + slopeResEffect * slopeResEffect); + if (pDCA > nSigmaPDCA * sigmaPDCAWithRes) { + return false; + } + + return true; + } + + 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 + if (!pDCACut(mchTrack, collision, nSigmaPdcaCut)) { + return false; + } + + return true; + } + + template + bool IsGoodMFT(const T& mftTrack, + double chi2Cut, + double pTCut, + std::array etaCut) + { + // chi2 cut + if (mftTrack.chi2() > chi2Cut) + return false; + + // transverse momentum cut + if (mftTrack.pt() < pTCut) { + return false; // skip low-momentum tracks + } + + // Eta cut + double eta = mftTrack.eta(); + if ((eta < etaCut[0] || eta > etaCut[1])) { + return false; + } + + return true; + } + + + 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(grpmagPath, metadata, ts); + o2::base::Propagator::initFieldFromGRP(grpmag); + LOGF(info, "Set field for muons"); + o2::mch::TrackExtrap::setField(); + if (!o2::base::GeometryManager::isGeometryLoaded()) { + ccdbManager->get(geoPath); + } + o2::mch::TrackExtrap::setField(); + auto* fieldB = static_cast(TGeoGlobalMagField::Instance()->GetField()); + if (fieldB) { + double centerMFT[3] = {0, 0, -61.4}; // Field at center of MFT + mBzAtMftCenter = fieldB->getBz(centerMFT); + // std::cout << "fieldB: " << (void*)fieldB << std::endl; + } + } + + void init(o2::framework::InitContext&) + { + // Load geometry + ccdbManager->setURL(ccdburl); + ccdbManager->setCaching(true); + ccdbManager->setLocalObjectValidityChecking(); + fCCDBApi.init(ccdburl); + mRunNumber = 0; + + if (!o2::base::GeometryManager::isGeometryLoaded()) { + LOGF(info, "Load geometry from CCDB"); + ccdbManager->get(geoPath); + } + } + + template + void skimBestMuonMatches(TMuons const& muons) + { + std::unordered_map> mCandidates; + for (const auto& muon : muons) { + if (static_cast(muon.trackType()) < 2) { + auto muonID = muon.matchMCHTrackId(); + auto chi2 = muon.chi2MatchMCHMFT(); + if (fConfigVariousOptions.fUseML.value) { + std::vector output; + std::vector inputML = matchingMlResponse.getInputFeaturesTest(muon); + matchingMlResponse.isSelectedMl(inputML, 0, output); + chi2 = output[0]; + } + if (mCandidates.find(muonID) == mCandidates.end()) { + mCandidates[muonID] = {chi2, muon.globalIndex()}; + } else { + if (chi2 < mCandidates[muonID].first) { + mCandidates[muonID] = {chi2, muon.globalIndex()}; + } + } + } + } + for (auto& pairCand : mCandidates) { + fBestMatch[pairCand.second.second] = true; + } + } + + void processMC(MyEvents const& collisions, + aod::BCsWithTimestamps const& bcs, + MyMuonsMC const& muonTracks, + MyMFTsMC const& mftTracks, + MyMFTCovariances const& mftCovs, + aod::McParticles const& /*mcParticles*/) + { + auto bc = bcs.begin(); + initCCDB(bc); + + mftCovIndexes.clear(); + for (auto& mftTrackCov : mftCovs) { + mftCovIndexes[mftTrackCov.matchMFTTrackId()] = mftTrackCov.globalIndex(); + } + + fwdMatchMLCandidates.reserve(muonTracks.size()); + + for (auto muon : muonTracks) { + // only consider global MFT-MCH-MID matches + if (static_cast(muon.trackType()) >= 2) { + continue; + } + + if (!muon.has_collision()) { + continue; + } + + if (fKeepBestMatch) { + if (fBestMatch.find(muon.globalIndex()) == fBestMatch.end()) { + continue; + } + } + + const auto& collision = collisions.rawIteratorAt(muon.collisionId()); + + auto muontrack = muon.template matchMCHTrack_as(); + auto mfttrack = muon.template matchMFTTrack_as(); + auto const& mfttrackcov = mfCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]); + + auto muonTime = 0.f; + auto mftTime = 0.f; + + o2::track::TrackParCovFwd mftprop = VarManager::FwdToTrackPar(mfttrack, mfttrackcov); + o2::dataformats::GlobalFwdTrack muonprop = VarManager::FwdToTrackPar(muontrack, muontrack); + if (fConfigVariousOptions.fzMatching.value < 0.) { + mftprop = VarManager::PropagateFwd(mfttrack, mfttrackcov, fConfigVariousOptions.fzMatching.value); + muonprop = VarManager::PropagateMuon(muontrack, collision, VarManager::kToMatching); + } + auto muonpropCov = muonprop.getCovariances(); + auto mftpropCov = mftprop.getCovariances(); + + if (!IsGoodMuon(muontrack, collision, fTrackChi2MchUp, fPMchLow, fPtMchLow, {fEtaMchLow, cfgEtaMchUp}, {cfgRabsLow, cfgRabsUp}, fSigmaPdcaUp)){ + continue; + } + + + if (!IsGoodMFT(mfttrack, fTrackChi2MFTUp, fPtMFTLow, {fEtaMFTLow, cfgEtaMFTUp})){ + continue; + } + + bool IsAmbig = (muon.compatibleCollIds().size() != 1); + int MFTMult = collision.mftNtracks(); + + fwdMatchMLCandidates( + muonprop.getX(), + muonprop.getY(), + muonprop.getPhi(), + muonprop.getTgl(), + muonprop.getSigned1Pt(), + muonTime, + muontrack.timeRes(), + muontrack.chi2(), + muontrack.pDca(), + muontrack.rAbs(), + muonpropCov(0, 0), + muonpropCov(1, 1), + muonpropCov(2, 2), + muonpropCov(3, 3), + muonpropCov(4, 4), + muonpropCov(1, 0), + muonpropCov(2, 0), + muonpropCov(2, 1), + muonpropCov(3, 0), + muonpropCov(3, 1), + muonpropCov(3, 2), + muonpropCov(4, 0), + muonpropCov(4, 1), + muonpropCov(4, 2), + muonpropCov(4, 3), + + mftprop.getX(), + mftprop.getY(), + mftprop.getPhi(), + mftprop.getTgl(), + mftprop.getSigned1Pt(), + mftTime, + mfttrack.timeRes(), + mfttrack.chi2(), + mftpropCov(0, 0), + mftpropCov(1, 1), + mftpropCov(2, 2), + mftpropCov(3, 3), + mftpropCov(4, 4), + mftpropCov(1, 0), + mftpropCov(2, 0), + mftpropCov(2, 1), + mftpropCov(3, 0), + mftpropCov(3, 1), + mftpropCov(3, 2), + mftpropCov(4, 0), + mftpropCov(4, 1), + mftpropCov(4, 2), + mftpropCov(4, 3), + + muon.chi2(), + muon.chi2MatchMCHMFT(), + muon.fwddcaX(), + muon.fwddcaY(), + IsAmbig, + MFTMult, + + muontrack.mcMask(), + mfttrack.mcMask(), + muon.mcMask(), + ); + + + + } + + } + + PROCESS_SWITCH(mftMchMatcher, processMC, "process_MC", true); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc)}; +}; + From 95a6ae2bcf7cb3c3707f1ca1ff16457016b4d69d Mon Sep 17 00:00:00 2001 From: Maurice Coquet Date: Fri, 2 Jan 2026 10:27:16 +0100 Subject: [PATCH 2/2] clang format --- PWGDQ/Core/VarManager.h | 2 +- PWGDQ/TableProducer/tableMaker_withAssoc.cxx | 2 +- PWGDQ/Tasks/mftMchMatcher.cxx | 536 ------------------- 3 files changed, 2 insertions(+), 538 deletions(-) delete mode 100644 PWGDQ/Tasks/mftMchMatcher.cxx diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index d79a19d3b36..afcbb382044 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -1480,7 +1480,7 @@ KFPVertex VarManager::createKFPVertexFromCollision(const T& collision) template o2::dataformats::GlobalFwdTrack VarManager::PropagateMuon(const T& muon, const C& collision, const int endPoint) { - o2::track::TrackParCovFwd fwdtrack = o2::aod::fwdtrackutils::getTrackParCovFwdShift(muon, fgzShiftFwd, muon); + o2::track::TrackParCovFwd fwdtrack = o2::aod::fwdtrackutils::getTrackParCovFwdShift(muon, fgzShiftFwd, muon); o2::dataformats::GlobalFwdTrack propmuon; if (static_cast(muon.trackType()) > 2) { o2::dataformats::GlobalFwdTrack track; diff --git a/PWGDQ/TableProducer/tableMaker_withAssoc.cxx b/PWGDQ/TableProducer/tableMaker_withAssoc.cxx index 019c85b821b..33b2ddfea1f 100644 --- a/PWGDQ/TableProducer/tableMaker_withAssoc.cxx +++ b/PWGDQ/TableProducer/tableMaker_withAssoc.cxx @@ -1481,7 +1481,7 @@ struct TableMaker { } } else { fGrpMag = fCCDB->getForTimeStamp(fConfigCCDB.fConfigGrpMagPath, bcs.begin().timestamp()); - auto* fZShift = fCCDB->getForTimeStamp>(fConfigCCDB.fZShiftPath, bcs.begin().timestamp()); + auto* fZShift = fCCDB->getForTimeStamp>(fConfigCCDB.fZShiftPath, bcs.begin().timestamp()); if (fGrpMag != nullptr) { o2::base::Propagator::initFieldFromGRP(fGrpMag); VarManager::SetMagneticField(fGrpMag->getNominalL3Field()); diff --git a/PWGDQ/Tasks/mftMchMatcher.cxx b/PWGDQ/Tasks/mftMchMatcher.cxx deleted file mode 100644 index efa3b297821..00000000000 --- a/PWGDQ/Tasks/mftMchMatcher.cxx +++ /dev/null @@ -1,536 +0,0 @@ -// Copyright 2019-2020 CERN and copyright holders of ALICE O2. -// See https://alice-o2.web.cern.ch/copdyright 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 mftMchMatcher.cxx -/// \brief MFT-MCH matching tool for data preparation - -#include "Common/DataModel/EventSelection.h" - -#include "CCDB/BasicCCDBManager.h" -#include "DataFormatsParameters/GRPMagField.h" -#include "Framework/ASoAHelpers.h" -#include "Framework/AnalysisTask.h" -#include "Framework/runDataProcessing.h" -#include "GlobalTracking/MatchGlobalFwd.h" -#include "MFTTracking/Constants.h" -#include "DetectorsBase/Propagator.h" -#include "Field/MagneticField.h" -#include - -#include - -#include -#include -#include -#include -#include -#include -#include -#include - -using namespace o2; -using namespace o2::framework; -using namespace o2::aod; - -using MyEvents = soa::Join; -using MyMuonsWithCov = soa::Join; -using MyMuonsMC = soa::Join; -using MyMFTs = aod::MFTTracks; -using MyMFTCovariances = aod::MFTTracksCov; -using MyMFTsMC = soa::Join; - -using MyMuon = MyMuonsWithCov::iterator; -using MyMuonMC = MyMuonsMC::iterator; -using MyMFT = MyMFTs::iterator; -using MyMFTCovariance = MyMFTCovariances::iterator; - - -using namespace std; - -using SMatrix55 = ROOT::Math::SMatrix>; -using SMatrix55Std = ROOT::Math::SMatrix; -using SMatrix5 = ROOT::Math::SVector; - -namespace o2::aod -{ -namespace fwdmatchcandidates -{ -DECLARE_SOA_COLUMN(XMCH, xMCH, float); -DECLARE_SOA_COLUMN(YMCH, yMCH, float); -DECLARE_SOA_COLUMN(PhiMCH, phiMCH, float); -DECLARE_SOA_COLUMN(TanlMCH, tanlMCH, float); -DECLARE_SOA_COLUMN(InvQPtMCH, invQPtMCH, float); -DECLARE_SOA_COLUMN(TimeMCH, timeMCH, float); -DECLARE_SOA_COLUMN(TimeResMCH, timeResMCH, float); -DECLARE_SOA_COLUMN(Chi2MCH, chi2MCH, float); -DECLARE_SOA_COLUMN(Rabs, rabs, float); -DECLARE_SOA_COLUMN(PDCA, pDCA, float); -DECLARE_SOA_COLUMN(McMaskMCH, mcMaskMCH, int); - -DECLARE_SOA_COLUMN(SigmaXMCH, sigmaXMCH, float); -DECLARE_SOA_COLUMN(SigmaYMCH, sigmaYMCH, float); -DECLARE_SOA_COLUMN(SigmaPhiMCH, sigmaPhiMCH, float); -DECLARE_SOA_COLUMN(SigmaTglMCH, sigmaTglMCH, float); -DECLARE_SOA_COLUMN(Sigma1PtMCH, sigma1PtMCH, float); -DECLARE_SOA_COLUMN(RhoXYMCH, rhoXYMCH, float); -DECLARE_SOA_COLUMN(RhoPhiXMCH, rhoPhiXMCH, float); -DECLARE_SOA_COLUMN(RhoPhiYMCH, rhoPhiYMCH, float); -DECLARE_SOA_COLUMN(RhoTglXMCH, rhoTglXMCH, float); -DECLARE_SOA_COLUMN(RhoTglYMCH, rhoTglYMCH, float); -DECLARE_SOA_COLUMN(RhoTglPhiMCH, rhoTglPhiMCH, float); -DECLARE_SOA_COLUMN(Rho1PtXMCH, rho1PtXMCH, float); -DECLARE_SOA_COLUMN(Rho1PtYMCH, rho1PtYMCH, float); -DECLARE_SOA_COLUMN(Rho1PtPhiMCH, rho1PtPhiMCH, float); -DECLARE_SOA_COLUMN(Rho1PtTglMCH, rho1PtTglMCH, float); - -DECLARE_SOA_COLUMN(XMFT, xMFT, float); -DECLARE_SOA_COLUMN(YMFT, yMFT, float); -DECLARE_SOA_COLUMN(PhiMFT, phiMFT, float); -DECLARE_SOA_COLUMN(TanlMFT, tanlMFT, float); -DECLARE_SOA_COLUMN(InvQPtMFT, invQPtMFT, float); -DECLARE_SOA_COLUMN(TimeMFT, timeMFT, float); -DECLARE_SOA_COLUMN(TimeResMFT, timeResMFT, float); -DECLARE_SOA_COLUMN(Chi2MFT, chi2MFT, float); -DECLARE_SOA_COLUMN(McMaskMFT, mcMaskMFT, int); - -DECLARE_SOA_COLUMN(SigmaXMFT, sigmaXMFT, float); -DECLARE_SOA_COLUMN(SigmaYMFT, sigmaYMFT, float); -DECLARE_SOA_COLUMN(SigmaPhiMFT, sigmaPhiMFT, float); -DECLARE_SOA_COLUMN(SigmaTglMFT, sigmaTglMFT, float); -DECLARE_SOA_COLUMN(Sigma1PtMFT, sigma1PtMFT, float); -DECLARE_SOA_COLUMN(RhoXYMFT, rhoXYMFT, float); -DECLARE_SOA_COLUMN(RhoPhiYMFT, rhoPhiYMFT, float); -DECLARE_SOA_COLUMN(RhoPhiXMFT, rhoPhiXMFT, float); -DECLARE_SOA_COLUMN(RhoTglXMFT, rhoTglXMFT, float); -DECLARE_SOA_COLUMN(RhoTglYMFT, rhoTglYMFT, float); -DECLARE_SOA_COLUMN(RhoTglPhiMFT, rhoTglPhiMFT, float); -DECLARE_SOA_COLUMN(Rho1PtXMFT, rho1PtXMFT, float); -DECLARE_SOA_COLUMN(Rho1PtYMFT, rho1PtYMFT, float); -DECLARE_SOA_COLUMN(Rho1PtPhiMFT, rho1PtPhiMFT, float); -DECLARE_SOA_COLUMN(Rho1PtTglMFT, rho1PtTglMFT, float); - -DECLARE_SOA_COLUMN(Chi2Glob, chi2Glob, float); -DECLARE_SOA_COLUMN(Chi2Match, chi2Match, float); -DECLARE_SOA_COLUMN(IsAmbig, isAmbig, bool); -DECLARE_SOA_COLUMN(MFTMult, mftMult, int); -DECLARE_SOA_COLUMN(DCAX, dcaX, float); -DECLARE_SOA_COLUMN(DCAY, dcaY, float); -DECLARE_SOA_COLUMN(McMaskGlob, mcMaskGlob, int); -} - -DECLARE_SOA_TABLE(FwdMatchMLCandidates, "AOD", "FWDMLCAND", - fwdmatchcandidates::XMCH, - fwdmatchcandidates::YMCH, - fwdmatchcandidates::PhiMCH, - fwdmatchcandidates::TanlMCH, - fwdmatchcandidates::InvQPtMCH, - fwdmatchcandidates::TimeMCH, - fwdmatchcandidates::TimeResMCH, - fwdmatchcandidates::Chi2MCH, - fwdmatchcandidates::PDCA, - fwdmatchcandidates::Rabs, - fwdmatchcandidates::CXXMCH, - fwdmatchcandidates::CYYMCH, - fwdmatchcandidates::CPhiPhiMCH, - fwdmatchcandidates::CTglTglMCH, - fwdmatchcandidates::C1Pt1PtMCH, - fwdmatchcandidates::CXYMCH, - fwdmatchcandidates::CPhiYMCH, - fwdmatchcandidates::CPhiXMCH, - fwdmatchcandidates::CTglXMCH, - fwdmatchcandidates::CTglYMCH, - fwdmatchcandidates::CTglPhiMCH, - fwdmatchcandidates::C1PtXMCH, - fwdmatchcandidates::C1PtYMCH, - fwdmatchcandidates::C1PtPhiMCH, - fwdmatchcandidates::C1PtTglMCH, - fwdmatchcandidates::XMFT, - fwdmatchcandidates::YMFT, - fwdmatchcandidates::PhiMFT, - fwdmatchcandidates::TanlMFT, - fwdmatchcandidates::InvQPtMFT, - fwdmatchcandidates::TimeMFT, - fwdmatchcandidates::TimeResMFT, - fwdmatchcandidates::Chi2MFT, - fwdmatchcandidates::CXXMFT, - fwdmatchcandidates::CYYMFT, - fwdmatchcandidates::CPhiPhiMFT, - fwdmatchcandidates::CTglTglMFT, - fwdmatchcandidates::C1Pt1PtMFT, - fwdmatchcandidates::CXYMFT, - fwdmatchcandidates::CPhiYMFT, - fwdmatchcandidates::CPhiXMFT, - fwdmatchcandidates::CTglXMFT, - fwdmatchcandidates::CTglYMFT, - fwdmatchcandidates::CTglPhiMFT, - fwdmatchcandidates::C1PtXMFT, - fwdmatchcandidates::C1PtYMFT, - fwdmatchcandidates::C1PtPhiMFT, - fwdmatchcandidates::C1PtTglMFT, - fwdmatchcandidates::Chi2Glob, - fwdmatchcandidates::Chi2Match, - fwdmatchcandidates::DCAX, - fwdmatchcandidates::DCAY, - fwdmatchcandidates::IsAmbig, - fwdmatchcandidates::MFTMult, - fwdmatchcandidates::McMaskMCH, - fwdmatchcandidates::McMaskMFT, - fwdmatchcandidates::McMaskGlob); -} - -struct mftMchMatcher { - Produces fwdMatchMLCandidates; - //// Variables for selecting muon tracks - Configurable fPMchLow{"cfgPMchLow", 0.0f, ""}; - Configurable fPtMchLow{"cfgPtMchLow", 0.7f, ""}; - Configurable fEtaMchLow{"cfgEtaMchLow", -4.0f, ""}; - Configurable fEtaMchUp{"cfgEtaMchUp", -2.5f, ""}; - Configurable fRabsLow{"cfgRabsLow", 17.6f, ""}; - Configurable fRabsUp{"cfgRabsUp", 89.5f, ""}; - Configurable fSigmaPdcaUp{"cfgPdcaUp", 6.f, ""}; - Configurable fTrackChi2MchUp{"cfgTrackChi2MchUp", 5.f, ""}; - Configurable fMatchingChi2MchMidUp{"cfgMatchingChi2MchMidUp", 999.f, ""}; - - //// Variables for selecting mft tracks - Configurable fEtaMftLow{"cfgEtaMftlow", -3.6f, ""}; - Configurable fEtaMftUp{"cfgEtaMftup", -2.5f, ""}; - Configurable fTrackChi2MFTUp{"cfgTrackChi2MchUp", 10.f, ""}; - Configurable fPtMFTLow{"cfgPtMchLow", 0.1f, ""}; - - //// Variables for matching configuration - Configurable fMatchingPlaneZ{"cfgMatchingPlaneZ", -77.5f, ""}; - Configurable fMaxCandidates{"cfgMaxCandidates", 0, ""}; - - //// Variables for ccdb - Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; - Configurable grpPath{"grpPath", "GLO/GRP/GRP", "Path of the grp file"}; - Configurable grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; - Configurable geoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"}; - - double mBzAtMftCenter{0}; - o2::globaltracking::MatchGlobalFwd mExtrap; - - int mRunNumber{0}; // needed to detect if the run changed and trigger update of magnetic field - Service ccdbManager; - o2::ccdb::CcdbApi fCCDBApi; - - std::unordered_map mftCovIndexes; - - HistogramRegistry registry{"registry", {}}; - - 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; - - double thetaAbs = TMath::ATan(mchTrack.rAtAbsorberEnd() / 505.) * TMath::RadToDeg(); - - // propagate muon track to vertex - auto mchTrackAtVertex = FwdtoMCH(FwdToTrackPar(mchTrack, mchTrack)); - o2::mch::TrackExtrap::extrapToVertex(mchTrackAtVertex, collision.posX(), collision.posY(), collision.posZ(), collision.covXX(), collision.covYY()); - - // double pUncorr = mchTrack.p(); - double p = mchTrackAtVertex.p(); - - double pDCA = mchTrack.pDca(); - double sigmaPDCA = (thetaAbs < 3) ? sigmaPDCA23 : sigmaPDCA310; - double nrp = nSigmaPDCA * relPRes * p; - double pResEffect = sigmaPDCA / (1. - nrp / (1. + nrp)); - double slopeResEffect = 535. * slopeRes * p; - double sigmaPDCAWithRes = TMath::Sqrt(pResEffect * pResEffect + slopeResEffect * slopeResEffect); - if (pDCA > nSigmaPDCA * sigmaPDCAWithRes) { - return false; - } - - return true; - } - - 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 - if (!pDCACut(mchTrack, collision, nSigmaPdcaCut)) { - return false; - } - - return true; - } - - template - bool IsGoodMFT(const T& mftTrack, - double chi2Cut, - double pTCut, - std::array etaCut) - { - // chi2 cut - if (mftTrack.chi2() > chi2Cut) - return false; - - // transverse momentum cut - if (mftTrack.pt() < pTCut) { - return false; // skip low-momentum tracks - } - - // Eta cut - double eta = mftTrack.eta(); - if ((eta < etaCut[0] || eta > etaCut[1])) { - return false; - } - - return true; - } - - - 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(grpmagPath, metadata, ts); - o2::base::Propagator::initFieldFromGRP(grpmag); - LOGF(info, "Set field for muons"); - o2::mch::TrackExtrap::setField(); - if (!o2::base::GeometryManager::isGeometryLoaded()) { - ccdbManager->get(geoPath); - } - o2::mch::TrackExtrap::setField(); - auto* fieldB = static_cast(TGeoGlobalMagField::Instance()->GetField()); - if (fieldB) { - double centerMFT[3] = {0, 0, -61.4}; // Field at center of MFT - mBzAtMftCenter = fieldB->getBz(centerMFT); - // std::cout << "fieldB: " << (void*)fieldB << std::endl; - } - } - - void init(o2::framework::InitContext&) - { - // Load geometry - ccdbManager->setURL(ccdburl); - ccdbManager->setCaching(true); - ccdbManager->setLocalObjectValidityChecking(); - fCCDBApi.init(ccdburl); - mRunNumber = 0; - - if (!o2::base::GeometryManager::isGeometryLoaded()) { - LOGF(info, "Load geometry from CCDB"); - ccdbManager->get(geoPath); - } - } - - template - void skimBestMuonMatches(TMuons const& muons) - { - std::unordered_map> mCandidates; - for (const auto& muon : muons) { - if (static_cast(muon.trackType()) < 2) { - auto muonID = muon.matchMCHTrackId(); - auto chi2 = muon.chi2MatchMCHMFT(); - if (fConfigVariousOptions.fUseML.value) { - std::vector output; - std::vector inputML = matchingMlResponse.getInputFeaturesTest(muon); - matchingMlResponse.isSelectedMl(inputML, 0, output); - chi2 = output[0]; - } - if (mCandidates.find(muonID) == mCandidates.end()) { - mCandidates[muonID] = {chi2, muon.globalIndex()}; - } else { - if (chi2 < mCandidates[muonID].first) { - mCandidates[muonID] = {chi2, muon.globalIndex()}; - } - } - } - } - for (auto& pairCand : mCandidates) { - fBestMatch[pairCand.second.second] = true; - } - } - - void processMC(MyEvents const& collisions, - aod::BCsWithTimestamps const& bcs, - MyMuonsMC const& muonTracks, - MyMFTsMC const& mftTracks, - MyMFTCovariances const& mftCovs, - aod::McParticles const& /*mcParticles*/) - { - auto bc = bcs.begin(); - initCCDB(bc); - - mftCovIndexes.clear(); - for (auto& mftTrackCov : mftCovs) { - mftCovIndexes[mftTrackCov.matchMFTTrackId()] = mftTrackCov.globalIndex(); - } - - fwdMatchMLCandidates.reserve(muonTracks.size()); - - for (auto muon : muonTracks) { - // only consider global MFT-MCH-MID matches - if (static_cast(muon.trackType()) >= 2) { - continue; - } - - if (!muon.has_collision()) { - continue; - } - - if (fKeepBestMatch) { - if (fBestMatch.find(muon.globalIndex()) == fBestMatch.end()) { - continue; - } - } - - const auto& collision = collisions.rawIteratorAt(muon.collisionId()); - - auto muontrack = muon.template matchMCHTrack_as(); - auto mfttrack = muon.template matchMFTTrack_as(); - auto const& mfttrackcov = mfCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]); - - auto muonTime = 0.f; - auto mftTime = 0.f; - - o2::track::TrackParCovFwd mftprop = VarManager::FwdToTrackPar(mfttrack, mfttrackcov); - o2::dataformats::GlobalFwdTrack muonprop = VarManager::FwdToTrackPar(muontrack, muontrack); - if (fConfigVariousOptions.fzMatching.value < 0.) { - mftprop = VarManager::PropagateFwd(mfttrack, mfttrackcov, fConfigVariousOptions.fzMatching.value); - muonprop = VarManager::PropagateMuon(muontrack, collision, VarManager::kToMatching); - } - auto muonpropCov = muonprop.getCovariances(); - auto mftpropCov = mftprop.getCovariances(); - - if (!IsGoodMuon(muontrack, collision, fTrackChi2MchUp, fPMchLow, fPtMchLow, {fEtaMchLow, cfgEtaMchUp}, {cfgRabsLow, cfgRabsUp}, fSigmaPdcaUp)){ - continue; - } - - - if (!IsGoodMFT(mfttrack, fTrackChi2MFTUp, fPtMFTLow, {fEtaMFTLow, cfgEtaMFTUp})){ - continue; - } - - bool IsAmbig = (muon.compatibleCollIds().size() != 1); - int MFTMult = collision.mftNtracks(); - - fwdMatchMLCandidates( - muonprop.getX(), - muonprop.getY(), - muonprop.getPhi(), - muonprop.getTgl(), - muonprop.getSigned1Pt(), - muonTime, - muontrack.timeRes(), - muontrack.chi2(), - muontrack.pDca(), - muontrack.rAbs(), - muonpropCov(0, 0), - muonpropCov(1, 1), - muonpropCov(2, 2), - muonpropCov(3, 3), - muonpropCov(4, 4), - muonpropCov(1, 0), - muonpropCov(2, 0), - muonpropCov(2, 1), - muonpropCov(3, 0), - muonpropCov(3, 1), - muonpropCov(3, 2), - muonpropCov(4, 0), - muonpropCov(4, 1), - muonpropCov(4, 2), - muonpropCov(4, 3), - - mftprop.getX(), - mftprop.getY(), - mftprop.getPhi(), - mftprop.getTgl(), - mftprop.getSigned1Pt(), - mftTime, - mfttrack.timeRes(), - mfttrack.chi2(), - mftpropCov(0, 0), - mftpropCov(1, 1), - mftpropCov(2, 2), - mftpropCov(3, 3), - mftpropCov(4, 4), - mftpropCov(1, 0), - mftpropCov(2, 0), - mftpropCov(2, 1), - mftpropCov(3, 0), - mftpropCov(3, 1), - mftpropCov(3, 2), - mftpropCov(4, 0), - mftpropCov(4, 1), - mftpropCov(4, 2), - mftpropCov(4, 3), - - muon.chi2(), - muon.chi2MatchMCHMFT(), - muon.fwddcaX(), - muon.fwddcaY(), - IsAmbig, - MFTMult, - - muontrack.mcMask(), - mfttrack.mcMask(), - muon.mcMask(), - ); - - - - } - - } - - PROCESS_SWITCH(mftMchMatcher, processMC, "process_MC", true); -}; - -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) -{ - return WorkflowSpec{ - adaptAnalysisTask(cfgc)}; -}; -