From e8ef4854d19bb151618e883529f4f2562cc99b68 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Thu, 16 Oct 2025 15:07:44 +0200 Subject: [PATCH 1/2] A3 OTF: Check available LUT before smearing --- ALICE3/Core/DelphesO2TrackSmearer.cxx | 31 ++--- ALICE3/Core/DelphesO2TrackSmearer.h | 25 ++-- ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx | 122 ++++++++++--------- ALICE3/TableProducer/OTF/onTheFlyTofPid.cxx | 46 ++++--- 4 files changed, 125 insertions(+), 99 deletions(-) diff --git a/ALICE3/Core/DelphesO2TrackSmearer.cxx b/ALICE3/Core/DelphesO2TrackSmearer.cxx index 0002160a49b..599d9cd9413 100644 --- a/ALICE3/Core/DelphesO2TrackSmearer.cxx +++ b/ALICE3/Core/DelphesO2TrackSmearer.cxx @@ -153,12 +153,13 @@ bool TrackSmearer::loadTable(int pdg, const char* filename, bool forceReload) /*****************************************************************/ -lutEntry_t* - TrackSmearer::getLUTEntry(int pdg, float nch, float radius, float eta, float pt, float& interpolatedEff) +lutEntry_t* TrackSmearer::getLUTEntry(const int pdg, const float nch, const float radius, const float eta, const float pt, float& interpolatedEff) { - auto ipdg = getIndexPDG(pdg); - if (!mLUTHeader[ipdg]) + const int ipdg = getIndexPDG(pdg); + if (!mLUTHeader[ipdg]) { + LOG(error) << " --- getLUTEntry: LUT header not loaded for pdg=" << pdg << ". Returning nullptr."; return nullptr; + } auto inch = mLUTHeader[ipdg]->nchmap.find(nch); auto irad = mLUTHeader[ipdg]->radmap.find(radius); auto ieta = mLUTHeader[ipdg]->etamap.find(eta); @@ -290,7 +291,7 @@ bool TrackSmearer::smearTrack(O2Track& o2track, int pdg, float nch) } auto eta = o2track.getEta(); float interpolatedEff = 0.0f; - auto lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, interpolatedEff); + lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, interpolatedEff); if (!lutEntry || !lutEntry->valid) return false; return smearTrack(o2track, lutEntry, interpolatedEff); @@ -298,20 +299,20 @@ bool TrackSmearer::smearTrack(O2Track& o2track, int pdg, float nch) /*****************************************************************/ // relative uncertainty on pt -double TrackSmearer::getPtRes(int pdg, float nch, float eta, float pt) +double TrackSmearer::getPtRes(const int pdg, const float nch, const float eta, const float pt) { float dummy = 0.0f; - auto lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy); + lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy); auto val = std::sqrt(lutEntry->covm[14]) * lutEntry->pt; return val; } /*****************************************************************/ // relative uncertainty on eta -double TrackSmearer::getEtaRes(int pdg, float nch, float eta, float pt) +double TrackSmearer::getEtaRes(const int pdg, const float nch, const float eta, const float pt) { float dummy = 0.0f; - auto lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy); + lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy); auto sigmatgl = std::sqrt(lutEntry->covm[9]); // sigmatgl2 auto etaRes = std::fabs(std::sin(2.0 * std::atan(std::exp(-eta)))) * sigmatgl; // propagate tgl to eta uncertainty etaRes /= lutEntry->eta; // relative uncertainty @@ -319,27 +320,27 @@ double TrackSmearer::getEtaRes(int pdg, float nch, float eta, float pt) } /*****************************************************************/ // absolute uncertainty on pt -double TrackSmearer::getAbsPtRes(int pdg, float nch, float eta, float pt) +double TrackSmearer::getAbsPtRes(const int pdg, const float nch, const float eta, const float pt) { float dummy = 0.0f; - auto lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy); + lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy); auto val = std::sqrt(lutEntry->covm[14]) * lutEntry->pt * lutEntry->pt; return val; } /*****************************************************************/ // absolute uncertainty on eta -double TrackSmearer::getAbsEtaRes(int pdg, float nch, float eta, float pt) +double TrackSmearer::getAbsEtaRes(const int pdg, const float nch, const float eta, const float pt) { float dummy = 0.0f; - auto lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy); + lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0., eta, pt, dummy); auto sigmatgl = std::sqrt(lutEntry->covm[9]); // sigmatgl2 auto etaRes = std::fabs(std::sin(2.0 * std::atan(std::exp(-eta)))) * sigmatgl; // propagate tgl to eta uncertainty return etaRes; } /*****************************************************************/ // efficiency -double TrackSmearer::getEfficiency(int pdg, float nch, float eta, float pt) +double TrackSmearer::getEfficiency(const int pdg, const float nch, const float eta, const float pt) { float efficiency = 0.0f; getLUTEntry(pdg, nch, 0., eta, pt, efficiency); @@ -360,7 +361,7 @@ double TrackSmearer::getEfficiency(int pdg, float nch, float eta, float pt) // return true; // #if 0 -// auto lutEntry = getLUTEntry(track.PID, 0., 0., track.Eta, track.PT); +// lutEntry_t* lutEntry = getLUTEntry(track.PID, 0., 0., track.Eta, track.PT); // if (!lutEntry) // return; diff --git a/ALICE3/Core/DelphesO2TrackSmearer.h b/ALICE3/Core/DelphesO2TrackSmearer.h index 770b51e75a6..027809fa004 100644 --- a/ALICE3/Core/DelphesO2TrackSmearer.h +++ b/ALICE3/Core/DelphesO2TrackSmearer.h @@ -180,23 +180,24 @@ class TrackSmearer /** LUT methods **/ bool loadTable(int pdg, const char* filename, bool forceReload = false); - void useEfficiency(bool val) { mUseEfficiency = val; } //; - void interpolateEfficiency(bool val) { mInterpolateEfficiency = val; } //; - void skipUnreconstructed(bool val) { mSkipUnreconstructed = val; } //; - void setWhatEfficiency(int val) { mWhatEfficiency = val; } //; - lutHeader_t* getLUTHeader(int pdg) { return mLUTHeader[getIndexPDG(pdg)]; } //; - lutEntry_t* getLUTEntry(int pdg, float nch, float radius, float eta, float pt, float& interpolatedEff); + bool hasTable(int pdg) { return (mLUTHeader[getIndexPDG(pdg)] != nullptr); } //; + void useEfficiency(bool val) { mUseEfficiency = val; } //; + void interpolateEfficiency(bool val) { mInterpolateEfficiency = val; } //; + void skipUnreconstructed(bool val) { mSkipUnreconstructed = val; } //; + void setWhatEfficiency(int val) { mWhatEfficiency = val; } //; + lutHeader_t* getLUTHeader(int pdg) { return mLUTHeader[getIndexPDG(pdg)]; } //; + lutEntry_t* getLUTEntry(const int pdg, const float nch, const float radius, const float eta, const float pt, float& interpolatedEff); bool smearTrack(O2Track& o2track, lutEntry_t* lutEntry, float interpolatedEff); bool smearTrack(O2Track& o2track, int pdg, float nch); // bool smearTrack(Track& track, bool atDCA = true); // Only in DelphesO2 - double getPtRes(int pdg, float nch, float eta, float pt); - double getEtaRes(int pdg, float nch, float eta, float pt); - double getAbsPtRes(int pdg, float nch, float eta, float pt); - double getAbsEtaRes(int pdg, float nch, float eta, float pt); - double getEfficiency(int pdg, float nch, float eta, float pt); + double getPtRes(const int pdg, const float nch, const float eta, const float pt); + double getEtaRes(const int pdg, const float nch, const float eta, const float pt); + double getAbsPtRes(const int pdg, const float nch, const float eta, const float pt); + double getAbsEtaRes(const int pdg, const float nch, const float eta, const float pt); + double getEfficiency(const int pdg, const float nch, const float eta, const float pt); - int getIndexPDG(int pdg) + int getIndexPDG(const int pdg) { switch (abs(pdg)) { case 11: diff --git a/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx index ad7dccc6748..961d293144e 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx @@ -48,7 +48,6 @@ #include #include #include -#include #include #include #include @@ -57,6 +56,7 @@ #include #include #include +#include #include #include @@ -84,7 +84,7 @@ struct OnTheFlyRichPid { Service ccdb; // master setting: magnetic field - Configurable magneticField{"magneticField", 0, "magnetic field (kilogauss) if 0, taken from the tracker task"}; + float magneticField = 0; // Always taken from the tracker task // add rich-specific configurables here Configurable bRichNumberOfSectors{"bRichNumberOfSectors", 21, "barrel RICH number of sectors"}; @@ -290,12 +290,8 @@ struct OnTheFlyRichPid { { pRandomNumberGenerator.SetSeed(0); // fully randomize - if (magneticField.value < o2::constants::math::Epsilon) { - LOG(info) << "Getting the magnetic field from the on-the-fly tracker task"; - if (!getTaskOptionValue(initContext, "on-the-fly-tracker", magneticField, false)) { - LOG(fatal) << "Could not get Bz from on-the-fly-tracker task"; - } - LOG(info) << "Bz = " << magneticField.value << " T"; + if (!getTaskOptionValue(initContext, "on-the-fly-tracker", "magneticField", magneticField, false)) { + LOG(fatal) << "Could not get Bz from on-the-fly-tracker task"; } // Load LUT for pt and eta smearing @@ -812,21 +808,29 @@ struct OnTheFlyRichPid { 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, 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::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]}; + static constexpr int kParticlePdgs[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 kParticleMasses[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::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 float hypothesisAngleBarrelRich = kErrorValue; - const bool hypothesisAngleBarrelRichOk = cherenkovAngle(recoTrack.getP(), kMasses[ii], aerogelRindex[iSecor], hypothesisAngleBarrelRich); + const bool hypothesisAngleBarrelRichOk = cherenkovAngle(recoTrack.getP(), kParticleMasses[ii], aerogelRindex[iSecor], hypothesisAngleBarrelRich); signalBarrelRich[ii] = hypothesisAngleBarrelRichOk; // Particle is above the threshold and enough photons // Evaluate total sigma (layer + tracking resolution) @@ -836,11 +840,13 @@ struct OnTheFlyRichPid { double ptResolution = transverseMomentum * transverseMomentum * std::sqrt(recoTrack.getSigma1Pt2()); double etaResolution = std::fabs(std::sin(2.0 * std::atan(std::exp(-recoTrack.getEta())))) * std::sqrt(recoTrack.getSigmaTgl2()); if (flagRICHLoadDelphesLUTs) { - ptResolution = mSmearer.getAbsPtRes(kPdgArray[ii], dNdEta, recoTrack.getEta(), transverseMomentum); - etaResolution = mSmearer.getAbsEtaRes(kPdgArray[ii], dNdEta, recoTrack.getEta(), transverseMomentum); + if (mSmearer.hasTable(kParticlePdgs[ii])) { + ptResolution = mSmearer.getAbsPtRes(kParticlePdgs[ii], dNdEta, recoTrack.getEta(), transverseMomentum); + etaResolution = mSmearer.getAbsEtaRes(kParticlePdgs[ii], dNdEta, recoTrack.getEta(), transverseMomentum); + } } // cout << endl << "Pt resolution: " << ptResolution << ", Eta resolution: " << etaResolution << endl << endl; - const float barrelTrackAngularReso = calculateTrackAngularResolutionAdvanced(recoTrack.getP() / std::cosh(recoTrack.getEta()), recoTrack.getEta(), ptResolution, etaResolution, kMasses[ii], aerogelRindex[iSecor]); + const float barrelTrackAngularReso = calculateTrackAngularResolutionAdvanced(recoTrack.getP() / std::cosh(recoTrack.getEta()), recoTrack.getEta(), ptResolution, etaResolution, kParticleMasses[ii], aerogelRindex[iSecor]); barrelTotalAngularReso = std::hypot(barrelRICHAngularResolution, barrelTrackAngularReso); if (doQAplots && hypothesisAngleBarrelRich > kErrorValue + 1. && @@ -848,64 +854,64 @@ struct OnTheFlyRichPid { barrelRICHAngularResolution > kErrorValue + 1. && flagReachesRadiator) { switch (mcParticle.pdgCode()) { - case kPdgArray[kEl]: // Electron - case -kPdgArray[kEl]: // Positron + case kParticlePdgs[kEl]: // Electron + case -kParticlePdgs[kEl]: // Positron if (ii == kEl) { histos.fill(HIST("h2dBarrelAngularResTrackElecVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); histos.fill(HIST("h2dBarrelAngularResTotalElecVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); } break; - case kPdgArray[kMu]: // Muon - case -kPdgArray[kMu]: // AntiMuon + case kParticlePdgs[kMu]: // Muon + case -kParticlePdgs[kMu]: // AntiMuon if (ii == kMu) { histos.fill(HIST("h2dBarrelAngularResTrackMuonVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); histos.fill(HIST("h2dBarrelAngularResTotalMuonVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); } break; - case kPdgArray[kPi]: // Pion - case -kPdgArray[kPi]: // AntiPion + case kParticlePdgs[kPi]: // Pion + case -kParticlePdgs[kPi]: // AntiPion if (ii == kPi) { histos.fill(HIST("h2dBarrelAngularResTrackPionVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); histos.fill(HIST("h2dBarrelAngularResTotalPionVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); } break; - case kPdgArray[kKa]: // Kaon - case -kPdgArray[kKa]: // AntiKaon + case kParticlePdgs[kKa]: // Kaon + case -kParticlePdgs[kKa]: // AntiKaon if (ii == kKa) { histos.fill(HIST("h2dBarrelAngularResTrackKaonVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); histos.fill(HIST("h2dBarrelAngularResTotalKaonVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); } break; - case kPdgArray[kPr]: // Proton - case -kPdgArray[kPr]: // AntiProton + case kParticlePdgs[kPr]: // Proton + case -kParticlePdgs[kPr]: // AntiProton if (ii == kPr) { histos.fill(HIST("h2dBarrelAngularResTrackProtVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); histos.fill(HIST("h2dBarrelAngularResTotalProtVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); } break; - case kPdgArray[kDe]: // Deuteron - case -kPdgArray[kDe]: // AntiDeuteron + case kParticlePdgs[kDe]: // Deuteron + case -kParticlePdgs[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 + case kParticlePdgs[kTr]: // Triton + case -kParticlePdgs[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 + case kParticlePdgs[kHe3]: // Helium3 + case -kParticlePdgs[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 + case kParticlePdgs[kAl]: // Alpha + case -kParticlePdgs[kAl]: // AntiAlpha if (ii == kAl) { histos.fill(HIST("h2dBarrelAngularResTrackAlVsP"), recoTrack.getP(), 1000.0 * barrelTrackAngularReso); histos.fill(HIST("h2dBarrelAngularResTotalAlVsP"), recoTrack.getP(), 1000.0 * barrelTotalAngularReso); @@ -941,40 +947,40 @@ struct OnTheFlyRichPid { histos.fill(HIST("hSectorID"), iSecor); switch (mcParticle.pdgCode()) { - case kPdgArray[kEl]: // Electron - case -kPdgArray[kEl]: // Positron + case kParticlePdgs[kEl]: // Electron + case -kParticlePdgs[kEl]: // Positron histos.fill(HIST("h2dBarrelNsigmaTrueElecVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]); histos.fill(HIST("h2dBarrelNsigmaTrueElecVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]); histos.fill(HIST("h2dBarrelNsigmaTrueElecVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]); histos.fill(HIST("h2dBarrelNsigmaTrueElecVsKaonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[3]); histos.fill(HIST("h2dBarrelNsigmaTrueElecVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]); break; - case kPdgArray[kMu]: // Muon - case -kPdgArray[kMu]: // AntiMuon + case kParticlePdgs[kMu]: // Muon + case -kParticlePdgs[kMu]: // AntiMuon histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]); histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]); histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]); histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsKaonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[3]); histos.fill(HIST("h2dBarrelNsigmaTrueMuonVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]); break; - case kPdgArray[kPi]: // Pion - case -kPdgArray[kPi]: // AntiPion + case kParticlePdgs[kPi]: // Pion + case -kParticlePdgs[kPi]: // AntiPion histos.fill(HIST("h2dBarrelNsigmaTruePionVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]); histos.fill(HIST("h2dBarrelNsigmaTruePionVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]); histos.fill(HIST("h2dBarrelNsigmaTruePionVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]); histos.fill(HIST("h2dBarrelNsigmaTruePionVsKaonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[3]); histos.fill(HIST("h2dBarrelNsigmaTruePionVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]); break; - case kPdgArray[kKa]: // Kaon - case -kPdgArray[kKa]: // AntiKaon + case kParticlePdgs[kKa]: // Kaon + case -kParticlePdgs[kKa]: // AntiKaon histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]); histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]); histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]); histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsKaonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[3]); histos.fill(HIST("h2dBarrelNsigmaTrueKaonVsProtHypothesis"), recoTrack.getP(), nSigmaBarrelRich[4]); break; - case kPdgArray[kPr]: // Proton - case -kPdgArray[kPr]: // AntiProton + case kParticlePdgs[kPr]: // Proton + case -kParticlePdgs[kPr]: // AntiProton histos.fill(HIST("h2dBarrelNsigmaTrueProtVsElecHypothesis"), recoTrack.getP(), nSigmaBarrelRich[0]); histos.fill(HIST("h2dBarrelNsigmaTrueProtVsMuonHypothesis"), recoTrack.getP(), nSigmaBarrelRich[1]); histos.fill(HIST("h2dBarrelNsigmaTrueProtVsPionHypothesis"), recoTrack.getP(), nSigmaBarrelRich[2]); @@ -982,31 +988,31 @@ struct OnTheFlyRichPid { 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 + case kParticlePdgs[kDe]: // Deuteron + case -kParticlePdgs[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 + case kParticlePdgs[kTr]: // Triton + case -kParticlePdgs[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 + case kParticlePdgs[kHe3]: // Helium3 + case -kParticlePdgs[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 + case kParticlePdgs[kAl]: // Alpha + case -kParticlePdgs[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]); diff --git a/ALICE3/TableProducer/OTF/onTheFlyTofPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyTofPid.cxx index f1755937566..2b091aae532 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTofPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTofPid.cxx @@ -42,7 +42,6 @@ #include #include #include -#include #include #include #include @@ -51,6 +50,7 @@ #include #include #include +#include #include #include @@ -583,8 +583,25 @@ 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, o2::constants::physics::kDeuteron, o2::constants::physics::kTriton, o2::constants::physics::kHelium3, o2::constants::physics::kAlpha}; - float masses[kParticles]; + 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}; + static constexpr float kParticleMasses[kParticles] = {o2::constants::physics::MassElectron, + o2::constants::physics::MassMuon, + o2::constants::physics::MassPionCharged, + o2::constants::physics::MassKaonCharged, + o2::constants::physics::MassProton, + o2::constants::physics::MassDeuteron, + o2::constants::physics::MassTriton, + o2::constants::physics::MassHelium3, + o2::constants::physics::MassAlpha}; + static constexpr float kParticleCharges[kParticles] = {1.f, 1.f, 1.f, 1.f, 1.f, 1.f, 1.f, 2.f, 2.f}; float momentumHypotheses[kParticles]; // Store momentum hypothesis for each particle auto truePdgInfo = pdg->GetParticle(mcParticle.pdgCode()); @@ -616,6 +633,7 @@ struct OnTheFlyTofPid { } } + // For every mass hypothesis compute the expected time, the delta with respect to it and the nsigma for (int ii = 0; ii < kParticles; ii++) { expectedTimeInnerTOF[ii] = -100; expectedTimeOuterTOF[ii] = -100; @@ -624,10 +642,8 @@ struct OnTheFlyTofPid { nSigmaInnerTOF[ii] = -100; nSigmaOuterTOF[ii] = -100; - auto pdgInfoThis = pdg->GetParticle(kParticlePdgs[ii]); - masses[ii] = pdgInfoThis->Mass(); - momentumHypotheses[ii] = rigidity * (std::abs(pdgInfoThis->Charge()) / 3.0f); // Total momentum for this hypothesis - const float v = computeParticleVelocity(momentumHypotheses[ii], masses[ii]); + momentumHypotheses[ii] = rigidity * kParticleCharges[ii]; // Total momentum for this hypothesis + const float v = computeParticleVelocity(momentumHypotheses[ii], kParticleMasses[ii]); expectedTimeInnerTOF[ii] = trackLengthInnerTOF / v; expectedTimeOuterTOF[ii] = trackLengthOuterTOF / v; @@ -639,25 +655,27 @@ struct OnTheFlyTofPid { float innerTotalTimeReso = simConfig.innerTOFTimeReso; float outerTotalTimeReso = simConfig.outerTOFTimeReso; if (simConfig.flagIncludeTrackTimeRes) { - double ptResolution = std::pow(momentumHypotheses[ii] / std::cosh(pseudorapidity), 2) * std::sqrt(trkWithTime.mMomentum.second); + const float transverseMomentum = momentumHypotheses[ii] / std::cosh(pseudorapidity); + double ptResolution = transverseMomentum * transverseMomentum * 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, momentumHypotheses[ii] / std::cosh(pseudorapidity)); - etaResolution = mSmearer.getAbsEtaRes(pdgInfoThis->PdgCode(), dNdEta, pseudorapidity, momentumHypotheses[ii] / std::cosh(pseudorapidity)); + if (mSmearer.hasTable(kParticlePdgs[ii])) { // Only if the LUT for this particle was loaded + ptResolution = mSmearer.getAbsPtRes(kParticlePdgs[ii], dNdEta, pseudorapidity, transverseMomentum); + etaResolution = mSmearer.getAbsEtaRes(kParticlePdgs[ii], dNdEta, pseudorapidity, transverseMomentum); + } } - 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); + const float innerTrackTimeReso = calculateTrackTimeResolutionAdvanced(transverseMomentum, pseudorapidity, ptResolution, etaResolution, kParticleMasses[ii], simConfig.innerTOFRadius, simConfig.magneticField); + const float outerTrackTimeReso = calculateTrackTimeResolutionAdvanced(transverseMomentum, pseudorapidity, ptResolution, etaResolution, kParticleMasses[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 (std::fabs(mcParticle.pdgCode()) == kParticlePdgs[ii]) { if (trackLengthRecoInnerTOF > 0) { h2dInnerTimeResTrack[ii]->Fill(momentumHypotheses[ii], innerTrackTimeReso); h2dInnerTimeResTotal[ii]->Fill(momentumHypotheses[ii], innerTotalTimeReso); } if (trackLengthRecoOuterTOF > 0) { - 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; From 25bc54c78087d34d3e4fb257aed0f336ec264118 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Thu, 16 Oct 2025 17:21:34 +0200 Subject: [PATCH 2/2] Refactor magneticField to use Configurable type --- ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx index 961d293144e..357e208f95f 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx @@ -84,7 +84,7 @@ struct OnTheFlyRichPid { Service ccdb; // master setting: magnetic field - float magneticField = 0; // Always taken from the tracker task + Configurable magneticField{"magneticField", 0, "magnetic field (kilogauss) if 0, taken from the tracker task"}; // add rich-specific configurables here Configurable bRichNumberOfSectors{"bRichNumberOfSectors", 21, "barrel RICH number of sectors"}; @@ -290,8 +290,12 @@ struct OnTheFlyRichPid { { pRandomNumberGenerator.SetSeed(0); // fully randomize - if (!getTaskOptionValue(initContext, "on-the-fly-tracker", "magneticField", magneticField, false)) { - LOG(fatal) << "Could not get Bz from on-the-fly-tracker task"; + if (magneticField.value < o2::constants::math::Epsilon) { + LOG(info) << "Getting the magnetic field from the on-the-fly tracker task"; + if (!getTaskOptionValue(initContext, "on-the-fly-tracker", magneticField, false)) { + LOG(fatal) << "Could not get Bz from on-the-fly-tracker task"; + } + LOG(info) << "Bz = " << magneticField.value << " T"; } // Load LUT for pt and eta smearing