From 444859dbbd249c6e8b0ae04fc332c8ba27c1bb12 Mon Sep 17 00:00:00 2001 From: Daiki Sekihata Date: Fri, 14 Nov 2025 17:25:59 +0100 Subject: [PATCH] PWGEM/Dilepton: update polarization for dd. acc. --- PWGEM/Dilepton/Core/Dilepton.h | 435 ++++++++++++++++++++++++++- PWGEM/Dilepton/Utils/EMTrack.h | 33 +- PWGEM/Dilepton/Utils/PairUtilities.h | 100 +++--- 3 files changed, 507 insertions(+), 61 deletions(-) diff --git a/PWGEM/Dilepton/Core/Dilepton.h b/PWGEM/Dilepton/Core/Dilepton.h index 222a690470a..9b0eb989a66 100644 --- a/PWGEM/Dilepton/Core/Dilepton.h +++ b/PWGEM/Dilepton/Core/Dilepton.h @@ -55,6 +55,8 @@ #include #include #include +#include +#include #include #include #include @@ -87,6 +89,7 @@ using FilteredMyMuon = FilteredMyMuons::iterator; using MyEMH_electron = o2::aod::pwgem::dilepton::utils::EventMixingHandler, std::pair, EMTrack>; using MyEMH_muon = o2::aod::pwgem::dilepton::utils::EventMixingHandler, std::pair, EMFwdTrack>; +using MyEMH_pair = o2::aod::pwgem::dilepton::utils::EventMixingHandler, std::pair, std::tuple>; template struct Dilepton { @@ -116,8 +119,8 @@ struct Dilepton { 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 cfgNumContribMin{"cfgNumContribMin", 0, "min. numContrib"}; - Configurable cfgNumContribMax{"cfgNumContribMax", 65000, "max. numContrib"}; + // Configurable cfgNumContribMin{"cfgNumContribMin", 0, "min. numContrib"}; + // Configurable cfgNumContribMax{"cfgNumContribMax", 65000, "max. numContrib"}; Configurable cfgApplyWeightTTCA{"cfgApplyWeightTTCA", false, "flag to apply weighting by 1/N"}; Configurable cfgDCAType{"cfgDCAType", 0, "type of DCA for output. 0:3D, 1:XY, 2:Z, else:3D"}; Configurable cfgUseSignedDCA{"cfgUseSignedDCA", false, "flag to use signs in the DCA calculation"}; @@ -292,24 +295,35 @@ struct Dilepton { Configurable> requiredMFTDisks{"requiredMFTDisks", std::vector{0}, "hit map on MFT disks [0,1,2,3,4]. logical-OR of each double-sided disk"}; } dimuoncuts; + struct : ConfigurableGroup { + std::string prefix = "accBins"; + ConfigurableAxis ConfMllAccBins{"ConfMllAccBins", {40, 0, 4}, "mll bins for acceptance for plarization"}; + ConfigurableAxis ConfPtllAccBins{"ConfPtllAccBins", {100, 0, 10}, "pTll bins for acceptance for plarization"}; + ConfigurableAxis ConfEtallAccBins{"ConEtallAccBins", {30, -1.5f, 1.5f}, "etall bins for acceptance for plarization"}; // pair pseudo-rapidity + ConfigurableAxis ConfPhillAccBins{"ConPhillAccBins", {36, 0.f, 2 * M_PI}, "phill bins for acceptance for plarization"}; // pair pseudo-rapidity + } accBins; + o2::aod::rctsel::RCTFlagsChecker rctChecker; - o2::ccdb::CcdbApi ccdbApi; + // o2::ccdb::CcdbApi ccdbApi; Service ccdb; int mRunNumber; float d_bz; - // o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE; HistogramRegistry fRegistry{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; - static constexpr std::string_view event_cut_types[2] = {"before/", "after/"}; + // static constexpr std::string_view event_cut_types[2] = {"before/", "after/"}; static constexpr std::string_view event_pair_types[2] = {"same/", "mix/"}; + std::mt19937 engine; std::vector cent_bin_edges; std::vector zvtx_bin_edges; std::vector ep_bin_edges; std::vector occ_bin_edges; - int nmod = -1; // this is for flow analysis - int subdet2 = -1; // this is for flow analysis - int subdet3 = -1; // this is for flow analysis + std::vector mll_bin_edges; + std::vector ptll_bin_edges; + std::vector etall_bin_edges; + std::vector phill_bin_edges; + + int nmod = -1; // this is for flow analysis float leptonM1 = 0.f; float leptonM2 = 0.f; @@ -403,6 +417,111 @@ struct Dilepton { emh_pos = new TEMH(ndepth); emh_neg = new TEMH(ndepth); + emh_pair_uls = new MyEMH_pair(ndepth); + emh_pair_lspp = new MyEMH_pair(ndepth); + emh_pair_lsmm = new MyEMH_pair(ndepth); + + if (accBins.ConfMllAccBins.value[0] == VARIABLE_WIDTH) { + mll_bin_edges = std::vector(accBins.ConfMllAccBins.value.begin(), accBins.ConfMllAccBins.value.end()); + mll_bin_edges.erase(mll_bin_edges.begin()); + for (const auto& edge : mll_bin_edges) { + LOGF(info, "VARIABLE_WIDTH: mll_bin_edges = %f", edge); + } + } else { + int nbins = static_cast(accBins.ConfMllAccBins.value[0]); + float xmin = static_cast(accBins.ConfMllAccBins.value[1]); + float xmax = static_cast(accBins.ConfMllAccBins.value[2]); + mll_bin_edges.resize(nbins + 1); + for (int i = 0; i < nbins + 1; i++) { + mll_bin_edges[i] = (xmax - xmin) / (nbins)*i + xmin; + LOGF(info, "FIXED_WIDTH: mll_bin_edges[%d] = %f", i, mll_bin_edges[i]); + } + } + + if (accBins.ConfPtllAccBins.value[0] == VARIABLE_WIDTH) { + ptll_bin_edges = std::vector(accBins.ConfPtllAccBins.value.begin(), accBins.ConfPtllAccBins.value.end()); + ptll_bin_edges.erase(ptll_bin_edges.begin()); + for (const auto& edge : ptll_bin_edges) { + LOGF(info, "VARIABLE_WIDTH: ptll_bin_edges = %f", edge); + } + } else { + int nbins = static_cast(accBins.ConfPtllAccBins.value[0]); + float xmin = static_cast(accBins.ConfPtllAccBins.value[1]); + float xmax = static_cast(accBins.ConfPtllAccBins.value[2]); + ptll_bin_edges.resize(nbins + 1); + for (int i = 0; i < nbins + 1; i++) { + ptll_bin_edges[i] = (xmax - xmin) / (nbins)*i + xmin; + LOGF(info, "FIXED_WIDTH: ptll_bin_edges[%d] = %f", i, ptll_bin_edges[i]); + } + } + + if (accBins.ConfEtallAccBins.value[0] == VARIABLE_WIDTH) { + etall_bin_edges = std::vector(accBins.ConfEtallAccBins.value.begin(), accBins.ConfEtallAccBins.value.end()); + etall_bin_edges.erase(etall_bin_edges.begin()); + for (const auto& edge : etall_bin_edges) { + LOGF(info, "VARIABLE_WIDTH: etall_bin_edges = %f", edge); + } + } else { + int nbins = static_cast(accBins.ConfEtallAccBins.value[0]); + float xmin = static_cast(accBins.ConfEtallAccBins.value[1]); + float xmax = static_cast(accBins.ConfEtallAccBins.value[2]); + etall_bin_edges.resize(nbins + 1); + for (int i = 0; i < nbins + 1; i++) { + etall_bin_edges[i] = (xmax - xmin) / (nbins)*i + xmin; + LOGF(info, "FIXED_WIDTH: etall_bin_edges[%d] = %f", i, etall_bin_edges[i]); + } + } + + if (accBins.ConfPhillAccBins.value[0] == VARIABLE_WIDTH) { + phill_bin_edges = std::vector(accBins.ConfPhillAccBins.value.begin(), accBins.ConfPhillAccBins.value.end()); + phill_bin_edges.erase(phill_bin_edges.begin()); + for (const auto& edge : phill_bin_edges) { + LOGF(info, "VARIABLE_WIDTH: phill_bin_edges = %f", edge); + } + } else { + int nbins = static_cast(accBins.ConfPhillAccBins.value[0]); + float xmin = static_cast(accBins.ConfPhillAccBins.value[1]); + float xmax = static_cast(accBins.ConfPhillAccBins.value[2]); + phill_bin_edges.resize(nbins + 1); + for (int i = 0; i < nbins + 1; i++) { + phill_bin_edges[i] = (xmax - xmin) / (nbins)*i + xmin; + LOGF(info, "FIXED_WIDTH: phill_bin_edges[%d] = %f", i, phill_bin_edges[i]); + } + } + + int nM = mll_bin_edges.size(); + int nPt = ptll_bin_edges.size(); + int nEta = etall_bin_edges.size(); + int nPhi = phill_bin_edges.size(); + + // emhs_pair_uls.resize(nM); + // emhs_pair_lspp.resize(nM); + // emhs_pair_lsmm.resize(nM); + // for (int im = 0;im(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kHFll)) { const AxisSpec axis_dphi_ee{36, -M_PI / 2., 3. / 2. * M_PI, "#Delta#varphi = #varphi_{l1} - #varphi_{l2} (rad.)"}; // for kHFll const AxisSpec axis_deta_ee{40, -2., 2., "#Delta#eta = #eta_{l1} - #eta_{l2}"}; @@ -943,11 +1099,15 @@ struct Dilepton { o2::math_utils::bringToPMPi(dphi_e_ee); float cos_thetaPol = 999, phiPol = 999.f; + auto arrM = std::array{static_cast(v12.Px()), static_cast(v12.Py()), static_cast(v12.Pz()), static_cast(v12.M())}; + auto arrD = t1.sign() * t1.pt() > t2.sign() * t2.pt() ? std::array{t1.px(), t1.py(), t1.pz(), leptonM1} : std::array{t2.px(), t2.py(), t2.pz(), leptonM2}; + if (cfgPolarizationFrame == 0) { - o2::aod::pwgem::dilepton::utils::pairutil::getAngleCS(std::array{t1.px(), t1.py(), t1.pz(), leptonM1}, std::array{t2.px(), t2.py(), t2.pz(), leptonM2}, beamE1, beamE2, beamP1, beamP2, t1.sign(), cos_thetaPol, phiPol); + o2::aod::pwgem::dilepton::utils::pairutil::getAngleCS(arrM, arrD, beamE1, beamE2, beamP1, beamP2, cos_thetaPol, phiPol); } else if (cfgPolarizationFrame == 1) { - o2::aod::pwgem::dilepton::utils::pairutil::getAngleHX(std::array{t1.px(), t1.py(), t1.pz(), leptonM1}, std::array{t2.px(), t2.py(), t2.pz(), leptonM2}, beamE1, beamE2, beamP1, beamP2, t1.sign(), cos_thetaPol, phiPol); + o2::aod::pwgem::dilepton::utils::pairutil::getAngleHX(arrM, arrD, beamE1, beamE2, beamP1, beamP2, cos_thetaPol, phiPol); } + o2::math_utils::bringToPMPi(phiPol); if (t1.sign() * t2.sign() < 0) { // ULS @@ -1000,11 +1160,16 @@ struct Dilepton { } } else if (cfgAnalysisType == static_cast(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kPolarization)) { float cos_thetaPol = 999, phiPol = 999.f; + auto arrM = std::array{static_cast(v12.Px()), static_cast(v12.Py()), static_cast(v12.Pz()), static_cast(v12.M())}; + auto random_sign = std::pow(-1, engine() % 2); // -1^0 = +1 or -1^1 = -1; + auto arrD = t1.sign() * t2.sign() < 0 ? (t1.sign() > 0 ? std::array{t1.px(), t1.py(), t1.pz(), leptonM1} : std::array{t2.px(), t2.py(), t2.pz(), leptonM2}) : (random_sign > 0 ? std::array{t1.px(), t1.py(), t1.pz(), leptonM1} : std::array{t2.px(), t2.py(), t2.pz(), leptonM2}); + if (cfgPolarizationFrame == 0) { - o2::aod::pwgem::dilepton::utils::pairutil::getAngleCS(std::array{t1.px(), t1.py(), t1.pz(), leptonM1}, std::array{t2.px(), t2.py(), t2.pz(), leptonM2}, beamE1, beamE2, beamP1, beamP2, t1.sign(), cos_thetaPol, phiPol); + o2::aod::pwgem::dilepton::utils::pairutil::getAngleCS(arrM, arrD, beamE1, beamE2, beamP1, beamP2, cos_thetaPol, phiPol); } else if (cfgPolarizationFrame == 1) { - o2::aod::pwgem::dilepton::utils::pairutil::getAngleHX(std::array{t1.px(), t1.py(), t1.pz(), leptonM1}, std::array{t2.px(), t2.py(), t2.pz(), leptonM2}, beamE1, beamE2, beamP1, beamP2, t1.sign(), cos_thetaPol, phiPol); + o2::aod::pwgem::dilepton::utils::pairutil::getAngleHX(arrM, arrD, beamE1, beamE2, beamP1, beamP2, cos_thetaPol, phiPol); } + o2::math_utils::bringToPMPi(phiPol); float quadmom = (3.f * std::pow(cos_thetaPol, 2) - 1.f) / 2.f; @@ -1015,6 +1180,55 @@ struct Dilepton { } else if (t1.sign() < 0 && t2.sign() < 0) { // LS-- fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("lsmm/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), cos_thetaPol, phiPol, quadmom, weight); } + + if constexpr (ev_id == 0) { // same event + int mbin = lower_bound(mll_bin_edges.begin(), mll_bin_edges.end(), v12.M()) - mll_bin_edges.begin() - 1; + if (mbin < 0) { + mbin = 0; + } else if (static_cast(mll_bin_edges.size()) - 2 < mbin) { + mbin = static_cast(mll_bin_edges.size()) - 2; + } + + int ptbin = lower_bound(ptll_bin_edges.begin(), ptll_bin_edges.end(), v12.Pt()) - ptll_bin_edges.begin() - 1; + if (ptbin < 0) { + ptbin = 0; + } else if (static_cast(ptll_bin_edges.size()) - 2 < ptbin) { + ptbin = static_cast(ptll_bin_edges.size()) - 2; + } + + int etabin = lower_bound(etall_bin_edges.begin(), etall_bin_edges.end(), v12.Eta()) - etall_bin_edges.begin() - 1; + if (etabin < 0) { + etabin = 0; + } else if (static_cast(etall_bin_edges.size()) - 2 < etabin) { + etabin = static_cast(etall_bin_edges.size()) - 2; + } + + float phi12 = v12.Phi(); + o2::math_utils::bringTo02Pi(phi12); + int phibin = lower_bound(phill_bin_edges.begin(), phill_bin_edges.end(), phi12) - phill_bin_edges.begin() - 1; + if (phibin < 0) { + phibin = 0; + } else if (static_cast(phill_bin_edges.size()) - 2 < phibin) { + phibin = static_cast(phill_bin_edges.size()) - 2; + } + + auto key_df_collision = std::make_pair(ndf, collision.globalIndex()); + float phi12_tmp = v12.Phi(); + o2::math_utils::bringTo02Pi(phi12_tmp); + EMPair empair = EMPair(v12.Pt(), v12.Eta(), phi12_tmp, v12.M(), t1.sign() + t2.sign()); + empair.setPositiveLegPxPyPzM(arrD[0], arrD[1], arrD[2], leptonM1); + // empair.setNegativeLegPtEtaPhiM(t2.pt(), t2.eta(), t2.phi(), leptonM2); + empair.setPairDCA(pair_dca); + auto pair_tmp = std::make_tuple(mbin, ptbin, etabin, phibin, empair); + if (t1.sign() * t2.sign() < 0) { // ULS + emh_pair_uls->AddTrackToEventPool(key_df_collision, pair_tmp); + } else if (t1.sign() > 0 && t2.sign() > 0) { // LS++ + emh_pair_lspp->AddTrackToEventPool(key_df_collision, pair_tmp); + } else if (t1.sign() < 0 && t2.sign() < 0) { // LS-- + emh_pair_lsmm->AddTrackToEventPool(key_df_collision, pair_tmp); + } + } + } else if (cfgAnalysisType == static_cast(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kHFll)) { float dphi = v1.Phi() - v2.Phi(); dphi = RecoDecay::constrainAngle(dphi, -o2::constants::math::PIHalf); @@ -1090,7 +1304,7 @@ struct Dilepton { } Filter collisionFilter_centrality = (cfgCentMin < o2::aod::cent::centFT0M && o2::aod::cent::centFT0M < cfgCentMax) || (cfgCentMin < o2::aod::cent::centFT0A && o2::aod::cent::centFT0A < cfgCentMax) || (cfgCentMin < o2::aod::cent::centFT0C && o2::aod::cent::centFT0C < cfgCentMax); - Filter collisionFilter_numContrib = cfgNumContribMin <= o2::aod::collision::numContrib && o2::aod::collision::numContrib < cfgNumContribMax; + // Filter collisionFilter_numContrib = cfgNumContribMin <= o2::aod::collision::numContrib && o2::aod::collision::numContrib < cfgNumContribMax; Filter collisionFilter_occupancy_track = eventcuts.cfgTrackOccupancyMin <= o2::aod::evsel::trackOccupancyInTimeRange && o2::aod::evsel::trackOccupancyInTimeRange < eventcuts.cfgTrackOccupancyMax; Filter collisionFilter_occupancy_ft0c = eventcuts.cfgFT0COccupancyMin <= o2::aod::evsel::ft0cOccupancyInTimeRange && o2::aod::evsel::ft0cOccupancyInTimeRange < eventcuts.cfgFT0COccupancyMax; using FilteredMyCollisions = soa::Filtered; @@ -1128,6 +1342,13 @@ struct Dilepton { TEMH* emh_pos = nullptr; TEMH* emh_neg = nullptr; + MyEMH_pair* emh_pair_uls = nullptr; + MyEMH_pair* emh_pair_lspp = nullptr; + MyEMH_pair* emh_pair_lsmm = nullptr; + + // std::vector>>> emhs_pair_uls; // 4D{m, pt, eta, phi} + // std::vector>>> emhs_pair_lspp; // 4D{m, pt, eta, phi} + // std::vector>>> emhs_pair_lsmm; // 4D{m, pt, eta, phi} std::map, uint64_t> map_mixed_eventId_to_globalBC; std::vector used_trackIds_per_col; @@ -1279,7 +1500,7 @@ struct Dilepton { // LOGF(info, "collision.globalIndex() = %d, collision.posZ() = %f, centrality = %f, ep2 = %f, collision.trackOccupancyInTimeRange() = %d, zbin = %d, centbin = %d, epbin = %d, occbin = %d", collision.globalIndex(), collision.posZ(), centrality, ep2, collision.trackOccupancyInTimeRange(), zbin, centbin, epbin, occbin); std::tuple key_bin = std::make_tuple(zbin, centbin, epbin, occbin); - std::pair key_df_collision = std::make_pair(ndf, collision.globalIndex()); + std::pair key_df_collision = std::make_pair(ndf, collision.globalIndex()); // this gives the current event. // make a vector of selected photons in this collision. auto selected_posTracks_in_this_event = emh_pos->GetTracksPerCollision(key_df_collision); @@ -1332,11 +1553,195 @@ struct Dilepton { } } // end of loop over mixed event pool + if (cfgAnalysisType == static_cast(o2::aod::pwgem::dilepton::utils::pairutil::DileptonAnalysisType::kPolarization)) { // only for polarization + auto selected_pairs_uls_in_this_event = emh_pair_uls->GetTracksPerCollision(key_df_collision); + auto selected_pairs_lspp_in_this_event = emh_pair_lspp->GetTracksPerCollision(key_df_collision); + auto selected_pairs_lsmm_in_this_event = emh_pair_lsmm->GetTracksPerCollision(key_df_collision); + auto collisionIds_in_mixing_pool = emh_pair_uls->GetCollisionIdsFromEventPool(key_bin); + float weight = 1.f; + + for (const auto& mix_dfId_collisionId : collisionIds_in_mixing_pool) { + auto pairs_uls_from_event_pool = emh_pair_uls->GetTracksPerCollision(mix_dfId_collisionId); + auto pairs_lspp_from_event_pool = emh_pair_lspp->GetTracksPerCollision(mix_dfId_collisionId); + auto pairs_lsmm_from_event_pool = emh_pair_lsmm->GetTracksPerCollision(mix_dfId_collisionId); + + for (const auto& pair1 : selected_pairs_uls_in_this_event) { // ULS mix + auto empair1 = std::get<4>(pair1); + auto v_pos = empair1.getPositiveLeg(); // pt, eta, phi, M + // auto v_neg = empair1.getNegativeLeg(); // pt, eta, phi, M + auto arrD = std::array{static_cast(v_pos.Px()), static_cast(v_pos.Py()), static_cast(v_pos.Pz()), leptonM1}; + + int mbin = lower_bound(mll_bin_edges.begin(), mll_bin_edges.end(), empair1.mass()) - mll_bin_edges.begin() - 1; + if (mbin < 0) { + mbin = 0; + } else if (static_cast(mll_bin_edges.size()) - 2 < mbin) { + mbin = static_cast(mll_bin_edges.size()) - 2; + } + + int ptbin = lower_bound(ptll_bin_edges.begin(), ptll_bin_edges.end(), empair1.pt()) - ptll_bin_edges.begin() - 1; + if (ptbin < 0) { + ptbin = 0; + } else if (static_cast(ptll_bin_edges.size()) - 2 < ptbin) { + ptbin = static_cast(ptll_bin_edges.size()) - 2; + } + + int etabin = lower_bound(etall_bin_edges.begin(), etall_bin_edges.end(), empair1.eta()) - etall_bin_edges.begin() - 1; + if (etabin < 0) { + etabin = 0; + } else if (static_cast(etall_bin_edges.size()) - 2 < etabin) { + etabin = static_cast(etall_bin_edges.size()) - 2; + } + + int phibin = lower_bound(phill_bin_edges.begin(), phill_bin_edges.end(), empair1.phi()) - phill_bin_edges.begin() - 1; + if (phibin < 0) { + phibin = 0; + } else if (static_cast(phill_bin_edges.size()) - 2 < phibin) { + phibin = static_cast(phill_bin_edges.size()) - 2; + } + + for (const auto& pair2 : std::views::filter(pairs_uls_from_event_pool, [&mbin, &ptbin, &etabin, &phibin](std::tuple t) { return std::get<0>(t) == mbin && std::get<1>(t) == ptbin && std::get<2>(t) == etabin && std::get<3>(t) == phibin; })) { + auto empair2 = std::get<4>(pair2); + auto arrM = std::array{static_cast(empair2.px()), static_cast(empair2.py()), static_cast(empair2.pz()), static_cast(empair2.mass())}; + + float cos_thetaPol = 999, phiPol = 999.f; + if (cfgPolarizationFrame == 0) { + o2::aod::pwgem::dilepton::utils::pairutil::getAngleCS(arrM, arrD, beamE1, beamE2, beamP1, beamP2, cos_thetaPol, phiPol); + } else if (cfgPolarizationFrame == 1) { + o2::aod::pwgem::dilepton::utils::pairutil::getAngleHX(arrM, arrD, beamE1, beamE2, beamP1, beamP2, cos_thetaPol, phiPol); + } + o2::math_utils::bringToPMPi(phiPol); + float quadmom = (3.f * std::pow(cos_thetaPol, 2) - 1.f) / 2.f; + fRegistry.fill(HIST("Pair/mix/uls/hsAcc"), empair1.mass(), empair1.pt(), empair1.getPairDCA(), empair1.rapidity(), cos_thetaPol, phiPol, quadmom, weight); + } + } // end of ULS + + for (const auto& pair1 : selected_pairs_lspp_in_this_event) { // LS++ + auto empair1 = std::get<4>(pair1); + auto v_pos = empair1.getPositiveLeg(); // pt, eta, phi, M + // auto v_neg = empair1.getNegativeLeg(); // pt, eta, phi, M + // auto arrD = +1 * v_pos.Pt() > +1 * v_neg.Pt() ? std::array{static_cast(v_pos.Px()), static_cast(v_pos.Py()), static_cast(v_pos.Pz()), leptonM1} : std::array{static_cast(v_neg.Px()), static_cast(v_neg.Py()), static_cast(v_neg.Pz()), leptonM2}; + auto arrD = std::array{static_cast(v_pos.Px()), static_cast(v_pos.Py()), static_cast(v_pos.Pz()), leptonM1}; + + int mbin = lower_bound(mll_bin_edges.begin(), mll_bin_edges.end(), empair1.mass()) - mll_bin_edges.begin() - 1; + if (mbin < 0) { + mbin = 0; + } else if (static_cast(mll_bin_edges.size()) - 2 < mbin) { + mbin = static_cast(mll_bin_edges.size()) - 2; + } + + int ptbin = lower_bound(ptll_bin_edges.begin(), ptll_bin_edges.end(), empair1.pt()) - ptll_bin_edges.begin() - 1; + if (ptbin < 0) { + ptbin = 0; + } else if (static_cast(ptll_bin_edges.size()) - 2 < ptbin) { + ptbin = static_cast(ptll_bin_edges.size()) - 2; + } + + int etabin = lower_bound(etall_bin_edges.begin(), etall_bin_edges.end(), empair1.eta()) - etall_bin_edges.begin() - 1; + if (etabin < 0) { + etabin = 0; + } else if (static_cast(etall_bin_edges.size()) - 2 < etabin) { + etabin = static_cast(etall_bin_edges.size()) - 2; + } + + int phibin = lower_bound(phill_bin_edges.begin(), phill_bin_edges.end(), empair1.phi()) - phill_bin_edges.begin() - 1; + if (phibin < 0) { + phibin = 0; + } else if (static_cast(phill_bin_edges.size()) - 2 < phibin) { + phibin = static_cast(phill_bin_edges.size()) - 2; + } + + for (const auto& pair2 : std::views::filter(pairs_lspp_from_event_pool, [&mbin, &ptbin, &etabin, &phibin](std::tuple t) { return std::get<0>(t) == mbin && std::get<1>(t) == ptbin && std::get<2>(t) == etabin && std::get<3>(t) == phibin; })) { + auto empair2 = std::get<4>(pair2); + auto arrM = std::array{static_cast(empair2.px()), static_cast(empair2.py()), static_cast(empair2.pz()), static_cast(empair2.mass())}; + + float cos_thetaPol = 999, phiPol = 999.f; + if (cfgPolarizationFrame == 0) { + o2::aod::pwgem::dilepton::utils::pairutil::getAngleCS(arrM, arrD, beamE1, beamE2, beamP1, beamP2, cos_thetaPol, phiPol); + } else if (cfgPolarizationFrame == 1) { + o2::aod::pwgem::dilepton::utils::pairutil::getAngleHX(arrM, arrD, beamE1, beamE2, beamP1, beamP2, cos_thetaPol, phiPol); + } + o2::math_utils::bringToPMPi(phiPol); + float quadmom = (3.f * std::pow(cos_thetaPol, 2) - 1.f) / 2.f; + fRegistry.fill(HIST("Pair/mix/lspp/hsAcc"), empair1.mass(), empair1.pt(), empair1.getPairDCA(), empair1.rapidity(), cos_thetaPol, phiPol, quadmom, weight); + } + } // end of LS++ + + for (const auto& pair1 : selected_pairs_lsmm_in_this_event) { // LS++ + auto empair1 = std::get<4>(pair1); + auto v_pos = empair1.getPositiveLeg(); // pt, eta, phi, M + // auto v_neg = empair1.getNegativeLeg(); // pt, eta, phi, M + // auto arrD = -1 * v_pos.Pt() > -1 * v_neg.Pt() ? std::array{static_cast(v_pos.Px()), static_cast(v_pos.Py()), static_cast(v_pos.Pz()), leptonM1} : std::array{static_cast(v_neg.Px()), static_cast(v_neg.Py()), static_cast(v_neg.Pz()), leptonM2}; + auto arrD = std::array{static_cast(v_pos.Px()), static_cast(v_pos.Py()), static_cast(v_pos.Pz()), leptonM1}; + + int mbin = lower_bound(mll_bin_edges.begin(), mll_bin_edges.end(), empair1.mass()) - mll_bin_edges.begin() - 1; + if (mbin < 0) { + mbin = 0; + } else if (static_cast(mll_bin_edges.size()) - 2 < mbin) { + mbin = static_cast(mll_bin_edges.size()) - 2; + } + + int ptbin = lower_bound(ptll_bin_edges.begin(), ptll_bin_edges.end(), empair1.pt()) - ptll_bin_edges.begin() - 1; + if (ptbin < 0) { + ptbin = 0; + } else if (static_cast(ptll_bin_edges.size()) - 2 < ptbin) { + ptbin = static_cast(ptll_bin_edges.size()) - 2; + } + + int etabin = lower_bound(etall_bin_edges.begin(), etall_bin_edges.end(), empair1.eta()) - etall_bin_edges.begin() - 1; + if (etabin < 0) { + etabin = 0; + } else if (static_cast(etall_bin_edges.size()) - 2 < etabin) { + etabin = static_cast(etall_bin_edges.size()) - 2; + } + + int phibin = lower_bound(phill_bin_edges.begin(), phill_bin_edges.end(), empair1.phi()) - phill_bin_edges.begin() - 1; + if (phibin < 0) { + phibin = 0; + } else if (static_cast(phill_bin_edges.size()) - 2 < phibin) { + phibin = static_cast(phill_bin_edges.size()) - 2; + } + + for (const auto& pair2 : std::views::filter(pairs_lsmm_from_event_pool, [&mbin, &ptbin, &etabin, &phibin](std::tuple t) { return std::get<0>(t) == mbin && std::get<1>(t) == ptbin && std::get<2>(t) == etabin && std::get<3>(t) == phibin; })) { + auto empair2 = std::get<4>(pair2); + auto arrM = std::array{static_cast(empair2.px()), static_cast(empair2.py()), static_cast(empair2.pz()), static_cast(empair2.mass())}; + + float cos_thetaPol = 999, phiPol = 999.f; + if (cfgPolarizationFrame == 0) { + o2::aod::pwgem::dilepton::utils::pairutil::getAngleCS(arrM, arrD, beamE1, beamE2, beamP1, beamP2, cos_thetaPol, phiPol); + } else if (cfgPolarizationFrame == 1) { + o2::aod::pwgem::dilepton::utils::pairutil::getAngleHX(arrM, arrD, beamE1, beamE2, beamP1, beamP2, cos_thetaPol, phiPol); + } + o2::math_utils::bringToPMPi(phiPol); + float quadmom = (3.f * std::pow(cos_thetaPol, 2) - 1.f) / 2.f; + fRegistry.fill(HIST("Pair/mix/lsmm/hsAcc"), empair1.mass(), empair1.pt(), empair1.getPairDCA(), empair1.rapidity(), cos_thetaPol, phiPol, quadmom, weight); + } + } // end of LS++ + + } // end of loop over mixed event pool + + } // end of if polarization + if (nuls > 0 || nlspp > 0 || nlsmm > 0) { map_mixed_eventId_to_globalBC[key_df_collision] = collision.globalBC(); emh_pos->AddCollisionIdAtLast(key_bin, key_df_collision); emh_neg->AddCollisionIdAtLast(key_bin, key_df_collision); - } + emh_pair_uls->AddCollisionIdAtLast(key_bin, key_df_collision); + emh_pair_lspp->AddCollisionIdAtLast(key_bin, key_df_collision); + emh_pair_lsmm->AddCollisionIdAtLast(key_bin, key_df_collision); + + // for (int im = 0;imAddCollisionIdAtLast(key_bin, key_df_collision); + // emhs_pair_lspp[im][ipt][ieta][iphi]->AddCollisionIdAtLast(key_bin, key_df_collision); + // emhs_pair_lsmm[im][ipt][ieta][iphi]->AddCollisionIdAtLast(key_bin, key_df_collision); + // } + // } + // } + // } + + } // end of if pair exist } // end of collision loop diff --git a/PWGEM/Dilepton/Utils/EMTrack.h b/PWGEM/Dilepton/Utils/EMTrack.h index 9c0a6f0691e..29bfdd94929 100644 --- a/PWGEM/Dilepton/Utils/EMTrack.h +++ b/PWGEM/Dilepton/Utils/EMTrack.h @@ -52,6 +52,7 @@ class EMTrack float cZY() const { return fCZY; } float cZZ() const { return fCZZ; } + float rapidity() const { return std::log((std::sqrt(std::pow(fMass, 2) + std::pow(fPt * std::cosh(fEta), 2)) + fPt * std::sinh(fEta)) / std::sqrt(std::pow(fMass, 2) + std::pow(fPt, 2))); } float p() const { return fPt * std::cosh(fEta); } float px() const { return fPt * std::cos(fPhi); } float py() const { return fPt * std::sin(fPhi); } @@ -150,8 +151,9 @@ class EMTrackWithCov : public EMTrack class EMPair : public EMTrack { public: - EMPair(float pt, float eta, float phi, float mass) : EMTrack(pt, eta, phi, mass, 0, 0, 0, 0, 0, 0) + EMPair(float pt, float eta, float phi, float mass, int8_t charge = 0) : EMTrack(pt, eta, phi, mass, charge, 0, 0, 0, 0, 0) { + fPairDCA = 999.f; fVPos = ROOT::Math::PtEtaPhiMVector(0, 0, 0, 0); fVNeg = ROOT::Math::PtEtaPhiMVector(0, 0, 0, 0); fVx = 0.f; @@ -161,6 +163,9 @@ class EMPair : public EMTrack ~EMPair() {} + void setPairDCA(float dca) { fPairDCA = dca; } + float getPairDCA() const { return fPairDCA; } + void setPositiveLegPtEtaPhiM(float pt, float eta, float phi, float m) { fVPos.SetPt(pt); @@ -176,6 +181,29 @@ class EMPair : public EMTrack fVNeg.SetM(m); } + void setPositiveLegPxPyPzM(float px, float py, float pz, float m) + { + float pt = std::sqrt(px * px + py * py); + float eta = std::atanh(pz / sqrt(std::pow(px, 2) + std::pow(py, 2) + std::pow(pz, 2))); + float phi = std::atan2(py, px); + + fVPos.SetPt(pt); + fVPos.SetEta(eta); + fVPos.SetPhi(phi); + fVPos.SetM(m); + } + void setNegativeLegPxPyPzM(float px, float py, float pz, float m) + { + float pt = std::sqrt(px * px + py * py); + float eta = std::atanh(pz / std::sqrt(pow(px, 2) + std::pow(py, 2) + std::pow(pz, 2))); + float phi = std::atan2(py, px); + + fVNeg.SetPt(pt); + fVNeg.SetEta(eta); + fVNeg.SetPhi(phi); + fVNeg.SetM(m); + } + ROOT::Math::PtEtaPhiMVector getPositiveLeg() const { return fVPos; } ROOT::Math::PtEtaPhiMVector getNegativeLeg() const { return fVNeg; } @@ -189,10 +217,11 @@ class EMPair : public EMTrack float vy() const { return fVy; } float vz() const { return fVz; } float v0radius() const { return std::sqrt(std::pow(fVx, 2) + std::pow(fVy, 2)); } - float eta_cp() const { return std::atanh(fVz / sqrt(pow(fVx, 2) + pow(fVy, 2) + pow(fVz, 2))); } + float eta_cp() const { return std::atanh(fVz / std::sqrt(std::pow(fVx, 2) + std::pow(fVy, 2) + std::pow(fVz, 2))); } float phi_cp() const { return std::atan2(fVy, fVx); } protected: + float fPairDCA; ROOT::Math::PtEtaPhiMVector fVPos; ROOT::Math::PtEtaPhiMVector fVNeg; diff --git a/PWGEM/Dilepton/Utils/PairUtilities.h b/PWGEM/Dilepton/Utils/PairUtilities.h index 9e59cd72704..0758289bbee 100644 --- a/PWGEM/Dilepton/Utils/PairUtilities.h +++ b/PWGEM/Dilepton/Utils/PairUtilities.h @@ -73,71 +73,86 @@ using SMatrix5 = ROOT::Math::SVector; //_______________________________________________________________________ template -void getAngleHX(std::array const& t1, std::array const& t2, const float beamE1, const float beamE2, const float beamP1, const float beamP2, const int8_t c1, float& cos_thetaHX, float& phiHX) +void getAngleHX(std::array const& tM, std::array const& tD, const float beamE1, const float beamE2, const float beamP1, const float beamP2, float& cos_thetaHX, float& phiHX) { - // t1[0] = px, t1[1] = py, t1[2] = pz, t1[3] = mass; - ROOT::Math::PxPyPzEVector v1(t1[0], t1[1], t1[2], std::sqrt(std::pow(t1[0], 2) + std::pow(t1[1], 2) + std::pow(t1[2], 2) + std::pow(t1[3], 2))); - ROOT::Math::PxPyPzEVector v2(t2[0], t2[1], t2[2], std::sqrt(std::pow(t2[0], 2) + std::pow(t2[1], 2) + std::pow(t2[2], 2) + std::pow(t2[3], 2))); - ROOT::Math::PxPyPzEVector v12 = v1 + v2; + // tM is mother momentum. tD is daugher momentum. Boost daugher to mother's rest frame. + ROOT::Math::PxPyPzEVector vM(tM[0], tM[1], tM[2], std::sqrt(std::pow(tM[0], 2) + std::pow(tM[1], 2) + std::pow(tM[2], 2) + std::pow(tM[3], 2))); + ROOT::Math::PxPyPzEVector vD(tD[0], tD[1], tD[2], std::sqrt(std::pow(tD[0], 2) + std::pow(tD[1], 2) + std::pow(tD[2], 2) + std::pow(tD[3], 2))); ROOT::Math::PxPyPzEVector Beam1(0., 0., -beamP1, beamE1); ROOT::Math::PxPyPzEVector Beam2(0., 0., beamP2, beamE2); - // Boost to center of mass frame. i.e. rest frame of pair - ROOT::Math::Boost boostv12{v12.BoostToCM()}; - ROOT::Math::XYZVectorF v1_CM{(boostv12(v1).Vect()).Unit()}; - ROOT::Math::XYZVectorF v2_CM{(boostv12(v2).Vect()).Unit()}; - ROOT::Math::XYZVectorF Beam1_CM{(boostv12(Beam1).Vect()).Unit()}; - ROOT::Math::XYZVectorF Beam2_CM{(boostv12(Beam2).Vect()).Unit()}; - // LOGF(info, "boostv12(v12).Vect().X() = %f, boostv12(v12).Vect().Y() = %f, boostv12(v12).Vect().Z() = %f", boostv12(v12).Vect().X(), boostv12(v12).Vect().Y(), boostv12(v12).Vect().Z()); // expected to be (0,0,0) + // Boost to pair rest frame + ROOT::Math::Boost boostvM{vM.BoostToCM()}; + ROOT::Math::XYZVectorF vD_PRF{(boostvM(vD).Vect()).Unit()}; + ROOT::Math::XYZVectorF Beam1_PRF{(boostvM(Beam1).Vect()).Unit()}; + ROOT::Math::XYZVectorF Beam2_PRF{(boostvM(Beam2).Vect()).Unit()}; // Helicity frame - ROOT::Math::XYZVectorF zaxis_HX{(v12.Vect()).Unit()}; - ROOT::Math::XYZVectorF yaxis_HX{(Beam1_CM.Cross(Beam2_CM)).Unit()}; + ROOT::Math::XYZVectorF zaxis_HX{(vM.Vect()).Unit()}; + ROOT::Math::XYZVectorF yaxis_HX{(Beam1_PRF.Cross(Beam2_PRF)).Unit()}; ROOT::Math::XYZVectorF xaxis_HX{(yaxis_HX.Cross(zaxis_HX)).Unit()}; - // pdgCode : 11 for electron, -11 for positron - // pdgCode : 13 for negative muon, -13 for positive muon - // LOGF(info, "zaxis_HX.Dot(v1_CM) = %f , zaxis_HX.Dot(v2_CM) = %f", zaxis_HX.Dot(v1_CM), zaxis_HX.Dot(v2_CM)); // absolute value is identical. only sign is opposite. - - cos_thetaHX = c1 > 0 ? zaxis_HX.Dot(v1_CM) : zaxis_HX.Dot(v2_CM); - phiHX = c1 > 0 ? std::atan2(yaxis_HX.Dot(v1_CM), xaxis_HX.Dot(v1_CM)) : std::atan2(yaxis_HX.Dot(v2_CM), xaxis_HX.Dot(v2_CM)); + cos_thetaHX = zaxis_HX.Dot(vD_PRF); + phiHX = std::atan2(yaxis_HX.Dot(vD_PRF), xaxis_HX.Dot(vD_PRF)); } - -//_______________________________________________________________________ //_______________________________________________________________________ template -void getAngleCS(std::array const& t1, std::array const& t2, const float beamE1, const float beamE2, const float beamP1, const float beamP2, const int8_t c1, float& cos_thetaCS, float& phiCS) +void getAngleCS(std::array const& tM, std::array const& tD, const float beamE1, const float beamE2, const float beamP1, const float beamP2, float& cos_thetaCS, float& phiCS) { - // t1[0] = px, t1[1] = py, t1[2] = pz, t1[3] = mass; - ROOT::Math::PxPyPzEVector v1(t1[0], t1[1], t1[2], std::sqrt(std::pow(t1[0], 2) + std::pow(t1[1], 2) + std::pow(t1[2], 2) + std::pow(t1[3], 2))); - ROOT::Math::PxPyPzEVector v2(t2[0], t2[1], t2[2], std::sqrt(std::pow(t2[0], 2) + std::pow(t2[1], 2) + std::pow(t2[2], 2) + std::pow(t2[3], 2))); - ROOT::Math::PxPyPzEVector v12 = v1 + v2; + // tM is mother momentum. tD is daugher momentum. Boost daugher to mother's rest frame. + ROOT::Math::PxPyPzEVector vM(tM[0], tM[1], tM[2], std::sqrt(std::pow(tM[0], 2) + std::pow(tM[1], 2) + std::pow(tM[2], 2) + std::pow(tM[3], 2))); + ROOT::Math::PxPyPzEVector vD(tD[0], tD[1], tD[2], std::sqrt(std::pow(tD[0], 2) + std::pow(tD[1], 2) + std::pow(tD[2], 2) + std::pow(tD[3], 2))); ROOT::Math::PxPyPzEVector Beam1(0., 0., -beamP1, beamE1); ROOT::Math::PxPyPzEVector Beam2(0., 0., beamP2, beamE2); - // Boost to center of mass frame. i.e. rest frame of pair - ROOT::Math::Boost boostv12{v12.BoostToCM()}; - ROOT::Math::XYZVectorF v1_CM{(boostv12(v1).Vect()).Unit()}; - ROOT::Math::XYZVectorF v2_CM{(boostv12(v2).Vect()).Unit()}; - ROOT::Math::XYZVectorF Beam1_CM{(boostv12(Beam1).Vect()).Unit()}; - ROOT::Math::XYZVectorF Beam2_CM{(boostv12(Beam2).Vect()).Unit()}; - // LOGF(info, "boostv12(v12).Vect().X() = %f, boostv12(v12).Vect().Y() = %f, boostv12(v12).Vect().Z() = %f", boostv12(v12).Vect().X(), boostv12(v12).Vect().Y(), boostv12(v12).Vect().Z()); // expected to be (0,0,0) + // Boost to pair rest frame + ROOT::Math::Boost boostvM{vM.BoostToCM()}; + ROOT::Math::XYZVectorF vD_PRF{(boostvM(vD).Vect()).Unit()}; + ROOT::Math::XYZVectorF Beam1_PRF{(boostvM(Beam1).Vect()).Unit()}; + ROOT::Math::XYZVectorF Beam2_PRF{(boostvM(Beam2).Vect()).Unit()}; // Collins-Soper frame - ROOT::Math::XYZVectorF zaxis_CS{((Beam1_CM.Unit() - Beam2_CM.Unit()).Unit())}; - ROOT::Math::XYZVectorF yaxis_CS{(Beam1_CM.Cross(Beam2_CM)).Unit()}; + ROOT::Math::XYZVectorF zaxis_CS{((Beam1_PRF.Unit() - Beam2_PRF.Unit()).Unit())}; + ROOT::Math::XYZVectorF yaxis_CS{(Beam1_PRF.Cross(Beam2_PRF)).Unit()}; ROOT::Math::XYZVectorF xaxis_CS{(yaxis_CS.Cross(zaxis_CS)).Unit()}; - // pdgCode : 11 for electron, -11 for positron - // pdgCode : 13 for negative muon, -13 for positive muon - // LOGF(info, "zaxis_CS.Dot(v1_CM) = %f , zaxis_CS.Dot(v2_CM) = %f", zaxis_CS.Dot(v1_CM), zaxis_CS.Dot(v2_CM)); // absolute value is identical. only sign is opposite. + cos_thetaCS = zaxis_CS.Dot(vD_PRF); + phiCS = std::atan2(yaxis_CS.Dot(vD_PRF), xaxis_CS.Dot(vD_PRF)); +} +//_______________________________________________________________________ +template +void getAngleHX(std::array const& t1, std::array const& t2, const float beamE1, const float beamE2, const float beamP1, const float beamP2, const int8_t c1, float& cos_thetaHX, float& phiHX) +{ + // t1[0] = px, t1[1] = py, t1[2] = pz, t1[3] = mass; + ROOT::Math::PxPyPzEVector v1(t1[0], t1[1], t1[2], std::sqrt(std::pow(t1[0], 2) + std::pow(t1[1], 2) + std::pow(t1[2], 2) + std::pow(t1[3], 2))); + ROOT::Math::PxPyPzEVector v2(t2[0], t2[1], t2[2], std::sqrt(std::pow(t2[0], 2) + std::pow(t2[1], 2) + std::pow(t2[2], 2) + std::pow(t2[3], 2))); + ROOT::Math::PxPyPzEVector v12 = v1 + v2; + std::array arrM{static_cast(v12.Px()), static_cast(v12.Py()), static_cast(v12.Pz()), static_cast(v12.M())}; - cos_thetaCS = c1 > 0 ? zaxis_CS.Dot(v1_CM) : zaxis_CS.Dot(v2_CM); - phiCS = c1 > 0 ? std::atan2(yaxis_CS.Dot(v1_CM), xaxis_CS.Dot(v1_CM)) : std::atan2(yaxis_CS.Dot(v2_CM), xaxis_CS.Dot(v2_CM)); + if (c1 > 0) { + getAngleHX(arrM, t1, beamE1, beamE2, beamP1, beamP2, cos_thetaHX, phiHX); + } else { + getAngleHX(arrM, t2, beamE1, beamE2, beamP1, beamP2, cos_thetaHX, phiHX); + } } +//_______________________________________________________________________ +template +void getAngleCS(std::array const& t1, std::array const& t2, const float beamE1, const float beamE2, const float beamP1, const float beamP2, const int8_t c1, float& cos_thetaCS, float& phiCS) +{ + // t1[0] = px, t1[1] = py, t1[2] = pz, t1[3] = mass; + ROOT::Math::PxPyPzEVector v1(t1[0], t1[1], t1[2], std::sqrt(std::pow(t1[0], 2) + std::pow(t1[1], 2) + std::pow(t1[2], 2) + std::pow(t1[3], 2))); + ROOT::Math::PxPyPzEVector v2(t2[0], t2[1], t2[2], std::sqrt(std::pow(t2[0], 2) + std::pow(t2[1], 2) + std::pow(t2[2], 2) + std::pow(t2[3], 2))); + ROOT::Math::PxPyPzEVector v12 = v1 + v2; + std::array arrM{static_cast(v12.Px()), static_cast(v12.Py()), static_cast(v12.Pz()), static_cast(v12.M())}; + if (c1 > 0) { + getAngleCS(arrM, t1, beamE1, beamE2, beamP1, beamP2, cos_thetaCS, phiCS); + } else { + getAngleCS(arrM, t2, beamE1, beamE2, beamP1, beamP2, cos_thetaCS, phiCS); + } +} //_______________________________________________________________________ template bool isSVFound(TDCAFitter fitter, TCollision const& collision, TTrack const& t1, TTrack const& t2, float& pca, float& lxy, float& cosPA) @@ -172,7 +187,6 @@ bool isSVFound(TDCAFitter fitter, TCollision const& collision, TTrack const& t1, return false; } } - //_______________________________________________________________________ // function call without cosPA template @@ -344,13 +358,11 @@ inline float pairDCAQuadSum(const float dca1, const float dca2) { return std::sqrt((dca1 * dca1 + dca2 * dca2) / 2.); } - //_______________________________________________________________________ inline float pairDCASignQuadSum(const float dca1, const float dca2, const int8_t charge1, const int8_t charge2) { return charge1 * charge2 * TMath::Sign(1., dca1) * TMath::Sign(1., dca2) * std::sqrt((dca1 * dca1 + dca2 * dca2) / 2.); } - //_______________________________________________________________________ } // namespace o2::aod::pwgem::dilepton::utils::pairutil #endif // PWGEM_DILEPTON_UTILS_PAIRUTILITIES_H_