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
43#ifdef GPUCA_CADEBUG_ENABLED
44#include "../utils/qconfig.h"
46#endif
47
48#ifndef GPUCA_GPUCODE_DEVICE
49#include <cmath>
50#include <cstdlib>
51#endif
52
53using namespace o2::gpu;
54using namespace o2::tpc;
55
56GPUd() 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)
57{
58 static constexpr float kDeg2Rad = M_PI / 180.f;
59 CADEBUG(static constexpr float kSectAngle = 2 * M_PI / 18.f);
60
61 const GPUParam& GPUrestrict() param = merger->Param();
62
63 GPUdEdx dEdx, dEdxAlt;
66 prop.SetMaterialTPC();
67 prop.SetPolynomialField(&param.polynomialField);
68 prop.SetMaxSinPhi(maxSinPhi);
69 prop.SetToyMCEventsFlag(param.par.toyMCEventsFlag);
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->OutputTracks()[iTrk].Looper()) {
98 storeOuter = 1;
99 }
100 } else if (iWay == nWays - 2 && merger->OutputTracks()[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) {
219 bool dodEdx = param.par.dodEdx && param.dodEdxDownscaled && param.rec.tpc.adddEdxSubThresholdClusters && iWay == nWays - 1 && CAMath::Abs(cluster.row - lastRow) == 2 && cluster.leg == clusters[maxN - 1].leg;
220 dodEdx = AttachClustersPropagate(merger, cluster.sector, lastRow, cluster.row, iTrk, cluster.leg == clusters[maxN - 1].leg, prop, inFlyDirection, GPUCA_MAX_SIN_PHI, dodEdx);
221 if (dodEdx) {
222 dEdx.fillSubThreshold(lastRow - wayDirection);
223 dEdxAlt.fillSubThreshold(lastRow - wayDirection);
224 }
225 }
226
227 int32_t err = prop.PropagateToXAlpha(xx, clAlpha, inFlyDirection);
228 // clang-format off
229 CADEBUG(if (!CheckCov()){printf("INVALID COV AFTER PROPAGATE!!!\n");});
230 // clang-format on
231 if (err == -2) // Rotation failed, try to bring to new x with old alpha first, rotate, and then propagate to x, alpha
232 {
233 CADEBUG(printf("REROTATE\n"));
234 if (prop.PropagateToXAlpha(xx, prop.GetAlpha(), inFlyDirection) == 0) {
235 err = prop.PropagateToXAlpha(xx, clAlpha, inFlyDirection);
236 }
237 }
238 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)) {
239 goodRows = 0;
240 } else {
241 goodRows++;
242 }
243 if (err == 0) {
244 lastRow = cluster.row;
245 lastSector = cluster.sector;
246 }
247 // clang-format off
248 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));
249 // clang-format on
250
251 if (crossCE) {
252 if (param.rec.tpc.addErrorsCECrossing) {
253 if (param.rec.tpc.addErrorsCECrossing >= 2) {
254 AddCovDiagErrorsWithCorrelations(param.rec.tpc.errorsCECrossing);
255 } else {
256 AddCovDiagErrors(param.rec.tpc.errorsCECrossing);
257 }
258 } else if (mC[2] < 0.5f) {
259 mC[2] = 0.5f;
260 }
261 }
262
263 if (err == 0 && changeDirection) {
264 const float mirrordY = prop.GetMirroredYTrack();
265 CADEBUG(printf(" -- MirroredY: %f --> %f", mP[0], mirrordY));
266 if (CAMath::Abs(yy - mP[0]) > CAMath::Abs(yy - mirrordY)) {
267 CADEBUG(printf(" - Mirroring!!!"));
268 if (allowModification) {
269 AttachClustersMirror<0>(merger, cluster.sector, cluster.row, iTrk, yy, prop); // TODO: Never true, will always call FollowCircle above, really???
270 }
271 MirrorTo(prop, yy, zz, inFlyDirection, param, cluster.row, clusterState, true, cluster.sector);
272 noFollowCircle = false;
273
274 lastUpdateX = mX;
275 lastLeg = cluster.leg;
276 lastRow = 255;
277 N++;
278 resetT0 = initResetT0();
279 // clang-format off
280 CADEBUG(printf("\n"));
281 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]));
282 // clang-format on
283 continue;
284 }
285 }
286
287 float uncorrectedY = -1e6f;
288 if (allowModification) {
289 uncorrectedY = AttachClusters(merger, cluster.sector, cluster.row, iTrk, cluster.leg == clusters[maxN - 1].leg, prop);
290 }
291
292 const int32_t err2 = mNDF > 0 && CAMath::Abs(prop.GetSinPhi0()) >= maxSinForUpdate;
293 if (err || err2) {
294 if (mC[0] > param.rec.tpc.trackFitCovLimit || mC[2] > param.rec.tpc.trackFitCovLimit) {
295 break;
296 }
297 MarkClusters(clusters, ihitMergeFirst, ihit, wayDirection, GPUTPCGMMergedTrackHit::flagNotFit);
298 nMissed2++;
299 NTolerated++;
300 CADEBUG(printf(" --- break (%d, %d)\n", err, err2));
301 continue;
302 }
303 CADEBUG(printf("\n"));
304
305 int32_t retVal;
306 float threshold = 3.f + (lastUpdateX >= 0 ? (CAMath::Abs(mX - lastUpdateX) / 2) : 0.f);
307 if (mNDF > 5 && (CAMath::Abs(yy - mP[0]) > threshold || CAMath::Abs(zz - mP[1]) > threshold)) {
309 } else {
310 int8_t rejectChi2 = attempt ? 0 : ((param.rec.tpc.mergerInterpolateErrors && CAMath::Abs(ihit - ihitMergeFirst) <= 1) ? (refit ? (GPUTPCGMPropagator::rejectInterFill + ((nWays - iWay) & 1)) : 0) : (allowModification && goodRows > 5));
311#if EXTRACT_RESIDUALS == 1
312 if (iWay == nWays - 1 && interpolation.hit[ihit].errorY > (GPUCA_MERGER_INTERPOLATION_ERROR_TYPE_A)0) {
313 const float Iz0 = interpolation.hit[ihit].posY - mP[0];
314 const float Iz1 = interpolation.hit[ihit].posZ - mP[1];
315 float Iw0 = mC[2] + (float)interpolation.hit[ihit].errorZ;
316 float Iw2 = mC[0] + (float)interpolation.hit[ihit].errorY;
317 float Idet1 = 1.f / CAMath::Max(1e-10f, Iw0 * Iw2 - mC[1] * mC[1]);
318 const float Ik00 = (mC[0] * Iw0 + mC[1] * mC[1]) * Idet1;
319 const float Ik01 = (mC[0] * mC[1] + mC[1] * Iw2) * Idet1;
320 const float Ik10 = (mC[1] * Iw0 + mC[2] * mC[1]) * Idet1;
321 const float Ik11 = (mC[1] * mC[1] + mC[2] * Iw2) * Idet1;
322 const float ImP0 = mP[0] + Ik00 * Iz0 + Ik01 * Iz1;
323 const float ImP1 = mP[1] + Ik10 * Iz0 + Ik11 * Iz1;
324 const float ImC0 = mC[0] - Ik00 * mC[0] + Ik01 * mC[1];
325 const float ImC2 = mC[2] - Ik10 * mC[1] + Ik11 * mC[2];
326 auto& tup = GPUROOTDump<TNtuple>::get("clusterres", "row:clX:clY:clZ:angle:trkX:trkY:trkZ:trkSinPhi:trkDzDs:trkQPt:trkSigmaY2:trkSigmaZ2trkSigmaQPt2");
327 tup.Fill((float)cluster.row, xx, yy, zz, clAlpha, mX, ImP0, ImP1, mP[2], mP[3], mP[4], ImC0, ImC2, mC[14]);
328 }
329#endif
331 if (param.rec.tpc.rejectEdgeClustersInTrackFit && uncorrectedY > -1e6f && param.rejectEdgeClusterByY(uncorrectedY, cluster.row, CAMath::Sqrt(mC[0]))) { // uncorrectedY > -1e6f implies allowModification
333 } else {
334 const float time = merger->GetConstantMem()->ioPtrs.clustersNative ? merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear[cluster.num].getTime() : -1.f;
335 const float invSqrtCharge = merger->GetConstantMem()->ioPtrs.clustersNative ? CAMath::InvSqrt(merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear[cluster.num].qMax) : 0.f;
336 const float invCharge = merger->GetConstantMem()->ioPtrs.clustersNative ? (1.f / merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear[cluster.num].qMax) : 0.f;
337 float invAvgCharge = (sumInvSqrtCharge += invSqrtCharge) / ++nAvgCharge;
338 invAvgCharge *= invAvgCharge;
339 retVal = prop.Update(yy, zz, cluster.row, param, clusterState, rejectChi2, &interpolation.hit[ihit], refit, cluster.sector, time, invAvgCharge, invCharge GPUCA_DEBUG_STREAMER_CHECK(, &debugVals));
340 }
341 GPUCA_DEBUG_STREAMER_CHECK(if (o2::utils::DebugStreamer::checkStream(o2::utils::StreamFlags::streamUpdateTrack, iTrk)) {
342 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);
343 });
344 }
345 // clang-format off
346 CADEBUG(if (!CheckCov()) GPUError("INVALID COV AFTER UPDATE!!!"));
347 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));
348 // clang-format on
349
350 ConstrainSinPhi();
351 if (retVal == 0) // track is updated
352 {
353 if (storeOuter == 1 && cluster.leg == clusters[maxN - 1].leg) {
354 StoreOuter(outerParam, prop, 2);
355 storeOuter = 255;
356 }
357 noFollowCircle2 = false;
358 lastUpdateX = mX;
359 covYYUpd = mC[0];
360 nMissed = nMissed2 = 0;
361 UnmarkClusters(clusters, ihitMergeFirst, ihit, wayDirection, GPUTPCGMMergedTrackHit::flagNotFit);
362 N++;
363 ihitStart = ihit;
364 float dy = mP[0] - prop.Model().Y();
365 float dz = mP[1] - prop.Model().Z();
366 if (CAMath::Abs(mP[4]) * param.qptB5Scaler > 10 && --resetT0 <= 0 && CAMath::Abs(mP[2]) < 0.15f && dy * dy + dz * dz > 1) {
367 CADEBUG(printf("Reinit linearization\n"));
368 prop.SetTrack(this, prop.GetAlpha());
369 }
370 if (param.par.dodEdx && param.dodEdxDownscaled && iWay == nWays - 1 && cluster.leg == clusters[maxN - 1].leg) { // TODO: Costimize flag to remove, and option to remove double-clusters
371 bool acc = (clusterState & param.rec.tpc.dEdxClusterRejectionFlagMask) == 0, accAlt = (clusterState & param.rec.tpc.dEdxClusterRejectionFlagMaskAlt) == 0;
372 if (acc || accAlt) {
373 float qtot = 0, qmax = 0, pad = 0, relTime = 0;
374 const int32_t clusterCount = (ihit - ihitMergeFirst) * wayDirection + 1;
375 for (int32_t iTmp = ihitMergeFirst; iTmp != ihit + wayDirection; iTmp += wayDirection) {
376 if (merger->GetConstantMem()->ioPtrs.clustersNative == nullptr) {
377 qtot += clustersXYZ[ihit].amp;
378 } else {
379 const ClusterNative& cl = merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear[cluster.num];
380 qtot += cl.qTot;
381 qmax = CAMath::Max<float>(qmax, cl.qMax);
382 pad += cl.getPad();
383 relTime += cl.getTime();
384 }
385 }
386 qtot /= clusterCount; // TODO: Weighted Average
387 pad /= clusterCount;
388 relTime /= clusterCount;
389 relTime = relTime - CAMath::Round(relTime);
390 if (acc) {
391 dEdx.fillCluster(qtot, qmax, cluster.row, cluster.sector, mP[2], mP[3], merger->GetConstantMem()->calibObjects, zz, pad, relTime);
392 }
393 if (accAlt) {
394 dEdxAlt.fillCluster(qtot, qmax, cluster.row, cluster.sector, mP[2], mP[3], merger->GetConstantMem()->calibObjects, zz, pad, relTime);
395 }
396 }
397 }
398 } else if (retVal >= GPUTPCGMPropagator::updateErrorClusterRejected) { // cluster far away form the track
399 if (allowModification) {
400 MarkClusters(clusters, ihitMergeFirst, ihit, wayDirection, GPUTPCGMMergedTrackHit::flagRejectDistance);
401 } else if (iWay == nWays - 1) {
402 MarkClusters(clusters, ihitMergeFirst, ihit, wayDirection, GPUTPCGMMergedTrackHit::flagRejectErr);
403 }
404 nMissed++;
405 nMissed2++;
406 } else {
407 break; // bad chi2 for the whole track, stop the fit
408 }
409 }
410 if (((nWays - iWay) & 1) && (clusters[0].sector < 18) == (clusters[maxN - 1].sector < 18)) {
411 ShiftZ2(clusters, clustersXYZ, merger, maxN);
412 }
413 }
414 ConstrainSinPhi();
415
416 GPUCA_DEBUG_STREAMER_CHECK(if (o2::utils::DebugStreamer::checkStream(o2::utils::StreamFlags::streamUpdateTrack, iTrk)) {
417 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";
418 })
419
420 if (!(N + NTolerated >= GPUCA_TRACKLET_SELECTOR_MIN_HITS_B5(mP[4] * param.qptB5Scaler) && 2 * NTolerated <= CAMath::Max(10, N) && CheckNumericalQuality(covYYUpd))) {
421 return false; // TODO: NTolerated should never become that large, check what is going wrong!
422 }
423 if (param.rec.tpc.minNClustersFinalTrack != -1 && N + NTolerated < param.rec.tpc.minNClustersFinalTrack) {
424 return false;
425 }
426
427 // 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.
428
429 if (param.par.dodEdx && param.dodEdxDownscaled) {
430 dEdx.computedEdx(merger->OutputTracksdEdx()[iTrk], param);
431 dEdxAlt.computedEdx(merger->OutputTracksdEdxAlt()[iTrk], param);
432 }
433 Alpha = prop.GetAlpha();
434 MoveToReference(prop, param, Alpha);
435 NormalizeAlpha(Alpha);
436
437 return true;
438}
439
440GPUdni() void GPUTPCGMTrackParam::MoveToReference(GPUTPCGMPropagator& prop, const GPUParam& param, float& Alpha)
441{
442 static constexpr float kDeg2Rad = M_PI / 180.f;
443 static constexpr float kSectAngle = 2 * M_PI / 18.f;
444
445 if (param.rec.tpc.trackReferenceX <= 500) {
446 GPUTPCGMTrackParam save = *this;
447 float saveAlpha = Alpha;
448 for (int32_t attempt = 0; attempt < 3; attempt++) {
449 float dAngle = CAMath::Round(CAMath::ATan2(mP[0], mX) / kDeg2Rad / 20.f) * kSectAngle;
450 Alpha += dAngle;
451 if (prop.PropagateToXAlpha(param.rec.tpc.trackReferenceX, Alpha, 0)) {
452 break;
453 }
454 ConstrainSinPhi();
455 if (CAMath::Abs(mP[0]) <= mX * CAMath::Tan(kSectAngle / 2.f)) {
456 return;
457 }
458 }
459 *this = save;
460 Alpha = saveAlpha;
461 }
462 if (CAMath::Abs(mP[0]) > mX * CAMath::Tan(kSectAngle / 2.f)) {
463 float dAngle = CAMath::Round(CAMath::ATan2(mP[0], mX) / kDeg2Rad / 20.f) * kSectAngle;
464 Rotate(dAngle);
465 ConstrainSinPhi();
466 Alpha += dAngle;
467 }
468}
469
470GPUd() 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)
471{
472 if (mirrorParameters) {
473 prop.Mirror(inFlyDirection);
474 }
475 float err2Y, err2Z;
476 prop.GetErr2(err2Y, err2Z, param, toZ, row, clusterState, sector, -1.f, 0.f, 0.f); // Use correct time / avgCharge
477 prop.Model().Y() = mP[0] = toY;
478 prop.Model().Z() = mP[1] = toZ;
479 if (mC[0] < err2Y) {
480 mC[0] = err2Y;
481 }
482 if (mC[2] < err2Z) {
483 mC[2] = err2Z;
484 }
485 if (CAMath::Abs(mC[5]) < 0.1f) {
486 mC[5] = mC[5] > 0 ? 0.1f : -0.1f;
487 }
488 if (mC[9] < 1.f) {
489 mC[9] = 1.f;
490 }
491 mC[1] = mC[4] = mC[6] = mC[8] = mC[11] = mC[13] = 0;
492 prop.SetTrack(this, prop.GetAlpha());
493 mNDF = -3;
494 mChi2 = 0;
495}
496
497GPUd() 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)
498{
499 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) {
500 float maxDistY, maxDistZ;
501 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
502 maxDistY = (maxDistY + mC[0]) * 20.f;
503 maxDistZ = (maxDistZ + mC[2]) * 20.f;
504 int32_t noReject = 0; // Cannot reject if simple estimation of y/z fails (extremely unlike case)
505 if (CAMath::Abs(clAlpha - prop.GetAlpha()) > 1.e-4f) {
506 noReject = prop.RotateToAlpha(clAlpha);
507 }
508 float projY = 0, projZ = 0;
509 if (noReject == 0) {
510 noReject |= prop.GetPropagatedYZ(xx, projY, projZ);
511 }
512 float count = 0.f;
513 xx = yy = zz = 0.f;
514 clusterState = 0;
515 while (true) {
516 float clx, cly, clz, clamp;
517 if (merger->Param().par.earlyTpcTransform) {
518 const float zOffset = (clusters[ihit].sector < 18) == (clusters[0].sector < 18) ? mTZOffset : -mTZOffset;
519 clx = clustersXYZ[ihit].x;
520 cly = clustersXYZ[ihit].y;
521 clz = clustersXYZ[ihit].z - zOffset;
522 clamp = clustersXYZ[ihit].amp;
523 } else {
524 const ClusterNative& GPUrestrict() cl = merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear[clusters[ihit].num];
525 clamp = cl.qTot;
526 merger->GetConstantMem()->calibObjects.fastTransformHelper->Transform(clusters[ihit].sector, clusters[ihit].row, cl.getPad(), cl.getTime(), clx, cly, clz, mTZOffset);
527 }
528 float dy = cly - projY;
529 float dz = clz - projZ;
530 if (noReject == 0 && (dy * dy > maxDistY || dz * dz > maxDistZ)) {
531 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));
532 if (rejectChi2) {
534 }
535 } else {
536 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)));
537 xx += clx * clamp; // TODO: Weight in pad/time instead of XYZ
538 yy += cly * clamp;
539 zz += clz * clamp;
540 clusterState |= clusters[ihit].state;
541 count += clamp;
542 }
543 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)) {
544 break;
545 }
546 ihit += wayDirection;
547 }
548 if (count < 0.1f) {
549 CADEBUG(printf("\t\tNo matching cluster in double-row, skipping\n"));
550 return -1;
551 }
552 xx /= count;
553 yy /= count;
554 zz /= count;
555 }
556 return 0;
557}
558
559GPUd() float GPUTPCGMTrackParam::AttachClusters(const GPUTPCGMMerger* GPUrestrict() Merger, int32_t sector, int32_t iRow, int32_t iTrack, bool goodLeg, GPUTPCGMPropagator& prop)
560{
561 float Y, Z;
562 if (Merger->Param().par.earlyTpcTransform) {
563 Y = mP[0];
564 Z = mP[1];
565 } else {
566 float X = 0;
567 Merger->GetConstantMem()->calibObjects.fastTransformHelper->InverseTransformYZtoX(sector, iRow, mP[0], mP[1], X);
568 if (prop.GetPropagatedYZ(X, Y, Z)) {
569 Y = mP[0];
570 Z = mP[1];
571 }
572 }
573 return AttachClusters(Merger, sector, iRow, iTrack, goodLeg, Y, Z);
574}
575
576GPUd() float GPUTPCGMTrackParam::AttachClusters(const GPUTPCGMMerger* GPUrestrict() Merger, int32_t sector, int32_t iRow, int32_t iTrack, bool goodLeg, float Y, float Z)
577{
578 if (Merger->Param().rec.tpc.disableRefitAttachment & 1) {
579 return -1e6f;
580 }
581 const GPUTPCTracker& GPUrestrict() tracker = *(Merger -> GetConstantMem()->tpcTrackers + sector);
582 const GPUTPCRow& GPUrestrict() row = tracker.Row(iRow);
583#ifndef GPUCA_TEXTURE_FETCH_CONSTRUCTOR
584 GPUglobalref() const cahit2* hits = tracker.HitData(row);
585 GPUglobalref() const calink* firsthit = tracker.FirstHitInBin(row);
586#endif
587 if (row.NHits() == 0) {
588 return -1e6f;
589 }
590
591 const float zOffset = Merger->Param().par.earlyTpcTransform ? ((Merger->OutputTracks()[iTrack].CSide() ^ (sector >= 18)) ? -mTZOffset : mTZOffset) : Merger->GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->convVertexTimeToZOffset(sector, mTZOffset, Merger->Param().continuousMaxTimeBin);
592 const float y0 = row.Grid().YMin();
593 const float stepY = row.HstepY();
594 const float z0 = row.Grid().ZMin() - zOffset; // We can use our own ZOffset, since this is only used temporarily anyway
595 const float stepZ = row.HstepZ();
596 int32_t bin, ny, nz;
597
598 float err2Y, err2Z;
599 Merger->Param().GetClusterErrors2(sector, iRow, Z, mP[2], mP[3], -1.f, 0.f, 0.f, err2Y, err2Z); // TODO: Use correct time/avgCharge
600 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
601 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
602 const float tubeY = CAMath::Sqrt(sy2);
603 const float tubeZ = CAMath::Sqrt(sz2);
604 const float sy21 = 1.f / sy2;
605 const float sz21 = 1.f / sz2;
606 float uncorrectedY, uncorrectedZ;
607 if (Merger->Param().par.earlyTpcTransform) {
608 uncorrectedY = Y;
609 uncorrectedZ = Z;
610 } else {
611 Merger->GetConstantMem()->calibObjects.fastTransformHelper->InverseTransformYZtoNominalYZ(sector, iRow, Y, Z, uncorrectedY, uncorrectedZ);
612 }
613
614 if (CAMath::Abs(uncorrectedY) > row.getTPCMaxY()) {
615 return uncorrectedY;
616 }
617 row.Grid().GetBinArea(uncorrectedY, uncorrectedZ + zOffset, tubeY, tubeZ, bin, ny, nz);
618
619 const int32_t nBinsY = row.Grid().Ny();
620 const int32_t idOffset = tracker.Data().ClusterIdOffset();
621 const int32_t* ids = &(tracker.Data().ClusterDataIndex()[row.HitNumberOffset()]);
622 uint32_t myWeight = Merger->TrackOrderAttach()[iTrack] | gputpcgmmergertypes::attachAttached | gputpcgmmergertypes::attachTube;
623 GPUAtomic(uint32_t)* const weights = Merger->ClusterAttachment();
624 if (goodLeg) {
626 }
627 for (int32_t k = 0; k <= nz; k++) {
628 const int32_t mybin = bin + k * nBinsY;
629 const uint32_t hitFst = CA_TEXTURE_FETCH(calink, gAliTexRefu, firsthit, mybin);
630 const uint32_t hitLst = CA_TEXTURE_FETCH(calink, gAliTexRefu, firsthit, mybin + ny + 1);
631 for (uint32_t ih = hitFst; ih < hitLst; ih++) {
632 int32_t id = idOffset + ids[ih];
633 GPUAtomic(uint32_t)* const weight = weights + id;
634#if !defined(GPUCA_NO_ATOMIC_PRECHECK) && GPUCA_NO_ATOMIC_PRECHECK < 1
635 if (myWeight <= *weight) {
636 continue;
637 }
638#endif
639 const cahit2 hh = CA_TEXTURE_FETCH(cahit2, gAliTexRefu2, hits, ih);
640 const float y = y0 + hh.x * stepY;
641 const float z = z0 + hh.y * stepZ;
642 const float dy = y - uncorrectedY;
643 const float dz = z - uncorrectedZ;
644 if (dy * dy * sy21 + dz * dz * sz21 <= CAMath::Sqrt(2.f)) {
645 // CADEBUG(printf("Found Y %f Z %f\n", y, z));
646 CAMath::AtomicMax(weight, myWeight);
647 }
648 }
649 }
650 return uncorrectedY;
651}
652
653GPUd() 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)
654{
655 static constexpr float kSectAngle = 2 * M_PI / 18.f;
656 if (Merger->Param().rec.tpc.disableRefitAttachment & 2) {
657 return dodEdx;
658 }
659 if (CAMath::Abs(lastRow - toRow) < 2) {
660 return dodEdx;
661 }
662 int32_t step = toRow > lastRow ? 1 : -1;
663 float xx = mX - GPUTPCGeometry::Row2X(lastRow);
664 for (int32_t iRow = lastRow + step; iRow != toRow; iRow += step) {
665 if (CAMath::Abs(mP[2]) > maxSinPhi) {
666 return dodEdx;
667 }
668 if (CAMath::Abs(mP[0]) > CAMath::Abs(mX) * CAMath::Tan(kSectAngle / 2.f)) {
669 return dodEdx;
670 }
671 int32_t err = prop.PropagateToXAlpha(xx + GPUTPCGeometry::Row2X(iRow), prop.GetAlpha(), inFlyDirection);
672 if (err) {
673 return dodEdx;
674 }
675 if (dodEdx && iRow + step == toRow) {
676 float yUncorrected, zUncorrected;
677 Merger->GetConstantMem()->calibObjects.fastTransformHelper->InverseTransformYZtoNominalYZ(sector, iRow, mP[0], mP[1], yUncorrected, zUncorrected);
678 uint32_t pad = CAMath::Float2UIntRn(GPUTPCGeometry::LinearY2Pad(sector, iRow, yUncorrected));
679 if (pad >= GPUTPCGeometry::NPads(iRow) || (Merger->GetConstantMem()->calibObjects.dEdxCalibContainer && Merger->GetConstantMem()->calibObjects.dEdxCalibContainer->isDead(sector, iRow, pad))) {
680 dodEdx = false;
681 }
682 }
683 CADEBUG(printf("Attaching in row %d\n", iRow));
684 AttachClusters(Merger, sector, iRow, iTrack, goodLeg, prop);
685 }
686 return dodEdx;
687}
688
689GPUd() bool GPUTPCGMTrackParam::FollowCircleChk(float lrFactor, float toY, float toX, bool up, bool right)
690{
691 return CAMath::Abs(mX * lrFactor - toY) > 1.f && // transport further in Y
692 CAMath::Abs(mP[2]) < 0.7f && // rotate back
693 (up ? (-mP[0] * lrFactor > toX || (right ^ (mP[2] > 0))) : (-mP[0] * lrFactor < toX || (right ^ (mP[2] < 0)))); // don't overshoot in X
694}
695
696GPUdii() void GPUTPCGMTrackParam::StoreOuter(gputpcgmmergertypes::GPUTPCOuterParam* outerParam, const GPUTPCGMPropagator& prop, int32_t phase)
697{
698 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])));
699 for (int32_t i = 0; i < 5; i++) {
700 outerParam->P[i] = mP[i];
701 }
702 for (int32_t i = 0; i < 15; i++) {
703 outerParam->C[i] = mC[i];
704 }
705 outerParam->X = mX;
706 outerParam->alpha = prop.GetAlpha();
707}
708
709GPUdic(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)
710{
711 uint32_t nLoopData = CAMath::AtomicAdd(&Merger->Memory()->nLoopData, 1u);
712 if (nLoopData >= Merger->NMaxTracks()) {
713 Merger->raiseError(GPUErrors::ERROR_MERGER_LOOPER_OVERFLOW, nLoopData, Merger->NMaxTracks());
714 CAMath::AtomicExch(&Merger->Memory()->nLoopData, Merger->NMaxTracks());
715 return;
716 }
718 data.param = *this;
719 data.alpha = alpha;
720 data.track = iTrack;
721 data.toAlpha = toAlpha;
722 data.toY = toY;
723 data.toX = toX;
724 data.sector = sector;
725 data.row = iRow;
726 data.toSector = toSector;
727 data.toRow = toRow;
728 data.inFlyDirection = inFlyDirection;
729 Merger->LoopData()[nLoopData] = data;
730}
731
732GPUdii() void GPUTPCGMTrackParam::RefitLoop(const GPUTPCGMMerger* GPUrestrict() Merger, int32_t loopIdx)
733{
735 prop.SetMaterialTPC();
736 prop.SetPolynomialField(&Merger->Param().polynomialField);
737 prop.SetMaxSinPhi(GPUCA_MAX_SIN_PHI);
738 prop.SetToyMCEventsFlag(Merger->Param().par.toyMCEventsFlag);
739 prop.SetMatLUT(Merger->Param().rec.useMatLUT ? Merger->GetConstantMem()->calibObjects.matLUT : nullptr);
740 prop.SetSeedingErrors(false);
741 prop.SetFitInProjections(true);
742 prop.SetPropagateBzOnly(false);
743
744 GPUTPCGMLoopData& data = Merger->LoopData()[loopIdx];
745 prop.SetTrack(&data.param, data.alpha);
746 if (data.toSector == -1) {
747 data.param.AttachClustersMirror<1>(Merger, data.sector, data.row, data.track, data.toY, prop, true);
748 } else {
749 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);
750 }
751}
752
753template <int32_t I>
754GPUdic(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)
755{
756 static constexpr float kSectAngle = 2 * M_PI / 18.f;
757 if (Merger->Param().rec.tpc.disableRefitAttachment & 4) {
758 return 1;
759 }
760 if (Merger->Param().rec.tpc.looperInterpolationInExtraPass && phase2 == false) {
761 StoreAttachMirror(Merger, sector, iRow, iTrack, toAlpha, toY, toX, toSector, toRow, inFlyDirection, prop.GetAlpha());
762 return 1;
763 }
764 const GPUParam& GPUrestrict() param = Merger->Param();
765 bool right;
766 float dAlpha = toAlpha - prop.GetAlpha();
767 int32_t sectorSide = sector >= (GPUCA_NSECTORS / 2) ? (GPUCA_NSECTORS / 2) : 0;
768 if (CAMath::Abs(dAlpha) > 0.001f) {
769 right = CAMath::Abs(dAlpha) < CAMath::Pi() ? (dAlpha > 0) : (dAlpha < 0);
770 } else {
771 right = toY > mP[0];
772 }
773 bool up = (mP[2] < 0) ^ right;
774 int32_t targetRow = up ? (GPUCA_ROW_COUNT - 1) : 0;
775 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...
776 // clang-format off
777 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));
778 // clang-format on
779
780 AttachClustersPropagate(Merger, sector, iRow, targetRow, iTrack, false, prop, inFlyDirection, 0.7f);
781 if (prop.RotateToAlpha(prop.GetAlpha() + (CAMath::Pi() / 2.f) * lrFactor)) {
782 return 1;
783 }
784 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));
785 while (sector != toSector || FollowCircleChk(lrFactor, toY, toX, up, right)) {
786 while ((sector != toSector) ? (CAMath::Abs(mX) <= CAMath::Abs(mP[0]) * CAMath::Tan(kSectAngle / 2.f)) : FollowCircleChk(lrFactor, toY, toX, up, right)) {
787 int32_t err = prop.PropagateToXAlpha(mX + 1.f, prop.GetAlpha(), inFlyDirection);
788 if (err) {
789 CADEBUG(printf("\t\tpropagation error (%d)\n", err));
790 prop.RotateToAlpha(prop.GetAlpha() - (CAMath::Pi() / 2.f) * lrFactor);
791 return 1;
792 }
793 CADEBUG(printf("\tPropagated to y = %f: X %f Z %f SinPhi %f\n", mX, mP[0], mP[1], mP[2]));
794 for (int32_t j = 0; j < GPUCA_ROW_COUNT; j++) {
795 float rowX = GPUTPCGeometry::Row2X(j);
796 if (CAMath::Abs(rowX - (-mP[0] * lrFactor)) < 1.5f) {
797 CADEBUG(printf("\t\tAttempt row %d (Y %f Z %f)\n", j, mX * lrFactor, mP[1]));
798 AttachClusters(Merger, sector, j, iTrack, false, mX * lrFactor, mP[1]);
799 }
800 }
801 }
802 if (sector != toSector) {
803 if (right) {
804 if (++sector >= sectorSide + 18) {
805 sector -= 18;
806 }
807 } else {
808 if (--sector < sectorSide) {
809 sector += 18;
810 }
811 }
812 CADEBUG(printf("\tRotating to sector %d\n", sector));
813 if (prop.RotateToAlpha(param.Alpha(sector) + (CAMath::Pi() / 2.f) * lrFactor)) {
814 CADEBUG(printf("\t\trotation error\n"));
815 prop.RotateToAlpha(prop.GetAlpha() - (CAMath::Pi() / 2.f) * lrFactor);
816 return 1;
817 }
818 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]));
819 }
820 }
821 CADEBUG(printf("\tRotating back\n"));
822 for (int32_t i = 0; i < 2; i++) {
823 if (prop.RotateToAlpha(prop.GetAlpha() + (CAMath::Pi() / 2.f) * lrFactor) == 0) {
824 break;
825 }
826 if (i) {
827 CADEBUG(printf("Final rotation failed\n"));
828 return 1;
829 }
830 CADEBUG(printf("\tresetting physical model\n"));
831 prop.SetTrack(this, prop.GetAlpha());
832 }
833 prop.Rotate180();
834 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]));
835 iRow = toRow;
836 float dx = toX - GPUTPCGeometry::Row2X(toRow);
837 if (up ^ (toX > mX)) {
838 if (up) {
839 while (iRow < GPUCA_ROW_COUNT - 2 && GPUTPCGeometry::Row2X(iRow + 1) + dx <= mX) {
840 iRow++;
841 }
842 } else {
843 while (iRow > 1 && GPUTPCGeometry::Row2X(iRow - 1) + dx >= mX) {
844 iRow--;
845 }
846 }
847 prop.PropagateToXAlpha(GPUTPCGeometry::Row2X(iRow) + dx, prop.GetAlpha(), inFlyDirection);
848 AttachClustersPropagate(Merger, sector, iRow, toRow, iTrack, false, prop, inFlyDirection);
849 }
850 if (prop.PropagateToXAlpha(toX, prop.GetAlpha(), inFlyDirection)) {
851 mX = toX;
852 }
853 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]));
854 return (0);
855}
856
857template <int32_t I>
858GPUdni() void GPUTPCGMTrackParam::AttachClustersMirror(const GPUTPCGMMerger* GPUrestrict() Merger, int32_t sector, int32_t iRow, int32_t iTrack, float toY, GPUTPCGMPropagator& GPUrestrict() prop, bool phase2)
859{
860 static constexpr float kSectAngle = 2 * M_PI / 18.f;
861
862 if (Merger->Param().rec.tpc.disableRefitAttachment & 8) {
863 return;
864 }
865 if (Merger->Param().rec.tpc.looperInterpolationInExtraPass && phase2 == false) {
866 StoreAttachMirror(Merger, sector, iRow, iTrack, 0, toY, 0, -1, 0, 0, prop.GetAlpha());
867 return;
868 }
869 // Note that the coordinate system is rotated by 90 degree swapping X and Y!
870 float X = mP[2] > 0 ? mP[0] : -mP[0];
871 float toX = mP[2] > 0 ? toY : -toY;
872 float Y = mP[2] > 0 ? -mX : mX;
873 float Z = mP[1];
874 if (CAMath::Abs(mP[2]) >= GPUCA_MAX_SIN_PHI_LOW) {
875 return;
876 }
877 float SinPhi = CAMath::Sqrt(1 - mP[2] * mP[2]) * (mP[2] > 0 ? -1 : 1);
878 if (CAMath::Abs(SinPhi) >= GPUCA_MAX_SIN_PHI_LOW) {
879 return;
880 }
881 float b = prop.GetBz(prop.GetAlpha(), mX, mP[0], mP[1]);
882
883 int32_t count = CAMath::Float2IntRn(CAMath::Abs((toX - X) * 2.f));
884 if (count == 0) {
885 return;
886 }
887 float dx = (toX - X) / count;
888 const float myRowX = GPUTPCGeometry::Row2X(iRow);
889 // printf("AttachMirror\n");
890 // printf("X %f Y %f Z %f SinPhi %f toY %f -->\n", mX, mP[0], mP[1], mP[2], toY);
891 // 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);
892 while (count--) {
893 float ex = CAMath::Sqrt(1 - SinPhi * SinPhi);
894 float exi = 1.f / ex;
895 float dxBzQ = dx * -b * mP[4];
896 float newSinPhi = SinPhi + dxBzQ;
897 if (CAMath::Abs(newSinPhi) > GPUCA_MAX_SIN_PHI_LOW) {
898 return;
899 }
900 float dS = dx * exi;
901 float h2 = dS * exi * exi;
902 float h4 = .5f * h2 * dxBzQ;
903
904 X += dx;
905 Y += dS * SinPhi + h4;
906 Z += dS * mP[3];
907 SinPhi = newSinPhi;
908 if (CAMath::Abs(X) > CAMath::Abs(Y) * CAMath::Tan(kSectAngle / 2.f)) {
909 continue;
910 }
911
912 // 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);
913
914 float paramX = mP[2] > 0 ? -Y : Y;
915 int32_t step = paramX >= mX ? 1 : -1;
916 int32_t found = 0;
917 for (int32_t j = iRow; j >= 0 && j < GPUCA_ROW_COUNT && found < 3; j += step) {
918 float rowX = mX + GPUTPCGeometry::Row2X(j) - myRowX;
919 if (CAMath::Abs(rowX - paramX) < 1.5f) {
920 // printf("Attempt row %d\n", j);
921 AttachClusters(Merger, sector, j, iTrack, false, mP[2] > 0 ? X : -X, Z);
922 }
923 }
924 }
925}
926
927GPUd() void GPUTPCGMTrackParam::ShiftZ2(const GPUTPCGMMergedTrackHit* clusters, GPUTPCGMMergedTrackHitXYZ* clustersXYZ, const GPUTPCGMMerger* merger, int32_t N)
928{
929 float tzInner, tzOuter;
930 float xInner, xOuter;
931 if (N == 0) {
932 N = 1;
933 }
934 if (merger->Param().par.earlyTpcTransform) {
935 tzInner = clustersXYZ[N - 1].z;
936 tzOuter = clustersXYZ[0].z;
937 xInner = clustersXYZ[N - 1].x;
938 xOuter = clustersXYZ[0].x;
939 } else {
940 const auto& GPUrestrict() cls = merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear;
941 tzInner = cls[clusters[N - 1].num].getTime();
942 tzOuter = cls[clusters[0].num].getTime();
943 xInner = GPUTPCGeometry::Row2X(clusters[N - 1].row);
944 xOuter = GPUTPCGeometry::Row2X(clusters[0].row);
945 }
946 ShiftZ(merger, clusters[0].sector, tzInner, tzOuter, xInner, xOuter);
947}
948
949GPUd() void GPUTPCGMTrackParam::ShiftZ(const GPUTPCGMMerger* GPUrestrict() merger, int32_t sector, float tz1, float tz2, float x1, float x2)
950{
951 if (!merger->Param().par.continuousTracking) {
952 return;
953 }
954 const float r1 = CAMath::Max(0.0001f, CAMath::Abs(mP[4] * merger->Param().polynomialField.GetNominalBz()));
955
956 const float dist2 = mX * mX + mP[0] * mP[0];
957 const float dist1r2 = dist2 * r1 * r1;
958 float deltaZ = 0.f;
959 bool beamlineReached = false;
960 if (dist1r2 < 4) {
961 const float alpha = CAMath::ACos(1 - 0.5f * dist1r2); // Angle of a circle, such that |(cosa, sina) - (1,0)| == dist
962 const float beta = CAMath::ATan2(mP[0], mX);
963 const int32_t comp = mP[2] > CAMath::Sin(beta);
964 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
965 const float res = CAMath::Abs(sinab - mP[2]);
966
967 if (res < 0.2) {
968 const float r = 1.f / r1;
969 const float dS = alpha * r;
970 float z0 = dS * mP[3];
971 if (CAMath::Abs(z0) > GPUTPCGeometry::TPCLength()) {
972 z0 = z0 > 0 ? GPUTPCGeometry::TPCLength() : -GPUTPCGeometry::TPCLength();
973 }
974 deltaZ = mP[1] - z0;
975 beamlineReached = true;
976
977 // 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);
978 }
979 }
980
981 if (!beamlineReached) {
982 if (merger->Param().par.earlyTpcTransform) {
983 float basez, basex;
984 if (CAMath::Abs(tz1) < CAMath::Abs(tz2)) {
985 basez = tz1;
986 basex = x1;
987 } else {
988 basez = tz2;
989 basex = x2;
990 }
991 float refZ = ((basez > 0) ? merger->Param().rec.tpc.defaultZOffsetOverR : -merger->Param().rec.tpc.defaultZOffsetOverR) * basex;
992 deltaZ = basez - refZ - mTZOffset;
993 } else {
994 float baset, basex;
995 if (CAMath::Abs(tz1) > CAMath::Abs(tz2)) {
996 baset = tz1;
997 basex = x1;
998 } else {
999 baset = tz2;
1000 basex = x2;
1001 }
1002 float refZ = ((sector < GPUCA_NSECTORS / 2) ? merger->Param().rec.tpc.defaultZOffsetOverR : -merger->Param().rec.tpc.defaultZOffsetOverR) * basex;
1003 float basez;
1004 merger->GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->TransformIdealZ(sector, baset, basez, mTZOffset);
1005 deltaZ = basez - refZ;
1006 }
1007 }
1008 if (merger->Param().par.earlyTpcTransform) {
1009 mTZOffset += deltaZ;
1010 mP[1] -= deltaZ;
1011 deltaZ = 0;
1012 float zMax = CAMath::Max(tz1, tz2);
1013 float zMin = CAMath::Min(tz1, tz2);
1014 // printf("Z Check: Clusters %f %f, min %f max %f vtx %f\n", tz1, tz2, zMin, zMax, mTZOffset);
1015 if (zMin < 0 && zMin - mTZOffset < -GPUTPCGeometry::TPCLength()) {
1016 deltaZ = zMin - mTZOffset + GPUTPCGeometry::TPCLength();
1017 } else if (zMax > 0 && zMax - mTZOffset > GPUTPCGeometry::TPCLength()) {
1018 deltaZ = zMax - mTZOffset - GPUTPCGeometry::TPCLength();
1019 }
1020 if (zMin < 0 && zMax - (mTZOffset + deltaZ) > 0) {
1021 deltaZ = zMax - mTZOffset;
1022 } else if (zMax > 0 && zMin - (mTZOffset + deltaZ) < 0) {
1023 deltaZ = zMin - mTZOffset;
1024 }
1025 // 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);
1026 mTZOffset += deltaZ;
1027 mP[1] -= deltaZ;
1028 } else {
1029 float deltaT = merger->GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->convDeltaZtoDeltaTimeInTimeFrame(sector, deltaZ);
1030 mTZOffset += deltaT;
1031 mP[1] -= deltaZ;
1032 const float maxT = CAMath::Min(tz1, tz2) - merger->GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->getT0();
1033 const float minT = CAMath::Max(tz1, tz2) - merger->GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->getMaxDriftTime(sector);
1034 // printf("T Check: Clusters %f %f, min %f max %f vtx %f\n", tz1, tz2, minT, maxT, mTZOffset);
1035 deltaT = 0.f;
1036 if (mTZOffset < minT) {
1037 deltaT = minT - mTZOffset;
1038 }
1039 if (mTZOffset + deltaT > maxT) {
1040 deltaT = maxT - mTZOffset;
1041 }
1042 if (deltaT != 0.f) {
1043 deltaZ = merger->GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->convDeltaTimeToDeltaZinTimeFrame(sector, deltaT);
1044 // 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);
1045 mTZOffset += deltaT;
1046 mP[1] -= deltaZ;
1047 }
1048 }
1049 // printf("\n");
1050}
1051
1052GPUd() bool GPUTPCGMTrackParam::CheckCov() const
1053{
1054 const float* c = mC;
1055 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]) &&
1056 (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]);
1057 return ok;
1058}
1059
1060GPUd() bool GPUTPCGMTrackParam::CheckNumericalQuality(float overrideCovYY) const
1061{
1062 //* Check that the track parameters and covariance matrix are reasonable
1063 bool ok = CAMath::Finite(mX) && CAMath::Finite(mChi2);
1064 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"));
1065 const float* c = mC;
1066 for (int32_t i = 0; i < 15; i++) {
1067 ok = ok && CAMath::Finite(c[i]);
1068 }
1069 CADEBUG(printf("OK1 %d\n", (int32_t)ok));
1070 for (int32_t i = 0; i < 5; i++) {
1071 ok = ok && CAMath::Finite(mP[i]);
1072 }
1073 CADEBUG(printf("OK2 %d\n", (int32_t)ok));
1074 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) {
1075 ok = 0;
1076 }
1077 CADEBUG(printf("OK3 %d\n", (int32_t)ok));
1078 if (CAMath::Abs(mP[2]) > GPUCA_MAX_SIN_PHI) {
1079 ok = 0;
1080 }
1081 CADEBUG(printf("OK4 %d\n", (int32_t)ok));
1082 if (!CheckCov()) {
1083 ok = false;
1084 }
1085 CADEBUG(printf("OK5 %d\n", (int32_t)ok));
1086 return ok;
1087}
1088
1089GPUd() void GPUTPCGMTrackParam::RefitTrack(GPUTPCGMMergedTrack& GPUrestrict() track, int32_t iTrk, GPUTPCGMMerger* GPUrestrict() merger, int32_t attempt) // TODO: Inline me, once __forceinline__ is fixed by HIP
1090{
1091 if (!track.OK()) {
1092 return;
1093 }
1094
1095 // clang-format off
1096 CADEBUG(if (DEBUG_SINGLE_TRACK >= 0 && iTrk != DEBUG_SINGLE_TRACK) { track.SetNClusters(0); track.SetOK(0); return; } );
1097 // clang-format on
1098
1099 int32_t nTrackHits = track.NClusters();
1100 int32_t NTolerated = 0; // Clusters not fit but tollerated for track length cut
1101 GPUTPCGMTrackParam t = track.Param();
1102 float Alpha = track.Alpha();
1103 CADEBUG(int32_t nTrackHitsOld = nTrackHits; float ptOld = t.QPt());
1104 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());
1105 CADEBUG(printf("Finished Fit Track %d\n", iTrk));
1106 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)));
1107
1108 if (!ok && attempt == 0 && merger->Param().rec.tpc.retryRefit) {
1109 for (uint32_t i = 0; i < track.NClusters(); i++) {
1110 merger->Clusters()[track.FirstClusterRef() + i].state &= GPUTPCGMMergedTrackHit::clustererAndSharedFlags;
1111 }
1112 CADEBUG(printf("Track rejected, marking for retry\n"));
1113 if (merger->Param().rec.tpc.retryRefit == 2) {
1114 nTrackHits = track.NClusters();
1115 NTolerated = 0; // Clusters not fit but tollerated for track length cut
1116 t = track.Param();
1117 Alpha = track.Alpha();
1118 ok = t.Fit(merger, iTrk, merger->Clusters() + track.FirstClusterRef(), merger->ClustersXYZ() + track.FirstClusterRef(), nTrackHits, NTolerated, Alpha, 1, GPUCA_MAX_SIN_PHI, &track.OuterParam());
1119 } else {
1120 uint32_t nRefit = CAMath::AtomicAdd(&merger->Memory()->nRetryRefit, 1u);
1121 merger->RetryRefitIds()[nRefit] = iTrk;
1122 return;
1123 }
1124 }
1125 if (CAMath::Abs(t.QPt()) < 1.e-4f) {
1126 t.QPt() = 1.e-4f;
1127 }
1128
1129 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->OutputTracks()[iTrk].Looper()); });
1130
1131 track.SetOK(ok);
1132 track.SetNClustersFitted(nTrackHits);
1133 track.Param() = t;
1134 track.Alpha() = Alpha;
1135
1136 if (track.OK()) {
1137 int32_t ind = track.FirstClusterRef();
1138 const GPUParam& GPUrestrict() param = merger->Param();
1139 float alphaa = param.Alpha(merger->Clusters()[ind].sector);
1140 float xx, yy, zz;
1141 if (merger->Param().par.earlyTpcTransform) {
1142 xx = merger->ClustersXYZ()[ind].x;
1143 yy = merger->ClustersXYZ()[ind].y;
1144 zz = merger->ClustersXYZ()[ind].z - track.Param().GetTZOffset();
1145 } else {
1146 const ClusterNative& GPUrestrict() cl = merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear[merger->Clusters()[ind].num];
1147 merger->GetConstantMem()->calibObjects.fastTransformHelper->Transform(merger->Clusters()[ind].sector, merger->Clusters()[ind].row, cl.getPad(), cl.getTime(), xx, yy, zz, track.Param().GetTZOffset());
1148 }
1149 float sinA, cosA;
1150 CAMath::SinCos(alphaa - track.Alpha(), sinA, cosA);
1151 track.SetLastX(xx * cosA - yy * sinA);
1152 track.SetLastY(xx * sinA + yy * cosA);
1153 track.SetLastZ(zz);
1154 // merger->DebugRefitMergedTrack(track);
1155 }
1156}
1157
1158GPUd() void GPUTPCGMTrackParam::Rotate(float alpha)
1159{
1160 float cA, sA;
1161 CAMath::SinCos(alpha, sA, cA);
1162 float x0 = mX;
1163 float sinPhi0 = mP[2], cosPhi0 = CAMath::Sqrt(1 - mP[2] * mP[2]);
1164 float cosPhi = cosPhi0 * cA + sinPhi0 * sA;
1165 float sinPhi = -cosPhi0 * sA + sinPhi0 * cA;
1166 float j0 = cosPhi0 / cosPhi;
1167 float j2 = cosPhi / cosPhi0;
1168 mX = x0 * cA + mP[0] * sA;
1169 mP[0] = -x0 * sA + mP[0] * cA;
1170 mP[2] = sinPhi;
1171 mC[0] *= j0 * j0;
1172 mC[1] *= j0;
1173 mC[3] *= j0;
1174 mC[6] *= j0;
1175 mC[10] *= j0;
1176
1177 mC[3] *= j2;
1178 mC[4] *= j2;
1179 mC[5] *= j2 * j2;
1180 mC[8] *= j2;
1181 mC[12] *= j2;
1182 if (cosPhi < 0) { // change direction ( t0 direction is already changed in t0.UpdateValues(); )
1183 SinPhi() = -SinPhi();
1184 DzDs() = -DzDs();
1185 QPt() = -QPt();
1186 mC[3] = -mC[3];
1187 mC[4] = -mC[4];
1188 mC[6] = -mC[6];
1189 mC[7] = -mC[7];
1190 mC[10] = -mC[10];
1191 mC[11] = -mC[11];
1192 }
1193}
1194
1195GPUd() void GPUTPCGMTrackParam::AddCovDiagErrors(const float* GPUrestrict() errors2)
1196{
1197 mC[0] += errors2[0];
1198 mC[2] += errors2[1];
1199 mC[5] += errors2[2];
1200 mC[9] += errors2[3];
1201 mC[14] += errors2[4];
1202}
1203
1204GPUd() void GPUTPCGMTrackParam::AddCovDiagErrorsWithCorrelations(const float* GPUrestrict() errors2)
1205{
1206 const int32_t diagMap[5] = {0, 2, 5, 9, 14};
1207 const float oldDiag[5] = {mC[0], mC[2], mC[5], mC[9], mC[14]};
1208 for (int32_t i = 0; i < 5; i++) {
1209 mC[diagMap[i]] += errors2[i];
1210 for (int32_t j = 0; j < i; j++) {
1211 mC[diagMap[i - 1] + j + 1] *= gpu::CAMath::Sqrt(mC[diagMap[i]] * mC[diagMap[j]] / (oldDiag[i] * oldDiag[j]));
1212 }
1213 }
1214}
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_MAX_SIN_PHI_LOW
#define GPUCA_TRACKLET_SELECTOR_MIN_HITS_B5(QPTB5)
#define GPUCA_MAX_SIN_PHI
#define GPUCA_MERGER_INTERPOLATION_ERROR_TYPE_A
#define CADEBUG(...)
Definition GPUDef.h:85
#define CA_TEXTURE_FETCH(type, texture, address, entry)
Definition GPUDef.h:64
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_MERGER_INTERPOLATION_ERROR_TYPE_A errorZ
GPUCA_MERGER_INTERPOLATION_ERROR_TYPE_A errorY
InterpolationErrorHit hit[GPUCA_MERGER_MAX_TRACK_CLUSTERS]
std::vector< Cluster > clusters
std::vector< int > row