AliceO2Group · davidrohr · Dec 3, 2020 · Sep 29, 2020 · Dec 2, 2020
@@ -352,6 +352,7 @@
 
 #define GPUCA_THREAD_COUNT_SCAN 512 // TODO: WARNING!!! Must not be GPUTYPE-dependent right now! // TODO: Fix!
 
+#define GPUCA_LB_GPUTPCCFCheckPadBaseline GPUCA_WARP_SIZE
 #define GPUCA_LB_GPUTPCCFChargeMapFiller_fillIndexMap GPUCA_LB_CLUSTER_FINDER
 #define GPUCA_LB_GPUTPCCFChargeMapFiller_fillFromDigits GPUCA_LB_CLUSTER_FINDER
 #define GPUCA_LB_GPUTPCCFChargeMapFiller_findFragmentStart GPUCA_LB_CLUSTER_FINDER

@@ -77,6 +77,7 @@ using namespace GPUCA_NAMESPACE::gpu;
 #include "GPUTPCCFClusterizer.cxx"
 #include "GPUTPCCFDeconvolution.cxx"
 #include "GPUTPCCFMCLabelFlattener.cxx"
+#include "GPUTPCCFCheckPadBaseline.cxx"
 #include "GPUTPCCFDecodeZS.cxx"
 #include "GPUTPCCFGather.cxx"
 

@@ -70,6 +70,7 @@ GPUCA_KRNL_LB((GPUTPCCompressionGatherKernels, buffered64   ), (simple), (), ())
 GPUCA_KRNL_LB((GPUTPCCompressionGatherKernels, buffered128  ), (simple), (), ())
 GPUCA_KRNL_LB((GPUTPCCompressionGatherKernels, multiBlock   ), (simple), (), ())
 
+GPUCA_KRNL_LB((GPUTPCCFCheckPadBaseline                     ), (single), (), ())
 GPUCA_KRNL_LB((GPUTPCCFChargeMapFiller,  fillIndexMap       ), (single), (), ())
 GPUCA_KRNL_LB((GPUTPCCFChargeMapFiller,  fillFromDigits     ), (single), (), ())
 GPUCA_KRNL_LB((GPUTPCCFChargeMapFiller,  findFragmentStart  ), (single), (, char setPositions), (, setPositions))

