From fc73bd91103f2b23327323d4dd88586621b31551 Mon Sep 17 00:00:00 2001 From: Anton Alkin Date: Wed, 8 Apr 2026 12:22:35 +0200 Subject: [PATCH 01/11] Flat LUT structure --- ALICE3/Core/FlatLutEntry.cxx | 174 +++++++++++++++++++++++++++++++++++ ALICE3/Core/FlatLutEntry.h | 158 +++++++++++++++++++++++++++++++ 2 files changed, 332 insertions(+) create mode 100644 ALICE3/Core/FlatLutEntry.cxx create mode 100644 ALICE3/Core/FlatLutEntry.h diff --git a/ALICE3/Core/FlatLutEntry.cxx b/ALICE3/Core/FlatLutEntry.cxx new file mode 100644 index 00000000000..f1b77d3c056 --- /dev/null +++ b/ALICE3/Core/FlatLutEntry.cxx @@ -0,0 +1,174 @@ +// 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. + +#include "FlatLutEntry.h" + +#include + +#include + +namespace o2::delphes +{ + +float map_t::fracPositionWithinBin(float val) const +{ + float width = (max - min) / nbins; + int bin; + float returnVal = 0.5f; + if (log) { + bin = static_cast((std::log10(val) - min) / width); + returnVal = ((std::log10(val) - min) / width) - bin; + } else { + bin = static_cast((val - min) / width); + returnVal = val / width - bin; + } + return returnVal; +} + +int map_t::find(float val) const +{ + float width = (max - min) / nbins; + int bin; + if (log) { + bin = static_cast((std::log10(val) - min) / width); + } else { + bin = static_cast((val - min) / width); + } + if (bin < 0) { + return 0; + } + if (bin > nbins - 1) { + return nbins - 1; + } + return bin; +} + +void map_t::print() const +{ + LOGF(info, "nbins = %d, min = %f, max = %f, log = %s \n", nbins, min, max, log ? "on" : "off"); +} + +bool lutHeader_t::check_version() const +{ + return (version == LUTCOVM_VERSION); +} + +void lutHeader_t::print() const +{ + LOGF(info, " version: %d \n", version); + LOGF(info, " pdg: %d \n", pdg); + LOGF(info, " field: %f \n", field); + LOGF(info, " nchmap: "); + nchmap.print(); + LOGF(info, " radmap: "); + radmap.print(); + LOGF(info, " etamap: "); + etamap.print(); + LOGF(info, " ptmap: "); + ptmap.print(); +} + +void FlatLutData::initialize(const lutHeader_t& header) +{ + mNchBins = header.nchmap.nbins; + mRadBins = header.radmap.nbins; + mEtaBins = header.etamap.nbins; + mPtBins = header.ptmap.nbins; + + size_t headerSize = sizeof(lutHeader_t); + size_t numEntries = static_cast(mNchBins) * mRadBins * mEtaBins * mPtBins; + size_t entriesSize = numEntries * sizeof(lutEntry_t); + size_t totalSize = headerSize + entriesSize; + + mData.resize(totalSize); + + // Write header at the beginning + std::memcpy(mData.data(), &header, headerSize); +} + +size_t FlatLutData::getEntryOffset(int nch_bin, int rad_bin, int eta_bin, int pt_bin) const +{ + size_t headerSize = sizeof(lutHeader_t); + + // Linear index: nch varies slowest, pt varies fastest + // idx = nch * (rad*eta*pt) + rad * (eta*pt) + eta * pt + pt + size_t linearIdx = static_cast(nch_bin) * (mRadBins * mEtaBins * mPtBins) + static_cast(rad_bin) * (mEtaBins * mPtBins) + static_cast(eta_bin) * mPtBins + static_cast(pt_bin); + + return headerSize + linearIdx * sizeof(lutEntry_t); +} + +lutEntry_t* FlatLutData::getEntry(int nch_bin, int rad_bin, int eta_bin, int pt_bin) +{ + size_t offset = getEntryOffset(nch_bin, rad_bin, eta_bin, pt_bin); + return reinterpret_cast(mData.data() + offset); +} + +const lutEntry_t* FlatLutData::getEntry(int nch_bin, int rad_bin, int eta_bin, int pt_bin) const +{ + size_t offset = getEntryOffset(nch_bin, rad_bin, eta_bin, pt_bin); + return reinterpret_cast(mData.data() + offset); +} + +FlatLutData FlatLutData::fromBuffer(const uint8_t* buffer, size_t size) +{ + FlatLutData data; + // Validate buffer + if (size < sizeof(lutHeader_t)) { + LOG(fatal) << "Buffer too small for LUT header"; + } + + const auto* header = reinterpret_cast(buffer); + data.mNchBins = header->nchmap.nbins; + data.mRadBins = header->radmap.nbins; + data.mEtaBins = header->etamap.nbins; + data.mPtBins = header->ptmap.nbins; + + size_t expectedSize = sizeof(lutHeader_t) + static_cast(data.mNchBins) * data.mRadBins * data.mEtaBins * data.mPtBins * sizeof(lutEntry_t); + + if (size < expectedSize) { + LOG(fatal) << "Buffer size mismatch: expected " << expectedSize << ", got " << size; + } + + // Copy buffer + data.mData.resize(size); + std::memcpy(data.mData.data(), buffer, size); + + return data; +} + +FlatLutData FlatLutData::fromExternalBuffer(uint8_t* buffer, size_t size) +{ + FlatLutData data; + // Validate buffer + if (size < sizeof(lutHeader_t)) { + LOG(fatal) << "Buffer too small for LUT header"; + } + + const auto* header = reinterpret_cast(buffer); + data.mNchBins = header->nchmap.nbins; + data.mRadBins = header->radmap.nbins; + data.mEtaBins = header->etamap.nbins; + data.mPtBins = header->ptmap.nbins; + + size_t expectedSize = sizeof(lutHeader_t) + static_cast(data.mNchBins) * data.mRadBins * data.mEtaBins * data.mPtBins * sizeof(lutEntry_t); + + if (size < expectedSize) { + LOG(fatal) << "Buffer size mismatch: expected " << expectedSize << " got " << size; + } + + // Store reference to external buffer (no copy) + // WARNING: Caller must ensure buffer lifetime exceeds FlatLutData usage + data.mData.assign(buffer, buffer + size); + + return data; +} + +} // namespace o2::delphes diff --git a/ALICE3/Core/FlatLutEntry.h b/ALICE3/Core/FlatLutEntry.h new file mode 100644 index 00000000000..18a2db63680 --- /dev/null +++ b/ALICE3/Core/FlatLutEntry.h @@ -0,0 +1,158 @@ +// 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. + +#ifndef ALICE3_CORE_FLATLUTENTRY_H_ +#define ALICE3_CORE_FLATLUTENTRY_H_ + +#include +#include +#include + +#define LUTCOVM_VERSION 20260408 + +namespace o2::delphes +{ + +/** + * @brief Flat LUT entry structure + */ +struct lutEntry_t { + float nch = 0.f; + float eta = 0.f; + float pt = 0.f; + bool valid = false; + float eff = 0.f; + float eff2 = 0.f; + float itof = 0.f; + float otof = 0.f; + float covm[15] = {0.f}; + float eigval[5] = {0.f}; + float eigvec[5][5] = {{0.f}}; + float eiginv[5][5] = {{0.f}}; + + void print() const; +}; + +/** + * @brief Binning map + */ +struct map_t { + int nbins = 1; + float min = 0.f; + float max = 1.e6f; + bool log = false; + + float eval(int bin) const + { + float width = (max - min) / nbins; + float val = min + (bin + 0.5f) * width; + if (log) + return std::pow(10.f, val); + return val; + } + + float fracPositionWithinBin(float val) const; + int find(float val) const; + void print() const; +}; + +/** + * @brief LUT header + */ +struct lutHeader_t { + int version = LUTCOVM_VERSION; + int pdg = 0; + float mass = 0.f; + float field = 0.f; + map_t nchmap; + map_t radmap; + map_t etamap; + map_t ptmap; + + bool check_version() const; + void print() const; +}; + +/** + * @brief Flat LUT data container - single contiguous buffer + * Memory layout: [header][entry_0][entry_1]...[entry_N] + * + * All entries stored sequentially in a single allocation. + * Can be directly mapped from file or shared memory without copying. + */ +class FlatLutData +{ + public: + FlatLutData() = default; + + /** + * @brief Initialize from binning information + * Pre-allocates contiguous memory for all entries + */ + void initialize(const lutHeader_t& header); + + /** + * @brief Get LUT entry by bin indices + * O(1) access via linear index calculation + */ + lutEntry_t* getEntry(int nch_bin, int rad_bin, int eta_bin, int pt_bin); + const lutEntry_t* getEntry(int nch_bin, int rad_bin, int eta_bin, int pt_bin) const; + + /** + * @brief Get raw data buffer (for serialization/shared memory) + */ + uint8_t* data() { return mData.data(); } + const uint8_t* data() const { return mData.data(); } + + /** + * @brief Total size in bytes + */ + size_t bytes() const { return mData.size(); } + + /** + * @brief Construct from external buffer (e.g., shared memory or file mapping) + */ + static FlatLutData fromBuffer(const uint8_t* buffer, size_t size); + + /** + * @brief Reference-based access without copying + * Useful when data is already in shared memory + */ + static FlatLutData fromExternalBuffer(uint8_t* buffer, size_t size); + + const lutHeader_t& header() const + { + return *reinterpret_cast(mData.data()); + } + + lutHeader_t& header() + { + return *reinterpret_cast(mData.data()); + } + + private: + /** + * @brief Linear index calculation for entry access + */ + size_t getEntryOffset(int nch_bin, int rad_bin, int eta_bin, int pt_bin) const; + + std::vector mData; + + // Cache dimensions for quick access + int mNchBins = 0; + int mRadBins = 0; + int mEtaBins = 0; + int mPtBins = 0; +}; + +} // namespace o2::delphes + +#endif // ALICE3_CORE_FLATLUTENTRY_H_ From ee2cc965cb771bbb0979359a2f14ee229a0d7764 Mon Sep 17 00:00:00 2001 From: Anton Alkin Date: Thu, 9 Apr 2026 11:11:43 +0200 Subject: [PATCH 02/11] update FlatLutEntry --- ALICE3/Core/FlatLutEntry.cxx | 99 +++++++++++++++++++++++------------- ALICE3/Core/FlatLutEntry.h | 68 +++++++++++++++++-------- 2 files changed, 110 insertions(+), 57 deletions(-) diff --git a/ALICE3/Core/FlatLutEntry.cxx b/ALICE3/Core/FlatLutEntry.cxx index f1b77d3c056..3234c6e41ce 100644 --- a/ALICE3/Core/FlatLutEntry.cxx +++ b/ALICE3/Core/FlatLutEntry.cxx @@ -12,6 +12,7 @@ #include "FlatLutEntry.h" #include +#include #include @@ -89,9 +90,9 @@ void FlatLutData::initialize(const lutHeader_t& header) size_t totalSize = headerSize + entriesSize; mData.resize(totalSize); - // Write header at the beginning std::memcpy(mData.data(), &header, headerSize); + updateRef(); } size_t FlatLutData::getEntryOffset(int nch_bin, int rad_bin, int eta_bin, int pt_bin) const @@ -105,69 +106,95 @@ size_t FlatLutData::getEntryOffset(int nch_bin, int rad_bin, int eta_bin, int pt return headerSize + linearIdx * sizeof(lutEntry_t); } -lutEntry_t* FlatLutData::getEntry(int nch_bin, int rad_bin, int eta_bin, int pt_bin) +const lutEntry_t* FlatLutData::getEntryRef(int nch_bin, int rad_bin, int eta_bin, int pt_bin) const { size_t offset = getEntryOffset(nch_bin, rad_bin, eta_bin, pt_bin); - return reinterpret_cast(mData.data() + offset); + return reinterpret_cast(mDataRef.data() + offset); } -const lutEntry_t* FlatLutData::getEntry(int nch_bin, int rad_bin, int eta_bin, int pt_bin) const +lutEntry_t* FlatLutData::getEntry(int nch_bin, int rad_bin, int eta_bin, int pt_bin) { size_t offset = getEntryOffset(nch_bin, rad_bin, eta_bin, pt_bin); - return reinterpret_cast(mData.data() + offset); + return reinterpret_cast(mData.data() + offset); } -FlatLutData FlatLutData::fromBuffer(const uint8_t* buffer, size_t size) +const lutHeader_t& FlatLutData::getHeaderRef() const { - FlatLutData data; - // Validate buffer - if (size < sizeof(lutHeader_t)) { - LOG(fatal) << "Buffer too small for LUT header"; - } + return *reinterpret_cast(mDataRef.data()); +} - const auto* header = reinterpret_cast(buffer); - data.mNchBins = header->nchmap.nbins; - data.mRadBins = header->radmap.nbins; - data.mEtaBins = header->etamap.nbins; - data.mPtBins = header->ptmap.nbins; +lutHeader_t& FlatLutData::getHeader() +{ + return *reinterpret_cast(mData.data()); +} - size_t expectedSize = sizeof(lutHeader_t) + static_cast(data.mNchBins) * data.mRadBins * data.mEtaBins * data.mPtBins * sizeof(lutEntry_t); +void FlatLutData::updateRef() +{ + mDataRef = std::span{mData.data(), mData.size()}; +} - if (size < expectedSize) { - LOG(fatal) << "Buffer size mismatch: expected " << expectedSize << ", got " << size; - } +void FlatLutData::cacheDimensions() +{ + auto const& header = getHeaderRef(); + mNchBins = header.nchmap.nbins; + mRadBins = header.radmap.nbins; + mEtaBins = header.etamap.nbins; + mPtBins = header.ptmap.nbins; +} - // Copy buffer - data.mData.resize(size); - std::memcpy(data.mData.data(), buffer, size); +void FlatLutData::adopt(const uint8_t* buffer, size_t size) +{ + mData.resize(size); + std::memcpy(mData.data(), buffer, size); + updateRef(); + cacheDimensions(); +} - return data; +void FlatLutData::view(const uint8_t* buffer, size_t size) +{ + mData.clear(); + mDataRef = std::span{buffer, size}; + cacheDimensions(); } -FlatLutData FlatLutData::fromExternalBuffer(uint8_t* buffer, size_t size) +void FlatLutData::validateBuffer(const uint8_t* buffer, size_t size) { - FlatLutData data; // Validate buffer if (size < sizeof(lutHeader_t)) { - LOG(fatal) << "Buffer too small for LUT header"; + throw framework::runtime_error_f("Buffer too small for LUT header: expected at least %zu, got %zu", sizeof(lutHeader_t), size); } const auto* header = reinterpret_cast(buffer); - data.mNchBins = header->nchmap.nbins; - data.mRadBins = header->radmap.nbins; - data.mEtaBins = header->etamap.nbins; - data.mPtBins = header->ptmap.nbins; + auto mNchBins = header->nchmap.nbins; + auto mRadBins = header->radmap.nbins; + auto mEtaBins = header->etamap.nbins; + auto mPtBins = header->ptmap.nbins; - size_t expectedSize = sizeof(lutHeader_t) + static_cast(data.mNchBins) * data.mRadBins * data.mEtaBins * data.mPtBins * sizeof(lutEntry_t); + size_t expectedSize = sizeof(lutHeader_t) + static_cast(mNchBins) * mRadBins * mEtaBins * mPtBins * sizeof(lutEntry_t); if (size < expectedSize) { - LOG(fatal) << "Buffer size mismatch: expected " << expectedSize << " got " << size; + throw framework::runtime_error_f("Buffer size mismatch: expected %zu, got %zu", expectedSize, size); } +} - // Store reference to external buffer (no copy) - // WARNING: Caller must ensure buffer lifetime exceeds FlatLutData usage - data.mData.assign(buffer, buffer + size); +FlatLutData FlatLutData::AdoptFromBuffer(const uint8_t* buffer, size_t size) +{ + validateBuffer(buffer, size); + FlatLutData data; + + // Copy buffer + data.adopt(buffer, size); + return data; +} +FlatLutData FlatLutData::ViewFromBuffer(const uint8_t* buffer, size_t size) +{ + validateBuffer(buffer, size); + FlatLutData data; + + // Store reference to external buffer + // WARNING: Caller must ensure buffer lifetime exceeds FlatLutData usage + data.view(buffer, size); return data; } diff --git a/ALICE3/Core/FlatLutEntry.h b/ALICE3/Core/FlatLutEntry.h index 18a2db63680..b0eb15350d8 100644 --- a/ALICE3/Core/FlatLutEntry.h +++ b/ALICE3/Core/FlatLutEntry.h @@ -14,6 +14,7 @@ #include #include +#include #include #define LUTCOVM_VERSION 20260408 @@ -100,43 +101,45 @@ class FlatLutData void initialize(const lutHeader_t& header); /** - * @brief Get LUT entry by bin indices - * O(1) access via linear index calculation + * @brief Get LUT entry by bin indices (view) + */ + const lutEntry_t* getEntryRef(int nch_bin, int rad_bin, int eta_bin, int pt_bin) const; + + /** + * @brief Get LUT entry by bin indices (owned) */ lutEntry_t* getEntry(int nch_bin, int rad_bin, int eta_bin, int pt_bin); - const lutEntry_t* getEntry(int nch_bin, int rad_bin, int eta_bin, int pt_bin) const; /** - * @brief Get raw data buffer (for serialization/shared memory) + * @brief Get LUT header (view) */ - uint8_t* data() { return mData.data(); } - const uint8_t* data() const { return mData.data(); } + const lutHeader_t& getHeaderRef() const; /** - * @brief Total size in bytes + * @brief Get LUT header (owned) */ - size_t bytes() const { return mData.size(); } + lutHeader_t& getHeader(); /** - * @brief Construct from external buffer (e.g., shared memory or file mapping) + * @brief Get raw data buffer */ - static FlatLutData fromBuffer(const uint8_t* buffer, size_t size); + uint8_t* data() { return mData.data(); } // owned + const uint8_t* data() const { return mDataRef.data(); } //view /** - * @brief Reference-based access without copying - * Useful when data is already in shared memory + * @brief Total size in bytes */ - static FlatLutData fromExternalBuffer(uint8_t* buffer, size_t size); + size_t bytes() const { return mDataRef.size(); } - const lutHeader_t& header() const - { - return *reinterpret_cast(mData.data()); - } + /** + * @brief Construct a new FlatLutData from external buffer as a copy + */ + static FlatLutData AdoptFromBuffer(const uint8_t* buffer, size_t size); - lutHeader_t& header() - { - return *reinterpret_cast(mData.data()); - } + /** + * @brief Construct a new FlatLutData from external buffer as a view + */ + static FlatLutData ViewFromBuffer(const uint8_t* buffer, size_t size); private: /** @@ -144,7 +147,30 @@ class FlatLutData */ size_t getEntryOffset(int nch_bin, int rad_bin, int eta_bin, int pt_bin) const; + /** + * @brief Update dimensions from the current header + */ + void cacheDimensions(); + + /** + * @brief Adopt a buffer by copying + */ + void adopt(const uint8_t* buffer, size_t size); + /** + * @brief Adopt a buffer as a view + */ + void view(const uint8_t* buffer, size_t size); + /** + * @brief Update mDataRef from mData + */ + void updateRef(); + /** + * @brief Validate external buffer + */ + static void validateBuffer(const uint8_t* buffer, size_t size); + std::vector mData; + std::span mDataRef; // Cache dimensions for quick access int mNchBins = 0; From 5e410506a622bffae45ae4aa670b9a39bf832155 Mon Sep 17 00:00:00 2001 From: Anton Alkin Date: Thu, 9 Apr 2026 11:12:18 +0200 Subject: [PATCH 03/11] start adding FlatTrackSmearer --- ALICE3/Core/FlatTrackSmearer.cxx | 66 +++++++++++++++++++++++++++ ALICE3/Core/FlatTrackSmearer.h | 76 ++++++++++++++++++++++++++++++++ 2 files changed, 142 insertions(+) create mode 100644 ALICE3/Core/FlatTrackSmearer.cxx create mode 100644 ALICE3/Core/FlatTrackSmearer.h diff --git a/ALICE3/Core/FlatTrackSmearer.cxx b/ALICE3/Core/FlatTrackSmearer.cxx new file mode 100644 index 00000000000..199ae940638 --- /dev/null +++ b/ALICE3/Core/FlatTrackSmearer.cxx @@ -0,0 +1,66 @@ +// 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. + +#include "FlatTrackSmearer.h" + +namespace o2::delphes { +int TrackSmearer::getIndexPDG(int pdg) const +{ + switch (std::abs(pdg)) { + case 11: + return 0; // Electron + case 13: + return 1; // Muon + case 211: + return 2; // Pion + case 321: + return 3; // Kaon + case 2212: + return 4; // Proton + case 1000010020: + return 5; // Deuteron + case 1000010030: + return 6; // Triton + case 1000020030: + return 7; // Helium3 + case 1000020040: + return 8; // Alphas + default: + return 2; // Default: pion + } +} + +const char* TrackSmearer::getParticleName(int pdg) const +{ + switch (std::abs(pdg)) { + case 11: + return "electron"; + case 13: + return "muon"; + case 211: + return "pion"; + case 321: + return "kaon"; + case 2212: + return "proton"; + case 1000010020: + return "deuteron"; + case 1000010030: + return "triton"; + case 1000020030: + return "helium3"; + case 1000020040: + return "alpha"; + default: + return "pion"; // Default: pion + } +} +} // namespace o2::delphes diff --git a/ALICE3/Core/FlatTrackSmearer.h b/ALICE3/Core/FlatTrackSmearer.h new file mode 100644 index 00000000000..078a0bda6db --- /dev/null +++ b/ALICE3/Core/FlatTrackSmearer.h @@ -0,0 +1,76 @@ +// 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. + +#ifndef ALICE3_CORE_FLATTRACKSMEARER_H_ +#define ALICE3_CORE_FLATTRACKSMEARER_H_ + +#include "FlatLutEntry.h" + +#include +#include + +namespace o2::delphes { +/** + * @brief Track smearing with flat LUT backend + */ +using O2Track = o2::track::TrackParCov; +class TrackSmearer +{ + public: + TrackSmearer() = default; + ~TrackSmearer() = default; + + /** LUT methods **/ + bool loadTable(int pdg, const char* filename, bool forceReload = false); + bool hasTable(int pdg) const; + + void useEfficiency(bool val) { mUseEfficiency = val; } + void interpolateEfficiency(bool val) { mInterpolateEfficiency = val; } + void skipUnreconstructed(bool val) { mSkipUnreconstructed = val; } + void setWhatEfficiency(int val) { mWhatEfficiency = val; } + + const lutHeader_t* getLUTHeader(int pdg) const; + const lutEntry_t* getLUTEntry(int pdg, float nch, float radius, float eta, float pt, float& interpolatedEff) const; + + bool smearTrack(O2Track& o2track, const lutEntry_t* lutEntry, float interpolatedEff); + bool smearTrack(O2Track& o2track, int pdg, float nch); + + double getPtRes(int pdg, float nch, float eta, float pt) const; + double getEtaRes(int pdg, float nch, float eta, float pt) const; + double getAbsPtRes(int pdg, float nch, float eta, float pt) const; + double getAbsEtaRes(int pdg, float nch, float eta, float pt) const; + double getEfficiency(int pdg, float nch, float eta, float pt) const; + + int getIndexPDG(int pdg) const; + + const char* getParticleName(int pdg) const; + + void setdNdEta(float val) { mdNdEta = val; } + void setCcdbManager(o2::ccdb::BasicCCDBManager* mgr) { mCcdbManager = mgr; } + + protected: + static constexpr unsigned int nLUTs = 9; // Number of LUT available + lutHeader_t const* mHeaders[nLUTs]; // header references for quick access + FlatLutData mLUTData[nLUTs]; // NEW: Flat data storage + + bool mUseEfficiency = true; + bool mInterpolateEfficiency = false; + bool mSkipUnreconstructed = true; // don't smear tracks that are not reco'ed + int mWhatEfficiency = 1; + float mdNdEta = 1600.f; + + private: + o2::ccdb::BasicCCDBManager* mCcdbManager = nullptr; +}; + +} + +#endif // ALICE3_CORE_FLATTRACKSMEARER_H_ From 0db84ccc83d707812597d8fec5553f1bd67c5c2c Mon Sep 17 00:00:00 2001 From: Anton Alkin Date: Thu, 9 Apr 2026 11:12:44 +0200 Subject: [PATCH 04/11] fixup! start adding FlatTrackSmearer --- ALICE3/Core/FlatTrackSmearer.cxx | 3 ++- ALICE3/Core/FlatTrackSmearer.h | 9 +++++---- 2 files changed, 7 insertions(+), 5 deletions(-) diff --git a/ALICE3/Core/FlatTrackSmearer.cxx b/ALICE3/Core/FlatTrackSmearer.cxx index 199ae940638..73479e1fb84 100644 --- a/ALICE3/Core/FlatTrackSmearer.cxx +++ b/ALICE3/Core/FlatTrackSmearer.cxx @@ -11,7 +11,8 @@ #include "FlatTrackSmearer.h" -namespace o2::delphes { +namespace o2::delphes +{ int TrackSmearer::getIndexPDG(int pdg) const { switch (std::abs(pdg)) { diff --git a/ALICE3/Core/FlatTrackSmearer.h b/ALICE3/Core/FlatTrackSmearer.h index 078a0bda6db..26122f7730c 100644 --- a/ALICE3/Core/FlatTrackSmearer.h +++ b/ALICE3/Core/FlatTrackSmearer.h @@ -17,7 +17,8 @@ #include #include -namespace o2::delphes { +namespace o2::delphes +{ /** * @brief Track smearing with flat LUT backend */ @@ -58,8 +59,8 @@ class TrackSmearer protected: static constexpr unsigned int nLUTs = 9; // Number of LUT available - lutHeader_t const* mHeaders[nLUTs]; // header references for quick access - FlatLutData mLUTData[nLUTs]; // NEW: Flat data storage + lutHeader_t const* mHeaders[nLUTs]; // header references for quick access + FlatLutData mLUTData[nLUTs]; // NEW: Flat data storage bool mUseEfficiency = true; bool mInterpolateEfficiency = false; @@ -71,6 +72,6 @@ class TrackSmearer o2::ccdb::BasicCCDBManager* mCcdbManager = nullptr; }; -} +} // namespace o2::delphes #endif // ALICE3_CORE_FLATTRACKSMEARER_H_ From df4f545d7b491e0aca5f989ce02e18130c171a77 Mon Sep 17 00:00:00 2001 From: Anton Alkin Date: Thu, 9 Apr 2026 13:26:17 +0200 Subject: [PATCH 05/11] make file reading make sense --- ALICE3/Core/FlatLutEntry.cxx | 58 +++++++++++++++++++++++++++++ ALICE3/Core/FlatLutEntry.h | 21 +++++++++++ ALICE3/Core/FlatTrackSmearer.cxx | 63 ++++++++++++++++++++++++++++++++ ALICE3/Core/FlatTrackSmearer.h | 3 +- 4 files changed, 143 insertions(+), 2 deletions(-) diff --git a/ALICE3/Core/FlatLutEntry.cxx b/ALICE3/Core/FlatLutEntry.cxx index 3234c6e41ce..f7d3f8cd427 100644 --- a/ALICE3/Core/FlatLutEntry.cxx +++ b/ALICE3/Core/FlatLutEntry.cxx @@ -142,6 +142,14 @@ void FlatLutData::cacheDimensions() mPtBins = header.ptmap.nbins; } +void FlatLutData::resetDimensions() +{ + mNchBins = 0; + mRadBins = 0; + mEtaBins = 0; + mPtBins = 0; +} + void FlatLutData::adopt(const uint8_t* buffer, size_t size) { mData.resize(size); @@ -198,4 +206,54 @@ FlatLutData FlatLutData::ViewFromBuffer(const uint8_t* buffer, size_t size) return data; } +bool FlatLutData::isLoaded() const +{ + return ((!mData.empty()) || (!mDataRef.empty())); +} + +FlatLutData FlatLutData::loadFromFile(const char* filename) +{ + std::ifstream lutFile(filename, std::ifstream::binary); + if (!lutFile.is_open()) { + throw framework::runtime_error_f("Cannot open LUT file: %s", filename); + } + + // Read header first + lutHeader_t tempHeader; + lutFile.read(reinterpret_cast(&tempHeader), sizeof(lutHeader_t)); + if (lutFile.gcount() != static_cast(sizeof(lutHeader_t))) { + throw framework::runtime_error_f("Failed to read LUT header from %s", filename); + } + + if (!tempHeader.check_version()) { + throw framework::runtime_error_f("LUT header version mismatch: expected %d, got %d", LUTCOVM_VERSION, tempHeader.version); + } + + FlatLutData data; + + // Initialize flat data structure + data.initialize(tempHeader); + + // Read all entries sequentially into flat buffer + size_t headerSize = sizeof(lutHeader_t); + size_t numEntries = static_cast(data.mNchBins) * data.mRadBins * data.mEtaBins * data.mPtBins; + size_t entriesSize = numEntries * sizeof(lutEntry_t); + + lutFile.read(reinterpret_cast(data.data() + headerSize), entriesSize); + if (lutFile.gcount() != static_cast(entriesSize)) { + throw framework::runtime_error_f("Failed to read LUT entries from %s: expected %zu bytes, got %zu", filename, entriesSize, static_cast(lutFile.gcount())); + } + + lutFile.close(); + LOGF(info, "Successfully loaded LUT from %s: %zu entries", filename, numEntries); + return data; +} + +void FlatLutData::reset() +{ + mData.clear(); + updateRef(); + resetDimensions(); +} + } // namespace o2::delphes diff --git a/ALICE3/Core/FlatLutEntry.h b/ALICE3/Core/FlatLutEntry.h index b0eb15350d8..5bdbda46e1f 100644 --- a/ALICE3/Core/FlatLutEntry.h +++ b/ALICE3/Core/FlatLutEntry.h @@ -93,6 +93,8 @@ class FlatLutData { public: FlatLutData() = default; + FlatLutData(FlatLutData&&) = default; + FlatLutData& operator=(FlatLutData&&) = default; /** * @brief Initialize from binning information @@ -141,6 +143,21 @@ class FlatLutData */ static FlatLutData ViewFromBuffer(const uint8_t* buffer, size_t size); + /** + * @brief Construct a new FlatLutData from a file + */ + static FlatLutData loadFromFile(const char* filename); + + /** + * @brief Check if the LUT is loaded + */ + bool isLoaded() const; + + /** + * @brief Reset LUT to empty + */ + void reset(); + private: /** * @brief Linear index calculation for entry access @@ -151,6 +168,10 @@ class FlatLutData * @brief Update dimensions from the current header */ void cacheDimensions(); + /** + * @brief Reset dimensions to 0 + */ + void resetDimensions(); /** * @brief Adopt a buffer by copying diff --git a/ALICE3/Core/FlatTrackSmearer.cxx b/ALICE3/Core/FlatTrackSmearer.cxx index 73479e1fb84..396b61285b5 100644 --- a/ALICE3/Core/FlatTrackSmearer.cxx +++ b/ALICE3/Core/FlatTrackSmearer.cxx @@ -11,6 +11,10 @@ #include "FlatTrackSmearer.h" +#include "ALICE3/Core/GeometryContainer.h" + +#include + namespace o2::delphes { int TrackSmearer::getIndexPDG(int pdg) const @@ -64,4 +68,63 @@ const char* TrackSmearer::getParticleName(int pdg) const return "pion"; // Default: pion } } + +bool TrackSmearer::loadTable(int pdg, const char* filename, bool forceReload) +{ + if (!filename || filename[0] == '\0') { + LOGF(info, "No LUT file provided for PDG %d. Skipping load.", pdg); + return false; + } + + const auto ipdg = getIndexPDG(pdg); + LOGF(info, "Loading %s LUT file: '%s'", getParticleName(pdg), filename); + + if (mLUTData[ipdg].isLoaded() && !forceReload) { + LOGF(info, "LUT table for PDG %d already loaded (index %d)", pdg, ipdg); + return false; + } + + const std::string localFilename = o2::fastsim::GeometryEntry::accessFile(filename, "./.ALICE3/LUTs/", mCcdbManager, 10); + + try { + mLUTData[ipdg] = FlatLutData::loadFromFile(localFilename.c_str()); + } catch (framework::RuntimeErrorRef ref) { + LOGF(error, "%s", framework::error_from_ref(ref).what); + return false; + } + + // Validate header + const auto& header = mLUTData[ipdg].getHeaderRef(); + + bool specialPdgCase = false; + switch (pdg) { + case o2::constants::physics::kAlpha: + // Special case: Allow Alpha particles to use He3 LUT + specialPdgCase = (header.pdg == o2::constants::physics::kHelium3); + if (specialPdgCase) { + LOGF(info, "Alpha particles (PDG %d) will use He3 LUT data (PDG %d)", pdg, header.pdg); + } + break; + default: + break; + } + + if (header.pdg != pdg && !specialPdgCase) { + LOGF(error, "LUT header PDG mismatch: expected %d, got %d", pdg, header.pdg); + mLUTData[ipdg].reset(); + return false; + } + + LOGF(info, "Successfully read LUT for PDG %d: %s", pdg, localFilename.c_str()); + header.print(); + + return true; +} + +bool TrackSmearer::hasTable(int pdg) const +{ + const int ipdg = getIndexPDG(pdg); + return mLUTData[ipdg].isLoaded(); +} + } // namespace o2::delphes diff --git a/ALICE3/Core/FlatTrackSmearer.h b/ALICE3/Core/FlatTrackSmearer.h index 26122f7730c..dac83005083 100644 --- a/ALICE3/Core/FlatTrackSmearer.h +++ b/ALICE3/Core/FlatTrackSmearer.h @@ -27,7 +27,6 @@ class TrackSmearer { public: TrackSmearer() = default; - ~TrackSmearer() = default; /** LUT methods **/ bool loadTable(int pdg, const char* filename, bool forceReload = false); @@ -60,7 +59,7 @@ class TrackSmearer protected: static constexpr unsigned int nLUTs = 9; // Number of LUT available lutHeader_t const* mHeaders[nLUTs]; // header references for quick access - FlatLutData mLUTData[nLUTs]; // NEW: Flat data storage + FlatLutData mLUTData[nLUTs]; // Flat data storage bool mUseEfficiency = true; bool mInterpolateEfficiency = false; From 63df43bc34458049cc2e174fcf269601cd62a87c Mon Sep 17 00:00:00 2001 From: Anton Alkin Date: Fri, 10 Apr 2026 10:52:19 +0200 Subject: [PATCH 06/11] rearrange code; preview header before reading full file or copying buffer --- ALICE3/Core/FlatLutEntry.cxx | 55 ++++++++++-------- ALICE3/Core/FlatLutEntry.h | 12 +++- ALICE3/Core/FlatTrackSmearer.cxx | 97 +++++++++++++++++++++++++------- ALICE3/Core/FlatTrackSmearer.h | 4 ++ 4 files changed, 123 insertions(+), 45 deletions(-) diff --git a/ALICE3/Core/FlatLutEntry.cxx b/ALICE3/Core/FlatLutEntry.cxx index f7d3f8cd427..2c97ed26838 100644 --- a/ALICE3/Core/FlatLutEntry.cxx +++ b/ALICE3/Core/FlatLutEntry.cxx @@ -167,16 +167,14 @@ void FlatLutData::view(const uint8_t* buffer, size_t size) void FlatLutData::validateBuffer(const uint8_t* buffer, size_t size) { - // Validate buffer - if (size < sizeof(lutHeader_t)) { - throw framework::runtime_error_f("Buffer too small for LUT header: expected at least %zu, got %zu", sizeof(lutHeader_t), size); + auto header = PreviewHeader(buffer, size); + if (!header.check_version()) { + throw framework::runtime_error_f("LUT header version mismatch: expected %d, got %d", LUTCOVM_VERSION, header.version); } - - const auto* header = reinterpret_cast(buffer); - auto mNchBins = header->nchmap.nbins; - auto mRadBins = header->radmap.nbins; - auto mEtaBins = header->etamap.nbins; - auto mPtBins = header->ptmap.nbins; + auto mNchBins = header.nchmap.nbins; + auto mRadBins = header.radmap.nbins; + auto mEtaBins = header.etamap.nbins; + auto mPtBins = header.ptmap.nbins; size_t expectedSize = sizeof(lutHeader_t) + static_cast(mNchBins) * mRadBins * mEtaBins * mPtBins * sizeof(lutEntry_t); @@ -185,6 +183,18 @@ void FlatLutData::validateBuffer(const uint8_t* buffer, size_t size) } } +lutHeader_t FlatLutData::PreviewHeader(const uint8_t* buffer, size_t size) +{ + if (size < sizeof(lutHeader_t)) { + throw framework::runtime_error_f("Buffer too small for LUT header: expected at least %zu, got %zu", sizeof(lutHeader_t), size); + } + const auto* header = reinterpret_cast(buffer); + if (!header->check_version()) { + throw framework::runtime_error_f("LUT header version mismatch: expected %d, got %d", LUTCOVM_VERSION, header->version); + } + return *header; +} + FlatLutData FlatLutData::AdoptFromBuffer(const uint8_t* buffer, size_t size) { validateBuffer(buffer, size); @@ -211,23 +221,23 @@ bool FlatLutData::isLoaded() const return ((!mData.empty()) || (!mDataRef.empty())); } -FlatLutData FlatLutData::loadFromFile(const char* filename) +lutHeader_t FlatLutData::PreviewHeader(std::ifstream& file, const char* filename) { - std::ifstream lutFile(filename, std::ifstream::binary); - if (!lutFile.is_open()) { - throw framework::runtime_error_f("Cannot open LUT file: %s", filename); - } - - // Read header first lutHeader_t tempHeader; - lutFile.read(reinterpret_cast(&tempHeader), sizeof(lutHeader_t)); - if (lutFile.gcount() != static_cast(sizeof(lutHeader_t))) { + file.read(reinterpret_cast(&tempHeader), sizeof(lutHeader_t)); + if (file.gcount() != static_cast(sizeof(lutHeader_t))) { throw framework::runtime_error_f("Failed to read LUT header from %s", filename); } - if (!tempHeader.check_version()) { throw framework::runtime_error_f("LUT header version mismatch: expected %d, got %d", LUTCOVM_VERSION, tempHeader.version); } + return tempHeader; +} + +FlatLutData FlatLutData::loadFromFile(std::ifstream& file, const char* filename) +{ + // Read header first + lutHeader_t tempHeader = PreviewHeader(file, filename); FlatLutData data; @@ -239,12 +249,11 @@ FlatLutData FlatLutData::loadFromFile(const char* filename) size_t numEntries = static_cast(data.mNchBins) * data.mRadBins * data.mEtaBins * data.mPtBins; size_t entriesSize = numEntries * sizeof(lutEntry_t); - lutFile.read(reinterpret_cast(data.data() + headerSize), entriesSize); - if (lutFile.gcount() != static_cast(entriesSize)) { - throw framework::runtime_error_f("Failed to read LUT entries from %s: expected %zu bytes, got %zu", filename, entriesSize, static_cast(lutFile.gcount())); + file.read(reinterpret_cast(data.data() + headerSize), entriesSize); + if (file.gcount() != static_cast(entriesSize)) { + throw framework::runtime_error_f("Failed to read LUT entries from %s: expected %zu bytes, got %zu", filename, entriesSize, static_cast(file.gcount())); } - lutFile.close(); LOGF(info, "Successfully loaded LUT from %s: %zu entries", filename, numEntries); return data; } diff --git a/ALICE3/Core/FlatLutEntry.h b/ALICE3/Core/FlatLutEntry.h index 5bdbda46e1f..c2508efa565 100644 --- a/ALICE3/Core/FlatLutEntry.h +++ b/ALICE3/Core/FlatLutEntry.h @@ -146,7 +146,17 @@ class FlatLutData /** * @brief Construct a new FlatLutData from a file */ - static FlatLutData loadFromFile(const char* filename); + static FlatLutData loadFromFile(std::ifstream& file, const char* filename); + + /** + * @brief Preview buffer header for version and other compatibility checks + */ + static lutHeader_t PreviewHeader(const uint8_t* buffer, size_t size); + + /** + * @brief Preview file-stored header for version and other compatibility checks + */ + static lutHeader_t PreviewHeader(std::ifstream& file, const char* filename); /** * @brief Check if the LUT is loaded diff --git a/ALICE3/Core/FlatTrackSmearer.cxx b/ALICE3/Core/FlatTrackSmearer.cxx index 396b61285b5..655d025155a 100644 --- a/ALICE3/Core/FlatTrackSmearer.cxx +++ b/ALICE3/Core/FlatTrackSmearer.cxx @@ -77,25 +77,95 @@ bool TrackSmearer::loadTable(int pdg, const char* filename, bool forceReload) } const auto ipdg = getIndexPDG(pdg); - LOGF(info, "Loading %s LUT file: '%s'", getParticleName(pdg), filename); - if (mLUTData[ipdg].isLoaded() && !forceReload) { LOGF(info, "LUT table for PDG %d already loaded (index %d)", pdg, ipdg); return false; } + std::ifstream lutFile(filename, std::ifstream::binary); + if (!lutFile.is_open()) { + throw framework::runtime_error_f("Cannot open LUT file: %s", filename); + } + + LOGF(info, "Loading %s LUT file: '%s'", getParticleName(pdg), filename); const std::string localFilename = o2::fastsim::GeometryEntry::accessFile(filename, "./.ALICE3/LUTs/", mCcdbManager, 10); try { - mLUTData[ipdg] = FlatLutData::loadFromFile(localFilename.c_str()); + auto header = FlatLutData::PreviewHeader(lutFile, localFilename.c_str()); + // Validate header + if (header.pdg != pdg && !checkSpecialCase(pdg, header)) { + LOGF(error, "LUT header PDG mismatch: expected %d, got %d; not loading", pdg, header.pdg); + return false; + } + mLUTData[ipdg] = FlatLutData::loadFromFile(lutFile, localFilename.c_str()); } catch (framework::RuntimeErrorRef ref) { LOGF(error, "%s", framework::error_from_ref(ref).what); return false; } + mHeaders[ipdg] = &mLUTData[ipdg].getHeaderRef(); - // Validate header - const auto& header = mLUTData[ipdg].getHeaderRef(); + LOGF(info, "Successfully read LUT for PDG %d: %s", pdg, localFilename.c_str()); + mHeaders[ipdg]->print(); + return true; +} +bool TrackSmearer::adoptTable(int pdg, const uint8_t* buffer, size_t size, bool forceReload) +{ + const auto ipdg = getIndexPDG(pdg); + if (mLUTData[ipdg].isLoaded() && !forceReload) { + LOGF(info, "LUT table for PDG %d already loaded (index %d)", pdg, ipdg); + return false; + } + try { + auto header = FlatLutData::PreviewHeader(buffer, size); + if (header.pdg != pdg && !checkSpecialCase(pdg, header)) { + LOGF(error, "LUT header PDG mismatch: expected %d, got %d", pdg, header.pdg); + return false; + } + mLUTData[ipdg] = FlatLutData::AdoptFromBuffer(buffer, size); + } catch (framework::RuntimeErrorRef ref) { + LOGF(error, "%s", framework::error_from_ref(ref).what); + } + mHeaders[ipdg] = &mLUTData[ipdg].getHeaderRef(); + + LOGF(info, "Successfully adopted LUT for PDG %d", pdg); + mHeaders[ipdg]->print(); + return true; +} + +bool TrackSmearer::viewTable(int pdg, const uint8_t* buffer, size_t size, bool forceReload) +{ + const auto ipdg = getIndexPDG(pdg); + if (mLUTData[ipdg].isLoaded() && !forceReload) { + LOGF(info, "LUT table for PDG %d already loaded (index %d)", pdg, ipdg); + return false; + } + try { + auto header = FlatLutData::PreviewHeader(buffer, size); + if (header.pdg != pdg && !checkSpecialCase(pdg, header)) { + LOGF(error, "LUT header PDG mismatch: expected %d, got %d", pdg, header.pdg); + return false; + } + mLUTData[ipdg] = FlatLutData::ViewFromBuffer(buffer, size); + } catch (framework::RuntimeErrorRef ref) { + LOGF(error, "%s", framework::error_from_ref(ref).what); + } + mHeaders[ipdg] = &mLUTData[ipdg].getHeaderRef(); + + LOGF(info, "Successfully adopted LUT for PDG %d", pdg); + mHeaders[ipdg]->print(); + return true; +} + +bool TrackSmearer::hasTable(int pdg) const +{ + const int ipdg = getIndexPDG(pdg); + return mLUTData[ipdg].isLoaded(); +} + +bool TrackSmearer::checkSpecialCase(int pdg, lutHeader_t const& header) +{ + // Validate header bool specialPdgCase = false; switch (pdg) { case o2::constants::physics::kAlpha: @@ -108,23 +178,8 @@ bool TrackSmearer::loadTable(int pdg, const char* filename, bool forceReload) default: break; } - - if (header.pdg != pdg && !specialPdgCase) { - LOGF(error, "LUT header PDG mismatch: expected %d, got %d", pdg, header.pdg); - mLUTData[ipdg].reset(); - return false; - } - - LOGF(info, "Successfully read LUT for PDG %d: %s", pdg, localFilename.c_str()); - header.print(); - - return true; + return specialPdgCase; } -bool TrackSmearer::hasTable(int pdg) const -{ - const int ipdg = getIndexPDG(pdg); - return mLUTData[ipdg].isLoaded(); -} } // namespace o2::delphes diff --git a/ALICE3/Core/FlatTrackSmearer.h b/ALICE3/Core/FlatTrackSmearer.h index dac83005083..2d4bb8428b3 100644 --- a/ALICE3/Core/FlatTrackSmearer.h +++ b/ALICE3/Core/FlatTrackSmearer.h @@ -30,6 +30,8 @@ class TrackSmearer /** LUT methods **/ bool loadTable(int pdg, const char* filename, bool forceReload = false); + bool adoptTable(int pdg, const uint8_t* buffer, size_t size, bool forceReload = false); + bool viewTable(int pdg, const uint8_t* buffer, size_t size, bool forceReload = false); bool hasTable(int pdg) const; void useEfficiency(bool val) { mUseEfficiency = val; } @@ -69,6 +71,8 @@ class TrackSmearer private: o2::ccdb::BasicCCDBManager* mCcdbManager = nullptr; + + bool checkSpecialCase(int pdg, lutHeader_t const& header); }; } // namespace o2::delphes From fb8d91bed84f735f233e68efb641b7eb0a100c2d Mon Sep 17 00:00:00 2001 From: Anton Alkin Date: Fri, 10 Apr 2026 12:23:28 +0200 Subject: [PATCH 07/11] add the rest of the methods --- ALICE3/Core/FlatTrackSmearer.cxx | 236 ++++++++++++++++++++++++++++++- ALICE3/Core/FlatTrackSmearer.h | 3 +- 2 files changed, 231 insertions(+), 8 deletions(-) diff --git a/ALICE3/Core/FlatTrackSmearer.cxx b/ALICE3/Core/FlatTrackSmearer.cxx index 655d025155a..67155d0a898 100644 --- a/ALICE3/Core/FlatTrackSmearer.cxx +++ b/ALICE3/Core/FlatTrackSmearer.cxx @@ -15,6 +15,8 @@ #include +#include + namespace o2::delphes { int TrackSmearer::getIndexPDG(int pdg) const @@ -69,6 +71,15 @@ const char* TrackSmearer::getParticleName(int pdg) const } } +void TrackSmearer::setWhatEfficiency(int val) +{ + // FIXME: this really should be an enum + if (val > 2) { + throw framework::runtime_error_f("getLUTEntry: unknown efficiency type %d", mWhatEfficiency); + } + mWhatEfficiency = val; +} + bool TrackSmearer::loadTable(int pdg, const char* filename, bool forceReload) { if (!filename || filename[0] == '\0') { @@ -102,10 +113,9 @@ bool TrackSmearer::loadTable(int pdg, const char* filename, bool forceReload) LOGF(error, "%s", framework::error_from_ref(ref).what); return false; } - mHeaders[ipdg] = &mLUTData[ipdg].getHeaderRef(); LOGF(info, "Successfully read LUT for PDG %d: %s", pdg, localFilename.c_str()); - mHeaders[ipdg]->print(); + mLUTData[ipdg].getHeaderRef().print(); return true; } @@ -126,10 +136,9 @@ bool TrackSmearer::adoptTable(int pdg, const uint8_t* buffer, size_t size, bool } catch (framework::RuntimeErrorRef ref) { LOGF(error, "%s", framework::error_from_ref(ref).what); } - mHeaders[ipdg] = &mLUTData[ipdg].getHeaderRef(); LOGF(info, "Successfully adopted LUT for PDG %d", pdg); - mHeaders[ipdg]->print(); + mLUTData[ipdg].getHeaderRef().print(); return true; } @@ -150,10 +159,9 @@ bool TrackSmearer::viewTable(int pdg, const uint8_t* buffer, size_t size, bool f } catch (framework::RuntimeErrorRef ref) { LOGF(error, "%s", framework::error_from_ref(ref).what); } - mHeaders[ipdg] = &mLUTData[ipdg].getHeaderRef(); LOGF(info, "Successfully adopted LUT for PDG %d", pdg); - mHeaders[ipdg]->print(); + mLUTData[ipdg].getHeaderRef().print(); return true; } @@ -181,5 +189,221 @@ bool TrackSmearer::checkSpecialCase(int pdg, lutHeader_t const& header) return specialPdgCase; } +const lutHeader_t* TrackSmearer::getLUTHeader(int pdg) const +{ + const int ipdg = getIndexPDG(pdg); + if (!mLUTData[ipdg].isLoaded()) { + return nullptr; + } + return &mLUTData[ipdg].getHeaderRef(); +} + +const lutEntry_t* TrackSmearer::getLUTEntry(const int pdg, const float nch, const float radius, const float eta, const float pt, float& interpolatedEff) const +{ + const int ipdg = getIndexPDG(pdg); + if (!mLUTData[ipdg].isLoaded()) { + return nullptr; + } + + const auto& header = mLUTData[ipdg].getHeaderRef(); + + auto inch = header.nchmap.find(nch); + auto irad = header.radmap.find(radius); + auto ieta = header.etamap.find(eta); + auto ipt = header.ptmap.find(pt); + + // Interpolate efficiency if requested + auto fraction = header.nchmap.fracPositionWithinBin(nch); + if (mInterpolateEfficiency) { + static constexpr float kFractionThreshold = 0.5f; + if (fraction > kFractionThreshold) { + switch (mWhatEfficiency) { + case 1: { + const auto* entry_curr = mLUTData[ipdg].getEntryRef(inch, irad, ieta, ipt); + if (inch < header.nchmap.nbins - 1) { + const auto* entry_next = mLUTData[ipdg].getEntryRef(inch + 1, irad, ieta, ipt); + interpolatedEff = (1.5f - fraction) * entry_curr->eff + (-0.5f + fraction) * entry_next->eff; + } else { + interpolatedEff = entry_curr->eff; + } + break; + } + case 2: { + const auto* entry_curr = mLUTData[ipdg].getEntryRef(inch, irad, ieta, ipt); + if (inch < header.nchmap.nbins - 1) { + const auto* entry_next = mLUTData[ipdg].getEntryRef(inch + 1, irad, ieta, ipt); + interpolatedEff = (1.5f - fraction) * entry_curr->eff2 + (-0.5f + fraction) * entry_next->eff2; + } else { + interpolatedEff = entry_curr->eff2; + } + break; + } + } + } else { + float comparisonValue = header.nchmap.log ? std::log10(nch) : nch; + switch (mWhatEfficiency) { + case 1: { + const auto* entry_curr = mLUTData[ipdg].getEntryRef(inch, irad, ieta, ipt); + if (inch > 0 && comparisonValue < header.nchmap.max) { + const auto* entry_prev = mLUTData[ipdg].getEntryRef(inch - 1, irad, ieta, ipt); + interpolatedEff = (0.5f + fraction) * entry_curr->eff + (0.5f - fraction) * entry_prev->eff; + } else { + interpolatedEff = entry_curr->eff; + } + break; + } + case 2: { + const auto* entry_curr = mLUTData[ipdg].getEntryRef(inch, irad, ieta, ipt); + if (inch > 0 && comparisonValue < header.nchmap.max) { + const auto* entry_prev = mLUTData[ipdg].getEntryRef(inch - 1, irad, ieta, ipt); + interpolatedEff = (0.5f + fraction) * entry_curr->eff2 + (0.5f - fraction) * entry_prev->eff2; + } else { + interpolatedEff = entry_curr->eff2; + } + break; + } + } + } + } else { + const auto* entry = mLUTData[ipdg].getEntryRef(inch, irad, ieta, ipt); + if (entry) { + switch (mWhatEfficiency) { + case 1: + interpolatedEff = entry->eff; + break; + case 2: + interpolatedEff = entry->eff2; + break; + } + } + } + + return mLUTData[ipdg].getEntryRef(inch, irad, ieta, ipt); +} + +bool TrackSmearer::smearTrack(O2Track& o2track, const lutEntry_t* lutEntry, float interpolatedEff) +{ + bool isReconstructed = true; + + // Generate efficiency + if (mUseEfficiency) { + auto eff = 0.f; + switch (mWhatEfficiency) { + case 1: + eff = lutEntry->eff; + break; + case 2: + eff = lutEntry->eff2; + break; + } + if (mInterpolateEfficiency) { + eff = interpolatedEff; + } + if (gRandom->Uniform() > eff) {//FIXME: use a fixed RNG instead of whatever ROOT has as a default + isReconstructed = false; + } + } + + // Return false already now in case not reco'ed + if (!isReconstructed && mSkipUnreconstructed) { + return false; + } + + // Transform params vector and smear + static constexpr int kParSize = 5; + double params[kParSize]; + for (int i = 0; i < kParSize; ++i) { + double val = 0.; + for (int j = 0; j < kParSize; ++j) { + val += lutEntry->eigvec[j][i] * o2track.getParam(j); + } + params[i] = gRandom->Gaus(val, std::sqrt(lutEntry->eigval[i])); + } + + // Transform back params vector + for (int i = 0; i < kParSize; ++i) { + double val = 0.; + for (int j = 0; j < kParSize; ++j) { + val += lutEntry->eiginv[j][i] * params[j]; + } + o2track.setParam(val, i); + } + + // Sanity check that par[2] sin(phi) is in [-1, 1] + if (std::fabs(o2track.getParam(2)) > 1.) { + LOGF(warn, "smearTrack failed sin(phi) sanity check: %f", o2track.getParam(2)); + } + + // Set covariance matrix + static constexpr int kCovMatSize = 15; + for (int i = 0; i < kCovMatSize; ++i) { + o2track.setCov(lutEntry->covm[i], i); + } + + return isReconstructed; +} + +bool TrackSmearer::smearTrack(O2Track& o2track, int pdg, float nch) +{ + auto pt = o2track.getPt(); + switch (pdg) { + case o2::constants::physics::kHelium3: + case -o2::constants::physics::kHelium3: + pt *= 2.f; + break; + } + + auto eta = o2track.getEta(); + float interpolatedEff = 0.0f; + const lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0.f, eta, pt, interpolatedEff); + + if (!lutEntry || !lutEntry->valid) { + return false; + } + + return smearTrack(o2track, lutEntry, interpolatedEff); +} + +double TrackSmearer::getPtRes(const int pdg, const float nch, const float eta, const float pt) const +{ + float dummy = 0.0f; + const lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0.f, eta, pt, dummy); + auto val = std::sqrt(lutEntry->covm[14]) * lutEntry->pt; + return val; +} + +double TrackSmearer::getEtaRes(const int pdg, const float nch, const float eta, const float pt) const +{ + float dummy = 0.0f; + const lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0.f, eta, pt, dummy); + auto sigmatgl = std::sqrt(lutEntry->covm[9]); // sigmatgl2 + auto etaRes = std::fabs(std::sin(2.0 * std::atan(std::exp(-eta)))) * sigmatgl; // propagate tgl to eta uncertainty + etaRes /= lutEntry->eta; // relative uncertainty + return etaRes; +} + +double TrackSmearer::getAbsPtRes(const int pdg, const float nch, const float eta, const float pt) const +{ + float dummy = 0.0f; + const lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0.f, eta, pt, dummy); + auto val = std::sqrt(lutEntry->covm[14]) * lutEntry->pt * lutEntry->pt; + return val; +} + +double TrackSmearer::getAbsEtaRes(const int pdg, const float nch, const float eta, const float pt) const +{ + float dummy = 0.0f; + const lutEntry_t* lutEntry = getLUTEntry(pdg, nch, 0.f, eta, pt, dummy); + auto sigmatgl = std::sqrt(lutEntry->covm[9]); // sigmatgl2 + auto etaRes = std::fabs(std::sin(2.0 * std::atan(std::exp(-eta)))) * sigmatgl; // propagate tgl to eta uncertainty + return etaRes; +} + +double TrackSmearer::getEfficiency(const int pdg, const float nch, const float eta, const float pt) const +{ + float efficiency = 0.0f; + (void)getLUTEntry(pdg, nch, 0.f, eta, pt, efficiency); + return efficiency; +} } // namespace o2::delphes diff --git a/ALICE3/Core/FlatTrackSmearer.h b/ALICE3/Core/FlatTrackSmearer.h index 2d4bb8428b3..1d56715961d 100644 --- a/ALICE3/Core/FlatTrackSmearer.h +++ b/ALICE3/Core/FlatTrackSmearer.h @@ -37,7 +37,7 @@ class TrackSmearer void useEfficiency(bool val) { mUseEfficiency = val; } void interpolateEfficiency(bool val) { mInterpolateEfficiency = val; } void skipUnreconstructed(bool val) { mSkipUnreconstructed = val; } - void setWhatEfficiency(int val) { mWhatEfficiency = val; } + void setWhatEfficiency(int val); const lutHeader_t* getLUTHeader(int pdg) const; const lutEntry_t* getLUTEntry(int pdg, float nch, float radius, float eta, float pt, float& interpolatedEff) const; @@ -60,7 +60,6 @@ class TrackSmearer protected: static constexpr unsigned int nLUTs = 9; // Number of LUT available - lutHeader_t const* mHeaders[nLUTs]; // header references for quick access FlatLutData mLUTData[nLUTs]; // Flat data storage bool mUseEfficiency = true; From 393cf92e63a459336d7cbd9137fa6e8493d306e5 Mon Sep 17 00:00:00 2001 From: Anton Alkin Date: Fri, 10 Apr 2026 12:25:34 +0200 Subject: [PATCH 08/11] update CMakeLists.txt --- ALICE3/Core/CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/ALICE3/Core/CMakeLists.txt b/ALICE3/Core/CMakeLists.txt index 7f8672c7d64..a00ac8d02dc 100644 --- a/ALICE3/Core/CMakeLists.txt +++ b/ALICE3/Core/CMakeLists.txt @@ -12,7 +12,7 @@ o2physics_add_library(ALICE3Core SOURCES TOFResoALICE3.cxx TrackUtilities.cxx - DelphesO2TrackSmearer.cxx + FlatTrackSmearer.cxx GeometryContainer.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore) @@ -20,7 +20,7 @@ o2physics_add_library(ALICE3Core o2physics_target_root_dictionary(ALICE3Core HEADERS TOFResoALICE3.h TrackUtilities.h - DelphesO2TrackSmearer.h + FlatLutEntry.h GeometryContainer.h LINKDEF ALICE3CoreLinkDef.h) From 6d2032368825f87633fa62096a77a037bb4d5571 Mon Sep 17 00:00:00 2001 From: Anton Alkin Date: Fri, 10 Apr 2026 13:00:40 +0200 Subject: [PATCH 09/11] add updated writer --- ALICE3/Core/FlatLutWriter.cxx | 588 ++++++++++++++++++++++++++++++++++ ALICE3/Core/FlatLutWriter.h | 91 ++++++ 2 files changed, 679 insertions(+) create mode 100644 ALICE3/Core/FlatLutWriter.cxx create mode 100644 ALICE3/Core/FlatLutWriter.h diff --git a/ALICE3/Core/FlatLutWriter.cxx b/ALICE3/Core/FlatLutWriter.cxx new file mode 100644 index 00000000000..3a4df00ff00 --- /dev/null +++ b/ALICE3/Core/FlatLutWriter.cxx @@ -0,0 +1,588 @@ +// 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. + +#include "ALICE3/Core/FlatLutWriter.h" + +#include "TrackUtilities.h" + +#include "ALICE3/Core/FlatTrackSmearer.h" + +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include + +#include + +using namespace o2::delphes; + +namespace o2::fastsim +{ +void FlatLutWriter::print() const +{ + LOG(info) << " --- Printing configuration of LUT writer --- "; + LOG(info) << " -> etaMaxBarrel = " << etaMaxBarrel; + LOG(info) << " -> usePara = " << usePara; + LOG(info) << " -> useDipole = " << useDipole; + LOG(info) << " -> useFlatDipole = " << useFlatDipole; + LOG(info) << " -> mAtLeastHits = " << mAtLeastHits; + LOG(info) << " -> mAtLeastCorr = " << mAtLeastCorr; + LOG(info) << " -> mAtLeastFake = " << mAtLeastFake; + LOG(info) << " -> Nch Binning: = " << mNchBinning.toString(); + LOG(info) << " -> Radius Binning: = " << mRadiusBinning.toString(); + LOG(info) << " -> Eta Binning: = " << mEtaBinning.toString(); + LOG(info) << " -> Pt Binning: = " << mPtBinning.toString(); + LOG(info) << " --- End of configuration --- "; +} + +std::string FlatLutWriter::LutBinning::toString() const +{ + std::string str = ""; + str.append(log ? "log" : "lin"); + str.append(" nbins: "); + str.append(std::to_string(nbins)); + str.append(" min: "); + str.append(std::to_string(min)); + str.append(" max: "); + str.append(std::to_string(max)); + return str; +} + +bool FlatLutWriter::fatSolve(lutEntry_t& lutEntry, + float pt, + float eta, + const float mass, + size_t itof, + size_t otof, + int q, + const float nch) +{ + lutEntry.valid = false; + + static TLorentzVector tlv; + tlv.SetPtEtaPhiM(pt, eta, 0.f, mass); + o2::track::TrackParCov trkIn; + o2::upgrade::convertTLorentzVectorToO2Track(q, tlv, {0.f, 0.f, 0.f}, trkIn); + + o2::track::TrackParCov trkOut; + const int status = fat.FastTrack(trkIn, trkOut, nch); + if (status <= mAtLeastHits) { + LOGF(debug, "fatSolve: FastTrack failed with status %d (threshold %d)", status, mAtLeastHits); + return false; + } + + LOGF(debug, "fatSolve: FastTrack succeeded with status %d", status); + + lutEntry.valid = true; + lutEntry.itof = fat.GetGoodHitProb(itof); + lutEntry.otof = fat.GetGoodHitProb(otof); + + static constexpr int nCov = 15; + for (int i = 0; i < nCov; ++i) + lutEntry.covm[i] = trkOut.getCov()[i]; + + // Define the efficiency + auto totfake = 0.f; + lutEntry.eff = 1.f; + for (size_t i = 1; i < fat.GetNLayers(); ++i) { + if (fat.IsLayerInert(i)) + continue; // skip inert layers + auto igoodhit = fat.GetGoodHitProb(i); + if (igoodhit <= 0.f || i == itof || i == otof) + continue; + lutEntry.eff *= igoodhit; + auto pairfake = 0.f; + for (size_t j = i + 1; j < fat.GetNLayers(); ++j) { + auto jgoodhit = fat.GetGoodHitProb(j); + if (jgoodhit <= 0.f || j == itof || j == otof) + continue; + pairfake = (1.f - igoodhit) * (1.f - jgoodhit); + break; + } + totfake += pairfake; + } + lutEntry.eff2 = (1.f - totfake); + + return true; +} + +#ifdef USE_FWD_PARAM +bool FlatLutWriter::fwdSolve(float* covm, float pt, float eta, float mass) +{ + if (fwdRes(covm, pt, eta, mass) < 0) + return false; + return true; +} +#else +bool FlatLutWriter::fwdSolve(float*, float, float, float) +{ + return false; +} +#endif + +bool FlatLutWriter::fwdPara(lutEntry_t& lutEntry, float pt, float eta, float mass, float Bfield) +{ + lutEntry.valid = false; + + // Parametrised forward response; interpolates between FAT at eta = 1.75 and a fixed parametrisation at eta = 4 + // Only diagonal elements + static constexpr float etaLimit = 4.0f; + if (std::fabs(eta) < etaMaxBarrel || std::fabs(eta) > etaLimit) { + return false; + } + + if (!fatSolve(lutEntry, pt, etaMaxBarrel, mass)) { + return false; + } + + static constexpr int nCov = 15; + float covmbarrel[nCov] = {0.f}; + for (int i = 0; i < nCov; ++i) { + covmbarrel[i] = lutEntry.covm[i]; + } + + // Parametrisation at eta = 4 + const double beta = 1. / std::sqrt(1.0 + mass * mass / pt / pt / std::cosh(eta) / std::cosh(eta)); + const float dcaPos = 2.5e-4f / std::sqrt(3.f); // 2.5 micron/sqrt(3) + const float r0 = 0.5f; // layer 0 radius [cm] + const float r1 = 1.3f; + const float r2 = 2.5f; + const float x0layer = 0.001f; // material budget (rad length) per layer + + const double sigmaAlpha = 0.0136 / beta / pt * std::sqrt(x0layer * std::cosh(eta)) * (1.0 + 0.038 * std::log(x0layer * std::cosh(eta))); + const double dcaxyMs = sigmaAlpha * r0 * std::sqrt(1.0 + r1 * r1 / (r2 - r0) / (r2 - r0)); + const double dcaxy2 = dcaPos * dcaPos + dcaxyMs * dcaxyMs; + + const double dcazMs = sigmaAlpha * r0 * std::cosh(eta); + const double dcaz2 = dcaPos * dcaPos + dcazMs * dcazMs; + + const float Leta = 2.8f / std::sinh(eta) - 0.01f * r0; // m + const double relmomresPos = 10e-6 * pt / 0.3 / Bfield / Leta / Leta * std::sqrt(720.0 / 15.0); + + const float relmomresBarrel = std::sqrt(covmbarrel[14]) * pt; + const float rOuter = 1.f; // m + const float relmomresPosBarrel = 10e-6f * pt / 0.3f / Bfield / rOuter / rOuter / std::sqrt(720.f / 15.f); + const float relmomresMSBarrel = std::sqrt(relmomresBarrel * relmomresBarrel - relmomresPosBarrel * relmomresPosBarrel); + + // Interpolate MS contribution (rel resolution 0.4 at eta = 4) + const float relmomresMSEta4 = 0.4f / beta * 0.5f / Bfield; + const float relmomresMS = relmomresMSEta4 * std::pow(relmomresMSEta4 / relmomresMSBarrel, (std::fabs(eta) - 4.f) / (4.f - etaMaxBarrel)); + const float momresTot = pt * std::sqrt(relmomresPos * relmomresPos + relmomresMS * relmomresMS); // total absolute mom reso + + // Fill cov matrix diag + for (int i = 0; i < 15; ++i) { + lutEntry.covm[i] = 0.f; + } + + lutEntry.covm[0] = covmbarrel[0]; + if (dcaxy2 > lutEntry.covm[0]) { + lutEntry.covm[0] = dcaxy2; + } + lutEntry.covm[2] = covmbarrel[2]; + if (dcaz2 > lutEntry.covm[2]) { + lutEntry.covm[2] = dcaz2; + } + lutEntry.covm[5] = covmbarrel[5]; // sigma^2 sin(phi) + lutEntry.covm[9] = covmbarrel[9]; // sigma^2 tanl + lutEntry.covm[14] = momresTot * momresTot / pt / pt / pt / pt; // sigma^2 1/pt + + // Check that all numbers are numbers + for (int i = 0; i < 15; ++i) { + if (std::isnan(lutEntry.covm[i])) { + LOGF(info, "lutEntry.covm[%d] is NaN", i); + return false; + } + } + return true; +} + +void FlatLutWriter::lutWrite(const char* filename, int pdg, float field, size_t itof, size_t otof) +{ + if (useFlatDipole && useDipole) { + LOGF(error, "Both dipole and dipole flat flags are on, please use only one of them"); + return; + } + + // Open output file for binary writing + std::ofstream lutFile(filename, std::ofstream::binary); + if (!lutFile.is_open()) { + LOGF(error, "Failed to open output file: %s", filename); + return; + } + + // Create and write header + lutHeader_t lutHeader; + lutHeader.pdg = pdg; + + const TParticlePDG* particle = TDatabasePDG::Instance()->GetParticle(pdg); + if (!particle) { + LOGF(fatal, "Cannot find particle with PDG code %d", pdg); + lutFile.close(); + return; + } + + lutHeader.mass = particle->Mass(); + const int q = std::abs(particle->Charge()) / 3; + if (q <= 0) { + LOGF(error, "Negative or null charge (%f) for pdg code %d", particle->Charge(), pdg); + lutFile.close(); + return; + } + + lutHeader.field = field; + + // Set binning maps in header + auto setMap = [](map_t& map, const LutBinning& b) { + map.log = b.log; + map.nbins = b.nbins; + map.min = b.min; + map.max = b.max; + }; + + setMap(lutHeader.nchmap, mNchBinning); + setMap(lutHeader.radmap, mRadiusBinning); + setMap(lutHeader.etamap, mEtaBinning); + setMap(lutHeader.ptmap, mPtBinning); + + // Write header to file + lutFile.write(reinterpret_cast(&lutHeader), sizeof(lutHeader_t)); + if (!lutFile.good()) { + LOGF(error, "Failed to write LUT header to %s", filename); + lutFile.close(); + return; + } + + // Get dimensions + const int nnch = lutHeader.nchmap.nbins; + const int nrad = lutHeader.radmap.nbins; + const int neta = lutHeader.etamap.nbins; + const int npt = lutHeader.ptmap.nbins; + + LOGF(info, "Writing LUT with dimensions: nch=%d rad=%d eta=%d pt=%d (total=%zu entries)", nnch, nrad, neta, npt, static_cast(nnch) * nrad * neta * npt); + + lutEntry_t lutEntry; + int nCalls = 0; + int successfulCalls = 0; + int failedCalls = 0; + + // Write all entries sequentially + for (int inch = 0; inch < nnch; ++inch) { + LOGF(info, "Writing nch bin %d/%d", inch, nnch); + auto nch = lutHeader.nchmap.eval(inch); + lutEntry.nch = nch; + fat.SetdNdEtaCent(nch); + + for (int irad = 0; irad < nrad; ++irad) { + for (int ieta = 0; ieta < neta; ++ieta) { + auto eta = lutHeader.etamap.eval(ieta); + lutEntry.eta = eta; + + for (int ipt = 0; ipt < npt; ++ipt) { + nCalls++; + lutEntry.pt = lutHeader.ptmap.eval(ipt); + lutEntry.valid = true; + + // Solve for this bin + if (std::fabs(eta) <= etaMaxBarrel) { + // Full lever arm region (barrel) + LOGF(debug, "Solving barrel: pt=%f eta=%f", lutEntry.pt, lutEntry.eta); + successfulCalls++; + + if (!fatSolve(lutEntry, lutEntry.pt, lutEntry.eta, lutHeader.mass, itof, otof, q, nch)) { + lutEntry.valid = false; + lutEntry.eff = 0.f; + lutEntry.eff2 = 0.f; + for (int i = 0; i < 15; ++i) { + lutEntry.covm[i] = 0.f; + } + successfulCalls--; + failedCalls++; + } + } else { + // Forward region + LOGF(debug, "Solving forward: pt=%f eta=%f", lutEntry.pt, lutEntry.eta); + lutEntry.eff = 1.f; + lutEntry.eff2 = 1.f; + bool retval = true; + successfulCalls++; + + if (useFlatDipole) { + // Using the parametrization at the border of the barrel + retval = fatSolve(lutEntry, lutEntry.pt, etaMaxBarrel, lutHeader.mass, itof, otof, q, nch); + } else if (usePara) { + retval = fwdPara(lutEntry, lutEntry.pt, lutEntry.eta, lutHeader.mass, field); + } else { + retval = fwdSolve(lutEntry.covm, lutEntry.pt, lutEntry.eta, lutHeader.mass); + } + + if (useDipole) { + // Using the parametrization at the border of the barrel only for efficiency and momentum resolution + lutEntry_t lutEntryBarrel; + retval = fatSolve(lutEntryBarrel, lutEntry.pt, etaMaxBarrel, lutHeader.mass, itof, otof, q, nch); + lutEntry.valid = lutEntryBarrel.valid; + lutEntry.covm[14] = lutEntryBarrel.covm[14]; + lutEntry.eff = lutEntryBarrel.eff; + lutEntry.eff2 = lutEntryBarrel.eff2; + } + + if (!retval) { + LOGF(debug, "Forward solve failed"); + lutEntry.valid = false; + for (int i = 0; i < 15; ++i) { + lutEntry.covm[i] = 0.f; + } + successfulCalls--; + failedCalls++; + } + } + + // Diagonalize covariance matrix + diagonalise(lutEntry); + + // Write entry to file + lutFile.write(reinterpret_cast(&lutEntry), sizeof(lutEntry_t)); + if (!lutFile.good()) { + LOGF(error, "Failed to write LUT entry at index [%d,%d,%d,%d]", inch, irad, ieta, ipt); + lutFile.close(); + return; + } + } + } + } + } + + lutFile.close(); + + LOGF(info, "Finished writing LUT file: %s", filename); + LOGF(info, "Total calls: %d, successful: %d, failed: %d", nCalls, successfulCalls, failedCalls); +} + +void FlatLutWriter::diagonalise(lutEntry_t& lutEntry) +{ + static constexpr int kEig = 5; + TMatrixDSym m(kEig); + + for (int i = 0, k = 0; i < kEig; ++i) { + for (int j = 0; j < i + 1; ++j, ++k) { + m(i, j) = lutEntry.covm[k]; + m(j, i) = lutEntry.covm[k]; + } + } + + TMatrixDSymEigen eigen(m); + + // Eigenvalues + TVectorD eigenVal = eigen.GetEigenValues(); + for (int i = 0; i < kEig; ++i) + lutEntry.eigval[i] = eigenVal[i]; + + // Eigenvectors + TMatrixD eigenVec = eigen.GetEigenVectors(); + for (int i = 0; i < kEig; ++i) + for (int j = 0; j < kEig; ++j) + lutEntry.eigvec[i][j] = eigenVec[i][j]; + + // Inverse eigenvectors + eigenVec.Invert(); + for (int i = 0; i < kEig; ++i) + for (int j = 0; j < kEig; ++j) + lutEntry.eiginv[i][j] = eigenVec[i][j]; +} + +TGraph* FlatLutWriter::lutRead(const char* filename, int pdg, int what, int vs, float nch, float radius, float eta, float pt) +{ + LOGF(info, "Reading LUT file: %s", filename); + + static const int kNch = 0; + static const int kEta = 1; + static const int kPt = 2; + + static const int kEfficiency = 0; + static const int kEfficiency2 = 1; + static const int kEfficiencyInnerTOF = 2; + static const int kEfficiencyOuterTOF = 3; + static const int kPtResolution = 4; + static const int kRPhiResolution = 5; + static const int kZResolution = 6; + + // Use TrackSmearer to load and access the LUT + o2::delphes::TrackSmearer smearer; + if (!smearer.loadTable(pdg, filename)) { + LOGF(error, "Failed to load LUT from %s", filename); + return nullptr; + } + + const lutHeader_t* lutHeader = smearer.getLUTHeader(pdg); + if (!lutHeader) { + LOGF(error, "No LUT header for PDG %d", pdg); + return nullptr; + } + + lutHeader->print(); + + map_t lutMap; + switch (vs) { + case kNch: + lutMap = lutHeader->nchmap; + break; + case kEta: + lutMap = lutHeader->etamap; + break; + case kPt: + lutMap = lutHeader->ptmap; + break; + default: + LOGF(error, "Unknown vs: %d", vs); + return nullptr; + } + + auto nbins = lutMap.nbins; + auto g = new TGraph(); + g->SetName(Form("lut_%s_%d_vs_%d_what_%d", filename, pdg, vs, what)); + g->SetTitle(Form("LUT for %s, pdg %d, vs %d, what %d", filename, pdg, vs, what)); + + // Set axis labels + switch (vs) { + case kNch: + LOGF(info, "Plot versus Nch"); + g->GetXaxis()->SetTitle("Nch"); + break; + case kEta: + LOGF(info, "Plot versus Eta"); + g->GetXaxis()->SetTitle("#eta"); + break; + case kPt: + LOGF(info, "Plot versus Pt"); + g->GetXaxis()->SetTitle("p_{T} (GeV/c)"); + break; + } + + switch (what) { + case kEfficiency: + LOGF(info, "Plot efficiency"); + g->GetYaxis()->SetTitle("Efficiency (%)"); + break; + case kEfficiency2: + LOGF(info, "Plot efficiency2"); + g->GetYaxis()->SetTitle("Efficiency2 (%)"); + break; + case kEfficiencyInnerTOF: + LOGF(info, "Plot inner TOF efficiency"); + g->GetYaxis()->SetTitle("Inner TOF Efficiency (%)"); + break; + case kEfficiencyOuterTOF: + LOGF(info, "Plot outer TOF efficiency"); + g->GetYaxis()->SetTitle("Outer TOF Efficiency (%)"); + break; + case kPtResolution: + LOGF(info, "Plot pt resolution"); + g->GetYaxis()->SetTitle("p_{T} Resolution (%)"); + break; + case kRPhiResolution: + LOGF(info, "Plot rphi resolution"); + g->GetYaxis()->SetTitle("R#phi Resolution (#mum)"); + break; + case kZResolution: + LOGF(info, "Plot z resolution"); + g->GetYaxis()->SetTitle("Z Resolution (#mum)"); + break; + default: + LOGF(error, "Unknown what: %d", what); + delete g; + return nullptr; + } + + // Fill graph by iterating over one dimension + bool canBeInvalid = true; + for (int i = 0; i < nbins; ++i) { + switch (vs) { + case kNch: + nch = lutMap.eval(i); + break; + case kEta: + eta = lutMap.eval(i); + break; + case kPt: + pt = lutMap.eval(i); + break; + } + + float eff = 0.f; + const lutEntry_t* lutEntry = smearer.getLUTEntry(pdg, nch, radius, eta, pt, eff); + + if (!lutEntry || !lutEntry->valid || lutEntry->eff == 0.f) { + if (!canBeInvalid) { + LOGF(warning, "Entry became invalid at bin %d", i); + } + continue; + } + canBeInvalid = false; + + double cen = 0.f; + switch (vs) { + case kNch: + cen = lutEntry->nch; + break; + case kEta: + cen = lutEntry->eta; + break; + case kPt: + cen = lutEntry->pt; + break; + } + + double val = 0.f; + switch (what) { + case kEfficiency: + val = lutEntry->eff * 100.f; // efficiency (%) + break; + case kEfficiency2: + val = lutEntry->eff2 * 100.f; // efficiency (%) + break; + case kEfficiencyInnerTOF: + val = lutEntry->itof * 100.f; // efficiency (%) + break; + case kEfficiencyOuterTOF: + val = lutEntry->otof * 100.f; // efficiency (%) + break; + case kPtResolution: + val = std::sqrt(lutEntry->covm[14]) * lutEntry->pt * 100.f; // pt resolution (%) + break; + case kRPhiResolution: + val = std::sqrt(lutEntry->covm[0]) * 1.e4f; // rphi resolution (um) + break; + case kZResolution: + val = std::sqrt(lutEntry->covm[1]) * 1.e4f; // z resolution (um) + break; + default: + LOGF(error, "Unknown what: %d", what); + break; + } + g->AddPoint(cen, val); + } + + LOGF(info, "Read %d points from LUT", g->GetN()); + return g; +} + +} // namespace o2::fastsim diff --git a/ALICE3/Core/FlatLutWriter.h b/ALICE3/Core/FlatLutWriter.h new file mode 100644 index 00000000000..f2240b160c4 --- /dev/null +++ b/ALICE3/Core/FlatLutWriter.h @@ -0,0 +1,91 @@ +// 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. + +#ifndef ALICE3_CORE_FLATLUTWRITER_H_ +#define ALICE3_CORE_FLATLUTWRITER_H_ + +#include "ALICE3/Core/FastTracker.h" +#include "ALICE3/Core/FlatLutEntry.h" +namespace o2::fastsim +{ +using lutEntry_t = o2::delphes::lutEntry_t; +/** + * @brief LUT writer using flat binary format + * + * Generates lookup tables in the optimized flat format compatible with FlatLutData. + * The output format is: + * [lutHeader_t][lutEntry_t_0][lutEntry_t_1]...[lutEntry_t_N] + * + * This enables zero-copy loading into the TrackSmearer and direct shared memory mapping. + */ +class FlatLutWriter +{ + public: + FlatLutWriter() = default; + + // Setters for binning configuration + void setBinningNch(bool log, int nbins, float min, float max) { mNchBinning = {log, nbins, min, max}; } + void setBinningRadius(bool log, int nbins, float min, float max) { mRadiusBinning = {log, nbins, min, max}; } + void setBinningEta(bool log, int nbins, float min, float max) { mEtaBinning = {log, nbins, min, max}; } + void setBinningPt(bool log, int nbins, float min, float max) { mPtBinning = {log, nbins, min, max}; } + + void setEtaMaxBarrel(float eta) { etaMaxBarrel = eta; } + void setAtLeastHits(int n) { mAtLeastHits = n; } + void setAtLeastCorr(int n) { mAtLeastCorr = n; } + void setAtLeastFake(int n) { mAtLeastFake = n; } + + bool fatSolve(lutEntry_t& lutEntry, + float pt = 0.1f, + float eta = 0.0f, + const float mass = o2::track::pid_constants::sMasses[o2::track::PID::Pion], + size_t itof = 0, + size_t otof = 0, + int q = 1, + const float nch = 1.0f); + + void print() const; + bool fwdSolve(float* covm, float pt = 0.1f, float eta = 0.0f, float mass = o2::track::pid_constants::sMasses[o2::track::PID::Pion]); + bool fwdPara(lutEntry_t& lutEntry, float pt = 0.1f, float eta = 0.0f, float mass = o2::track::pid_constants::sMasses[o2::track::PID::Pion], float Bfield = 0.5f); + void lutWrite(const char* filename = "lutCovm.dat", int pdg = 211, float field = 0.2f, size_t itof = 0, size_t otof = 0); + TGraph* lutRead(const char* filename, int pdg, int what, int vs, float nch = 0.f, float radius = 0.f, float eta = 0.f, float pt = 0.f); + + o2::fastsim::FastTracker fat; + + private: + void diagonalise(lutEntry_t& lutEntry); + + float etaMaxBarrel = 1.75f; + bool usePara = true; // use fwd parametrisation + bool useDipole = false; // use dipole i.e. flat parametrization for efficiency and momentum resolution + bool useFlatDipole = false; // use dipole i.e. flat parametrization outside of the barrel + + int mAtLeastHits = 4; + int mAtLeastCorr = 4; + int mAtLeastFake = 0; + + // Binning of the LUT to make + struct LutBinning { + bool log; + int nbins; + float min; + float max; + std::string toString() const; + }; + + LutBinning mNchBinning = {true, 20, 0.5f, 3.5f}; + LutBinning mRadiusBinning = {false, 1, 0.0f, 100.0f}; + LutBinning mEtaBinning = {false, 80, -4.0f, 4.0f}; + LutBinning mPtBinning = {true, 200, -2.0f, 2.0f}; +}; + +} // namespace o2::fastsim + +#endif // ALICE3_CORE_FLATLUTWRITER_H_ From 52307467ae80a540b73d21bce02859122f5278e3 Mon Sep 17 00:00:00 2001 From: Anton Alkin Date: Fri, 10 Apr 2026 13:02:12 +0200 Subject: [PATCH 10/11] update CMakeLists.txt and link def --- ALICE3/Core/CMakeLists.txt | 4 ++-- ALICE3/Core/FastTrackerLinkDef.h | 2 +- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/ALICE3/Core/CMakeLists.txt b/ALICE3/Core/CMakeLists.txt index a00ac8d02dc..37c722e12ff 100644 --- a/ALICE3/Core/CMakeLists.txt +++ b/ALICE3/Core/CMakeLists.txt @@ -27,7 +27,7 @@ o2physics_target_root_dictionary(ALICE3Core o2physics_add_library(FastTracker SOURCES FastTracker.cxx DetLayer.cxx - DelphesO2LutWriter.cxx + FlatLutWriter.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::ALICE3Core) @@ -35,5 +35,5 @@ o2physics_add_library(FastTracker o2physics_target_root_dictionary(FastTracker HEADERS FastTracker.h DetLayer.h - DelphesO2LutWriter.h + FlatLutWriter.h LINKDEF FastTrackerLinkDef.h) diff --git a/ALICE3/Core/FastTrackerLinkDef.h b/ALICE3/Core/FastTrackerLinkDef.h index 4986a1879f7..9885d281a4d 100644 --- a/ALICE3/Core/FastTrackerLinkDef.h +++ b/ALICE3/Core/FastTrackerLinkDef.h @@ -18,6 +18,6 @@ #pragma link C++ class o2::fastsim::GeometryContainer + ; #pragma link C++ class o2::fastsim::FastTracker + ; -#pragma link C++ class o2::fastsim::DelphesO2LutWriter + ; +#pragma link C++ class o2::fastsim::FlatLutWriter + ; #endif // ALICE3_CORE_FASTTRACKERLINKDEF_H_ From af0bd6aa4f2a7035aaedd911ed3b65f7c89bba18 Mon Sep 17 00:00:00 2001 From: Anton Alkin Date: Fri, 10 Apr 2026 13:05:01 +0200 Subject: [PATCH 11/11] make it compile --- ALICE3/Core/CMakeLists.txt | 2 ++ 1 file changed, 2 insertions(+) diff --git a/ALICE3/Core/CMakeLists.txt b/ALICE3/Core/CMakeLists.txt index 37c722e12ff..3e9358af98c 100644 --- a/ALICE3/Core/CMakeLists.txt +++ b/ALICE3/Core/CMakeLists.txt @@ -12,6 +12,7 @@ o2physics_add_library(ALICE3Core SOURCES TOFResoALICE3.cxx TrackUtilities.cxx + FlatLutEntry.cxx FlatTrackSmearer.cxx GeometryContainer.cxx PUBLIC_LINK_LIBRARIES O2::Framework @@ -27,6 +28,7 @@ o2physics_target_root_dictionary(ALICE3Core o2physics_add_library(FastTracker SOURCES FastTracker.cxx DetLayer.cxx + FlatLutEntry.cxx FlatLutWriter.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore