diff --git a/ALICE3/Core/FlatLutEntry.cxx b/ALICE3/Core/FlatLutEntry.cxx index 5110e9bff4c..0415f7630f4 100644 --- a/ALICE3/Core/FlatLutEntry.cxx +++ b/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 @@ -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((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(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(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(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); + 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(mNchBins) * mRadBins * mEtaBins * mPtBins * sizeof(lutEntry_t); + const size_t expectedSize = sizeof(lutHeader_t) + static_cast(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(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(&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()) { + 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; diff --git a/ALICE3/Core/FlatLutEntry.h b/ALICE3/Core/FlatLutEntry.h index 7786523170b..08ab38dcab5 100644 --- a/ALICE3/Core/FlatLutEntry.h +++ b/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(nch_bin) * (mRadBins * mEtaBins * mPtBins) + static_cast(rad_bin) * (mEtaBins * mPtBins) + static_cast(eta_bin) * mPtBins + static_cast(pt_bin); + } + /** * @brief Get LUT entry by bin indices (view) */ diff --git a/ALICE3/Core/FlatTrackSmearer.cxx b/ALICE3/Core/FlatTrackSmearer.cxx index d90c8a8997b..668473369be 100644 --- a/ALICE3/Core/FlatTrackSmearer.cxx +++ b/ALICE3/Core/FlatTrackSmearer.cxx @@ -27,6 +27,7 @@ #include #include #include +#include namespace o2::delphes { @@ -139,7 +140,7 @@ bool TrackSmearer::adoptTable(int pdg, const uint8_t* buffer, size_t size, bool 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 @@ bool TrackSmearer::viewTable(int pdg, const uint8_t* buffer, size_t size, bool f 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; diff --git a/ALICE3/Macros/drawTrackSmearer.C b/ALICE3/Macros/drawTrackSmearer.C new file mode 100644 index 00000000000..3e43a1842c7 --- /dev/null +++ b/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 +#include +#include +#include +#include + +namespace +{ +constexpr int kElectronPdg = static_cast(PDG_t::kElectron); +constexpr int kMuonPdg = static_cast(PDG_t::kMuonMinus); +constexpr int kPionPdg = static_cast(PDG_t::kPiPlus); +constexpr int kProtonPdg = static_cast(PDG_t::kProton); +constexpr int kKaonPdg = static_cast(PDG_t::kKPlus); + +const std::vector> 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>& 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"); + } +}