Skip to content

Commit 4ec9cc8

Browse files
alcalivaalcaliva
committed
[PWGLF] Add histograms for efficiency factorization
1 parent 92456b8 commit 4ec9cc8

File tree

1 file changed

+176
-89
lines changed

1 file changed

+176
-89
lines changed

‎PWGLF/Tasks/Strangeness/strangenessInJets.cxx‎

Lines changed: 176 additions & 89 deletions
Original file line numberDiff line numberDiff line change
@@ -302,6 +302,19 @@ struct StrangenessInJets {
302302
registryMC.add("Proton_generated_in_jet", "Proton_generated_in_jet", HistType::kTH2F, {multBinning, ptAxisLongLived});
303303
registryMC.add("Proton_generated_in_ue", "Proton_generated_in_ue", HistType::kTH2F, {multBinning, ptAxisLongLived});
304304
}
305+
306+
// Histograms to calculate probability of hyperons to be found within jets
307+
if (enabledSignals.value[ParticleOfInterest::kV0Particles]) {
308+
registryMC.add("K0s_generated_fullevent", "K0s_generated_fullevent", HistType::kTH2F, {multBinning, ptAxis});
309+
registryMC.add("Lambda_generated_fullevent", "Lambda_generated_fullevent", HistType::kTH2F, {multBinning, ptAxis});
310+
registryMC.add("AntiLambda_generated_fullevent", "AntiLambda_generated_fullevent", HistType::kTH2F, {multBinning, ptAxis});
311+
}
312+
if (enabledSignals.value[ParticleOfInterest::kCascades]) {
313+
registryMC.add("XiPos_generated_fullevent", "XiPos_generated_fullevent", HistType::kTH2F, {multBinning, ptAxis});
314+
registryMC.add("XiNeg_generated_fullevent", "XiNeg_generated_fullevent", HistType::kTH2F, {multBinning, ptAxis});
315+
registryMC.add("OmegaPos_generated_fullevent", "OmegaPos_generated_fullevent", HistType::kTH2F, {multBinning, ptAxis});
316+
registryMC.add("OmegaNeg_generated_fullevent", "OmegaNeg_generated_fullevent", HistType::kTH2F, {multBinning, ptAxis});
317+
}
305318
}
306319

307320
// Histograms for mc reconstructed
@@ -350,58 +363,21 @@ struct StrangenessInJets {
350363
registryMC.add("Proton_reconstructed_in_jet", "Proton_reconstructed_in_jet", HistType::kTH2F, {multBinning, ptAxisLongLived});
351364
registryMC.add("Proton_reconstructed_in_ue", "Proton_reconstructed_in_ue", HistType::kTH2F, {multBinning, ptAxisLongLived});
352365
}
353-
}
354-
}
355-
356-
/*
357-
// Calculation of perpendicular axes
358-
void getPerpendicularAxis(TVector3 p, TVector3& u, double sign)
359-
{
360-
// initialization
361-
double ux(0), uy(0), uz(0);
362-
363-
// components of vector p
364-
const double px = p.X();
365-
const double py = p.Y();
366-
const double pz = p.Z();
367-
368-
// protection 1
369-
if (px == 0 && py != 0) {
370-
uy = -(pz * pz) / py;
371-
ux = sign * std::sqrt(py * py - (pz * pz * pz * pz) / (py * py));
372-
uz = pz;
373-
u.SetXYZ(ux, uy, uz);
374-
return;
375-
}
376366

377-
// protection 2
378-
if (py == 0 && px != 0) {
379-
ux = -(pz * pz) / px;
380-
uy = sign * std::sqrt(px * px - (pz * pz * pz * pz) / (px * px));
381-
uz = pz;
382-
u.SetXYZ(ux, uy, uz);
383-
return;
384-
}
385-
386-
// equation parameters
387-
const double a = px * px + py * py;
388-
const double b = 2.0 * px * pz * pz;
389-
const double c = pz * pz * pz * pz - py * py * py * py - px * px * py * py;
390-
const double delta = b * b - 4.0 * a * c;
391-
392-
// protection agains delta<0
393-
if (delta < 0) {
394-
return;
367+
// Histograms to calculate probability of hyperons to be found within jets
368+
if (enabledSignals.value[ParticleOfInterest::kV0Particles]) {
369+
registryMC.add("K0s_reconstructed_fullevent", "K0s_reconstructed_fullevent", HistType::kTH2F, {multBinning, ptAxis});
370+
registryMC.add("Lambda_reconstructed_fullevent", "Lambda_reconstructed_fullevent", HistType::kTH2F, {multBinning, ptAxis});
371+
registryMC.add("AntiLambda_reconstructed_fullevent", "AntiLambda_reconstructed_fullevent", HistType::kTH2F, {multBinning, ptAxis});
372+
}
373+
if (enabledSignals.value[ParticleOfInterest::kCascades]) {
374+
registryMC.add("XiPos_reconstructed_fullevent", "XiPos_reconstructed_fullevent", HistType::kTH2F, {multBinning, ptAxis});
375+
registryMC.add("XiNeg_reconstructed_fullevent", "XiNeg_reconstructed_fullevent", HistType::kTH2F, {multBinning, ptAxis});
376+
registryMC.add("OmegaPos_reconstructed_fullevent", "OmegaPos_reconstructed_fullevent", HistType::kTH2F, {multBinning, ptAxis});
377+
registryMC.add("OmegaNeg_reconstructed_fullevent", "OmegaNeg_reconstructed_fullevent", HistType::kTH2F, {multBinning, ptAxis});
378+
}
395379
}
396-
397-
// solutions
398-
ux = (-b + sign * std::sqrt(delta)) / (2.0 * a);
399-
uy = (-pz * pz - px * ux) / py;
400-
uz = pz;
401-
u.SetXYZ(ux, uy, uz);
402-
return;
403380
}
404-
*/
405381

