From a85e282bee0e29473634915f19e4c66cc82032d2 Mon Sep 17 00:00:00 2001 From: shahoian Date: Wed, 21 May 2025 16:35:02 +0200 Subject: [PATCH] Add TPC cluster selector helper for tracking studies --- .../study/CMakeLists.txt | 7 ++ .../GlobalTrackingStudy/TPCClusSelector.h | 92 ++++++++++++++ .../study/src/GlobalTrackingStudyLinkDef.h | 1 + .../study/src/TPCClusSelector.cxx | 117 ++++++++++++++++++ 4 files changed, 217 insertions(+) create mode 100644 Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/TPCClusSelector.h create mode 100644 Detectors/GlobalTrackingWorkflow/study/src/TPCClusSelector.cxx diff --git a/Detectors/GlobalTrackingWorkflow/study/CMakeLists.txt b/Detectors/GlobalTrackingWorkflow/study/CMakeLists.txt index 398e7eb215f2e..776d3946283c3 100644 --- a/Detectors/GlobalTrackingWorkflow/study/CMakeLists.txt +++ b/Detectors/GlobalTrackingWorkflow/study/CMakeLists.txt @@ -12,6 +12,7 @@ #add_compile_options(-O0 -g -fPIC) o2_add_library(GlobalTrackingStudy + TARGETVARNAME targetName SOURCES src/TPCTrackStudy.cxx src/TrackingStudy.cxx src/SVStudy.cxx @@ -23,6 +24,7 @@ o2_add_library(GlobalTrackingStudy src/TrackInfoExt.cxx src/TrackMCStudyConfig.cxx src/TrackMCStudyTypes.cxx + src/TPCClusSelector.cxx PUBLIC_LINK_LIBRARIES O2::GlobalTracking O2::GlobalTrackingWorkflowReaders O2::GlobalTrackingWorkflowHelpers @@ -73,3 +75,8 @@ o2_add_executable(dump-workfow COMPONENT_NAME bc-tracks SOURCES src/track-dump-workflow.cxx PUBLIC_LINK_LIBRARIES O2::GlobalTrackingStudy) + +if (OpenMP_CXX_FOUND) + target_compile_definitions(${targetName} PRIVATE WITH_OPENMP) + target_link_libraries(${targetName} PRIVATE OpenMP::OpenMP_CXX) +endif() diff --git a/Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/TPCClusSelector.h b/Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/TPCClusSelector.h new file mode 100644 index 0000000000000..c1765558458c2 --- /dev/null +++ b/Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/TPCClusSelector.h @@ -0,0 +1,92 @@ +// 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. + +// helper class for TPC clusters selection + +#ifndef ALICEO2_TPCCLUSSELECTOR_H +#define ALICEO2_TPCCLUSSELECTOR_H + +#include +#include +#include + +namespace o2::tpc +{ +class ClusterNativeAccess; + +class TPCClusSelector +{ + // helper to select TPC cluster matching to certain timebin and optionally pads range + // example of usage: + /* + TPCClusSelector clSel; + o2::tpc::ClusterNativeHelper::Reader tcpClusterReader; + tcpClusterReader.init(native_clusters_file.c_str()); + o2::tpc::ClusterNativeAccess tpcClusterIdxStruct; + std::unique_ptr tpcClusterBuffer; ///< buffer for clusters in tpcClusterIdxStruct + o2::tpc::ClusterNativeHelper::ConstMCLabelContainerViewWithBuffer tpcClusterMCBuffer; ///< buffer for mc labels + + tcpClusterReader.read(iTF); + tcpClusterReader.fillIndex(tpcClusterIdxStruct, tpcClusterBuffer, tpcClusterMCBuffer); + + clSel.fill(tpcClusterIdxStruct); // Create sorted index + // to get i-th cluster in orderer timebins: + const auto& clus = tpcClusterIdxStruct.clusters[sector][row][ clSel.getIndex(sector, row, i)]; + + // to get sorted indices range of clusters in the tbmin:tbmax range + auto rng = clSel.findClustersRange(sector, row, tbmin, tbmax, tpcClusterIdxStruct); + if (rng.first>rng.second) { // nothing is found } + const auto& cln = tpcClusterIdxStruct.clusters[sector][row][clSel.getIndex(sector, row, rng.first )]; /... + + // to get number of clusters in tbmin:tbmax, padmin:padmax range (and optionally get the list) + std::vector cllist; // optional list + int nfnd = clSel.findClustersEntries(sector, row, tbmin, tbmax, padmin, padmax, tpcClusterIdxStruct, &cllist); + for (int i=0;i findClustersRange(int sec, int row, float tbmin, float tbmax, const o2::tpc::ClusterNativeAccess& tpcClusterIdxStruct); + int findClustersEntries(int sec, int row, float tbmin, float tbmax, float padmin, float padmax, const o2::tpc::ClusterNativeAccess& tpcClusterIdxStruct, std::vector* clIDDirect = nullptr); + void fill(const o2::tpc::ClusterNativeAccess& tpcClusterIdxStruct); + + int getNThreads() const { return mNThreads; } + void setNThreads(int n); + + private: + struct Sector { + static constexpr int NRows = 152; + std::array, NRows> rows; + void clear() + { + for (auto& r : rows) + r.clear(); + } + }; + + static constexpr int NSectors = 36; + std::array mSectors{}; + int mNThreads = 1; + + ClassDefNV(TPCClusSelector, 1); +}; + +} // namespace o2::tpc + +#endif diff --git a/Detectors/GlobalTrackingWorkflow/study/src/GlobalTrackingStudyLinkDef.h b/Detectors/GlobalTrackingWorkflow/study/src/GlobalTrackingStudyLinkDef.h index f666132c9c1cf..f0d3e7d4d0b4e 100644 --- a/Detectors/GlobalTrackingWorkflow/study/src/GlobalTrackingStudyLinkDef.h +++ b/Detectors/GlobalTrackingWorkflow/study/src/GlobalTrackingStudyLinkDef.h @@ -38,5 +38,6 @@ #pragma link C++ class std::vector < o2::trackstudy::ClResTPCCont> + ; #pragma link C++ class o2::trackstudy::TrackPairInfo + ; #pragma link C++ class std::vector < o2::trackstudy::TrackPairInfo> + ; +#pragma ling C++ class o2::tpc::TPCClusSelector + ; #endif diff --git a/Detectors/GlobalTrackingWorkflow/study/src/TPCClusSelector.cxx b/Detectors/GlobalTrackingWorkflow/study/src/TPCClusSelector.cxx new file mode 100644 index 0000000000000..e5b28fb0fd62b --- /dev/null +++ b/Detectors/GlobalTrackingWorkflow/study/src/TPCClusSelector.cxx @@ -0,0 +1,117 @@ +// 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. + +// helper class for TPC clusters selection +#include "GlobalTrackingStudy/TPCClusSelector.h" +#include "DataFormatsTPC/ClusterNativeHelper.h" +#include "Framework/Logger.h" +#include +#ifdef WITH_OPENMP +#include +#endif + +using namespace o2::tpc; + +void TPCClusSelector::setNThreads(int n) +{ +#ifndef WITH_OPENMP + if (n > 1) { + LOGP(warn, "No OpenMP"); + } + n = 1; +#endif + mNThreads = n; +} + +std::pair TPCClusSelector::findClustersRange(int sec, int row, float tbmin, float tbmax, const o2::tpc::ClusterNativeAccess& tpcClusterIdxStruct) +{ + // find sorted indices of clusters in the [tbmin:tbmax] range, if not found, return {-1,-2} + const auto& vidx = mSectors[sec].rows[row]; + const auto* clarr = tpcClusterIdxStruct.clusters[sec][row]; + // use binary search to find 1st cluster with time >= tb + int ncl = vidx.size(), left = 0, right = ncl; + while (left < right) { + int mid = left + (right - left) / 2; + if (clarr[vidx[mid]].getTime() < tbmin) { + left = mid + 1; + } else { + right = mid; + } + } + if (left == ncl || clarr[vidx[left]].getTime() > tbmax) { + return {-1, -2}; // all clusters have time < tbmin or no clusters in the range [tbmin:tbmax] + } + int idmin = left, idmax = left, idtst = idmin; + // look at smaller times + while (++idtst < ncl && clarr[vidx[idtst]].getTime() <= tbmax) { + idmax = idtst; + } + return {idmin, idmax}; +} + +int TPCClusSelector::findClustersEntries(int sec, int row, float tbmin, float tbmax, float padmin, float padmax, const o2::tpc::ClusterNativeAccess& tpcClusterIdxStruct, std::vector* clIDDirect) +{ + // find direct cluster indices for tbmin:tbmas / padmin/padmax range, fill clIDDirect vector if provided + const auto& vidx = mSectors[sec].rows[row]; + const auto* clarr = tpcClusterIdxStruct.clusters[sec][row]; + // use binary search to find 1st cluster with time >= tb + int ncl = vidx.size(), left = 0, right = ncl; + if (clIDDirect) { + clIDDirect->clear(); + } + while (left < right) { + int mid = left + (right - left) / 2; + if (clarr[vidx[mid]].getTime() < tbmin) { + left = mid + 1; + } else { + right = mid; + } + } + if (left == ncl || clarr[vidx[left]].getTime() > tbmax) { + return 0; // all clusters have time < tbmin or no clusters in the range [tbmin:tbmax] + } + int nclf = 0; + while (left < ncl) { + const auto& cl = clarr[vidx[left]]; + if (cl.getTime() > tbmax) { + break; + } + if (cl.getPad() >= padmin && cl.getPad() <= padmax) { + nclf++; + if (clIDDirect) { + clIDDirect->push_back(vidx[left]); + } + } + } + return nclf; +} + +void TPCClusSelector::fill(const o2::tpc::ClusterNativeAccess& tpcClusterIdxStruct) +{ + for (int is = 0; is < NSectors; is++) { + auto& sect = mSectors[is]; +#ifdef WITH_OPENMP +#pragma omp parallel for schedule(dynamic) num_threads(mNThreads) +#endif + for (int ir = 0; ir < Sector::NRows; ir++) { + size_t ncl = tpcClusterIdxStruct.nClusters[is][ir]; + if (ncl >= 0xffff) { + LOGP(error, "Row {} of sector {} has {} clusters, truncating to {}", ir, is, ncl, int(0xffff)); + ncl = 0xffff; + } + auto& rowidx = sect.rows[ir]; + rowidx.resize(ncl); + std::iota(rowidx.begin(), rowidx.end(), 0); + const auto* clus = tpcClusterIdxStruct.clusters[is][ir]; // C array of clusters + std::sort(rowidx.begin(), rowidx.end(), [&](size_t a, size_t b) { return clus[a].getTime() < clus[b].getTime(); }); + } + } +}