From 15b0f6a3ecd333698ae3bfa7ef739725d62d69be Mon Sep 17 00:00:00 2001 From: Henrik Fribert Date: Tue, 23 Sep 2025 16:53:14 +0200 Subject: [PATCH 01/23] Tracker PID initial implementation --- ALICE3/DataModel/OTFPIDTrk.h | 39 +- ALICE3/TableProducer/OTF/onTheFlyTracker.cxx | 86 +-- .../TableProducer/OTF/onTheFlyTrackerPid.cxx | 573 +++++++++++++++--- 3 files changed, 568 insertions(+), 130 deletions(-) diff --git a/ALICE3/DataModel/OTFPIDTrk.h b/ALICE3/DataModel/OTFPIDTrk.h index b10923a892b..f0a66496ee5 100644 --- a/ALICE3/DataModel/OTFPIDTrk.h +++ b/ALICE3/DataModel/OTFPIDTrk.h @@ -11,7 +11,6 @@ /// /// \file OTFPIDTrk.h -/// \author Berkin Ulukutlu TUM /// \author Henrik Fribert TUM /// \author Nicolò Jacazio Università del Piemonte Orientale /// \since May 22, 2025 @@ -29,16 +28,18 @@ namespace o2::aod namespace upgrade::trk { -DECLARE_SOA_COLUMN(TimeOverThresholdBarrel, timeOverThresholdBarrel, float); //! Time over threshold for the barrel layers -DECLARE_SOA_COLUMN(ClusterSizeBarrel, clusterSizeBarrel, float); //! Cluster size for the barrel layers +DECLARE_SOA_COLUMN(TimeOverThresholdBarrel, timeOverThresholdBarrel, float); //! Time over threshold for the Barrel layers DECLARE_SOA_COLUMN(TimeOverThresholdForward, timeOverThresholdForward, float); //! Time over threshold for the Forward layers -DECLARE_SOA_COLUMN(ClusterSizeForward, clusterSizeForward, float); //! Cluster size for the barrel layers DECLARE_SOA_COLUMN(NSigmaTrkEl, nSigmaEl, float); //! NSigma electron from the tracker layers DECLARE_SOA_COLUMN(NSigmaTrkMu, nSigmaMu, float); //! NSigma muon from the tracker layers DECLARE_SOA_COLUMN(NSigmaTrkPi, nSigmaPi, float); //! NSigma pion from the tracker layers DECLARE_SOA_COLUMN(NSigmaTrkKa, nSigmaKa, float); //! NSigma kaon from the tracker layers DECLARE_SOA_COLUMN(NSigmaTrkPr, nSigmaPr, float); //! NSigma proton from the tracker layers +DECLARE_SOA_COLUMN(NSigmaTrkDe, nSigmaDe, float); //! NSigma deuteron from the tracker layers +DECLARE_SOA_COLUMN(NSigmaTrkTr, nSigmaTr, float); //! NSigma triton from the tracker layers +DECLARE_SOA_COLUMN(NSigmaTrkHe, nSigmaHe, float); //! NSigma helium-3 from the tracker layers +DECLARE_SOA_COLUMN(NSigmaTrkAl, nSigmaAl, float); //! NSigma alpha from the tracker layers DECLARE_SOA_DYNAMIC_COLUMN(NSigmaTrk, nSigmaTrk, //! General function to get the nSigma for the tracker layers [](const float el, @@ -46,6 +47,10 @@ DECLARE_SOA_DYNAMIC_COLUMN(NSigmaTrk, nSigmaTrk, //! General function to get the const float pi, const float ka, const float pr, + const float de, + const float tr, + const float he, + const float al, const int id) -> float { switch (std::abs(id)) { case 0: @@ -58,8 +63,16 @@ DECLARE_SOA_DYNAMIC_COLUMN(NSigmaTrk, nSigmaTrk, //! General function to get the return ka; case 4: return pr; + case 5: + return de; + case 6: + return tr; + case 7: + return he; + case 8: + return al; default: - LOG(fatal) << "Unrecognized PDG code for InnerTOF"; + LOG(fatal) << "Unrecognized PDG code"; return 999.f; } }); @@ -67,20 +80,30 @@ DECLARE_SOA_DYNAMIC_COLUMN(NSigmaTrk, nSigmaTrk, //! General function to get the } // namespace upgrade::trk DECLARE_SOA_TABLE(UpgradeTrkPidSignals, "AOD", "UPGRADETRKSIG", - upgrade::trk::TimeOverThresholdBarrel, - upgrade::trk::ClusterSizeBarrel); + o2::soa::Index<>, + upgrade::trk::TimeOverThresholdBarrel); DECLARE_SOA_TABLE(UpgradeTrkPids, "AOD", "UPGRADETRKPID", + o2::soa::Index<>, upgrade::trk::NSigmaTrkEl, upgrade::trk::NSigmaTrkMu, upgrade::trk::NSigmaTrkPi, upgrade::trk::NSigmaTrkKa, upgrade::trk::NSigmaTrkPr, + upgrade::trk::NSigmaTrkDe, + upgrade::trk::NSigmaTrkTr, + upgrade::trk::NSigmaTrkHe, + upgrade::trk::NSigmaTrkAl, upgrade::trk::NSigmaTrk); + upgrade::trk::NSigmaTrkPr, + upgrade::trk::NSigmaTrkDe, + upgrade::trk::NSigmaTrkTr, + upgrade::trk::NSigmaTrkHe, + upgrade::trk::NSigmaTrkAl + >); using UpgradeTrkPidSignal = UpgradeTrkPidSignals::iterator; using UpgradeTrkPid = UpgradeTrkPids::iterator; diff --git a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx index 0f785534159..cbb187b8359 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx @@ -113,6 +113,7 @@ struct OnTheFlyTracker { Configurable lutDe{"lutDe", "lutCovm.de.dat", "LUT for deuterons"}; Configurable lutTr{"lutTr", "lutCovm.tr.dat", "LUT for tritons"}; Configurable lutHe3{"lutHe3", "lutCovm.he3.dat", "LUT for Helium-3"}; + Configurable lutAl{"lutAl", "lutCovm.he3.dat", "LUT for Alphas"}; // To be created, for now propagating as He-3 struct : ConfigurableGroup { ConfigurableAxis axisMomentum{"axisMomentum", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f, 1.7f, 1.8f, 1.9f, 2.0f, 2.2f, 2.4f, 2.6f, 2.8f, 3.0f, 3.2f, 3.4f, 3.6f, 3.8f, 4.0f, 4.4f, 4.8f, 5.2f, 5.6f, 6.0f, 6.5f, 7.0f, 7.5f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 17.0f, 19.0f, 21.0f, 23.0f, 25.0f, 30.0f, 35.0f, 40.0f, 50.0f}, "#it{p} (GeV/#it{c})"}; @@ -138,12 +139,11 @@ struct OnTheFlyTracker { Configurable minSiliconHits{"minSiliconHits", 6, "minimum number of silicon hits to accept track"}; Configurable minSiliconHitsIfTPCUsed{"minSiliconHitsIfTPCUsed", 2, "minimum number of silicon hits to accept track in case TPC info is present"}; Configurable minTPCClusters{"minTPCClusters", 70, "minimum number of TPC hits necessary to consider minSiliconHitsIfTPCUsed"}; - Configurable alice3detector{"alice3detector", 2, "0: ALICE 3 v1, 1: ALICE 3 v4, 2: ALICE 3 Sep 2025"}; + Configurable alice3detector{"alice3detector", 0, "0: ALICE 3 v1, 1: ALICE 3 v4"}; Configurable applyZacceptance{"applyZacceptance", false, "apply z limits to detector layers or not"}; Configurable applyMSCorrection{"applyMSCorrection", true, "apply ms corrections for secondaries or not"}; Configurable applyElossCorrection{"applyElossCorrection", true, "apply eloss corrections for secondaries or not"}; Configurable applyEffCorrection{"applyEffCorrection", true, "apply efficiency correction or not"}; - Configurable scaleVD{"scaleVD", 1, "scale x0 and xrho in VD layers"}; Configurable> pixelRes{"pixelRes", {0.00025, 0.00025, 0.001, 0.001}, "RPhiIT, ZIT, RPhiOT, ZOT"}; } fastTrackerSettings; // allows for gap between peak and bg in case someone wants to @@ -240,32 +240,44 @@ struct OnTheFlyTracker { // For TGenPhaseSpace seed TRandom3 rand; - Service ccdb; void init(o2::framework::InitContext&) { - - ccdb->setURL("http://alice-ccdb.cern.ch"); - ccdb->setTimestamp(-1); - if (enableLUT) { - mSmearer.setCcdbManager(ccdb.operator->()); - - auto loadLUT = [&](int pdg, const std::string& lutFile) { - bool success = mSmearer.loadTable(pdg, lutFile.c_str()); - if (!success && !lutFile.empty()) { - LOG(fatal) << "Having issue with loading the LUT " << pdg << " " << lutFile; + std::map mapPdgLut; + const char* lutElChar = lutEl->c_str(); + const char* lutMuChar = lutMu->c_str(); + const char* lutPiChar = lutPi->c_str(); + const char* lutKaChar = lutKa->c_str(); + const char* lutPrChar = lutPr->c_str(); + + LOGF(info, "Will load electron lut file ..: %s", lutElChar); + LOGF(info, "Will load muon lut file ......: %s", lutMuChar); + LOGF(info, "Will load pion lut file ......: %s", lutPiChar); + LOGF(info, "Will load kaon lut file ......: %s", lutKaChar); + LOGF(info, "Will load proton lut file ....: %s", lutPrChar); + + mapPdgLut.insert(std::make_pair(11, lutElChar)); + mapPdgLut.insert(std::make_pair(13, lutMuChar)); + mapPdgLut.insert(std::make_pair(211, lutPiChar)); + mapPdgLut.insert(std::make_pair(321, lutKaChar)); + mapPdgLut.insert(std::make_pair(2212, lutPrChar)); + + if (enableNucleiSmearing) { + const char* lutDeChar = lutDe->c_str(); + const char* lutTrChar = lutTr->c_str(); + const char* lutHe3Char = lutHe3->c_str(); + const char* lutAlChar = lutAl->c_str(); + mapPdgLut.insert(std::make_pair(1000010020, lutDeChar)); + mapPdgLut.insert(std::make_pair(1000010030, lutTrChar)); + mapPdgLut.insert(std::make_pair(1000020030, lutHe3Char)); + mapPdgLut.insert(std::make_pair(1000020030, lutAlChar)); + } + for (const auto& e : mapPdgLut) { + if (!mSmearer.loadTable(e.first, e.second)) { + LOG(fatal) << "Having issue with loading the LUT " << e.first << " " << e.second; } - }; - loadLUT(11, lutEl.value); - loadLUT(13, lutMu.value); - loadLUT(211, lutPi.value); - loadLUT(321, lutKa.value); - loadLUT(2212, lutPr.value); - loadLUT(1000010020, lutDe.value); - loadLUT(1000010030, lutTr.value); - loadLUT(1000020030, lutHe3.value); - + } // interpolate efficiencies if requested to do so mSmearer.interpolateEfficiency(static_cast(interpolateLutEfficiencyVsNch)); @@ -413,22 +425,12 @@ struct OnTheFlyTracker { fastTracker.SetApplyMSCorrection(fastTrackerSettings.applyMSCorrection); fastTracker.SetApplyElossCorrection(fastTrackerSettings.applyElossCorrection); - switch (fastTrackerSettings.alice3detector) { - case 0: - fastTracker.AddSiliconALICE3v2(fastTrackerSettings.pixelRes); - break; - - case 1: - fastTracker.AddSiliconALICE3v4(fastTrackerSettings.pixelRes); - fastTracker.AddTPC(0.1, 0.1); - break; - - case 2: - fastTracker.AddSiliconALICE3(fastTrackerSettings.scaleVD, fastTrackerSettings.pixelRes); - break; - - default: - break; + if (fastTrackerSettings.alice3detector == 0) { + fastTracker.AddSiliconALICE3v2(fastTrackerSettings.pixelRes); + } + if (fastTrackerSettings.alice3detector == 1) { + fastTracker.AddSiliconALICE3v4(fastTrackerSettings.pixelRes); + fastTracker.AddTPC(0.1, 0.1); } // print fastTracker settings @@ -516,7 +518,8 @@ struct OnTheFlyTracker { continue; } const auto pdg = std::abs(mcParticle.pdgCode()); - if (pdg != kElectron && pdg != kMuonMinus && pdg != kPiPlus && pdg != kKPlus && pdg != kProton) { + if (pdg != kElectron && pdg != kMuonMinus && pdg != kPiPlus && pdg != kKPlus && pdg != kProton + && pdg != 1000010020 && pdg != 1000010030 && pdg != 1000020030 && pdg != 1000020040) { if (!cascadeDecaySettings.decayXi) { continue; } else if (pdg != 3312) { @@ -563,7 +566,8 @@ struct OnTheFlyTracker { continue; } } - if (pdg != kElectron && pdg != kMuonMinus && pdg != kPiPlus && pdg != kKPlus && pdg != kProton) { + if (pdg != kElectron && pdg != kMuonMinus && pdg != kPiPlus && pdg != kKPlus && pdg != kProton + && pdg != 1000010020 && pdg != 1000010030 && pdg != 1000020030 && pdg != 1000020040) { if (!cascadeDecaySettings.decayXi) { continue; } else if (pdg != 3312) { diff --git a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx index 42d2bcc5252..6c914890066 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx @@ -11,10 +11,10 @@ /// /// \file onTheFlyTrackerPid.cxx /// -/// \brief This task produces the PID information that can be obtained from the tracker layers (i.e. cluster size and ToT). -/// It currently contemplates 5 particle types: electrons, muons, pions, kaons and protons. +/// \brief This task produces the PID information that can be obtained from the tracker layers (i.e. ToT and possibly cluster size). +/// So far only ToT implemented. It currently contemplates 5 (9) particle types: electrons, muons, pions, kaons and +/// protons (as well as deuterons, tritons, helium-3 and alphas/helium-4 if added via the event generator). /// -/// \author Berkin Ulukutlu TUM /// \author Henrik Fribert TUM /// \author Nicolò Jacazio Università del Piemonte Orientale /// \since May 22, 2025 @@ -22,9 +22,15 @@ #include #include +#include #include #include #include +#include +#include +#include +#include +#include #include "Framework/AnalysisDataModel.h" #include "Framework/AnalysisTask.h" @@ -46,12 +52,10 @@ #include "DataFormatsCalibration/MeanVertexObject.h" #include "CommonConstants/GeomConstants.h" #include "CommonConstants/PhysicsConstants.h" -#include "TRandom3.h" #include "TF1.h" #include "TH2F.h" #include "TVector3.h" #include "TString.h" -#include "ALICE3/DataModel/OTFRICH.h" #include "DetectorsVertexing/HelixHelper.h" #include "TableHelper.h" #include "ALICE3/Core/DelphesO2TrackSmearer.h" @@ -60,116 +64,523 @@ using namespace o2; using namespace o2::framework; +static constexpr std::array kTrackerRadii = {0.5f, 1.2f, 2.5f, 3.75f, 7.0f, 12.0f, 20.0f, 30.0f, 45.0f, 60.0f, 80.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) { + mPdgToIndexMap.reserve(10); + mIndexToPdgMap.reserve(10); + } + ToTLUT() = delete; + + ~ToTLUT() { + for (auto hist_ptr : mLUTHistogramFlat) { + if (hist_ptr) { + delete hist_ptr; + } + } + } + + bool load(int pdg, const std::string& filename) { + TFile* f = TFile::Open(filename.c_str()); + if (!f || f->IsZombie()) { + std::cerr << "ERROR: Failed to open LUT file: " << filename << std::endl; + return false; + } + + int current_pdg_idx; + auto it = mPdgToIndexMap.find(pdg); + if (it == mPdgToIndexMap.end()) { + current_pdg_idx = mIndexToPdgMap.size(); + mPdgToIndexMap[pdg] = current_pdg_idx; + mIndexToPdgMap.push_back(pdg); + + size_t totalSize = (current_pdg_idx + 1) * mMaxLayers * mEtaBins * mPtBins; + if (mLUTHistogramFlat.size() < totalSize) { + mLUTHistogramFlat.resize(totalSize, nullptr); + } + } else { + current_pdg_idx = it->second; + } + + bool success = true; + for (int layer = 0; layer < mMaxLayers; ++layer) { + for (int etaBin = 0; etaBin < mEtaBins; ++etaBin) { + float etaMin_bin = mEtaMin + etaBin * mEtaBinWidth; + float etaMax_bin = etaMin_bin + mEtaBinWidth; + for (int ptBin = 0; ptBin < mPtBins; ++ptBin) { + float ptMin_bin = mPtMin + ptBin * mPtBinWidth; + float ptMax_bin = ptMin_bin + mPtBinWidth; + TString histName = Form("tot_%d_barrel%d_eta%.2f-%.2f_pt%.2f-%.2f", pdg, layer, etaMin_bin, etaMax_bin, ptMin_bin, ptMax_bin); + + TH1F* hist_from_file = dynamic_cast(f->Get(histName)); + if (hist_from_file) { + TH1F* cloned_hist = static_cast(hist_from_file->Clone()); + cloned_hist->SetDirectory(nullptr); + + size_t flatIdx = getFlatIndex(current_pdg_idx, layer, etaBin, ptBin); + mLUTHistogramFlat[flatIdx] = cloned_hist; + } else { + size_t flatIdx = getFlatIndex(current_pdg_idx, layer, etaBin, ptBin); + mLUTHistogramFlat[flatIdx] = nullptr; + success = false; + } + } + } + } + + f->Close(); + delete f; + return success; + } + + TH1F* getHistogramForSampling(int pdg_idx, int layer, int etaBin, int ptBin) const { + if (pdg_idx < 0 || static_cast(pdg_idx) >= getNumPdgTypes() || + layer < 0 || layer >= mMaxLayers || + etaBin < 0 || etaBin >= mEtaBins || ptBin < 0 || ptBin >= mPtBins) { + return nullptr; + } + size_t flatIdx = getFlatIndex(pdg_idx, layer, etaBin, ptBin); + return (flatIdx < mLUTHistogramFlat.size()) ? mLUTHistogramFlat[flatIdx] : nullptr; + } + + int getPdgIndex(int pdgCode) const { + auto it = mPdgToIndexMap.find(pdgCode); + if (it != mPdgToIndexMap.end()) { + return it->second; + } + return -1; + } + + inline int getEtaBin(float eta) const { + const float clampedEta = std::max(mEtaMin, std::min(eta, mEtaMax - 1e-6f)); + return std::min(static_cast((clampedEta - mEtaMin) / mEtaBinWidth), mEtaBins - 1); + } + + inline int getPtBin(float pt) const { + const float clampedPt = std::max(mPtMin, std::min(pt, mPtMax - 1e-6f)); + return std::min(static_cast((clampedPt - mPtMin) / mPtBinWidth), mPtBins - 1); + } + + inline size_t getFlatIndex(int pdgIdx, int layer, int etaBin, int ptBin) const { + return ((pdgIdx * mMaxLayers + layer) * mEtaBins + etaBin) * mPtBins + ptBin; + } + + size_t getNumPdgTypes() const { + return mIndexToPdgMap.size(); + } + + std::vector mLUTHistogramFlat; + + std::unordered_map mPdgToIndexMap; + std::vector mIndexToPdgMap; + + float mTruncationFraction; + int mMaxLayers; + int mEtaBins; + float mEtaMin; + float mEtaMax; + int mPtBins; + float mPtMin; + float mPtMax; + + float mEtaBinWidth; + float mPtBinWidth; +}; + +static constexpr int kNumHypothesisParticles = 9; +std::array, kNumHypothesisParticles>, kNumHypothesisParticles> h2dBarrelNsigmaTrue; +std::array, kNumHypothesisParticles> h2dHitsPerTrackVsP; +std::array, kNumHypothesisParticles> h2dToTvsP_PerParticle; + struct OnTheFlyTrackerPid { + + float calculateNsigma(float measuredToT, float expectedToT, float resolution) { + if (resolution <= 0) return 999.f; + return (measuredToT - expectedToT) / resolution; + } + + float getToTMeanFromMomentumSlice(std::shared_ptr hist, float momentum) { + if (!hist) return -1.f; + int binX = hist->GetXaxis()->FindBin(momentum); + TH1D* proj = hist->ProjectionY("temp", binX, binX); + if (proj->GetEntries() < 10) { + delete proj; + return -1.f; + } + float mean = proj->GetMean(); + delete proj; + return mean; + } + + float getToTResolutionFromMomentumSlice(std::shared_ptr hist, float momentum) { + if (!hist) return -1.f; + int binX = hist->GetXaxis()->FindBin(momentum); + TH1D* proj = hist->ProjectionY("temp", binX, binX); + if (proj->GetEntries() < 10) { + delete proj; + return -1.f; + } + float stddev = proj->GetStdDev(); + delete proj; + return stddev; + } + + float computeTrackLength(o2::track::TrackParCov track, float radius, float magneticField) + { + float length = -100; + o2::math_utils::CircleXYf_t trcCircle; + float sna, csa; + track.getCircleParams(magneticField, trcCircle, sna, csa); + + const float centerDistance = std::hypot(trcCircle.xC, trcCircle.yC); + + if (centerDistance < trcCircle.rC + radius && centerDistance > std::fabs(trcCircle.rC - radius)) { + length = 0.0f; + const float ux = trcCircle.xC / centerDistance; + const float uy = trcCircle.yC / centerDistance; + const float vx = -uy; + const float vy = +ux; + const float radical = (centerDistance * centerDistance - trcCircle.rC * trcCircle.rC + radius * radius) / (2.0f * centerDistance); + const float displace = (0.5f / centerDistance) * std::sqrt( + (-centerDistance + trcCircle.rC - radius) * + (-centerDistance - trcCircle.rC + radius) * + (-centerDistance + trcCircle.rC + radius) * + (centerDistance + trcCircle.rC + radius)); + + const float point1[2] = {radical * ux + displace * vx, radical * uy + displace * vy}; + const float point2[2] = {radical * ux - displace * vx, radical * uy - displace * vy}; + + std::array mom; + track.getPxPyPzGlo(mom); + const float scalarProduct1 = point1[0] * mom[0] + point1[1] * mom[1]; + const float scalarProduct2 = point2[0] * mom[0] + point2[1] * mom[1]; + + std::array startPoint; + track.getXYZGlo(startPoint); + + float cosAngle = -1000, modulus = -1000; + + if (scalarProduct1 > scalarProduct2) { + modulus = std::hypot(point1[0] - trcCircle.xC, point1[1] - trcCircle.yC) * std::hypot(startPoint[0] - trcCircle.xC, startPoint[1] - trcCircle.yC); + cosAngle = (point1[0] - trcCircle.xC) * (startPoint[0] - trcCircle.xC) + (point1[1] - trcCircle.yC) * (startPoint[1] - trcCircle.yC); + } else { + modulus = std::hypot(point2[0] - trcCircle.xC, point2[1] - trcCircle.yC) * std::hypot(startPoint[0] - trcCircle.xC, startPoint[1] - trcCircle.yC); + cosAngle = (point2[0] - trcCircle.xC) * (startPoint[0] - trcCircle.xC) + (point2[1] - trcCircle.yC) * (startPoint[1] - trcCircle.yC); + } + cosAngle /= modulus; + length = trcCircle.rC * std::acos(cosAngle); + length *= std::sqrt(1.0f + track.getTgl() * track.getTgl()); + } + return length; + } + Produces tableUpgradeTrkPidSignals; Produces tableUpgradeTrkPids; - // necessary for particle charges Service pdg; + std::unique_ptr mToTLUT; - static constexpr int kMaxBarrelLayers = 8; - static constexpr int kMaxForwardLayers = 9; + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; - struct : ConfigurableGroup { - Configurable efficiencyFormula{"efficiencyFormula", "1.0/(1.0+exp(-(x-0.01)/0.2))", "ROOT TF1 formula for efficiency"}; - Configurable landauFormula{"landauFormula", "TMath::Landau(x, 1, 1, true)", "ROOT TF1 formula for Landau distribution (e.g. ToT response)"}; - Configurable averageMethod{"averageMethod", 0, "Method to average the ToT and cluster size. 0: truncated mean"}; - } simConfig; + Configurable lutTotEl{"lutTotEl", "lut_tot_11.root", "ToT LUT for electrons"}; + Configurable lutTotMu{"lutTotMu", "lut_tot_13.root", "ToT LUT for muons"}; + Configurable lutTotPi{"lutTotPi", "lut_tot_211.root", "ToT LUT for pions"}; + Configurable lutTotKa{"lutTotKa", "lut_tot_321.root", "ToT LUT for kaons"}; + Configurable lutTotPr{"lutTotPr", "lut_tot_2212.root", "ToT LUT for protons"}; + Configurable lutTotDe{"lutTotDe", "lut_tot_1000010020.root", "ToT LUT for deuteron"}; + Configurable lutTotTr{"lutTotTr", "lut_tot_1000010030.root", "ToT LUT for triton"}; + Configurable lutTotHe{"lutTotHe", "lut_tot_1000020030.root", "ToT LUT for helium-3"}; + Configurable lutTotAl{"lutTotAl", "lut_tot_1000020040.root", "ToT LUT for alphas"}; - TF1* mEfficiency = nullptr; - static constexpr int kEtaBins = 50; - static constexpr float kEtaMin = -2.5; - static constexpr float kEtaMax = 2.5; - static constexpr int kPtBins = 200; - static constexpr float kPtMin = 0.0; - static constexpr float kPtMax = 20.0; + 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 mMaxBarrelLayers{"maxBarrelLayers", 11, "Maximum number of barrel layers"}; + Configurable mEtaBins{"etaBins", 50, "Number of eta bins for LUTs and histograms"}; + Configurable mEtaMin{"etaMin", -2.5f, "Minimum eta value"}; + Configurable mEtaMax{"etaMax", 2.5f, "Maximum eta value"}; + Configurable mPtBins{"ptBins", 500, "Number of pT bins for LUTs (LUTs are pT-based)"}; + Configurable mPtMin{"ptMin", 0.0f, "Minimum pT value for LUT binning"}; + Configurable mPtMax{"ptMax", 10.0f, "Maximum pT value for LUT binning"}; + Configurable mNumLogBins{"numLogBins", 200, "Number of logarithmic momentum bins"}; - std::array, kEtaBins> mElossPi; + std::vector mLogBins; - void init(o2::framework::InitContext&) - { + std::array mHypothesisPdgCodes = { + 11, // Electron + 13, // Muon + 211, // Pion + 321, // Kaon + 2212, // Proton + 1000010020, // Deuteron + 1000010030, // Triton + 1000020030, // Helium-3 + 1000020040 // Alpha + }; - for (int i = 0; i < kEtaBins; i++) { - for (int j = 0; j < kPtBins; j++) { - mElossPi[i][j] = new TF1(Form("mElossPi_%d_%d", i, j), simConfig.landauFormula.value.c_str(), 0, 20); - } + void init(o2::framework::InitContext&) { + if (static_cast(mMaxBarrelLayers.value) > kTrackerRadii.size()) { + LOG(fatal) << "Configured maxBarrelLayers (" << mMaxBarrelLayers.value + << ") exceeds the size of kTrackerRadii (" << kTrackerRadii.size() + << "). Please adjust maxBarrelLayers."; + } + + mToTLUT = std::make_unique(truncationFraction.value, mMaxBarrelLayers.value, + mEtaBins.value, mEtaMin.value, mEtaMax.value, + mPtBins.value, mPtMin.value, mPtMax.value); + + bool loaded = true; + loaded &= mToTLUT->load(11, lutTotEl.value); + loaded &= mToTLUT->load(13, lutTotMu.value); + loaded &= mToTLUT->load(211, lutTotPi.value); + loaded &= mToTLUT->load(321, lutTotKa.value); + loaded &= mToTLUT->load(2212, lutTotPr.value); + loaded &= mToTLUT->load(1000010020, lutTotDe.value); + loaded &= mToTLUT->load(1000010030, lutTotTr.value); + loaded &= mToTLUT->load(1000020030, lutTotHe.value); + loaded &= mToTLUT->load(1000020040, lutTotAl.value); + + if (!loaded) { + std::cerr << "WARNING: Failed to load one or more ToT LUTs. PID results might be incomplete." << std::endl; + } + + // Logarithmic momentum bins + mLogBins.clear(); + double pMin = 0.05; + double pMax = 10; + double logMin = std::log10(pMin); + double logMax = std::log10(pMax); + double dLog = (logMax - logMin) / mNumLogBins.value; + for (int i = 0; i <= mNumLogBins.value; ++i) { + mLogBins.push_back(std::pow(10, logMin + i * dLog)); + } + + const AxisSpec axisMomentum{mLogBins, "#it{p/z} (GeV/#it{c})"}; + const AxisSpec axisToT{600, 0., 300., "ToT (#mus/10#mum)"}; + const AxisSpec axisNsigma{200, -10., 10., "N#sigma"}; + const AxisSpec axisLayer{mMaxBarrelLayers.value, -0.5, static_cast(mMaxBarrelLayers.value) - 0.5, "Layer"}; + const AxisSpec axisHitsPerTrack{mMaxBarrelLayers.value + 1, -0.5, static_cast(mMaxBarrelLayers.value) + 0.5, "# Hits per track"}; + + histos.add("hToTvsP", "ToT vs #it{p/z}; #it{p/z} (GeV/#it{c}); ToT (#mus/10#mum)", kTH2F, {axisMomentum, axisToT}); + histos.add("hToTvsPt", "ToT vs #it{p}; #it{p} (GeV/#it{c}); ToT (#mus/10#mum)", kTH2F, {mLogBins, axisToT}); + histos.add("hHitLayers", "Number of hits on each detector layer;Layer;Counts", kTH1F, {axisLayer}); + histos.add("hHitMultiplicity", "Hit multiplicity along the track; # Hits per track;Counts", kTH1F, {axisHitsPerTrack}); + + std::vector> particleInfo = { + {11, "Elec"}, {13, "Muon"}, {211, "Pion"}, {321, "Kaon"}, {2212, "Prot"}, + {1000010020, "Deut"}, {1000010030, "Trit"}, {1000020030, "He3"}, {1000020040, "Al"} + }; + + for (size_t i_true = 0; i_true < particleInfo.size(); ++i_true) { + std::string trueName = particleInfo[i_true].second; + std::string trueNamePretty = trueName; // Fallback + if (trueName == "Elec") trueNamePretty = "#it{e}"; + else if (trueName == "Muon") trueNamePretty = "#it{#mu}"; + else if (trueName == "Pion") trueNamePretty = "#it{#pi}"; + else if (trueName == "Kaon") trueNamePretty = "#it{K}"; + else if (trueName == "Prot") trueNamePretty = "#it{p}"; + else if (trueName == "Deut") trueNamePretty = "#it{d}"; + else if (trueName == "Trit") trueNamePretty = "#it{t}"; + else if (trueName == "He3") trueNamePretty = "#it{^{3}He}"; + else if (trueName == "Al") trueNamePretty = "#it{^{4}He}"; + + std::string hitsVsPName = "HitsPerTrack/hHitsPerTrackVsP_" + trueName; + std::string hitsVsPTitle = "N_hits vs #it{p/z} for " + trueNamePretty + "; #it{p/z} (GeV/#it{c}); N_hits"; + h2dHitsPerTrackVsP[i_true] = histos.add(hitsVsPName.c_str(), hitsVsPTitle.c_str(), kTH2F, {axisMomentum, axisHitsPerTrack}); + + std::string totVsPName = "ToTvsP/hToTvsP_" + trueName; + std::string totVsPTitle = "ToT vs #it{p/z} for " + trueNamePretty + "; #it{p/z} (GeV/#it{c}); ToT (#mus/10#mum)"; + h2dToTvsP_PerParticle[i_true] = histos.add(totVsPName.c_str(), totVsPTitle.c_str(), kTH2F, {axisMomentum, axisToT}); + + + for (size_t i_hyp = 0; i_hyp < particleInfo.size(); ++i_hyp) { + std::string hypName = particleInfo[i_hyp].second; + std::string hypNamePretty = hypName; // Fallback + if (hypName == "Elec") hypNamePretty = "#it{e}"; + else if (hypName == "Muon") hypNamePretty = "#it{#mu}"; + else if (hypName == "Pion") hypNamePretty = "#it{#pi}"; + else if (hypName == "Kaon") hypNamePretty = "#it{K}"; + else if (hypName == "Prot") hypNamePretty = "#it{p}"; + else if (hypName == "Deut") hypNamePretty = "#it{d}"; + else if (hypName == "Trit") hypNamePretty = "#it{t}"; + else if (hypName == "He3") hypNamePretty = "#it{^{3}He}"; + else if (hypName == "Al") hypNamePretty = "#it{^{4}He}"; + + std::string histName = "NSigma/BarrelNsigmaTrue" + trueName + "Vs" + hypName + "Hypothesis"; + std::string histTitle = "Nsigma (True " + trueNamePretty + " vs Hyp " + hypNamePretty + "); #it{p/z} (GeV/#it{c}); N#sigma"; + h2dBarrelNsigmaTrue[i_true][i_hyp] = histos.add(histName.c_str(), histTitle.c_str(), kTH2F, {axisMomentum, axisNsigma}); + } } - mEfficiency = new TF1("mEfficiency", simConfig.efficiencyFormula.value.c_str(), 0, 20); } - void process(soa::Join::iterator const&, + void process(soa::Join::iterator const& collision, soa::Join const& tracks, - aod::McParticles const&, - aod::McCollisions const&) - { - std::array timeOverThresholdBarrel; - std::array clusterSizeBarrel; - // std::array timeOverThresholdForward; - // std::array clusterSizeForward; - - auto noSignalTrack = [&]() { - tableUpgradeTrkPidSignals(0.f, 0.f); // no PID information - tableUpgradeTrkPids(0.f, 0.f, 0.f, 0.f, 0.f); // no PID information - }; + aod::McParticles const& /*mcParticles*/, + aod::McCollisions const& /*mcCollisions*/) { + + o2::dataformats::VertexBase mcPvVtx({0.0f, 0.0f, 0.0f}, {0.}); + + if (collision.has_mcCollision()) { + const auto& mcCollisionObject = collision.mcCollision(); + mcPvVtx.setX(mcCollisionObject.posX()); + mcPvVtx.setY(mcCollisionObject.posY()); + mcPvVtx.setZ(mcCollisionObject.posZ()); + } + + int nTracksProcessed = 0; + int nValidToT = 0; + + for (const auto& track : tracks) { + nTracksProcessed++; + float truncatedMeanToT = -1.0f; + std::array nSigmaValues; + nSigmaValues.fill(999.f); - for (const auto& track : tracks) { if (!track.has_mcParticle()) { - noSignalTrack(); + tableUpgradeTrkPidSignals(truncatedMeanToT); + tableUpgradeTrkPids(nSigmaValues[0], nSigmaValues[1], nSigmaValues[2], nSigmaValues[3], + nSigmaValues[4], nSigmaValues[5], nSigmaValues[6], nSigmaValues[7], nSigmaValues[8]); continue; } + const auto& mcParticle = track.mcParticle(); + const auto& pdgInfo = pdg->GetParticle(mcParticle.pdgCode()); if (!pdgInfo) { - LOG(warning) << "PDG code " << mcParticle.pdgCode() << " not found in the database"; - noSignalTrack(); + tableUpgradeTrkPidSignals(truncatedMeanToT); + tableUpgradeTrkPids(nSigmaValues[0], nSigmaValues[1], nSigmaValues[2], nSigmaValues[3], + nSigmaValues[4], nSigmaValues[5], nSigmaValues[6], nSigmaValues[7], nSigmaValues[8]); continue; } + const float pt = mcParticle.pt(); + const float p = mcParticle.p(); const float eta = mcParticle.eta(); + const int truePdgCode = std::abs(mcParticle.pdgCode()); - const int binnedPt = static_cast((pt - kPtMin) / kPtBins); - const int binnedEta = static_cast((eta - kEtaMin) / kEtaBins); - if (binnedPt < 0 || binnedPt >= kPtBins || binnedEta < 0 || binnedEta >= kEtaBins) { - noSignalTrack(); - continue; + float rigidity = p; + if (pdgInfo) { + const float trueCharge = std::abs(pdgInfo->Charge()) / 3.0f; + if (trueCharge > 0) { + rigidity = p / trueCharge; + } + } + + int true_pdg_idx = mToTLUT->getPdgIndex(truePdgCode); + if (true_pdg_idx == -1) { + tableUpgradeTrkPidSignals(truncatedMeanToT); + tableUpgradeTrkPids(nSigmaValues[0], nSigmaValues[1], nSigmaValues[2], nSigmaValues[3], + nSigmaValues[4], nSigmaValues[5], nSigmaValues[6], nSigmaValues[7], nSigmaValues[8]); + continue; + } + + const int binnedPt = mToTLUT->getPtBin(pt); + const int binnedEta = mToTLUT->getEtaBin(std::abs(eta)); + + uint16_t hitMap = 0; + int nHitLayers = 0; + o2::track::TrackParCov o2track = o2::upgrade::convertMCParticleToO2Track(mcParticle, pdg); + + float xPv = -100.f; + static constexpr float kTrkXThreshold = -99.f; + if (o2track.propagateToDCA(mcPvVtx, dBz)) { + xPv = o2track.getX(); + } + + if (xPv > kTrkXThreshold) { + for (int layer = 0; layer < mMaxBarrelLayers.value; ++layer) { + float layerRadius = kTrackerRadii[layer]; + float trackLength = computeTrackLength(o2track, layerRadius, dBz); + + if (trackLength > 0) { + hitMap |= (1 << layer); + histos.fill(HIST("hHitLayers"), layer); + nHitLayers++; + } + } + } + + histos.fill(HIST("hHitMultiplicity"), nHitLayers); + h2dHitsPerTrackVsP[true_pdg_idx]->Fill(rigidity, nHitLayers); + + std::vector validToTs; + + for (int layer = 3; layer < mMaxBarrelLayers.value; ++layer) { + if ((hitMap >> layer) & 0x1) { + TH1F* totHist = mToTLUT->getHistogramForSampling(true_pdg_idx, layer, binnedEta, binnedPt); + + if (totHist && totHist->GetEntries() > 1) { + float sampledToT = totHist->GetRandom(); + validToTs.push_back(sampledToT); + } + } } - for (int i = 0; i < kMaxBarrelLayers; i++) { - timeOverThresholdBarrel[i] = -1; - clusterSizeBarrel[i] = -1; - - // Check if layer is efficient - if (mEfficiency->Eval(pt) > gRandom->Uniform(0, 1)) { - timeOverThresholdBarrel[i] = mElossPi[binnedEta][binnedPt]->GetRandom(); // Simulate ToT - clusterSizeBarrel[i] = mElossPi[binnedEta][binnedPt]->GetRandom(); // Simulate cluster size + + truncatedMeanToT = -1.0f; + const size_t nValid = validToTs.size(); + size_t nUse = 0; + + if (nValid == 3 || nValid == 4) nUse = 2; + else if (nValid == 1 || nValid == 2) nUse = 1; + else if (nValid == 5 || nValid == 6) nUse = 3; + else if (nValid >= 7) nUse = 4; + + if (nUse > 0 && nValid >= nUse) { + nValidToT++; + std::sort(validToTs.begin(), validToTs.end()); + float sum = 0.0f; + for (size_t i = 0; i < nUse; ++i) { + sum += validToTs[i]; } + truncatedMeanToT = sum / static_cast(nUse); + + histos.fill(HIST("hToTvsPt"), p, truncatedMeanToT); + histos.fill(HIST("hToTvsP"), rigidity, truncatedMeanToT); + h2dToTvsP_PerParticle[true_pdg_idx]->Fill(rigidity, truncatedMeanToT); } - // Now we do the average - switch (simConfig.averageMethod) { - case 0: { // truncated mean - float meanToT = 0; - float meanClusterSize = 0; - // Order them by ToT - std::sort(timeOverThresholdBarrel.begin(), timeOverThresholdBarrel.end()); - std::sort(clusterSizeBarrel.begin(), clusterSizeBarrel.end()); - static constexpr int kTruncatedMean = 5; - // Take the mean of the first 5 values - for (int i = 0; i < kTruncatedMean; i++) { - meanToT += timeOverThresholdBarrel[i]; - meanClusterSize += clusterSizeBarrel[i]; + nSigmaValues.fill(999.f); + + if (truncatedMeanToT > 0) { + for (size_t i_hyp = 0; i_hyp < mHypothesisPdgCodes.size(); ++i_hyp) { + int hypPdgCode = mHypothesisPdgCodes[i_hyp]; + int hyp_pdg_idx = mToTLUT->getPdgIndex(hypPdgCode); + + if (hyp_pdg_idx == -1) { + nSigmaValues[i_hyp] = 999.f; + continue; } - meanToT /= kTruncatedMean; - meanClusterSize /= kTruncatedMean; - // Fill the table - tableUpgradeTrkPidSignals(meanToT, meanClusterSize); - } break; - - default: - LOG(fatal) << "Unknown average method " << simConfig.averageMethod; - break; + + float expectedToT = getToTMeanFromMomentumSlice(h2dToTvsP_PerParticle[hyp_pdg_idx], rigidity); + float resolution = getToTResolutionFromMomentumSlice(h2dToTvsP_PerParticle[hyp_pdg_idx], rigidity); + + if (expectedToT > 0 && resolution > 0) { + nSigmaValues[i_hyp] = calculateNsigma(truncatedMeanToT, expectedToT, resolution); + h2dBarrelNsigmaTrue[true_pdg_idx][i_hyp]->Fill(rigidity, nSigmaValues[i_hyp]); + } else { + nSigmaValues[i_hyp] = 999.f; + } + } } + + tableUpgradeTrkPidSignals(truncatedMeanToT); + tableUpgradeTrkPids(nSigmaValues[0], nSigmaValues[1], nSigmaValues[2], nSigmaValues[3], + nSigmaValues[4], nSigmaValues[5], nSigmaValues[6], nSigmaValues[7], nSigmaValues[8]); } } }; -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask(cfgc)}; } +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +} \ No newline at end of file From 2018ef3e61d2270db722217c4455ea89a23d0508 Mon Sep 17 00:00:00 2001 From: Henrik Fribert Date: Tue, 23 Sep 2025 17:11:34 +0200 Subject: [PATCH 02/23] onTheFlyTracker changes --- ALICE3/TableProducer/OTF/onTheFlyTracker.cxx | 86 ++++++++++---------- 1 file changed, 41 insertions(+), 45 deletions(-) diff --git a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx index cbb187b8359..0f785534159 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx @@ -113,7 +113,6 @@ struct OnTheFlyTracker { Configurable lutDe{"lutDe", "lutCovm.de.dat", "LUT for deuterons"}; Configurable lutTr{"lutTr", "lutCovm.tr.dat", "LUT for tritons"}; Configurable lutHe3{"lutHe3", "lutCovm.he3.dat", "LUT for Helium-3"}; - Configurable lutAl{"lutAl", "lutCovm.he3.dat", "LUT for Alphas"}; // To be created, for now propagating as He-3 struct : ConfigurableGroup { ConfigurableAxis axisMomentum{"axisMomentum", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f, 1.7f, 1.8f, 1.9f, 2.0f, 2.2f, 2.4f, 2.6f, 2.8f, 3.0f, 3.2f, 3.4f, 3.6f, 3.8f, 4.0f, 4.4f, 4.8f, 5.2f, 5.6f, 6.0f, 6.5f, 7.0f, 7.5f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 17.0f, 19.0f, 21.0f, 23.0f, 25.0f, 30.0f, 35.0f, 40.0f, 50.0f}, "#it{p} (GeV/#it{c})"}; @@ -139,11 +138,12 @@ struct OnTheFlyTracker { Configurable minSiliconHits{"minSiliconHits", 6, "minimum number of silicon hits to accept track"}; Configurable minSiliconHitsIfTPCUsed{"minSiliconHitsIfTPCUsed", 2, "minimum number of silicon hits to accept track in case TPC info is present"}; Configurable minTPCClusters{"minTPCClusters", 70, "minimum number of TPC hits necessary to consider minSiliconHitsIfTPCUsed"}; - Configurable alice3detector{"alice3detector", 0, "0: ALICE 3 v1, 1: ALICE 3 v4"}; + Configurable alice3detector{"alice3detector", 2, "0: ALICE 3 v1, 1: ALICE 3 v4, 2: ALICE 3 Sep 2025"}; Configurable applyZacceptance{"applyZacceptance", false, "apply z limits to detector layers or not"}; Configurable applyMSCorrection{"applyMSCorrection", true, "apply ms corrections for secondaries or not"}; Configurable applyElossCorrection{"applyElossCorrection", true, "apply eloss corrections for secondaries or not"}; Configurable applyEffCorrection{"applyEffCorrection", true, "apply efficiency correction or not"}; + Configurable scaleVD{"scaleVD", 1, "scale x0 and xrho in VD layers"}; Configurable> pixelRes{"pixelRes", {0.00025, 0.00025, 0.001, 0.001}, "RPhiIT, ZIT, RPhiOT, ZOT"}; } fastTrackerSettings; // allows for gap between peak and bg in case someone wants to @@ -240,44 +240,32 @@ struct OnTheFlyTracker { // For TGenPhaseSpace seed TRandom3 rand; + Service ccdb; void init(o2::framework::InitContext&) { + + ccdb->setURL("http://alice-ccdb.cern.ch"); + ccdb->setTimestamp(-1); + if (enableLUT) { - std::map mapPdgLut; - const char* lutElChar = lutEl->c_str(); - const char* lutMuChar = lutMu->c_str(); - const char* lutPiChar = lutPi->c_str(); - const char* lutKaChar = lutKa->c_str(); - const char* lutPrChar = lutPr->c_str(); - - LOGF(info, "Will load electron lut file ..: %s", lutElChar); - LOGF(info, "Will load muon lut file ......: %s", lutMuChar); - LOGF(info, "Will load pion lut file ......: %s", lutPiChar); - LOGF(info, "Will load kaon lut file ......: %s", lutKaChar); - LOGF(info, "Will load proton lut file ....: %s", lutPrChar); - - mapPdgLut.insert(std::make_pair(11, lutElChar)); - mapPdgLut.insert(std::make_pair(13, lutMuChar)); - mapPdgLut.insert(std::make_pair(211, lutPiChar)); - mapPdgLut.insert(std::make_pair(321, lutKaChar)); - mapPdgLut.insert(std::make_pair(2212, lutPrChar)); - - if (enableNucleiSmearing) { - const char* lutDeChar = lutDe->c_str(); - const char* lutTrChar = lutTr->c_str(); - const char* lutHe3Char = lutHe3->c_str(); - const char* lutAlChar = lutAl->c_str(); - mapPdgLut.insert(std::make_pair(1000010020, lutDeChar)); - mapPdgLut.insert(std::make_pair(1000010030, lutTrChar)); - mapPdgLut.insert(std::make_pair(1000020030, lutHe3Char)); - mapPdgLut.insert(std::make_pair(1000020030, lutAlChar)); - } - for (const auto& e : mapPdgLut) { - if (!mSmearer.loadTable(e.first, e.second)) { - LOG(fatal) << "Having issue with loading the LUT " << e.first << " " << e.second; + mSmearer.setCcdbManager(ccdb.operator->()); + + auto loadLUT = [&](int pdg, const std::string& lutFile) { + bool success = mSmearer.loadTable(pdg, lutFile.c_str()); + if (!success && !lutFile.empty()) { + LOG(fatal) << "Having issue with loading the LUT " << pdg << " " << lutFile; } - } + }; + loadLUT(11, lutEl.value); + loadLUT(13, lutMu.value); + loadLUT(211, lutPi.value); + loadLUT(321, lutKa.value); + loadLUT(2212, lutPr.value); + loadLUT(1000010020, lutDe.value); + loadLUT(1000010030, lutTr.value); + loadLUT(1000020030, lutHe3.value); + // interpolate efficiencies if requested to do so mSmearer.interpolateEfficiency(static_cast(interpolateLutEfficiencyVsNch)); @@ -425,12 +413,22 @@ struct OnTheFlyTracker { fastTracker.SetApplyMSCorrection(fastTrackerSettings.applyMSCorrection); fastTracker.SetApplyElossCorrection(fastTrackerSettings.applyElossCorrection); - if (fastTrackerSettings.alice3detector == 0) { - fastTracker.AddSiliconALICE3v2(fastTrackerSettings.pixelRes); - } - if (fastTrackerSettings.alice3detector == 1) { - fastTracker.AddSiliconALICE3v4(fastTrackerSettings.pixelRes); - fastTracker.AddTPC(0.1, 0.1); + switch (fastTrackerSettings.alice3detector) { + case 0: + fastTracker.AddSiliconALICE3v2(fastTrackerSettings.pixelRes); + break; + + case 1: + fastTracker.AddSiliconALICE3v4(fastTrackerSettings.pixelRes); + fastTracker.AddTPC(0.1, 0.1); + break; + + case 2: + fastTracker.AddSiliconALICE3(fastTrackerSettings.scaleVD, fastTrackerSettings.pixelRes); + break; + + default: + break; } // print fastTracker settings @@ -518,8 +516,7 @@ struct OnTheFlyTracker { continue; } const auto pdg = std::abs(mcParticle.pdgCode()); - if (pdg != kElectron && pdg != kMuonMinus && pdg != kPiPlus && pdg != kKPlus && pdg != kProton - && pdg != 1000010020 && pdg != 1000010030 && pdg != 1000020030 && pdg != 1000020040) { + if (pdg != kElectron && pdg != kMuonMinus && pdg != kPiPlus && pdg != kKPlus && pdg != kProton) { if (!cascadeDecaySettings.decayXi) { continue; } else if (pdg != 3312) { @@ -566,8 +563,7 @@ struct OnTheFlyTracker { continue; } } - if (pdg != kElectron && pdg != kMuonMinus && pdg != kPiPlus && pdg != kKPlus && pdg != kProton - && pdg != 1000010020 && pdg != 1000010030 && pdg != 1000020030 && pdg != 1000020040) { + if (pdg != kElectron && pdg != kMuonMinus && pdg != kPiPlus && pdg != kKPlus && pdg != kProton) { if (!cascadeDecaySettings.decayXi) { continue; } else if (pdg != 3312) { From 3fdc82ed067f188643c2cf8d2b2142a99037415e Mon Sep 17 00:00:00 2001 From: Henrik Fribert Date: Tue, 23 Sep 2025 17:39:23 +0200 Subject: [PATCH 03/23] Fix: Linter and Clang-formatting --- ALICE3/DataModel/OTFPIDTrk.h | 21 +- ALICE3/TableProducer/OTF/onTheFlyTracker.cxx | 6 +- .../TableProducer/OTF/onTheFlyTrackerPid.cxx | 489 ++++++++++-------- 3 files changed, 289 insertions(+), 227 deletions(-) diff --git a/ALICE3/DataModel/OTFPIDTrk.h b/ALICE3/DataModel/OTFPIDTrk.h index f0a66496ee5..ce8a5d07a4f 100644 --- a/ALICE3/DataModel/OTFPIDTrk.h +++ b/ALICE3/DataModel/OTFPIDTrk.h @@ -31,15 +31,15 @@ namespace upgrade::trk DECLARE_SOA_COLUMN(TimeOverThresholdBarrel, timeOverThresholdBarrel, float); //! Time over threshold for the Barrel layers DECLARE_SOA_COLUMN(TimeOverThresholdForward, timeOverThresholdForward, float); //! Time over threshold for the Forward layers -DECLARE_SOA_COLUMN(NSigmaTrkEl, nSigmaEl, float); //! NSigma electron from the tracker layers -DECLARE_SOA_COLUMN(NSigmaTrkMu, nSigmaMu, float); //! NSigma muon from the tracker layers -DECLARE_SOA_COLUMN(NSigmaTrkPi, nSigmaPi, float); //! NSigma pion from the tracker layers -DECLARE_SOA_COLUMN(NSigmaTrkKa, nSigmaKa, float); //! NSigma kaon from the tracker layers -DECLARE_SOA_COLUMN(NSigmaTrkPr, nSigmaPr, float); //! NSigma proton from the tracker layers -DECLARE_SOA_COLUMN(NSigmaTrkDe, nSigmaDe, float); //! NSigma deuteron from the tracker layers -DECLARE_SOA_COLUMN(NSigmaTrkTr, nSigmaTr, float); //! NSigma triton from the tracker layers -DECLARE_SOA_COLUMN(NSigmaTrkHe, nSigmaHe, float); //! NSigma helium-3 from the tracker layers -DECLARE_SOA_COLUMN(NSigmaTrkAl, nSigmaAl, float); //! NSigma alpha from the tracker layers +DECLARE_SOA_COLUMN(NSigmaTrkEl, nSigmaTrkEl, float); //! NSigma electron from the tracker layers +DECLARE_SOA_COLUMN(NSigmaTrkMu, nSigmaTrkMu, float); //! NSigma muon from the tracker layers +DECLARE_SOA_COLUMN(NSigmaTrkPi, nSigmaTrkPi, float); //! NSigma pion from the tracker layers +DECLARE_SOA_COLUMN(NSigmaTrkKa, nSigmaTrkKa, float); //! NSigma kaon from the tracker layers +DECLARE_SOA_COLUMN(NSigmaTrkPr, nSigmaTrkPr, float); //! NSigma proton from the tracker layers +DECLARE_SOA_COLUMN(NSigmaTrkDe, nSigmaTrkDe, float); //! NSigma deuteron from the tracker layers +DECLARE_SOA_COLUMN(NSigmaTrkTr, nSigmaTrkTr, float); //! NSigma triton from the tracker layers +DECLARE_SOA_COLUMN(NSigmaTrkHe, nSigmaTrkHe, float); //! NSigma helium-3 from the tracker layers +DECLARE_SOA_COLUMN(NSigmaTrkAl, nSigmaTrkAl, float); //! NSigma alpha from the tracker layers DECLARE_SOA_DYNAMIC_COLUMN(NSigmaTrk, nSigmaTrk, //! General function to get the nSigma for the tracker layers [](const float el, @@ -102,8 +102,7 @@ DECLARE_SOA_TABLE(UpgradeTrkPids, "AOD", "UPGRADETRKPID", upgrade::trk::NSigmaTrkDe, upgrade::trk::NSigmaTrkTr, upgrade::trk::NSigmaTrkHe, - upgrade::trk::NSigmaTrkAl - >); + upgrade::trk::NSigmaTrkAl>); using UpgradeTrkPidSignal = UpgradeTrkPidSignals::iterator; using UpgradeTrkPid = UpgradeTrkPids::iterator; diff --git a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx index 0f785534159..461456409e7 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx @@ -113,6 +113,7 @@ struct OnTheFlyTracker { Configurable lutDe{"lutDe", "lutCovm.de.dat", "LUT for deuterons"}; Configurable lutTr{"lutTr", "lutCovm.tr.dat", "LUT for tritons"}; Configurable lutHe3{"lutHe3", "lutCovm.he3.dat", "LUT for Helium-3"}; + Configurable lutAl{"lutAl", "lutCovm.he3.dat", "LUT for Alphas"}; // To be created, for now propagating as He-3 struct : ConfigurableGroup { ConfigurableAxis axisMomentum{"axisMomentum", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f, 1.7f, 1.8f, 1.9f, 2.0f, 2.2f, 2.4f, 2.6f, 2.8f, 3.0f, 3.2f, 3.4f, 3.6f, 3.8f, 4.0f, 4.4f, 4.8f, 5.2f, 5.6f, 6.0f, 6.5f, 7.0f, 7.5f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 17.0f, 19.0f, 21.0f, 23.0f, 25.0f, 30.0f, 35.0f, 40.0f, 50.0f}, "#it{p} (GeV/#it{c})"}; @@ -265,6 +266,7 @@ struct OnTheFlyTracker { loadLUT(1000010020, lutDe.value); loadLUT(1000010030, lutTr.value); loadLUT(1000020030, lutHe3.value); + loadLUT(1000020040, lutHe3.value); // interpolate efficiencies if requested to do so mSmearer.interpolateEfficiency(static_cast(interpolateLutEfficiencyVsNch)); @@ -516,7 +518,7 @@ struct OnTheFlyTracker { continue; } const auto pdg = std::abs(mcParticle.pdgCode()); - if (pdg != kElectron && pdg != kMuonMinus && pdg != kPiPlus && pdg != kKPlus && pdg != kProton) { + if (pdg != kElectron && pdg != kMuonMinus && pdg != kPiPlus && pdg != kKPlus && pdg != kProton && pdg != 1000010020 && pdg != 1000010030 && pdg != 1000020030 && pdg != 1000020040) { if (!cascadeDecaySettings.decayXi) { continue; } else if (pdg != 3312) { @@ -563,7 +565,7 @@ struct OnTheFlyTracker { continue; } } - if (pdg != kElectron && pdg != kMuonMinus && pdg != kPiPlus && pdg != kKPlus && pdg != kProton) { + if (pdg != kElectron && pdg != kMuonMinus && pdg != kPiPlus && pdg != kKPlus && pdg != kProton && pdg != 1000010020 && pdg != 1000010030 && pdg != 1000020030 && pdg != 1000020040) { if (!cascadeDecaySettings.decayXi) { continue; } else if (pdg != 3312) { diff --git a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx index 6c914890066..14e05efd242 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx @@ -20,119 +20,146 @@ /// \since May 22, 2025 /// -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include +#include "TableHelper.h" +#include "ALICE3/Core/DelphesO2TrackSmearer.h" +#include "ALICE3/Core/TrackUtilities.h" +#include "ALICE3/DataModel/OTFPIDTrk.h" +#include "Common/Core/trackUtilities.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include "CCDB/BasicCCDBManager.h" +#include "CCDB/CcdbApi.h" +#include "CommonConstants/GeomConstants.h" +#include "CommonConstants/PhysicsConstants.h" +#include "CommonUtils/NameConf.h" +#include "DataFormatsCalibration/MeanVertexObject.h" +#include "DataFormatsParameters/GRPMagField.h" +#include "DetectorsBase/GeometryManager.h" +#include "DetectorsBase/Propagator.h" +#include "DetectorsVertexing/HelixHelper.h" +#include "Framework/ASoAHelpers.h" #include "Framework/AnalysisDataModel.h" #include "Framework/AnalysisTask.h" -#include "Framework/runDataProcessing.h" -#include "Framework/RunningWorkflowInfo.h" #include "Framework/HistogramRegistry.h" #include "Framework/O2DatabasePDGPlugin.h" -#include "Framework/ASoAHelpers.h" -#include "Common/DataModel/TrackSelectionTables.h" -#include "Common/Core/trackUtilities.h" -#include "ALICE3/Core/TrackUtilities.h" +#include "Framework/RunningWorkflowInfo.h" +#include "Framework/runDataProcessing.h" #include "ReconstructionDataFormats/DCA.h" -#include "DetectorsBase/Propagator.h" -#include "DetectorsBase/GeometryManager.h" -#include "CommonUtils/NameConf.h" -#include "CCDB/CcdbApi.h" -#include "CCDB/BasicCCDBManager.h" -#include "DataFormatsParameters/GRPMagField.h" -#include "DataFormatsCalibration/MeanVertexObject.h" -#include "CommonConstants/GeomConstants.h" -#include "CommonConstants/PhysicsConstants.h" + #include "TF1.h" #include "TH2F.h" -#include "TVector3.h" #include "TString.h" -#include "DetectorsVertexing/HelixHelper.h" -#include "TableHelper.h" -#include "ALICE3/Core/DelphesO2TrackSmearer.h" -#include "ALICE3/DataModel/OTFPIDTrk.h" +#include "TVector3.h" +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include using namespace o2; using namespace o2::framework; static constexpr std::array kTrackerRadii = {0.5f, 1.2f, 2.5f, 3.75f, 7.0f, 12.0f, 20.0f, 30.0f, 45.0f, 60.0f, 80.0f}; -class ToTLUT { -public: +// Constants for magic numbers +static constexpr int kMinEntriesForProjection = 10; +static constexpr double kMinMomentumForLogBins = 0.05; +static constexpr double kMaxMomentumForLogBins = 10.0; +static constexpr int kMinLayerForTruncation = 3; + +// Constants for validation conditions +static constexpr size_t kMinValidHits = 1; +static constexpr size_t kLowValidHits = 2; +static constexpr size_t kMidLowValidHits = 3; +static constexpr size_t kMidHighValidHits = 4; +static constexpr size_t kHighValidHits1 = 5; +static constexpr size_t kHighValidHits2 = 6; +static constexpr size_t kMaxValidHits = 7; + +// Constants for truncated mean calculation +static constexpr size_t kMaxValidHitsForTruncation7Plus = 4; +static constexpr size_t kMaxValidHitsForTruncation56 = 3; +static constexpr size_t kMaxValidHitsForTruncation34 = 2; +static constexpr size_t kMaxValidHitsForTruncation12 = 1; + +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) { + : mTruncationFraction(truncationFractionVal), + mMaxLayers(maxLayers), + mEtaBins(etaBins), + mEtaMin(etaMin), + mEtaMax(etaMax), + mPtBins(ptBins), + mPtMin(ptMin), + mPtMax(ptMax), + mEtaBinWidth((etaMax - etaMin) / etaBins), + mPtBinWidth((ptMax - ptMin) / ptBins) + { mPdgToIndexMap.reserve(10); mIndexToPdgMap.reserve(10); } ToTLUT() = delete; - ~ToTLUT() { - for (auto hist_ptr : mLUTHistogramFlat) { + ~ToTLUT() + { + for (const auto& hist_ptr : mLUTHistogramFlat) { if (hist_ptr) { delete hist_ptr; } } } - bool load(int pdg, const std::string& filename) { + bool load(int pdg, const std::string& filename) + { TFile* f = TFile::Open(filename.c_str()); if (!f || f->IsZombie()) { - std::cerr << "ERROR: Failed to open LUT file: " << filename << std::endl; + LOG(error) << "Failed to open LUT file: " << filename; return false; } - int current_pdg_idx; + int currentPdgIdx; auto it = mPdgToIndexMap.find(pdg); if (it == mPdgToIndexMap.end()) { - current_pdg_idx = mIndexToPdgMap.size(); - mPdgToIndexMap[pdg] = current_pdg_idx; - mIndexToPdgMap.push_back(pdg); - - size_t totalSize = (current_pdg_idx + 1) * mMaxLayers * mEtaBins * mPtBins; - if (mLUTHistogramFlat.size() < totalSize) { - mLUTHistogramFlat.resize(totalSize, nullptr); - } + currentPdgIdx = mIndexToPdgMap.size(); + mPdgToIndexMap[pdg] = currentPdgIdx; + mIndexToPdgMap.push_back(pdg); + + size_t totalSize = (currentPdgIdx + 1) * mMaxLayers * mEtaBins * mPtBins; + if (mLUTHistogramFlat.size() < totalSize) { + mLUTHistogramFlat.resize(totalSize, nullptr); + } } else { - current_pdg_idx = it->second; + currentPdgIdx = it->second; } bool success = true; for (int layer = 0; layer < mMaxLayers; ++layer) { for (int etaBin = 0; etaBin < mEtaBins; ++etaBin) { - float etaMin_bin = mEtaMin + etaBin * mEtaBinWidth; - float etaMax_bin = etaMin_bin + mEtaBinWidth; + float etaMinBin = mEtaMin + etaBin * mEtaBinWidth; + float etaMaxBin = etaMinBin + mEtaBinWidth; for (int ptBin = 0; ptBin < mPtBins; ++ptBin) { - float ptMin_bin = mPtMin + ptBin * mPtBinWidth; - float ptMax_bin = ptMin_bin + mPtBinWidth; - TString histName = Form("tot_%d_barrel%d_eta%.2f-%.2f_pt%.2f-%.2f", pdg, layer, etaMin_bin, etaMax_bin, ptMin_bin, ptMax_bin); - - TH1F* hist_from_file = dynamic_cast(f->Get(histName)); - if (hist_from_file) { - TH1F* cloned_hist = static_cast(hist_from_file->Clone()); - cloned_hist->SetDirectory(nullptr); - - size_t flatIdx = getFlatIndex(current_pdg_idx, layer, etaBin, ptBin); - mLUTHistogramFlat[flatIdx] = cloned_hist; + float ptMinBin = mPtMin + ptBin * mPtBinWidth; + float ptMaxBin = ptMinBin + mPtBinWidth; + 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)); + if (histFromFile) { + TH1F* clonedHist = static_cast(histFromFile->Clone()); + clonedHist->SetDirectory(nullptr); + + size_t flatIdx = getFlatIndex(currentPdgIdx, layer, etaBin, ptBin); + mLUTHistogramFlat[flatIdx] = clonedHist; } else { - size_t flatIdx = getFlatIndex(current_pdg_idx, layer, etaBin, ptBin); + size_t flatIdx = getFlatIndex(currentPdgIdx, layer, etaBin, ptBin); mLUTHistogramFlat[flatIdx] = nullptr; success = false; } @@ -145,47 +172,53 @@ class ToTLUT { return success; } - TH1F* getHistogramForSampling(int pdg_idx, int layer, int etaBin, int ptBin) const { - if (pdg_idx < 0 || static_cast(pdg_idx) >= getNumPdgTypes() || + TH1F* getHistogramForSampling(int pdgIdx, int layer, int etaBin, int ptBin) const + { + if (pdgIdx < 0 || static_cast(pdgIdx) >= getNumPdgTypes() || layer < 0 || layer >= mMaxLayers || etaBin < 0 || etaBin >= mEtaBins || ptBin < 0 || ptBin >= mPtBins) { - return nullptr; + return nullptr; } - size_t flatIdx = getFlatIndex(pdg_idx, layer, etaBin, ptBin); + size_t flatIdx = getFlatIndex(pdgIdx, layer, etaBin, ptBin); return (flatIdx < mLUTHistogramFlat.size()) ? mLUTHistogramFlat[flatIdx] : nullptr; } - int getPdgIndex(int pdgCode) const { - auto it = mPdgToIndexMap.find(pdgCode); - if (it != mPdgToIndexMap.end()) { - return it->second; - } - return -1; + int getPdgIndex(int pdgCode) const + { + auto it = mPdgToIndexMap.find(pdgCode); + if (it != mPdgToIndexMap.end()) { + return it->second; + } + return -1; } - inline int getEtaBin(float eta) const { + inline int getEtaBin(float eta) const + { const float clampedEta = std::max(mEtaMin, std::min(eta, mEtaMax - 1e-6f)); return std::min(static_cast((clampedEta - mEtaMin) / mEtaBinWidth), mEtaBins - 1); } - inline int getPtBin(float pt) const { + inline int getPtBin(float pt) const + { const float clampedPt = std::max(mPtMin, std::min(pt, mPtMax - 1e-6f)); return std::min(static_cast((clampedPt - mPtMin) / mPtBinWidth), mPtBins - 1); } - inline size_t getFlatIndex(int pdgIdx, int layer, int etaBin, int ptBin) const { + inline size_t getFlatIndex(int pdgIdx, int layer, int etaBin, int ptBin) const + { return ((pdgIdx * mMaxLayers + layer) * mEtaBins + etaBin) * mPtBins + ptBin; } - size_t getNumPdgTypes() const { + size_t getNumPdgTypes() const + { return mIndexToPdgMap.size(); } std::vector mLUTHistogramFlat; - + std::unordered_map mPdgToIndexMap; std::vector mIndexToPdgMap; - + float mTruncationFraction; int mMaxLayers; int mEtaBins; @@ -194,7 +227,7 @@ class ToTLUT { int mPtBins; float mPtMin; float mPtMax; - + float mEtaBinWidth; float mPtBinWidth; }; @@ -202,20 +235,24 @@ class ToTLUT { static constexpr int kNumHypothesisParticles = 9; std::array, kNumHypothesisParticles>, kNumHypothesisParticles> h2dBarrelNsigmaTrue; std::array, kNumHypothesisParticles> h2dHitsPerTrackVsP; -std::array, kNumHypothesisParticles> h2dToTvsP_PerParticle; +std::array, kNumHypothesisParticles> h2dToTvsPperParticle; struct OnTheFlyTrackerPid { - float calculateNsigma(float measuredToT, float expectedToT, float resolution) { - if (resolution <= 0) return 999.f; + float calculateNsigma(float measuredToT, float expectedToT, float resolution) + { + if (resolution <= 0) + return 999.f; return (measuredToT - expectedToT) / resolution; } - float getToTMeanFromMomentumSlice(std::shared_ptr hist, float momentum) { - if (!hist) return -1.f; + float getToTMeanFromMomentumSlice(std::shared_ptr hist, float momentum) + { + if (!hist) + return -1.f; int binX = hist->GetXaxis()->FindBin(momentum); TH1D* proj = hist->ProjectionY("temp", binX, binX); - if (proj->GetEntries() < 10) { + if (proj->GetEntries() < kMinEntriesForProjection) { delete proj; return -1.f; } @@ -224,11 +261,13 @@ struct OnTheFlyTrackerPid { return mean; } - float getToTResolutionFromMomentumSlice(std::shared_ptr hist, float momentum) { - if (!hist) return -1.f; + float getToTResolutionFromMomentumSlice(std::shared_ptr hist, float momentum) + { + if (!hist) + return -1.f; int binX = hist->GetXaxis()->FindBin(momentum); TH1D* proj = hist->ProjectionY("temp", binX, binX); - if (proj->GetEntries() < 10) { + if (proj->GetEntries() < kMinEntriesForProjection) { delete proj; return -1.f; } @@ -306,39 +345,40 @@ struct OnTheFlyTrackerPid { 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 mMaxBarrelLayers{"maxBarrelLayers", 11, "Maximum number of barrel layers"}; - Configurable mEtaBins{"etaBins", 50, "Number of eta bins for LUTs and histograms"}; - Configurable mEtaMin{"etaMin", -2.5f, "Minimum eta value"}; - Configurable mEtaMax{"etaMax", 2.5f, "Maximum eta value"}; - Configurable mPtBins{"ptBins", 500, "Number of pT bins for LUTs (LUTs are pT-based)"}; - Configurable mPtMin{"ptMin", 0.0f, "Minimum pT value for LUT binning"}; - Configurable mPtMax{"ptMax", 10.0f, "Maximum pT value for LUT binning"}; - Configurable mNumLogBins{"numLogBins", 200, "Number of logarithmic momentum bins"}; + 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"}; std::vector mLogBins; std::array mHypothesisPdgCodes = { - 11, // Electron - 13, // Muon - 211, // Pion - 321, // Kaon - 2212, // Proton - 1000010020, // Deuteron - 1000010030, // Triton - 1000020030, // Helium-3 - 1000020040 // Alpha + 11, // Electron + 13, // Muon + 211, // Pion + 321, // Kaon + 2212, // Proton + 1000010020, // Deuteron + 1000010030, // Triton + 1000020030, // Helium-3 + 1000020040 // Alpha }; - void init(o2::framework::InitContext&) { - if (static_cast(mMaxBarrelLayers.value) > kTrackerRadii.size()) { - LOG(fatal) << "Configured maxBarrelLayers (" << mMaxBarrelLayers.value - << ") exceeds the size of kTrackerRadii (" << kTrackerRadii.size() - << "). Please adjust maxBarrelLayers."; + void init(o2::framework::InitContext&) + { + if (static_cast(maxBarrelLayers.value) > kTrackerRadii.size()) { + LOG(fatal) << "Configured maxBarrelLayers (" << maxBarrelLayers.value + << ") exceeds the size of kTrackerRadii (" << kTrackerRadii.size() + << "). Please adjust maxBarrelLayers."; } - mToTLUT = std::make_unique(truncationFraction.value, mMaxBarrelLayers.value, - mEtaBins.value, mEtaMin.value, mEtaMax.value, - mPtBins.value, mPtMin.value, mPtMax.value); + mToTLUT = std::make_unique(truncationFraction.value, maxBarrelLayers.value, + etaBins.value, etaMin.value, etaMax.value, + ptBins.value, ptMin.value, ptMax.value); bool loaded = true; loaded &= mToTLUT->load(11, lutTotEl.value); @@ -352,25 +392,25 @@ struct OnTheFlyTrackerPid { loaded &= mToTLUT->load(1000020040, lutTotAl.value); if (!loaded) { - std::cerr << "WARNING: Failed to load one or more ToT LUTs. PID results might be incomplete." << std::endl; + LOG(warning) << "Failed to load one or more ToT LUTs. PID results might be incomplete."; } // Logarithmic momentum bins mLogBins.clear(); - double pMin = 0.05; - double pMax = 10; + double pMin = kMinMomentumForLogBins; + double pMax = kMaxMomentumForLogBins; double logMin = std::log10(pMin); double logMax = std::log10(pMax); - double dLog = (logMax - logMin) / mNumLogBins.value; - for (int i = 0; i <= mNumLogBins.value; ++i) { + 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/z} (GeV/#it{c})"}; const AxisSpec axisToT{600, 0., 300., "ToT (#mus/10#mum)"}; const AxisSpec axisNsigma{200, -10., 10., "N#sigma"}; - const AxisSpec axisLayer{mMaxBarrelLayers.value, -0.5, static_cast(mMaxBarrelLayers.value) - 0.5, "Layer"}; - const AxisSpec axisHitsPerTrack{mMaxBarrelLayers.value + 1, -0.5, static_cast(mMaxBarrelLayers.value) + 0.5, "# Hits per track"}; + const AxisSpec axisLayer{maxBarrelLayers.value, -0.5, static_cast(maxBarrelLayers.value) - 0.5, "Layer"}; + const AxisSpec axisHitsPerTrack{maxBarrelLayers.value + 1, -0.5, static_cast(maxBarrelLayers.value) + 0.5, "# Hits per track"}; histos.add("hToTvsP", "ToT vs #it{p/z}; #it{p/z} (GeV/#it{c}); ToT (#mus/10#mum)", kTH2F, {axisMomentum, axisToT}); histos.add("hToTvsPt", "ToT vs #it{p}; #it{p} (GeV/#it{c}); ToT (#mus/10#mum)", kTH2F, {mLogBins, axisToT}); @@ -378,70 +418,86 @@ struct OnTheFlyTrackerPid { histos.add("hHitMultiplicity", "Hit multiplicity along the track; # Hits per track;Counts", kTH1F, {axisHitsPerTrack}); std::vector> particleInfo = { - {11, "Elec"}, {13, "Muon"}, {211, "Pion"}, {321, "Kaon"}, {2212, "Prot"}, - {1000010020, "Deut"}, {1000010030, "Trit"}, {1000020030, "He3"}, {1000020040, "Al"} - }; - - for (size_t i_true = 0; i_true < particleInfo.size(); ++i_true) { - std::string trueName = particleInfo[i_true].second; - std::string trueNamePretty = trueName; // Fallback - if (trueName == "Elec") trueNamePretty = "#it{e}"; - else if (trueName == "Muon") trueNamePretty = "#it{#mu}"; - else if (trueName == "Pion") trueNamePretty = "#it{#pi}"; - else if (trueName == "Kaon") trueNamePretty = "#it{K}"; - else if (trueName == "Prot") trueNamePretty = "#it{p}"; - else if (trueName == "Deut") trueNamePretty = "#it{d}"; - else if (trueName == "Trit") trueNamePretty = "#it{t}"; - else if (trueName == "He3") trueNamePretty = "#it{^{3}He}"; - else if (trueName == "Al") trueNamePretty = "#it{^{4}He}"; - - std::string hitsVsPName = "HitsPerTrack/hHitsPerTrackVsP_" + trueName; - std::string hitsVsPTitle = "N_hits vs #it{p/z} for " + trueNamePretty + "; #it{p/z} (GeV/#it{c}); N_hits"; - h2dHitsPerTrackVsP[i_true] = histos.add(hitsVsPName.c_str(), hitsVsPTitle.c_str(), kTH2F, {axisMomentum, axisHitsPerTrack}); - - std::string totVsPName = "ToTvsP/hToTvsP_" + trueName; - std::string totVsPTitle = "ToT vs #it{p/z} for " + trueNamePretty + "; #it{p/z} (GeV/#it{c}); ToT (#mus/10#mum)"; - h2dToTvsP_PerParticle[i_true] = histos.add(totVsPName.c_str(), totVsPTitle.c_str(), kTH2F, {axisMomentum, axisToT}); - - - for (size_t i_hyp = 0; i_hyp < particleInfo.size(); ++i_hyp) { - std::string hypName = particleInfo[i_hyp].second; - std::string hypNamePretty = hypName; // Fallback - if (hypName == "Elec") hypNamePretty = "#it{e}"; - else if (hypName == "Muon") hypNamePretty = "#it{#mu}"; - else if (hypName == "Pion") hypNamePretty = "#it{#pi}"; - else if (hypName == "Kaon") hypNamePretty = "#it{K}"; - else if (hypName == "Prot") hypNamePretty = "#it{p}"; - else if (hypName == "Deut") hypNamePretty = "#it{d}"; - else if (hypName == "Trit") hypNamePretty = "#it{t}"; - else if (hypName == "He3") hypNamePretty = "#it{^{3}He}"; - else if (hypName == "Al") hypNamePretty = "#it{^{4}He}"; - - std::string histName = "NSigma/BarrelNsigmaTrue" + trueName + "Vs" + hypName + "Hypothesis"; - std::string histTitle = "Nsigma (True " + trueNamePretty + " vs Hyp " + hypNamePretty + "); #it{p/z} (GeV/#it{c}); N#sigma"; - h2dBarrelNsigmaTrue[i_true][i_hyp] = histos.add(histName.c_str(), histTitle.c_str(), kTH2F, {axisMomentum, axisNsigma}); - } + {11, "Elec"}, {13, "Muon"}, {211, "Pion"}, {321, "Kaon"}, {2212, "Prot"}, {1000010020, "Deut"}, {1000010030, "Trit"}, {1000020030, "He3"}, {1000020040, "Al"}}; + + for (size_t iTrue = 0; iTrue < particleInfo.size(); ++iTrue) { + std::string trueName = particleInfo[iTrue].second; + std::string trueNamePretty = trueName; // Fallback + if (trueName == "Elec") + trueNamePretty = "#it{e}"; + else if (trueName == "Muon") + trueNamePretty = "#it{#mu}"; + else if (trueName == "Pion") + trueNamePretty = "#it{#pi}"; + else if (trueName == "Kaon") + trueNamePretty = "#it{K}"; + else if (trueName == "Prot") + trueNamePretty = "#it{p}"; + else if (trueName == "Deut") + trueNamePretty = "#it{d}"; + else if (trueName == "Trit") + trueNamePretty = "#it{t}"; + else if (trueName == "He3") + trueNamePretty = "#it{^{3}He}"; + else if (trueName == "Al") + trueNamePretty = "#it{^{4}He}"; + + std::string hitsVsPName = "HitsPerTrack/hHitsPerTrackVsP_" + trueName; + std::string hitsVsPTitle = "N_hits vs #it{p/z} for " + trueNamePretty + "; #it{p/z} (GeV/#it{c}); N_hits"; + h2dHitsPerTrackVsP[iTrue] = histos.add(hitsVsPName.c_str(), hitsVsPTitle.c_str(), kTH2F, {axisMomentum, axisHitsPerTrack}); + + std::string totVsPName = "ToTvsP/hToTvsP_" + trueName; + std::string totVsPTitle = "ToT vs #it{p/z} for " + trueNamePretty + "; #it{p/z} (GeV/#it{c}); ToT (#mus/10#mum)"; + h2dToTvsPperParticle[iTrue] = histos.add(totVsPName.c_str(), totVsPTitle.c_str(), kTH2F, {axisMomentum, axisToT}); + + for (size_t iHyp = 0; iHyp < particleInfo.size(); ++iHyp) { + std::string hypName = particleInfo[iHyp].second; + std::string hypNamePretty = hypName; // Fallback + if (hypName == "Elec") + hypNamePretty = "#it{e}"; + else if (hypName == "Muon") + hypNamePretty = "#it{#mu}"; + else if (hypName == "Pion") + hypNamePretty = "#it{#pi}"; + else if (hypName == "Kaon") + hypNamePretty = "#it{K}"; + else if (hypName == "Prot") + hypNamePretty = "#it{p}"; + else if (hypName == "Deut") + hypNamePretty = "#it{d}"; + else if (hypName == "Trit") + hypNamePretty = "#it{t}"; + else if (hypName == "He3") + hypNamePretty = "#it{^{3}He}"; + else if (hypName == "Al") + hypNamePretty = "#it{^{4}He}"; + + std::string histName = "NSigma/BarrelNsigmaTrue" + trueName + "Vs" + hypName + "Hypothesis"; + std::string histTitle = "Nsigma (True " + trueNamePretty + " vs Hyp " + hypNamePretty + "); #it{p/z} (GeV/#it{c}); N#sigma"; + h2dBarrelNsigmaTrue[iTrue][iHyp] = histos.add(histName.c_str(), histTitle.c_str(), kTH2F, {axisMomentum, axisNsigma}); + } } } void process(soa::Join::iterator const& collision, soa::Join const& tracks, aod::McParticles const& /*mcParticles*/, - aod::McCollisions const& /*mcCollisions*/) { + aod::McCollisions const& /*mcCollisions*/) + { o2::dataformats::VertexBase mcPvVtx({0.0f, 0.0f, 0.0f}, {0.}); if (collision.has_mcCollision()) { - const auto& mcCollisionObject = collision.mcCollision(); - mcPvVtx.setX(mcCollisionObject.posX()); - mcPvVtx.setY(mcCollisionObject.posY()); - mcPvVtx.setZ(mcCollisionObject.posZ()); + const auto& mcCollisionObject = collision.mcCollision(); + mcPvVtx.setX(mcCollisionObject.posX()); + mcPvVtx.setY(mcCollisionObject.posY()); + mcPvVtx.setZ(mcCollisionObject.posZ()); } int nTracksProcessed = 0; int nValidToT = 0; - for (const auto& track : tracks) { + for (const auto& track : tracks) { nTracksProcessed++; float truncatedMeanToT = -1.0f; std::array nSigmaValues; @@ -450,17 +506,17 @@ struct OnTheFlyTrackerPid { if (!track.has_mcParticle()) { tableUpgradeTrkPidSignals(truncatedMeanToT); tableUpgradeTrkPids(nSigmaValues[0], nSigmaValues[1], nSigmaValues[2], nSigmaValues[3], - nSigmaValues[4], nSigmaValues[5], nSigmaValues[6], nSigmaValues[7], nSigmaValues[8]); + nSigmaValues[4], nSigmaValues[5], nSigmaValues[6], nSigmaValues[7], nSigmaValues[8]); continue; } const auto& mcParticle = track.mcParticle(); - + const auto& pdgInfo = pdg->GetParticle(mcParticle.pdgCode()); if (!pdgInfo) { tableUpgradeTrkPidSignals(truncatedMeanToT); tableUpgradeTrkPids(nSigmaValues[0], nSigmaValues[1], nSigmaValues[2], nSigmaValues[3], - nSigmaValues[4], nSigmaValues[5], nSigmaValues[6], nSigmaValues[7], nSigmaValues[8]); + nSigmaValues[4], nSigmaValues[5], nSigmaValues[6], nSigmaValues[7], nSigmaValues[8]); continue; } @@ -477,12 +533,12 @@ struct OnTheFlyTrackerPid { } } - int true_pdg_idx = mToTLUT->getPdgIndex(truePdgCode); - if (true_pdg_idx == -1) { - tableUpgradeTrkPidSignals(truncatedMeanToT); - tableUpgradeTrkPids(nSigmaValues[0], nSigmaValues[1], nSigmaValues[2], nSigmaValues[3], - nSigmaValues[4], nSigmaValues[5], nSigmaValues[6], nSigmaValues[7], nSigmaValues[8]); - continue; + int truePdgIdx = mToTLUT->getPdgIndex(truePdgCode); + if (truePdgIdx == -1) { + tableUpgradeTrkPidSignals(truncatedMeanToT); + tableUpgradeTrkPids(nSigmaValues[0], nSigmaValues[1], nSigmaValues[2], nSigmaValues[3], + nSigmaValues[4], nSigmaValues[5], nSigmaValues[6], nSigmaValues[7], nSigmaValues[8]); + continue; } const int binnedPt = mToTLUT->getPtBin(pt); @@ -491,7 +547,7 @@ struct OnTheFlyTrackerPid { uint16_t hitMap = 0; int nHitLayers = 0; o2::track::TrackParCov o2track = o2::upgrade::convertMCParticleToO2Track(mcParticle, pdg); - + float xPv = -100.f; static constexpr float kTrkXThreshold = -99.f; if (o2track.propagateToDCA(mcPvVtx, dBz)) { @@ -499,10 +555,10 @@ struct OnTheFlyTrackerPid { } if (xPv > kTrkXThreshold) { - for (int layer = 0; layer < mMaxBarrelLayers.value; ++layer) { + for (int layer = 0; layer < maxBarrelLayers.value; ++layer) { float layerRadius = kTrackerRadii[layer]; float trackLength = computeTrackLength(o2track, layerRadius, dBz); - + if (trackLength > 0) { hitMap |= (1 << layer); histos.fill(HIST("hHitLayers"), layer); @@ -512,14 +568,14 @@ struct OnTheFlyTrackerPid { } histos.fill(HIST("hHitMultiplicity"), nHitLayers); - h2dHitsPerTrackVsP[true_pdg_idx]->Fill(rigidity, nHitLayers); + h2dHitsPerTrackVsP[truePdgIdx]->Fill(rigidity, nHitLayers); std::vector validToTs; - - for (int layer = 3; layer < mMaxBarrelLayers.value; ++layer) { + + for (int layer = kMinLayerForTruncation; layer < maxBarrelLayers.value; ++layer) { if ((hitMap >> layer) & 0x1) { - TH1F* totHist = mToTLUT->getHistogramForSampling(true_pdg_idx, layer, binnedEta, binnedPt); - + TH1F* totHist = mToTLUT->getHistogramForSampling(truePdgIdx, layer, binnedEta, binnedPt); + if (totHist && totHist->GetEntries() > 1) { float sampledToT = totHist->GetRandom(); validToTs.push_back(sampledToT); @@ -531,10 +587,14 @@ struct OnTheFlyTrackerPid { const size_t nValid = validToTs.size(); size_t nUse = 0; - if (nValid == 3 || nValid == 4) nUse = 2; - else if (nValid == 1 || nValid == 2) nUse = 1; - else if (nValid == 5 || nValid == 6) nUse = 3; - else if (nValid >= 7) nUse = 4; + if (nValid == kMidLowValidHits || nValid == kMidHighValidHits) + nUse = kMaxValidHitsForTruncation34; + else if (nValid == kMinValidHits || nValid == kLowValidHits) + nUse = kMaxValidHitsForTruncation12; + else if (nValid == kHighValidHits1 || nValid == kHighValidHits2) + nUse = kMaxValidHitsForTruncation56; + else if (nValid >= kMaxValidHits) + nUse = kMaxValidHitsForTruncation7Plus; if (nUse > 0 && nValid >= nUse) { nValidToT++; @@ -547,40 +607,41 @@ struct OnTheFlyTrackerPid { histos.fill(HIST("hToTvsPt"), p, truncatedMeanToT); histos.fill(HIST("hToTvsP"), rigidity, truncatedMeanToT); - h2dToTvsP_PerParticle[true_pdg_idx]->Fill(rigidity, truncatedMeanToT); + h2dToTvsPperParticle[truePdgIdx]->Fill(rigidity, truncatedMeanToT); } nSigmaValues.fill(999.f); if (truncatedMeanToT > 0) { - for (size_t i_hyp = 0; i_hyp < mHypothesisPdgCodes.size(); ++i_hyp) { - int hypPdgCode = mHypothesisPdgCodes[i_hyp]; - int hyp_pdg_idx = mToTLUT->getPdgIndex(hypPdgCode); + for (size_t iHyp = 0; iHyp < mHypothesisPdgCodes.size(); ++iHyp) { + int hypPdgCode = mHypothesisPdgCodes[iHyp]; + int hypPdgIdx = mToTLUT->getPdgIndex(hypPdgCode); - if (hyp_pdg_idx == -1) { - nSigmaValues[i_hyp] = 999.f; - continue; + if (hypPdgIdx == -1) { + nSigmaValues[iHyp] = 999.f; + continue; } - float expectedToT = getToTMeanFromMomentumSlice(h2dToTvsP_PerParticle[hyp_pdg_idx], rigidity); - float resolution = getToTResolutionFromMomentumSlice(h2dToTvsP_PerParticle[hyp_pdg_idx], rigidity); + float expectedToT = getToTMeanFromMomentumSlice(h2dToTvsPperParticle[hypPdgIdx], rigidity); + float resolution = getToTResolutionFromMomentumSlice(h2dToTvsPperParticle[hypPdgIdx], rigidity); if (expectedToT > 0 && resolution > 0) { - nSigmaValues[i_hyp] = calculateNsigma(truncatedMeanToT, expectedToT, resolution); - h2dBarrelNsigmaTrue[true_pdg_idx][i_hyp]->Fill(rigidity, nSigmaValues[i_hyp]); + nSigmaValues[iHyp] = calculateNsigma(truncatedMeanToT, expectedToT, resolution); + h2dBarrelNsigmaTrue[truePdgIdx][iHyp]->Fill(rigidity, nSigmaValues[iHyp]); } else { - nSigmaValues[i_hyp] = 999.f; + nSigmaValues[iHyp] = 999.f; } } } tableUpgradeTrkPidSignals(truncatedMeanToT); tableUpgradeTrkPids(nSigmaValues[0], nSigmaValues[1], nSigmaValues[2], nSigmaValues[3], - nSigmaValues[4], nSigmaValues[5], nSigmaValues[6], nSigmaValues[7], nSigmaValues[8]); + nSigmaValues[4], nSigmaValues[5], nSigmaValues[6], nSigmaValues[7], nSigmaValues[8]); } } }; -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ return WorkflowSpec{adaptAnalysisTask(cfgc)}; -} \ No newline at end of file +} From 99a9efa91c7d14df9273038be66e2235b180d7be Mon Sep 17 00:00:00 2001 From: Henrik Fribert Date: Wed, 24 Sep 2025 12:41:46 +0200 Subject: [PATCH 04/23] Fix: Remove Alpha LUT from onTheFlyTracker.cxx --- ALICE3/TableProducer/OTF/onTheFlyTracker.cxx | 1 - 1 file changed, 1 deletion(-) diff --git a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx index 461456409e7..74b7e05247f 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx @@ -266,7 +266,6 @@ struct OnTheFlyTracker { loadLUT(1000010020, lutDe.value); loadLUT(1000010030, lutTr.value); loadLUT(1000020030, lutHe3.value); - loadLUT(1000020040, lutHe3.value); // interpolate efficiencies if requested to do so mSmearer.interpolateEfficiency(static_cast(interpolateLutEfficiencyVsNch)); From e77349634674ab6bc3a41ea281936b75c9c77371 Mon Sep 17 00:00:00 2001 From: Henrik Fribert Date: Wed, 24 Sep 2025 13:06:57 +0200 Subject: [PATCH 05/23] Fix: Missing include --- ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx | 1 + 1 file changed, 1 insertion(+) diff --git a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx index 14e05efd242..01f65c8583d 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx @@ -62,6 +62,7 @@ #include #include #include +#include using namespace o2; using namespace o2::framework; From 42701cea5f79df9af6711611315e804898daeb82 Mon Sep 17 00:00:00 2001 From: Henrik Fribert Date: Wed, 24 Sep 2025 13:10:15 +0200 Subject: [PATCH 06/23] Fix: Formatting --- ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx index 01f65c8583d..2eb3716f6a8 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx @@ -58,11 +58,11 @@ #include #include #include +#include #include #include #include #include -#include using namespace o2; using namespace o2::framework; From 69175cd461a30d5dc69a08eff8acc90a076cbfdb Mon Sep 17 00:00:00 2001 From: Henrik Fribert Date: Wed, 24 Sep 2025 14:22:24 +0200 Subject: [PATCH 07/23] Fix: Removed unused variables --- ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx | 5 ----- 1 file changed, 5 deletions(-) diff --git a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx index 2eb3716f6a8..99c6a31733e 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx @@ -495,11 +495,7 @@ struct OnTheFlyTrackerPid { mcPvVtx.setZ(mcCollisionObject.posZ()); } - int nTracksProcessed = 0; - int nValidToT = 0; - for (const auto& track : tracks) { - nTracksProcessed++; float truncatedMeanToT = -1.0f; std::array nSigmaValues; nSigmaValues.fill(999.f); @@ -598,7 +594,6 @@ struct OnTheFlyTrackerPid { nUse = kMaxValidHitsForTruncation7Plus; if (nUse > 0 && nValid >= nUse) { - nValidToT++; std::sort(validToTs.begin(), validToTs.end()); float sum = 0.0f; for (size_t i = 0; i < nUse; ++i) { From 0ecbf20d5007b0443f3db112acedb835c485b553 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Wed, 24 Sep 2025 15:48:56 +0200 Subject: [PATCH 08/23] Refactor long-lived particle handling in tracker --- ALICE3/TableProducer/OTF/onTheFlyTracker.cxx | 20 ++++++++++++++++---- 1 file changed, 16 insertions(+), 4 deletions(-) diff --git a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx index 74b7e05247f..57d10fddac3 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx @@ -113,7 +113,7 @@ struct OnTheFlyTracker { Configurable lutDe{"lutDe", "lutCovm.de.dat", "LUT for deuterons"}; Configurable lutTr{"lutTr", "lutCovm.tr.dat", "LUT for tritons"}; Configurable lutHe3{"lutHe3", "lutCovm.he3.dat", "LUT for Helium-3"}; - Configurable lutAl{"lutAl", "lutCovm.he3.dat", "LUT for Alphas"}; // To be created, for now propagating as He-3 + Configurable lutAl{"lutAl", "lutCovm.he3.dat", "LUT for Alphas"}; struct : ConfigurableGroup { ConfigurableAxis axisMomentum{"axisMomentum", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f, 1.7f, 1.8f, 1.9f, 2.0f, 2.2f, 2.4f, 2.6f, 2.8f, 3.0f, 3.2f, 3.4f, 3.6f, 3.8f, 4.0f, 4.4f, 4.8f, 5.2f, 5.6f, 6.0f, 6.5f, 7.0f, 7.5f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 17.0f, 19.0f, 21.0f, 23.0f, 25.0f, 30.0f, 35.0f, 40.0f, 50.0f}, "#it{p} (GeV/#it{c})"}; @@ -507,6 +507,16 @@ struct OnTheFlyTracker { auto ir = irSampler.generateCollisionTime(); const float eventCollisionTime = ir.timeInBCNS; + constexpr std::array longLivedHandledPDGs = {kElectron, + kMuonMinus, + kPiPlus, + kKPlus, + kProton, + o2::constants::physics::kDeuteron, + o2::constants::physics::kTriton, + o2::constants::physics::kHelium3, + o2::constants::physics::kAlpha}; + // First we compute the number of charged particles in the event dNdEta = 0.f; for (const auto& mcParticle : mcParticles) { @@ -517,7 +527,8 @@ struct OnTheFlyTracker { continue; } const auto pdg = std::abs(mcParticle.pdgCode()); - if (pdg != kElectron && pdg != kMuonMinus && pdg != kPiPlus && pdg != kKPlus && pdg != kProton && pdg != 1000010020 && pdg != 1000010030 && pdg != 1000020030 && pdg != 1000020040) { + const bool longLivedToBeHandled = std::find(longLivedHandledPDGs.begin(), longLivedHandledPDGs.end(), pdg) != longLivedHandledPDGs.end(); + if (!longLivedToBeHandled) { if (!cascadeDecaySettings.decayXi) { continue; } else if (pdg != 3312) { @@ -564,8 +575,9 @@ struct OnTheFlyTracker { continue; } } - if (pdg != kElectron && pdg != kMuonMinus && pdg != kPiPlus && pdg != kKPlus && pdg != kProton && pdg != 1000010020 && pdg != 1000010030 && pdg != 1000020030 && pdg != 1000020040) { - if (!cascadeDecaySettings.decayXi) { + const bool longLivedToBeHandled = std::find(longLivedHandledPDGs.begin(), longLivedHandledPDGs.end(), pdg) != longLivedHandledPDGs.end(); + if (!longLivedToBeHandled) { + if (!(cascadeDecaySettings.decayXi)) { continue; } else if (pdg != 3312) { continue; From 4c4d522f9d02998e14d60c13e9a50edd578e3e01 Mon Sep 17 00:00:00 2001 From: Henrik Fribert Date: Wed, 24 Sep 2025 19:25:23 +0200 Subject: [PATCH 09/23] Fix: Missing alpha propagation --- ALICE3/Core/DelphesO2TrackSmearer.cxx | 13 +++++++++---- ALICE3/Core/DelphesO2TrackSmearer.h | 6 +++++- ALICE3/TableProducer/OTF/onTheFlyTracker.cxx | 1 + 3 files changed, 15 insertions(+), 5 deletions(-) diff --git a/ALICE3/Core/DelphesO2TrackSmearer.cxx b/ALICE3/Core/DelphesO2TrackSmearer.cxx index af72d0cf8b0..fd06a0c9251 100644 --- a/ALICE3/Core/DelphesO2TrackSmearer.cxx +++ b/ALICE3/Core/DelphesO2TrackSmearer.cxx @@ -106,10 +106,15 @@ bool TrackSmearer::loadTable(int pdg, const char* filename, bool forceReload) return false; } if (mLUTHeader[ipdg]->pdg != pdg) { - LOG(info) << " --- LUT header PDG mismatch: expected/detected = " << pdg << "/" << mLUTHeader[ipdg]->pdg << std::endl; - delete mLUTHeader[ipdg]; - mLUTHeader[ipdg] = nullptr; - return false; + // Special case: Allow Alpha particles to use He3 LUT + if (pdg == o2::constants::physics::kAlpha && mLUTHeader[ipdg]->pdg == o2::constants::physics::kHelium3) { + LOG(info) << " --- Alpha particles (PDG " << pdg << ") will use He3 LUT data (PDG " << mLUTHeader[ipdg]->pdg << ")" << std::endl; + } else { + LOG(info) << " --- LUT header PDG mismatch: expected/detected = " << pdg << "/" << mLUTHeader[ipdg]->pdg << std::endl; + delete mLUTHeader[ipdg]; + mLUTHeader[ipdg] = nullptr; + return false; + } } const int nnch = mLUTHeader[ipdg]->nchmap.nbins; const int nrad = mLUTHeader[ipdg]->radmap.nbins; diff --git a/ALICE3/Core/DelphesO2TrackSmearer.h b/ALICE3/Core/DelphesO2TrackSmearer.h index 2e2198bd68f..48b05d3169e 100644 --- a/ALICE3/Core/DelphesO2TrackSmearer.h +++ b/ALICE3/Core/DelphesO2TrackSmearer.h @@ -214,6 +214,8 @@ class TrackSmearer return 6; // Triton case 1000020030: return 7; // Helium3 + case 1000020040: + return 8; // Alphas default: return 2; // Default: pion } @@ -238,6 +240,8 @@ class TrackSmearer return "triton"; case 1000020030: return "helium3"; + case 1000020040: + return "alpha"; default: return "pion"; // Default: pion } @@ -246,7 +250,7 @@ class TrackSmearer void setCcdbManager(o2::ccdb::BasicCCDBManager* mgr) { mCcdbManager = mgr; } //; protected: - static constexpr unsigned int nLUTs = 8; // Number of LUT available + static constexpr unsigned int nLUTs = 9; // Number of LUT available lutHeader_t* mLUTHeader[nLUTs] = {nullptr}; lutEntry_t***** mLUTEntry[nLUTs] = {nullptr}; bool mUseEfficiency = true; diff --git a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx index 57d10fddac3..ff5f931dc9e 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx @@ -266,6 +266,7 @@ struct OnTheFlyTracker { loadLUT(1000010020, lutDe.value); loadLUT(1000010030, lutTr.value); loadLUT(1000020030, lutHe3.value); + loadLUT(1000020040, lutAl.value); // interpolate efficiencies if requested to do so mSmearer.interpolateEfficiency(static_cast(interpolateLutEfficiencyVsNch)); From 7dd97f6c8ae4e0406e06c6196ac6d904825e83a6 Mon Sep 17 00:00:00 2001 From: Henrik Fribert Date: Wed, 24 Sep 2025 19:32:41 +0200 Subject: [PATCH 10/23] Fix: MegaLinter --- ALICE3/Core/DelphesO2TrackSmearer.h | 1 + 1 file changed, 1 insertion(+) diff --git a/ALICE3/Core/DelphesO2TrackSmearer.h b/ALICE3/Core/DelphesO2TrackSmearer.h index 48b05d3169e..770b51e75a6 100644 --- a/ALICE3/Core/DelphesO2TrackSmearer.h +++ b/ALICE3/Core/DelphesO2TrackSmearer.h @@ -29,6 +29,7 @@ #include +#include #include #include #include From 1eca32249de4533d02fe9292d7e90f5e3c3db1fa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Thu, 25 Sep 2025 05:26:53 +0200 Subject: [PATCH 11/23] Refactor PDG mismatch handling for special cases --- ALICE3/Core/DelphesO2TrackSmearer.cxx | 22 +++++++++++++--------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/ALICE3/Core/DelphesO2TrackSmearer.cxx b/ALICE3/Core/DelphesO2TrackSmearer.cxx index fd06a0c9251..33f87ec4549 100644 --- a/ALICE3/Core/DelphesO2TrackSmearer.cxx +++ b/ALICE3/Core/DelphesO2TrackSmearer.cxx @@ -105,16 +105,20 @@ bool TrackSmearer::loadTable(int pdg, const char* filename, bool forceReload) mLUTHeader[ipdg] = nullptr; return false; } - if (mLUTHeader[ipdg]->pdg != pdg) { - // Special case: Allow Alpha particles to use He3 LUT - if (pdg == o2::constants::physics::kAlpha && mLUTHeader[ipdg]->pdg == o2::constants::physics::kHelium3) { + bool specialPdgCase = false; + switch(mLUTHeader[ipdg]->pdg) { // Handle special cases + case o2::constants::physics::kAlpha; // Special case: Allow Alpha particles to use He3 LUT + specialPdgCase = (mLUTHeader[ipdg]->pdg == o2::constants::physics::kHelium3); + if (specialPdgCase) LOG(info) << " --- Alpha particles (PDG " << pdg << ") will use He3 LUT data (PDG " << mLUTHeader[ipdg]->pdg << ")" << std::endl; - } else { - LOG(info) << " --- LUT header PDG mismatch: expected/detected = " << pdg << "/" << mLUTHeader[ipdg]->pdg << std::endl; - delete mLUTHeader[ipdg]; - mLUTHeader[ipdg] = nullptr; - return false; - } + break; + default: + } + if (mLUTHeader[ipdg]->pdg != pdg && !specialPdgCase) { + LOG(info) << " --- LUT header PDG mismatch: expected/detected = " << pdg << "/" << mLUTHeader[ipdg]->pdg << std::endl; + delete mLUTHeader[ipdg]; + mLUTHeader[ipdg] = nullptr; + return false; } const int nnch = mLUTHeader[ipdg]->nchmap.nbins; const int nrad = mLUTHeader[ipdg]->radmap.nbins; From 14a09675c781c2b2c0f7822ac46b2bf2206300f8 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Thu, 25 Sep 2025 05:30:19 +0200 Subject: [PATCH 12/23] Refactor condition for cascadeDecaySettings check --- ALICE3/TableProducer/OTF/onTheFlyTracker.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx index ff5f931dc9e..37212633f61 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTracker.cxx @@ -578,7 +578,7 @@ struct OnTheFlyTracker { } const bool longLivedToBeHandled = std::find(longLivedHandledPDGs.begin(), longLivedHandledPDGs.end(), pdg) != longLivedHandledPDGs.end(); if (!longLivedToBeHandled) { - if (!(cascadeDecaySettings.decayXi)) { + if (!cascadeDecaySettings.decayXi) { continue; } else if (pdg != 3312) { continue; From 2c9a8ea5c72d5c1c26394200b65caf41257b087d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Thu, 25 Sep 2025 05:45:34 +0200 Subject: [PATCH 13/23] Support CCDB --- .../TableProducer/OTF/onTheFlyTrackerPid.cxx | 70 ++++++++++++------- 1 file changed, 46 insertions(+), 24 deletions(-) diff --git a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx index 99c6a31733e..592acb628b6 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx @@ -28,29 +28,29 @@ #include "Common/Core/trackUtilities.h" #include "Common/DataModel/TrackSelectionTables.h" -#include "CCDB/BasicCCDBManager.h" -#include "CCDB/CcdbApi.h" -#include "CommonConstants/GeomConstants.h" -#include "CommonConstants/PhysicsConstants.h" -#include "CommonUtils/NameConf.h" -#include "DataFormatsCalibration/MeanVertexObject.h" -#include "DataFormatsParameters/GRPMagField.h" -#include "DetectorsBase/GeometryManager.h" -#include "DetectorsBase/Propagator.h" -#include "DetectorsVertexing/HelixHelper.h" -#include "Framework/ASoAHelpers.h" -#include "Framework/AnalysisDataModel.h" -#include "Framework/AnalysisTask.h" -#include "Framework/HistogramRegistry.h" -#include "Framework/O2DatabasePDGPlugin.h" -#include "Framework/RunningWorkflowInfo.h" -#include "Framework/runDataProcessing.h" -#include "ReconstructionDataFormats/DCA.h" - -#include "TF1.h" -#include "TH2F.h" -#include "TString.h" -#include "TVector3.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include #include #include @@ -119,8 +119,30 @@ class ToTLUT } } - bool load(int pdg, const std::string& filename) + bool load(int pdg, const std::string& filename, o2::ccdb::BasicCCDBManager *ccdb=nullptr) { + if (strncmp(filename, "ccdb:", 5) == 0) { // Check if filename starts with "ccdb:" + // NOTE: this eventually will directly fetch the object from CCDB instead of accessing it via the root file + LOG(info) << " --- ToT LUT file source identified as CCDB."; + std::string path = std::string(filename).substr(5); // Remove "ccdb:" prefix + const std::string outPath = "/tmp/ToTLUTs/"; + filename = Form("%s/%s/snapshot.root", outPath.c_str(), path.c_str()); + TFile checkFile(filename.c_str(), "READ"); + if (!checkFile.IsOpen()) { // File does not exist, retrieve from CCDB + LOG(info) << " --- CCDB source detected for PDG " << pdg << ": " << path; + if (!ccdb) { + LOG(fatal) << " --- CCDB manager not set. Please provide it."; + } + std::map metadata; + ccdb->getCCDBAccessor().retrieveBlob(path, outPath, metadata, 1); + // Add CCDB handling logic here if needed + LOG(info) << " --- Now retrieving ToT LUT file from CCDB to: " << filename; + } else { // File exists, proceed to load + LOG(info) << " --- ToT LUT file already exists: " << filename << ". Skipping download."; + checkFile.Close(); + } + return load(pdg, filename); + } TFile* f = TFile::Open(filename.c_str()); if (!f || f->IsZombie()) { LOG(error) << "Failed to open LUT file: " << filename; From eca467c412e9404d49efec868c40f37cef88776c Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Thu, 25 Sep 2025 03:46:22 +0000 Subject: [PATCH 14/23] Please consider the following formatting changes --- ALICE3/Core/DelphesO2TrackSmearer.cxx | 13 +++--- .../TableProducer/OTF/onTheFlyTrackerPid.cxx | 46 +++++++++---------- 2 files changed, 30 insertions(+), 29 deletions(-) diff --git a/ALICE3/Core/DelphesO2TrackSmearer.cxx b/ALICE3/Core/DelphesO2TrackSmearer.cxx index 33f87ec4549..b9ff70804ac 100644 --- a/ALICE3/Core/DelphesO2TrackSmearer.cxx +++ b/ALICE3/Core/DelphesO2TrackSmearer.cxx @@ -106,13 +106,14 @@ bool TrackSmearer::loadTable(int pdg, const char* filename, bool forceReload) return false; } bool specialPdgCase = false; - switch(mLUTHeader[ipdg]->pdg) { // Handle special cases + switch (mLUTHeader[ipdg]->pdg) { // Handle special cases case o2::constants::physics::kAlpha; // Special case: Allow Alpha particles to use He3 LUT - specialPdgCase = (mLUTHeader[ipdg]->pdg == o2::constants::physics::kHelium3); - if (specialPdgCase) - LOG(info) << " --- Alpha particles (PDG " << pdg << ") will use He3 LUT data (PDG " << mLUTHeader[ipdg]->pdg << ")" << std::endl; - break; - default: + specialPdgCase = (mLUTHeader[ipdg]->pdg == o2::constants::physics::kHelium3); + if (specialPdgCase) + LOG(info) + << " --- Alpha particles (PDG " << pdg << ") will use He3 LUT data (PDG " << mLUTHeader[ipdg]->pdg << ")" << std::endl; + break; + default: } if (mLUTHeader[ipdg]->pdg != pdg && !specialPdgCase) { LOG(info) << " --- LUT header PDG mismatch: expected/detected = " << pdg << "/" << mLUTHeader[ipdg]->pdg << std::endl; diff --git a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx index 592acb628b6..628e48f5cf5 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx @@ -48,11 +48,11 @@ #include #include +#include +#include #include #include #include -#include -#include #include #include @@ -119,30 +119,30 @@ class ToTLUT } } - bool load(int pdg, const std::string& filename, o2::ccdb::BasicCCDBManager *ccdb=nullptr) + bool load(int pdg, const std::string& filename, o2::ccdb::BasicCCDBManager* ccdb = nullptr) { - if (strncmp(filename, "ccdb:", 5) == 0) { // Check if filename starts with "ccdb:" - // NOTE: this eventually will directly fetch the object from CCDB instead of accessing it via the root file - LOG(info) << " --- ToT LUT file source identified as CCDB."; - std::string path = std::string(filename).substr(5); // Remove "ccdb:" prefix - const std::string outPath = "/tmp/ToTLUTs/"; - filename = Form("%s/%s/snapshot.root", outPath.c_str(), path.c_str()); - TFile checkFile(filename.c_str(), "READ"); - if (!checkFile.IsOpen()) { // File does not exist, retrieve from CCDB - LOG(info) << " --- CCDB source detected for PDG " << pdg << ": " << path; - if (!ccdb) { - LOG(fatal) << " --- CCDB manager not set. Please provide it."; + if (strncmp(filename, "ccdb:", 5) == 0) { // Check if filename starts with "ccdb:" + // NOTE: this eventually will directly fetch the object from CCDB instead of accessing it via the root file + LOG(info) << " --- ToT LUT file source identified as CCDB."; + std::string path = std::string(filename).substr(5); // Remove "ccdb:" prefix + const std::string outPath = "/tmp/ToTLUTs/"; + filename = Form("%s/%s/snapshot.root", outPath.c_str(), path.c_str()); + TFile checkFile(filename.c_str(), "READ"); + if (!checkFile.IsOpen()) { // File does not exist, retrieve from CCDB + LOG(info) << " --- CCDB source detected for PDG " << pdg << ": " << path; + if (!ccdb) { + LOG(fatal) << " --- CCDB manager not set. Please provide it."; + } + std::map metadata; + ccdb->getCCDBAccessor().retrieveBlob(path, outPath, metadata, 1); + // Add CCDB handling logic here if needed + LOG(info) << " --- Now retrieving ToT LUT file from CCDB to: " << filename; + } else { // File exists, proceed to load + LOG(info) << " --- ToT LUT file already exists: " << filename << ". Skipping download."; + checkFile.Close(); } - std::map metadata; - ccdb->getCCDBAccessor().retrieveBlob(path, outPath, metadata, 1); - // Add CCDB handling logic here if needed - LOG(info) << " --- Now retrieving ToT LUT file from CCDB to: " << filename; - } else { // File exists, proceed to load - LOG(info) << " --- ToT LUT file already exists: " << filename << ". Skipping download."; - checkFile.Close(); + return load(pdg, filename); } - return load(pdg, filename); - } TFile* f = TFile::Open(filename.c_str()); if (!f || f->IsZombie()) { LOG(error) << "Failed to open LUT file: " << filename; From 3434287f65dddc60041d03ec52c4f90e39c118f1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Nicol=C3=B2=20Jacazio?= Date: Thu, 25 Sep 2025 05:53:56 +0200 Subject: [PATCH 15/23] Add service --- .../TableProducer/OTF/onTheFlyTrackerPid.cxx | 22 +++++++++++-------- 1 file changed, 13 insertions(+), 9 deletions(-) diff --git a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx index 628e48f5cf5..6fdd1ccf7d9 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx @@ -390,9 +390,13 @@ struct OnTheFlyTrackerPid { 1000020030, // Helium-3 1000020040 // Alpha }; + Service ccdb; void init(o2::framework::InitContext&) { + ccdb->setURL("http://alice-ccdb.cern.ch"); + ccdb->setTimestamp(-1); + if (static_cast(maxBarrelLayers.value) > kTrackerRadii.size()) { LOG(fatal) << "Configured maxBarrelLayers (" << maxBarrelLayers.value << ") exceeds the size of kTrackerRadii (" << kTrackerRadii.size() @@ -404,15 +408,15 @@ struct OnTheFlyTrackerPid { ptBins.value, ptMin.value, ptMax.value); bool loaded = true; - loaded &= mToTLUT->load(11, lutTotEl.value); - loaded &= mToTLUT->load(13, lutTotMu.value); - loaded &= mToTLUT->load(211, lutTotPi.value); - loaded &= mToTLUT->load(321, lutTotKa.value); - loaded &= mToTLUT->load(2212, lutTotPr.value); - loaded &= mToTLUT->load(1000010020, lutTotDe.value); - loaded &= mToTLUT->load(1000010030, lutTotTr.value); - loaded &= mToTLUT->load(1000020030, lutTotHe.value); - loaded &= mToTLUT->load(1000020040, lutTotAl.value); + loaded &= mToTLUT->load(11, lutTotEl.value, ccdb); + loaded &= mToTLUT->load(13, lutTotMu.value, ccdb); + loaded &= mToTLUT->load(211, lutTotPi.value, ccdb); + loaded &= mToTLUT->load(321, lutTotKa.value, ccdb); + loaded &= mToTLUT->load(2212, lutTotPr.value, ccdb); + loaded &= mToTLUT->load(1000010020, lutTotDe.value, ccdb); + loaded &= mToTLUT->load(1000010030, lutTotTr.value, ccdb); + loaded &= mToTLUT->load(1000020030, lutTotHe.value, ccdb); + loaded &= mToTLUT->load(1000020040, lutTotAl.value, ccdb); if (!loaded) { LOG(warning) << "Failed to load one or more ToT LUTs. PID results might be incomplete."; From ad09ab0541702eb85c90ccd053dfd2c1b54b0369 Mon Sep 17 00:00:00 2001 From: Henrik Fribert Date: Thu, 25 Sep 2025 15:22:59 +0200 Subject: [PATCH 16/23] Fix: Ccdb handling & syntax fixes --- ALICE3/Core/DelphesO2TrackSmearer.cxx | 7 +- .../TableProducer/OTF/onTheFlyTrackerPid.cxx | 89 +++++++++++-------- 2 files changed, 56 insertions(+), 40 deletions(-) diff --git a/ALICE3/Core/DelphesO2TrackSmearer.cxx b/ALICE3/Core/DelphesO2TrackSmearer.cxx index b9ff70804ac..3ac880ac320 100644 --- a/ALICE3/Core/DelphesO2TrackSmearer.cxx +++ b/ALICE3/Core/DelphesO2TrackSmearer.cxx @@ -106,14 +106,15 @@ bool TrackSmearer::loadTable(int pdg, const char* filename, bool forceReload) return false; } bool specialPdgCase = false; - switch (mLUTHeader[ipdg]->pdg) { // Handle special cases - case o2::constants::physics::kAlpha; // Special case: Allow Alpha particles to use He3 LUT + switch (pdg) { // Handle special cases + case o2::constants::physics::kAlpha: // Special case: Allow Alpha particles to use He3 LUT specialPdgCase = (mLUTHeader[ipdg]->pdg == o2::constants::physics::kHelium3); if (specialPdgCase) LOG(info) << " --- Alpha particles (PDG " << pdg << ") will use He3 LUT data (PDG " << mLUTHeader[ipdg]->pdg << ")" << std::endl; break; - default: + default: + break; } if (mLUTHeader[ipdg]->pdg != pdg && !specialPdgCase) { LOG(info) << " --- LUT header PDG mismatch: expected/detected = " << pdg << "/" << mLUTHeader[ipdg]->pdg << std::endl; diff --git a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx index 6fdd1ccf7d9..6ab9df2443e 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx @@ -57,6 +57,7 @@ #include #include #include +#include #include #include #include @@ -119,30 +120,39 @@ class ToTLUT } } - bool load(int pdg, const std::string& filename, o2::ccdb::BasicCCDBManager* ccdb = nullptr) + void setCcdbManager(o2::ccdb::BasicCCDBManager* mgr) { mCcdbManager = mgr; } + + bool load(int pdg, const std::string& filename) { - if (strncmp(filename, "ccdb:", 5) == 0) { // Check if filename starts with "ccdb:" - // NOTE: this eventually will directly fetch the object from CCDB instead of accessing it via the root file - LOG(info) << " --- ToT LUT file source identified as CCDB."; - std::string path = std::string(filename).substr(5); // Remove "ccdb:" prefix + if (!filename.empty() && strncmp(filename.c_str(), "ccdb:", 5) == 0) { + std::string basePath = std::string(filename).substr(5); + std::string path = basePath + "/PDG_" + std::to_string(pdg); const std::string outPath = "/tmp/ToTLUTs/"; - filename = Form("%s/%s/snapshot.root", outPath.c_str(), path.c_str()); - TFile checkFile(filename.c_str(), "READ"); - if (!checkFile.IsOpen()) { // File does not exist, retrieve from CCDB - LOG(info) << " --- CCDB source detected for PDG " << pdg << ": " << path; - if (!ccdb) { - LOG(fatal) << " --- CCDB manager not set. Please provide it."; + + std::string localFilename = Form("%s/lut_tot_%d.root", outPath.c_str(), pdg); + std::ifstream checkFile(localFilename); + if (!checkFile.is_open()) { + if (!mCcdbManager) { + LOG(fatal) << "CCDB manager not set. Please set it before loading LUT from CCDB."; } std::map metadata; - ccdb->getCCDBAccessor().retrieveBlob(path, outPath, metadata, 1); - // Add CCDB handling logic here if needed - LOG(info) << " --- Now retrieving ToT LUT file from CCDB to: " << filename; - } else { // File exists, proceed to load - LOG(info) << " --- ToT LUT file already exists: " << filename << ". Skipping download."; - checkFile.Close(); + mCcdbManager->getCCDBAccessor().retrieveBlob(path, outPath, metadata, 1); + + std::string foundFile = Form("%s/%s/snapshot.root", outPath.c_str(), path.c_str()); + std::ifstream testFile(foundFile); + if (!testFile.is_open()) { + LOG(error) << "Could not find downloaded CCDB file for PDG " << pdg; + return false; + } + testFile.close(); + + return load(pdg, foundFile); + } else { + checkFile.close(); + return load(pdg, localFilename); } - return load(pdg, filename); } + TFile* f = TFile::Open(filename.c_str()); if (!f || f->IsZombie()) { LOG(error) << "Failed to open LUT file: " << filename; @@ -253,6 +263,9 @@ class ToTLUT float mEtaBinWidth; float mPtBinWidth; + + private: + o2::ccdb::BasicCCDBManager* mCcdbManager = nullptr; }; static constexpr int kNumHypothesisParticles = 9; @@ -352,19 +365,20 @@ struct OnTheFlyTrackerPid { Produces tableUpgradeTrkPids; Service pdg; + Service ccdb; std::unique_ptr mToTLUT; HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; - Configurable lutTotEl{"lutTotEl", "lut_tot_11.root", "ToT LUT for electrons"}; - Configurable lutTotMu{"lutTotMu", "lut_tot_13.root", "ToT LUT for muons"}; - Configurable lutTotPi{"lutTotPi", "lut_tot_211.root", "ToT LUT for pions"}; - Configurable lutTotKa{"lutTotKa", "lut_tot_321.root", "ToT LUT for kaons"}; - Configurable lutTotPr{"lutTotPr", "lut_tot_2212.root", "ToT LUT for protons"}; - Configurable lutTotDe{"lutTotDe", "lut_tot_1000010020.root", "ToT LUT for deuteron"}; - Configurable lutTotTr{"lutTotTr", "lut_tot_1000010030.root", "ToT LUT for triton"}; - Configurable lutTotHe{"lutTotHe", "lut_tot_1000020030.root", "ToT LUT for helium-3"}; - Configurable lutTotAl{"lutTotAl", "lut_tot_1000020040.root", "ToT LUT for alphas"}; + Configurable lutTotEl{"lutTotEl", "ccdb:Users/h/hfribert/ToT_LUTs", "ToT LUT for electrons"}; + Configurable lutTotMu{"lutTotMu", "ccdb:Users/h/hfribert/ToT_LUTs", "ToT LUT for muons"}; + Configurable lutTotPi{"lutTotPi", "ccdb:Users/h/hfribert/ToT_LUTs", "ToT LUT for pions"}; + Configurable lutTotKa{"lutTotKa", "ccdb:Users/h/hfribert/ToT_LUTs", "ToT LUT for kaons"}; + Configurable lutTotPr{"lutTotPr", "ccdb:Users/h/hfribert/ToT_LUTs", "ToT LUT for protons"}; + Configurable lutTotDe{"lutTotDe", "ccdb:Users/h/hfribert/ToT_LUTs", "ToT LUT for deuteron"}; + Configurable lutTotTr{"lutTotTr", "ccdb:Users/h/hfribert/ToT_LUTs", "ToT LUT for triton"}; + 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"}; @@ -390,7 +404,6 @@ struct OnTheFlyTrackerPid { 1000020030, // Helium-3 1000020040 // Alpha }; - Service ccdb; void init(o2::framework::InitContext&) { @@ -407,16 +420,18 @@ struct OnTheFlyTrackerPid { etaBins.value, etaMin.value, etaMax.value, ptBins.value, ptMin.value, ptMax.value); + mToTLUT->setCcdbManager(ccdb.operator->()); + bool loaded = true; - loaded &= mToTLUT->load(11, lutTotEl.value, ccdb); - loaded &= mToTLUT->load(13, lutTotMu.value, ccdb); - loaded &= mToTLUT->load(211, lutTotPi.value, ccdb); - loaded &= mToTLUT->load(321, lutTotKa.value, ccdb); - loaded &= mToTLUT->load(2212, lutTotPr.value, ccdb); - loaded &= mToTLUT->load(1000010020, lutTotDe.value, ccdb); - loaded &= mToTLUT->load(1000010030, lutTotTr.value, ccdb); - loaded &= mToTLUT->load(1000020030, lutTotHe.value, ccdb); - loaded &= mToTLUT->load(1000020040, lutTotAl.value, ccdb); + loaded &= mToTLUT->load(11, lutTotEl.value); + loaded &= mToTLUT->load(13, lutTotMu.value); + loaded &= mToTLUT->load(211, lutTotPi.value); + loaded &= mToTLUT->load(321, lutTotKa.value); + loaded &= mToTLUT->load(2212, lutTotPr.value); + loaded &= mToTLUT->load(1000010020, lutTotDe.value); + loaded &= mToTLUT->load(1000010030, lutTotTr.value); + loaded &= mToTLUT->load(1000020030, lutTotHe.value); + loaded &= mToTLUT->load(1000020040, lutTotAl.value); if (!loaded) { LOG(warning) << "Failed to load one or more ToT LUTs. PID results might be incomplete."; From d00847e15c14742d13b2fd23750f3e481d3a8651 Mon Sep 17 00:00:00 2001 From: Henrik Fribert Date: Thu, 25 Sep 2025 15:25:56 +0200 Subject: [PATCH 17/23] Fix: Formatting --- ALICE3/Core/DelphesO2TrackSmearer.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/ALICE3/Core/DelphesO2TrackSmearer.cxx b/ALICE3/Core/DelphesO2TrackSmearer.cxx index 3ac880ac320..0002160a49b 100644 --- a/ALICE3/Core/DelphesO2TrackSmearer.cxx +++ b/ALICE3/Core/DelphesO2TrackSmearer.cxx @@ -106,12 +106,12 @@ bool TrackSmearer::loadTable(int pdg, const char* filename, bool forceReload) return false; } bool specialPdgCase = false; - switch (pdg) { // Handle special cases - case o2::constants::physics::kAlpha: // Special case: Allow Alpha particles to use He3 LUT + switch (pdg) { // Handle special cases + case o2::constants::physics::kAlpha: // Special case: Allow Alpha particles to use He3 LUT specialPdgCase = (mLUTHeader[ipdg]->pdg == o2::constants::physics::kHelium3); if (specialPdgCase) LOG(info) - << " --- Alpha particles (PDG " << pdg << ") will use He3 LUT data (PDG " << mLUTHeader[ipdg]->pdg << ")" << std::endl; + << " --- Alpha particles (PDG " << pdg << ") will use He3 LUT data (PDG " << mLUTHeader[ipdg]->pdg << ")" << std::endl; break; default: break; From 39ed52aadc6c889a9a50fc4e87ef473482c69173 Mon Sep 17 00:00:00 2001 From: Henrik Fribert Date: Mon, 29 Sep 2025 10:31:47 +0200 Subject: [PATCH 18/23] Feat: Add nuclei to TOF and RICH --- ALICE3/DataModel/OTFRICH.h | 34 +++++++- ALICE3/DataModel/OTFTOF.h | 92 +++++++++++++++++--- ALICE3/TableProducer/OTF/onTheFlyRichPid.cxx | 92 +++++++++++++++++--- ALICE3/TableProducer/OTF/onTheFlyTofPid.cxx | 81 ++++++++++------- 4 files changed, 241 insertions(+), 58 deletions(-) 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()) { From 1689bad66a58a2a8a3478dd9f5468d715fc08fcb Mon Sep 17 00:00:00 2001 From: Henrik Fribert Date: Mon, 29 Sep 2025 10:33:24 +0200 Subject: [PATCH 19/23] Feat: Task for evaluating combined PID --- ALICE3/Tasks/CMakeLists.txt | 5 +- ALICE3/Tasks/alice3PidEvaluation.cxx | 375 +++++++++++++++++++++++++++ 2 files changed, 379 insertions(+), 1 deletion(-) create mode 100644 ALICE3/Tasks/alice3PidEvaluation.cxx 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..3c405856a8d --- /dev/null +++ b/ALICE3/Tasks/alice3PidEvaluation.cxx @@ -0,0 +1,375 @@ +// 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 + +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*/) + { + int totalTracks = 0; + int analyzedTracks = 0; + + 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) { + totalTracks++; + + 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(); + + analyzedTracks++; + + 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)}; +} From d72611bfa932246bb026d394060066b2931c31e4 Mon Sep 17 00:00:00 2001 From: Henrik Fribert Date: Tue, 30 Sep 2025 14:55:40 +0200 Subject: [PATCH 20/23] Fix: Missing includes --- ALICE3/Tasks/alice3PidEvaluation.cxx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/ALICE3/Tasks/alice3PidEvaluation.cxx b/ALICE3/Tasks/alice3PidEvaluation.cxx index 3c405856a8d..f8dc2b73093 100644 --- a/ALICE3/Tasks/alice3PidEvaluation.cxx +++ b/ALICE3/Tasks/alice3PidEvaluation.cxx @@ -40,6 +40,8 @@ #include "TVector3.h" #include +#include +#include #include using namespace o2; From 356d989cc92cbd062fd5debf91ca5da4ccd7e917 Mon Sep 17 00:00:00 2001 From: Henrik Fribert Date: Wed, 1 Oct 2025 14:17:01 +0200 Subject: [PATCH 21/23] Feat: Added configurables to only load range of LUT bins in trackerPID --- .../TableProducer/OTF/onTheFlyTrackerPid.cxx | 67 ++++++++++++------- 1 file changed, 44 insertions(+), 23 deletions(-) diff --git a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx index 6ab9df2443e..810946389ee 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) + 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()) { From a3d0bfd9f41c2e3b6e3479d9be201cb2db73fbc6 Mon Sep 17 00:00:00 2001 From: Henrik Fribert Date: Wed, 1 Oct 2025 14:24:25 +0200 Subject: [PATCH 22/23] Fix: MegaLinter --- ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx index 810946389ee..b9e05d0ec75 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx @@ -103,7 +103,7 @@ static constexpr float kLUTPtMax = 10.0f; class ToTLUT { public: - ToTLUT(int maxLayers, float analysisEtaMinVal = 0.0f, float analysisEtaMaxVal = 1.0f, float analysisPtMinVal = 0.0f, float analysisPtMaxVal = 10.0f) + 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), From f5e9829ca09da324137438a3541c9ae5715e55f7 Mon Sep 17 00:00:00 2001 From: Henrik Fribert Date: Wed, 1 Oct 2025 18:16:52 +0200 Subject: [PATCH 23/23] Fix: Unused variables --- ALICE3/Tasks/alice3PidEvaluation.cxx | 6 ------ 1 file changed, 6 deletions(-) diff --git a/ALICE3/Tasks/alice3PidEvaluation.cxx b/ALICE3/Tasks/alice3PidEvaluation.cxx index f8dc2b73093..2ee64ea7623 100644 --- a/ALICE3/Tasks/alice3PidEvaluation.cxx +++ b/ALICE3/Tasks/alice3PidEvaluation.cxx @@ -149,8 +149,6 @@ struct Alice3PidEvaluation { void process(soa::Join const& tracks, aod::McParticles const& /*mcParticles*/) { - int totalTracks = 0; - int analyzedTracks = 0; auto isValidNSigma = [](float nSigma) -> bool { return (nSigma < kInvalidNSigmaValue && nSigma > -kInvalidNSigmaValue); @@ -258,8 +256,6 @@ struct Alice3PidEvaluation { }; for (const auto& track : tracks) { - totalTracks++; - if (!track.has_mcParticle()) { continue; } @@ -290,8 +286,6 @@ struct Alice3PidEvaluation { trackerNSigma[kHelium3] = track.nSigmaTrkHe(); trackerNSigma[kAlpha] = track.nSigmaTrkAl(); - analyzedTracks++; - analyzeDetector(trackerHists, trackerNSigma, trueParticleIndex, p); std::array tofInnerNSigma;