From 8c785379a7a70596882ed18cbfc0503db96ecbbc Mon Sep 17 00:00:00 2001 From: Rafael Pezzi Date: Tue, 30 Nov 2021 16:23:29 +0100 Subject: [PATCH 1/2] Improve FwdMatcher configuration --- .../include/GlobalTracking/MatchGlobalFwd.h | 84 ++---- .../GlobalTracking/MatchGlobalFwdParam.h | 3 +- .../GlobalTracking/src/MatchGlobalFwd.cxx | 278 +++++++++--------- 3 files changed, 166 insertions(+), 199 deletions(-) diff --git a/Detectors/GlobalTracking/include/GlobalTracking/MatchGlobalFwd.h b/Detectors/GlobalTracking/include/GlobalTracking/MatchGlobalFwd.h index 8d96c60c28d12..83071eb71ecc5 100644 --- a/Detectors/GlobalTracking/include/GlobalTracking/MatchGlobalFwd.h +++ b/Detectors/GlobalTracking/include/GlobalTracking/MatchGlobalFwd.h @@ -78,6 +78,11 @@ struct TrackLocMCH : public o2::dataformats::GlobalFwdTrack { ClassDefNV(TrackLocMCH, 0); }; +using o2::dataformats::GlobalFwdTrack; +using o2::track::TrackParCovFwd; +typedef std::function MatchingFunc_t; +typedef std::function CutFunc_t; + class MatchGlobalFwd { public: @@ -99,13 +104,13 @@ class MatchGlobalFwd enum MatchingType : uint8_t { ///< MFT-MCH matching modes MATCHINGFUNC, ///< Matching function-based MFT-MCH track matching - MATCHINGMCLABEL, ///< Monte Carlo label-based MFT-MCH track matching - MATCHINGUPSTREAM ///< MFT-MCH track matching loaded from input file + MATCHINGUPSTREAM, ///< MFT-MCH track matching loaded from input file + MATCHINGUNDEFINED }; static constexpr Double_t sLastMFTPlaneZ = o2::mft::constants::mft::LayerZCoordinate()[9]; - MatchGlobalFwd() = default; + MatchGlobalFwd(); ~MatchGlobalFwd() = default; void run(const o2::globaltracking::RecoContainer& inp); @@ -133,6 +138,7 @@ class MatchGlobalFwd private: void updateTimeDependentParams(); + void fillBuiltinFunctions(); bool prepareMCHData(); bool prepareMFTData(); @@ -148,8 +154,9 @@ class MatchGlobalFwd ///< Matches MFT tracks in one MFT ROFrame with all MCH tracks in the overlapping MCH ROFrames template void ROFMatch(int MFTROFId, int firstMCHROFId, int lastMCHROFId); - void fitTracks(); // Fit all matched tracks - void fitGlobalMuonTrack(o2::dataformats::GlobalFwdTrack&); // Kalman filter fit global Forward track by attaching MFT clusters + + void fitTracks(); ///< Fit all matched tracks + void fitGlobalMuonTrack(o2::dataformats::GlobalFwdTrack&); ///< Kalman filter fit global Forward track by attaching MFT clusters bool computeCluster(o2::dataformats::GlobalFwdTrack& track, const MFTCluster& cluster, int& startingLayerID); template @@ -223,77 +230,28 @@ class MatchGlobalFwd return true; } - void configMatching(const std::string& matchingFcn, const std::string& cutFcn); - double matchingEval(const TrackLocMCH& mchTrack, const TrackLocMFT& mftTrack); - bool matchingCut(const TrackLocMCH&, const TrackLocMFT&); - double (MatchGlobalFwd::*mMatchFunc)(const TrackLocMCH& mchTrack, const TrackLocMFT& mftTrack) = &MatchGlobalFwd::noMatchFcn; - void setMatchingFunction(double (MatchGlobalFwd::*func)(const TrackLocMCH&, const TrackLocMFT&)) - { - mMatchFunc = func; - } - bool (MatchGlobalFwd::*mCutFunc)(const TrackLocMCH& mchTrack, const TrackLocMFT& mftTrack); - - void setCutFunction(bool (MatchGlobalFwd::*func)(const TrackLocMCH&, const TrackLocMFT&)) - { - mCutFunc = func; - } - /// Matching methods - /// Position - double matchMFT_MCH_TracksXY(const TrackLocMCH& mchTrack, const TrackLocMFT& mftTrack); - /// Position & Angles - double matchMFT_MCH_TracksXYPhiTanl(const TrackLocMCH& mchTrack, const TrackLocMFT& mftTrack); - /// Position, Angles & Charged Momentum - double matchMFT_MCH_TracksAllParam(const TrackLocMCH& mchTrack, const TrackLocMFT& mftTrack); - /// Hiroshima's Matching - double matchHiroshima(const TrackLocMCH& mchTrack, const TrackLocMFT& mftTrack); - /// MC label matching - double noMatchFcn(const TrackLocMCH& mchTrack, const TrackLocMFT& mftTrack) - { - throw std::runtime_error("noMatchFcn function must not be called!"); - } - - /// Cut functions - bool cutDisabled(const TrackLocMCH& mchTrack, const TrackLocMFT& mftTrack) { return true; }; - - bool matchCut3Sigma(const TrackLocMCH& mchTrack, const TrackLocMFT& mftTrack) - { - auto dx = mchTrack.getX() - mftTrack.getX(); - auto dy = mchTrack.getY() - mftTrack.getY(); - auto dPhi = mchTrack.getPhi() - mftTrack.getPhi(); - auto dTanl = TMath::Abs(mchTrack.getTanl() - mftTrack.getTanl()); - auto dInvQPt = TMath::Abs(mchTrack.getInvQPt() - mftTrack.getInvQPt()); - auto distanceSq = dx * dx + dy * dy; - auto cutDistanceSq = 9 * (mchTrack.getSigma2X() + mchTrack.getSigma2Y()); - auto cutPhiSq = 9 * (mchTrack.getSigma2Phi() + mftTrack.getSigma2Phi()); - auto cutTanlSq = 9 * (mchTrack.getSigma2Tanl() + mftTrack.getSigma2Tanl()); - auto cutInvQPtSq = 9 * (mchTrack.getSigma2InvQPt() + mftTrack.getSigma2InvQPt()); - return (distanceSq < cutDistanceSq) and (dPhi * dPhi < cutPhiSq) and (dTanl * dTanl < cutTanlSq) and (dInvQPt * dInvQPt < cutInvQPtSq); + MatchingFunc_t mMatchFunc = [](const GlobalFwdTrack& mchtrack, const TrackParCovFwd& mfttrack) -> double { + throw std::runtime_error("MatchGlobalFwd: matching function not configured!"); }; - bool matchCut3SigmaXYAngles(const TrackLocMCH& mchTrack, const TrackLocMFT& mftTrack) - { - auto dx = mchTrack.getX() - mftTrack.getX(); - auto dy = mchTrack.getY() - mftTrack.getY(); - auto dPhi = mchTrack.getPhi() - mftTrack.getPhi(); - auto dTanl = TMath::Abs(mchTrack.getTanl() - mftTrack.getTanl()); - auto distanceSq = dx * dx + dy * dy; - auto cutDistanceSq = 9 * (mchTrack.getSigma2X() + mchTrack.getSigma2Y()); - auto cutPhiSq = 9 * (mchTrack.getSigma2Phi() + mftTrack.getSigma2Phi()); - auto cutTanlSq = 9 * (mchTrack.getSigma2Tanl() + mftTrack.getSigma2Tanl()); - return (distanceSq < cutDistanceSq) and (dPhi * dPhi < cutPhiSq) and (dTanl * dTanl < cutTanlSq); + CutFunc_t mCutFunc = [](const GlobalFwdTrack& mchtrack, const TrackParCovFwd& mfttrack) -> bool { + throw std::runtime_error("MatchGlobalFwd: track pair candidate cut function not configured!"); }; /// Converts mchTrack parameters to Forward coordinate system o2::dataformats::GlobalFwdTrack MCHtoFwd(const o2::mch::TrackParam& mchTrack); float mBz = -5.f; ///< nominal Bz in kGauss - float mMatchingPlaneZ = sLastMFTPlaneZ; ///< MCH-MFT matching plane Z position TODO: Make configurable + float mMatchingPlaneZ = sLastMFTPlaneZ; ///< MCH-MFT matching plane Z position Float_t mMFTDiskThicknessInX0 = 0.042 / 5; o2::InteractionRecord mStartIR{0, 0}; ///< IR corresponding to the start of the TF int mMFTROFrameLengthInBC = 0; ///< MFT RO frame in BC (for MFT cont. mode only) float mMFTROFrameLengthMUS = -1.; ///< MFT RO frame in \mus float mMFTROFrameLengthMUSInv = -1.; ///< MFT RO frame in \mus inverse + std::map mMatchingFunctionMap; ///< MFT-MCH Matching function + std::map mCutFunctionMap; ///< MFT-MCH Candidate cut function + o2::BunchFilling mBunchFilling; std::array mClosestBunchAbove; // closest filled bunch from above std::array mClosestBunchBelow; // closest filled bunch from below @@ -329,7 +287,7 @@ class MatchGlobalFwd bool mMCTruthON = false; ///< Flag availability of MC truth bool mUseMIDMCHMatch = false; ///< Flag for using MCHMID matches (TrackMCHMID) bool mSaveAll = false; ///< Flag to save all MCH-MFT matching pairs - MatchingType mMatchingType = MATCHINGFUNC; + MatchingType mMatchingType = MATCHINGUNDEFINED; TGeoManager* mGeoManager; }; diff --git a/Detectors/GlobalTracking/include/GlobalTracking/MatchGlobalFwdParam.h b/Detectors/GlobalTracking/include/GlobalTracking/MatchGlobalFwdParam.h index e7b3bb8807cf2..e0b1cb856a7d4 100644 --- a/Detectors/GlobalTracking/include/GlobalTracking/MatchGlobalFwdParam.h +++ b/Detectors/GlobalTracking/include/GlobalTracking/MatchGlobalFwdParam.h @@ -24,9 +24,10 @@ namespace globaltracking // ** Parameters for Global forward matching // ** struct GlobalFwdMatchingParam : public o2::conf::ConfigurableParamHelper { - bool useMIDMatch = false; + bool useMIDMatch = false; ///< Use input from MCH-MID mathing std::string matchFcn = "matchALL"; ///< MFT-MCH matching score evaluation std::string cutFcn = "cutDisabled"; ///< MFT-MCH candicate cut function + bool MCMatching = false; ///< MFT-MCH matching computed from MCLabels double matchPlaneZ = -77.5; ///< MFT-MCH matching plane z coordinate bool isMatchUpstream() const { return matchFcn.find("matchUpstream") < matchFcn.length(); } diff --git a/Detectors/GlobalTracking/src/MatchGlobalFwd.cxx b/Detectors/GlobalTracking/src/MatchGlobalFwd.cxx index 64387c5b7dbbd..ee91e7ebee866 100644 --- a/Detectors/GlobalTracking/src/MatchGlobalFwd.cxx +++ b/Detectors/GlobalTracking/src/MatchGlobalFwd.cxx @@ -24,49 +24,28 @@ void MatchGlobalFwd::init() mMatchingPlaneZ = matchingParam.matchPlaneZ; LOG(info) << "MFTMCH matchingPlaneZ = " << mMatchingPlaneZ; - auto& matchingFcn = matchingParam.matchFcn; - LOG(info) << "Match function string = " << matchingFcn; - - if (matchingFcn.find("matchMC") < matchingFcn.length()) { - LOG(info) << " ==> Setting MC Label matching: " << matchingFcn; - LOG(info) << " MC Label matching: Matching scores defaults to matchMFT_MCH_TracksAllParam, unless another function is explictly specified."; - mMatchingType = MATCHINGMCLABEL; - setMatchingFunction(&MatchGlobalFwd::matchMFT_MCH_TracksAllParam); - if (!mMCTruthON) { - LOG(fatal) << "Label matching requries MC Labels!"; - } - } + auto& matchingFcnStr = matchingParam.matchFcn; + LOG(info) << "Match function string = " << matchingFcnStr; - if (matchingFcn.find("matchALL") < matchingFcn.length()) { - LOG(info) << " ==> Setting MatchingFunction matchALL."; - setMatchingFunction(&MatchGlobalFwd::matchMFT_MCH_TracksAllParam); - } else if (matchingFcn.find("matchPhiTanlXY") < matchingFcn.length()) { - LOG(info) << " ==> Setting MatchingFunction matchPhiTanlXY."; - setMatchingFunction(&MatchGlobalFwd::matchMFT_MCH_TracksXYPhiTanl); - } else if (matchingFcn.find("matchXY") < matchingFcn.length()) { - LOG(info) << " ==> Setting MatchingFunction matchXY."; - setMatchingFunction(&MatchGlobalFwd::matchMFT_MCH_TracksXY); - } else if (matchingFcn.find("matchHiroshima") < matchingFcn.length()) { - LOG(info) << " ==> Setting MatchingFunction Hiroshima."; - setMatchingFunction(&MatchGlobalFwd::matchHiroshima); - } else if (matchingParam.isMatchUpstream()) { + if (matchingParam.isMatchUpstream()) { LOG(info) << " ==> Setting Upstream matching."; - setMatchingFunction(&MatchGlobalFwd::noMatchFcn); mMatchingType = MATCHINGUPSTREAM; + } else { + if (mMatchingFunctionMap.find(matchingFcnStr) != mMatchingFunctionMap.end()) { + mMatchFunc = mMatchingFunctionMap[matchingFcnStr]; + mMatchingType = MATCHINGFUNC; + LOG(info) << " Found built-in matching function " << matchingFcnStr; + } else { + throw std::invalid_argument("Invalid matching function! Aborting..."); + } } - auto& cutFcn = matchingParam.cutFcn; - LOG(info) << "MFTMCH pair candidate cut function string = " << cutFcn; - - if (cutFcn.find("cutDisabled") < cutFcn.length()) { - LOG(info) << " ==> Setting CutFunction: cutDisabled"; - setCutFunction(&MatchGlobalFwd::cutDisabled); - } else if (cutFcn.find("matchCut3Sigma") < cutFcn.length()) { - LOG(info) << " ==> Setting MatchingFunction matchCut3Sigma"; - setCutFunction(&MatchGlobalFwd::matchCut3Sigma); - } else if (cutFcn.find("matchCut3SigmaXYAngles") < cutFcn.length()) { - LOG(info) << " ==> Setting MatchingFunction matchCut3SigmaXYAngles: "; - setCutFunction(&MatchGlobalFwd::matchCut3SigmaXYAngles); + auto& cutFcnStr = matchingParam.cutFcn; + LOG(info) << "MFTMCH pair candidate cut function string = " << cutFcnStr; + + if (mCutFunctionMap.find(cutFcnStr) != mCutFunctionMap.end()) { + mCutFunc = mCutFunctionMap[cutFcnStr]; + LOG(info) << " Found built-in cut function " << cutFcnStr; } else { throw std::invalid_argument("Invalid cut function! Aborting..."); } @@ -81,6 +60,9 @@ void MatchGlobalFwd::init() //_________________________________________________________ void MatchGlobalFwd::run(const o2::globaltracking::RecoContainer& inp) { + + auto& matchingParam = GlobalFwdMatchingParam::Instance(); + mRecoCont = &inp; mStartIR = inp.startIR; @@ -90,18 +72,19 @@ void MatchGlobalFwd::run(const o2::globaltracking::RecoContainer& inp) return; } - switch (mMatchingType) { - case MATCHINGFUNC: - mSaveAll ? doMatching() : doMatching(); - break; - case MATCHINGUPSTREAM: - loadMatches(); - break; - case MATCHINGMCLABEL: - doMCMatching(); - break; - default: - LOG(fatal) << "Invalid MFTMCH matching mode"; + if (matchingParam.MCMatching) { // MC Label matching + mMCTruthON ? doMCMatching() : throw std::runtime_error("Label matching requries MC Labels!"); + } else { + switch (mMatchingType) { + case MATCHINGFUNC: + mSaveAll ? doMatching() : doMatching(); + break; + case MATCHINGUPSTREAM: + loadMatches(); + break; + default: + LOG(fatal) << "Invalid MFTMCH matching mode"; + } } fitTracks(); @@ -375,14 +358,14 @@ void MatchGlobalFwd::ROFMatch(int MFTROFId, int firstMCHROFId, int lastMCHROFId) if (mMCTruthON) { matchLabel = computeLabel(MCHId, MFTId); } - if (matchingCut(thisMCHTrack, thisMFTTrack)) { + if (mCutFunc(thisMCHTrack, thisMFTTrack)) { thisMCHTrack.countMFTCandidate(); if (mMCTruthON) { if (matchLabel.isCorrect()) { thisMCHTrack.setCloseMatch(); } } - auto chi2 = matchingEval(thisMCHTrack, thisMFTTrack); + auto chi2 = mMatchFunc(thisMCHTrack, thisMFTTrack); if (chi2 < thisMCHTrack.getMFTMCHMatchingChi2()) { thisMCHTrack.setMFTTrackID(MFTId); thisMCHTrack.setMFTMCHMatchingChi2(chi2); @@ -462,7 +445,7 @@ void MatchGlobalFwd::doMCMatching() if (matchLabel.isCorrect()) { nTrue++; thisMCHTrack.setCloseMatch(); - auto chi2 = matchingEval(thisMCHTrack, thisMFTTrack); + auto chi2 = mMatchFunc(thisMCHTrack, thisMFTTrack); thisMCHTrack.setMFTTrackID(MFTId); thisMCHTrack.setMFTMCHMatchingChi2(chi2); thisMCHTrack.setTimeMUS(thisMCHTrack.tBracket.getMin(), thisMCHTrack.tBracket.delta()); @@ -603,18 +586,6 @@ bool MatchGlobalFwd::computeCluster(o2::dataformats::GlobalFwdTrack& track, cons return false; } -//_________________________________________________________________________________________________ -double MatchGlobalFwd::matchingEval(const TrackLocMCH& mchTrack, const TrackLocMFT& mftTrack) -{ - return (this->*mMatchFunc)(mchTrack, mftTrack); -} - -//_________________________________________________________________________________________________ -bool MatchGlobalFwd::matchingCut(const TrackLocMCH& mchTrack, const TrackLocMFT& mftTrack) -{ - return (this->*mCutFunc)(mchTrack, mftTrack); -} - //_________________________________________________________ void MatchGlobalFwd::setMFTROFrameLengthMUS(float fums) { @@ -656,43 +627,6 @@ void MatchGlobalFwd::setBunchFilling(const o2::BunchFilling& bf) } } -//_________________________________________________________________________________________________ -double MatchGlobalFwd::matchMFT_MCH_TracksAllParam(const TrackLocMCH& mchTrack, const TrackLocMFT& mftTrack) -{ - // Match two tracks evaluating all parameters: X,Y, phi, tanl & q/pt - - SMatrix55Sym I = ROOT::Math::SMatrixIdentity(), H_k, V_k; - SVector5 m_k(mftTrack.getX(), mftTrack.getY(), mftTrack.getPhi(), - mftTrack.getTanl(), mftTrack.getInvQPt()), - r_k_kminus1; - SVector5 GlobalMuonTrackParameters = mchTrack.getParameters(); - SMatrix55Sym GlobalMuonTrackCovariances = mchTrack.getCovariances(); - V_k(0, 0) = mftTrack.getCovariances()(0, 0); - V_k(1, 1) = mftTrack.getCovariances()(1, 1); - V_k(2, 2) = mftTrack.getCovariances()(2, 2); - V_k(3, 3) = mftTrack.getCovariances()(3, 3); - V_k(4, 4) = mftTrack.getCovariances()(4, 4); - H_k(0, 0) = 1.0; - H_k(1, 1) = 1.0; - H_k(2, 2) = 1.0; - H_k(3, 3) = 1.0; - H_k(4, 4) = 1.0; - - // Covariance of residuals - SMatrix55Std invResCov = (V_k + ROOT::Math::Similarity(H_k, GlobalMuonTrackCovariances)); - invResCov.Invert(); - - // Kalman Gain Matrix - SMatrix55Std K_k = GlobalMuonTrackCovariances * ROOT::Math::Transpose(H_k) * invResCov; - - // Update Parameters - r_k_kminus1 = m_k - H_k * GlobalMuonTrackParameters; // Residuals of prediction - - auto matchChi2Track = ROOT::Math::Similarity(r_k_kminus1, invResCov); - - return matchChi2Track; -} - //_________________________________________________________________________________________________ o2::dataformats::GlobalFwdTrack MatchGlobalFwd::MCHtoFwd(const o2::mch::TrackParam& mchParam) { @@ -779,38 +713,50 @@ o2::dataformats::GlobalFwdTrack MatchGlobalFwd::MCHtoFwd(const o2::mch::TrackPar } //_________________________________________________________________________________________________ -double MatchGlobalFwd::matchMFT_MCH_TracksXY(const TrackLocMCH& mchTrack, const TrackLocMFT& mftTrack) +MatchGlobalFwd::MatchGlobalFwd() { - // Calculate Matching Chi2 - X and Y positions - SMatrix55Sym I = ROOT::Math::SMatrixIdentity(); - SMatrix25 H_k; - SMatrix22 V_k; - SVector2 m_k(mftTrack.getX(), mftTrack.getY()), r_k_kminus1; - SVector5 GlobalMuonTrackParameters = mchTrack.getParameters(); - SMatrix55Sym GlobalMuonTrackCovariances = mchTrack.getCovariances(); - V_k(0, 0) = mftTrack.getCovariances()(0, 0); - V_k(1, 1) = mftTrack.getCovariances()(1, 1); - H_k(0, 0) = 1.0; - H_k(1, 1) = 1.0; + // Define built-in matching functions - // Covariance of residuals - SMatrix22 invResCov = (V_k + ROOT::Math::Similarity(H_k, GlobalMuonTrackCovariances)); - invResCov.Invert(); + //________________________________________________________________________________ + mMatchingFunctionMap["matchALL"] = [](const GlobalFwdTrack& mchTrack, const TrackParCovFwd& mftTrack) -> double { + // Match two tracks evaluating all parameters: X,Y, phi, tanl & q/pt - // Kalman Gain Matrix - SMatrix52 K_k = GlobalMuonTrackCovariances * ROOT::Math::Transpose(H_k) * invResCov; + SMatrix55Sym I = ROOT::Math::SMatrixIdentity(), H_k, V_k; + SVector5 m_k(mftTrack.getX(), mftTrack.getY(), mftTrack.getPhi(), + mftTrack.getTanl(), mftTrack.getInvQPt()), + r_k_kminus1; + SVector5 GlobalMuonTrackParameters = mchTrack.getParameters(); + SMatrix55Sym GlobalMuonTrackCovariances = mchTrack.getCovariances(); + V_k(0, 0) = mftTrack.getCovariances()(0, 0); + V_k(1, 1) = mftTrack.getCovariances()(1, 1); + V_k(2, 2) = mftTrack.getCovariances()(2, 2); + V_k(3, 3) = mftTrack.getCovariances()(3, 3); + V_k(4, 4) = mftTrack.getCovariances()(4, 4); + H_k(0, 0) = 1.0; + H_k(1, 1) = 1.0; + H_k(2, 2) = 1.0; + H_k(3, 3) = 1.0; + H_k(4, 4) = 1.0; - // Residuals of prediction - r_k_kminus1 = m_k - H_k * GlobalMuonTrackParameters; - auto matchChi2Track = ROOT::Math::Similarity(r_k_kminus1, invResCov); + // Covariance of residuals + SMatrix55Std invResCov = (V_k + ROOT::Math::Similarity(H_k, GlobalMuonTrackCovariances)); + invResCov.Invert(); - return matchChi2Track; -} + // Kalman Gain Matrix + SMatrix55Std K_k = GlobalMuonTrackCovariances * ROOT::Math::Transpose(H_k) * invResCov; + + // Update Parameters + r_k_kminus1 = m_k - H_k * GlobalMuonTrackParameters; // Residuals of prediction + + auto matchChi2Track = ROOT::Math::Similarity(r_k_kminus1, invResCov); + + return matchChi2Track; + }; + + //________________________________________________________________________________ + mMatchingFunctionMap["matchsXYPhiTanl"] = [](const GlobalFwdTrack& mchTrack, const TrackParCovFwd& mftTrack) -> double { -//_________________________________________________________________________________________________ -double MatchGlobalFwd::matchMFT_MCH_TracksXYPhiTanl(const TrackLocMCH& mchTrack, const TrackLocMFT& mftTrack) -{ // Match two tracks evaluating positions & angles SMatrix55Sym I = ROOT::Math::SMatrixIdentity(); @@ -842,14 +788,41 @@ double MatchGlobalFwd::matchMFT_MCH_TracksXYPhiTanl(const TrackLocMCH& mchTrack, auto matchChi2Track = ROOT::Math::Similarity(r_k_kminus1, invResCov); - return matchChi2Track; -} + return matchChi2Track; }; -//_________________________________________________________________________________________________ -double MatchGlobalFwd::matchHiroshima(const TrackLocMCH& mchTrack, const TrackLocMFT& mftTrack) -{ + //________________________________________________________________________________ + mMatchingFunctionMap["matchXY"] = [](const GlobalFwdTrack& mchTrack, const TrackParCovFwd& mftTrack) -> double { + + // Calculate Matching Chi2 - X and Y positions + + SMatrix55Sym I = ROOT::Math::SMatrixIdentity(); + SMatrix25 H_k; + SMatrix22 V_k; + SVector2 m_k(mftTrack.getX(), mftTrack.getY()), r_k_kminus1; + SVector5 GlobalMuonTrackParameters = mchTrack.getParameters(); + SMatrix55Sym GlobalMuonTrackCovariances = mchTrack.getCovariances(); + V_k(0, 0) = mftTrack.getCovariances()(0, 0); + V_k(1, 1) = mftTrack.getCovariances()(1, 1); + H_k(0, 0) = 1.0; + H_k(1, 1) = 1.0; + + // Covariance of residuals + SMatrix22 invResCov = (V_k + ROOT::Math::Similarity(H_k, GlobalMuonTrackCovariances)); + invResCov.Invert(); - //Hiroshima's Matching function + // Kalman Gain Matrix + SMatrix52 K_k = GlobalMuonTrackCovariances * ROOT::Math::Transpose(H_k) * invResCov; + + // Residuals of prediction + r_k_kminus1 = m_k - H_k * GlobalMuonTrackParameters; + auto matchChi2Track = ROOT::Math::Similarity(r_k_kminus1, invResCov); + + return matchChi2Track; }; + + //________________________________________________________________________________ + mMatchingFunctionMap["matchNeedsName"] = [this](const GlobalFwdTrack& mchTrack, const TrackParCovFwd& mftTrack) -> double { + + //Hiroshima's Matching function needs a physics-based name //Matching constants Double_t LAbs = 415.; //Absorber Length[cm] @@ -906,5 +879,40 @@ double MatchGlobalFwd::matchHiroshima(const TrackLocMCH& mchTrack, const TrackLo auto scoreY = TMath::Sqrt(dycircle * dycircle + dthetaycircle * dthetaycircle); auto score = TMath::Sqrt(scoreX * scoreX + scoreY * scoreY); - return score; -}; + return score; }; + + // Define built-in candidate cut functions + + //________________________________________________________________________________ + mCutFunctionMap["cutDisabled"] = [](const GlobalFwdTrack& mchTrack, const TrackParCovFwd& mftTrack) -> bool { + return true; + }; + + //________________________________________________________________________________ + mCutFunctionMap["cut3Sigma"] = [](const GlobalFwdTrack& mchTrack, const TrackParCovFwd& mftTrack) -> bool { + auto dx = mchTrack.getX() - mftTrack.getX(); + auto dy = mchTrack.getY() - mftTrack.getY(); + auto dPhi = mchTrack.getPhi() - mftTrack.getPhi(); + auto dTanl = TMath::Abs(mchTrack.getTanl() - mftTrack.getTanl()); + auto dInvQPt = TMath::Abs(mchTrack.getInvQPt() - mftTrack.getInvQPt()); + auto distanceSq = dx * dx + dy * dy; + auto cutDistanceSq = 9 * (mchTrack.getSigma2X() + mchTrack.getSigma2Y()); + auto cutPhiSq = 9 * (mchTrack.getSigma2Phi() + mftTrack.getSigma2Phi()); + auto cutTanlSq = 9 * (mchTrack.getSigma2Tanl() + mftTrack.getSigma2Tanl()); + auto cutInvQPtSq = 9 * (mchTrack.getSigma2InvQPt() + mftTrack.getSigma2InvQPt()); + return (distanceSq < cutDistanceSq) and (dPhi * dPhi < cutPhiSq) and (dTanl * dTanl < cutTanlSq) and (dInvQPt * dInvQPt < cutInvQPtSq); + }; + + //________________________________________________________________________________ + mCutFunctionMap["cut3SigmaXYAngles"] = [](const GlobalFwdTrack& mchTrack, const TrackParCovFwd& mftTrack) -> bool { + auto dx = mchTrack.getX() - mftTrack.getX(); + auto dy = mchTrack.getY() - mftTrack.getY(); + auto dPhi = mchTrack.getPhi() - mftTrack.getPhi(); + auto dTanl = TMath::Abs(mchTrack.getTanl() - mftTrack.getTanl()); + auto distanceSq = dx * dx + dy * dy; + auto cutDistanceSq = 9 * (mchTrack.getSigma2X() + mchTrack.getSigma2Y()); + auto cutPhiSq = 9 * (mchTrack.getSigma2Phi() + mftTrack.getSigma2Phi()); + auto cutTanlSq = 9 * (mchTrack.getSigma2Tanl() + mftTrack.getSigma2Tanl()); + return (distanceSq < cutDistanceSq) and (dPhi * dPhi < cutPhiSq) and (dTanl * dTanl < cutTanlSq); + }; +} From 5990194a08460f6e174a5ab916fb092d7bae87f2 Mon Sep 17 00:00:00 2001 From: Rafael Pezzi Date: Tue, 30 Nov 2021 19:18:56 +0100 Subject: [PATCH 2/2] Option to load MFTMCH matching function from external file --- .../include/GlobalTracking/MatchGlobalFwd.h | 51 ++++++++++++------- .../GlobalTracking/MatchGlobalFwdParam.h | 13 +++-- .../GlobalTracking/src/MatchGlobalFwd.cxx | 3 ++ 3 files changed, 45 insertions(+), 22 deletions(-) diff --git a/Detectors/GlobalTracking/include/GlobalTracking/MatchGlobalFwd.h b/Detectors/GlobalTracking/include/GlobalTracking/MatchGlobalFwd.h index 83071eb71ecc5..5c45ac0360aa7 100644 --- a/Detectors/GlobalTracking/include/GlobalTracking/MatchGlobalFwd.h +++ b/Detectors/GlobalTracking/include/GlobalTracking/MatchGlobalFwd.h @@ -10,7 +10,7 @@ // or submit itself to any jurisdiction. /// \file MatchGlobalFwd.h -/// \brief Class to perform MCH MFT matching +/// \brief Class to perform MFT MCH (and MID) matching /// \author rafael.pezzi@cern.ch #ifndef ALICEO2_GLOBTRACKING_MATCHGLOBALFWD_ @@ -23,6 +23,7 @@ #include #include #include "CommonConstants/LHCConstants.h" +#include "CommonUtils/ConfigurationMacroHelper.h" #include "CommonDataFormat/BunchFilling.h" #include "ITSMFTReconstruction/ChipMappingMFT.h" #include "MCHTracking/TrackExtrap.h" @@ -83,25 +84,25 @@ using o2::track::TrackParCovFwd; typedef std::function MatchingFunc_t; typedef std::function CutFunc_t; +using MFTCluster = o2::BaseCluster; +using BracketF = o2::math_utils::Bracket; +using SMatrix55Std = ROOT::Math::SMatrix; +using SMatrix55Sym = ROOT::Math::SMatrix>; + +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 SMatrix54 = ROOT::Math::SMatrix; +using SMatrix22 = ROOT::Math::SMatrix; +using SMatrix25 = ROOT::Math::SMatrix; +using SMatrix52 = ROOT::Math::SMatrix; + class MatchGlobalFwd { public: - using MFTCluster = o2::BaseCluster; - using BracketF = o2::math_utils::Bracket; - using SMatrix55Std = ROOT::Math::SMatrix; - using SMatrix55Sym = ROOT::Math::SMatrix>; - - 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 SMatrix54 = ROOT::Math::SMatrix; - using SMatrix22 = ROOT::Math::SMatrix; - using SMatrix25 = ROOT::Math::SMatrix; - using SMatrix52 = ROOT::Math::SMatrix; - enum MatchingType : uint8_t { ///< MFT-MCH matching modes MATCHINGFUNC, ///< Matching function-based MFT-MCH track matching MATCHINGUPSTREAM, ///< MFT-MCH track matching loaded from input file @@ -238,6 +239,22 @@ class MatchGlobalFwd throw std::runtime_error("MatchGlobalFwd: track pair candidate cut function not configured!"); }; + bool loadExternalMatchingFunction() + { + // Loads MFTMCH Matching function from external file + + auto& matchingParam = GlobalFwdMatchingParam::Instance(); + + const auto& extFuncMacroFile = matchingParam.extMatchFuncFile; + const auto& extFuncName = matchingParam.extMatchFuncName; + + LOG(info) << "Loading external MFTMCH matching function: function name = " << extFuncName << " ; Filename = " << extFuncMacroFile; + + auto func = o2::conf::GetFromMacro(extFuncMacroFile.c_str(), extFuncName.c_str(), "o2::globaltracking::MatchingFunc_t*", "mtcFcn"); + mMatchFunc = (*func); + return true; + } + /// Converts mchTrack parameters to Forward coordinate system o2::dataformats::GlobalFwdTrack MCHtoFwd(const o2::mch::TrackParam& mchTrack); diff --git a/Detectors/GlobalTracking/include/GlobalTracking/MatchGlobalFwdParam.h b/Detectors/GlobalTracking/include/GlobalTracking/MatchGlobalFwdParam.h index e0b1cb856a7d4..cf1cd99d0cc91 100644 --- a/Detectors/GlobalTracking/include/GlobalTracking/MatchGlobalFwdParam.h +++ b/Detectors/GlobalTracking/include/GlobalTracking/MatchGlobalFwdParam.h @@ -24,13 +24,16 @@ namespace globaltracking // ** Parameters for Global forward matching // ** struct GlobalFwdMatchingParam : public o2::conf::ConfigurableParamHelper { - bool useMIDMatch = false; ///< Use input from MCH-MID mathing - std::string matchFcn = "matchALL"; ///< MFT-MCH matching score evaluation - std::string cutFcn = "cutDisabled"; ///< MFT-MCH candicate cut function - bool MCMatching = false; ///< MFT-MCH matching computed from MCLabels - double matchPlaneZ = -77.5; ///< MFT-MCH matching plane z coordinate + bool useMIDMatch = false; ///< Use input from MCH-MID mathing + std::string matchFcn = "matchALL"; ///< MFT-MCH matching score evaluation + std::string extMatchFuncFile = "FwdMatchFunc.C"; ///< File name for external input matching function + std::string extMatchFuncName = "getMatchingFunction()"; ///< Name of external matching function getter + std::string cutFcn = "cutDisabled"; ///< MFT-MCH candicate cut function + bool MCMatching = false; ///< MFT-MCH matching computed from MCLabels + double matchPlaneZ = -77.5; ///< MFT-MCH matching plane z coordinate bool isMatchUpstream() const { return matchFcn.find("matchUpstream") < matchFcn.length(); } + bool matchingExternalFunction() const { return matchFcn.find("matchExternal") < matchFcn.length(); } bool saveAll = false; ///< Option to save all MFTMCH pair combinations. O2ParamDef(GlobalFwdMatchingParam, "FwdMatching"); diff --git a/Detectors/GlobalTracking/src/MatchGlobalFwd.cxx b/Detectors/GlobalTracking/src/MatchGlobalFwd.cxx index ee91e7ebee866..896496db2e2f0 100644 --- a/Detectors/GlobalTracking/src/MatchGlobalFwd.cxx +++ b/Detectors/GlobalTracking/src/MatchGlobalFwd.cxx @@ -30,6 +30,9 @@ void MatchGlobalFwd::init() if (matchingParam.isMatchUpstream()) { LOG(info) << " ==> Setting Upstream matching."; mMatchingType = MATCHINGUPSTREAM; + } else if (matchingParam.matchingExternalFunction()) { + loadExternalMatchingFunction(); + mMatchingType = MATCHINGFUNC; } else { if (mMatchingFunctionMap.find(matchingFcnStr) != mMatchingFunctionMap.end()) { mMatchFunc = mMatchingFunctionMap[matchingFcnStr];