@@ -79,6 +79,7 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger* GPUrestrict() merger, int32_
7979 float lastUpdateX = -1 .f ;
8080 uint8_t lastRow = 255 ;
8181 uint8_t lastSector = 255 ;
82+ float deltaZ = 0 .f ;
8283
8384 for (int32_t iWay = 0 ; iWay < nWays; iWay++) {
8485 int32_t nMissed = 0 , nMissed2 = 0 ;
@@ -117,7 +118,7 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger* GPUrestrict() merger, int32_
117118 }
118119
119120 if ((param.rec .tpc .trackFitRejectMode > 0 && nMissed >= param.rec .tpc .trackFitRejectMode ) || nMissed2 >= param.rec .tpc .trackFitMaxRowMissedHard || clusters[ihit].state & GPUTPCGMMergedTrackHit::flagReject) {
120- CADEBUG (printf (" \t Skipping hit, %d hits rejected, flag %X\n " , nMissed, (int32_t )clusters[ihit].state ));
121+ CADEBUG (printf (" \t Skipping hit %d , %d hits rejected, flag %X\n " , ihit , nMissed, (int32_t )clusters[ihit].state ));
121122 if (finalOutInFit && !(clusters[ihit].state & GPUTPCGMMergedTrackHit::flagReject)) {
122123 clusters[ihit].state |= GPUTPCGMMergedTrackHit::flagRejectErr;
123124 }
@@ -225,14 +226,25 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger* GPUrestrict() merger, int32_
225226 }
226227 CADEBUG (printf (" \n " ));
227228
228- int32_t retValUpd;
229+ int32_t retValUpd = 0 , retValInt = 0 ;
229230 float threshold = 3 .f + (lastUpdateX >= 0 ? (CAMath::Abs (mX - lastUpdateX) / 2 ) : 0 .f );
230231 if (mNDF > (int32_t )param.rec .tpc .mergerNonInterpolateRejectMinNDF && (CAMath::Abs (yy - mP [0 ]) > threshold || CAMath::Abs (zz - mP [1 ]) > threshold)) {
231232 retValUpd = GPUTPCGMPropagator::updateErrorClusterRejectedDistance;
232233 } else {
233- int8_t rejectChi2 = attempt ? 0 // In second attempt, we do not reject
234- : (param.rec .tpc .mergerInterpolateErrors && CAMath::Abs (ihit - ihitMergeFirst) <= 1 ) ? (finalOutInFit ? (GPUTPCGMPropagator::rejectInterFill + !(iWay & 1 )) : 0 ) // reject via interpolation
235- : (allowChangeClusters && goodRows > 5 ); // normal rejection during the fit
234+ int8_t rejectChi2 = 0 ;
235+ if (attempt == 0 ) {
236+ if (param.rec .tpc .mergerInterpolateErrors && CAMath::Abs (ihit - ihitMergeFirst) <= 1 ) {
237+ if (iWay == nWays - 3 ) {
238+ rejectChi2 = GPUTPCGMPropagator::rejectInterFill;
239+ } else if (iWay == nWays - 2 ) {
240+ rejectChi2 = GPUTPCGMPropagator::rejectInterReject;
241+ } else if (iWay == nWays - 1 ) {
242+ rejectChi2 = (param.rec .tpc .mergerInterpolateRejectAlsoOnCurrentPosition && GetNDF () > (int32_t )param.rec .tpc .mergerNonInterpolateRejectMinNDF ) ? GPUTPCGMPropagator::rejectDirect : 0 ;
243+ }
244+ } else {
245+ rejectChi2 = allowChangeClusters && goodRows > 5 ;
246+ }
247+ }
236248
237249 float err2Y, err2Z;
238250 const float time = merger->GetConstantMem ()->ioPtrs .clustersNative ? merger->GetConstantMem ()->ioPtrs .clustersNative ->clustersLinear [cluster.num ].getTime () : -1 .f ;
@@ -243,18 +255,15 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger* GPUrestrict() merger, int32_
243255
244256 prop.GetErr2 (err2Y, err2Z, param, zz, cluster.row , clusterState, cluster.sector , time, invAvgCharge, invCharge);
245257
246- int retValInt = 0 ;
247258 if (rejectChi2 >= GPUTPCGMPropagator::rejectInterFill) {
248259 if (rejectChi2 == GPUTPCGMPropagator::rejectInterReject && interpolation.hit [ihit].errorY < (GPUCA_PAR_MERGER_INTERPOLATION_ERROR_TYPE_A)0 ) {
249260 rejectChi2 = GPUTPCGMPropagator::rejectDirect;
250261 } else {
251- retValInt = prop.InterpolateReject (param, yy, zz, clusterState, rejectChi2, &interpolation.hit [ihit], err2Y, err2Z);
262+ retValInt = prop.InterpolateReject (param, yy, zz, clusterState, rejectChi2, &interpolation.hit [ihit], err2Y, err2Z, deltaZ );
252263 }
253264 }
254265
255- if (retValInt) {
256- retValUpd = retValInt;
257- } else if (param.rec .tpc .rejectEdgeClustersInTrackFit && uncorrectedY > -1e6f && param.rejectEdgeClusterByY (uncorrectedY, cluster.row , CAMath::Sqrt (mC [0 ]))) { // uncorrectedY > -1e6f implies allowChangeClusters
266+ if (param.rec .tpc .rejectEdgeClustersInTrackFit && uncorrectedY > -1e6f && param.rejectEdgeClusterByY (uncorrectedY, cluster.row , CAMath::Sqrt (mC [0 ]))) { // uncorrectedY > -1e6f implies allowChangeClusters
258267 retValUpd = GPUTPCGMPropagator::updateErrorClusterRejectedEdge;
259268 } else {
260269 retValUpd = prop.Update (yy, zz, cluster.row , param, clusterState, rejectChi2, refit, err2Y, err2Z);
@@ -265,11 +274,11 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger* GPUrestrict() merger, int32_
265274 }
266275 // clang-format off
267276 CADEBUG (if (!CheckCov ()) GPUError (" INVALID COV AFTER UPDATE!!!" ));
268- CADEBUG (printf (" \t %21sFit Alpha %8.3f , X %8.3f - Y %8.3f, Z %8.3f - QPt %7.2f (%7.2f), SP %5.2f (%5.2f), DzDs %5.2f %16s --- Cov sY %8.3f sZ %8.3f sSP %8.3f sPt %8.3f - YPt %8.3f - FErr %d\n " , " " , prop.GetAlpha (), mX , mP [0 ], mP [1 ], mP [4 ], prop.GetQPt0 (), mP [2 ], prop.GetSinPhi0 (), mP [3 ], " " , sqrtf (mC [0 ]), sqrtf (mC [2 ]), sqrtf (mC [5 ]), sqrtf (mC [14 ]), mC [10 ], retValUpd));
277+ CADEBUG (printf (" \t %21sFit Alpha %8.3f , X %8.3f - Y %8.3f, Z %8.3f - QPt %7.2f (%7.2f), SP %5.2f (%5.2f), DzDs %5.2f %16s --- Cov sY %8.3f sZ %8.3f sSP %8.3f sPt %8.3f - YPt %8.3f - FErr %d %d \n " , " " , prop.GetAlpha (), mX , mP [0 ], mP [1 ], mP [4 ], prop.GetQPt0 (), mP [2 ], prop.GetSinPhi0 (), mP [3 ], " " , sqrtf (mC [0 ]), sqrtf (mC [2 ]), sqrtf (mC [5 ]), sqrtf (mC [14 ]), mC [10 ], retValUpd, retValInt ));
269278 // clang-format on
270279
271- ConstrainSinPhi (); // TODO: Limit using ConstrainSinPhi everywhere!
272- if (retValUpd == 0 ) // track is updated
280+ ConstrainSinPhi (); // TODO: Limit using ConstrainSinPhi everywhere!
281+ if (! retValUpd && !retValInt ) // track is updated
273282 {
274283 lastUpdateX = mX ;
275284 covYYUpd = mC [0 ];
@@ -311,14 +320,16 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger* GPUrestrict() merger, int32_
311320 }
312321 }
313322 }
314- } else if (retValUpd >= GPUTPCGMPropagator::updateErrorClusterRejected) { // cluster far away form the track
315- if (allowChangeClusters) {
323+ } else if (retValInt || retValUpd >= GPUTPCGMPropagator::updateErrorClusterRejected) { // cluster far away form the track
324+ if (retValInt || allowChangeClusters) {
316325 MarkClusters (clusters, ihitMergeFirst, ihit, wayDirection, GPUTPCGMMergedTrackHit::flagRejectDistance);
317326 } else if (finalFit) {
318327 MarkClusters (clusters, ihitMergeFirst, ihit, wayDirection, GPUTPCGMMergedTrackHit::flagRejectErr);
319328 }
320- nMissed++;
321- nMissed2++;
329+ if (!retValInt) {
330+ nMissed++;
331+ nMissed2++;
332+ }
322333 } else {
323334 break ; // bad chi2 for the whole track, stop the fit
324335 }
@@ -328,7 +339,9 @@ GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger* GPUrestrict() merger, int32_
328339 CADEBUG (printf (" \t\t STORING %d lastRow %d row %d out %d\n " , iTrk, (int )lastRow, (int )clusters[(iWay & 1 ) ? (maxN - 1 ) : 0 ].row , lastRow > clusters[(iWay & 1 ) ? (maxN - 1 ) : 0 ].row ));
329340 }
330341 if (!(iWay & 1 ) && !finalFit && !track.CCE () && !track.Looper ()) {
331- ShiftZ (clusters, merger, maxN);
342+ deltaZ = ShiftZ (clusters, merger, maxN);
343+ } else {
344+ deltaZ = 0 .f ;
332345 }
333346 }
334347 ConstrainSinPhi ();
@@ -775,21 +788,21 @@ GPUdi() void GPUTPCGMTrackParam::AttachClustersLooper(const GPUTPCGMMerger* GPUr
775788 }
776789}
777790
778- GPUd () void GPUTPCGMTrackParam::ShiftZ(const GPUTPCGMMergedTrackHit* clusters, const GPUTPCGMMerger* merger, int32_t N)
791+ GPUd () float GPUTPCGMTrackParam::ShiftZ(const GPUTPCGMMergedTrackHit* clusters, const GPUTPCGMMerger* merger, int32_t N)
779792{
780793 if (N == 0 ) {
781794 N = 1 ;
782795 }
783796 const auto & GPUrestrict () cls = merger->GetConstantMem ()->ioPtrs .clustersNative ->clustersLinear ;
784797 float z0 = cls[clusters[0 ].num ].getTime (), zn = cls[clusters[N - 1 ].num ].getTime ();
785798 const auto tmp = zn > z0 ? std::array<float , 3 >{zn, z0, GPUTPCGeometry::Row2X (clusters[N - 1 ].row )} : std::array<float , 3 >{z0, zn, GPUTPCGeometry::Row2X (clusters[0 ].row )};
786- ShiftZ (merger, clusters[0 ].sector , tmp[0 ], tmp[1 ], tmp[2 ]);
799+ return ShiftZ (merger, clusters[0 ].sector , tmp[0 ], tmp[1 ], tmp[2 ]);
787800}
788801
789- GPUd () void GPUTPCGMTrackParam::ShiftZ(const GPUTPCGMMerger* GPUrestrict () merger, int32_t sector, float cltmax, float cltmin, float clx)
802+ GPUd () float GPUTPCGMTrackParam::ShiftZ(const GPUTPCGMMerger* GPUrestrict () merger, int32_t sector, float cltmax, float cltmin, float clx)
790803{
791804 if (!merger->Param ().par .continuousTracking ) {
792- return ;
805+ return 0 . f ;
793806 }
794807 float deltaZ = 0 .f ;
795808 bool beamlineReached = false ;
@@ -828,7 +841,6 @@ GPUd() void GPUTPCGMTrackParam::ShiftZ(const GPUTPCGMMerger* GPUrestrict() merge
828841 {
829842 float deltaT = merger->GetConstantMem ()->calibObjects .fastTransformHelper ->getCorrMap ()->convDeltaZtoDeltaTimeInTimeFrame (sector, deltaZ);
830843 mTOffset += deltaT;
831- mP [1 ] -= deltaZ;
832844 const float maxT = cltmin - merger->GetConstantMem ()->calibObjects .fastTransformHelper ->getCorrMap ()->getT0 ();
833845 const float minT = cltmax - merger->GetConstantMem ()->calibObjects .fastTransformHelper ->getCorrMap ()->getMaxDriftTime (sector);
834846 // printf("T Check: Clusters %f %f, min %f max %f vtx %f\n", tz1, tz2, minT, maxT, mTOffset);
@@ -840,13 +852,14 @@ GPUd() void GPUTPCGMTrackParam::ShiftZ(const GPUTPCGMMerger* GPUrestrict() merge
840852 deltaT = maxT - mTOffset ;
841853 }
842854 if (deltaT != 0 .f ) {
843- deltaZ = merger->GetConstantMem ()->calibObjects .fastTransformHelper ->getCorrMap ()->convDeltaTimeToDeltaZinTimeFrame (sector, deltaT);
855+ deltaZ + = merger->GetConstantMem ()->calibObjects .fastTransformHelper ->getCorrMap ()->convDeltaTimeToDeltaZinTimeFrame (sector, deltaT);
844856 // 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);
845857 mTOffset += deltaT;
846- mP [1 ] -= deltaZ;
847858 }
859+ mP [1 ] -= deltaZ;
848860 }
849861 // printf("\n");
862+ return -deltaZ;
850863}
851864
852865GPUd () bool GPUTPCGMTrackParam::CheckCov() const
@@ -861,28 +874,23 @@ GPUd() bool GPUTPCGMTrackParam::CheckNumericalQuality(float overrideCovYY) const
861874{
862875 // * Check that the track parameters and covariance matrix are reasonable
863876 bool ok = CAMath::Finite (mX ) && CAMath::Finite (mChi2 );
864- CADEBUG (printf (" OK %d - %f - " , (int32_t )ok, mX ); for (int32_t i = 0 ; i < 5 ; i++) { printf (" %f " , mP [i]); } printf (" - " ); for (int32_t i = 0 ; i < 15 ; i++) { printf (" %f " , mC [i]); } printf (" \n " ));
877+ // CADEBUG(printf("OK %d - %f - ", (int32_t)ok, mX); for (int32_t i = 0; i < 5; i++) { printf("%f ", mP[i]); } printf(" - "); for (int32_t i = 0; i < 15; i++) { printf("%f ", mC[i]); } printf("\n"));
865878 const float * c = mC ;
866879 for (int32_t i = 0 ; i < 15 ; i++) {
867880 ok = ok && CAMath::Finite (c[i]);
868881 }
869- CADEBUG (printf (" OK1 %d\n " , (int32_t )ok));
870882 for (int32_t i = 0 ; i < 5 ; i++) {
871883 ok = ok && CAMath::Finite (mP [i]);
872884 }
873- CADEBUG (printf (" OK2 %d\n " , (int32_t )ok));
874885 if ((overrideCovYY > 0 ? overrideCovYY : c[0 ]) > 4 .f * 4 .f || c[2 ] > 4 .f * 4 .f || c[5 ] > 2 .f * 2 .f || c[9 ] > 2 .f * 2 .f ) {
875886 ok = 0 ;
876887 }
877- CADEBUG (printf (" OK3 %d\n " , (int32_t )ok));
878888 if (CAMath::Abs (mP [2 ]) > GPUCA_MAX_SIN_PHI) {
879889 ok = 0 ;
880890 }
881- CADEBUG (printf (" OK4 %d\n " , (int32_t )ok));
882891 if (!CheckCov ()) {
883892 ok = false ;
884893 }
885- CADEBUG (printf (" OK5 %d\n " , (int32_t )ok));
886894 return ok;
887895}
888896
@@ -903,7 +911,7 @@ GPUdii() void GPUTPCGMTrackParam::RefitTrack(GPUTPCGMMergedTrack& GPUrestrict()
903911 CADEBUG (int32_t nTrackHitsOld = nTrackHits; float ptOld = t.QPt ());
904912 bool ok = t.Fit (merger, iTrk, merger->Clusters () + track.FirstClusterRef (), nTrackHits, NTolerated, Alpha, attempt, GPUCA_MAX_SIN_PHI, track);
905913 CADEBUG (printf (" Finished Fit Track %d\n " , iTrk));
906- CADEBUG (printf (" OUTPUT hits %d -> %d+%d = %d, QPt %f -> %f, SP %f, ok %d chi2 %f chi2ndf %f\n " , nTrackHitsOld, nTrackHits, NTolerated, nTrackHits + NTolerated, ptOld, t.QPt (), t.SinPhi (), (int32_t )ok, t.Chi2 (), t.Chi2 () / CAMath::Max (1 , nTrackHits)));
914+ CADEBUG (printf (" OUTPUT hits %d -> %d+%d = %d, QPt %f -> %f, SP %f, OK %d chi2 %f chi2ndf %f\n " , nTrackHitsOld, nTrackHits, NTolerated, nTrackHits + NTolerated, ptOld, t.QPt (), t.SinPhi (), (int32_t )ok, t.Chi2 (), t.Chi2 () / CAMath::Max (1 , nTrackHits)));
907915
908916 if (!ok && attempt == 0 && merger->Param ().rec .tpc .retryRefit ) {
909917 for (uint32_t i = 0 ; i < track.NClusters (); i++) {
0 commit comments