Add TPC cluster selector helper for tracking studies

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()

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 <vector>
+#include <array>
+#include <Rtypes.h>
+
+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<o2::tpc::ClusterNative[]> 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<int> cllist; // optional list
+    int nfnd = clSel.findClustersEntries(sector, row, tbmin, tbmax, padmin, padmax, tpcClusterIdxStruct, &cllist);
+    for (int i=0;i<nfnd;i++) {
+      const auto& cln = tpcClusterIdxStruct.clusters[sector][row][cllist[i]]; /...  direct indices!
+    }
+   */
+
+ public:
+  void clear()
+  {
+    for (auto& s : mSectors)
+      s.clear();
+  }
+  size_t getIndex(int sec, int row, uint32_t icl) const { return mSectors[sec].rows[row][icl]; }
+
+  std::pair<int, int> 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<int>* 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<std::vector<uint16_t>, NRows> rows;
+    void clear()
+    {
+      for (auto& r : rows)
+        r.clear();
+    }
+  };
+
+  static constexpr int NSectors = 36;
+  std::array<Sector, NSectors> mSectors{};
+  int mNThreads = 1;
+
+  ClassDefNV(TPCClusSelector, 1);
+};
+
+} // namespace o2::tpc
+
+#endif

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

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 <numeric>
+#ifdef WITH_OPENMP
+#include <omp.h>
+#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<int, int> 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<int>* 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(); });
+    }
+  }
+}