@@ -50,6 +50,9 @@ AddOptionRTC(tpcTubeMaxSize2, float, 2.5f * 2.5f, "", 0, "Square of max tube siz
 AddOptionRTC(trdMinTrackPt, float, .5f, "", 0, "Min Pt for tracks to be propagated through the TRD")
 AddOptionRTC(trdMaxChi2, float, 15.f, "", 0, "Max chi2 for TRD tracklets to be matched to a track")
 AddOptionRTC(trdPenaltyChi2, float, 12.f, "", 0, "Chi2 penalty for no available TRD tracklet (effective chi2 cut value)")
+AddOptionRTC(noisyPadsQuickCheck, unsigned char, 1, "", 0, "Only check first fragment for noisy pads instead of all fragments (when test is enabled).")
+AddOptionRTC(maxTimeBinAboveThresholdIn1000Bin, unsigned short, 500, "", 0, "Except pad from cluster finding if total number of charges in a fragment is above this baseline (disable = 0)")
+AddOptionRTC(maxConsecTimeBinAboveThreshold, unsigned short, 200, "", 0, "Except pad from cluster finding if number of consecutive charges in a fragment is above this baseline (disable = 0)")
 AddOptionRTC(tpcCFqmaxCutoff, unsigned char, 3, "", 0, "Cluster Finder rejects cluster with qmax below this threshold")
 AddOptionRTC(tpcCFqtotCutoff, unsigned char, 0, "", 0, "Cluster Finder rejects cluster with qtot below this threshold")
 AddOptionRTC(tpcCFinnerThreshold, unsigned char, 0, "", 0, "Cluster Finder extends cluster if inner charge above this threshold")

@@ -163,6 +163,7 @@ if(ALIGPU_BUILD_TYPE STREQUAL "O2" OR CONFIG_O2_EXTENSIONS)
         TPCClusterFinder/GPUTPCClusterFinder.cxx
         TPCClusterFinder/ClusterAccumulator.cxx
         TPCClusterFinder/MCLabelAccumulator.cxx
+        TPCClusterFinder/GPUTPCCFCheckPadBaseline.cxx
         TPCClusterFinder/GPUTPCCFStreamCompaction.cxx
         TPCClusterFinder/GPUTPCCFChargeMapFiller.cxx
         TPCClusterFinder/GPUTPCCFPeakFinder.cxx

@@ -1162,6 +1162,10 @@ int GPUChainTracking::RunTPCClusterizer(bool synchronizeOutput)
         using PeakMapType = decltype(*clustererShadow.mPpeakMap);
         runKernel<GPUMemClean16>(GetGridAutoStep(lane, RecoStep::TPCClusterFinding), krnlRunRangeNone, {}, clustererShadow.mPchargeMap, TPCMapMemoryLayout<ChargeMapType>::items() * sizeof(ChargeMapType));
         runKernel<GPUMemClean16>(GetGridAutoStep(lane, RecoStep::TPCClusterFinding), krnlRunRangeNone, {}, clustererShadow.mPpeakMap, TPCMapMemoryLayout<PeakMapType>::items() * sizeof(PeakMapType));
+        if (fragment.index == 0) {
+          using HasLostBaselineType = decltype(*clustererShadow.mPpadHasLostBaseline);
+          runKernel<GPUMemClean16>(GetGridAutoStep(lane, RecoStep::TPCClusterFinding), krnlRunRangeNone, {}, clustererShadow.mPpadHasLostBaseline, TPC_PADS_IN_SECTOR * sizeof(HasLostBaselineType));
+        }
         DoDebugAndDump(RecoStep::TPCClusterFinding, 0, clusterer, &GPUTPCClusterFinder::DumpChargeMap, *mDebugFile, "Zeroed Charges");
 
         if (mIOPtrs.tpcZS && mCFContext->nPagesSector[iSlice]) {
@@ -1221,6 +1225,14 @@ int GPUChainTracking::RunTPCClusterizer(bool synchronizeOutput)
           runKernel<GPUTPCCFChargeMapFiller, GPUTPCCFChargeMapFiller::fillIndexMap>(GetGrid(clusterer.mPmemory->counters.nDigitsInFragment, lane, GPUReconstruction::krnlDeviceType::CPU), {iSlice}, {});
         }
 
+        bool checkForNoisyPads = (rec()->GetParam().rec.maxTimeBinAboveThresholdIn1000Bin > 0) || (rec()->GetParam().rec.maxConsecTimeBinAboveThreshold > 0);
+        checkForNoisyPads &= (rec()->GetParam().rec.noisyPadsQuickCheck ? fragment.index == 0 : true);
+
+        if (checkForNoisyPads) {
+          int nBlocks = TPC_PADS_IN_SECTOR / GPUTPCCFCheckPadBaseline::getPadsPerBlock(doGPU);
+          runKernel<GPUTPCCFCheckPadBaseline>(GetGridBlk(nBlocks, lane), {iSlice}, {});
+        }
+
         runKernel<GPUTPCCFPeakFinder>(GetGrid(clusterer.mPmemory->counters.nPositions, lane), {iSlice}, {});
         DoDebugAndDump(RecoStep::TPCClusterFinding, 0, clusterer, &GPUTPCClusterFinder::DumpPeaks, *mDebugFile);
 
@@ -1310,7 +1322,7 @@ int GPUChainTracking::RunTPCClusterizer(bool synchronizeOutput)
         transferRunning[lane] = 1;
       }
 
