GPU TPC: Dynamically increase the protection / decrease the removal tube at inner pad rows, edge pads, high local occupancy

Author: davidrohr

GPU/GPUTracking/Base/GPUReconstruction.cxx

@@ -748,7 +748,7 @@ void* GPUReconstruction::AllocateDirectMemory(size_t size, int32_t type)
   void*& poolend = (type & GPUMemoryResource::MEMORY_GPU) ? mDeviceMemoryPoolEnd : mHostMemoryPoolEnd;
   char* retVal;
   if ((type & GPUMemoryResource::MEMORY_STACK)) {
-    poolend = (char*)poolend - size;
+    poolend = (char*)poolend - size; // TODO: Implement overflow check
     poolend = (char*)poolend - GPUProcessor::getAlignmentMod<GPUCA_MEMALIGN>(poolend);
     retVal = (char*)poolend;
   } else {

GPU/GPUTracking/Definitions/GPUSettingsList.h

@@ -68,11 +68,12 @@ AddOptionRTC(minNClustersFinalTrack, int32_t, -1, "", 0, "required min number of
 AddOptionRTC(searchWindowDZDR, float, 2.5f, "", 0, "Use DZDR window for seeding instead of neighboursSearchArea")
 AddOptionRTC(trackReferenceX, float, 1000.f, "", 0, "Transport all tracks to this X after tracking (disabled if > 500, auto = 1000)")
 AddOptionRTC(zsThreshold, float, 2.0f, "", 0, "Zero-Suppression threshold")
-AddOptionRTC(tubeProtectSigma2, float, 5.f * 5.f, "", 0, "Max sigma2 to mark adjacent cluster for protection")
-AddOptionRTC(tubeProtectMaxSize2, float, 3.5f * 3.5f, "", 0, "Square of max tube size (if smaller than tubeProtectChi2)")
-AddOptionRTC(tubeProtectMinSize2, float, 1.0f * 1.0f, "", 0, "Square of min tube size (if larger than tubeProtectChi2)")
-AddOptionRTC(tubeRemoveSigma2, float, 1.f * 1.f, "", 0, "Max sigma2 to mark adjacent cluster for removal")
-AddOptionRTC(tubeRemoveMaxSize2, float, 1.5f * 1.5f, "", 0, "Square of max tube size (if smaller than tubeRejectChi2)")
+AddOptionRTC(tubeProtectSigma2, float, 4.f * 4.f, "", 0, "Max sigma2 to mark adjacent cluster for protection")
+AddOptionRTC(tubeProtectMaxSize2, float, 2.f * 2.f, "", 0, "Square of max tube size (if smaller than tubeProtectChi2)")
+AddOptionRTC(tubeProtectMinSize2, float, 0.5f * 0.5f, "", 0, "Square of min tube size (if larger than tubeProtectChi2)")
+AddOptionRTC(tubeRemoveSigma2, float, 1.25f * 1.25f, "", 0, "Max sigma2 to mark adjacent cluster for removal")
+AddOptionRTC(tubeRemoveMaxSize2, float, 2.5f * 2.5f, "", 0, "Square of max tube size (if smaller than tubeRejectChi2)")
+AddOptionRTC(tubeExtraProtectMinOccupancy, uint32_t, 1500, "", 0, "Increase Protection, decrease removal by factor 2, when above this lokal occupancy / rowx")
 AddOptionRTC(clustersShiftTimebins, float, 0, "", 0, "Shift of TPC clusters (applied during CTF cluster decoding)")
 AddOptionRTC(clustersShiftTimebinsClusterizer, float, 0, "", 0, "Shift of TPC clusters (applied during CTF clusterization)")
 AddOptionRTC(clustersEdgeFixDistance, float, 0.f, "", 0, "If >0, revert cluster.flag edge bit distance to edge exceeds this parameter (fixed during CTF decoding)")
@@ -162,6 +163,9 @@ AddOptionRTC(dEdxClusterRejectionFlagMask, int8_t, o2::gpu::GPUTPCGMMergedTrackH
 AddOptionRTC(dEdxClusterRejectionFlagMaskAlt, int8_t, o2::gpu::GPUTPCGMMergedTrackHit::flagEdge, "", 0, "OR mask of TPC flags that will reject the cluster in alternative dEdx")
 AddOptionRTC(rejectEdgeClustersInSeeding, int8_t, 0, "", 0, "Reject edge clusters based on uncorrected track Y during seeding")
 AddOptionRTC(rejectEdgeClustersInTrackFit, int8_t, 0, "", 0, "Reject edge clusters based on uncorrected track Y during track fit")
+AddOptionRTC(tubeExtraProtectMinRow, uint8_t, 20, "", 0, "Increase Protection, decrease removal by factor 2, when below this row")
+AddOptionRTC(tubeExtraProtectEdgePads, uint8_t, 2, "", 0, "Increase Protection, decrease removal by factor 2, when on this number of pads from the edge")
+
 AddOptionArray(PID_remap, int8_t, 9, (0, 1, 2, 3, 4, 5, 6, 7, 8), "", 0, "Remap Ipid to PID_reamp[Ipid] (no remap if<0)") // BUG: CUDA cannot yet hand AddOptionArrayRTC
 AddHelp("help", 'h')
 EndConfig()

GPU/GPUTracking/Merger/GPUTPCGMTrackParam.cxx

@@ -248,7 +248,6 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger* GPUrestrict() merger, int32_
         const float invCharge = merger->GetConstantMem()->ioPtrs.clustersNative ? (1.f / merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear[cluster.num].qMax) : 0.f;
         float invAvgCharge = (sumInvSqrtCharge += invSqrtCharge) / ++nAvgCharge;
         invAvgCharge *= invAvgCharge;
-
         prop.GetErr2(err2Y, err2Z, param, zz, cluster.row, clusterState, cluster.sector, time, invAvgCharge, invCharge);
 
         if (rejectChi2 >= GPUTPCGMPropagator::rejectInterFill) {
@@ -491,7 +490,8 @@ GPUd() float GPUTPCGMTrackParam::AttachClusters(const GPUTPCGMMerger* GPUrestric
 
 GPUd() float GPUTPCGMTrackParam::AttachClusters(const GPUTPCGMMerger* GPUrestrict() Merger, int32_t sector, int32_t iRow, int32_t iTrack, bool goodLeg, float Y, float Z)
 {
-  if (Merger->Param().rec.tpc.disableRefitAttachment & 1) {
+  const auto& param = Merger->Param();
+  if (param.rec.tpc.disableRefitAttachment & 1) {
     return -1e6f;
   }
   const GPUTPCTracker& GPUrestrict() tracker = *(Merger->GetConstantMem()->tpcTrackers + sector);
@@ -502,34 +502,40 @@ GPUd() float GPUTPCGMTrackParam::AttachClusters(const GPUTPCGMMerger* GPUrestric
     return -1e6f;
   }
 
-  const float zOffset = Merger->Param().par.continuousTracking ? Merger->GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->convVertexTimeToZOffset(sector, mTOffset, Merger->Param().continuousMaxTimeBin) : 0;
+  const float zOffset = param.par.continuousTracking ? Merger->GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->convVertexTimeToZOffset(sector, mTOffset, param.continuousMaxTimeBin) : 0; // TODO: do some validatiomns for the transform conv functions...
   const float y0 = row.Grid().YMin();
   const float stepY = row.HstepY();
   const float z0 = row.Grid().ZMin() - zOffset; // We can use our own ZOffset, since this is only used temporarily anyway
   const float stepZ = row.HstepZ();
   int32_t bin, ny, nz;
 
-  bool protect = CAMath::Abs(GetQPt() * Merger->Param().qptB5Scaler) <= Merger->Param().rec.tpc.rejectQPtB5 && goodLeg;
+  float uncorrectedY, uncorrectedZ;
+  Merger->GetConstantMem()->calibObjects.fastTransformHelper->InverseTransformYZtoNominalYZ(sector, iRow, Y, Z, uncorrectedY, uncorrectedZ);
+  if (CAMath::Abs(uncorrectedY) > row.getTPCMaxY()) {
+    return uncorrectedY;
+  }
 
+  bool protect = CAMath::Abs(GetQPt() * param.qptB5Scaler) <= param.rec.tpc.rejectQPtB5 && goodLeg;
   float err2Y, err2Z;
-  Merger->Param().GetClusterErrors2(sector, iRow, Z, mP[2], mP[3], -1.f, 0.f, 0.f, err2Y, err2Z);                                       // TODO: Use correct time/avgCharge
-  const float tubeMaxSize2 = protect ? Merger->Param().rec.tpc.tubeProtectMaxSize2 : Merger->Param().rec.tpc.tubeRemoveMaxSize2;
-  const float tubeMinSize2 = protect ? Merger->Param().rec.tpc.tubeProtectMinSize2 : 0.f;
-  const float tubeSigma2 = protect ? Merger->Param().rec.tpc.tubeProtectSigma2 : Merger->Param().rec.tpc.tubeRemoveSigma2;
+  param.GetClusterErrors2(sector, iRow, Z, mP[2], mP[3], -1.f, 0.f, 0.f, err2Y, err2Z); // TODO: Use correct time/avgCharge
+  const float tubeMaxSize2 = protect ? param.rec.tpc.tubeProtectMaxSize2 : param.rec.tpc.tubeRemoveMaxSize2;
+  const float tubeMinSize2 = protect ? param.rec.tpc.tubeProtectMinSize2 : 0.f;
+  float tubeSigma2 = protect ? param.rec.tpc.tubeProtectSigma2 : param.rec.tpc.tubeRemoveSigma2;
+  uint32_t pad = CAMath::Float2UIntRn(GPUTPCGeometry::LinearY2Pad(sector, iRow, uncorrectedY));
+  float time = Merger->GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->InverseTransformInTimeFrame(sector, mP[1] + (param.par.continuousTracking ? Merger->GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->convVertexTimeToZOffset(sector, mTOffset, param.continuousMaxTimeBin) : 0), param.continuousMaxTimeBin); // TODO: Simplify this call in TPCFastTransform
+  if (iRow < param.rec.tpc.tubeExtraProtectMinRow ||
+      pad < param.rec.tpc.tubeExtraProtectEdgePads || pad >= (uint32_t)(GPUTPCGeometry::NPads(iRow) - param.rec.tpc.tubeExtraProtectEdgePads) ||
+      param.GetUnscaledMult(time) / GPUTPCGeometry::Row2X(iRow) > param.rec.tpc.tubeExtraProtectMinOccupancy) {
+    tubeSigma2 *= protect ? 2 : 0.5;
+  }
   const float sy2 = CAMath::Max(tubeMinSize2, CAMath::Min(tubeMaxSize2, tubeSigma2 * (err2Y + CAMath::Abs(mC[0])))); // Cov can be bogus when following circle
   const float sz2 = CAMath::Max(tubeMinSize2, CAMath::Min(tubeMaxSize2, tubeSigma2 * (err2Z + CAMath::Abs(mC[2])))); // In that case we should provide the track error externally
   const float tubeY = CAMath::Sqrt(sy2);
   const float tubeZ = CAMath::Sqrt(sz2);
   const float sy21 = 1.f / sy2;
   const float sz21 = 1.f / sz2;
-  float uncorrectedY, uncorrectedZ;
-  Merger->GetConstantMem()->calibObjects.fastTransformHelper->InverseTransformYZtoNominalYZ(sector, iRow, Y, Z, uncorrectedY, uncorrectedZ);
 
-  if (CAMath::Abs(uncorrectedY) > row.getTPCMaxY()) {
-    return uncorrectedY;
-  }
   row.Grid().GetBinArea(uncorrectedY, uncorrectedZ + zOffset, tubeY, tubeZ, bin, ny, nz);
-
   const int32_t nBinsY = row.Grid().Ny();
   const int32_t idOffset = tracker.Data().ClusterIdOffset();
   const int32_t* ids = &(tracker.Data().ClusterDataIndex()[row.HitNumberOffset()]);

GPU/GPUTracking/TPCConvert/GPUTPCConvertImpl.h

@@ -28,7 +28,7 @@ class GPUTPCConvertImpl
  public:
   GPUd() static void convert(const GPUConstantMem& GPUrestrict() cm, int32_t sector, int32_t row, float pad, float time, float& GPUrestrict() x, float& GPUrestrict() y, float& GPUrestrict() z)
   {
-    if (cm.param.par.continuousTracking) {
+    if (cm.param.par.continuousTracking) { // TODO: This might be wrong, don't we just need to do TransformInTimeframe always
       cm.calibObjects.fastTransformHelper->getCorrMap()->TransformInTimeFrame(sector, row, pad, time, x, y, z, cm.param.continuousMaxTimeBin);
     } else {
       cm.calibObjects.fastTransformHelper->Transform(sector, row, pad, time, x, y, z);

GPU/TPCFastTransformation/TPCFastTransform.h

@@ -194,6 +194,7 @@ class TPCFastTransform : public FlatObject
 
   /// Inverse transformation
   GPUd() void InverseTransformInTimeFrame(int32_t slice, int32_t row, float /*x*/, float y, float z, float& pad, float& time, float maxTimeBin) const;
+  GPUd() float InverseTransformInTimeFrame(int32_t slice, float z, float maxTimeBin) const;
 
   /// Inverse transformation: Transformed Y and Z -> transformed X
   GPUd() void InverseTransformYZtoX(int32_t slice, int32_t row, float y, float z, float& x, const TPCFastTransform* ref = nullptr, const TPCFastTransform* ref2 = nullptr, float scale = 0.f, float scale2 = 0.f, int32_t scaleMode = 0) const;
@@ -667,6 +668,13 @@ GPUdi() void TPCFastTransform::InverseTransformInTimeFrame(int32_t slice, int32_
   convUVtoPadTimeInTimeFrame(slice, row, u, v, pad, time, maxTimeBin);
 }
 
+GPUdi() float TPCFastTransform::InverseTransformInTimeFrame(int32_t slice, float z, float maxTimeBin) const
+{
+  float pad, time;
+  InverseTransformInTimeFrame(slice, 0, 0, 0, z, pad, time, maxTimeBin);
+  return time;
+}
+
 GPUdi() void TPCFastTransform::TransformIdealZ(int32_t slice, float time, float& z, float vertexTime) const
 {
   /// _______________ The main method: cluster transformation _______________________