TPC: Cluster filter for PbPb 2023

wiechula · davidrohr · commit b6432870f650 · 2025-03-06T18:14:25.000+01:00
diff --git a/GPU/GPUTracking/Debug/GPUTPCClusterFilter.cxx b/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;
 }
diff --git a/GPU/GPUTracking/Debug/GPUTPCClusterFilter.h b/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
 
diff --git a/GPU/GPUTracking/Definitions/GPUSettingsList.h b/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")
diff --git a/GPU/GPUTracking/Global/GPUChainTracking.cxx b/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) {
diff --git a/GPU/GPUTracking/Global/GPUChainTrackingClusterizer.cxx b/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)]);
diff --git a/GPU/GPUTracking/Global/GPUChainTrackingCompression.cxx b/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)) {
diff --git a/GPU/GPUTracking/Global/GPUChainTrackingDebugAndProfiling.cxx b/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;
           }

Original file line number	Diff line number	Diff line change
`@@ -294,16 +294,12 @@ bool GPUChainTracking::ValidateSettings()`
`294`	`294`	`return false;`
`295`	`295`	`}`
`296`	`296`	`}`
`297`		`- if ((GetRecoSteps() & GPUDataTypes::RecoStep::TPCDecompression) && GetProcessingSettings().tpcApplyCFCutsAtDecoding && !GetProcessingSettings().tpcUseOldCPUDecoding) {`
`298`		`- GPUError("tpcApplyCFCutsAtDecoding currently requires tpcUseOldCPUDecoding");`
`299`		`- return false;`
`300`		`- }`
`301`	`297`	`if ((GetRecoStepsGPU() & GPUDataTypes::RecoStep::TPCCompression) && !(GetRecoStepsGPU() & GPUDataTypes::RecoStep::TPCCompression) && (GetProcessingSettings().tpcCompressionGatherMode == 1 \|\| GetProcessingSettings().tpcCompressionGatherMode == 3)) {`
`302`	`298`	`GPUError("Invalid tpcCompressionGatherMode for compression on CPU");`
`303`	`299`	`return false;`
`304`	`300`	`}`
`305`		`- if (GetProcessingSettings().tpcApplyDebugClusterFilter == 1 && (GetRecoStepsGPU() & GPUDataTypes::RecoStep::TPCClusterFinding \|\| GetProcessingSettings().delayedOutput \|\| GetProcessingSettings().runMC)) {`
`306`		`- GPUError("tpcApplyDebugClusterFilter cannot be used with GPU clusterization or with delayedOutput for GPU or with MC labels");`
	`301`	`+ if (GetProcessingSettings().tpcApplyClusterFilterOnCPU > 0 && (GetRecoStepsGPU() & GPUDataTypes::RecoStep::TPCClusterFinding \|\| GetProcessingSettings().delayedOutput \|\| GetProcessingSettings().runMC)) {`
	`302`	`+ GPUError("tpcApplyClusterFilterOnCPU cannot be used with GPU clusterization or with delayedOutput for GPU or with MC labels");`
`307`	`303`	`return false;`
`308`	`304`	`}`
`309`	`305`	`if (GetRecoSteps() & RecoStep::TRDTracking) {`
`@@ -815,7 +811,7 @@ int32_t GPUChainTracking::RunChainFinalize()`
`815`	`811`
`816`	`812`	`PrintDebugOutput();`
`817`	`813`
`818`		`- //PrintMemoryRelations();`
	`814`	`+ // PrintMemoryRelations();`
`819`	`815`
`820`	`816`	`if (GetProcessingSettings().eventDisplay) {`
`821`	`817`	`if (!mDisplayRunning) {`
Original file line number	Diff line number	Diff line change
`@@ -295,7 +295,8 @@ void GPUChainTracking::SanityCheck()`
`295`	`295`
`296`	`296`	`void GPUChainTracking::RunTPCClusterFilter(o2::tpc::ClusterNativeAccess* clusters, std::function<o2::tpc::ClusterNative*(size_t)> allocator, bool applyClusterCuts)`
`297`	`297`	`{`
`298`		`- GPUTPCClusterFilter clusterFilter(*clusters);`
	`298`	`+ const uint8_t filterType = GetProcessingSettings().tpcApplyClusterFilterOnCPU;`
	`299`	`+ GPUTPCClusterFilter clusterFilter(*clusters, filterType);`
`299`	`300`	`o2::tpc::ClusterNative* outputBuffer = nullptr;`
`300`	`301`	`for (int32_t iPhase = 0; iPhase < 2; iPhase++) {`
`301`	`302`	`uint32_t countTotal = 0;`
`@@ -312,7 +313,7 @@ void GPUChainTracking::RunTPCClusterFilter(o2::tpc::ClusterNativeAccess* cluster`
`312`	`313`	`if (param().tpcCutTimeBin > 0) {`
`313`	`314`	`keep = keep && cl.getTime() < param().tpcCutTimeBin;`
`314`	`315`	`}`
`315`		`- keep = keep && (!GetProcessingSettings().tpcApplyDebugClusterFilter \|\| clusterFilter.filter(iSector, iRow, cl));`
	`316`	`+ keep = keep && (!filterType \|\| clusterFilter.filter(iSector, iRow, cl));`
`316`	`317`	`if (iPhase && keep) {`
`317`	`318`	`outputBuffer[countTotal] = cl;`
`318`	`319`	`}`