Skip to content

Commit 2dd7aca

Browse files
SuJeong-Jialibuild
SuJeong-Ji
and
alibuild
authored
[PWGLF] Add histogram for centrality distribution for correction factors (#16512)
Co-authored-by: ALICE Action Bot <alibuild@cern.ch>
1 parent 3ea8e55 commit 2dd7aca

1 file changed

Lines changed: 140 additions & 34 deletions

File tree

‎PWGLF/Tasks/Resonances/chk892LI.cxx‎

Lines changed: 140 additions & 34 deletions
Original file line numberDiff line numberDiff line change
@@ -145,6 +145,9 @@ struct Chk892LI {
145145
ConfigurableAxis cfgBinsVtxZ{"cfgBinsVtxZ", {VARIABLE_WIDTH, -10.0, -9.0, -8.0, -7.0, -6.0, -5.0, -4.0, -3.0, -2.0, -1.0, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0}, "Binning of the z-vertex axis"};
146146
Configurable<int> cNbinsDiv{"cNbinsDiv", 1, "Integer to divide the number of bins"};
147147
Configurable<int> cNbinsDivQA{"cNbinsDivQA", 1, "Integer to divide the number of bins for QA"};
148+
149+
Configurable<std::vector<float>> cfgGenMultCuts{"cfgGenMultCuts", std::vector<float>{500, 300, 200, 120, 80, 50, 30, 10, 0}, "Generated multiplicity lower cuts corresponding to reco centrality bins"};
150+
Configurable<std::vector<float>> cfgCentBinCentres{"cfgCentBinCentres", std::vector<float>{2.5, 7.5, 15.0, 25.0, 35.0, 45.0, 55.0, 75.0, 95.0}, "Reco centrality bin centres"};
148151
} AxisConfig;
149152

150153
/// Event cuts
@@ -253,7 +256,10 @@ struct Chk892LI {
253256

254257
Configurable<bool> cfgRecoUseInelGt0{"cfgRecoUseInelGt0", false, "Data/Reco require INEL>0"};
255258
Configurable<bool> cfgTruthUseInelGt0{"cfgTruthUseInelGt0", false, "Truth denominator: require INEL>0"};
256-
Configurable<bool> cfgTruthIncludeZvtx{"cfgTruthIncludeZvtx", true, "Truth denominator: also require |vtxz|<cfgEvtZvtx"};
259+
Configurable<bool> cfgTruthIncludeZvtx{"cfgTruthIncludeZvtx", false, "Truth denominator: also require |vtxz|<cfgEvtZvtx"};
260+
261+
Configurable<float> cfgGenMultEtaMax{"cfgGenMultEtaMax", 0.5, "Max Eta for generated mid-rapidity multiplicity"};
262+
Configurable<float> cfgGenMultEtaMin{"cfgGenMultEtaMin", -0.5, "Min Eta for generated mid-rapidity multiplicity"};
257263

258264
float lCentrality;
259265

@@ -422,6 +428,8 @@ struct Chk892LI {
422428
// MC
423429
if (doprocessMC) {
424430

431+
AxisSpec genMultAxis{500, -0.5, 499.5, "N_{ch}^{gen} (|#eta|<0.5)"};
432+
425433
histos.add("QACent_woCut", "Centrality without cut", HistType::kTH1F, {centAxis});
426434
histos.add("QACent_woCentCut", "Centrality without cent cut", HistType::kTH1F, {centAxis});
427435
histos.add("QACent_wCentCut", "Centrality with cent cut", HistType::kTH1F, {centAxis});
@@ -445,7 +453,18 @@ struct Chk892LI {
445453
histos.add("Correction/sigLoss_num_pri_pos", "Gen primary Kstar selected by vertex position (|y|<0.5, selected events) in reco class", HistType::kTH2F, {ptAxis, centAxis});
446454
histos.add("Correction/EF_den", "Gen events (truth class)", HistType::kTH1F, {centAxis});
447455
histos.add("Correction/EF_num", "Reco events (selected events)", HistType::kTH1F, {centAxis});
456+
histos.add("Correction/RecoCentVsGenMult", "Reco centrality vs generated mid-rapidity multiplicity", HistType::kTH2F, {centAxis, genMultAxis});
457+
histos.add("Correction/sigLoss_num_vsGenMult", "Generated Kstar in selected events vs gen mult", HistType::kTH2F, {ptAxis, genMultAxis});
458+
histos.add("Correction/sigLoss_num_vsGenMult_pri", "Generated Kstar in selected events vs gen mult", HistType::kTH2F, {ptAxis, genMultAxis});
459+
histos.add("Correction/sigLoss_num_vsGenMult_pri_pos", "Generated Kstar in selected events vs gen mult", HistType::kTH2F, {ptAxis, genMultAxis});
460+
histos.add("Correction/sigLoss_den_vsGenMult", "Generated Kstar vs generated mid-rapidity multiplicity", HistType::kTH2F, {ptAxis, genMultAxis});
461+
histos.add("Correction/sigLoss_den_vsGenMult_pri", "Generated primary Kstar vs generated mid-rapidity multiplicity", HistType::kTH2F, {ptAxis, genMultAxis});
462+
histos.add("Correction/sigLoss_den_vsGenMult_pri_pos", "Generated primary Kstar selected by vertex position vs generated mid-rapidity multiplicity", HistType::kTH2F, {ptAxis, genMultAxis});
463+
histos.add("Correction/EF_num_vsGenMult", "Selected reco-associated events vs gen mult", HistType::kTH1F, {genMultAxis});
464+
histos.add("Correction/EF_den_vsGenMult", "Truth selected generated events vs generated multiplicity", HistType::kTH1F, {genMultAxis});
465+
448466
histos.add("Correction/MCTruthCent_all", "MC truth FT0M centrality (all mcCollisions)", HistType::kTH1F, {centAxis});
467+
histos.add("Correction/MCTruthCent_allowed", "MC truth FT0M centrality (allowed mcCollisions)", HistType::kTH1F, {centAxis});
449468
histos.add("Correction/MCTruthCent_cut", "MC truth FT0M centrality (truth selection applied)", HistType::kTH1F, {centAxis});
450469

451470
histos.add("Correction/setSizes", "Sizes of sets", HistType::kTH1F, {{4, -0.5, 3.5}});
@@ -496,6 +515,53 @@ struct Chk892LI {
496515
}
497516
}
498517

518+
template <typename McPartsT>
519+
int getGenMidRapMultiplicity(McPartsT const& partsThisMc)
520+
{
521+
int nCh = 0;
522+
523+
for (auto const& part : partsThisMc) {
524+
if (!part.isPhysicalPrimary()) {
525+
continue;
526+
}
527+
if (part.eta() > cfgGenMultEtaMax || part.eta() < cfgGenMultEtaMin) {
528+
continue;
529+
}
530+
531+
auto pdgParticle = pdg->GetParticle(part.pdgCode());
532+
if (!pdgParticle) {
533+
continue;
534+
}
535+
if (pdgParticle->Charge() == 0) {
536+
continue;
537+
}
538+
539+
nCh++;
540+
}
541+
542+
return nCh;
543+
}
544+
545+
float getCentClassFromGenMult(int nCh)
546+
{
547+
const auto& cuts = AxisConfig.cfgGenMultCuts.value;
548+
const auto& centres = AxisConfig.cfgCentBinCentres.value;
549+
550+
if (cuts.size() != centres.size()) {
551+
LOGF(fatal,
552+
"cfgGenMultCuts size (%zu) and cfgCentBinCentres size (%zu) must be same",
553+
cuts.size(), centres.size());
554+
}
555+
556+
for (size_t i = 0; i < cuts.size(); ++i) {
557+
if (nCh >= cuts[i]) {
558+
return centres[i];
559+
}
560+
}
561+
562+
return kInvalidCentrality;
563+
}
564+
499565
// Track selection
500566
template <typename TrackType>
501567
bool trackCut(TrackType const& track)
@@ -560,8 +626,6 @@ struct Chk892LI {
560626
return false;
561627
if (PIDCuts.cfgTPConly)
562628
return true;
563-
// if (candidate.pt() <= PIDCuts.cfgTOFMinPt)
564-
// return true;
565629

566630
if (candidate.hasTOF()) {
567631
const bool tofPIDPassed = std::abs(candidate.tofNSigmaPi()) < PIDCuts.cfgMaxTOFnSigmaPion;
@@ -720,10 +784,10 @@ struct Chk892LI {
720784
std::unordered_map<int64_t, float> centTruthByAllowed;
721785
std::unordered_set<int64_t> refClassIds;
722786
std::unordered_map<int64_t, float> refCentByMcId;
787+
std::unordered_map<int64_t, int> genMultByMcId;
723788

724-
template <typename McCollsT>
725-
// void buildAllowedMcIds(McCollsT const& mcCollisions, RecoEventsT const& events)
726-
void buildAllowedMcIds(McCollsT const& mcCollisions, MCEventCandidates const& events)
789+
template <typename McCollsT, typename McPartsT>
790+
void buildAllowedMcIds(McCollsT const& mcCollisions, MCEventCandidates const& events, McPartsT const& mcparts)
727791
{
728792
allowedMcIds.clear();
729793
centTruthByAllowed.clear();
@@ -770,14 +834,6 @@ struct Chk892LI {
770834
histos.fill(HIST("QACent_woCentCut"), lCentrality);
771835
}
772836

773-
if (lCentrality < EventCuts.cfgEventCentralityMin || lCentrality > EventCuts.cfgEventCentralityMax) {
774-
continue;
775-
}
776-
777-
if (doprocessMC) {
778-
histos.fill(HIST("QACent_wCentCut"), lCentrality);
779-
}
780-
781837
atLeastOneMatch = true;
782838

783839
// choose best collision by largest number of PV contributors
@@ -792,11 +848,29 @@ struct Chk892LI {
792848
continue;
793849
}
794850

851+
auto partsThisMc = mcparts.sliceBy(perMCCollision, mcid);
852+
const int genMult = getGenMidRapMultiplicity(partsThisMc);
853+
const float genCentClass = getCentClassFromGenMult(genMult);
854+
855+
if (genCentClass == kInvalidCentrality)
856+
continue;
857+
858+
if (genCentClass < EventCuts.cfgEventCentralityMin || genCentClass > EventCuts.cfgEventCentralityMax) {
859+
continue;
860+
}
861+
862+
if (doprocessMC) {
863+
histos.fill(HIST("QACent_wCentCut"), bestCent);
864+
}
865+
795866
allowedMcIds.insert(mcid);
796-
centTruthByAllowed.emplace(mcid, bestCent);
867+
centTruthByAllowed.emplace(mcid, genCentClass);
868+
genMultByMcId[mcid] = genMult;
797869

798870
if (doprocessMC) {
799-
histos.fill(HIST("QAMCCent_allowed"), bestCent);
871+
histos.fill(HIST("Correction/EF_num_vsGenMult"), genMult);
872+
histos.fill(HIST("Correction/RecoCentVsGenMult"), bestCent, genMult);
873+
histos.fill(HIST("QAMCCent_allowed"), genCentClass);
800874
}
801875
}
802876
} // buildAllowedMcIds
@@ -808,31 +882,33 @@ struct Chk892LI {
808882
refCentByMcId.clear();
809883

810884
for (const auto& coll : mccolls) {
811-
bool pass = true;
812-
813-
if (cfgTruthIncludeZvtx && std::abs(coll.posZ()) >= EventCuts.cfgEvtZvtx)
814-
pass = false;
885+
const auto mcid = coll.globalIndex();
815886

816-
if (pass && cfgTruthUseInelGt0) {
817-
auto partsThisMc = mcparts.sliceBy(perMCCollision, coll.globalIndex());
818-
if (!pwglf::isINELgtNmc(partsThisMc, 0, pdg))
819-
pass = false;
887+
if (cfgTruthIncludeZvtx && std::abs(coll.posZ()) >= EventCuts.cfgEvtZvtx) {
888+
continue;
820889
}
821890

822-
if (!pass)
891+
auto partsThisMc = mcparts.sliceBy(perMCCollision, mcid);
892+
893+
if (cfgTruthUseInelGt0 && !pwglf::isINELgtNmc(partsThisMc, 0, pdg)) {
823894
continue;
895+
}
824896

825-
const auto mcid = coll.globalIndex();
897+
const int genMult = getGenMidRapMultiplicity(partsThisMc);
898+
const float genCentClass = getCentClassFromGenMult(genMult);
826899

827-
auto it = centTruthByAllowed.find(mcid);
828-
if (it == centTruthByAllowed.end()) {
900+
if (genCentClass == kInvalidCentrality) {
901+
continue;
902+
}
903+
if (genCentClass < EventCuts.cfgEventCentralityMin || genCentClass > EventCuts.cfgEventCentralityMax) {
829904
continue;
830905
}
831906

832907
refClassIds.insert(mcid);
908+
refCentByMcId.emplace(mcid, genCentClass);
909+
genMultByMcId[mcid] = genMult;
833910

834-
const float lCentrality = it->second;
835-
refCentByMcId.emplace(mcid, lCentrality);
911+
histos.fill(HIST("Correction/EF_den_vsGenMult"), genMult);
836912

837913
} // for
838914
} // buildReferenceMcIds
@@ -1115,8 +1191,16 @@ struct Chk892LI {
11151191

11161192
const float lCentrality = iter->second;
11171193

1194+
auto itMult = genMultByMcId.find(mcid);
1195+
if (itMult != genMultByMcId.end()) {
1196+
histos.fill(HIST("Correction/sigLoss_num_vsGenMult"), part.pt(), itMult->second);
1197+
}
11181198
histos.fill(HIST("Correction/sigLoss_num"), part.pt(), lCentrality);
1199+
11191200
if (part.vt() == 0) {
1201+
if (itMult != genMultByMcId.end()) {
1202+
histos.fill(HIST("Correction/sigLoss_num_vsGenMult_pri"), part.pt(), itMult->second);
1203+
}
11201204
histos.fill(HIST("Correction/sigLoss_num_pri"), part.pt(), lCentrality);
11211205
}
11221206

@@ -1129,6 +1213,9 @@ struct Chk892LI {
11291213
const float distanceFromPV = std::sqrt(dx * dx + dy * dy + dz * dz);
11301214

11311215
if (distanceFromPV < fMaxPosPV) {
1216+
if (itMult != genMultByMcId.end()) {
1217+
histos.fill(HIST("Correction/sigLoss_num_vsGenMult_pri_pos"), part.pt(), itMult->second);
1218+
}
11321219
histos.fill(HIST("Correction/sigLoss_num_pri_pos"), part.pt(), lCentrality);
11331220
}
11341221
}
@@ -1155,8 +1242,16 @@ struct Chk892LI {
11551242
const float lCentrality = iter->second;
11561243

11571244
histos.fill(HIST("Correction/sigLoss_den"), part.pt(), lCentrality);
1245+
1246+
auto itMult = genMultByMcId.find(mcid);
1247+
if (itMult != genMultByMcId.end()) {
1248+
histos.fill(HIST("Correction/sigLoss_den_vsGenMult"), part.pt(), itMult->second);
1249+
}
11581250
if (part.vt() == 0) {
11591251
histos.fill(HIST("Correction/sigLoss_den_pri"), part.pt(), lCentrality);
1252+
if (itMult != genMultByMcId.end()) {
1253+
histos.fill(HIST("Correction/sigLoss_den_vsGenMult_pri"), part.pt(), itMult->second);
1254+
}
11601255
}
11611256

11621257
const auto mcc = part.mcCollision_as<MCTrueEventCandidates>();
@@ -1169,6 +1264,9 @@ struct Chk892LI {
11691264

11701265
if (distanceFromPV < fMaxPosPV) {
11711266
histos.fill(HIST("Correction/sigLoss_den_pri_pos"), part.pt(), lCentrality);
1267+
if (itMult != genMultByMcId.end()) {
1268+
histos.fill(HIST("Correction/sigLoss_den_vsGenMult_pri_pos"), part.pt(), itMult->second);
1269+
}
11721270
}
11731271
}
11741272
} // fillSigLossDen
@@ -1428,7 +1526,8 @@ struct Chk892LI {
14281526
MCEventCandidates const& events,
14291527
MCTrueEventCandidates const& mccolls)
14301528
{
1431-
buildAllowedMcIds(mccolls, events);
1529+
genMultByMcId.clear();
1530+
buildAllowedMcIds(mccolls, events, mcpart);
14321531
buildReferenceMcIds(mccolls, mcpart);
14331532
effK0sProcessGen(mcpart);
14341533
effK0sProcessReco(v0s);
@@ -1445,6 +1544,7 @@ struct Chk892LI {
14451544
const float lCentrality = iter->second;
14461545
histos.fill(HIST("Correction/EF_den"), lCentrality);
14471546
}
1547+
14481548
for (const auto& mcid : allowedMcIds) {
14491549
auto iter = centTruthByAllowed.find(mcid);
14501550
if (iter == centTruthByAllowed.end())
@@ -1463,8 +1563,14 @@ struct Chk892LI {
14631563
histos.fill(HIST("Correction/setSizes"), 2.0, nIntersect);
14641564
histos.fill(HIST("Correction/setSizes"), 3.0, allowedMcIds.size() - nIntersect);
14651565

1466-
for (auto const& [mcid, lCentrality] : centTruthByAllowed) {
1467-
histos.fill(HIST("Correction/MCTruthCent_all"), lCentrality);
1566+
for (auto const& mcc : mccolls) {
1567+
auto partsThisMc = mcpart.sliceBy(perMCCollision, mcc.globalIndex());
1568+
const int genMult = getGenMidRapMultiplicity(partsThisMc);
1569+
const float genCentClass = getCentClassFromGenMult(genMult);
1570+
if (genCentClass == kInvalidCentrality) {
1571+
continue;
1572+
}
1573+
histos.fill(HIST("Correction/MCTruthCent_all"), genCentClass);
14681574
}
14691575

14701576
for (const auto& mcid : refClassIds) {
@@ -1488,7 +1594,7 @@ struct Chk892LI {
14881594
histos.fill(HIST("Correction/hNEventsMCTruth"), 2.0);
14891595

14901596
auto partsThisMc = mcpart.sliceBy(perMCCollision, mcid);
1491-
if (pwglf::isINELgtNmc(partsThisMc, 0, pdg)) {
1597+
if (!cfgTruthUseInelGt0 || pwglf::isINELgtNmc(partsThisMc, 0, pdg)) {
14921598
histos.fill(HIST("Correction/hNEventsMCTruth"), 3.0);
14931599
}
14941600
}

0 commit comments

Comments
 (0)