From e2e70a183083a9a8affe4b3b40f21afe1965b5d5 Mon Sep 17 00:00:00 2001 From: Henrik Fribert Date: Wed, 14 Oct 2026 18:05:46 +0200 Subject: [PATCH 1/5] Changed LUT approach for ToT --- .../TableProducer/OTF/onTheFlyTrackerPid.cxx | 846 +++++------------- 1 file changed, 221 insertions(+), 625 deletions(-) diff --git a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx index bf7a1f1350e..0d65ef5f642 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx @@ -11,23 +11,18 @@ /// /// \file onTheFlyTrackerPid.cxx /// -/// \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). +/// \brief This task simulates the PID signal that can be obtained from the tracker layers as ToT (Time-over-Threshold). +/// The dE/dx of every track is sampled from a per-track LUT made from the full simulation. +/// It contemplates 9 particle types: electrons, muons, pions, kaons, protons, deuterons, tritons, +/// helium-3 and alphas/helium-4. /// /// \author Henrik Fribert TUM /// \author Nicolò Jacazio Università del Piemonte Orientale /// \since May 22, 2025 /// -#include "GeometryContainer.h" - -#include "ALICE3/Core/TrackUtilities.h" -#include "ALICE3/DataModel/OTFPIDTrk.h" - #include #include -#include #include #include #include @@ -38,732 +33,333 @@ #include #include #include -#include -#include -#include -#include #include #include #include +#include +#include #include #include #include #include #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}; - -// 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 + +/// Local path of a LUT file; "ccdb:" is downloaded once to /tmp/ToTLUTs//snapshot.root +std::string getLocalLutFile(const std::string& filename, o2::ccdb::BasicCCDBManager* ccdbManager) { - public: - explicit ToTLUT(int maxLayers, - int etaBins, float etaMin, float etaMax, - int pBins, float pMin, float pMax, - float analysisEtaMinVal = 0.0f, float analysisEtaMaxVal = 1.0f, - float analysisPMinVal = 0.0f, float analysisPMaxVal = 10.0f) - : mMaxLayers(maxLayers), - mAnalysisEtaMin(analysisEtaMinVal), - mAnalysisEtaMax(analysisEtaMaxVal), - mAnalysisPMin(analysisPMinVal), - mAnalysisPMax(analysisPMaxVal), - mEtaBins(etaBins), - mEtaMin(etaMin), - mEtaMax(etaMax), - mPBins(pBins), - mPMin(pMin), - mPMax(pMax), - mEtaBinWidth((etaMax - etaMin) / etaBins), - mPBinWidth((pMax - pMin) / pBins) - { - mPdgToIndexMap.reserve(10); - mIndexToPdgMap.reserve(10); + if (!filename.starts_with("ccdb:")) { + return filename; } - ToTLUT() = delete; - - ~ToTLUT() - { - for (const auto& hist_ptr : mLUTHistogramFlat) { - - delete hist_ptr; + const std::string ccdbPath = filename.substr(5); + const std::string outPath = "/tmp/ToTLUTs/"; + const std::string localFilename = outPath + ccdbPath + "/snapshot.root"; + if (!std::ifstream(localFilename).is_open()) { + if (!ccdbManager) { + LOG(fatal) << "CCDB manager not set. Please set it before loading LUT from CCDB."; + } + std::map metadata; + LOG(info) << "Retrieving " << localFilename << " from CCDB path: " << ccdbPath; + ccdbManager->getCCDBAccessor().retrieveBlob(ccdbPath, outPath, metadata, 1); + if (!std::ifstream(localFilename).is_open()) { + LOG(fatal) << "Could not find downloaded CCDB file " << localFilename; } } + return localFilename; +} - void setCcdbManager(o2::ccdb::BasicCCDBManager* mgr) { mCcdbManager = mgr; } - - bool load(int pdg, const std::string& filename) +/// Per-track dE/dx LUT; for every species, |eta| bin and rigidity bin the +/// distribution of the track dE/dx in e-, as a histogram. The dE/dx of a track is sampled from it, +/// interpolated linearly in log(rigidity) between the rigidity bin centres with the same random number. +class TrackDedxLut +{ + public: + template + bool load(const std::string& filename, o2::ccdb::BasicCCDBManager* ccdbManager, const std::array& pdgCodes) { - if (filename.empty()) { - LOG(warning) << "Provided filename is empty for PDG " << pdg; + std::unique_ptr f(TFile::Open(getLocalLutFile(filename, ccdbManager).c_str())); + if (!f || f->IsZombie()) { + LOG(fatal) << "Failed to open dE/dx LUT file: " << filename; return false; } - if (filename.starts_with("ccdb:")) { // Check if filename starts with "ccdb:" - std::string ccdbPath = filename.substr(5); // remove "ccdb:" prefix - const std::string outPath = "/tmp/ToTLUTs/"; - const std::string localFilename = outPath + ccdbPath + "/snapshot.root"; - std::ifstream checkFile(localFilename); - if (!checkFile.is_open()) { // File is not found, need to download it from CCDB - if (!mCcdbManager) { - LOG(fatal) << "CCDB manager not set. Please set it before loading LUT from CCDB."; - } - std::map metadata; - LOG(info) << "Retrieving " << localFilename << " from CCDB path: " << ccdbPath; - mCcdbManager->getCCDBAccessor().retrieveBlob(ccdbPath, outPath, metadata, 1); - std::ifstream testFile(localFilename); - if (!testFile.is_open()) { - LOG(fatal) << "Could not find downloaded CCDB file for PDG " << pdg << " in file " << localFilename; - return false; - } - testFile.close(); - return load(pdg, localFilename); - } // File is found, proceed to load it - checkFile.close(); - return load(pdg, localFilename); - } - // In case the file is already available locally - TFile* f = TFile::Open(filename.c_str()); - if (!f || f->IsZombie()) { - LOG(fatal) << "Failed to open LUT file: " << filename; + const auto* etaBins = f->Get("etaBins"); + const auto* rigBins = f->Get("rigBins"); + const auto* dedxBins = f->Get("dedxBins"); + if (!etaBins || !rigBins || !dedxBins) { + LOG(fatal) << "No etaBins/rigBins/dedxBins in dE/dx LUT file: " << filename; return false; } - - int currentPdgIdx = 0; - auto it = mPdgToIndexMap.find(pdg); - if (it == mPdgToIndexMap.end()) { - currentPdgIdx = mIndexToPdgMap.size(); - mPdgToIndexMap[pdg] = currentPdgIdx; - mIndexToPdgMap.push_back(pdg); - - size_t totalSize = (currentPdgIdx + 1) * mMaxLayers * mEtaBins * mPBins; - if (mLUTHistogramFlat.size() < totalSize) { - mLUTHistogramFlat.resize(totalSize, nullptr); - } - } else { - currentPdgIdx = it->second; + mEtaEdges = binEdges(etaBins); + mRigEdges = binEdges(rigBins); + mDedxEdges = binEdges(dedxBins); + mLogRigCentres.clear(); + for (size_t i = 0; i + 1 < mRigEdges.size(); ++i) { + mLogRigCentres.push_back(0.5 * (std::log(mRigEdges[i]) + std::log(mRigEdges[i + 1]))); } - - bool success = true; - for (int layer = 0; layer < mMaxLayers; ++layer) { - 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) { + const size_t nBins = (mEtaEdges.size() - 1) * (mRigEdges.size() - 1); + const size_t nDedx = mDedxEdges.size() - 1; + + bool allLoaded = true; + for (const int pdg : pdgCodes) { + std::unique_ptr> h(f->Get>(Form("h_%d", pdg))); + if (!h || h->size() != nBins * nDedx) { + LOG(warning) << "No dE/dx LUT for PDG " << pdg << " in " << filename; + allLoaded = false; + continue; + } + // cumulative distributions, only for the filled bins + auto& offset = mOffset[pdg]; + auto& cdf = mCdf[pdg]; + offset.assign(nBins, -1); + for (size_t b = 0; b < nBins; ++b) { + const float* hb = h->data() + b * nDedx; + double sum = 0.; + for (size_t k = 0; k < nDedx; ++k) { + sum += hb[k]; + } + if (sum <= 0.) { continue; } - - for (int pBin = 0; pBin < mPBins; ++pBin) { - float pMinBin = mPMin + pBin * mPBinWidth; - float pMaxBin = pMinBin + mPBinWidth; - float pCenter = (pMinBin + pMaxBin) / 2.0f; - - if (pCenter < mAnalysisPMin || pCenter >= mAnalysisPMax) { - continue; - } - - TString histName = Form("tot_%d_barrel%d_eta%.2f-%.2f_p%.2f-%.2f", pdg, layer, etaMinBin, etaMaxBin, pMinBin, pMaxBin); - - TH1F* histFromFile = dynamic_cast(f->Get(histName)); - if (histFromFile) { - TH1F* clonedHist = dynamic_cast(histFromFile->Clone()); - clonedHist->SetDirectory(nullptr); - - size_t flatIdx = getFlatIndex(currentPdgIdx, layer, etaBin, pBin); - mLUTHistogramFlat[flatIdx] = clonedHist; - } else { - size_t flatIdx = getFlatIndex(currentPdgIdx, layer, etaBin, pBin); - mLUTHistogramFlat[flatIdx] = nullptr; - if (layer >= kMinLayerForTruncation) { - success = false; - } - } + offset[b] = cdf.size(); + double cum = 0.; + for (size_t k = 0; k < nDedx; ++k) { + cum += hb[k]; + cdf.push_back(cum / sum); } } } - - f->Close(); - delete f; - return success; + LOG(info) << "Loaded dE/dx LUT " << filename << " (" << mEtaEdges.size() - 1 << " |eta| x " << mRigEdges.size() - 1 + << " rigidity x " << mDedxEdges.size() - 1 << " dE/dx bins)"; + return allLoaded; } - TH1F* getHistogramForSampling(int pdgIdx, int layer, int etaBin, int pBin) const + /// Sampled dE/dx for uniform random number u, or -1 if there is no LUT entry + float sample(int pdg, float eta, float rigidity, double u) const { - if (pdgIdx < 0 || static_cast(pdgIdx) >= getNumPdgTypes() || - layer < 0 || layer >= mMaxLayers || - etaBin < 0 || etaBin >= mEtaBins || pBin < 0 || pBin >= mPBins) { - return nullptr; + const auto it = mOffset.find(pdg); + const float absEta = std::abs(eta); + if (it == mOffset.end() || mEtaEdges.empty() || absEta < mEtaEdges.front() || absEta >= mEtaEdges.back() || + rigidity < mRigEdges.front() || rigidity >= mRigEdges.back()) { + return -1.f; } - size_t flatIdx = getFlatIndex(pdgIdx, layer, etaBin, pBin); - return (flatIdx < mLUTHistogramFlat.size()) ? mLUTHistogramFlat[flatIdx] : nullptr; - } - - int getPdgIndex(int pdgCode) const - { - auto it = mPdgToIndexMap.find(pdgCode); - if (it != mPdgToIndexMap.end()) { - return it->second; + const int nRig = mRigEdges.size() - 1; + const int iEta = std::upper_bound(mEtaEdges.begin(), mEtaEdges.end(), absEta) - mEtaEdges.begin() - 1; + const int iRig = std::upper_bound(mRigEdges.begin(), mRigEdges.end(), rigidity) - mRigEdges.begin() - 1; + const int offset = it->second[iEta * nRig + iRig]; + if (offset < 0) { + return -1.f; } - 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 getPBin(float p) const - { - const float clampedP = std::max(mPMin, std::min(p, mPMax - 1e-6f)); - return std::min(static_cast((clampedP - mPMin) / mPBinWidth), mPBins - 1); - } - - inline size_t getFlatIndex(int pdgIdx, int layer, int etaBin, int pBin) const - { - return ((pdgIdx * mMaxLayers + layer) * mEtaBins + etaBin) * mPBins + pBin; - } - - size_t getNumPdgTypes() const - { - return mIndexToPdgMap.size(); + const float* cdf = mCdf.at(pdg).data(); + float dedx = inverseCdf(cdf + offset, u); + + // neighbouring bin on the side of the track, no interpolation if it is empty + const double logRig = std::log(rigidity); + const int iNext = logRig < mLogRigCentres[iRig] ? iRig - 1 : iRig + 1; + if (iNext >= 0 && iNext < nRig && it->second[iEta * nRig + iNext] >= 0) { + const float w = (logRig - mLogRigCentres[iRig]) / (mLogRigCentres[iNext] - mLogRigCentres[iRig]); + dedx += w * (inverseCdf(cdf + it->second[iEta * nRig + iNext], u) - dedx); + } + return dedx; } - std::vector mLUTHistogramFlat; - - std::unordered_map mPdgToIndexMap; - std::vector mIndexToPdgMap; - - int mMaxLayers; - float mAnalysisEtaMin; - float mAnalysisEtaMax; - float mAnalysisPMin; - float mAnalysisPMax; - int mEtaBins; - float mEtaMin; - float mEtaMax; - int mPBins; - float mPMin; - float mPMax; - - float mEtaBinWidth; - float mPBinWidth; - private: - o2::ccdb::BasicCCDBManager* mCcdbManager = nullptr; -}; - -struct OnTheFlyTrackerPid { - static constexpr int kNumHypothesisParticles = 9; - std::array, kNumHypothesisParticles>, kNumHypothesisParticles> h2dBarrelNsigmaTrue; - std::array, kNumHypothesisParticles> h2dHitsPerTrackVsP; - std::array, kNumHypothesisParticles> h2dToTvsPperParticle; - std::array, kNumHypothesisParticles> h2dHitsPerTrackVsPLin; - std::array, kNumHypothesisParticles> h2dToTvsPperParticleLin; - - float calculateNsigma(float measuredToT, float expectedToT, float resolution) + /// dE/dx at which the cumulative distribution reaches u, uniform within the dE/dx bin + float inverseCdf(const float* cdf, double u) const { - if (resolution <= 0) { - return 999.f; - } - return (measuredToT - expectedToT) / resolution; + const int nDedx = mDedxEdges.size() - 1; + const int k = std::min(static_cast(std::upper_bound(cdf, cdf + nDedx, u) - cdf), nDedx - 1); + const double lo = k == 0 ? 0. : cdf[k - 1]; + const double frac = cdf[k] > lo ? (u - lo) / (cdf[k] - lo) : 0.5; + return mDedxEdges[k] + frac * (mDedxEdges[k + 1] - mDedxEdges[k]); } - float getToTMeanFromMomentumSlice(const std::shared_ptr& hist, float momentum) + static std::vector binEdges(const TH1* h) { - if (!hist) { - return -1.f; - } - int binX = hist->GetXaxis()->FindBin(momentum); - TH1D* proj = hist->ProjectionY("temp", binX, binX); - if (proj->GetEntries() < kMinEntriesForProjection) { - delete proj; - return -1.f; + std::vector edges; + for (int i = 1; i <= h->GetNbinsX() + 1; ++i) { + edges.push_back(h->GetXaxis()->GetBinLowEdge(i)); } - float mean = proj->GetMean(); - delete proj; - return mean; + return edges; } - float getToTResolutionFromMomentumSlice(const 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() < kMinEntriesForProjection) { - delete proj; - return -1.f; - } - float stddev = proj->GetStdDev(); - delete proj; - return stddev; - } + std::vector mEtaEdges; + std::vector mRigEdges; + std::vector mLogRigCentres; + std::vector mDedxEdges; + std::unordered_map> mOffset; // |pdg| -> start of the cdf of bin etaBin * nRig + rigBin, -1 if empty + std::unordered_map> mCdf; // |pdg| -> cumulative dE/dx distributions of the filled bins +}; - float computeTrackLength(o2::track::TrackParCov track, float radius, float magneticField) - { - o2::math_utils::CircleXYf_t trcCircle; - float sna = NAN, csa = NAN; - 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)) { - float 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 std::array point1 = {radical * ux + displace * vx, radical * uy + displace * vy}; - const std::array point2 = {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; - } - return -100.f; - } +struct OnTheFlyTrackerPid { + static constexpr int kNumSpecies = 9; + static constexpr std::array kPdgCodes = {11, 13, 211, 321, 2212, 1000010020, 1000010030, 1000020030, 1000020040}; + static constexpr std::array kNames = {"Elec", "Muon", "Pion", "Kaon", "Prot", "Deut", "Trit", "He3", "Al"}; + static constexpr std::array kPrettyNames = {"#it{e}", "#it{#mu}", "#it{#pi}", "#it{K}", "#it{p}", "#it{d}", "#it{t}", "#it{^{3}He}", "#it{^{4}He}"}; - Produces tableUpgradeTrkPidSignals; - Produces tableUpgradeTrkPids; + std::array, kNumSpecies> h2dToTvsPperParticle; + std::array, kNumSpecies> h2dToTvsPperParticleLin; + std::array, kNumSpecies> h3dToTvsPEtaPerParticle; + std::array, kNumSpecies> h2dTracks; + std::array, kNumSpecies> h2dTracksWithDedx; Service pdgDatabase{}; Service ccdb{}; - std::unique_ptr mToTLUT; + TrackDedxLut mTrackLut; + TRandom3 mRandom; HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; - Configurable lutTotEl{"lutTotEl", "ccdb:Users/h/hfribert/ToT_LUTs_ALMIRA/Scale1/PDG_11", "ToT LUT for electrons"}; - Configurable lutTotMu{"lutTotMu", "ccdb:Users/h/hfribert/ToT_LUTs_ALMIRA/Scale1/PDG_13", "ToT LUT for muons"}; - Configurable lutTotPi{"lutTotPi", "ccdb:Users/h/hfribert/ToT_LUTs_ALMIRA/Scale1/PDG_211", "ToT LUT for pions"}; - Configurable lutTotKa{"lutTotKa", "ccdb:Users/h/hfribert/ToT_LUTs_ALMIRA/Scale1/PDG_321", "ToT LUT for kaons"}; - Configurable lutTotPr{"lutTotPr", "ccdb:Users/h/hfribert/ToT_LUTs_ALMIRA/Scale1/PDG_2212", "ToT LUT for protons"}; - Configurable lutTotDe{"lutTotDe", "ccdb:Users/h/hfribert/ToT_LUTs_ALMIRA/Scale1/PDG_1000010020", "ToT LUT for deuteron"}; - Configurable lutTotTr{"lutTotTr", "ccdb:Users/h/hfribert/ToT_LUTs_ALMIRA/Scale1/PDG_1000010030", "ToT LUT for triton"}; - Configurable lutTotHe{"lutTotHe", "ccdb:Users/h/hfribert/ToT_LUTs_ALMIRA/Scale1/PDG_1000020030", "ToT LUT for helium-3"}; - Configurable lutTotAl{"lutTotAl", "ccdb:Users/h/hfribert/ToT_LUTs_ALMIRA/Scale1/PDG_1000020040", "ToT LUT for alphas"}; - - Configurable maxBarrelLayers{"maxBarrelLayers", 11, "Maximum number of barrel layers"}; - 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 analysisPMin{"analysisPMin", 0.0f, "Minimum |p| (GeV/c) for LUT loading optimization"}; - Configurable analysisPMax{"analysisPMax", 10.0f, "Maximum |p| (GeV/c) for LUT loading optimization"}; - - // LUT binning (has to match LUTs) - Configurable lutEtaBins{"lutEtaBins", 50, "Number of eta bins in the LUT"}; - Configurable lutEtaMin{"lutEtaMin", -2.5f, "Minimum eta in the LUT"}; - Configurable lutEtaMax{"lutEtaMax", 2.5f, "Maximum eta in the LUT"}; - Configurable lutPBins{"lutPBins", 500, "Number of momentum bins in the LUT"}; - Configurable lutPMin{"lutPMin", 0.0f, "Minimum momentum in the LUT (GeV/c)"}; - Configurable lutPMax{"lutPMax", 10.0f, "Maximum momentum in the LUT (GeV/c)"}; - - // ToT histogram axis - Configurable totHistBins{"totHistBins", 600, "Number of bins on the ToT axis"}; - Configurable totHistMin{"totHistMin", 0.0f, "Minimum of the ToT axis (#mus/10#mum)"}; - Configurable totHistMax{"totHistMax", 10.0f, "Maximum of the ToT axis (#mus/10#mum)"}; + Configurable trackLut{"trackLut", "", "Per-track dE/dx LUT from MakeToTLUT.C (local file or ccdb:)"}; + Configurable seed{"seed", 0, "Seed for the dE/dx sampling (0: different in every job)"}; + Configurable primariesOnly{"primariesOnly", false, "Only physical primaries (otherwise all tracks, including decay products)"}; + Configurable useTrueKinematics{"useTrueKinematics", false, "Fill the histograms vs the true instead of the reconstructed p/z and |eta| (the dE/dx is always sampled at the true ones)"}; + + // dE/dx axis, logarithmic + Configurable numDedxBins{"numDedxBins", 300, "Number of logarithmic dE/dx bins"}; + Configurable dedxMin{"dedxMin", 1e2f, "Minimum of the dE/dx axis in e-"}; + Configurable dedxMax{"dedxMax", 1e5f, "Maximum of the dE/dx axis in e-"}; + + ConfigurableAxis axisEta{"axisEta", {VARIABLE_WIDTH, 0.0f, 0.5f, 1.0f, 1.5f, 2.0f, 2.5f, 3.0f, 4.0f}, "|#eta| binning of the eta-differential histograms"}; // Momentum axis binning choice - Configurable enableLogPBins{"enableLogPBins", true, "Enable logarithmic momentum binning for ToT vs p histograms"}; - Configurable enableLinearPBins{"enableLinearPBins", false, "Enable linear momentum binning for ToT vs p histograms"}; + Configurable numLogBins{"numLogBins", 200, "Number of logarithmic momentum bins"}; + Configurable enableLogPBins{"enableLogPBins", true, "Enable logarithmic momentum binning for dE/dx vs p histograms"}; + Configurable enableLinearPBins{"enableLinearPBins", false, "Enable linear momentum binning for dE/dx vs p histograms"}; Configurable numLinearBins{"numLinearBins", 500, "Number of linear momentum bins"}; Configurable linearPMin{"linearPMin", 0.0f, "Minimum momentum for linear-binned histograms (GeV/c)"}; Configurable linearPMax{"linearPMax", 10.0f, "Maximum momentum for linear-binned histograms (GeV/c)"}; - 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 - }; - - // Configuration defined at init time - o2::fastsim::GeometryContainer mGeoContainer; - float mMagneticField = 0.0f; - void init(o2::framework::InitContext& initContext) + void init(o2::framework::InitContext&) { ccdb->setURL("http://alice-ccdb.cern.ch"); ccdb->setTimestamp(-1); - mGeoContainer.setCcdbManager(ccdb.operator->()); - mGeoContainer.init(initContext); - mMagneticField = mGeoContainer.getFloatValue(0, "global", "magneticfield"); - - 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(maxBarrelLayers.value, - lutEtaBins.value, lutEtaMin.value, lutEtaMax.value, - lutPBins.value, lutPMin.value, lutPMax.value, - analysisEtaMin.value, analysisEtaMax.value, - analysisPMin.value, analysisPMax.value); - - mToTLUT->setCcdbManager(ccdb.operator->()); - - 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) { - LOG(warning) << "Failed to load one or more ToT LUTs. PID results might be incomplete."; + mRandom.SetSeed(seed.value); + if (!mTrackLut.load(trackLut.value, ccdb.operator->(), kPdgCodes)) { + LOG(warning) << "dE/dx LUT incomplete, species without LUT get no tracker PID."; } + bookHistograms(); + } - // Logarithmic momentum bins - mLogBins.clear(); - double pMin = kMinMomentumForLogBins; - double pMax = kMaxMomentumForLogBins; - 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)); + static std::vector logBinning(int nBins, double min, double max) + { + std::vector edges; + const double logMin = std::log10(min); + const double dLog = (std::log10(max) - logMin) / nBins; + for (int i = 0; i <= nBins; ++i) { + edges.push_back(std::pow(10, logMin + i * dLog)); } + return edges; + } - const AxisSpec axisToT{totHistBins.value, static_cast(totHistMin.value), static_cast(totHistMax.value), "ToT (#mus/10#mum)"}; - const AxisSpec axisNsigma{200, -10., 10., "N#sigma"}; - 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"}; - + void bookHistograms() + { if (!enableLogPBins.value && !enableLinearPBins.value) { LOG(fatal) << "At least one of enableLogPBins or enableLinearPBins must be true."; } - 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}); + const std::vector logBins = logBinning(numLogBins.value, kMinMomentumForLogBins, kMaxMomentumForLogBins); + const std::string signalTitle = "d#it{E}/d#it{x} (e^{-})"; + const AxisSpec axisDedx{logBinning(numDedxBins.value, dedxMin.value, dedxMax.value), signalTitle}; + const AxisSpec axisMomentumLog{logBins, "#it{p/z} (GeV/#it{c})"}; + const AxisSpec axisMomentumLin{numLinearBins.value, static_cast(linearPMin.value), static_cast(linearPMax.value), "#it{p/z} (GeV/#it{c})"}; + const AxisSpec axisAbsEta{axisEta, "|#eta|"}; if (enableLogPBins.value) { - const AxisSpec axisMomentumLog{mLogBins, "#it{p/z} (GeV/#it{c})"}; - histos.add("hToTvsP", "ToT vs #it{p/z}; #it{p/z} (GeV/#it{c}); ToT (#mus/10#mum)", kTH2F, {axisMomentumLog, axisToT}); - histos.add("hToTvsMom", "ToT vs #it{p}; #it{p} (GeV/#it{c}); ToT (#mus/10#mum)", kTH2F, {mLogBins, axisToT}); + histos.add("hToTvsP", ("Signal vs #it{p/z}; #it{p/z} (GeV/#it{c}); " + signalTitle).c_str(), kTH2F, {axisMomentumLog, axisDedx}); + histos.add("hToTvsMom", ("Signal vs #it{p}; #it{p} (GeV/#it{c}); " + signalTitle).c_str(), kTH2F, {logBins, axisDedx}); } - if (enableLinearPBins.value) { - const AxisSpec axisMomentumLin{numLinearBins.value, static_cast(linearPMin.value), static_cast(linearPMax.value), "#it{p/z} (GeV/#it{c})"}; - histos.add("hToTvsPLin", "ToT vs #it{p/z} (linear); #it{p/z} (GeV/#it{c}); ToT (#mus/10#mum)", kTH2F, {axisMomentumLin, axisToT}); - histos.add("hToTvsMomLin", "ToT vs #it{p} (linear); #it{p} (GeV/#it{c}); ToT (#mus/10#mum)", kTH2F, {axisMomentumLin, axisToT}); + histos.add("hToTvsPLin", ("Signal vs #it{p/z} (linear); #it{p/z} (GeV/#it{c}); " + signalTitle).c_str(), kTH2F, {axisMomentumLin, axisDedx}); + histos.add("hToTvsMomLin", ("Signal vs #it{p} (linear); #it{p} (GeV/#it{c}); " + signalTitle).c_str(), kTH2F, {axisMomentumLin, axisDedx}); } - std::vector> particleInfo = { - {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}"; - } - + for (int i = 0; i < kNumSpecies; ++i) { if (enableLogPBins.value) { - const AxisSpec axisMomentumLog{mLogBins, "#it{p/z} (GeV/#it{c})"}; - - 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, {axisMomentumLog, 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, {axisMomentumLog, 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}"; - } - - const std::string histName = std::format("NSigma/BarrelNsigmaTrue{}Vs{}Hypothesis", trueName, hypName); - const std::string histTitle = std::format("Nsigma (True {} vs Hyp {}); #it{{p/z}} (GeV/#it{{c}}); N#sigma", trueNamePretty, hypNamePretty); - h2dBarrelNsigmaTrue[iTrue][iHyp] = histos.add(histName.c_str(), histTitle.c_str(), kTH2F, {axisMomentumLog, axisNsigma}); - } + const std::string name = std::format("ToTvsP/hToTvsP_{}", kNames[i]); + const std::string title = std::format("Signal vs #it{{p/z}} for {}; #it{{p/z}} (GeV/#it{{c}}); {}", kPrettyNames[i], signalTitle); + h2dToTvsPperParticle[i] = histos.add(name.c_str(), title.c_str(), kTH2F, {axisMomentumLog, axisDedx}); } - if (enableLinearPBins.value) { - const AxisSpec axisMomentumLin{numLinearBins.value, static_cast(linearPMin.value), static_cast(linearPMax.value), "#it{p/z} (GeV/#it{c})"}; + const std::string name = std::format("ToTvsPLin/hToTvsPLin_{}", kNames[i]); + const std::string title = std::format("Signal vs #it{{p/z}} for {} (linear); #it{{p/z}} (GeV/#it{{c}}); {}", kPrettyNames[i], signalTitle); + h2dToTvsPperParticleLin[i] = histos.add(name.c_str(), title.c_str(), kTH2F, {axisMomentumLin, axisDedx}); + } - const std::string hitsVsPLinName = std::format("HitsPerTrackLin/hHitsPerTrackVsPLin_{}", trueName); - const std::string hitsVsPLinTitle = std::format("N_hits vs #it{{p/z}} for {} (linear); #it{{p/z}} (GeV/#it{{c}}); N_hits", trueNamePretty); - h2dHitsPerTrackVsPLin[iTrue] = histos.add(hitsVsPLinName.c_str(), hitsVsPLinTitle.c_str(), kTH2F, {axisMomentumLin, axisHitsPerTrack}); + const std::string nameEta = std::format("ToTvsPEta/hToTvsPEta_{}", kNames[i]); + const std::string titleEta = std::format("Signal vs #it{{p/z}} and |#eta| for {}; #it{{p/z}} (GeV/#it{{c}}); |#eta|; {}", kPrettyNames[i], signalTitle); + h3dToTvsPEtaPerParticle[i] = histos.add(nameEta.c_str(), titleEta.c_str(), kTH3F, {axisMomentumLog, axisAbsEta, axisDedx}); - const std::string totVsPLinName = std::format("ToTvsPLin/hToTvsPLin_{}", trueName); - const std::string totVsPLinTitle = std::format("ToT vs #it{{p/z}} for {} (linear); #it{{p/z}} (GeV/#it{{c}}); ToT (#mus/10#mum)", trueNamePretty); - h2dToTvsPperParticleLin[iTrue] = histos.add(totVsPLinName.c_str(), totVsPLinTitle.c_str(), kTH2F, {axisMomentumLin, axisToT}); - } + const std::string nameTracks = std::format("Bookkeeping/hTracks_{}", kNames[i]); + const std::string titleTracks = std::format("Reconstructed {}; #it{{p/z}} (GeV/#it{{c}}); |#eta|", kPrettyNames[i]); + h2dTracks[i] = histos.add(nameTracks.c_str(), titleTracks.c_str(), kTH2F, {axisMomentumLog, axisAbsEta}); + + const std::string nameWithDedx = std::format("Bookkeeping/hTracksWithDedx_{}", kNames[i]); + const std::string titleWithDedx = std::format("Reconstructed {} with dE/dx; #it{{p/z}} (GeV/#it{{c}}); |#eta|", kPrettyNames[i]); + h2dTracksWithDedx[i] = histos.add(nameWithDedx.c_str(), titleWithDedx.c_str(), kTH2F, {axisMomentumLog, axisAbsEta}); } } - void process(soa::Join::iterator const& collision, - soa::Join const& tracks, - aod::McParticles const& /*mcParticles*/, - aod::McCollisions const& /*mcCollisions*/) + void process(soa::Join const& tracks, + aod::McParticles const&) { - 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()); - } - for (const auto& track : tracks) { - float truncatedMeanToT = -1.0f; - std::array nSigmaValues{}; - nSigmaValues.fill(999.f); - 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]); continue; } - const auto& mcParticle = track.mcParticle(); - - const auto& pdgInfo = pdgDatabase->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]); + if (primariesOnly.value && !mcParticle.isPhysicalPrimary()) { continue; } - - const float p = mcParticle.p(); - const float eta = mcParticle.eta(); - const int truePdgCode = std::abs(mcParticle.pdgCode()); - - float rigidity = p; - if (pdgInfo) { - const float trueCharge = std::abs(pdgInfo->Charge()) / 3.0f; - if (trueCharge > 0) { - rigidity = p / trueCharge; - } - } - - 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]); + const int pdg = std::abs(mcParticle.pdgCode()); + const auto itSpecies = std::find(kPdgCodes.begin(), kPdgCodes.end(), pdg); + if (itSpecies == kPdgCodes.end()) { continue; } - - const int binnedP = mToTLUT->getPBin(p); - const int binnedEta = mToTLUT->getEtaBin(std::abs(eta)); - - uint16_t hitMap = 0; - int nHitLayers = 0; - o2::track::TrackParCov o2track = o2::upgrade::convertMCParticleToO2Track(mcParticle, pdgDatabase); - - float xPv = -100.f; - static constexpr float kTrkXThreshold = -99.f; - if (o2track.propagateToDCA(mcPvVtx, mMagneticField)) { - xPv = o2track.getX(); + const auto* pdgInfo = pdgDatabase->GetParticle(mcParticle.pdgCode()); + if (!pdgInfo || pdgInfo->Charge() == 0) { + continue; } - - if (xPv > kTrkXThreshold) { - for (int layer = 0; layer < maxBarrelLayers.value; ++layer) { - float layerRadius = kTrackerRadii[layer]; - float trackLength = computeTrackLength(o2track, layerRadius, mMagneticField); - - if (trackLength > 0) { - hitMap |= (1 << layer); - histos.fill(HIST("hHitLayers"), layer); - nHitLayers++; - } - } + const int iSpecies = std::distance(kPdgCodes.begin(), itSpecies); + const float absCharge = std::abs(pdgInfo->Charge()) / 3.f; + const float rigidityTrue = mcParticle.p() / absCharge; + + // the OTF tracks carry p/z + const float rigidity = useTrueKinematics.value ? rigidityTrue : track.p(); + const float absEta = std::abs(useTrueKinematics.value ? mcParticle.eta() : track.eta()); + const float p = rigidity * absCharge; + h2dTracks[iSpecies]->Fill(rigidity, absEta); + + const float dedx = mTrackLut.sample(pdg, mcParticle.eta(), rigidityTrue, mRandom.Rndm()); + if (dedx < 0.f) { + continue; } - - histos.fill(HIST("hHitMultiplicity"), nHitLayers); + h2dTracksWithDedx[iSpecies]->Fill(rigidity, absEta); + h3dToTvsPEtaPerParticle[iSpecies]->Fill(rigidity, absEta, dedx); if (enableLogPBins.value) { - h2dHitsPerTrackVsP[truePdgIdx]->Fill(rigidity, nHitLayers); + histos.fill(HIST("hToTvsMom"), p, dedx); + histos.fill(HIST("hToTvsP"), rigidity, dedx); + h2dToTvsPperParticle[iSpecies]->Fill(rigidity, dedx); } if (enableLinearPBins.value) { - h2dHitsPerTrackVsPLin[truePdgIdx]->Fill(rigidity, nHitLayers); - } - - std::vector validToTs; - - for (int layer = kMinLayerForTruncation; layer < maxBarrelLayers.value; ++layer) { - if (((hitMap >> layer) & 0x1) != 0) { - TH1F* totHist = mToTLUT->getHistogramForSampling(truePdgIdx, layer, binnedEta, binnedP); - - if (totHist && totHist->GetEntries() > 1) { - float sampledToT = totHist->GetRandom(); - validToTs.push_back(sampledToT); - } - } + histos.fill(HIST("hToTvsMomLin"), p, dedx); + histos.fill(HIST("hToTvsPLin"), rigidity, dedx); + h2dToTvsPperParticleLin[iSpecies]->Fill(rigidity, dedx); } - - truncatedMeanToT = -1.0f; - const size_t nValid = validToTs.size(); - size_t nUse = 0; - - 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) { - 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); - - if (enableLogPBins.value) { - histos.fill(HIST("hToTvsMom"), p, truncatedMeanToT); - histos.fill(HIST("hToTvsP"), rigidity, truncatedMeanToT); - h2dToTvsPperParticle[truePdgIdx]->Fill(rigidity, truncatedMeanToT); - } - if (enableLinearPBins.value) { - histos.fill(HIST("hToTvsMomLin"), p, truncatedMeanToT); - histos.fill(HIST("hToTvsPLin"), rigidity, truncatedMeanToT); - h2dToTvsPperParticleLin[truePdgIdx]->Fill(rigidity, truncatedMeanToT); - } - } - - nSigmaValues.fill(999.f); - - if (enableLogPBins.value && truncatedMeanToT > 0) { - for (size_t iHyp = 0; iHyp < mHypothesisPdgCodes.size(); ++iHyp) { - int hypPdgCode = mHypothesisPdgCodes[iHyp]; - int hypPdgIdx = mToTLUT->getPdgIndex(hypPdgCode); - - if (hypPdgIdx == -1) { - nSigmaValues[iHyp] = 999.f; - continue; - } - - float expectedToT = getToTMeanFromMomentumSlice(h2dToTvsPperParticle[hypPdgIdx], rigidity); - float resolution = getToTResolutionFromMomentumSlice(h2dToTvsPperParticle[hypPdgIdx], rigidity); - - if (expectedToT > 0 && resolution > 0) { - nSigmaValues[iHyp] = calculateNsigma(truncatedMeanToT, expectedToT, resolution); - h2dBarrelNsigmaTrue[truePdgIdx][iHyp]->Fill(rigidity, nSigmaValues[iHyp]); - } else { - 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]); } } }; From 8fa339e27e663bd0921eb3675b952d4490decedf Mon Sep 17 00:00:00 2001 From: Henrik Fribert Date: Wed, 14 Oct 2026 18:07:44 +0200 Subject: [PATCH 2/5] Fix linter --- .../TableProducer/OTF/onTheFlyTrackerPid.cxx | 56 +++++++++---------- 1 file changed, 28 insertions(+), 28 deletions(-) diff --git a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx index 0d65ef5f642..b449e533178 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx @@ -56,8 +56,8 @@ using namespace o2; using namespace o2::framework; -static constexpr double kMinMomentumForLogBins = 0.05; -static constexpr double kMaxMomentumForLogBins = 10.0; +static constexpr double MinMomentumForLogBins = 0.05; +static constexpr double MaxMomentumForLogBins = 10.0; /// Local path of a LUT file; "ccdb:" is downloaded once to /tmp/ToTLUTs//snapshot.root std::string getLocalLutFile(const std::string& filename, o2::ccdb::BasicCCDBManager* ccdbManager) @@ -114,7 +114,7 @@ class TrackDedxLut const size_t nDedx = mDedxEdges.size() - 1; bool allLoaded = true; - for (const int pdg : pdgCodes) { + for (const auto& pdg : pdgCodes) { std::unique_ptr> h(f->Get>(Form("h_%d", pdg))); if (!h || h->size() != nBins * nDedx) { LOG(warning) << "No dE/dx LUT for PDG " << pdg << " in " << filename; @@ -205,16 +205,16 @@ class TrackDedxLut }; struct OnTheFlyTrackerPid { - static constexpr int kNumSpecies = 9; - static constexpr std::array kPdgCodes = {11, 13, 211, 321, 2212, 1000010020, 1000010030, 1000020030, 1000020040}; - static constexpr std::array kNames = {"Elec", "Muon", "Pion", "Kaon", "Prot", "Deut", "Trit", "He3", "Al"}; - static constexpr std::array kPrettyNames = {"#it{e}", "#it{#mu}", "#it{#pi}", "#it{K}", "#it{p}", "#it{d}", "#it{t}", "#it{^{3}He}", "#it{^{4}He}"}; + static constexpr int NumSpecies = 9; + static constexpr std::array PdgCodes = {11, 13, 211, 321, 2212, 1000010020, 1000010030, 1000020030, 1000020040}; + static constexpr std::array SpeciesNames = {"Elec", "Muon", "Pion", "Kaon", "Prot", "Deut", "Trit", "He3", "Al"}; + static constexpr std::array SpeciesLabels = {"#it{e}", "#it{#mu}", "#it{#pi}", "#it{K}", "#it{p}", "#it{d}", "#it{t}", "#it{^{3}He}", "#it{^{4}He}"}; - std::array, kNumSpecies> h2dToTvsPperParticle; - std::array, kNumSpecies> h2dToTvsPperParticleLin; - std::array, kNumSpecies> h3dToTvsPEtaPerParticle; - std::array, kNumSpecies> h2dTracks; - std::array, kNumSpecies> h2dTracksWithDedx; + std::array, NumSpecies> h2dToTvsPperParticle; + std::array, NumSpecies> h2dToTvsPperParticleLin; + std::array, NumSpecies> h3dToTvsPEtaPerParticle; + std::array, NumSpecies> h2dTracks; + std::array, NumSpecies> h2dTracksWithDedx; Service pdgDatabase{}; Service ccdb{}; @@ -249,7 +249,7 @@ struct OnTheFlyTrackerPid { ccdb->setTimestamp(-1); mRandom.SetSeed(seed.value); - if (!mTrackLut.load(trackLut.value, ccdb.operator->(), kPdgCodes)) { + if (!mTrackLut.load(trackLut.value, ccdb.operator->(), PdgCodes)) { LOG(warning) << "dE/dx LUT incomplete, species without LUT get no tracker PID."; } bookHistograms(); @@ -272,7 +272,7 @@ struct OnTheFlyTrackerPid { LOG(fatal) << "At least one of enableLogPBins or enableLinearPBins must be true."; } - const std::vector logBins = logBinning(numLogBins.value, kMinMomentumForLogBins, kMaxMomentumForLogBins); + const std::vector logBins = logBinning(numLogBins.value, MinMomentumForLogBins, MaxMomentumForLogBins); const std::string signalTitle = "d#it{E}/d#it{x} (e^{-})"; const AxisSpec axisDedx{logBinning(numDedxBins.value, dedxMin.value, dedxMax.value), signalTitle}; const AxisSpec axisMomentumLog{logBins, "#it{p/z} (GeV/#it{c})"}; @@ -288,28 +288,28 @@ struct OnTheFlyTrackerPid { histos.add("hToTvsMomLin", ("Signal vs #it{p} (linear); #it{p} (GeV/#it{c}); " + signalTitle).c_str(), kTH2F, {axisMomentumLin, axisDedx}); } - for (int i = 0; i < kNumSpecies; ++i) { + for (int i = 0; i < NumSpecies; ++i) { if (enableLogPBins.value) { - const std::string name = std::format("ToTvsP/hToTvsP_{}", kNames[i]); - const std::string title = std::format("Signal vs #it{{p/z}} for {}; #it{{p/z}} (GeV/#it{{c}}); {}", kPrettyNames[i], signalTitle); + const std::string name = std::format("ToTvsP/hToTvsP_{}", SpeciesNames[i]); + const std::string title = std::format("Signal vs #it{{p/z}} for {}; #it{{p/z}} (GeV/#it{{c}}); {}", SpeciesLabels[i], signalTitle); h2dToTvsPperParticle[i] = histos.add(name.c_str(), title.c_str(), kTH2F, {axisMomentumLog, axisDedx}); } if (enableLinearPBins.value) { - const std::string name = std::format("ToTvsPLin/hToTvsPLin_{}", kNames[i]); - const std::string title = std::format("Signal vs #it{{p/z}} for {} (linear); #it{{p/z}} (GeV/#it{{c}}); {}", kPrettyNames[i], signalTitle); + const std::string name = std::format("ToTvsPLin/hToTvsPLin_{}", SpeciesNames[i]); + const std::string title = std::format("Signal vs #it{{p/z}} for {} (linear); #it{{p/z}} (GeV/#it{{c}}); {}", SpeciesLabels[i], signalTitle); h2dToTvsPperParticleLin[i] = histos.add(name.c_str(), title.c_str(), kTH2F, {axisMomentumLin, axisDedx}); } - const std::string nameEta = std::format("ToTvsPEta/hToTvsPEta_{}", kNames[i]); - const std::string titleEta = std::format("Signal vs #it{{p/z}} and |#eta| for {}; #it{{p/z}} (GeV/#it{{c}}); |#eta|; {}", kPrettyNames[i], signalTitle); + const std::string nameEta = std::format("ToTvsPEta/hToTvsPEta_{}", SpeciesNames[i]); + const std::string titleEta = std::format("Signal vs #it{{p/z}} and |#eta| for {}; #it{{p/z}} (GeV/#it{{c}}); |#eta|; {}", SpeciesLabels[i], signalTitle); h3dToTvsPEtaPerParticle[i] = histos.add(nameEta.c_str(), titleEta.c_str(), kTH3F, {axisMomentumLog, axisAbsEta, axisDedx}); - const std::string nameTracks = std::format("Bookkeeping/hTracks_{}", kNames[i]); - const std::string titleTracks = std::format("Reconstructed {}; #it{{p/z}} (GeV/#it{{c}}); |#eta|", kPrettyNames[i]); + const std::string nameTracks = std::format("Bookkeeping/hTracks_{}", SpeciesNames[i]); + const std::string titleTracks = std::format("Reconstructed {}; #it{{p/z}} (GeV/#it{{c}}); |#eta|", SpeciesLabels[i]); h2dTracks[i] = histos.add(nameTracks.c_str(), titleTracks.c_str(), kTH2F, {axisMomentumLog, axisAbsEta}); - const std::string nameWithDedx = std::format("Bookkeeping/hTracksWithDedx_{}", kNames[i]); - const std::string titleWithDedx = std::format("Reconstructed {} with dE/dx; #it{{p/z}} (GeV/#it{{c}}); |#eta|", kPrettyNames[i]); + const std::string nameWithDedx = std::format("Bookkeeping/hTracksWithDedx_{}", SpeciesNames[i]); + const std::string titleWithDedx = std::format("Reconstructed {} with dE/dx; #it{{p/z}} (GeV/#it{{c}}); |#eta|", SpeciesLabels[i]); h2dTracksWithDedx[i] = histos.add(nameWithDedx.c_str(), titleWithDedx.c_str(), kTH2F, {axisMomentumLog, axisAbsEta}); } } @@ -326,15 +326,15 @@ struct OnTheFlyTrackerPid { continue; } const int pdg = std::abs(mcParticle.pdgCode()); - const auto itSpecies = std::find(kPdgCodes.begin(), kPdgCodes.end(), pdg); - if (itSpecies == kPdgCodes.end()) { + const auto itSpecies = std::find(PdgCodes.begin(), PdgCodes.end(), pdg); + if (itSpecies == PdgCodes.end()) { continue; } const auto* pdgInfo = pdgDatabase->GetParticle(mcParticle.pdgCode()); if (!pdgInfo || pdgInfo->Charge() == 0) { continue; } - const int iSpecies = std::distance(kPdgCodes.begin(), itSpecies); + const int iSpecies = std::distance(PdgCodes.begin(), itSpecies); const float absCharge = std::abs(pdgInfo->Charge()) / 3.f; const float rigidityTrue = mcParticle.p() / absCharge; From 80942b6cb81c947faab2ce82de3a74cde8a0a046 Mon Sep 17 00:00:00 2001 From: Henrik Fribert Date: Thu, 15 Oct 2026 09:09:08 +0200 Subject: [PATCH 3/5] Retrigger CI From 85cee1ba03663980256de95c421802270bcc0b03 Mon Sep 17 00:00:00 2001 From: Henrik Fribert Date: Thu, 15 Oct 2026 10:43:05 +0200 Subject: [PATCH 4/5] Include for std::distance --- 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 b449e533178..2dd3d72ef33 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx @@ -47,6 +47,7 @@ #include #include #include +#include #include #include #include From d77c23c587c1125505a2e6712997eca62ead8e53 Mon Sep 17 00:00:00 2001 From: Henrik Fribert Date: Thu, 15 Oct 2026 13:33:16 +0200 Subject: [PATCH 5/5] Add back tracker PID tables --- .../TableProducer/OTF/onTheFlyTrackerPid.cxx | 184 ++++++++++++------ 1 file changed, 127 insertions(+), 57 deletions(-) diff --git a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx index 2dd3d72ef33..51a8dcb95db 100644 --- a/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx +++ b/ALICE3/TableProducer/OTF/onTheFlyTrackerPid.cxx @@ -21,6 +21,8 @@ /// \since May 22, 2025 /// +#include "ALICE3/DataModel/OTFPIDTrk.h" + #include #include #include @@ -122,15 +124,22 @@ class TrackDedxLut allLoaded = false; continue; } - // cumulative distributions, only for the filled bins + // cumulative distributions and raw moments, only for the filled bins auto& offset = mOffset[pdg]; auto& cdf = mCdf[pdg]; + auto& mean = mMean[pdg]; + auto& sigma = mSigma[pdg]; offset.assign(nBins, -1); + mean.assign(nBins, -1.f); + sigma.assign(nBins, -1.f); for (size_t b = 0; b < nBins; ++b) { const float* hb = h->data() + b * nDedx; - double sum = 0.; + double sum = 0., m1 = 0., m2 = 0.; for (size_t k = 0; k < nDedx; ++k) { + const double x = 0.5 * (mDedxEdges[k] + mDedxEdges[k + 1]); sum += hb[k]; + m1 += hb[k] * x; + m2 += hb[k] * x * x; } if (sum <= 0.) { continue; @@ -141,6 +150,9 @@ class TrackDedxLut cum += hb[k]; cdf.push_back(cum / sum); } + const double mu = m1 / sum; + mean[b] = mu; + sigma[b] = std::sqrt(std::max(0., m2 / sum - mu * mu)); } } LOG(info) << "Loaded dE/dx LUT " << filename << " (" << mEtaEdges.size() - 1 << " |eta| x " << mRigEdges.size() - 1 @@ -150,34 +162,67 @@ class TrackDedxLut /// Sampled dE/dx for uniform random number u, or -1 if there is no LUT entry float sample(int pdg, float eta, float rigidity, double u) const + { + int bin = -1, binNext = -1; + float weight = 0.f; + if (!locate(pdg, eta, rigidity, bin, binNext, weight)) { + return -1.f; + } + const auto& offset = mOffset.at(pdg); + const float* cdf = mCdf.at(pdg).data(); + float dedx = inverseCdf(cdf + offset[bin], u); + if (binNext >= 0) { + dedx += weight * (inverseCdf(cdf + offset[binNext], u) - dedx); + } + return dedx; + } + + bool getMeanSigma(int pdg, float eta, float rigidity, float& mean, float& sigma) const + { + int bin = -1, binNext = -1; + float weight = 0.f; + if (!locate(pdg, eta, rigidity, bin, binNext, weight)) { + return false; + } + const auto& m = mMean.at(pdg); + const auto& s = mSigma.at(pdg); + mean = m[bin]; + sigma = s[bin]; + if (binNext >= 0) { + mean += weight * (m[binNext] - mean); + sigma += weight * (s[binNext] - sigma); + } + return sigma > 0.f; + } + + private: + /// LUT bin of the track and the neighbouring rigidity bin on its side (-1 if empty) + bool locate(int pdg, float eta, float rigidity, int& bin, int& binNext, float& weight) const { const auto it = mOffset.find(pdg); const float absEta = std::abs(eta); if (it == mOffset.end() || mEtaEdges.empty() || absEta < mEtaEdges.front() || absEta >= mEtaEdges.back() || rigidity < mRigEdges.front() || rigidity >= mRigEdges.back()) { - return -1.f; + return false; } const int nRig = mRigEdges.size() - 1; const int iEta = std::upper_bound(mEtaEdges.begin(), mEtaEdges.end(), absEta) - mEtaEdges.begin() - 1; const int iRig = std::upper_bound(mRigEdges.begin(), mRigEdges.end(), rigidity) - mRigEdges.begin() - 1; - const int offset = it->second[iEta * nRig + iRig]; - if (offset < 0) { - return -1.f; + bin = iEta * nRig + iRig; + if (it->second[bin] < 0) { + return false; } - const float* cdf = mCdf.at(pdg).data(); - float dedx = inverseCdf(cdf + offset, u); - - // neighbouring bin on the side of the track, no interpolation if it is empty + binNext = -1; + weight = 0.f; const double logRig = std::log(rigidity); const int iNext = logRig < mLogRigCentres[iRig] ? iRig - 1 : iRig + 1; if (iNext >= 0 && iNext < nRig && it->second[iEta * nRig + iNext] >= 0) { - const float w = (logRig - mLogRigCentres[iRig]) / (mLogRigCentres[iNext] - mLogRigCentres[iRig]); - dedx += w * (inverseCdf(cdf + it->second[iEta * nRig + iNext], u) - dedx); + binNext = iEta * nRig + iNext; + weight = (logRig - mLogRigCentres[iRig]) / (mLogRigCentres[iNext] - mLogRigCentres[iRig]); } - return dedx; + return true; } - private: /// dE/dx at which the cumulative distribution reaches u, uniform within the dE/dx bin float inverseCdf(const float* cdf, double u) const { @@ -203,6 +248,8 @@ class TrackDedxLut std::vector mDedxEdges; std::unordered_map> mOffset; // |pdg| -> start of the cdf of bin etaBin * nRig + rigBin, -1 if empty std::unordered_map> mCdf; // |pdg| -> cumulative dE/dx distributions of the filled bins + std::unordered_map> mMean; // |pdg| -> mean dE/dx per bin, -1 if empty + std::unordered_map> mSigma; }; struct OnTheFlyTrackerPid { @@ -210,6 +257,7 @@ struct OnTheFlyTrackerPid { static constexpr std::array PdgCodes = {11, 13, 211, 321, 2212, 1000010020, 1000010030, 1000020030, 1000020040}; static constexpr std::array SpeciesNames = {"Elec", "Muon", "Pion", "Kaon", "Prot", "Deut", "Trit", "He3", "Al"}; static constexpr std::array SpeciesLabels = {"#it{e}", "#it{#mu}", "#it{#pi}", "#it{K}", "#it{p}", "#it{d}", "#it{t}", "#it{^{3}He}", "#it{^{4}He}"}; + static constexpr float NoNsigma = 999.f; std::array, NumSpecies> h2dToTvsPperParticle; std::array, NumSpecies> h2dToTvsPperParticleLin; @@ -217,6 +265,9 @@ struct OnTheFlyTrackerPid { std::array, NumSpecies> h2dTracks; std::array, NumSpecies> h2dTracksWithDedx; + Produces tableUpgradeTrkPidSignals; + Produces tableUpgradeTrkPids; + Service pdgDatabase{}; Service ccdb{}; TrackDedxLut mTrackLut; @@ -226,7 +277,7 @@ struct OnTheFlyTrackerPid { Configurable trackLut{"trackLut", "", "Per-track dE/dx LUT from MakeToTLUT.C (local file or ccdb:)"}; Configurable seed{"seed", 0, "Seed for the dE/dx sampling (0: different in every job)"}; - Configurable primariesOnly{"primariesOnly", false, "Only physical primaries (otherwise all tracks, including decay products)"}; + Configurable primariesOnly{"primariesOnly", false, "Fill the histograms only with physical primaries (the tables are filled for all tracks)"}; Configurable useTrueKinematics{"useTrueKinematics", false, "Fill the histograms vs the true instead of the reconstructed p/z and |eta| (the dE/dx is always sampled at the true ones)"}; // dE/dx axis, logarithmic @@ -315,52 +366,71 @@ struct OnTheFlyTrackerPid { } } + /// dE/dx of the track (-1 if there is none) and its nsigma for every species + template + float processTrack(TTrack const& track, std::array& nSigma) + { + if (!track.has_mcParticle()) { + return -1.f; + } + const auto& mcParticle = track.mcParticle(); + const int pdg = std::abs(mcParticle.pdgCode()); + const auto itSpecies = std::find(PdgCodes.begin(), PdgCodes.end(), pdg); + if (itSpecies == PdgCodes.end()) { + return -1.f; + } + const auto* pdgInfo = pdgDatabase->GetParticle(mcParticle.pdgCode()); + if (!pdgInfo || pdgInfo->Charge() == 0) { + return -1.f; + } + const int iSpecies = std::distance(PdgCodes.begin(), itSpecies); + const float absCharge = std::abs(pdgInfo->Charge()) / 3.f; + const float rigidityTrue = mcParticle.p() / absCharge; + const float dedx = mTrackLut.sample(pdg, mcParticle.eta(), rigidityTrue, mRandom.Rndm()); + + if (dedx > 0.f) { + for (int i = 0; i < NumSpecies; ++i) { + float mean = 0.f, sigma = 0.f; + if (mTrackLut.getMeanSigma(PdgCodes[i], track.eta(), track.p(), mean, sigma)) { + nSigma[i] = (dedx - mean) / sigma; + } + } + } + + if (primariesOnly.value && !mcParticle.isPhysicalPrimary()) { + return dedx; + } + const float rigidity = useTrueKinematics.value ? rigidityTrue : track.p(); + const float absEta = std::abs(useTrueKinematics.value ? mcParticle.eta() : track.eta()); + const float p = rigidity * absCharge; + h2dTracks[iSpecies]->Fill(rigidity, absEta); + if (dedx < 0.f) { + return dedx; + } + h2dTracksWithDedx[iSpecies]->Fill(rigidity, absEta); + h3dToTvsPEtaPerParticle[iSpecies]->Fill(rigidity, absEta, dedx); + if (enableLogPBins.value) { + histos.fill(HIST("hToTvsMom"), p, dedx); + histos.fill(HIST("hToTvsP"), rigidity, dedx); + h2dToTvsPperParticle[iSpecies]->Fill(rigidity, dedx); + } + if (enableLinearPBins.value) { + histos.fill(HIST("hToTvsMomLin"), p, dedx); + histos.fill(HIST("hToTvsPLin"), rigidity, dedx); + h2dToTvsPperParticleLin[iSpecies]->Fill(rigidity, dedx); + } + return dedx; + } + void process(soa::Join const& tracks, aod::McParticles const&) { + std::array nSigma{}; for (const auto& track : tracks) { - if (!track.has_mcParticle()) { - continue; - } - const auto& mcParticle = track.mcParticle(); - if (primariesOnly.value && !mcParticle.isPhysicalPrimary()) { - continue; - } - const int pdg = std::abs(mcParticle.pdgCode()); - const auto itSpecies = std::find(PdgCodes.begin(), PdgCodes.end(), pdg); - if (itSpecies == PdgCodes.end()) { - continue; - } - const auto* pdgInfo = pdgDatabase->GetParticle(mcParticle.pdgCode()); - if (!pdgInfo || pdgInfo->Charge() == 0) { - continue; - } - const int iSpecies = std::distance(PdgCodes.begin(), itSpecies); - const float absCharge = std::abs(pdgInfo->Charge()) / 3.f; - const float rigidityTrue = mcParticle.p() / absCharge; - - // the OTF tracks carry p/z - const float rigidity = useTrueKinematics.value ? rigidityTrue : track.p(); - const float absEta = std::abs(useTrueKinematics.value ? mcParticle.eta() : track.eta()); - const float p = rigidity * absCharge; - h2dTracks[iSpecies]->Fill(rigidity, absEta); - - const float dedx = mTrackLut.sample(pdg, mcParticle.eta(), rigidityTrue, mRandom.Rndm()); - if (dedx < 0.f) { - continue; - } - h2dTracksWithDedx[iSpecies]->Fill(rigidity, absEta); - h3dToTvsPEtaPerParticle[iSpecies]->Fill(rigidity, absEta, dedx); - if (enableLogPBins.value) { - histos.fill(HIST("hToTvsMom"), p, dedx); - histos.fill(HIST("hToTvsP"), rigidity, dedx); - h2dToTvsPperParticle[iSpecies]->Fill(rigidity, dedx); - } - if (enableLinearPBins.value) { - histos.fill(HIST("hToTvsMomLin"), p, dedx); - histos.fill(HIST("hToTvsPLin"), rigidity, dedx); - h2dToTvsPperParticleLin[iSpecies]->Fill(rigidity, dedx); - } + nSigma.fill(NoNsigma); + const float dedx = processTrack(track, nSigma); + tableUpgradeTrkPidSignals(dedx); + tableUpgradeTrkPids(nSigma[0], nSigma[1], nSigma[2], nSigma[3], nSigma[4], nSigma[5], nSigma[6], nSigma[7], nSigma[8]); } } };