GPU TPC: Fix deltaTimeFwd/deltaTimeBwd after track rebuild

Author: davidrohr

GPU/GPUTracking/Merger/GPUTPCGMO2Output.cxx

@@ -215,9 +215,11 @@ GPUdii() void GPUTPCGMO2Output::Thread<GPUTPCGMO2Output::output>(int32_t nBlocks
 
     if (track.PrevSegment() >= 0) {
       const GPUTPCGMMergedTrack* chkTrk = track.GetFirstSegment(tracks, merger.Param().rec.enableCyclicGraphWorkarounds);
-      const auto& firstPrevCluster = trackClusters[chkTrk->FirstClusterRef()];
-      t1 = clusters->clustersLinear[firstPrevCluster.num].getTime();
-      sector1 = firstPrevCluster.sector;
+      if (chkTrk->NClusters()) { // First segment can have lost all its clusters in the rebuild, then its cluster reference is invalid
+        const auto& firstPrevCluster = trackClusters[chkTrk->FirstClusterRef()];
+        t1 = clusters->clustersLinear[firstPrevCluster.num].getTime();
+        sector1 = firstPrevCluster.sector;
+      }
     }
 
     bool cce = track.CCE() && ((sector1 < MAXSECTOR / 2) ^ (sector2 < MAXSECTOR / 2));

GPU/GPUTracking/Merger/GPUTPCGMTrackParam.cxx

@@ -81,6 +81,10 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger& GPUrestrict() merger, int32_
 
   const int32_t nWays = param.rec.tpc.nWays;
   const int32_t maxN = N;
+  if (rebuilt && maxN > 0 && param.par.continuousTracking && !track.CCE() && !track.Looper()) { // The rebuilt track can have other first / last clusters than the track the Z shift was computed for
+    const float t0 = merger.GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[0].num].getTime(), tN = merger.GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[maxN - 1].num].getTime();
+    ConstrainZ(merger, clusters[0].sector, CAMath::Max(t0, tN), CAMath::Min(t0, tN));
+  }
   int32_t ihitStart = 0;
   int32_t interpolatedStart = 0;
   float covYYUpd = 0.f;
@@ -1327,28 +1331,29 @@ GPUd() float GPUTPCGMTrackParam::ShiftZ(const GPUTPCGMMerger& GPUrestrict() merg
     merger.GetConstantMem()->calibObjects.fastTransform->TransformIdealZ(sector, cltmax, basez, mTOffset);
     deltaZ = basez - refZ;
   }
-  {
-    float deltaT = merger.GetConstantMem()->calibObjects.fastTransform->convDeltaZtoDeltaTimeInTimeFrame(sector, deltaZ);
+  mTOffset += merger.GetConstantMem()->calibObjects.fastTransform->convDeltaZtoDeltaTimeInTimeFrame(sector, deltaZ);
+  return -ConstrainZ(merger, sector, cltmax, cltmin, deltaZ);
+}
+
+GPUd() float GPUTPCGMTrackParam::ConstrainZ(const GPUTPCGMMerger& GPUrestrict() merger, uint32_t sector, float cltmax, float cltmin, float deltaZ)
+{
+  const float maxT = cltmin - merger.GetConstantMem()->calibObjects.fastTransform->getT0();
+  const float minT = cltmax - merger.GetConstantMem()->calibObjects.fastTransform->getMaxDriftTime(sector);
+  // printf("T Check: Clusters %f %f, min %f max %f vtx %f\n", tz1, tz2, minT, maxT, mTOffset);
+  float deltaT = 0.f;
+  if (mTOffset < minT) {
+    deltaT = minT - mTOffset;
+  }
+  if (mTOffset + deltaT > maxT) {
+    deltaT = maxT - mTOffset;
+  }
+  if (deltaT != 0.f) {
+    deltaZ += merger.GetConstantMem()->calibObjects.fastTransform->convDeltaTimeToDeltaZinTimeFrame(sector, deltaT);
+    // printf("Moving clusters to TPC Range: QPt %f, New mTOffset %f, t1 %f, t2 %f, Shift %f in Z: %f to %f --> %f to %f in T\n", mP[4], mTOffset + deltaT, tz1, tz2, deltaZ, tz2 - mTOffset, tz1 - mTOffset, tz2 - mTOffset - deltaT, tz1 - mTOffset - deltaT);
     mTOffset += deltaT;
-    const float maxT = cltmin - merger.GetConstantMem()->calibObjects.fastTransform->getT0();
-    const float minT = cltmax - merger.GetConstantMem()->calibObjects.fastTransform->getMaxDriftTime(sector);
-    // printf("T Check: Clusters %f %f, min %f max %f vtx %f\n", tz1, tz2, minT, maxT, mTOffset);
-    deltaT = 0.f;
-    if (mTOffset < minT) {
-      deltaT = minT - mTOffset;
-    }
-    if (mTOffset + deltaT > maxT) {
-      deltaT = maxT - mTOffset;
-    }
-    if (deltaT != 0.f) {
-      deltaZ += merger.GetConstantMem()->calibObjects.fastTransform->convDeltaTimeToDeltaZinTimeFrame(sector, deltaT);
-      // printf("Moving clusters to TPC Range: QPt %f, New mTOffset %f, t1 %f, t2 %f, Shift %f in Z: %f to %f --> %f to %f in T\n", mP[4], mTOffset + deltaT, tz1, tz2, deltaZ, tz2 - mTOffset, tz1 - mTOffset, tz2 - mTOffset - deltaT, tz1 - mTOffset - deltaT);
-      mTOffset += deltaT;
-    }
-    mP[1] -= deltaZ;
   }
-  // printf("\n");
-  return -deltaZ;
+  mP[1] -= deltaZ;
+  return deltaZ;
 }
 
 GPUd() bool GPUTPCGMTrackParam::CheckCov() const

GPU/GPUTracking/Merger/GPUTPCGMTrackParam.h

@@ -203,6 +203,7 @@ class GPUTPCGMTrackParam
   GPUd() void Rotate(float alpha);
   GPUd() float ShiftZ(const GPUTPCGMMerger& merger, uint32_t sector, float cltmax, float cltmin, float clx);
   GPUd() float ShiftZ(const GPUTPCGMMergedTrackHit* clusters, const GPUTPCGMMerger& merger, int32_t N);
+  GPUd() float ConstrainZ(const GPUTPCGMMerger& merger, uint32_t sector, float cltmax, float cltmin, float deltaZ = 0.f); // Shift Z / TOffset such that the clusters are within the TPC drift volume, applies in addition the pending shift deltaZ to Z, returns the total shift
 
   GPUd() static float Reciprocal(float x) { return 1.f / x; }
   GPUdi() static void Assign(float& x, bool mask, float v)