From df6849ad8d9830682420517605f0bfb27cae5753 Mon Sep 17 00:00:00 2001 From: Daiki Sekihata Date: Thu, 26 Dec 2024 23:57:22 +0900 Subject: [PATCH] [PWGEM/Dilepton] update 3D smearing (#9137) --- .../TableProducer/skimmerPrimaryElectron.cxx | 303 +++++++++++++++--- PWGEM/Dilepton/Tasks/smearing.cxx | 18 +- PWGEM/Dilepton/Utils/MomentumSmearer.h | 203 ++++++++++-- 3 files changed, 449 insertions(+), 75 deletions(-) diff --git a/PWGEM/Dilepton/TableProducer/skimmerPrimaryElectron.cxx b/PWGEM/Dilepton/TableProducer/skimmerPrimaryElectron.cxx index 18a0cb8c988..83efda8c7c3 100644 --- a/PWGEM/Dilepton/TableProducer/skimmerPrimaryElectron.cxx +++ b/PWGEM/Dilepton/TableProducer/skimmerPrimaryElectron.cxx @@ -589,7 +589,10 @@ struct skimmerPrimaryElectron { struct prefilterPrimaryElectron { enum class EM_Electron_PF : int { - kElFromPC = 0, // electron from photon conversion + kElFromPC = 0, // electron from photon conversion + kElFromPi0_1 = 1, // electron from pi0 dalitz decay, threshold 1 + kElFromPi0_2 = 2, // electron from pi0 dalitz decay, threshold 2 + kElFromPi0_3 = 3, // electron from pi0 dalitz decay, threshold 3 }; Produces ele_pfb; @@ -610,7 +613,7 @@ struct prefilterPrimaryElectron { Configurable max_dcaxy{"max_dcaxy", 0.3, "DCAxy To PV for loose track sample"}; Configurable max_dcaz{"max_dcaz", 0.3, "DCAz To PV for loose track sample"}; Configurable minpt{"minpt", 0.1, "min pt for track for loose track sample"}; - Configurable maxeta{"maxeta", 1.2, "eta acceptance for loose track sample"}; + Configurable maxeta{"maxeta", 0.9, "eta acceptance for loose track sample"}; Configurable min_ncluster_tpc{"min_ncluster_tpc", 0, "min ncluster tpc"}; Configurable mincrossedrows{"mincrossedrows", 70, "min crossed rows"}; Configurable max_frac_shared_clusters_tpc{"max_frac_shared_clusters_tpc", 999.f, "max fraction of shared clusters in TPC"}; @@ -624,12 +627,9 @@ struct prefilterPrimaryElectron { Configurable slope{"slope", 0.0185, "slope for m vs. phiv"}; Configurable intercept{"intercept", -0.0280, "intercept for m vs. phiv"}; - HistogramRegistry fRegistry{ - "fRegistry", - { - {"hMvsPhiV_PV", "mass vs. phiv;#varphi_{V} (rad.);m_{ee} (GeV/c^{2})", {HistType::kTH2F, {{90, .0f, TMath::Pi()}, {100, 0, 0.1}}}}, - }, - }; + Configurable> max_mee_vec{"max_mee_vec", std::vector{0.08, 0.10, 0.12}, "vector fo max mee for prefilter in ULS. Please sort this by increasing order."}; // currently, 3 thoresholds are allowed. + + HistogramRegistry fRegistry{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; int mRunNumber; float d_bz; @@ -645,6 +645,22 @@ struct prefilterPrimaryElectron { ccdb->setCaching(true); ccdb->setLocalObjectValidityChecking(); ccdb->setFatalWhenNull(false); + + if (!doprocessDummy) { + addHistograms(); + } + } + + void addHistograms() + { + fRegistry.add("Track/hPt", "pT;p_{T} (GeV/c)", kTH1F, {{1000, 0.0f, 10}}, false); + fRegistry.add("Track/hEtaPhi", "#eta vs. #varphi;#varphi (rad.);#eta", kTH2F, {{180, 0, 2 * M_PI}, {40, -1.0f, 1.0f}}, false); + fRegistry.add("Track/hTPCNsigmaEl", "loose track TPC PID", kTH2F, {{1000, 0.f, 10}, {100, -5, +5}}); + fRegistry.add("Pair/before/uls/hMvsPt", "mass vs. pT;m_{ee} (GeV/c^{2});p_{T,ee} (GeV/c)", kTH2F, {{400, 0, 4}, {100, 0, 10}}); + fRegistry.add("Pair/before/uls/hMvsPhiV", "mass vs. phiv;#varphi_{V} (rad.);m_{ee} (GeV/c^{2})", kTH2F, {{90, 0.f, M_PI}, {100, 0, 1.f}}); + fRegistry.addClone("Pair/before/uls/", "Pair/before/lspp/"); + fRegistry.addClone("Pair/before/uls/", "Pair/before/lsmm/"); + fRegistry.addClone("Pair/before/", "Pair/after/"); } void initCCDB(aod::BCsWithTimestamps::iterator const& bc) @@ -748,11 +764,7 @@ struct prefilterPrimaryElectron { template bool reconstructPC(TCollision const& collision, TTrack1 const& ele, TTrack2 const& pos) { - if (pos.sign() * ele.sign() > 0) { // reject same sign pair - return false; - } - - float mee, phiv = 0; + float mee = 0, phiv = 0; gpu::gpustd::array dcaInfo; std::array pVec_recalc = {0, 0, 0}; // px, py, pz @@ -778,8 +790,6 @@ struct prefilterPrimaryElectron { phiv = o2::aod::pwgem::dilepton::utils::pairutil::getPhivPair(pos.px(), pos.py(), pos.pz(), pVec_recalc[0], pVec_recalc[1], pVec_recalc[2], pos.sign(), ele.sign(), d_bz); } - fRegistry.fill(HIST("hMvsPhiV_PV"), phiv, mee); - if (mee < slope * phiv + intercept) { return true; } else { @@ -822,7 +832,8 @@ struct prefilterPrimaryElectron { if (!checkTrack(collision, track)) { continue; } - + fRegistry.fill(HIST("Track/hPt"), track.pt()); + fRegistry.fill(HIST("Track/hEtaPhi"), track.phi(), track.eta()); if (track.sign() > 0) { posTracks_per_coll.emplace_back(track); } else { @@ -830,37 +841,125 @@ struct prefilterPrimaryElectron { } } - for (auto& empos : positrons_per_coll) { - // auto pos = tracks.rawIteratorAt(empos.trackId()); // use rawIterator, if the table is filtered. - for (auto& ele : negTracks_per_coll) { - if (!checkTrack(collision, ele)) { - continue; - } + for (auto& ele : negTracks_per_coll) { + if (!checkTrack(collision, ele)) { + continue; + } + gpu::gpustd::array dcaInfo; + std::array pVec_recalc = {0, 0, 0}; // px, py, pz + auto track_par_cov_recalc = getTrackParCov(ele); + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, track_par_cov_recalc, 2.f, matCorr, &dcaInfo); + getPxPyPz(track_par_cov_recalc, pVec_recalc); + + for (auto& empos : positrons_per_coll) { if (empos.trackId() == ele.globalIndex()) { continue; } - bool isPC = reconstructPC<-1>(collision, ele, empos); - if (isPC) { - pfb_map[empos.globalIndex()] |= (uint8_t(1) << static_cast(EM_Electron_PF::kElFromPC)); + + ROOT::Math::PtEtaPhiMVector v1(track_par_cov_recalc.getPt(), track_par_cov_recalc.getEta(), track_par_cov_recalc.getPhi(), o2::constants::physics::MassElectron); // loose track + ROOT::Math::PtEtaPhiMVector v2(empos.pt(), empos.eta(), empos.phi(), o2::constants::physics::MassElectron); // signal track + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + float phiv = o2::aod::pwgem::dilepton::utils::pairutil::getPhivPair(empos.px(), empos.py(), empos.pz(), pVec_recalc[0], pVec_recalc[1], pVec_recalc[2], empos.sign(), ele.sign(), d_bz); + fRegistry.fill(HIST("Pair/before/uls/hMvsPhiV"), phiv, v12.M()); + fRegistry.fill(HIST("Pair/before/uls/hMvsPt"), v12.M(), v12.Pt()); + if (v12.M() < max_mee_vec->at(static_cast(max_mee_vec->size()) - 1)) { + fRegistry.fill(HIST("Track/hTPCNsigmaEl"), ele.tpcInnerParam(), ele.tpcNSigmaEl()); + } + for (int i = 0; i < static_cast(max_mee_vec->size()); i++) { + if (v12.M() < max_mee_vec->at(i)) { + pfb_map[empos.globalIndex()] |= (uint8_t(1) << (static_cast(EM_Electron_PF::kElFromPi0_1) + i)); + } } - } // end of loose electron loop - } // end of signal positon loop - for (auto& emele : electrons_per_coll) { - // auto ele = tracks.rawIteratorAt(emele.trackId()); // use rawIterator, if the table is filtered. - for (auto& pos : posTracks_per_coll) { - if (!checkTrack(collision, pos)) { // track cut is applied to loose sample - continue; + if (v12.M() < slope * phiv + intercept) { + pfb_map[empos.globalIndex()] |= (uint8_t(1) << static_cast(EM_Electron_PF::kElFromPC)); } + + } // end of signal positon loop + } // end of loose electron loop + + for (auto& pos : posTracks_per_coll) { + if (!checkTrack(collision, pos)) { // track cut is applied to loose sample + continue; + } + gpu::gpustd::array dcaInfo; + std::array pVec_recalc = {0, 0, 0}; // px, py, pz + auto track_par_cov_recalc = getTrackParCov(pos); + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, track_par_cov_recalc, 2.f, matCorr, &dcaInfo); + getPxPyPz(track_par_cov_recalc, pVec_recalc); + for (auto& emele : electrons_per_coll) { if (emele.trackId() == pos.globalIndex()) { continue; } - bool isPC = reconstructPC<+1>(collision, emele, pos); - if (isPC) { + + ROOT::Math::PtEtaPhiMVector v1(emele.pt(), emele.eta(), emele.phi(), o2::constants::physics::MassElectron); // signal track + ROOT::Math::PtEtaPhiMVector v2(track_par_cov_recalc.getPt(), track_par_cov_recalc.getEta(), track_par_cov_recalc.getPhi(), o2::constants::physics::MassElectron); // loose track + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + float phiv = o2::aod::pwgem::dilepton::utils::pairutil::getPhivPair(pVec_recalc[0], pVec_recalc[1], pVec_recalc[2], emele.px(), emele.py(), emele.pz(), pos.sign(), emele.sign(), d_bz); + fRegistry.fill(HIST("Pair/before/uls/hMvsPhiV"), phiv, v12.M()); + fRegistry.fill(HIST("Pair/before/uls/hMvsPt"), v12.M(), v12.Pt()); + if (v12.M() < max_mee_vec->at(static_cast(max_mee_vec->size()) - 1)) { + fRegistry.fill(HIST("Track/hTPCNsigmaEl"), pos.tpcInnerParam(), pos.tpcNSigmaEl()); + } + for (int i = 0; i < static_cast(max_mee_vec->size()); i++) { + if (v12.M() < max_mee_vec->at(i)) { + pfb_map[emele.globalIndex()] |= (uint8_t(1) << (static_cast(EM_Electron_PF::kElFromPi0_1) + i)); + } + } + + if (v12.M() < slope * phiv + intercept) { pfb_map[emele.globalIndex()] |= (uint8_t(1) << static_cast(EM_Electron_PF::kElFromPC)); } - } // end of loose positon loop - } // end of signal electron loop + } // end of signal electron loop + } // end of loose positon loop + + for (auto& pos : posTracks_per_coll) { + if (!checkTrack(collision, pos)) { // track cut is applied to loose sample + continue; + } + gpu::gpustd::array dcaInfo; + std::array pVec_recalc = {0, 0, 0}; // px, py, pz + auto track_par_cov_recalc = getTrackParCov(pos); + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, track_par_cov_recalc, 2.f, matCorr, &dcaInfo); + getPxPyPz(track_par_cov_recalc, pVec_recalc); + for (auto& empos : positrons_per_coll) { + if (empos.trackId() == pos.globalIndex()) { + continue; + } + + ROOT::Math::PtEtaPhiMVector v1(empos.pt(), empos.eta(), empos.phi(), o2::constants::physics::MassElectron); // signal track + ROOT::Math::PtEtaPhiMVector v2(track_par_cov_recalc.getPt(), track_par_cov_recalc.getEta(), track_par_cov_recalc.getPhi(), o2::constants::physics::MassElectron); // loose track + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + float phiv = o2::aod::pwgem::dilepton::utils::pairutil::getPhivPair(pVec_recalc[0], pVec_recalc[1], pVec_recalc[2], empos.px(), empos.py(), empos.pz(), pos.sign(), empos.sign(), d_bz); + fRegistry.fill(HIST("Pair/before/lspp/hMvsPhiV"), phiv, v12.M()); + fRegistry.fill(HIST("Pair/before/lspp/hMvsPt"), v12.M(), v12.Pt()); + } // end of signal positron loop + } // end of loose positon loop + + for (auto& ele : negTracks_per_coll) { + if (!checkTrack(collision, ele)) { + continue; + } + gpu::gpustd::array dcaInfo; + std::array pVec_recalc = {0, 0, 0}; // px, py, pz + auto track_par_cov_recalc = getTrackParCov(ele); + o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, track_par_cov_recalc, 2.f, matCorr, &dcaInfo); + getPxPyPz(track_par_cov_recalc, pVec_recalc); + + for (auto& emele : electrons_per_coll) { + if (emele.trackId() == ele.globalIndex()) { + continue; + } + + ROOT::Math::PtEtaPhiMVector v1(track_par_cov_recalc.getPt(), track_par_cov_recalc.getEta(), track_par_cov_recalc.getPhi(), o2::constants::physics::MassElectron); // loose track + ROOT::Math::PtEtaPhiMVector v2(emele.pt(), emele.eta(), emele.phi(), o2::constants::physics::MassElectron); // signal track + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + float phiv = o2::aod::pwgem::dilepton::utils::pairutil::getPhivPair(emele.px(), emele.py(), emele.pz(), pVec_recalc[0], pVec_recalc[1], pVec_recalc[2], emele.sign(), ele.sign(), d_bz); + fRegistry.fill(HIST("Pair/before/lsmm/hMvsPhiV"), phiv, v12.M()); + fRegistry.fill(HIST("Pair/before/lsmm/hMvsPt"), v12.M(), v12.Pt()); + + } // end of signal electron loop + } // end of loose electron loop posTracks_per_coll.clear(); negTracks_per_coll.clear(); @@ -872,6 +971,26 @@ struct prefilterPrimaryElectron { ele_pfb(pfb_map[ele.globalIndex()]); } + // check prefilter + for (auto& collision : collisions) { + auto positrons_per_coll = positrons->sliceByCachedUnsorted(o2::aod::emprimaryelectron::collisionId, collision.globalIndex(), cache); // signal sample + auto electrons_per_coll = electrons->sliceByCachedUnsorted(o2::aod::emprimaryelectron::collisionId, collision.globalIndex(), cache); // signal sample + + for (auto& [ele, pos] : combinations(CombinationsFullIndexPolicy(electrons_per_coll, positrons_per_coll))) { + if (pfb_map[ele.globalIndex()] != 0 || pfb_map[pos.globalIndex()] != 0) { + continue; + } + + ROOT::Math::PtEtaPhiMVector v1(ele.pt(), ele.eta(), ele.phi(), o2::constants::physics::MassElectron); + ROOT::Math::PtEtaPhiMVector v2(pos.pt(), pos.eta(), pos.phi(), o2::constants::physics::MassElectron); + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + float phiv = o2::aod::pwgem::dilepton::utils::pairutil::getPhivPair(pos.px(), pos.py(), pos.pz(), ele.px(), ele.py(), ele.pz(), pos.sign(), ele.sign(), d_bz); + fRegistry.fill(HIST("Pair/after/uls/hMvsPhiV"), phiv, v12.M()); + fRegistry.fill(HIST("Pair/after/uls/hMvsPt"), v12.M(), v12.Pt()); + + } // end of ULS pairing + } // end of collision loop + pfb_map.clear(); } PROCESS_SWITCH(prefilterPrimaryElectron, processPrefilter_TTCA, "process prefilter with TTCA", false); @@ -893,6 +1012,21 @@ struct prefilterPrimaryElectron { auto positrons_per_coll = positrons->sliceByCachedUnsorted(o2::aod::emprimaryelectron::collisionId, collision.globalIndex(), cache); // signal sample auto electrons_per_coll = electrons->sliceByCachedUnsorted(o2::aod::emprimaryelectron::collisionId, collision.globalIndex(), cache); // signal sample + for (auto& pos : posTracks_per_coll) { + if (!checkTrack(collision, pos)) { // track cut is applied to loose sample + continue; + } + fRegistry.fill(HIST("Track/hPt"), pos.pt()); + fRegistry.fill(HIST("Track/hEtaPhi"), pos.phi(), pos.eta()); + } + for (auto& neg : negTracks_per_coll) { + if (!checkTrack(collision, neg)) { // track cut is applied to loose sample + continue; + } + fRegistry.fill(HIST("Track/hPt"), neg.pt()); + fRegistry.fill(HIST("Track/hEtaPhi"), neg.phi(), neg.eta()); + } + for (auto& [ele, empos] : combinations(CombinationsFullIndexPolicy(negTracks_per_coll, positrons_per_coll))) { // auto pos = tracks.rawIteratorAt(empos.trackId()); // use rawIterator, if the table is filtered. if (!checkTrack(collision, ele)) { // track cut is applied to loose sample @@ -902,11 +1036,26 @@ struct prefilterPrimaryElectron { continue; } - bool isPC = reconstructPC<-1>(collision, ele, empos); - if (isPC) { + ROOT::Math::PtEtaPhiMVector v1(ele.pt(), ele.eta(), ele.phi(), o2::constants::physics::MassElectron); // loose track + ROOT::Math::PtEtaPhiMVector v2(empos.pt(), empos.eta(), empos.phi(), o2::constants::physics::MassElectron); // signal track + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + float phiv = o2::aod::pwgem::dilepton::utils::pairutil::getPhivPair(empos.px(), empos.py(), empos.pz(), ele.px(), ele.py(), ele.pz(), empos.sign(), ele.sign(), d_bz); + fRegistry.fill(HIST("Pair/before/uls/hMvsPhiV"), phiv, v12.M()); + fRegistry.fill(HIST("Pair/before/uls/hMvsPt"), v12.M(), v12.Pt()); + if (v12.M() < max_mee_vec->at(static_cast(max_mee_vec->size()) - 1)) { + fRegistry.fill(HIST("Track/hTPCNsigmaEl"), ele.tpcInnerParam(), ele.tpcNSigmaEl()); + } + for (int i = 0; i < static_cast(max_mee_vec->size()); i++) { + if (v12.M() < max_mee_vec->at(i)) { + pfb_map[empos.globalIndex()] |= (uint8_t(1) << (static_cast(EM_Electron_PF::kElFromPi0_1) + i)); + } + } + + if (v12.M() < slope * phiv + intercept) { pfb_map[empos.globalIndex()] |= (uint8_t(1) << static_cast(EM_Electron_PF::kElFromPC)); } - } + + } // end of ULS pairing for (auto& [pos, emele] : combinations(CombinationsFullIndexPolicy(posTracks_per_coll, electrons_per_coll))) { // auto ele = tracks.rawIteratorAt(emele.trackId()); // use rawIterator, if the table is filtered. @@ -916,17 +1065,87 @@ struct prefilterPrimaryElectron { if (emele.trackId() == pos.globalIndex()) { continue; } - bool isPC = reconstructPC<+1>(collision, emele, pos); - if (isPC) { + + ROOT::Math::PtEtaPhiMVector v1(emele.pt(), emele.eta(), emele.phi(), o2::constants::physics::MassElectron); // signal track + ROOT::Math::PtEtaPhiMVector v2(pos.pt(), pos.eta(), pos.phi(), o2::constants::physics::MassElectron); // loose track + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + float phiv = o2::aod::pwgem::dilepton::utils::pairutil::getPhivPair(pos.px(), pos.py(), pos.pz(), emele.px(), emele.py(), emele.pz(), pos.sign(), emele.sign(), d_bz); + fRegistry.fill(HIST("Pair/before/uls/hMvsPhiV"), phiv, v12.M()); + fRegistry.fill(HIST("Pair/before/uls/hMvsPt"), v12.M(), v12.Pt()); + if (v12.M() < max_mee_vec->at(static_cast(max_mee_vec->size()) - 1)) { + fRegistry.fill(HIST("Track/hTPCNsigmaEl"), pos.tpcInnerParam(), pos.tpcNSigmaEl()); + } + for (int i = 0; i < static_cast(max_mee_vec->size()); i++) { + if (v12.M() < max_mee_vec->at(i)) { + pfb_map[emele.globalIndex()] |= (uint8_t(1) << (static_cast(EM_Electron_PF::kElFromPi0_1) + i)); + } + } + + if (v12.M() < slope * phiv + intercept) { pfb_map[emele.globalIndex()] |= (uint8_t(1) << static_cast(EM_Electron_PF::kElFromPC)); } - } + + } // end of ULS pairing + + for (auto& [pos, empos] : combinations(CombinationsFullIndexPolicy(posTracks_per_coll, positrons_per_coll))) { + // auto pos = tracks.rawIteratorAt(empos.trackId()); // use rawIterator, if the table is filtered. + if (!checkTrack(collision, pos)) { // track cut is applied to loose sample + continue; + } + if (empos.trackId() == pos.globalIndex()) { + continue; + } + + ROOT::Math::PtEtaPhiMVector v1(pos.pt(), pos.eta(), pos.phi(), o2::constants::physics::MassElectron); // loose track + ROOT::Math::PtEtaPhiMVector v2(empos.pt(), empos.eta(), empos.phi(), o2::constants::physics::MassElectron); // signal track + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + float phiv = o2::aod::pwgem::dilepton::utils::pairutil::getPhivPair(empos.px(), empos.py(), empos.pz(), pos.px(), pos.py(), pos.pz(), empos.sign(), pos.sign(), d_bz); + fRegistry.fill(HIST("Pair/before/lspp/hMvsPhiV"), phiv, v12.M()); + fRegistry.fill(HIST("Pair/before/lspp/hMvsPt"), v12.M(), v12.Pt()); + } // end of LS++ pairing + + for (auto& [ele, emele] : combinations(CombinationsFullIndexPolicy(negTracks_per_coll, electrons_per_coll))) { + // auto ele = tracks.rawIteratorAt(emele.trackId()); // use rawIterator, if the table is filtered. + if (!checkTrack(collision, ele)) { // track cut is applied to loose sample + continue; + } + if (emele.trackId() == ele.globalIndex()) { + continue; + } + + ROOT::Math::PtEtaPhiMVector v1(ele.pt(), ele.eta(), ele.phi(), o2::constants::physics::MassElectron); // loose track + ROOT::Math::PtEtaPhiMVector v2(emele.pt(), emele.eta(), emele.phi(), o2::constants::physics::MassElectron); // signal track + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + float phiv = o2::aod::pwgem::dilepton::utils::pairutil::getPhivPair(emele.px(), emele.py(), emele.pz(), ele.px(), ele.py(), ele.pz(), emele.sign(), ele.sign(), d_bz); + fRegistry.fill(HIST("Pair/before/lsmm/hMvsPhiV"), phiv, v12.M()); + fRegistry.fill(HIST("Pair/before/lsmm/hMvsPt"), v12.M(), v12.Pt()); + } // end of LS-- pairing + } // end of collision loop for (auto& ele : primaryelectrons) { ele_pfb(pfb_map[ele.globalIndex()]); } + // check prefilter + for (auto& collision : collisions) { + auto positrons_per_coll = positrons->sliceByCachedUnsorted(o2::aod::emprimaryelectron::collisionId, collision.globalIndex(), cache); // signal sample + auto electrons_per_coll = electrons->sliceByCachedUnsorted(o2::aod::emprimaryelectron::collisionId, collision.globalIndex(), cache); // signal sample + + for (auto& [ele, pos] : combinations(CombinationsFullIndexPolicy(electrons_per_coll, positrons_per_coll))) { + if (pfb_map[ele.globalIndex()] != 0 || pfb_map[pos.globalIndex()] != 0) { + continue; + } + + ROOT::Math::PtEtaPhiMVector v1(ele.pt(), ele.eta(), ele.phi(), o2::constants::physics::MassElectron); + ROOT::Math::PtEtaPhiMVector v2(pos.pt(), pos.eta(), pos.phi(), o2::constants::physics::MassElectron); + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + float phiv = o2::aod::pwgem::dilepton::utils::pairutil::getPhivPair(pos.px(), pos.py(), pos.pz(), ele.px(), ele.py(), ele.pz(), pos.sign(), ele.sign(), d_bz); + fRegistry.fill(HIST("Pair/after/uls/hMvsPhiV"), phiv, v12.M()); + fRegistry.fill(HIST("Pair/after/uls/hMvsPt"), v12.M(), v12.Pt()); + } // end of ULS pairing + } // end of collision loop + pfb_map.clear(); } PROCESS_SWITCH(prefilterPrimaryElectron, processPrefilter_SA, "process prefilter standalone", false); diff --git a/PWGEM/Dilepton/Tasks/smearing.cxx b/PWGEM/Dilepton/Tasks/smearing.cxx index 82944dc48a0..6fe3399d1d4 100644 --- a/PWGEM/Dilepton/Tasks/smearing.cxx +++ b/PWGEM/Dilepton/Tasks/smearing.cxx @@ -43,7 +43,9 @@ struct ApplySmearing { struct : ConfigurableGroup { std::string prefix = "electron_filename_group"; + Configurable fConfigNDSmearing{"cfgNDSmearing", false, "apply ND-correlated smearing"}; Configurable fConfigResFileName{"cfgResFileName", "", "name of resolution file"}; + Configurable fConfigResNDHistName{"cfgResNDHistName", "hs_reso", "name of ND resolution file"}; Configurable fConfigResPtHistName{"cfgResPtHistName", "RelPtResArrCocktail", "histogram name for pt in resolution file"}; Configurable fConfigResEtaHistName{"cfgResEtaHistName", "EtaResArr", "histogram name for eta in resolution file"}; Configurable fConfigResPhiPosHistName{"cfgResPhiPosHistName", "PhiPosResArr", "histogram name for phi pos in resolution file"}; @@ -59,7 +61,9 @@ struct ApplySmearing { struct : ConfigurableGroup { std::string prefix = "sa_muon_filename_group"; + Configurable fConfigNDSmearing{"cfgNDSmearing", false, "apply ND-correlated smearing"}; Configurable fConfigResFileName{"cfgResFileName", "", "name of resolution file"}; + Configurable fConfigResNDHistName{"cfgResNDHistName", "hs_reso", "name of ND resolution file"}; Configurable fConfigResPtHistName{"cfgResPtHistName", "RelPtResArrCocktail", "histogram name for pt in resolution file"}; Configurable fConfigResEtaHistName{"cfgResEtaHistName", "EtaResArr", "histogram name for eta in resolution file"}; Configurable fConfigResPhiPosHistName{"cfgResPhiPosHistName", "PhiPosResArr", "histogram name for phi pos in resolution file"}; @@ -75,7 +79,9 @@ struct ApplySmearing { struct : ConfigurableGroup { std::string prefix = "gl_muon_filename_group"; + Configurable fConfigNDSmearing{"cfgNDSmearing", false, "apply ND-correlated smearing"}; Configurable fConfigResFileName{"cfgResFileName", "", "name of resolution file"}; + Configurable fConfigResNDHistName{"cfgResNDHistName", "hs_reso", "name of ND resolution file"}; Configurable fConfigResPtHistName{"cfgResPtHistName", "RelPtResArrCocktail", "histogram name for pt in resolution file"}; Configurable fConfigResEtaHistName{"cfgResEtaHistName", "EtaResArr", "histogram name for eta in resolution file"}; Configurable fConfigResPhiPosHistName{"cfgResPhiPosHistName", "PhiPosResArr", "histogram name for phi pos in resolution file"}; @@ -96,7 +102,9 @@ struct ApplySmearing { void init(InitContext&) { + smearer_Electron.setNDSmearing(electron_filenames.fConfigNDSmearing.value); smearer_Electron.setResFileName(TString(electron_filenames.fConfigResFileName)); + smearer_Electron.setResNDHistName(TString(electron_filenames.fConfigResNDHistName)); smearer_Electron.setResPtHistName(TString(electron_filenames.fConfigResPtHistName)); smearer_Electron.setResEtaHistName(TString(electron_filenames.fConfigResEtaHistName)); smearer_Electron.setResPhiPosHistName(TString(electron_filenames.fConfigResPhiPosHistName)); @@ -106,7 +114,9 @@ struct ApplySmearing { smearer_Electron.setDCAFileName(TString(electron_filenames.fConfigDCAFileName)); smearer_Electron.setDCAHistName(TString(electron_filenames.fConfigDCAHistName)); + smearer_StandaloneMuon.setNDSmearing(sa_muon_filenames.fConfigNDSmearing.value); smearer_StandaloneMuon.setResFileName(TString(sa_muon_filenames.fConfigResFileName)); + smearer_StandaloneMuon.setResNDHistName(TString(sa_muon_filenames.fConfigResNDHistName)); smearer_StandaloneMuon.setResPtHistName(TString(sa_muon_filenames.fConfigResPtHistName)); smearer_StandaloneMuon.setResEtaHistName(TString(sa_muon_filenames.fConfigResEtaHistName)); smearer_StandaloneMuon.setResPhiPosHistName(TString(sa_muon_filenames.fConfigResPhiPosHistName)); @@ -116,7 +126,9 @@ struct ApplySmearing { smearer_StandaloneMuon.setDCAFileName(TString(sa_muon_filenames.fConfigDCAFileName)); smearer_StandaloneMuon.setDCAHistName(TString(sa_muon_filenames.fConfigDCAHistName)); + smearer_GlobalMuon.setNDSmearing(gl_muon_filenames.fConfigNDSmearing.value); smearer_GlobalMuon.setResFileName(TString(gl_muon_filenames.fConfigResFileName)); + smearer_GlobalMuon.setResNDHistName(TString(gl_muon_filenames.fConfigResNDHistName)); smearer_GlobalMuon.setResPtHistName(TString(gl_muon_filenames.fConfigResPtHistName)); smearer_GlobalMuon.setResEtaHistName(TString(gl_muon_filenames.fConfigResEtaHistName)); smearer_GlobalMuon.setResPhiPosHistName(TString(gl_muon_filenames.fConfigResPhiPosHistName)); @@ -239,12 +251,14 @@ struct ApplySmearing { } } - void processDummyMCanalysis(ReducedMCTracks const&) {} + void processDummyMCanalysisEM(aod::EMMCParticles const&) {} + void processDummyMCanalysisDQ(ReducedMCTracks const&) {} PROCESS_SWITCH(ApplySmearing, processMCanalysisEM, "Run for MC analysis which uses skimmed EM data format", false); PROCESS_SWITCH(ApplySmearing, processMCanalysisDQ, "Run for MC analysis which uses skimmed DQ data format", false); PROCESS_SWITCH(ApplySmearing, processCocktail, "Run for cocktail analysis", false); - PROCESS_SWITCH(ApplySmearing, processDummyMCanalysis, "Dummy process function", false); + PROCESS_SWITCH(ApplySmearing, processDummyMCanalysisEM, "Dummy process function", false); + PROCESS_SWITCH(ApplySmearing, processDummyMCanalysisDQ, "Dummy process function", false); PROCESS_SWITCH(ApplySmearing, processDummyCocktail, "Dummy process function", true); }; diff --git a/PWGEM/Dilepton/Utils/MomentumSmearer.h b/PWGEM/Dilepton/Utils/MomentumSmearer.h index 002e7d6d247..d49a8c074b6 100644 --- a/PWGEM/Dilepton/Utils/MomentumSmearer.h +++ b/PWGEM/Dilepton/Utils/MomentumSmearer.h @@ -15,13 +15,18 @@ #ifndef PWGEM_DILEPTON_UTILS_MOMENTUMSMEARER_H_ #define PWGEM_DILEPTON_UTILS_MOMENTUMSMEARER_H_ -#include +#include + +#include +#include +#include +#include #include #include #include #include -#include -#include + +#include "CCDB/BasicCCDBManager.h" #include "Framework/Logger.h" using namespace o2::framework; @@ -48,6 +53,23 @@ class MomentumSmearer init(); } + /// Constructor with resolution ND sparse histogram + MomentumSmearer(TString resFileName, TString resNDHistName) + { + setResFileName(resFileName); + setResNDHistName(resNDHistName); + setResPtHistName(""); + setResEtaHistName(""); + setResPhiPosHistName(""); + setResPhiNegHistName(""); + setEffFileName(""); + setEffHistName(""); + setDCAFileName(""); + setDCAHistName(""); + fDoNDSmearing = true; + init(); + } + /// Constructor with resolution histograms and efficiency MomentumSmearer(TString resFileName, TString resPtHistName, TString resEtaHistName, TString resPhiPosHistName, TString resPhiNegHistName, TString effFileName, TString effHistName) { @@ -115,10 +137,46 @@ class MomentumSmearer } } + void fillVecResoND(THnSparseF* hs_reso) + { + LOGP(info, "prepare TH3D"); + const int npt = hs_reso->GetAxis(0)->GetNbins(); + const int neta = hs_reso->GetAxis(1)->GetNbins(); + const int nphi = hs_reso->GetAxis(2)->GetNbins(); + const int nch = hs_reso->GetAxis(3)->GetNbins(); + LOGF(info, "npt = %d, neta = %d, nphi = %d, nch = %d without under- and overflow bins", npt, neta, nphi, nch); + // fVecResoND.reserve(npt * neta * nphi * nch); + + fVecResoND.resize(npt, std::vector>>(neta, std::vector>(nphi, std::vector(nch)))); + // auto h3 = reinterpret_cast(hs_reso->Projection(4, 5, 6)); + // h3->SetName(Form("h3reso_pt%d_eta%d_phi%d_ch%d", 0, 0, 0, 0)); + // fVecResoND[0][0][0][0] = h3; + + for (int ipt = 0; ipt < npt; ipt++) { + hs_reso->GetAxis(0)->SetRange(ipt + 1, ipt + 1); + for (int ieta = 0; ieta < neta; ieta++) { + hs_reso->GetAxis(1)->SetRange(ieta + 1, ieta + 1); + for (int iphi = 0; iphi < nphi; iphi++) { + hs_reso->GetAxis(2)->SetRange(iphi + 1, iphi + 1); + for (int ich = 0; ich < nch; ich++) { + if (-0.5 < hs_reso->GetAxis(3)->GetBinCenter(ich + 1) && hs_reso->GetAxis(3)->GetBinCenter(ich + 1) < 0.5) { + continue; + } + hs_reso->GetAxis(3)->SetRange(ich + 1, ich + 1); + auto h3 = reinterpret_cast(hs_reso->Projection(4, 5, 6)); + h3->SetName(Form("h3reso_pt%d_eta%d_phi%d_ch%d", ipt, ieta, iphi, ich)); + fVecResoND[ipt][ieta][iphi][ich] = h3; + } // end of charge loop + } // end of phi loop + } // end of eta loop + } // end of pt loop + } + void init() { - if (fInitialized) + if (fInitialized) { return; + } if ((fResFileName.BeginsWith("alien://") || fEffFileName.BeginsWith("alien://") || fDCAFileName.BeginsWith("alien://")) && (!fFromCcdb)) { TGrid::Connect("alien://"); @@ -156,34 +214,48 @@ class MomentumSmearer fFile->Close(); } } - if (fResType != 0) { - fResoPt = reinterpret_cast(listRes->FindObject(fResPtHistName)); - if (!fResoPt) { - LOGP(fatal, "Could not open {} from file {}", fResPtHistName.Data(), fResFileName.Data()); - } - fResoEta = reinterpret_cast(listRes->FindObject(fResEtaHistName)); - if (!fResoEta) { - LOGP(fatal, "Could not open {} from file {}", fResEtaHistName.Data(), fResFileName.Data()); - } + LOGF(info, "Apply ND-correlated smearing: fDoNDSmearing = %d , fResNDHistName = %s", fDoNDSmearing, fResNDHistName.Data()); - fResoPhi_Pos = reinterpret_cast(listRes->FindObject(fResPhiPosHistName)); - if (!fResoPhi_Pos) { - LOGP(fatal, "Could not open {} from file {}", fResPhiPosHistName.Data(), fResFileName.Data()); + if (fDoNDSmearing) { + if (fResType != 0) { + fResoND = reinterpret_cast(listRes->FindObject(fResNDHistName)); + if (!fResoND) { + LOGP(fatal, "Could not open {} from file {}", fResNDHistName.Data(), fResFileName.Data()); + } + fillVecResoND(fResoND); } + } else { + if (fResType != 0) { + fResoPt = reinterpret_cast(listRes->FindObject(fResPtHistName)); + if (!fResoPt) { + LOGP(fatal, "Could not open {} from file {}", fResPtHistName.Data(), fResFileName.Data()); + } + + fResoEta = reinterpret_cast(listRes->FindObject(fResEtaHistName)); + if (!fResoEta) { + LOGP(fatal, "Could not open {} from file {}", fResEtaHistName.Data(), fResFileName.Data()); + } - fResoPhi_Neg = reinterpret_cast(listRes->FindObject(fResPhiNegHistName)); - if (!fResoPhi_Neg) { - LOGP(fatal, "Could not open {} from file {}", fResPhiNegHistName.Data(), fResFileName.Data()); + fResoPhi_Pos = reinterpret_cast(listRes->FindObject(fResPhiPosHistName)); + if (!fResoPhi_Pos) { + LOGP(fatal, "Could not open {} from file {}", fResPhiPosHistName.Data(), fResFileName.Data()); + } + + fResoPhi_Neg = reinterpret_cast(listRes->FindObject(fResPhiNegHistName)); + if (!fResoPhi_Neg) { + LOGP(fatal, "Could not open {} from file {}", fResPhiNegHistName.Data(), fResFileName.Data()); + } + fillVecReso(fResoPt, fVecResoPt); + fillVecReso(fResoEta, fVecResoEta); + fillVecReso(fResoPhi_Pos, fVecResoPhi_Pos); + fillVecReso(fResoPhi_Neg, fVecResoPhi_Neg); } - fillVecReso(fResoPt, fVecResoPt); - fillVecReso(fResoEta, fVecResoEta); - fillVecReso(fResoPhi_Pos, fVecResoPhi_Pos); - fillVecReso(fResoPhi_Neg, fVecResoPhi_Neg); } - if (!fFromCcdb) + if (!fFromCcdb) { delete listRes; + } LOGP(info, "Set efficiency histos"); TList* listEff = new TList(); @@ -237,8 +309,9 @@ class MomentumSmearer } } - if (!fFromCcdb) + if (!fFromCcdb) { delete listEff; + } LOGP(info, "Set DCA histos"); TList* listDCA = new TList(); @@ -280,8 +353,9 @@ class MomentumSmearer fillVecReso(fDCA, fVecDCA); } - if (!fFromCcdb) + if (!fFromCcdb) { delete listDCA; + } fInitialized = true; } @@ -312,13 +386,72 @@ class MomentumSmearer phismeared = phigen; return; } - applySmearing(ptgen, ptgen, ptgen, ptsmeared, fResoPt, fVecResoPt); - applySmearing(ptgen, etagen, 1., etasmeared, fResoEta, fVecResoEta); - if (ch > 0) { - applySmearing(ptgen, phigen, 1., phismeared, fResoPhi_Pos, fVecResoPhi_Pos); + + if (fDoNDSmearing) { + applySmearingND(ch, ptgen, etagen, phigen, ptsmeared, etasmeared, phismeared); } else { - applySmearing(ptgen, phigen, 1., phismeared, fResoPhi_Neg, fVecResoPhi_Neg); + applySmearing(ptgen, ptgen, ptgen, ptsmeared, fResoPt, fVecResoPt); + applySmearing(ptgen, etagen, 1., etasmeared, fResoEta, fVecResoEta); + if (ch > 0) { + applySmearing(ptgen, phigen, 1., phismeared, fResoPhi_Pos, fVecResoPhi_Pos); + } else { + applySmearing(ptgen, phigen, 1., phismeared, fResoPhi_Neg, fVecResoPhi_Neg); + } + } + } + + void applySmearingND(const int ch, const float ptgen, const float etagen, const float phigen, float& ptsmeared, float& etasmeared, float& phismeared) + { + const int npt = fResoND->GetAxis(0)->GetNbins(); + const int neta = fResoND->GetAxis(1)->GetNbins(); + const int nphi = fResoND->GetAxis(2)->GetNbins(); + const int nch = fResoND->GetAxis(3)->GetNbins(); + int ptbin = fResoND->GetAxis(0)->FindBin(ptgen); + int etabin = fResoND->GetAxis(1)->FindBin(etagen); + int phibin = fResoND->GetAxis(2)->FindBin(phigen); + int chbin = fResoND->GetAxis(3)->FindBin(ch); + + // protection + if (ptbin < 1) { + ptbin = 1; + } else if (ptbin > npt) { + ptbin = npt; + } + + // protection + if (etabin < 1) { + etabin = 1; + } else if (etabin > neta) { + etabin = neta; + } + + // protection + if (phibin < 1) { + phibin = 1; + } else if (phibin > nphi) { + phibin = nphi; + } + + // protection + if (chbin < 1) { + chbin = 1; + } else if (chbin > nch) { + chbin = nch; + } + + // ptbin = 1; + // etabin = 1; + // phibin = 1; + // chbin = 1; + + double dpt_rel = 0, deta = 0, dphi = 0; + if (fVecResoND[ptbin - 1][etabin - 1][phibin - 1][chbin - 1]->GetEntries() > 0) { + fVecResoND[ptbin - 1][etabin - 1][phibin - 1][chbin - 1]->GetRandom3(dpt_rel, deta, dphi); } + ptsmeared = ptgen - dpt_rel * ptgen; + etasmeared = etagen - deta; + phismeared = phigen - dphi; + // LOGF(info, "ptgen = %f (GeV/c), etagen = %f, phigen = %f (rad.), ptsmeared = %f (GeV/c), etasmeared = %f, phismeared = %f (rad.)", ptgen, etagen, phigen, ptsmeared, etasmeared, phismeared); } float getEfficiency(float pt, float eta, float phi) @@ -408,7 +541,9 @@ class MomentumSmearer } // setters + void setNDSmearing(bool flag) { fDoNDSmearing = flag; } void setResFileName(TString resFileName) { fResFileName = resFileName; } + void setResNDHistName(TString resNDHistName) { fResNDHistName = resNDHistName; } void setResPtHistName(TString resPtHistName) { fResPtHistName = resPtHistName; } void setResEtaHistName(TString resEtaHistName) { fResEtaHistName = resEtaHistName; } void setResPhiPosHistName(TString resPhiPosHistName) { fResPhiPosHistName = resPhiPosHistName; } @@ -428,7 +563,9 @@ class MomentumSmearer void setTimestamp(int64_t timestamp) { fTimestamp = timestamp; } // getters + bool getNDSmearing() { return fDoNDSmearing; } TString getResFileName() { return fResFileName; } + TString getResNDHistName() { return fResNDHistName; } TString getResPtHistName() { return fResPtHistName; } TString getResEtaHistName() { return fResEtaHistName; } TString getResPhiPosHistName() { return fResPhiPosHistName; } @@ -448,7 +585,9 @@ class MomentumSmearer private: bool fInitialized = false; + bool fDoNDSmearing = false; TString fResFileName; + TString fResNDHistName; TString fResPtHistName; TString fResEtaHistName; TString fResPhiPosHistName; @@ -463,10 +602,12 @@ class MomentumSmearer int fEffType = 0; int fResType = 0; int fDCAType = 0; + THnSparseF* fResoND; TH2F* fResoPt; TH2F* fResoEta; TH2F* fResoPhi_Pos; TH2F* fResoPhi_Neg; + std::vector>>> fVecResoND; std::vector fVecResoPt; std::vector fVecResoEta; std::vector fVecResoPhi_Pos;