GPUTPCCFCheckPadBaseline.cxx

Go to the documentation of this file.

// 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.
 
 
#include "GPUTPCCFCheckPadBaseline.h"
#include "CfArray2D.h"
#include "PackedCharge.h"
#include "clusterFinderDefs.h"
#include "DataFormatsTPC/ClusterNative.h"
#include "GPUCommonAlgorithm.h"
 
#ifndef GPUCA_GPUCODE
#include "MCLabelAccumulator.h"
#include "utils/VcShim.h"
#include <vector>
#endif
 
#if 0
#define DPRINT(...) printf(__VA_ARGS__)
#define DPRINTB(...) \
  if (iThread == 0)  \
  printf(__VA_ARGS__)
#define DPRINTB_IF(test, ...) \
  if (iThread == 0 && (test)) \
  printf(__VA_ARGS__)
#else
#define DPRINT(...) ((void)0)
#define DPRINTB(...) ((void)0)
#define DPRINTB_IF(test, ...) ((void)0)
#endif
 
using namespace o2::gpu;
using namespace o2::gpu::tpccf;
 
using Kernel = GPUTPCCFCheckPadBaseline;
 
static GPUdi() HIPTailDescriptor* GetHIPTails(GPUTPCClusterFinder& clusterer, int32_t row)
{
  // HIP TAILS: indexing starts at 1, so 0 index indicates no connection
  return clusterer.mPhipTailsByRow + row * GPUTPCCFHIPClusterizer::MaxHIPTailsPerRow;
}
 
static GPUdi() Charge UpdateHIPTailFilter(Charge filteredCharge, Charge charge, Charge alpha)
{
  return filteredCharge + alpha * (charge - filteredCharge);
}
 
static GPUdi() float HIPTailTimeVariance(const HIPTailDescriptor& tail)
{
  const float length = tail.tailEnd > tail.tailStart ? float(tail.tailEnd - tail.tailStart) : 1.f;
  return (length * length - 1.f) * (1.f / 12.f);
}
 
// Collect tails marked for closing across the workgroup using a prefix scan,
// then cooperatively zero the charge map entries for each closed tail.
// Caller must set acc.activeHIPTail.end before calling if the tail is open.
static GPUdi() uint16_t CloseHIPTails(
  Kernel::GPUSharedMemory& smem,
  GPUTPCClusterFinder& clusterer,
  int32_t iThread, int32_t nThreads,
  int16_t iPadHandle,
  CfChargePos basePos,
  CfArray2D<PackedCharge>& chargeMap,
  Kernel::PadChargeAccu& acc,
  bool shouldCloseTail)
{
  const uint32_t row = basePos.row();
  const uint16_t nClosedTails = work_group_count(shouldCloseTail);
 
  auto* nHIPTails = clusterer.mPnHIPTails;
  auto* hipTails = GetHIPTails(clusterer, row);
 
  if (nClosedTails > 0) {
    int16_t iClosedTail = work_group_scan_inclusive_add((int16_t)shouldCloseTail) - 1;
    const bool shouldStoreTail = shouldCloseTail && acc.activeHIPTail.Length() > 0;
    uint16_t nStoredTails = work_group_count(shouldStoreTail);
    int16_t iStoredTail = work_group_scan_inclusive_add((int16_t)shouldStoreTail) - 1;
 
    // Use exactly one atomic add per closing call to reduce differences in
    // tail ordering between runs.
    if (nStoredTails > 0) {
      if (iThread == 0) {
        smem.tailStoreBase = CAMath::AtomicAdd(&nHIPTails[row], (uint32_t)nStoredTails);
      }
      GPUbarrier();
    }
    if (shouldCloseTail) {
      smem.tailsClosedPad[iClosedTail] = iPadHandle;
      smem.tailsClosed[iClosedTail] = acc.activeHIPTail;
      smem.tailsClosedStoreIdx[iClosedTail] = GPUTPCCFHIPTailConnector::MaxHIPTailsPerRow;
 
      if (shouldStoreTail) {
        const uint32_t idx = smem.tailStoreBase + iStoredTail + 1;
        smem.tailsClosedStoreIdx[iClosedTail] = idx;
        if (idx < GPUTPCCFHIPTailConnector::MaxHIPTailsPerRow) {
          hipTails[idx] = {0, 0, (uint16_t)iPadHandle,
                           (uint16_t)acc.activeHIPTail.start, (uint16_t)acc.activeHIPTail.end,
                           acc.activeSatStart, acc.activeSatEnd,
                           0.f, 0.f};
        }
      }
 
      acc.tailFilterCharge = 0;
      acc.activeHIPTail.Reset();
      acc.activeSatStart = acc.activeSatEnd = -1;
    }
 
    GPUbarrier();
  }
 
  // TODO: performance improvement -> parallelize this loop across tails
  for (uint16_t iTail = 0; iTail < nClosedTails; iTail++) {
    const auto tailPad = smem.tailsClosedPad[iTail];
    const auto tail = smem.tailsClosed[iTail];
    const uint32_t tailStoreIdx = smem.tailsClosedStoreIdx[iTail];
 
    Charge qTot = 0.f;
    Charge qMax = 0.f;
    for (uint16_t iTime = iThread; iTime < tail.Length(); iTime += nThreads) {
      const int16_t time = tail.start + iTime;
      auto pos = basePos.delta({tailPad, time});
      const Charge q = chargeMap[pos].unpack();
      qTot += q;
      qMax = CAMath::Max(qMax, q);
      chargeMap[pos] = PackedCharge{0};
    }
 
    smem.tailQTotScratch[iThread] = qTot;
    smem.tailQMaxScratch[iThread] = qMax;
    GPUbarrier();
    for (uint16_t active = nThreads; active > 1;) {
      const uint16_t stride = (active + 1) / 2;
      if (iThread < active - stride) {
        smem.tailQTotScratch[iThread] += smem.tailQTotScratch[iThread + stride];
        smem.tailQMaxScratch[iThread] = CAMath::Max(smem.tailQMaxScratch[iThread], smem.tailQMaxScratch[iThread + stride]);
      }
      active = stride;
      GPUbarrier();
    }
 
    if (iThread == 0 && tailStoreIdx < GPUTPCCFHIPTailConnector::MaxHIPTailsPerRow) {
      HIPTailDescriptor& tailDescriptor = hipTails[tailStoreIdx];
      tailDescriptor.qTot = smem.tailQTotScratch[0];
      tailDescriptor.qMax = smem.tailQMaxScratch[0];
    }
  }
 
  return nClosedTails;
}
 
