From 589e749004c8f08e12d2366a68041c3f1308e356 Mon Sep 17 00:00:00 2001 From: saganatt <8majak8@gmail.com> Date: Sun, 23 Jan 2022 00:03:51 +0100 Subject: [PATCH 1/4] Added Bayes task, uncertain about the results --- Tools/PIDML/qaPid.cxx | 137 +++++++++++++++++++----------------------- 1 file changed, 63 insertions(+), 74 deletions(-) diff --git a/Tools/PIDML/qaPid.cxx b/Tools/PIDML/qaPid.cxx index d50319204a5..6dcea38b07c 100644 --- a/Tools/PIDML/qaPid.cxx +++ b/Tools/PIDML/qaPid.cxx @@ -153,6 +153,27 @@ struct pid { } } + template + void fillPidHistos(const T& track, const int pdgCode, bool isPidTrue) + { + if (isPidTrue) { + histReg.fill(HIST(pidTrueRegistryNames[i]), track.pt()); + histReg.fill(HIST(TPCPidTrueRegistryNames[i]), track.p(), track.tpcSignal()); + histReg.fill(HIST("TPCSignalPidTrue"), track.p(), track.tpcSignal()); + histReg.fill(HIST("TOFSignalPidTrue"), track.p(), track.beta()); + histReg.fill(HIST(TOFPidTrueRegistryNames[i]), track.p(), track.beta()); + } else { + histReg.fill(HIST(pidFalseRegistryNames[i]), track.pt()); + histReg.fill(HIST(TPCPidFalseRegistryNames[i]), track.p(), track.tpcSignal()); + histReg.fill(HIST("TPCSignalPidFalse"), track.p(), track.tpcSignal()); + histReg.fill(HIST("TOFSignalPidFalse"), track.p(), track.beta()); + histReg.fill(HIST(TOFPidFalseRegistryNames[i]), track.p(), track.beta()); + + double pt = track.pt(); + fillContaminationRegistry(i, pdgCode, pt); + } + } + // nb of particles (5 particles - Pi, Pr, Ka, e, mu, and 5 antiparticles) static const int numParticles = 10; @@ -166,6 +187,7 @@ struct pid { static constexpr std::string_view TOFPidFalseRegistryNames[numParticles] = {"TOFPidFalse/211", "TOFPidFalse/2212", "TOFPidFalse/321", "TOFPidFalse/11", "TOFPidFalse/13", "TOFPidFalse/0211", "TOFPidFalse/02212", "TOFPidFalse/0321", "TOFPidFalse/011", "TOFPidFalse/013"}; static constexpr int pdgCodes[numParticles] = {211, 2212, 321, 11, 13, -211, -2212, -321, -11, -13}; + static constexpr int pidToPdg[numParticles] = {11, 13, 211, 321, 2212, -11, -13, -211, -321, -2212}; // charge with index i corresponds to i-th particle from pdgCodes array static constexpr int particleCharge[numParticles] = {1, 1, 1, -1, -1, -1, -1, -1, 1, 1}; // momentum value when to switch from pid based only on TPC (below the value) to combination of TPC and TOF (above the value) @@ -310,44 +332,23 @@ struct pid { */ const float p = track.p(); - bool isPidFalse = false; - if ((p < pSwitch[i]) & (track.sign() == particleCharge[i])) { if (abs(tpcNSigmas[i]) < nsigmacut.value) { if (pdgCode == pdgCodes[i]) { - histReg.fill(HIST(pidTrueRegistryNames[i]), track.pt()); - histReg.fill(HIST(TPCPidTrueRegistryNames[i]), track.p(), track.tpcSignal()); - histReg.fill(HIST("TPCSignalPidTrue"), track.p(), track.tpcSignal()); - histReg.fill(HIST("TOFSignalPidTrue"), track.p(), track.beta()); - histReg.fill(HIST(TOFPidTrueRegistryNames[i]), track.p(), track.beta()); + fillPidHistos(track, pdgCode, true); } else { - isPidFalse = true; + fillPidHistos(track, pdgCode, false); } } } else if ((p >= pSwitch[i]) & (track.sign() == particleCharge[i])) { if (sqrt(pow(tpcNSigmas[i], 2) + pow(tofNSigmas[i], 2)) < nsigmacut.value) { if (pdgCode == pdgCodes[i]) { - histReg.fill(HIST(pidTrueRegistryNames[i]), track.pt()); - histReg.fill(HIST(TPCPidTrueRegistryNames[i]), track.p(), track.tpcSignal()); - histReg.fill(HIST("TPCSignalPidTrue"), track.p(), track.tpcSignal()); - histReg.fill(HIST("TOFSignalPidTrue"), track.p(), track.beta()); - histReg.fill(HIST(TOFPidTrueRegistryNames[i]), track.p(), track.beta()); + fillPidHistos(track, pdgCode, true); } else { - isPidFalse = true; + fillPidHistos(track, pdgCode, false); } } } - - if (isPidFalse) { - histReg.fill(HIST(pidFalseRegistryNames[i]), track.pt()); - histReg.fill(HIST(TPCPidFalseRegistryNames[i]), track.p(), track.tpcSignal()); - histReg.fill(HIST("TPCSignalPidFalse"), track.p(), track.tpcSignal()); - histReg.fill(HIST("TOFSignalPidFalse"), track.p(), track.beta()); - histReg.fill(HIST(TOFPidFalseRegistryNames[i]), track.p(), track.beta()); - - double pt = track.pt(); - fillContaminationRegistry(i, pdgCode, pt); - } } template @@ -358,8 +359,6 @@ struct pid { // list of Nsigmas for particles float particleNSigma[arrLen]; - bool isPidFalse = false; - // calculate Nsigmas for every particle for (int j = 0; j < arrLen; ++j) { if (p < pSwitch[j]) { @@ -373,40 +372,21 @@ struct pid { float tmp_NSigma = abs(tpcNSigmas[i]); if ((tmp_NSigma < nsigmacut.value) & (indexOfSmallestElement(particleNSigma, arrLen) == i)) { if (pdgCode == pdgCodes[i]) { - histReg.fill(HIST(pidTrueRegistryNames[i]), track.pt()); - histReg.fill(HIST(TPCPidTrueRegistryNames[i]), track.p(), track.tpcSignal()); - histReg.fill(HIST("TPCSignalPidTrue"), track.p(), track.tpcSignal()); - histReg.fill(HIST("TOFSignalPidTrue"), track.p(), track.beta()); - histReg.fill(HIST(TOFPidTrueRegistryNames[i]), track.p(), track.beta()); + fillPidHistos(track, pdgCode, true); } else { - isPidFalse = true; + fillPidHistos(track, pdgCode, false); } } } else if ((p >= pSwitch[i]) & (track.sign() == particleCharge[i])) { float tmp_NSigma = combinedSignal(tpcNSigmas[i], tofNSigmas[i]); if ((tmp_NSigma < nsigmacut.value) & (indexOfSmallestElement(particleNSigma, arrLen) == i)) { if (pdgCode == pdgCodes[i]) { - histReg.fill(HIST(pidTrueRegistryNames[i]), track.pt()); - histReg.fill(HIST(TPCPidTrueRegistryNames[i]), track.p(), track.tpcSignal()); - histReg.fill(HIST("TPCSignalPidTrue"), track.p(), track.tpcSignal()); - histReg.fill(HIST("TOFSignalPidTrue"), track.p(), track.beta()); - histReg.fill(HIST(TOFPidTrueRegistryNames[i]), track.p(), track.beta()); + fillPidHistos(track, pdgCode, true); } else { - isPidFalse = true; + fillPidHistos(track, pdgCode, false); } } } - - if (isPidFalse) { - histReg.fill(HIST(pidFalseRegistryNames[i]), track.pt()); - histReg.fill(HIST(TPCPidFalseRegistryNames[i]), track.p(), track.tpcSignal()); - histReg.fill(HIST("TPCSignalPidFalse"), track.p(), track.tpcSignal()); - histReg.fill(HIST("TOFSignalPidFalse"), track.p(), track.beta()); - histReg.fill(HIST(TOFPidFalseRegistryNames[i]), track.p(), track.beta()); - - double pt = track.pt(); - fillContaminationRegistry(i, pdgCode, pt); - } } template @@ -414,8 +394,6 @@ struct pid { { const float p = track.p(); - bool isPidFalse = false; - // list of Nsigmas for particles float particleNSigma[arrLen]; @@ -441,47 +419,53 @@ struct pid { float tmp_NSigma = abs(tpcNSigmas[i]); if ((tmp_NSigma < nsigmacut.value) & (indexOfSmallestElement(particleNSigma, arrLen) == i)) { if (pdgCode == pdgCodes[i]) { - histReg.fill(HIST(pidTrueRegistryNames[i]), track.pt()); - histReg.fill(HIST(TPCPidTrueRegistryNames[i]), track.p(), track.tpcSignal()); - histReg.fill(HIST("TPCSignalPidTrue"), track.p(), track.tpcSignal()); - histReg.fill(HIST("TOFSignalPidTrue"), track.p(), track.beta()); - histReg.fill(HIST(TOFPidTrueRegistryNames[i]), track.p(), track.beta()); + fillPidHistos(track, pdgCode, true); } else { - isPidFalse = true; + fillPidHistos(track, pdgCode, false); } } } else if ((p >= pSwitch[i]) & (track.sign() == particleCharge[i])) { float tmp_NSigma = combinedSignal(tpcNSigmas[i], tofNSigmas[i]); if ((tmp_NSigma < nsigmacut.value) & (indexOfSmallestElement(particleNSigma, arrLen) == i)) { if (pdgCode == pdgCodes[i]) { - histReg.fill(HIST(pidTrueRegistryNames[i]), track.pt()); - histReg.fill(HIST(TPCPidTrueRegistryNames[i]), track.p(), track.tpcSignal()); - histReg.fill(HIST("TPCSignalPidTrue"), track.p(), track.tpcSignal()); - histReg.fill(HIST("TOFSignalPidTrue"), track.p(), track.beta()); - histReg.fill(HIST(TOFPidTrueRegistryNames[i]), track.p(), track.beta()); + fillPidHistos(track, pdgCode, true); } else { - isPidFalse = true; + fillPidHistos(track, pdgCode, false); } } } - if (isPidFalse) { - histReg.fill(HIST(pidFalseRegistryNames[i]), track.pt()); - histReg.fill(HIST(TPCPidFalseRegistryNames[i]), track.p(), track.tpcSignal()); - histReg.fill(HIST("TPCSignalPidFalse"), track.p(), track.tpcSignal()); - histReg.fill(HIST("TOFSignalPidFalse"), track.p(), track.beta()); - histReg.fill(HIST(TOFPidFalseRegistryNames[i]), track.p(), track.beta()); + } + } + + template + int getPdgFromPid(const T& track) + { + // Convert PID to PDG code + // We don't identify particles with PID bigger than Proton, putting dummy code so they will be skipped. + int bayesPdg = track.bayesID() > o2::track::PID::Proton ? 999 : pidToPdg[track.bayesID()]; + if (track.bayesID() <= o2::track::PID::Proton && track.sign() == -1 * particleCharge[track.bayesID()]) { // Check if antiparticle + bayesPdg += numParticles / 2; + } + return bayesPdg; + } - double pt = track.pt(); - fillContaminationRegistry(i, pdgCode, pt); + template + void pidBayes(const T& track, const int pdgCode, const int bayesPdg) + { + if (bayesPdg == pdgCode) { + if (pdgCode == pdgCodes[i]) { + fillPidHistos(track, pdgCode, true); + } else { + fillPidHistos(track, pdgCode, false); } } } Configurable nsigmacut{"nsigmacut", 2.5, "Value of the NSigma cut"}; - Configurable strategy{"Strategy", 1, "1-PID with Nsigma method, 2-PID with NSigma and condition for minimal Nsigma value for particle, 3-Exlcusive condition for NSigma"}; + Configurable strategy{"strategy", 1, "1-PID with Nsigma method, 2-PID with NSigma and condition for minimal Nsigma value for particle, 3-Exlcusive condition for NSigma, 4-Bayesian PID"}; Filter trackFilter = aod::track::isGlobalTrack == static_cast(true); - using pidTracks = soa::Filtered>; + using pidTracks = soa::Filtered>; void process(pidTracks const& tracks, aod::McParticles const& mcParticles) { for (auto& track : tracks) { @@ -542,6 +526,11 @@ struct pid { static_for<0, 9>([&](auto i) { pidExclusiveStrategy(track, pdgCode, tpcNSigmas, tofNSigmas, 5); }); + } else if (strategy.value == 4) { + int bayesPdg = getPdgFromPid(track); + static_for<0, 9>([&](auto i) { + pidBayes(track, pdgCode, bayesPdg); + }); } } } From c3a5752aa8a94c810e5034723af6a2d03ce1cd2c Mon Sep 17 00:00:00 2001 From: saganatt <8majak8@gmail.com> Date: Mon, 24 Jan 2022 10:05:14 +0100 Subject: [PATCH 2/4] Fully working Bayes PID, uncertain results --- Tools/PIDML/qaPid.cxx | 102 +++++++++++++++++++++++++----------------- 1 file changed, 61 insertions(+), 41 deletions(-) diff --git a/Tools/PIDML/qaPid.cxx b/Tools/PIDML/qaPid.cxx index 6dcea38b07c..c256148f311 100644 --- a/Tools/PIDML/qaPid.cxx +++ b/Tools/PIDML/qaPid.cxx @@ -25,7 +25,7 @@ using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; -struct pid { +struct EvaluatePid { double combinedSignal(float val1, float val2) { return sqrt(pow(val1, 2) + pow(val2, 2)); @@ -461,50 +461,58 @@ struct pid { } } + template + void fillMcHistos(const T& track, const int pdgCode) + { + // pions + if (pdgCode == 211) { + histReg.fill(HIST("MC/211"), track.pt()); + } else if (pdgCode == -211) { + histReg.fill(HIST("MC/0211"), track.pt()); + } + // protons + else if (pdgCode == 2212) { + histReg.fill(HIST("MC/2212"), track.pt()); + } else if (pdgCode == -2212) { + histReg.fill(HIST("MC/02212"), track.pt()); + } + // kaons + else if (pdgCode == 321) { + histReg.fill(HIST("MC/321"), track.pt()); + } else if (pdgCode == -321) { + histReg.fill(HIST("MC/0321"), track.pt()); + } + // electrons + else if (pdgCode == 11) { + histReg.fill(HIST("MC/11"), track.pt()); + ; + } else if (pdgCode == -11) { + histReg.fill(HIST("MC/011"), track.pt()); + } + // muons + else if (pdgCode == 13) { + histReg.fill(HIST("MC/13"), track.pt()); + } else if (pdgCode == -13) { + histReg.fill(HIST("MC/013"), track.pt()); + } else { + histReg.fill(HIST("MC/else"), track.pt()); + } + } + Configurable nsigmacut{"nsigmacut", 2.5, "Value of the NSigma cut"}; Configurable strategy{"strategy", 1, "1-PID with Nsigma method, 2-PID with NSigma and condition for minimal Nsigma value for particle, 3-Exlcusive condition for NSigma, 4-Bayesian PID"}; Filter trackFilter = aod::track::isGlobalTrack == static_cast(true); - using pidTracks = soa::Filtered>; + using pidTracks = soa::Filtered>; + using pidBayesTracks = soa::Filtered>; + void process(pidTracks const& tracks, aod::McParticles const& mcParticles) { for (auto& track : tracks) { auto particle = track.mcParticle(); int pdgCode = particle.pdgCode(); - // pions - if (pdgCode == 211) { - histReg.fill(HIST("MC/211"), track.pt()); - } else if (pdgCode == -211) { - histReg.fill(HIST("MC/0211"), track.pt()); - } - // protons - else if (pdgCode == 2212) { - histReg.fill(HIST("MC/2212"), track.pt()); - } else if (pdgCode == -2212) { - histReg.fill(HIST("MC/02212"), track.pt()); - } - // kaons - else if (pdgCode == 321) { - histReg.fill(HIST("MC/321"), track.pt()); - } else if (pdgCode == -321) { - histReg.fill(HIST("MC/0321"), track.pt()); - } - // electrons - else if (pdgCode == 11) { - histReg.fill(HIST("MC/11"), track.pt()); - ; - } else if (pdgCode == -11) { - histReg.fill(HIST("MC/011"), track.pt()); - } - // muons - else if (pdgCode == 13) { - histReg.fill(HIST("MC/13"), track.pt()); - } else if (pdgCode == -13) { - histReg.fill(HIST("MC/013"), track.pt()); - } else { - histReg.fill(HIST("MC/else"), track.pt()); - } + fillMcHistos(track, pdgCode); // PID // indexes for particles: 0-Pi, 1-Pr, 2-Ka, 3-El, 4-Mu @@ -526,19 +534,31 @@ struct pid { static_for<0, 9>([&](auto i) { pidExclusiveStrategy(track, pdgCode, tpcNSigmas, tofNSigmas, 5); }); - } else if (strategy.value == 4) { - int bayesPdg = getPdgFromPid(track); - static_for<0, 9>([&](auto i) { - pidBayes(track, pdgCode, bayesPdg); - }); } } } + PROCESS_SWITCH(EvaluatePid, process, "Process standard strategies", true); + + void processBayes(pidBayesTracks const& tracks, aod::McParticles const& mcParticles) + { + for (auto& track : tracks) { + auto particle = track.mcParticle(); + int pdgCode = particle.pdgCode(); + + fillMcHistos(track, pdgCode); + + int bayesPdg = getPdgFromPid(track); + static_for<0, 9>([&](auto i) { + pidBayes(track, pdgCode, bayesPdg); + }); + } + } + PROCESS_SWITCH(EvaluatePid, processBayes, "Process Bayesian strategy", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ - adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc), }; } From 4c8dd799da4618a1e7fe05ec8de4482f92711571 Mon Sep 17 00:00:00 2001 From: saganatt <8majak8@gmail.com> Date: Mon, 24 Jan 2022 18:03:42 +0100 Subject: [PATCH 3/4] Renamed process --- Tools/PIDML/qaPid.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Tools/PIDML/qaPid.cxx b/Tools/PIDML/qaPid.cxx index c256148f311..54e40b373fd 100644 --- a/Tools/PIDML/qaPid.cxx +++ b/Tools/PIDML/qaPid.cxx @@ -506,7 +506,7 @@ struct EvaluatePid { using pidTracks = soa::Filtered>; using pidBayesTracks = soa::Filtered>; - void process(pidTracks const& tracks, aod::McParticles const& mcParticles) + void processStandard(pidTracks const& tracks, aod::McParticles const& mcParticles) { for (auto& track : tracks) { auto particle = track.mcParticle(); @@ -537,7 +537,7 @@ struct EvaluatePid { } } } - PROCESS_SWITCH(EvaluatePid, process, "Process standard strategies", true); + PROCESS_SWITCH(EvaluatePid, processStandard, "Process standard strategies", true); void processBayes(pidBayesTracks const& tracks, aod::McParticles const& mcParticles) { From d6c8b459285fc6e96f551915c7494e5fa3a47728 Mon Sep 17 00:00:00 2001 From: saganatt <8majak8@gmail.com> Date: Tue, 25 Jan 2022 16:56:50 +0100 Subject: [PATCH 4/4] Buggy process switches, back to single process --- Tools/PIDML/qaPid.cxx | 27 +++++++-------------------- 1 file changed, 7 insertions(+), 20 deletions(-) diff --git a/Tools/PIDML/qaPid.cxx b/Tools/PIDML/qaPid.cxx index 54e40b373fd..c152927dcd1 100644 --- a/Tools/PIDML/qaPid.cxx +++ b/Tools/PIDML/qaPid.cxx @@ -503,10 +503,9 @@ struct EvaluatePid { Configurable strategy{"strategy", 1, "1-PID with Nsigma method, 2-PID with NSigma and condition for minimal Nsigma value for particle, 3-Exlcusive condition for NSigma, 4-Bayesian PID"}; Filter trackFilter = aod::track::isGlobalTrack == static_cast(true); - using pidTracks = soa::Filtered>; - using pidBayesTracks = soa::Filtered>; + using pidTracks = soa::Filtered>; - void processStandard(pidTracks const& tracks, aod::McParticles const& mcParticles) + void process(pidTracks const& tracks, aod::McParticles const& mcParticles) { for (auto& track : tracks) { auto particle = track.mcParticle(); @@ -534,26 +533,14 @@ struct EvaluatePid { static_for<0, 9>([&](auto i) { pidExclusiveStrategy(track, pdgCode, tpcNSigmas, tofNSigmas, 5); }); + } else if (strategy.value == 4) { + int bayesPdg = getPdgFromPid(track); + static_for<0, 9>([&](auto i) { + pidBayes(track, pdgCode, bayesPdg); + }); } } } - PROCESS_SWITCH(EvaluatePid, processStandard, "Process standard strategies", true); - - void processBayes(pidBayesTracks const& tracks, aod::McParticles const& mcParticles) - { - for (auto& track : tracks) { - auto particle = track.mcParticle(); - int pdgCode = particle.pdgCode(); - - fillMcHistos(track, pdgCode); - - int bayesPdg = getPdgFromPid(track); - static_for<0, 9>([&](auto i) { - pidBayes(track, pdgCode, bayesPdg); - }); - } - } - PROCESS_SWITCH(EvaluatePid, processBayes, "Process Bayesian strategy", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)