Skip to content

Commit 1592596

Browse files
abmodakabmodak
authored
[PWGLF] Add updated track efficiency estimation (#13820)
1 parent bfc91e5 commit 1592596

File tree

1 file changed

+173
-9
lines changed

1 file changed

+173
-9
lines changed

‎PWGLF/Tasks/GlobalEventProperties/heavyionMultiplicity.cxx‎

Lines changed: 173 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -92,6 +92,31 @@ enum {
9292
kSpeciesend
9393
};
9494

95+
enum {
96+
kGenTrkTypebegin = 0,
97+
kGenAll = 1,
98+
kGenPion,
99+
kGenKaon,
100+
kGenProton,
101+
kGenOther,
102+
kGenTrkTypeend
103+
};
104+
105+
enum {
106+
kRecTrkTypebegin = 0,
107+
kRecoAll = 1,
108+
kRecoPrimary,
109+
kRecoPion,
110+
kRecoKaon,
111+
kRecoProton,
112+
kRecoOther,
113+
kRecoSecondary,
114+
kRecoWeakDecay,
115+
kRecoFake,
116+
kRecoBkg,
117+
kRecTrkTypeend
118+
};
119+
95120
static constexpr TrackSelectionFlags::flagtype TrackSelectionIts =
96121
TrackSelectionFlags::kITSNCls | TrackSelectionFlags::kITSChi2NDF |
97122
TrackSelectionFlags::kITSHits;
@@ -113,12 +138,15 @@ AxisSpec axisCent{100, 0, 100, "#Cent"};
113138
AxisSpec axisTrackType = {kTrackTypeend - 1, +kTrackTypebegin + 0.5, +kTrackTypeend - 0.5, "", "TrackTypeAxis"};
114139
AxisSpec axisGenPtVary = {kGenpTend - 1, +kGenpTbegin + 0.5, +kGenpTend - 0.5, "", "GenpTVaryAxis"};
115140
AxisSpec axisSpecies = {kSpeciesend - 1, +kSpeciesbegin + 0.5, +kSpeciesend - 0.5, "", "SpeciesAxis"};
141+
AxisSpec axisGenTrkType = {kGenTrkTypeend - 1, +kGenTrkTypebegin + 0.5, +kGenTrkTypeend - 0.5, "", "GenTrackTypeAxis"};
142+
AxisSpec axisRecTrkType = {kRecTrkTypeend - 1, +kRecTrkTypebegin + 0.5, +kRecTrkTypeend - 0.5, "", "RecTrackTypeAxis"};
116143
AxisSpec axisMassK0s = {200, 0.4, 0.6, "K0sMass", "K0sMass"};
117144
AxisSpec axisMassLambda = {200, 1.07, 1.17, "Lambda/AntiLamda Mass", "Lambda/AntiLamda Mass"};
118145
AxisSpec axisTracks{9, 0.5, 9.5, "#tracks", "TrackAxis"};
119-
auto static constexpr kMinCharge = 3.f;
120-
auto static constexpr kMinpTcut = 0.1f;
121-
auto static constexpr kNItslayers = 7;
146+
AxisSpec axisDeltaEta{50, -1.0, +1.0, "#Delta(#eta)"};
147+
auto static constexpr KminCharge = 3.f;
148+
auto static constexpr KminPtCut = 0.1f;
149+
auto static constexpr KnItsLayers = 7;
122150

123151
struct HeavyionMultiplicity {
124152

@@ -229,7 +257,7 @@ struct HeavyionMultiplicity {
229257
auto* x2 = htrack->GetAxis(1);
230258
x2->SetBinLabel(1, "All tracks");
231259
x2->SetBinLabel(2, "Non-fake tracks");
232-
for (int i = 0; i < kNItslayers; i++) {
260+
for (int i = 0; i < KnItsLayers; i++) {
233261
x2->SetBinLabel(i + 3, Form("layer %d", i));
234262
}
235263
}
@@ -293,6 +321,22 @@ struct HeavyionMultiplicity {
293321
histos.add("hgendndetaVscentBeforeEvtSel", "hgendndetaBeforeEvtSel vs centrality", kTH2F, {axisEta, impactParAxis});
294322
histos.add("hgendndetaVscentAfterEvtSel", "hgendndetaAfterEvtSel vs centrality", kTH2F, {axisEta, impactParAxis});
295323
}
324+
325+
if (doprocessMCeff) {
326+
histos.add("hGenMCvertexZ", "hGenMCvertexZ", kTH1D, {axisVtxZ}, false);
327+
histos.add("hGenMCvtxzcent", "hGenMCvtxzcent", kTH3D, {axisVtxZ, centAxis, axisOccupancy}, false);
328+
histos.add("hGenMCAssoRecvertexZ", "hGenMCAssoRecvertexZ", kTH1D, {axisVtxZ}, false);
329+
histos.add("hGenMCAssoRecvtxzcent", "hGenMCAssoRecvtxzcent", kTH3D, {axisVtxZ, centAxis, axisOccupancy}, false);
330+
histos.add("hGenMCdndeta", "hGenMCdndeta", kTHnSparseD, {axisVtxZ, centAxis, axisOccupancy, axisEta, axisPhi}, false);
331+
histos.add("hGenMCAssoRecdndeta", "hGenMCAssoRecdndeta", kTHnSparseD, {axisVtxZ, centAxis, axisOccupancy, axisEta, axisPhi, axisGenTrkType, axisGenPtVary}, false);
332+
333+
histos.add("hRecMCvertexZ", "hRecMCvertexZ", kTH1D, {axisVtxZ}, false);
334+
histos.add("hRecMCvtxzcent", "hRecMCvtxzcent", kTH3D, {axisVtxZ, centAxis, axisOccupancy}, false);
335+
histos.add("hRecMCcentrality", "hRecMCcentrality", kTH1D, {axisCent}, false);
336+
histos.add("hRecMCphivseta", "hRecMCphivseta", kTH2D, {axisPhi2, axisEta}, false);
337+
histos.add("hRecMCdndeta", "hRecMCdndeta", kTHnSparseD, {axisVtxZ, centAxis, axisOccupancy, axisEta, axisPhi, axisRecTrkType}, false);
338+
histos.add("etaResolution", "etaResolution", kTH2D, {axisEta, axisDeltaEta});
339+
}
296340
}
297341

298342
template <typename CheckCol>
@@ -402,7 +446,7 @@ struct HeavyionMultiplicity {
402446
if (pdgTrack == nullptr) {
403447
return false;
404448
}
405-
if (std::abs(pdgTrack->Charge()) < kMinCharge) {
449+
if (std::abs(pdgTrack->Charge()) < KminCharge) {
406450
return false;
407451
}
408452
if (std::abs(track.eta()) >= etaRange) {
@@ -543,7 +587,7 @@ struct HeavyionMultiplicity {
543587
continue;
544588
}
545589
histos.fill(HIST("hmcgendndeta"), RecCol.posZ(), selColCent(RecCol), particle.eta(), particle.phi(), static_cast<double>(kSpAll), kNoGenpTVar);
546-
if (particle.pt() < kMinpTcut) {
590+
if (particle.pt() < KminPtCut) {
547591
histos.fill(HIST("hmcgendndeta"), RecCol.posZ(), selColCent(RecCol), particle.eta(), particle.phi(), static_cast<double>(kSpAll), kGenpTup, -10.0 * particle.pt() + 2);
548592
histos.fill(HIST("hmcgendndeta"), RecCol.posZ(), selColCent(RecCol), particle.eta(), particle.phi(), static_cast<double>(kSpAll), kGenpTdown, 5.0 * particle.pt() + 0.5);
549593
} else {
@@ -607,7 +651,7 @@ struct HeavyionMultiplicity {
607651
continue;
608652
}
609653
histos.fill(HIST("hmcgendndpt"), selColCent(RecCol), particle.pt(), kNoGenpTVar);
610-
if (particle.pt() < kMinpTcut) {
654+
if (particle.pt() < KminPtCut) {
611655
histos.fill(HIST("hmcgendndpt"), selColCent(RecCol), particle.pt(), kGenpTup, -10.0 * particle.pt() + 2);
612656
histos.fill(HIST("hmcgendndpt"), selColCent(RecCol), particle.pt(), kGenpTdown, 5.0 * particle.pt() + 0.5);
613657
} else {
@@ -641,7 +685,7 @@ struct HeavyionMultiplicity {
641685
}
642686
histos.fill(HIST("hTracksCount"), selColCent(RecCol), 1);
643687
bool isFakeItsTracks = false;
644-
for (int i = 0; i < kNItslayers; i++) {
688+
for (int i = 0; i < KnItsLayers; i++) {
645689
if (Rectrack.mcMask() & 1 << i) {
646690
isFakeItsTracks = true;
647691
histos.fill(HIST("hTracksCount"), selColCent(RecCol), i + 3);
@@ -793,7 +837,7 @@ struct HeavyionMultiplicity {
793837
continue;
794838
}
795839
histos.fill(HIST("hmcgendndetapp"), RecCol.posZ(), RecCol.centFT0M(), particle.eta(), particle.phi(), static_cast<double>(kSpAll), kNoGenpTVar);
796-
if (particle.pt() < kMinpTcut) {
840+
if (particle.pt() < KminPtCut) {
797841
histos.fill(HIST("hmcgendndetapp"), RecCol.posZ(), RecCol.centFT0M(), particle.eta(), particle.phi(), static_cast<double>(kSpAll), kGenpTup, -10.0 * particle.pt() + 2);
798842
histos.fill(HIST("hmcgendndetapp"), RecCol.posZ(), RecCol.centFT0M(), particle.eta(), particle.phi(), static_cast<double>(kSpAll), kGenpTdown, 5.0 * particle.pt() + 0.5);
799843
} else {
@@ -916,6 +960,125 @@ struct HeavyionMultiplicity {
916960
}
917961
}
918962

963+
void processMCeff(soa::Join<aod::McCollisions, aod::McCollsExtra>::iterator const& mcCollision, CollisionMCRecTable const& RecCols, TrackMCTrueTable const& GenParticles, FilTrackMCRecTable const& RecTracks)
964+
{
965+
auto gencent = -999;
966+
auto genoccu = -999;
967+
bool atLeastOne = false;
968+
969+
for (const auto& RecCol : RecCols) {
970+
if (!isEventSelected(RecCol)) {
971+
continue;
972+
}
973+
if (RecCol.globalIndex() != mcCollision.bestCollisionIndex()) {
974+
continue;
975+
}
976+
atLeastOne = true;
977+
gencent = selColCent(RecCol);
978+
genoccu = selColOccu(RecCol);
979+
}
980+
981+
histos.fill(HIST("hGenMCvertexZ"), mcCollision.posZ());
982+
histos.fill(HIST("hGenMCvtxzcent"), mcCollision.posZ(), gencent, genoccu);
983+
984+
if (atLeastOne) {
985+
histos.fill(HIST("hGenMCAssoRecvertexZ"), mcCollision.posZ());
986+
histos.fill(HIST("hGenMCAssoRecvtxzcent"), mcCollision.posZ(), gencent, genoccu);
987+
}
988+
989+
for (const auto& particle : GenParticles) {
990+
if (!isGenTrackSelected(particle)) {
991+
continue;
992+
}
993+
histos.fill(HIST("hGenMCdndeta"), mcCollision.posZ(), gencent, genoccu, particle.eta(), particle.phi());
994+
if (atLeastOne) {
995+
histos.fill(HIST("hGenMCAssoRecdndeta"), mcCollision.posZ(), gencent, genoccu, particle.eta(), particle.phi(), static_cast<double>(kGenAll), kNoGenpTVar);
996+
if (particle.pt() < KminPtCut) {
997+
histos.fill(HIST("hGenMCAssoRecdndeta"), mcCollision.posZ(), gencent, genoccu, particle.eta(), particle.phi(), static_cast<double>(kGenAll), kGenpTup, -10.0 * particle.pt() + 2);
998+
histos.fill(HIST("hGenMCAssoRecdndeta"), mcCollision.posZ(), gencent, genoccu, particle.eta(), particle.phi(), static_cast<double>(kGenAll), kGenpTdown, 5.0 * particle.pt() + 0.5);
999+
} else {
1000+
histos.fill(HIST("hGenMCAssoRecdndeta"), mcCollision.posZ(), gencent, genoccu, particle.eta(), particle.phi(), static_cast<double>(kSpAll), kGenpTup);
1001+
histos.fill(HIST("hGenMCAssoRecdndeta"), mcCollision.posZ(), gencent, genoccu, particle.eta(), particle.phi(), static_cast<double>(kGenAll), kGenpTdown);
1002+
}
1003+
int pid = 0;
1004+
switch (std::abs(particle.pdgCode())) {
1005+
case PDG_t::kPiPlus:
1006+
pid = kGenPion;
1007+
break;
1008+
case PDG_t::kKPlus:
1009+
pid = kGenKaon;
1010+
break;
1011+
case PDG_t::kProton:
1012+
pid = kGenProton;
1013+
break;
1014+
default:
1015+
pid = kGenOther;
1016+
break;
1017+
}
1018+
histos.fill(HIST("hGenMCAssoRecdndeta"), mcCollision.posZ(), gencent, genoccu, particle.eta(), particle.phi(), static_cast<double>(pid), kNoGenpTVar);
1019+
} // Associated with reco col
1020+
} // track (mcgen) loop
1021+
1022+
for (const auto& RecCol : RecCols) {
1023+
if (!isEventSelected(RecCol)) {
1024+
continue;
1025+
}
1026+
if (RecCol.globalIndex() != mcCollision.bestCollisionIndex()) {
1027+
continue;
1028+
}
1029+
histos.fill(HIST("hRecMCvertexZ"), RecCol.posZ());
1030+
histos.fill(HIST("hRecMCcentrality"), selColCent(RecCol));
1031+
histos.fill(HIST("hRecMCvtxzcent"), RecCol.posZ(), selColCent(RecCol), selColOccu(RecCol));
1032+
1033+
auto recTracksPart = RecTracks.sliceBy(perCollision, RecCol.globalIndex());
1034+
std::vector<int> mclabels;
1035+
for (const auto& Rectrack : recTracksPart) {
1036+
if (!isTrackSelected(Rectrack)) {
1037+
continue;
1038+
}
1039+
histos.fill(HIST("hRecMCphivseta"), Rectrack.phi(), Rectrack.eta());
1040+
histos.fill(HIST("hRecMCdndeta"), RecCol.posZ(), selColCent(RecCol), selColOccu(RecCol), Rectrack.eta(), Rectrack.phi(), static_cast<double>(kRecoAll));
1041+
if (Rectrack.has_mcParticle()) {
1042+
int pid = 0;
1043+
auto mcpart = Rectrack.mcParticle();
1044+
histos.fill(HIST("etaResolution"), Rectrack.eta(), Rectrack.eta() - mcpart.eta());
1045+
if (mcpart.isPhysicalPrimary()) {
1046+
pid = kRecoPrimary;
1047+
switch (std::abs(mcpart.pdgCode())) {
1048+
case PDG_t::kPiPlus:
1049+
pid = kRecoPion;
1050+
break;
1051+
case PDG_t::kKPlus:
1052+
pid = kRecoKaon;
1053+
break;
1054+
case PDG_t::kProton:
1055+
pid = kRecoProton;
1056+
break;
1057+
default:
1058+
pid = kRecoOther;
1059+
break;
1060+
}
1061+
} else {
1062+
pid = kRecoSecondary;
1063+
}
1064+
if (mcpart.has_mothers()) {
1065+
auto mcpartMother = mcpart.template mothers_as<aod::McParticles>().front();
1066+
if (mcpartMother.pdgCode() == PDG_t::kK0Short || std::abs(mcpartMother.pdgCode()) == PDG_t::kLambda0) {
1067+
pid = kRecoWeakDecay;
1068+
}
1069+
}
1070+
if (find(mclabels.begin(), mclabels.end(), Rectrack.mcParticleId()) != mclabels.end()) {
1071+
pid = kRecoFake;
1072+
}
1073+
mclabels.push_back(Rectrack.mcParticleId());
1074+
histos.fill(HIST("hRecMCdndeta"), RecCol.posZ(), selColCent(RecCol), selColOccu(RecCol), mcpart.eta(), mcpart.phi(), static_cast<double>(pid));
1075+
} else {
1076+
histos.fill(HIST("hRecMCdndeta"), RecCol.posZ(), selColCent(RecCol), selColOccu(RecCol), Rectrack.eta(), Rectrack.phi(), static_cast<double>(kRecoBkg));
1077+
}
1078+
} // track (mcrec) loop
1079+
} // collision loop
1080+
}
1081+
9191082
PROCESS_SWITCH(HeavyionMultiplicity, processData, "process data CentFT0C", false);
9201083
PROCESS_SWITCH(HeavyionMultiplicity, processCorrelation, "do correlation study in data", false);
9211084
PROCESS_SWITCH(HeavyionMultiplicity, processMonteCarlo, "process MC CentFT0C", false);
@@ -926,6 +1089,7 @@ struct HeavyionMultiplicity {
9261089
PROCESS_SWITCH(HeavyionMultiplicity, processppMonteCarlo, "process pp MC", false);
9271090
PROCESS_SWITCH(HeavyionMultiplicity, processGen, "process pure MC gen", false);
9281091
PROCESS_SWITCH(HeavyionMultiplicity, processEvtLossSigLossMC, "process Signal Loss, Event Loss", false);
1092+
PROCESS_SWITCH(HeavyionMultiplicity, processMCeff, "process extra efficiency function", false);
9291093
};
9301094

9311095
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

0 commit comments

Comments
 (0)