template <bool CheckHIPTrigger, bool CheckHIPTailEnd>
static GPUdi() void ScanCachedCharges(Kernel::GPUSharedMemory& smem, uint16_t timeOffset, uint16_t pad, Charge hipTailThreshold, Charge hipTailFilterAlpha, Kernel::PadChargeAccu& acc)
{
  for (int32_t i = 0; i < Kernel::NumOfCachedTBs; i++) {
    const Charge qs = smem.charges[i][pad];
    const int16_t curTB = timeOffset + i;
 
    acc.totalCharges += qs > 0;
    acc.consecCharges = qs > 0 ? acc.consecCharges + 1 : 0;
    acc.maxConsecCharges = CAMath::Max(acc.consecCharges, acc.maxConsecCharges);
    acc.maxCharge = CAMath::Max<Charge>(qs, acc.maxCharge);
 
    if (qs >= hipTailThreshold) {
      if (acc.aboveThresholdStart < 0) {
        acc.aboveThresholdStart = curTB;
      }
    } else {
      acc.aboveThresholdStart = -1;
    }
 
    if constexpr (CheckHIPTrigger) {
      // Track saturated plateaus. A plateau continuing from the previous chunk always triggers there as well,
      // so tracking them only in chunks with a trigger is sufficient.
      const bool isSaturated = qs >= Charge(Kernel::MaxADC);
      if (isSaturated && !acc.plateauOpen) {
        acc.plateauStart = curTB;
      }
      acc.plateauOpen = isSaturated;
      // Plateau of the active tail continues from a previous chunk
      if (isSaturated && acc.activeSatStart > -1 && acc.activeSatStart == acc.plateauStart) {
        acc.activeSatEnd = curTB;
      }
 
      if (acc.HIPtb < 0 && isSaturated) {
        acc.HIPtb = acc.aboveThresholdStart; // start of rising edge, not first sat TB
        acc.satStart = acc.plateauStart;     // Plateau may have started in the previous chunk if it crosses the chunk boundary
        smem.tails[pad] = {acc.HIPtb, 0};    // Broadcast HIP start TB to neighboring pads / threads
      }
      if (isSaturated && acc.satStart > -1 && acc.satStart == acc.plateauStart) {
        acc.satEnd = curTB;
      }
    }
 
    if constexpr (CheckHIPTailEnd) {
      if (acc.activeHIPTail.IsOpen()) {
        acc.tailFilterCharge = UpdateHIPTailFilter(acc.tailFilterCharge, qs, hipTailFilterAlpha);
        if (acc.tailFilterCharge < hipTailThreshold) {
          acc.activeHIPTail.end = curTB;
        }
      }
    }
  }
}
 
template <>
GPUd() void GPUTPCCFCheckPadBaseline::Thread<0>(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, GPUSharedMemory& smem, processorType& clusterer)
{
#ifdef GPUCA_GPUCODE
  CheckBaselineGPU(nBlocks, nThreads, iBlock, iThread, smem, clusterer);
#else
  CheckBaselineCPU(nBlocks, nThreads, iBlock, iThread, smem, clusterer);
#endif
}
 
