Skip to content

Commit ae2f69a

Browse files
MaximVirtaMaxim Virtaalibuild
authored
[PWGCF] PID QA improved & PID efficiency added (#15504)
Co-authored-by: Maxim Virta <maximus@Maxims-MacBook-Pro.local> Co-authored-by: ALICE Action Bot <alibuild@cern.ch>
1 parent 4799092 commit ae2f69a

File tree

1 file changed

+111
-67
lines changed

1 file changed

+111
-67
lines changed

‎PWGCF/GenericFramework/Tasks/flowGfwV02.cxx‎

Lines changed: 111 additions & 67 deletions
Original file line numberDiff line numberDiff line change
@@ -105,6 +105,7 @@ struct FlowGfwV02 {
105105
O2_DEFINE_CONFIGURABLE(cfgUseAdditionalEventCut, bool, false, "Use additional event cut on mult correlations")
106106
O2_DEFINE_CONFIGURABLE(cfgEfficiency, std::string, "", "CCDB path to efficiency object")
107107
O2_DEFINE_CONFIGURABLE(cfgAcceptance, std::string, "", "CCDB path to acceptance object")
108+
O2_DEFINE_CONFIGURABLE(cfgPIDEfficiency, bool, false, "Use PID efficiency for efficiency calculation")
108109
O2_DEFINE_CONFIGURABLE(cfgFixedMultMin, int, 1, "Minimum for fixed nch range");
109110
O2_DEFINE_CONFIGURABLE(cfgFixedMultMax, int, 3000, "Maximum for fixed nch range");
110111
O2_DEFINE_CONFIGURABLE(cfgTofPtCut, float, 0.5f, "Minimum pt to use TOF N-sigma")
@@ -186,7 +187,7 @@ struct FlowGfwV02 {
186187
Service<ccdb::BasicCCDBManager> ccdb;
187188

188189
struct Config {
189-
TH1D* mEfficiency = nullptr;
190+
TH1D* mEfficiency[4] = {nullptr, nullptr, nullptr, nullptr};
190191
GFWWeights* mAcceptance;
191192
bool correctionsLoaded = false;
192193
} cfg;
@@ -306,17 +307,24 @@ struct FlowGfwV02 {
306307
pidStates.itsNsigmaCut[IndProtonLow] = nSigmas->getData()[IndProtonLow][kITS];
307308

308309
if (cfgGetNsigmaQA) {
309-
if (!cfgUseItsPID) {
310-
registry.add("TofTpcNsigma_before", "", {HistType::kTHnSparseD, {cfgPIDHistograms.cfgAxisNsigmaTPC, cfgPIDHistograms.cfgAxisNsigmaTOF, cfgPIDHistograms.cfgAxisPt}});
311-
registry.add("TofTpcNsigma_after", "", {HistType::kTHnSparseD, {cfgPIDHistograms.cfgAxisNsigmaTPC, cfgPIDHistograms.cfgAxisNsigmaTOF, cfgPIDHistograms.cfgAxisPt}});
312-
}
313310
if (cfgUseItsPID) {
314-
registry.add("TofItsNsigma_before", "", {HistType::kTHnSparseD, {cfgPIDHistograms.cfgAxisNsigmaITS, cfgPIDHistograms.cfgAxisNsigmaTOF, cfgPIDHistograms.cfgAxisPt}});
315-
registry.add("TofItsNsigma_after", "", {HistType::kTHnSparseD, {cfgPIDHistograms.cfgAxisNsigmaITS, cfgPIDHistograms.cfgAxisNsigmaTOF, cfgPIDHistograms.cfgAxisPt}});
311+
registry.add("QA_PID/before/TofItsNsigma", "", {HistType::kTHnSparseD, {cfgPIDHistograms.cfgAxisNsigmaITS, cfgPIDHistograms.cfgAxisNsigmaTOF, cfgPIDHistograms.cfgAxisPt}});
312+
} else {
313+
registry.add("QA_PID/before/TofTpcNsigma_pions", "", {HistType::kTHnSparseD, {cfgPIDHistograms.cfgAxisNsigmaTPC, cfgPIDHistograms.cfgAxisNsigmaTOF, cfgPIDHistograms.cfgAxisPt}});
314+
registry.add("QA_PID/before/TofTpcNsigma_kaons", "", {HistType::kTHnSparseD, {cfgPIDHistograms.cfgAxisNsigmaTPC, cfgPIDHistograms.cfgAxisNsigmaTOF, cfgPIDHistograms.cfgAxisPt}});
315+
registry.add("QA_PID/before/TofTpcNsigma_protons", "", {HistType::kTHnSparseD, {cfgPIDHistograms.cfgAxisNsigmaTPC, cfgPIDHistograms.cfgAxisNsigmaTOF, cfgPIDHistograms.cfgAxisPt}});
316316
}
317317

318-
registry.add("TpcdEdx_ptwise", "", {HistType::kTH2D, {cfgPIDHistograms.cfgAxisTpcSignal, cfgPIDHistograms.cfgAxisPt}});
319-
registry.add("TpcdEdx_ptwise_afterCut", "", {HistType::kTH2D, {cfgPIDHistograms.cfgAxisTpcSignal, cfgPIDHistograms.cfgAxisPt}});
318+
registry.add("QA_PID/before/TpcdEdx_ptwise_pions", "", {HistType::kTH2D, {cfgPIDHistograms.cfgAxisTpcSignal, cfgPIDHistograms.cfgAxisPt}});
319+
registry.add("QA_PID/before/ExpTpcdEdx_ptwise_pions", "", {HistType::kTH2D, {cfgPIDHistograms.cfgAxisTpcSignal, cfgPIDHistograms.cfgAxisPt}});
320+
registry.add("QA_PID/before/ExpSigma_ptwise_pions", "", {HistType::kTH2D, {cfgPIDHistograms.cfgAxisTpcSignal, cfgPIDHistograms.cfgAxisPt}});
321+
registry.add("QA_PID/before/TpcdEdx_ptwise_kaons", "", {HistType::kTH2D, {cfgPIDHistograms.cfgAxisTpcSignal, cfgPIDHistograms.cfgAxisPt}});
322+
registry.add("QA_PID/before/ExpTpcdEdx_ptwise_kaons", "", {HistType::kTH2D, {cfgPIDHistograms.cfgAxisTpcSignal, cfgPIDHistograms.cfgAxisPt}});
323+
registry.add("QA_PID/before/ExpSigma_ptwise_kaons", "", {HistType::kTH2D, {cfgPIDHistograms.cfgAxisTpcSignal, cfgPIDHistograms.cfgAxisPt}});
324+
registry.add("QA_PID/before/TpcdEdx_ptwise_protons", "", {HistType::kTH2D, {cfgPIDHistograms.cfgAxisTpcSignal, cfgPIDHistograms.cfgAxisPt}});
325+
registry.add("QA_PID/before/ExpTpcdEdx_ptwise_protons", "", {HistType::kTH2D, {cfgPIDHistograms.cfgAxisTpcSignal, cfgPIDHistograms.cfgAxisPt}});
326+
registry.add("QA_PID/before/ExpSigma_ptwise_protons", "", {HistType::kTH2D, {cfgPIDHistograms.cfgAxisTpcSignal, cfgPIDHistograms.cfgAxisPt}});
327+
registry.addClone("QA_PID/before/", "QA_PID/after/");
320328
}
321329

322330
o2::analysis::gfw::regions.SetNames(cfgRegions->GetNames());
@@ -588,11 +596,22 @@ struct FlowGfwV02 {
588596
cfg.mAcceptance = ccdb->getForTimeStamp<GFWWeights>(cfgAcceptance.value, timestamp);
589597
}
590598
if (!cfgEfficiency.value.empty()) {
591-
cfg.mEfficiency = ccdb->getForTimeStamp<TH1D>(cfgEfficiency, timestamp);
592-
if (cfg.mEfficiency == nullptr) {
599+
cfg.mEfficiency[PidCharged] = ccdb->getForTimeStamp<TH1D>(cfgEfficiency, timestamp);
600+
if (cfg.mEfficiency[PidCharged] == nullptr) {
593601
LOGF(fatal, "Could not load efficiency histogram from %s", cfgEfficiency.value.c_str());
594602
}
595-
LOGF(info, "Loaded efficiency histogram from %s (%p)", cfgEfficiency.value.c_str(), (void*)cfg.mEfficiency);
603+
LOGF(info, "Loaded efficiency histogram from %s (%p)", cfgEfficiency.value.c_str(), (void*)cfg.mEfficiency[PidCharged]);
604+
}
605+
if (cfgPIDEfficiency) {
606+
const std::array<std::string, 4> pidStrings = {"ch", "pi", "ka", "pr"};
607+
for (int i = 1; i < 4; i++) {
608+
609+
cfg.mEfficiency[i] = ccdb->getForTimeStamp<TH1D>(cfgEfficiency.value + pidStrings[i], timestamp);
610+
if (cfg.mEfficiency[i] == nullptr) {
611+
LOGF(fatal, "Could not load PID efficiency histogram from %s", cfgEfficiency.value + pidStrings[i].c_str());
612+
}
613+
LOGF(info, "Loaded PID efficiency histogram from %s (%p)", cfgEfficiency.value + pidStrings[i].c_str(), (void*)cfg.mEfficiency[i]);
614+
}
596615
}
597616
cfg.correctionsLoaded = true;
598617
}
@@ -605,7 +624,7 @@ struct FlowGfwV02 {
605624
cfg.mAcceptance = ccdb->getForRun<GFWWeights>(cfgAcceptance.value, runnumber);
606625
}
607626
if (!cfgEfficiency.value.empty()) {
608-
cfg.mEfficiency = ccdb->getForRun<TH1D>(cfgEfficiency.value, runnumber);
627+
cfg.mEfficiency[PidCharged] = ccdb->getForRun<TH1D>(cfgEfficiency.value, runnumber);
609628
}
610629
cfg.correctionsLoaded = true;
611630
}
@@ -638,66 +657,17 @@ struct FlowGfwV02 {
638657
}
639658

640659
template <typename TTrack>
641-
double getEfficiency(TTrack track)
660+
double getEfficiency(TTrack track, const int& pid = PidCharged)
642661
{
643662
double eff = 1.;
644-
if (cfg.mEfficiency)
645-
eff = cfg.mEfficiency->GetBinContent(cfg.mEfficiency->FindBin(track.pt()));
663+
if (cfg.mEfficiency[pid])
664+
eff = cfg.mEfficiency[pid]->GetBinContent(cfg.mEfficiency[pid]->FindBin(track.pt()));
646665
if (eff == 0)
647666
return -1.;
648667
else
649668
return 1. / eff;
650669
}
651670

652-
template <typename TTrack>
653-
void fillNsigmaAfterCut(TTrack track1, Int_t pid) // function to fill the QA after Nsigma selection
654-
{
655-
switch (pid) {
656-
case 1: // For Pions
657-
if (!cfgUseItsPID) {
658-
if (cfgGetdEdx) {
659-
double tpcExpSignalPi = track1.tpcSignal() - (track1.tpcNSigmaPi() * track1.tpcExpSigmaPi());
660-
661-
registry.fill(HIST("TpcdEdx_ptwise_afterCut"), track1.pt(), track1.tpcSignal(), track1.tofNSigmaPi());
662-
registry.fill(HIST("ExpTpcdEdx_ptwise_afterCut"), track1.pt(), tpcExpSignalPi, track1.tofNSigmaPi());
663-
registry.fill(HIST("ExpSigma_ptwise_afterCut"), track1.pt(), track1.tpcExpSigmaPi(), track1.tofNSigmaPi());
664-
}
665-
registry.fill(HIST("TofTpcNsigma_after"), track1.tpcNSigmaPi(), track1.tofNSigmaPi(), track1.pt());
666-
}
667-
if (cfgUseItsPID)
668-
registry.fill(HIST("TofItsNsigma_after"), pidStates.itsResponse.nSigmaITS<o2::track::PID::Pion>(track1), track1.tofNSigmaPi(), track1.pt());
669-
break;
670-
case 2: // For Kaons
671-
if (!cfgUseItsPID) {
672-
if (cfgGetdEdx) {
673-
double tpcExpSignalKa = track1.tpcSignal() - (track1.tpcNSigmaKa() * track1.tpcExpSigmaKa());
674-
675-
registry.fill(HIST("TpcdEdx_ptwise_afterCut"), track1.pt(), track1.tpcSignal(), track1.tofNSigmaKa());
676-
registry.fill(HIST("ExpTpcdEdx_ptwise_afterCut"), track1.pt(), tpcExpSignalKa, track1.tofNSigmaKa());
677-
registry.fill(HIST("ExpSigma_ptwise_afterCut"), track1.pt(), track1.tpcExpSigmaKa(), track1.tofNSigmaKa());
678-
}
679-
registry.fill(HIST("TofTpcNsigma_after"), track1.tpcNSigmaKa(), track1.tofNSigmaKa(), track1.pt());
680-
}
681-
if (cfgUseItsPID)
682-
registry.fill(HIST("TofItsNsigma_after"), pidStates.itsResponse.nSigmaITS<o2::track::PID::Kaon>(track1), track1.tofNSigmaKa(), track1.pt());
683-
break;
684-
case 3: // For Protons
685-
if (!cfgUseItsPID) {
686-
if (cfgGetdEdx) {
687-
double tpcExpSignalPr = track1.tpcSignal() - (track1.tpcNSigmaPr() * track1.tpcExpSigmaPr());
688-
689-
registry.fill(HIST("TpcdEdx_ptwise_afterCut"), track1.pt(), track1.tpcSignal(), track1.tofNSigmaPr());
690-
registry.fill(HIST("ExpTpcdEdx_ptwise_afterCut"), track1.pt(), tpcExpSignalPr, track1.tofNSigmaPr());
691-
registry.fill(HIST("ExpSigma_ptwise_afterCut"), track1.pt(), track1.tpcExpSigmaPr(), track1.tofNSigmaPr());
692-
}
693-
registry.fill(HIST("TofTpcNsigma_after"), track1.tpcNSigmaPr(), track1.tofNSigmaPr(), track1.pt());
694-
}
695-
if (cfgUseItsPID)
696-
registry.fill(HIST("TofItsNsigma_after"), pidStates.itsResponse.nSigmaITS<o2::track::PID::Proton>(track1), track1.tofNSigmaPr(), track1.pt());
697-
break;
698-
} // end of switch
699-
}
700-
701671
template <typename TCollision>
702672
bool eventSelected(TCollision collision, const int& multTrk, const float& centrality)
703673
{
@@ -930,11 +900,14 @@ struct FlowGfwV02 {
930900
for (const auto& track : tracks) {
931901
processTrack(track, vtxz, xaxis.multiplicity, run, acceptedTracks);
932902
if (track.eta() > -0.4 && track.eta() < 0.4)
933-
pidStates.hPtMid[PidCharged]->Fill(track.pt(), getEfficiency(track));
903+
pidStates.hPtMid[PidCharged]->Fill(track.pt(), getEfficiency(track, PidCharged));
934904
// If PID is identified, fill pt spectrum for the corresponding particle
935905
int pidInd = getNsigmaPID(track);
936906
if (pidInd != -1 && track.eta() > -0.4 && track.eta() < 0.4) {
937-
pidStates.hPtMid[pidInd]->Fill(track.pt(), getEfficiency(track)); // TODO: Add PID index to the efficiency histogram
907+
if (cfgPIDEfficiency)
908+
pidStates.hPtMid[pidInd]->Fill(track.pt(), getEfficiency(track, pidInd));
909+
else
910+
pidStates.hPtMid[pidInd]->Fill(track.pt(), getEfficiency(track, PidCharged)); // Default to charged particles if PID efficiency is not used
938911
}
939912
}
940913
if (cfgConsistentEventFlag & 1)
@@ -1005,6 +978,10 @@ struct FlowGfwV02 {
1005978
fillTrackQA<kBefore>(track, vtxz);
1006979
registry.fill(HIST("trackQA/before/nch_pt"), multiplicity, track.pt());
1007980
}
981+
982+
if (cfgGetNsigmaQA)
983+
fillPidQA<kBefore>(track, getNsigmaPID(track));
984+
1008985
if (!trackSelected(track))
1009986
return;
1010987

@@ -1014,6 +991,9 @@ struct FlowGfwV02 {
1014991
fillTrackQA<kAfter>(track, vtxz);
1015992
registry.fill(HIST("trackQA/after/nch_pt"), multiplicity, track.pt());
1016993
}
994+
995+
if (cfgGetNsigmaQA)
996+
fillPidQA<kAfter>(track, getNsigmaPID(track));
1017997
}
1018998

1019999
template <DataType dt, typename TTrack>
@@ -1045,6 +1025,70 @@ struct FlowGfwV02 {
10451025
return;
10461026
}
10471027

1028+
template <QAFillTime ft, typename TTrack>
1029+
inline void fillPidQA(TTrack track, const int& pid)
1030+
{
1031+
// Fill Nsigma QA
1032+
if (!cfgUseItsPID) {
1033+
if (ft == kBefore) {
1034+
registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TofTpcNsigma_pions"), track.tpcNSigmaPi(), track.tofNSigmaPi(), track.pt());
1035+
if (cfgGetdEdx) {
1036+
double tpcExpSignalPi = track.tpcSignal() - (track.tpcNSigmaPi() * track.tpcExpSigmaPi());
1037+
1038+
registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TpcdEdx_ptwise_pions"), track.pt(), track.tpcSignal(), track.tofNSigmaPi());
1039+
registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpTpcdEdx_ptwise_pions"), track.pt(), tpcExpSignalPi, track.tofNSigmaPi());
1040+
registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpSigma_ptwise_pions"), track.pt(), track.tpcExpSigmaPi(), track.tofNSigmaPi());
1041+
}
1042+
registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TofTpcNsigma_kaons"), track.tpcNSigmaKa(), track.tofNSigmaKa(), track.pt());
1043+
if (cfgGetdEdx) {
1044+
double tpcExpSignalKa = track.tpcSignal() - (track.tpcNSigmaKa() * track.tpcExpSigmaKa());
1045+
1046+
registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TpcdEdx_ptwise_kaons"), track.pt(), track.tpcSignal(), track.tofNSigmaKa());
1047+
registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpTpcdEdx_ptwise_kaons"), track.pt(), tpcExpSignalKa, track.tofNSigmaKa());
1048+
registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpSigma_ptwise_kaons"), track.pt(), track.tpcExpSigmaKa(), track.tofNSigmaKa());
1049+
}
1050+
registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TofTpcNsigma_protons"), track.tpcNSigmaPr(), track.tofNSigmaPr(), track.pt());
1051+
if (cfgGetdEdx) {
1052+
double tpcExpSignalPr = track.tpcSignal() - (track.tpcNSigmaPr() * track.tpcExpSigmaPr());
1053+
1054+
registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TpcdEdx_ptwise_protons"), track.pt(), track.tpcSignal(), track.tofNSigmaPr());
1055+
registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpTpcdEdx_ptwise_protons"), track.pt(), tpcExpSignalPr, track.tofNSigmaPr());
1056+
registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpSigma_ptwise_protons"), track.pt(), track.tpcExpSigmaPr(), track.tofNSigmaPr());
1057+
}
1058+
} else if (ft == kAfter) {
1059+
if (pid == PidPions) {
1060+
registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TofTpcNsigma_pions"), track.tpcNSigmaPi(), track.tofNSigmaPi(), track.pt());
1061+
if (cfgGetdEdx) {
1062+
double tpcExpSignalPi = track.tpcSignal() - (track.tpcNSigmaPi() * track.tpcExpSigmaPi());
1063+
1064+
registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TpcdEdx_ptwise_pions"), track.pt(), track.tpcSignal(), track.tofNSigmaPi());
1065+
registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpTpcdEdx_ptwise_pions"), track.pt(), tpcExpSignalPi, track.tofNSigmaPi());
1066+
registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpSigma_ptwise_pions"), track.pt(), track.tpcExpSigmaPi(), track.tofNSigmaPi());
1067+
}
1068+
}
1069+
if (pid == PidKaons) {
1070+
registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TofTpcNsigma_kaons"), track.tpcNSigmaKa(), track.tofNSigmaKa(), track.pt());
1071+
if (cfgGetdEdx) {
1072+
double tpcExpSignalKa = track.tpcSignal() - (track.tpcNSigmaKa() * track.tpcExpSigmaKa());
1073+
1074+
registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TpcdEdx_ptwise_kaons"), track.pt(), track.tpcSignal(), track.tofNSigmaKa());
1075+
registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpTpcdEdx_ptwise_kaons"), track.pt(), tpcExpSignalKa, track.tofNSigmaKa());
1076+
registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpSigma_ptwise_kaons"), track.pt(), track.tpcExpSigmaKa(), track.tofNSigmaKa());
1077+
}
1078+
}
1079+
if (pid == PidProtons) {
1080+
registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TofTpcNsigma_protons"), track.tpcNSigmaPr(), track.tofNSigmaPr(), track.pt());
1081+
if (cfgGetdEdx) {
1082+
double tpcExpSignalPr = track.tpcSignal() - (track.tpcNSigmaPr() * track.tpcExpSigmaPr());
1083+
registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("TpcdEdx_ptwise_protons"), track.pt(), track.tpcSignal(), track.tofNSigmaPr());
1084+
registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpTpcdEdx_ptwise_protons"), track.pt(), tpcExpSignalPr, track.tofNSigmaPr());
1085+
registry.fill(HIST("QA_PID/") + HIST(FillTimeName[ft]) + HIST("ExpSigma_ptwise_protons"), track.pt(), track.tpcExpSigmaPr(), track.tofNSigmaPr());
1086+
}
1087+
}
1088+
}
1089+
}
1090+
}
1091+
10481092
template <QAFillTime ft, typename TTrack>
10491093
inline void fillTrackQA(TTrack track, const float vtxz)
10501094
{

0 commit comments

Comments
 (0)