TPC: Cluster filter for PbPb 2023

GPU/GPUTracking/Debug/GPUTPCClusterFilter.cxx

@@ -12,20 +12,62 @@
 /// \file GPUTPCClusterFilter.cxx
 /// \author David Rohr
 
+#include "GPUCommonLogger.h"
 #include "GPUTPCClusterFilter.h"
 #include "DataFormatsTPC/ClusterNative.h"
 
 using namespace o2::gpu;
 
-GPUTPCClusterFilter::GPUTPCClusterFilter(const o2::tpc::ClusterNativeAccess& clusters)
+GPUTPCClusterFilter::GPUTPCClusterFilter(const o2::tpc::ClusterNativeAccess& clusters, uint8_t filterType)
+  : mFilterType(filterType)
 {
-  // Could initialize private variables based on the clusters here
+  if (filterType == 1) {
+    // Custom filter settings go here
+
+  } else if (filterType == 2) {
+    // PbPb23 filter
+    mClusterStats = std::make_unique<std::vector<int>[]>(MaxStacks);
+    static bool called = false;
+    if (!called) {
+      LOGP(info, "GPUTPCClusterFilter called for PbPb 2023 settings");
+      called = true;
+    }
+
+    for (uint32_t iSector = 0; iSector < GPUCA_NSECTORS; iSector++) {
+      for (uint32_t iRow = 0; iRow < GPUCA_ROW_COUNT; iRow++) {
+        const uint32_t globalStack = getGlobalStack(iSector, iRow);
+        mClusterStats[globalStack].resize(MaxTimeBin);
+
+        for (uint32_t k = 0; k < clusters.nClusters[iSector][iRow]; k++) {
+          const o2::tpc::ClusterNative& cl = clusters.clusters[iSector][iRow][k];
+          const int clTime = static_cast<int>(cl.getTime());
+          const float clQmax = cl.getQmax();
+
+          if (clQmax < 12) {
+            if (clTime >= static_cast<int>(mClusterStats[globalStack].size())) {
+              mClusterStats[globalStack].resize(mClusterStats[globalStack].size() + 445);
+            }
+            ++mClusterStats[globalStack][clTime];
+          }
+        }
+      }
+    }
+  }
 }
 
 bool GPUTPCClusterFilter::filter(uint32_t sector, uint32_t row, o2::tpc::ClusterNative& cl)
 {
   // Return true to keep the cluster, false to drop it.
   // May change cluster properties by modifying the cl reference.
   // Note that this function might be called multiple times for the same cluster, in which case the final modified cl reference goes into the output clusters.
+  if (mFilterType == 2) {
+    const uint32_t globalStack = getGlobalStack(sector, row);
+    const int clTime = static_cast<int>(cl.getTime());
+    const float clQmax = cl.getQmax();
+    if ((mClusterStats[globalStack][clTime] > 40 && clQmax < 12) || (mClusterStats[globalStack][clTime] > 200)) {
+      return false;
+    }
+  }
+
   return true;
 }

GPU/GPUTracking/Debug/GPUTPCClusterFilter.h

@@ -15,7 +15,10 @@
 #ifndef GPUTPCCLUSTERFILTER_H
 #define GPUTPCCLUSTERFILTER_H
 
+#include <memory>
 #include <cstdint>
+#include <vector>
+#include "GPUDefConstantsAndSettings.h"
 
 namespace o2::tpc
 {
@@ -28,8 +31,29 @@ namespace o2::gpu
 class GPUTPCClusterFilter
 {
  public:
-  GPUTPCClusterFilter(const o2::tpc::ClusterNativeAccess& clusters);
+  GPUTPCClusterFilter(const o2::tpc::ClusterNativeAccess& clusters, uint8_t filterType);
   bool filter(uint32_t sector, uint32_t row, o2::tpc::ClusterNative& cl);
+
+ private:
+  static constexpr uint32_t MaxTimeBin = 14256;
+  static constexpr uint32_t MaxStacks = GPUCA_NSECTORS * 4;
+  uint8_t mFilterType = 0; //< 0: off, 1: custom, 2: PbPb23
+
+  std::unique_ptr<std::vector<int>[]> mClusterStats; //< Number of clusters per stack and time bin
+
+  uint32_t getGlobalStack(uint32_t sector, uint32_t row) const
+  {
+    int stack = 3;
+    if (row < 63) {
+      stack = 0;
+    } else if (row < 97) {
+      stack = 1;
+    } else if (row < 127) {
+      stack = 2;
+    }
+
+    return sector * 4 + stack;
+  };
 };
 } // namespace o2::gpu
 