// Charges are stored in a 2D array (pad and time) using a tiling layout.
// Tiles are 8 pads x 4 timebins large stored in time-major layout and make up a single cacheline.
//
// This kernel processes one row per block. Threads cooperatively load chunks
// of 4 consecutive time bins for all pads into shared memory. Thread `i` then processes charges for pad `i` in shared memory.
// Blocks require `nextMultipleOf<64>(138 * 4) = 576` threads to process the largest TPC rows with 138 pads correctly.
GPUd() void GPUTPCCFCheckPadBaseline::CheckBaselineGPU(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, GPUSharedMemory& smem, processorType& clusterer)
{
#if defined(GPUCA_GPUCODE) && !defined(GPUCA_DEVICE_SINGLE_THREAD) // Needs 576 threads per block, CPUDevice RTC uses the host-compiled kernel
  static_assert(GPUCA_GET_THREAD_COUNT(GPUCA_LB_GPUTPCCFCheckPadBaseline) == 576);
  if (iBlock >= (int32_t)GPUTPCGeometry::NROWS) {
    return;
  }
 
  const CfFragment& frag = clusterer.mPmemory->frag;
  const bool hipFilterOn = clusterer.Param().rec.tpc.hipTailFilter;
  const Charge hipTailThreshold = clusterer.Param().rec.tpc.hipTailFilterThreshold;
  const Charge hipTailFilterAlpha = clusterer.Param().rec.tpc.hipTailFilterAlpha;
  CfArray2D<PackedCharge> chargeMap(reinterpret_cast<PackedCharge*>(clusterer.mPchargeMap));
 
  constexpr GPUTPCGeometry geo;
 
  const auto iRow = iBlock;
  const auto nPads = geo.NPads(iRow);
  const CfChargePos basePos{(Row)iRow, 0, 0};
 
  PadChargeAccu acc;
 
  const int16_t iPadOffset = iThread % MaxNPadsPerRow;
  const int16_t iTimeOffset = iThread / MaxNPadsPerRow;
  const int16_t iPadHandle = iThread;
  const bool handlePad = iPadHandle < nPads;
 
  if (iPadHandle < MaxNPadsPerRow) {
    smem.tails[iPadHandle] = {-1, -1};
  }
  GPUbarrier();
 
  // Pad filter scans the entire fragments including overlap.
  // Minimal runtime overhead and prevents headaches later on as
  // saturated signal in overlap region can create tails in the next fragment
  // even when cleared in current fragment as they're decoded twice
  const TPCFragmentTime firstTB = 0;
  const TPCFragmentTime lastTB = frag.length;
 
  for (uint16_t t = firstTB; t < lastTB; t += NumOfCachedTBs) {
 
    bool thisThreadHasTrigger = false;
    for (uint16_t tt = 0; tt < NumOfCachedTBs; tt += TimebinsPerCacheline) {
      const TPCFragmentTime iTimeLoad = t + tt + iTimeOffset;
 
      const CfChargePos pos = basePos.delta({iPadOffset, iTimeLoad});
 
      const Charge ql = iTimeLoad < lastTB && iPadOffset < nPads ? chargeMap[pos].unpack() : 0;
      smem.charges[tt + iTimeOffset][iPadOffset] = ql;
 
      thisThreadHasTrigger |= ql >= Charge(MaxADC);
    }
 
    bool hasHIPTrigger = false;
    if (hipFilterOn) {
      hasHIPTrigger = work_group_any(thisThreadHasTrigger);
    } else {
      // Need a barrier here even if HIP filter is disabled
      GPUbarrier();
    }
 
    acc.HIPtb = -1;
    acc.satStart = acc.satEnd = -1;
    if (!hasHIPTrigger) {
      acc.plateauOpen = false; // No saturated TB in this chunk
    }
 
    if (handlePad) {
 
      // TODO: is this really necessary?
      // Why is the old version so much slower, when we just add short branches to the loop???
      if (!hasHIPTrigger) [[likely]] {
        if (!acc.activeHIPTail.IsOpen()) {
          ScanCachedCharges<false, false>(smem, t, iPadHandle, hipTailThreshold, hipTailFilterAlpha, acc);
        } else {
          ScanCachedCharges<false, true>(smem, t, iPadHandle, hipTailThreshold, hipTailFilterAlpha, acc);
        }
      } else {
        if (!acc.activeHIPTail.IsOpen()) {
          ScanCachedCharges<true, false>(smem, t, iPadHandle, hipTailThreshold, hipTailFilterAlpha, acc);
        } else {
          ScanCachedCharges<true, true>(smem, t, iPadHandle, hipTailThreshold, hipTailFilterAlpha, acc);
        }
      }
    }
 
    GPUbarrier();
 
    if (hasHIPTrigger) [[unlikely]] {
 
      DPRINTB("%d: Trigger!\n", iBlock);
 
      if (handlePad && acc.HIPtb < 0) {
 
        // Search neighboring pads for trigger
        for (int16_t i = -SSClusterPadWidth; i < 0; i++) {
          const auto p = iPadHandle + i;
          if (p > -1) {
            acc.HIPtb = CAMath::Max(smem.tails[p].start, acc.HIPtb);
          }
        }
 
        for (int16_t i = 1; i <= SSClusterPadWidth; i++) {
          const auto p = iPadHandle + i;
          if (p < MaxNPadsPerRow) {
            acc.HIPtb = CAMath::Max(smem.tails[p].start, acc.HIPtb);
          }
        }
      }
 
      bool shouldCloseTail = acc.HIPtb > -1 && acc.activeHIPTail.HasValue();
      // End the old tail at the rising edge of the new trigger, also if the tail filter already ended it later in this chunk.
      // Otherwise the old tail would take the saturated samples of the new trigger.
      if (shouldCloseTail && (acc.activeHIPTail.IsOpen() || acc.activeHIPTail.end > acc.HIPtb)) {
        DPRINT("%d: end = %d\n", iThread, acc.HIPtb);
        acc.activeHIPTail.end = acc.HIPtb;
      }
 
      // Tails are closed at the rising edge of the new trigger, which can lie before the saturated plateau of the old tail.
      // The new tail inherits the plateau if it now contains saturated samples, the old one keeps it only if it still does.
      int16_t newSatStart = acc.satStart;
      int16_t newSatEnd = acc.satEnd;
      if (shouldCloseTail && acc.activeSatStart > -1) {
        if (acc.satStart < 0 && acc.activeSatEnd >= acc.activeHIPTail.end) {
          newSatStart = acc.activeSatStart;
          newSatEnd = acc.activeSatEnd;
        }
        if (acc.activeSatStart >= acc.activeHIPTail.end) {
          acc.activeSatStart = acc.activeSatEnd = -1;
        }
      }
 
      CloseHIPTails(smem, clusterer, iThread, nThreads, iPadHandle, basePos, chargeMap, acc, shouldCloseTail);
 
      GPUbarrier();
 
      if (acc.HIPtb > -1) {
        DPRINT("%d: start = %d\n", iThread, acc.HIPtb);
        acc.activeHIPTail.SetOpen(acc.HIPtb);
        acc.tailFilterCharge = Charge(MaxADC);
        acc.activeSatStart = newSatStart;
        acc.activeSatEnd = newSatEnd;
      }
 
      // Clear smem between iterations to prevent stale entries
      if (handlePad) {
        smem.tails[iPadHandle].Reset();
      }
 
      GPUbarrier();
 
    } // if (hasHIPTrigger)
 
  } // for (uint16_t t = firstTB; t < lastTB; t += NumOfCachedTBs)
 
  if (handlePad) {
    updatePadBaseline(basePos.gpad + iPadHandle, clusterer, acc.totalCharges, acc.maxConsecCharges, acc.maxCharge);
  }
 
  // --- Close remaining tails
  const bool shouldCloseTail = acc.activeHIPTail.HasValue();
 
  // Call `work_group_any` here, instead of always counting.
  // This is important as `work_group_count` is a lot slower
  // and has a lot of overhead if no HIPs were found.
  if (work_group_any(shouldCloseTail)) {
    if (shouldCloseTail && acc.activeHIPTail.IsOpen()) {
      acc.activeHIPTail.end = lastTB;
    }
 
    [[maybe_unused]] const uint16_t nClosedTails = CloseHIPTails(smem, clusterer, iThread, nThreads, iPadHandle, basePos, chargeMap, acc, shouldCloseTail);
 
    DPRINTB_IF(nClosedTails > 0, "%d: Close remaining tails (%d)\n", iBlock, nClosedTails);
  }
 
#endif
}
 
