217 CADEBUG(int32_t ii; printf(
"\nRefitting track\n"));
218 typename internal::refitTrackTypes<S>::propagator prop;
220 float TrackParCovChi2 = 0.f;
221 convertTrack<S, T, typename internal::refitTrackTypes<S>::propagator>(trk,
trkX, prop, &TrackParCovChi2);
224 [[maybe_unused]] int32_t eLossSign = 0;
225 if constexpr (std::is_same_v<T, GPUTPCGMMergedTrack>) {
227 tOffset =
trkX.GetParam().GetTOffset();
228 }
else if constexpr (std::is_same_v<T, TrackTPC>) {
230 tOffset =
trkX.getTime0();
231 }
else if constexpr (std::is_same_v<T, TrackParCovWithArgs>) {
233 tOffset =
trkX.time0;
234 eLossSign =
trkX.eLossSign;
236 static_assert(
"Invalid template");
238 if constexpr (std::is_same_v<S, GPUTPCGMTrackParam>) {
239 CADEBUG(printf(
"\t%21sInit Alpha %8.3f , X %8.3f - Y %8.3f, Z %8.3f - QPt %7.2f (%7.2f), SP %5.2f (%5.2f) --- Cov sY %8.3f sZ %8.3f sSP %8.3f sPt %8.3f\n",
"", prop.GetAlpha(), trk.GetX(), trk.Par()[0], trk.Par()[1], trk.Par()[4], prop.GetQPt0(), trk.Par()[2], prop.GetSinPhi0(), sqrtf(trk.Cov()[0]), sqrtf(trk.Cov()[2]), sqrtf(trk.Cov()[5]), sqrtf(trk.Cov()[14])));
242 int32_t direction = outward ? -1 : 1;
246 uint8_t sector = 255,
row = 255;
247 int32_t lastSector = -1, currentSector = -1, currentRow = -1;
248 int16_t clusterState = 0, nextState = 0;
250 float sumInvSqrtCharge = 0.f;
251 int32_t nAvgCharge = 0;
253 for (int32_t
i =
start;
i != stop;
i += cl ? 0 : direction) {
255 float time = 0.f, invCharge = 0.f, invSqrtCharge = 0.f;
260 if constexpr (std::is_same_v<T, GPUTPCGMMergedTrack>) {
261 const auto& hit = mPtrackHits[
trkX.FirstClusterRef() +
i];
262 cl = &mPclusterNative->clustersLinear[hit.num];
265 if (
i + direction != stop) {
273 nextState = mPclusterState[hit.num];
274 }
else if constexpr (std::is_same_v<T, TrackTPC>) {
275 cl = &
trkX.getCluster(mPtrackHitReferences,
i, *mPclusterNative, sector,
row);
276 nextState = mPclusterState[cl - mPclusterNative->clustersLinear];
277 }
else if constexpr (std::is_same_v<T, TrackParCovWithArgs>) {
278 cl = &TrackTPC::getCluster(mPtrackHitReferences,
i, *mPclusterNative, sector,
row,
trkX.clusRef);
279 nextState = mPclusterState[cl - mPclusterNative->clustersLinear];
281 static_assert(
"Invalid template");
284 if (
clusters == 0 || (
row == currentRow && sector == currentSector)) {
291 mPfastTransform->Transform(sector,
row, cl->getPad(), cl->getTime(),
x,
y,
z, tOffset);
292 CADEBUG(printf(
"\tHit %3d/%3d Row %3d: Cluster Alpha %8.3f %3d, X %8.3f - Y %8.3f, Z %8.3f - State %d\n", ii,
count,
row, mPparam->Alpha(sector), (int32_t)sector,
x,
y,
z, (int32_t)nextState));
294 currentSector = sector;
296 clusterState = nextState;
297 time = cl->getTime();
298 invSqrtCharge = CAMath::InvSqrt(cl->
qMax);
299 invCharge = (1.f / cl->
qMax);
302 mPfastTransform->Transform(sector,
row, cl->getPad(), cl->getTime(), xx, yy, zz, tOffset);
303 CADEBUG(printf(
"\tHit %3d/%3d Row %3d: Cluster Alpha %8.3f %3d, X %8.3f - Y %8.3f, Z %8.3f - State %d\n", ii,
count,
row, mPparam->Alpha(sector), (int32_t)sector, xx, yy, zz, (int32_t)nextState));
304 x += xx * cl->getQtot();
305 y += yy * cl->getQtot();
306 z += zz * cl->getQtot();
308 clusterState |= nextState;
312 if (
i + direction != stop) {
325 CADEBUG(printf(
"\tMerged Hit Row %3d: Cluster Alpha %8.3f %3d, X %8.3f - Y %8.3f, Z %8.3f - State %d\n",
row, mPparam->Alpha(sector), (int32_t)sector,
x,
y,
z, (int32_t)clusterState));
328 float invAvgCharge = (sumInvSqrtCharge += invSqrtCharge) / ++nAvgCharge;
329 invAvgCharge *= invAvgCharge;
331 if constexpr (std::is_same_v<S, GPUTPCGMTrackParam>) {
332 if (prop.PropagateToXAlpha(
x, mPparam->Alpha(currentSector), !outward)) {
337 trk.ResetCovariance();
338 }
else if (lastSector != -1 && (lastSector < 18) != (sector < 18)) {
339 if (mPparam->rec.tpc.addErrorsCECrossing) {
340 if (mPparam->rec.tpc.addErrorsCECrossing >= 2) {
341 trk.AddCovDiagErrorsWithCorrelations(mPparam->rec.tpc.errorsCECrossing);
343 trk.AddCovDiagErrors(mPparam->rec.tpc.errorsCECrossing);
345 }
else if (trk.Cov()[2] < 0.5f) {
349 CADEBUG(printf(
"\t%21sPropaga Alpha %8.3f , X %8.3f - Y %8.3f, Z %8.3f - QPt %7.2f (%7.2f), SP %5.2f (%5.2f) --- Res %8.3f %8.3f --- Cov sY %8.3f sZ %8.3f sSP %8.3f sPt %8.3f - YPt %8.3f\n",
"", prop.GetAlpha(),
x, trk.Par()[0], trk.Par()[1], trk.Par()[4], prop.GetQPt0(), trk.Par()[2], prop.GetSinPhi0(), trk.Par()[0] -
y, trk.Par()[1] -
z, sqrtf(trk.Cov()[0]), sqrtf(trk.Cov()[2]), sqrtf(trk.Cov()[5]), sqrtf(trk.Cov()[14]), trk.Cov()[10]));
351 if (prop.Update(
y,
z,
row, *mPparam, clusterState, 0,
true, sector,
time, invAvgCharge, invCharge)) {
355 trk.ConstrainSinPhi();
356 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\n",
"", prop.GetAlpha(),
x, trk.Par()[0], trk.Par()[1], trk.Par()[4], prop.GetQPt0(), trk.Par()[2], prop.GetSinPhi0(), trk.Par()[3],
"", sqrtf(trk.Cov()[0]), sqrtf(trk.Cov()[2]), sqrtf(trk.Cov()[5]), sqrtf(trk.Cov()[14]), trk.Cov()[10]));
357 }
else if constexpr (std::is_same_v<S, TrackParCov>) {
358 if (!trk.rotate(mPparam->Alpha(currentSector))) {
362 if (!prop->PropagateToXBxByBz(trk,
x, constants::MAX_SIN_PHI_LOW, Propagator::MAX_STEP, Propagator::MatCorrType::USEMatCorrLUT,
nullptr, eLossSign)) {
366 if (lastSector != -1 && (lastSector < 18) != (sector < 18)) {
367 if (mPparam->rec.tpc.addErrorsCECrossing) {
368 trk.updateCov(mPparam->rec.tpc.errorsCECrossing, mPparam->rec.tpc.addErrorsCECrossing >= 2);
369 }
else if (trk.getCov()[2] < 0.5f) {
374 trk.resetCovariance();
375 float bzkG = prop->getNominalBz(), qptB5Scale = CAMath::Abs(bzkG) > 0.1f ? CAMath::Abs(bzkG) / 5.006680f : 1.f;
376 float q2pt2 = trk.getQ2Pt() * trk.getQ2Pt(), q2pt2Wgh = q2pt2 * qptB5Scale * qptB5Scale;
377 float err2 = (100.f + q2pt2Wgh) / (1.f + q2pt2Wgh) * q2pt2;
378 trk.setCov(err2, 14);
379 TrackParCovChi2 = 0.f;
381 CADEBUG(printf(
"\t%21sPropaga Alpha %8.3f , X %8.3f - Y %8.3f, Z %8.3f - QPt %7.2f (%7.2f), SP %5.2f (%5.2f) --- Res %8.3f %8.3f --- Cov sY %8.3f sZ %8.3f sSP %8.3f sPt %8.3f - YPt %8.3f\n",
"", trk.getAlpha(),
x, trk.getParams()[0], trk.getParams()[1], trk.getParams()[4], trk.getParams()[4], trk.getParams()[2], trk.getParams()[2], trk.getParams()[0] -
y, trk.getParams()[1] -
z, sqrtf(trk.getCov()[0]), sqrtf(trk.getCov()[2]), sqrtf(trk.getCov()[5]), sqrtf(trk.getCov()[14]), trk.getCov()[10]));
382 std::array<float, 2> p = {
y,
z};
383 std::array<float, 3>
c = {0, 0, 0};
384 GPUTPCGMPropagator::GetErr2(
c[0],
c[2], *mPparam, getPar(trk)[2], getPar(trk)[3],
z,
x,
y, currentRow, clusterState, sector,
time, invAvgCharge, invCharge,
false);
385 TrackParCovChi2 += trk.getPredictedChi2(p,
c);
386 if (!trk.update(p,
c)) {
390 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\n",
"", trk.getAlpha(),
x, trk.getParams()[0], trk.getParams()[1], trk.getParams()[4], trk.getParams()[4], trk.getParams()[2], trk.getParams()[2], trk.getParams()[3],
"", sqrtf(trk.getCov()[0]), sqrtf(trk.getCov()[2]), sqrtf(trk.getCov()[5]), sqrtf(trk.getCov()[14]), trk.getCov()[10]));
392 static_assert(
"Invalid template");
397 if constexpr (std::is_same_v<S, GPUTPCGMTrackParam>) {
398 float alpha = prop.GetAlpha();
399 trk.MoveToReference(prop, *mPparam,
alpha);
400 trk.NormalizeAlpha(
alpha);
401 prop.SetAlpha(
alpha);
402 }
else if constexpr (std::is_same_v<S, TrackParCov>) {
403 constexpr float kDeg2Rad = M_PI / 180.f;
404 constexpr float kSectAngle = 2 * M_PI / 18.f;
405 if (mPparam->rec.tpc.trackReferenceX <= 500) {
407 const int32_t eLossSignRef = (mPparam->rec.tpc.trackReferenceX > trk.getX()) == outward ? eLossSign : -eLossSign;
408 if (prop->PropagateToXBxByBz(trk, mPparam->rec.tpc.trackReferenceX, Propagator::MAX_SIN_PHI, Propagator::MAX_STEP, Propagator::MatCorrType::USEMatCorrLUT,
nullptr, eLossSignRef)) {
409 if (CAMath::Abs(trk.getY()) > trk.getX() * CAMath::Tan(kSectAngle / 2.f)) {
410 float newAlpha = trk.getAlpha() + CAMath::Round(CAMath::ATan2(trk.getY(), trk.getX()) / kDeg2Rad / 20.f) * kSectAngle;
411 GPUTPCGMTrackParam::NormalizeAlpha(newAlpha);
412 trk.rotate(newAlpha) && prop->PropagateToXBxByBz(trk, mPparam->rec.tpc.trackReferenceX, Propagator::MAX_SIN_PHI, Propagator::MAX_STEP, Propagator::MatCorrType::USEMatCorrLUT,
nullptr, eLossSignRef);
417 static_assert(
"Invalid template");
420 convertTrack<T, S, typename internal::refitTrackTypes<S>::propagator>(
trkX, trk, prop, &TrackParCovChi2);