diff --git a/PWGEM/Dilepton/Core/DileptonHadronMPC.h b/PWGEM/Dilepton/Core/DileptonHadronMPC.h index 56b2c837cc8..e6d7eb0960e 100644 --- a/PWGEM/Dilepton/Core/DileptonHadronMPC.h +++ b/PWGEM/Dilepton/Core/DileptonHadronMPC.h @@ -35,7 +35,6 @@ #include "CCDB/BasicCCDBManager.h" #include "CommonConstants/LHCConstants.h" -#include "DataFormatsParameters/GRPECSObject.h" #include "DataFormatsParameters/GRPLHCIFData.h" #include "DataFormatsParameters/GRPMagField.h" #include "DataFormatsParameters/GRPObject.h" @@ -56,6 +55,7 @@ #include #include #include +#include #include #include @@ -99,15 +99,14 @@ struct DileptonHadronMPC { Configurable cfgEP2Estimator_for_Mix{"cfgEP2Estimator_for_Mix", 3, "FT0M:0, FT0A:1, FT0C:2, BTot:3, BPos:4, BNeg:5, FV0A:6"}; Configurable cfgCentEstimator{"cfgCentEstimator", 2, "FT0M:0, FT0A:1, FT0C:2"}; Configurable cfgOccupancyEstimator{"cfgOccupancyEstimator", 0, "FT0C:0, Track:1"}; - Configurable cfgDoMix{"cfgDoMix", true, "flag for event mixing"}; + Configurable cfgDoMix{"cfgDoMix", true, "flag for event mixing between 2 leptons (i.e. R factor)"}; Configurable ndepth_lepton{"ndepth_lepton", 100, "depth for event mixing between lepton-lepton"}; - Configurable ndepth_hadron{"ndepth_hadron", 1, "depth for event mixing between hadron-hadron"}; + Configurable ndepth_hadron{"ndepth_hadron", 10, "depth for event mixing between hadron-hadron"}; Configurable ndiff_bc_mix{"ndiff_bc_mix", 594, "difference in global BC required in mixed events"}; ConfigurableAxis ConfVtxBins{"ConfVtxBins", {VARIABLE_WIDTH, -10.0f, -8.f, -6.f, -4.f, -2.f, 0.f, 2.f, 4.f, 6.f, 8.f, 10.f}, "Mixing bins - z-vertex"}; ConfigurableAxis ConfCentBins{"ConfCentBins", {VARIABLE_WIDTH, 0.0f, 5.0f, 10.0f, 20.0f, 30.0f, 40.0f, 50.0f, 60.0f, 70.0f, 80.0f, 90.0f, 100.f, 999.f}, "Mixing bins - centrality"}; ConfigurableAxis ConfEPBins{"ConfEPBins", {16, -M_PI / 2, +M_PI / 2}, "Mixing bins - event plane angle"}; ConfigurableAxis ConfOccupancyBins{"ConfOccupancyBins", {VARIABLE_WIDTH, -1, 1e+10}, "Mixing bins - occupancy"}; - Configurable cfg_swt_name{"cfg_swt_name", "fHighTrackMult", "desired software trigger name"}; // 1 trigger name per 1 task. fHighTrackMult, fHighFt0Mult // Configurable cfgNtracksPV08Min{"cfgNtracksPV08Min", -1, "min. multNTracksPV"}; // Configurable cfgNtracksPV08Max{"cfgNtracksPV08Max", static_cast(1e+9), "max. multNTracksPV"}; Configurable cfgApplyWeightTTCA{"cfgApplyWeightTTCA", false, "flag to apply weighting by 1/N"}; @@ -129,9 +128,9 @@ struct DileptonHadronMPC { std::string prefix = "eventcut_group"; Configurable cfgZvtxMin{"cfgZvtxMin", -10.f, "min. Zvtx"}; Configurable cfgZvtxMax{"cfgZvtxMax", +10.f, "max. Zvtx"}; - Configurable cfgRequireSel8{"cfgRequireSel8", true, "require sel8 in event cut"}; + Configurable cfgRequireSel8{"cfgRequireSel8", false, "require sel8 in event cut"}; Configurable cfgRequireFT0AND{"cfgRequireFT0AND", true, "require FT0AND in event cut"}; - Configurable cfgRequireNoTFB{"cfgRequireNoTFB", false, "require No time frame border in event cut"}; + Configurable cfgRequireNoTFB{"cfgRequireNoTFB", true, "require No time frame border in event cut"}; Configurable cfgRequireNoITSROFB{"cfgRequireNoITSROFB", false, "require no ITS readout frame border in event cut"}; Configurable cfgRequireNoSameBunchPileup{"cfgRequireNoSameBunchPileup", false, "require no same bunch pileup in event cut"}; Configurable cfgRequireVertexITSTPC{"cfgRequireVertexITSTPC", false, "require Vertex ITSTPC in event cut"}; // ITS-TPC matched track contributes PV. @@ -259,8 +258,8 @@ struct DileptonHadronMPC { Configurable cfg_max_eta_track{"cfg_max_eta_track", -2.5, "max eta for single track"}; Configurable cfg_min_phi_track{"cfg_min_phi_track", 0.f, "min phi for single track"}; Configurable cfg_max_phi_track{"cfg_max_phi_track", 6.3, "max phi for single track"}; - Configurable cfg_min_ncluster_mft{"cfg_min_ncluster_mft", 6, "min ncluster MFT"}; - Configurable cfg_min_ncluster_mch{"cfg_min_ncluster_mch", 8, "min ncluster MCH"}; + Configurable cfg_min_ncluster_mft{"cfg_min_ncluster_mft", 5, "min ncluster MFT"}; + Configurable cfg_min_ncluster_mch{"cfg_min_ncluster_mch", 5, "min ncluster MCH"}; Configurable cfg_max_chi2{"cfg_max_chi2", 1e+6, "max chi2/ndf"}; Configurable cfg_max_chi2mft{"cfg_max_chi2mft", 1e+6, "max chi2/ndf"}; // Configurable cfg_max_matching_chi2_mftmch{"cfg_max_matching_chi2_mftmch", 40, "max chi2 for MFT-MCH matching"}; @@ -288,9 +287,8 @@ struct DileptonHadronMPC { Configurable cfg_max_eta_track{"cfg_max_eta_track", +0.8, "max eta for ref. track"}; // Configurable cfg_min_phi_track{"cfg_min_phi_track", 0.0, "min phi for ref. track"}; // Configurable cfg_max_phi_track{"cfg_max_phi_track", 6.3, "max phi for ref. track"}; - // Configurable cfg_max_dcaxy{"cfg_max_dcaxy", 0.5, "max dca XY for single track in cm"}; - // Configurable cfg_max_dcaz{"cfg_max_dcaz", 0.5, "max dca Z for single track in cm"}; - Configurable cfg_track_bits{"cfg_track_bits", 5765, "required track bits"}; // default:645, loose:0, tight:778 + Configurable cfg_track_bits{"cfg_track_bits", 0, "required track bits"}; + Configurable cfg_reftrack_type{"cfg_reftrack_type", 0, "ref. track type for 2PC. 0:kCB, 1:kMFTsa"}; // this is unique in the derived data. Check which track type you store as aod::EMPrimaryTracks. } trackcuts; struct : ConfigurableGroup { @@ -417,11 +415,11 @@ struct DileptonHadronMPC { } if (doprocessTriggerAnalysis) { - LOGF(info, "Trigger analysis is enabled. Desired trigger name = %s", cfg_swt_name.value); + LOGF(info, "Trigger analysis is enabled. Desired trigger name = %s", zorroGroup.cfg_swt_name.value); fRegistry.add("NormTrigger/hInspectedTVX", "inspected TVX;run number;N_{TVX}", kTProfile, {{80000, 520000.5, 600000.5}}, true); fRegistry.add("NormTrigger/hScalers", "trigger counter before DS;run number;counter", kTProfile, {{80000, 520000.5, 600000.5}}, true); fRegistry.add("NormTrigger/hSelections", "trigger counter after DS;run number;counter", kTProfile, {{80000, 520000.5, 600000.5}}, true); - auto hTriggerCounter = fRegistry.add("NormTrigger/hTriggerCounter", Form("trigger counter of %s;run number;", cfg_swt_name.value.data()), kTH2D, {{80000, 520000.5, 600000.5}, {2, -0.5, 1.5}}, false); + auto hTriggerCounter = fRegistry.add("NormTrigger/hTriggerCounter", Form("trigger counter of %s;run number;", zorroGroup.cfg_swt_name.value.data()), kTH2D, {{80000, 520000.5, 600000.5}, {2, -0.5, 1.5}}, false); hTriggerCounter->GetYaxis()->SetBinLabel(1, "Analyzed Trigger"); hTriggerCounter->GetYaxis()->SetBinLabel(2, "Analyzed TOI"); } @@ -496,6 +494,7 @@ struct DileptonHadronMPC { static constexpr std::string_view qvec_det_names[7] = {"FT0M", "FT0A", "FT0C", "BTot", "BPos", "BNeg", "FV0A"}; // event info o2::aod::pwgem::dilepton::utils::eventhistogram::addEventHistograms<-1>(&fRegistry); + fRegistry.add("Event/after/hMultNTracksPV_RefTracks", "Track mult. correlation;N_{track}^{PV};N_{track} for ref. flow", kTH2D, {{501, -0.5, 500.5}, {501, -0.5, 500.5}}, false); std::string mass_axis_title = "m_{ll} (GeV/c^{2})"; std::string pair_pt_axis_title = "p_{T,ll} (GeV/c)"; @@ -557,7 +556,8 @@ struct DileptonHadronMPC { fRegistry.add("DileptonHadron/same/uls/hs", "dilepton-hadron 2PC", kTHnSparseD, {axis_mass, axis_pt, axis_dca, axis_y, axis_deta, axis_dphi}, true); fRegistry.addClone("DileptonHadron/same/uls/", "DileptonHadron/same/lspp/"); fRegistry.addClone("DileptonHadron/same/uls/", "DileptonHadron/same/lsmm/"); - // fRegistry.addClone("DileptonHadron/same/", "DileptonHadron/mix/"); + fRegistry.addClone("DileptonHadron/same/", "DileptonHadron/mix/"); + fRegistry.add("DileptonHadron/mix/hDiffBC", "diff. global BC in mixed event;|BC_{current} - BC_{mixed}|", kTH1D, {{10001, -0.5, 10000.5}}, true); // hadron-hadron const AxisSpec axis_dphi_hh{90, -M_PI / 2, 3 * M_PI / 2, "#Delta#varphi = #varphi_{h}^{ref1} - #varphi_{h}^{ref2} (rad.)"}; @@ -701,13 +701,11 @@ struct DileptonHadronMPC { fEMTrackCut.SetTrackPtRange(trackcuts.cfg_min_pt_track, trackcuts.cfg_max_pt_track); fEMTrackCut.SetTrackEtaRange(trackcuts.cfg_min_eta_track, trackcuts.cfg_max_eta_track); // fEMTrackCut.SetTrackPhiRange(trackcuts.cfg_min_phi_track, trackcuts.cfg_max_phi_track); - // fEMTrackCut.SetTrackMaxDcaXY(trackcuts.cfg_max_dcaxy); - // fEMTrackCut.SetTrackMaxDcaZ(trackcuts.cfg_max_dcaz); fEMTrackCut.SetTrackBit(trackcuts.cfg_track_bits); } - template - bool fillDilepton(TCollision const& collision, TTrack1 const& t1, TTrack2 const& t2, TCut const& cut, TAllTracks const& tracks) + template + bool fillDilepton(TCollision const& collision, TTrack1 const& t1, TTrack2 const& t2, TCut const& cut) { if constexpr (ev_id == 0) { if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { @@ -725,12 +723,19 @@ struct DileptonHadronMPC { return false; } - if (!o2::aod::pwgem::dilepton::utils::emtrackutil::isBestMatch(t1, cut, tracks)) { + if (!map_best_match_globalmuon[t1.globalIndex()]) { return false; } - if (!o2::aod::pwgem::dilepton::utils::emtrackutil::isBestMatch(t2, cut, tracks)) { + if (!map_best_match_globalmuon[t2.globalIndex()]) { 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; + // } } } @@ -826,8 +831,8 @@ struct DileptonHadronMPC { return true; } - template - bool fillDileptonHadron(TTrack1 const& t1, TTrack2 const& t2, TCut const& cut, TAllTracks const& tracks, TRefTrack const& t3) + template + bool fillDileptonHadron(TTrack1 const& t1, TTrack2 const& t2, TCut const& cut, TRefTrack const& t3) { // this function must be called, if dilepton passes the cut. @@ -850,12 +855,19 @@ struct DileptonHadronMPC { return false; } - if (!o2::aod::pwgem::dilepton::utils::emtrackutil::isBestMatch(t1, cut, tracks)) { + if (!map_best_match_globalmuon[t1.globalIndex()]) { return false; } - if (!o2::aod::pwgem::dilepton::utils::emtrackutil::isBestMatch(t2, cut, tracks)) { + if (!map_best_match_globalmuon[t2.globalIndex()]) { 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 (!fEMTrackCut.IsSelected(t3)) { // for charged track @@ -874,11 +886,10 @@ struct DileptonHadronMPC { } float weight = 1.f; - if (cfgApplyWeightTTCA) { - weight = map_weight[std::make_pair(t1.globalIndex(), t2.globalIndex())]; - } - if (ev_id == 1) { - weight = 1.f; + if constexpr (ev_id == 0) { + if (cfgApplyWeightTTCA) { + weight = map_weight[std::make_pair(t1.globalIndex(), t2.globalIndex())]; + } } ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), leptonM1); @@ -926,7 +937,7 @@ struct DileptonHadronMPC { } template - bool fillHadronHadron(TRefTrack const& t1, TRefTrack const& t2, TLeptons const& posLeptons, TLeptons const& negLeptons, TLeptonCut const& cut) + bool fillHadronHadron(TRefTrack const& t1, TRefTrack const& t2, TLeptons const&, TLeptons const&, TLeptonCut const&) { if constexpr (ev_id == 0) { if (!fEMTrackCut.IsSelected(t1) || !fEMTrackCut.IsSelected(t2)) { // for charged track @@ -935,36 +946,37 @@ struct DileptonHadronMPC { // Leptons should not be in reference track sample. if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { - for (const auto& pos : posLeptons) { // leptons per collision - if (dielectroncuts.cfg_pid_scheme == static_cast(DielectronCut::PIDSchemes::kPIDML)) { - if (!cut.template IsSelectedTrack(pos)) { - continue; - } - } else { // cut based - if (!cut.template IsSelectedTrack(pos)) { - continue; - } - } - // if (t1.trackId() == pos.trackId() || t2.trackId() == pos.trackId()) { - // return false; - // } - } // end of pos lepton loop - - for (const auto& neg : negLeptons) { // leptons per collision - if (dielectroncuts.cfg_pid_scheme == static_cast(DielectronCut::PIDSchemes::kPIDML)) { - if (!cut.template IsSelectedTrack(neg)) { - continue; - } - } else { // cut based - if (!cut.template IsSelectedTrack(neg)) { - continue; - } - } - // if (t1.trackId() == neg.trackId() || t2.trackId() == neg.trackId()) { - // return false; - // } - } // end of neg lepton lopp - + if (trackcuts.cfg_reftrack_type == static_cast(o2::aod::pwgem::dilepton::utils::emtrackutil::RefTrackType::kCB)) { + // for (const auto& pos : posLeptons) { // leptons per collision + // if (dielectroncuts.cfg_pid_scheme == static_cast(DielectronCut::PIDSchemes::kPIDML)) { + // if (!cut.template IsSelectedTrack(pos)) { + // continue; + // } + // } else { // cut based + // if (!cut.template IsSelectedTrack(pos)) { + // continue; + // } + // } + // // if (t1.trackId() == pos.trackId() || t2.trackId() == pos.trackId()) { + // // return false; + // // } + // } // end of pos lepton loop + + // for (const auto& neg : negLeptons) { // leptons per collision + // if (dielectroncuts.cfg_pid_scheme == static_cast(DielectronCut::PIDSchemes::kPIDML)) { + // if (!cut.template IsSelectedTrack(neg)) { + // continue; + // } + // } else { // cut based + // if (!cut.template IsSelectedTrack(neg)) { + // continue; + // } + // } + // // if (t1.trackId() == neg.trackId() || t2.trackId() == neg.trackId()) { + // // return false; + // // } + // } // end of neg lepton lopp + } } // end of if kDielectron } // end of if same event @@ -1040,8 +1052,8 @@ struct DileptonHadronMPC { std::vector used_trackIds_per_col; int ndf = 0; - template - void run2PC(TCollisions const& collisions, TLeptons const& posTracks, TLeptons const& negTracks, TPresilce const& perCollision, TCut const& cut, TAllTracks const& tracks, TRefTracks const& refTracks) + template + void run2PC(TCollisions const& collisions, TLeptons const& posTracks, TLeptons const& negTracks, TPresilce const& perCollision, TCut const& cut, TRefTracks const& refTracks) { for (const auto& collision : collisions) { initCCDB(collision); @@ -1077,6 +1089,7 @@ struct DileptonHadronMPC { fRegistry.fill(HIST("Event/after/hEP2_CentFT0C_forMix"), collision.centFT0C(), ep2); auto refTracks_per_coll = refTracks.sliceBy(perCollision_track, collision.globalIndex()); + fRegistry.fill(HIST("Event/after/hMultNTracksPV_RefTracks"), collision.multNTracksPV(), refTracks_per_coll.size()); auto posTracks_per_coll = posTracks.sliceByCached(perCollision, collision.globalIndex(), cache); auto negTracks_per_coll = negTracks.sliceByCached(perCollision, collision.globalIndex(), cache); @@ -1084,29 +1097,29 @@ struct DileptonHadronMPC { int nuls = 0, nlspp = 0, nlsmm = 0; for (const auto& [pos, neg] : combinations(CombinationsFullIndexPolicy(posTracks_per_coll, negTracks_per_coll))) { // ULS - bool is_pair_ok = fillDilepton<0>(collision, pos, neg, cut, tracks); + bool is_pair_ok = fillDilepton<0>(collision, pos, neg, cut); if (is_pair_ok) { nuls++; for (const auto& refTrack : refTracks_per_coll) { - fillDileptonHadron<0>(pos, neg, cut, tracks, refTrack); + fillDileptonHadron<0>(pos, neg, cut, refTrack); } } } for (const auto& [pos1, pos2] : combinations(CombinationsStrictlyUpperIndexPolicy(posTracks_per_coll, posTracks_per_coll))) { // LS++ - bool is_pair_ok = fillDilepton<0>(collision, pos1, pos2, cut, tracks); + bool is_pair_ok = fillDilepton<0>(collision, pos1, pos2, cut); if (is_pair_ok) { nlspp++; for (const auto& refTrack : refTracks_per_coll) { - fillDileptonHadron<0>(pos1, pos2, cut, tracks, refTrack); + fillDileptonHadron<0>(pos1, pos2, cut, refTrack); } } } for (const auto& [neg1, neg2] : combinations(CombinationsStrictlyUpperIndexPolicy(negTracks_per_coll, negTracks_per_coll))) { // LS-- - bool is_pair_ok = fillDilepton<0>(collision, neg1, neg2, cut, tracks); + bool is_pair_ok = fillDilepton<0>(collision, neg1, neg2, cut); if (is_pair_ok) { nlsmm++; for (const auto& refTrack : refTracks_per_coll) { - fillDileptonHadron<0>(neg1, neg2, cut, tracks, refTrack); + fillDileptonHadron<0>(neg1, neg2, cut, refTrack); } } } @@ -1128,7 +1141,10 @@ struct DileptonHadronMPC { } // store ref tracks } } + // LOGF(info, "collision.globalIndex() = %d, collision.centFT0M() = %f, refTracks_per_coll.size() = %d", collision.globalIndex(), collision.centFT0M(), refTracks_per_coll.size()); + for (const auto& [ref1, ref2] : combinations(CombinationsStrictlyUpperIndexPolicy(refTracks_per_coll, refTracks_per_coll))) { + // TODO: remove lepton candidates from reference track sample in case of kCB. fillHadronHadron<0>(ref1, ref2, posTracks_per_coll, negTracks_per_coll, cut); } } @@ -1206,25 +1222,25 @@ struct DileptonHadronMPC { for (const auto& pos : selected_posTracks_in_this_event) { // ULS mix for (const auto& neg : negTracks_from_event_pool) { - fillDilepton<1>(collision, pos, neg, cut, nullptr); + fillDilepton<1>(collision, pos, neg, cut); } } for (const auto& neg : selected_negTracks_in_this_event) { // ULS mix for (const auto& pos : posTracks_from_event_pool) { - fillDilepton<1>(collision, neg, pos, cut, nullptr); + fillDilepton<1>(collision, neg, pos, cut); } } for (const auto& pos1 : selected_posTracks_in_this_event) { // LS++ mix for (const auto& pos2 : posTracks_from_event_pool) { - fillDilepton<1>(collision, pos1, pos2, cut, nullptr); + fillDilepton<1>(collision, pos1, pos2, cut); } } for (const auto& neg1 : selected_negTracks_in_this_event) { // LS-- mix for (const auto& neg2 : negTracks_from_event_pool) { - fillDilepton<1>(collision, neg1, neg2, cut, nullptr); + fillDilepton<1>(collision, neg1, neg2, cut); } } } // end of loop over mixed event pool for lepton-lepton @@ -1233,6 +1249,89 @@ struct DileptonHadronMPC { auto selected_refTracks_in_this_event = emh_ref->GetTracksPerCollision(key_df_collision); auto collisionIds_in_mixing_pool_hadron = emh_ref->GetCollisionIdsFromEventPool(key_bin); + // for ULS and hadron mix + for (const auto& pos : selected_posTracks_in_this_event) { + for (const auto& neg : selected_negTracks_in_this_event) { + + for (const auto& mix_dfId_collisionId : collisionIds_in_mixing_pool_hadron) { + int mix_dfId = mix_dfId_collisionId.first; + int mix_collisionId = mix_dfId_collisionId.second; + if (collision.globalIndex() == mix_collisionId && ndf == mix_dfId) { // this never happens. only protection. + continue; + } + auto globalBC_mix = map_mixed_eventId_to_globalBC[mix_dfId_collisionId]; + uint64_t diffBC = std::max(collision.globalBC(), globalBC_mix) - std::min(collision.globalBC(), globalBC_mix); + fRegistry.fill(HIST("DileptonHadron/mix/hDiffBC"), diffBC); + if (diffBC < ndiff_bc_mix) { + continue; + } + + auto refTracks_from_event_pool = emh_ref->GetTracksPerCollision(mix_dfId_collisionId); + for (const auto& ref : refTracks_from_event_pool) { + fillDileptonHadron<1>(pos, neg, cut, ref); + } + + } // end of loop over mixed event pool for dilepton-hadron + } + } + + // for LS++ and hadron mix + for (size_t i1 = 0; i1 < selected_posTracks_in_this_event.size(); i1++) { + auto pos1 = selected_posTracks_in_this_event[i1]; + for (size_t i2 = i1 + 1; i2 < selected_posTracks_in_this_event.size(); i2++) { + auto pos2 = selected_posTracks_in_this_event[i2]; + + for (const auto& mix_dfId_collisionId : collisionIds_in_mixing_pool_hadron) { + int mix_dfId = mix_dfId_collisionId.first; + int mix_collisionId = mix_dfId_collisionId.second; + if (collision.globalIndex() == mix_collisionId && ndf == mix_dfId) { // this never happens. only protection. + continue; + } + + auto globalBC_mix = map_mixed_eventId_to_globalBC[mix_dfId_collisionId]; + uint64_t diffBC = std::max(collision.globalBC(), globalBC_mix) - std::min(collision.globalBC(), globalBC_mix); + fRegistry.fill(HIST("DileptonHadron/mix/hDiffBC"), diffBC); + if (diffBC < ndiff_bc_mix) { + continue; + } + + auto refTracks_from_event_pool = emh_ref->GetTracksPerCollision(mix_dfId_collisionId); + for (const auto& ref : refTracks_from_event_pool) { + fillDileptonHadron<1>(pos1, pos2, cut, ref); + } + + } // end of loop over mixed event pool for dilepton-hadron + } + } + + // for LS-- and hadron mix + for (size_t i1 = 0; i1 < selected_negTracks_in_this_event.size(); i1++) { + auto neg1 = selected_negTracks_in_this_event[i1]; + for (size_t i2 = i1 + 1; i2 < selected_negTracks_in_this_event.size(); i2++) { + auto neg2 = selected_negTracks_in_this_event[i2]; + + for (const auto& mix_dfId_collisionId : collisionIds_in_mixing_pool_hadron) { + int mix_dfId = mix_dfId_collisionId.first; + int mix_collisionId = mix_dfId_collisionId.second; + if (collision.globalIndex() == mix_collisionId && ndf == mix_dfId) { // this never happens. only protection. + continue; + } + auto globalBC_mix = map_mixed_eventId_to_globalBC[mix_dfId_collisionId]; + uint64_t diffBC = std::max(collision.globalBC(), globalBC_mix) - std::min(collision.globalBC(), globalBC_mix); + fRegistry.fill(HIST("DileptonHadron/mix/hDiffBC"), diffBC); + if (diffBC < ndiff_bc_mix) { + continue; + } + + auto refTracks_from_event_pool = emh_ref->GetTracksPerCollision(mix_dfId_collisionId); + for (const auto& ref : refTracks_from_event_pool) { + fillDileptonHadron<1>(neg1, neg2, cut, ref); + } + + } // end of loop over mixed event pool for dilepton-hadron + } + } + for (const auto& mix_dfId_collisionId : collisionIds_in_mixing_pool_hadron) { int mix_dfId = mix_dfId_collisionId.first; int mix_collisionId = mix_dfId_collisionId.second; @@ -1269,8 +1368,8 @@ struct DileptonHadronMPC { } // end of DF - template - bool isPairOK(TTrack1 const& t1, TTrack2 const& t2, TCut const& cut, TAllTracks const& tracks) + template + bool isPairOK(TTrack1 const& t1, TTrack2 const& t2, TCut const& cut) { if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { if (dielectroncuts.cfg_pid_scheme == static_cast(DielectronCut::PIDSchemes::kPIDML)) { @@ -1287,12 +1386,19 @@ struct DileptonHadronMPC { return false; } - if (!o2::aod::pwgem::dilepton::utils::emtrackutil::isBestMatch(t1, cut, tracks)) { + if (!map_best_match_globalmuon[t1.globalIndex()]) { return false; } - if (!o2::aod::pwgem::dilepton::utils::emtrackutil::isBestMatch(t2, cut, tracks)) { + if (!map_best_match_globalmuon[t2.globalIndex()]) { 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) { @@ -1332,17 +1438,17 @@ struct DileptonHadronMPC { auto negTracks_per_coll = negTracks.sliceByCached(perCollision, collision.globalIndex(), cache); for (const auto& [pos, neg] : combinations(CombinationsFullIndexPolicy(posTracks_per_coll, negTracks_per_coll))) { // ULS - if (isPairOK(pos, neg, cut, tracks)) { + if (isPairOK(pos, neg, cut)) { passed_pairIds.emplace_back(std::make_pair(pos.globalIndex(), neg.globalIndex())); } } for (const auto& [pos1, pos2] : combinations(CombinationsStrictlyUpperIndexPolicy(posTracks_per_coll, posTracks_per_coll))) { // LS++ - if (isPairOK(pos1, pos2, cut, tracks)) { + if (isPairOK(pos1, pos2, cut)) { passed_pairIds.emplace_back(std::make_pair(pos1.globalIndex(), pos2.globalIndex())); } } for (const auto& [neg1, neg2] : combinations(CombinationsStrictlyUpperIndexPolicy(negTracks_per_coll, negTracks_per_coll))) { // LS-- - if (isPairOK(neg1, neg2, cut, tracks)) { + if (isPairOK(neg1, neg2, cut)) { passed_pairIds.emplace_back(std::make_pair(neg1.globalIndex(), neg2.globalIndex())); } } @@ -1383,6 +1489,8 @@ struct DileptonHadronMPC { passed_pairIds.shrink_to_fit(); } + std::unordered_map map_best_match_globalmuon; + void processAnalysis(FilteredMyCollisions const& collisions, FilteredRefTracks const& refTracks, Types const&... args) { if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) { @@ -1390,15 +1498,17 @@ struct DileptonHadronMPC { if (cfgApplyWeightTTCA) { fillPairWeightMap(collisions, positive_electrons, negative_electrons, o2::aod::emprimaryelectron::emeventId, fDielectronCut, electrons); } - run2PC(collisions, positive_electrons, negative_electrons, o2::aod::emprimaryelectron::emeventId, fDielectronCut, electrons, refTracks); + run2PC(collisions, positive_electrons, negative_electrons, o2::aod::emprimaryelectron::emeventId, fDielectronCut, refTracks); } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { auto muons = std::get<0>(std::tie(args...)); + map_best_match_globalmuon = findBestMatchMap(muons, fDimuonCut); if (cfgApplyWeightTTCA) { fillPairWeightMap(collisions, positive_muons, negative_muons, o2::aod::emprimarymuon::emeventId, fDimuonCut, muons); } - run2PC(collisions, positive_muons, negative_muons, o2::aod::emprimarymuon::emeventId, fDimuonCut, muons, refTracks); + run2PC(collisions, positive_muons, negative_muons, o2::aod::emprimarymuon::emeventId, fDimuonCut, refTracks); } map_weight.clear(); + map_best_match_globalmuon.clear(); ndf++; } PROCESS_SWITCH(DileptonHadronMPC, processAnalysis, "run dilepton analysis", true); @@ -1410,15 +1520,17 @@ struct DileptonHadronMPC { if (cfgApplyWeightTTCA) { fillPairWeightMap(collisions, positive_electrons, negative_electrons, o2::aod::emprimaryelectron::emeventId, fDielectronCut, electrons); } - run2PC(collisions, positive_electrons, negative_electrons, o2::aod::emprimaryelectron::emeventId, fDielectronCut, electrons, refTracks); + run2PC(collisions, positive_electrons, negative_electrons, o2::aod::emprimaryelectron::emeventId, fDielectronCut, refTracks); } else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) { auto muons = std::get<0>(std::tie(args...)); + map_best_match_globalmuon = findBestMatchMap(muons, fDimuonCut); if (cfgApplyWeightTTCA) { fillPairWeightMap(collisions, positive_muons, negative_muons, o2::aod::emprimarymuon::emeventId, fDimuonCut, muons); } - run2PC(collisions, positive_muons, negative_muons, o2::aod::emprimarymuon::emeventId, fDimuonCut, muons, refTracks); + run2PC(collisions, positive_muons, negative_muons, o2::aod::emprimarymuon::emeventId, fDimuonCut, refTracks); } map_weight.clear(); + map_best_match_globalmuon.clear(); ndf++; } PROCESS_SWITCH(DileptonHadronMPC, processTriggerAnalysis, "run dilepton analysis on triggered data", false); diff --git a/PWGEM/Dilepton/Core/EMTrackCut.cxx b/PWGEM/Dilepton/Core/EMTrackCut.cxx index ed9e345a0f9..cd995e08836 100644 --- a/PWGEM/Dilepton/Core/EMTrackCut.cxx +++ b/PWGEM/Dilepton/Core/EMTrackCut.cxx @@ -38,12 +38,6 @@ void EMTrackCut::SetTrackPhiRange(float minPhi, float maxPhi) LOG(info) << "EMTrack Cut, set track phi range (rad.): " << mMinTrackPhi << " - " << mMaxTrackPhi; } -// void EMTrackCut::SetTrackMaxDcaXYPtDep(std::function ptDepCut) -// { -// mMaxDcaXYPtDep = ptDepCut; -// LOG(info) << "EMTrack Cut, set max DCA xy pt dep: " << mMaxDcaXYPtDep(1.0); -// } - void EMTrackCut::SetTrackBit(uint16_t bit) { mTrackBit = bit; diff --git a/PWGEM/Dilepton/Core/EMTrackCut.h b/PWGEM/Dilepton/Core/EMTrackCut.h index 2194aabb7cf..118c462cbf7 100644 --- a/PWGEM/Dilepton/Core/EMTrackCut.h +++ b/PWGEM/Dilepton/Core/EMTrackCut.h @@ -42,9 +42,11 @@ class EMTrackCut : public TNamed if (!IsSelectedTrack(track, EMTrackCuts::kTrackPtRange)) { return false; } + if (!IsSelectedTrack(track, EMTrackCuts::kTrackEtaRange)) { return false; } + if (!IsSelectedTrack(track, EMTrackCuts::kTrackPhiRange)) { return false; } @@ -69,15 +71,8 @@ class EMTrackCut : public TNamed case EMTrackCuts::kTrackPhiRange: return track.phi() > mMinTrackPhi && track.phi() < mMaxTrackPhi; - case EMTrackCuts::kTrackBit: { - // for (int i = 0; i < 10; i++) { - // if ((mTrackBit & (1 << i)) > 0 && !((track.trackBit() & (1 << i)) > 0)) { - // return false; - // } - // } - // return true; + case EMTrackCuts::kTrackBit: return (track.trackBit() & mTrackBit) >= mTrackBit; - } default: return false; @@ -91,14 +86,10 @@ class EMTrackCut : public TNamed void SetTrackBit(uint16_t bits); private: - // kinematic cuts float mMinTrackPt{0.f}, mMaxTrackPt{1e10f}; // range in pT float mMinTrackEta{-1e10f}, mMaxTrackEta{1e10f}; // range in eta float mMinTrackPhi{0.f}, mMaxTrackPhi{6.3}; // range in phi - - // track quality cuts uint16_t mTrackBit{0}; - // std::function mMaxDcaXYPtDep{}; // max dca in xy plane as function of pT ClassDef(EMTrackCut, 1); }; diff --git a/PWGEM/Dilepton/TableProducer/skimmerPrimaryMFTTrack.cxx b/PWGEM/Dilepton/TableProducer/skimmerPrimaryMFTTrack.cxx index 83835874d87..ad27348ba7d 100644 --- a/PWGEM/Dilepton/TableProducer/skimmerPrimaryMFTTrack.cxx +++ b/PWGEM/Dilepton/TableProducer/skimmerPrimaryMFTTrack.cxx @@ -15,6 +15,7 @@ #include "PWGEM/Dilepton/DataModel/dileptonTables.h" #include "PWGEM/Dilepton/Utils/EMTrackUtilities.h" +#include "Common/Core/RecoDecay.h" #include "Common/Core/TableHelper.h" #include "Common/Core/fwdtrackUtilities.h" @@ -183,7 +184,7 @@ struct skimmerPrimaryMFTTrack { float pt = trackPar.getPt(); float eta = trackPar.getEta(); float phi = trackPar.getPhi(); - o2::math_utils::bringTo02Pi(phi); + phi = RecoDecay::constrainAngle(phi, 0, 1U); if (pt < cfgPtMin || cfgPtMax < pt) { return; @@ -205,36 +206,36 @@ struct skimmerPrimaryMFTTrack { } if (mfttrack.nClusters() >= 6) { - trackBit |= static_cast(RefMFTTrackBit::kNclsMFT6); + SETBIT(trackBit, static_cast(o2::aod::pwgem::dilepton::utils::emtrackutil::RefMFTTrackBit::kNclsMFT6)); } if (mfttrack.nClusters() >= 7) { - trackBit |= static_cast(RefMFTTrackBit::kNclsMFT7); + SETBIT(trackBit, static_cast(o2::aod::pwgem::dilepton::utils::emtrackutil::RefMFTTrackBit::kNclsMFT7)); } if (mfttrack.nClusters() >= 8) { - trackBit |= static_cast(RefMFTTrackBit::kNclsMFT8); + SETBIT(trackBit, static_cast(o2::aod::pwgem::dilepton::utils::emtrackutil::RefMFTTrackBit::kNclsMFT8)); } if (mfttrack.chi2() / ndf < 3.f) { - trackBit |= static_cast(RefMFTTrackBit::kChi2MFT3); + SETBIT(trackBit, static_cast(o2::aod::pwgem::dilepton::utils::emtrackutil::RefMFTTrackBit::kChi2MFT3)); } if (mfttrack.chi2() / ndf < 2.f) { - trackBit |= static_cast(RefMFTTrackBit::kChi2MFT2); + SETBIT(trackBit, static_cast(o2::aod::pwgem::dilepton::utils::emtrackutil::RefMFTTrackBit::kChi2MFT2)); } if (std::fabs(dcaXY) < 0.05) { - trackBit |= static_cast(RefMFTTrackBit::kDCAxy005cm); + SETBIT(trackBit, static_cast(o2::aod::pwgem::dilepton::utils::emtrackutil::RefMFTTrackBit::kDCAxy005cm)); } if (std::fabs(dcaXY) < 0.04) { - trackBit |= static_cast(RefMFTTrackBit::kDCAxy004cm); + SETBIT(trackBit, static_cast(o2::aod::pwgem::dilepton::utils::emtrackutil::RefMFTTrackBit::kDCAxy004cm)); } if (std::fabs(dcaXY) < 0.03) { - trackBit |= static_cast(RefMFTTrackBit::kDCAxy003cm); + SETBIT(trackBit, static_cast(o2::aod::pwgem::dilepton::utils::emtrackutil::RefMFTTrackBit::kDCAxy003cm)); } if (std::fabs(dcaXY) < 0.02) { - trackBit |= static_cast(RefMFTTrackBit::kDCAxy002cm); + SETBIT(trackBit, static_cast(o2::aod::pwgem::dilepton::utils::emtrackutil::RefMFTTrackBit::kDCAxy002cm)); } if (std::fabs(dcaXY) < 0.01) { - trackBit |= static_cast(RefMFTTrackBit::kDCAxy001cm); + SETBIT(trackBit, static_cast(o2::aod::pwgem::dilepton::utils::emtrackutil::RefMFTTrackBit::kDCAxy001cm)); } emprimarytracks(mfttrack.sign() / pt, eta, phi, trackBit); diff --git a/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuon.cxx b/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuon.cxx index bba65c5e5d7..a64ae751244 100644 --- a/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuon.cxx +++ b/PWGEM/Dilepton/TableProducer/skimmerPrimaryMuon.cxx @@ -77,18 +77,18 @@ struct skimmerPrimaryMuon { 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"}; - 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 minEtaGL{"minEtaGL", -1e+10, "min. eta acceptance for MFT-MCH-MID"}; + Configurable maxEtaGL{"maxEtaGL", 0, "max. eta acceptance for MFT-MCH-MID"}; Configurable minRabsGL{"minRabsGL", 17.6, "min. R at absorber end for global muon (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"}; - Configurable maxDCAxy{"maxDCAxy", 0.2, "max. DCAxy for global muons"}; + Configurable maxDCAxy{"maxDCAxy", 1e+10, "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 maxMatchingChi2MCHMFT{"maxMatchingChi2MCHMFT", 1e+10, "max. chi2 for MCH-MFT matching"}; Configurable maxChi2SA{"maxChi2SA", 1e+10, "max. chi2 for standalone muon"}; - Configurable maxChi2GL{"maxChi2GL", 12, "max. chi2 for global muon"}; + Configurable maxChi2GL{"maxChi2GL", 1e+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"}; Configurable minNmuon{"minNmuon", 0, "min number of muon candidates per collision"}; @@ -182,10 +182,10 @@ struct skimmerPrimaryMuon { fRegistry.add("MFTMCHMID/hNclustersMFT", "NclustersMFT;Nclusters MFT", kTH1F, {{11, -0.5f, 10.5}}, false); fRegistry.add("MFTMCHMID/hRatAbsorberEnd", "R at absorber end;R at absorber end (cm)", kTH1F, {{100, 0.0f, 100}}, false); fRegistry.add("MFTMCHMID/hPDCA_Rabs", "pDCA vs. Rabs;R at absorber end (cm);p #times DCA (GeV/c #upoint cm)", kTH2F, {{100, 0, 100}, {100, 0.0f, 1000}}, false); - fRegistry.add("MFTMCHMID/hChi2", "chi2;chi2/ndf", kTH1F, {{200, 0.0f, 20}}, false); - fRegistry.add("MFTMCHMID/hChi2MFT", "chi2 MFT;chi2 MFT/ndf", kTH1F, {{200, 0.0f, 20}}, false); - fRegistry.add("MFTMCHMID/hChi2MatchMCHMID", "chi2 match MCH-MID;chi2", kTH1F, {{200, 0.0f, 20}}, false); - fRegistry.add("MFTMCHMID/hChi2MatchMCHMFT", "chi2 match MCH-MFT;chi2", kTH1F, {{200, 0.0f, 100}}, false); + fRegistry.add("MFTMCHMID/hChi2_Pt", "chi2;p_{T} (GeV/c);chi2/ndf", kTH2F, {{100, 0, 10}, {200, 0.0f, 20}}, false); + fRegistry.add("MFTMCHMID/hChi2MFT_Pt", "chi2 MFT;p_{T} (GeV/c);chi2 MFT/ndf", kTH2F, {{100, 0, 10}, {200, 0.0f, 20}}, false); + fRegistry.add("MFTMCHMID/hChi2MatchMCHMID_Pt", "chi2 match MCH-MID;p_{T} (GeV/c);matching chi2 MCH-MID", kTH2F, {{100, 0, 10}, {200, 0.0f, 20}}, false); + fRegistry.add("MFTMCHMID/hChi2MatchMCHMFT_Pt", "chi2 match MCH-MFT;p_{T} (GeV/c);matching chi2 MCH-MFT", kTH2F, {{100, 0, 10}, {200, 0.0f, 100}}, false); fRegistry.add("MFTMCHMID/hDCAxy2D", "DCA x vs. y;DCA_{x} (cm);DCA_{y} (cm)", kTH2F, {{200, -1, 1}, {200, -1, +1}}, false); fRegistry.add("MFTMCHMID/hDCAxy2DinSigma", "DCA x vs. y in sigma;DCA_{x} (#sigma);DCA_{y} (#sigma)", kTH2F, {{200, -10, 10}, {200, -10, +10}}, false); fRegistry.add("MFTMCHMID/hDCAxy", "DCAxy;DCA_{xy} (cm);", kTH1F, {{100, 0, 1}}, false); @@ -491,10 +491,10 @@ struct skimmerPrimaryMuon { 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 / ndf_mft); - fRegistry.fill(HIST("MFTMCHMID/hChi2MatchMCHMID"), fwdtrack.chi2MatchMCHMID()); - fRegistry.fill(HIST("MFTMCHMID/hChi2MatchMCHMFT"), fwdtrack.chi2MatchMCHMFT()); + fRegistry.fill(HIST("MFTMCHMID/hChi2_Pt"), pt, fwdtrack.chi2() / ndf_mchmft); + fRegistry.fill(HIST("MFTMCHMID/hChi2MFT_Pt"), pt, chi2mft / ndf_mft); + fRegistry.fill(HIST("MFTMCHMID/hChi2MatchMCHMID_Pt"), pt, fwdtrack.chi2MatchMCHMID()); + fRegistry.fill(HIST("MFTMCHMID/hChi2MatchMCHMFT_Pt"), pt, fwdtrack.chi2MatchMCHMFT()); fRegistry.fill(HIST("MFTMCHMID/hDCAxy2D"), dcaX, dcaY); fRegistry.fill(HIST("MFTMCHMID/hDCAxy2DinSigma"), dcaX / std::sqrt(cXX), dcaY / std::sqrt(cYY)); fRegistry.fill(HIST("MFTMCHMID/hDCAxy"), dcaXY); @@ -519,10 +519,10 @@ struct skimmerPrimaryMuon { 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 / ndf_mft); - fRegistry.fill(HIST("MCHMID/hChi2MatchMCHMID"), fwdtrack.chi2MatchMCHMID()); - fRegistry.fill(HIST("MCHMID/hChi2MatchMCHMFT"), fwdtrack.chi2MatchMCHMFT()); + fRegistry.fill(HIST("MCHMID/hChi2_Pt"), pt, fwdtrack.chi2()); + fRegistry.fill(HIST("MCHMID/hChi2MFT_Pt"), pt, chi2mft / ndf_mft); + fRegistry.fill(HIST("MCHMID/hChi2MatchMCHMID_Pt"), pt, fwdtrack.chi2MatchMCHMID()); + fRegistry.fill(HIST("MCHMID/hChi2MatchMCHMFT_Pt"), pt, fwdtrack.chi2MatchMCHMFT()); fRegistry.fill(HIST("MCHMID/hDCAxy2D"), dcaX, dcaY); fRegistry.fill(HIST("MCHMID/hDCAxy2DinSigma"), dcaX / std::sqrt(cXX), dcaY / std::sqrt(cYY)); fRegistry.fill(HIST("MCHMID/hDCAxy"), dcaXY); @@ -622,38 +622,6 @@ struct skimmerPrimaryMuon { // 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); } - // // PresliceUnsorted perMFTTrack = o2::aod::fwdtrack::matchMFTTrackId; - // template - // bool isBestMatch(TFwdTrack const& fwdtrack, TFwdTracks const& fwdtracks, TMFTTracks const&) - // { - // if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { - // std::map map_chi2MCHMFT; - // map_chi2MCHMFT[fwdtrack.globalIndex()] = fwdtrack.chi2MatchMCHMFT(); // add myself - // // LOGF(info, "add myself: fwdtrack.globalIndex() = %d, fwdtrack.chi2MatchMCHMFT() = %f", fwdtrack.globalIndex(), fwdtrack.chi2MatchMCHMFT()); - - // auto glMuonsPerMFT = std::views::filter(vec_min_chi2MatchMCHMFT, [&](std::tuple t) { return std::get<2>(t) == fwdtrack.matchMFTTrackId() && std::get<1>(t) != fwdtrack.matchMCHTrackId() && std::get<0>(t) != fwdtrack.globalIndex(); }); - // for (const auto& candidate : glMuonsPerMFT) { - // map_chi2MCHMFT[std::get<0>(candidate)] = fwdtracks.rawIteratorAt(std::get<0>(candidate)).chi2MatchMCHMFT(); - // // LOGF(info, "same MFT found: candidate.globalIndex() = %d, candidate.chi2MatchMCHMFT() = %f", std::get<0>(candidate), fwdtracks.rawIteratorAt(std::get<0>(candidate)).chi2MatchMCHMFT()); - // } - - // auto it = std::min_element(map_chi2MCHMFT.begin(), map_chi2MCHMFT.end(), [](decltype(map_chi2MCHMFT)::value_type& l, decltype(map_chi2MCHMFT)::value_type& r) -> bool { return l.second < r.second; }); - - // // LOGF(info, "min: globalIndex = %d, chi2 = %f", it->first, it->second); - // // LOGF(info, "bool = %d", it->first == fwdtrack.globalIndex()); - - // if (it->first == fwdtrack.globalIndex()) { // search for minimum matching-chi2 - // map_chi2MCHMFT.clear(); - // return true; - // } else { - // map_chi2MCHMFT.clear(); - // return false; - // } - // } else { - // return true; - // } - // } - SliceCache cache; Preslice perCollision = o2::aod::fwdtrack::collisionId; Preslice fwdtrackIndicesPerCollision = aod::track_association::collisionId; @@ -693,10 +661,6 @@ struct skimmerPrimaryMuon { // continue; // } - // if (!isBestMatch(fwdtrack, fwdtracks, mfttracks)) { - // continue; - // } - if (!fillFwdTrackTable(collision, fwdtrack, nullptr, false)) { continue; } @@ -783,10 +747,6 @@ struct skimmerPrimaryMuon { // continue; // } - // if (!isBestMatch(fwdtrack, fwdtracks, mfttracks)) { - // continue; - // } - if (!fillFwdTrackTable(collision, fwdtrack, nullptr, mapAmb[fwdtrack.globalIndex()])) { continue; } @@ -833,98 +793,94 @@ struct skimmerPrimaryMuon { } PROCESS_SWITCH(skimmerPrimaryMuon, processRec_TTCA, "process reconstructed info", false); - 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(); - } - - vec_min_chi2MatchMCHMFT.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) { - 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; - // } - - // if (!isBestMatch(fwdtrack, fwdtracks, mfttracks)) { - // continue; - // } - - if (!fillFwdTrackTable(collision, fwdtrack, mftCovs, mapAmb[fwdtrack.globalIndex()])) { - continue; - } - - if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { - multiMapGLMuonsPerCollision.insert(std::make_pair(collision.globalIndex(), fwdtrack.globalIndex())); - } - if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { - multiMapSAMuonsPerCollision.insert(std::make_pair(collision.globalIndex(), fwdtrack.globalIndex())); - } - - } // end of fwdtrack loop - } // end of collision loop + // 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(); + // } - for (const auto& collision : collisions) { - int count_samuons = multiMapSAMuonsPerCollision.count(collision.globalIndex()); - int count_glmuons = multiMapGLMuonsPerCollision.count(collision.globalIndex()); - if (fillQAHistograms) { - fRegistry.fill(HIST("MCHMID/hNmu"), count_samuons); - fRegistry.fill(HIST("MFTMCHMID/hNmu"), count_glmuons); - } - if (count_samuons >= minNmuon) { - auto range_samuons = multiMapSAMuonsPerCollision.equal_range(collision.globalIndex()); - for (auto it = range_samuons.first; it != range_samuons.second; it++) { - auto fwdtrack = fwdtracks.rawIteratorAt(it->second); - fillFwdTrackTable(collision, fwdtrack, mftCovs, mapAmb[fwdtrack.globalIndex()]); - } - } - if (count_glmuons >= minNmuon) { - auto range_glmuons = multiMapGLMuonsPerCollision.equal_range(collision.globalIndex()); - for (auto it = range_glmuons.first; it != range_glmuons.second; it++) { - auto fwdtrack = fwdtracks.rawIteratorAt(it->second); - fillFwdTrackTable(collision, fwdtrack, mftCovs, mapAmb[fwdtrack.globalIndex()]); - } - } - } // end of collision loop + // vec_min_chi2MatchMCHMFT.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) { + // 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; + // } - multiMapSAMuonsPerCollision.clear(); - multiMapGLMuonsPerCollision.clear(); - 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); + // 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; + // // } + + // if (!fillFwdTrackTable(collision, fwdtrack, mftCovs, mapAmb[fwdtrack.globalIndex()])) { + // continue; + // } + + // if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { + // multiMapGLMuonsPerCollision.insert(std::make_pair(collision.globalIndex(), fwdtrack.globalIndex())); + // } + // if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { + // multiMapSAMuonsPerCollision.insert(std::make_pair(collision.globalIndex(), fwdtrack.globalIndex())); + // } + + // } // end of fwdtrack loop + // } // end of collision loop + + // for (const auto& collision : collisions) { + // int count_samuons = multiMapSAMuonsPerCollision.count(collision.globalIndex()); + // int count_glmuons = multiMapGLMuonsPerCollision.count(collision.globalIndex()); + // if (fillQAHistograms) { + // fRegistry.fill(HIST("MCHMID/hNmu"), count_samuons); + // fRegistry.fill(HIST("MFTMCHMID/hNmu"), count_glmuons); + // } + // if (count_samuons >= minNmuon) { + // auto range_samuons = multiMapSAMuonsPerCollision.equal_range(collision.globalIndex()); + // for (auto it = range_samuons.first; it != range_samuons.second; it++) { + // auto fwdtrack = fwdtracks.rawIteratorAt(it->second); + // fillFwdTrackTable(collision, fwdtrack, mftCovs, mapAmb[fwdtrack.globalIndex()]); + // } + // } + // if (count_glmuons >= minNmuon) { + // auto range_glmuons = multiMapGLMuonsPerCollision.equal_range(collision.globalIndex()); + // for (auto it = range_glmuons.first; it != range_glmuons.second; it++) { + // auto fwdtrack = fwdtracks.rawIteratorAt(it->second); + // fillFwdTrackTable(collision, fwdtrack, mftCovs, mapAmb[fwdtrack.globalIndex()]); + // } + // } + // } // end of collision loop + + // multiMapSAMuonsPerCollision.clear(); + // multiMapGLMuonsPerCollision.clear(); + // 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& mfttracks, aod::BCsWithTimestamps const&) { @@ -959,10 +915,6 @@ struct skimmerPrimaryMuon { // continue; // } - // if (!isBestMatch(fwdtrack, fwdtracks, mfttracks)) { - // continue; - // } - if (!fillFwdTrackTable(collision, fwdtrack, nullptr, false)) { continue; } @@ -1048,10 +1000,6 @@ struct skimmerPrimaryMuon { // continue; // } - // if (!isBestMatch(fwdtrack, fwdtracks, mfttracks)) { - // continue; - // } - if (!fillFwdTrackTable(collision, fwdtrack, nullptr, mapAmb[fwdtrack.globalIndex()])) { continue; } @@ -1098,97 +1046,94 @@ struct skimmerPrimaryMuon { } PROCESS_SWITCH(skimmerPrimaryMuon, processRec_TTCA_SWT, "process reconstructed info", false); - 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(); - } - vec_min_chi2MatchMCHMFT.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) { - 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) { - auto bc = collision.template bc_as(); - initCCDB(bc); - if (!collision.isSelected()) { - continue; - } - if (collision.swtaliastmp_raw() == 0) { - continue; - } - - auto fwdtrackIdsThisCollision = fwdtrackIndices.sliceBy(fwdtrackIndicesPerCollision, collision.globalIndex()); - for (const auto& fwdtrackId : fwdtrackIdsThisCollision) { - 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; - // } - // if (!isBestMatch(fwdtrack, fwdtracks, mfttracks)) { - // continue; - // } - - if (!fillFwdTrackTable(collision, fwdtrack, mftCovs, mapAmb[fwdtrack.globalIndex()])) { - continue; - } - - if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { - multiMapGLMuonsPerCollision.insert(std::make_pair(collision.globalIndex(), fwdtrack.globalIndex())); - } - if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { - multiMapSAMuonsPerCollision.insert(std::make_pair(collision.globalIndex(), fwdtrack.globalIndex())); - } - - } // end of fwdtrack loop - } // end of collision loop - - for (const auto& collision : collisions) { - int count_samuons = multiMapSAMuonsPerCollision.count(collision.globalIndex()); - int count_glmuons = multiMapGLMuonsPerCollision.count(collision.globalIndex()); - if (fillQAHistograms) { - fRegistry.fill(HIST("MCHMID/hNmu"), count_samuons); - fRegistry.fill(HIST("MFTMCHMID/hNmu"), count_glmuons); - } - if (count_samuons >= minNmuon) { - auto range_samuons = multiMapSAMuonsPerCollision.equal_range(collision.globalIndex()); - for (auto it = range_samuons.first; it != range_samuons.second; it++) { - auto fwdtrack = fwdtracks.rawIteratorAt(it->second); - fillFwdTrackTable(collision, fwdtrack, mftCovs, mapAmb[fwdtrack.globalIndex()]); - } - } - if (count_glmuons >= minNmuon) { - auto range_glmuons = multiMapGLMuonsPerCollision.equal_range(collision.globalIndex()); - for (auto it = range_glmuons.first; it != range_glmuons.second; it++) { - auto fwdtrack = fwdtracks.rawIteratorAt(it->second); - fillFwdTrackTable(collision, fwdtrack, mftCovs, mapAmb[fwdtrack.globalIndex()]); - } - } - } // end of collision loop + // 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(); + // } + // vec_min_chi2MatchMCHMFT.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) { + // 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) { + // auto bc = collision.template bc_as(); + // initCCDB(bc); + // if (!collision.isSelected()) { + // continue; + // } + // if (collision.swtaliastmp_raw() == 0) { + // continue; + // } - multiMapSAMuonsPerCollision.clear(); - multiMapGLMuonsPerCollision.clear(); - 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); + // auto fwdtrackIdsThisCollision = fwdtrackIndices.sliceBy(fwdtrackIndicesPerCollision, collision.globalIndex()); + // for (const auto& fwdtrackId : fwdtrackIdsThisCollision) { + // 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; + // // } + + // if (!fillFwdTrackTable(collision, fwdtrack, mftCovs, mapAmb[fwdtrack.globalIndex()])) { + // continue; + // } + + // if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { + // multiMapGLMuonsPerCollision.insert(std::make_pair(collision.globalIndex(), fwdtrack.globalIndex())); + // } + // if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { + // multiMapSAMuonsPerCollision.insert(std::make_pair(collision.globalIndex(), fwdtrack.globalIndex())); + // } + + // } // end of fwdtrack loop + // } // end of collision loop + + // for (const auto& collision : collisions) { + // int count_samuons = multiMapSAMuonsPerCollision.count(collision.globalIndex()); + // int count_glmuons = multiMapGLMuonsPerCollision.count(collision.globalIndex()); + // if (fillQAHistograms) { + // fRegistry.fill(HIST("MCHMID/hNmu"), count_samuons); + // fRegistry.fill(HIST("MFTMCHMID/hNmu"), count_glmuons); + // } + // if (count_samuons >= minNmuon) { + // auto range_samuons = multiMapSAMuonsPerCollision.equal_range(collision.globalIndex()); + // for (auto it = range_samuons.first; it != range_samuons.second; it++) { + // auto fwdtrack = fwdtracks.rawIteratorAt(it->second); + // fillFwdTrackTable(collision, fwdtrack, mftCovs, mapAmb[fwdtrack.globalIndex()]); + // } + // } + // if (count_glmuons >= minNmuon) { + // auto range_glmuons = multiMapGLMuonsPerCollision.equal_range(collision.globalIndex()); + // for (auto it = range_glmuons.first; it != range_glmuons.second; it++) { + // auto fwdtrack = fwdtracks.rawIteratorAt(it->second); + // fillFwdTrackTable(collision, fwdtrack, mftCovs, mapAmb[fwdtrack.globalIndex()]); + // } + // } + // } // end of collision loop + + // multiMapSAMuonsPerCollision.clear(); + // multiMapGLMuonsPerCollision.clear(); + // 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& mfttracks, aod::BCsWithTimestamps const&, aod::McParticles const&) { @@ -1223,9 +1168,6 @@ struct skimmerPrimaryMuon { // 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; // } - // if (!isBestMatch(fwdtrack, fwdtracks, mfttracks)) { - // continue; - // } if (!fillFwdTrackTable(collision, fwdtrack, nullptr, false)) { continue; @@ -1314,9 +1256,6 @@ struct skimmerPrimaryMuon { // 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; // } - // if (!isBestMatch(fwdtrack, fwdtracks, mfttracks)) { - // continue; - // } if (!fillFwdTrackTable(collision, fwdtrack, nullptr, mapAmb[fwdtrack.globalIndex()])) { continue; @@ -1364,104 +1303,102 @@ struct skimmerPrimaryMuon { } PROCESS_SWITCH(skimmerPrimaryMuon, processMC_TTCA, "process reconstructed and MC info", false); - 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, aod::McParticles const&) - { - for (const auto& mfttrackConv : mftCovs) { - map_mfttrackcovs[mfttrackConv.matchMFTTrackId()] = mfttrackConv.globalIndex(); - } - vec_min_chi2MatchMCHMFT.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) { - 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) { - auto bc = collision.template bc_as(); - initCCDB(bc); - if (!collision.isSelected()) { - continue; - } - if (!collision.has_mcCollision()) { - continue; - } - - auto fwdtrackIdsThisCollision = fwdtrackIndices.sliceBy(fwdtrackIndicesPerCollision, collision.globalIndex()); - for (const auto& fwdtrackId : fwdtrackIdsThisCollision) { - 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; - // } - // if (!isBestMatch(fwdtrack, fwdtracks, mfttracks)) { - // continue; - // } - - if (!fillFwdTrackTable(collision, fwdtrack, mftCovs, mapAmb[fwdtrack.globalIndex()])) { - continue; - } - - if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { - multiMapGLMuonsPerCollision.insert(std::make_pair(collision.globalIndex(), fwdtrack.globalIndex())); - } - if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { - multiMapSAMuonsPerCollision.insert(std::make_pair(collision.globalIndex(), fwdtrack.globalIndex())); - } - - } // end of fwdtrack loop - } // end of collision loop - - for (const auto& collision : collisions) { - int count_samuons = multiMapSAMuonsPerCollision.count(collision.globalIndex()); - int count_glmuons = multiMapGLMuonsPerCollision.count(collision.globalIndex()); - if (fillQAHistograms) { - fRegistry.fill(HIST("MCHMID/hNmu"), count_samuons); - fRegistry.fill(HIST("MFTMCHMID/hNmu"), count_glmuons); - } - if (count_samuons >= minNmuon) { - auto range_samuons = multiMapSAMuonsPerCollision.equal_range(collision.globalIndex()); - for (auto it = range_samuons.first; it != range_samuons.second; it++) { - auto fwdtrack = fwdtracks.rawIteratorAt(it->second); - fillFwdTrackTable(collision, fwdtrack, mftCovs, mapAmb[fwdtrack.globalIndex()]); - } - } - if (count_glmuons >= minNmuon) { - auto range_glmuons = multiMapGLMuonsPerCollision.equal_range(collision.globalIndex()); - for (auto it = range_glmuons.first; it != range_glmuons.second; it++) { - auto fwdtrack = fwdtracks.rawIteratorAt(it->second); - fillFwdTrackTable(collision, fwdtrack, mftCovs, mapAmb[fwdtrack.globalIndex()]); - } - } - } // end of collision loop + // 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, aod::McParticles const&) + // { + // for (const auto& mfttrackConv : mftCovs) { + // map_mfttrackcovs[mfttrackConv.matchMFTTrackId()] = mfttrackConv.globalIndex(); + // } + // vec_min_chi2MatchMCHMFT.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) { + // 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) { + // auto bc = collision.template bc_as(); + // initCCDB(bc); + // if (!collision.isSelected()) { + // continue; + // } + // if (!collision.has_mcCollision()) { + // continue; + // } - multiMapSAMuonsPerCollision.clear(); - multiMapGLMuonsPerCollision.clear(); - 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); + // auto fwdtrackIdsThisCollision = fwdtrackIndices.sliceBy(fwdtrackIndicesPerCollision, collision.globalIndex()); + // for (const auto& fwdtrackId : fwdtrackIdsThisCollision) { + // 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; + // // } + + // if (!fillFwdTrackTable(collision, fwdtrack, mftCovs, mapAmb[fwdtrack.globalIndex()])) { + // continue; + // } + + // if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::GlobalMuonTrack) { + // multiMapGLMuonsPerCollision.insert(std::make_pair(collision.globalIndex(), fwdtrack.globalIndex())); + // } + // if (fwdtrack.trackType() == o2::aod::fwdtrack::ForwardTrackTypeEnum::MuonStandaloneTrack) { + // multiMapSAMuonsPerCollision.insert(std::make_pair(collision.globalIndex(), fwdtrack.globalIndex())); + // } + + // } // end of fwdtrack loop + // } // end of collision loop + + // for (const auto& collision : collisions) { + // int count_samuons = multiMapSAMuonsPerCollision.count(collision.globalIndex()); + // int count_glmuons = multiMapGLMuonsPerCollision.count(collision.globalIndex()); + // if (fillQAHistograms) { + // fRegistry.fill(HIST("MCHMID/hNmu"), count_samuons); + // fRegistry.fill(HIST("MFTMCHMID/hNmu"), count_glmuons); + // } + // if (count_samuons >= minNmuon) { + // auto range_samuons = multiMapSAMuonsPerCollision.equal_range(collision.globalIndex()); + // for (auto it = range_samuons.first; it != range_samuons.second; it++) { + // auto fwdtrack = fwdtracks.rawIteratorAt(it->second); + // fillFwdTrackTable(collision, fwdtrack, mftCovs, mapAmb[fwdtrack.globalIndex()]); + // } + // } + // if (count_glmuons >= minNmuon) { + // auto range_glmuons = multiMapGLMuonsPerCollision.equal_range(collision.globalIndex()); + // for (auto it = range_glmuons.first; it != range_glmuons.second; it++) { + // auto fwdtrack = fwdtracks.rawIteratorAt(it->second); + // fillFwdTrackTable(collision, fwdtrack, mftCovs, mapAmb[fwdtrack.globalIndex()]); + // } + // } + // } // end of collision loop + + // multiMapSAMuonsPerCollision.clear(); + // multiMapGLMuonsPerCollision.clear(); + // 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); }; + struct associateAmbiguousMuon { Produces em_amb_muon_ids; diff --git a/PWGEM/Dilepton/TableProducer/skimmerPrimaryTrack.cxx b/PWGEM/Dilepton/TableProducer/skimmerPrimaryTrack.cxx index bc88fd7fef0..b7c6b98824f 100644 --- a/PWGEM/Dilepton/TableProducer/skimmerPrimaryTrack.cxx +++ b/PWGEM/Dilepton/TableProducer/skimmerPrimaryTrack.cxx @@ -15,6 +15,7 @@ #include "PWGEM/Dilepton/DataModel/dileptonTables.h" #include "PWGEM/Dilepton/Utils/EMTrackUtilities.h" +#include "Common/Core/RecoDecay.h" #include "Common/Core/TableHelper.h" #include "Common/Core/trackUtilities.h" #include "Common/DataModel/CollisionAssociationTables.h" @@ -260,7 +261,7 @@ struct skimmerPrimaryTrack { float pt = trackParCov.getPt(); float eta = trackParCov.getEta(); float phi = trackParCov.getPhi(); - o2::math_utils::bringTo02Pi(phi); + phi = RecoDecay::constrainAngle(phi, 0, 1U); uint16_t trackBit = 0; // As minimal cuts, following cuts are applied. The cut values are hardcoded on the purpose for consistent bit operation. @@ -273,49 +274,50 @@ struct skimmerPrimaryTrack { // Ncr/Nf ratio in TPC > 0.8 if (track.itsNCls() >= 5) { - trackBit |= static_cast(RefTrackBit::kNclsITS5); + SETBIT(trackBit, static_cast(o2::aod::pwgem::dilepton::utils::emtrackutil::RefTrackBit::kNclsITS5)); } if (track.itsNCls() >= 6) { - trackBit |= static_cast(RefTrackBit::kNclsITS6); + SETBIT(trackBit, static_cast(o2::aod::pwgem::dilepton::utils::emtrackutil::RefTrackBit::kNclsITS6)); } if (track.tpcNClsCrossedRows() >= 70) { - trackBit |= static_cast(RefTrackBit::kNcrTPC70); + SETBIT(trackBit, static_cast(o2::aod::pwgem::dilepton::utils::emtrackutil::RefTrackBit::kNcrTPC70)); } if (track.tpcNClsCrossedRows() >= 90) { - trackBit |= static_cast(RefTrackBit::kNcrTPC90); + SETBIT(trackBit, static_cast(o2::aod::pwgem::dilepton::utils::emtrackutil::RefTrackBit::kNcrTPC90)); } if (track.tpcNClsFound() >= 50) { - trackBit |= static_cast(RefTrackBit::kNclsTPC50); + SETBIT(trackBit, static_cast(o2::aod::pwgem::dilepton::utils::emtrackutil::RefTrackBit::kNclsTPC50)); } if (track.tpcNClsFound() >= 70) { - trackBit |= static_cast(RefTrackBit::kNclsTPC70); + SETBIT(trackBit, static_cast(o2::aod::pwgem::dilepton::utils::emtrackutil::RefTrackBit::kNclsTPC70)); } if (track.tpcNClsFound() >= 90) { - trackBit |= static_cast(RefTrackBit::kNclsTPC90); + SETBIT(trackBit, static_cast(o2::aod::pwgem::dilepton::utils::emtrackutil::RefTrackBit::kNclsTPC90)); } if (track.tpcChi2NCl() < 4.f) { - trackBit |= static_cast(RefTrackBit::kChi2TPC4); + SETBIT(trackBit, static_cast(o2::aod::pwgem::dilepton::utils::emtrackutil::RefTrackBit::kChi2TPC4)); } if (track.tpcChi2NCl() < 3.f) { - trackBit |= static_cast(RefTrackBit::kChi2TPC3); + SETBIT(trackBit, static_cast(o2::aod::pwgem::dilepton::utils::emtrackutil::RefTrackBit::kChi2TPC3)); } + if (track.tpcFractionSharedCls() < 0.7) { - trackBit |= static_cast(RefTrackBit::kFracSharedTPC07); + SETBIT(trackBit, static_cast(o2::aod::pwgem::dilepton::utils::emtrackutil::RefTrackBit::kFracSharedTPC07)); } if (std::fabs(dcaZ) < 0.5) { - trackBit |= static_cast(RefTrackBit::kDCAz05cm); + SETBIT(trackBit, static_cast(o2::aod::pwgem::dilepton::utils::emtrackutil::RefTrackBit::kDCAz05cm)); } if (std::fabs(dcaZ) < 0.3) { - trackBit |= static_cast(RefTrackBit::kDCAz03cm); + SETBIT(trackBit, static_cast(o2::aod::pwgem::dilepton::utils::emtrackutil::RefTrackBit::kDCAz03cm)); } if (std::fabs(dcaXY) < 0.5) { - trackBit |= static_cast(RefTrackBit::kDCAxy05cm); + SETBIT(trackBit, static_cast(o2::aod::pwgem::dilepton::utils::emtrackutil::RefTrackBit::kDCAxy05cm)); } if (std::fabs(dcaXY) < 0.3) { - trackBit |= static_cast(RefTrackBit::kDCAxy03cm); + SETBIT(trackBit, static_cast(o2::aod::pwgem::dilepton::utils::emtrackutil::RefTrackBit::kDCAxy03cm)); } emprimarytracks(/*collision.globalIndex(),*/ /*track.globalIndex(),*/ track.sign() / pt, eta, phi, trackBit); diff --git a/PWGEM/Dilepton/Utils/EMTrackUtilities.h b/PWGEM/Dilepton/Utils/EMTrackUtilities.h index 88daecccb62..47e73e31f96 100644 --- a/PWGEM/Dilepton/Utils/EMTrackUtilities.h +++ b/PWGEM/Dilepton/Utils/EMTrackUtilities.h @@ -29,34 +29,43 @@ namespace o2::aod::pwgem::dilepton::utils::emtrackutil { -enum class RefTrackBit : uint16_t { // This is not for leptons, but charged tracks for reference flow. - kNclsITS5 = 1, - kNclsITS6 = 2, - kNcrTPC70 = 4, - kNcrTPC90 = 8, - kNclsTPC50 = 16, // (not necessary, if ncr is used.) - kNclsTPC70 = 32, // (not necessary, if ncr is used.) - kNclsTPC90 = 64, // (not necessary, if ncr is used.) - kChi2TPC4 = 128, - kChi2TPC3 = 256, - kFracSharedTPC07 = 512, - kDCAxy05cm = 1024, // default is 1 cm - kDCAxy03cm = 2048, - kDCAz05cm = 4096, // default is 1cm - kDCAz03cm = 8192, +enum class RefTrackType : int { // charged tracks for reference flow. + kCB = 0, + kMFTsa = 1, }; -enum class RefMFTTrackBit : uint16_t { // This is not for leptons, but charged tracks for reference flow. - kNclsMFT6 = 1, // default is 5 - kNclsMFT7 = 2, - kNclsMFT8 = 4, - kChi2MFT3 = 8, // default is 4 - kChi2MFT2 = 16, - kDCAxy005cm = 32, // default is 0.06 cm - kDCAxy004cm = 64, - kDCAxy003cm = 128, - kDCAxy002cm = 256, - kDCAxy001cm = 512, +// This is not for leptons, but charged tracks for reference flow. +enum class RefTrackBit : int { + kNclsITS5 = 0, + kNclsITS6, + kNcrTPC70, + kNcrTPC90, + kNclsTPC50, // (not necessary, if ncr is used.) + kNclsTPC70, // (not necessary, if ncr is used.) + kNclsTPC90, // (not necessary, if ncr is used.) + kChi2TPC4, + kChi2TPC3, + kFracSharedTPC07, + kDCAxy05cm, // default is 1 cm + kDCAxy03cm, + kDCAz05cm, // default is 1cm + kDCAz03cm, + kNCuts, +}; + +// This is not for leptons, but charged tracks for reference flow. +enum class RefMFTTrackBit : int { + kNclsMFT6 = 0, // default is 5 + kNclsMFT7, + kNclsMFT8, + kChi2MFT3, // default is 4 + kChi2MFT2, + kDCAxy005cm, // default is 0.06 cm + kDCAxy004cm, + kDCAxy003cm, + kDCAxy002cm, + kDCAxy001cm, + kNCuts, }; //_______________________________________________________________________