GPU/GPUTracking/Definitions/GPUSettingsList.h

@@ -291,7 +291,7 @@ AddOption(tpcDownscaledEdx, uint8_t, 0, "", 0, "If != 0, downscale dEdx processi
 AddOption(tpcMaxAttachedClustersPerSectorRow, uint32_t, 51000, "", 0, "Maximum number of TPC attached clusters which can be decoded per SectorRow")
 AddOption(tpcUseOldCPUDecoding, bool, false, "", 0, "Enable old CPU-based TPC decoding")
 AddOption(tpcApplyCFCutsAtDecoding, bool, false, "", 0, "Apply cluster cuts from clusterization during decoding of compressed clusters")
-AddOption(tpcApplyDebugClusterFilter, bool, false, "", 0, "Apply custom cluster filter of GPUTPCClusterFilter class")
+AddOption(tpcApplyClusterFilterOnCPU, uint8_t, 0, "", 0, "Apply custom cluster filter of GPUTPCClusterFilter class, 0: off, 1: debug, 2: PbPb23")
 AddOption(RTCcacheFolder, std::string, "./rtccache/", "", 0, "Folder in which the cache file is stored")
 AddOption(RTCprependCommand, std::string, "", "", 0, "Prepend RTC compilation commands by this string")
 AddOption(RTCoverrideArchitecture, std::string, "", "", 0, "Override arhcitecture part of RTC compilation command line")

GPU/GPUTracking/Global/GPUChainTracking.cxx

@@ -294,16 +294,12 @@ bool GPUChainTracking::ValidateSettings()
       return false;
     }
   }
-  if ((GetRecoSteps() & GPUDataTypes::RecoStep::TPCDecompression) && GetProcessingSettings().tpcApplyCFCutsAtDecoding && !GetProcessingSettings().tpcUseOldCPUDecoding) {
-    GPUError("tpcApplyCFCutsAtDecoding currently requires tpcUseOldCPUDecoding");
-    return false;
-  }
   if ((GetRecoStepsGPU() & GPUDataTypes::RecoStep::TPCCompression) && !(GetRecoStepsGPU() & GPUDataTypes::RecoStep::TPCCompression) && (GetProcessingSettings().tpcCompressionGatherMode == 1 || GetProcessingSettings().tpcCompressionGatherMode == 3)) {
     GPUError("Invalid tpcCompressionGatherMode for compression on CPU");
     return false;
   }