GPUd() void GPUTPCCFCheckPadBaseline::CheckBaselineCPU(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, GPUSharedMemory& smem, processorType& clusterer)
{
#ifndef GPUCA_GPUCODE
  if (iBlock >= (int32_t)GPUTPCGeometry::NROWS) {
    return;
  }
 
  constexpr GPUTPCGeometry geo;
  const int32_t row = iBlock;
  const int32_t nPads = geo.NPads(row);
  const int32_t nVecPads = (nPads + PadsPerCacheline - 1) / PadsPerCacheline;
 
  const CfFragment& frag = clusterer.mPmemory->frag;
  const bool hipFilterOn = clusterer.Param().rec.tpc.hipTailFilter;
  const Charge hipTailThreshold = clusterer.Param().rec.tpc.hipTailFilterThreshold;
  const Charge hipTailFilterAlpha = clusterer.Param().rec.tpc.hipTailFilterAlpha;
  auto* nHIPTails = clusterer.mPnHIPTails;
  auto* hipTails = GetHIPTails(clusterer, row);
 
  CfArray2D<PackedCharge> chargeMap(reinterpret_cast<PackedCharge*>(clusterer.mPchargeMap));
 
  using UShort8 = Vc::fixed_size_simd<uint16_t, PadsPerCacheline>;
  using Short8 = Vc::fixed_size_simd<int16_t, PadsPerCacheline>;
  using Charge8 = Vc::fixed_size_simd<float, PadsPerCacheline>;
 
  std::vector<UShort8> totalChargesV(nVecPads, UShort8{Vc::Zero});
  std::vector<UShort8> consecChargesV(nVecPads, UShort8{Vc::Zero});
  std::vector<UShort8> maxConsecChargesV(nVecPads, UShort8{Vc::Zero});
  std::vector<Charge8> maxChargeV(nVecPads, Charge8{Vc::Zero});
 
  std::vector<Short8> localHipTbV(nVecPads, -1);
  std::vector<Short8> broadcastHipTbV(nVecPads, -1);
  std::vector<Short8> localSatStartV(nVecPads, -1); // start of the saturated plateau that triggered in the current chunk, only set for pads with a local trigger
  std::vector<Short8> localSatEndV(nVecPads, -1);   // end of that plateau as far as seen in the current chunk
  std::vector<Short8> plateauStartV(nVecPads, -1);  // first TB of the current / last saturated plateau on this pad
  std::vector<Short8> plateauOpenV(nVecPads, 0);    // 1 while the previous TB was saturated
  std::vector<Short8> aboveThresholdStartV(nVecPads, -1);
  std::vector<Short8> activeHIPTailStartV(nVecPads, -1);
  std::vector<Short8> activeHIPTailEndV(nVecPads, -1);
  std::vector<Short8> activeHIPTailSatStartV(nVecPads, -1); // start of the saturated plateau that triggered the active tail, -1 if inherited from a neighbor
  std::vector<Short8> activeHIPTailSatEndV(nVecPads, -1);   // end of that plateau, extended while the plateau continues into later chunks
  std::vector<Charge8> tailFilterChargeV(nVecPads, Charge8{Vc::Zero});
 
  for (int16_t t = 0; t < frag.length; t += NumOfCachedTBs) {
 
    bool hasAnyTrigger = false;
 
    // Run actual noisy pad filter and look for HIP trigger
    for (int16_t iVecPad = 0; iVecPad < nVecPads; iVecPad++) {
 
      auto totalCharges = totalChargesV[iVecPad];
      auto consecCharges = consecChargesV[iVecPad];
      auto maxConsecCharges = maxConsecChargesV[iVecPad];
      auto maxCharge = maxChargeV[iVecPad];
 
      auto hipTb = Short8(-1);
      auto satStart = Short8(-1);
      auto satEnd = Short8(-1);
      auto plateauStart = plateauStartV[iVecPad];
      auto plateauOpen = plateauOpenV[iVecPad];
      const auto activeHIPTailSatStart = activeHIPTailSatStartV[iVecPad];
      auto activeHIPTailSatEnd = activeHIPTailSatEndV[iVecPad];
      auto aboveThresholdStart = aboveThresholdStartV[iVecPad];
      auto activeHIPTailStart = activeHIPTailStartV[iVecPad];
      auto activeHIPTailEnd = activeHIPTailEndV[iVecPad];
      auto tailFilterCharge = tailFilterChargeV[iVecPad];
 
      const CfChargePos basePos(row, iVecPad * PadsPerCacheline, t);
 
      for (tpccf::TPCFragmentTime localtime = 0; localtime < NumOfCachedTBs; localtime++) {
 
        const uint16_t* packedChargeStart = reinterpret_cast<uint16_t*>(&chargeMap[basePos.delta({0, localtime})]);
        const UShort8 packedCharges = t + localtime < frag.length
                                        ? UShort8{packedChargeStart, Vc::Aligned}
                                        : UShort8{Vc::Zero};
        const auto isCharge = packedCharges != 0;
 
        const auto unpackedCharges = Charge8(packedCharges) / Charge(1 << PackedCharge::DecimalBits);
 
        if (isCharge.isNotEmpty()) {
          totalCharges(isCharge)++;
          consecCharges += 1;
          consecCharges(not isCharge) = 0;
          maxConsecCharges = Vc::max(consecCharges, maxConsecCharges);
 
          // Manually unpack charges to float.
          // Duplicated from PackedCharge::unpack to generate vectorized code:
          //   Charge unpack() const { return Charge(mVal & ChargeMask) / Charge(1 << DecimalBits); }
          // Note that PackedCharge has to cut off the highest 2 bits via ChargeMask as they are used for flags by the cluster finder
          // and are not part of the charge value. We can skip this step because the cluster finder hasn't run yet
          // and thus these bits are guarenteed to be zero.
          maxCharge = Vc::max(maxCharge, unpackedCharges);
 
          const auto aboveRisingEdge = unpackedCharges >= hipTailThreshold;
          const auto startRisingEdge = aboveRisingEdge && aboveThresholdStart < 0;
          aboveThresholdStart(startRisingEdge) = t + localtime;
          aboveThresholdStart(!aboveRisingEdge) = -1;
 
          // Track saturated plateaus across chunk boundaries
          const auto isSaturated = unpackedCharges >= Charge(MaxADC);
          plateauStart(isSaturated && plateauOpen < 1) = t + localtime;
          plateauOpen = 0;
          plateauOpen(isSaturated) = 1;
          // Plateau of the active tail continues from a previous chunk
          activeHIPTailSatEnd(isSaturated && activeHIPTailSatStart > -1 && (activeHIPTailSatStart - plateauStart) == 0) = t + localtime;
 
          const auto hasNewTrigger = hipTb < 0 && isSaturated;
          hipTb(hasNewTrigger) = aboveThresholdStart;
          satStart(hasNewTrigger) = plateauStart; // Plateau may have started in the previous chunk if it crosses the chunk boundary
          satEnd(isSaturated && satStart > -1 && (satStart - plateauStart) == 0) = t + localtime;
          hasAnyTrigger |= hasNewTrigger.isNotEmpty();
        } else {
          consecCharges = 0;
          aboveThresholdStart = -1;
          plateauOpen = 0;
        }
 
        const auto tailOpen = activeHIPTailStart > -1 && activeHIPTailEnd < 0;
        tailFilterCharge(tailOpen) = tailFilterCharge + hipTailFilterAlpha * (unpackedCharges - tailFilterCharge);
        activeHIPTailEnd(tailOpen && tailFilterCharge < hipTailThreshold) = t + localtime;
      } // for (tpccf::TPCFragmentTime localtime = 0; localtime < TimebinsPerCacheline; localtime++)
 
      totalChargesV[iVecPad] = totalCharges;
      consecChargesV[iVecPad] = consecCharges;
      maxConsecChargesV[iVecPad] = maxConsecCharges;
      maxChargeV[iVecPad] = maxCharge;
 
      localHipTbV[iVecPad] = hipTb;
      localSatStartV[iVecPad] = satStart;
      localSatEndV[iVecPad] = satEnd;
      plateauStartV[iVecPad] = plateauStart;
      plateauOpenV[iVecPad] = plateauOpen;
      activeHIPTailSatEndV[iVecPad] = activeHIPTailSatEnd;
      aboveThresholdStartV[iVecPad] = aboveThresholdStart;
      activeHIPTailStartV[iVecPad] = activeHIPTailStart;
      activeHIPTailEndV[iVecPad] = activeHIPTailEnd;
      tailFilterChargeV[iVecPad] = tailFilterCharge;
 
    } // for (int16_t iVecPad = 0; iVecPad < nVecPads; iVecPad++)
 
    if (hasAnyTrigger) {
      broadcastHipTbV = localHipTbV;
    }
 
    // Broadcast trigger times to neighboring pads across the whole
    for (int16_t iVecPad = 0; iVecPad < nVecPads && hasAnyTrigger; iVecPad++) {
 
      const auto hipTb = localHipTbV[iVecPad];
 
      const auto hasHipTrigger = hipTb > -1;
      if (hasHipTrigger.isNotEmpty()) [[unlikely]] {
 
        // TODO: This could be vectorised, but doesn't seem necessary
        for (uint16_t p = 0; p < PadsPerCacheline; p++) {
          if (hasHipTrigger[p]) {
            const int16_t pad = iVecPad * PadsPerCacheline + p;
            const int16_t neighborSt = CAMath::Max(0, pad - SSClusterPadWidth);
            const int16_t neighborEnd = CAMath::Min(nPads, pad + SSClusterPadWidth + 1);
            for (int16_t np = neighborSt; np < neighborEnd; np++) {
              if (np == pad) {
                continue;
              }
              const auto pv = np / PadsPerCacheline;
              const auto pi = np % PadsPerCacheline;
              // GPU keeps a pad's own trigger time; only pads without a local trigger inherit from neighbors.
              if (localHipTbV[pv][pi] < 0) {
                broadcastHipTbV[pv][pi] = CAMath::Max<int16_t>(hipTb[p], broadcastHipTbV[pv][pi]);
              }
            } // for (int16_t np = neighborSt; np < neighborEnd; np++)
          } // if (hasHipTrigger[p]) {
        } // for (uint16_t p = 0; p < PadsPerCacheline; p++)
      } // if (hasHipTrigger.isNotEmpty())
    } // for (int16_t iVecPad = 0; iVecPad < nVecPads; iVecPad++)
 
    // Close old tails for all pads, open new tails in case of overlap
    for (int16_t iVecPad = 0; iVecPad < nVecPads && hasAnyTrigger; iVecPad++) {
 
      auto hipTb = broadcastHipTbV[iVecPad];
      const auto satStart = localSatStartV[iVecPad];
      const auto satEnd = localSatEndV[iVecPad];
      auto aboveThresholdStart = aboveThresholdStartV[iVecPad];
      auto activeHIPTailStart = activeHIPTailStartV[iVecPad];
      auto activeHIPTailEnd = activeHIPTailEndV[iVecPad];
      auto activeHIPTailSatStart = activeHIPTailSatStartV[iVecPad];
      auto activeHIPTailSatEnd = activeHIPTailSatEndV[iVecPad];
      auto tailFilterCharge = tailFilterChargeV[iVecPad];
 
      const auto shouldCloseTail = hipTb > -1 && activeHIPTailStart > -1;
      // End the old tail at the rising edge of the new trigger, also if the tail filter already ended it later in this chunk.
      // Otherwise the old tail would take the saturated samples of the new trigger.
      activeHIPTailEnd(shouldCloseTail && !(activeHIPTailEnd >= 0 && (activeHIPTailEnd - hipTb) < 1)) = hipTb;
 
      // Tails are closed at the rising edge of the new trigger, which can lie before the saturated plateau of the old tail.
      // The new tail inherits the plateau if it now contains saturated samples, the old one keeps it only if it still does.
      auto newSatStart = satStart;
      auto newSatEnd = satEnd;
      const auto inheritSat = shouldCloseTail && satStart < 0 && activeHIPTailSatStart > -1 && (activeHIPTailSatEnd - activeHIPTailEnd) >= 0;
      newSatStart(inheritSat) = activeHIPTailSatStart;
      newSatEnd(inheritSat) = activeHIPTailSatEnd;
      const auto oldLosesSat = shouldCloseTail && activeHIPTailSatStart > -1 && (activeHIPTailSatStart - activeHIPTailEnd) >= 0;
      activeHIPTailSatStart(oldLosesSat) = -1;
      activeHIPTailSatEnd(oldLosesSat) = -1;
 
      // Closing tails will store them to global memory and zero the range
      // So it's enough to disable this part to fully disable the tail filter
      if (hipFilterOn && shouldCloseTail.isNotEmpty()) {
        for (int16_t p = 0; p < PadsPerCacheline; p++) {
          const int16_t pad = iVecPad * PadsPerCacheline + p;
          if (shouldCloseTail[p] && pad < nPads) {
            Charge tailQtot = 0;
            Charge tailQMax = 0;
            for (int16_t tt = activeHIPTailStart[p]; tt < activeHIPTailEnd[p]; tt++) {
              const CfChargePos basePos(row, iVecPad * PadsPerCacheline, 0);
              const auto pos = basePos.delta({p, tt});
              const auto q = chargeMap[pos].unpack();
              tailQtot += q;
              tailQMax = CAMath::Max(tailQMax, q);
              chargeMap[pos] = PackedCharge{0};
            }
 
            if (activeHIPTailEnd[p] > activeHIPTailStart[p]) { // Prune empty tails
              const auto tailIdx = CAMath::AtomicAdd<uint32_t>(&nHIPTails[row], 1) + 1;
              if (tailIdx < GPUTPCCFHIPTailConnector::MaxHIPTailsPerRow) {
                hipTails[tailIdx] = {
                  .iPrev = 0,
                  .iNext = 0,
                  .pad = uint16_t(pad),
                  .tailStart = uint16_t(activeHIPTailStart[p]),
                  .tailEnd = uint16_t(activeHIPTailEnd[p]),
                  .satStart = int16_t(activeHIPTailSatStart[p]),
                  .satEnd = int16_t(activeHIPTailSatEnd[p]),
                  .qTot = tailQtot,
                  .qMax = tailQMax,
                };
              }
            }
 
          } // if (shouldCloseTail[p] && pad < nPads)
        } // for (uint16_t p = 0; p < PadsPerCacheline; p++)
      } // if (shouldCloseThipFilterOn && shouldCloseTail.isNotEmpty())
 
      activeHIPTailStart(hipTb > -1) = hipTb;
      activeHIPTailEnd(hipTb > -1) = -1;
      activeHIPTailSatStart(hipTb > -1) = newSatStart;
      activeHIPTailSatEnd(hipTb > -1) = newSatEnd;
      tailFilterCharge(hipTb > -1) = MaxADC;
 
      aboveThresholdStartV[iVecPad] = aboveThresholdStart;
      activeHIPTailStartV[iVecPad] = activeHIPTailStart;
      activeHIPTailEndV[iVecPad] = activeHIPTailEnd;
      activeHIPTailSatStartV[iVecPad] = activeHIPTailSatStart;
      activeHIPTailSatEndV[iVecPad] = activeHIPTailSatEnd;
      tailFilterChargeV[iVecPad] = tailFilterCharge;
 
    } // for (int32_t iVecPad = 0; iVecPad < nVecPads; iVecPad++)
  } // for (auto t = 0; t < frag.length; t += TimebinsPerCacheline)
 
  // Close old tails for all pads, open new tails in case of overlap
  for (int16_t iVecPad = 0; iVecPad < nVecPads; iVecPad++) {
 
    auto activeHIPTailStart = activeHIPTailStartV[iVecPad];
    auto activeHIPTailEnd = activeHIPTailEndV[iVecPad];
    const auto activeHIPTailSatStart = activeHIPTailSatStartV[iVecPad];
    const auto activeHIPTailSatEnd = activeHIPTailSatEndV[iVecPad];
 
    const auto shouldCloseTail = activeHIPTailStart > -1;
    activeHIPTailEnd(shouldCloseTail && activeHIPTailEnd < 0) = frag.length;
 
    if (hipFilterOn && shouldCloseTail.isNotEmpty()) {
      for (int16_t p = 0; p < PadsPerCacheline; p++) {
        const int16_t pad = iVecPad * PadsPerCacheline + p;
        if (shouldCloseTail[p] && pad < nPads) {
          Charge tailQtot = 0;
          Charge tailQMax = 0;
          for (int16_t tt = activeHIPTailStart[p]; tt < activeHIPTailEnd[p]; tt++) {
            const CfChargePos basePos(row, iVecPad * PadsPerCacheline, 0);
            const auto pos = basePos.delta({p, tt});
            const auto q = chargeMap[pos].unpack();
            tailQtot += q;
            tailQMax = CAMath::Max(tailQMax, q);
            chargeMap[pos] = PackedCharge{0};
          }
 
          if (activeHIPTailEnd[p] > activeHIPTailStart[p]) { // Prune empty tails
            const auto tailIdx = CAMath::AtomicAdd<uint32_t>(&nHIPTails[row], 1) + 1;
            if (tailIdx < GPUTPCCFHIPTailConnector::MaxHIPTailsPerRow) {
              hipTails[tailIdx] = {
                .iPrev = 0,
                .iNext = 0,
                .pad = uint16_t(pad),
                .tailStart = uint16_t(activeHIPTailStart[p]),
                .tailEnd = uint16_t(activeHIPTailEnd[p]),
                .satStart = int16_t(activeHIPTailSatStart[p]),
                .satEnd = int16_t(activeHIPTailSatEnd[p]),
                .qTot = tailQtot,
                .qMax = tailQMax,
              };
            }
          }
 
        } // if (shouldCloseTail[p] && pad < nPads)
      } // for (uint16_t p = 0; p < PadsPerCacheline; p++)
    } // if (hipFilterOn && shouldCloseTail.isNotEmpty())
  } // for (int16_t iVecPad = 0; iVecPad < nVecPads; iVecPad++)
 
  for (int32_t iVecPad = 0; iVecPad < nVecPads; iVecPad++) {
 
    const UShort8 totalCharges = totalChargesV[iVecPad];
    const UShort8 maxConsecCharges = maxConsecChargesV[iVecPad];
    const Charge8 maxCharge = maxChargeV[iVecPad];
 
    const CfChargePos basePos(row, iVecPad * PadsPerCacheline, 0);
 
    for (tpccf::Pad localpad = 0; localpad < PadsPerCacheline; localpad++) {
      updatePadBaseline(basePos.gpad + localpad, clusterer, totalCharges[localpad], maxConsecCharges[localpad], maxCharge[localpad]);
    }
  }
#endif
}
 
