diff --git a/PWGDQ/Tasks/qaMatching.cxx b/PWGDQ/Tasks/qaMatching.cxx index 9611f690500..20277a78696 100644 --- a/PWGDQ/Tasks/qaMatching.cxx +++ b/PWGDQ/Tasks/qaMatching.cxx @@ -69,7 +69,6 @@ #include #include #include -#include #include #include #include @@ -85,6 +84,77 @@ using namespace o2; using namespace o2::framework; using namespace o2::aod; +namespace qamatching +{ +DECLARE_SOA_COLUMN(P, p, float); +DECLARE_SOA_COLUMN(Pt, pt, float); +DECLARE_SOA_COLUMN(Eta, eta, float); +DECLARE_SOA_COLUMN(Phi, phi, float); +DECLARE_SOA_COLUMN(MatchLabel, matchLabel, int8_t); +DECLARE_SOA_COLUMN(TrackId, trackId, int64_t); +DECLARE_SOA_COLUMN(MatchType, matchType, int8_t); +DECLARE_SOA_COLUMN(MatchScore, matchScore, float); +DECLARE_SOA_COLUMN(MatchRanking, matchRanking, int32_t); +DECLARE_SOA_COLUMN(MftMultiplicity, mftMultiplicity, int32_t); +DECLARE_SOA_COLUMN(TrackType, trackType, int8_t); +DECLARE_SOA_COLUMN(MftMatchAttempts, mftMatchAttempts, int32_t); +DECLARE_SOA_COLUMN(XAtVtx, xAtVtx, float); +DECLARE_SOA_COLUMN(YAtVtx, yAtVtx, float); +DECLARE_SOA_COLUMN(ZAtVtx, zAtVtx, float); +DECLARE_SOA_COLUMN(PxAtVtx, pxAtVtx, float); +DECLARE_SOA_COLUMN(PyAtVtx, pyAtVtx, float); +DECLARE_SOA_COLUMN(PzAtVtx, pzAtVtx, float); +DECLARE_SOA_COLUMN(ColX, colX, float); +DECLARE_SOA_COLUMN(ColY, colY, float); +DECLARE_SOA_COLUMN(ColZ, colZ, float); +} // namespace qamatching + +namespace o2::aod +{ +DECLARE_SOA_TABLE(QaMatchingEvents, "AOD", "QAMEVT", + o2::soa::Index<>, + qamatching::MftMultiplicity, + qamatching::ColX, + qamatching::ColY, + qamatching::ColZ); +} // namespace o2::aod + +namespace qamatching +{ +DECLARE_SOA_INDEX_COLUMN_FULL(ReducedEvent, reducedEvent, int32_t, o2::aod::QaMatchingEvents, ""); +} // namespace qamatching + +namespace o2::aod +{ +DECLARE_SOA_TABLE(QaMatchingMCHTrack, "AOD", "QAMCHTRK", + qamatching::ReducedEventId, + qamatching::TrackId, + qamatching::TrackType, + qamatching::P, + qamatching::Pt, + qamatching::Eta, + qamatching::Phi, + qamatching::MftMatchAttempts, + qamatching::XAtVtx, + qamatching::YAtVtx, + qamatching::ZAtVtx, + qamatching::PxAtVtx, + qamatching::PyAtVtx, + qamatching::PzAtVtx); +DECLARE_SOA_TABLE(QaMatchingCandidates, "AOD", "QAMCAND", + qamatching::ReducedEventId, + qamatching::MatchLabel, + qamatching::TrackId, + qamatching::P, qamatching::Pt, qamatching::Eta, qamatching::Phi, + qamatching::MatchType, qamatching::MatchScore, qamatching::MatchRanking, + qamatching::XAtVtx, + qamatching::YAtVtx, + qamatching::ZAtVtx, + qamatching::PxAtVtx, + qamatching::PyAtVtx, + qamatching::PzAtVtx); +} // namespace o2::aod + using MyEvents = soa::Join; using MyMuons = soa::Join; using MyMuonsMC = soa::Join; @@ -98,20 +168,20 @@ using MyMFT = MyMFTs::iterator; using MyMFTCovariance = MyMFTCovariances::iterator; using SMatrix55 = ROOT::Math::SMatrix>; -using SMatrix5 = ROOT::Math::SVector; +using SMatrix5 = ROOT::Math::SVector; static float chi2ToScore(float chi2, int ndf, float chi2max) { - double p = -TMath::Log10(ROOT::Math::chisquared_cdf_c(chi2, ndf)); - double pnorm = -TMath::Log10(ROOT::Math::chisquared_cdf_c(chi2max, ndf)); + double p = -std::log10(ROOT::Math::chisquared_cdf_c(chi2, ndf)); + double pnorm = -std::log10(ROOT::Math::chisquared_cdf_c(chi2max, ndf)); double result = (1.f / (p / pnorm + 1.f)); return static_cast(result); } -struct qaMatching { +struct QaMatching { template - using matrix = std::array, nr>; + using Matrix = std::array, nr>; enum MuonMatchType { kMatchTypeTrueLeading = 0, @@ -125,6 +195,22 @@ struct qaMatching { kMatchTypeUndefined }; + static constexpr int GlobalTrackTypeMax = 2; + static constexpr int MchMidTrackType = 3; + static constexpr int FirstDecayMotherRank = 2; + static constexpr int MftTrackTypeStandard = 0; + static constexpr int MftTrackTypeCA = 1; + static constexpr int ThetaAbsBoundaryDeg = 3; + static constexpr double SlopeResolutionZ = 535.; + static constexpr int MatchingDegreesOfFreedom = 5; + static constexpr size_t MinCandidatesForDeltaChi2 = 2; + static constexpr float MatchingScoreChi2Max = 50.f; + static constexpr float InvalidDeltaChi2 = -1.f; + static constexpr int ExtrapolationMethodStandard = 0; + static constexpr int ExtrapolationMethodMftFirstPoint = 2; + static constexpr int ExtrapolationMethodVertex = 3; + static constexpr int ExtrapolationMethodMftDca = 4; + struct MatchingCandidate { int64_t collisionId{-1}; int64_t globalTrackId{-1}; @@ -141,170 +227,171 @@ struct qaMatching { }; //// 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, ""}; + Configurable cfgPMchLow{"cfgPMchLow", 0.0f, ""}; + Configurable cfgPtMchLow{"cfgPtMchLow", 0.7f, ""}; + Configurable cfgEtaMchLow{"cfgEtaMchLow", -4.0f, ""}; + Configurable cfgEtaMchUp{"cfgEtaMchUp", -2.5f, ""}; + Configurable cfgRabsLow{"cfgRabsLow", 17.6f, ""}; + Configurable cfgRabsUp{"cfgRabsUp", 89.5f, ""}; + Configurable cfgPdcaUp{"cfgPdcaUp", 6.f, ""}; + Configurable cfgTrackChi2MchUp{"cfgTrackChi2MchUp", 5.f, ""}; + Configurable cfgMatchingChi2MchMidUp{"cfgMatchingChi2MchMidUp", 999.f, ""}; //// Variables for selecting mft tracks - Configurable fEtaMftLow{"cfgEtaMftlow", -3.6f, ""}; - Configurable fEtaMftUp{"cfgEtaMftup", -2.5f, ""}; - Configurable fTrackNClustMftLow{"cfgTrackNClustMftLow", 7, ""}; - Configurable fTrackChi2MftUp{"cfgTrackChi2MftUp", 999.f, ""}; + Configurable cfgEtaMftLow{"cfgEtaMftLow", -3.6f, ""}; + Configurable cfgEtaMftUp{"cfgEtaMftUp", -2.5f, ""}; + Configurable cfgTrackNClustMftLow{"cfgTrackNClustMftLow", 7, ""}; + Configurable cfgTrackChi2MftUp{"cfgTrackChi2MftUp", 999.f, ""}; //// Variables for selecting global tracks - Configurable fMatchingChi2ScoreMftMchLow{"cfgMatchingChi2ScoreMftMchLow", chi2ToScore(50.f, 5, 50.f), ""}; + Configurable cfgMatchingChi2ScoreMftMchLow{"cfgMatchingChi2ScoreMftMchLow", chi2ToScore(50.f, 5, 50.f), ""}; //// Variables for selecting tagged muons - Configurable fMuonTaggingNCrossedMftPlanesLow{"cfgMuonTaggingNCrossedMftPlanesLow", 5, ""}; - Configurable fMuonTaggingTrackChi2MchUp{"cfgMuonTaggingTrackChi2MchUp", 5.f, ""}; - Configurable fMuonTaggingPMchLow{"cfgMuonTaggingPMchLow", 0.0f, ""}; - Configurable fMuonTaggingPtMchLow{"cfgMuonTaggingPtMchLow", 0.7f, ""}; - Configurable fMuonTaggingEtaMchLow{"cfgMuonTaggingEtaMchLow", -3.6f, ""}; - Configurable fMuonTaggingEtaMchUp{"cfgMuonTaggingEtaMchUp", -2.5f, ""}; - Configurable fMuonTaggingRabsLow{"cfgMuonTaggingRabsLow", 17.6f, ""}; - Configurable fMuonTaggingRabsUp{"cfgMuonTaggingRabsUp", 89.5f, ""}; - Configurable fMuonTaggingSigmaPdcaUp{"cfgMuonTaggingPdcaUp", 4.f, ""}; - Configurable fMuonTaggingChi2DiffLow{"cfgMuonTaggingChi2DiffLow", 100.f, ""}; + 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 cfgMuonTaggingChi2DiffLow{"cfgMuonTaggingChi2DiffLow", 100.f, ""}; //// Variables for ccdb - Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable ccdbUrl{"ccdbUrl", "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 grpMagPath{"grpMagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; Configurable geoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"}; // CCDB connection configurables struct : ConfigurableGroup { - Configurable fConfigCcdbUrl{"ccdb-url-", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; - 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 fConfigGrpPath{"grpPath-", "GLO/GRP/GRP", "Path of the grp file"}; - Configurable fConfigGeoPath{"geoPath-", "GLO/Config/GeometryAligned", "Path of the geometry file"}; - Configurable fConfigGrpMagPath{"grpmagPath-", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; - } fConfigCCDB; + 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; struct : ConfigurableGroup { - Configurable fCreatePdgMomHistograms{"cfgCreatePdgMomHistograms", false, "create matching characteristics plots with particle mom PDG codes"}; - } fConfigQAs; + Configurable cfgCreatePdgMomHistograms{"cfgCreatePdgMomHistograms", false, "create matching characteristics plots with particle mom PDG codes"}; + } configQas; /// Variables for histograms configuration - Configurable fNCandidatesMax{"cfgNCandidatesMax", 5, "Number of matching candidates stored for each muon track"}; - Configurable fMftTrackMultiplicityMax{"cfgMftTrackMultiplicityMax", 1000, "Maximum number of MFT tracks per collision"}; + Configurable cfgNCandidatesMax{"cfgNCandidatesMax", 5, "Number of matching candidates stored for each muon track"}; + Configurable cfgMftTrackMultiplicityMax{"cfgMftTrackMultiplicityMax", 1000, "Maximum number of MFT tracks per collision"}; + Configurable cfgQaMatchingAodDebug{"cfgQaMatchingAodDebug", 0, "If >0, print AO2D filling debug (0=off, N=max collisions)"}; double mBzAtMftCenter{0}; o2::globaltracking::MatchGlobalFwd mExtrap; - using MatchingFunc_t = std::function(const o2::dataformats::GlobalFwdTrack& mchtrack, const o2::track::TrackParCovFwd& mfttrack)>; - std::map mMatchingFunctionMap; ///< MFT-MCH Matching function + using MatchingFunc = std::function(const o2::dataformats::GlobalFwdTrack& mchtrack, const o2::track::TrackParCovFwd& mfttrack)>; + std::map mMatchingFunctionMap; ///< MFT-MCH Matching function // Chi2 matching interface - static constexpr int sChi2FunctionsNum = 5; + static constexpr int Chi2FunctionsNum = 5; struct : ConfigurableGroup { - Configurable fFunctionLabel_1{"cfgChi2FunctionLabel_1", std::string{"ProdAll"}, "Text label identifying this chi2 matching method"}; - Configurable fFunctionLabel_2{"cfgChi2FunctionLabel_2", std::string{"MatchXYPhiTanlMom"}, "Text label identifying this chi2 matching method"}; - Configurable fFunctionLabel_3{"cfgChi2FunctionLabel_3", std::string{"MatchXYPhiTanl"}, "Text label identifying this chi2 matching method"}; - Configurable fFunctionLabel_4{"cfgChi2FunctionLabel_4", std::string{""}, "Text label identifying this chi2 matching method"}; - Configurable fFunctionLabel_5{"cfgChi2FunctionLabel_5", std::string{""}, "Text label identifying this chi2 matching method"}; - std::array*, sChi2FunctionsNum> fFunctionLabel{ - &fFunctionLabel_1, &fFunctionLabel_2, &fFunctionLabel_3, &fFunctionLabel_4, &fFunctionLabel_5}; - - Configurable fFunctionNames_1{"cfgChi2FunctionNames_1", std::string{"prod"}, "Name of the chi2 matching function"}; - Configurable fFunctionNames_2{"cfgChi2FunctionNames_2", std::string{"matchALL"}, "Name of the chi2 matching function"}; - Configurable fFunctionNames_3{"cfgChi2FunctionNames_3", std::string{"matchXYPhiTanl"}, "Name of the chi2 matching function"}; - Configurable fFunctionNames_4{"cfgChi2FunctionNames_4", std::string{""}, "Name of the chi2 matching function"}; - Configurable fFunctionNames_5{"cfgChi2FunctionNames_5", std::string{""}, "Name of the chi2 matching function"}; - std::array*, sChi2FunctionsNum> fFunctionName{ - &fFunctionNames_1, &fFunctionNames_2, &fFunctionNames_3, &fFunctionNames_4, &fFunctionNames_5}; - - Configurable fMatchingScoreCut_1{"cfgChi2FunctionMatchingScoreCut_1", 0.f, "Minimum score value for selecting good matches"}; - Configurable fMatchingScoreCut_2{"cfgChi2FunctionMatchingScoreCut_2", 0.5f, "Minimum score value for selecting good matches"}; - Configurable fMatchingScoreCut_3{"cfgChi2FunctionMatchingScoreCut_3", 0.5f, "Minimum score value for selecting good matches"}; - Configurable fMatchingScoreCut_4{"cfgChi2FunctionMatchingScoreCut_4", 0.5f, "Minimum score value for selecting good matches"}; - Configurable fMatchingScoreCut_5{"cfgChi2FunctionMatchingScoreCut_5", 0.5f, "Minimum score value for selecting good matches"}; - std::array*, sChi2FunctionsNum> fMatchingScoreCut{ - &fMatchingScoreCut_1, &fMatchingScoreCut_2, &fMatchingScoreCut_3, &fMatchingScoreCut_4, &fMatchingScoreCut_5}; - - Configurable fMatchingPlaneZ_1{"cfgChi2FunctionMatchingPlaneZ_1", static_cast(o2::mft::constants::mft::LayerZCoordinate()[9]), "Z position of the matching plane"}; - Configurable fMatchingPlaneZ_2{"cfgChi2FunctionMatchingPlaneZ_2", static_cast(o2::mft::constants::mft::LayerZCoordinate()[9]), "Z position of the matching plane"}; - Configurable fMatchingPlaneZ_3{"cfgChi2FunctionMatchingPlaneZ_3", static_cast(o2::mft::constants::mft::LayerZCoordinate()[9]), "Z position of the matching plane"}; - Configurable fMatchingPlaneZ_4{"cfgChi2FunctionMatchingPlaneZ_4", static_cast(o2::mft::constants::mft::LayerZCoordinate()[9]), "Z position of the matching plane"}; - Configurable fMatchingPlaneZ_5{"cfgChi2FunctionMatchingPlaneZ_5", static_cast(o2::mft::constants::mft::LayerZCoordinate()[9]), "Z position of the matching plane"}; - std::array*, sChi2FunctionsNum> fMatchingPlaneZ{ - &fMatchingPlaneZ_1, &fMatchingPlaneZ_2, &fMatchingPlaneZ_3, &fMatchingPlaneZ_4, &fMatchingPlaneZ_5}; - - Configurable fMatchingExtrapMethod_1{"cfgChi2MatchingExtrapMethod_1", static_cast(0), "Method for MCH track extrapolation to maching plane"}; - Configurable fMatchingExtrapMethod_2{"cfgChi2MatchingExtrapMethod_2", static_cast(0), "Method for MCH track extrapolation to maching plane"}; - Configurable fMatchingExtrapMethod_3{"cfgChi2MatchingExtrapMethod_3", static_cast(0), "Method for MCH track extrapolation to maching plane"}; - Configurable fMatchingExtrapMethod_4{"cfgChi2MatchingExtrapMethod_4", static_cast(0), "Method for MCH track extrapolation to maching plane"}; - Configurable fMatchingExtrapMethod_5{"cfgChi2MatchingExtrapMethod_5", static_cast(0), "Method for MCH track extrapolation to maching plane"}; - std::array*, sChi2FunctionsNum> fMatchingExtrapMethod{ - &fMatchingExtrapMethod_1, &fMatchingExtrapMethod_2, &fMatchingExtrapMethod_3, &fMatchingExtrapMethod_4, &fMatchingExtrapMethod_5}; - } fConfigChi2MatchingOptions; + Configurable cfgChi2FunctionLabel1{"cfgChi2FunctionLabel1", std::string{"ProdAll"}, "Text label identifying this chi2 matching method"}; + Configurable cfgChi2FunctionLabel2{"cfgChi2FunctionLabel2", std::string{"MatchXYPhiTanlMom"}, "Text label identifying this chi2 matching method"}; + Configurable cfgChi2FunctionLabel3{"cfgChi2FunctionLabel3", std::string{"MatchXYPhiTanl"}, "Text label identifying this chi2 matching method"}; + Configurable cfgChi2FunctionLabel4{"cfgChi2FunctionLabel4", std::string{""}, "Text label identifying this chi2 matching method"}; + Configurable cfgChi2FunctionLabel5{"cfgChi2FunctionLabel5", std::string{""}, "Text label identifying this chi2 matching method"}; + std::array*, Chi2FunctionsNum> functionLabels{ + &cfgChi2FunctionLabel1, &cfgChi2FunctionLabel2, &cfgChi2FunctionLabel3, &cfgChi2FunctionLabel4, &cfgChi2FunctionLabel5}; + + Configurable cfgChi2FunctionName1{"cfgChi2FunctionName1", std::string{"prod"}, "Name of the chi2 matching function"}; + Configurable cfgChi2FunctionName2{"cfgChi2FunctionName2", std::string{"matchALL"}, "Name of the chi2 matching function"}; + Configurable cfgChi2FunctionName3{"cfgChi2FunctionName3", std::string{"matchXYPhiTanl"}, "Name of the chi2 matching function"}; + Configurable cfgChi2FunctionName4{"cfgChi2FunctionName4", std::string{""}, "Name of the chi2 matching function"}; + Configurable cfgChi2FunctionName5{"cfgChi2FunctionName5", std::string{""}, "Name of the chi2 matching function"}; + std::array*, Chi2FunctionsNum> functionNames{ + &cfgChi2FunctionName1, &cfgChi2FunctionName2, &cfgChi2FunctionName3, &cfgChi2FunctionName4, &cfgChi2FunctionName5}; + + Configurable cfgChi2FunctionMatchingScoreCut1{"cfgChi2FunctionMatchingScoreCut1", 0.f, "Minimum score value for selecting good matches"}; + Configurable cfgChi2FunctionMatchingScoreCut2{"cfgChi2FunctionMatchingScoreCut2", 0.5f, "Minimum score value for selecting good matches"}; + Configurable cfgChi2FunctionMatchingScoreCut3{"cfgChi2FunctionMatchingScoreCut3", 0.5f, "Minimum score value for selecting good matches"}; + Configurable cfgChi2FunctionMatchingScoreCut4{"cfgChi2FunctionMatchingScoreCut4", 0.5f, "Minimum score value for selecting good matches"}; + Configurable cfgChi2FunctionMatchingScoreCut5{"cfgChi2FunctionMatchingScoreCut5", 0.5f, "Minimum score value for selecting good matches"}; + std::array*, Chi2FunctionsNum> matchingScoreCuts{ + &cfgChi2FunctionMatchingScoreCut1, &cfgChi2FunctionMatchingScoreCut2, &cfgChi2FunctionMatchingScoreCut3, &cfgChi2FunctionMatchingScoreCut4, &cfgChi2FunctionMatchingScoreCut5}; + + Configurable cfgChi2FunctionMatchingPlaneZ1{"cfgChi2FunctionMatchingPlaneZ1", static_cast(o2::mft::constants::mft::LayerZCoordinate()[9]), "Z position of the matching plane"}; + Configurable cfgChi2FunctionMatchingPlaneZ2{"cfgChi2FunctionMatchingPlaneZ2", static_cast(o2::mft::constants::mft::LayerZCoordinate()[9]), "Z position of the matching plane"}; + Configurable cfgChi2FunctionMatchingPlaneZ3{"cfgChi2FunctionMatchingPlaneZ3", static_cast(o2::mft::constants::mft::LayerZCoordinate()[9]), "Z position of the matching plane"}; + Configurable cfgChi2FunctionMatchingPlaneZ4{"cfgChi2FunctionMatchingPlaneZ4", static_cast(o2::mft::constants::mft::LayerZCoordinate()[9]), "Z position of the matching plane"}; + Configurable cfgChi2FunctionMatchingPlaneZ5{"cfgChi2FunctionMatchingPlaneZ5", static_cast(o2::mft::constants::mft::LayerZCoordinate()[9]), "Z position of the matching plane"}; + std::array*, Chi2FunctionsNum> matchingPlaneZs{ + &cfgChi2FunctionMatchingPlaneZ1, &cfgChi2FunctionMatchingPlaneZ2, &cfgChi2FunctionMatchingPlaneZ3, &cfgChi2FunctionMatchingPlaneZ4, &cfgChi2FunctionMatchingPlaneZ5}; + + Configurable cfgChi2MatchingExtrapMethod1{"cfgChi2MatchingExtrapMethod1", 0, "Method for MCH track extrapolation to maching plane"}; + Configurable cfgChi2MatchingExtrapMethod2{"cfgChi2MatchingExtrapMethod2", 0, "Method for MCH track extrapolation to maching plane"}; + Configurable cfgChi2MatchingExtrapMethod3{"cfgChi2MatchingExtrapMethod3", 0, "Method for MCH track extrapolation to maching plane"}; + Configurable cfgChi2MatchingExtrapMethod4{"cfgChi2MatchingExtrapMethod4", 0, "Method for MCH track extrapolation to maching plane"}; + Configurable cfgChi2MatchingExtrapMethod5{"cfgChi2MatchingExtrapMethod5", 0, "Method for MCH track extrapolation to maching plane"}; + std::array*, Chi2FunctionsNum> matchingExtrapMethods{ + &cfgChi2MatchingExtrapMethod1, &cfgChi2MatchingExtrapMethod2, &cfgChi2MatchingExtrapMethod3, &cfgChi2MatchingExtrapMethod4, &cfgChi2MatchingExtrapMethod5}; + } configChi2MatchingOptions; // ML interface - static constexpr int sMLModelsNum = 5; + static constexpr int MlModelsNum = 5; struct : ConfigurableGroup { - Configurable fModelLabel_1{"cfgMLModelLabel_1", std::string{""}, "Text label identifying this group of ML models"}; - Configurable fModelLabel_2{"cfgMLModelLabel_2", std::string{""}, "Text label identifying this group of ML models"}; - Configurable fModelLabel_3{"cfgMLModelLabel_3", std::string{""}, "Text label identifying this group of ML models"}; - Configurable fModelLabel_4{"cfgMLModelLabel_4", std::string{""}, "Text label identifying this group of ML models"}; - Configurable fModelLabel_5{"cfgMLModelLabel_5", std::string{""}, "Text label identifying this group of ML models"}; - std::array*, sMLModelsNum> fModelLabel{ - &fModelLabel_1, &fModelLabel_2, &fModelLabel_3, &fModelLabel_4, &fModelLabel_5}; - - Configurable fModelPathCCDB_1{"cfgMLModelPathCCDB_1", "Users/m/mcoquet/MLTest", "Paths of models on CCDB"}; - Configurable fModelPathCCDB_2{"cfgMLModelPathsCCDB_2", std::string{""}, "Paths of models on CCDB"}; - Configurable fModelPathCCDB_3{"cfgMLModelPathsCCDB_3", std::string{""}, "Paths of models on CCDB"}; - Configurable fModelPathCCDB_4{"cfgMLModelPathsCCDB_4", std::string{""}, "Paths of models on CCDB"}; - Configurable fModelPathCCDB_5{"cfgMLModelPathsCCDB_5", std::string{""}, "Paths of models on CCDB"}; - std::array*, sMLModelsNum> fModelPathCCDB{ - &fModelPathCCDB_1, &fModelPathCCDB_2, &fModelPathCCDB_3, &fModelPathCCDB_4, &fModelPathCCDB_5}; - - Configurable fModelName_1{"cfgMLModelName_1", "model.onnx", "ONNX file names for each pT bin (if not from CCDB full path)"}; - Configurable fModelName_2{"cfgMLModelNames_2", std::string{""}, "ONNX file names for each pT bin (if not from CCDB full path)"}; - Configurable fModelName_3{"cfgMLModelNames_3", std::string{""}, "ONNX file names for each pT bin (if not from CCDB full path)"}; - Configurable fModelName_4{"cfgMLModelNames_4", std::string{""}, "ONNX file names for each pT bin (if not from CCDB full path)"}; - Configurable fModelName_5{"cfgMLModelNames_5", std::string{""}, "ONNX file names for each pT bin (if not from CCDB full path)"}; - std::array*, sMLModelsNum> fModelName{ - &fModelName_1, &fModelName_2, &fModelName_3, &fModelName_4, &fModelName_5}; - - Configurable fInputFeatures_1{"cfgMLInputFeatures_1", "chi2MCHMFT", "Names of ML model input features"}; - Configurable fInputFeatures_2{"cfgMLInputFeatures_2", std::string{""}, "Names of ML model input features"}; - Configurable fInputFeatures_3{"cfgMLInputFeatures_3", std::string{""}, "Names of ML model input features"}; - Configurable fInputFeatures_4{"cfgMLInputFeatures_4", std::string{""}, "Names of ML model input features"}; - Configurable fInputFeatures_5{"cfgMLInputFeatures_5", std::string{""}, "Names of ML model input features"}; - std::array*, sMLModelsNum> fInputFeatures{ - &fInputFeatures_1, &fInputFeatures_2, &fInputFeatures_3, &fInputFeatures_4, &fInputFeatures_5}; - - Configurable fMatchingScoreCut_1{"cfgMLModelMatchingScoreCut_1", 0.f, "Minimum score value for selecting good matches"}; - Configurable fMatchingScoreCut_2{"cfgMLModelMatchingScoreCut_2", 0.f, "Minimum score value for selecting good matches"}; - Configurable fMatchingScoreCut_3{"cfgMLModelMatchingScoreCut_3", 0.f, "Minimum score value for selecting good matches"}; - Configurable fMatchingScoreCut_4{"cfgMLModelMatchingScoreCut_4", 0.f, "Minimum score value for selecting good matches"}; - Configurable fMatchingScoreCut_5{"cfgMLModelMatchingScoreCut_5", 0.f, "Minimum score value for selecting good matches"}; - std::array*, sMLModelsNum> fMatchingScoreCut{ - &fMatchingScoreCut_1, &fMatchingScoreCut_2, &fMatchingScoreCut_3, &fMatchingScoreCut_4, &fMatchingScoreCut_5}; - - Configurable fMatchingPlaneZ_1{"cfgMLModelMatchingPlaneZ_1", static_cast(o2::mft::constants::mft::LayerZCoordinate()[9]), "Z position of the matching plane"}; - Configurable fMatchingPlaneZ_2{"cfgMLModelMatchingPlaneZ_2", 0.f, "Z position of the matching plane"}; - Configurable fMatchingPlaneZ_3{"cfgMLModelMatchingPlaneZ_3", 0.f, "Z position of the matching plane"}; - Configurable fMatchingPlaneZ_4{"cfgMLModelMatchingPlaneZ_4", 0.f, "Z position of the matching plane"}; - Configurable fMatchingPlaneZ_5{"cfgMLModelMatchingPlaneZ_5", 0.f, "Z position of the matching plane"}; - std::array*, sMLModelsNum> fMatchingPlaneZ{ - &fMatchingPlaneZ_1, &fMatchingPlaneZ_2, &fMatchingPlaneZ_3, &fMatchingPlaneZ_4, &fMatchingPlaneZ_5}; - - Configurable fMatchingExtrapMethod_1{"cfgMLMatchingExtrapMethod_1", static_cast(0), "Method for MCH track extrapolation to maching plane"}; - Configurable fMatchingExtrapMethod_2{"cfgMLMatchingExtrapMethod_2", static_cast(0), "Method for MCH track extrapolation to maching plane"}; - Configurable fMatchingExtrapMethod_3{"cfgMLMatchingExtrapMethod_3", static_cast(0), "Method for MCH track extrapolation to maching plane"}; - Configurable fMatchingExtrapMethod_4{"cfgMLMatchingExtrapMethod_4", static_cast(0), "Method for MCH track extrapolation to maching plane"}; - Configurable fMatchingExtrapMethod_5{"cfgMLMatchingExtrapMethod_5", static_cast(0), "Method for MCH track extrapolation to maching plane"}; - std::array*, sMLModelsNum> fMatchingExtrapMethod{ - &fMatchingExtrapMethod_1, &fMatchingExtrapMethod_2, &fMatchingExtrapMethod_3, &fMatchingExtrapMethod_4, &fMatchingExtrapMethod_5}; - } fConfigMlOptions; + Configurable cfgMlModelLabel1{"cfgMlModelLabel1", std::string{""}, "Text label identifying this group of ML models"}; + Configurable cfgMlModelLabel2{"cfgMlModelLabel2", std::string{""}, "Text label identifying this group of ML models"}; + Configurable cfgMlModelLabel3{"cfgMlModelLabel3", std::string{""}, "Text label identifying this group of ML models"}; + Configurable cfgMlModelLabel4{"cfgMlModelLabel4", std::string{""}, "Text label identifying this group of ML models"}; + Configurable cfgMlModelLabel5{"cfgMlModelLabel5", std::string{""}, "Text label identifying this group of ML models"}; + std::array*, MlModelsNum> modelLabels{ + &cfgMlModelLabel1, &cfgMlModelLabel2, &cfgMlModelLabel3, &cfgMlModelLabel4, &cfgMlModelLabel5}; + + Configurable cfgMlModelPathCcdb1{"cfgMlModelPathCcdb1", "Users/m/mcoquet/MLTest", "Paths of models on CCDB"}; + Configurable cfgMlModelPathCcdb2{"cfgMlModelPathCcdb2", std::string{""}, "Paths of models on CCDB"}; + Configurable cfgMlModelPathCcdb3{"cfgMlModelPathCcdb3", std::string{""}, "Paths of models on CCDB"}; + Configurable cfgMlModelPathCcdb4{"cfgMlModelPathCcdb4", std::string{""}, "Paths of models on CCDB"}; + Configurable cfgMlModelPathCcdb5{"cfgMlModelPathCcdb5", std::string{""}, "Paths of models on CCDB"}; + std::array*, MlModelsNum> modelPathCcds{ + &cfgMlModelPathCcdb1, &cfgMlModelPathCcdb2, &cfgMlModelPathCcdb3, &cfgMlModelPathCcdb4, &cfgMlModelPathCcdb5}; + + Configurable cfgMlModelName1{"cfgMlModelName1", "model.onnx", "ONNX file names for each pT bin (if not from CCDB full path)"}; + Configurable cfgMlModelName2{"cfgMlModelName2", std::string{""}, "ONNX file names for each pT bin (if not from CCDB full path)"}; + Configurable cfgMlModelName3{"cfgMlModelName3", std::string{""}, "ONNX file names for each pT bin (if not from CCDB full path)"}; + Configurable cfgMlModelName4{"cfgMlModelName4", std::string{""}, "ONNX file names for each pT bin (if not from CCDB full path)"}; + Configurable cfgMlModelName5{"cfgMlModelName5", std::string{""}, "ONNX file names for each pT bin (if not from CCDB full path)"}; + std::array*, MlModelsNum> modelNames{ + &cfgMlModelName1, &cfgMlModelName2, &cfgMlModelName3, &cfgMlModelName4, &cfgMlModelName5}; + + Configurable cfgMlInputFeatures1{"cfgMlInputFeatures1", "chi2MCHMFT", "Names of ML model input features"}; + Configurable cfgMlInputFeatures2{"cfgMlInputFeatures2", std::string{""}, "Names of ML model input features"}; + Configurable cfgMlInputFeatures3{"cfgMlInputFeatures3", std::string{""}, "Names of ML model input features"}; + Configurable cfgMlInputFeatures4{"cfgMlInputFeatures4", std::string{""}, "Names of ML model input features"}; + Configurable cfgMlInputFeatures5{"cfgMlInputFeatures5", std::string{""}, "Names of ML model input features"}; + std::array*, MlModelsNum> inputFeatures{ + &cfgMlInputFeatures1, &cfgMlInputFeatures2, &cfgMlInputFeatures3, &cfgMlInputFeatures4, &cfgMlInputFeatures5}; + + Configurable cfgMlModelMatchingScoreCut1{"cfgMlModelMatchingScoreCut1", 0.f, "Minimum score value for selecting good matches"}; + Configurable cfgMlModelMatchingScoreCut2{"cfgMlModelMatchingScoreCut2", 0.f, "Minimum score value for selecting good matches"}; + Configurable cfgMlModelMatchingScoreCut3{"cfgMlModelMatchingScoreCut3", 0.f, "Minimum score value for selecting good matches"}; + Configurable cfgMlModelMatchingScoreCut4{"cfgMlModelMatchingScoreCut4", 0.f, "Minimum score value for selecting good matches"}; + Configurable cfgMlModelMatchingScoreCut5{"cfgMlModelMatchingScoreCut5", 0.f, "Minimum score value for selecting good matches"}; + std::array*, MlModelsNum> matchingScoreCuts{ + &cfgMlModelMatchingScoreCut1, &cfgMlModelMatchingScoreCut2, &cfgMlModelMatchingScoreCut3, &cfgMlModelMatchingScoreCut4, &cfgMlModelMatchingScoreCut5}; + + Configurable cfgMlModelMatchingPlaneZ1{"cfgMlModelMatchingPlaneZ1", static_cast(o2::mft::constants::mft::LayerZCoordinate()[9]), "Z position of the matching plane"}; + Configurable cfgMlModelMatchingPlaneZ2{"cfgMlModelMatchingPlaneZ2", 0.f, "Z position of the matching plane"}; + Configurable cfgMlModelMatchingPlaneZ3{"cfgMlModelMatchingPlaneZ3", 0.f, "Z position of the matching plane"}; + Configurable cfgMlModelMatchingPlaneZ4{"cfgMlModelMatchingPlaneZ4", 0.f, "Z position of the matching plane"}; + Configurable cfgMlModelMatchingPlaneZ5{"cfgMlModelMatchingPlaneZ5", 0.f, "Z position of the matching plane"}; + std::array*, MlModelsNum> matchingPlaneZs{ + &cfgMlModelMatchingPlaneZ1, &cfgMlModelMatchingPlaneZ2, &cfgMlModelMatchingPlaneZ3, &cfgMlModelMatchingPlaneZ4, &cfgMlModelMatchingPlaneZ5}; + + Configurable cfgMlMatchingExtrapMethod1{"cfgMlMatchingExtrapMethod1", 0, "Method for MCH track extrapolation to maching plane"}; + Configurable cfgMlMatchingExtrapMethod2{"cfgMlMatchingExtrapMethod2", 0, "Method for MCH track extrapolation to maching plane"}; + Configurable cfgMlMatchingExtrapMethod3{"cfgMlMatchingExtrapMethod3", 0, "Method for MCH track extrapolation to maching plane"}; + Configurable cfgMlMatchingExtrapMethod4{"cfgMlMatchingExtrapMethod4", 0, "Method for MCH track extrapolation to maching plane"}; + Configurable cfgMlMatchingExtrapMethod5{"cfgMlMatchingExtrapMethod5", 0, "Method for MCH track extrapolation to maching plane"}; + std::array*, MlModelsNum> matchingExtrapMethods{ + &cfgMlMatchingExtrapMethod1, &cfgMlMatchingExtrapMethod2, &cfgMlMatchingExtrapMethod3, &cfgMlMatchingExtrapMethod4, &cfgMlMatchingExtrapMethod5}; + } configMlOptions; std::vector binsPtMl; std::array cutValues; @@ -382,25 +469,29 @@ struct qaMatching { HistogramRegistry registryDimuon{"registryDimuon", {}}; std::unordered_map matchingHistos; - matrix dimuonHistos; + Matrix dimuonHistos; + + Produces qaMatchingEvents; + Produces qaMatchingMCHTrack; + Produces qaMatchingCandidates; struct EfficiencyPlotter { - o2::framework::HistPtr p_num; - o2::framework::HistPtr p_den; - o2::framework::HistPtr p_pdg_num; - o2::framework::HistPtr p_pdg_den; - o2::framework::HistPtr pt_num; - o2::framework::HistPtr pt_den; - o2::framework::HistPtr pt_pdg_num; - o2::framework::HistPtr pt_pdg_den; - o2::framework::HistPtr phi_num; - o2::framework::HistPtr phi_den; - o2::framework::HistPtr phi_pdg_num; - o2::framework::HistPtr phi_pdg_den; - o2::framework::HistPtr eta_num; - o2::framework::HistPtr eta_den; - o2::framework::HistPtr eta_pdg_num; - o2::framework::HistPtr eta_pdg_den; + o2::framework::HistPtr pNum; + o2::framework::HistPtr pDen; + o2::framework::HistPtr pPdgNum; + o2::framework::HistPtr pPdgDen; + o2::framework::HistPtr ptNum; + o2::framework::HistPtr ptDen; + o2::framework::HistPtr ptPdgNum; + o2::framework::HistPtr ptPdgDen; + o2::framework::HistPtr phiNum; + o2::framework::HistPtr phiDen; + o2::framework::HistPtr phiPdgNum; + o2::framework::HistPtr phiPdgDen; + o2::framework::HistPtr etaNum; + o2::framework::HistPtr etaDen; + o2::framework::HistPtr etaPdgNum; + o2::framework::HistPtr etaPdgDen; EfficiencyPlotter(std::string path, std::string title, HistogramRegistry& registry, bool createPdgMomHistograms) @@ -417,112 +508,114 @@ struct qaMatching { // momentum dependence histName = path + "p_num"; histTitle = title + " vs. p - num"; - p_num = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {pAxis}}); + pNum = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {pAxis}}); histName = path + "p_den"; histTitle = title + " vs. p - den"; - p_den = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {pAxis}}); + pDen = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {pAxis}}); if (createPdgMomHistograms) { histName = path + "p_pdg_num"; histTitle = title + " vs. p vs pdg ID - num"; - p_pdg_num = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {pAxis, motherPDGAxis}}); + pPdgNum = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {pAxis, motherPDGAxis}}); histName = path + "p_pdg_den"; histTitle = title + " vs. p vs pdg ID - den"; - p_pdg_den = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {pAxis, motherPDGAxis}}); + pPdgDen = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {pAxis, motherPDGAxis}}); } // pT dependence histName = path + "pt_num"; histTitle = title + " vs. p_{T} - num"; - pt_num = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {pTAxis}}); + ptNum = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {pTAxis}}); histName = path + "pt_den"; histTitle = title + " vs. p_{T} - den"; - pt_den = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {pTAxis}}); + ptDen = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {pTAxis}}); if (createPdgMomHistograms) { histName = path + "pt_pdg_num"; histTitle = title + " vs. p_{T} vs pdg ID - num"; - pt_pdg_num = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {pTAxis, motherPDGAxis}}); + ptPdgNum = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {pTAxis, motherPDGAxis}}); histName = path + "pt_pdg_den"; histTitle = title + " vs. p_{T} vs pdg ID - den"; - pt_pdg_den = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {pTAxis, motherPDGAxis}}); + ptPdgDen = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {pTAxis, motherPDGAxis}}); } // eta dependence histName = path + "eta_num"; histTitle = title + " vs. #eta - num"; - eta_num = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {etaAxis}}); + etaNum = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {etaAxis}}); histName = path + "eta_den"; histTitle = title + " vs. #eta - den"; - eta_den = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {etaAxis}}); + etaDen = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {etaAxis}}); if (createPdgMomHistograms) { histName = path + "eta_pdg_num"; histTitle = title + " vs. #eta vs pdg ID - num"; - eta_pdg_num = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {etaAxis, motherPDGAxis}}); + etaPdgNum = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {etaAxis, motherPDGAxis}}); histName = path + "eta_pdg_den"; histTitle = title + " vs. #eta vs pdg ID - den"; - eta_pdg_den = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {etaAxis, motherPDGAxis}}); + etaPdgDen = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {etaAxis, motherPDGAxis}}); } // phi dependence histName = path + "phi_num"; histTitle = title + " vs. #phi - num"; - phi_num = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {phiAxis}}); + phiNum = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {phiAxis}}); histName = path + "phi_den"; histTitle = title + " vs. #phi - den"; - phi_den = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {phiAxis}}); + phiDen = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH1F, {phiAxis}}); if (createPdgMomHistograms) { histName = path + "phi_pdg_num"; histTitle = title + " vs. #phi vs pdg ID - num"; - phi_pdg_num = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {phiAxis, motherPDGAxis}}); + phiPdgNum = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {phiAxis, motherPDGAxis}}); histName = path + "phi_pdg_den"; histTitle = title + " vs. #phi vs pdg ID - den"; - phi_pdg_den = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {phiAxis, motherPDGAxis}}); + phiPdgDen = registry.add(histName.c_str(), histTitle.c_str(), {HistType::kTH2F, {phiAxis, motherPDGAxis}}); } } template - void Fill(const T& track, bool passed) + void fill(const T& track, bool passed) { - double phi = track.phi() * 180 / TMath::Pi(); - std::get>(p_den)->Fill(track.p()); - std::get>(pt_den)->Fill(track.pt()); - std::get>(eta_den)->Fill(track.eta()); - std::get>(phi_den)->Fill(phi); + constexpr double RadToDeg = 180. / o2::constants::math::PI; + double phi = track.phi() * RadToDeg; + std::get>(pDen)->Fill(track.p()); + std::get>(ptDen)->Fill(track.pt()); + std::get>(etaDen)->Fill(track.eta()); + std::get>(phiDen)->Fill(phi); if (passed) { - std::get>(p_num)->Fill(track.p()); - std::get>(pt_num)->Fill(track.pt()); - std::get>(eta_num)->Fill(track.eta()); - std::get>(phi_num)->Fill(phi); + std::get>(pNum)->Fill(track.p()); + std::get>(ptNum)->Fill(track.pt()); + std::get>(etaNum)->Fill(track.eta()); + std::get>(phiNum)->Fill(phi); } } // Study the PDG origin of particles and their effect on the purity score template - void Fill(const T& track, int pdgCode, bool passed) + void fill(const T& track, int pdgCode, bool passed) { - double phi = track.phi() * 180 / TMath::Pi(); - std::get>(p_pdg_den)->Fill(track.p(), pdgCode); - std::get>(pt_pdg_den)->Fill(track.pt(), pdgCode); - std::get>(eta_pdg_den)->Fill(track.eta(), pdgCode); - std::get>(phi_pdg_den)->Fill(phi, pdgCode); + constexpr double RadToDeg = 180. / o2::constants::math::PI; + double phi = track.phi() * RadToDeg; + std::get>(pPdgDen)->Fill(track.p(), pdgCode); + std::get>(ptPdgDen)->Fill(track.pt(), pdgCode); + std::get>(etaPdgDen)->Fill(track.eta(), pdgCode); + std::get>(phiPdgDen)->Fill(phi, pdgCode); if (passed) { - std::get>(p_pdg_num)->Fill(track.p(), pdgCode); - std::get>(pt_pdg_num)->Fill(track.pt(), pdgCode); - std::get>(eta_pdg_num)->Fill(track.eta(), pdgCode); - std::get>(phi_pdg_num)->Fill(phi, pdgCode); + std::get>(pPdgNum)->Fill(track.p(), pdgCode); + std::get>(ptPdgNum)->Fill(track.pt(), pdgCode); + std::get>(etaPdgNum)->Fill(track.eta(), pdgCode); + std::get>(phiPdgNum)->Fill(phi, pdgCode); } } }; @@ -808,7 +901,7 @@ struct qaMatching { CollisionInfos fCollisionInfos; template - void initCCDB(BC const& bc) + void initCcdb(BC const& bc) { if (mRunNumber == bc.runNumber()) return; @@ -817,7 +910,7 @@ struct qaMatching { std::map metadata; auto soreor = o2::ccdb::BasicCCDBManager::getRunDuration(fCCDBApi, mRunNumber); auto ts = soreor.first; - auto grpmag = fCCDBApi.retrieveFromTFileAny(grpmagPath, metadata, ts); + auto grpmag = fCCDBApi.retrieveFromTFileAny(grpMagPath, metadata, ts); o2::base::Propagator::initFieldFromGRP(grpmag); LOGF(info, "Set field for muons"); VarManager::SetupMuonMagField(); @@ -827,13 +920,13 @@ struct qaMatching { 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); + double centerMft[3] = {0, 0, -61.4}; // Field at center of MFT + mBzAtMftCenter = fieldB->getBz(centerMft); // std::cout << "fieldB: " << (void*)fieldB << std::endl; } } - void CreateMatchingHistosMC() + void createMatchingHistosMc() { AxisSpec chi2Axis = {1000, 0, 1000, "chi^{2}"}; AxisSpec chi2AxisSmall = {200, 0, 100, "chi^{2}"}; @@ -843,30 +936,30 @@ struct qaMatching { AxisSpec phiAxis = {90, -180, 180, "#phi (degrees)"}; std::string histPath = "matching/MC/"; - AxisSpec trackPositionXAtMFTAxis = {100, -15, 15, "MFT x (cm)"}; - AxisSpec trackPositionYAtMFTAxis = {100, -15, 15, "MFT y (cm)"}; - registry.add((histPath + "pairedMCHTracksAtMFT").c_str(), "Paired MCH tracks position at MFT end", {HistType::kTH2F, {trackPositionXAtMFTAxis, trackPositionYAtMFTAxis}}); - registry.add((histPath + "pairedMFTTracksAtMFT").c_str(), "Paired MFT tracks position at MFT end", {HistType::kTH2F, {trackPositionXAtMFTAxis, trackPositionYAtMFTAxis}}); - registry.add((histPath + "selectedMCHTracksAtMFT").c_str(), "Selected MCH tracks position at MFT end", {HistType::kTH2F, {trackPositionXAtMFTAxis, trackPositionYAtMFTAxis}}); - registry.add((histPath + "selectedMCHTracksAtMFTTrue").c_str(), "Selected MCH tracks position at MFT end - true", {HistType::kTH2F, {trackPositionXAtMFTAxis, trackPositionYAtMFTAxis}}); - registry.add((histPath + "selectedMCHTracksAtMFTFake").c_str(), "Selected MCH tracks position at MFT end - fake", {HistType::kTH2F, {trackPositionXAtMFTAxis, trackPositionYAtMFTAxis}}); + AxisSpec trackPositionXAtMftAxis = {100, -15, 15, "MFT x (cm)"}; + AxisSpec trackPositionYAtMftAxis = {100, -15, 15, "MFT y (cm)"}; + registry.add((histPath + "pairedMCHTracksAtMFT").c_str(), "Paired MCH tracks position at MFT end", {HistType::kTH2F, {trackPositionXAtMftAxis, trackPositionYAtMftAxis}}); + registry.add((histPath + "pairedMFTTracksAtMFT").c_str(), "Paired MFT tracks position at MFT end", {HistType::kTH2F, {trackPositionXAtMftAxis, trackPositionYAtMftAxis}}); + registry.add((histPath + "selectedMCHTracksAtMFT").c_str(), "Selected MCH tracks position at MFT end", {HistType::kTH2F, {trackPositionXAtMftAxis, trackPositionYAtMftAxis}}); + registry.add((histPath + "selectedMCHTracksAtMFTTrue").c_str(), "Selected MCH tracks position at MFT end - true", {HistType::kTH2F, {trackPositionXAtMftAxis, trackPositionYAtMftAxis}}); + registry.add((histPath + "selectedMCHTracksAtMFTFake").c_str(), "Selected MCH tracks position at MFT end - fake", {HistType::kTH2F, {trackPositionXAtMftAxis, trackPositionYAtMftAxis}}); - fChi2MatchingPlotter = std::make_unique(histPath + "Prod/", ®istryMatching, fConfigQAs.fCreatePdgMomHistograms, fMftTrackMultiplicityMax, fNCandidatesMax); + fChi2MatchingPlotter = std::make_unique(histPath + "Prod/", ®istryMatching, configQas.cfgCreatePdgMomHistograms, cfgMftTrackMultiplicityMax, cfgNCandidatesMax); int registryIndex = 0; for (const auto& [label, func] : matchingChi2Functions) { - fMatchingPlotters[label] = std::make_unique(histPath + label + "/", registryMatchingVec[registryIndex], fConfigQAs.fCreatePdgMomHistograms, fMftTrackMultiplicityMax, fNCandidatesMax); + fMatchingPlotters[label] = std::make_unique(histPath + label + "/", registryMatchingVec[registryIndex], configQas.cfgCreatePdgMomHistograms, cfgMftTrackMultiplicityMax, cfgNCandidatesMax); registryIndex += 1; } for (const auto& [label, response] : matchingMlResponses) { - fMatchingPlotters[label] = std::make_unique(histPath + label + "/", (registryMatchingVec[registryIndex]), fConfigQAs.fCreatePdgMomHistograms, fMftTrackMultiplicityMax, fNCandidatesMax); + fMatchingPlotters[label] = std::make_unique(histPath + label + "/", (registryMatchingVec[registryIndex]), configQas.cfgCreatePdgMomHistograms, cfgMftTrackMultiplicityMax, cfgNCandidatesMax); registryIndex += 1; } - fTaggedMuonsMatchingPlotter = std::make_unique(histPath + "Tagged/", ®istryMatching, fConfigQAs.fCreatePdgMomHistograms, fMftTrackMultiplicityMax, fNCandidatesMax); - fSelectedMuonsMatchingPlotter = std::make_unique(histPath + "Selected/", ®istryMatching, fConfigQAs.fCreatePdgMomHistograms, fMftTrackMultiplicityMax, fNCandidatesMax); + fTaggedMuonsMatchingPlotter = std::make_unique(histPath + "Tagged/", ®istryMatching, configQas.cfgCreatePdgMomHistograms, cfgMftTrackMultiplicityMax, cfgNCandidatesMax); + fSelectedMuonsMatchingPlotter = std::make_unique(histPath + "Selected/", ®istryMatching, configQas.cfgCreatePdgMomHistograms, cfgMftTrackMultiplicityMax, cfgNCandidatesMax); } - void CreateDimuonHistos() + void createDimuonHistos() { AxisSpec invMassAxis = {400, 1, 5, "M_{#mu^{+}#mu^{-}} (GeV/c^{2})"}; AxisSpec invMassCorrelationAxis = {400, 0, 8, "M_{#mu^{+}#mu^{-}} (GeV/c^{2})"}; @@ -896,7 +989,7 @@ struct qaMatching { registryDimuon.add("dimuon/MC/invariantMass_ScaledMftKine_GlobalMuonCuts_GoodMatches_vs_match_type", "#mu^{+}#mu^{-} invariant mass vs. match tye (global muon cuts, rescaled MFT, good matches)", {HistType::kTH2F, {invMassAxis, matchTypeAxis}}); } - void InitMatchingFunctions() + void initMatchingFunctions() { using SMatrix55Std = ROOT::Math::SMatrix; using SMatrix55Sym = ROOT::Math::SMatrix>; @@ -915,31 +1008,31 @@ struct qaMatching { mMatchingFunctionMap["matchALL"] = [](const o2::dataformats::GlobalFwdTrack& mchTrack, const o2::track::TrackParCovFwd& mftTrack) -> std::tuple { // Match two tracks evaluating all parameters: X,Y, phi, tanl & q/pt - SMatrix55Sym 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; + SMatrix55Sym hK, vK; + SVector5 mK(mftTrack.getX(), mftTrack.getY(), mftTrack.getPhi(), + mftTrack.getTanl(), mftTrack.getInvQPt()), + rKKminus1; + SVector5 globalMuonTrackParameters = mchTrack.getParameters(); + SMatrix55Sym 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 = (V_k + ROOT::Math::Similarity(H_k, GlobalMuonTrackCovariances)); + SMatrix55Std invResCov = (vK + ROOT::Math::Similarity(hK, globalMuonTrackCovariances)); invResCov.Invert(); // Update Parameters - r_k_kminus1 = m_k - H_k * GlobalMuonTrackParameters; // Residuals of prediction + rKKminus1 = mK - hK * globalMuonTrackParameters; // Residuals of prediction - auto matchChi2Track = ROOT::Math::Similarity(r_k_kminus1, invResCov); + auto matchChi2Track = ROOT::Math::Similarity(rKKminus1, invResCov); // return chi2 and NDF return {matchChi2Track, 5}; @@ -949,30 +1042,30 @@ struct qaMatching { mMatchingFunctionMap["matchXYPhiTanl"] = [](const o2::dataformats::GlobalFwdTrack& mchTrack, const o2::track::TrackParCovFwd& mftTrack) -> std::tuple { // Match two tracks evaluating positions & angles - SMatrix45 H_k; - SMatrix44 V_k; - SVector4 m_k(mftTrack.getX(), mftTrack.getY(), mftTrack.getPhi(), - mftTrack.getTanl()), - 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); - H_k(0, 0) = 1.0; - H_k(1, 1) = 1.0; - H_k(2, 2) = 1.0; - H_k(3, 3) = 1.0; + SMatrix45 hK; + SMatrix44 vK; + SVector4 mK(mftTrack.getX(), mftTrack.getY(), mftTrack.getPhi(), + mftTrack.getTanl()), + rKKminus1; + SVector5 globalMuonTrackParameters = mchTrack.getParameters(); + SMatrix55Sym 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 = (V_k + ROOT::Math::Similarity(H_k, GlobalMuonTrackCovariances)); + SMatrix44 invResCov = (vK + ROOT::Math::Similarity(hK, globalMuonTrackCovariances)); invResCov.Invert(); // Residuals of prediction - r_k_kminus1 = m_k - H_k * GlobalMuonTrackParameters; + rKKminus1 = mK - hK * globalMuonTrackParameters; - auto matchChi2Track = ROOT::Math::Similarity(r_k_kminus1, invResCov); + auto matchChi2Track = ROOT::Math::Similarity(rKKminus1, invResCov); // return chi2 and NDF return {matchChi2Track, 4}; @@ -982,23 +1075,23 @@ struct qaMatching { mMatchingFunctionMap["matchXY"] = [](const o2::dataformats::GlobalFwdTrack& mchTrack, const o2::track::TrackParCovFwd& mftTrack) -> std::tuple { // Calculate Matching Chi2 - X and Y positions - 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; + SMatrix25 hK; + SMatrix22 vK; + SVector2 mK(mftTrack.getX(), mftTrack.getY()), rKKminus1; + SVector5 globalMuonTrackParameters = mchTrack.getParameters(); + SMatrix55Sym 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 = (V_k + ROOT::Math::Similarity(H_k, GlobalMuonTrackCovariances)); + SMatrix22 invResCov = (vK + ROOT::Math::Similarity(hK, globalMuonTrackCovariances)); invResCov.Invert(); // Residuals of prediction - r_k_kminus1 = m_k - H_k * GlobalMuonTrackParameters; - auto matchChi2Track = ROOT::Math::Similarity(r_k_kminus1, invResCov); + rKKminus1 = mK - hK * globalMuonTrackParameters; + auto matchChi2Track = ROOT::Math::Similarity(rKKminus1, invResCov); // return reduced chi2 return {matchChi2Track, 2}; @@ -1008,10 +1101,10 @@ struct qaMatching { void init(o2::framework::InitContext&) { // Load geometry - ccdbManager->setURL(ccdburl); + ccdbManager->setURL(ccdbUrl); ccdbManager->setCaching(true); ccdbManager->setLocalObjectValidityChecking(); - fCCDBApi.init(ccdburl); + fCCDBApi.init(ccdbUrl); mRunNumber = 0; if (!o2::base::GeometryManager::isGeometryLoaded()) { @@ -1020,13 +1113,13 @@ struct qaMatching { } // Matching functions - InitMatchingFunctions(); - for (size_t funcId = 0; funcId < sChi2FunctionsNum; funcId++) { - auto label = fConfigChi2MatchingOptions.fFunctionLabel[funcId]->value; - auto funcName = fConfigChi2MatchingOptions.fFunctionName[funcId]->value; - auto scoreMin = fConfigChi2MatchingOptions.fMatchingScoreCut[funcId]->value; - auto matchingPlaneZ = fConfigChi2MatchingOptions.fMatchingPlaneZ[funcId]->value; - auto extrapMethod = fConfigChi2MatchingOptions.fMatchingExtrapMethod[funcId]->value; + initMatchingFunctions(); + for (size_t funcId = 0; funcId < Chi2FunctionsNum; funcId++) { + auto label = configChi2MatchingOptions.functionLabels[funcId]->value; + auto funcName = configChi2MatchingOptions.functionNames[funcId]->value; + auto scoreMin = configChi2MatchingOptions.matchingScoreCuts[funcId]->value; + auto matchingPlaneZ = configChi2MatchingOptions.matchingPlaneZs[funcId]->value; + auto extrapMethod = configChi2MatchingOptions.matchingExtrapMethods[funcId]->value; if (label == "" || funcName == "") break; @@ -1045,20 +1138,20 @@ struct qaMatching { cutDirMl = {cuts_ml::CutNot}; o2::framework::LabeledArray mycutsMl(cutValues.data(), 1, 1, std::vector{"pT bin 0"}, std::vector{"score"}); - for (size_t modelId = 0; modelId < sMLModelsNum; modelId++) { - auto label = fConfigMlOptions.fModelLabel[modelId]->value; - auto modelPath = fConfigMlOptions.fModelPathCCDB[modelId]->value; - auto inputFeatures = fConfigMlOptions.fInputFeatures[modelId]->value; - auto modelName = fConfigMlOptions.fModelName[modelId]->value; - auto scoreMin = fConfigMlOptions.fMatchingScoreCut[modelId]->value; - auto matchingPlaneZ = fConfigMlOptions.fMatchingPlaneZ[modelId]->value; - auto extrapMethod = fConfigMlOptions.fMatchingExtrapMethod[modelId]->value; + for (size_t modelId = 0; modelId < MlModelsNum; modelId++) { + auto label = configMlOptions.modelLabels[modelId]->value; + auto modelPath = configMlOptions.modelPathCcds[modelId]->value; + auto inputFeatures = configMlOptions.inputFeatures[modelId]->value; + auto modelName = configMlOptions.modelNames[modelId]->value; + auto scoreMin = configMlOptions.matchingScoreCuts[modelId]->value; + auto matchingPlaneZ = configMlOptions.matchingPlaneZs[modelId]->value; + auto extrapMethod = configMlOptions.matchingExtrapMethods[modelId]->value; if (label == "" || modelPath == "" || inputFeatures == "" || modelName == "") break; matchingMlResponses[label].configure(binsPtMl, mycutsMl, cutDirMl, 1); - matchingMlResponses[label].setModelPathsCCDB(std::vector{modelName}, fCCDBApi, std::vector{modelPath}, fConfigCCDB.fConfigNoLaterThan.value); + matchingMlResponses[label].setModelPathsCCDB(std::vector{modelName}, fCCDBApi, std::vector{modelPath}, configCcdb.cfgCcdbNoLaterThan.value); matchingMlResponses[label].cacheInputFeaturesIndices(std::vector{inputFeatures}); matchingMlResponses[label].init(); @@ -1071,23 +1164,25 @@ struct qaMatching { AxisSpec trackTypeAxis = {static_cast(nTrackTypes), 0.0, static_cast(nTrackTypes), "track type"}; registry.add("nTracksPerType", "Number of tracks per type", {HistType::kTH1F, {trackTypeAxis}}); - AxisSpec tracksMultiplicityAxis = {fMftTrackMultiplicityMax, 0, static_cast(fMftTrackMultiplicityMax), "tracks multiplicity"}; + AxisSpec tracksMultiplicityAxis = {cfgMftTrackMultiplicityMax, 0, static_cast(cfgMftTrackMultiplicityMax), "tracks multiplicity"}; registry.add("tracksMultiplicityMFT", "MFT tracks multiplicity", {HistType::kTH1F, {tracksMultiplicityAxis}}); registry.add("tracksMultiplicityMCH", "MCH tracks multiplicity", {HistType::kTH1F, {tracksMultiplicityAxis}}); - CreateMatchingHistosMC(); - CreateDimuonHistos(); + createMatchingHistosMc(); + createDimuonHistos(); } template - bool pDCACut(const T& mchTrack, const C& collision, double nSigmaPDCA) + 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(); + 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); @@ -1096,11 +1191,11 @@ struct qaMatching { double p = mchTrackAtVertex.getP(); double pDCA = mchTrack.pDca(); - double sigmaPDCA = (thetaAbs < 3) ? sigmaPDCA23 : sigmaPDCA310; + double sigmaPDCA = (thetaAbs < ThetaAbsBoundaryDeg) ? 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); + double slopeResEffect = SlopeResolutionZ * slopeRes * p; + double sigmaPDCAWithRes = std::sqrt(pResEffect * pResEffect + slopeResEffect * slopeResEffect); if (pDCA > nSigmaPDCA * sigmaPDCAWithRes) { return false; } @@ -1109,7 +1204,7 @@ struct qaMatching { } template - bool IsGoodMuon(const T& mchTrack, const C& collision, + bool isGoodMuon(const T& mchTrack, const C& collision, double chi2Cut, double pCut, double pTCut, @@ -1144,7 +1239,7 @@ struct qaMatching { } // pDCA cut - if (!pDCACut(mchTrack, collision, nSigmaPdcaCut)) { + if (!pDcaCut(mchTrack, collision, nSigmaPdcaCut)) { return false; } @@ -1152,19 +1247,19 @@ struct qaMatching { } template - bool IsGoodMuon(const T& muonTrack, const C& collision) + bool isGoodMuon(const T& muonTrack, const C& collision) { - return IsGoodMuon(muonTrack, collision, fTrackChi2MchUp, fPMchLow, fPtMchLow, {fEtaMchLow, fEtaMchUp}, {fRabsLow, fRabsUp}, fSigmaPdcaUp); + return isGoodMuon(muonTrack, collision, cfgTrackChi2MchUp, cfgPMchLow, cfgPtMchLow, {cfgEtaMchLow, cfgEtaMchUp}, {cfgRabsLow, cfgRabsUp}, cfgPdcaUp); } template - bool IsGoodGlobalMuon(const T& muonTrack, const C& collision) + bool isGoodGlobalMuon(const T& muonTrack, const C& collision) { - return IsGoodMuon(muonTrack, collision, fTrackChi2MchUp, fPMchLow, fPtMchLow, {fEtaMftLow, fEtaMftUp}, {fRabsLow, fRabsUp}, fSigmaPdcaUp); + return isGoodMuon(muonTrack, collision, cfgTrackChi2MchUp, cfgPMchLow, cfgPtMchLow, {cfgEtaMftLow, cfgEtaMftUp}, {cfgRabsLow, cfgRabsUp}, cfgPdcaUp); } template - bool IsGoodMFT(const T& mftTrack, + bool isGoodMft(const T& mftTrack, double chi2Cut, int nClustersCut) { @@ -1183,17 +1278,17 @@ struct qaMatching { } template - bool IsGoodMFT(const T& mftTrack) + bool isGoodMft(const T& mftTrack) { - return IsGoodMFT(mftTrack, fTrackChi2MftUp, fTrackNClustMftLow); + return isGoodMft(mftTrack, cfgTrackChi2MftUp, cfgTrackNClustMftLow); } template - bool IsGoodGlobalMatching(const TMUON& muonTrack, + bool isGoodGlobalMatching(const TMUON& muonTrack, double matchingScore, double matchingScoreCut) { - if (static_cast(muonTrack.trackType()) > 2) + if (static_cast(muonTrack.trackType()) > GlobalTrackTypeMax) return false; // MFT-MCH matching score cut @@ -1204,15 +1299,15 @@ struct qaMatching { } template - bool IsGoodGlobalMatching(const TMUON& muonTrack, double matchingScore) + bool isGoodGlobalMatching(const TMUON& muonTrack, double matchingScore) { - return IsGoodGlobalMatching(muonTrack, matchingScore, fMatchingChi2ScoreMftMchLow); + return isGoodGlobalMatching(muonTrack, matchingScore, cfgMatchingChi2ScoreMftMchLow); } template - bool IsTrueGlobalMatching(const TMUON& muonTrack, const std::vector>& matchablePairs) + bool isTrueGlobalMatching(const TMUON& muonTrack, const std::vector>& matchablePairs) { - if (static_cast(muonTrack.trackType()) > 2) + if (static_cast(muonTrack.trackType()) > GlobalTrackTypeMax) return false; int64_t mchTrackId = static_cast(muonTrack.matchMCHTrackId()); @@ -1223,18 +1318,18 @@ struct qaMatching { return (std::find(matchablePairs.begin(), matchablePairs.end(), trackIndexes) != matchablePairs.end()); } - bool IsMatchableMCH(int64_t mchTrackId, const std::vector>& matchablePairs) + bool isMatchableMch(int64_t mchTrackId, const std::vector>& matchablePairs) { - for (auto [id1, id2] : matchablePairs) { + for (const auto& [id1, id2] : matchablePairs) { if (mchTrackId == id1) return true; } return false; } - std::optional> GetMatchablePairForMCH(int64_t mchTrackId, const std::vector>& matchablePairs) + std::optional> getMatchablePairForMch(int64_t mchTrackId, const std::vector>& matchablePairs) { - for (auto pair : matchablePairs) { + for (const auto& pair : matchablePairs) { if (mchTrackId == pair.first) return pair; } @@ -1242,7 +1337,7 @@ struct qaMatching { } template - o2::dataformats::GlobalFwdTrack FwdToTrackPar(const T& track) + o2::dataformats::GlobalFwdTrack fwdToTrackPar(const T& track) { double chi2 = track.chi2(); SMatrix5 tpars(track.x(), track.y(), track.phi(), track.tgl(), track.signed1Pt()); @@ -1259,7 +1354,7 @@ struct qaMatching { } template - o2::dataformats::GlobalFwdTrack FwdToTrackPar(const T& track, const C& cov) + o2::dataformats::GlobalFwdTrack fwdToTrackPar(const T& track, const C& cov) { double chi2 = track.chi2(); SMatrix5 tpars(track.x(), track.y(), track.phi(), track.tgl(), track.signed1Pt()); @@ -1275,7 +1370,7 @@ struct qaMatching { return fwdtrack; } - o2::dataformats::GlobalFwdTrack PropagateToZMCH(const o2::dataformats::GlobalFwdTrack& muon, const double z) + o2::dataformats::GlobalFwdTrack propagateToZMch(const o2::dataformats::GlobalFwdTrack& muon, const double z) { auto mchTrack = mExtrap.FwdtoMCH(muon); @@ -1300,7 +1395,7 @@ struct qaMatching { } template - o2::dataformats::GlobalFwdTrack PropagateToZMCH(const T& muon, const double z) + o2::dataformats::GlobalFwdTrack propagateToZMch(const T& muon, const double z) { double chi2 = muon.chi2(); SMatrix5 tpars(muon.x(), muon.y(), muon.phi(), muon.tgl(), muon.signed1Pt()); @@ -1314,10 +1409,10 @@ struct qaMatching { track.setZ(fwdtrack.getZ()); track.setCovariances(fwdtrack.getCovariances()); - return PropagateToZMCH(track, z); + return propagateToZMch(track, z); } - o2::dataformats::GlobalFwdTrack PropagateToZMFT(const o2::dataformats::GlobalFwdTrack& mftTrack, const double z) + o2::dataformats::GlobalFwdTrack propagateToZMft(const o2::dataformats::GlobalFwdTrack& mftTrack, const double z) { o2::dataformats::GlobalFwdTrack fwdtrack{mftTrack}; fwdtrack.propagateToZ(z, mBzAtMftCenter); @@ -1325,10 +1420,10 @@ struct qaMatching { } template - o2::dataformats::GlobalFwdTrack PropagateToZMFT(const TMFT& mftTrack, const CMFT& mftCov, const double z) + o2::dataformats::GlobalFwdTrack propagateToZMft(const TMFT& mftTrack, const CMFT& mftCov, const double z) { - o2::dataformats::GlobalFwdTrack fwdtrack = FwdToTrackPar(mftTrack, mftCov); - return PropagateToZMFT(fwdtrack, z); + o2::dataformats::GlobalFwdTrack fwdtrack = fwdToTrackPar(mftTrack, mftCov); + return propagateToZMft(fwdtrack, z); } // method 0: standard extrapolation @@ -1336,60 +1431,60 @@ struct qaMatching { // method 2: MCH track extrapolation constrained to the first MFT track point, MFT extrapolation using MCH momentum // method 3: MCH track extrapolation constrained to the collision point, MFT extrapolation using MCH momentum template - o2::dataformats::GlobalFwdTrack PropagateToMatchingPlaneMCH(const TMCH& mchTrack, const TMFT& mftTrack, const CMFT& mftTrackCov, const C& collision, const double z, int method) + o2::dataformats::GlobalFwdTrack propagateToMatchingPlaneMch(const TMCH& mchTrack, const TMFT& mftTrack, const CMFT& mftTrackCov, const C& collision, const double z, int method) { - if (method == 0 || method == 1) { + if (method == ExtrapolationMethodStandard || method == 1) { // simple extrapolation upstream through the absorber - return PropagateToZMCH(mchTrack, z); + return propagateToZMch(mchTrack, z); } - if (method == 2) { + if (method == ExtrapolationMethodMftFirstPoint) { // extrapolation to the first MFT point and then back to the matching plane - auto mftTrackPar = FwdToTrackPar(mftTrack, mftTrackCov); - // std::cout << std::format("[PropagateToMatchingPlaneMCH] extrapolating to MFT: x={:0.3f} y={:0.3f} z={:0.3f}", mftTrackPar.getX(), mftTrackPar.getY(), mftTrackPar.getZ()) << std::endl; - auto mchTrackAtMFT = PropagateToVertexMCH(FwdToTrackPar(mchTrack, mchTrack), + auto mftTrackPar = fwdToTrackPar(mftTrack, mftTrackCov); + // std::cout << std::format("[propagateToMatchingPlaneMch] extrapolating to MFT: x={:0.3f} y={:0.3f} z={:0.3f}", mftTrackPar.getX(), mftTrackPar.getY(), mftTrackPar.getZ()) << std::endl; + auto mchTrackAtMFT = propagateToVertexMch(fwdToTrackPar(mchTrack, mchTrack), mftTrackPar.getX(), mftTrackPar.getY(), mftTrackPar.getZ(), mftTrackPar.getSigma2X(), mftTrackPar.getSigma2Y()); - // std::cout << std::format("[PropagateToMatchingPlaneMCH] extrapolating to z={:0.3f}", z) << std::endl; - return PropagateToZMCH(mchTrackAtMFT, z); + // std::cout << std::format("[propagateToMatchingPlaneMch] extrapolating to z={:0.3f}", z) << std::endl; + return propagateToZMch(mchTrackAtMFT, z); } - if (method == 3) { + if (method == ExtrapolationMethodVertex) { // extrapolation to the vertex and then back to the matching plane auto mchTrackAtVertex = VarManager::PropagateMuon(mchTrack, collision, VarManager::kToVertex); - return PropagateToZMCH(mchTrackAtVertex, z); + return propagateToZMch(mchTrackAtVertex, z); } - if (method == 4) { + if (method == ExtrapolationMethodMftDca) { // extrapolation to the MFT DCA and then back to the matching plane - auto mftTrackDCA = PropagateToZMFT(FwdToTrackPar(mftTrack, mftTrackCov), collision.posZ()); - auto mchTrackAtDCA = PropagateToVertexMCH(FwdToTrackPar(mchTrack, mchTrack), + auto mftTrackDCA = propagateToZMft(fwdToTrackPar(mftTrack, mftTrackCov), collision.posZ()); + auto mchTrackAtDCA = propagateToVertexMch(fwdToTrackPar(mchTrack, mchTrack), mftTrackDCA.getX(), mftTrackDCA.getY(), mftTrackDCA.getZ(), mftTrackDCA.getSigma2X(), mftTrackDCA.getSigma2Y()); - return PropagateToZMCH(mchTrackAtDCA, z); + return propagateToZMch(mchTrackAtDCA, z); } return {}; } template - o2::dataformats::GlobalFwdTrack PropagateToMatchingPlaneMFT(const TMCH& mchTrack, const TMFT& mftTrack, const CMFT& mftTrackCov, const C& collision, const double z, int method) + o2::dataformats::GlobalFwdTrack propagateToMatchingPlaneMft(const TMCH& mchTrack, const TMFT& mftTrack, const CMFT& mftTrackCov, const C& collision, const double z, int method) { - if (method == 0) { + if (method == ExtrapolationMethodStandard) { // extrapolation with MFT tools - return PropagateToZMFT(mftTrack, mftTrackCov, z); + return propagateToZMft(mftTrack, mftTrackCov, z); } if (method > 0) { // extrapolation with MCH tools auto mchTrackAtVertex = VarManager::PropagateMuon(mchTrack, collision, VarManager::kToVertex); double pMCH = mchTrackAtVertex.getP(); - double px = pMCH * sin(M_PI / 2 - atan(mftTrack.tgl())) * cos(mftTrack.phi()); - double py = pMCH * sin(M_PI / 2 - atan(mftTrack.tgl())) * sin(mftTrack.phi()); - double pt = std::sqrt(std::pow(px, 2) + std::pow(py, 2)); + double px = pMCH * std::sin(o2::constants::math::PIHalf - std::atan(mftTrack.tgl())) * std::cos(mftTrack.phi()); + double py = pMCH * std::sin(o2::constants::math::PIHalf - std::atan(mftTrack.tgl())) * std::sin(mftTrack.phi()); + double pt = std::hypot(px, py); double sign = mchTrack.sign(); - o2::dataformats::GlobalFwdTrack track = FwdToTrackPar(mftTrack, mftTrackCov); + o2::dataformats::GlobalFwdTrack track = fwdToTrackPar(mftTrack, mftTrackCov); // update momentum in track parameters and errors auto newCov = track.getCovariances(); @@ -1412,7 +1507,7 @@ struct qaMatching { return {}; } - o2::dataformats::GlobalFwdTrack PropagateToVertexMCH(const o2::dataformats::GlobalFwdTrack& muon, + o2::dataformats::GlobalFwdTrack propagateToVertexMch(const o2::dataformats::GlobalFwdTrack& muon, const double vx, const double vy, const double vz, const double covVx, const double covVy) { @@ -1430,10 +1525,10 @@ struct qaMatching { } template - o2::dataformats::GlobalFwdTrack PropagateToVertexMCH(const TMCH& muon, + o2::dataformats::GlobalFwdTrack propagateToVertexMch(const TMCH& muon, const C& collision) { - return PropagateToVertexMCH(FwdToTrackPar(muon, muon), + return propagateToVertexMch(fwdToTrackPar(muon, muon), collision.posX(), collision.posY(), collision.posZ(), @@ -1441,7 +1536,7 @@ struct qaMatching { collision.covYY()); } - o2::dataformats::GlobalFwdTrack PropagateToVertexMFT(o2::dataformats::GlobalFwdTrack muon, + o2::dataformats::GlobalFwdTrack propagateToVertexMft(o2::dataformats::GlobalFwdTrack muon, const float vx, const float vy, const float vz, const float covVx, const float covVy) { @@ -1457,10 +1552,10 @@ struct qaMatching { } template - o2::dataformats::GlobalFwdTrack PropagateToVertexMFT(const TMFT& muon, + o2::dataformats::GlobalFwdTrack propagateToVertexMft(const TMFT& muon, const C& collision) { - return PropagateToVertexMFT(FwdToTrackPar(muon), + return propagateToVertexMft(fwdToTrackPar(muon), collision.posX(), collision.posY(), collision.posZ(), @@ -1469,15 +1564,15 @@ struct qaMatching { } template - o2::dataformats::GlobalFwdTrack PropagateToVertexMFT(const TMFT& muon, + o2::dataformats::GlobalFwdTrack propagateToVertexMft(const TMFT& muon, const TMCH& mchTrack, const C& collision) { // extrapolation with MCH tools - auto mchTrackAtMFT = mExtrap.FwdtoMCH(FwdToTrackPar(mchTrack)); + auto mchTrackAtMFT = mExtrap.FwdtoMCH(fwdToTrackPar(mchTrack)); o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchTrackAtMFT, muon.z()); - auto muonTrackProp = mExtrap.FwdtoMCH(FwdToTrackPar(muon)); + auto muonTrackProp = mExtrap.FwdtoMCH(fwdToTrackPar(muon)); // update global track momentum from the MCH track double pRatio = muonTrackProp.p() / mchTrackAtMFT.p(); @@ -1496,7 +1591,7 @@ struct qaMatching { } template - void GetMotherParticles(MCP const& mcParticle, std::vector>& motherParticlesVec) + void getMotherParticles(MCP const& mcParticle, std::vector>& motherParticlesVec) { const auto& motherParticles = mcParticle.template mothers_as(); if (motherParticles.empty()) { @@ -1505,11 +1600,11 @@ struct qaMatching { const auto& motherParticle = motherParticles[0]; motherParticlesVec.emplace_back(std::make_pair(static_cast(motherParticle.pdgCode()), static_cast(motherParticle.globalIndex()))); - GetMotherParticles(motherParticle, motherParticlesVec); + getMotherParticles(motherParticle, motherParticlesVec); } template - std::vector> GetMotherParticles(T const& track) + std::vector> getMotherParticles(T const& track) { std::vector> result; if (!track.has_mcParticle()) @@ -1518,15 +1613,15 @@ struct qaMatching { const auto& mcParticle = track.mcParticle(); result.emplace_back(std::make_pair(static_cast(mcParticle.pdgCode()), static_cast(mcParticle.globalIndex()))); - GetMotherParticles(mcParticle, result); + getMotherParticles(mcParticle, result); return result; } template - int GetDecayRanking(TMCH const& mchTrack, TMFTs const& mftTracks) + int getDecayRanking(TMCH const& mchTrack, TMFTs const& mftTracks) { - auto mchMotherParticles = GetMotherParticles(mchTrack); + auto mchMotherParticles = getMotherParticles(mchTrack); int decayRanking = 0; // search for an MFT track that is associated to one of the MCH mother particles @@ -1555,14 +1650,14 @@ struct qaMatching { } template - void FillMatchablePairs(CollisionInfo& collisionInfo, + void fillMatchablePairs(CollisionInfo& collisionInfo, TMUON const& muonTracks, TMFT const& mftTracks) { collisionInfo.matchablePairs.clear(); for (const auto& muonTrack : muonTracks) { // only consider MCH standalone or MCH-MID matches - if (static_cast(muonTrack.trackType()) <= 2) + if (static_cast(muonTrack.trackType()) <= GlobalTrackTypeMax) continue; // only consider tracks associated to the current collision @@ -1576,7 +1671,7 @@ struct qaMatching { continue; // get the index associated to the MC particle auto muonMcParticle = muonTrack.mcParticle(); - if (std::abs(muonMcParticle.pdgCode()) != 13) + if (std::abs(muonMcParticle.pdgCode()) != kMuonMinus) continue; int64_t muonMcTrackIndex = muonMcParticle.globalIndex(); @@ -1598,7 +1693,7 @@ struct qaMatching { } template - int GetTrueMatchIndex(TMUON const& muonTracks, + int getTrueMatchIndex(TMUON const& muonTracks, const std::vector& matchCandidatesVector, const std::vector>& matchablePairs) { @@ -1609,7 +1704,7 @@ struct qaMatching { for (size_t i = 0; i < matchCandidatesVector.size(); i++) { auto const& muonTrack = muonTracks.rawIteratorAt(matchCandidatesVector[i].globalTrackId); - if (IsTrueGlobalMatching(muonTrack, matchablePairs)) { + if (isTrueGlobalMatching(muonTrack, matchablePairs)) { trueMatchIndex = i + 1; break; } @@ -1617,92 +1712,43 @@ struct qaMatching { return trueMatchIndex; } - template - bool IsMuon(const TMCH& mchTrack, - const TMFT& mftTrack) - { - // skip tracks that do not have an associated MC particle - if (!mchTrack.has_mcParticle()) - return false; - if (!mftTrack.has_mcParticle()) - return false; - - // get the index associated to the MC particles - auto mchMcParticle = mchTrack.mcParticle(); - auto mftMcParticle = mftTrack.mcParticle(); - if (mchMcParticle.globalIndex() != mftMcParticle.globalIndex()) - return false; - - if (std::abs(mchMcParticle.pdgCode()) != 13) - return false; - - return true; - } - - template - bool IsMuon(const TMUON& muonTrack, - TMUONS const& /*muonTracks*/, - TMFTS const& /*mftTracks*/) - { - if (static_cast(muonTrack.trackType()) >= 2) - return false; - - auto const& mchTrack = muonTrack.template matchMCHTrack_as(); - auto const& mftTrack = muonTrack.template matchMFTTrack_as(); - - return IsMuon(mchTrack, mftTrack); - } - template - MuonMatchType GetMatchType(const TMUON& muonTrack, - TMUONS const& muonTracks, + MuonMatchType getMatchType(const TMUON& muonTrack, + TMUONS const& /*muonTracks*/, TMFTS const& mftTracks, const std::vector>& matchablePairs, int ranking) { - if (static_cast(muonTrack.trackType()) > 2) + if (static_cast(muonTrack.trackType()) > GlobalTrackTypeMax) return kMatchTypeUndefined; auto const& mchTrack = muonTrack.template matchMCHTrack_as(); - - bool isPaired = IsMatchableMCH(mchTrack.globalIndex(), matchablePairs); - bool isMuon = IsMuon(muonTrack, muonTracks, mftTracks); - int decayRanking = GetDecayRanking(mchTrack, mftTracks); + bool isPaired = isMatchableMch(mchTrack.globalIndex(), matchablePairs); + int decayRanking = getDecayRanking(mchTrack, mftTracks); MuonMatchType result{kMatchTypeUndefined}; if (isPaired) { - if (isMuon) { - result = (ranking == 1) ? kMatchTypeTrueLeading : kMatchTypeTrueNonLeading; - } else { - result = (ranking == 1) ? kMatchTypeWrongLeading : kMatchTypeWrongNonLeading; - } - } else if (decayRanking == 2) { + result = (ranking == 1) ? kMatchTypeTrueLeading : kMatchTypeTrueNonLeading; + } else if (decayRanking == FirstDecayMotherRank) { result = (ranking == 1) ? kMatchTypeDecayLeading : kMatchTypeDecayNonLeading; } else { result = (ranking == 1) ? kMatchTypeFakeLeading : kMatchTypeFakeNonLeading; } - - if (result == kMatchTypeUndefined) { - std::cout << std::format("[GetMatchType] isPaired={} isMuon={} decayRanking={} result={}", - isPaired, isMuon, decayRanking, static_cast(result)) - << std::endl; - } - return result; } // for each MCH standalone track, collect the associated matching candidates template - void GetSelectedMuons(const CollisionInfo& collisionInfo, + void getSelectedMuons(const CollisionInfo& collisionInfo, C const& collisions, TMUON const& muonTracks, std::vector& selectedMuons) { selectedMuons.clear(); - for (auto muonTrack : muonTracks) { + for (const auto& muonTrack : muonTracks) { // only consider MCH-MID matches - if (static_cast(muonTrack.trackType()) != 3) { + if (static_cast(muonTrack.trackType()) != MchMidTrackType) { continue; } @@ -1715,13 +1761,13 @@ struct qaMatching { const auto& collision = collisions.rawIteratorAt(muonTrack.collisionId()); // select MCH tracks with strict quality cuts - if (!IsGoodMuon(muonTrack, collision, - fMuonTaggingTrackChi2MchUp, - fMuonTaggingPMchLow, - fMuonTaggingPtMchLow, - {fMuonTaggingEtaMchLow, fMuonTaggingEtaMchUp}, - {fMuonTaggingRabsLow, fMuonTaggingRabsUp}, - fMuonTaggingSigmaPdcaUp)) { + if (!isGoodMuon(muonTrack, collision, + cfgMuonTaggingTrackChi2MchUp, + cfgMuonTaggingPMchLow, + cfgMuonTaggingPtMchLow, + {cfgMuonTaggingEtaMchLow, cfgMuonTaggingEtaMchUp}, + {cfgMuonTaggingRabsLow, cfgMuonTaggingRabsUp}, + cfgMuonTaggingPdcaUp)) { continue; } @@ -1729,7 +1775,7 @@ struct qaMatching { auto mchTrackAtVertex = VarManager::PropagateMuon(muonTrack, collision, VarManager::kToVertex); // propagate the track from the vertex to the first MFT plane - const auto& extrapToMFTfirst = PropagateToZMCH(mchTrackAtVertex, o2::mft::constants::mft::LayerZCoordinate()[0]); + const auto& extrapToMFTfirst = propagateToZMch(mchTrackAtVertex, o2::mft::constants::mft::LayerZCoordinate()[0]); double rFront = std::sqrt(extrapToMFTfirst.getX() * extrapToMFTfirst.getX() + extrapToMFTfirst.getY() * extrapToMFTfirst.getY()); double rMinFront = 3.f; double rMaxFront = 9.f; @@ -1737,7 +1783,7 @@ struct qaMatching { continue; // propagate the track from the vertex to the last MFT plane - const auto& extrapToMFTlast = PropagateToZMCH(mchTrackAtVertex, o2::mft::constants::mft::LayerZCoordinate()[9]); + const auto& extrapToMFTlast = propagateToZMch(mchTrackAtVertex, o2::mft::constants::mft::LayerZCoordinate()[9]); double rBack = std::sqrt(extrapToMFTlast.getX() * extrapToMFTlast.getX() + extrapToMFTlast.getY() * extrapToMFTlast.getY()); double rMinBack = 5.f; double rMaxBack = 12.f; @@ -1751,13 +1797,13 @@ struct qaMatching { // for each MCH standalone track, collect the associated matching candidates template - void GetTaggedMuons(const CollisionInfo& collisionInfo, + void getTaggedMuons(const CollisionInfo& collisionInfo, TMUON const& muonTracks, const std::vector& selectedMuons, std::vector& taggedMuons) { taggedMuons.clear(); - for (auto [mchIndex, globalTracksVector] : collisionInfo.matchingCandidates) { + for (const auto& [mchIndex, globalTracksVector] : collisionInfo.matchingCandidates) { // check if the current muon is selected if (std::find(selectedMuons.begin(), selectedMuons.end(), mchIndex) == selectedMuons.end()) @@ -1773,22 +1819,22 @@ struct qaMatching { auto const& muonTrack1 = muonTracks.rawIteratorAt(globalTracksVector[1].globalTrackId); double chi2diff = muonTrack1.chi2MatchMCHMFT() - muonTrack0.chi2MatchMCHMFT(); - if (chi2diff < fMuonTaggingChi2DiffLow) + if (chi2diff < cfgMuonTaggingChi2DiffLow) continue; taggedMuons.emplace_back(mchIndex); } } - void GetMuonPairs(const CollisionInfo& collisionInfo, + void getMuonPairs(const CollisionInfo& collisionInfo, std::vector& muonPairs, std::vector& globalMuonPairs) { // outer loop over muon tracks - for (auto mchIndex1 : collisionInfo.mchTracks) { + for (const auto& mchIndex1 : collisionInfo.mchTracks) { // inner loop over muon tracks - for (auto mchIndex2 : collisionInfo.mchTracks) { + for (const auto& mchIndex2 : collisionInfo.mchTracks) { // avoid double-counting of muon pairs if (mchIndex2 <= mchIndex1) continue; @@ -1799,10 +1845,10 @@ struct qaMatching { } // outer loop over global muon tracks - for (auto& [mchIndex1, matchingCandidates1] : collisionInfo.matchingCandidates) { + for (const auto& [mchIndex1, matchingCandidates1] : collisionInfo.matchingCandidates) { // inner loop over global muon tracks - for (auto& [mchIndex2, matchingCandidates2] : collisionInfo.matchingCandidates) { + for (const auto& [mchIndex2, matchingCandidates2] : collisionInfo.matchingCandidates) { // avoid double-counting of muon pairs if (mchIndex2 <= mchIndex1) continue; @@ -1813,7 +1859,7 @@ struct qaMatching { } } - double GetMuMuInvariantMass(const o2::mch::TrackParam& track1, const o2::mch::TrackParam& track2) + double getMuMuInvariantMass(const o2::mch::TrackParam& track1, const o2::mch::TrackParam& track2) { ROOT::Math::PxPyPzMVector muon1{ track1.px(), @@ -1832,7 +1878,7 @@ struct qaMatching { return dimuon.M(); } - double GetMuMuInvariantMass(const o2::dataformats::GlobalFwdTrack& track1, const o2::dataformats::GlobalFwdTrack& track2) + double getMuMuInvariantMass(const o2::dataformats::GlobalFwdTrack& track1, const o2::dataformats::GlobalFwdTrack& track2) { ROOT::Math::PxPyPzMVector muon1{ track1.getPx(), @@ -1852,7 +1898,7 @@ struct qaMatching { } template - int GetMftMchMatchAttempts(EVT const& collisions, + int getMftMchMatchAttempts(EVT const& collisions, BC const& bcs, TMUON const& mchTrack, TMFTS const& mftTracks) @@ -1887,7 +1933,7 @@ struct qaMatching { } template - void FillCollisions(EVT const& collisions, + void fillCollisions(EVT const& collisions, BC const& bcs, TMUON const& muonTracks, TMFT const& mftTracks, @@ -1908,19 +1954,8 @@ struct qaMatching { int64_t collisionIndex = collision.globalIndex(); auto bc = bcs.rawIteratorAt(collision.bcId()); - /*// require a minimum BC gap between the current conllision and the previous/next ones - const auto& collisionPrev = collisions.rawIteratorAt(collisionIds[cid-1]); - const auto& collisionNext = collisions.rawIteratorAt(collisionIds[cid+1]); - auto bcPrev = bcs.rawIteratorAt(collisionPrev.bcId()); - auto bcNext = bcs.rawIteratorAt(collisionNext.bcId()); - int64_t deltaPrev = bc.globalBC() - bcPrev.globalBC(); - int64_t deltaNext = bcNext.globalBC() - bc.globalBC(); - if (deltaPrev < 50 || deltaNext < 50) { - continue; - }*/ - // fill collision information for global muon tracks (MFT-MCH-MID matches) - for (auto muonTrack : muonTracks) { + for (const auto& muonTrack : muonTracks) { if (!muonTrack.has_collision()) continue; @@ -1934,18 +1969,18 @@ struct qaMatching { collisionInfo.zVertex = collision.posZ(); if (collisionInfo.matchablePairs.empty()) { - FillMatchablePairs(collisionInfo, muonTracks, mftTracks); + fillMatchablePairs(collisionInfo, muonTracks, mftTracks); } - if (static_cast(muonTrack.trackType()) > 2) { + if (static_cast(muonTrack.trackType()) > GlobalTrackTypeMax) { // standalone MCH or MCH-MID tracks int64_t mchTrackIndex = muonTrack.globalIndex(); collisionInfo.mchTracks.push_back(mchTrackIndex); } else { // global muon tracks (MFT-MCH or MFT-MCH-MID) int64_t muonTrackIndex = muonTrack.globalIndex(); - double matchChi2 = muonTrack.chi2MatchMCHMFT() / 5.f; - double matchScore = chi2ToScore(muonTrack.chi2MatchMCHMFT(), 5, 50.f); + double matchChi2 = muonTrack.chi2MatchMCHMFT() / MatchingDegreesOfFreedom; + double matchScore = chi2ToScore(muonTrack.chi2MatchMCHMFT(), MatchingDegreesOfFreedom, MatchingScoreChi2Max); auto const& mchTrack = muonTrack.template matchMCHTrack_as(); int64_t mchTrackIndex = mchTrack.globalIndex(); auto const& mftTrack = muonTrack.template matchMFTTrack_as(); @@ -1988,7 +2023,7 @@ struct qaMatching { } // fill collision information for MFT standalone tracks - for (auto mftTrack : mftTracks) { + for (const auto& mftTrack : mftTracks) { if (!mftTrack.has_collision()) continue; @@ -2012,19 +2047,23 @@ struct qaMatching { return (track1.matchScore > track2.matchScore); }; - for (auto& [collisionIndex, collisionInfo] : collisionInfos) { - for (auto& [mchIndex, globalTracksVector] : collisionInfo.matchingCandidates) { + for (auto collisionInfoIt = collisionInfos.begin(); collisionInfoIt != collisionInfos.end(); ++collisionInfoIt) { + auto& collisionInfo = collisionInfoIt->second; + for (auto matchingCandidatesIt = collisionInfo.matchingCandidates.begin(); matchingCandidatesIt != collisionInfo.matchingCandidates.end(); ++matchingCandidatesIt) { + auto& mchIndex = matchingCandidatesIt->first; + auto& globalTracksVector = matchingCandidatesIt->second; std::sort(globalTracksVector.begin(), globalTracksVector.end(), compareMatchingScore); const auto& mchTrack = muonTracks.rawIteratorAt(mchIndex); - auto mftMchMatchAttempts = GetMftMchMatchAttempts(collisions, bcs, mchTrack, mftTracks); + auto mftMchMatchAttempts = getMftMchMatchAttempts(collisions, bcs, mchTrack, mftTracks); int ranking = 1; - for (auto& candidate : globalTracksVector) { + for (auto candidateIt = globalTracksVector.begin(); candidateIt != globalTracksVector.end(); ++candidateIt) { + auto& candidate = *candidateIt; const auto& muonTrack = muonTracks.rawIteratorAt(candidate.globalTrackId); candidate.matchRanking = ranking; candidate.matchRankingProd = ranking; - candidate.matchType = GetMatchType(muonTrack, muonTracks, mftTracks, collisionInfo.matchablePairs, ranking); + candidate.matchType = getMatchType(muonTrack, muonTracks, mftTracks, collisionInfo.matchablePairs, ranking); candidate.mftMchMatchAttempts = mftMchMatchAttempts; ranking += 1; } @@ -2033,7 +2072,7 @@ struct qaMatching { } template - void FillMatchingPlotsMC(C const& collision, + void fillMatchingPlotsMc(C const& collision, const CollisionInfo& collisionInfo, TMUON const& muonTracks, TMFT const& mftTracks, @@ -2042,7 +2081,7 @@ struct qaMatching { const std::vector>& matchablePairs, double matchingScoreCut, MatchingPlotter* plotter, - bool verbose = false) + bool /*verbose*/ = false) { int mftTrackMult = collisionInfo.mftTracks.size(); @@ -2051,7 +2090,7 @@ struct qaMatching { for (const auto& [mchIndex, globalTracksVector] : matchingCandidates) { // check if the MCH track belongs to a matchable pair - bool isPairedMCH = IsMatchableMCH(static_cast(mchIndex), matchablePairs); + bool isPairedMCH = isMatchableMch(static_cast(mchIndex), matchablePairs); // get the standalone MCH track auto const& mchTrack = muonTracks.rawIteratorAt(mchIndex); @@ -2059,24 +2098,24 @@ struct qaMatching { double mchPt = mchTrack.pt(); // MCH track quality flag - bool isGoodMCH = IsGoodGlobalMuon(mchTrack, collision); + bool isGoodMCH = isGoodGlobalMuon(mchTrack, collision); - auto matchablePair = GetMatchablePairForMCH(static_cast(mchIndex), matchablePairs); + auto matchablePair = getMatchablePairForMch(static_cast(mchIndex), matchablePairs); bool hasMatchablePair = matchablePair.has_value(); int decayRanking = 0; int mftTrackType = -1; - float dchi2 = (globalTracksVector.size() >= 2) ? globalTracksVector[1].matchChi2 - globalTracksVector[0].matchChi2 : -1; + float dchi2 = (globalTracksVector.size() >= MinCandidatesForDeltaChi2) ? globalTracksVector[1].matchChi2 - globalTracksVector[0].matchChi2 : InvalidDeltaChi2; if (hasMatchablePair) { auto const& pairedMftTrack = mftTracks.rawIteratorAt(matchablePair.value().second); - mftTrackType = pairedMftTrack.isCA() ? 1 : 0; - decayRanking = GetDecayRanking(mchTrack, mftTracks); + mftTrackType = pairedMftTrack.isCA() ? MftTrackTypeCA : MftTrackTypeStandard; + decayRanking = getDecayRanking(mchTrack, mftTracks); } // find the index of the matching candidate that corresponds to the true match // index=1 corresponds to the leading candidate // index=0 means no candidate was found that corresponds to the true match - int trueMatchIndex = GetTrueMatchIndex(muonTracks, globalTracksVector, matchablePairs); - int trueMatchIndexProd = GetTrueMatchIndex(muonTracks, matchingCandidatesProd.at(mchIndex), matchablePairs); + int trueMatchIndex = getTrueMatchIndex(muonTracks, globalTracksVector, matchablePairs); + int trueMatchIndexProd = getTrueMatchIndex(muonTracks, matchingCandidatesProd.at(mchIndex), matchablePairs); float mcParticleDz = -1000; if (mchTrack.has_mcParticle()) { @@ -2189,7 +2228,7 @@ struct qaMatching { // ==================================== // Matching properties - for (auto [mchIndex, globalTracksVector] : matchingCandidates) { + for (const auto& [mchIndex, globalTracksVector] : matchingCandidates) { if (globalTracksVector.size() < 1) continue; @@ -2197,7 +2236,7 @@ struct qaMatching { auto const& mchTrack = muonTracks.rawIteratorAt(mchIndex); // skip global muon tracks that do not pass the MCH and MFT quality cuts - if (!IsGoodGlobalMuon(mchTrack, collision)) + if (!isGoodGlobalMuon(mchTrack, collision)) continue; double mchMom = mchTrack.p(); @@ -2220,11 +2259,11 @@ struct qaMatching { } } - for (auto [mchIndex, globalTracksVector] : matchingCandidates) { + for (const auto& [mchIndex, globalTracksVector] : matchingCandidates) { if (globalTracksVector.size() < 1) continue; - int trueMatchIndex = GetTrueMatchIndex(muonTracks, globalTracksVector, matchablePairs); + int trueMatchIndex = getTrueMatchIndex(muonTracks, globalTracksVector, matchablePairs); // loop over candidates int candidateIndex = 1; @@ -2251,9 +2290,7 @@ struct qaMatching { // ==================================== // Matching purity - if (verbose) - std::cout << std::format(" Filling matching purity plots with score cut {}", matchingScoreCut) << std::endl; - for (auto [mchIndex, globalTracksVector] : matchingCandidates) { + for (const auto& [mchIndex, globalTracksVector] : matchingCandidates) { if (globalTracksVector.size() < 1) continue; @@ -2266,37 +2303,28 @@ struct qaMatching { auto const& mftTrack = muonTrack.template matchMFTTrack_as(); // skip global muon tracks that do not pass the MCH and MFT quality cuts - if (!IsGoodGlobalMuon(mchTrack, collision)) + if (!isGoodGlobalMuon(mchTrack, collision)) continue; - if (!IsGoodMFT(mftTrack)) + if (!isGoodMft(mftTrack)) continue; // skip candidates that do not pass the matching quality cuts - if (!IsGoodGlobalMatching(muonTrack, matchingScore, matchingScoreCut)) + if (!isGoodGlobalMatching(muonTrack, matchingScore, matchingScoreCut)) continue; // check if the matching candidate is a true one - bool isTrueMatch = IsTrueGlobalMatching(muonTrack, matchablePairs); + bool isTrueMatch = isTrueGlobalMatching(muonTrack, matchablePairs); // ---- MC ancestry ---- - auto motherParticles = GetMotherParticles(muonTrack); + auto motherParticles = getMotherParticles(muonTrack); int motherPDG = 0; if (motherParticles.size() > 1) { motherPDG = motherParticles[1].first; } - - if (verbose) { - std::cout << std::format(" MCH track #{} -> Muon track #{}, isTrueMatch={}", mchIndex, globalTracksVector[0].globalTrackId, isTrueMatch) << std::endl; - std::cout << " MC ancestry (pdg): "; - for (auto const& [pdg, idx] : motherParticles) { - std::cout << "(" << pdg << ") "; - } - std::cout << std::endl; - } // fill matching purity plots - plotter->fMatchingPurityPlotter.Fill(mchTrack, isTrueMatch); - if (fConfigQAs.fCreatePdgMomHistograms) { - plotter->fMatchingPurityPlotter.Fill(mchTrack, motherPDG, isTrueMatch); + plotter->fMatchingPurityPlotter.fill(mchTrack, isTrueMatch); + if (configQas.cfgCreatePdgMomHistograms) { + plotter->fMatchingPurityPlotter.fill(mchTrack, motherPDG, isTrueMatch); } } @@ -2304,13 +2332,13 @@ struct qaMatching { // Matching efficiencies // outer loop on matchable pairs - for (auto [matchableMchIndex, matchableMftIndex] : matchablePairs) { + for (const auto& [matchableMchIndex, matchableMftIndex] : matchablePairs) { // get the standalone MCH track auto const& mchTrack = muonTracks.rawIteratorAt(matchableMchIndex); // skip track pairs that do not pass the MCH and MFT quality cuts // we only consider matchable pairs that fulfill the track quality requirements - if (!IsGoodGlobalMuon(mchTrack, collision)) + if (!isGoodGlobalMuon(mchTrack, collision)) continue; bool goodMatchFound = false; @@ -2328,44 +2356,47 @@ struct qaMatching { auto const& mftTrack = muonTrack.template matchMFTTrack_as(); auto mftIndex = mftTrack.globalIndex(); - goodMatchFound = IsGoodGlobalMatching(muonTrack, matchingScore, matchingScoreCut); + goodMatchFound = isGoodGlobalMatching(muonTrack, matchingScore, matchingScoreCut); isTrueMatch = (mftIndex == matchableMftIndex); } } + if (cfgQaMatchingAodDebug > 0 && goodMatchFound && isTrueMatch) { + LOGF(info, + "[good&true] mchId=%lld trackType=%d p=%.3f pt=%.3f eta=%.3f phi=%.3f", + static_cast(mchTrack.globalIndex()), + static_cast(mchTrack.trackType()), + mchTrack.p(), + mchTrack.pt(), + mchTrack.eta(), + mchTrack.phi()); + } + // ---- MC ancestry ---- - auto motherParticles = GetMotherParticles(mchTrack); + auto motherParticles = getMotherParticles(mchTrack); int motherPDG = 0; if (motherParticles.size() > 1) { motherPDG = motherParticles[1].first; } - if (verbose) { - std::cout << " MC ancestry (pdg): "; - for (auto const& [pdg, idx] : motherParticles) { - std::cout << "(" << pdg << ") "; - } - std::cout << std::endl; - } - // fill matching efficiency plots - plotter->fPairingEfficiencyPlotter.Fill(mchTrack, goodMatchFound); - if (fConfigQAs.fCreatePdgMomHistograms) { - plotter->fPairingEfficiencyPlotter.Fill(mchTrack, motherPDG, goodMatchFound); + plotter->fPairingEfficiencyPlotter.fill(mchTrack, goodMatchFound); + if (configQas.cfgCreatePdgMomHistograms) { + plotter->fPairingEfficiencyPlotter.fill(mchTrack, motherPDG, goodMatchFound); } - plotter->fMatchingEfficiencyPlotter.Fill(mchTrack, (goodMatchFound && isTrueMatch)); - if (fConfigQAs.fCreatePdgMomHistograms) { - plotter->fMatchingEfficiencyPlotter.Fill(mchTrack, motherPDG, (goodMatchFound && isTrueMatch)); + plotter->fMatchingEfficiencyPlotter.fill(mchTrack, (goodMatchFound && isTrueMatch)); + if (configQas.cfgCreatePdgMomHistograms) { + plotter->fMatchingEfficiencyPlotter.fill(mchTrack, motherPDG, (goodMatchFound && isTrueMatch)); } - plotter->fFakeMatchingEfficiencyPlotter.Fill(mchTrack, (goodMatchFound && !isTrueMatch)); - if (fConfigQAs.fCreatePdgMomHistograms) { - plotter->fFakeMatchingEfficiencyPlotter.Fill(mchTrack, motherPDG, (goodMatchFound && !isTrueMatch)); + plotter->fFakeMatchingEfficiencyPlotter.fill(mchTrack, (goodMatchFound && !isTrueMatch)); + if (configQas.cfgCreatePdgMomHistograms) { + plotter->fFakeMatchingEfficiencyPlotter.fill(mchTrack, motherPDG, (goodMatchFound && !isTrueMatch)); } } } template - void FillDimuonPlotsMC(const CollisionInfo& collisionInfo, + void fillDimuonPlotsMc(const CollisionInfo& collisionInfo, C const& collisions, TMUON const& muonTracks, TMFT const& /*mftTracks*/) @@ -2373,9 +2404,9 @@ struct qaMatching { std::vector muonPairs; std::vector globalMuonPairs; - GetMuonPairs(collisionInfo, muonPairs, globalMuonPairs); + getMuonPairs(collisionInfo, muonPairs, globalMuonPairs); - for (auto& [muon1, muon2] : muonPairs) { + for (const auto& [muon1, muon2] : muonPairs) { auto const& collision = collisions.rawIteratorAt(muon1.first); auto mchIndex1 = muon1.second; @@ -2389,16 +2420,16 @@ struct qaMatching { if ((sign1 * sign2) >= 0) continue; - bool goodMuonTracks = (IsGoodMuon(muonTrack1, collision) && IsGoodMuon(muonTrack2, collision)); + bool goodMuonTracks = (isGoodMuon(muonTrack1, collision) && isGoodMuon(muonTrack2, collision)); if (goodMuonTracks) { - double mass = GetMuMuInvariantMass(PropagateToVertexMCH(muonTrack1, collision), - PropagateToVertexMCH(muonTrack2, collision)); + double mass = getMuMuInvariantMass(propagateToVertexMch(muonTrack1, collision), + propagateToVertexMch(muonTrack2, collision)); registryDimuon.get(HIST("dimuon/invariantMass_MuonKine_MuonCuts"))->Fill(mass); } } - for (auto& [muon1, muon2] : globalMuonPairs) { + for (const auto& [muon1, muon2] : globalMuonPairs) { auto& candidates1 = muon1.second; auto& candidates2 = muon2.second; @@ -2426,17 +2457,17 @@ struct qaMatching { matchType = candidates2[0].matchType * 10 + candidates1[0].matchType; } - bool goodGlobalMuonTracks = (IsGoodGlobalMuon(mchTrack1, collision) && IsGoodGlobalMuon(mchTrack2, collision)); + bool goodGlobalMuonTracks = (isGoodGlobalMuon(mchTrack1, collision) && isGoodGlobalMuon(mchTrack2, collision)); if (!goodGlobalMuonTracks) { continue; } - bool goodGlobalMuonMatches = (IsGoodGlobalMatching(muonTrack1, matchScore1) && IsGoodGlobalMatching(muonTrack2, matchScore2)); + bool goodGlobalMuonMatches = (isGoodGlobalMatching(muonTrack1, matchScore1) && isGoodGlobalMatching(muonTrack2, matchScore2)); - double massMCH = GetMuMuInvariantMass(PropagateToVertexMCH(mchTrack1, collision), - PropagateToVertexMCH(mchTrack2, collision)); - double mass = GetMuMuInvariantMass(PropagateToVertexMCH(muonTrack1, collision), - PropagateToVertexMCH(muonTrack2, collision)); + double massMCH = getMuMuInvariantMass(propagateToVertexMch(mchTrack1, collision), + propagateToVertexMch(mchTrack2, collision)); + double mass = getMuMuInvariantMass(propagateToVertexMch(muonTrack1, collision), + propagateToVertexMch(muonTrack2, collision)); registryDimuon.get(HIST("dimuon/invariantMass_MuonKine_GlobalMuonCuts"))->Fill(massMCH); registryDimuon.get(HIST("dimuon/invariantMass_ScaledMftKine_GlobalMuonCuts"))->Fill(mass); registryDimuon.get(HIST("dimuon/MC/invariantMass_MuonKine_GlobalMuonCuts_vs_match_type"))->Fill(massMCH, matchType); @@ -2452,7 +2483,7 @@ struct qaMatching { } template - void RunChi2Matching(C const& collisions, + void runChi2Matching(C const& collisions, BC const& bcs, TMUON const& muonTracks, TMFT const& mftTracks, @@ -2475,7 +2506,7 @@ struct qaMatching { return; auto matchingFunc = mMatchingFunctionMap.at(funcName); - for (auto& [mchIndex, globalTracksVector] : matchingCandidates) { + for (const auto& [mchIndex, globalTracksVector] : matchingCandidates) { auto const& mchTrack = muonTracks.rawIteratorAt(mchIndex); for (const auto& candidate : globalTracksVector) { @@ -2495,12 +2526,12 @@ struct qaMatching { auto const& mftTrackCov = mftCovs.rawIteratorAt(mftTrackCovs[mftTrack.globalIndex()]); // get tracks parameters in O2 format - auto mftTrackProp = FwdToTrackPar(mftTrack, mftTrackCov); - auto mchTrackProp = FwdToTrackPar(mchTrack, mchTrack); + auto mftTrackProp = fwdToTrackPar(mftTrack, mftTrackCov); + auto mchTrackProp = fwdToTrackPar(mchTrack, mchTrack); if (matchingPlaneZ < 0.) { - mftTrackProp = PropagateToMatchingPlaneMFT(mchTrack, mftTrack, mftTrackCov, collision, matchingPlaneZ, extrapMethod); - mchTrackProp = PropagateToMatchingPlaneMCH(mchTrack, mftTrack, mftTrackCov, collision, matchingPlaneZ, extrapMethod); + mftTrackProp = propagateToMatchingPlaneMft(mchTrack, mftTrack, mftTrackCov, collision, matchingPlaneZ, extrapMethod); + mchTrackProp = propagateToMatchingPlaneMch(mchTrack, mftTrack, mftTrackCov, collision, matchingPlaneZ, extrapMethod); } // run the chi2 matching function @@ -2548,17 +2579,20 @@ struct qaMatching { return (track1.matchScore > track2.matchScore); }; - for (auto& [mchIndex, globalTracksVector] : newMatchingCandidates) { + for (auto matchingCandidatesIt = newMatchingCandidates.begin(); matchingCandidatesIt != newMatchingCandidates.end(); ++matchingCandidatesIt) { + auto& mchIndex = matchingCandidatesIt->first; + auto& globalTracksVector = matchingCandidatesIt->second; std::sort(globalTracksVector.begin(), globalTracksVector.end(), compareMatchingScore); const auto& mchTrack = muonTracks.rawIteratorAt(mchIndex); - auto mftMchMatchAttempts = GetMftMchMatchAttempts(collisions, bcs, mchTrack, mftTracks); + auto mftMchMatchAttempts = getMftMchMatchAttempts(collisions, bcs, mchTrack, mftTracks); int ranking = 1; - for (auto& candidate : globalTracksVector) { + for (auto candidateIt = globalTracksVector.begin(); candidateIt != globalTracksVector.end(); ++candidateIt) { + auto& candidate = *candidateIt; const auto& muonTrack = muonTracks.rawIteratorAt(candidate.globalTrackId); candidate.matchRanking = ranking; - candidate.matchType = GetMatchType(muonTrack, muonTracks, mftTracks, matchablePairs, ranking); + candidate.matchType = getMatchType(muonTrack, muonTracks, mftTracks, matchablePairs, ranking); candidate.mftMchMatchAttempts = mftMchMatchAttempts; ranking += 1; } @@ -2566,7 +2600,7 @@ struct qaMatching { } template - void RunChi2Matching(C const& collisions, + void runChi2Matching(C const& collisions, BC const& bcs, TMUON const& muonTracks, TMFT const& mftTracks, @@ -2596,11 +2630,11 @@ struct qaMatching { auto matchingPlaneZ = matchingPlanesZ[label]; auto extrapMethod = matchingExtrapMethod[label]; - RunChi2Matching(collisions, bcs, muonTracks, mftTracks, mftCovs, funcName, matchingPlaneZ, extrapMethod, matchablePairs, matchingCandidates, newMatchingCandidates); + runChi2Matching(collisions, bcs, muonTracks, mftTracks, mftCovs, funcName, matchingPlaneZ, extrapMethod, matchablePairs, matchingCandidates, newMatchingCandidates); } template - void RunMLMatching(C const& collisions, + void runMlMatching(C const& collisions, BC const& bcs, TMUON const& muonTracks, TMFT const& mftTracks, @@ -2616,7 +2650,7 @@ struct qaMatching { return; auto& mlResponse = mlIter->second; - for (auto& [mchIndex, globalTracksVector] : matchingCandidates) { + for (const auto& [mchIndex, globalTracksVector] : matchingCandidates) { auto const& mchTrack = muonTracks.rawIteratorAt(mchIndex); for (const auto& candidate : globalTracksVector) { auto const& muonTrack = muonTracks.rawIteratorAt(candidate.globalTrackId); @@ -2636,14 +2670,14 @@ struct qaMatching { // std::cout << fmt::format("Covariance matrix for MFT track #{} retrieved", mftTrack.globalIndex()) << std::endl; // get tracks parameters in O2 format - auto mftTrackProp = FwdToTrackPar(mftTrack, mftTrackCov); - auto mchTrackProp = FwdToTrackPar(mchTrack, mchTrack); + auto mftTrackProp = fwdToTrackPar(mftTrack, mftTrackCov); + auto mchTrackProp = fwdToTrackPar(mchTrack, mchTrack); // extrapolate to the matching plane auto matchingPlaneZ = matchingPlanesZ[label]; if (matchingPlaneZ < 0.) { - mftTrackProp = PropagateToZMFT(mftTrackProp, matchingPlaneZ); - mchTrackProp = PropagateToZMCH(mchTrackProp, matchingPlaneZ); + mftTrackProp = propagateToZMft(mftTrackProp, matchingPlaneZ); + mchTrackProp = propagateToZMch(mchTrackProp, matchingPlaneZ); } // run the ML model @@ -2651,8 +2685,8 @@ struct qaMatching { std::vector inputML = mlResponse.getInputFeaturesGlob(muonTrack, mchTrackProp, mftTrackProp, collision); mlResponse.isSelectedMl(inputML, 0, output); float matchScore = output[0]; - float matchChi2Prod = muonTrack.chi2MatchMCHMFT() / 5.f; - float matchScoreProd = chi2ToScore(muonTrack.chi2MatchMCHMFT(), 5, 50.f); + float matchChi2Prod = muonTrack.chi2MatchMCHMFT() / MatchingDegreesOfFreedom; + float matchScoreProd = chi2ToScore(muonTrack.chi2MatchMCHMFT(), MatchingDegreesOfFreedom, MatchingScoreChi2Max); // std::cout << std::format("Matching score: {}, Chi2: {}", matchingScore, muonTrack.chi2MatchMCHMFT()) << std::endl; // check if a vector of global muon candidates is already available for the current MCH index @@ -2693,17 +2727,20 @@ struct qaMatching { return (track1.matchScore > track2.matchScore); }; - for (auto& [mchIndex, globalTracksVector] : newMatchingCandidates) { + for (auto matchingCandidatesIt = newMatchingCandidates.begin(); matchingCandidatesIt != newMatchingCandidates.end(); ++matchingCandidatesIt) { + auto& mchIndex = matchingCandidatesIt->first; + auto& globalTracksVector = matchingCandidatesIt->second; std::sort(globalTracksVector.begin(), globalTracksVector.end(), compareMatchingScore); const auto& mchTrack = muonTracks.rawIteratorAt(mchIndex); - auto mftMchMatchAttempts = GetMftMchMatchAttempts(collisions, bcs, mchTrack, mftTracks); + auto mftMchMatchAttempts = getMftMchMatchAttempts(collisions, bcs, mchTrack, mftTracks); int ranking = 1; - for (auto& candidate : globalTracksVector) { + for (auto candidateIt = globalTracksVector.begin(); candidateIt != globalTracksVector.end(); ++candidateIt) { + auto& candidate = *candidateIt; const auto& muonTrack = muonTracks.rawIteratorAt(candidate.globalTrackId); candidate.matchRanking = ranking; - candidate.matchType = GetMatchType(muonTrack, muonTracks, mftTracks, matchablePairs, ranking); + candidate.matchType = getMatchType(muonTrack, muonTracks, mftTracks, matchablePairs, ranking); candidate.mftMchMatchAttempts = mftMchMatchAttempts; ranking += 1; } @@ -2711,7 +2748,7 @@ struct qaMatching { } template - void ProcessCollisionMC(const CollisionInfo& collisionInfo, + void processCollisionMc(const CollisionInfo& collisionInfo, C const& collisions, BC const& bcs, TMUON const& muonTracks, @@ -2724,30 +2761,30 @@ struct qaMatching { registry.get(HIST("tracksMultiplicityMCH"))->Fill(collisionInfo.mchTracks.size()); // Chi2-based matching analysis - FillMatchingPlotsMC(collision, collisionInfo, muonTracks, mftTracks, collisionInfo.matchingCandidates, collisionInfo.matchingCandidates, collisionInfo.matchablePairs, fMatchingChi2ScoreMftMchLow, fChi2MatchingPlotter.get(), false); - for (auto& [label, func] : matchingChi2Functions) { + fillMatchingPlotsMc(collision, collisionInfo, muonTracks, mftTracks, collisionInfo.matchingCandidates, collisionInfo.matchingCandidates, collisionInfo.matchablePairs, cfgMatchingChi2ScoreMftMchLow, fChi2MatchingPlotter.get(), false); + for (const auto& [label, func] : matchingChi2Functions) { MatchingCandidates matchingCandidates; - RunChi2Matching(collisions, bcs, muonTracks, mftTracks, mftCovs, label, collisionInfo.matchablePairs, collisionInfo.matchingCandidates, matchingCandidates); + runChi2Matching(collisions, bcs, muonTracks, mftTracks, mftCovs, label, collisionInfo.matchablePairs, collisionInfo.matchingCandidates, matchingCandidates); auto* plotter = fMatchingPlotters.at(label).get(); double matchingScoreCut = matchingScoreCuts.at(label); - FillMatchingPlotsMC(collision, collisionInfo, muonTracks, mftTracks, matchingCandidates, collisionInfo.matchingCandidates, collisionInfo.matchablePairs, matchingScoreCut, plotter, false); + fillMatchingPlotsMc(collision, collisionInfo, muonTracks, mftTracks, matchingCandidates, collisionInfo.matchingCandidates, collisionInfo.matchablePairs, matchingScoreCut, plotter, false); } // ML-based matching analysis - for (auto& [label, mlResponse] : matchingMlResponses) { + for (const auto& [label, mlResponse] : matchingMlResponses) { MatchingCandidates matchingCandidates; - RunMLMatching(collisions, bcs, muonTracks, mftTracks, mftCovs, label, collisionInfo.matchablePairs, collisionInfo.matchingCandidates, matchingCandidates); + runMlMatching(collisions, bcs, muonTracks, mftTracks, mftCovs, label, collisionInfo.matchablePairs, collisionInfo.matchingCandidates, matchingCandidates); auto* plotter = fMatchingPlotters.at(label).get(); double matchingScoreCut = matchingScoreCuts.at(label); - FillMatchingPlotsMC(collision, collisionInfo, muonTracks, mftTracks, matchingCandidates, collisionInfo.matchingCandidates, collisionInfo.matchablePairs, matchingScoreCut, plotter); + fillMatchingPlotsMc(collision, collisionInfo, muonTracks, mftTracks, matchingCandidates, collisionInfo.matchingCandidates, collisionInfo.matchablePairs, matchingScoreCut, plotter); } // Muons tagging - for (auto [mchIndex, mftIndex] : collisionInfo.matchablePairs) { + for (const auto& [mchIndex, mftIndex] : collisionInfo.matchablePairs) { auto const& mchTrack = muonTracks.rawIteratorAt(mchIndex); if (!mchTrack.has_collision()) continue; @@ -2763,37 +2800,141 @@ struct qaMatching { // extrapolate to the matching plane auto z = o2::mft::constants::mft::LayerZCoordinate()[9]; - auto mchTrackProp = PropagateToZMCH(mchTrackAtVertex, z); - auto mftTrackProp = PropagateToZMFT(FwdToTrackPar(mftTrack, mftTrackCov), z); + auto mchTrackProp = propagateToZMch(mchTrackAtVertex, z); + auto mftTrackProp = propagateToZMft(fwdToTrackPar(mftTrack, mftTrackCov), z); registry.get(HIST("matching/MC/pairedMCHTracksAtMFT"))->Fill(mchTrackProp.getX(), mchTrackProp.getY()); registry.get(HIST("matching/MC/pairedMFTTracksAtMFT"))->Fill(mftTrackProp.getX(), mftTrackProp.getY()); } std::vector selectedMuons; - GetSelectedMuons(collisionInfo, collisions, muonTracks, selectedMuons); + getSelectedMuons(collisionInfo, collisions, muonTracks, selectedMuons); MatchingCandidates selectedMatchingCandidates; - for (auto [mchIndex, globalTracksVector] : collisionInfo.matchingCandidates) { + for (const auto& [mchIndex, globalTracksVector] : collisionInfo.matchingCandidates) { if (std::find(selectedMuons.begin(), selectedMuons.end(), mchIndex) != selectedMuons.end()) { selectedMatchingCandidates[mchIndex] = globalTracksVector; } } - FillMatchingPlotsMC(collision, collisionInfo, muonTracks, mftTracks, selectedMatchingCandidates, collisionInfo.matchingCandidates, collisionInfo.matchablePairs, fMatchingChi2ScoreMftMchLow, fSelectedMuonsMatchingPlotter.get()); + fillMatchingPlotsMc(collision, collisionInfo, muonTracks, mftTracks, selectedMatchingCandidates, collisionInfo.matchingCandidates, collisionInfo.matchablePairs, cfgMatchingChi2ScoreMftMchLow, fSelectedMuonsMatchingPlotter.get()); std::vector taggedMuons; - GetTaggedMuons(collisionInfo, muonTracks, selectedMuons, taggedMuons); + getTaggedMuons(collisionInfo, muonTracks, selectedMuons, taggedMuons); MatchingCandidates taggedMatchingCandidates; - for (auto [mchIndex, globalTracksVector] : collisionInfo.matchingCandidates) { + for (const auto& [mchIndex, globalTracksVector] : collisionInfo.matchingCandidates) { if (std::find(taggedMuons.begin(), taggedMuons.end(), mchIndex) != taggedMuons.end()) { taggedMatchingCandidates[mchIndex] = globalTracksVector; } } - FillMatchingPlotsMC(collision, collisionInfo, muonTracks, mftTracks, taggedMatchingCandidates, collisionInfo.matchingCandidates, collisionInfo.matchablePairs, fMatchingChi2ScoreMftMchLow, fTaggedMuonsMatchingPlotter.get()); + fillMatchingPlotsMc(collision, collisionInfo, muonTracks, mftTracks, taggedMatchingCandidates, collisionInfo.matchingCandidates, collisionInfo.matchablePairs, cfgMatchingChi2ScoreMftMchLow, fTaggedMuonsMatchingPlotter.get()); // Di-muon analysis - FillDimuonPlotsMC(collisionInfo, collisions, muonTracks, mftTracks); + fillDimuonPlotsMc(collisionInfo, collisions, muonTracks, mftTracks); + } + + template + void fillQaMatchingAodTablesForCollision(TCOLLISION const& collision, + TMUON const& muonTracks, + const MatchingCandidates& matchingCandidates, + int8_t matchLabel, + int32_t reducedEventId) + { + for (const auto& [mchIndex, candidates] : matchingCandidates) { + if (candidates.empty()) { + continue; + } + + const auto& mchTrack = muonTracks.rawIteratorAt(mchIndex); + if (!isGoodGlobalMuon(mchTrack, collision)) { + continue; + } + + for (const auto& candidate : candidates) { + const auto& candidateTrack = muonTracks.rawIteratorAt(candidate.globalTrackId); + auto candidateTrackAtVertex = VarManager::PropagateMuon(candidateTrack, collision, VarManager::kToVertex); + qaMatchingCandidates( + reducedEventId, + matchLabel, + mchIndex, + static_cast(candidateTrack.p()), + static_cast(candidateTrack.pt()), + static_cast(candidateTrack.eta()), + static_cast(candidateTrack.phi()), + static_cast(candidate.matchType), + static_cast(candidate.matchScore), + static_cast(candidate.matchRanking), + static_cast(candidateTrackAtVertex.getX()), + static_cast(candidateTrackAtVertex.getY()), + static_cast(candidateTrackAtVertex.getZ()), + static_cast(candidateTrackAtVertex.getPx()), + static_cast(candidateTrackAtVertex.getPy()), + static_cast(candidateTrackAtVertex.getPz())); + } + } + } + + template + void fillQaMatchingAodEventForCollision(const CollisionInfo& collisionInfo, + TCOLLISION const& collision, + int32_t reducedEventId, + int& debugCounter) + { + int32_t mftMultiplicity = static_cast(collisionInfo.mftTracks.size()); + qaMatchingEvents( + mftMultiplicity, + static_cast(collision.posX()), + static_cast(collision.posY()), + static_cast(collision.posZ())); + + if (cfgQaMatchingAodDebug > 0 && debugCounter < cfgQaMatchingAodDebug) { + LOGF(info, "[AO2D] reducedEvent=%", reducedEventId); + debugCounter += 1; + } + } + + template + void fillQaMatchingMchTracksForCollision(const CollisionInfo& collisionInfo, + TCOLLISIONS const& collisions, + TCOLLISION const& collision, + TMUON const& muonTracks, + TMFT const& mftTracks, + TBC const& bcs, + int32_t reducedEventId) + { + std::vector mchIds; + for (const auto& mchIndex : collisionInfo.mchTracks) { + if (std::find(mchIds.begin(), mchIds.end(), mchIndex) == mchIds.end()) { + mchIds.emplace_back(mchIndex); + } + } + for (const auto& [mchIndex, candidates] : collisionInfo.matchingCandidates) { + (void)candidates; + if (std::find(mchIds.begin(), mchIds.end(), mchIndex) == mchIds.end()) { + mchIds.emplace_back(mchIndex); + } + } + + for (const auto& mchIndex : mchIds) { + auto const& mchTrack = muonTracks.rawIteratorAt(mchIndex); + int mftMchMatchAttempts = getMftMchMatchAttempts(collisions, bcs, mchTrack, mftTracks); + auto mchTrackAtVertex = VarManager::PropagateMuon(mchTrack, collision, VarManager::kToVertex); + qaMatchingMCHTrack( + reducedEventId, + mchIndex, + static_cast(mchTrack.trackType()), + static_cast(mchTrack.p()), + static_cast(mchTrack.pt()), + static_cast(mchTrack.eta()), + static_cast(mchTrack.phi()), + static_cast(mftMchMatchAttempts), + static_cast(mchTrackAtVertex.getX()), + static_cast(mchTrackAtVertex.getY()), + static_cast(mchTrackAtVertex.getZ()), + static_cast(mchTrackAtVertex.getPx()), + static_cast(mchTrackAtVertex.getPy()), + static_cast(mchTrackAtVertex.getPz())); + } } void processQAMC(MyEvents const& collisions, @@ -2804,29 +2945,72 @@ struct qaMatching { aod::McParticles const& /*mcParticles*/) { auto bc = bcs.begin(); - initCCDB(bc); + initCcdb(bc); - for (auto& muon : muonTracks) { + for (const auto& muon : muonTracks) { registry.get(HIST("nTracksPerType"))->Fill(static_cast(muon.trackType())); } - FillCollisions(collisions, bcs, muonTracks, mftTracks, fCollisionInfos); + fillCollisions(collisions, bcs, muonTracks, mftTracks, fCollisionInfos); mftTrackCovs.clear(); - for (auto& mftTrackCov : mftCovs) { + for (const auto& mftTrackCov : mftCovs) { mftTrackCovs[mftTrackCov.matchMFTTrackId()] = mftTrackCov.globalIndex(); } + std::unordered_map reducedEventIds; + int32_t reducedEventCounter = 0; + for (auto const& [collisionIndex, collisionInfo] : fCollisionInfos) { + reducedEventIds.emplace(collisionInfo.index, reducedEventCounter); + reducedEventCounter += 1; + } + + int debugCounter = 0; + for (auto const& [collisionIndex, collisionInfo] : fCollisionInfos) { + auto it = reducedEventIds.find(collisionInfo.index); + if (it == reducedEventIds.end()) { + continue; + } + int32_t reducedEventId = it->second; + auto collision = collisions.rawIteratorAt(collisionInfo.index); + fillQaMatchingAodEventForCollision(collisionInfo, collision, reducedEventId, debugCounter); + fillQaMatchingMchTracksForCollision(collisionInfo, collisions, collision, muonTracks, mftTracks, bcs, reducedEventId); + } + + struct AodLabel { + const char* name; + int8_t id; + }; + std::array aodLabels{{{"ProdAll", 0}, {"MatchXYPhiTanl", 1}, {"MatchXYPhiTanlMom", 2}}}; + for (const auto& aodLabel : aodLabels) { + if (matchingChi2Functions.find(aodLabel.name) == matchingChi2Functions.end()) { + LOGF(warn, "[AO2D] Chi2 label not found: %s", aodLabel.name); + continue; + } + debugCounter = 0; + for (auto const& [collisionIndex, collisionInfo] : fCollisionInfos) { + auto it = reducedEventIds.find(collisionInfo.index); + if (it == reducedEventIds.end()) { + continue; + } + int32_t reducedEventId = it->second; + MatchingCandidates matchingCandidates; + runChi2Matching(collisions, bcs, muonTracks, mftTracks, mftCovs, aodLabel.name, collisionInfo.matchablePairs, collisionInfo.matchingCandidates, matchingCandidates); + auto collision = collisions.rawIteratorAt(collisionInfo.index); + fillQaMatchingAodTablesForCollision(collision, muonTracks, matchingCandidates, aodLabel.id, reducedEventId); + } + } + for (auto const& [collisionIndex, collisionInfo] : fCollisionInfos) { - ProcessCollisionMC(collisionInfo, collisions, bcs, muonTracks, mftTracks, mftCovs); + processCollisionMc(collisionInfo, collisions, bcs, muonTracks, mftTracks, mftCovs); } } - PROCESS_SWITCH(qaMatching, processQAMC, "processQAMC", true); + PROCESS_SWITCH(QaMatching, processQAMC, "processQAMC", true); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ - adaptAnalysisTask(cfgc)}; + adaptAnalysisTask(cfgc)}; };