-  if (GetProcessingSettings().tpcApplyDebugClusterFilter == 1 && (GetRecoStepsGPU() & GPUDataTypes::RecoStep::TPCClusterFinding || GetProcessingSettings().delayedOutput || GetProcessingSettings().runMC)) {
-    GPUError("tpcApplyDebugClusterFilter cannot be used with GPU clusterization or with delayedOutput for GPU or with MC labels");
+  if (GetProcessingSettings().tpcApplyClusterFilterOnCPU > 0 && (GetRecoStepsGPU() & GPUDataTypes::RecoStep::TPCClusterFinding || GetProcessingSettings().delayedOutput || GetProcessingSettings().runMC)) {
+    GPUError("tpcApplyClusterFilterOnCPU cannot be used with GPU clusterization or with delayedOutput for GPU or with MC labels");
     return false;
   }
   if (GetRecoSteps() & RecoStep::TRDTracking) {
@@ -815,7 +811,7 @@ int32_t GPUChainTracking::RunChainFinalize()
 
   PrintDebugOutput();
 
-  //PrintMemoryRelations();
+  // PrintMemoryRelations();
 
   if (GetProcessingSettings().eventDisplay) {
     if (!mDisplayRunning) {

GPU/GPUTracking/Global/GPUChainTrackingClusterizer.cxx

@@ -635,7 +635,7 @@ int32_t GPUChainTracking::RunTPCClusterizer(bool synchronizeOutput)
     if (mWaitForFinalInputs) {
       GPUFatal("Cannot use waitForFinalInput callback without delayed output");
     }
-    if (!GetProcessingSettings().tpcApplyDebugClusterFilter) {
+    if (!GetProcessingSettings().tpcApplyClusterFilterOnCPU) {
       AllocateRegisteredMemory(mInputsHost->mResourceClusterNativeOutput, mSubOutputControls[GPUTrackingOutputs::getIndex(&GPUTrackingOutputs::clustersNative)]);
       tmpNativeClusters = mInputsHost->mPclusterNativeOutput;
     } else {
@@ -1021,7 +1021,7 @@ int32_t GPUChainTracking::RunTPCClusterizer(bool synchronizeOutput)
     tmpNativeAccess->clustersMCTruth = mcLabelsConstView;
     tmpNativeAccess->setOffsetPtrs();
     mIOPtrs.clustersNative = tmpNativeAccess;
-    if (GetProcessingSettings().tpcApplyDebugClusterFilter) {
+    if (GetProcessingSettings().tpcApplyClusterFilterOnCPU) {
       auto allocator = [this, &tmpNativeClusters](size_t size) {
         this->mInputsHost->mNClusterNative = size;
         this->AllocateRegisteredMemory(this->mInputsHost->mResourceClusterNativeOutput, this->mSubOutputControls[GPUTrackingOutputs::getIndex(&GPUTrackingOutputs::clustersNative)]);

GPU/GPUTracking/Global/GPUChainTrackingCompression.cxx

@@ -201,6 +201,11 @@ int32_t GPUChainTracking::RunTPCCompression()
 
 int32_t GPUChainTracking::RunTPCDecompression()
 {
+  const bool runFiltering = GetProcessingSettings().tpcApplyCFCutsAtDecoding || (GetProcessingSettings().tpcApplyClusterFilterOnCPU > 0) || (param().tpcCutTimeBin > 0);
+  if (runFiltering && !GetProcessingSettings().tpcUseOldCPUDecoding) {
+    GPUFatal("tpcApplyCFCutsAtDecoding, tpcApplyClusterFilterOnCPU and tpcCutTimeBin currently require tpcUseOldCPUDecoding");
+  }
+
   if (GetProcessingSettings().tpcUseOldCPUDecoding) {
     const auto& threadContext = GetThreadContext();
     TPCClusterDecompressor decomp;
@@ -214,7 +219,6 @@ int32_t GPUChainTracking::RunTPCDecompression()
       return ((tmpBuffer = std::make_unique<ClusterNative[]>(size))).get();
     };
     auto& decompressTimer = getTimer<TPCClusterDecompressor>("TPCDecompression", 0);
-    bool runFiltering = GetProcessingSettings().tpcApplyCFCutsAtDecoding;
     auto allocatorUse = runFiltering ? std::function<ClusterNative*(size_t)>{allocatorTmp} : std::function<ClusterNative*(size_t)>{allocatorFinal};
     decompressTimer.Start();
     if (decomp.decompress(mIOPtrs.tpcCompressedClusters, *mClusterNativeAccess, allocatorUse, param(), GetProcessingSettings().deterministicGPUReconstruction)) {

GPU/GPUTracking/Global/GPUChainTrackingDebugAndProfiling.cxx

@@ -295,7 +295,8 @@ void GPUChainTracking::SanityCheck()
 
 void GPUChainTracking::RunTPCClusterFilter(o2::tpc::ClusterNativeAccess* clusters, std::function<o2::tpc::ClusterNative*(size_t)> allocator, bool applyClusterCuts)
 {
-  GPUTPCClusterFilter clusterFilter(*clusters);
+  const uint8_t filterType = GetProcessingSettings().tpcApplyClusterFilterOnCPU;
+  GPUTPCClusterFilter clusterFilter(*clusters, filterType);
   o2::tpc::ClusterNative* outputBuffer = nullptr;
   for (int32_t iPhase = 0; iPhase < 2; iPhase++) {
     uint32_t countTotal = 0;
@@ -312,7 +313,7 @@ void GPUChainTracking::RunTPCClusterFilter(o2::tpc::ClusterNativeAccess* cluster
           if (param().tpcCutTimeBin > 0) {
             keep = keep && cl.getTime() < param().tpcCutTimeBin;
           }
-          keep = keep && (!GetProcessingSettings().tpcApplyDebugClusterFilter || clusterFilter.filter(iSector, iRow, cl));
+          keep = keep && (!filterType || clusterFilter.filter(iSector, iRow, cl));
           if (iPhase && keep) {
             outputBuffer[countTotal] = cl;
           }