-      if (not propagateMCLabels) {
+      if (not propagateMCLabels || clusterer.mPmemory->counters.nClusters == 0) {
         continue;
       }
 

@@ -72,6 +72,21 @@ struct CfFragment {
     return (hasBacklog ? t < OverlapTimebins : false) || (hasFuture ? t >= (length - OverlapTimebins) : false);
   }
 
+  GPUdi() tpccf::TPCFragmentTime lengthWithoutOverlap() const
+  {
+    return length - (hasBacklog ? OverlapTimebins : 0) - (hasFuture ? OverlapTimebins : 0);
+  }
+
+  GPUdi() tpccf::TPCFragmentTime firstNonOverlapTimeBin() const
+  {
+    return (hasBacklog ? OverlapTimebins : 0);
+  }
+
+  GPUdi() tpccf::TPCFragmentTime lastNonOverlapTimeBin() const
+  {
+    return length - (hasFuture ? OverlapTimebins : 0);
+  }
+
   GPUdi() tpccf::TPCFragmentTime toLocal(tpccf::TPCTime t) const
   {
     return t - first();

@@ -0,0 +1,80 @@
+// Copyright CERN and copyright holders of ALICE O2. This software is
+// distributed under the terms of the GNU General Public License v3 (GPL
+// Version 3), copied verbatim in the file "COPYING".
+//
+// See http://alice-o2.web.cern.ch/license for full licensing information.
+//
+// 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.
+
+/// \file GPUTPCCFCheckPadBaseline.h
+/// \author Felix Weiglhofer
+
+#include "GPUTPCCFCheckPadBaseline.h"
+#include "Array2D.h"
+#include "PackedCharge.h"
+#include "clusterFinderDefs.h"
+
+using namespace GPUCA_NAMESPACE::gpu;
+using namespace GPUCA_NAMESPACE::gpu::tpccf;
+
+template <>
+GPUd() void GPUTPCCFCheckPadBaseline::Thread<0>(int nBlocks, int nThreads, int iBlock, int iThread, GPUSharedMemory& smem, processorType& clusterer)
+{
+  static_assert(TPC_MAX_FRAGMENT_LEN % NumOfCachedTimebins == 0);
+
+  Array2D<PackedCharge> chargeMap(reinterpret_cast<PackedCharge*>(clusterer.mPchargeMap));
+
+  int totalCharges = 0;
+  int consecCharges = 0;
+  int maxConsecCharges = 0;
+
+  int localPadId = iThread / NumOfCachedTimebins;
+  int localTimeBin = iThread % NumOfCachedTimebins;
+  bool handlePad = localTimeBin == 0;
+  int basePad = iBlock * PadsPerBlock;
+
+  CfFragment& fragment = clusterer.mPmemory->fragment;
+
+  ChargePos basePos = padToChargePos(basePad + localPadId, clusterer);
+
+  for (tpccf::TPCFragmentTime t = localTimeBin + fragment.firstNonOverlapTimeBin(); t < fragment.lastNonOverlapTimeBin(); t += NumOfCachedTimebins) {
+    ChargePos pos = basePos.delta({0, t});
+    smem.charges[localPadId][localTimeBin] = (pos.valid()) ? chargeMap[pos].unpack() : 0;
+    GPUbarrierWarp();
+    if (handlePad) {
+      for (int i = 0; i < NumOfCachedTimebins; i++) {
+        Charge q = smem.charges[localPadId][i];
+        totalCharges += (q > 0);
+        consecCharges = (q > 0) ? consecCharges + 1 : 0;
+        maxConsecCharges = CAMath::Max(consecCharges, maxConsecCharges);
+      }
+    }
+  }
+
+  GPUbarrierWarp();
+
+  if (handlePad) {
+    int totalChargesBaseline = clusterer.Param().rec.maxTimeBinAboveThresholdIn1000Bin * fragment.lengthWithoutOverlap() / 1000;
+    int consecChargesBaseline = clusterer.Param().rec.maxConsecTimeBinAboveThreshold;
+    bool hasLostBaseline = (totalChargesBaseline > 0 && totalCharges >= totalChargesBaseline) || (consecChargesBaseline > 0 && maxConsecCharges >= consecChargesBaseline);
+    clusterer.mPpadHasLostBaseline[basePad + localPadId] |= hasLostBaseline;
+  }
+}
+
+GPUd() ChargePos GPUTPCCFCheckPadBaseline::padToChargePos(int pad, const GPUTPCClusterFinder& clusterer)
+{
+  const GPUTPCGeometry& geo = clusterer.Param().tpcGeometry;
+
+  int padOffset = 0;
+  for (Row r = 0; r < TPC_NUM_OF_ROWS; r++) {
+    int padInRow = pad - padOffset;
+    if (0 <= padInRow && padInRow < geo.NPads(r)) {
+      return ChargePos{r, Pad(padInRow), 0};
+    }
+    padOffset += geo.NPads(r);
+  }
+
+  return ChargePos{0, 0, INVALID_TIME_BIN};
+}
@@ -0,0 +1,75 @@
+// Copyright CERN and copyright holders of ALICE O2. This software is
+// distributed under the terms of the GNU General Public License v3 (GPL
+// Version 3), copied verbatim in the file "COPYING".
+//
+// See http://alice-o2.web.cern.ch/license for full licensing information.
+//
+// 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.
+
+/// \file GPUTPCCFCheckPadBaseline.h
+/// \author Felix Weiglhofer
+
+#ifndef O2_GPU_GPU_TPC_CF_CHECK_PAD_BASELINE_H
+#define O2_GPU_GPU_TPC_CF_CHECK_PAD_BASELINE_H
+
+#include "GPUGeneralKernels.h"
+#include "GPUConstantMem.h"
+
+#include "clusterFinderDefs.h"
+
+namespace GPUCA_NAMESPACE::gpu
+{
+
+class GPUTPCCFCheckPadBaseline : public GPUKernelTemplate
+{
+
+ private:
+  // Only use these constants on device side...
+  // Use getPadsPerBlock() for host side
+  enum {
+    PadsPerBlockGPU = 4, // Number of pads in a single cache line
+    PadsPerBlockCPU = 1,
+#ifdef GPUCA_GPUCODE
+    PadsPerBlock = PadsPerBlockGPU,
+#else
+    PadsPerBlock = PadsPerBlockCPU,
+#endif
+    NumOfCachedTimebins = GPUCA_GET_THREAD_COUNT(GPUCA_LB_GPUTPCCFCheckPadBaseline) / PadsPerBlock,
+  };
+
+ public:
+  struct GPUSharedMemory {
+    tpccf::Charge charges[PadsPerBlock][NumOfCachedTimebins];
+  };
+
+#ifdef HAVE_O2HEADERS
+  typedef GPUTPCClusterFinder processorType;
+  GPUhdi() static processorType* Processor(GPUConstantMem& processors)
+  {
+    return processors.tpcClusterer;
+  }
+#endif
+
+  GPUhdi() CONSTEXPR static GPUDataTypes::RecoStep GetRecoStep()
+  {
+    return GPUDataTypes::RecoStep::TPCClusterFinding;
+  }
+
+  // Use this to get num of pads per block on host side. Can't use constant there.
+  static int getPadsPerBlock(bool isGPU)
+  {
+    return (isGPU) ? PadsPerBlockGPU : PadsPerBlockCPU;
+  }
+
+  template <int iKernel = defaultKernel>
+  GPUd() static void Thread(int nBlocks, int nThreads, int iBlock, int iThread, GPUSharedMemory& smem, processorType& clusterer);
+
+ private:
+  GPUd() static ChargePos padToChargePos(int pad, const GPUTPCClusterFinder&);
+};
+
+} // namespace GPUCA_NAMESPACE::gpu
+
+#endif
@@ -26,7 +26,7 @@ GPUdii() void GPUTPCCFPeakFinder::Thread<0>(int nBlocks, int nThreads, int iBloc
 {
   Array2D<PackedCharge> chargeMap(reinterpret_cast<PackedCharge*>(clusterer.mPchargeMap));
   Array2D<uchar> isPeakMap(clusterer.mPpeakMap);
-  findPeaksImpl(get_num_groups(0), get_local_size(0), get_group_id(0), get_local_id(0), smem, chargeMap, clusterer.mPpositions, clusterer.mPmemory->counters.nPositions, clusterer.Param().rec, clusterer.mPisPeak, isPeakMap);
+  findPeaksImpl(get_num_groups(0), get_local_size(0), get_group_id(0), get_local_id(0), smem, chargeMap, clusterer.mPpadHasLostBaseline, clusterer.mPpositions, clusterer.mPmemory->counters.nPositions, clusterer.Param().rec, *clusterer.GetConstantMem()->calibObjects.tpcPadGain, clusterer.mPisPeak, isPeakMap);
 }
 
 GPUdii() bool GPUTPCCFPeakFinder::isPeak(
@@ -91,9 +91,11 @@ GPUdii() bool GPUTPCCFPeakFinder::isPeak(
 
 GPUd() void GPUTPCCFPeakFinder::findPeaksImpl(int nBlocks, int nThreads, int iBlock, int iThread, GPUSharedMemory& smem,
                                               const Array2D<PackedCharge>& chargeMap,
+                                              const uchar* padHasLostBaseline,
                                               const ChargePos* positions,
                                               SizeT digitnum,
                                               const GPUSettingsRec& calib,
+                                              const TPCPadGainCalib& gainCorrection, // Only used for globalPad() function
                                               uchar* isPeakPredicate,
                                               Array2D<uchar>& peakMap)
 {
@@ -105,8 +107,10 @@ GPUd() void GPUTPCCFPeakFinder::findPeaksImpl(int nBlocks, int nThreads, int iBl
   ChargePos pos = positions[CAMath::Min(idx, (SizeT)(digitnum - 1))];
   Charge charge = pos.valid() ? chargeMap[pos].unpack() : Charge(0);
 
-  uchar peak;
-  peak = isPeak(smem, charge, pos, SCRATCH_PAD_SEARCH_N, chargeMap, calib, smem.posBcast, smem.buf);
+  bool hasLostBaseline = padHasLostBaseline[gainCorrection.globalPad(pos.row(), pos.pad())];
+  charge = (hasLostBaseline) ? 0.f : charge;
+
+  uchar peak = isPeak(smem, charge, pos, SCRATCH_PAD_SEARCH_N, chargeMap, calib, smem.posBcast, smem.buf);
 
   // Exit early if dummy. See comment above.
   bool iamDummy = (idx >= digitnum);

@@ -52,7 +52,7 @@ class GPUTPCCFPeakFinder : public GPUKernelTemplate
   GPUd() static void Thread(int nBlocks, int nThreads, int iBlock, int iThread, GPUSharedMemory& smem, processorType& clusterer, Args... args);
 
  private:
-  static GPUd() void findPeaksImpl(int, int, int, int, GPUSharedMemory&, const Array2D<PackedCharge>&, const ChargePos*, tpccf::SizeT, const GPUSettingsRec&, uchar*, Array2D<uchar>&);
+  static GPUd() void findPeaksImpl(int, int, int, int, GPUSharedMemory&, const Array2D<PackedCharge>&, const uchar*, const ChargePos*, tpccf::SizeT, const GPUSettingsRec&, const TPCPadGainCalib&, uchar*, Array2D<uchar>&);
 
   static GPUd() bool isPeak(GPUSharedMemory&, tpccf::Charge, const ChargePos&, ushort, const Array2D<PackedCharge>&, const GPUSettingsRec&, ChargePos*, PackedCharge*);
 };

@@ -75,6 +75,7 @@ void* GPUTPCClusterFinder::SetPointersOutput(void* mem)
 
 void* GPUTPCClusterFinder::SetPointersScratch(void* mem)
 {
+  computePointerWithAlignment(mem, mPpadHasLostBaseline, TPC_PADS_IN_SECTOR);
   computePointerWithAlignment(mem, mPpositions, mNMaxDigitsFragment);
   computePointerWithAlignment(mem, mPpeakPositions, mNMaxPeaks);
   computePointerWithAlignment(mem, mPfilteredPeakPositions, mNMaxClusters);

@@ -48,6 +48,8 @@ struct TPCPadGainCalib;
 
 struct ChargePos;
 
+class GPUTPCGeometry;
+
 class GPUTPCClusterFinder : public GPUProcessor
 {
  public:
@@ -96,6 +98,7 @@ class GPUTPCClusterFinder : public GPUProcessor
   unsigned char* mPzs = nullptr;
   ZSOffset* mPzsOffsets = nullptr;
   MinMaxCN* mMinMaxCN = nullptr;
+  unsigned char* mPpadHasLostBaseline = nullptr;
   tpc::Digit* mPdigits = nullptr; // input digits, only set if ZS is skipped
   ChargePos* mPpositions = nullptr;
   ChargePos* mPpeakPositions = nullptr;
@@ -136,6 +139,7 @@ class GPUTPCClusterFinder : public GPUProcessor
   short mOutputId = -1;
 
   GPUdi() float getGainCorrection(tpccf::Row, tpccf::Pad) const;
+  GPUdi() const GPUTPCGeometry* getGeometry() const;
 
 #ifndef GPUCA_GPUCODE
   void DumpDigits(std::ostream& out);

@@ -22,6 +22,7 @@
 #include "GPUTPCCFStreamCompaction.h"
 #include "GPUTPCCFClusterizer.h"
 #include "GPUTPCCFMCLabelFlattener.h"
+#include "GPUTPCCFCheckPadBaseline.h"
 #include "GPUTPCCFDecodeZS.h"
 #include "GPUTPCCFGather.h"
 

@@ -58,6 +58,11 @@ struct TPCPadGainCalib {
     return mGainCorrection[sector].get(globalPad(row, pad));
   }
 
+  GPUdi() unsigned short globalPad(tpccf::Row row, tpccf::Pad pad) const
+  {
+    return mPadOffsetPerRow[row] + pad;
+  }
+
  private:
   template <typename T = unsigned short>
   class SectorPadGainCorrection
@@ -121,10 +126,6 @@ struct TPCPadGainCalib {
   unsigned short mPadOffsetPerRow[TPC_NUM_OF_ROWS];
   SectorPadGainCorrection<unsigned short> mGainCorrection[TPC_SECTORS];
 
-  GPUdi() unsigned short globalPad(tpccf::Row row, tpccf::Pad pad) const
-  {
-    return mPadOffsetPerRow[row] + pad;
-  }
 };
 
 } // namespace GPUCA_NAMESPACE::gpu