[ALICE3] Add LUT utilities

ALICE3/Core/FlatLutEntry.cxx

@@ -9,6 +9,9 @@
 // granted to it by virtue of its status as an Intergovernmental Organization
 // or submit itself to any jurisdiction.
 
+/// \file FlatLutEntry.cxx
+/// \brief Flat LUT implementation for compact helper tables used by the ALICE3 track smearing workflow.
+
 #include "FlatLutEntry.h"
 
 #include <Framework/Logger.h>
@@ -24,12 +27,48 @@
 
 namespace o2::delphes
 {
+namespace
+{
+constexpr int kNumCovarianceTerms = 15;
+constexpr int kNumEigenModes = 5;
+constexpr float kBinCenterOffset = 0.5f;
+} // namespace
+
+void lutEntry_t::print() const
+{
+  LOGF(info, "  nch = %f, eta = %f, pt = %f, valid = %s\n", nch, eta, pt, valid ? "true" : "false");
+  LOGF(info, "  eff = %f, eff2 = %f, itof = %f, otof = %f\n", eff, eff2, itof, otof);
+  LOGF(info, "  covm: ");
+  for (int i = 0; i < kNumCovarianceTerms; ++i) {
+    LOGF(info, "%f ", covm[i]);
+  }
+  LOGF(info, "\n");
+  LOGF(info, "  eigval: ");
+  for (int i = 0; i < kNumEigenModes; ++i) {
+    LOGF(info, "%f ", eigval[i]);
+  }
+  LOGF(info, "\n");
+  LOGF(info, "  eigvec:\n");
+  for (int i = 0; i < kNumEigenModes; ++i) {
+    for (int j = 0; j < kNumEigenModes; ++j) {
+      LOGF(info, "%f ", eigvec[i][j]);
+    }
+    LOGF(info, "\n");
+  }
+  LOGF(info, "  eiginv:\n");
+  for (int i = 0; i < kNumEigenModes; ++i) {
+    for (int j = 0; j < kNumEigenModes; ++j) {
+      LOGF(info, "%f ", eiginv[i][j]);
+    }
+    LOGF(info, "\n");
+  }
+}
 
 float map_t::fracPositionWithinBin(float val) const
 {
   float width = (max - min) / nbins;
   int bin;
-  float returnVal = 0.5f;
+  float returnVal = kBinCenterOffset;
   if (log) {
     bin = static_cast<int>((std::log10(val) - min) / width);
     returnVal = ((std::log10(val) - min) / width) - bin;
@@ -63,7 +102,7 @@ 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
+bool lutHeader_t::checkVersion() const
 {
   return (version == LUTCOVM_VERSION);
 }
@@ -90,10 +129,10 @@ void FlatLutData::initialize(const lutHeader_t& header)
   mEtaBins = header.etamap.nbins;
   mPtBins = header.ptmap.nbins;
 
-  size_t headerSize = sizeof(lutHeader_t);
-  size_t numEntries = static_cast<size_t>(mNchBins) * mRadBins * mEtaBins * mPtBins;
-  size_t entriesSize = numEntries * sizeof(lutEntry_t);
-  size_t totalSize = headerSize + entriesSize;
+  constexpr size_t headerSize = sizeof(lutHeader_t);
+  const size_t numEntries = static_cast<size_t>(mNchBins) * mRadBins * mEtaBins * mPtBins;
+  const size_t entriesSize = numEntries * sizeof(lutEntry_t);
+  const size_t totalSize = headerSize + entriesSize;
 
   mData.resize(totalSize);
   // Write header at the beginning
@@ -103,13 +142,10 @@ void FlatLutData::initialize(const lutHeader_t& header)
 
 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<size_t>(nch_bin) * (mRadBins * mEtaBins * mPtBins) + static_cast<size_t>(rad_bin) * (mEtaBins * mPtBins) + static_cast<size_t>(eta_bin) * mPtBins + static_cast<size_t>(pt_bin);
-
-  return headerSize + linearIdx * sizeof(lutEntry_t);
+  static constexpr size_t HeaderSize = sizeof(lutHeader_t);
+  const size_t linearIdx = getEntryIndex(nch_bin, rad_bin, eta_bin, pt_bin);
+  static constexpr size_t EntrySize = sizeof(lutEntry_t);
+  return HeaderSize + linearIdx * EntrySize;
 }
 
 const lutEntry_t* FlatLutData::getEntryRef(int nch_bin, int rad_bin, int eta_bin, int pt_bin) const
@@ -173,26 +209,26 @@ void FlatLutData::view(const uint8_t* buffer, size_t size)
 
 void FlatLutData::validateBuffer(const uint8_t* buffer, size_t size)
 {
-  auto header = PreviewHeader(buffer, size);
-  auto mNchBins = header.nchmap.nbins;
-  auto mRadBins = header.radmap.nbins;
-  auto mEtaBins = header.etamap.nbins;
-  auto mPtBins = header.ptmap.nbins;
+  auto header = previewHeader(buffer, size);
+  const auto nchBins = header.nchmap.nbins;
+  const auto radBins = header.radmap.nbins;
+  const auto etaBins = header.etamap.nbins;
+  const auto ptBins = header.ptmap.nbins;
 
-  size_t expectedSize = sizeof(lutHeader_t) + static_cast<size_t>(mNchBins) * mRadBins * mEtaBins * mPtBins * sizeof(lutEntry_t);
+  const size_t expectedSize = sizeof(lutHeader_t) + static_cast<size_t>(nchBins) * radBins * etaBins * ptBins * sizeof(lutEntry_t);
 
   if (size < expectedSize) {
     throw framework::runtime_error_f("Buffer size mismatch: expected %zu, got %zu", expectedSize, size);
   }
 }
 
-lutHeader_t FlatLutData::PreviewHeader(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<const lutHeader_t*>(buffer);
-  if (!header->check_version()) {
+  if (!header->checkVersion()) {
     throw framework::runtime_error_f("LUT header version mismatch: expected %d, got %d", LUTCOVM_VERSION, header->version);
   }
   return *header;
@@ -229,14 +265,14 @@ bool FlatLutData::isLoaded() const
   return ((!mData.empty()) || (!mDataRef.empty()));
 }
 
-lutHeader_t FlatLutData::PreviewHeader(std::ifstream& file, const char* filename)
+lutHeader_t FlatLutData::previewHeader(std::ifstream& file, const char* filename)
 {
   lutHeader_t tempHeader;
   file.read(reinterpret_cast<char*>(&tempHeader), sizeof(lutHeader_t));
   if (file.gcount() != static_cast<std::streamsize>(sizeof(lutHeader_t))) {
     throw framework::runtime_error_f("Failed to read LUT header from %s", filename);
   }
-  if (!tempHeader.check_version()) {
+  if (!tempHeader.checkVersion()) {
     throw framework::runtime_error_f("LUT header version mismatch: expected %d, got %d", LUTCOVM_VERSION, tempHeader.version);
   }
   return tempHeader;
@@ -245,7 +281,7 @@ lutHeader_t FlatLutData::PreviewHeader(std::ifstream& file, const char* filename
 FlatLutData FlatLutData::loadFromFile(std::ifstream& file, const char* filename)
 {
   // Read header first
-  lutHeader_t tempHeader = PreviewHeader(file, filename);
+  lutHeader_t tempHeader = previewHeader(file, filename);
 
   FlatLutData data;
 

ALICE3/Core/FlatLutEntry.h

@@ -9,6 +9,9 @@
 // granted to it by virtue of its status as an Intergovernmental Organization
 // or submit itself to any jurisdiction.
 
+/// \file FlatLutEntry.h
+/// \brief Flat LUT data structures and buffer handling for the ALICE3 fast smearing backend.
+
 #ifndef ALICE3_CORE_FLATLUTENTRY_H_
 #define ALICE3_CORE_FLATLUTENTRY_H_
 
@@ -25,9 +28,9 @@ namespace o2::delphes
 {
 
 /**
- * @brief Flat LUT entry structure
+ * @brief Flat LUT entry structure.
  */
-struct lutEntry_t {
+struct LutEntry {
   float nch = 0.f;
   float eta = 0.f;
   float pt = 0.f;
@@ -44,10 +47,12 @@ struct lutEntry_t {
   void print() const;
 };
 
+using lutEntry_t = LutEntry;
+
 /**
- * @brief Binning map
+ * @brief Binning map.
  */
-struct map_t {
+struct Map {
   int nbins = 1;
   float min = 0.f;
   float max = 1.e6f;
@@ -67,10 +72,12 @@ struct map_t {
   void print() const;
 };
 
+using map_t = Map;
+
 /**
- * @brief LUT header
+ * @brief LUT header.
  */
-struct lutHeader_t {
+struct LutHeader {
   int version = LUTCOVM_VERSION;
   int pdg = 0;
   float mass = 0.f;
@@ -80,12 +87,14 @@ struct lutHeader_t {
   map_t etamap;
   map_t ptmap;
 
-  bool check_version() const;
+  bool checkVersion() const;
   void print() const;
 };
 
+using lutHeader_t = LutHeader;
+
 /**
- * @brief Flat LUT data container - single contiguous buffer
+ * @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.
@@ -104,6 +113,13 @@ class FlatLutData
    */
   void initialize(const lutHeader_t& header);
 
+  size_t getEntryIndex(int nch_bin, int rad_bin, int eta_bin, int pt_bin) const
+  {
+    // Linear index: nch varies slowest, pt varies fastest
+    // idx = nch * (rad*eta*pt) + rad * (eta*pt) + eta * pt + pt
+    return static_cast<size_t>(nch_bin) * (mRadBins * mEtaBins * mPtBins) + static_cast<size_t>(rad_bin) * (mEtaBins * mPtBins) + static_cast<size_t>(eta_bin) * mPtBins + static_cast<size_t>(pt_bin);
+  }
+
   /**
    * @brief Get LUT entry by bin indices (view)
    */

ALICE3/Core/FlatTrackSmearer.cxx

@@ -27,6 +27,7 @@
 #include <fstream>
 #include <span>
 #include <string>
+#include <vector>
 
 namespace o2::delphes
 {
@@ -85,7 +86,7 @@
 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;
@@ -139,7 +140,7 @@
     return false;
   }
   try {
-    auto header = FlatLutData::PreviewHeader(buffer, size);
+    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;
@@ -162,7 +163,7 @@
     return false;
   }
   try {
-    auto header = FlatLutData::PreviewHeader(buffer, size);
+    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;

ALICE3/Macros/drawTrackSmearer.C

@@ -0,0 +1,96 @@
+// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
+// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
+// All rights not expressly granted are reserved.
+//
+// This software is distributed under the terms of the GNU General Public
+// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
+//
+// In applying this license CERN does not waive the privileges and immunities
+// granted to it by virtue of its status as an Intergovernmental Organization
+// or submit itself to any jurisdiction.
+
+/// \file drawTrackSmearer.C
+/// \brief Draw the ALICE3 track-smearing performance curves for key particle species.
+
+#include "FlatTrackSmearer.h"
+#include "TrackUtilities.h"
+
+#include <TAxis.h>
+#include <TCanvas.h>
+#include <TColor.h>
+#include <TGraph.h>
+#include <TPDGCode.h>
+
+namespace
+{
+constexpr int kElectronPdg = static_cast<int>(PDG_t::kElectron);
+constexpr int kMuonPdg = static_cast<int>(PDG_t::kMuonMinus);
+constexpr int kPionPdg = static_cast<int>(PDG_t::kPiPlus);
+constexpr int kProtonPdg = static_cast<int>(PDG_t::kProton);
+constexpr int kKaonPdg = static_cast<int>(PDG_t::kKPlus);
+
+const std::vector<std::pair<int, std::string>> kDefaultLutFiles = {{kElectronPdg, "/tmp/lut/lutCov.acts.11.2T.dNdEta5.dat"},
+                                                                   {kMuonPdg, "/tmp/lut/lutCov.acts.13.2T.dNdEta5.dat"},
+                                                                   {kPionPdg, "/tmp/lut/lutCov.acts.211.2T.dNdEta5.dat"},
+                                                                   {kProtonPdg, "/tmp/lut/lutCov.acts.2212.2T.dNdEta5.dat"},
+                                                                   {kKaonPdg, "/tmp/lut/lutCov.acts.321.2T.dNdEta5.dat"}};
+} // namespace
+
+void drawTrackSmearer(const std::vector<std::pair<int, std::string>>& filenames = kDefaultLutFiles)
+{
+
+  o2::delphes::TrackSmearer trackSmearer;
+
+  TCanvas* cPtReso = new TCanvas("cPtReso", "cPtReso", 800, 600);
+  TCanvas* cPtEff = new TCanvas("cPtEff", "cPtEff", 800, 600);
+
+  for (const auto& [pdg, filename] : filenames) {
+    trackSmearer.loadTable(pdg, filename.c_str(), true);
+    TGraph* gPt = new TGraph();
+    gPt->SetName(Form("gPt_%d", pdg));
+    gPt->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
+    gPt->GetYaxis()->SetTitle("pt resolution");
+
+    TGraph* gPtEff = new TGraph();
+    gPtEff->SetName(Form("gPtEff_%d", pdg));
+    gPtEff->GetXaxis()->SetTitle("#it{p}_{T} (GeV/#it{c})");
+    gPtEff->GetYaxis()->SetTitle("efficiency");
+
+    for (int i = 0; i < trackSmearer.getLUTHeader(pdg)->ptmap.nbins; i++) {
+      const float pt = trackSmearer.getLUTHeader(pdg)->ptmap.eval(i);
+      const float res = trackSmearer.getPtRes(pdg, 0, 0., pt);
+      gPt->AddPoint(pt, res / pt);
+      const float eff = trackSmearer.getEfficiency(pdg, 0, 0., pt);
+      gPtEff->AddPoint(pt, eff);
+    }
+
+    int color = 0;
+    if (pdg == kElectronPdg) {
+      color = TColor::GetColor("#e41a1c");
+    } else if (pdg == kMuonPdg) {
+      color = TColor::GetColor("#377eb8");
+    } else if (pdg == kPionPdg) {
+      color = TColor::GetColor("#4daf4a");
+    } else if (pdg == kProtonPdg) {
+      color = TColor::GetColor("#984ea3");
+    } else if (pdg == kKaonPdg) {
+      color = TColor::GetColor("#ff7f00");
+    }
+    gPt->SetLineColor(color);
+    gPtEff->SetLineColor(color);
+
+    cPtReso->cd();
+    if (cPtReso->GetListOfPrimitives()->GetEntries() == 0) {
+      gPt->Draw("ALP");
+    } else {
+      gPt->Draw("LP SAME");
+    }
+    cPtEff->cd();
+    if (cPtEff->GetListOfPrimitives()->GetEntries() == 0) {
+      gPtEff->Draw("ALP");
+    } else {
+      gPtEff->Draw("LP SAME");
+    }
+    gPt->SaveAs("/tmp/gPt.root");
+  }
+}