406382
// Delta phi calculation
407383
double getDeltaPhi(double a1, double a2)
@@ -1246,14 +1222,58 @@ struct StrangenessInJets {
12461222
strHadronMomentum.emplace_back(particle.px(), particle.py(), particle.pz());
12471223
}
12481224

1249-
// Select physical primary particles or HF decay products
1250-
if (!isPhysicalPrimaryOrFromHF(particle, mcParticles))
1251-
continue;
1252-
1225+
// Eta and pt selection
12531226
double minPtParticle = 0.1;
12541227
if (particle.eta() < etaMin || particle.eta() > etaMax || particle.pt() < minPtParticle)
12551228
continue;
12561229

1230+
// Fill generated histograms for full event
1231+
if (particle.isPhysicalPrimary()) {
1232+
switch (particle.pdgCode()) {
1233+
case kK0Short:
1234+
if (enabledSignals.value[ParticleOfInterest::kV0Particles]) {
1235+
registryMC.fill(HIST("K0s_generated_fullevent"), genMultiplicity, particle.pt());
1236+
}
1237+
break;
1238+
case kLambda0:
1239+
if (enabledSignals.value[ParticleOfInterest::kV0Particles]) {
1240+
registryMC.fill(HIST("Lambda_generated_fullevent"), genMultiplicity, particle.pt());
1241+
}
1242+
break;
1243+
case kLambda0Bar:
1244+
if (enabledSignals.value[ParticleOfInterest::kV0Particles]) {
1245+
registryMC.fill(HIST("AntiLambda_generated_fullevent"), genMultiplicity, particle.pt());
1246+
}
1247+
break;
1248+
case kXiMinus:
1249+
if (enabledSignals.value[ParticleOfInterest::kCascades]) {
1250+
registryMC.fill(HIST("XiNeg_generated_fullevent"), genMultiplicity, particle.pt());
1251+
}
1252+
break;
1253+
case kXiPlusBar:
1254+
if (enabledSignals.value[ParticleOfInterest::kCascades]) {
1255+
registryMC.fill(HIST("XiPos_generated_fullevent"), genMultiplicity, particle.pt());
1256+
}
1257+
break;
1258+
case kOmegaMinus:
1259+
if (enabledSignals.value[ParticleOfInterest::kCascades]) {
1260+
registryMC.fill(HIST("OmegaNeg_generated_fullevent"), genMultiplicity, particle.pt());
1261+
}
1262+
break;
1263+
case kOmegaPlusBar:
1264+
if (enabledSignals.value[ParticleOfInterest::kCascades]) {
1265+
registryMC.fill(HIST("OmegaPos_generated_fullevent"), genMultiplicity, particle.pt());
1266+
}
1267+
break;
1268+
default:
1269+
break;
1270+
}
1271+
}
1272+
1273+
// Select physical primary particles or HF decay products
1274+
if (!isPhysicalPrimaryOrFromHF(particle, mcParticles))
1275+
continue;
1276+
12571277
// Build 4-momentum assuming charged pion mass
12581278
static constexpr float MassPionChargedSquared = o2::constants::physics::MassPionCharged * o2::constants::physics::MassPionCharged;
12591279
const double energy = std::sqrt(particle.p() * particle.p() + MassPionChargedSquared);
@@ -1472,7 +1492,7 @@ struct StrangenessInJets {
14721492
// Reconstructed MC events
14731493
void processMCreconstructed(SimCollisions const& collisions, soa::Join<aod::McCollisions, aod::McCentFT0Ms> const&,
14741494
DaughterTracksMC const& mcTracks, aod::V0Datas const& fullV0s,
1475-
aod::CascDataExt const& Cascades, const aod::McParticles&)
1495+
aod::CascDataExt const& Cascades, aod::McParticles const& mcParticles)
14761496
{
14771497
// Define per-event containers
14781498
std::vector<fastjet::PseudoJet> fjParticles;
@@ -1520,6 +1540,89 @@ struct StrangenessInJets {
15201540
auto cascPerColl = Cascades.sliceBy(perCollisionCasc, collision.globalIndex());
15211541
auto tracksPerColl = mcTracks.sliceBy(perCollisionTrk, collision.globalIndex());
15221542

1543+
// V0 particles
1544+
if (enabledSignals.value[ParticleOfInterest::kV0Particles]) {
1545+
for (const auto& v0 : v0sPerColl) {
1546+
const auto& pos = v0.posTrack_as<DaughterTracksMC>();
1547+
const auto& neg = v0.negTrack_as<DaughterTracksMC>();
1548+
1549+
// Get MC particles
1550+
if (!pos.has_mcParticle() || !neg.has_mcParticle())
1551+
continue;
1552+
auto posParticle = pos.mcParticle_as<aod::McParticles>();
1553+
auto negParticle = neg.mcParticle_as<aod::McParticles>();
1554+
if (!posParticle.has_mothers() || !negParticle.has_mothers())
1555+
continue;
1556+
1557+
// Get Mothers
1558+
auto motherPos = mcParticles.iteratorAt(posParticle.mothersIds()[0]);
1559+
auto motherNeg = mcParticles.iteratorAt(negParticle.mothersIds()[0]);
1560+
if (motherPos != motherNeg)
1561+
continue;
1562+
if (!motherPos.isPhysicalPrimary())
1563+
continue;
1564+
1565+
// K0s
1566+
if (passedK0ShortSelection(v0, pos, neg) && motherPos.pdgCode() == kK0Short) {
1567+
registryMC.fill(HIST("K0s_reconstructed_fullevent"), multiplicity, v0.pt());
1568+
}
1569+
// Lambda
1570+
if (passedLambdaSelection(v0, pos, neg) && motherPos.pdgCode() == kLambda0) {
1571+
registryMC.fill(HIST("Lambda_reconstructed_fullevent"), multiplicity, v0.pt());
1572+
}
1573+
// AntiLambda
1574+
if (passedAntiLambdaSelection(v0, pos, neg) && motherPos.pdgCode() == kLambda0Bar) {
1575+
registryMC.fill(HIST("AntiLambda_reconstructed_fullevent"), multiplicity, v0.pt());
1576+
}
1577+
}
1578+
}
1579+
1580+
// Cascades
1581+
if (enabledSignals.value[ParticleOfInterest::kCascades]) {
1582+
for (const auto& casc : cascPerColl) {
1583+
auto bach = casc.bachelor_as<DaughterTracksMC>();
1584+
auto pos = casc.posTrack_as<DaughterTracksMC>();
1585+
auto neg = casc.negTrack_as<DaughterTracksMC>();
1586+
1587+
// Get MC particles
1588+
if (!bach.has_mcParticle() || !pos.has_mcParticle() || !neg.has_mcParticle())
1589+
continue;
1590+
auto posParticle = pos.mcParticle_as<aod::McParticles>();
1591+
auto negParticle = neg.mcParticle_as<aod::McParticles>();
1592+
auto bachParticle = bach.mcParticle_as<aod::McParticles>();
1593+
if (!posParticle.has_mothers() || !negParticle.has_mothers() || !bachParticle.has_mothers())
1594+
continue;
1595+
1596+
// Select particles originating from the same parent
1597+
auto motherPos = mcParticles.iteratorAt(posParticle.mothersIds()[0]);
1598+
auto motherNeg = mcParticles.iteratorAt(negParticle.mothersIds()[0]);
1599+
auto motherBach = mcParticles.iteratorAt(bachParticle.mothersIds()[0]);
1600+
if (motherPos != motherNeg)
1601+
continue;
1602+
if (std::abs(motherPos.pdgCode()) != kLambda0)
1603+
continue;
1604+
if (!motherBach.isPhysicalPrimary())
1605+
continue;
1606+
1607+
// Xi+
1608+
if (passedXiSelection(casc, pos, neg, bach, collision) && bach.sign() > 0 && motherBach.pdgCode() == kXiPlusBar) {
1609+
registryMC.fill(HIST("XiPos_reconstructed_fullevent"), multiplicity, casc.pt());
1610+
}
1611+
// Xi-
1612+
if (passedXiSelection(casc, pos, neg, bach, collision) && bach.sign() < 0 && motherBach.pdgCode() == kXiMinus) {
1613+
registryMC.fill(HIST("XiNeg_reconstructed_fullevent"), multiplicity, casc.pt());
1614+
}
1615+
// Omega+
1616+
if (passedOmegaSelection(casc, pos, neg, bach, collision) && bach.sign() > 0 && motherBach.pdgCode() == kOmegaPlusBar) {
1617+
registryMC.fill(HIST("OmegaPos_reconstructed_fullevent"), multiplicity, casc.pt());
1618+
}
1619+
// Omega-
1620+
if (passedOmegaSelection(casc, pos, neg, bach, collision) && bach.sign() < 0 && motherBach.pdgCode() == kOmegaMinus) {
1621+
registryMC.fill(HIST("OmegaNeg_reconstructed_fullevent"), multiplicity, casc.pt());
1622+
}
1623+
}
1624+
}
1625+
15231626
// Loop over reconstructed tracks
15241627
for (auto const& track : tracksPerColl) {
15251628
if (!passedTrackSelectionForJetReconstruction(track))
@@ -1596,19 +1699,12 @@ struct StrangenessInJets {
15961699
continue;
15971700

15981701
// Select particles originating from the same parent
1599-
int pdgParent(0);
1600-
bool isPhysPrim = false;
1601-
for (const auto& particleMotherOfNeg : negParticle.mothers_as<aod::McParticles>()) {
1602-
for (const auto& particleMotherOfPos : posParticle.mothers_as<aod::McParticles>()) {
1603-
if (particleMotherOfNeg == particleMotherOfPos) {
1604-
pdgParent = particleMotherOfNeg.pdgCode();
1605-
isPhysPrim = particleMotherOfNeg.isPhysicalPrimary();
1606-
}
1607-
}
1608-
}
1609-
if (pdgParent == 0)
1702+
auto motherPos = mcParticles.iteratorAt(posParticle.mothersIds()[0]);
1703+
auto motherNeg = mcParticles.iteratorAt(negParticle.mothersIds()[0]);
1704+
if (motherPos != motherNeg)
16101705
continue;
1611-
1706+
bool isPhysPrim = motherPos.isPhysicalPrimary());
1707+
16121708
// Compute distance from jet and UE axes
16131709
double deltaEtaJet = v0dir.Eta() - selectedJet[i].Eta();
16141710
double deltaPhiJet = getDeltaPhi(v0dir.Phi(), selectedJet[i].Phi());
@@ -1621,7 +1717,7 @@ struct StrangenessInJets {
16211717
double deltaRue2 = std::sqrt(deltaEtaUe2 * deltaEtaUe2 + deltaPhiUe2 * deltaPhiUe2);
16221718

16231719
// K0s
1624-
if (passedK0ShortSelection(v0, pos, neg) && pdgParent == kK0Short && isPhysPrim) {
1720+
if (passedK0ShortSelection(v0, pos, neg) && motherPos.pdgCode() == kK0Short && isPhysPrim) {
16251721
if (deltaRjet < rJet) {
16261722
registryMC.fill(HIST("K0s_reconstructed_jet"), multiplicity, v0.pt());
16271723
}
@@ -1630,7 +1726,7 @@ struct StrangenessInJets {
16301726
}
16311727
}
16321728
// Lambda
1633-
if (passedLambdaSelection(v0, pos, neg) && pdgParent == kLambda0 && isPhysPrim) {
1729+
if (passedLambdaSelection(v0, pos, neg) && motherPos.pdgCode() == kLambda0 && isPhysPrim) {
16341730
if (deltaRjet < rJet) {
16351731
registryMC.fill(HIST("Lambda_reconstructed_jet"), multiplicity, v0.pt());
16361732
}
@@ -1639,7 +1735,7 @@ struct StrangenessInJets {
16391735
}
16401736
}
16411737
// AntiLambda
1642-
if (passedAntiLambdaSelection(v0, pos, neg) && pdgParent == kLambda0Bar && isPhysPrim) {
1738+
if (passedAntiLambdaSelection(v0, pos, neg) && motherPos.pdgCode() == kLambda0Bar && isPhysPrim) {
16431739
if (deltaRjet < rJet) {
16441740
registryMC.fill(HIST("AntiLambda_reconstructed_jet"), multiplicity, v0.pt());
16451741
}
@@ -1696,23 +1792,14 @@ struct StrangenessInJets {
16961792
continue;
16971793

16981794
// Select particles originating from the same parent
1699-
int pdgParent(0);
1700-
bool isPhysPrim = false;
1701-
for (const auto& particleMotherOfNeg : negParticle.mothers_as<aod::McParticles>()) {
1702-
for (const auto& particleMotherOfPos : posParticle.mothers_as<aod::McParticles>()) {
1703-
for (const auto& particleMotherOfBach : bachParticle.mothers_as<aod::McParticles>()) {
1704-
if (particleMotherOfNeg != particleMotherOfPos)
1705-
continue;
1706-
if (std::abs(particleMotherOfNeg.pdgCode()) != kLambda0)
1707-
continue;
1708-
isPhysPrim = particleMotherOfBach.isPhysicalPrimary();
1709-
pdgParent = particleMotherOfBach.pdgCode();
1710-
}
1711-
}
1712-
}
1713-
if (pdgParent == 0)
1795+
auto motherPos = mcParticles.iteratorAt(posParticle.mothersIds()[0]);
1796+
auto motherNeg = mcParticles.iteratorAt(negParticle.mothersIds()[0]);
1797+
auto motherBach = mcParticles.iteratorAt(bachParticle.mothersIds()[0]);
1798+
if (motherPos != motherNeg)
1799+
continue;
1800+
if (std::abs(motherPos.pdgCode()) != kLambda0)
17141801
continue;
1715-
if (!isPhysPrim)
1802+
if (!motherBach.isPhysicalPrimary())
17161803
continue;
17171804

17181805
// Compute distances from jet and UE axes
@@ -1728,7 +1815,7 @@ struct StrangenessInJets {
17281815
double deltaRue2 = std::sqrt(deltaEtaUe2 * deltaEtaUe2 + deltaPhiUe2 * deltaPhiUe2);
17291816

17301817
// Xi+
1731-
if (passedXiSelection(casc, pos, neg, bach, collision) && bach.sign() > 0 && pdgParent == kXiPlusBar) {
1818+
if (passedXiSelection(casc, pos, neg, bach, collision) && bach.sign() > 0 && motherBach.pdgCode() == kXiPlusBar) {
17321819
if (deltaRjet < rJet) {
17331820
registryMC.fill(HIST("XiPos_reconstructed_jet"), multiplicity, casc.pt());
17341821
}
@@ -1737,7 +1824,7 @@ struct StrangenessInJets {
17371824
}
17381825
}
17391826
// Xi-
1740-
if (passedXiSelection(casc, pos, neg, bach, collision) && bach.sign() < 0 && pdgParent == kXiMinus) {
1827+
if (passedXiSelection(casc, pos, neg, bach, collision) && bach.sign() < 0 && motherBach.pdgCode() == kXiMinus) {
17411828
if (deltaRjet < rJet) {
17421829
registryMC.fill(HIST("XiNeg_reconstructed_jet"), multiplicity, casc.pt());
17431830
}
@@ -1746,7 +1833,7 @@ struct StrangenessInJets {
17461833
}
17471834
}
17481835
// Omega+
1749-
if (passedOmegaSelection(casc, pos, neg, bach, collision) && bach.sign() > 0 && pdgParent == kOmegaPlusBar) {
1836+
if (passedOmegaSelection(casc, pos, neg, bach, collision) && bach.sign() > 0 && motherBach.pdgCode() == kOmegaPlusBar) {
17501837
if (deltaRjet < rJet) {
17511838
registryMC.fill(HIST("OmegaPos_reconstructed_jet"), multiplicity, casc.pt());
17521839
}
@@ -1755,7 +1842,7 @@ struct StrangenessInJets {
17551842
}
17561843
}
17571844
// Omega-
1758-
if (passedOmegaSelection(casc, pos, neg, bach, collision) && bach.sign() < 0 && pdgParent == kOmegaMinus) {
1845+
if (passedOmegaSelection(casc, pos, neg, bach, collision) && bach.sign() < 0 && motherBach.pdgCode() == kOmegaMinus) {
17591846
if (deltaRjet < rJet) {
17601847
registryMC.fill(HIST("OmegaNeg_reconstructed_jet"), multiplicity, casc.pt());
17611848
}

0 commit comments

Comments
 (0)