Project
Loading...
Searching...
No Matches
GPUTPCGMTrackParam.cxx
Go to the documentation of this file.
1// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
2// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
3// All rights not expressly granted are reserved.
4//
5// This software is distributed under the terms of the GNU General Public
6// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
7//
8// In applying this license CERN does not waive the privileges and immunities
9// granted to it by virtue of its status as an Intergovernmental Organization
10// or submit itself to any jurisdiction.
11
14
15#define GPUCA_CADEBUG 0
16#define DEBUG_SINGLE_TRACK -1
17#define EXTRACT_RESIDUALS 0
18
19#if EXTRACT_RESIDUALS == 1
20#include "GPUROOTDump.h"
21#endif
22
23#include "GPUTPCDef.h"
24#include "GPUTPCGMTrackParam.h"
26#include "GPUTPCGMPropagator.h"
27#include "GPUTPCGMBorderTrack.h"
28#include "GPUTPCGMMergedTrack.h"
30#include "GPUTPCGMMerger.h"
31#include "GPUTPCTracker.h"
32#include "GPUTPCClusterData.h"
33#include "GPUdEdx.h"
34#include "GPUParam.h"
35#include "GPUO2DataTypes.h"
36#include "GPUConstantMem.h"
37#include "TPCFastTransform.h"
39#include "GPUTPCConvertImpl.h"
40#include "GPUTPCGMMergerTypes.h"
41#include "GPUParam.inc"
42#include "GPUGetConstexpr.h"
43
44#ifdef GPUCA_CADEBUG_ENABLED
45#include "GPUSettings.h"
47#endif
48
49#ifndef GPUCA_GPUCODE_DEVICE
50#include <cmath>
51#include <cstdlib>
52#endif
53
54using namespace o2::gpu;
55using namespace o2::tpc;
56
57GPUd() bool GPUTPCGMTrackParam::Fit(GPUTPCGMMerger* GPUrestrict() merger, int32_t iTrk, GPUTPCGMMergedTrackHit* GPUrestrict() clusters, GPUTPCGMMergedTrackHitXYZ* GPUrestrict() clustersXYZ, int32_t& GPUrestrict() N, int32_t& GPUrestrict() NTolerated, float& GPUrestrict() Alpha, int32_t attempt, float maxSinPhi, gputpcgmmergertypes::GPUTPCOuterParam* GPUrestrict() outerParam)
58{
59 static constexpr float kDeg2Rad = M_PI / 180.f;
60 CADEBUG(static constexpr float kSectAngle = 2 * M_PI / 18.f);
61
62 const GPUParam& GPUrestrict() param = merger->Param();
63
64 GPUdEdx dEdx, dEdxAlt;
67 prop.SetMaterialTPC();
68 prop.SetPolynomialField(&param.polynomialField);
69 prop.SetMaxSinPhi(maxSinPhi);
70 if ((clusters[0].sector < 18) == (clusters[N - 1].sector < 18)) {
71 ShiftZ2(clusters, clustersXYZ, merger, N);
72 }
73 if (param.rec.tpc.mergerInterpolateErrors) {
74 for (int32_t i = 0; i < N; i++) {
75 interpolation.hit[i].errorY = -1;
76 }
77 }
78
79 const int32_t nWays = param.rec.tpc.nWays;
80 const int32_t maxN = N;
81 int32_t ihitStart = 0;
82 float covYYUpd = 0.f;
83 float lastUpdateX = -1.f;
84 uint8_t lastRow = 255;
85 uint8_t lastSector = 255;
86 uint8_t storeOuter = 0;
87
88 for (int32_t iWay = 0; iWay < nWays; iWay++) {
89 int32_t nMissed = 0, nMissed2 = 0;
90 float sumInvSqrtCharge = 0.f;
91 int32_t nAvgCharge = 0;
92
93 if (iWay && storeOuter != 255 && param.rec.tpc.nWaysOuter && outerParam) {
94 storeOuter = 0;
95 if (iWay == nWays - 1) {
96 StoreOuter(outerParam, prop, 0);
97 if (merger->MergedTracks()[iTrk].Looper()) {
98 storeOuter = 1;
99 }
100 } else if (iWay == nWays - 2 && merger->MergedTracks()[iTrk].Looper()) {
101 storeOuter = 2;
102 }
103 }
104
105 int32_t resetT0 = initResetT0();
106 const bool refit = (nWays == 1 || iWay >= 1);
107 const float maxSinForUpdate = CAMath::Sin(70.f * kDeg2Rad);
108
109 ResetCovariance();
110 prop.SetSeedingErrors(!(refit && attempt == 0));
111 prop.SetFitInProjections(param.rec.fitInProjections == -1 ? (iWay != 0) : param.rec.fitInProjections);
112 prop.SetPropagateBzOnly(param.rec.fitPropagateBzOnly > iWay);
113 prop.SetMatLUT((param.rec.useMatLUT && iWay == nWays - 1) ? merger->GetConstantMem()->calibObjects.matLUT : nullptr);
114 prop.SetTrack(this, iWay ? prop.GetAlpha() : Alpha);
115 ConstrainSinPhi(prop.GetFitInProjections() ? 0.95f : GPUCA_MAX_SIN_PHI_LOW);
116 CADEBUG(printf("Fitting track %d way %d (sector %d, alpha %f)\n", iTrk, iWay, CAMath::Float2IntRn(prop.GetAlpha() / kSectAngle) + (mP[1] < 0 ? 18 : 0), prop.GetAlpha()));
117
118 N = 0;
119 lastUpdateX = -1;
120 const bool inFlyDirection = iWay & 1;
121 uint8_t lastLeg = clusters[ihitStart].leg;
122 const int32_t wayDirection = (iWay & 1) ? -1 : 1;
123
124 bool noFollowCircle = false, noFollowCircle2 = false;
125 int32_t goodRows = 0;
126 for (int32_t ihit = ihitStart; ihit >= 0 && ihit < maxN; ihit += wayDirection) {
127 const bool crossCE = lastSector != 255 && ((lastSector < 18) ^ (clusters[ihit].sector < 18));
128 if (crossCE) {
129 lastSector = clusters[ihit].sector;
130 noFollowCircle2 = true;
131 }
132
133 if (storeOuter == 2 && clusters[ihit].leg == clusters[maxN - 1].leg - 1) {
134 if (lastLeg == clusters[maxN - 1].leg) {
135 StoreOuter(outerParam, prop, 1);
136 storeOuter = 255;
137 } else {
138 storeOuter = 0;
139 }
140 }
141
142 if ((param.rec.tpc.trackFitRejectMode > 0 && nMissed >= param.rec.tpc.trackFitRejectMode) || nMissed2 >= param.rec.tpc.trackFitMaxRowMissedHard || clusters[ihit].state & GPUTPCGMMergedTrackHit::flagReject) {
143 CADEBUG(printf("\tSkipping hit, %d hits rejected, flag %X\n", nMissed, (int32_t)clusters[ihit].state));
144 if (iWay + 2 >= nWays && !(clusters[ihit].state & GPUTPCGMMergedTrackHit::flagReject)) {
146 }
147 continue;
148 }
149
150 const bool allowModification = refit && (iWay == 0 || (((nWays - iWay) & 1) ? (ihit >= CAMath::Min(maxN / 2, 30)) : (ihit <= CAMath::Max(maxN / 2, maxN - 30))));
151
152 int32_t ihitMergeFirst = ihit;
153 uint8_t clusterState = clusters[ihit].state;
154 const float clAlpha = param.Alpha(clusters[ihit].sector);
155 float xx, yy, zz;
156 if (param.par.earlyTpcTransform) {
157 const float zOffset = (clusters[ihit].sector < 18) == (clusters[0].sector < 18) ? mTZOffset : -mTZOffset;
158 xx = clustersXYZ[ihit].x;
159 yy = clustersXYZ[ihit].y;
160 zz = clustersXYZ[ihit].z - zOffset;
161 } else {
162 const ClusterNative& GPUrestrict() cl = merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[ihit].num];
163 merger->GetConstantMem()->calibObjects.fastTransformHelper->Transform(clusters[ihit].sector, clusters[ihit].row, cl.getPad(), cl.getTime(), xx, yy, zz, mTZOffset);
164 }
165 // clang-format off
166 CADEBUG(printf("\tHit %3d/%3d Row %3d: Cluster Alpha %8.3f %3d, X %8.3f - Y %8.3f, Z %8.3f (Missed %d)\n", ihit, maxN, (int32_t)clusters[ihit].row, clAlpha, (int32_t)clusters[ihit].sector, xx, yy, zz, nMissed));
167 // CADEBUG(if ((uint32_t)merger->GetTrackingChain()->mIOPtrs.nMCLabelsTPC > clusters[ihit].num))
168 // CADEBUG({printf(" MC:"); for (int32_t i = 0; i < 3; i++) {int32_t mcId = merger->GetTrackingChain()->mIOPtrs.mcLabelsTPC[clusters[ihit].num].fClusterID[i].fMCID; if (mcId >= 0) printf(" %d", mcId); } } printf("\n"));
169 // clang-format on
170 if (MergeDoubleRowClusters(ihit, wayDirection, clusters, clustersXYZ, merger, prop, xx, yy, zz, maxN, clAlpha, clusterState, allowModification) == -1) {
171 nMissed++;
172 nMissed2++;
173 continue;
174 }
175
176 if (param.rec.tpc.rejectIFCLowRadiusCluster) {
177 const float r2 = xx * xx + yy * yy;
178 const float rmax = (83.5f + param.rec.tpc.sysClusErrorMinDist);
179 if (r2 < rmax * rmax) {
180 MarkClusters(clusters, ihitMergeFirst, ihit, wayDirection, GPUTPCGMMergedTrackHit::flagRejectErr);
181 }
182 }
183
184 const auto& cluster = clusters[ihit];
185
186 bool changeDirection = (cluster.leg - lastLeg) & 1;
187 // clang-format off
188 CADEBUG(if (changeDirection) printf("\t\tChange direction\n"));
189 CADEBUG(printf("\tLeg %3d Sector %2d %4sTrack Alpha %8.3f %s, X %8.3f - Y %8.3f, Z %8.3f - QPt %7.2f (%7.2f), SP %5.2f (%5.2f) %28s --- Cov sY %8.3f sZ %8.3f sSP %8.3f sPt %8.3f - YPt %8.3f\n", (int32_t)cluster.leg, (int32_t)cluster.sector, "", prop.GetAlpha(), (CAMath::Abs(prop.GetAlpha() - clAlpha) < 0.01 ? " " : " R!"), mX, mP[0], mP[1], mP[4], prop.GetQPt0(), mP[2], prop.GetSinPhi0(), "", sqrtf(mC[0]), sqrtf(mC[2]), sqrtf(mC[5]), sqrtf(mC[14]), mC[10]));
190 // clang-format on
191 if (allowModification && changeDirection && !noFollowCircle && !noFollowCircle2) {
192 bool tryFollow = lastRow != 255;
193 if (tryFollow) {
194 const GPUTPCGMTrackParam backup = *this;
195 const float backupAlpha = prop.GetAlpha();
196 if (FollowCircle<0>(merger, prop, lastSector, lastRow, iTrk, clAlpha, xx, yy, cluster.sector, cluster.row, inFlyDirection)) {
197 CADEBUG(printf("Error during follow circle, resetting track!\n"));
198 *this = backup;
199 prop.SetTrack(this, backupAlpha);
200 noFollowCircle = true;
201 tryFollow = false;
202 }
203 }
204 if (tryFollow) {
205 MirrorTo(prop, yy, zz, inFlyDirection, param, cluster.row, clusterState, false, cluster.sector);
206 lastUpdateX = mX;
207 lastLeg = cluster.leg;
208 lastSector = cluster.sector;
209 lastRow = 255;
210 N++;
211 resetT0 = initResetT0();
212 // clang-format off
213 CADEBUG(printf("\n"));
214 CADEBUG(printf("\t%21sMirror Alpha %8.3f , X %8.3f - Y %8.3f, Z %8.3f - QPt %7.2f (%7.2f), SP %5.2f (%5.2f) %28s --- Cov sY %8.3f sZ %8.3f sSP %8.3f sPt %8.3f - YPt %8.3f\n", "", prop.GetAlpha(), mX, mP[0], mP[1], mP[4], prop.GetQPt0(), mP[2], prop.GetSinPhi0(), "", sqrtf(mC[0]), sqrtf(mC[2]), sqrtf(mC[5]), sqrtf(mC[14]), mC[10]));
215 // clang-format on
216 continue;
217 }
218 } else if (allowModification && lastRow != 255 && CAMath::Abs(cluster.row - lastRow) > 1) {
220 bool dodEdx = param.dodEdxEnabled && param.rec.tpc.adddEdxSubThresholdClusters && iWay == nWays - 1 && CAMath::Abs(cluster.row - lastRow) == 2 && cluster.leg == clusters[maxN - 1].leg;
221 dodEdx = AttachClustersPropagate(merger, cluster.sector, lastRow, cluster.row, iTrk, cluster.leg == clusters[maxN - 1].leg, prop, inFlyDirection, GPUCA_MAX_SIN_PHI, dodEdx);
222 if (dodEdx) {
223 dEdx.fillSubThreshold(lastRow - wayDirection);
224 if GPUCA_RTC_CONSTEXPR (GPUCA_GET_CONSTEXPR(param.rec.tpc, dEdxClusterRejectionFlagMask) != GPUCA_GET_CONSTEXPR(param.rec.tpc, dEdxClusterRejectionFlagMaskAlt)) {
225 dEdxAlt.fillSubThreshold(lastRow - wayDirection);
226 }
227 }
228 }
229 }
230
231 int32_t err = prop.PropagateToXAlpha(xx, clAlpha, inFlyDirection);
232 // clang-format off
233 CADEBUG(if (!CheckCov()){printf("INVALID COV AFTER PROPAGATE!!!\n");});
234 // clang-format on
235 if (err == -2) // Rotation failed, try to bring to new x with old alpha first, rotate, and then propagate to x, alpha
236 {
237 CADEBUG(printf("REROTATE\n"));
238 if (prop.PropagateToXAlpha(xx, prop.GetAlpha(), inFlyDirection) == 0) {
239 err = prop.PropagateToXAlpha(xx, clAlpha, inFlyDirection);
240 }
241 }
242 if (lastRow == 255 || CAMath::Abs((int32_t)lastRow - (int32_t)cluster.row) > 5 || lastSector != cluster.sector || (param.rec.tpc.trackFitRejectMode < 0 && -nMissed <= param.rec.tpc.trackFitRejectMode)) {
243 goodRows = 0;
244 } else {
245 goodRows++;
246 }
247 if (err == 0) {
248 lastRow = cluster.row;
249 lastSector = cluster.sector;
250 }
251 // clang-format off
252 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 - Err %d", "", prop.GetAlpha(), mX, mP[0], mP[1], mP[4], prop.GetQPt0(), mP[2], prop.GetSinPhi0(), mP[0] - yy, mP[1] - zz, sqrtf(mC[0]), sqrtf(mC[2]), sqrtf(mC[5]), sqrtf(mC[14]), mC[10], err));
253 // clang-format on
254
255 if (crossCE) {
256 if (param.rec.tpc.addErrorsCECrossing) {
257 if (param.rec.tpc.addErrorsCECrossing >= 2) {
258 AddCovDiagErrorsWithCorrelations(param.rec.tpc.errorsCECrossing);
259 } else {
260 AddCovDiagErrors(param.rec.tpc.errorsCECrossing);
261 }
262 } else if (mC[2] < 0.5f) {
263 mC[2] = 0.5f;
264 }
265 }
266
267 if (err == 0 && changeDirection) {
268 const float mirrordY = prop.GetMirroredYTrack();
269 CADEBUG(printf(" -- MirroredY: %f --> %f", mP[0], mirrordY));
270 if (CAMath::Abs(yy - mP[0]) > CAMath::Abs(yy - mirrordY)) {
271 CADEBUG(printf(" - Mirroring!!!"));
272 if (allowModification) {
273 AttachClustersMirror<0>(merger, cluster.sector, cluster.row, iTrk, yy, prop); // TODO: Never true, will always call FollowCircle above, really???
274 }
275 MirrorTo(prop, yy, zz, inFlyDirection, param, cluster.row, clusterState, true, cluster.sector);
276 noFollowCircle = false;
277
278 lastUpdateX = mX;
279 lastLeg = cluster.leg;
280 lastRow = 255;
281 N++;
282 resetT0 = initResetT0();
283 // clang-format off
284 CADEBUG(printf("\n"));
285 CADEBUG(printf("\t%21sMirror Alpha %8.3f , X %8.3f - Y %8.3f, Z %8.3f - QPt %7.2f (%7.2f), SP %5.2f (%5.2f) %28s --- Cov sY %8.3f sZ %8.3f sSP %8.3f sPt %8.3f - YPt %8.3f\n", "", prop.GetAlpha(), mX, mP[0], mP[1], mP[4], prop.GetQPt0(), mP[2], prop.GetSinPhi0(), "", sqrtf(mC[0]), sqrtf(mC[2]), sqrtf(mC[5]), sqrtf(mC[14]), mC[10]));
286 // clang-format on
287 continue;
288 }
289 }
290
291 float uncorrectedY = -1e6f;
292 if (allowModification) {
293 uncorrectedY = AttachClusters(merger, cluster.sector, cluster.row, iTrk, cluster.leg == clusters[maxN - 1].leg, prop);
294 }
295
296 const int32_t err2 = mNDF > 0 && CAMath::Abs(prop.GetSinPhi0()) >= maxSinForUpdate;
297 if (err || err2) {
298 if (mC[0] > param.rec.tpc.trackFitCovLimit || mC[2] > param.rec.tpc.trackFitCovLimit) {
299 break;
300 }
301 MarkClusters(clusters, ihitMergeFirst, ihit, wayDirection, GPUTPCGMMergedTrackHit::flagHighIncl);
302 nMissed2++;
303 NTolerated++;
304 CADEBUG(printf(" --- break (%d, %d)\n", err, err2));
305 continue;
306 }
307 CADEBUG(printf("\n"));
308
309 int32_t retVal;
310 float threshold = 3.f + (lastUpdateX >= 0 ? (CAMath::Abs(mX - lastUpdateX) / 2) : 0.f);
311 if (mNDF > 5 && (CAMath::Abs(yy - mP[0]) > threshold || CAMath::Abs(zz - mP[1]) > threshold)) {
313 } else {
314 int8_t rejectChi2 = attempt ? 0 : ((param.rec.tpc.mergerInterpolateErrors && CAMath::Abs(ihit - ihitMergeFirst) <= 1) ? (refit ? (GPUTPCGMPropagator::rejectInterFill + ((nWays - iWay) & 1)) : 0) : (allowModification && goodRows > 5));
315#if EXTRACT_RESIDUALS == 1
316 if (iWay == nWays - 1 && interpolation.hit[ihit].errorY > (GPUCA_PAR_MERGER_INTERPOLATION_ERROR_TYPE_A)0) {
317 const float Iz0 = interpolation.hit[ihit].posY - mP[0];
318 const float Iz1 = interpolation.hit[ihit].posZ - mP[1];
319 float Iw0 = mC[2] + (float)interpolation.hit[ihit].errorZ;
320 float Iw2 = mC[0] + (float)interpolation.hit[ihit].errorY;
321 float Idet1 = 1.f / CAMath::Max(1e-10f, Iw0 * Iw2 - mC[1] * mC[1]);
322 const float Ik00 = (mC[0] * Iw0 + mC[1] * mC[1]) * Idet1;
323 const float Ik01 = (mC[0] * mC[1] + mC[1] * Iw2) * Idet1;
324 const float Ik10 = (mC[1] * Iw0 + mC[2] * mC[1]) * Idet1;
325 const float Ik11 = (mC[1] * mC[1] + mC[2] * Iw2) * Idet1;
326 const float ImP0 = mP[0] + Ik00 * Iz0 + Ik01 * Iz1;
327 const float ImP1 = mP[1] + Ik10 * Iz0 + Ik11 * Iz1;
328 const float ImC0 = mC[0] - Ik00 * mC[0] + Ik01 * mC[1];
329 const float ImC2 = mC[2] - Ik10 * mC[1] + Ik11 * mC[2];
330 auto& tup = GPUROOTDump<TNtuple>::get("clusterres", "row:clX:clY:clZ:angle:trkX:trkY:trkZ:trkSinPhi:trkDzDs:trkQPt:trkSigmaY2:trkSigmaZ2trkSigmaQPt2");
331 tup.Fill((float)cluster.row, xx, yy, zz, clAlpha, mX, ImP0, ImP1, mP[2], mP[3], mP[4], ImC0, ImC2, mC[14]);
332 }
333#endif
335 if (param.rec.tpc.rejectEdgeClustersInTrackFit && uncorrectedY > -1e6f && param.rejectEdgeClusterByY(uncorrectedY, cluster.row, CAMath::Sqrt(mC[0]))) { // uncorrectedY > -1e6f implies allowModification
337 } else {
338 const float time = merger->GetConstantMem()->ioPtrs.clustersNative ? merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear[cluster.num].getTime() : -1.f;
339 const float invSqrtCharge = merger->GetConstantMem()->ioPtrs.clustersNative ? CAMath::InvSqrt(merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear[cluster.num].qMax) : 0.f;
340 const float invCharge = merger->GetConstantMem()->ioPtrs.clustersNative ? (1.f / merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear[cluster.num].qMax) : 0.f;
341 float invAvgCharge = (sumInvSqrtCharge += invSqrtCharge) / ++nAvgCharge;
342 invAvgCharge *= invAvgCharge;
343 retVal = prop.Update(yy, zz, cluster.row, param, clusterState, rejectChi2, &interpolation.hit[ihit], refit, cluster.sector, time, invAvgCharge, invCharge GPUCA_DEBUG_STREAMER_CHECK(, &debugVals));
344 }
345 GPUCA_DEBUG_STREAMER_CHECK(if (o2::utils::DebugStreamer::checkStream(o2::utils::StreamFlags::streamUpdateTrack, iTrk)) {
346 merger->DebugStreamerUpdate(iTrk, ihit, xx, yy, zz, cluster, merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear[cluster.num], *this, prop, interpolation.hit[ihit], rejectChi2, refit, retVal, sumInvSqrtCharge / nAvgCharge * sumInvSqrtCharge / nAvgCharge, yy, zz, clusterState, debugVals.retVal, debugVals.err2Y, debugVals.err2Z);
347 });
348 }
349 // clang-format off
350 CADEBUG(if (!CheckCov()) GPUError("INVALID COV AFTER UPDATE!!!"));
351 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 - Err %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], retVal));
352 // clang-format on
353
354 ConstrainSinPhi();
355 if (retVal == 0) // track is updated
356 {
357 if (storeOuter == 1 && cluster.leg == clusters[maxN - 1].leg) {
358 StoreOuter(outerParam, prop, 2);
359 storeOuter = 255;
360 }
361 noFollowCircle2 = false;
362 lastUpdateX = mX;
363 covYYUpd = mC[0];
364 nMissed = nMissed2 = 0;
365 UnmarkClusters(clusters, ihitMergeFirst, ihit, wayDirection, GPUTPCGMMergedTrackHit::flagHighIncl);
366 N++;
367 ihitStart = ihit;
368 float dy = mP[0] - prop.Model().Y();
369 float dz = mP[1] - prop.Model().Z();
370 if (CAMath::Abs(mP[4]) * param.qptB5Scaler > 10 && --resetT0 <= 0 && CAMath::Abs(mP[2]) < 0.15f && dy * dy + dz * dz > 1) {
371 CADEBUG(printf("Reinit linearization\n"));
372 prop.SetTrack(this, prop.GetAlpha());
373 }
375 if (param.dodEdxEnabled && iWay == nWays - 1 && cluster.leg == clusters[maxN - 1].leg) { // TODO: Costimize flag to remove, and option to remove double-clusters
376 bool acc = (clusterState & param.rec.tpc.dEdxClusterRejectionFlagMask) == 0, accAlt = (clusterState & param.rec.tpc.dEdxClusterRejectionFlagMaskAlt) == 0;
377 if (acc || accAlt) {
378 float qtot = 0, qmax = 0, pad = 0, relTime = 0;
379 const int32_t clusterCount = (ihit - ihitMergeFirst) * wayDirection + 1;
380 for (int32_t iTmp = ihitMergeFirst; iTmp != ihit + wayDirection; iTmp += wayDirection) {
381 if (merger->GetConstantMem()->ioPtrs.clustersNative == nullptr) {
382 qtot += clustersXYZ[ihit].amp;
383 } else {
384 const ClusterNative& cl = merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear[cluster.num];
385 qtot += cl.qTot;
386 qmax = CAMath::Max<float>(qmax, cl.qMax);
387 pad += cl.getPad();
388 relTime += cl.getTime();
389 }
390 }
391 qtot /= clusterCount; // TODO: Weighted Average
392 pad /= clusterCount;
393 relTime /= clusterCount;
394 relTime = relTime - CAMath::Round(relTime);
395 if (acc) {
396 dEdx.fillCluster(qtot, qmax, cluster.row, cluster.sector, mP[2], mP[3], merger->GetConstantMem()->calibObjects, zz, pad, relTime);
397 }
398 if GPUCA_RTC_CONSTEXPR (GPUCA_GET_CONSTEXPR(param.rec.tpc, dEdxClusterRejectionFlagMask) != GPUCA_GET_CONSTEXPR(param.rec.tpc, dEdxClusterRejectionFlagMaskAlt)) {
399 if (accAlt) {
400 dEdxAlt.fillCluster(qtot, qmax, cluster.row, cluster.sector, mP[2], mP[3], merger->GetConstantMem()->calibObjects, zz, pad, relTime);
401 }
402 }
403 }
404 }
405 }
406 } else if (retVal >= GPUTPCGMPropagator::updateErrorClusterRejected) { // cluster far away form the track
407 if (allowModification) {
408 MarkClusters(clusters, ihitMergeFirst, ihit, wayDirection, GPUTPCGMMergedTrackHit::flagRejectDistance);
409 } else if (iWay == nWays - 1) {
410 MarkClusters(clusters, ihitMergeFirst, ihit, wayDirection, GPUTPCGMMergedTrackHit::flagRejectErr);
411 }
412 nMissed++;
413 nMissed2++;
414 } else {
415 break; // bad chi2 for the whole track, stop the fit
416 }
417 }
418 if (((nWays - iWay) & 1) && (iWay != nWays - 1) && (clusters[0].sector < 18) == (clusters[maxN - 1].sector < 18)) {
419 ShiftZ2(clusters, clustersXYZ, merger, maxN);
420 }
421 }
422 ConstrainSinPhi();
423
424 GPUCA_DEBUG_STREAMER_CHECK(if (o2::utils::DebugStreamer::checkStream(o2::utils::StreamFlags::streamUpdateTrack, iTrk)) {
425 o2::utils::DebugStreamer::instance()->getStreamer("debug_accept_track", "UPDATE") << o2::utils::DebugStreamer::instance()->getUniqueTreeName("debug_accept_track").data() << "iTrk=" << iTrk << "outerParam=" << *outerParam << "track=" << this << "ihitStart=" << ihitStart << "\n";
426 })
427
428 if (!(N + NTolerated >= GPUCA_TRACKLET_SELECTOR_MIN_HITS_B5(mP[4] * param.qptB5Scaler) && 2 * NTolerated <= CAMath::Max(10, N) && CheckNumericalQuality(covYYUpd))) {
429 return false; // TODO: NTolerated should never become that large, check what is going wrong!
430 }
431 if (param.rec.tpc.minNClustersFinalTrack != -1 && N + NTolerated < param.rec.tpc.minNClustersFinalTrack) {
432 return false;
433 }
434
435 // TODO: we have looping tracks here with 0 accepted clusters in the primary leg. In that case we should refit the track using only the primary leg.
436
437 if (param.par.dodEdx && param.dodEdxEnabled) {
438 dEdx.computedEdx(merger->MergedTracksdEdx()[iTrk], param);
439 if GPUCA_RTC_CONSTEXPR (GPUCA_GET_CONSTEXPR(param.rec.tpc, dEdxClusterRejectionFlagMask) != GPUCA_GET_CONSTEXPR(param.rec.tpc, dEdxClusterRejectionFlagMaskAlt)) {
440 dEdxAlt.computedEdx(merger->MergedTracksdEdxAlt()[iTrk], param);
441 }
442 }
443 Alpha = prop.GetAlpha();
444 MoveToReference(prop, param, Alpha);
445 NormalizeAlpha(Alpha);
446
447 return true;
448}
449
450GPUdni() void GPUTPCGMTrackParam::MoveToReference(GPUTPCGMPropagator& prop, const GPUParam& param, float& Alpha)
451{
452 static constexpr float kDeg2Rad = M_PI / 180.f;
453 static constexpr float kSectAngle = 2 * M_PI / 18.f;
454
455 if (param.rec.tpc.trackReferenceX <= 500) {
456 GPUTPCGMTrackParam save = *this;
457 float saveAlpha = Alpha;
458 for (int32_t attempt = 0; attempt < 3; attempt++) {
459 float dAngle = CAMath::Round(CAMath::ATan2(mP[0], mX) / kDeg2Rad / 20.f) * kSectAngle;
460 Alpha += dAngle;
461 if (prop.PropagateToXAlpha(param.rec.tpc.trackReferenceX, Alpha, 0)) {
462 break;
463 }
464 ConstrainSinPhi();
465 if (CAMath::Abs(mP[0]) <= mX * CAMath::Tan(kSectAngle / 2.f)) {
466 return;
467 }
468 }
469 *this = save;
470 Alpha = saveAlpha;
471 }
472 if (CAMath::Abs(mP[0]) > mX * CAMath::Tan(kSectAngle / 2.f)) {
473 float dAngle = CAMath::Round(CAMath::ATan2(mP[0], mX) / kDeg2Rad / 20.f) * kSectAngle;
474 Rotate(dAngle);
475 ConstrainSinPhi();
476 Alpha += dAngle;
477 }
478}
479
480GPUd() void GPUTPCGMTrackParam::MirrorTo(GPUTPCGMPropagator& GPUrestrict() prop, float toY, float toZ, bool inFlyDirection, const GPUParam& param, uint8_t row, uint8_t clusterState, bool mirrorParameters, int8_t sector)
481{
482 if (mirrorParameters) {
483 prop.Mirror(inFlyDirection);
484 }
485 float err2Y, err2Z;
486 prop.GetErr2(err2Y, err2Z, param, toZ, row, clusterState, sector, -1.f, 0.f, 0.f); // Use correct time / avgCharge
487 prop.Model().Y() = mP[0] = toY;
488 prop.Model().Z() = mP[1] = toZ;
489 if (mC[0] < err2Y) {
490 mC[0] = err2Y;
491 }
492 if (mC[2] < err2Z) {
493 mC[2] = err2Z;
494 }
495 if (CAMath::Abs(mC[5]) < 0.1f) {
496 mC[5] = mC[5] > 0 ? 0.1f : -0.1f;
497 }
498 if (mC[9] < 1.f) {
499 mC[9] = 1.f;
500 }
501 mC[1] = mC[4] = mC[6] = mC[8] = mC[11] = mC[13] = 0;
502 prop.SetTrack(this, prop.GetAlpha());
503 mNDF = -3;
504 mChi2 = 0;
505}
506
507GPUd() int32_t GPUTPCGMTrackParam::MergeDoubleRowClusters(int32_t& ihit, int32_t wayDirection, GPUTPCGMMergedTrackHit* GPUrestrict() clusters, GPUTPCGMMergedTrackHitXYZ* clustersXYZ, const GPUTPCGMMerger* GPUrestrict() merger, GPUTPCGMPropagator& GPUrestrict() prop, float& GPUrestrict() xx, float& GPUrestrict() yy, float& GPUrestrict() zz, int32_t maxN, float clAlpha, uint8_t& GPUrestrict() clusterState, bool rejectChi2)
508{
509 if (ihit + wayDirection >= 0 && ihit + wayDirection < maxN && clusters[ihit].row == clusters[ihit + wayDirection].row && clusters[ihit].sector == clusters[ihit + wayDirection].sector && clusters[ihit].leg == clusters[ihit + wayDirection].leg) {
510 float maxDistY, maxDistZ;
511 prop.GetErr2(maxDistY, maxDistZ, merger->Param(), zz, clusters[ihit].row, 0, clusters[ihit].sector, -1.f, 0.f, 0.f); // TODO: Use correct time, avgCharge
512 maxDistY = (maxDistY + mC[0]) * 20.f;
513 maxDistZ = (maxDistZ + mC[2]) * 20.f;
514 int32_t noReject = 0; // Cannot reject if simple estimation of y/z fails (extremely unlike case)
515 if (CAMath::Abs(clAlpha - prop.GetAlpha()) > 1.e-4f) {
516 noReject = prop.RotateToAlpha(clAlpha);
517 }
518 float projY = 0, projZ = 0;
519 if (noReject == 0) {
520 noReject |= prop.GetPropagatedYZ(xx, projY, projZ);
521 }
522 float count = 0.f;
523 xx = yy = zz = 0.f;
524 clusterState = 0;
525 while (true) {
526 float clx, cly, clz, clamp;
527 if (merger->Param().par.earlyTpcTransform) {
528 const float zOffset = (clusters[ihit].sector < 18) == (clusters[0].sector < 18) ? mTZOffset : -mTZOffset;
529 clx = clustersXYZ[ihit].x;
530 cly = clustersXYZ[ihit].y;
531 clz = clustersXYZ[ihit].z - zOffset;
532 clamp = clustersXYZ[ihit].amp;
533 } else {
534 const ClusterNative& GPUrestrict() cl = merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[ihit].num];
535 clamp = cl.qTot;
536 merger->GetConstantMem()->calibObjects.fastTransformHelper->Transform(clusters[ihit].sector, clusters[ihit].row, cl.getPad(), cl.getTime(), clx, cly, clz, mTZOffset);
537 }
538 float dy = cly - projY;
539 float dz = clz - projZ;
540 if (noReject == 0 && (dy * dy > maxDistY || dz * dz > maxDistZ)) {
541 CADEBUG(printf("\t\tRejecting double-row cluster: dy %f, dz %f, chiY %f, chiZ %f (Y: trk %f prj %f cl %f - Z: trk %f prj %f cl %f)\n", dy, dz, sqrtf(maxDistY), sqrtf(maxDistZ), mP[0], projY, cly, mP[1], projZ, clz));
542 if (rejectChi2) {
544 }
545 } else {
546 CADEBUG(printf("\t\tMerging hit row %d X %f Y %f Z %f (dy %f, dz %f, chiY %f, chiZ %f)\n", clusters[ihit].row, clx, cly, clz, dy, dz, sqrtf(maxDistY), sqrtf(maxDistZ)));
547 xx += clx * clamp; // TODO: Weight in pad/time instead of XYZ
548 yy += cly * clamp;
549 zz += clz * clamp;
550 clusterState |= clusters[ihit].state;
551 count += clamp;
552 }
553 if (!(ihit + wayDirection >= 0 && ihit + wayDirection < maxN && clusters[ihit].row == clusters[ihit + wayDirection].row && clusters[ihit].sector == clusters[ihit + wayDirection].sector && clusters[ihit].leg == clusters[ihit + wayDirection].leg)) {
554 break;
555 }
556 ihit += wayDirection;
557 }
558 if (count < 0.1f) {
559 CADEBUG(printf("\t\tNo matching cluster in double-row, skipping\n"));
560 return -1;
561 }
562 xx /= count;
563 yy /= count;
564 zz /= count;
565 }
566 return 0;
567}
568
569GPUd() float GPUTPCGMTrackParam::AttachClusters(const GPUTPCGMMerger* GPUrestrict() Merger, int32_t sector, int32_t iRow, int32_t iTrack, bool goodLeg, GPUTPCGMPropagator& prop)
570{
571 float Y, Z;
572 if (Merger->Param().par.earlyTpcTransform) {
573 Y = mP[0];
574 Z = mP[1];
575 } else {
576 float X = 0;
577 Merger->GetConstantMem()->calibObjects.fastTransformHelper->InverseTransformYZtoX(sector, iRow, mP[0], mP[1], X);
578 if (prop.GetPropagatedYZ(X, Y, Z)) {
579 Y = mP[0];
580 Z = mP[1];
581 }
582 }
583 return AttachClusters(Merger, sector, iRow, iTrack, goodLeg, Y, Z);
584}
585
586GPUd() float GPUTPCGMTrackParam::AttachClusters(const GPUTPCGMMerger* GPUrestrict() Merger, int32_t sector, int32_t iRow, int32_t iTrack, bool goodLeg, float Y, float Z)
587{
588 if (Merger->Param().rec.tpc.disableRefitAttachment & 1) {
589 return -1e6f;
590 }
591 const GPUTPCTracker& GPUrestrict() tracker = *(Merger -> GetConstantMem()->tpcTrackers + sector);
592 const GPUTPCRow& GPUrestrict() row = tracker.Row(iRow);
593 GPUglobalref() const cahit2* hits = tracker.HitData(row);
594 GPUglobalref() const calink* firsthit = tracker.FirstHitInBin(row);
595 if (row.NHits() == 0) {
596 return -1e6f;
597 }
598
599 const float zOffset = Merger->Param().par.earlyTpcTransform ? ((Merger->MergedTracks()[iTrack].CSide() ^ (sector >= 18)) ? -mTZOffset : mTZOffset) : Merger->GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->convVertexTimeToZOffset(sector, mTZOffset, Merger->Param().continuousMaxTimeBin);
600 const float y0 = row.Grid().YMin();
601 const float stepY = row.HstepY();
602 const float z0 = row.Grid().ZMin() - zOffset; // We can use our own ZOffset, since this is only used temporarily anyway
603 const float stepZ = row.HstepZ();
604 int32_t bin, ny, nz;
605
606 float err2Y, err2Z;
607 Merger->Param().GetClusterErrors2(sector, iRow, Z, mP[2], mP[3], -1.f, 0.f, 0.f, err2Y, err2Z); // TODO: Use correct time/avgCharge
608 const float sy2 = CAMath::Min(Merger->Param().rec.tpc.tubeMaxSize2, Merger->Param().rec.tpc.tubeChi2 * (err2Y + CAMath::Abs(mC[0]))); // Cov can be bogus when following circle
609 const float sz2 = CAMath::Min(Merger->Param().rec.tpc.tubeMaxSize2, Merger->Param().rec.tpc.tubeChi2 * (err2Z + CAMath::Abs(mC[2]))); // In that case we should provide the track error externally
610 const float tubeY = CAMath::Sqrt(sy2);
611 const float tubeZ = CAMath::Sqrt(sz2);
612 const float sy21 = 1.f / sy2;
613 const float sz21 = 1.f / sz2;
614 float uncorrectedY, uncorrectedZ;
615 if (Merger->Param().par.earlyTpcTransform) {
616 uncorrectedY = Y;
617 uncorrectedZ = Z;
618 } else {
619 Merger->GetConstantMem()->calibObjects.fastTransformHelper->InverseTransformYZtoNominalYZ(sector, iRow, Y, Z, uncorrectedY, uncorrectedZ);
620 }
621
622 if (CAMath::Abs(uncorrectedY) > row.getTPCMaxY()) {
623 return uncorrectedY;
624 }
625 row.Grid().GetBinArea(uncorrectedY, uncorrectedZ + zOffset, tubeY, tubeZ, bin, ny, nz);
626
627 const int32_t nBinsY = row.Grid().Ny();
628 const int32_t idOffset = tracker.Data().ClusterIdOffset();
629 const int32_t* ids = &(tracker.Data().ClusterDataIndex()[row.HitNumberOffset()]);
630 uint32_t myWeight = Merger->TrackOrderAttach()[iTrack] | gputpcgmmergertypes::attachAttached | gputpcgmmergertypes::attachTube;
631 GPUAtomic(uint32_t)* const weights = Merger->ClusterAttachment();
632 if (goodLeg) {
634 }
635 for (int32_t k = 0; k <= nz; k++) {
636 const int32_t mybin = bin + k * nBinsY;
637 const uint32_t hitFst = firsthit[mybin];
638 const uint32_t hitLst = firsthit[mybin + ny + 1];
639 for (uint32_t ih = hitFst; ih < hitLst; ih++) {
640 int32_t id = idOffset + ids[ih];
641 GPUAtomic(uint32_t)* const weight = weights + id;
642 if constexpr (GPUCA_PAR_NO_ATOMIC_PRECHECK == 0) {
643 if (myWeight <= *weight) {
644 continue;
645 }
646 }
647 const cahit2 hh = hits[ih];
648 const float y = y0 + hh.x * stepY;
649 const float z = z0 + hh.y * stepZ;
650 const float dy = y - uncorrectedY;
651 const float dz = z - uncorrectedZ;
652 if (dy * dy * sy21 + dz * dz * sz21 <= CAMath::Sqrt(2.f)) {
653 // CADEBUG(printf("Found Y %f Z %f\n", y, z));
654 CAMath::AtomicMax(weight, myWeight);
655 }
656 }
657 }
658 return uncorrectedY;
659}
660
661GPUd() bool GPUTPCGMTrackParam::AttachClustersPropagate(const GPUTPCGMMerger* GPUrestrict() Merger, int32_t sector, int32_t lastRow, int32_t toRow, int32_t iTrack, bool goodLeg, GPUTPCGMPropagator& GPUrestrict() prop, bool inFlyDirection, float maxSinPhi, bool dodEdx)
662{
663 static constexpr float kSectAngle = 2 * M_PI / 18.f;
664 if (Merger->Param().rec.tpc.disableRefitAttachment & 2) {
665 return dodEdx;
666 }
667 if (CAMath::Abs(lastRow - toRow) < 2) {
668 return dodEdx;
669 }
670 int32_t step = toRow > lastRow ? 1 : -1;
671 float xx = mX - GPUTPCGeometry::Row2X(lastRow);
672 for (int32_t iRow = lastRow + step; iRow != toRow; iRow += step) {
673 if (CAMath::Abs(mP[2]) > maxSinPhi) {
674 return dodEdx;
675 }
676 if (CAMath::Abs(mP[0]) > CAMath::Abs(mX) * CAMath::Tan(kSectAngle / 2.f)) {
677 return dodEdx;
678 }
679 int32_t err = prop.PropagateToXAlpha(xx + GPUTPCGeometry::Row2X(iRow), prop.GetAlpha(), inFlyDirection);
680 if (err) {
681 return dodEdx;
682 }
683 if (dodEdx && iRow + step == toRow) {
684 float yUncorrected, zUncorrected;
685 Merger->GetConstantMem()->calibObjects.fastTransformHelper->InverseTransformYZtoNominalYZ(sector, iRow, mP[0], mP[1], yUncorrected, zUncorrected);
686 uint32_t pad = CAMath::Float2UIntRn(GPUTPCGeometry::LinearY2Pad(sector, iRow, yUncorrected));
687 if (pad >= GPUTPCGeometry::NPads(iRow) || (Merger->GetConstantMem()->calibObjects.dEdxCalibContainer && Merger->GetConstantMem()->calibObjects.dEdxCalibContainer->isDead(sector, iRow, pad))) {
688 dodEdx = false;
689 }
690 }
691 CADEBUG(printf("Attaching in row %d\n", iRow));
692 AttachClusters(Merger, sector, iRow, iTrack, goodLeg, prop);
693 }
694 return dodEdx;
695}
696
697GPUd() bool GPUTPCGMTrackParam::FollowCircleChk(float lrFactor, float toY, float toX, bool up, bool right)
698{
699 return CAMath::Abs(mX * lrFactor - toY) > 1.f && // transport further in Y
700 CAMath::Abs(mP[2]) < 0.7f && // rotate back
701 (up ? (-mP[0] * lrFactor > toX || (right ^ (mP[2] > 0))) : (-mP[0] * lrFactor < toX || (right ^ (mP[2] < 0)))); // don't overshoot in X
702}
703
704GPUdii() void GPUTPCGMTrackParam::StoreOuter(gputpcgmmergertypes::GPUTPCOuterParam* outerParam, const GPUTPCGMPropagator& prop, int32_t phase)
705{
706 CADEBUG(printf("\t%21sStorO%d 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", "", phase, prop.GetAlpha(), mX, mP[0], mP[1], mP[4], prop.GetQPt0(), mP[2], prop.GetSinPhi0(), sqrtf(mC[0]), sqrtf(mC[2]), sqrtf(mC[5]), sqrtf(mC[14])));
707 for (int32_t i = 0; i < 5; i++) {
708 outerParam->P[i] = mP[i];
709 }
710 for (int32_t i = 0; i < 15; i++) {
711 outerParam->C[i] = mC[i];
712 }
713 outerParam->X = mX;
714 outerParam->alpha = prop.GetAlpha();
715}
716
717GPUdic(0, 1) void GPUTPCGMTrackParam::StoreAttachMirror(const GPUTPCGMMerger* GPUrestrict() Merger, int32_t sector, int32_t iRow, int32_t iTrack, float toAlpha, float toY, float toX, int32_t toSector, int32_t toRow, bool inFlyDirection, float alpha)
718{
719 uint32_t nLoopData = CAMath::AtomicAdd(&Merger->Memory()->nLoopData, 1u);
720 if (nLoopData >= Merger->NMaxTracks()) {
721 Merger->raiseError(GPUErrors::ERROR_MERGER_LOOPER_OVERFLOW, nLoopData, Merger->NMaxTracks());
722 CAMath::AtomicExch(&Merger->Memory()->nLoopData, Merger->NMaxTracks());
723 return;
724 }
726 data.param = *this;
727 data.alpha = alpha;
728 data.track = iTrack;
729 data.toAlpha = toAlpha;
730 data.toY = toY;
731 data.toX = toX;
732 data.sector = sector;
733 data.row = iRow;
734 data.toSector = toSector;
735 data.toRow = toRow;
736 data.inFlyDirection = inFlyDirection;
737 Merger->LoopData()[nLoopData] = data;
738}
739
740GPUdii() void GPUTPCGMTrackParam::RefitLoop(const GPUTPCGMMerger* GPUrestrict() Merger, int32_t loopIdx)
741{
743 prop.SetMaterialTPC();
744 prop.SetPolynomialField(&Merger->Param().polynomialField);
745 prop.SetMaxSinPhi(GPUCA_MAX_SIN_PHI);
746 prop.SetMatLUT(Merger->Param().rec.useMatLUT ? Merger->GetConstantMem()->calibObjects.matLUT : nullptr);
747 prop.SetSeedingErrors(false);
748 prop.SetFitInProjections(true);
749 prop.SetPropagateBzOnly(false);
750
751 GPUTPCGMLoopData& data = Merger->LoopData()[loopIdx];
752 prop.SetTrack(&data.param, data.alpha);
753 if (data.toSector == -1) {
754 data.param.AttachClustersMirror<1>(Merger, data.sector, data.row, data.track, data.toY, prop, true);
755 } else {
756 data.param.FollowCircle<1>(Merger, prop, data.sector, data.row, data.track, data.toAlpha, data.toX, data.toY, data.toSector, data.toRow, data.inFlyDirection, true);
757 }
758}
759
760template <int32_t I>
761GPUdic(0, 1) int32_t GPUTPCGMTrackParam::FollowCircle(const GPUTPCGMMerger* GPUrestrict() Merger, GPUTPCGMPropagator& GPUrestrict() prop, int32_t sector, int32_t iRow, int32_t iTrack, float toAlpha, float toX, float toY, int32_t toSector, int32_t toRow, bool inFlyDirection, bool phase2)
762{
763 static constexpr float kSectAngle = 2 * M_PI / 18.f;
764 if (Merger->Param().rec.tpc.disableRefitAttachment & 4) {
765 return 1;
766 }
767 const bool inExtraPass = Merger->Param().rec.tpc.looperInterpolationInExtraPass == -1 ? GPUCA_PAR_MERGER_SPLIT_LOOP_INTERPOLATION : Merger->Param().rec.tpc.looperInterpolationInExtraPass;
768 if (inExtraPass && phase2 == false) {
769 StoreAttachMirror(Merger, sector, iRow, iTrack, toAlpha, toY, toX, toSector, toRow, inFlyDirection, prop.GetAlpha());
770 return 1;
771 }
772 const GPUParam& GPUrestrict() param = Merger->Param();
773 bool right;
774 float dAlpha = toAlpha - prop.GetAlpha();
775 int32_t sectorSide = sector >= (GPUCA_NSECTORS / 2) ? (GPUCA_NSECTORS / 2) : 0;
776 if (CAMath::Abs(dAlpha) > 0.001f) {
777 right = CAMath::Abs(dAlpha) < CAMath::Pi() ? (dAlpha > 0) : (dAlpha < 0);
778 } else {
779 right = toY > mP[0];
780 }
781 bool up = (mP[2] < 0) ^ right;
782 int32_t targetRow = up ? (GPUCA_ROW_COUNT - 1) : 0;
783 float lrFactor = mP[2] < 0 ? -1.f : 1.f; // !(right ^ down) // TODO: shouldn't it be "right ? 1.f : -1.f", but that gives worse results...
784 // clang-format off
785 CADEBUG(printf("CIRCLE Track %d: Sector %d Alpha %f X %f Y %f Z %f SinPhi %f DzDs %f - Next hit: Sector %d Alpha %f X %f Y %f - Right %d Up %d dAlpha %f lrFactor %f\n", iTrack, sector, prop.GetAlpha(), mX, mP[0], mP[1], mP[2], mP[3], toSector, toAlpha, toX, toY, (int32_t)right, (int32_t)up, dAlpha, lrFactor));
786 // clang-format on
787
788 AttachClustersPropagate(Merger, sector, iRow, targetRow, iTrack, false, prop, inFlyDirection, 0.7f);
789 if (prop.RotateToAlpha(prop.GetAlpha() + (CAMath::Pi() / 2.f) * lrFactor)) {
790 return 1;
791 }
792 CADEBUG(printf("\tRotated: X %f Y %f Z %f SinPhi %f (Alpha %f / %f)\n", mP[0], mX, mP[1], mP[2], prop.GetAlpha(), prop.GetAlpha() + CAMath::Pi() / 2.f));
793 while (sector != toSector || FollowCircleChk(lrFactor, toY, toX, up, right)) {
794 while ((sector != toSector) ? (CAMath::Abs(mX) <= CAMath::Abs(mP[0]) * CAMath::Tan(kSectAngle / 2.f)) : FollowCircleChk(lrFactor, toY, toX, up, right)) {
795 int32_t err = prop.PropagateToXAlpha(mX + 1.f, prop.GetAlpha(), inFlyDirection);
796 if (err) {
797 CADEBUG(printf("\t\tpropagation error (%d)\n", err));
798 prop.RotateToAlpha(prop.GetAlpha() - (CAMath::Pi() / 2.f) * lrFactor);
799 return 1;
800 }
801 CADEBUG(printf("\tPropagated to y = %f: X %f Z %f SinPhi %f\n", mX, mP[0], mP[1], mP[2]));
802 for (int32_t j = 0; j < GPUCA_ROW_COUNT; j++) {
803 float rowX = GPUTPCGeometry::Row2X(j);
804 if (CAMath::Abs(rowX - (-mP[0] * lrFactor)) < 1.5f) {
805 CADEBUG(printf("\t\tAttempt row %d (Y %f Z %f)\n", j, mX * lrFactor, mP[1]));
806 AttachClusters(Merger, sector, j, iTrack, false, mX * lrFactor, mP[1]);
807 }
808 }
809 }
810 if (sector != toSector) {
811 if (right) {
812 if (++sector >= sectorSide + 18) {
813 sector -= 18;
814 }
815 } else {
816 if (--sector < sectorSide) {
817 sector += 18;
818 }
819 }
820 CADEBUG(printf("\tRotating to sector %d\n", sector));
821 if (prop.RotateToAlpha(param.Alpha(sector) + (CAMath::Pi() / 2.f) * lrFactor)) {
822 CADEBUG(printf("\t\trotation error\n"));
823 prop.RotateToAlpha(prop.GetAlpha() - (CAMath::Pi() / 2.f) * lrFactor);
824 return 1;
825 }
826 CADEBUG(printf("\tAfter Rotatin Alpha %f Position X %f Y %f Z %f SinPhi %f\n", prop.GetAlpha(), mP[0], mX, mP[1], mP[2]));
827 }
828 }
829 CADEBUG(printf("\tRotating back\n"));
830 for (int32_t i = 0; i < 2; i++) {
831 if (prop.RotateToAlpha(prop.GetAlpha() + (CAMath::Pi() / 2.f) * lrFactor) == 0) {
832 break;
833 }
834 if (i) {
835 CADEBUG(printf("Final rotation failed\n"));
836 return 1;
837 }
838 CADEBUG(printf("\tresetting physical model\n"));
839 prop.SetTrack(this, prop.GetAlpha());
840 }
841 prop.Rotate180();
842 CADEBUG(printf("\tMirrored position: Alpha %f X %f Y %f Z %f SinPhi %f DzDs %f\n", prop.GetAlpha(), mX, mP[0], mP[1], mP[2], mP[3]));
843 iRow = toRow;
844 float dx = toX - GPUTPCGeometry::Row2X(toRow);
845 if (up ^ (toX > mX)) {
846 if (up) {
847 while (iRow < GPUCA_ROW_COUNT - 2 && GPUTPCGeometry::Row2X(iRow + 1) + dx <= mX) {
848 iRow++;
849 }
850 } else {
851 while (iRow > 1 && GPUTPCGeometry::Row2X(iRow - 1) + dx >= mX) {
852 iRow--;
853 }
854 }
855 prop.PropagateToXAlpha(GPUTPCGeometry::Row2X(iRow) + dx, prop.GetAlpha(), inFlyDirection);
856 AttachClustersPropagate(Merger, sector, iRow, toRow, iTrack, false, prop, inFlyDirection);
857 }
858 if (prop.PropagateToXAlpha(toX, prop.GetAlpha(), inFlyDirection)) {
859 mX = toX;
860 }
861 CADEBUG(printf("Final position: Alpha %f X %f Y %f Z %f SinPhi %f DzDs %f\n", prop.GetAlpha(), mX, mP[0], mP[1], mP[2], mP[3]));
862 return (0);
863}
864
865template <int32_t I>
866GPUdni() void GPUTPCGMTrackParam::AttachClustersMirror(const GPUTPCGMMerger* GPUrestrict() Merger, int32_t sector, int32_t iRow, int32_t iTrack, float toY, GPUTPCGMPropagator& GPUrestrict() prop, bool phase2)
867{
868 static constexpr float kSectAngle = 2 * M_PI / 18.f;
869
870 if (Merger->Param().rec.tpc.disableRefitAttachment & 8) {
871 return;
872 }
873 const bool inExtraPass = Merger->Param().rec.tpc.looperInterpolationInExtraPass == -1 ? GPUCA_PAR_MERGER_SPLIT_LOOP_INTERPOLATION : Merger->Param().rec.tpc.looperInterpolationInExtraPass;
874 if (inExtraPass && phase2 == false) {
875 StoreAttachMirror(Merger, sector, iRow, iTrack, 0, toY, 0, -1, 0, 0, prop.GetAlpha());
876 return;
877 }
878 // Note that the coordinate system is rotated by 90 degree swapping X and Y!
879 float X = mP[2] > 0 ? mP[0] : -mP[0];
880 float toX = mP[2] > 0 ? toY : -toY;
881 float Y = mP[2] > 0 ? -mX : mX;
882 float Z = mP[1];
883 if (CAMath::Abs(mP[2]) >= GPUCA_MAX_SIN_PHI_LOW) {
884 return;
885 }
886 float SinPhi = CAMath::Sqrt(1 - mP[2] * mP[2]) * (mP[2] > 0 ? -1 : 1);
887 if (CAMath::Abs(SinPhi) >= GPUCA_MAX_SIN_PHI_LOW) {
888 return;
889 }
890 float b = prop.GetBz(prop.GetAlpha(), mX, mP[0], mP[1]);
891
892 int32_t count = CAMath::Float2IntRn(CAMath::Abs((toX - X) * 2.f));
893 if (count == 0) {
894 return;
895 }
896 float dx = (toX - X) / count;
897 const float myRowX = GPUTPCGeometry::Row2X(iRow);
898 // printf("AttachMirror\n");
899 // printf("X %f Y %f Z %f SinPhi %f toY %f -->\n", mX, mP[0], mP[1], mP[2], toY);
900 // printf("X %f Y %f Z %f SinPhi %f, count %d dx %f (to: %f)\n", X, Y, Z, SinPhi, count, dx, X + count * dx);
901 while (count--) {
902 float ex = CAMath::Sqrt(1 - SinPhi * SinPhi);
903 float exi = 1.f / ex;
904 float dxBzQ = dx * -b * mP[4];
905 float newSinPhi = SinPhi + dxBzQ;
906 if (CAMath::Abs(newSinPhi) > GPUCA_MAX_SIN_PHI_LOW) {
907 return;
908 }
909 float dS = dx * exi;
910 float h2 = dS * exi * exi;
911 float h4 = .5f * h2 * dxBzQ;
912
913 X += dx;
914 Y += dS * SinPhi + h4;
915 Z += dS * mP[3];
916 SinPhi = newSinPhi;
917 if (CAMath::Abs(X) > CAMath::Abs(Y) * CAMath::Tan(kSectAngle / 2.f)) {
918 continue;
919 }
920
921 // printf("count %d: At X %f Y %f Z %f SinPhi %f\n", count, mP[2] > 0 ? -Y : Y, mP[2] > 0 ? X : -X, Z, SinPhi);
922
923 float paramX = mP[2] > 0 ? -Y : Y;
924 int32_t step = paramX >= mX ? 1 : -1;
925 int32_t found = 0;
926 for (int32_t j = iRow; j >= 0 && j < GPUCA_ROW_COUNT && found < 3; j += step) {
927 float rowX = mX + GPUTPCGeometry::Row2X(j) - myRowX;
928 if (CAMath::Abs(rowX - paramX) < 1.5f) {
929 // printf("Attempt row %d\n", j);
930 AttachClusters(Merger, sector, j, iTrack, false, mP[2] > 0 ? X : -X, Z);
931 }
932 }
933 }
934}
935
936GPUd() void GPUTPCGMTrackParam::ShiftZ2(const GPUTPCGMMergedTrackHit* clusters, GPUTPCGMMergedTrackHitXYZ* clustersXYZ, const GPUTPCGMMerger* merger, int32_t N)
937{
938 float tzInner, tzOuter;
939 float xInner, xOuter;
940 if (N == 0) {
941 N = 1;
942 }
943 if (merger->Param().par.earlyTpcTransform) {
944 tzInner = clustersXYZ[N - 1].z;
945 tzOuter = clustersXYZ[0].z;
946 xInner = clustersXYZ[N - 1].x;
947 xOuter = clustersXYZ[0].x;
948 } else {
949 const auto& GPUrestrict() cls = merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear;
950 tzInner = cls[clusters[N - 1].num].getTime();
951 tzOuter = cls[clusters[0].num].getTime();
952 xInner = GPUTPCGeometry::Row2X(clusters[N - 1].row);
953 xOuter = GPUTPCGeometry::Row2X(clusters[0].row);
954 }
955 ShiftZ(merger, clusters[0].sector, tzInner, tzOuter, xInner, xOuter);
956}
957
958GPUd() void GPUTPCGMTrackParam::ShiftZ(const GPUTPCGMMerger* GPUrestrict() merger, int32_t sector, float tz1, float tz2, float x1, float x2)
959{
960 if (!merger->Param().par.continuousTracking) {
961 return;
962 }
963 const float r1 = CAMath::Max(0.0001f, CAMath::Abs(mP[4] * merger->Param().polynomialField.GetNominalBz()));
964
965 const float dist2 = mX * mX + mP[0] * mP[0];
966 const float dist1r2 = dist2 * r1 * r1;
967 float deltaZ = 0.f;
968 bool beamlineReached = false;
969 if (dist1r2 < 4) {
970 const float alpha = CAMath::ACos(1 - 0.5f * dist1r2); // Angle of a circle, such that |(cosa, sina) - (1,0)| == dist
971 const float beta = CAMath::ATan2(mP[0], mX);
972 const int32_t comp = mP[2] > CAMath::Sin(beta);
973 const float sinab = CAMath::Sin((comp ? 0.5f : -0.5f) * alpha + beta); // Angle of circle through origin and track position, to be compared to Snp
974 const float res = CAMath::Abs(sinab - mP[2]);
975
976 if (res < 0.2) {
977 const float r = 1.f / r1;
978 const float dS = alpha * r;
979 float z0 = dS * mP[3];
980 if (CAMath::Abs(z0) > GPUTPCGeometry::TPCLength()) {
981 z0 = z0 > 0 ? GPUTPCGeometry::TPCLength() : -GPUTPCGeometry::TPCLength();
982 }
983 deltaZ = mP[1] - z0;
984 beamlineReached = true;
985
986 // printf("X %9.3f Y %9.3f QPt %9.3f R %9.3f --> Alpha %9.3f Snp %9.3f Snab %9.3f Res %9.3f dS %9.3f z0 %9.3f\n", mX, mP[0], mP[4], r, alpha / 3.1415 * 180, mP[2], sinab, res, dS, z0);
987 }
988 }
989
990 if (!beamlineReached) {
991 if (merger->Param().par.earlyTpcTransform) {
992 float basez, basex;
993 if (CAMath::Abs(tz1) < CAMath::Abs(tz2)) {
994 basez = tz1;
995 basex = x1;
996 } else {
997 basez = tz2;
998 basex = x2;
999 }
1000 float refZ = ((basez > 0) ? merger->Param().rec.tpc.defaultZOffsetOverR : -merger->Param().rec.tpc.defaultZOffsetOverR) * basex;
1001 deltaZ = basez - refZ - mTZOffset;
1002 } else {
1003 float baset, basex;
1004 if (CAMath::Abs(tz1) > CAMath::Abs(tz2)) {
1005 baset = tz1;
1006 basex = x1;
1007 } else {
1008 baset = tz2;
1009 basex = x2;
1010 }
1011 float refZ = ((sector < GPUCA_NSECTORS / 2) ? merger->Param().rec.tpc.defaultZOffsetOverR : -merger->Param().rec.tpc.defaultZOffsetOverR) * basex;
1012 float basez;
1013 merger->GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->TransformIdealZ(sector, baset, basez, mTZOffset);
1014 deltaZ = basez - refZ;
1015 }
1016 }
1017 if (merger->Param().par.earlyTpcTransform) {
1018 mTZOffset += deltaZ;
1019 mP[1] -= deltaZ;
1020 deltaZ = 0;
1021 float zMax = CAMath::Max(tz1, tz2);
1022 float zMin = CAMath::Min(tz1, tz2);
1023 // printf("Z Check: Clusters %f %f, min %f max %f vtx %f\n", tz1, tz2, zMin, zMax, mTZOffset);
1024 if (zMin < 0 && zMin - mTZOffset < -GPUTPCGeometry::TPCLength()) {
1025 deltaZ = zMin - mTZOffset + GPUTPCGeometry::TPCLength();
1026 } else if (zMax > 0 && zMax - mTZOffset > GPUTPCGeometry::TPCLength()) {
1027 deltaZ = zMax - mTZOffset - GPUTPCGeometry::TPCLength();
1028 }
1029 if (zMin < 0 && zMax - (mTZOffset + deltaZ) > 0) {
1030 deltaZ = zMax - mTZOffset;
1031 } else if (zMax > 0 && zMin - (mTZOffset + deltaZ) < 0) {
1032 deltaZ = zMin - mTZOffset;
1033 }
1034 // if (deltaZ != 0) printf("Moving clusters to TPC Range: Shift %f in Z: %f to %f --> %f to %f in Z\n", deltaZ, tz2 - mTZOffset, tz1 - mTZOffset, tz2 - mTZOffset - deltaZ, tz1 - mTZOffset - deltaZ);
1035 mTZOffset += deltaZ;
1036 mP[1] -= deltaZ;
1037 } else {
1038 float deltaT = merger->GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->convDeltaZtoDeltaTimeInTimeFrame(sector, deltaZ);
1039 mTZOffset += deltaT;
1040 mP[1] -= deltaZ;
1041 const float maxT = CAMath::Min(tz1, tz2) - merger->GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->getT0();
1042 const float minT = CAMath::Max(tz1, tz2) - merger->GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->getMaxDriftTime(sector);
1043 // printf("T Check: Clusters %f %f, min %f max %f vtx %f\n", tz1, tz2, minT, maxT, mTZOffset);
1044 deltaT = 0.f;
1045 if (mTZOffset < minT) {
1046 deltaT = minT - mTZOffset;
1047 }
1048 if (mTZOffset + deltaT > maxT) {
1049 deltaT = maxT - mTZOffset;
1050 }
1051 if (deltaT != 0.f) {
1052 deltaZ = merger->GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->convDeltaTimeToDeltaZinTimeFrame(sector, deltaT);
1053 // printf("Moving clusters to TPC Range: QPt %f, New mTZOffset %f, t1 %f, t2 %f, Shift %f in Z: %f to %f --> %f to %f in T\n", mP[4], mTZOffset + deltaT, tz1, tz2, deltaZ, tz2 - mTZOffset, tz1 - mTZOffset, tz2 - mTZOffset - deltaT, tz1 - mTZOffset - deltaT);
1054 mTZOffset += deltaT;
1055 mP[1] -= deltaZ;
1056 }
1057 }
1058 // printf("\n");
1059}
1060
1061GPUd() bool GPUTPCGMTrackParam::CheckCov() const
1062{
1063 const float* c = mC;
1064 bool ok = c[0] >= 0 && c[2] >= 0 && c[5] >= 0 && c[9] >= 0 && c[14] >= 0 && (c[1] * c[1] <= c[2] * c[0]) && (c[3] * c[3] <= c[5] * c[0]) && (c[4] * c[4] <= c[5] * c[2]) && (c[6] * c[6] <= c[9] * c[0]) && (c[7] * c[7] <= c[9] * c[2]) && (c[8] * c[8] <= c[9] * c[5]) &&
1065 (c[10] * c[10] <= c[14] * c[0]) && (c[11] * c[11] <= c[14] * c[2]) && (c[12] * c[12] <= c[14] * c[5]) && (c[13] * c[13] <= c[14] * c[9]);
1066 return ok;
1067}
1068
1069GPUd() bool GPUTPCGMTrackParam::CheckNumericalQuality(float overrideCovYY) const
1070{
1071 //* Check that the track parameters and covariance matrix are reasonable
1072 bool ok = CAMath::Finite(mX) && CAMath::Finite(mChi2);
1073 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"));
1074 const float* c = mC;
1075 for (int32_t i = 0; i < 15; i++) {
1076 ok = ok && CAMath::Finite(c[i]);
1077 }
1078 CADEBUG(printf("OK1 %d\n", (int32_t)ok));
1079 for (int32_t i = 0; i < 5; i++) {
1080 ok = ok && CAMath::Finite(mP[i]);
1081 }
1082 CADEBUG(printf("OK2 %d\n", (int32_t)ok));
1083 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) {
1084 ok = 0;
1085 }
1086 CADEBUG(printf("OK3 %d\n", (int32_t)ok));
1087 if (CAMath::Abs(mP[2]) > GPUCA_MAX_SIN_PHI) {
1088 ok = 0;
1089 }
1090 CADEBUG(printf("OK4 %d\n", (int32_t)ok));
1091 if (!CheckCov()) {
1092 ok = false;
1093 }
1094 CADEBUG(printf("OK5 %d\n", (int32_t)ok));
1095 return ok;
1096}
1097
1098GPUd() void GPUTPCGMTrackParam::RefitTrack(GPUTPCGMMergedTrack& GPUrestrict() track, int32_t iTrk, GPUTPCGMMerger* GPUrestrict() merger, int32_t attempt) // TODO: Inline me, once __forceinline__ is fixed by HIP
1099{
1100 if (!track.OK()) {
1101 return;
1102 }
1103
1104 // clang-format off
1105 CADEBUG(if (DEBUG_SINGLE_TRACK >= 0 && iTrk != DEBUG_SINGLE_TRACK) { track.SetNClusters(0); track.SetOK(0); return; } );
1106 // clang-format on
1107
1108 int32_t nTrackHits = track.NClusters();
1109 int32_t NTolerated = 0; // Clusters not fit but tollerated for track length cut
1110 GPUTPCGMTrackParam t = track.Param();
1111 float Alpha = track.Alpha();
1112 CADEBUG(int32_t nTrackHitsOld = nTrackHits; float ptOld = t.QPt());
1113 bool ok = t.Fit(merger, iTrk, merger->Clusters() + track.FirstClusterRef(), merger->Param().par.earlyTpcTransform ? merger->ClustersXYZ() + track.FirstClusterRef() : nullptr, nTrackHits, NTolerated, Alpha, attempt, GPUCA_MAX_SIN_PHI, &track.OuterParam());
1114 CADEBUG(printf("Finished Fit Track %d\n", iTrk));
1115 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)));
1116
1117 if (!ok && attempt == 0 && merger->Param().rec.tpc.retryRefit) {
1118 for (uint32_t i = 0; i < track.NClusters(); i++) {
1119 merger->Clusters()[track.FirstClusterRef() + i].state &= GPUTPCGMMergedTrackHit::clustererAndSharedFlags;
1120 }
1121 CADEBUG(printf("Track rejected, marking for retry\n"));
1122 if (merger->Param().rec.tpc.retryRefit == 2) {
1123 nTrackHits = track.NClusters();
1124 NTolerated = 0; // Clusters not fit but tollerated for track length cut
1125 t = track.Param();
1126 Alpha = track.Alpha();
1127 ok = t.Fit(merger, iTrk, merger->Clusters() + track.FirstClusterRef(), merger->ClustersXYZ() + track.FirstClusterRef(), nTrackHits, NTolerated, Alpha, 1, GPUCA_MAX_SIN_PHI, &track.OuterParam());
1128 } else {
1129 uint32_t nRefit = CAMath::AtomicAdd(&merger->Memory()->nRetryRefit, 1u);
1130 merger->RetryRefitIds()[nRefit] = iTrk;
1131 return;
1132 }
1133 }
1134 if (CAMath::Abs(t.QPt()) < 1.e-4f) {
1135 t.QPt() = 1.e-4f;
1136 }
1137
1138 CADEBUG(if (t.GetX() > 250) { printf("ERROR, Track %d at impossible X %f, Pt %f, Looper %d\n", iTrk, t.GetX(), CAMath::Abs(1.f / t.QPt()), (int32_t)merger->MergedTracks()[iTrk].Looper()); });
1139
1140 track.SetOK(ok);
1141 track.SetNClustersFitted(nTrackHits);
1142 track.Param() = t;
1143 track.Alpha() = Alpha;
1144
1145 // if (track.OK()) merger->DebugRefitMergedTrack(track);
1146}
1147
1148GPUd() void GPUTPCGMTrackParam::Rotate(float alpha)
1149{
1150 float cA, sA;
1151 CAMath::SinCos(alpha, sA, cA);
1152 float x0 = mX;
1153 float sinPhi0 = mP[2], cosPhi0 = CAMath::Sqrt(1 - mP[2] * mP[2]);
1154 float cosPhi = cosPhi0 * cA + sinPhi0 * sA;
1155 float sinPhi = -cosPhi0 * sA + sinPhi0 * cA;
1156 float j0 = cosPhi0 / cosPhi;
1157 float j2 = cosPhi / cosPhi0;
1158 mX = x0 * cA + mP[0] * sA;
1159 mP[0] = -x0 * sA + mP[0] * cA;
1160 mP[2] = sinPhi;
1161 mC[0] *= j0 * j0;
1162 mC[1] *= j0;
1163 mC[3] *= j0;
1164 mC[6] *= j0;
1165 mC[10] *= j0;
1166
1167 mC[3] *= j2;
1168 mC[4] *= j2;
1169 mC[5] *= j2 * j2;
1170 mC[8] *= j2;
1171 mC[12] *= j2;
1172 if (cosPhi < 0) { // change direction ( t0 direction is already changed in t0.UpdateValues(); )
1173 SinPhi() = -SinPhi();
1174 DzDs() = -DzDs();
1175 QPt() = -QPt();
1176 mC[3] = -mC[3];
1177 mC[4] = -mC[4];
1178 mC[6] = -mC[6];
1179 mC[7] = -mC[7];
1180 mC[10] = -mC[10];
1181 mC[11] = -mC[11];
1182 }
1183}
1184
1185GPUd() void GPUTPCGMTrackParam::AddCovDiagErrors(const float* GPUrestrict() errors2)
1186{
1187 mC[0] += errors2[0];
1188 mC[2] += errors2[1];
1189 mC[5] += errors2[2];
1190 mC[9] += errors2[3];
1191 mC[14] += errors2[4];
1192}
1193
1194GPUd() void GPUTPCGMTrackParam::AddCovDiagErrorsWithCorrelations(const float* GPUrestrict() errors2)
1195{
1196 const int32_t diagMap[5] = {0, 2, 5, 9, 14};
1197 const float oldDiag[5] = {mC[0], mC[2], mC[5], mC[9], mC[14]};
1198 for (int32_t i = 0; i < 5; i++) {
1199 mC[diagMap[i]] += errors2[i];
1200 for (int32_t j = 0; j < i; j++) {
1201 mC[diagMap[i - 1] + j + 1] *= gpu::CAMath::Sqrt(mC[diagMap[i]] * mC[diagMap[j]] / (oldDiag[i] * oldDiag[j]));
1202 }
1203 }
1204}
benchmark::State & state
Helper class to access correction maps.
uint64_t phase
Definition RawEventData.h:7
int16_t time
Definition RawEventData.h:4
int32_t i
#define GPUdic(...)
#define GPUAtomic(type)
#define GPUrestrict()
#define GPUglobalref()
#define GPUCA_DEBUG_STREAMER_CHECK(...)
#define GPUCA_RTC_CONSTEXPR
#define GPUCA_MAX_SIN_PHI_LOW
#define GPUCA_TRACKLET_SELECTOR_MIN_HITS_B5(QPTB5)
#define GPUCA_MAX_SIN_PHI
#define GPUCA_PAR_NO_ATOMIC_PRECHECK
#define GPUCA_PAR_MERGER_SPLIT_LOOP_INTERPOLATION
#define GPUCA_PAR_MERGER_INTERPOLATION_ERROR_TYPE_A
#define CADEBUG(...)
Definition GPUDef.h:79
#define GPUCA_GET_CONSTEXPR(obj, val)
int32_t retVal
uint8_t leg
GPUdii() void GPUTPCGMTrackParam
GPUdni() void GPUTPCGMTrackParam
#define DEBUG_SINGLE_TRACK
#define GPUCA_NSECTORS
#define GPUCA_ROW_COUNT
float float float & zMax
float float & zMin
uint32_t j
Definition RawData.h:0
uint32_t res
Definition RawData.h:0
uint32_t c
Definition RawData.h:2
Definition of TPCFastTransform class.
double num
static GPUROOTDump< T, Args... > & get(const char *name1, Names... names)
Definition GPUROOTDump.h:58
float const GPUParam float int32_t int16_t int8_t sector
float int32_t const GPUParam & param
int32_t GPUTPCGMMergedTrackHit GPUTPCGMMergedTrackHitXYZ int32_t int32_t float & Alpha
int32_t int32_t int32_t bool float Y
const GPUParam float & alpha
GLfloat GLfloat GLfloat alpha
Definition glcorearb.h:279
GLsizei const GLuint const GLfloat * weights
Definition glcorearb.h:5475
GLint GLsizei count
Definition glcorearb.h:399
GLuint * ids
Definition glcorearb.h:647
GLuint GLfloat GLfloat GLfloat x1
Definition glcorearb.h:5034
GLdouble GLdouble right
Definition glcorearb.h:4077
GLuint GLuint GLfloat weight
Definition glcorearb.h:5477
GLboolean GLboolean GLboolean b
Definition glcorearb.h:1233
GLboolean * data
Definition glcorearb.h:298
typedef void(APIENTRYP PFNGLCULLFACEPROC)(GLenum mode)
GLuint GLfloat x0
Definition glcorearb.h:5034
GLenum clamp
Definition glcorearb.h:1245
GLboolean r
Definition glcorearb.h:1233
GLenum GLfloat param
Definition glcorearb.h:271
GLuint id
Definition glcorearb.h:650
GLuint GLfloat GLfloat y0
Definition glcorearb.h:5034
GLdouble GLdouble GLdouble z
Definition glcorearb.h:843
uint32_t calink
Definition GPUTPCDef.h:30
const float3 float float x2
Definition MathUtils.h:42
Global TPC definitions and constants.
Definition SimTraits.h:167
GPUd() void PIDResponse
Definition PIDResponse.h:71
value_T step
Definition TrackUtils.h:42
value_T maxT
Definition TrackUtils.h:132
@ streamUpdateTrack
stream update track informations
std::string getTime(uint64_t ts)
GPUCA_PAR_MERGER_INTERPOLATION_ERROR_TYPE_A errorZ
GPUCA_PAR_MERGER_INTERPOLATION_ERROR_TYPE_A errorY
InterpolationErrorHit hit[GPUCA_MERGER_MAX_TRACK_CLUSTERS]
std::vector< Cluster > clusters
std::vector< int > row