diff --git a/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuon.cxx b/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuon.cxx index b595c111543..df73ca2391c 100644 --- a/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuon.cxx +++ b/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuon.cxx @@ -92,6 +92,8 @@ struct skimmerPrimaryMuon { Configurable minNmuon{"minNmuon", 0, "min number of muon candidates per collision"}; 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 cfgSlope_dr_chi2MatchMFTMCH{"cfgSlope_dr_chi2MatchMFTMCH", -0.15 / 30, "slope of chiMatchMCHMFT vs. dR"}; + Configurable cfgIntercept_dr_chi2MatchMFTMCH{"cfgIntercept_dr_chi2MatchMFTMCH", 1e+10f, "intercept of chiMatchMCHMFT vs. dR"}; // for z shift for propagation Configurable cfgApplyZShiftFromCCDB{"cfgApplyZShiftFromCCDB", false, "flag to apply z shift"}; @@ -194,8 +196,11 @@ struct skimmerPrimaryMuon { fRegistry.add("MFTMCHMID/hDCAy_PosZ", "DCAy vs. posZ;Z_{vtx} (cm);DCA_{y} (cm)", kTH2F, {{200, -10, +10}, {400, -0.2, +0.2}}, false); fRegistry.add("MFTMCHMID/hDCAx_Phi", "DCAx vs. #varphi;#varphi (rad.);DCA_{x} (cm)", kTH2F, {{90, 0, 2 * M_PI}, {400, -0.2, +0.2}}, false); fRegistry.add("MFTMCHMID/hDCAy_Phi", "DCAy vs. #varphi;#varphi (rad.);DCA_{y} (cm)", kTH2F, {{90, 0, 2 * M_PI}, {400, -0.2, +0.2}}, false); - fRegistry.add("MFTMCHMID/hNmu", "#mu multiplicity;N_{#mu} per collision", kTH1F, {{21, -0.5, 20.5}}, false); + fRegistry.add("MFTMCHMID/hdR_Chi2MatchMCHMFT", "dr vs. matching chi2 MCH-MFT;chi2 match MCH-MFT;#DeltaR", kTH2F, {{200, 0, 50}, {200, 0, 0.5}}, false); + fRegistry.add("MFTMCHMID/hdR_Chi2", "dr vs. chi2;global chi2/ndf;#DeltaR", kTH2F, {{100, 0, 10}, {200, 0, 0.5}}, false); + fRegistry.add("MFTMCHMID/hChi2_Chi2MatchMCHMFT", "chi2 vs. matching chi2 MCH-MFT;chi2 match MCH-MFT;global chi2/ndf", kTH2F, {{200, 0, 50}, {100, 0, 10}}, false); + fRegistry.addClone("MFTMCHMID/", "MCHMID/"); fRegistry.add("MFTMCHMID/hDCAxResolutionvsPt", "DCA_{x} vs. p_{T};p_{T} (GeV/c);DCA_{x} resolution (#mum);", kTH2F, {{100, 0, 10.f}, {500, 0, 500}}, false); fRegistry.add("MFTMCHMID/hDCAyResolutionvsPt", "DCA_{y} vs. p_{T};p_{T} (GeV/c);DCA_{y} resolution (#mum);", kTH2F, {{100, 0, 10.f}, {500, 0, 500}}, false); @@ -525,31 +530,72 @@ struct skimmerPrimaryMuon { std::unordered_map map_mfttrackcovs; std::vector> vec_min_chi2MatchMCHMFT; // std::pair -> chi2MatchMCHMFT; - template - void findBestMatchPerMCHMID(TMuons const& muons) + std::vector> vec_min_dr; // std::pair -> dr; + + template + void findBestMatchPerMCHMID(TCollision const& collision, TFwdTrack const& fwdtrack, TFwdTracks const& fwdtracks, TMFTTracks const&) { - 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 constexpr (isMC) { + if (!fwdtrack.has_mcParticle()) { + return; + } + } + + const auto& muons_per_MCHMID = fwdtracks.sliceBy(fwdtracksPerMCHTrack, fwdtrack.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()); + + o2::dataformats::GlobalFwdTrack propmuonAtPV_Matched = propagateMuon(fwdtrack, fwdtrack, collision, propagationPoint::kToVertex, matchingZ, mBz, mZShift); + float etaMatchedMCHMID = propmuonAtPV_Matched.getEta(); + float phiMatchedMCHMID = propmuonAtPV_Matched.getPhi(); + o2::math_utils::bringTo02Pi(phiMatchedMCHMID); + + float min_chi2MatchMCHMFT = 1e+10, min_dr = 1e+10; + std::tuple tupleIds_at_min_chi2mftmch; + std::tuple tupleIds_at_min_dr; + 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(); // MFTsa + + o2::dataformats::GlobalFwdTrack propmuonAtPV = propagateMuon(muon_tmp, muon_tmp, collision, propagationPoint::kToVertex, matchingZ, mBz, mZShift); + float eta = propmuonAtPV.getEta(); + float phi = propmuonAtPV.getPhi(); + o2::math_utils::bringTo02Pi(phi); + + float deta = etaMatchedMCHMID - eta; + float dphi = phiMatchedMCHMID - phi; + o2::math_utils::bringToPMPi(dphi); + float dr = std::sqrt(deta * deta + dphi * dphi); + + int ndf = 2 * (mchtrack.nClusters() + mfttrack.nClusters()) - 5; + + fRegistry.fill(HIST("MFTMCHMID/hdR_Chi2MatchMCHMFT"), muon_tmp.chi2MatchMCHMFT(), dr); + fRegistry.fill(HIST("MFTMCHMID/hdR_Chi2"), muon_tmp.chi2() / ndf, dr); + fRegistry.fill(HIST("MFTMCHMID/hChi2_Chi2MatchMCHMFT"), muon_tmp.chi2MatchMCHMFT(), muon_tmp.chi2() / ndf); + if (cfgSlope_dr_chi2MatchMFTMCH * muon_tmp.chi2MatchMCHMFT() + cfgIntercept_dr_chi2MatchMFTMCH < dr) { + continue; + } + + // 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_chi2mftmch = std::make_tuple(muon_tmp.globalIndex(), muon_tmp.matchMCHTrackId(), muon_tmp.matchMFTTrackId()); + } + + if (dr < min_dr) { + min_dr = dr; + tupleIds_at_min_dr = tupleId; } - 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 + } + vec_min_chi2MatchMCHMFT.emplace_back(tupleIds_at_min_chi2mftmch); + vec_min_dr.emplace_back(tupleIds_at_min_dr); + // 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); } SliceCache cache; @@ -560,9 +606,19 @@ struct skimmerPrimaryMuon { std::unordered_multimap multiMapSAMuonsPerCollision; // collisionId -> trackIds std::unordered_multimap multiMapGLMuonsPerCollision; // collisionId -> trackIds - void processRec_SA(MyCollisions const& collisions, MyFwdTracks const& fwdtracks, aod::MFTTracks const&, aod::BCsWithTimestamps const&) + void processRec_SA(MyCollisions const& collisions, MyFwdTracks const& fwdtracks, aod::MFTTracks const& mfttracks, aod::BCsWithTimestamps const&) { - findBestMatchPerMCHMID(fwdtracks); + 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); + const auto& fwdtracks_per_coll = fwdtracks.sliceBy(perCollision, collision.globalIndex()); + for (const auto& fwdtrack : fwdtracks_per_coll) { + findBestMatchPerMCHMID(collision, fwdtrack, fwdtracks, mfttracks); + } // end of fwdtrack loop + } // end of collision loop for (const auto& collision : collisions) { const auto& bc = collision.template bc_as(); @@ -624,12 +680,24 @@ struct skimmerPrimaryMuon { 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(skimmerPrimaryMuon, processRec_SA, "process reconstructed info", false); - void processRec_TTCA(MyCollisions const& collisions, MyFwdTracks const& fwdtracks, aod::MFTTracks const&, aod::BCsWithTimestamps const&, aod::FwdTrackAssoc const& fwdtrackIndices) + void processRec_TTCA(MyCollisions const& collisions, MyFwdTracks const& fwdtracks, aod::MFTTracks const& mfttracks, aod::BCsWithTimestamps const&, aod::FwdTrackAssoc const& fwdtrackIndices) { - findBestMatchPerMCHMID(fwdtracks); + 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); + } // end of fwdtrack loop + } // end of collision loop std::unordered_map mapAmb; // fwdtrack.globalIndex() -> bool isAmb; for (const auto& fwdtrack : fwdtracks) { @@ -699,15 +767,28 @@ struct skimmerPrimaryMuon { 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(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) + void processRec_TTCA_withMFTCov(MyCollisions const& collisions, MyFwdTracks const& fwdtracks, aod::MFTTracks const& mfttracks, aod::BCsWithTimestamps const&, aod::FwdTrackAssoc const& fwdtrackIndices, aod::MFTTracksCov const& mftCovs) { for (const auto& mfttrackConv : mftCovs) { map_mfttrackcovs[mfttrackConv.matchMFTTrackId()] = mfttrackConv.globalIndex(); } - findBestMatchPerMCHMID(fwdtracks); + + 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); + } // end of fwdtrack loop + } // end of collision loop std::unordered_map mapAmb; // fwdtrack.globalIndex() -> bool isAmb; for (const auto& fwdtrack : fwdtracks) { @@ -777,12 +858,23 @@ struct skimmerPrimaryMuon { 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(skimmerPrimaryMuon, processRec_TTCA_withMFTCov, "process reconstructed info", false); - void processRec_SA_SWT(MyCollisionsWithSWT const& collisions, MyFwdTracks const& fwdtracks, aod::MFTTracks const&, aod::BCsWithTimestamps const&) + void processRec_SA_SWT(MyCollisionsWithSWT const& collisions, MyFwdTracks const& fwdtracks, aod::MFTTracks const& mfttracks, aod::BCsWithTimestamps const&) { - findBestMatchPerMCHMID(fwdtracks); + 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); + const auto& fwdtracks_per_coll = fwdtracks.sliceBy(perCollision, collision.globalIndex()); + for (const auto& fwdtrack : fwdtracks_per_coll) { + findBestMatchPerMCHMID(collision, fwdtrack, fwdtracks, mfttracks); + } // end of fwdtrack loop + } // end of collision loop for (const auto& collision : collisions) { const auto& bc = collision.template bc_as(); @@ -847,12 +939,24 @@ struct skimmerPrimaryMuon { 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(skimmerPrimaryMuon, processRec_SA_SWT, "process reconstructed info only with standalone", false); - void processRec_TTCA_SWT(MyCollisionsWithSWT const& collisions, MyFwdTracks const& fwdtracks, aod::MFTTracks const&, aod::BCsWithTimestamps const&, aod::FwdTrackAssoc const& fwdtrackIndices) + void processRec_TTCA_SWT(MyCollisionsWithSWT const& collisions, MyFwdTracks const& fwdtracks, aod::MFTTracks const& mfttracks, aod::BCsWithTimestamps const&, aod::FwdTrackAssoc const& fwdtrackIndices) { - findBestMatchPerMCHMID(fwdtracks); + 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); + } // end of fwdtrack loop + } // end of collision loop std::unordered_map mapAmb; // fwdtrack.globalIndex() -> bool isAmb; for (const auto& fwdtrack : fwdtracks) { @@ -924,15 +1028,27 @@ struct skimmerPrimaryMuon { 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(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) + void processRec_TTCA_SWT_withMFTCov(MyCollisionsWithSWT const& collisions, MyFwdTracks const& fwdtracks, aod::MFTTracks const& mfttracks, aod::BCsWithTimestamps const&, aod::FwdTrackAssoc const& fwdtrackIndices, aod::MFTTracksCov const& mftCovs) { for (const auto& mfttrackConv : mftCovs) { map_mfttrackcovs[mfttrackConv.matchMFTTrackId()] = mfttrackConv.globalIndex(); } - findBestMatchPerMCHMID(fwdtracks); + 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); + } // end of fwdtrack loop + } // end of collision loop std::unordered_map mapAmb; // fwdtrack.globalIndex() -> bool isAmb; for (const auto& fwdtrack : fwdtracks) { @@ -1004,12 +1120,23 @@ struct skimmerPrimaryMuon { 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(skimmerPrimaryMuon, processRec_TTCA_SWT_withMFTCov, "process reconstructed info", false); - void processMC_SA(soa::Join const& collisions, MyFwdTracksMC const& fwdtracks, MFTTracksMC const&, aod::BCsWithTimestamps const&) + void processMC_SA(soa::Join const& collisions, MyFwdTracksMC const& fwdtracks, MFTTracksMC const& mfttracks, aod::BCsWithTimestamps const&) { - findBestMatchPerMCHMID(fwdtracks); + 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 fwdtracks_per_coll = fwdtracks.sliceBy(perCollision, collision.globalIndex()); + for (const auto& fwdtrack : fwdtracks_per_coll) { + findBestMatchPerMCHMID(collision, fwdtrack, fwdtracks, mfttracks); + } // end of fwdtrack loop + } // end of collision loop for (const auto& collision : collisions) { auto bc = collision.template bc_as(); @@ -1075,12 +1202,24 @@ struct skimmerPrimaryMuon { 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(skimmerPrimaryMuon, processMC_SA, "process reconstructed and MC info", false); - void processMC_TTCA(soa::Join const& collisions, MyFwdTracksMC const& fwdtracks, MFTTracksMC const&, aod::BCsWithTimestamps const&, aod::FwdTrackAssoc const& fwdtrackIndices) + void processMC_TTCA(soa::Join const& collisions, MyFwdTracksMC const& fwdtracks, MFTTracksMC const& mfttracks, aod::BCsWithTimestamps const&, aod::FwdTrackAssoc const& fwdtrackIndices) { - findBestMatchPerMCHMID(fwdtracks); + 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); + } // end of fwdtrack loop + } // end of collision loop std::unordered_map mapAmb; // fwdtrack.globalIndex() -> bool isAmb; for (const auto& fwdtrack : fwdtracks) { @@ -1155,15 +1294,27 @@ struct skimmerPrimaryMuon { 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(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) + void processMC_TTCA_withMFTCov(soa::Join const& collisions, MyFwdTracksMC const& fwdtracks, MFTTracksMC const& mfttracks, aod::BCsWithTimestamps const&, aod::FwdTrackAssoc const& fwdtrackIndices, aod::MFTTracksCov const& mftCovs) { for (const auto& mfttrackConv : mftCovs) { map_mfttrackcovs[mfttrackConv.matchMFTTrackId()] = mfttrackConv.globalIndex(); } - findBestMatchPerMCHMID(fwdtracks); + 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); + } // end of fwdtrack loop + } // end of collision loop std::unordered_map mapAmb; // fwdtrack.globalIndex() -> bool isAmb; for (const auto& fwdtrack : fwdtracks) { @@ -1238,6 +1389,8 @@ struct skimmerPrimaryMuon { 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(skimmerPrimaryMuon, processMC_TTCA_withMFTCov, "process reconstructed and MC with MFTCov info", false); diff --git a/PWGEM/Dilepton/Tasks/matchingMFT.cxx b/PWGEM/Dilepton/Tasks/matchingMFT.cxx index 10ceff463c6..de8dfc1d722 100644 --- a/PWGEM/Dilepton/Tasks/matchingMFT.cxx +++ b/PWGEM/Dilepton/Tasks/matchingMFT.cxx @@ -97,6 +97,8 @@ struct matchingMFT { 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 cfgSlope_dr_chi2MatchMFTMCH{"cfgSlope_dr_chi2MatchMFTMCH", -0.15 / 30, "slope of chiMatchMCHMFT vs. dR"}; + Configurable cfgIntercept_dr_chi2MatchMFTMCH{"cfgIntercept_dr_chi2MatchMFTMCH", 1e+10f, "intercept of chiMatchMCHMFT vs. dR"}; Configurable cfgCentEstimator{"cfgCentEstimator", 2, "FT0M:0, FT0A:1, FT0C:2"}; Configurable cfgCentMin{"cfgCentMin", -1.f, "min. centrality"}; @@ -227,6 +229,8 @@ struct matchingMFT { 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/hdR_Chi2MatchMCHMFT", "dr vs. matching chi2 MCH-MFT;chi2 match MCH-MFT;#DeltaR", kTH2F, {{200, 0, 50}, {200, 0, 0.5}}, false); + fRegistry.add("MFTMCHMID/primary/correct/hdR_Chi2", "dr vs. chi2;global chi2/ndf;#DeltaR", kTH2F, {{100, 0, 10}, {200, 0, 0.5}}, false); + fRegistry.add("MFTMCHMID/primary/correct/hChi2_Chi2MatchMCHMFT", "chi2 vs. matching chi2 MCH-MFT;chi2 match MCH-MFT;global chi2/ndf", kTH2F, {{200, 0, 50}, {100, 0, 10}}, false); fRegistry.add("MFTMCHMID/primary/correct/hCorrectAsocc", "correct fwdtrack-to-collision association", kTH1F, {{2, -0.5, +1.5}}, false); fRegistry.add("MFTMCHMID/primary/correct/hIsCA", "cellular automaton;isCA", kTH1F, {{2, -0.5, 1.5}}, false); fRegistry.add("MFTMCHMID/primary/correct/hProdVtxZ", "prod. vtx Z of muon;V_{z} (cm)", kTH1F, {{200, -100, 100}}, false); @@ -744,8 +748,8 @@ struct matchingMFT { 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&) + template + void findBestMatchPerMCHMID(TCollision const& collision, TFwdTrack const& fwdtrack, TFwdTracks const& fwdtracks, TMFTTracks const&) { if (fwdtrack.trackType() != o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { return; @@ -777,12 +781,6 @@ struct matchingMFT { continue; } - // float dx = 999.f, dy = 999.f; - // if constexpr (withMFTCov) { - // getDxDyAtMatchingPlane(collision, muon_tmp, mftCovs, dx, dy); - // } - // float dr = std::sqrt(dx * dx + dy * dy); - o2::dataformats::GlobalFwdTrack propmuonAtPV = propagateMuon(muon_tmp, muon_tmp, collision, propagationPoint::kToVertex, matchingZ, mBz, mZShift); float eta = propmuonAtPV.getEta(); float phi = propmuonAtPV.getPhi(); @@ -793,10 +791,14 @@ struct matchingMFT { o2::math_utils::bringToPMPi(dphi); float dr = std::sqrt(deta * deta + dphi * dphi); + if (cfgSlope_dr_chi2MatchMFTMCH * muon_tmp.chi2MatchMCHMFT() + cfgIntercept_dr_chi2MatchMFTMCH < dr) { + continue; + } + // 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); + float chi2ndf = muon_tmp.chi2() / (2.f * (mchtrack.nClusters() + mfttrack.nClusters()) - 5.f); bool isPrimary = mcParticle_MCHMID.isPhysicalPrimary() || mcParticle_MCHMID.producedByGenerator(); bool isMatched = (mcParticle_MFT.globalIndex() == mcParticle_MCHMID.globalIndex()) && (mcParticle_MFT.mcCollisionId() == mcParticle_MCHMID.mcCollisionId()); @@ -804,14 +806,22 @@ struct matchingMFT { if (isPrimary) { if (isMatched) { fRegistry.fill(HIST("MFTMCHMID/primary/correct/hdR_Chi2MatchMCHMFT"), muon_tmp.chi2MatchMCHMFT(), dr); + fRegistry.fill(HIST("MFTMCHMID/primary/correct/hdR_Chi2"), chi2ndf, dr); + fRegistry.fill(HIST("MFTMCHMID/primary/correct/hChi2_Chi2MatchMCHMFT"), muon_tmp.chi2MatchMCHMFT(), chi2ndf); } else { fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hdR_Chi2MatchMCHMFT"), muon_tmp.chi2MatchMCHMFT(), dr); + fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hdR_Chi2"), chi2ndf, dr); + fRegistry.fill(HIST("MFTMCHMID/primary/wrong/hChi2_Chi2MatchMCHMFT"), muon_tmp.chi2MatchMCHMFT(), chi2ndf); } } else { if (isMatched) { fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hdR_Chi2MatchMCHMFT"), muon_tmp.chi2MatchMCHMFT(), dr); + fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hdR_Chi2"), chi2ndf, dr); + fRegistry.fill(HIST("MFTMCHMID/secondary/correct/hChi2_Chi2MatchMCHMFT"), muon_tmp.chi2MatchMCHMFT(), chi2ndf); } else { fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hdR_Chi2MatchMCHMFT"), muon_tmp.chi2MatchMCHMFT(), dr); + fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hdR_Chi2"), chi2ndf, dr); + fRegistry.fill(HIST("MFTMCHMID/secondary/wrong/hChi2_Chi2MatchMCHMFT"), muon_tmp.chi2MatchMCHMFT(), chi2ndf); } } @@ -1002,7 +1012,7 @@ struct matchingMFT { initCCDB(bc); auto fwdtracks_per_coll = fwdtracks.sliceBy(perCollision, collision.globalIndex()); for (const auto& fwdtrack : fwdtracks_per_coll) { - findBestMatchPerMCHMID(collision, fwdtrack, fwdtracks, mfttracks, nullptr); + findBestMatchPerMCHMID(collision, fwdtrack, fwdtracks, mfttracks); } // end of fwdtrack loop } // end of collision loop @@ -1064,7 +1074,7 @@ struct matchingMFT { auto fwdtrackIdsThisCollision = fwdtrackIndices.sliceBy(fwdtrackIndicesPerCollision, collision.globalIndex()); for (const auto& fwdtrackId : fwdtrackIdsThisCollision) { auto fwdtrack = fwdtrackId.template fwdtrack_as(); - findBestMatchPerMCHMID(collision, fwdtrack, fwdtracks, mfttracks, nullptr); + findBestMatchPerMCHMID(collision, fwdtrack, fwdtracks, mfttracks); } // end of fwdtrack loop } // end of collision loop @@ -1133,7 +1143,7 @@ struct matchingMFT { auto fwdtrackIdsThisCollision = fwdtrackIndices.sliceBy(fwdtrackIndicesPerCollision, collision.globalIndex()); for (const auto& fwdtrackId : fwdtrackIdsThisCollision) { auto fwdtrack = fwdtrackId.template fwdtrack_as(); - findBestMatchPerMCHMID(collision, fwdtrack, fwdtracks, mfttracks, mftCovs); + findBestMatchPerMCHMID(collision, fwdtrack, fwdtracks, mfttracks); } // end of fwdtrack loop } // end of collision loop