Skip to content

Commit efff2f7

Browse files
shahor02shahor02
authored
Add TPC cluster selector helper for tracking studies (#14308)
1 parent 0386f65 commit efff2f7

File tree

4 files changed

+217
-0
lines changed

4 files changed

+217
-0
lines changed

‎Detectors/GlobalTrackingWorkflow/study/CMakeLists.txt‎

Lines changed: 7 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -12,6 +12,7 @@
1212
#add_compile_options(-O0 -g -fPIC)
1313

1414
o2_add_library(GlobalTrackingStudy
15+
TARGETVARNAME targetName
1516
SOURCES src/TPCTrackStudy.cxx
1617
src/TrackingStudy.cxx
1718
src/SVStudy.cxx
@@ -23,6 +24,7 @@ o2_add_library(GlobalTrackingStudy
2324
src/TrackInfoExt.cxx
2425
src/TrackMCStudyConfig.cxx
2526
src/TrackMCStudyTypes.cxx
27+
src/TPCClusSelector.cxx
2628
PUBLIC_LINK_LIBRARIES O2::GlobalTracking
2729
O2::GlobalTrackingWorkflowReaders
2830
O2::GlobalTrackingWorkflowHelpers
@@ -73,3 +75,8 @@ o2_add_executable(dump-workfow
7375
COMPONENT_NAME bc-tracks
7476
SOURCES src/track-dump-workflow.cxx
7577
PUBLIC_LINK_LIBRARIES O2::GlobalTrackingStudy)
78+
79+
if (OpenMP_CXX_FOUND)
80+
target_compile_definitions(${targetName} PRIVATE WITH_OPENMP)
81+
target_link_libraries(${targetName} PRIVATE OpenMP::OpenMP_CXX)
82+
endif()

‎Detectors/GlobalTrackingWorkflow/study/include/GlobalTrackingStudy/TPCClusSelector.h‎

Lines changed: 92 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,92 @@
1+
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2+
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3+
// All rights not expressly granted are reserved.
4+
//
5+
// This software is distributed under the terms of the GNU General Public
6+
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7+
//
8+
// In applying this license CERN does not waive the privileges and immunities
9+
// granted to it by virtue of its status as an Intergovernmental Organization
10+
// or submit itself to any jurisdiction.
11+
12+
// helper class for TPC clusters selection
13+
14+
#ifndef ALICEO2_TPCCLUSSELECTOR_H
15+
#define ALICEO2_TPCCLUSSELECTOR_H
16+
17+
#include <vector>
18+
#include <array>
19+
#include <Rtypes.h>
20+
21+
namespace o2::tpc
22+
{
23+
class ClusterNativeAccess;
24+
25+
class TPCClusSelector
26+
{
27+
// helper to select TPC cluster matching to certain timebin and optionally pads range
28+
// example of usage:
29+
/*
30+
TPCClusSelector clSel;
31+
o2::tpc::ClusterNativeHelper::Reader tcpClusterReader;
32+
tcpClusterReader.init(native_clusters_file.c_str());
33+
o2::tpc::ClusterNativeAccess tpcClusterIdxStruct;
34+
std::unique_ptr<o2::tpc::ClusterNative[]> tpcClusterBuffer; ///< buffer for clusters in tpcClusterIdxStruct
35+
o2::tpc::ClusterNativeHelper::ConstMCLabelContainerViewWithBuffer tpcClusterMCBuffer; ///< buffer for mc labels
36+
37+
tcpClusterReader.read(iTF);
38+
tcpClusterReader.fillIndex(tpcClusterIdxStruct, tpcClusterBuffer, tpcClusterMCBuffer);
39+
40+
clSel.fill(tpcClusterIdxStruct); // Create sorted index
41+
// to get i-th cluster in orderer timebins:
42+
const auto& clus = tpcClusterIdxStruct.clusters[sector][row][ clSel.getIndex(sector, row, i)];
43+
44+
// to get sorted indices range of clusters in the tbmin:tbmax range
45+
auto rng = clSel.findClustersRange(sector, row, tbmin, tbmax, tpcClusterIdxStruct);
46+
if (rng.first>rng.second) { // nothing is found }
47+
const auto& cln = tpcClusterIdxStruct.clusters[sector][row][clSel.getIndex(sector, row, rng.first )]; /...
48+
49+
// to get number of clusters in tbmin:tbmax, padmin:padmax range (and optionally get the list)
50+
std::vector<int> cllist; // optional list
51+
int nfnd = clSel.findClustersEntries(sector, row, tbmin, tbmax, padmin, padmax, tpcClusterIdxStruct, &cllist);
52+
for (int i=0;i<nfnd;i++) {
53+
const auto& cln = tpcClusterIdxStruct.clusters[sector][row][cllist[i]]; /... direct indices!
54+
}
55+
*/
56+
57+
public:
58+
void clear()
59+
{
60+
for (auto& s : mSectors)
61+
s.clear();
62+
}
63+
size_t getIndex(int sec, int row, uint32_t icl) const { return mSectors[sec].rows[row][icl]; }
64+
65+
std::pair<int, int> findClustersRange(int sec, int row, float tbmin, float tbmax, const o2::tpc::ClusterNativeAccess& tpcClusterIdxStruct);
66+
int findClustersEntries(int sec, int row, float tbmin, float tbmax, float padmin, float padmax, const o2::tpc::ClusterNativeAccess& tpcClusterIdxStruct, std::vector<int>* clIDDirect = nullptr);
67+
void fill(const o2::tpc::ClusterNativeAccess& tpcClusterIdxStruct);
68+
69+
int getNThreads() const { return mNThreads; }
70+
void setNThreads(int n);
71+
72+
private:
73+
struct Sector {
74+
static constexpr int NRows = 152;
75+
std::array<std::vector<uint16_t>, NRows> rows;
76+
void clear()
77+
{
78+
for (auto& r : rows)
79+
r.clear();
80+
}
81+
};
82+
83+
static constexpr int NSectors = 36;
84+
std::array<Sector, NSectors> mSectors{};
85+
int mNThreads = 1;
86+
87+
ClassDefNV(TPCClusSelector, 1);
88+
};
89+
90+
} // namespace o2::tpc
91+
92+
#endif

‎Detectors/GlobalTrackingWorkflow/study/src/GlobalTrackingStudyLinkDef.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -38,5 +38,6 @@
3838
#pragma link C++ class std::vector < o2::trackstudy::ClResTPCCont> + ;
3939
#pragma link C++ class o2::trackstudy::TrackPairInfo + ;
4040
#pragma link C++ class std::vector < o2::trackstudy::TrackPairInfo> + ;
41+
#pragma ling C++ class o2::tpc::TPCClusSelector + ;
4142

4243
#endif

‎Detectors/GlobalTrackingWorkflow/study/src/TPCClusSelector.cxx‎

Lines changed: 117 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,117 @@
1+
// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2+
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3+
// All rights not expressly granted are reserved.
4+
//
5+
// This software is distributed under the terms of the GNU General Public
6+
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7+
//
8+
// In applying this license CERN does not waive the privileges and immunities
9+
// granted to it by virtue of its status as an Intergovernmental Organization
10+
// or submit itself to any jurisdiction.
11+
12+
// helper class for TPC clusters selection
13+
#include "GlobalTrackingStudy/TPCClusSelector.h"
14+
#include "DataFormatsTPC/ClusterNativeHelper.h"
15+
#include "Framework/Logger.h"
16+
#include <numeric>
17+
#ifdef WITH_OPENMP
18+
#include <omp.h>
19+
#endif
20+
21+
using namespace o2::tpc;
22+
23+
void TPCClusSelector::setNThreads(int n)
24+
{
25+
#ifndef WITH_OPENMP
26+
if (n > 1) {
27+
LOGP(warn, "No OpenMP");
28+
}
29+
n = 1;
30+
#endif
31+
mNThreads = n;
32+
}
33+
34+
std::pair<int, int> TPCClusSelector::findClustersRange(int sec, int row, float tbmin, float tbmax, const o2::tpc::ClusterNativeAccess& tpcClusterIdxStruct)
35+
{
36+
// find sorted indices of clusters in the [tbmin:tbmax] range, if not found, return {-1,-2}
37+
const auto& vidx = mSectors[sec].rows[row];
38+
const auto* clarr = tpcClusterIdxStruct.clusters[sec][row];
39+
// use binary search to find 1st cluster with time >= tb
40+
int ncl = vidx.size(), left = 0, right = ncl;
41+
while (left < right) {
42+
int mid = left + (right - left) / 2;
43+
if (clarr[vidx[mid]].getTime() < tbmin) {
44+
left = mid + 1;
45+
} else {
46+
right = mid;
47+
}
48+
}
49+
if (left == ncl || clarr[vidx[left]].getTime() > tbmax) {
50+
return {-1, -2}; // all clusters have time < tbmin or no clusters in the range [tbmin:tbmax]
51+
}
52+
int idmin = left, idmax = left, idtst = idmin;
53+
// look at smaller times
54+
while (++idtst < ncl && clarr[vidx[idtst]].getTime() <= tbmax) {
55+
idmax = idtst;
56+
}
57+
return {idmin, idmax};
58+
}
59+
60+
int TPCClusSelector::findClustersEntries(int sec, int row, float tbmin, float tbmax, float padmin, float padmax, const o2::tpc::ClusterNativeAccess& tpcClusterIdxStruct, std::vector<int>* clIDDirect)
61+
{
62+
// find direct cluster indices for tbmin:tbmas / padmin/padmax range, fill clIDDirect vector if provided
63+
const auto& vidx = mSectors[sec].rows[row];
64+
const auto* clarr = tpcClusterIdxStruct.clusters[sec][row];
65+
// use binary search to find 1st cluster with time >= tb
66+
int ncl = vidx.size(), left = 0, right = ncl;
67+
if (clIDDirect) {
68+
clIDDirect->clear();
69+
}
70+
while (left < right) {
71+
int mid = left + (right - left) / 2;
72+
if (clarr[vidx[mid]].getTime() < tbmin) {
73+
left = mid + 1;
74+
} else {
75+
right = mid;
76+
}
77+
}
78+
if (left == ncl || clarr[vidx[left]].getTime() > tbmax) {
79+
return 0; // all clusters have time < tbmin or no clusters in the range [tbmin:tbmax]
80+
}
81+
int nclf = 0;
82+
while (left < ncl) {
83+
const auto& cl = clarr[vidx[left]];
84+
if (cl.getTime() > tbmax) {
85+
break;
86+
}
87+
if (cl.getPad() >= padmin && cl.getPad() <= padmax) {
88+
nclf++;
89+
if (clIDDirect) {
90+
clIDDirect->push_back(vidx[left]);
91+
}
92+
}
93+
}
94+
return nclf;
95+
}
96+
97+
void TPCClusSelector::fill(const o2::tpc::ClusterNativeAccess& tpcClusterIdxStruct)
98+
{
99+
for (int is = 0; is < NSectors; is++) {
100+
auto& sect = mSectors[is];
101+
#ifdef WITH_OPENMP
102+
#pragma omp parallel for schedule(dynamic) num_threads(mNThreads)
103+
#endif
104+
for (int ir = 0; ir < Sector::NRows; ir++) {
105+
size_t ncl = tpcClusterIdxStruct.nClusters[is][ir];
106+
if (ncl >= 0xffff) {
107+
LOGP(error, "Row {} of sector {} has {} clusters, truncating to {}", ir, is, ncl, int(0xffff));
108+
ncl = 0xffff;
109+
}
110+
auto& rowidx = sect.rows[ir];
111+
rowidx.resize(ncl);
112+
std::iota(rowidx.begin(), rowidx.end(), 0);
113+
const auto* clus = tpcClusterIdxStruct.clusters[is][ir]; // C array of clusters
114+
std::sort(rowidx.begin(), rowidx.end(), [&](size_t a, size_t b) { return clus[a].getTime() < clus[b].getTime(); });
115+
}
116+
}
117+
}

0 commit comments

Comments
 (0)