GPU TPC: Improve track extrapolation for better efficiency / less fake rate

Author: davidrohr

GPU/GPUTracking/Definitions/GPUSettingsList.h

@@ -104,6 +104,13 @@ AddOptionRTC(trackletMaxSharedFraction, float, 0.1f, "", 0, "Max fraction of sha
 AddOptionRTC(trackletMinSharedNormFactor, float, 0.f, "", 0, "Max shared defined as trackletMinSharedNormFactor*max(current_nhits,trackletMinSharedNormFactor*minHits,1)")
 AddOptionRTC(rebuildTrackMaxSharedFraction, float, 0.1f, "", 0, "Max fraction of shared clusters for rebuilt tracks")
 AddOptionRTC(rebuildTrackMaxNonIntCov, float, 4.f, "", 0, "Max Err2 allowed for non-interpolated cluster attachment during rebuild")
+AddOptionRTC(rebuildErrorWeightFactor, float, 5.f, "", 0, "Factor for the error (distance^2) of a cluster candidate in its hit weight 1 / (1 + factor * error) during rebuild")
+AddOptionRTC(rebuildTrackExtrMaxAvgChi2, float, 4.f, "", 0, "Max average chi2 per cluster to accept an extrapolated segment (0 = disabled)")
+AddOptionRTC(rebuildTrackAltCandidateMaxErrFactor, float, 2.f, "", 0, "Alternative cluster candidates of a row are used only if their distance^2 is below this factor times the distance^2 of the best candidate")
+AddOptionRTC(rebuildTrackAltCandidateGoodErr, float, 0.25f, "", 0, "Alternative cluster candidates of a row with distance^2 below this value are always used")
+AddOptionRTC(rebuildTrackErrFactorOneSided, float, 1.5f, "", 0, "Factor applied to the error (distance^2) of cluster candidates found with the track position from only one side during rebuild")
+AddOptionRTC(rebuildTrackErrFactorExtrapolation, float, 1.f, "", 0, "Factor applied to the error (distance^2) of cluster candidates found by track extrapolation during rebuild")
+AddOptionRTC(rebuildTrackMaxInterpolationDistY, float, 4.f, "", 0, "Max distance in Y between the track positions from both sides to use their interpolation for the cluster search in rows without cluster during rebuild")
 AddOptionRTC(maxTimeBinAboveThresholdIn1000Bin, uint16_t, 500, "", 0, "Except pad from cluster finding if total number of charges in a fragment is above this baseline (disable = 0)")
 AddOptionRTC(maxConsecTimeBinAboveThreshold, uint16_t, 200, "", 0, "Except pad from cluster finding if number of consecutive charges in a fragment is above this baseline (disable = 0)")
 AddOptionRTC(noisyPadSaturationThreshold, uint16_t, 700, "", 0, "Threshold where a timebin is considered saturated, disabling the noisy pad check for that pad")
@@ -168,9 +175,12 @@ AddOptionRTC(rejectEdgeClustersInSeeding, int8_t, 0, "", 0, "Reject edge cluster
 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")
-AddOptionRTC(rebuildTrackExtrMaxMissingRows, uint8_t, 15, "", 0, "Maximum total number of rows allowed to be missing during track extrapolation")
-AddOptionRTC(rebuildTrackExtrMinConsecGoodRows, uint8_t, 6, "", 0, "Minimum number of consecutive rows required to accept extrapolated segment")
+AddOptionRTC(rebuildTrackExtrMaxMissingRows, uint8_t, 20, "", 0, "Maximum total number of rows allowed to be missing during track extrapolation")
+AddOptionRTC(rebuildTrackExtrMinConsecGoodRows, uint8_t, 8, "", 0, "Minimum number of consecutive rows required to accept extrapolated segment")
+AddOptionRTC(rebuildTrackExtrMinConsecGoodRowsEdge, uint8_t, 3, "", 0, "Minimum number of consecutive rows required to accept extrapolated segment")
 AddOptionRTC(rebuildTrackExtrConsecGoodRowsMaxGap, uint8_t, 1, "", 0, "Max row gap size to consider the segment as consecutive rows")
+AddOptionRTC(rebuildTrackCandidateErrChi2, uint8_t, 0, "", 0, "Rank cluster candidates during rebuild by chi2 instead of distance^2 (rebuildTrackAltCandidateGoodErr is then in units of chi2)")
+AddOptionRTC(rebuildTrackNonIntMinNDF, uint8_t, 10, "", 0, "NDF the track must exceed for non-interpolated cluster attachment (one-sided search and extrapolation) during rebuild")
 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 handle AddOptionArrayRTC
 AddHelp("help", 'h')
 EndConfig()

GPU/GPUTracking/Global/GPUChainTrackingMerger.cxx

@@ -263,10 +263,12 @@ int32_t GPUChainTracking::RunTPCTrackingMerger(bool synchronizeOutput)
   runKernel<GPUTPCGMMergerFinalize, 2>(GetGridAuto(0, deviceType));
   if (param().rec.tpc.mergeLoopersAfterburner) {
     runKernel<GPUTPCGMMergerMergeLoopers, 0>(doGPU ? GetGrid(Merger.NMergedTracks(), 0, deviceType) : GetGridAuto(0, deviceType));
-    if (doGPU) {
-      TransferMemoryResourceLinkToHost(RecoStep::TPCMerging, Merger.MemoryResMemory(), 0);
-      SynchronizeStream(0); // TODO: could probably synchronize on an event after runKernel<GPUTPCGMMergerMergeLoopers, 1>
-    }
+  }
+  if (doGPU) {
+    TransferMemoryResourceLinkToHost(RecoStep::TPCMerging, Merger.MemoryResMemory(), 0);
+    SynchronizeStream(0); // TODO: could probably synchronize on an event after runKernel<GPUTPCGMMergerMergeLoopers, 1>
+  }
+  if (param().rec.tpc.mergeLoopersAfterburner) {
     runKernel<GPUTPCGMMergerMergeLoopers, 1>(GetGridAuto(0, deviceType));
     runKernel<GPUTPCGMMergerMergeLoopers, 2>(doGPU ? GetGrid(Merger.Memory()->nLooperMatchCandidates, 0, deviceType) : GetGridAuto(0, deviceType));
     DoDebugAndDump(RecoStep::TPCMerging, GPUChainTrackingDebugFlags::TPCMergingLoopers, Merger, &GPUTPCGMMerger::DumpLoopers, *mDebugFile);

GPU/GPUTracking/Merger/GPUTPCGMMerger.cxx

@@ -2184,7 +2184,7 @@ GPUd() void GPUTPCGMMerger::PrepareHitWeights(int32_t nBlocks, int32_t nThreads,
         auto* candidates = &mClusterCandidates[(i * GPUTPCGeometry::NROWS + cl.row) * Param().rec.tpc.rebuildTrackInFitClusterCandidates];
         if (cl.row != lastRow && candidates[0].id == 0 && (!(cl.state & GPUTPCGMMergedTrackHit::flagReject) || trk.GetParam().GetNDF() <= 0)) {
           candidates[0].id = cl.num + 2;
-          candidates[0].best = 128; // TODO: Needs better handling in the future, particularly should be marked as shared correctly
+          candidates[0].best = 128; // TODO: Needs better handling in the future, particularly should be marked as shared correctly. Tracks with NDF <= 0 get no fallback in rows with candidates, and get those back only via the shared recovery
           candidates[0].weight = cl.state;
           candidates[0].sector = cl.sector;
           lastRow = cl.row;
@@ -2204,15 +2204,14 @@ GPUd() void GPUTPCGMMerger::ComputeHitWeights(int32_t nBlocks, int32_t nThreads,
     for (uint32_t j = 0; j < GPUTPCGeometry::NROWS; j++) {
       auto* candidates = &mClusterCandidates[(i * GPUTPCGeometry::NROWS + j) * Param().rec.tpc.rebuildTrackInFitClusterCandidates];
       if (candidates[0].best == 0) {
-        if (candidates[iteration].id >= 2) {
-          const int32_t nCl = trk.NClusters();
-          float weightf = nCl * (nCl * 2 - 5) * 128 / trk.GetParam().GetChi2() / (1.f + candidates[iteration].error);
+        if (candidates[iteration].id >= 2 && (iteration == 0 || candidates[iteration].error <= CAMath::Max(Param().rec.tpc.rebuildTrackAltCandidateGoodErr, candidates[0].error * Param().rec.tpc.rebuildTrackAltCandidateMaxErrFactor))) { // Don't fall back to a much worse candidate, the best one is then rather a shared cluster
+          float weightf = (trk.GetParam().GetNDF() + 5) * trk.GetParam().GetNDF() * 64 / trk.GetParam().GetChi2() / (1.f + Param().rec.tpc.rebuildErrorWeightFactor * candidates[iteration].error);
           int32_t weight = weightf < 0.f ? 0 : (weightf > 2e9f ? 2e9f : (int32_t)weightf);
           candidates[iteration].weight = weight;
           CAMath::AtomicMax(&mHitWeights[candidates[iteration].id - 2], weight);
           CADEBUG(printf("REBUILD: iTrk %d Iteration %d Row %d Cluster %d Weight %d\n", i, iteration, j, candidates[iteration].id - 2, weight));
         } else {
-          candidates[0].best = -1;
+          candidates[0].best = -2;
         }
       }
     }
@@ -2259,7 +2258,7 @@ GPUd() void GPUTPCGMMerger::ResolveHitWeights2(int32_t nBlocks, int32_t nThreads
         if (candidates[0].best != 128 && candidates[candidates[0].best - 1].weight != mHitWeights[candidates[candidates[0].best - 1].id - 2]) {
           mSharedCount2[candidates[candidates[0].best - 1].id - 2] = 1; // shared between tracks of different iterations
         }
-      } else if (candidates[0].best == 0 && candidates[0].id >= 2) {
+      } else if (candidates[0].best != -1 && candidates[0].id >= 2) {
         sharedRowCandidates[nSharedRowCandidates++] = j;
       }
     }