Skip to content

Commit 5b64739

Browse files
dsekihatdsekihat
committed
PWGEM/Dilepton: update DileptonMC for acc
1 parent 992ab80 commit 5b64739

File tree

1 file changed

+149
-10
lines changed

1 file changed

+149
-10
lines changed

‎PWGEM/Dilepton/Core/DileptonMC.h‎

Lines changed: 149 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -369,7 +369,8 @@ struct DileptonMC {
369369
const AxisSpec axis_pt{ConfPtllBins, pair_pt_axis_title};
370370
const AxisSpec axis_y{ConfYllBins, pair_y_axis_title};
371371
const AxisSpec axis_dca{ConfDCAllBins, pair_dca_axis_title};
372-
const AxisSpec axis_pt_meson{ConfPtllBins, "p_{T} (GeV/c)"}; // for omega, phi meson pT spectra
372+
const AxisSpec axis_pt_meson{ConfPtllBins, "p_{T}^{VM} (GeV/c)"}; // for omega, phi meson pT spectra
373+
const AxisSpec axis_y_meson{ConfYllBins, "y^{VM}"}; // for omega, phi meson pT spectra
373374

374375
const AxisSpec axis_dca_narrow{ConfDCAllNarrowBins, pair_dca_axis_title};
375376
const AxisSpec axis_dpt{ConfDPtBins, "#Delta p_{T,1}^{gen-rec} + #Delta p_{T,2}^{gen-rec} (GeV/c)"};
@@ -415,9 +416,6 @@ struct DileptonMC {
415416
// fRegistry.addClone("Generated/sm/PromptPi0/", "Generated/sm/Upsilon2S/");
416417
// fRegistry.addClone("Generated/sm/PromptPi0/", "Generated/sm/Upsilon3S/");
417418

418-
fRegistry.add("Generated/sm/Omega2ll/uls/hPtY", "pT of #omega meson", kTH2D, {axis_y, axis_pt_meson}, true);
419-
fRegistry.add("Generated/sm/Phi2ll/uls/hPtY", "pT of #phi meson", kTH2D, {axis_y, axis_pt_meson}, true);
420-
421419
fRegistry.add("Generated/ccbar/c2l_c2l/uls/hs", "generated dilepton", kTHnSparseD, {axis_mass, axis_pt, axis_y, axis_dphi_ee, axis_deta_ee, axis_cos_theta_pol, axis_phi_pol, axis_quadmom, axis_aco, axis_asym_pt, axis_dphi_e_ee}, true);
422420
fRegistry.addClone("Generated/ccbar/c2l_c2l/uls/", "Generated/ccbar/c2l_c2l/lspp/");
423421
fRegistry.addClone("Generated/ccbar/c2l_c2l/uls/", "Generated/ccbar/c2l_c2l/lsmm/");
@@ -461,6 +459,14 @@ struct DileptonMC {
461459
}
462460
}
463461

462+
// evaluate acceptance for polarization
463+
fRegistry.add("Generated/VM/All/Phi/hs", "gen. VM #rightarrow ll", kTHnSparseD, {axis_mass, axis_pt, axis_y, axis_cos_theta_pol, axis_phi_pol, axis_quadmom}, true);
464+
fRegistry.addClone("Generated/VM/All/Phi/", "Generated/VM/All/Rho/");
465+
fRegistry.addClone("Generated/VM/All/Phi/", "Generated/VM/All/Omega/");
466+
fRegistry.addClone("Generated/VM/All/Phi/", "Generated/VM/All/PromptJPsi/");
467+
fRegistry.addClone("Generated/VM/All/Phi/", "Generated/VM/All/NonPromptJPsi/");
468+
fRegistry.addClone("Generated/VM/All/", "Generated/VM/Acc/");
469+
464470
// reconstructed pair info
465471
fRegistry.add("Pair/sm/Photon/uls/hs", "rec. dilepton", kTHnSparseD, {axis_mass, axis_pt, axis_y, axis_dphi_ee, axis_deta_ee, axis_cos_theta_pol, axis_phi_pol, axis_quadmom, axis_aco, axis_asym_pt, axis_dphi_e_ee, axis_dca}, true);
466472

@@ -542,6 +548,11 @@ struct DileptonMC {
542548
fRegistry.addClone("Pair/corr_bkg_lh/", "Pair/corr_bkg_hh/");
543549
fRegistry.addClone("Pair/corr_bkg_lh/", "Pair/comb_bkg/");
544550

551+
if (doprocessGen_VM) {
552+
fRegistry.add("Generated/VM/Omega/hPtY", "pT vs. y of #omega", kTH2D, {axis_y_meson, axis_pt_meson}, true); // for pT spectrum
553+
fRegistry.add("Generated/VM/Phi/hPtY", "pT vs. y of #phi", kTH2D, {axis_y_meson, axis_pt_meson}, true); // for pT spectrum
554+
}
555+
545556
if (cfgFillUnfolding) {
546557
// for 2D unfolding
547558
const AxisSpec axis_mass_gen{ConfMllBins, "m_{ll}^{gen} (GeV/c^{2})"};
@@ -880,6 +891,17 @@ struct DileptonMC {
880891
eta = lepton.eta();
881892
}
882893

894+
return isInAcceptance(pt, eta);
895+
896+
// if ((mctrackcuts.min_mcPt < pt && pt < mctrackcuts.max_mcPt) && (mctrackcuts.min_mcEta < eta && eta < mctrackcuts.max_mcEta)) {
897+
// return true;
898+
// } else {
899+
// return false;
900+
// }
901+
}
902+
903+
bool isInAcceptance(const float pt, const float eta)
904+
{
883905
if ((mctrackcuts.min_mcPt < pt && pt < mctrackcuts.max_mcPt) && (mctrackcuts.min_mcEta < eta && eta < mctrackcuts.max_mcEta)) {
884906
return true;
885907
} else {
@@ -1825,7 +1847,6 @@ struct DileptonMC {
18251847
o2::math_utils::bringToPMPi(phiPol);
18261848
float quadmom = (3.f * std::pow(cos_thetaPol, 2) - 1.f) / 2.f;
18271849

1828-
// bool isInAcc = isInAcceptance<isSmeared>(t1) && isInAcceptance<isSmeared>(t2);
18291850
if (!isInAcceptance<isSmeared>(t1) || !isInAcceptance<isSmeared>(t2)) {
18301851
return false;
18311852
}
@@ -1897,7 +1918,7 @@ struct DileptonMC {
18971918
case static_cast<int>(EM_HFeeType::kBCe_Be_SameB):
18981919
fillGenHistograms<18>(sign1, sign2, 0, 0, v12.M(), v12.Pt(), v12.Rapidity(), dphi, deta, cos_thetaPol, phiPol, quadmom, aco, asym, std::fabs(dphi_e_ee), weight); // b2c2l_b2l_sameb
18991920
break;
1900-
case static_cast<int>(EM_HFeeType::kBCe_Be_DiffB): // LS
1921+
case static_cast<int>(EM_HFeeType::kBCe_Be_DiffB):
19011922
fillGenHistograms<19>(sign1, sign2, 0, 0, v12.M(), v12.Pt(), v12.Rapidity(), dphi, deta, cos_thetaPol, phiPol, quadmom, aco, asym, std::fabs(dphi_e_ee), weight); // b2c2l_b2l_diffb
19021923
break;
19031924
default:
@@ -1907,6 +1928,113 @@ struct DileptonMC {
19071928
return false;
19081929
}
19091930

1931+
template <typename TMCParticle, typename TMCParticles>
1932+
bool fillGenParticleAcc(TMCParticle const& mcParticle, TMCParticles const& mcParticles)
1933+
{
1934+
if (!mcParticle.isPhysicalPrimary() && !mcParticle.producedByGenerator()) {
1935+
return false;
1936+
}
1937+
if (mcParticle.daughtersIds().size() < 2) {
1938+
return false;
1939+
}
1940+
1941+
if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDielectron) {
1942+
if (mcParticle.y() < dielectroncuts.cfg_min_pair_y || dielectroncuts.cfg_max_pair_y < mcParticle.y()) {
1943+
return false;
1944+
}
1945+
} else if constexpr (pairtype == o2::aod::pwgem::dilepton::utils::pairutil::DileptonPairType::kDimuon) {
1946+
if (mcParticle.y() < dimuoncuts.cfg_min_pair_y || dimuoncuts.cfg_max_pair_y < mcParticle.y()) {
1947+
return false;
1948+
}
1949+
}
1950+
1951+
int pdg = mcParticle.pdgCode();
1952+
if (std::abs(pdg) == 113 && mcParticle.daughtersIds().size() != 2) { // reject dalitz decay
1953+
return false;
1954+
}
1955+
if (std::abs(pdg) == 223 && mcParticle.daughtersIds().size() != 2) { // reject dalitz decay
1956+
return false;
1957+
}
1958+
if (std::abs(pdg) == 333 && mcParticle.daughtersIds().size() != 2) { // reject dalitz decay
1959+
return false;
1960+
}
1961+
// accept radiative decay of charmonia (ee + multiple gamma).
1962+
1963+
float pt1 = 0.f, eta1 = 0.f, phi1 = 0.f, sign1 = 0.f;
1964+
float pt2 = 0.f, eta2 = 0.f, phi2 = 0.f, sign2 = 0.f;
1965+
std::vector<std::array<float, 4>> vDau;
1966+
vDau.reserve(mcParticle.daughtersIds().size());
1967+
for (const auto& daughterId : mcParticle.daughtersIds()) {
1968+
auto dau = mcParticles.rawIteratorAt(daughterId);
1969+
if (std::abs(dau.pdgCode()) == pdg_lepton) {
1970+
vDau.emplace_back(std::array<float, 4>{dau.pt(), dau.eta(), dau.phi(), dau.pdgCode() > 0 ? -1.f : +1.f});
1971+
}
1972+
}
1973+
if (vDau.size() != 2 || vDau[0][3] * vDau[1][3] > 0.f) { // decay objects must be ULS 2 leptons.
1974+
vDau.clear();
1975+
vDau.shrink_to_fit();
1976+
return false;
1977+
}
1978+
1979+
// LOGF(info, "mcParticle.globalIndex() = %d, mcParticle.pdgCode() = %d, mcParticle.daughtersIds().size() = %d", mcParticle.globalIndex(), mcParticle.pdgCode(), mcParticle.daughtersIds().size());
1980+
1981+
ROOT::Math::PtEtaPhiMVector v1(vDau[0][0], vDau[0][1], vDau[0][2], leptonM1);
1982+
ROOT::Math::PtEtaPhiMVector v2(vDau[1][0], vDau[1][1], vDau[1][2], leptonM2);
1983+
ROOT::Math::PtEtaPhiMVector v12 = v1 + v2;
1984+
1985+
std::array<float, 4> arrP1 = {static_cast<float>(v1.Px()), static_cast<float>(v1.Py()), static_cast<float>(v1.Pz()), leptonM1};
1986+
std::array<float, 4> arrP2 = {static_cast<float>(v2.Px()), static_cast<float>(v2.Py()), static_cast<float>(v2.Pz()), leptonM2};
1987+
float cos_thetaPol = 999, phiPol = 999.f;
1988+
if (cfgPolarizationFrame == 0) {
1989+
o2::aod::pwgem::dilepton::utils::pairutil::getAngleCS(arrP1, arrP2, beamE1, beamE2, beamP1, beamP2, vDau[0][3] > 0 ? 1 : -1, cos_thetaPol, phiPol);
1990+
} else if (cfgPolarizationFrame == 1) {
1991+
o2::aod::pwgem::dilepton::utils::pairutil::getAngleHX(arrP1, arrP2, beamE1, beamE2, beamP1, beamP2, vDau[0][3] > 0 ? 1 : -1, cos_thetaPol, phiPol);
1992+
}
1993+
o2::math_utils::bringToPMPi(phiPol);
1994+
float quadmom = (3.f * std::pow(cos_thetaPol, 2) - 1.f) / 2.f;
1995+
1996+
float weight = 1.f;
1997+
switch (std::abs(mcParticle.pdgCode())) {
1998+
case 113:
1999+
fRegistry.fill(HIST("Generated/VM/All/Rho/hs"), v12.M(), mcParticle.pt(), mcParticle.y(), cos_thetaPol, phiPol, quadmom, weight);
2000+
if (isInAcceptance(v1.Pt(), v1.Eta()) && isInAcceptance(v2.Pt(), v2.Eta())) {
2001+
fRegistry.fill(HIST("Generated/VM/Acc/Rho/hs"), v12.M(), mcParticle.pt(), mcParticle.y(), cos_thetaPol, phiPol, quadmom, weight);
2002+
}
2003+
break;
2004+
case 223:
2005+
fRegistry.fill(HIST("Generated/VM/All/Omega/hs"), v12.M(), mcParticle.pt(), mcParticle.y(), cos_thetaPol, phiPol, quadmom, weight);
2006+
if (isInAcceptance(v1.Pt(), v1.Eta()) && isInAcceptance(v2.Pt(), v2.Eta())) {
2007+
fRegistry.fill(HIST("Generated/VM/Acc/Omega/hs"), v12.M(), mcParticle.pt(), mcParticle.y(), cos_thetaPol, phiPol, quadmom, weight);
2008+
}
2009+
break;
2010+
case 333:
2011+
fRegistry.fill(HIST("Generated/VM/All/Phi/hs"), v12.M(), mcParticle.pt(), mcParticle.y(), cos_thetaPol, phiPol, quadmom, weight);
2012+
if (isInAcceptance(v1.Pt(), v1.Eta()) && isInAcceptance(v2.Pt(), v2.Eta())) {
2013+
fRegistry.fill(HIST("Generated/VM/Acc/Phi/hs"), v12.M(), mcParticle.pt(), mcParticle.y(), cos_thetaPol, phiPol, quadmom, weight);
2014+
}
2015+
break;
2016+
case 443:
2017+
if (IsFromBeauty(mcParticle, mcParticles) > 0) {
2018+
fRegistry.fill(HIST("Generated/VM/All/NonPromptJPsi/hs"), v12.M(), mcParticle.pt(), mcParticle.y(), cos_thetaPol, phiPol, quadmom, weight);
2019+
if (isInAcceptance(v1.Pt(), v1.Eta()) && isInAcceptance(v2.Pt(), v2.Eta())) {
2020+
fRegistry.fill(HIST("Generated/VM/Acc/NonPromptJPsi/hs"), v12.M(), mcParticle.pt(), mcParticle.y(), cos_thetaPol, phiPol, quadmom, weight);
2021+
}
2022+
} else {
2023+
fRegistry.fill(HIST("Generated/VM/All/PromptJPsi/hs"), v12.M(), mcParticle.pt(), mcParticle.y(), cos_thetaPol, phiPol, quadmom, weight);
2024+
if (isInAcceptance(v1.Pt(), v1.Eta()) && isInAcceptance(v2.Pt(), v2.Eta())) {
2025+
fRegistry.fill(HIST("Generated/VM/Acc/PromptJPsi/hs"), v12.M(), mcParticle.pt(), mcParticle.y(), cos_thetaPol, phiPol, quadmom, weight);
2026+
}
2027+
}
2028+
break;
2029+
default:
2030+
break;
2031+
}
2032+
2033+
vDau.clear();
2034+
vDau.shrink_to_fit();
2035+
return false;
2036+
}
2037+
19102038
template <bool isSmeared, typename TCollisions, typename TMCLeptons, typename TPreslice, typename TCut, typename TAllTracks, typename TMCCollisions, typename TMCParticles>
19112039
void runTruePairing(TCollisions const& collisions, TMCLeptons const& posTracks, TMCLeptons const& negTracks, TPreslice const& perCollision, TCut const& cut, TAllTracks const& tracks, TMCCollisions const& mccollisions, TMCParticles const& mcparticles)
19122040
{
@@ -1986,7 +2114,19 @@ struct DileptonMC {
19862114

19872115
for (const auto& [t1, t2] : combinations(CombinationsStrictlyUpperIndexPolicy(negTracks_per_coll, negTracks_per_coll))) { // LS--
19882116
fillGenPairInfo<isSmeared>(t1, t2, mcparticles);
1989-
} // end of true LS++ pair loop
2117+
} // end of true LS-- pair loop
2118+
2119+
// acceptance for polarization of vector mesons
2120+
auto mcParticles_per_coll = mcparticles.sliceBy(perMcCollision, mccollision.globalIndex());
2121+
for (const auto& mcParticle : mcParticles_per_coll) {
2122+
if (!mcParticle.isPhysicalPrimary() && !mcParticle.producedByGenerator()) {
2123+
continue;
2124+
}
2125+
int pdg = std::abs(mcParticle.pdgCode());
2126+
if (pdg == 113 || pdg == 223 || pdg == 333 || pdg == 443) { // select only VMs
2127+
fillGenParticleAcc(mcParticle, mcparticles); // VMs that decay into dilepton are stored in derived data. This is sufficient for polarization.
2128+
}
2129+
} // end of mc particle loop
19902130
} // end of collision loop
19912131
}
19922132

@@ -2462,7 +2602,6 @@ struct DileptonMC {
24622602
auto mctracks_per_coll = mcparticles.sliceBy(perMcCollision_vm, mccollision.globalIndex());
24632603

24642604
for (const auto& mctrack : mctracks_per_coll) {
2465-
24662605
if (!(mctrack.isPhysicalPrimary() || mctrack.producedByGenerator())) {
24672606
continue;
24682607
}
@@ -2479,10 +2618,10 @@ struct DileptonMC {
24792618

24802619
switch (std::abs(mctrack.pdgCode())) {
24812620
case 223:
2482-
fRegistry.fill(HIST("Generated/sm/Omega2ll/uls/hPtY"), mctrack.y(), mctrack.pt(), 1.f / mctrack.dsf());
2621+
fRegistry.fill(HIST("Generated/VM/Omega/hPtY"), mctrack.y(), mctrack.pt(), 1.f / mctrack.dsf());
24832622
break;
24842623
case 333:
2485-
fRegistry.fill(HIST("Generated/sm/Phi2ll/uls/hPtY"), mctrack.y(), mctrack.pt(), 1.f / mctrack.dsf());
2624+
fRegistry.fill(HIST("Generated/VM/Phi/hPtY"), mctrack.y(), mctrack.pt(), 1.f / mctrack.dsf());
24862625
break;
24872626
default:
24882627
break;

0 commit comments

Comments
 (0)