diff --git a/ALICE3/DataModel/OTFRICH.h b/ALICE3/DataModel/OTFRICH.h index d4d9c5257ce..05771dda57b 100644 --- a/ALICE3/DataModel/OTFRICH.h +++ b/ALICE3/DataModel/OTFRICH.h @@ -31,12 +31,20 @@ DECLARE_SOA_COLUMN(NSigmaMuonRich, nSigmaMuonRich, float); //! NSigma mu DECLARE_SOA_COLUMN(NSigmaPionRich, nSigmaPionRich, float); //! NSigma pion BarrelRich DECLARE_SOA_COLUMN(NSigmaKaonRich, nSigmaKaonRich, float); //! NSigma kaon BarrelRich DECLARE_SOA_COLUMN(NSigmaProtonRich, nSigmaProtonRich, float); //! NSigma proton BarrelRich +DECLARE_SOA_COLUMN(NSigmaDeuteronRich, nSigmaDeuteronRich, float); //! NSigma deuteron BarrelRich +DECLARE_SOA_COLUMN(NSigmaTritonRich, nSigmaTritonRich, float); //! NSigma triton BarrelRich +DECLARE_SOA_COLUMN(NSigmaHelium3Rich, nSigmaHelium3Rich, float); //! NSigma helium3 BarrelRich +DECLARE_SOA_COLUMN(NSigmaAlphaRich, nSigmaAlphaRich, float); //! NSigma alpha BarrelRich DECLARE_SOA_DYNAMIC_COLUMN(NSigmaRich, nSigmaRich, //! General function to get the nSigma for the RICH [](const float el, const float mu, const float pi, const float ka, const float pr, + const float de, + const float tr, + const float he3, + const float al, const int id) -> float { switch (std::abs(id)) { case 0: @@ -49,6 +57,14 @@ DECLARE_SOA_DYNAMIC_COLUMN(NSigmaRich, nSigmaRich, //! General f return ka; case 4: return pr; + case 5: + return de; + case 6: + return tr; + case 7: + return he3; + case 8: + return al; default: LOG(fatal) << "Unrecognized PDG code for RICH"; return 999.f; @@ -62,6 +78,10 @@ DECLARE_SOA_COLUMN(HasSigMu, hasSigMu, bool); //! Has nSigma muon BarrelRi DECLARE_SOA_COLUMN(HasSigPi, hasSigPi, bool); //! Has nSigma pion BarrelRich (is pion over threshold) DECLARE_SOA_COLUMN(HasSigKa, hasSigKa, bool); //! Has nSigma kaon BarrelRich (is kaon over threshold) DECLARE_SOA_COLUMN(HasSigPr, hasSigPr, bool); //! Has nSigma proton BarrelRich (is proton over threshold) +DECLARE_SOA_COLUMN(HasSigDe, hasSigDe, bool); //! Has nSigma deuteron BarrelRich (is deuteron over threshold) +DECLARE_SOA_COLUMN(HasSigTr, hasSigTr, bool); //! Has nSigma triton BarrelRich (is triton over threshold) +DECLARE_SOA_COLUMN(HasSigHe3, hasSigHe3, bool); //! Has nSigma helium3 BarrelRich (is helium3 over threshold) +DECLARE_SOA_COLUMN(HasSigAl, hasSigAl, bool); //! Has nSigma alpha BarrelRich (is alpha over threshold) } // namespace upgrade_rich DECLARE_SOA_TABLE(UpgradeRichs, "AOD", "UPGRADERICH", @@ -70,11 +90,19 @@ DECLARE_SOA_TABLE(UpgradeRichs, "AOD", "UPGRADERICH", upgrade_rich::NSigmaPionRich, upgrade_rich::NSigmaKaonRich, upgrade_rich::NSigmaProtonRich, + upgrade_rich::NSigmaDeuteronRich, + upgrade_rich::NSigmaTritonRich, + upgrade_rich::NSigmaHelium3Rich, + upgrade_rich::NSigmaAlphaRich, upgrade_rich::NSigmaRich); + upgrade_rich::NSigmaProtonRich, + upgrade_rich::NSigmaDeuteronRich, + upgrade_rich::NSigmaTritonRich, + upgrade_rich::NSigmaHelium3Rich, + upgrade_rich::NSigmaAlphaRich>); using UpgradeRich = UpgradeRichs::iterator; @@ -85,6 +113,10 @@ DECLARE_SOA_TABLE(UpgradeRichSignals, "AOD", "UPGRADERICHSIG", upgrade_rich::HasSigPi, upgrade_rich::HasSigKa, upgrade_rich::HasSigPr, + upgrade_rich::HasSigDe, + upgrade_rich::HasSigTr, + upgrade_rich::HasSigHe3, + upgrade_rich::HasSigAl, upgrade_rich::HasSigInGas); using UpgradeRichSignal = UpgradeRichSignals::iterator; diff --git a/ALICE3/DataModel/OTFTOF.h b/ALICE3/DataModel/OTFTOF.h index f6f7c39234b..150efc91f5e 100644 --- a/ALICE3/DataModel/OTFTOF.h +++ b/ALICE3/DataModel/OTFTOF.h @@ -39,34 +39,54 @@ DECLARE_SOA_COLUMN(NSigmaMuonInnerTOF, nSigmaMuonInnerTOF, float); //! DECLARE_SOA_COLUMN(NSigmaPionInnerTOF, nSigmaPionInnerTOF, float); //! NSigma pion InnerTOF DECLARE_SOA_COLUMN(NSigmaKaonInnerTOF, nSigmaKaonInnerTOF, float); //! NSigma kaon InnerTOF DECLARE_SOA_COLUMN(NSigmaProtonInnerTOF, nSigmaProtonInnerTOF, float); //! NSigma proton InnerTOF +DECLARE_SOA_COLUMN(NSigmaDeuteronInnerTOF, nSigmaDeuteronInnerTOF, float); //! NSigma deuteron InnerTOF +DECLARE_SOA_COLUMN(NSigmaTritonInnerTOF, nSigmaTritonInnerTOF, float); //! NSigma triton InnerTOF +DECLARE_SOA_COLUMN(NSigmaHelium3InnerTOF, nSigmaHelium3InnerTOF, float); //! NSigma helium3 InnerTOF +DECLARE_SOA_COLUMN(NSigmaAlphaInnerTOF, nSigmaAlphaInnerTOF, float); //! NSigma alpha InnerTOF DECLARE_SOA_COLUMN(InnerTOFTrackTimeReco, innerTOFTrackTimeReco, float); //! Track time measured at the InnerTOF DECLARE_SOA_COLUMN(InnerTOFTrackLengthReco, innerTOFTrackLengthReco, float); //! track length for calculation of InnerTOF (reconstructed) -DECLARE_SOA_COLUMN(InnerTOFExpectedTimeEl, innerTOFExpectedTimeEl, float); //! Reconstructed expected time at the InnerTOF for the Electron mass hypotheses -DECLARE_SOA_COLUMN(InnerTOFExpectedTimeMu, innerTOFExpectedTimeMu, float); //! Reconstructed expected time at the InnerTOF for the Muon mass hypotheses -DECLARE_SOA_COLUMN(InnerTOFExpectedTimePi, innerTOFExpectedTimePi, float); //! Reconstructed expected time at the InnerTOF for the Pion mass hypotheses -DECLARE_SOA_COLUMN(InnerTOFExpectedTimeKa, innerTOFExpectedTimeKa, float); //! Reconstructed expected time at the InnerTOF for the Kaon mass hypotheses -DECLARE_SOA_COLUMN(InnerTOFExpectedTimePr, innerTOFExpectedTimePr, float); //! Reconstructed expected time at the InnerTOF for the Proton mass hypotheses +DECLARE_SOA_COLUMN(InnerTOFExpectedTimeEl, innerTOFExpectedTimeEl, float); //! Reconstructed expected time at the InnerTOF for the Electron mass hypotheses +DECLARE_SOA_COLUMN(InnerTOFExpectedTimeMu, innerTOFExpectedTimeMu, float); //! Reconstructed expected time at the InnerTOF for the Muon mass hypotheses +DECLARE_SOA_COLUMN(InnerTOFExpectedTimePi, innerTOFExpectedTimePi, float); //! Reconstructed expected time at the InnerTOF for the Pion mass hypotheses +DECLARE_SOA_COLUMN(InnerTOFExpectedTimeKa, innerTOFExpectedTimeKa, float); //! Reconstructed expected time at the InnerTOF for the Kaon mass hypotheses +DECLARE_SOA_COLUMN(InnerTOFExpectedTimePr, innerTOFExpectedTimePr, float); //! Reconstructed expected time at the InnerTOF for the Proton mass hypotheses +DECLARE_SOA_COLUMN(InnerTOFExpectedTimeDe, innerTOFExpectedTimeDe, float); //! Reconstructed expected time at the InnerTOF for the Deuteron mass hypotheses +DECLARE_SOA_COLUMN(InnerTOFExpectedTimeTr, innerTOFExpectedTimeTr, float); //! Reconstructed expected time at the InnerTOF for the Triton mass hypotheses +DECLARE_SOA_COLUMN(InnerTOFExpectedTimeHe3, innerTOFExpectedTimeHe3, float); //! Reconstructed expected time at the InnerTOF for the Helium3 mass hypotheses +DECLARE_SOA_COLUMN(InnerTOFExpectedTimeAl, innerTOFExpectedTimeAl, float); //! Reconstructed expected time at the InnerTOF for the Alpha mass hypotheses DECLARE_SOA_COLUMN(NSigmaElectronOuterTOF, nSigmaElectronOuterTOF, float); //! NSigma electron OuterTOF DECLARE_SOA_COLUMN(NSigmaMuonOuterTOF, nSigmaMuonOuterTOF, float); //! NSigma muon OuterTOF DECLARE_SOA_COLUMN(NSigmaPionOuterTOF, nSigmaPionOuterTOF, float); //! NSigma pion OuterTOF DECLARE_SOA_COLUMN(NSigmaKaonOuterTOF, nSigmaKaonOuterTOF, float); //! NSigma kaon OuterTOF DECLARE_SOA_COLUMN(NSigmaProtonOuterTOF, nSigmaProtonOuterTOF, float); //! NSigma proton OuterTOF +DECLARE_SOA_COLUMN(NSigmaDeuteronOuterTOF, nSigmaDeuteronOuterTOF, float); //! NSigma deuteron OuterTOF +DECLARE_SOA_COLUMN(NSigmaTritonOuterTOF, nSigmaTritonOuterTOF, float); //! NSigma triton OuterTOF +DECLARE_SOA_COLUMN(NSigmaHelium3OuterTOF, nSigmaHelium3OuterTOF, float); //! NSigma helium3 OuterTOF +DECLARE_SOA_COLUMN(NSigmaAlphaOuterTOF, nSigmaAlphaOuterTOF, float); //! NSigma alpha OuterTOF DECLARE_SOA_COLUMN(OuterTOFTrackTimeReco, outerTOFTrackTimeReco, float); //! Track time measured at the OuterTOF DECLARE_SOA_COLUMN(OuterTOFTrackLengthReco, outerTOFTrackLengthReco, float); //! track length for calculation of OuterTOF (reconstructed) -DECLARE_SOA_COLUMN(OuterTOFExpectedTimeEl, outerTOFExpectedTimeEl, float); //! Reconstructed expected time at the OuterTOF for the Electron mass hypotheses -DECLARE_SOA_COLUMN(OuterTOFExpectedTimeMu, outerTOFExpectedTimeMu, float); //! Reconstructed expected time at the OuterTOF for the Muon mass hypotheses -DECLARE_SOA_COLUMN(OuterTOFExpectedTimePi, outerTOFExpectedTimePi, float); //! Reconstructed expected time at the OuterTOF for the Pion mass hypotheses -DECLARE_SOA_COLUMN(OuterTOFExpectedTimeKa, outerTOFExpectedTimeKa, float); //! Reconstructed expected time at the OuterTOF for the Kaon mass hypotheses -DECLARE_SOA_COLUMN(OuterTOFExpectedTimePr, outerTOFExpectedTimePr, float); //! Reconstructed expected time at the OuterTOF for the Proton mass hypotheses -DECLARE_SOA_DYNAMIC_COLUMN(NSigmaInnerTOF, nSigmaInnerTOF, //! General function to get the nSigma for the InnerTOF +DECLARE_SOA_COLUMN(OuterTOFExpectedTimeEl, outerTOFExpectedTimeEl, float); //! Reconstructed expected time at the OuterTOF for the Electron mass hypotheses +DECLARE_SOA_COLUMN(OuterTOFExpectedTimeMu, outerTOFExpectedTimeMu, float); //! Reconstructed expected time at the OuterTOF for the Muon mass hypotheses +DECLARE_SOA_COLUMN(OuterTOFExpectedTimePi, outerTOFExpectedTimePi, float); //! Reconstructed expected time at the OuterTOF for the Pion mass hypotheses +DECLARE_SOA_COLUMN(OuterTOFExpectedTimeKa, outerTOFExpectedTimeKa, float); //! Reconstructed expected time at the OuterTOF for the Kaon mass hypotheses +DECLARE_SOA_COLUMN(OuterTOFExpectedTimePr, outerTOFExpectedTimePr, float); //! Reconstructed expected time at the OuterTOF for the Proton mass hypotheses +DECLARE_SOA_COLUMN(OuterTOFExpectedTimeDe, outerTOFExpectedTimeDe, float); //! Reconstructed expected time at the OuterTOF for the Deuteron mass hypotheses +DECLARE_SOA_COLUMN(OuterTOFExpectedTimeTr, outerTOFExpectedTimeTr, float); //! Reconstructed expected time at the OuterTOF for the Triton mass hypotheses +DECLARE_SOA_COLUMN(OuterTOFExpectedTimeHe3, outerTOFExpectedTimeHe3, float); //! Reconstructed expected time at the OuterTOF for the Helium3 mass hypotheses +DECLARE_SOA_COLUMN(OuterTOFExpectedTimeAl, outerTOFExpectedTimeAl, float); //! Reconstructed expected time at the OuterTOF for the Alpha mass hypotheses +DECLARE_SOA_DYNAMIC_COLUMN(NSigmaInnerTOF, nSigmaInnerTOF, //! General function to get the nSigma for the InnerTOF [](const float el, const float mu, const float pi, const float ka, const float pr, + const float de, + const float tr, + const float he3, + const float al, const int id) -> float { switch (std::abs(id)) { case 0: @@ -79,6 +99,14 @@ DECLARE_SOA_DYNAMIC_COLUMN(NSigmaInnerTOF, nSigmaInnerTOF, //! G return ka; case 4: return pr; + case 5: + return de; + case 6: + return tr; + case 7: + return he3; + case 8: + return al; default: LOG(fatal) << "Unrecognized PDG code for InnerTOF"; return 999.f; @@ -90,6 +118,10 @@ DECLARE_SOA_DYNAMIC_COLUMN(NSigmaOuterTOF, nSigmaOuterTOF, //! General function const float pi, const float ka, const float pr, + const float de, + const float tr, + const float he3, + const float al, const int id) -> float { switch (std::abs(id)) { case 0: @@ -102,6 +134,14 @@ DECLARE_SOA_DYNAMIC_COLUMN(NSigmaOuterTOF, nSigmaOuterTOF, //! General function return ka; case 4: return pr; + case 5: + return de; + case 6: + return tr; + case 7: + return he3; + case 8: + return al; default: LOG(fatal) << "Unrecognized PDG code for InnerTOF"; return 999.f; @@ -124,6 +164,10 @@ DECLARE_SOA_TABLE(UpgradeTofs, "AOD", "UPGRADETOF", upgrade_tof::NSigmaPionInnerTOF, upgrade_tof::NSigmaKaonInnerTOF, upgrade_tof::NSigmaProtonInnerTOF, + upgrade_tof::NSigmaDeuteronInnerTOF, + upgrade_tof::NSigmaTritonInnerTOF, + upgrade_tof::NSigmaHelium3InnerTOF, + upgrade_tof::NSigmaAlphaInnerTOF, upgrade_tof::InnerTOFTrackTimeReco, upgrade_tof::InnerTOFTrackLengthReco, upgrade_tof::NSigmaElectronOuterTOF, @@ -131,18 +175,30 @@ DECLARE_SOA_TABLE(UpgradeTofs, "AOD", "UPGRADETOF", upgrade_tof::NSigmaPionOuterTOF, upgrade_tof::NSigmaKaonOuterTOF, upgrade_tof::NSigmaProtonOuterTOF, + upgrade_tof::NSigmaDeuteronOuterTOF, + upgrade_tof::NSigmaTritonOuterTOF, + upgrade_tof::NSigmaHelium3OuterTOF, + upgrade_tof::NSigmaAlphaOuterTOF, upgrade_tof::OuterTOFTrackTimeReco, upgrade_tof::OuterTOFTrackLengthReco, upgrade_tof::NSigmaInnerTOF, + upgrade_tof::NSigmaProtonInnerTOF, + upgrade_tof::NSigmaDeuteronInnerTOF, + upgrade_tof::NSigmaTritonInnerTOF, + upgrade_tof::NSigmaHelium3InnerTOF, + upgrade_tof::NSigmaAlphaInnerTOF>, upgrade_tof::NSigmaOuterTOF); + upgrade_tof::NSigmaProtonOuterTOF, + upgrade_tof::NSigmaDeuteronOuterTOF, + upgrade_tof::NSigmaTritonOuterTOF, + upgrade_tof::NSigmaHelium3OuterTOF, + upgrade_tof::NSigmaAlphaOuterTOF>); DECLARE_SOA_TABLE(UpgradeTofExpectedTimes, "AOD", "UPGRADETOFEXPT", upgrade_tof::InnerTOFExpectedTimeEl, @@ -150,11 +206,19 @@ DECLARE_SOA_TABLE(UpgradeTofExpectedTimes, "AOD", "UPGRADETOFEXPT", upgrade_tof::InnerTOFExpectedTimePi, upgrade_tof::InnerTOFExpectedTimeKa, upgrade_tof::InnerTOFExpectedTimePr, + upgrade_tof::InnerTOFExpectedTimeDe, + upgrade_tof::InnerTOFExpectedTimeTr, + upgrade_tof::InnerTOFExpectedTimeHe3, + upgrade_tof::InnerTOFExpectedTimeAl, upgrade_tof::OuterTOFExpectedTimeEl, upgrade_tof::OuterTOFExpectedTimeMu, upgrade_tof::OuterTOFExpectedTimePi, upgrade_tof::OuterTOFExpectedTimeKa, - upgrade_tof::OuterTOFExpectedTimePr); + upgrade_tof::OuterTOFExpectedTimePr, + upgrade_tof::OuterTOFExpectedTimeDe, + upgrade_tof::OuterTOFExpectedTimeTr, + upgrade_tof::OuterTOFExpectedTimeHe3, + upgrade_tof::OuterTOFExpectedTimeAl); using UpgradeTofMC = UpgradeTofMCs::iterator; using UpgradeTof = UpgradeTofs::iterator; diff --git a/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx index 9ccc09b3a5d..bbb903b307c 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx @@ -319,6 +319,7 @@ struct OnTheFlyRichPid { loadLUT(1000010020, "lutDe"); loadLUT(1000010030, "lutTr"); loadLUT(1000020030, "lutHe3"); + loadLUT(1000020040, "lutAl"); } if (doQAplots) { @@ -333,9 +334,9 @@ struct OnTheFlyRichPid { histos.add("h2dAngularResolutionVsEtaBarrelRICH", "h2dAngularResolutionVsEtaBarrelRICH", kTH2F, {axisEta, axisRingAngularResolution}); histos.add("hSectorID", "hSectorID", kTH1F, {axisSector}); - const int kNspec = 5; // electron, muon, pion, kaon, proton - std::string particleNames1[kNspec] = {"#it{e}", "#it{#mu}", "#it{#pi}", "#it{K}", "#it{p}"}; - std::string particleNames2[kNspec] = {"Elec", "Muon", "Pion", "Kaon", "Prot"}; + const int kNspec = 9; // electron, muon, pion, kaon, proton, deuteron, triton, helium3, alpha + std::string particleNames1[kNspec] = {"#it{e}", "#it{#mu}", "#it{#pi}", "#it{K}", "#it{p}", "#it{d}", "#it{t}", "^{3}He", "#it{#alpha}"}; + std::string particleNames2[kNspec] = {"Elec", "Muon", "Pion", "Kaon", "Prot", "Deut", "Trit", "He3", "Al"}; for (int iTrue = 0; iTrue < kNspec; iTrue++) { std::string nameTitleBarrelTrackRes = "h2dBarrelAngularResTrack" + particleNames2[iTrue] + "VsP"; std::string nameTitleBarrelTotalRes = "h2dBarrelAngularResTotal" + particleNames2[iTrue] + "VsP"; @@ -735,8 +736,8 @@ struct OnTheFlyRichPid { for (const auto& track : tracks) { auto fillDummyValues = [&](bool gasRich = false) { - upgradeRich(kErrorValue, kErrorValue, kErrorValue, kErrorValue, kErrorValue); - upgradeRichSignal(false, false, false, false, false, false, gasRich); + upgradeRich(kErrorValue, kErrorValue, kErrorValue, kErrorValue, kErrorValue, kErrorValue, kErrorValue, kErrorValue, kErrorValue); + upgradeRichSignal(false, false, false, false, false, false, false, false, false, false, gasRich); }; // first step: find precise arrival time (if any) @@ -798,21 +799,29 @@ struct OnTheFlyRichPid { } // Straight to Nsigma - static constexpr int kNspecies = 5; + static constexpr int kNspecies = 9; static constexpr int kEl = 0; static constexpr int kMu = 1; static constexpr int kPi = 2; static constexpr int kKa = 3; static constexpr int kPr = 4; - float nSigmaBarrelRich[kNspecies] = {kErrorValue, kErrorValue, kErrorValue, kErrorValue, kErrorValue}; - bool signalBarrelRich[kNspecies] = {false, false, false, false, false}; + static constexpr int kDe = 5; + static constexpr int kTr = 6; + static constexpr int kHe3 = 7; + static constexpr int kAl = 8; + float nSigmaBarrelRich[kNspecies] = {kErrorValue, kErrorValue, kErrorValue, kErrorValue, kErrorValue, kErrorValue, kErrorValue, kErrorValue, kErrorValue}; + bool signalBarrelRich[kNspecies] = {false, false, false, false, false, false, false, false, false}; float deltaThetaBarrelRich[kNspecies]; //, nSigmaBarrelRich[kNspecies]; - static constexpr int kPdgArray[kNspecies] = {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton}; + static constexpr int kPdgArray[kNspecies] = {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton, o2::constants::physics::kDeuteron, o2::constants::physics::kTriton, o2::constants::physics::kHelium3, o2::constants::physics::kAlpha}; static constexpr float kMasses[kNspecies] = {o2::track::pid_constants::sMasses[o2::track::PID::Electron], o2::track::pid_constants::sMasses[o2::track::PID::Muon], o2::track::pid_constants::sMasses[o2::track::PID::Pion], o2::track::pid_constants::sMasses[o2::track::PID::Kaon], - o2::track::pid_constants::sMasses[o2::track::PID::Proton]}; + o2::track::pid_constants::sMasses[o2::track::PID::Proton], + o2::track::pid_constants::sMasses[o2::track::PID::Deuteron], + o2::track::pid_constants::sMasses[o2::track::PID::Triton], + o2::track::pid_constants::sMasses[o2::track::PID::Helium3], + o2::track::pid_constants::sMasses[o2::track::PID::Alpha]}; for (int ii = 0; ii < kNspecies; ii++) { // Loop on the particle hypotheses @@ -874,6 +883,34 @@ struct OnTheFlyRichPid { histos.fill(HIST("h2dBarrelAngularResTotalProtVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); } break; + case kPdgArray[kDe]: // Deuteron + case -kPdgArray[kDe]: // AntiDeuteron + if (ii == kDe) { + histos.fill(HIST("h2dBarrelAngularResTrackDeutVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); + histos.fill(HIST("h2dBarrelAngularResTotalDeutVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); + } + break; + case kPdgArray[kTr]: // Triton + case -kPdgArray[kTr]: // AntiTriton + if (ii == kTr) { + histos.fill(HIST("h2dBarrelAngularResTrackTritVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); + histos.fill(HIST("h2dBarrelAngularResTotalTritVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); + } + break; + case kPdgArray[kHe3]: // Helium3 + case -kPdgArray[kHe3]: // AntiHelium3 + if (ii == kHe3) { + histos.fill(HIST("h2dBarrelAngularResTrackHe3VsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); + histos.fill(HIST("h2dBarrelAngularResTotalHe3VsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); + } + break; + case kPdgArray[kAl]: // Alpha + case -kPdgArray[kAl]: // AntiAlpha + if (ii == kAl) { + histos.fill(HIST("h2dBarrelAngularResTrackAlVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); + histos.fill(HIST("h2dBarrelAngularResTotalAlVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); + } + break; default: break; } @@ -943,6 +980,37 @@ struct OnTheFlyRichPid { histos.fill(HIST("h2dBarrelNsigmaTrueProtVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]); histos.fill(HIST("h2dBarrelNsigmaTrueProtVsKaonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[3]); histos.fill(HIST("h2dBarrelNsigmaTrueProtVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]); + histos.fill(HIST("h2dBarrelNsigmaTrueProtVsDeutHypothesis"), recoTrack.getP(), nSigmaBarrelRich[5]); + break; + case kPdgArray[kDe]: // Deuteron + case -kPdgArray[kDe]: // AntiDeuteron + histos.fill(HIST("h2dBarrelNsigmaTrueDeutVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]); + histos.fill(HIST("h2dBarrelNsigmaTrueDeutVsDeutHypothesis"), recoTrack.getP(), nSigmaBarrelRich[5]); + histos.fill(HIST("h2dBarrelNsigmaTrueDeutVsTritHypothesis"), recoTrack.getP(), nSigmaBarrelRich[6]); + histos.fill(HIST("h2dBarrelNsigmaTrueDeutVsHe3Hypothesis"), recoTrack.getP(), nSigmaBarrelRich[7]); + histos.fill(HIST("h2dBarrelNsigmaTrueDeutVsAlHypothesis"), recoTrack.getP(), nSigmaBarrelRich[8]); + break; + case kPdgArray[kTr]: // Triton + case -kPdgArray[kTr]: // AntiTriton + histos.fill(HIST("h2dBarrelNsigmaTrueTritVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]); + histos.fill(HIST("h2dBarrelNsigmaTrueTritVsDeutHypothesis"), recoTrack.getP(), nSigmaBarrelRich[5]); + histos.fill(HIST("h2dBarrelNsigmaTrueTritVsTritHypothesis"), recoTrack.getP(), nSigmaBarrelRich[6]); + histos.fill(HIST("h2dBarrelNsigmaTrueTritVsHe3Hypothesis"), recoTrack.getP(), nSigmaBarrelRich[7]); + histos.fill(HIST("h2dBarrelNsigmaTrueTritVsAlHypothesis"), recoTrack.getP(), nSigmaBarrelRich[8]); + break; + case kPdgArray[kHe3]: // Helium3 + case -kPdgArray[kHe3]: // AntiHelium3 + histos.fill(HIST("h2dBarrelNsigmaTrueHe3VsDeutHypothesis"), recoTrack.getP(), nSigmaBarrelRich[5]); + histos.fill(HIST("h2dBarrelNsigmaTrueHe3VsTritHypothesis"), recoTrack.getP(), nSigmaBarrelRich[6]); + histos.fill(HIST("h2dBarrelNsigmaTrueHe3VsHe3Hypothesis"), recoTrack.getP(), nSigmaBarrelRich[7]); + histos.fill(HIST("h2dBarrelNsigmaTrueHe3VsAlHypothesis"), recoTrack.getP(), nSigmaBarrelRich[8]); + break; + case kPdgArray[kAl]: // Alpha + case -kPdgArray[kAl]: // AntiAlpha + histos.fill(HIST("h2dBarrelNsigmaTrueAlVsDeutHypothesis"), recoTrack.getP(), nSigmaBarrelRich[5]); + histos.fill(HIST("h2dBarrelNsigmaTrueAlVsTritHypothesis"), recoTrack.getP(), nSigmaBarrelRich[6]); + histos.fill(HIST("h2dBarrelNsigmaTrueAlVsHe3Hypothesis"), recoTrack.getP(), nSigmaBarrelRich[7]); + histos.fill(HIST("h2dBarrelNsigmaTrueAlVsAlHypothesis"), recoTrack.getP(), nSigmaBarrelRich[8]); break; default: break; @@ -951,8 +1019,8 @@ struct OnTheFlyRichPid { } // Sigmas have been fully calculated. Please populate the NSigma helper table (once per track) - upgradeRich(nSigmaBarrelRich[0], nSigmaBarrelRich[1], nSigmaBarrelRich[2], nSigmaBarrelRich[3], nSigmaBarrelRich[4]); - upgradeRichSignal(expectedAngleBarrelRichOk, signalBarrelRich[0], signalBarrelRich[1], signalBarrelRich[2], signalBarrelRich[3], signalBarrelRich[4], expectedAngleBarrelGasRichOk); + upgradeRich(nSigmaBarrelRich[0], nSigmaBarrelRich[1], nSigmaBarrelRich[2], nSigmaBarrelRich[3], nSigmaBarrelRich[4], nSigmaBarrelRich[5], nSigmaBarrelRich[6], nSigmaBarrelRich[7], nSigmaBarrelRich[8]); + upgradeRichSignal(expectedAngleBarrelRichOk, signalBarrelRich[0], signalBarrelRich[1], signalBarrelRich[2], signalBarrelRich[3], signalBarrelRich[4], signalBarrelRich[5], signalBarrelRich[6], signalBarrelRich[7], signalBarrelRich[8], expectedAngleBarrelGasRichOk); } } }; diff --git a/ALICE3/TableProducer/OTF/onTheFlyTofPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyTofPid.cxx index b28d48d92ff..2e4fbc5d05d 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTofPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTofPid.cxx @@ -65,14 +65,14 @@ using namespace o2; using namespace o2::framework; -std::array, 5> h2dInnerTimeResTrack; -std::array, 5> h2dInnerTimeResTotal; -std::array, 5> h2dOuterTimeResTrack; -std::array, 5> h2dOuterTimeResTotal; -std::array, 5>, 5> h2dInnerNsigmaTrue; -std::array, 5>, 5> h2dOuterNsigmaTrue; -std::array, 5>, 5> h2dInnerDeltaTrue; -std::array, 5>, 5> h2dOuterDeltaTrue; +std::array, 9> h2dInnerTimeResTrack; +std::array, 9> h2dInnerTimeResTotal; +std::array, 9> h2dOuterTimeResTrack; +std::array, 9> h2dOuterTimeResTotal; +std::array, 9>, 9> h2dInnerNsigmaTrue; +std::array, 9>, 9> h2dOuterNsigmaTrue; +std::array, 9>, 9> h2dInnerDeltaTrue; +std::array, 9>, 9> h2dOuterDeltaTrue; struct OnTheFlyTofPid { Produces upgradeTofMC; @@ -135,7 +135,7 @@ struct OnTheFlyTofPid { // for handling basic QA histograms if requested HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; OutputObj listEfficiency{"efficiency"}; - static constexpr int kParticles = 5; + static constexpr int kParticles = 9; void init(o2::framework::InitContext& initContext) { @@ -169,6 +169,7 @@ struct OnTheFlyTofPid { loadLUT(1000010020, "lutDe"); loadLUT(1000010030, "lutTr"); loadLUT(1000020030, "lutHe3"); + loadLUT(1000020040, "lutAl"); } if (plotsConfig.doQAplots) { @@ -184,17 +185,20 @@ struct OnTheFlyTofPid { listEfficiency->Add(new TEfficiency("effEventTime", "effEventTime", plotsConfig.nBinsMult, 0.0f, plotsConfig.maxMultRange)); const AxisSpec axisMomentum{static_cast(plotsConfig.nBinsP), 0.0f, +10.0f, "#it{p} (GeV/#it{c})"}; + const AxisSpec axisRigidity{static_cast(plotsConfig.nBinsP), 0.0f, +10.0f, "#it{p} / |#it{z}| (GeV/#it{c})"}; const AxisSpec axisMomentumSmall{static_cast(plotsConfig.nBinsP), 0.0f, +1.0f, "#it{p} (GeV/#it{c})"}; const AxisSpec axisVelocity{static_cast(plotsConfig.nBinsBeta), 0.0f, +1.1f, "Measured #beta"}; const AxisSpec axisTrackLengthInner{static_cast(plotsConfig.nBinsTrackLengthInner), 0.0f, 60.0f, "Track length (cm)"}; const AxisSpec axisTrackLengthOuter{static_cast(plotsConfig.nBinsTrackLengthOuter), 0.0f, 300.0f, "Track length (cm)"}; const AxisSpec axisTrackDeltaLength{static_cast(plotsConfig.nBinsTrackDeltaLength), 0.0f, 30.0f, "Delta Track length (cm)"}; histos.add("iTOF/h2dVelocityVsMomentumInner", "h2dVelocityVsMomentumInner", kTH2F, {axisMomentum, axisVelocity}); + histos.add("iTOF/h2dVelocityVsRigidityInner", "h2dVelocityVsRigidityInner", kTH2F, {axisRigidity, axisVelocity}); histos.add("iTOF/h2dTrackLengthInnerVsPt", "h2dTrackLengthInnerVsPt", kTH2F, {axisMomentumSmall, axisTrackLengthInner}); histos.add("iTOF/h2dTrackLengthInnerRecoVsPt", "h2dTrackLengthInnerRecoVsPt", kTH2F, {axisMomentumSmall, axisTrackLengthInner}); histos.add("iTOF/h2dDeltaTrackLengthInnerVsPt", "h2dDeltaTrackLengthInnerVsPt", kTH2F, {axisMomentumSmall, axisTrackDeltaLength}); histos.add("oTOF/h2dVelocityVsMomentumOuter", "h2dVelocityVsMomentumOuter", kTH2F, {axisMomentum, axisVelocity}); + histos.add("oTOF/h2dVelocityVsRigidityOuter", "h2dVelocityVsRigidityOuter", kTH2F, {axisRigidity, axisVelocity}); histos.add("oTOF/h2dTrackLengthOuterVsPt", "h2dTrackLengthOuterVsPt", kTH2F, {axisMomentumSmall, axisTrackLengthOuter}); histos.add("oTOF/h2dTrackLengthOuterRecoVsPt", "h2dTrackLengthOuterRecoVsPt", kTH2F, {axisMomentumSmall, axisTrackLengthOuter}); histos.add("oTOF/h2dDeltaTrackLengthOuterVsPt", "h2dDeltaTrackLengthOuterVsPt", kTH2F, {axisMomentumSmall, axisTrackDeltaLength}); @@ -206,8 +210,8 @@ struct OnTheFlyTofPid { histos.add("h2dRelativePtResolution", "h2dRelativePtResolution", kTH2F, {axisPt, axisRelativePt}); histos.add("h2dRelativeEtaResolution", "h2dRelativeEtaResolution", kTH2F, {axisEta, axisRelativeEta}); - std::string particleNames[kParticles] = {"#it{e}", "#it{#mu}", "#it{#pi}", "#it{K}", "#it{p}"}; - std::string particleNames2[kParticles] = {"Elec", "Muon", "Pion", "Kaon", "Prot"}; + std::string particleNames[kParticles] = {"#it{e}", "#it{#mu}", "#it{#pi}", "#it{K}", "#it{p}", "#it{d}", "#it{t}", "^{3}He", "#it{#alpha}"}; + std::string particleNames2[kParticles] = {"Elec", "Muon", "Pion", "Kaon", "Prot", "Deut", "Trit", "He3", "Al"}; for (int iTrue = 0; iTrue < kParticles; iTrue++) { auto addHistogram = [&](const std::string& name, const AxisSpec& axis) { return histos.add(name, "", kTH2F, {axisMomentum, axis}); @@ -579,8 +583,20 @@ struct OnTheFlyTofPid { static std::array expectedTimeInnerTOF, expectedTimeOuterTOF; static std::array deltaTimeInnerTOF, deltaTimeOuterTOF; static std::array nSigmaInnerTOF, nSigmaOuterTOF; - static constexpr int kParticlePdgs[kParticles] = {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton}; + static constexpr int kParticlePdgs[kParticles] = {kElectron, kMuonMinus, kPiPlus, kKPlus, kProton, o2::constants::physics::kDeuteron, o2::constants::physics::kTriton, o2::constants::physics::kHelium3, o2::constants::physics::kAlpha}; float masses[kParticles]; + float momentumHypotheses[kParticles]; // Store momentum hypothesis for each particle + + auto truePdgInfo = pdg->GetParticle(mcParticle.pdgCode()); + float rigidity = momentum; // fallback to momentum if charge unknown + + // Use MC truth charge for rigidity calculation + if (truePdgInfo) { + const float trueCharge = std::abs(truePdgInfo->Charge()) / 3.0f; + if (trueCharge > 0) { + rigidity = momentum / trueCharge; + } + } if (plotsConfig.doQAplots) { // unit conversion: length in cm, time in ps @@ -588,11 +604,13 @@ struct OnTheFlyTofPid { const float outerBeta = (trackLengthOuterTOF / measuredTimeOuterTOF) / o2::constants::physics::LightSpeedCm2PS; if (trackLengthRecoInnerTOF > 0) { histos.fill(HIST("iTOF/h2dVelocityVsMomentumInner"), momentum, innerBeta); + histos.fill(HIST("iTOF/h2dVelocityVsRigidityInner"), rigidity, innerBeta); histos.fill(HIST("iTOF/h2dTrackLengthInnerVsPt"), noSmearingPt, trackLengthInnerTOF); histos.fill(HIST("iTOF/h2dTrackLengthInnerRecoVsPt"), noSmearingPt, trackLengthRecoInnerTOF); } if (trackLengthRecoOuterTOF > 0) { histos.fill(HIST("oTOF/h2dVelocityVsMomentumOuter"), momentum, outerBeta); + histos.fill(HIST("oTOF/h2dVelocityVsRigidityOuter"), rigidity, outerBeta); histos.fill(HIST("oTOF/h2dTrackLengthOuterVsPt"), noSmearingPt, trackLengthOuterTOF); histos.fill(HIST("oTOF/h2dTrackLengthOuterRecoVsPt"), noSmearingPt, trackLengthRecoOuterTOF); } @@ -608,7 +626,8 @@ struct OnTheFlyTofPid { auto pdgInfoThis = pdg->GetParticle(kParticlePdgs[ii]); masses[ii] = pdgInfoThis->Mass(); - const float v = computeParticleVelocity(momentum, masses[ii]); + momentumHypotheses[ii] = rigidity * (std::abs(pdgInfoThis->Charge()) / 3.0f); // Total momentum for this hypothesis + const float v = computeParticleVelocity(momentumHypotheses[ii], masses[ii]); expectedTimeInnerTOF[ii] = trackLengthInnerTOF / v; expectedTimeOuterTOF[ii] = trackLengthOuterTOF / v; @@ -620,27 +639,27 @@ struct OnTheFlyTofPid { float innerTotalTimeReso = simConfig.innerTOFTimeReso; float outerTotalTimeReso = simConfig.outerTOFTimeReso; if (simConfig.flagIncludeTrackTimeRes) { - double ptResolution = std::pow(momentum / std::cosh(pseudorapidity), 2) * std::sqrt(trkWithTime.mMomentum.second); + double ptResolution = std::pow(momentumHypotheses[ii] / std::cosh(pseudorapidity), 2) * std::sqrt(trkWithTime.mMomentum.second); double etaResolution = std::fabs(std::sin(2.0 * std::atan(std::exp(-pseudorapidity)))) * std::sqrt(trkWithTime.mPseudorapidity.second); if (simConfig.flagTOFLoadDelphesLUTs) { - ptResolution = mSmearer.getAbsPtRes(pdgInfoThis->PdgCode(), dNdEta, pseudorapidity, momentum / std::cosh(pseudorapidity)); - etaResolution = mSmearer.getAbsEtaRes(pdgInfoThis->PdgCode(), dNdEta, pseudorapidity, momentum / std::cosh(pseudorapidity)); + ptResolution = mSmearer.getAbsPtRes(pdgInfoThis->PdgCode(), dNdEta, pseudorapidity, momentumHypotheses[ii] / std::cosh(pseudorapidity)); + etaResolution = mSmearer.getAbsEtaRes(pdgInfoThis->PdgCode(), dNdEta, pseudorapidity, momentumHypotheses[ii] / std::cosh(pseudorapidity)); } - float innerTrackTimeReso = calculateTrackTimeResolutionAdvanced(momentum / std::cosh(pseudorapidity), pseudorapidity, ptResolution, etaResolution, masses[ii], simConfig.innerTOFRadius, simConfig.magneticField); - float outerTrackTimeReso = calculateTrackTimeResolutionAdvanced(momentum / std::cosh(pseudorapidity), pseudorapidity, ptResolution, etaResolution, masses[ii], simConfig.outerTOFRadius, simConfig.magneticField); + float innerTrackTimeReso = calculateTrackTimeResolutionAdvanced(momentumHypotheses[ii] / std::cosh(pseudorapidity), pseudorapidity, ptResolution, etaResolution, masses[ii], simConfig.innerTOFRadius, simConfig.magneticField); + float outerTrackTimeReso = calculateTrackTimeResolutionAdvanced(momentumHypotheses[ii] / std::cosh(pseudorapidity), pseudorapidity, ptResolution, etaResolution, masses[ii], simConfig.outerTOFRadius, simConfig.magneticField); innerTotalTimeReso = std::hypot(simConfig.innerTOFTimeReso, innerTrackTimeReso); outerTotalTimeReso = std::hypot(simConfig.outerTOFTimeReso, outerTrackTimeReso); if (plotsConfig.doQAplots) { if (std::fabs(mcParticle.pdgCode()) == pdg->GetParticle(kParticlePdgs[ii])->PdgCode()) { if (trackLengthRecoInnerTOF > 0) { - h2dInnerTimeResTrack[ii]->Fill(momentum, innerTrackTimeReso); - h2dInnerTimeResTotal[ii]->Fill(momentum, innerTotalTimeReso); + h2dInnerTimeResTrack[ii]->Fill(momentumHypotheses[ii], innerTrackTimeReso); + h2dInnerTimeResTotal[ii]->Fill(momentumHypotheses[ii], innerTotalTimeReso); } if (trackLengthRecoOuterTOF > 0) { - const float transverseMomentum = momentum / std::cosh(pseudorapidity); - h2dOuterTimeResTrack[ii]->Fill(momentum, outerTrackTimeReso); - h2dOuterTimeResTotal[ii]->Fill(momentum, outerTotalTimeReso); + const float transverseMomentum = momentumHypotheses[ii] / std::cosh(pseudorapidity); + h2dOuterTimeResTrack[ii]->Fill(momentumHypotheses[ii], outerTrackTimeReso); + h2dOuterTimeResTotal[ii]->Fill(momentumHypotheses[ii], outerTotalTimeReso); static constexpr int kIdPion = 2; if (ii == kIdPion) { histos.fill(HIST("h2dRelativePtResolution"), transverseMomentum, 100.0 * ptResolution / transverseMomentum); @@ -667,14 +686,14 @@ struct OnTheFlyTofPid { } if (trackLengthRecoInnerTOF > 0) { for (int iii = 0; iii < kParticles; iii++) { - h2dInnerNsigmaTrue[ii][iii]->Fill(momentum, nSigmaInnerTOF[iii]); - h2dInnerDeltaTrue[ii][iii]->Fill(momentum, deltaTimeInnerTOF[iii]); + h2dInnerNsigmaTrue[ii][iii]->Fill(momentumHypotheses[ii], nSigmaInnerTOF[iii]); + h2dInnerDeltaTrue[ii][iii]->Fill(momentumHypotheses[ii], deltaTimeInnerTOF[iii]); } } if (trackLengthRecoOuterTOF > 0) { for (int iii = 0; iii < kParticles; iii++) { - h2dOuterNsigmaTrue[ii][iii]->Fill(momentum, nSigmaOuterTOF[iii]); - h2dOuterDeltaTrue[ii][iii]->Fill(momentum, deltaTimeOuterTOF[iii]); + h2dOuterNsigmaTrue[ii][iii]->Fill(momentumHypotheses[ii], nSigmaOuterTOF[iii]); + h2dOuterDeltaTrue[ii][iii]->Fill(momentumHypotheses[ii], deltaTimeOuterTOF[iii]); } } } @@ -691,12 +710,12 @@ struct OnTheFlyTofPid { // Sigmas have been fully calculated. Please populate the NSigma helper table (once per track) upgradeTof(tzero[0], tzero[1], - nSigmaInnerTOF[0], nSigmaInnerTOF[1], nSigmaInnerTOF[2], nSigmaInnerTOF[3], nSigmaInnerTOF[4], + nSigmaInnerTOF[0], nSigmaInnerTOF[1], nSigmaInnerTOF[2], nSigmaInnerTOF[3], nSigmaInnerTOF[4], nSigmaInnerTOF[5], nSigmaInnerTOF[6], nSigmaInnerTOF[7], nSigmaInnerTOF[8], measuredTimeInnerTOF, trackLengthRecoInnerTOF, - nSigmaOuterTOF[0], nSigmaOuterTOF[1], nSigmaOuterTOF[2], nSigmaOuterTOF[3], nSigmaOuterTOF[4], + nSigmaOuterTOF[0], nSigmaOuterTOF[1], nSigmaOuterTOF[2], nSigmaOuterTOF[3], nSigmaOuterTOF[4], nSigmaOuterTOF[5], nSigmaOuterTOF[6], nSigmaOuterTOF[7], nSigmaOuterTOF[8], measuredTimeOuterTOF, trackLengthRecoOuterTOF); - upgradeTofExpectedTime(expectedTimeInnerTOF[0], expectedTimeInnerTOF[1], expectedTimeInnerTOF[2], expectedTimeInnerTOF[3], expectedTimeInnerTOF[4], - expectedTimeOuterTOF[0], expectedTimeOuterTOF[1], expectedTimeOuterTOF[2], expectedTimeOuterTOF[3], expectedTimeOuterTOF[4]); + upgradeTofExpectedTime(expectedTimeInnerTOF[0], expectedTimeInnerTOF[1], expectedTimeInnerTOF[2], expectedTimeInnerTOF[3], expectedTimeInnerTOF[4], expectedTimeInnerTOF[5], expectedTimeInnerTOF[6], expectedTimeInnerTOF[7], expectedTimeInnerTOF[8], + expectedTimeOuterTOF[0], expectedTimeOuterTOF[1], expectedTimeOuterTOF[2], expectedTimeOuterTOF[3], expectedTimeOuterTOF[4], expectedTimeOuterTOF[5], expectedTimeOuterTOF[6], expectedTimeOuterTOF[7], expectedTimeOuterTOF[8]); } if (trackWithTimeIndex != tracks.size()) { diff --git a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx index 6ab9df2443e..b9e05d0ec75 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx @@ -91,20 +91,32 @@ static constexpr size_t kMaxValidHitsForTruncation56 = 3; static constexpr size_t kMaxValidHitsForTruncation34 = 2; static constexpr size_t kMaxValidHitsForTruncation12 = 1; +// Constants for LUT binning +// To do: Include in LUT header or similar +static constexpr int kLUTEtaBins = 50; +static constexpr float kLUTEtaMin = -2.5f; +static constexpr float kLUTEtaMax = 2.5f; +static constexpr int kLUTPtBins = 500; +static constexpr float kLUTPtMin = 0.0f; +static constexpr float kLUTPtMax = 10.0f; + class ToTLUT { public: - ToTLUT(float truncationFractionVal, int maxLayers, int etaBins, float etaMin, float etaMax, int ptBins, float ptMin, float ptMax) - : mTruncationFraction(truncationFractionVal), - mMaxLayers(maxLayers), - mEtaBins(etaBins), - mEtaMin(etaMin), - mEtaMax(etaMax), - mPtBins(ptBins), - mPtMin(ptMin), - mPtMax(ptMax), - mEtaBinWidth((etaMax - etaMin) / etaBins), - mPtBinWidth((ptMax - ptMin) / ptBins) + explicit ToTLUT(int maxLayers, float analysisEtaMinVal = 0.0f, float analysisEtaMaxVal = 1.0f, float analysisPtMinVal = 0.0f, float analysisPtMaxVal = 10.0f) + : mMaxLayers(maxLayers), + mAnalysisEtaMin(analysisEtaMinVal), + mAnalysisEtaMax(analysisEtaMaxVal), + mAnalysisPtMin(analysisPtMinVal), + mAnalysisPtMax(analysisPtMaxVal), + mEtaBins(kLUTEtaBins), + mEtaMin(kLUTEtaMin), + mEtaMax(kLUTEtaMax), + mPtBins(kLUTPtBins), + mPtMin(kLUTPtMin), + mPtMax(kLUTPtMax), + mEtaBinWidth((kLUTEtaMax - kLUTEtaMin) / kLUTEtaBins), + mPtBinWidth((kLUTPtMax - kLUTPtMin) / kLUTPtBins) { mPdgToIndexMap.reserve(10); mIndexToPdgMap.reserve(10); @@ -179,9 +191,21 @@ class ToTLUT for (int etaBin = 0; etaBin < mEtaBins; ++etaBin) { float etaMinBin = mEtaMin + etaBin * mEtaBinWidth; float etaMaxBin = etaMinBin + mEtaBinWidth; + + float etaCenter = (etaMinBin + etaMaxBin) / 2.0f; + if (std::abs(etaCenter) < mAnalysisEtaMin || std::abs(etaCenter) > mAnalysisEtaMax) { + continue; + } + for (int ptBin = 0; ptBin < mPtBins; ++ptBin) { float ptMinBin = mPtMin + ptBin * mPtBinWidth; float ptMaxBin = ptMinBin + mPtBinWidth; + float ptCenter = (ptMinBin + ptMaxBin) / 2.0f; + + if (ptCenter < mAnalysisPtMin || ptCenter >= mAnalysisPtMax) { + continue; + } + TString histName = Form("tot_%d_barrel%d_eta%.2f-%.2f_pt%.2f-%.2f", pdg, layer, etaMinBin, etaMaxBin, ptMinBin, ptMaxBin); TH1F* histFromFile = dynamic_cast(f->Get(histName)); @@ -252,8 +276,11 @@ class ToTLUT std::unordered_map mPdgToIndexMap; std::vector mIndexToPdgMap; - float mTruncationFraction; int mMaxLayers; + float mAnalysisEtaMin; + float mAnalysisEtaMax; + float mAnalysisPtMin; + float mAnalysisPtMax; int mEtaBins; float mEtaMin; float mEtaMax; @@ -380,16 +407,13 @@ struct OnTheFlyTrackerPid { Configurable lutTotHe{"lutTotHe", "ccdb:Users/h/hfribert/ToT_LUTs", "ToT LUT for helium-3"}; Configurable lutTotAl{"lutTotAl", "ccdb:Users/h/hfribert/ToT_LUTs", "ToT LUT for alphas"}; - Configurable truncationFraction{"truncationFraction", 0.80f, "Fraction of lower entries to consider for truncated standard deviation"}; Configurable dBz{"dBz", 20, "magnetic field (kilogauss) for track propagation"}; Configurable maxBarrelLayers{"maxBarrelLayers", 11, "Maximum number of barrel layers"}; - Configurable etaBins{"etaBins", 50, "Number of eta bins for LUTs and histograms"}; - Configurable etaMin{"etaMin", -2.5f, "Minimum eta value"}; - Configurable etaMax{"etaMax", 2.5f, "Maximum eta value"}; - Configurable ptBins{"ptBins", 500, "Number of pT bins for LUTs (LUTs are pT-based)"}; - Configurable ptMin{"ptMin", 0.0f, "Minimum pT value for LUT binning"}; - Configurable ptMax{"ptMax", 10.0f, "Maximum pT value for LUT binning"}; Configurable numLogBins{"numLogBins", 200, "Number of logarithmic momentum bins"}; + Configurable analysisEtaMin{"analysisEtaMin", 0.0f, "Minimum |eta| for LUT loading optimization"}; + Configurable analysisEtaMax{"analysisEtaMax", 1.0f, "Maximum |eta| for LUT loading optimization"}; + Configurable analysisPtMin{"analysisPtMin", 0.0f, "Minimum pT (GeV/c) for LUT loading optimization"}; + Configurable analysisPtMax{"analysisPtMax", 10.0f, "Maximum pT (GeV/c) for LUT loading optimization"}; std::vector mLogBins; @@ -416,9 +440,7 @@ struct OnTheFlyTrackerPid { << "). Please adjust maxBarrelLayers."; } - mToTLUT = std::make_unique(truncationFraction.value, maxBarrelLayers.value, - etaBins.value, etaMin.value, etaMax.value, - ptBins.value, ptMin.value, ptMax.value); + mToTLUT = std::make_unique(maxBarrelLayers.value, analysisEtaMin.value, analysisEtaMax.value, analysisPtMin.value, analysisPtMax.value); mToTLUT->setCcdbManager(ccdb.operator->()); @@ -526,7 +548,6 @@ struct OnTheFlyTrackerPid { aod::McParticles const& /*mcParticles*/, aod::McCollisions const& /*mcCollisions*/) { - o2::dataformats::VertexBase mcPvVtx({0.0f, 0.0f, 0.0f}, {0.}); if (collision.has_mcCollision()) { diff --git a/ALICE3/Tasks/CMakeLists.txt b/ALICE3/Tasks/CMakeLists.txt index 06864096cf5..42fb53a0a25 100644 --- a/ALICE3/Tasks/CMakeLists.txt +++ b/ALICE3/Tasks/CMakeLists.txt @@ -74,4 +74,7 @@ o2physics_add_dpl_workflow(alice3-tracking-performance PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) - +o2physics_add_dpl_workflow(alice3-pid-evaluation + SOURCES alice3PidEvaluation.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) diff --git a/ALICE3/Tasks/alice3PidEvaluation.cxx b/ALICE3/Tasks/alice3PidEvaluation.cxx new file mode 100644 index 00000000000..2ee64ea7623 --- /dev/null +++ b/ALICE3/Tasks/alice3PidEvaluation.cxx @@ -0,0 +1,371 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +/// +/// \file alice3PidEvaluation.cxx +/// +/// \brief This task computes purity and efficiency from the OTF PID tables for multiple detectors. +/// Analyzes individual detectors (Tracker, TOF Inner, TOF Outer, RICH), +/// as well as combined detector performance using quadrature combination +/// of nSigma values. +/// +/// \author Henrik Fribert TUM +/// \since August 14, 2025 +/// + +#include "ALICE3/DataModel/OTFPIDTrk.h" +#include "ALICE3/DataModel/OTFRICH.h" +#include "ALICE3/DataModel/OTFTOF.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include "CommonUtils/NameConf.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/RunningWorkflowInfo.h" +#include "Framework/runDataProcessing.h" +#include "ReconstructionDataFormats/DCA.h" + +#include "TH1F.h" +#include "TH2F.h" +#include "TProfile.h" +#include "TVector3.h" + +#include +#include +#include +#include + +using namespace o2; +using namespace o2::framework; + +struct Alice3PidEvaluation { + + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + + static constexpr float kInvalidNSigmaValue = 999.0f; + + Configurable maxNSigmaForIdentification{"maxNSigmaForIdentification", 3.0f, "Maximum |nSigma| allowed for particle identification (closest-hypothesis rule)"}; + Configurable numLogBins{"numLogBins", 200, "Number of logarithmic momentum bins"}; + Configurable useClosestHypothesisRule{"useClosestHypothesisRule", true, "Use closest-hypothesis rule: assign track to hypothesis with smallest |nSigma|"}; + Configurable useMinimalIdentification{"useMinimalIdentification", false, "Require that only one hypothesis is within the cutoff"}; + Configurable includeTrackerInCombined{"includeTrackerInCombined", true, "Include Tracker in combined analysis"}; + Configurable includeTofInnerInCombined{"includeTofInnerInCombined", true, "Include TOF Inner in combined analysis"}; + Configurable includeTofOuterInCombined{"includeTofOuterInCombined", true, "Include TOF Outer in combined analysis"}; + Configurable includeRichInCombined{"includeRichInCombined", false, "Include RICH in combined analysis"}; + + std::vector mLogBins; + + enum PidHypothesis { kElectron, + kMuon, + kPion, + kKaon, + kProton, + kDeuteron, + kTriton, + kHelium3, + kAlpha, + kCount }; + static constexpr std::array kHypothesisPdg = {11, 13, 211, 321, 2212, 1000010020, 1000010030, 1000020030, 1000020040}; + static constexpr std::array kHypothesisNames = {"Electron", "Muon", "Pion", "Kaon", "Proton", "Deuteron", "Triton", "Helium3", "Alpha"}; + + struct DetectorHistograms { + std::array, PidHypothesis::kCount> hTotalTrue; + std::array, PidHypothesis::kCount> hEfficiency; + std::array, PidHypothesis::kCount> hPurityAsHypothesis; + }; + + std::shared_ptr hDetectorParticipation2D; + + DetectorHistograms trackerHists; + DetectorHistograms tofInnerHists; + DetectorHistograms tofOuterHists; + DetectorHistograms richHists; + DetectorHistograms combinedHists; + DetectorHistograms combinedNoTrackerHists; + void init(o2::framework::InitContext&) + { + LOG(info) << "Initializing multi-detector PID evaluation using closest-hypothesis rule"; + LOG(info) << "Maximum |nSigma| for identification: " << maxNSigmaForIdentification.value; + LOG(info) << "Closest-hypothesis rule: " << (useClosestHypothesisRule.value ? "ENABLED" : "DISABLED"); + LOG(info) << "Require unique identification: " << (useMinimalIdentification.value ? "YES" : "NO"); + LOG(info) << "Combined analysis includes: " + << (includeTrackerInCombined.value ? "Tracker " : "") + << (includeTofInnerInCombined.value ? "TOF_Inner " : "") + << (includeTofOuterInCombined.value ? "TOF_Outer " : "") + << (includeRichInCombined.value ? "RICH " : ""); + + mLogBins.clear(); + double pMin = 0.05; + double pMax = 10; + double logMin = std::log10(pMin); + double logMax = std::log10(pMax); + double dLog = (logMax - logMin) / numLogBins.value; + for (int i = 0; i <= numLogBins.value; ++i) { + mLogBins.push_back(std::pow(10, logMin + i * dLog)); + } + const AxisSpec axisMomentum{mLogBins, "#it{p} (GeV/#it{c})"}; + + auto createDetectorHistograms = [&](DetectorHistograms& detHists, const std::string& detectorName) { + for (int trueIdx = 0; trueIdx < PidHypothesis::kCount; ++trueIdx) { + const auto& trueName = kHypothesisNames[trueIdx]; + detHists.hTotalTrue[trueIdx] = histos.add(Form("%s/hTotalTrue%s", detectorName.c_str(), trueName), + Form("%s: Total True %s; #it{p} (GeV/#it{c})", detectorName.c_str(), trueName), + kTH1F, {axisMomentum}); + detHists.hEfficiency[trueIdx] = histos.add(Form("%s/hEfficiency%s", detectorName.c_str(), trueName), + Form("%s: PID Efficiency for %s; #it{p} (GeV/#it{c}); Efficiency", detectorName.c_str(), trueName), + kTProfile, {axisMomentum}); + } + + for (int hypIdx = 0; hypIdx < PidHypothesis::kCount; ++hypIdx) { + const auto& hypName = kHypothesisNames[hypIdx]; + detHists.hPurityAsHypothesis[hypIdx] = histos.add(Form("%s/hPurityAs%s", detectorName.c_str(), hypName), + Form("%s: Purity when selecting as %s; #it{p} (GeV/#it{c}); Purity", detectorName.c_str(), hypName), + kTProfile, {axisMomentum}); + } + }; + + createDetectorHistograms(trackerHists, "Tracker"); + createDetectorHistograms(tofInnerHists, "TOF_Inner"); + createDetectorHistograms(tofOuterHists, "TOF_Outer"); + createDetectorHistograms(richHists, "RICH"); + createDetectorHistograms(combinedHists, "Combined"); + createDetectorHistograms(combinedNoTrackerHists, "Combined_NoTracker"); + + const AxisSpec axisDetectorCount{5, -0.5, 4.5, "Number of detectors"}; + hDetectorParticipation2D = histos.add("Combined/hDetectorParticipation2D", + "Detector participation vs momentum; #it{p} (GeV/#it{c}); Number of detectors", + kTH2F, {axisMomentum, axisDetectorCount}); + } + + void process(soa::Join const& tracks, + aod::McParticles const& /*mcParticles*/) + { + + auto isValidNSigma = [](float nSigma) -> bool { + return (nSigma < kInvalidNSigmaValue && nSigma > -kInvalidNSigmaValue); + }; + + auto computeCombinedNSigma = [&](const std::vector>& detectorNSigmas, float p) -> std::array { + std::array combinedNSigma; + int totalValidDetectors = 0; + for (const auto& detNSigma : detectorNSigmas) { + bool detectorHasValidMeasurement = false; + for (int hypIdx = 0; hypIdx < PidHypothesis::kCount; ++hypIdx) { + if (isValidNSigma(detNSigma[hypIdx])) { + detectorHasValidMeasurement = true; + break; + } + } + if (detectorHasValidMeasurement) { + totalValidDetectors++; + } + } + + for (int hypIdx = 0; hypIdx < PidHypothesis::kCount; ++hypIdx) { + float sumSquares = 0.0f; + int validDetectors = 0; + + for (const auto& detNSigma : detectorNSigmas) { + if (isValidNSigma(detNSigma[hypIdx])) { + sumSquares += detNSigma[hypIdx] * detNSigma[hypIdx]; + validDetectors++; + } + } + + if (validDetectors > 0) { + combinedNSigma[hypIdx] = std::sqrt(sumSquares); + } else { + combinedNSigma[hypIdx] = kInvalidNSigmaValue; + } + } + if (totalValidDetectors > 0) { + hDetectorParticipation2D->Fill(p, totalValidDetectors); + } + + return combinedNSigma; + }; + + auto analyzeDetector = [&](DetectorHistograms& detHists, const std::array& nSigmaValues, + int trueParticleIndex, float p) { + detHists.hTotalTrue[trueParticleIndex]->Fill(p); + bool hasValidNSigma = false; + for (int i = 0; i < PidHypothesis::kCount; ++i) { + if (isValidNSigma(nSigmaValues[i])) { + hasValidNSigma = true; + break; + } + } + if (!hasValidNSigma) { + return; + } + + bool correctlyIdentified = false; + int selectedHypothesis = -1; + + if (useClosestHypothesisRule.value) { + float minAbsNSigma = kInvalidNSigmaValue; + int bestHypothesis = -1; + int validHypothesesCount = 0; + + for (int hypIdx = 0; hypIdx < PidHypothesis::kCount; ++hypIdx) { + if (isValidNSigma(nSigmaValues[hypIdx])) { + float absNSigma = std::fabs(nSigmaValues[hypIdx]); + if (absNSigma < minAbsNSigma) { + minAbsNSigma = absNSigma; + bestHypothesis = hypIdx; + } + if (absNSigma < maxNSigmaForIdentification.value) { + validHypothesesCount++; + } + } + } + + if (bestHypothesis >= 0 && minAbsNSigma < maxNSigmaForIdentification.value) { + if (useMinimalIdentification.value && validHypothesesCount > 1) { + selectedHypothesis = -1; + } else { + selectedHypothesis = bestHypothesis; + } + } + correctlyIdentified = (selectedHypothesis == trueParticleIndex); + } else { + correctlyIdentified = (std::fabs(nSigmaValues[trueParticleIndex]) < maxNSigmaForIdentification.value); + for (int hypIdx = 0; hypIdx < PidHypothesis::kCount; ++hypIdx) { + if (std::fabs(nSigmaValues[hypIdx]) < maxNSigmaForIdentification.value) { + bool isCorrect = (hypIdx == trueParticleIndex); + detHists.hPurityAsHypothesis[hypIdx]->Fill(p, isCorrect ? 1.0 : 0.0); + } + } + } + + detHists.hEfficiency[trueParticleIndex]->Fill(p, correctlyIdentified ? 1.0 : 0.0); + + if (useClosestHypothesisRule.value && selectedHypothesis >= 0) { + bool isCorrect = (selectedHypothesis == trueParticleIndex); + detHists.hPurityAsHypothesis[selectedHypothesis]->Fill(p, isCorrect ? 1.0 : 0.0); + } + }; + + for (const auto& track : tracks) { + if (!track.has_mcParticle()) { + continue; + } + + const auto& mcParticle = track.mcParticle(); + const float p = mcParticle.p(); + const int truePdg = std::abs(mcParticle.pdgCode()); + + int trueParticleIndex = -1; + for (int i = 0; i < PidHypothesis::kCount; ++i) { + if (kHypothesisPdg[i] == truePdg) { + trueParticleIndex = i; + break; + } + } + if (trueParticleIndex == -1) { + continue; + } + + std::array trackerNSigma; + trackerNSigma[kElectron] = track.nSigmaTrkEl(); + trackerNSigma[kMuon] = track.nSigmaTrkMu(); + trackerNSigma[kPion] = track.nSigmaTrkPi(); + trackerNSigma[kKaon] = track.nSigmaTrkKa(); + trackerNSigma[kProton] = track.nSigmaTrkPr(); + trackerNSigma[kDeuteron] = track.nSigmaTrkDe(); + trackerNSigma[kTriton] = track.nSigmaTrkTr(); + trackerNSigma[kHelium3] = track.nSigmaTrkHe(); + trackerNSigma[kAlpha] = track.nSigmaTrkAl(); + + analyzeDetector(trackerHists, trackerNSigma, trueParticleIndex, p); + + std::array tofInnerNSigma; + tofInnerNSigma[kElectron] = track.nSigmaElectronInnerTOF(); + tofInnerNSigma[kMuon] = track.nSigmaMuonInnerTOF(); + tofInnerNSigma[kPion] = track.nSigmaPionInnerTOF(); + tofInnerNSigma[kKaon] = track.nSigmaKaonInnerTOF(); + tofInnerNSigma[kProton] = track.nSigmaProtonInnerTOF(); + tofInnerNSigma[kDeuteron] = track.nSigmaDeuteronInnerTOF(); + tofInnerNSigma[kTriton] = track.nSigmaTritonInnerTOF(); + tofInnerNSigma[kHelium3] = track.nSigmaHelium3InnerTOF(); + tofInnerNSigma[kAlpha] = track.nSigmaAlphaInnerTOF(); + + analyzeDetector(tofInnerHists, tofInnerNSigma, trueParticleIndex, p); + + std::array tofOuterNSigma; + tofOuterNSigma[kElectron] = track.nSigmaElectronOuterTOF(); + tofOuterNSigma[kMuon] = track.nSigmaMuonOuterTOF(); + tofOuterNSigma[kPion] = track.nSigmaPionOuterTOF(); + tofOuterNSigma[kKaon] = track.nSigmaKaonOuterTOF(); + tofOuterNSigma[kProton] = track.nSigmaProtonOuterTOF(); + tofOuterNSigma[kDeuteron] = track.nSigmaDeuteronOuterTOF(); + tofOuterNSigma[kTriton] = track.nSigmaTritonOuterTOF(); + tofOuterNSigma[kHelium3] = track.nSigmaHelium3OuterTOF(); + tofOuterNSigma[kAlpha] = track.nSigmaAlphaOuterTOF(); + + analyzeDetector(tofOuterHists, tofOuterNSigma, trueParticleIndex, p); + + std::array richNSigma; + richNSigma[kElectron] = track.nSigmaElectronRich(); + richNSigma[kMuon] = track.nSigmaMuonRich(); + richNSigma[kPion] = track.nSigmaPionRich(); + richNSigma[kKaon] = track.nSigmaKaonRich(); + richNSigma[kProton] = track.nSigmaProtonRich(); + richNSigma[kDeuteron] = track.nSigmaDeuteronRich(); + richNSigma[kTriton] = track.nSigmaTritonRich(); + richNSigma[kHelium3] = track.nSigmaHelium3Rich(); + richNSigma[kAlpha] = track.nSigmaAlphaRich(); + + analyzeDetector(richHists, richNSigma, trueParticleIndex, p); + + std::vector> allDetectorNSigmas; + + if (includeTrackerInCombined.value) { + allDetectorNSigmas.push_back(trackerNSigma); + } + if (includeTofInnerInCombined.value) { + allDetectorNSigmas.push_back(tofInnerNSigma); + } + if (includeTofOuterInCombined.value) { + allDetectorNSigmas.push_back(tofOuterNSigma); + } + if (includeRichInCombined.value) { + allDetectorNSigmas.push_back(richNSigma); + } + if (!allDetectorNSigmas.empty()) { + std::array combinedNSigma = computeCombinedNSigma(allDetectorNSigmas, p); + analyzeDetector(combinedHists, combinedNSigma, trueParticleIndex, p); + } + + std::vector> noTrackerDetectorNSigmas; + + if (includeTofInnerInCombined.value) { + noTrackerDetectorNSigmas.push_back(tofInnerNSigma); + } + if (includeTofOuterInCombined.value) { + noTrackerDetectorNSigmas.push_back(tofOuterNSigma); + } + if (includeRichInCombined.value) { + noTrackerDetectorNSigmas.push_back(richNSigma); + } + if (!noTrackerDetectorNSigmas.empty()) { + std::array combinedNoTrackerNSigma = computeCombinedNSigma(noTrackerDetectorNSigmas, p); + analyzeDetector(combinedNoTrackerHists, combinedNoTrackerNSigma, trueParticleIndex, p); + } + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +}