GPU: simplifaction, refactoring and some minor improvements

Author: davidrohr

GPU/GPUTracking/Definitions/GPUSettingsList.h

@@ -137,7 +137,7 @@ AddOptionRTC(nWays, uint8_t, 3, "", 0, "Do N fit passes in final fit of merger (
 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, <0: reject at max -n hits")
+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")
 AddOptionRTC(dEdxTruncHigh, uint8_t, 77, "", 0, "High truncation threshold, fraction of 128")

GPU/GPUTracking/Merger/GPUTPCGMPropagator.cxx

@@ -608,15 +608,15 @@ GPUd() float GPUTPCGMPropagator::PredictChi2(float posY, float posZ, float err2Y
   }
 }
 
-GPUd() int32_t GPUTPCGMPropagator::Update(float posY, float posZ, int32_t iRow, const GPUParam& GPUrestrict() param, int16_t clusterState, int8_t rejectChi2, bool refit, int8_t sector, float time, float avgInvCharge, float invCharge)
+GPUd() int32_t GPUTPCGMPropagator::Update(float posY, float posZ, int32_t iRow, const GPUParam& GPUrestrict() param, int16_t clusterState, bool rejectChi2, bool refit, int8_t sector, float time, float avgInvCharge, float invCharge)
 {
   float err2Y, err2Z;
   GetErr2(err2Y, err2Z, param, posZ, iRow, clusterState, sector, time, avgInvCharge, invCharge);
 
   return Update(posY, posZ, iRow, param, clusterState, rejectChi2, refit, err2Y, err2Z);
 }
 
-GPUd() int32_t GPUTPCGMPropagator::Update(float posY, float posZ, int32_t iRow, const GPUParam& GPUrestrict() param, int16_t clusterState, int8_t rejectChi2, bool refit, float err2Y, float err2Z)
+GPUd() int32_t GPUTPCGMPropagator::Update(float posY, float posZ, int32_t iRow, const GPUParam& GPUrestrict() param, int16_t clusterState, bool rejectChi2, bool refit, float err2Y, float err2Z)
 {
   if (mT->NDF() == -5) { // first measurement: no need to filter, as the result is known in advance. just set it.
     mT->ResetCovariance();
@@ -635,80 +635,83 @@ GPUd() int32_t GPUTPCGMPropagator::Update(float posY, float posZ, int32_t iRow,
     return 0;
   }
 
-  return Update(posY, posZ, clusterState, rejectChi2 == rejectDirect, err2Y, err2Z, &param);
+  return Update(posY, posZ, clusterState, rejectChi2, err2Y, err2Z, &param);
 }
 
-GPUd() int32_t GPUTPCGMPropagator::InterpolateReject(const GPUParam& GPUrestrict() param, float posY, float posZ, int16_t clusterState, int8_t rejectChi2, gputpcgmmergertypes::InterpolationErrorHit* inter, float err2Y, float err2Z, float deltaZ)
+GPUd() void GPUTPCGMPropagator::InterpolateFill(float posY, float posZ, gputpcgmmergertypes::InterpolationErrorHit* inter)
 {
   float* GPUrestrict() mC = mT->Cov();
   float* GPUrestrict() mP = mT->Par();
-  if (rejectChi2 == rejectInterFill) {
-    inter->posY = mP[0];
-    inter->posZ = mP[1];
-    if (mT->NDF() <= 0) {
-      inter->errorY = inter->errorZ = 100.f;
-    } else {
-      inter->errorY = mC[0];
-      inter->errorZ = mC[2];
-    }
-  } else if (rejectChi2 == rejectInterReject) {
-    float chi2Y, chi2Z;
-    if (mT->NDF() <= 0) {
-      chi2Y = CAMath::Square((float)inter->posY - posY) / ((float)inter->errorY + err2Y);
-      chi2Z = CAMath::Square((float)inter->posZ + deltaZ - posZ) / ((float)inter->errorZ + err2Z);
-    } else if (mFitInProjections) {
-      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];
-      // printf("\t%21sInterpo ----- abde artaf%16s Y %8.3f, Z %8.3f (Errors %f <-- (%f, %f) %f <-- (%f, %f))\n", "", "", ImP0, ImP1, sqrtf(ImC0), sqrtf(mC[0]), sqrtf(inter->errorY), sqrtf(ImC2), sqrtf(mC[2]), sqrtf(inter->errorZ));
-      const float Jz0 = posY - ImP0;
-      const float Jz1 = posZ - ImP1;
-      const float Jw0 = 1.f / (ImC0 + err2Y);
-      const float Jw2 = 1.f / (ImC2 + err2Z);
-      chi2Y = Jw0 * Jz0 * Jz0;
-      chi2Z = Jw2 * Jz1 * Jz1;
-    } else {
-      const float Iz0 = inter->posY - mP[0];
-      const float Iz1 = inter->posZ + deltaZ - mP[1];
-      float Iw0 = mC[2] + (float)inter->errorZ;
-      float Iw2 = mC[0] + (float)inter->errorY;
-      float Idet = CAMath::Max(1e-10f, Iw0 * Iw2 - mC[1] * mC[1]);
-      Idet = 1.f / Idet;
-      Iw0 *= Idet;
-      const float Iw1 = mC[1] * Idet;
-      Iw2 *= Idet;
-      const float Ik00 = mC[0] * Iw0 + mC[1] * Iw1;
-      const float Ik01 = mC[0] * Iw1 + mC[1] * Iw2;
-      const float Ik10 = mC[1] * Iw0 + mC[2] * Iw1;
-      const float Ik11 = mC[1] * Iw1 + mC[2] * Iw2;
-      const float ImP0 = mP[0] + Ik00 * Iz0 + Ik01 * Iz1;
-      const float ImP1 = mP[1] + Ik10 * Iz0 + Ik11 * Iz1;
-      const float ImC0 = mC[0] - Ik00 * mC[0] + Ik01 * mC[1];
-      const float ImC1 = mC[1] - Ik10 * mC[0] + Ik11 * mC[1];
-      const float ImC2 = mC[2] - Ik10 * mC[1] + Ik11 * mC[2];
-      const float Jz0 = posY - ImP0;
-      const float Jz1 = posZ - ImP1;
-      float Jw0 = ImC2 + err2Z;
-      float Jw2 = ImC0 + err2Y;
-      float Jdet = CAMath::Max(1e-10f, Jw0 * Jw2 - ImC1 * ImC1);
-      Jdet = 1.f / Jdet;
-      Jw0 *= Jdet;
-      const float Jw1 = ImC1 * Jdet;
-      Jw2 *= Jdet;
-      chi2Y = CAMath::Abs((Jw0 * Jz0 + Jw1 * Jz1) * Jz0);
-      chi2Z = CAMath::Abs((Jw1 * Jz0 + Jw2 * Jz1) * Jz1);
-    }
-    if (RejectCluster(chi2Y * param.rec.tpc.clusterRejectChi2TolleranceY, chi2Z * param.rec.tpc.clusterRejectChi2TolleranceZ, clusterState)) { // TODO: Relative Pt resolution decreases slightly, why?
-      // printf("Reject Cluster chiy2 %f chiz2 %f (Pos Y: %f - %f %f ; Pos Z: %f - %f %f)\n", chi2Y, chi2Z, posY, mP[0], (float)inter->posY, posZ, mP[1], (float)inter->posZ + deltaZ);
-      return updateErrorClusterRejectedInInterpolation;
-    }
+  inter->posY = mP[0];
+  inter->posZ = mP[1];
+  if (mT->NDF() <= 0) {
+    inter->errorY = inter->errorZ = 100.f;
+  } else {
+    inter->errorY = mC[0];
+    inter->errorZ = mC[2];
+  }
+}
+
+GPUd() int32_t GPUTPCGMPropagator::InterpolateReject(const GPUParam& GPUrestrict() param, float posY, float posZ, int16_t clusterState, gputpcgmmergertypes::InterpolationErrorHit* inter, float err2Y, float err2Z, float deltaZ)
+{
+  float* GPUrestrict() mC = mT->Cov();
+  float* GPUrestrict() mP = mT->Par();
+  float chi2Y, chi2Z;
+  if (mT->NDF() <= 0) {
+    chi2Y = CAMath::Square((float)inter->posY - posY) / ((float)inter->errorY + err2Y);
+    chi2Z = CAMath::Square((float)inter->posZ + deltaZ - posZ) / ((float)inter->errorZ + err2Z);
+  } else if (mFitInProjections) {
+    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];
+    // printf("\t%21sInterpo ----- abde artaf%16s Y %8.3f, Z %8.3f (Errors %f <-- (%f, %f) %f <-- (%f, %f))\n", "", "", ImP0, ImP1, sqrtf(ImC0), sqrtf(mC[0]), sqrtf(inter->errorY), sqrtf(ImC2), sqrtf(mC[2]), sqrtf(inter->errorZ));
+    const float Jz0 = posY - ImP0;
+    const float Jz1 = posZ - ImP1;
+    const float Jw0 = 1.f / (ImC0 + err2Y);
+    const float Jw2 = 1.f / (ImC2 + err2Z);
+    chi2Y = Jw0 * Jz0 * Jz0;
+    chi2Z = Jw2 * Jz1 * Jz1;
+  } else {
+    const float Iz0 = inter->posY - mP[0];
+    const float Iz1 = inter->posZ + deltaZ - mP[1];
+    float Iw0 = mC[2] + (float)inter->errorZ;
+    float Iw2 = mC[0] + (float)inter->errorY;
+    float Idet = CAMath::Max(1e-10f, Iw0 * Iw2 - mC[1] * mC[1]);
+    Idet = 1.f / Idet;
+    Iw0 *= Idet;
+    const float Iw1 = mC[1] * Idet;
+    Iw2 *= Idet;
+    const float Ik00 = mC[0] * Iw0 + mC[1] * Iw1;
+    const float Ik01 = mC[0] * Iw1 + mC[1] * Iw2;
+    const float Ik10 = mC[1] * Iw0 + mC[2] * Iw1;
+    const float Ik11 = mC[1] * Iw1 + mC[2] * Iw2;
+    const float ImP0 = mP[0] + Ik00 * Iz0 + Ik01 * Iz1;
+    const float ImP1 = mP[1] + Ik10 * Iz0 + Ik11 * Iz1;
+    const float ImC0 = mC[0] - Ik00 * mC[0] + Ik01 * mC[1];
+    const float ImC1 = mC[1] - Ik10 * mC[0] + Ik11 * mC[1];
+    const float ImC2 = mC[2] - Ik10 * mC[1] + Ik11 * mC[2];
+    const float Jz0 = posY - ImP0;
+    const float Jz1 = posZ - ImP1;
+    float Jw0 = ImC2 + err2Z;
+    float Jw2 = ImC0 + err2Y;
+    float Jdet = CAMath::Max(1e-10f, Jw0 * Jw2 - ImC1 * ImC1);
+    Jdet = 1.f / Jdet;
+    Jw0 *= Jdet;
+    const float Jw1 = ImC1 * Jdet;
+    Jw2 *= Jdet;
+    chi2Y = CAMath::Abs((Jw0 * Jz0 + Jw1 * Jz1) * Jz0);
+    chi2Z = CAMath::Abs((Jw1 * Jz0 + Jw2 * Jz1) * Jz1);
+  }
+  if (RejectCluster(chi2Y * param.rec.tpc.clusterRejectChi2TolleranceY, chi2Z * param.rec.tpc.clusterRejectChi2TolleranceZ, clusterState)) { // TODO: Relative Pt resolution decreases slightly, why?
+    // printf("Reject Cluster chiy2 %f chiz2 %f (Pos Y: %f - %f %f ; Pos Z: %f - %f %f)\n", chi2Y, chi2Z, posY, mP[0], (float)inter->posY, posZ, mP[1], (float)inter->posZ + deltaZ);
+    return updateErrorClusterRejectedInInterpolation;
   }
   return 0;
 }

GPU/GPUTracking/Merger/GPUTPCGMPropagator.h

@@ -57,11 +57,6 @@ class GPUTPCGMPropagator
     updateErrorClusterRejectedInUpdate = 4,
     updateErrorClusterRejectedEdge = 5
   };
-  enum RejectChi2Mode {
-    rejectDirect = 1,
-    rejectInterFill = 2,
-    rejectInterReject = 3
-  };
 
   GPUdDefault() GPUTPCGMPropagator() = default;
 
@@ -104,10 +99,11 @@ class GPUTPCGMPropagator
 
   GPUd() int32_t PropagateToXAlphaBz(float posX, float posAlpha, bool inFlyDirection);
 
-  GPUd() int32_t Update(float posY, float posZ, int32_t iRow, const GPUParam& param, int16_t clusterState, int8_t rejectChi2, bool refit, int8_t sector, float time, float avgInvCharge, float invCharge);
-  GPUd() int32_t Update(float posY, float posZ, int32_t iRow, const GPUParam& param, int16_t clusterState, int8_t rejectChi2, bool refit, float err2Y, float err2Z);
+  GPUd() int32_t Update(float posY, float posZ, int32_t iRow, const GPUParam& param, int16_t clusterState, bool rejectChi2, bool refit, int8_t sector, float time, float avgInvCharge, float invCharge);
+  GPUd() int32_t Update(float posY, float posZ, int32_t iRow, const GPUParam& param, int16_t clusterState, bool rejectChi2, bool refit, float err2Y, float err2Z);
   GPUd() int32_t Update(float posY, float posZ, int16_t clusterState, bool rejectChi2, float err2Y, float err2Z, const GPUParam* param = nullptr);
-  GPUd() int32_t InterpolateReject(const GPUParam& param, float posY, float posZ, int16_t clusterState, int8_t rejectChi2, gputpcgmmergertypes::InterpolationErrorHit* inter, float err2Y, float err2Z, float deltaZ);
+  GPUd() void InterpolateFill(float posY, float posZ, gputpcgmmergertypes::InterpolationErrorHit* inter);
+  GPUd() int32_t InterpolateReject(const GPUParam& param, float posY, float posZ, int16_t clusterState, gputpcgmmergertypes::InterpolationErrorHit* inter, float err2Y, float err2Z, float deltaZ);
   GPUd() float PredictChi2(float posY, float posZ, int32_t iRow, const GPUParam& param, int16_t clusterState, int8_t sideC, float time, float avgCharge, float charge) const;
   GPUd() float PredictChi2(float posY, float posZ, float err2Y, float err2Z) const;
   GPUd() static int32_t RejectCluster(float chiY, float chiZ, uint8_t clusterState)