GPU TPC: If value from other side of interpolation not stored, use only current value

Author: davidrohr

GPU/GPUTracking/Definitions/GPUSettingsList.h

@@ -109,6 +109,8 @@ AddOptionRTC(rejectEdgeClustersMargin, float, 0.f, "", 0, "Margin in cm of Y pos
 AddOptionRTC(rejectEdgeClustersSigmaMargin, float, 0.f, "", 0, "Margin factor for trackSigmaY when rejecting edge clusters based on uncorrected track Y")
 AddOptionRTC(trackletMaxSharedFraction, float, 0.1f, "", 0, "Max fraction of shared clusters for tracklet")
 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(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")
@@ -140,7 +142,6 @@ AddOptionRTC(cfEdgeTwoPads, uint8_t, 0, "", 0, "Flag clusters with peak on the 2
 AddOptionRTC(nWays, uint8_t, 3, "", 0, "Do N fit passes in final fit of merger (must be odd to end with inward fit)")
 AddOptionRTC(rebuildTrackInFit, uint8_t, 1, "", 0, "Rebuild track completely during fit based on clusters closed to interpolated track positions")
 AddOptionRTC(rebuildTrackInFitClusterCandidates, uint8_t, 3, "", 0, "Number of cluster candidates per row for rebuilt track")
-AddOptionRTC(rebuildTrackMaxSharedFraction, float, 0.1f, "", 0, "Max fraction of shared clusters for rebuilt tracks")
 AddOptionRTC(trackFitRejectMode, int8_t, 5, "", 0, "0: no limit on rejection or missed hits, >0: break after n rejected hits")
 AddOptionRTC(rejectIFCLowRadiusCluster, uint8_t, 1, "", 0, "Reject clusters that get the IFC mask error during refit")
 AddOptionRTC(dEdxTruncLow, uint8_t, 2, "", 0, "Low truncation threshold, fraction of 128")

GPU/GPUTracking/Global/GPUChainTracking.cxx

@@ -1002,6 +1002,7 @@ void GPUChainTracking::ApplySyncSettings(GPUSettingsProcessing& proc, GPUSetting
 {
   if (syncMode) {
     rec.useMatLUT = false;
+    rec.tpc.rebuildTrackMaxNonIntCov = 0.f;
   }
   if (proc.rtc.optSpecialCode == -1) {
     proc.rtc.optSpecialCode = syncMode;

GPU/GPUTracking/Merger/GPUTPCGMTrackParam.cxx

@@ -406,27 +406,32 @@ GPUdii() float GPUTPCGMTrackParam::FindBestInterpolatedHit(GPUTPCGMMerger& GPUre
   GPUglobalref() const cahit2* hits = tracker.HitData(rowData);
   GPUglobalref() const calink* firsthit = tracker.FirstHitInBin(rowData);
   float uncorrectedY = -1e6f, uncorrectedZ;
-  if (rowData.NHits() && inter.errorY >= (GPUCA_PAR_MERGER_INTERPOLATION_ERROR_TYPE_A)0) {
+  if (rowData.NHits() && (inter.errorY >= (GPUCA_PAR_MERGER_INTERPOLATION_ERROR_TYPE_A)0 || (param.rec.tpc.rebuildTrackMaxNonIntCov > 0 && mC[0] < param.rec.tpc.rebuildTrackMaxNonIntCov && mC[2] < param.rec.tpc.rebuildTrackMaxNonIntCov))) {
     const float zOffset = param.par.continuousTracking ? merger.GetConstantMem()->calibObjects.fastTransform->convVertexTimeToZOffset(sector, mTOffset, param.continuousMaxTimeBin) : 0;
     const float y0 = rowData.Grid().YMin();
     const float stepY = rowData.HstepY();
     const float z0 = rowData.Grid().ZMin() - zOffset; // We can use our own ZOffset, since this is only used temporarily anyway
     const float stepZ = rowData.HstepZ();
     int32_t bin, ny, nz;
 
-    float err2Y, err2Z;
-    param.GetClusterErrors2(sector, row, mP[1], mP[2], mP[3], -1.f, 0.f, 0.f, err2Y, err2Z); // TODO: Use correct time/avgCharge
-
-    const float Iz0 = inter.posY - mP[0];
-    const float Iz1 = inter.posZ + deltaZ - mP[1];
-    const float Iw0 = 1.f / (mC[0] + (float)inter.errorY);
-    const float Iw2 = 1.f / (mC[2] + (float)inter.errorZ);
-    const float Ik00 = mC[0] * Iw0;
-    const float Ik11 = mC[2] * Iw2;
-    const float ImP0 = mP[0] + Ik00 * Iz0;
-    const float ImP1 = mP[1] + Ik11 * Iz1;
-    const float ImC0 = mC[0] - Ik00 * mC[0];
-    const float ImC2 = mC[2] - Ik11 * mC[2];
+    float ImP0, ImP1, ImC0, ImC2;
+    if (inter.errorY >= (GPUCA_PAR_MERGER_INTERPOLATION_ERROR_TYPE_A)0) {
+      const float Iz0 = inter.posY - mP[0];
+      const float Iz1 = inter.posZ + deltaZ - mP[1];
+      const float Iw0 = 1.f / (mC[0] + (float)inter.errorY);
+      const float Iw2 = 1.f / (mC[2] + (float)inter.errorZ);
+      const float Ik00 = mC[0] * Iw0;
+      const float Ik11 = mC[2] * Iw2;
+      ImP0 = mP[0] + Ik00 * Iz0;
+      ImP1 = mP[1] + Ik11 * Iz1;
+      ImC0 = mC[0] - Ik00 * mC[0];
+      ImC2 = mC[2] - Ik11 * mC[2];
+    } else {
+      ImP0 = mP[0];
+      ImP1 = mP[1];
+      ImC0 = mC[0];
+      ImC2 = mC[2];
+    }
 
     merger.GetConstantMem()->calibObjects.fastTransform->InverseTransformYZtoNominalYZ(sector, row, ImP0, ImP1, uncorrectedY, uncorrectedZ);
 
@@ -435,6 +440,9 @@ GPUdii() float GPUTPCGMTrackParam::FindBestInterpolatedHit(GPUTPCGMMerger& GPUre
       nCandidates++;
     }
     if (CAMath::Abs(uncorrectedY) <= rowData.getTPCMaxY()) {
+      float err2Y, err2Z;
+      param.GetClusterErrors2(sector, row, mP[1], mP[2], mP[3], -1.f, 0.f, 0.f, err2Y, err2Z); // TODO: Use correct time/avgCharge
+
       const float kFactor = tracker.GetChiSeedFactor();
       const float sy2 = 4 * CAMath::Min(param.rec.tpc.hitSearchArea2, kFactor * (err2Y + CAMath::Abs(mC[0]))); // TODO: is 4 a good factor??
       const float sz2 = 4 * CAMath::Min(param.rec.tpc.hitSearchArea2, kFactor * (err2Z + CAMath::Abs(mC[2])));