GPUd() void GPUTPCCFCheckPadBaseline::updatePadBaseline(int32_t pad, const GPUTPCClusterFinder& clusterer, int32_t totalCharges, int32_t consecCharges, Charge maxCharge)
{
  const CfFragment& frag = clusterer.mPmemory->frag;
  const int32_t totalChargesBaseline = clusterer.Param().rec.tpc.maxTimeBinAboveThresholdIn1000Bin * frag.lengthWithoutOverlap() / 1000;
  const int32_t consecChargesBaseline = clusterer.Param().rec.tpc.maxConsecTimeBinAboveThreshold;
  const uint16_t saturationThreshold = clusterer.Param().rec.tpc.noisyPadSaturationThreshold;
  const bool isNoisy = (!saturationThreshold || maxCharge < saturationThreshold) && ((totalChargesBaseline > 0 && totalCharges >= totalChargesBaseline) || (consecChargesBaseline > 0 && consecCharges >= consecChargesBaseline));
 
  if (isNoisy) {
    clusterer.mPpadIsNoisy[pad] = true;
  }
}
 
// ======== HIP Tail Connector Kernel ========
 
template <>
GPUd() void GPUTPCCFHIPTailConnector::Thread<0>(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, GPUSharedMemory& smem, processorType& clusterer)
{
  if (iBlock >= (int32_t)GPUTPCGeometry::NROWS) {
    return;
  }
  const uint32_t row = iBlock;
 
  const uint32_t nTails = CAMath::Min(clusterer.mPnHIPTails[row], (uint32_t)MaxHIPTailsPerRow - 1);
 
  // HIP TAILS: indexing starts at 1, so 0 index indicates no connection
  HIPTailDescriptor* tails = GetHIPTails(clusterer, row);
 
#ifdef GPUCA_DETERMINISTIC_MODE
  // Races in tail comparisons and atomic swap can lead to slightly different clusters.
  // So need a sequential fallback for deterministic mode
  GPUCommonAlgorithm::sortInBlock(tails + 1, tails + nTails + 1, [](auto&& t1, auto&& t2) {
    if (t1.pad != t2.pad) {
      return t1.pad < t2.pad;
    } else if (t1.tailStart != t2.tailStart) {
      return t1.tailStart < t2.tailStart;
    } else if (t1.tailEnd != t2.tailEnd) {
      return t1.tailEnd < t2.tailEnd;
    } else if (t1.satStart != t2.satStart) {
      return t1.satStart < t2.satStart;
    } else if (t1.satEnd != t2.satEnd) {
      return t1.satEnd < t2.satEnd;
    } else if (t1.qTot != t2.qTot) {
      return t1.qTot < t2.qTot;
    } else {
      return t1.qMax < t2.qMax;
    }
  } GPUCA_THREAD_INFO_PROVIDE);
  if (iThread > 0) {
    return;
  }
  nThreads = 1;
#endif
 
  for (uint32_t iTail = iThread + 1; iTail <= nTails; iTail += nThreads) {
    auto* tail = &tails[iTail];
 
    // TODO: this is needed because tailStarts may vary due to rising edge
    // Better approach would be to also track the triggered timebin and match that instead
    uint16_t overlapWindowStart = tail->tailStart >= 5 ? tail->tailStart - 5 : 0;
    uint16_t overlapWindowEnd = tail->tailStart + 5;
 
    for (uint32_t jTail = iTail + 1; jTail <= nTails; jTail++) {
      auto* tailNext = &tails[jTail];
      if (tailNext->iPrev > 0) {
        continue;
      }
 
      const bool overlapPad = tailNext->pad >= tail->pad - GPUTPCCFCheckPadBaseline::SSClusterPadWidth && tailNext->pad <= tail->pad + GPUTPCCFCheckPadBaseline::SSClusterPadWidth;
      const bool overlapTime = tailNext->tailStart >= overlapWindowStart && tailNext->tailStart < overlapWindowEnd;
 
      if (overlapPad && overlapTime) {
        if (CAMath::AtomicCAS(&tailNext->iPrev, 0u, iTail)) {
          tail->iNext = jTail;
          break;
        }
      }
    }
  }
}
 
