From de994c03751790a0b62ee9851313f36fc139b096 Mon Sep 17 00:00:00 2001 From: Daiki Sekihata Date: Sat, 11 Oct 2025 17:37:02 +0200 Subject: [PATCH 1/2] PWGEM/Dilepton, Common: update muon analyses --- Common/Core/fwdtrackUtilities.h | 93 +++-- PWGEM/Dilepton/DataModel/dileptonTables.h | 44 ++- .../TableProducer/skimmerPrimaryMuon.cxx | 235 ++++++++++-- .../Dilepton/Tasks/Converters/CMakeLists.txt | 5 + .../Tasks/Converters/muonConverter1.cxx | 53 +++ .../Tasks/Converters/trackConverter1.cxx | 2 +- PWGEM/Dilepton/Tasks/createResolutionMap.cxx | 166 +++++++- PWGEM/Dilepton/Tasks/matchingMFT.cxx | 362 ++++++++++++++---- 8 files changed, 825 insertions(+), 135 deletions(-) create mode 100644 PWGEM/Dilepton/Tasks/Converters/muonConverter1.cxx diff --git a/Common/Core/fwdtrackUtilities.h b/Common/Core/fwdtrackUtilities.h index 0533bfea89c..a07f461394d 100644 --- a/Common/Core/fwdtrackUtilities.h +++ b/Common/Core/fwdtrackUtilities.h @@ -18,6 +18,7 @@ #ifndef COMMON_CORE_FWDTRACKUTILITIES_H_ #define COMMON_CORE_FWDTRACKUTILITIES_H_ +#include "Framework/AnalysisDataModel.h" #include #include #include @@ -29,6 +30,7 @@ #include #include +#include #include namespace o2::aod @@ -40,30 +42,71 @@ enum class propagationPoint : int { kToVertex = 0, kToDCA = 1, kToRabs = 2, + kToMatchingPlane = 3, }; using SMatrix55 = ROOT::Math::SMatrix>; using SMatrix55Std = ROOT::Math::SMatrix; using SMatrix5 = ROOT::Math::SVector; -/// propagate fwdtrack to a certain point. -template -o2::dataformats::GlobalFwdTrack propagateMuon(TFwdTrack const& muon, TCollision const& collision, const propagationPoint endPoint) +template +o2::track::TrackParCovFwd getTrackParCovFwd(TFwdTrack const& track, TFwdTrackCov const& cov) { - double chi2 = muon.chi2(); - SMatrix5 tpars(muon.x(), muon.y(), muon.phi(), muon.tgl(), muon.signed1Pt()); - std::vector v1{muon.cXX(), muon.cXY(), muon.cYY(), muon.cPhiX(), muon.cPhiY(), - muon.cPhiPhi(), muon.cTglX(), muon.cTglY(), muon.cTglPhi(), muon.cTglTgl(), - muon.c1PtX(), muon.c1PtY(), muon.c1PtPhi(), muon.c1PtTgl(), muon.c1Pt21Pt2()}; + // This function works for both (saMuon, saMuon) and (MFTTrack, MFTTrackCov). + // Don't use covariant matrix of global muons stored in AO2D.root. + + double chi2 = track.chi2(); + if constexpr (std::is_same_v, aod::MFTTracksCov::iterator>) { + chi2 = track.chi2(); + } else { + if (track.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { + chi2 = track.chi2(); + } else if (track.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { + chi2 = track.chi2() * (2.f * track.nClusters() - 5.f); + } + } + + SMatrix5 tpars(track.x(), track.y(), track.phi(), track.tgl(), track.signed1Pt()); + std::vector v1{cov.cXX(), cov.cXY(), cov.cYY(), cov.cPhiX(), cov.cPhiY(), + cov.cPhiPhi(), cov.cTglX(), cov.cTglY(), cov.cTglPhi(), cov.cTglTgl(), + cov.c1PtX(), cov.c1PtY(), cov.c1PtPhi(), cov.c1PtTgl(), cov.c1Pt21Pt2()}; SMatrix55 tcovs(v1.begin(), v1.end()); - o2::track::TrackParCovFwd fwdtrack{muon.z(), tpars, tcovs, chi2}; + o2::track::TrackParCovFwd trackparCov{track.z(), tpars, tcovs, chi2}; // this is chi2! Not chi2/ndf. + v1.clear(); + v1.shrink_to_fit(); + return trackparCov; +} + +/// propagate fwdtrack to a certain point. +template +o2::dataformats::GlobalFwdTrack propagateMuon(TFwdTrack const& muon, TFwdTrackCov const& cov, TCollision const& collision, const propagationPoint endPoint, const float matchingZ, const float bzkG) +{ + o2::track::TrackParCovFwd trackParCovFwd; + if (muon.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { + trackParCovFwd = getTrackParCovFwd(muon, cov); + } else if (muon.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { + trackParCovFwd = getTrackParCovFwd(muon, muon); + } else { + trackParCovFwd = getTrackParCovFwd(muon, muon); + } + + o2::dataformats::GlobalFwdTrack propmuon = propagateTrackParCovFwd(trackParCovFwd, muon.trackType(), collision, endPoint, matchingZ, bzkG); + return propmuon; +} + +template +o2::dataformats::GlobalFwdTrack propagateTrackParCovFwd(TFwdTrackParCov const& fwdtrackORG, uint8_t trackType, TCollision const& collision, const propagationPoint endPoint, const float matchingZ, const float bzkG) +{ + // TFwdTrackParCov is o2::track::TrackParCovFwd + + o2::track::TrackParCovFwd fwdtrack(fwdtrackORG); o2::dataformats::GlobalFwdTrack propmuon; o2::globaltracking::MatchGlobalFwd mMatching; - if (static_cast(muon.trackType()) > 2) { // MCH-MID or MCH standalone + if (trackType > 2) { // MCH-MID or MCH standalone o2::dataformats::GlobalFwdTrack track; - track.setParameters(tpars); + track.setParameters(fwdtrack.getParameters()); track.setZ(fwdtrack.getZ()); - track.setCovariances(tcovs); + track.setCovariances(fwdtrack.getCovariances()); auto mchTrack = mMatching.FwdtoMCH(track); if (endPoint == propagationPoint::kToVertex) { @@ -72,37 +115,42 @@ o2::dataformats::GlobalFwdTrack propagateMuon(TFwdTrack const& muon, TCollision o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchTrack, collision.posZ()); } else if (endPoint == propagationPoint::kToRabs) { o2::mch::TrackExtrap::extrapToZ(mchTrack, -505.); + } else if (endPoint == propagationPoint::kToMatchingPlane) { + o2::mch::TrackExtrap::extrapToVertexWithoutBranson(mchTrack, matchingZ); } auto proptrack = mMatching.MCHtoFwd(mchTrack); propmuon.setParameters(proptrack.getParameters()); propmuon.setZ(proptrack.getZ()); propmuon.setCovariances(proptrack.getCovariances()); - } else if (static_cast(muon.trackType()) < 2) { // MFT-MCH-MID - const double centerMFT[3] = {0, 0, -61.4}; - o2::field::MagneticField* field = static_cast(TGeoGlobalMagField::Instance()->GetField()); - auto Bz = field->getBz(centerMFT); // Get field at centre of MFT - auto geoMan = o2::base::GeometryManager::meanMaterialBudget(muon.x(), muon.y(), muon.z(), collision.posX(), collision.posY(), collision.posZ()); + } else if (trackType < 2) { // MFT-MCH-MID + // const double centerMFT[3] = {0, 0, -61.4}; + // o2::field::MagneticField* field = static_cast(TGeoGlobalMagField::Instance()->GetField()); + // auto Bz = field->getBz(centerMFT); // Get field at centre of MFT in kG. + + auto geoMan = o2::base::GeometryManager::meanMaterialBudget(fwdtrack.getX(), fwdtrack.getY(), fwdtrack.getZ(), collision.posX(), collision.posY(), collision.posZ()); auto x2x0 = static_cast(geoMan.meanX2X0); + if (endPoint == propagationPoint::kToVertex) { - fwdtrack.propagateToVtxhelixWithMCS(collision.posZ(), {collision.posX(), collision.posY()}, {collision.covXX(), collision.covYY()}, Bz, x2x0); + fwdtrack.propagateToVtxhelixWithMCS(collision.posZ(), {collision.posX(), collision.posY()}, {collision.covXX(), collision.covYY()}, bzkG, x2x0); } else if (endPoint == propagationPoint::kToDCA) { - fwdtrack.propagateToZhelix(collision.posZ(), Bz); + fwdtrack.propagateToZhelix(collision.posZ(), bzkG); + } else if (endPoint == propagationPoint::kToMatchingPlane) { + fwdtrack.propagateToZhelix(matchingZ, bzkG); } propmuon.setParameters(fwdtrack.getParameters()); propmuon.setZ(fwdtrack.getZ()); propmuon.setCovariances(fwdtrack.getCovariances()); } - v1.clear(); - v1.shrink_to_fit(); - return propmuon; } template o2::dataformats::GlobalFwdTrack refitGlobalMuonCov(TFwdTrack const& muon, TMFTTrack const& mft) { + // TFwdTrack and TMFTTrack are o2::track::TrackParCovFwd. + auto muonCov = muon.getCovariances(); auto mftCov = mft.getCovariances(); @@ -135,6 +183,7 @@ o2::dataformats::GlobalFwdTrack refitGlobalMuonCov(TFwdTrack const& muon, TMFTTr o2::dataformats::GlobalFwdTrack globalTrack; globalTrack.setParameters(mft.getParameters()); + globalTrack.setZ(mft.getZ()); globalTrack.setInvQPt(invQPtGlob); globalTrack.setCovariances(globalCov); diff --git a/PWGEM/Dilepton/DataModel/dileptonTables.h b/PWGEM/Dilepton/DataModel/dileptonTables.h index 9d03b354d77..ce9a36977b3 100644 --- a/PWGEM/Dilepton/DataModel/dileptonTables.h +++ b/PWGEM/Dilepton/DataModel/dileptonTables.h @@ -706,12 +706,17 @@ DECLARE_SOA_COLUMN(MFTTrackId, mfttrackId, int); DECLARE_SOA_COLUMN(MCHTrackId, mchtrackId, int); //! DECLARE_SOA_SELF_ARRAY_INDEX_COLUMN(GlobalMuonsWithSameMFT, globalMuonsWithSameMFT); //! self indices to global muons that have the same MFTTrackId DECLARE_SOA_SELF_ARRAY_INDEX_COLUMN(AmbiguousMuons, ambiguousMuons); -DECLARE_SOA_COLUMN(CXXatDCA, cXXatDCA, float); //! DCAx resolution squared at DCA -DECLARE_SOA_COLUMN(CYYatDCA, cYYatDCA, float); //! DCAy resolution squared at DCA -DECLARE_SOA_COLUMN(CXYatDCA, cXYatDCA, float); //! correlation term of DCAx,y resolution at DCA -DECLARE_SOA_COLUMN(PtMatchedMCHMID, ptMatchedMCHMID, float); //! pt of MCH-MID track in MFT-MCH-MID track at PV -DECLARE_SOA_COLUMN(EtaMatchedMCHMID, etaMatchedMCHMID, float); //! eta of MCH-MID track in MFT-MCH-MID track at PV -DECLARE_SOA_COLUMN(PhiMatchedMCHMID, phiMatchedMCHMID, float); //! phi of MCH-MID track in MFT-MCH-MID track at PV +DECLARE_SOA_COLUMN(CXXatDCA, cXXatDCA, float); //! DCAx resolution squared at DCA +DECLARE_SOA_COLUMN(CYYatDCA, cYYatDCA, float); //! DCAy resolution squared at DCA +DECLARE_SOA_COLUMN(CXYatDCA, cXYatDCA, float); //! correlation term of DCAx,y resolution at DCA +DECLARE_SOA_COLUMN(PtMatchedMCHMID, ptMatchedMCHMID, float); //! pt of MCH-MID track in MFT-MCH-MID track at PV +DECLARE_SOA_COLUMN(EtaMatchedMCHMID, etaMatchedMCHMID, float); //! eta of MCH-MID track in MFT-MCH-MID track at PV +DECLARE_SOA_COLUMN(PhiMatchedMCHMID, phiMatchedMCHMID, float); //! phi of MCH-MID track in MFT-MCH-MID track at PV +DECLARE_SOA_COLUMN(EtaMatchedMCHMIDatMP, etaMatchedMCHMIDatMP, float); //! eta of MCH-MID track in MFT-MCH-MID track at matching plane +DECLARE_SOA_COLUMN(PhiMatchedMCHMIDatMP, phiMatchedMCHMIDatMP, float); //! phi of MCH-MID track in MFT-MCH-MID track at matching plane +DECLARE_SOA_COLUMN(EtaMatchedMFTatMP, etaMatchedMFTatMP, float); //! eta of MFT track in MFT-MCH-MID track at matching plane +DECLARE_SOA_COLUMN(PhiMatchedMFTatMP, phiMatchedMFTatMP, float); //! phi of MFT track in MFT-MCH-MID track at matching plane + DECLARE_SOA_COLUMN(IsAssociatedToMPC, isAssociatedToMPC, bool); //! is associated to most probable collision DECLARE_SOA_COLUMN(IsAmbiguous, isAmbiguous, bool); //! is ambiguous DECLARE_SOA_COLUMN(Sign, sign, int8_t); //! @@ -742,7 +747,7 @@ DECLARE_SOA_DYNAMIC_COLUMN(MFTClusterMap, mftClusterMap, //! MFT cluster map, on return clmap; }); } // namespace emprimarymuon -DECLARE_SOA_TABLE(EMPrimaryMuons, "AOD", "EMPRIMARYMU", //! +DECLARE_SOA_TABLE(EMPrimaryMuons_000, "AOD", "EMPRIMARYMU", //! o2::soa::Index<>, emprimarymuon::CollisionId, emprimarymuon::FwdTrackId, emprimarymuon::MFTTrackId, emprimarymuon::MCHTrackId, fwdtrack::TrackType, fwdtrack::Pt, fwdtrack::Eta, fwdtrack::Phi, emprimarymuon::Sign, @@ -763,6 +768,31 @@ DECLARE_SOA_TABLE(EMPrimaryMuons, "AOD", "EMPRIMARYMU", //! emprimarymuon::Px, emprimarymuon::Py, emprimarymuon::Pz); + +DECLARE_SOA_TABLE_VERSIONED(EMPrimaryMuons_001, "AOD", "EMPRIMARYMU", 1, //! + o2::soa::Index<>, emprimarymuon::CollisionId, + emprimarymuon::FwdTrackId, emprimarymuon::MFTTrackId, emprimarymuon::MCHTrackId, fwdtrack::TrackType, + fwdtrack::Pt, fwdtrack::Eta, fwdtrack::Phi, emprimarymuon::Sign, + fwdtrack::FwdDcaX, fwdtrack::FwdDcaY, emprimarymuon::CXXatDCA, emprimarymuon::CYYatDCA, emprimarymuon::CXYatDCA, + emprimarymuon::PtMatchedMCHMID, emprimarymuon::EtaMatchedMCHMID, emprimarymuon::PhiMatchedMCHMID, + emprimarymuon::EtaMatchedMCHMIDatMP, emprimarymuon::PhiMatchedMCHMIDatMP, + emprimarymuon::EtaMatchedMFTatMP, emprimarymuon::PhiMatchedMFTatMP, + + fwdtrack::NClusters, fwdtrack::PDca, fwdtrack::RAtAbsorberEnd, + fwdtrack::Chi2, fwdtrack::Chi2MatchMCHMID, fwdtrack::Chi2MatchMCHMFT, + fwdtrack::MCHBitMap, fwdtrack::MIDBitMap, fwdtrack::MIDBoards, + fwdtrack::MFTClusterSizesAndTrackFlags, emprimarymuon::Chi2MFT, emprimarymuon::IsAssociatedToMPC, emprimarymuon::IsAmbiguous, + + // dynamic column + emprimarymuon::Signed1Pt, + emprimarymuon::NClustersMFT, + emprimarymuon::MFTClusterMap, + emprimarymuon::P, + emprimarymuon::Px, + emprimarymuon::Py, + emprimarymuon::Pz); + +using EMPrimaryMuons = EMPrimaryMuons_001; // iterators using EMPrimaryMuon = EMPrimaryMuons::iterator; diff --git a/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuon.cxx b/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuon.cxx index 280c80b8320..656f9136e12 100644 --- a/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuon.cxx +++ b/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuon.cxx @@ -63,7 +63,7 @@ struct skimmerPrimaryMuon { using MFTTracksMC = soa::Join; using MFTTrackMC = MFTTracksMC::iterator; - Produces emprimarymuons; + Produces emprimarymuons; Produces emprimarymuonscov; // Configurables @@ -71,7 +71,7 @@ struct skimmerPrimaryMuon { Configurable grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; Configurable geoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"}; Configurable fillQAHistograms{"fillQAHistograms", false, "flag to fill QA histograms"}; - Configurable minPt{"minPt", 0.2, "min pt for muon"}; + Configurable minPt{"minPt", 0.1, "min pt for muon"}; Configurable maxPt{"maxPt", 1e+10, "max pt for muon"}; Configurable minEtaSA{"minEtaSA", -4.0, "min. eta acceptance for MCH-MID"}; Configurable maxEtaSA{"maxEtaSA", -2.5, "max. eta acceptance for MCH-MID"}; @@ -81,17 +81,19 @@ struct skimmerPrimaryMuon { Configurable minRabs{"minRabs", 17.6, "min. R at absorber end"}; Configurable midRabs{"midRabs", 26.5, "middle R at absorber end for pDCA cut"}; Configurable maxRabs{"maxRabs", 89.5, "max. R at absorber end"}; - Configurable maxDCAxy{"maxDCAxy", 1e+10, "max. DCAxy for global muons"}; + Configurable maxDCAxy{"maxDCAxy", 1, "max. DCAxy for global muons"}; Configurable maxPDCAforLargeR{"maxPDCAforLargeR", 324.f, "max. pDCA for large R at absorber end"}; Configurable maxPDCAforSmallR{"maxPDCAforSmallR", 594.f, "max. pDCA for small R at absorber end"}; Configurable maxMatchingChi2MCHMFT{"maxMatchingChi2MCHMFT", 50.f, "max. chi2 for MCH-MFT matching"}; Configurable maxChi2SA{"maxChi2SA", 1e+6, "max. chi2 for standalone muon"}; - Configurable maxChi2GL{"maxChi2GL", 1e+6, "max. chi2 for global muon"}; + Configurable maxChi2GL{"maxChi2GL", 10, "max. chi2 for global muon"}; Configurable refitGlobalMuon{"refitGlobalMuon", true, "flag to refit global muon"}; + Configurable matchingZ{"matchingZ", -77.5, "z position where matching is performed"}; o2::ccdb::CcdbApi ccdbApi; Service ccdb; - int mRunNumber; + int mRunNumber = 0; + float mBz = 0; HistogramRegistry fRegistry{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; static constexpr std::string_view muon_types[5] = {"MFTMCHMID/", "MFTMCHMIDOtherMatch/", "MFTMCH/", "MCHMID/", "MCH/"}; @@ -108,6 +110,7 @@ struct skimmerPrimaryMuon { addHistograms(); } mRunNumber = 0; + mBz = 0; } void initCCDB(aod::BCsWithTimestamps::iterator const& bc) @@ -126,6 +129,10 @@ struct skimmerPrimaryMuon { ccdb->get(geoPath); } o2::mch::TrackExtrap::setField(); + const double centerMFT[3] = {0, 0, -61.4}; + o2::field::MagneticField* field = static_cast(TGeoGlobalMagField::Instance()->GetField()); + mBz = field->getBz(centerMFT); // Get field at centre of MFT + LOGF(info, "Bz at center of MFT = %f kZG", mBz); } void addHistograms() @@ -143,6 +150,8 @@ struct skimmerPrimaryMuon { fRegistry.add("MFTMCHMID/hDeltaPt_Pt", "#Deltap_{T}/p_{T} vs. p_{T};p_{T}^{gl} (GeV/c);(p_{T}^{sa} - p_{T}^{gl})/p_{T}^{gl}", kTH2F, {{100, 0, 10}, {200, -0.5, +0.5}}, false); fRegistry.add("MFTMCHMID/hDeltaEta_Pt", "#Delta#eta vs. p_{T};p_{T}^{gl} (GeV/c);#Delta#eta", kTH2F, {{100, 0, 10}, {200, -0.5, +0.5}}, false); fRegistry.add("MFTMCHMID/hDeltaPhi_Pt", "#Delta#varphi vs. p_{T};p_{T}^{gl} (GeV/c);#Delta#varphi (rad.)", kTH2F, {{100, 0, 10}, {200, -0.5, +0.5}}, false); + fRegistry.add("MFTMCHMID/hDeltaEtaAtMP_Pt", "#Delta#eta vs. p_{T} at MP;p_{T}^{gl} (GeV/c);#Delta#eta", kTH2F, {{100, 0, 10}, {200, -0.5, +0.5}}, false); + fRegistry.add("MFTMCHMID/hDeltaPhiAtMP_Pt", "#Delta#varphi vs. p_{T} at MP;p_{T}^{gl} (GeV/c);#Delta#varphi (rad.)", kTH2F, {{100, 0, 10}, {200, -0.5, +0.5}}, false); fRegistry.add("MFTMCHMID/hSign", "sign;sign", kTH1F, {{3, -1.5, +1.5}}, false); fRegistry.add("MFTMCHMID/hNclusters", "Nclusters;Nclusters", kTH1F, {{21, -0.5f, 20.5}}, false); fRegistry.add("MFTMCHMID/hNclustersMFT", "NclustersMFT;Nclusters MFT", kTH1F, {{11, -0.5f, 10.5}}, false); @@ -204,8 +213,8 @@ struct skimmerPrimaryMuon { return true; } - template - void fillFwdTrackTable(TCollision const& collision, TFwdTrack fwdtrack, const bool isAmbiguous) + template + void fillFwdTrackTable(TCollision const& collision, TFwdTrack fwdtrack, TMFTTracksCov const& mftCovs, const bool isAmbiguous) { if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack && fwdtrack.chi2MatchMCHMFT() > maxMatchingChi2MCHMFT) { return; @@ -219,13 +228,13 @@ struct skimmerPrimaryMuon { return; } - o2::dataformats::GlobalFwdTrack propmuonAtPV = propagateMuon(fwdtrack, collision, propagationPoint::kToVertex); + o2::dataformats::GlobalFwdTrack propmuonAtPV = propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToVertex, matchingZ, mBz); float pt = propmuonAtPV.getPt(); float eta = propmuonAtPV.getEta(); float phi = propmuonAtPV.getPhi(); o2::math_utils::bringTo02Pi(phi); - o2::dataformats::GlobalFwdTrack propmuonAtDCA = propagateMuon(fwdtrack, collision, propagationPoint::kToDCA); + o2::dataformats::GlobalFwdTrack propmuonAtDCA = propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToDCA, matchingZ, mBz); float cXXatDCA = propmuonAtDCA.getSigma2X(); float cYYatDCA = propmuonAtDCA.getSigma2Y(); float cXYatDCA = propmuonAtDCA.getSigmaXY(); @@ -257,6 +266,12 @@ struct skimmerPrimaryMuon { float chi2mft = 0.f; uint64_t mftClusterSizesAndTrackFlags = 0; int ndf_mchmft = 1; + int ndf_mft = 1; + + float etaMatchedMCHMIDatMP = 999.f; + float phiMatchedMCHMIDatMP = 999.f; + float etaMatchedMFTatMP = 999.f; + float phiMatchedMFTatMP = 999.f; if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { // apply r-absorber cut here to minimize the number of calling propagateMuon. @@ -274,6 +289,7 @@ struct skimmerPrimaryMuon { nClustersMFT = mfttrack.nClusters(); mftClusterSizesAndTrackFlags = mfttrack.mftClusterSizesAndTrackFlags(); ndf_mchmft = 2.f * (mchtrack.nClusters() + nClustersMFT) - 5.f; + ndf_mft = 2.f * nClustersMFT - 5.f; chi2mft = mfttrack.chi2(); // chi2mft = mfttrack.chi2() / (2.f * nClustersMFT - 5.f); @@ -282,13 +298,13 @@ struct skimmerPrimaryMuon { return; } - o2::dataformats::GlobalFwdTrack propmuonAtPV_Matched = propagateMuon(mchtrack, collision, propagationPoint::kToVertex); + o2::dataformats::GlobalFwdTrack propmuonAtPV_Matched = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToVertex, matchingZ, mBz); ptMatchedMCHMID = propmuonAtPV_Matched.getPt(); etaMatchedMCHMID = propmuonAtPV_Matched.getEta(); phiMatchedMCHMID = propmuonAtPV_Matched.getPhi(); o2::math_utils::bringTo02Pi(phiMatchedMCHMID); - o2::dataformats::GlobalFwdTrack propmuonAtDCA_Matched = propagateMuon(mchtrack, collision, propagationPoint::kToDCA); + o2::dataformats::GlobalFwdTrack propmuonAtDCA_Matched = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToDCA, matchingZ, mBz); float dcaX_Matched = propmuonAtDCA_Matched.getX() - collision.posX(); float dcaY_Matched = propmuonAtDCA_Matched.getY() - collision.posY(); float dcaXY_Matched = std::sqrt(dcaX_Matched * dcaX_Matched + dcaY_Matched * dcaY_Matched); @@ -306,9 +322,35 @@ struct skimmerPrimaryMuon { // y = mfttrack.y(); // z = mfttrack.z(); // tgl = mfttrack.tgl(); + + if constexpr (withMFTCov) { + auto mfttrackcov = mftCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]); + auto muonAtMP = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToMatchingPlane, matchingZ, mBz); // propagated to matching plane + o2::track::TrackParCovFwd mftsaAtMP = getTrackParCovFwd(mfttrack, mfttrackcov); // values at innermost update + mftsaAtMP.propagateToZhelix(matchingZ, mBz); // propagated to matching plane + etaMatchedMFTatMP = mftsaAtMP.getEta(); + phiMatchedMFTatMP = mftsaAtMP.getPhi(); + etaMatchedMCHMIDatMP = muonAtMP.getEta(); + phiMatchedMCHMIDatMP = muonAtMP.getPhi(); + o2::math_utils::bringTo02Pi(phiMatchedMCHMIDatMP); + o2::math_utils::bringTo02Pi(phiMatchedMFTatMP); + + // o2::track::TrackParCovFwd mftsa = getTrackParCovFwd(mfttrack, mfttrackcov); // values at innermost update + // o2::dataformats::GlobalFwdTrack globalMuonRefit = o2::aod::fwdtrackutils::refitGlobalMuonCov(propmuonAtPV_Matched, mftsa); // this is track at IU. + // auto globalMuonAtDCA = o2::aod::fwdtrackutils::propagateTrackParCovFwd(globalMuonRefit, fwdtrack.trackType(), collision, propagationPoint::kToDCA, matchingZ, mBz); + // pt = globalMuonAtDCA.getPt(); + // eta = globalMuonAtDCA.getEta(); + // phi = globalMuonAtDCA.getPhi(); + // o2::math_utils::bringTo02Pi(phi); + // cXXatDCA = globalMuonAtDCA.getSigma2X(); + // cYYatDCA = globalMuonAtDCA.getSigma2Y(); + // cXYatDCA = globalMuonAtDCA.getSigmaXY(); + // dcaX = globalMuonAtDCA.getX() - collision.posX(); + // dcaY = globalMuonAtDCA.getY() - collision.posY(); + } } } else if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { - o2::dataformats::GlobalFwdTrack propmuonAtRabs = propagateMuon(fwdtrack, collision, propagationPoint::kToRabs); // this is necessary only for MuonStandaloneTrack + o2::dataformats::GlobalFwdTrack propmuonAtRabs = propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToRabs, matchingZ, mBz); // this is necessary only for MuonStandaloneTrack float xAbs = propmuonAtRabs.getX(); float yAbs = propmuonAtRabs.getY(); rAtAbsorberEnd = std::sqrt(xAbs * xAbs + yAbs * yAbs); // Redo propagation only for muon tracks // propagation of MFT tracks alredy done in reconstruction @@ -325,12 +367,16 @@ struct skimmerPrimaryMuon { float dphi = phiMatchedMCHMID - phi; o2::math_utils::bringToPMPi(dphi); + float detaMP = etaMatchedMCHMIDatMP - etaMatchedMFTatMP; + float dphiMP = phiMatchedMCHMIDatMP - phiMatchedMFTatMP; + o2::math_utils::bringToPMPi(dphiMP); + bool isAssociatedToMPC = fwdtrack.collisionId() == collision.globalIndex(); // LOGF(info, "isAmbiguous = %d, isAssociatedToMPC = %d, fwdtrack.globalIndex() = %d, fwdtrack.collisionId() = %d, collision.globalIndex() = %d", isAmbiguous, isAssociatedToMPC, fwdtrack.globalIndex(), fwdtrack.collisionId(), collision.globalIndex()); emprimarymuons(collision.globalIndex(), fwdtrack.globalIndex(), fwdtrack.matchMFTTrackId(), fwdtrack.matchMCHTrackId(), fwdtrack.trackType(), pt, eta, phi, fwdtrack.sign(), dcaX, dcaY, cXXatDCA, cYYatDCA, cXYatDCA, ptMatchedMCHMID, etaMatchedMCHMID, phiMatchedMCHMID, - // x, y, z, tgl, + etaMatchedMCHMIDatMP, phiMatchedMCHMIDatMP, etaMatchedMFTatMP, phiMatchedMFTatMP, fwdtrack.nClusters(), pDCA, rAtAbsorberEnd, fwdtrack.chi2(), fwdtrack.chi2MatchMCHMID(), fwdtrack.chi2MatchMCHMFT(), fwdtrack.mchBitMap(), fwdtrack.midBitMap(), fwdtrack.midBoards(), mftClusterSizesAndTrackFlags, chi2mft, isAssociatedToMPC, isAmbiguous); @@ -359,13 +405,15 @@ struct skimmerPrimaryMuon { fRegistry.fill(HIST("MFTMCHMID/hDeltaPt_Pt"), pt, dpt); fRegistry.fill(HIST("MFTMCHMID/hDeltaEta_Pt"), pt, deta); fRegistry.fill(HIST("MFTMCHMID/hDeltaPhi_Pt"), pt, dphi); + fRegistry.fill(HIST("MFTMCHMID/hDeltaEtaAtMP_Pt"), pt, detaMP); + fRegistry.fill(HIST("MFTMCHMID/hDeltaPhiAtMP_Pt"), pt, dphiMP); fRegistry.fill(HIST("MFTMCHMID/hSign"), fwdtrack.sign()); fRegistry.fill(HIST("MFTMCHMID/hNclusters"), fwdtrack.nClusters()); fRegistry.fill(HIST("MFTMCHMID/hNclustersMFT"), nClustersMFT); fRegistry.fill(HIST("MFTMCHMID/hPDCA_Rabs"), rAtAbsorberEnd, pDCA); fRegistry.fill(HIST("MFTMCHMID/hRatAbsorberEnd"), rAtAbsorberEnd); fRegistry.fill(HIST("MFTMCHMID/hChi2"), fwdtrack.chi2() / ndf_mchmft); - fRegistry.fill(HIST("MFTMCHMID/hChi2MFT"), chi2mft); + fRegistry.fill(HIST("MFTMCHMID/hChi2MFT"), chi2mft / ndf_mft); fRegistry.fill(HIST("MFTMCHMID/hChi2MatchMCHMID"), fwdtrack.chi2MatchMCHMID()); fRegistry.fill(HIST("MFTMCHMID/hChi2MatchMCHMFT"), fwdtrack.chi2MatchMCHMFT()); fRegistry.fill(HIST("MFTMCHMID/hDCAxy2D"), dcaX, dcaY); @@ -382,13 +430,15 @@ struct skimmerPrimaryMuon { fRegistry.fill(HIST("MCHMID/hDeltaPt_Pt"), pt, dpt); fRegistry.fill(HIST("MCHMID/hDeltaEta_Pt"), pt, deta); fRegistry.fill(HIST("MCHMID/hDeltaPhi_Pt"), pt, dphi); + fRegistry.fill(HIST("MCHMID/hDeltaEtaAtMP_Pt"), pt, detaMP); + fRegistry.fill(HIST("MCHMID/hDeltaPhiAtMP_Pt"), pt, dphiMP); fRegistry.fill(HIST("MCHMID/hSign"), fwdtrack.sign()); fRegistry.fill(HIST("MCHMID/hNclusters"), fwdtrack.nClusters()); fRegistry.fill(HIST("MCHMID/hNclustersMFT"), nClustersMFT); fRegistry.fill(HIST("MCHMID/hPDCA_Rabs"), rAtAbsorberEnd, pDCA); fRegistry.fill(HIST("MCHMID/hRatAbsorberEnd"), rAtAbsorberEnd); fRegistry.fill(HIST("MCHMID/hChi2"), fwdtrack.chi2()); - fRegistry.fill(HIST("MCHMID/hChi2MFT"), chi2mft); + fRegistry.fill(HIST("MCHMID/hChi2MFT"), chi2mft / ndf_mft); fRegistry.fill(HIST("MCHMID/hChi2MatchMCHMID"), fwdtrack.chi2MatchMCHMID()); fRegistry.fill(HIST("MCHMID/hChi2MatchMCHMFT"), fwdtrack.chi2MatchMCHMFT()); fRegistry.fill(HIST("MCHMID/hDCAxy2D"), dcaX, dcaY); @@ -402,7 +452,7 @@ struct skimmerPrimaryMuon { } } - // std::map, float> mCandidates; // std::pair -> chi2MatchMCHMFT; + std::unordered_map map_mfttrackcovs; std::vector> vec_min_chi2MatchMCHMFT; // std::pair -> chi2MatchMCHMFT; template void findBestMatchPerMCHMID(TMuons const& muons) @@ -426,7 +476,6 @@ struct skimmerPrimaryMuon { } } vec_min_chi2MatchMCHMFT.emplace_back(tupleIds_at_min); - // mCandidates[tupleIds_at_min] = min_chi2MatchMCHMFT; // LOGF(info, "min: muon_tmp.globalIndex() = %d, muon_tmp.matchMCHTrackId() = %d, muon_tmp.matchMFTTrackId() = %d, muon_tmp.chi2MatchMCHMFT() = %f", std::get<0>(tupleIds_at_min), std::get<1>(tupleIds_at_min), std::get<2>(tupleIds_at_min), min_chi2MatchMCHMFT); } } // end of muon loop @@ -460,10 +509,11 @@ struct skimmerPrimaryMuon { continue; } - fillFwdTrackTable(collision, fwdtrack, false); + fillFwdTrackTable(collision, fwdtrack, nullptr, false); } // end of fwdtrack loop } // end of collision loop + map_mfttrackcovs.clear(); vec_min_chi2MatchMCHMFT.clear(); vec_min_chi2MatchMCHMFT.shrink_to_fit(); } @@ -498,15 +548,59 @@ struct skimmerPrimaryMuon { continue; } - fillFwdTrackTable(collision, fwdtrack, mapAmb[fwdtrack.globalIndex()]); + fillFwdTrackTable(collision, fwdtrack, nullptr, mapAmb[fwdtrack.globalIndex()]); } // end of fwdtrack loop } // end of collision loop + mapAmb.clear(); + map_mfttrackcovs.clear(); vec_min_chi2MatchMCHMFT.clear(); vec_min_chi2MatchMCHMFT.shrink_to_fit(); } PROCESS_SWITCH(skimmerPrimaryMuon, processRec_TTCA, "process reconstructed info", false); + void processRec_TTCA_withMFTCov(MyCollisions const& collisions, MyFwdTracks const& fwdtracks, aod::MFTTracks const&, aod::BCsWithTimestamps const&, aod::FwdTrackAssoc const& fwdtrackIndices, aod::MFTTracksCov const& mftCovs) + { + for (const auto& mfttrackConv : mftCovs) { + map_mfttrackcovs[mfttrackConv.matchMFTTrackId()] = mfttrackConv.globalIndex(); + } + findBestMatchPerMCHMID(fwdtracks); + + std::unordered_map mapAmb; // fwdtrack.globalIndex() -> bool isAmb; + for (const auto& fwdtrack : fwdtracks) { + const auto& fwdtrackIdsPerFwdTrack = fwdtrackIndices.sliceBy(fwdtrackIndicesPerFwdTrack, fwdtrack.globalIndex()); + mapAmb[fwdtrack.globalIndex()] = fwdtrackIdsPerFwdTrack.size() > 1; + // LOGF(info, "fwdtrack.globalIndex() = %d, ntimes = %d, isAmbiguous = %d", fwdtrack.globalIndex(), fwdtrackIdsPerFwdTrack.size(), mapAmb[fwdtrack.globalIndex()]); + } // end of fwdtrack loop + + for (const auto& collision : collisions) { + const auto& bc = collision.template bc_as(); + initCCDB(bc); + + if (!collision.isSelected()) { + continue; + } + + const auto& fwdtrackIdsThisCollision = fwdtrackIndices.sliceBy(fwdtrackIndicesPerCollision, collision.globalIndex()); + for (const auto& fwdtrackId : fwdtrackIdsThisCollision) { + const auto& fwdtrack = fwdtrackId.template fwdtrack_as(); + if (fwdtrack.trackType() != o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack && fwdtrack.trackType() != o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { + continue; + } + if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack && std::find(vec_min_chi2MatchMCHMFT.begin(), vec_min_chi2MatchMCHMFT.end(), std::make_tuple(fwdtrack.globalIndex(), fwdtrack.matchMCHTrackId(), fwdtrack.matchMFTTrackId())) == vec_min_chi2MatchMCHMFT.end()) { + continue; + } + + fillFwdTrackTable(collision, fwdtrack, mftCovs, mapAmb[fwdtrack.globalIndex()]); + } // end of fwdtrack loop + } // end of collision loop + mapAmb.clear(); + map_mfttrackcovs.clear(); + vec_min_chi2MatchMCHMFT.clear(); + vec_min_chi2MatchMCHMFT.shrink_to_fit(); + } + PROCESS_SWITCH(skimmerPrimaryMuon, processRec_TTCA_withMFTCov, "process reconstructed info", false); + void processRec_SA_SWT(MyCollisionsWithSWT const& collisions, MyFwdTracks const& fwdtracks, aod::MFTTracks const&, aod::BCsWithTimestamps const&) { findBestMatchPerMCHMID(fwdtracks); @@ -532,9 +626,10 @@ struct skimmerPrimaryMuon { continue; } - fillFwdTrackTable(collision, fwdtrack, false); + fillFwdTrackTable(collision, fwdtrack, nullptr, false); } // end of fwdtrack loop } // end of collision loop + map_mfttrackcovs.clear(); vec_min_chi2MatchMCHMFT.clear(); vec_min_chi2MatchMCHMFT.shrink_to_fit(); } @@ -571,15 +666,60 @@ struct skimmerPrimaryMuon { continue; } - fillFwdTrackTable(collision, fwdtrack, mapAmb[fwdtrack.globalIndex()]); + fillFwdTrackTable(collision, fwdtrack, nullptr, mapAmb[fwdtrack.globalIndex()]); } // end of fwdtrack loop } // end of collision loop mapAmb.clear(); + map_mfttrackcovs.clear(); vec_min_chi2MatchMCHMFT.clear(); vec_min_chi2MatchMCHMFT.shrink_to_fit(); } PROCESS_SWITCH(skimmerPrimaryMuon, processRec_TTCA_SWT, "process reconstructed info", false); + void processRec_TTCA_SWT_withMFTCov(MyCollisionsWithSWT const& collisions, MyFwdTracks const& fwdtracks, aod::MFTTracks const&, aod::BCsWithTimestamps const&, aod::FwdTrackAssoc const& fwdtrackIndices, aod::MFTTracksCov const& mftCovs) + { + for (const auto& mfttrackConv : mftCovs) { + map_mfttrackcovs[mfttrackConv.matchMFTTrackId()] = mfttrackConv.globalIndex(); + } + findBestMatchPerMCHMID(fwdtracks); + + std::unordered_map mapAmb; // fwdtrack.globalIndex() -> bool isAmb; + for (const auto& fwdtrack : fwdtracks) { + const auto& fwdtrackIdsPerFwdTrack = fwdtrackIndices.sliceBy(fwdtrackIndicesPerFwdTrack, fwdtrack.globalIndex()); + mapAmb[fwdtrack.globalIndex()] = fwdtrackIdsPerFwdTrack.size() > 1; + // LOGF(info, "fwdtrack.globalIndex() = %d, ntimes = %d, isAmbiguous = %d", fwdtrack.globalIndex(), fwdtrackIdsPerFwdTrack.size(), mapAmb[fwdtrack.globalIndex()]); + } // end of fwdtrack loop + + for (const auto& collision : collisions) { + const auto& bc = collision.template bc_as(); + initCCDB(bc); + if (!collision.isSelected()) { + continue; + } + if (collision.swtaliastmp_raw() == 0) { + continue; + } + + const auto& fwdtrackIdsThisCollision = fwdtrackIndices.sliceBy(fwdtrackIndicesPerCollision, collision.globalIndex()); + for (const auto& fwdtrackId : fwdtrackIdsThisCollision) { + const auto& fwdtrack = fwdtrackId.template fwdtrack_as(); + if (fwdtrack.trackType() != o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack && fwdtrack.trackType() != o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { + continue; + } + if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack && std::find(vec_min_chi2MatchMCHMFT.begin(), vec_min_chi2MatchMCHMFT.end(), std::make_tuple(fwdtrack.globalIndex(), fwdtrack.matchMCHTrackId(), fwdtrack.matchMFTTrackId())) == vec_min_chi2MatchMCHMFT.end()) { + continue; + } + + fillFwdTrackTable(collision, fwdtrack, mftCovs, mapAmb[fwdtrack.globalIndex()]); + } // end of fwdtrack loop + } // end of collision loop + mapAmb.clear(); + map_mfttrackcovs.clear(); + vec_min_chi2MatchMCHMFT.clear(); + vec_min_chi2MatchMCHMFT.shrink_to_fit(); + } + PROCESS_SWITCH(skimmerPrimaryMuon, processRec_TTCA_SWT_withMFTCov, "process reconstructed info", false); + void processMC_SA(soa::Join const& collisions, MyFwdTracksMC const& fwdtracks, MFTTracksMC const&, aod::BCsWithTimestamps const&) { findBestMatchPerMCHMID(fwdtracks); @@ -606,9 +746,10 @@ struct skimmerPrimaryMuon { continue; } - fillFwdTrackTable(collision, fwdtrack, false); + fillFwdTrackTable(collision, fwdtrack, nullptr, false); } // end of fwdtrack loop } // end of collision loop + map_mfttrackcovs.clear(); vec_min_chi2MatchMCHMFT.clear(); vec_min_chi2MatchMCHMFT.shrink_to_fit(); } @@ -648,15 +789,63 @@ struct skimmerPrimaryMuon { continue; } - fillFwdTrackTable(collision, fwdtrack, mapAmb[fwdtrack.globalIndex()]); + fillFwdTrackTable(collision, fwdtrack, nullptr, mapAmb[fwdtrack.globalIndex()]); } // end of fwdtrack loop } // end of collision loop mapAmb.clear(); + map_mfttrackcovs.clear(); vec_min_chi2MatchMCHMFT.clear(); vec_min_chi2MatchMCHMFT.shrink_to_fit(); } PROCESS_SWITCH(skimmerPrimaryMuon, processMC_TTCA, "process reconstructed and MC info", false); + void processMC_TTCA_withMFTCov(soa::Join const& collisions, MyFwdTracksMC const& fwdtracks, MFTTracksMC const&, aod::BCsWithTimestamps const&, aod::FwdTrackAssoc const& fwdtrackIndices, aod::MFTTracksCov const& mftCovs) + { + for (const auto& mfttrackConv : mftCovs) { + map_mfttrackcovs[mfttrackConv.matchMFTTrackId()] = mfttrackConv.globalIndex(); + } + findBestMatchPerMCHMID(fwdtracks); + + std::unordered_map mapAmb; // fwdtrack.globalIndex() -> bool isAmb; + for (const auto& fwdtrack : fwdtracks) { + const auto& fwdtrackIdsPerFwdTrack = fwdtrackIndices.sliceBy(fwdtrackIndicesPerFwdTrack, fwdtrack.globalIndex()); + mapAmb[fwdtrack.globalIndex()] = fwdtrackIdsPerFwdTrack.size() > 1; + // LOGF(info, "fwdtrack.globalIndex() = %d, ntimes = %d, isAmbiguous = %d", fwdtrack.globalIndex(), fwdtrackIdsPerFwdTrack.size(), mapAmb[fwdtrack.globalIndex()]); + } // end of fwdtrack loop + + for (const auto& collision : collisions) { + const auto& bc = collision.template bc_as(); + initCCDB(bc); + if (!collision.isSelected()) { + continue; + } + if (!collision.has_mcCollision()) { + continue; + } + + const auto& fwdtrackIdsThisCollision = fwdtrackIndices.sliceBy(fwdtrackIndicesPerCollision, collision.globalIndex()); + for (const auto& fwdtrackId : fwdtrackIdsThisCollision) { + const auto& fwdtrack = fwdtrackId.template fwdtrack_as(); + if (!fwdtrack.has_mcParticle()) { + continue; + } + if (fwdtrack.trackType() != o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack && fwdtrack.trackType() != o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { + continue; + } + if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack && std::find(vec_min_chi2MatchMCHMFT.begin(), vec_min_chi2MatchMCHMFT.end(), std::make_tuple(fwdtrack.globalIndex(), fwdtrack.matchMCHTrackId(), fwdtrack.matchMFTTrackId())) == vec_min_chi2MatchMCHMFT.end()) { + continue; + } + + fillFwdTrackTable(collision, fwdtrack, mftCovs, mapAmb[fwdtrack.globalIndex()]); + } // end of fwdtrack loop + } // end of collision loop + mapAmb.clear(); + map_mfttrackcovs.clear(); + vec_min_chi2MatchMCHMFT.clear(); + vec_min_chi2MatchMCHMFT.shrink_to_fit(); + } + PROCESS_SWITCH(skimmerPrimaryMuon, processMC_TTCA_withMFTCov, "process reconstructed and MC with MFTCov info", false); + void processDummy(aod::Collisions const&) {} PROCESS_SWITCH(skimmerPrimaryMuon, processDummy, "process dummy", true); }; diff --git a/PWGEM/Dilepton/Tasks/Converters/CMakeLists.txt b/PWGEM/Dilepton/Tasks/Converters/CMakeLists.txt index b41943ef7ff..489839883ec 100644 --- a/PWGEM/Dilepton/Tasks/Converters/CMakeLists.txt +++ b/PWGEM/Dilepton/Tasks/Converters/CMakeLists.txt @@ -45,3 +45,8 @@ o2physics_add_dpl_workflow(track-converter1 PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(muon-converter1 + SOURCES muonConverter1.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + diff --git a/PWGEM/Dilepton/Tasks/Converters/muonConverter1.cxx b/PWGEM/Dilepton/Tasks/Converters/muonConverter1.cxx new file mode 100644 index 00000000000..527faf7c723 --- /dev/null +++ b/PWGEM/Dilepton/Tasks/Converters/muonConverter1.cxx @@ -0,0 +1,53 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +// +// ======================== +// +// This code produces muon table 001 from 000. +// Please write to: daiki.sekihata@cern.ch + +#include "PWGEM/Dilepton/DataModel/dileptonTables.h" + +#include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisTask.h" +#include "Framework/runDataProcessing.h" + +using namespace o2; +using namespace o2::aod; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::soa; + +struct muonConverter1 { + Produces muon_001; + + void process(aod::EMPrimaryMuons_000 const& muons) + { + for (const auto& muon : muons) { + muon_001( + muon.collisionId(), + muon.fwdtrackId(), muon.mfttrackId(), muon.mchtrackId(), muon.trackType(), + muon.pt(), muon.eta(), muon.phi(), muon.sign(), + muon.fwdDcaX(), muon.fwdDcaY(), muon.cXXatDCA(), muon.cYYatDCA(), muon.cXYatDCA(), + muon.ptMatchedMCHMID(), muon.etaMatchedMCHMID(), muon.phiMatchedMCHMID(), + 0, 0, 0, 0, + muon.nClusters(), muon.pDca(), muon.rAtAbsorberEnd(), + muon.chi2(), muon.chi2MatchMCHMID(), muon.chi2MatchMCHMFT(), + muon.mchBitMap(), muon.midBitMap(), muon.midBoards(), + muon.mftClusterSizesAndTrackFlags(), muon.chi2MFT(), muon.isAssociatedToMPC(), muon.isAmbiguous()); + } // end of muon loop + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc, TaskName{"muon-converter1"})}; +} diff --git a/PWGEM/Dilepton/Tasks/Converters/trackConverter1.cxx b/PWGEM/Dilepton/Tasks/Converters/trackConverter1.cxx index fda87265b3c..c6f7bfce634 100644 --- a/PWGEM/Dilepton/Tasks/Converters/trackConverter1.cxx +++ b/PWGEM/Dilepton/Tasks/Converters/trackConverter1.cxx @@ -31,7 +31,7 @@ struct trackConverter1 { void process(aod::EMPrimaryTracks_000 const& tracks) { - for (auto& track : tracks) { + for (const auto& track : tracks) { track_001(track.collisionId(), track.trackId(), track.sign() / track.pt(), diff --git a/PWGEM/Dilepton/Tasks/createResolutionMap.cxx b/PWGEM/Dilepton/Tasks/createResolutionMap.cxx index 09dc09f37a0..cc93bbccde9 100644 --- a/PWGEM/Dilepton/Tasks/createResolutionMap.cxx +++ b/PWGEM/Dilepton/Tasks/createResolutionMap.cxx @@ -16,6 +16,7 @@ #include "PWGEM/Dilepton/Utils/MCUtilities.h" #include "Common/CCDB/RCTSelectionFlags.h" +#include "Common/Core/RecoDecay.h" #include "Common/Core/fwdtrackUtilities.h" #include "Common/Core/trackUtilities.h" #include "Common/DataModel/Centrality.h" @@ -46,6 +47,8 @@ #include #include #include +#include +#include #include #include @@ -167,9 +170,12 @@ struct CreateResolutionMap { Configurable cfg_max_reldpt{"cfg_max_reldpt", 1e+10f, "max. relative dpt between MFT-MCH-MID and MCH-MID"}; Configurable cfg_max_deta{"cfg_max_deta", 1e+10f, "max. deta between MFT-MCH-MID and MCH-MID"}; Configurable cfg_max_dphi{"cfg_max_dphi", 1e+10f, "max. dphi between MFT-MCH-MID and MCH-MID"}; + Configurable cfg_max_detaMP{"cfg_max_detaMP", 1e+10f, "max. deta between MFT and MCH-MID at matching plane"}; + Configurable cfg_max_dphiMP{"cfg_max_dphiMP", 1e+10f, "max. dphi between MFT and MCH-MID at matching plane"}; Configurable refitGlobalMuon{"refitGlobalMuon", true, "flag to refit global muon"}; Configurable requireMFTHitMap{"requireMFTHitMap", false, "flag to require MFT hit map"}; Configurable> requiredMFTDisks{"requiredMFTDisks", std::vector{4}, "hit map on MFT disks [0,1,2,3,4]. logical-OR of each double-sided disk"}; + Configurable matchingZ{"matchingZ", -77.5, "z position where matching is performed"}; } muoncuts; HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject}; @@ -180,7 +186,8 @@ struct CreateResolutionMap { o2::ccdb::CcdbApi ccdbApi; Service ccdb; int mRunNumber = 0; - float d_bz; + float mBzMFT = 0; + float d_bz = 0; // o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE; o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrLUT; o2::dataformats::VertexBase mVtx; @@ -215,6 +222,7 @@ struct CreateResolutionMap { mRunNumber = 0; d_bz = 0; + mBzMFT = 0; if (ConfPhiPositionCBGenBins.value[0] == VARIABLE_WIDTH) { phiPosition_bin_edges = std::vector(ConfPhiPositionCBGenBins.value.begin(), ConfPhiPositionCBGenBins.value.end()); @@ -305,7 +313,10 @@ struct CreateResolutionMap { ccdb->get(geoPath); } o2::mch::TrackExtrap::setField(); - return; + const double centerMFT[3] = {0, 0, -61.4}; + o2::field::MagneticField* field = static_cast(TGeoGlobalMagField::Instance()->GetField()); + mBzMFT = field->getBz(centerMFT); // Get field at centre of MFT + LOGF(info, "Bz at center of MFT = %f kZG", mBzMFT); } auto run3grp_timestamp = bc.timestamp(); @@ -518,8 +529,8 @@ struct CreateResolutionMap { return true; } - template - void fillMuon(TCollision const& collision, TMuon const& muon, const float centrality) + template + void fillMuon(TCollision const& collision, TMuon const& muon, TMFTTracksCov const& mftCovs, const float centrality) { auto mcparticle = muon.template mcParticle_as(); if (std::abs(mcparticle.pdgCode()) != 13 || !(mcparticle.isPhysicalPrimary() || mcparticle.producedByGenerator())) { @@ -532,11 +543,15 @@ struct CreateResolutionMap { return; } + if (muon.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack && std::find(vec_min_chi2MatchMCHMFT.begin(), vec_min_chi2MatchMCHMFT.end(), std::make_tuple(muon.globalIndex(), muon.matchMCHTrackId(), muon.matchMFTTrackId())) == vec_min_chi2MatchMCHMFT.end()) { + return; + } + if (muon.chi2MatchMCHMID() < 0.f) { // this should never happen. only for protection. return; } - o2::dataformats::GlobalFwdTrack propmuonAtPV = propagateMuon(muon, collision, propagationPoint::kToVertex); - o2::dataformats::GlobalFwdTrack propmuonAtDCA = propagateMuon(muon, collision, propagationPoint::kToDCA); + o2::dataformats::GlobalFwdTrack propmuonAtPV = propagateMuon(muon, muon, collision, propagationPoint::kToVertex, muoncuts.matchingZ, mBzMFT); + o2::dataformats::GlobalFwdTrack propmuonAtDCA = propagateMuon(muon, muon, collision, propagationPoint::kToDCA, muoncuts.matchingZ, mBzMFT); float pt = propmuonAtPV.getPt(); float eta = propmuonAtPV.getEta(); @@ -555,10 +570,15 @@ struct CreateResolutionMap { float phiMatchedMCHMID = propmuonAtPV.getPhi(); o2::math_utils::bringTo02Pi(phiMatchedMCHMID); + float etaMatchedMCHMIDatMP = 999.f; + float phiMatchedMCHMIDatMP = 999.f; + float etaMatchedMFTatMP = 999.f; + float phiMatchedMFTatMP = 999.f; + if (muon.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { // mcparticle for global MFT-MCH-MID is identical to mcparticle of MCH-MID track. If not, mismatch. - const auto& mchtrack = muon.template matchMCHTrack_as(); // MCH-MID - const auto& mfttrack = muon.template matchMFTTrack_as(); // MFTsa + auto mchtrack = muon.template matchMCHTrack_as(); // MCH-MID + auto mfttrack = muon.template matchMFTTrack_as(); // MFTsa if (!mchtrack.has_mcParticle() || !mfttrack.has_mcParticle()) { return; } @@ -571,12 +591,12 @@ struct CreateResolutionMap { return; } - o2::dataformats::GlobalFwdTrack propmuonAtPV_Matched = propagateMuon(mchtrack, collision, propagationPoint::kToVertex); + o2::dataformats::GlobalFwdTrack propmuonAtPV_Matched = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToVertex, muoncuts.matchingZ, mBzMFT); ptMatchedMCHMID = propmuonAtPV_Matched.getPt(); etaMatchedMCHMID = propmuonAtPV_Matched.getEta(); phiMatchedMCHMID = propmuonAtPV_Matched.getPhi(); o2::math_utils::bringTo02Pi(phiMatchedMCHMID); - o2::dataformats::GlobalFwdTrack propmuonAtDCA_Matched = propagateMuon(mchtrack, collision, propagationPoint::kToDCA); + o2::dataformats::GlobalFwdTrack propmuonAtDCA_Matched = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToDCA, muoncuts.matchingZ, mBzMFT); float dcaX_Matched = propmuonAtDCA_Matched.getX() - collision.posX(); float dcaY_Matched = propmuonAtDCA_Matched.getY() - collision.posY(); float dcaXY_Matched = std::sqrt(dcaX_Matched * dcaX_Matched + dcaY_Matched * dcaY_Matched); @@ -589,6 +609,7 @@ struct CreateResolutionMap { if (muon.chi2MatchMCHMFT() > muoncuts.cfg_max_matching_chi2_mftmch) { return; } + if (muoncuts.refitGlobalMuon) { eta = mfttrack.eta(); // phi = mfttrack.phi(); @@ -596,13 +617,49 @@ struct CreateResolutionMap { pt = propmuonAtPV_Matched.getP() * std::sin(2.f * std::atan(std::exp(-eta))); phi = propmuonAtDCA.getPhi(); o2::math_utils::bringTo02Pi(phi); + + if constexpr (withMFTCov) { + auto mfttrackcov = mftCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]); + auto muonAtMP = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToMatchingPlane, muoncuts.matchingZ, mBzMFT); // propagated to matching plane + o2::track::TrackParCovFwd mftsaAtMP = getTrackParCovFwd(mfttrack, mfttrackcov); // values at innermost update + mftsaAtMP.propagateToZhelix(muoncuts.matchingZ, mBzMFT); // propagated to matching plane + etaMatchedMFTatMP = mftsaAtMP.getEta(); + phiMatchedMFTatMP = mftsaAtMP.getPhi(); + etaMatchedMCHMIDatMP = muonAtMP.getEta(); + phiMatchedMCHMIDatMP = muonAtMP.getPhi(); + o2::math_utils::bringTo02Pi(phiMatchedMCHMIDatMP); + o2::math_utils::bringTo02Pi(phiMatchedMFTatMP); + + // o2::track::TrackParCovFwd mftsa = getTrackParCovFwd(mfttrack, mfttrackcov); // values at innermost update + // o2::dataformats::GlobalFwdTrack globalMuonRefit = o2::aod::fwdtrackutils::refitGlobalMuonCov(propmuonAtPV_Matched, mftsa); // this is track at IU. + // auto globalMuonAtDCA = o2::aod::fwdtrackutils::propagateTrackParCovFwd(globalMuonRefit, fwdtrack.trackType(), collision, propagationPoint::kToDCA, matchingZ, mBz); + // pt = globalMuonAtDCA.getPt(); + // eta = globalMuonAtDCA.getEta(); + // phi = globalMuonAtDCA.getPhi(); + // o2::math_utils::bringTo02Pi(phi); + // cXXatDCA = globalMuonAtDCA.getSigma2X(); + // cYYatDCA = globalMuonAtDCA.getSigma2Y(); + // cXYatDCA = globalMuonAtDCA.getSigmaXY(); + // dcaX = globalMuonAtDCA.getX() - collision.posX(); + // dcaY = globalMuonAtDCA.getY() - collision.posY(); + } } float dpt = (ptMatchedMCHMID - pt) / pt; + if (std::fabs(dpt) > muoncuts.cfg_max_reldpt) { + return; + } float deta = etaMatchedMCHMID - eta; float dphi = phiMatchedMCHMID - phi; o2::math_utils::bringToPMPi(dphi); - if (std::sqrt(std::pow(deta / muoncuts.cfg_max_deta, 2) + std::pow(dphi / muoncuts.cfg_max_dphi, 2)) > 1.f || std::fabs(dpt) > muoncuts.cfg_max_reldpt) { + if (std::sqrt(std::pow(deta / muoncuts.cfg_max_deta, 2) + std::pow(dphi / muoncuts.cfg_max_dphi, 2)) > 1.f) { + return; + } + + float detaMP = etaMatchedMCHMIDatMP - etaMatchedMFTatMP; + float dphiMP = phiMatchedMCHMIDatMP - phiMatchedMFTatMP; + o2::math_utils::bringToPMPi(dphiMP); + if (std::sqrt(std::pow(detaMP / muoncuts.cfg_max_detaMP, 2) + std::pow(dphiMP / muoncuts.cfg_max_dphiMP, 2)) > 1.f) { return; } @@ -616,7 +673,7 @@ struct CreateResolutionMap { } } else if (muon.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { - o2::dataformats::GlobalFwdTrack propmuonAtRabs = propagateMuon(muon, collision, propagationPoint::kToRabs); // this is necessary only for MuonStandaloneTrack + o2::dataformats::GlobalFwdTrack propmuonAtRabs = propagateMuon(muon, muon, collision, propagationPoint::kToRabs, muoncuts.matchingZ, mBzMFT); // this is necessary only for MuonStandaloneTrack float xAbs = propmuonAtRabs.getX(); float yAbs = propmuonAtRabs.getY(); rAtAbsorberEnd = std::sqrt(xAbs * xAbs + yAbs * yAbs); // Redo propagation only for muon tracks // propagation of MFT tracks alredy done in reconstruction @@ -729,6 +786,7 @@ struct CreateResolutionMap { SliceCache cache; Preslice perCollision_mid = o2::aod::track::collisionId; Preslice perCollision_fwd = o2::aod::fwdtrack::collisionId; + PresliceUnsorted fwdtracksPerMCHTrack = aod::fwdtrack::matchMCHTrackId; using MyCollisions = Join; using MyCollision = MyCollisions::iterator; @@ -802,6 +860,34 @@ struct CreateResolutionMap { } } + std::vector> vec_min_chi2MatchMCHMFT; // std::pair -> chi2MatchMCHMFT; + template + void findBestMatchPerMCHMID(TMuons const& muons) + { + vec_min_chi2MatchMCHMFT.reserve(muons.size()); + for (const auto& muon : muons) { + if (muon.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { + const auto& muons_per_MCHMID = muons.sliceBy(fwdtracksPerMCHTrack, muon.globalIndex()); + // LOGF(info, "stanadalone: muon.globalIndex() = %d, muon.chi2MatchMCHMFT() = %f", muon.globalIndex(), muon.chi2MatchMCHMFT()); + // LOGF(info, "muons_per_MCHMID.size() = %d", muons_per_MCHMID.size()); + + float min_chi2MatchMCHMFT = 1e+10; + std::tuple tupleIds_at_min; + for (const auto& muon_tmp : muons_per_MCHMID) { + if (muon_tmp.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { + // LOGF(info, "muon_tmp.globalIndex() = %d, muon_tmp.matchMCHTrackId() = %d, muon_tmp.matchMFTTrackId() = %d, muon_tmp.chi2MatchMCHMFT() = %f", muon_tmp.globalIndex(), muon_tmp.matchMCHTrackId(), muon_tmp.matchMFTTrackId(), muon_tmp.chi2MatchMCHMFT()); + if (0.f < muon_tmp.chi2MatchMCHMFT() && muon_tmp.chi2MatchMCHMFT() < min_chi2MatchMCHMFT) { + min_chi2MatchMCHMFT = muon_tmp.chi2MatchMCHMFT(); + tupleIds_at_min = std::make_tuple(muon_tmp.globalIndex(), muon_tmp.matchMCHTrackId(), muon_tmp.matchMFTTrackId()); + } + } + } + vec_min_chi2MatchMCHMFT.emplace_back(tupleIds_at_min); + // LOGF(info, "min: muon_tmp.globalIndex() = %d, muon_tmp.matchMCHTrackId() = %d, muon_tmp.matchMFTTrackId() = %d, muon_tmp.chi2MatchMCHMFT() = %f", std::get<0>(tupleIds_at_min), std::get<1>(tupleIds_at_min), std::get<2>(tupleIds_at_min), min_chi2MatchMCHMFT); + } + } // end of muon loop + } + void processElectronSA(MyCollisions const& collisions, aod::BCsWithTimestamps const&, MyTracks const& tracks, aod::McCollisions const&, aod::McParticles const&) { for (const auto& collision : collisions) { @@ -879,6 +965,7 @@ struct CreateResolutionMap { void processMuonSA(MyCollisions const& collisions, aod::BCsWithTimestamps const&, MyFwdTracks const& fwdtracks, MyMFTTracks const&, aod::McCollisions const&, aod::McParticles const&) { + findBestMatchPerMCHMID(fwdtracks); for (const auto& collision : collisions) { auto bc = collision.template foundBC_as(); initCCDB(bc); @@ -906,16 +993,19 @@ struct CreateResolutionMap { if (cfgEventGeneratorType >= 0 && mccollision_from_mctrack.getSubGeneratorId() != cfgEventGeneratorType) { continue; } - fillMuon(collision, muon, centrality); + fillMuon(collision, muon, nullptr, centrality); } // end of standalone muon loop } // end of collision loop + vec_min_chi2MatchMCHMFT.clear(); + vec_min_chi2MatchMCHMFT.shrink_to_fit(); } PROCESS_SWITCH(CreateResolutionMap, processMuonSA, "create resolution map for muon at forward rapidity", true); Preslice fwdtrackIndicesPerCollision = aod::track_association::collisionId; - void processMuonTTCA(MyCollisions const& collisions, aod::BCsWithTimestamps const&, MyFwdTracks const&, MyMFTTracks const&, aod::FwdTrackAssoc const& fwdtrackIndices, aod::McCollisions const&, aod::McParticles const&) + void processMuonTTCA(MyCollisions const& collisions, aod::BCsWithTimestamps const&, MyFwdTracks const& fwdtracks, MyMFTTracks const&, aod::FwdTrackAssoc const& fwdtrackIndices, aod::McCollisions const&, aod::McParticles const&) { + findBestMatchPerMCHMID(fwdtracks); for (const auto& collision : collisions) { auto bc = collision.template foundBC_as(); initCCDB(bc); @@ -943,12 +1033,58 @@ struct CreateResolutionMap { if (cfgEventGeneratorType >= 0 && mccollision_from_mctrack.getSubGeneratorId() != cfgEventGeneratorType) { continue; } - fillMuon(collision, muon, centrality); + fillMuon(collision, muon, nullptr, centrality); } // end of fwdtrack loop } // end of collision loop + vec_min_chi2MatchMCHMFT.clear(); + vec_min_chi2MatchMCHMFT.shrink_to_fit(); } PROCESS_SWITCH(CreateResolutionMap, processMuonTTCA, "create resolution map for muon at forward rapidity", false); + std::unordered_map map_mfttrackcovs; + void processMuonTTCA_withMFTCov(MyCollisions const& collisions, aod::BCsWithTimestamps const&, MyFwdTracks const& fwdtracks, MyMFTTracks const&, aod::MFTTracksCov const& mftCovs, aod::FwdTrackAssoc const& fwdtrackIndices, aod::McCollisions const&, aod::McParticles const&) + { + for (const auto& mfttrackConv : mftCovs) { + map_mfttrackcovs[mfttrackConv.matchMFTTrackId()] = mfttrackConv.globalIndex(); + } + findBestMatchPerMCHMID(fwdtracks); + + for (const auto& collision : collisions) { + auto bc = collision.template foundBC_as(); + initCCDB(bc); + + if (!isSelectedEvent(collision)) { + continue; + } + + if (!collision.has_mcCollision()) { + continue; + } + + float centrality = std::array{collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}[cfgCentEstimator]; + // auto mccollision = collision.template mcCollision_as(); + // registry.fill(HIST("Event/Muon/hImpPar_Centrality"), mccollision.impactParameter(), centrality); + + auto fwdtrackIdsThisCollision = fwdtrackIndices.sliceBy(fwdtrackIndicesPerCollision, collision.globalIndex()); + for (const auto& fwdtrackId : fwdtrackIdsThisCollision) { + auto muon = fwdtrackId.template fwdtrack_as(); + if (!muon.has_mcParticle()) { + continue; + } + auto mctrack = muon.template mcParticle_as(); + auto mccollision_from_mctrack = mctrack.template mcCollision_as(); + if (cfgEventGeneratorType >= 0 && mccollision_from_mctrack.getSubGeneratorId() != cfgEventGeneratorType) { + continue; + } + fillMuon(collision, muon, mftCovs, centrality); + } // end of fwdtrack loop + } // end of collision loop + map_mfttrackcovs.clear(); + vec_min_chi2MatchMCHMFT.clear(); + vec_min_chi2MatchMCHMFT.shrink_to_fit(); + } + PROCESS_SWITCH(CreateResolutionMap, processMuonTTCA_withMFTCov, "create resolution map for muon at forward rapidity", false); + void processGen(aod::McCollisions const& mcCollisions) { for (const auto& mccollision : mcCollisions) { diff --git a/PWGEM/Dilepton/Tasks/matchingMFT.cxx b/PWGEM/Dilepton/Tasks/matchingMFT.cxx index e9f6dadc28f..9a8d492996f 100644 --- a/PWGEM/Dilepton/Tasks/matchingMFT.cxx +++ b/PWGEM/Dilepton/Tasks/matchingMFT.cxx @@ -16,6 +16,7 @@ #include "TableHelper.h" #include "Common/CCDB/RCTSelectionFlags.h" +#include "Common/Core/RecoDecay.h" #include "Common/Core/fwdtrackUtilities.h" #include "Common/DataModel/Centrality.h" #include "Common/DataModel/CollisionAssociationTables.h" @@ -36,6 +37,7 @@ #include "TGeoGlobalMagField.h" +#include #include #include #include @@ -79,8 +81,12 @@ struct matchingMFT { Configurable maxRelDPt{"maxRelDPt", 1e+10f, "max. relative dpt between MFT-MCH-MID and MCH-MID"}; Configurable maxDEta{"maxDEta", 1e+10f, "max. deta between MFT-MCH-MID and MCH-MID"}; Configurable maxDPhi{"maxDPhi", 1e+10f, "max. dphi between MFT-MCH-MID and MCH-MID"}; + Configurable maxDEtaMP{"maxDEtaMP", 1e+10f, "max. deta between MFT and MCH-MID at matching plane Z"}; + Configurable maxDPhiMP{"maxDPhiMP", 1e+10f, "max. dphi between MFT and MCH-MID at matching plane Z"}; Configurable requireMFTHitMap{"requireMFTHitMap", false, "flag to require MFT hit map"}; Configurable> requiredMFTDisks{"requiredMFTDisks", std::vector{0}, "hit map on MFT disks [0,1,2,3,4]. logical-OR of each double-sided disk"}; + Configurable matchingZ{"matchingZ", -77.5, "z position where matching is performed"}; + Configurable cfgBestMatchFinder{"cfgBestMatchFinder", 0, "matching chi2:0, dr:1"}; Configurable cfgCentEstimator{"cfgCentEstimator", 2, "FT0M:0, FT0A:1, FT0C:2"}; Configurable cfgCentMin{"cfgCentMin", -1.f, "min. centrality"}; @@ -118,6 +124,7 @@ struct matchingMFT { o2::ccdb::CcdbApi ccdbApi; Service ccdb; int mRunNumber = -1; + float mBz = 0; template void initCCDB(TBC const& bc) @@ -136,6 +143,10 @@ struct matchingMFT { ccdb->get(geoPath); } o2::mch::TrackExtrap::setField(); + const double centerMFT[3] = {0, 0, -61.4}; + o2::field::MagneticField* field = static_cast(TGeoGlobalMagField::Instance()->GetField()); + mBz = field->getBz(centerMFT); // Get field at centre of MFT + LOGF(info, "Bz at center of MFT = %f kZG", mBz); } void addHistograms() @@ -161,12 +172,13 @@ struct matchingMFT { hMuonType->GetXaxis()->SetBinLabel(4, "MCH-MID"); hMuonType->GetXaxis()->SetBinLabel(5, "MCH standalone"); - const AxisSpec axis_pt{{0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 3, 4, 5, 6, 7, 8, 9, 10}, "p_{T}^{gl} (GeV/c)"}; + const AxisSpec axis_pt{{0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2, 2.5, 3, 3.5, 4, 4.5, 5, 6, 7, 8, 9, 10}, "p_{T}^{gl} (GeV/c)"}; fRegistry.add("MFTMCHMID/primary/correct/hPt", "pT;p_{T} (GeV/c)", kTH1F, {{100, 0.0f, 10}}, false); fRegistry.add("MFTMCHMID/primary/correct/hEtaPhi", "#eta vs. #varphi;#varphi (rad.);#eta", kTH2F, {{180, 0, 2 * M_PI}, {80, -5.f, -1.f}}, false); fRegistry.add("MFTMCHMID/primary/correct/hEtaPhi_MatchedMCHMID", "#eta vs. #varphi;#varphi (rad.);#eta", kTH2F, {{180, 0, 2 * M_PI}, {80, -5.f, -1.f}}, false); - fRegistry.add("MFTMCHMID/primary/correct/hsDelta", "diff. between GL and associated SA;p_{T}^{gl} (GeV/c);(p_{T}^{sa} - p_{T}^{gl})/p_{T}^{gl};#Delta#eta;#Delta#varphi (rad.);", kTHnSparseF, {axis_pt, {100, -0.5, +0.5}, {100, -0.5, +0.5}, {90, -M_PI / 4, M_PI / 4}}, false); + fRegistry.add("MFTMCHMID/primary/correct/hsDelta", "diff. between GL and associated SA;p_{T}^{gl} (GeV/c);(p_{T}^{sa} - p_{T}^{gl})/p_{T}^{gl};#eta^{sa} - #eta^{gl};#varphi^{sa} - #varphi^{gl} (rad.);", kTHnSparseF, {axis_pt, {200, -0.5, +0.5}, {200, -1, +1}, {90, -M_PI / 4, M_PI / 4}}, false); + fRegistry.add("MFTMCHMID/primary/correct/hsDeltaAtMP", "diff. between MFT and MCH-MID;p_{T}^{gl} (GeV/c);#varphi^{MCH-MID} - #varphi^{MFT} (rad.);#eta^{MCH-MID} - #eta^{MFT};", kTHnSparseF, {axis_pt, {90, -M_PI / 4, M_PI / 4}, {200, -1, +1}}, false); fRegistry.add("MFTMCHMID/primary/correct/hDiffCollId", "difference in collision index;collisionId_{TTCA} - collisionId_{MP}", kTH1F, {{41, -20.5, +20.5}}, false); fRegistry.add("MFTMCHMID/primary/correct/hSign", "sign;sign", kTH1F, {{3, -1.5, +1.5}}, false); fRegistry.add("MFTMCHMID/primary/correct/hNclusters", "Nclusters;Nclusters", kTH1F, {{21, -0.5f, 20.5}}, false); @@ -179,7 +191,6 @@ struct matchingMFT { fRegistry.add("MFTMCHMID/primary/correct/hChi2MatchMCHMID", "chi2 match MCH-MID;chi2", kTH1F, {{100, 0.0f, 100}}, false); fRegistry.add("MFTMCHMID/primary/correct/hChi2MatchMCHMFT", "chi2 match MCH-MFT;chi2", kTH1F, {{100, 0.0f, 100}}, false); fRegistry.add("MFTMCHMID/primary/correct/hDCAxy2D", "DCA x vs. y;DCA_{x} (cm);DCA_{y} (cm)", kTH2F, {{200, -0.5, 0.5}, {200, -0.5, +0.5}}, false); - fRegistry.add("MFTMCHMID/primary/correct/hDCAz", "DCA z;DCA_{z} (cm);", kTH1F, {{1000, 0, 10}}, false); fRegistry.add("MFTMCHMID/primary/correct/hDCAxy2DinSigma", "DCA x vs. y in sigma;DCA_{x} (#sigma);DCA_{y} (#sigma)", kTH2F, {{200, -10, 10}, {200, -10, +10}}, false); fRegistry.add("MFTMCHMID/primary/correct/hDCAxy", "DCAxy;DCA_{xy} (cm);", kTH1F, {{100, 0, 1}}, false); fRegistry.add("MFTMCHMID/primary/correct/hDCAxyinSigma", "DCAxy in sigma;DCA_{xy} (#sigma);", kTH1F, {{100, 0, 10}}, false); @@ -188,7 +199,7 @@ struct matchingMFT { fRegistry.add("MFTMCHMID/primary/correct/hDCAxyResolutionvsPt", "DCA_{xy} resolution vs. p_{T};p_{T} (GeV/c);DCA_{xy} resolution (#mum);", kTH2F, {{100, 0, 10.f}, {500, 0, 500}}, false); fRegistry.add("MFTMCHMID/primary/correct/hMCHBitMap", "MCH bit map;MCH bit map", kTH1F, {{1024, -0.5, 1023.5}}, false); fRegistry.add("MFTMCHMID/primary/correct/hMIDBitMap", "MID bit map;MID bit map", kTH1F, {{256, -0.5, 255.5}}, false); - fRegistry.add("MFTMCHMID/primary/correct/hMeanMFTClusterSize", "mean MFT cluster size vs. p;p (GeV/c); #times cos(#lambda)", kTH2F, {{100, 0, 100}, {100, 0, 1}}, false); + fRegistry.add("MFTMCHMID/primary/correct/hMeanMFTClusterSize", "mean MFT cluster size vs. p;p (GeV/c); #times sin(#lambda)", kTH2F, {{100, 0, 100}, {100, 0, 10}}, false); fRegistry.add("MFTMCHMID/primary/correct/hProdVtxZ", "prod. vtx Z of muon;V_{z} (cm)", kTH1F, {{200, -100, 100}}, false); fRegistry.add("MFTMCHMID/primary/correct/hRelDeltaPt", "pT resolution;p_{T}^{gen} (GeV/c);(p_{T}^{rec} - p_{T}^{gen})/p_{T}^{gen}", kTH2F, {{100, 0, 10}, {400, -1, +1}}, false); @@ -198,6 +209,7 @@ struct matchingMFT { fRegistry.add("MFTMCHMID/primary/correct/hDeltaPhi_Neg", "#varphi resolution;p_{T}^{gen} (GeV/c);#varphi^{rec} - #varphi^{gen} (rad.)", kTH2F, {{100, 0, 10}, {400, -0.2, +0.2}}, false); fRegistry.addClone("MFTMCHMID/primary/correct/", "MFTMCHMID/primary/wrong/"); fRegistry.addClone("MFTMCHMID/primary/", "MFTMCHMID/secondary/"); + fRegistry.addClone("MFTMCHMID/", "MCHMID/"); } bool isSelected(const float pt, const float eta, const float rAtAbsorberEnd, const float pDCA, const float chi2_per_ndf, const uint8_t trackType, const float dcaXY) @@ -279,25 +291,67 @@ struct matchingMFT { } if (n > 0) { - return static_cast(clsSize) / static_cast(n) * std::cos(std::atan(track.tgl())); + return static_cast(clsSize) / static_cast(n) * std::fabs(std::sin(std::atan(track.tgl()))); } else { return 0.f; } } - template - void fillHistograms(TCollision const& collision, TFwdTrack fwdtrack, TFwdTracks const&, TMFTTracks const&) + template + void getDeltaEtaDeltaPhiAtMatchingPlane(TCollision const& collision, TFwdTrack const& fwdtrack, TMFTrackCov const& mftCovs, float& deta, float& dphi) { - const auto& mchtrack = fwdtrack.template matchMCHTrack_as(); // MCH-MID - const auto& mfttrack = fwdtrack.template matchMFTTrack_as(); + if (fwdtrack.trackType() != o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { + deta = 999.f; + dphi = 999.f; + return; // do nothing + } + + auto mchtrack = fwdtrack.template matchMCHTrack_as(); // MCH-MID + auto mfttrack = fwdtrack.template matchMFTTrack_as(); + + if (mchtrack.trackType() != o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { + deta = 999.f; + dphi = 999.f; + return; // do nothing + } + + o2::dataformats::GlobalFwdTrack propmuonAtPV = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToMatchingPlane, matchingZ, mBz); + auto mfttrackcov = mftCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]); + + auto muonAtMP = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToMatchingPlane, matchingZ, mBz); // propagated to matching plane + o2::track::TrackParCovFwd mftsaAtMP = getTrackParCovFwd(mfttrack, mfttrackcov); // values at innermost update + mftsaAtMP.propagateToZhelix(matchingZ, mBz); // propagated to matching plane + deta = muonAtMP.getEta() - mftsaAtMP.getEta(); + dphi = muonAtMP.getPhi() - mftsaAtMP.getPhi(); + o2::math_utils::bringToPMPi(dphi); + // reldpt = (muonAtMP.getPt() - mftsaAtMP.getPt()) / muonAtMP.getPt(); + + // o2::track::TrackParCovFwd mftsa = getTrackParCovFwd(mfttrack, mfttrackcov); // values at innermost update + // o2::dataformats::GlobalFwdTrack globalMuonRefit = o2::aod::fwdtrackutils::refitGlobalMuonCov(propmuonAtPV, mftsa); // this is track at IU. + // auto globalMuonAtDCA = o2::aod::fwdtrackutils::propagateTrackParCovFwd(globalMuonRefit, fwdtrack.trackType(), collision, propagationPoint::kToMatchingPlane, matchingZ); + // deta = propmuonAtPV.getEta() - globalMuonAtDCA.getEta(); + // dphi = propmuonAtPV.getPhi() - globalMuonAtDCA.getPhi(); + // o2::math_utils::bringToPMPi(dphi); + // reldpt = (globalMuonAtDCA.getPt() - propmuonAtPV.getPt()) / propmuonAtPV.getPt(); + } + + template + void fillHistograms(TCollision const& collision, TFwdTrack fwdtrack, TFwdTracks const&, TMFTTracks const&, TMFTTracksCov const& mftCovs) + { + if (fwdtrack.trackType() != o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { // only for protection + return; + } + + auto mchtrack = fwdtrack.template matchMCHTrack_as(); // MCH-MID + auto mfttrack = fwdtrack.template matchMFTTrack_as(); if (!fwdtrack.has_mcParticle() || !mchtrack.has_mcParticle() || !mfttrack.has_mcParticle()) { return; } - const auto& mcParticle_MFTMCHMID = fwdtrack.template mcParticle_as(); // this is identical to mcParticle_MCHMID - const auto& mcParticle_MCHMID = mchtrack.template mcParticle_as(); // this is identical to mcParticle_MFTMCHMID - const auto& mcParticle_MFT = mfttrack.template mcParticle_as(); + auto mcParticle_MFTMCHMID = fwdtrack.template mcParticle_as(); // this is identical to mcParticle_MCHMID + auto mcParticle_MCHMID = mchtrack.template mcParticle_as(); // this is identical to mcParticle_MFTMCHMID + auto mcParticle_MFT = mfttrack.template mcParticle_as(); // LOGF(info, "mcParticle_MFTMCHMID.pdgCode() = %d, mcParticle_MCHMID.pdgCode() = %d, mcParticle_MFT.pdgCode() = %d", mcParticle_MFTMCHMID.pdgCode(), mcParticle_MCHMID.pdgCode(), mcParticle_MFT.pdgCode()); // LOGF(info, "mcParticle_MFTMCHMID.globalIndex() = %d, mcParticle_MCHMID.globalIndex() = %d, mcParticle_MFT.globalIndex() = %d", mcParticle_MFTMCHMID.globalIndex(), mcParticle_MCHMID.globalIndex(), mcParticle_MFT.globalIndex()); @@ -334,8 +388,10 @@ struct matchingMFT { bool isPrimary = mcParticle_MCHMID.isPhysicalPrimary() || mcParticle_MCHMID.producedByGenerator(); bool isMatched = (mcParticle_MFT.globalIndex() == mcParticle_MCHMID.globalIndex()) && (mcParticle_MFT.mcCollisionId() == mcParticle_MCHMID.mcCollisionId()); - o2::dataformats::GlobalFwdTrack propmuonAtPV = propagateMuon(fwdtrack, collision, propagationPoint::kToVertex); - o2::dataformats::GlobalFwdTrack propmuonAtDCA = propagateMuon(fwdtrack, collision, propagationPoint::kToDCA); + o2::dataformats::GlobalFwdTrack propmuonAtPV = propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToVertex, matchingZ, mBz); + o2::dataformats::GlobalFwdTrack propmuonAtDCA = propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToDCA, matchingZ, mBz); + o2::dataformats::GlobalFwdTrack propmuonAtDCA_Matched = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToDCA, matchingZ, mBz); + o2::dataformats::GlobalFwdTrack propmuonAtPV_Matched = propagateMuon(mchtrack, mchtrack, collision, propagationPoint::kToVertex, matchingZ, mBz); float p = propmuonAtPV.getP(); float pt = propmuonAtPV.getPt(); @@ -351,7 +407,6 @@ struct matchingMFT { float dcaY = propmuonAtDCA.getY() - collision.posY(); float rAtAbsorberEnd = fwdtrack.rAtAbsorberEnd(); // this works only for GlobalMuonTrack float dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY); - float dcaZ = -dcaXY * std::sinh(eta); float det = cXXatDCA * cYYatDCA - cXYatDCA * cXYatDCA; // determinanat float dcaXYinSigma = 999.f; @@ -362,19 +417,72 @@ struct matchingMFT { } float sigma_dcaXY = dcaXY / dcaXYinSigma; - o2::dataformats::GlobalFwdTrack propmuonAtDCA_Matched = propagateMuon(mchtrack, collision, propagationPoint::kToDCA); float dcaX_Matched = propmuonAtDCA_Matched.getX() - collision.posX(); float dcaY_Matched = propmuonAtDCA_Matched.getY() - collision.posY(); float dcaXY_Matched = std::sqrt(dcaX_Matched * dcaX_Matched + dcaY_Matched * dcaY_Matched); float pDCA = mchtrack.p() * dcaXY_Matched; + float dphiMP = 999.f, detaMP = 999.f; - o2::dataformats::GlobalFwdTrack propmuonAtPV_Matched = propagateMuon(mchtrack, collision, propagationPoint::kToVertex); if (refitGlobalMuon) { eta = mfttrack.eta(); - phi = mfttrack.phi(); - o2::math_utils::bringTo02Pi(phi); + // phi = mfttrack.phi(); + // o2::math_utils::bringTo02Pi(phi); pt = propmuonAtPV_Matched.getP() * std::sin(2.f * std::atan(std::exp(-eta))); p = propmuonAtPV_Matched.getP(); + phi = propmuonAtDCA.getPhi(); + o2::math_utils::bringTo02Pi(phi); + + if constexpr (withMFTCov) { + // auto mfttrackcov = mftCovs.rawIteratorAt(map_mfttrackcovs[mfttrack.globalIndex()]); + // o2::track::TrackParCovFwd mftsa = getTrackParCovFwd(mfttrack, mfttrackcov); // values at innermost update + // o2::dataformats::GlobalFwdTrack globalMuonRefit = o2::aod::fwdtrackutils::refitGlobalMuonCov(propmuonAtPV_Matched, mftsa); // this is track at IU. + // auto globalMuonAtDCA = o2::aod::fwdtrackutils::propagateTrackParCovFwd(globalMuonRefit, fwdtrack.trackType(), collision, propagationPoint::kToDCA, matchingZ, mBz); + + // eta = globalMuonAtDCA.getEta(); + // phi = globalMuonAtDCA.getPhi(); + // o2::math_utils::bringTo02Pi(phi); + // pt = globalMuonAtDCA.getPt(); + // p = globalMuonAtDCA.getP(); + + // eta = globalMuonRefit.getEta(); + // phi = globalMuonRefit.getPhi(); + // o2::math_utils::bringTo02Pi(phi); + // pt = globalMuonRefit.getPt(); + // p = globalMuonRefit.getP(); + + // cXXatDCA = globalMuonAtDCA.getSigma2X(); + // cYYatDCA = globalMuonAtDCA.getSigma2Y(); + // cXYatDCA = globalMuonAtDCA.getSigmaXY(); + // dcaX = globalMuonAtDCA.getX() - collision.posX(); + // dcaY = globalMuonAtDCA.getY() - collision.posY(); + // dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY); + // det = cXXatDCA * cYYatDCA - cXYatDCA * cXYatDCA; // determinanat + // if (det < 0) { + // dcaXYinSigma = 999.f; + // } else { + // dcaXYinSigma = std::sqrt(std::fabs((dcaX * dcaX * cYYatDCA + dcaY * dcaY * cXXatDCA - 2. * dcaX * dcaY * cXYatDCA) / det / 2.)); // dca xy in sigma + // } + // sigma_dcaXY = dcaXY / dcaXYinSigma; + + // o2::track::TrackParCovFwd mftsaAtDCA = getTrackParCovFwd(mfttrack, mfttrackcov); // values at innermost update + // mftsaAtDCA.propagateToZhelix(collision.posZ(), mBz); // propagated to matching plane + // cXXatDCA = mftsaAtDCA.getSigma2X(); + // cYYatDCA = mftsaAtDCA.getSigma2Y(); + // cXYatDCA = mftsaAtDCA.getSigmaXY(); + // dcaX = mftsaAtDCA.getX() - collision.posX(); + // dcaY = mftsaAtDCA.getY() - collision.posY(); + // dcaXY = std::sqrt(dcaX * dcaX + dcaY * dcaY); + // det = cXXatDCA * cYYatDCA - cXYatDCA * cXYatDCA; // determinanat + // if (det < 0) { + // dcaXYinSigma = 999.f; + // } else { + // dcaXYinSigma = std::sqrt(std::fabs((dcaX * dcaX * cYYatDCA + dcaY * dcaY * cXXatDCA - 2. * dcaX * dcaY * cXYatDCA) / det / 2.)); // dca xy in sigma + // } + // sigma_dcaXY = dcaXY / dcaXYinSigma; + + getDeltaEtaDeltaPhiAtMatchingPlane(collision, fwdtrack, mftCovs, detaMP, dphiMP); + o2::math_utils::bringToPMPi(dphiMP); + } } float ptMatchedMCHMID = propmuonAtPV_Matched.getPt(); @@ -385,7 +493,13 @@ struct matchingMFT { float deta = etaMatchedMCHMID - eta; float dphi = phiMatchedMCHMID - phi; o2::math_utils::bringToPMPi(dphi); - if (std::sqrt(std::pow(deta / maxDEta, 2) + std::pow(dphi / maxDPhi, 2)) > 1.f || std::fabs(dpt) > maxRelDPt) { + if (std::sqrt(std::pow(deta / maxDEta, 2) + std::pow(dphi / maxDPhi, 2)) > 1.f) { + return; + } + if (std::sqrt(std::pow(detaMP / maxDEtaMP, 2) + std::pow(dphiMP / maxDPhiMP, 2)) > 1.f) { + return; + } + if (std::fabs(dpt) > maxRelDPt) { return; } @@ -409,6 +523,7 @@ struct matchingMFT { fRegistry.fill(HIST("MFTMCHMID/primary/correct/hEtaPhi"), phi, eta); fRegistry.fill(HIST("MFTMCHMID/primary/correct/hEtaPhi_MatchedMCHMID"), phiMatchedMCHMID, etaMatchedMCHMID); fRegistry.fill(HIST("MFTMCHMID/primary/correct/hsDelta"), pt, dpt, deta, dphi); + fRegistry.fill(HIST("MFTMCHMID/primary/correct/hsDeltaAtMP"), pt, dphiMP, detaMP); fRegistry.fill(HIST("MFTMCHMID/primary/correct/hDiffCollId"), collision.globalIndex() - fwdtrack.collisionId()); fRegistry.fill(HIST("MFTMCHMID/primary/correct/hSign"), fwdtrack.sign()); fRegistry.fill(HIST("MFTMCHMID/primary/correct/hNclusters"), fwdtrack.nClusters()); @@ -421,7 +536,6 @@ struct matchingMFT { fRegistry.fill(HIST("MFTMCHMID/primary/correct/hChi2MatchMCHMID"), fwdtrack.chi2MatchMCHMID()); fRegistry.fill(HIST("MFTMCHMID/primary/correct/hChi2MatchMCHMFT"), fwdtrack.chi2MatchMCHMFT()); fRegistry.fill(HIST("MFTMCHMID/primary/correct/hDCAxy2D"), dcaX, dcaY); - fRegistry.fill(HIST("MFTMCHMID/primary/correct/hDCAz"), dcaZ); fRegistry.fill(HIST("MFTMCHMID/primary/correct/hDCAxy2DinSigma"), dcaX / std::sqrt(cXXatDCA), dcaY / std::sqrt(cYYatDCA)); fRegistry.fill(HIST("MFTMCHMID/primary/correct/hDCAxy"), dcaXY); fRegistry.fill(HIST("MFTMCHMID/primary/correct/hDCAxyinSigma"), dcaXYinSigma); @@ -446,6 +560,7 @@ struct matchingMFT { fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hEtaPhi"), phi, eta); fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hEtaPhi_MatchedMCHMID"), phiMatchedMCHMID, etaMatchedMCHMID); fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hsDelta"), pt, dpt, deta, dphi); + fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hsDeltaAtMP"), pt, dphiMP, detaMP); fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hDiffCollId"), collision.globalIndex() - fwdtrack.collisionId()); fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hSign"), fwdtrack.sign()); fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hNclusters"), fwdtrack.nClusters()); @@ -458,7 +573,6 @@ struct matchingMFT { fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hChi2MatchMCHMID"), fwdtrack.chi2MatchMCHMID()); fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hChi2MatchMCHMFT"), fwdtrack.chi2MatchMCHMFT()); fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hDCAxy2D"), dcaX, dcaY); - fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hDCAz"), dcaZ); fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hDCAxy2DinSigma"), dcaX / std::sqrt(cXXatDCA), dcaY / std::sqrt(cYYatDCA)); fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hDCAxy"), dcaXY); fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hDCAxyinSigma"), dcaXYinSigma); @@ -469,6 +583,7 @@ struct matchingMFT { fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hDCAyResolutionvsPt"), pt, std::sqrt(cYYatDCA) * 1e+4); // convert cm to um fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hDCAxyResolutionvsPt"), pt, sigma_dcaXY * 1e+4); // convert cm to um fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hProdVtxZ"), mcParticle_MFTMCHMID.vz()); + fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hRelDeltaPt"), mcParticle_MFTMCHMID.pt(), (pt - mcParticle_MFTMCHMID.pt()) / mcParticle_MFTMCHMID.pt()); if (mcParticle_MFTMCHMID.pdgCode() > 0) { fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hDeltaEta_Neg"), mcParticle_MFTMCHMID.pt(), eta - mcParticle_MFTMCHMID.eta()); @@ -484,6 +599,7 @@ struct matchingMFT { fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hEtaPhi"), phi, eta); fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hEtaPhi_MatchedMCHMID"), phiMatchedMCHMID, etaMatchedMCHMID); fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hsDelta"), pt, dpt, deta, dphi); + fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hsDeltaAtMP"), pt, dphiMP, detaMP); fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hDiffCollId"), collision.globalIndex() - fwdtrack.collisionId()); fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hSign"), fwdtrack.sign()); fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hNclusters"), fwdtrack.nClusters()); @@ -496,7 +612,6 @@ struct matchingMFT { fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hChi2MatchMCHMID"), fwdtrack.chi2MatchMCHMID()); fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hChi2MatchMCHMFT"), fwdtrack.chi2MatchMCHMFT()); fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hDCAxy2D"), dcaX, dcaY); - fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hDCAz"), dcaZ); fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hDCAxy2DinSigma"), dcaX / std::sqrt(cXXatDCA), dcaY / std::sqrt(cYYatDCA)); fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hDCAxy"), dcaXY); fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hDCAxyinSigma"), dcaXYinSigma); @@ -507,6 +622,7 @@ struct matchingMFT { fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hDCAyResolutionvsPt"), pt, std::sqrt(cYYatDCA) * 1e+4); // convert cm to um fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hDCAxyResolutionvsPt"), pt, sigma_dcaXY * 1e+4); // convert cm to um fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hProdVtxZ"), mcParticle_MFTMCHMID.vz()); + fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hRelDeltaPt"), mcParticle_MFTMCHMID.pt(), (pt - mcParticle_MFTMCHMID.pt()) / mcParticle_MFTMCHMID.pt()); if (mcParticle_MFTMCHMID.pdgCode() > 0) { fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hDeltaEta_Neg"), mcParticle_MFTMCHMID.pt(), eta - mcParticle_MFTMCHMID.eta()); @@ -520,6 +636,7 @@ struct matchingMFT { fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hEtaPhi"), phi, eta); fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hEtaPhi_MatchedMCHMID"), phiMatchedMCHMID, etaMatchedMCHMID); fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hsDelta"), pt, dpt, deta, dphi); + fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hsDeltaAtMP"), pt, dphiMP, detaMP); fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hDiffCollId"), collision.globalIndex() - fwdtrack.collisionId()); fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hSign"), fwdtrack.sign()); fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hNclusters"), fwdtrack.nClusters()); @@ -532,7 +649,6 @@ struct matchingMFT { fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hChi2MatchMCHMID"), fwdtrack.chi2MatchMCHMID()); fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hChi2MatchMCHMFT"), fwdtrack.chi2MatchMCHMFT()); fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hDCAxy2D"), dcaX, dcaY); - fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hDCAz"), dcaZ); fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hDCAxy2DinSigma"), dcaX / std::sqrt(cXXatDCA), dcaY / std::sqrt(cYYatDCA)); fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hDCAxy"), dcaXY); fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hDCAxyinSigma"), dcaXYinSigma); @@ -543,6 +659,7 @@ struct matchingMFT { fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hDCAyResolutionvsPt"), pt, std::sqrt(cYYatDCA) * 1e+4); // convert cm to um fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hDCAxyResolutionvsPt"), pt, sigma_dcaXY * 1e+4); // convert cm to um fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hProdVtxZ"), mcParticle_MFTMCHMID.vz()); + fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hRelDeltaPt"), mcParticle_MFTMCHMID.pt(), (pt - mcParticle_MFTMCHMID.pt()) / mcParticle_MFTMCHMID.pt()); if (mcParticle_MFTMCHMID.pdgCode() > 0) { fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hDeltaEta_Neg"), mcParticle_MFTMCHMID.pt(), eta - mcParticle_MFTMCHMID.eta()); @@ -570,32 +687,81 @@ struct matchingMFT { } std::vector> vec_min_chi2MatchMCHMFT; // std::pair -> chi2MatchMCHMFT; - template - void findBestMatchPerMCHMID(TMuons const& muons) + std::vector> vec_min_dr; // std::pair -> deta + dphi; + std::map, bool> mapCorrectMatch; + + template + void findBestMatchPerMCHMID(TCollision const& collision, TFwdTrack const& fwdtrack, TFwdTracks const& fwdtracks, TMFTTracks const&, TMFTTracksCov const& mftCovs) { - vec_min_chi2MatchMCHMFT.reserve(muons.size()); - for (const auto& muon : muons) { - if (muon.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { - const auto& muons_per_MCHMID = muons.sliceBy(fwdtracksPerMCHTrack, muon.globalIndex()); - // LOGF(info, "stanadalone: muon.globalIndex() = %d, muon.chi2MatchMCHMFT() = %f", muon.globalIndex(), muon.chi2MatchMCHMFT()); - // LOGF(info, "muons_per_MCHMID.size() = %d", muons_per_MCHMID.size()); - - float min_chi2MatchMCHMFT = 1e+10; - std::tuple tupleIds_at_min; - for (const auto& muon_tmp : muons_per_MCHMID) { - if (muon_tmp.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { - // LOGF(info, "muon_tmp.globalIndex() = %d, muon_tmp.matchMCHTrackId() = %d, muon_tmp.matchMFTTrackId() = %d, muon_tmp.chi2MatchMCHMFT() = %f", muon_tmp.globalIndex(), muon_tmp.matchMCHTrackId(), muon_tmp.matchMFTTrackId(), muon_tmp.chi2MatchMCHMFT()); - if (0.f < muon_tmp.chi2MatchMCHMFT() && muon_tmp.chi2MatchMCHMFT() < min_chi2MatchMCHMFT) { - min_chi2MatchMCHMFT = muon_tmp.chi2MatchMCHMFT(); - tupleIds_at_min = std::make_tuple(muon_tmp.globalIndex(), muon_tmp.matchMCHTrackId(), muon_tmp.matchMFTTrackId()); - } - } + if (fwdtrack.trackType() != o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { + return; + } + if (!fwdtrack.has_mcParticle()) { + return; + } + + std::tuple tupleIds_at_min_chi2mftmch; + std::tuple tupleIds_at_min_dr; + float min_chi2MatchMCHMFT = 1e+10; + float min_dr = 1e+10; + auto muons_per_MCHMID = fwdtracks.sliceBy(fwdtracksPerMCHTrack, fwdtrack.globalIndex()); + // LOGF(info, "muons_per_MCHMID.size() = %d", muons_per_MCHMID.size()); + + for (const auto& muon_tmp : muons_per_MCHMID) { + if (muon_tmp.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { + + auto tupleId = std::make_tuple(muon_tmp.globalIndex(), muon_tmp.matchMCHTrackId(), muon_tmp.matchMFTTrackId()); + auto mchtrack = muon_tmp.template matchMCHTrack_as(); // MCH-MID + auto mfttrack = muon_tmp.template matchMFTTrack_as(); + + if (!muon_tmp.has_mcParticle() || !mchtrack.has_mcParticle() || !mfttrack.has_mcParticle()) { + continue; } - vec_min_chi2MatchMCHMFT.emplace_back(tupleIds_at_min); - // mCandidates[tupleIds_at_min] = min_chi2MatchMCHMFT; - // LOGF(info, "min: muon_tmp.globalIndex() = %d, muon_tmp.matchMCHTrackId() = %d, muon_tmp.matchMFTTrackId() = %d, muon_tmp.chi2MatchMCHMFT() = %f", std::get<0>(tupleIds_at_min), std::get<1>(tupleIds_at_min), std::get<2>(tupleIds_at_min), min_chi2MatchMCHMFT); - } - } // end of muon loop + + float deta = 999.f, dphi = 999.f; + getDeltaEtaDeltaPhiAtMatchingPlane(collision, muon_tmp, mftCovs, deta, dphi); + float dr = std::sqrt(deta * deta + dphi * dphi); + + // auto mcParticle_MFTMCHMID = muon_tmp.template mcParticle_as(); // this is identical to mcParticle_MCHMID + auto mcParticle_MCHMID = mchtrack.template mcParticle_as(); // this is identical to mcParticle_MFTMCHMID + auto mcParticle_MFT = mfttrack.template mcParticle_as(); + // float chi2ndf = muon_tmp.chi2() / (2.f * (mchtrack.nClusters() + mfttrack.nClusters()) - 5.f); + + if (mcParticle_MFT.globalIndex() == mcParticle_MCHMID.globalIndex()) { + mapCorrectMatch[tupleId] = true; + } else { + mapCorrectMatch[tupleId] = false; + } + + // if (std::abs(mcParticle_MCHMID.pdgCode()) == 13 && mcParticle_MCHMID.isPhysicalPrimary()) { + // if (mcParticle_MFT.globalIndex() == mcParticle_MCHMID.globalIndex()) { + // LOGF(info, "This is correct match between MFT and MCH-MID: muon_tmp.globalIndex() = %d, chi2/ndf = %f, matching chi2/ndf = %f, mcParticle.pt() = %f, mcParticle.eta() = %f, mcParticle.phi() = %f, reldpt = %f, deta = %f, dphi = %f, dr = %f", muon_tmp.globalIndex(), chi2ndf, muon_tmp.chi2MatchMCHMFT(), mcParticle_MCHMID.pt(), mcParticle_MCHMID.eta(), mcParticle_MCHMID.phi(), reldpt, deta, dphi, dr); + // } else { + // LOGF(info, "This is wrong match between MFT and MCH-MID: muon_tmp.globalIndex() = %d, chi2/ndf = %f, matching chi2/ndf = %f , mcParticle.pt() = %f, mcParticle.eta() = %f, mcParticle.phi() = %f, reldpt = %f, deta = %f, dphi = %f, dr = %f" , muon_tmp.globalIndex(), chi2ndf, muon_tmp.chi2MatchMCHMFT(), mcParticle_MCHMID.pt(), mcParticle_MCHMID.eta(), mcParticle_MCHMID.phi(), reldpt, deta, dphi, dr); + // } + // } + + if (0.f < muon_tmp.chi2MatchMCHMFT() && muon_tmp.chi2MatchMCHMFT() < min_chi2MatchMCHMFT) { + min_chi2MatchMCHMFT = muon_tmp.chi2MatchMCHMFT(); + tupleIds_at_min_chi2mftmch = tupleId; + } + + if (dr < min_dr) { + min_dr = dr; + tupleIds_at_min_dr = tupleId; + } + + } // end of if global muon + } // end of candidates loop + + vec_min_chi2MatchMCHMFT.emplace_back(tupleIds_at_min_chi2mftmch); + vec_min_dr.emplace_back(tupleIds_at_min_dr); + + // auto mcParticleTMP = fwdtrack.template mcParticle_as(); // this is identical to mcParticle_MFTMCHMID + // if (std::abs(mcParticleTMP.pdgCode()) == 13 && mcParticleTMP.isPhysicalPrimary()) { + // LOGF(info, "min chi2: muon_tmp.globalIndex() = %d, muon_tmp.matchMCHTrackId() = %d, muon_tmp.matchMFTTrackId() = %d, muon_tmp.chi2MatchMCHMFT() = %f, correct match = %d", std::get<0>(tupleIds_at_min_chi2mftmch), std::get<1>(tupleIds_at_min_chi2mftmch), std::get<2>(tupleIds_at_min_chi2mftmch), min_chi2MatchMCHMFT, mapCorrectMatch[tupleIds_at_min_chi2mftmch]); + // LOGF(info, "min dr: muon_tmp.globalIndex() = %d, muon_tmp.matchMCHTrackId() = %d, muon_tmp.matchMFTTrackId() = %d, dr = %f, correct match = %d", std::get<0>(tupleIds_at_min_dr), std::get<1>(tupleIds_at_min_dr), std::get<2>(tupleIds_at_min_dr), min_dr, mapCorrectMatch[tupleIds_at_min_dr]); + // } } SliceCache cache; @@ -611,10 +777,8 @@ struct matchingMFT { void processWithoutFTTCA(FilteredMyCollisions const& collisions, MyFwdTracks const& fwdtracks, MyMFTTracks const& mfttracks, aod::BCsWithTimestamps const&, aod::McParticles const&) { - findBestMatchPerMCHMID(fwdtracks); - for (const auto& collision : collisions) { - const auto& bc = collision.template bc_as(); + auto bc = collision.template bc_as(); initCCDB(bc); fRegistry.fill(HIST("Event/hCollisionCounter"), 0); if (!collision.has_mcCollision()) { @@ -623,36 +787,34 @@ struct matchingMFT { if (cfgRequireGoodRCT && !rctChecker.checkTable(collision)) { continue; } - const float centralities[3] = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}; + float centralities[3] = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}; if (centralities[cfgCentEstimator] < cfgCentMin || cfgCentMax < centralities[cfgCentEstimator]) { continue; } fRegistry.fill(HIST("Event/hCollisionCounter"), 1); fillEventHistograms(collision); - const auto& fwdtracks_per_coll = fwdtracks.sliceBy(perCollision, collision.globalIndex()); + auto fwdtracks_per_coll = fwdtracks.sliceBy(perCollision, collision.globalIndex()); + for (const auto& fwdtrack : fwdtracks_per_coll) { if (fwdtrack.trackType() != o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { continue; } - if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack && std::find(vec_min_chi2MatchMCHMFT.begin(), vec_min_chi2MatchMCHMFT.end(), std::make_tuple(fwdtrack.globalIndex(), fwdtrack.matchMCHTrackId(), fwdtrack.matchMFTTrackId())) == vec_min_chi2MatchMCHMFT.end()) { - continue; - } - fillHistograms(collision, fwdtrack, fwdtracks, mfttracks); + fillHistograms(collision, fwdtrack, fwdtracks, mfttracks, nullptr); } // end of fwdtrack loop } // end of collision loop vec_min_chi2MatchMCHMFT.clear(); vec_min_chi2MatchMCHMFT.shrink_to_fit(); + vec_min_dr.clear(); + vec_min_dr.shrink_to_fit(); } PROCESS_SWITCH(matchingMFT, processWithoutFTTCA, "process without FTTCA", false); void processWithFTTCA(FilteredMyCollisions const& collisions, MyFwdTracks const& fwdtracks, MyMFTTracks const& mfttracks, aod::BCsWithTimestamps const&, aod::FwdTrackAssoc const& fwdtrackIndices, aod::McParticles const&) { - findBestMatchPerMCHMID(fwdtracks); - for (const auto& collision : collisions) { - const auto& bc = collision.template bc_as(); + auto bc = collision.template bc_as(); initCCDB(bc); fRegistry.fill(HIST("Event/hCollisionCounter"), 0); if (!collision.has_mcCollision()) { @@ -661,30 +823,96 @@ struct matchingMFT { if (cfgRequireGoodRCT && !rctChecker.checkTable(collision)) { continue; } - const float centralities[3] = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}; + float centralities[3] = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}; if (centralities[cfgCentEstimator] < cfgCentMin || cfgCentMax < centralities[cfgCentEstimator]) { continue; } fRegistry.fill(HIST("Event/hCollisionCounter"), 1); fillEventHistograms(collision); - const auto& fwdtrackIdsThisCollision = fwdtrackIndices.sliceBy(fwdtrackIndicesPerCollision, collision.globalIndex()); + auto fwdtrackIdsThisCollision = fwdtrackIndices.sliceBy(fwdtrackIndicesPerCollision, collision.globalIndex()); for (const auto& fwdtrackId : fwdtrackIdsThisCollision) { - const auto& fwdtrack = fwdtrackId.template fwdtrack_as(); + auto fwdtrack = fwdtrackId.template fwdtrack_as(); if (fwdtrack.trackType() != o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { continue; } - if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack && std::find(vec_min_chi2MatchMCHMFT.begin(), vec_min_chi2MatchMCHMFT.end(), std::make_tuple(fwdtrack.globalIndex(), fwdtrack.matchMCHTrackId(), fwdtrack.matchMFTTrackId())) == vec_min_chi2MatchMCHMFT.end()) { - continue; - } - fillHistograms(collision, fwdtrack, fwdtracks, mfttracks); + fillHistograms(collision, fwdtrack, fwdtracks, mfttracks, nullptr); } // end of fwdtrack loop } // end of collision loop vec_min_chi2MatchMCHMFT.clear(); vec_min_chi2MatchMCHMFT.shrink_to_fit(); + vec_min_dr.clear(); + vec_min_dr.shrink_to_fit(); } PROCESS_SWITCH(matchingMFT, processWithFTTCA, "process with FTTCA", true); + + std::unordered_map map_mfttrackcovs; + + void processWithFTTCA_withMFTCov(FilteredMyCollisions const& collisions, MyFwdTracks const& fwdtracks, MyMFTTracks const& mfttracks, aod::MFTTracksCov const& mftCovs, aod::BCsWithTimestamps const&, aod::FwdTrackAssoc const& fwdtrackIndices, aod::McParticles const&) + { + for (const auto& mfttrackConv : mftCovs) { + map_mfttrackcovs[mfttrackConv.matchMFTTrackId()] = mfttrackConv.globalIndex(); + } + + vec_min_chi2MatchMCHMFT.reserve(fwdtracks.size()); + vec_min_dr.reserve(fwdtracks.size()); + + for (const auto& collision : collisions) { + auto bc = collision.template bc_as(); + initCCDB(bc); + auto fwdtrackIdsThisCollision = fwdtrackIndices.sliceBy(fwdtrackIndicesPerCollision, collision.globalIndex()); + for (const auto& fwdtrackId : fwdtrackIdsThisCollision) { + auto fwdtrack = fwdtrackId.template fwdtrack_as(); + findBestMatchPerMCHMID(collision, fwdtrack, fwdtracks, mfttracks, mftCovs); + } // end of fwdtrack loop + } // end of collision loop + + for (const auto& collision : collisions) { + auto bc = collision.template bc_as(); + initCCDB(bc); + fRegistry.fill(HIST("Event/hCollisionCounter"), 0); + if (!collision.has_mcCollision()) { + continue; + } + if (cfgRequireGoodRCT && !rctChecker.checkTable(collision)) { + continue; + } + float centralities[3] = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}; + if (centralities[cfgCentEstimator] < cfgCentMin || cfgCentMax < centralities[cfgCentEstimator]) { + continue; + } + fRegistry.fill(HIST("Event/hCollisionCounter"), 1); + fillEventHistograms(collision); + + auto fwdtrackIdsThisCollision = fwdtrackIndices.sliceBy(fwdtrackIndicesPerCollision, collision.globalIndex()); + for (const auto& fwdtrackId : fwdtrackIdsThisCollision) { + auto fwdtrack = fwdtrackId.template fwdtrack_as(); + + if (cfgBestMatchFinder == 0) { // chi2 + if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack && std::find(vec_min_chi2MatchMCHMFT.begin(), vec_min_chi2MatchMCHMFT.end(), std::make_tuple(fwdtrack.globalIndex(), fwdtrack.matchMCHTrackId(), fwdtrack.matchMFTTrackId())) == vec_min_chi2MatchMCHMFT.end()) { + continue; + } + } else if (cfgBestMatchFinder == 1) { // dr + if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack && std::find(vec_min_dr.begin(), vec_min_dr.end(), std::make_tuple(fwdtrack.globalIndex(), fwdtrack.matchMCHTrackId(), fwdtrack.matchMFTTrackId())) == vec_min_dr.end()) { + continue; + } + } else { // best match is not selected. Histograms are filled with all global muons. + if (fwdtrack.trackType() != o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { + continue; + } + } + fillHistograms(collision, fwdtrack, fwdtracks, mfttracks, mftCovs); + } // end of fwdtrack loop + } // end of collision loop + + map_mfttrackcovs.clear(); + vec_min_chi2MatchMCHMFT.clear(); + vec_min_chi2MatchMCHMFT.shrink_to_fit(); + vec_min_dr.clear(); + vec_min_dr.shrink_to_fit(); + } + PROCESS_SWITCH(matchingMFT, processWithFTTCA_withMFTCov, "process with FTTCA with MFTCov", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { From cae76d2f729bdcab31c6dd543bea4e6e7549c245 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Sun, 12 Oct 2025 10:18:09 +0000 Subject: [PATCH 2/2] Please consider the following formatting changes --- PWGEM/Dilepton/Tasks/matchingMFT.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGEM/Dilepton/Tasks/matchingMFT.cxx b/PWGEM/Dilepton/Tasks/matchingMFT.cxx index 9a8d492996f..9c2307f417c 100644 --- a/PWGEM/Dilepton/Tasks/matchingMFT.cxx +++ b/PWGEM/Dilepton/Tasks/matchingMFT.cxx @@ -723,7 +723,7 @@ struct matchingMFT { float dr = std::sqrt(deta * deta + dphi * dphi); // auto mcParticle_MFTMCHMID = muon_tmp.template mcParticle_as(); // this is identical to mcParticle_MCHMID - auto mcParticle_MCHMID = mchtrack.template mcParticle_as(); // this is identical to mcParticle_MFTMCHMID + auto mcParticle_MCHMID = mchtrack.template mcParticle_as(); // this is identical to mcParticle_MFTMCHMID auto mcParticle_MFT = mfttrack.template mcParticle_as(); // float chi2ndf = muon_tmp.chi2() / (2.f * (mchtrack.nClusters() + mfttrack.nClusters()) - 5.f);