From fb866e0e736f6bf31dc9d58d8a803e8cc71c7675 Mon Sep 17 00:00:00 2001 From: Daiki Sekihata Date: Sat, 5 Apr 2025 13:39:44 +0200 Subject: [PATCH 1/2] PWGEM/Dilepton: update table producer of muon --- PWGEM/Dilepton/Core/Dilepton.h | 56 +++++++++------ PWGEM/Dilepton/Core/DileptonMC.h | 68 +++++++++++-------- PWGEM/Dilepton/Core/SingleTrackQC.h | 13 +++- PWGEM/Dilepton/Core/SingleTrackQCMC.h | 20 ++++-- PWGEM/Dilepton/DataModel/dileptonTables.h | 15 ++-- .../TableProducer/skimmerPrimaryMuon.cxx | 39 +++++++++-- PWGEM/Dilepton/Tasks/matchingMFT.cxx | 13 ++-- PWGEM/Dilepton/Utils/EMTrackUtilities.h | 27 ++++++++ 8 files changed, 176 insertions(+), 75 deletions(-) diff --git a/PWGEM/Dilepton/Core/Dilepton.h b/PWGEM/Dilepton/Core/Dilepton.h index fec5a3fcbad..19180bc2292 100644 --- a/PWGEM/Dilepton/Core/Dilepton.h +++ b/PWGEM/Dilepton/Core/Dilepton.h @@ -78,7 +78,7 @@ using MyElectron = MyElectrons::iterator; using FilteredMyElectrons = soa::Filtered; using FilteredMyElectron = FilteredMyElectrons::iterator; -using MyMuons = soa::Join; +using MyMuons = soa::Join; using MyMuon = MyMuons::iterator; using FilteredMyMuons = soa::Filtered; using FilteredMyMuon = FilteredMyMuons::iterator; @@ -771,8 +771,8 @@ struct Dilepton { } } - template - bool fillPairInfo(TCollision const& collision, TTrack1 const& t1, TTrack2 const& t2, TCut const& cut) + template + bool fillPairInfo(TCollision const& collision, TTrack1 const& t1, TTrack2 const& t2, TCut const& cut, TAllTracks const& tracks) { if constexpr (ev_id == 1) { if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { @@ -813,6 +813,13 @@ struct Dilepton { if (!cut.template IsSelectedTrack(t1) || !cut.template IsSelectedTrack(t2)) { return false; } + + if (!o2::aod::pwgem::dilepton::utils::emtrackutil::isBestMatch(t1, cut, tracks)) { + return false; + } + if (!o2::aod::pwgem::dilepton::utils::emtrackutil::isBestMatch(t2, cut, tracks)) { + return false; + } } } @@ -1112,8 +1119,8 @@ struct Dilepton { std::vector> used_trackIds; int ndf = 0; - template - void runPairing(TCollisions const& collisions, TLeptons const& posTracks, TLeptons const& negTracks, TPresilce const& perCollision, TCut const& cut) + template + void runPairing(TCollisions const& collisions, TLeptons const& posTracks, TLeptons const& negTracks, TPresilce const& perCollision, TCut const& cut, TAllTracks const& tracks) { for (auto& collision : collisions) { initCCDB(collision); @@ -1190,19 +1197,19 @@ struct Dilepton { int nuls = 0, nlspp = 0, nlsmm = 0; for (auto& [pos, neg] : combinations(CombinationsFullIndexPolicy(posTracks_per_coll, negTracks_per_coll))) { // ULS - bool is_pair_ok = fillPairInfo<0>(collision, pos, neg, cut); + bool is_pair_ok = fillPairInfo<0>(collision, pos, neg, cut, tracks); if (is_pair_ok) { nuls++; } } for (auto& [pos1, pos2] : combinations(CombinationsStrictlyUpperIndexPolicy(posTracks_per_coll, posTracks_per_coll))) { // LS++ - bool is_pair_ok = fillPairInfo<0>(collision, pos1, pos2, cut); + bool is_pair_ok = fillPairInfo<0>(collision, pos1, pos2, cut, tracks); if (is_pair_ok) { nlspp++; } } for (auto& [neg1, neg2] : combinations(CombinationsStrictlyUpperIndexPolicy(negTracks_per_coll, negTracks_per_coll))) { // LS-- - bool is_pair_ok = fillPairInfo<0>(collision, neg1, neg2, cut); + bool is_pair_ok = fillPairInfo<0>(collision, neg1, neg2, cut, tracks); if (is_pair_ok) { nlsmm++; } @@ -1282,25 +1289,25 @@ struct Dilepton { for (auto& pos : selected_posTracks_in_this_event) { // ULS mix for (auto& neg : negTracks_from_event_pool) { - fillPairInfo<1>(collision, pos, neg, cut); + fillPairInfo<1>(collision, pos, neg, cut, tracks); } } for (auto& neg : selected_negTracks_in_this_event) { // ULS mix for (auto& pos : posTracks_from_event_pool) { - fillPairInfo<1>(collision, neg, pos, cut); + fillPairInfo<1>(collision, neg, pos, cut, tracks); } } for (auto& pos1 : selected_posTracks_in_this_event) { // LS++ mix for (auto& pos2 : posTracks_from_event_pool) { - fillPairInfo<1>(collision, pos1, pos2, cut); + fillPairInfo<1>(collision, pos1, pos2, cut, tracks); } } for (auto& neg1 : selected_negTracks_in_this_event) { // LS-- mix for (auto& neg2 : negTracks_from_event_pool) { - fillPairInfo<1>(collision, neg1, neg2, cut); + fillPairInfo<1>(collision, neg1, neg2, cut, tracks); } } } // end of loop over mixed event pool @@ -1315,8 +1322,8 @@ struct Dilepton { } // end of DF - template - bool isPairOK(TCollision const& collision, TTrack1 const& t1, TTrack2 const& t2, TCut const& cut) + template + bool isPairOK(TCollision const& collision, TTrack1 const& t1, TTrack2 const& t2, TCut const& cut, TAllTracks const& tracks) { if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { if (dielectroncuts.cfg_pid_scheme == static_cast(DielectronCut::PIDSchemes::kPIDML)) { @@ -1332,6 +1339,13 @@ struct Dilepton { if (!cut.IsSelectedTrack(t1) || !cut.IsSelectedTrack(t2)) { return false; } + + if (!o2::aod::pwgem::dilepton::utils::emtrackutil::isBestMatch(t1, cut, tracks)) { + return false; + } + if (!o2::aod::pwgem::dilepton::utils::emtrackutil::isBestMatch(t2, cut, tracks)) { + return false; + } } if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { @@ -1398,17 +1412,17 @@ struct Dilepton { auto negTracks_per_coll = negTracks.sliceByCached(perCollision, collision.globalIndex(), cache); for (auto& [pos, neg] : combinations(CombinationsFullIndexPolicy(posTracks_per_coll, negTracks_per_coll))) { // ULS - if (isPairOK(collision, pos, neg, cut)) { + if (isPairOK(collision, pos, neg, cut, tracks)) { passed_pairIds.emplace_back(std::make_pair(pos.globalIndex(), neg.globalIndex())); } } for (auto& [pos1, pos2] : combinations(CombinationsStrictlyUpperIndexPolicy(posTracks_per_coll, posTracks_per_coll))) { // LS++ - if (isPairOK(collision, pos1, pos2, cut)) { + if (isPairOK(collision, pos1, pos2, cut, tracks)) { passed_pairIds.emplace_back(std::make_pair(pos1.globalIndex(), pos2.globalIndex())); } } for (auto& [neg1, neg2] : combinations(CombinationsStrictlyUpperIndexPolicy(negTracks_per_coll, negTracks_per_coll))) { // LS-- - if (isPairOK(collision, neg1, neg2, cut)) { + if (isPairOK(collision, neg1, neg2, cut, tracks)) { passed_pairIds.emplace_back(std::make_pair(neg1.globalIndex(), neg2.globalIndex())); } } @@ -1457,13 +1471,13 @@ struct Dilepton { if (cfgApplyWeightTTCA) { fillPairWeightMap(collisions, positive_electrons, negative_electrons, o2::aod::emprimaryelectron::emeventId, fDielectronCut, electrons); } - runPairing(collisions, positive_electrons, negative_electrons, o2::aod::emprimaryelectron::emeventId, fDielectronCut); + runPairing(collisions, positive_electrons, negative_electrons, o2::aod::emprimaryelectron::emeventId, fDielectronCut, electrons); } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { auto muons = std::get<0>(std::tie(args...)); if (cfgApplyWeightTTCA) { fillPairWeightMap(collisions, positive_muons, negative_muons, o2::aod::emprimarymuon::emeventId, fDimuonCut, muons); } - runPairing(collisions, positive_muons, negative_muons, o2::aod::emprimarymuon::emeventId, fDimuonCut); + runPairing(collisions, positive_muons, negative_muons, o2::aod::emprimarymuon::emeventId, fDimuonCut, muons); } map_weight.clear(); ndf++; @@ -1478,13 +1492,13 @@ struct Dilepton { if (cfgApplyWeightTTCA) { fillPairWeightMap(collisions, positive_electrons, negative_electrons, o2::aod::emprimaryelectron::emeventId, fDielectronCut, electrons); } - runPairing(collisions, positive_electrons, negative_electrons, o2::aod::emprimaryelectron::emeventId, fDielectronCut); + runPairing(collisions, positive_electrons, negative_electrons, o2::aod::emprimaryelectron::emeventId, fDielectronCut, electrons); } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { auto muons = std::get<0>(std::tie(args...)); if (cfgApplyWeightTTCA) { fillPairWeightMap(collisions, positive_muons, negative_muons, o2::aod::emprimarymuon::emeventId, fDimuonCut, muons); } - runPairing(collisions, positive_muons, negative_muons, o2::aod::emprimarymuon::emeventId, fDimuonCut); + runPairing(collisions, positive_muons, negative_muons, o2::aod::emprimarymuon::emeventId, fDimuonCut, muons); } map_weight.clear(); ndf++; diff --git a/PWGEM/Dilepton/Core/DileptonMC.h b/PWGEM/Dilepton/Core/DileptonMC.h index 90c88d7ab08..731a2493e7f 100644 --- a/PWGEM/Dilepton/Core/DileptonMC.h +++ b/PWGEM/Dilepton/Core/DileptonMC.h @@ -71,7 +71,7 @@ using MyMCElectron = MyMCElectrons::iterator; using FilteredMyMCElectrons = soa::Filtered; using FilteredMyMCElectron = FilteredMyMCElectrons::iterator; -using MyMCMuons = soa::Join; +using MyMCMuons = soa::Join; using MyMCMuon = MyMCMuons::iterator; using FilteredMyMCMuons = soa::Filtered; using FilteredMyMCMuon = FilteredMyMCMuons::iterator; @@ -753,8 +753,8 @@ struct DileptonMC { } } - template - bool fillTruePairInfo(TCollision const& collision, TMCCollisions const&, TTrack1 const& t1, TTrack2 const& t2, TCut const& cut, TMCParticles const& mcparticles) + template + bool fillTruePairInfo(TCollision const& collision, TMCCollisions const&, TTrack1 const& t1, TTrack2 const& t2, TCut const& cut, TAllTracks const& tracks, TMCParticles const& mcparticles) { auto t1mc = mcparticles.iteratorAt(t1.emmcparticleId()); auto t2mc = mcparticles.iteratorAt(t2.emmcparticleId()); @@ -789,6 +789,14 @@ struct DileptonMC { if (!cut.template IsSelectedTrack(t1) || !cut.template IsSelectedTrack(t2)) { return false; } + + if (!o2::aod::pwgem::dilepton::utils::emtrackutil::isBestMatch(t1, cut, tracks)) { + return false; + } + if (!o2::aod::pwgem::dilepton::utils::emtrackutil::isBestMatch(t2, cut, tracks)) { + return false; + } + if (!cut.IsSelectedPair(t1, t2)) { return false; } @@ -1184,8 +1192,8 @@ struct DileptonMC { Filter collisionFilter_occupancy_ft0c = eventcuts.cfgFT0COccupancyMin <= o2::aod::evsel::ft0cOccupancyInTimeRange && o2::aod::evsel::ft0cOccupancyInTimeRange < eventcuts.cfgFT0COccupancyMax; using FilteredMyCollisions = soa::Filtered; - template - void runTruePairing(TCollisions const& collisions, TMCLeptons const& posTracks, TMCLeptons const& negTracks, TPreslice const& perCollision, TCut const& cut, TMCCollisions const& mccollisions, TMCParticles const& mcparticles) + template + void runTruePairing(TCollisions const& collisions, TMCLeptons const& posTracks, TMCLeptons const& negTracks, TPreslice const& perCollision, TCut const& cut, TAllTracks const& tracks, TMCCollisions const& mccollisions, TMCParticles const& mcparticles) { for (auto& collision : collisions) { initCCDB(collision); @@ -1207,15 +1215,15 @@ struct DileptonMC { // LOGF(info, "centrality = %f , posTracks_per_coll.size() = %d, negTracks_per_coll.size() = %d", centralities[cfgCentEstimator], posTracks_per_coll.size(), negTracks_per_coll.size()); for (auto& [pos, neg] : combinations(CombinationsFullIndexPolicy(posTracks_per_coll, negTracks_per_coll))) { // ULS - fillTruePairInfo(collision, mccollisions, pos, neg, cut, mcparticles); + fillTruePairInfo(collision, mccollisions, pos, neg, cut, tracks, mcparticles); } // end of ULS pair loop for (auto& [pos1, pos2] : combinations(CombinationsStrictlyUpperIndexPolicy(posTracks_per_coll, posTracks_per_coll))) { // LS++ - fillTruePairInfo(collision, mccollisions, pos1, pos2, cut, mcparticles); + fillTruePairInfo(collision, mccollisions, pos1, pos2, cut, tracks, mcparticles); } // end of LS++ pair loop for (auto& [neg1, neg2] : combinations(CombinationsStrictlyUpperIndexPolicy(negTracks_per_coll, negTracks_per_coll))) { // LS-- - fillTruePairInfo(collision, mccollisions, neg1, neg2, cut, mcparticles); + fillTruePairInfo(collision, mccollisions, neg1, neg2, cut, tracks, mcparticles); } // end of LS-- pair loop } // end of collision loop @@ -1790,8 +1798,8 @@ struct DileptonMC { } // end of collision loop } - template - bool isPairOK(TCollision const& collision, TTrack1 const& t1, TTrack2 const& t2, TCut const& cut) + template + bool isPairOK(TCollision const& collision, TTrack1 const& t1, TTrack2 const& t2, TCut const& cut, TAllTracks const& tracks) { if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { if (dielectroncuts.cfg_pid_scheme == static_cast(DielectronCut::PIDSchemes::kPIDML)) { @@ -1807,6 +1815,12 @@ struct DileptonMC { if (!cut.template IsSelectedTrack(t1) || !cut.template IsSelectedTrack(t2)) { return false; } + if (!o2::aod::pwgem::dilepton::utils::emtrackutil::isBestMatch(t1, cut, tracks)) { + return false; + } + if (!o2::aod::pwgem::dilepton::utils::emtrackutil::isBestMatch(t2, cut, tracks)) { + return false; + } } if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { @@ -1859,7 +1873,7 @@ struct DileptonMC { continue; } - if (isPairOK(collision, pos, neg, cut)) { + if (isPairOK(collision, pos, neg, cut, tracks)) { passed_pairIds.emplace_back(std::make_pair(pos.globalIndex(), neg.globalIndex())); } } @@ -1875,7 +1889,7 @@ struct DileptonMC { continue; } - if (isPairOK(collision, pos1, pos2, cut)) { + if (isPairOK(collision, pos1, pos2, cut, tracks)) { passed_pairIds.emplace_back(std::make_pair(pos1.globalIndex(), pos2.globalIndex())); } } @@ -1890,7 +1904,7 @@ struct DileptonMC { if (cfgEventGeneratorType >= 0 && mccollision_from_neg2.getSubGeneratorId() != cfgEventGeneratorType) { continue; } - if (isPairOK(collision, neg1, neg2, cut)) { + if (isPairOK(collision, neg1, neg2, cut, tracks)) { passed_pairIds.emplace_back(std::make_pair(neg1.globalIndex(), neg2.globalIndex())); } } @@ -1966,8 +1980,8 @@ struct DileptonMC { return true; } - template - void fillUnfolding(TCollisions const& collisions, TTracks1 const& posTracks, TTracks2 const& negTracks, TPresilce const& perCollision, TCut const& cut, TMCCollisions const&, TMCParticles const& mcparticles) + template + void fillUnfolding(TCollisions const& collisions, TTracks1 const& posTracks, TTracks2 const& negTracks, TPresilce const& perCollision, TCut const& cut, TAllTracks const& tracks, TMCCollisions const&, TMCParticles const& mcparticles) { for (auto& collision : collisions) { initCCDB(collision); @@ -2011,7 +2025,7 @@ struct DileptonMC { continue; } - if (!isPairOK(collision, pos, neg, cut)) { // without acceptance + if (!isPairOK(collision, pos, neg, cut, tracks)) { // without acceptance continue; } @@ -2149,7 +2163,7 @@ struct DileptonMC { continue; } - if (!isPairOK(collision, pos1, pos2, cut)) { // without acceptance + if (!isPairOK(collision, pos1, pos2, cut, tracks)) { // without acceptance continue; } @@ -2206,7 +2220,7 @@ struct DileptonMC { if (cfgEventGeneratorType >= 0 && mccollision_from_neg2.getSubGeneratorId() != cfgEventGeneratorType) { continue; } - if (!isPairOK(collision, neg1, neg2, cut)) { // without acceptance + if (!isPairOK(collision, neg1, neg2, cut, tracks)) { // without acceptance continue; } if ((std::abs(mcneg1.pdgCode()) != pdg_lepton || std::abs(mcneg2.pdgCode()) != pdg_lepton) || (mcneg1.emmceventId() != mcneg2.emmceventId())) { @@ -2223,7 +2237,7 @@ struct DileptonMC { continue; } - if (!isPairOK(collision, neg1, neg2, cut)) { // without acceptance + if (!isPairOK(collision, neg1, neg2, cut, tracks)) { // without acceptance continue; } @@ -2291,19 +2305,19 @@ struct DileptonMC { if (cfgApplyWeightTTCA) { fillPairWeightMap(collisions, positive_electrons, negative_electrons, o2::aod::emprimaryelectron::emeventId, fDielectronCut, leptons, mccollisions, mcparticles); } - runTruePairing(collisions, positive_electrons, negative_electrons, o2::aod::emprimaryelectron::emeventId, fDielectronCut, mccollisions, mcparticles); + runTruePairing(collisions, positive_electrons, negative_electrons, o2::aod::emprimaryelectron::emeventId, fDielectronCut, leptons, mccollisions, mcparticles); runGenInfo(collisions, mccollisions, positive_electronsMC, negative_electronsMC, mcparticles); if (cfgFillUnfolding) { - fillUnfolding(collisions, positive_electrons, negative_electrons, o2::aod::emprimaryelectron::emeventId, fDielectronCut, mccollisions, mcparticles); + fillUnfolding(collisions, positive_electrons, negative_electrons, o2::aod::emprimaryelectron::emeventId, fDielectronCut, leptons, mccollisions, mcparticles); } } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { if (cfgApplyWeightTTCA) { fillPairWeightMap(collisions, positive_muons, negative_muons, o2::aod::emprimarymuon::emeventId, fDimuonCut, leptons, mccollisions, mcparticles); } - runTruePairing(collisions, positive_muons, negative_muons, o2::aod::emprimarymuon::emeventId, fDimuonCut, mccollisions, mcparticles); + runTruePairing(collisions, positive_muons, negative_muons, o2::aod::emprimarymuon::emeventId, fDimuonCut, leptons, mccollisions, mcparticles); runGenInfo(collisions, mccollisions, positive_muonsMC, negative_muonsMC, mcparticles); if (cfgFillUnfolding) { - fillUnfolding(collisions, positive_muons, negative_muons, o2::aod::emprimarymuon::emeventId, fDimuonCut, mccollisions, mcparticles); + fillUnfolding(collisions, positive_muons, negative_muons, o2::aod::emprimarymuon::emeventId, fDimuonCut, leptons, mccollisions, mcparticles); } } map_weight.clear(); @@ -2322,19 +2336,19 @@ struct DileptonMC { if (cfgApplyWeightTTCA) { fillPairWeightMap(collisions, positive_electrons, negative_electrons, o2::aod::emprimaryelectron::emeventId, fDielectronCut, leptons, mccollisions, mcparticles_smeared); } - runTruePairing(collisions, positive_electrons, negative_electrons, o2::aod::emprimaryelectron::emeventId, fDielectronCut, mccollisions, mcparticles_smeared); + runTruePairing(collisions, positive_electrons, negative_electrons, o2::aod::emprimaryelectron::emeventId, fDielectronCut, leptons, mccollisions, mcparticles_smeared); runGenInfo(collisions, mccollisions, positive_electronsMC_smeared, negative_electronsMC_smeared, mcparticles_smeared); if (cfgFillUnfolding) { - fillUnfolding(collisions, positive_electrons, negative_electrons, o2::aod::emprimaryelectron::emeventId, fDielectronCut, mccollisions, mcparticles_smeared); + fillUnfolding(collisions, positive_electrons, negative_electrons, o2::aod::emprimaryelectron::emeventId, fDielectronCut, leptons, mccollisions, mcparticles_smeared); } } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { if (cfgApplyWeightTTCA) { fillPairWeightMap(collisions, positive_muons, negative_muons, o2::aod::emprimarymuon::emeventId, fDimuonCut, leptons, mccollisions, mcparticles_smeared); } - runTruePairing(collisions, positive_muons, negative_muons, o2::aod::emprimarymuon::emeventId, fDimuonCut, mccollisions, mcparticles_smeared); + runTruePairing(collisions, positive_muons, negative_muons, o2::aod::emprimarymuon::emeventId, fDimuonCut, leptons, mccollisions, mcparticles_smeared); runGenInfo(collisions, mccollisions, positive_muonsMC_smeared, negative_muonsMC_smeared, mcparticles_smeared); if (cfgFillUnfolding) { - fillUnfolding(collisions, positive_muons, negative_muons, o2::aod::emprimarymuon::emeventId, fDimuonCut, mccollisions, mcparticles_smeared); + fillUnfolding(collisions, positive_muons, negative_muons, o2::aod::emprimarymuon::emeventId, fDimuonCut, leptons, mccollisions, mcparticles_smeared); } } map_weight.clear(); diff --git a/PWGEM/Dilepton/Core/SingleTrackQC.h b/PWGEM/Dilepton/Core/SingleTrackQC.h index 2075a1ba0cf..1d7a2c696cb 100644 --- a/PWGEM/Dilepton/Core/SingleTrackQC.h +++ b/PWGEM/Dilepton/Core/SingleTrackQC.h @@ -61,7 +61,7 @@ using MyElectron = MyElectrons::iterator; using FilteredMyElectrons = soa::Filtered; using FilteredMyElectron = FilteredMyElectrons::iterator; -using MyMuons = soa::Join; +using MyMuons = soa::Join; using MyMuon = MyMuons::iterator; using FilteredMyMuons = soa::Filtered; using FilteredMyMuon = FilteredMyMuons::iterator; @@ -267,7 +267,7 @@ struct SingleTrackQC { } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { const AxisSpec axis_pt{ConfPtlBins, "p_{T,#mu} (GeV/c)"}; const AxisSpec axis_eta{50, -6, -1, "#eta_{#mu}"}; - const AxisSpec axis_phi{90, 0, 2 * M_PI, "#varphi_{#mu} (rad.)"}; + const AxisSpec axis_phi{36, 0, 2 * M_PI, "#varphi_{#mu} (rad.)"}; const AxisSpec axis_dca{ConfDCABins, "DCA_{#mu}^{XY} (#sigma)"}; // track info @@ -281,7 +281,7 @@ struct SingleTrackQC { fRegistry.add("Track/positive/hDCAyRes_Pt", "DCA_{y} resolution vs. pT;p_{T} (GeV/c);DCA_{y} resolution (#mum)", kTH2F, {{100, 0, 10}, {500, 0, 500}}, false); fRegistry.add("Track/positive/hNclsMCH", "number of MCH clusters", kTH1F, {{21, -0.5, 20.5}}, false); fRegistry.add("Track/positive/hNclsMFT", "number of MFT clusters", kTH1F, {{11, -0.5, 10.5}}, false); - fRegistry.add("Track/positive/hPDCA", "pDCA;r at absorber end (cm);p #times DCA (GeV/c #upoint cm)", kTH2F, {{100, 0, 100}, {100, 0.0f, 1000}}, false); + fRegistry.add("Track/positive/hPDCA", "pDCA;R at absorber end (cm);p #times DCA (GeV/c #upoint cm)", kTH2F, {{100, 0, 100}, {100, 0.0f, 1000}}, false); fRegistry.add("Track/positive/hChi2", "chi2;chi2", kTH1F, {{100, 0.0f, 100}}, false); fRegistry.add("Track/positive/hChi2MatchMCHMID", "chi2 match MCH-MID;chi2", kTH1F, {{100, 0.0f, 100}}, false); fRegistry.add("Track/positive/hChi2MatchMCHMFT", "chi2 match MCH-MFT;chi2", kTH1F, {{100, 0.0f, 100}}, false); @@ -616,6 +616,10 @@ struct SingleTrackQC { if (!cut.template IsSelectedTrack(track)) { continue; } + if (!o2::aod::pwgem::dilepton::utils::emtrackutil::isBestMatch(track, cut, tracks)) { + continue; + } + fillMuonInfo(track); } // end of track loop } @@ -664,6 +668,9 @@ struct SingleTrackQC { if (!cut.template IsSelectedTrack(track)) { continue; } + if (!o2::aod::pwgem::dilepton::utils::emtrackutil::isBestMatch(track, cut, tracks)) { + continue; + } passed_trackIds.emplace_back(track.globalIndex()); } // end of track loop } diff --git a/PWGEM/Dilepton/Core/SingleTrackQCMC.h b/PWGEM/Dilepton/Core/SingleTrackQCMC.h index eaaa78d2687..c3aeedeb791 100644 --- a/PWGEM/Dilepton/Core/SingleTrackQCMC.h +++ b/PWGEM/Dilepton/Core/SingleTrackQCMC.h @@ -62,7 +62,7 @@ using MyMCElectrons = soa::Join; -using MyMCMuons = soa::Join; +using MyMCMuons = soa::Join; using MyMCMuon = MyMCMuons::iterator; using FilteredMyMCMuons = soa::Filtered; @@ -311,7 +311,7 @@ struct SingleTrackQCMC { } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { const AxisSpec axis_pt{ConfPtlBins, "p_{T,#mu} (GeV/c)"}; const AxisSpec axis_eta{50, -6, -1, "#eta_{#mu}"}; - const AxisSpec axis_phi{90, 0.0, 2 * M_PI, "#varphi_{#mu} (rad.)"}; + const AxisSpec axis_phi{36, 0.0, 2 * M_PI, "#varphi_{#mu} (rad.)"}; const AxisSpec axis_dca{ConfDCABins, "DCA_{#mu}^{XY} (#sigma)"}; const AxisSpec axis_charge_gen{3, -1.5, +1.5, "true charge"}; @@ -338,7 +338,7 @@ struct SingleTrackQCMC { fRegistry.add("Track/lf/positive/hDCAyRes_Pt", "DCA_{y} resolution vs. pT;p_{T} (GeV/c);DCA_{y} resolution (#mum)", kTH2F, {{100, 0, 10}, {500, 0, 500}}, false); fRegistry.add("Track/lf/positive/hNclsMCH", "number of MCH clusters", kTH1F, {{21, -0.5, 20.5}}, false); fRegistry.add("Track/lf/positive/hNclsMFT", "number of MFT clusters", kTH1F, {{11, -0.5, 10.5}}, false); - fRegistry.add("Track/lf/positive/hPDCA", "pDCA;r at absorber (cm);p #times DCA (GeV/c #upoint cm)", kTH2F, {{100, 0, 100}, {100, 0.0f, 1000}}, false); + fRegistry.add("Track/lf/positive/hPDCA", "pDCA;R at absorber (cm);p #times DCA (GeV/c #upoint cm)", kTH2F, {{100, 0, 100}, {100, 0.0f, 1000}}, false); fRegistry.add("Track/lf/positive/hChi2", "chi2;chi2", kTH1F, {{100, 0.0f, 100}}, false); fRegistry.add("Track/lf/positive/hChi2MatchMCHMID", "chi2 match MCH-MID;chi2", kTH1F, {{100, 0.0f, 100}}, false); fRegistry.add("Track/lf/positive/hChi2MatchMCHMFT", "chi2 match MCH-MFT;chi2", kTH1F, {{100, 0.0f, 100}}, false); @@ -728,6 +728,9 @@ struct SingleTrackQCMC { if (!isInAcceptance(mctrack)) { continue; } + if (!mctrack.has_mothers()) { + continue; + } auto mccollision_from_track = mctrack.template emmcevent_as(); if (cfgEventGeneratorType >= 0 && mccollision_from_track.getSubGeneratorId() != cfgEventGeneratorType) { @@ -752,11 +755,11 @@ struct SingleTrackQCMC { if (!cut.template IsSelectedTrack(track)) { continue; } + if (!o2::aod::pwgem::dilepton::utils::emtrackutil::isBestMatch(track, cut, tracks)) { + continue; + } } - if (!mctrack.has_mothers()) { - continue; - } auto mcmother = mcparticles.iteratorAt(mctrack.mothersIds()[0]); int pdg_mother = std::abs(mcmother.pdgCode()); @@ -941,10 +944,13 @@ struct SingleTrackQCMC { if (cfgEventGeneratorType >= 0 && mccollision_from_track.getSubGeneratorId() != cfgEventGeneratorType) { continue; } - if (!cut.template IsSelectedTrack(track)) { continue; } + + if (!o2::aod::pwgem::dilepton::utils::emtrackutil::isBestMatch(track, cut, tracks)) { + continue; + } passed_trackIds.emplace_back(track.globalIndex()); } // end of track loop } diff --git a/PWGEM/Dilepton/DataModel/dileptonTables.h b/PWGEM/Dilepton/DataModel/dileptonTables.h index 2084b0dbf4d..7ae349d91cc 100644 --- a/PWGEM/Dilepton/DataModel/dileptonTables.h +++ b/PWGEM/Dilepton/DataModel/dileptonTables.h @@ -529,11 +529,12 @@ using EMPrimaryElectronPrefilterBitDerived = EMPrimaryElectronsPrefilterBitDeriv namespace emprimarymuon { -DECLARE_SOA_INDEX_COLUMN(EMEvent, emevent); //! -DECLARE_SOA_COLUMN(CollisionId, collisionId, int); //! -DECLARE_SOA_COLUMN(FwdTrackId, fwdtrackId, int); //! -DECLARE_SOA_COLUMN(MFTTrackId, mfttrackId, int); //! -DECLARE_SOA_COLUMN(MCHTrackId, mchtrackId, int); //! +DECLARE_SOA_INDEX_COLUMN(EMEvent, emevent); //! +DECLARE_SOA_COLUMN(CollisionId, collisionId, int); //! +DECLARE_SOA_COLUMN(FwdTrackId, fwdtrackId, int); //! +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 @@ -621,6 +622,10 @@ DECLARE_SOA_TABLE(EMAmbiguousMuonSelfIds, "AOD", "EMAMBMUSELFID", emprimarymuon: // iterators using EMAmbiguousMuonSelfId = EMAmbiguousMuonSelfIds::iterator; +DECLARE_SOA_TABLE(EMGlobalMuonSelfIds, "AOD", "EMGLMUSELFID", emprimarymuon::GlobalMuonsWithSameMFTIds); // To be joined with EMPrimaryMuons table at analysis level. +// iterators +using EMGlobalMuonSelfId = EMGlobalMuonSelfIds::iterator; + } // namespace o2::aod #endif // PWGEM_DILEPTON_DATAMODEL_DILEPTONTABLES_H_ diff --git a/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuon.cxx b/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuon.cxx index 3e17b8845cc..b64682ec120 100644 --- a/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuon.cxx +++ b/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuon.cxx @@ -304,14 +304,14 @@ struct associateAmbiguousMuon { void process(aod::EMPrimaryMuons const& muons) { - for (auto& muon : muons) { + for (const auto& muon : muons) { auto muons_with_same_trackId = muons.sliceBy(perTrack, muon.fwdtrackId()); ambmuon_self_Ids.reserve(muons_with_same_trackId.size()); - for (auto& amp_muon : muons_with_same_trackId) { - if (amp_muon.globalIndex() == muon.globalIndex()) { // don't store myself. + for (const auto& amb_muon : muons_with_same_trackId) { + if (amb_muon.globalIndex() == muon.globalIndex()) { // don't store myself. continue; } - ambmuon_self_Ids.emplace_back(amp_muon.globalIndex()); + ambmuon_self_Ids.emplace_back(amb_muon.globalIndex()); } em_amb_muon_ids(ambmuon_self_Ids); ambmuon_self_Ids.clear(); @@ -319,9 +319,38 @@ struct associateAmbiguousMuon { } } }; +struct associateSameMFT { + Produces em_same_mft_ids; + + SliceCache cache; + PresliceUnsorted perMFTTrack = o2::aod::emprimarymuon::mfttrackId; + std::vector self_Ids; + + void process(aod::EMPrimaryMuons const& muons) + { + for (const auto& muon : muons) { + if (muon.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { + auto muons_with_same_mfttrackId = muons.sliceBy(perMFTTrack, muon.mfttrackId()); + self_Ids.reserve(muons_with_same_mfttrackId.size()); + for (const auto& global_muon : muons_with_same_mfttrackId) { + if (global_muon.globalIndex() == muon.globalIndex()) { // don't store myself. + continue; + } + self_Ids.emplace_back(global_muon.globalIndex()); + } + em_same_mft_ids(self_Ids); + self_Ids.clear(); + self_Ids.shrink_to_fit(); + } else { // for standalone muons + em_same_mft_ids(std::vector{}); // empty + } + } // end of muon loop + } +}; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ adaptAnalysisTask(cfgc, TaskName{"skimmer-primary-muon"}), - adaptAnalysisTask(cfgc, TaskName{"associate-ambiguous-muon"})}; + adaptAnalysisTask(cfgc, TaskName{"associate-ambiguous-muon"}), + adaptAnalysisTask(cfgc, TaskName{"associate-same-mft"})}; } diff --git a/PWGEM/Dilepton/Tasks/matchingMFT.cxx b/PWGEM/Dilepton/Tasks/matchingMFT.cxx index 81f3a91669f..b5e71cf7837 100644 --- a/PWGEM/Dilepton/Tasks/matchingMFT.cxx +++ b/PWGEM/Dilepton/Tasks/matchingMFT.cxx @@ -60,6 +60,7 @@ struct matchingMFT { Configurable maxEtaSA{"maxEtaSA", -2.5, "max. eta acceptance for MCH-MID"}; Configurable minEtaGL{"minEtaGL", -3.6, "min. eta acceptance for MFT-MCH-MID"}; Configurable maxEtaGL{"maxEtaGL", -2.5, "max. eta acceptance for MFT-MCH-MID"}; + Configurable minRabsGL{"minRabsGL", 27.6, "min. R at absorber end for global muons (min. eta = -3.6)"}; // std::tan(2.f * std::atan(std::exp(- -3.6)) ) * -505. 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"}; @@ -70,7 +71,6 @@ struct matchingMFT { Configurable maxChi2SA{"maxChi2SA", 1e+6, "max. chi2 for standalone muon"}; Configurable maxChi2GL{"maxChi2GL", 1e+6, "max. chi2 for global muon"}; Configurable refitGlobalMuon{"refitGlobalMuon", true, "flag to refit global muon"}; - Configurable applyEtaCutToSAinGL{"applyEtaCutToSAinGL", false, "flag to apply eta cut to samuon in global muon"}; Configurable requireTrueAssociation{"requireTrueAssociation", false, "flag to require true mc collision association"}; Configurable cfgCentEstimator{"cfgCentEstimator", 2, "FT0M:0, FT0A:1, FT0C:2"}; @@ -172,7 +172,7 @@ struct matchingMFT { fRegistry.addClone("MFTMCHMID/primary/", "MFTMCHMID/secondary/"); } - bool isSelected(const float pt, const float eta, const float rAtAbsorberEnd, const float pDCA, const float chi2, const uint8_t trackType, const float etaMatchedMCHMID, const float dcaXY) + bool isSelected(const float pt, const float eta, const float rAtAbsorberEnd, const float pDCA, const float chi2, const uint8_t trackType, const float dcaXY) { if (pt < minPt || maxPt < pt) { return false; @@ -188,15 +188,15 @@ struct matchingMFT { if (eta < minEtaGL || maxEtaGL < eta) { return false; } - if (applyEtaCutToSAinGL && (etaMatchedMCHMID < minEtaGL || maxEtaGL < etaMatchedMCHMID)) { - return false; - } if (maxDCAxy < dcaXY) { return false; } if (chi2 < 0.f || maxChi2GL < chi2) { return false; } + if (rAtAbsorberEnd < minRabsGL || maxRabs < rAtAbsorberEnd) { + return false; + } } else if (trackType == static_cast(o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack)) { if (eta < minEtaSA || maxEtaSA < eta) { return false; @@ -269,7 +269,6 @@ struct matchingMFT { 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; - int nClustersMFT = mfttrack.nClusters(); float chi2mft = mfttrack.chi2(); @@ -280,7 +279,7 @@ struct matchingMFT { pt = propmuonAtPV_Matched.getP() * std::sin(2.f * std::atan(std::exp(-eta))); } - if (!isSelected(pt, eta, rAtAbsorberEnd, pDCA, fwdtrack.chi2(), fwdtrack.trackType(), etaMatchedMCHMID, dcaXY)) { + if (!isSelected(pt, eta, rAtAbsorberEnd, pDCA, fwdtrack.chi2(), fwdtrack.trackType(), dcaXY)) { return; } diff --git a/PWGEM/Dilepton/Utils/EMTrackUtilities.h b/PWGEM/Dilepton/Utils/EMTrackUtilities.h index d71211c1321..cddcb08633b 100644 --- a/PWGEM/Dilepton/Utils/EMTrackUtilities.h +++ b/PWGEM/Dilepton/Utils/EMTrackUtilities.h @@ -73,6 +73,33 @@ float fwdDcaXYinSigma(T const& track) // } } //_______________________________________________________________________ +template +bool isBestMatch(TTrack const& track, TCut const& cut, TTracks const& tracks) +{ + // this is only for muon at forward rapidity + if (track.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { + std::map map_chi2MCHMFT; + map_chi2MCHMFT[track.globalIndex()] = track.chi2MatchMCHMFT(); // add myself + for (const auto& glmuonId : track.globalMuonsWithSameMFTIds()) { + const auto& candidate = tracks.rawIteratorAt(glmuonId); + if (candidate.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { + if (cut.template IsSelectedTrack(candidate)) { + map_chi2MCHMFT[candidate.globalIndex()] = candidate.chi2MatchMCHMFT(); + } + } + } + if (map_chi2MCHMFT.begin()->first == track.globalIndex()) { // search for minimum matching-chi2 + map_chi2MCHMFT.clear(); + return true; + } else { + map_chi2MCHMFT.clear(); + return false; + } + } else { + return true; + } +} +//_______________________________________________________________________ template float sigmaPt(T const& track) { From 238348049c6be97aac9a9e70c26b77af9310967f Mon Sep 17 00:00:00 2001 From: Daiki Sekihata Date: Sat, 5 Apr 2025 18:31:33 +0200 Subject: [PATCH 2/2] Update EMTrackUtilities.h --- PWGEM/Dilepton/Utils/EMTrackUtilities.h | 1 + 1 file changed, 1 insertion(+) diff --git a/PWGEM/Dilepton/Utils/EMTrackUtilities.h b/PWGEM/Dilepton/Utils/EMTrackUtilities.h index cddcb08633b..77a1c8ca7ed 100644 --- a/PWGEM/Dilepton/Utils/EMTrackUtilities.h +++ b/PWGEM/Dilepton/Utils/EMTrackUtilities.h @@ -18,6 +18,7 @@ #include #include #include +#include //_______________________________________________________________________ namespace o2::aod::pwgem::dilepton::utils::emtrackutil