// ======== HIP Clusterizer Kernel ========
 
template <>
GPUd() void GPUTPCCFHIPClusterizer::Thread<0>(int32_t nBlocks, int32_t nThreads, int32_t iBlock, int32_t iThread, GPUSharedMemory& smem, processorType& clusterer, uint8_t onlyMC)
{
  if (iBlock >= (int32_t)GPUTPCGeometry::NROWS) {
    return;
  }
 
  const uint32_t row = iBlock;
  uint32_t nTails = clusterer.mPnHIPTails[row];
  nTails = CAMath::Min(nTails, (uint32_t)MaxHIPTailsPerRow - 1);
 
  const auto* tails = GetHIPTails(clusterer, row);
  const auto& frag = clusterer.mPmemory->frag;
 
  auto* clusterPosInRow = clusterer.mPhipClusterPosInRow
                            ? clusterer.mPhipClusterPosInRow + row * MaxHIPTailsPerRow
                            : nullptr;
 
  for (uint32_t iTail = iThread + 1; iTail <= nTails; iTail += nThreads) {
 
    const auto* tail = &tails[iTail];
    if (tail->iPrev != 0) {
      continue;
    }
 
    CPU_ONLY(auto labelAcc = MCLabelAccumulator{clusterer});
 
    float qTot = 0;
    float qMax = 0;
    float padSum = 0;
    float padSqSum = 0;
    float satTimeSum = 0;
    uint32_t nSatTails = 0;
    uint32_t tailStart = (uint32_t)-1;
    uint32_t tailEnd = 0;
 
    // Zero-th element is empty tail
    for (; tail != tails; tail = &tails[tail->iNext]) {
      const float tailWeight = tail->qTot;
      const float tailPad = tail->pad;
      qMax = CAMath::Max(qMax, tail->qMax);
      qTot += tail->qTot;
      padSum += tailWeight * tailPad;
      padSqSum += tailWeight * tailPad * tailPad;
      if (tail->satStart >= 0 && tail->satEnd >= 0) {
        satTimeSum += 0.5f * (tail->satStart + tail->satEnd);
        nSatTails++;
      }
      tailStart = CAMath::Min<uint32_t>(tailStart, tail->tailStart);
      tailEnd = CAMath::Max<uint32_t>(tailEnd, tail->tailEnd);
 
      CPU_ONLY(labelAcc.collectTail(row, tail->pad, tail->tailStart, tail->tailEnd));
    }
 
    const float weightSum = CAMath::Max(qTot, 1.f);
    const float padMean = padSum / weightSum;
    const float padSigma = CAMath::Sqrt(CAMath::Max(0.f, padSqSum / weightSum - padMean * padMean));
 
    tpc::ClusterNative cn;
    cn.qMax = qMax;
    cn.setSaturatedQtot(qTot);
    cn.setSaturatedTailLength(tailEnd - tailStart);
    cn.setPad(padMean);
    cn.setSigmaPad(padSigma);
 
    if (cn.qMax >= 1023) {
 
      // Use the middle of the saturated plateau, averaged over all tails of the cluster that were triggered by saturation on their own pad.
      // Computed only here: chains consisting only of tails inherited from neighboring pads have no saturated plateau,
      // but these never contain a saturated charge and are dropped by the qMax cut.
      const float clusterTime = frag.start + satTimeSum / nSatTails - clusterer.Param().rec.tpc.clustersShiftTimebinsClusterizer;
      assert(!CAMath::IsNaN(clusterTime));
      cn.setTimeFlags(clusterTime, 0);
 
      uint32_t index;
 
      if (!onlyMC) {
        // Cut off clusters where the tail connection failed for some reason
        // TODO: Deduplicate with GPUTPCCFClusterizer::sortIntoBuckets (can't call cross-kernel).
        // TODO: Add error reporting for row cluster overflow.
        index = CAMath::AtomicAdd(&clusterer.mPclusterInRow[row], 1u);
        if (index < clusterer.mNMaxClusterPerRow) {
          clusterer.mPclusterByRow[clusterer.mNMaxClusterPerRow * row + index] = cn;
        }
        if (clusterPosInRow) {
          clusterPosInRow[iTail] = index;
        }
      } else {
        index = clusterPosInRow[iTail];
      }
 
      CPU_ONLY(labelAcc.commit(row, index, clusterer.mNMaxClusterPerRow));
    }
 
  } // for (uint32_t iTail = iThread + 1; iTail <= nTails; iTail += nThreads)
}