Project
Loading...
Searching...
No Matches
GPUTPCGMSectorTrack.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#include "GPUParam.h"
16#include "GPUTPCGMBorderTrack.h"
17#include "GPUTPCGMSectorTrack.h"
18#include "GPUO2DataTypes.h"
19#include "GPUTPCGMMerger.h"
20#include "GPUTPCConvertImpl.h"
21#include "GPUParam.inc"
22
23using namespace o2::gpu;
24using namespace o2::tpc;
25
26GPUd() void GPUTPCGMSectorTrack::Set(const GPUTPCGMMerger* merger, const GPUTPCTrack* sectorTr, float alpha, int32_t sector)
27{
28 const GPUTPCBaseTrackParam& t = sectorTr->Param();
29 mOrigTrack = sectorTr;
30 mParam.mX = t.GetX();
31 mParam.mY = t.GetY();
32 mParam.mZ = t.GetZ();
33 mParam.mDzDs = t.GetDzDs();
34 mParam.mSinPhi = t.GetSinPhi();
35 mParam.mQPt = t.GetQPt();
36 mParam.mCosPhi = CAMath::Sqrt(1.f - mParam.mSinPhi * mParam.mSinPhi);
37 mParam.mSecPhi = 1.f / mParam.mCosPhi;
38 mAlpha = alpha;
39 mSector = sector;
40 mTOffset = merger->GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->convZOffsetToVertexTime(sector, t.GetZOffset(), merger->Param().continuousMaxTimeBin);
41 mNClusters = sectorTr->NHits();
42}
43
44GPUd() void GPUTPCGMSectorTrack::Set(const GPUTPCGMTrackParam& trk, const GPUTPCTrack* sectorTr, float alpha, int32_t sector)
45{
46 mOrigTrack = sectorTr;
47 mParam.mX = trk.GetX();
48 mParam.mY = trk.GetY();
49 mParam.mZ = trk.GetZ();
50 mParam.mDzDs = trk.GetDzDs();
51 mParam.mSinPhi = trk.GetSinPhi();
52 mParam.mQPt = trk.GetQPt();
53 mParam.mCosPhi = CAMath::Sqrt(1.f - mParam.mSinPhi * mParam.mSinPhi);
54 mParam.mSecPhi = 1.f / mParam.mCosPhi;
55 mAlpha = alpha;
56 mSector = sector;
57 mTOffset = trk.GetTOffset();
58 mNClusters = sectorTr->NHits();
59 mParam.mC0 = trk.GetCov(0);
60 mParam.mC2 = trk.GetCov(2);
61 mParam.mC3 = trk.GetCov(3);
62 mParam.mC5 = trk.GetCov(5);
63 mParam.mC7 = trk.GetCov(7);
64 mParam.mC9 = trk.GetCov(9);
65 mParam.mC10 = trk.GetCov(10);
66 mParam.mC12 = trk.GetCov(12);
67 mParam.mC14 = trk.GetCov(14);
68}
69
70GPUd() void GPUTPCGMSectorTrack::SetParam2(const GPUTPCGMTrackParam& trk)
71{
72 mParam2.mX = trk.GetX();
73 mParam2.mY = trk.GetY();
74 mParam2.mZ = trk.GetZ();
75 mParam2.mDzDs = trk.GetDzDs();
76 mParam2.mSinPhi = trk.GetSinPhi();
77 mParam2.mQPt = trk.GetQPt();
78 mParam2.mCosPhi = CAMath::Sqrt(1.f - mParam2.mSinPhi * mParam2.mSinPhi);
79 mParam2.mSecPhi = 1.f / mParam2.mCosPhi;
80 mParam2.mC0 = trk.GetCov(0);
81 mParam2.mC2 = trk.GetCov(2);
82 mParam2.mC3 = trk.GetCov(3);
83 mParam2.mC5 = trk.GetCov(5);
84 mParam2.mC7 = trk.GetCov(7);
85 mParam2.mC9 = trk.GetCov(9);
86 mParam2.mC10 = trk.GetCov(10);
87 mParam2.mC12 = trk.GetCov(12);
88 mParam2.mC14 = trk.GetCov(14);
89}
90
91GPUd() bool GPUTPCGMSectorTrack::FilterErrors(const GPUTPCGMMerger* merger, int32_t iSector, float maxSinPhi, float sinPhiMargin)
92{
93 float lastX;
94 // float lastX = GPUTPCGeometry::Row2X(mOrigTrack->Cluster(mOrigTrack->NClusters() - 1).GetRow()); // TODO: Why is this needed to be set below, Row2X should work, but looses some tracks
95 float y, z;
96 int32_t row, index;
97 const GPUTPCTracker& trk = merger->GetConstantMem()->tpcTrackers[iSector];
98 const GPUTPCHitId& ic = trk.TrackHits()[mOrigTrack->FirstHitID() + mOrigTrack->NHits() - 1];
99 index = trk.Data().ClusterDataIndex(trk.Data().Row(ic.RowIndex()), ic.HitIndex()) + merger->GetConstantMem()->ioPtrs.clustersNative->clusterOffset[iSector][0];
100 row = ic.RowIndex();
101 const ClusterNative& cl = merger->GetConstantMem()->ioPtrs.clustersNative->clustersLinear[index];
102 GPUTPCConvertImpl::convert(*merger->GetConstantMem(), iSector, row, cl.getPad(), cl.getTime(), lastX, y, z);
103
104 const int32_t N = 3;
105
106 float bz = -merger->Param().bzCLight;
107
108 float k = mParam.mQPt * bz;
109 float dx = (1.f / N) * (lastX - mParam.mX);
110 float kdx = k * dx;
111 float dxBz = dx * bz;
112 float kdx205 = 2.f + kdx * kdx * 0.5f;
113
114 {
115 merger->Param().GetClusterErrors2(iSector, 0, mParam.mZ, mParam.mSinPhi, mParam.mDzDs, -1.f, 0.f, 0.f, mParam.mC0, mParam.mC2); // TODO: provide correct time and row
116#ifndef GPUCA_TPC_GEOMETRY_O2
117 float C0a, C2a;
118 merger->Param().GetClusterErrorsSeeding2(iSector, 0, mParam.mZ, mParam.mSinPhi, mParam.mDzDs, -1.f, C0a, C2a);
119 if (C0a > mParam.mC0) {
120 mParam.mC0 = C0a;
121 }
122 if (C2a > mParam.mC2) {
123 mParam.mC2 = C2a;
124 }
125#endif
126
127 mParam.mC3 = 0;
128 mParam.mC5 = 1;
129 mParam.mC7 = 0;
130 mParam.mC9 = 1;
131 mParam.mC10 = 0;
132 mParam.mC12 = 0;
133 mParam.mC14 = 10;
134 }
135
136 for (int32_t iStep = 0; iStep < N; iStep++) {
137 float err2Y, err2Z;
138
139 { // transport block
140 float ex = mParam.mCosPhi;
141 float ey = mParam.mSinPhi;
142 float ey1 = kdx + ey;
143 if (CAMath::Abs(ey1) > maxSinPhi) {
144 if (ey1 > maxSinPhi && ey1 < maxSinPhi + sinPhiMargin) {
145 ey1 = maxSinPhi - 0.01f;
146 } else if (ey1 > -maxSinPhi - sinPhiMargin) {
147 ey1 = -maxSinPhi + 0.01f;
148 } else {
149 return 0;
150 }
151 }
152
153 float ss = ey + ey1;
154 float ex1 = CAMath::Sqrt(1.f - ey1 * ey1);
155
156 float cc = ex + ex1;
157 float dxcci = dx / cc;
158
159 float dy = dxcci * ss;
160 float norm2 = 1.f + ey * ey1 + ex * ex1;
161 float dl = dxcci * CAMath::Sqrt(norm2 + norm2);
162
163 float dS;
164 {
165 float dSin = 0.5f * k * dl;
166 float a = dSin * dSin;
167 const float k2 = 1.f / 6.f;
168 const float k4 = 3.f / 40.f;
169 dS = dl + dl * a * (k2 + a * (k4)); //+ k6*a) );
170 }
171
172 float dz = dS * mParam.mDzDs;
173 float ex1i = 1.f / ex1;
174 {
175 merger->Param().GetClusterErrors2(iSector, 0, mParam.mZ, mParam.mSinPhi, mParam.mDzDs, -1.f, 0.f, 0.f, err2Y, err2Z); // TODO: Provide correct time / row
176#ifndef GPUCA_TPC_GEOMETRY_O2
177 float C0a, C2a;
178 merger->Param().GetClusterErrorsSeeding2(iSector, 0, mParam.mZ, mParam.mSinPhi, mParam.mDzDs, -1.f, C0a, C2a);
179 if (C0a > err2Y) {
180 err2Y = C0a;
181 }
182 if (C2a > err2Z) {
183 err2Z = C2a;
184 }
185#endif
186 }
187
188 float hh = kdx205 * dxcci * ex1i;
189 float h2 = hh * mParam.mSecPhi;
190
191 mParam.mX += dx;
192 mParam.mY += dy;
193 mParam.mZ += dz;
194 mParam.mSinPhi = ey1;
195 mParam.mCosPhi = ex1;
196 mParam.mSecPhi = ex1i;
197
198 float h4 = bz * dxcci * hh;
199
200 float c20 = mParam.mC3;
201 float c22 = mParam.mC5;
202 float c31 = mParam.mC7;
203 float c33 = mParam.mC9;
204 float c40 = mParam.mC10;
205 float c42 = mParam.mC12;
206 float c44 = mParam.mC14;
207
208 float c20ph4c42 = c20 + h4 * c42;
209 float h2c22 = h2 * c22;
210 float h4c44 = h4 * c44;
211 float n7 = c31 + dS * c33;
212 float n10 = c40 + h2 * c42 + h4c44;
213 float n12 = c42 + dxBz * c44;
214
215 mParam.mC0 += h2 * h2c22 + h4 * h4c44 + 2.f * (h2 * c20ph4c42 + h4 * c40);
216
217 mParam.mC3 = c20ph4c42 + h2c22 + dxBz * n10;
218 mParam.mC10 = n10;
219
220 mParam.mC5 = c22 + dxBz * (c42 + n12);
221 mParam.mC12 = n12;
222
223 mParam.mC2 += dS * (c31 + n7);
224 mParam.mC7 = n7;
225 } // end transport block
226
227 // Filter block
228
229 float c00 = mParam.mC0, c11 = mParam.mC2, c20 = mParam.mC3, c31 = mParam.mC7, c40 = mParam.mC10;
230
231 float mS0 = 1.f / (err2Y + c00);
232 float mS2 = 1.f / (err2Z + c11);
233
234 // K = CHtS
235
236 float k00, k11, k20, k31, k40;
237
238 k00 = c00 * mS0;
239 k20 = c20 * mS0;
240 k40 = c40 * mS0;
241
242 mParam.mC0 -= k00 * c00;
243 mParam.mC5 -= k20 * c20;
244 mParam.mC10 -= k00 * c40;
245 mParam.mC12 -= k40 * c20;
246 mParam.mC3 -= k20 * c00;
247 mParam.mC14 -= k40 * c40;
248
249 k11 = c11 * mS2;
250 k31 = c31 * mS2;
251
252 mParam.mC7 -= k31 * c11;
253 mParam.mC2 -= k11 * c11;
254 mParam.mC9 -= k31 * c31;
255 }
256
257 //* Check that the track parameters and covariance matrix are reasonable
258
259 bool ok = CAMath::Finite(mParam.mX) && CAMath::Finite(mParam.mY) && CAMath::Finite(mParam.mZ) && CAMath::Finite(mParam.mSinPhi) && CAMath::Finite(mParam.mDzDs) && CAMath::Finite(mParam.mQPt) && CAMath::Finite(mParam.mCosPhi) && CAMath::Finite(mParam.mSecPhi) && CAMath::Finite(mTOffset) && CAMath::Finite(mParam.mC0) && CAMath::Finite(mParam.mC2) &&
260 CAMath::Finite(mParam.mC3) && CAMath::Finite(mParam.mC5) && CAMath::Finite(mParam.mC7) && CAMath::Finite(mParam.mC9) && CAMath::Finite(mParam.mC10) && CAMath::Finite(mParam.mC12) && CAMath::Finite(mParam.mC14);
261
262 if (mParam.mC0 <= 0.f || mParam.mC2 <= 0.f || mParam.mC5 <= 0.f || mParam.mC9 <= 0.f || mParam.mC14 <= 0.f || mParam.mC0 > 5.f || mParam.mC2 > 5.f || mParam.mC5 > 2.f || mParam.mC9 > 2.f) {
263 ok = 0;
264 }
265
266 if (ok) {
267 ok = ok && (mParam.mC3 * mParam.mC3 <= mParam.mC5 * mParam.mC0) && (mParam.mC7 * mParam.mC7 <= mParam.mC9 * mParam.mC2) && (mParam.mC10 * mParam.mC10 <= mParam.mC14 * mParam.mC0) && (mParam.mC12 * mParam.mC12 <= mParam.mC14 * mParam.mC5);
268 }
269
270 return ok;
271}
272
273template <>
274GPUd() bool GPUTPCGMSectorTrack::TransportToX<2>(GPUTPCGMMerger* merger, float x, float Bz, GPUTPCGMBorderTrack& b, float maxSinPhi, bool doCov) const
275{
276 if (CAMath::Abs(x - mParam2.mX) < CAMath::Abs(x - mParam.mX) && mParam2.mX > 0) {
277 return TransportToX<1>(merger, x, Bz, b, maxSinPhi, doCov);
278 } else {
279 return TransportToX<0>(merger, x, Bz, b, maxSinPhi, doCov);
280 }
281}
282
283template <int I>
284GPUd() bool GPUTPCGMSectorTrack::TransportToX(GPUTPCGMMerger* merger, float x, float Bz, GPUTPCGMBorderTrack& b, float maxSinPhi, bool doCov) const
285{
286 const auto& param = I ? mParam2 : mParam;
287 Bz = -Bz;
288 float ex = param.mCosPhi;
289 float ey = param.mSinPhi;
290 float k = param.mQPt * Bz;
291 float dx = x - param.mX;
292 float ey1 = k * dx + ey;
293
294 if (CAMath::Abs(ey1) > maxSinPhi) {
295 return 0;
296 }
297
298 float ex1 = CAMath::Sqrt(1.f - ey1 * ey1);
299 float dxBz = dx * Bz;
300
301 float ss = ey + ey1;
302 float cc = ex + ex1;
303 float dxcci = dx / cc;
304 float norm2 = 1.f + ey * ey1 + ex * ex1;
305
306 float dy = dxcci * ss;
307
308 float dS;
309 {
310 float dl = dxcci * CAMath::Sqrt(norm2 + norm2);
311 float dSin = 0.5f * k * dl;
312 float a = dSin * dSin;
313 const float k2 = 1.f / 6.f;
314 const float k4 = 3.f / 40.f;
315 // const float k6 = 5.f/112.f;
316 dS = dl + dl * a * (k2 + a * (k4)); //+ k6*a) );
317 }
318
319 float dz = dS * param.mDzDs;
320
321 b.SetPar(0, param.mY + dy);
322 b.SetPar(1, param.mZ + dz);
323 b.SetPar(2, ey1);
324 b.SetPar(3, param.mDzDs);
325 b.SetPar(4, param.mQPt);
326 b.SetZOffsetLinear(merger->GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->convVertexTimeToZOffset(mSector, mTOffset, merger->Param().continuousMaxTimeBin));
327
328 if (!doCov) {
329 return (1);
330 }
331
332 float ex1i = 1.f / ex1;
333 float hh = dxcci * ex1i * norm2;
334 float h2 = hh * param.mSecPhi;
335 float h4 = Bz * dxcci * hh;
336
337 float c20 = param.mC3;
338 float c22 = param.mC5;
339 float c31 = param.mC7;
340 float c33 = param.mC9;
341 float c40 = param.mC10;
342 float c42 = param.mC12;
343 float c44 = param.mC14;
344
345 float c20ph4c42 = c20 + h4 * c42;
346 float h2c22 = h2 * c22;
347 float h4c44 = h4 * c44;
348 float n7 = c31 + dS * c33;
349
350 if (CAMath::Abs(param.mQPt) > 6.66f) // Special treatment for low Pt
351 {
352 b.SetCov(0, CAMath::Max(param.mC0, param.mC0 + h2 * h2c22 + h4 * h4c44 + 2.f * (h2 * c20ph4c42 + h4 * c40))); // Do not decrease Y cov for matching!
353 float C2tmp = dS * 2.f * c31;
354 if (C2tmp < 0) {
355 C2tmp = 0;
356 }
357 b.SetCov(1, param.mC2 + C2tmp + dS * dS * c33); // Incorrect formula, correct would be "dS * (c31 + n7)", but we need to make sure cov(Z) increases regardless of the direction of the propagation
358 } else {
359 b.SetCov(0, param.mC0 + h2 * h2c22 + h4 * h4c44 + 2.f * (h2 * c20ph4c42 + h4 * c40));
360 b.SetCov(1, param.mC2 + dS * (c31 + n7));
361 }
362 b.SetCov(2, c22 + dxBz * (c42 + c42 + dxBz * c44));
363 b.SetCov(3, c33);
364 b.SetCov(4, c44);
365 b.SetCovD(0, c20ph4c42 + h2c22 + dxBz * (c40 + h2 * c42 + h4c44));
366 b.SetCovD(1, n7);
367
368 b.LimitCov();
369
370 return 1;
371}
372
373template GPUdni() bool GPUTPCGMSectorTrack::TransportToX<0>(GPUTPCGMMerger* merger, float x, float Bz, GPUTPCGMBorderTrack& b, float maxSinPhi, bool doCov) const;
374template GPUdni() bool GPUTPCGMSectorTrack::TransportToX<1>(GPUTPCGMMerger* merger, float x, float Bz, GPUTPCGMBorderTrack& b, float maxSinPhi, bool doCov) const;
375
376GPUd() bool GPUTPCGMSectorTrack::TransportToXAlpha(GPUTPCGMMerger* merger, float newX, float sinAlpha, float cosAlpha, float Bz, GPUTPCGMBorderTrack& b, float maxSinPhi) const
377{
378 //*
379
380 float c00 = mParam.mC0;
381 float c11 = mParam.mC2;
382 float c20 = mParam.mC3;
383 float c22 = mParam.mC5;
384 float c31 = mParam.mC7;
385 float c33 = mParam.mC9;
386 float c40 = mParam.mC10;
387 float c42 = mParam.mC12;
388 float c44 = mParam.mC14;
389
390 float x, y;
391 float z = mParam.mZ;
392 float sinPhi = mParam.mSinPhi;
393 float cosPhi = mParam.mCosPhi;
394 float secPhi = mParam.mSecPhi;
395 float dzds = mParam.mDzDs;
396 float qpt = mParam.mQPt;
397
398 // Rotate the coordinate system in XY on the angle alpha
399 {
400 float sP = sinPhi, cP = cosPhi;
401 cosPhi = cP * cosAlpha + sP * sinAlpha;
402 sinPhi = -cP * sinAlpha + sP * cosAlpha;
403
404 if (CAMath::Abs(sinPhi) > GPUCA_MAX_SIN_PHI || CAMath::Abs(cP) < 1.e-2f) {
405 return 0;
406 }
407
408 secPhi = 1.f / cosPhi;
409 float j0 = cP * secPhi;
410 float j2 = cosPhi / cP;
411 x = mParam.mX * cosAlpha + mParam.mY * sinAlpha;
412 y = -mParam.mX * sinAlpha + mParam.mY * cosAlpha;
413
414 c00 *= j0 * j0;
415 c40 *= j0;
416
417 c22 *= j2 * j2;
418 c42 *= j2;
419 if (cosPhi < 0.f) { // rotate to 180'
420 cosPhi = -cosPhi;
421 secPhi = -secPhi;
422 sinPhi = -sinPhi;
423 dzds = -dzds;
424 qpt = -qpt;
425 c20 = -c20;
426 c31 = -c31;
427 c40 = -c40;
428 }
429 }
430
431 Bz = -Bz;
432 float ex = cosPhi;
433 float ey = sinPhi;
434 float k = qpt * Bz;
435 float dx = newX - x;
436 float ey1 = k * dx + ey;
437
438 if (CAMath::Abs(ey1) > maxSinPhi) {
439 return 0;
440 }
441
442 float ex1 = CAMath::Sqrt(1.f - ey1 * ey1);
443
444 float dxBz = dx * Bz;
445
446 float ss = ey + ey1;
447 float cc = ex + ex1;
448 float dxcci = dx / cc;
449 float norm2 = 1.f + ey * ey1 + ex * ex1;
450
451 float dy = dxcci * ss;
452
453 float dS;
454 {
455 float dl = dxcci * CAMath::Sqrt(norm2 + norm2);
456 float dSin = 0.5f * k * dl;
457 float a = dSin * dSin;
458 const float k2 = 1.f / 6.f;
459 const float k4 = 3.f / 40.f;
460 // const float k6 = 5.f/112.f;
461 dS = dl + dl * a * (k2 + a * (k4)); //+ k6*a) );
462 }
463
464 float ex1i = 1.f / ex1;
465 float dz = dS * dzds;
466
467 float hh = dxcci * ex1i * norm2;
468 float h2 = hh * secPhi;
469 float h4 = Bz * dxcci * hh;
470
471 float c20ph4c42 = c20 + h4 * c42;
472 float h2c22 = h2 * c22;
473 float h4c44 = h4 * c44;
474 float n7 = c31 + dS * c33;
475
476 b.SetPar(0, y + dy);
477 b.SetPar(1, z + dz);
478 b.SetPar(2, ey1);
479 b.SetPar(3, dzds);
480 b.SetPar(4, qpt);
481 b.SetZOffsetLinear(merger->GetConstantMem()->calibObjects.fastTransformHelper->getCorrMap()->convVertexTimeToZOffset(mSector, mTOffset, merger->Param().continuousMaxTimeBin));
482
483 b.SetCov(0, c00 + h2 * h2c22 + h4 * h4c44 + 2.f * (h2 * c20ph4c42 + h4 * c40));
484 b.SetCov(1, c11 + dS * (c31 + n7));
485 b.SetCov(2, c22 + dxBz * (c42 + c42 + dxBz * c44));
486 b.SetCov(3, c33);
487 b.SetCov(4, c44);
488 b.SetCovD(0, c20ph4c42 + h2c22 + dxBz * (c40 + h2 * c42 + h4c44));
489 b.SetCovD(1, n7);
490
491 b.LimitCov();
492
493 return 1;
494}
495
496GPUd() void GPUTPCGMSectorTrack::CopyBaseTrackCov()
497{
498 const float* GPUrestrict() cov = mOrigTrack -> Param().mC;
499 mParam.mC0 = cov[0];
500 mParam.mC2 = cov[2];
501 mParam.mC3 = cov[3];
502 mParam.mC5 = cov[5];
503 mParam.mC7 = cov[7];
504 mParam.mC9 = cov[9];
505 mParam.mC10 = cov[10];
506 mParam.mC12 = cov[12];
507 mParam.mC14 = cov[14];
508}
#define GPUdni()
#define GPUrestrict()
#define GPUCA_MAX_SIN_PHI
GLfloat GLfloat GLfloat alpha
Definition glcorearb.h:279
GLint GLenum GLint x
Definition glcorearb.h:403
GLuint index
Definition glcorearb.h:781
GLboolean GLboolean GLboolean b
Definition glcorearb.h:1233
GLint y
Definition glcorearb.h:270
typedef void(APIENTRYP PFNGLCULLFACEPROC)(GLenum mode)
GLenum GLfloat param
Definition glcorearb.h:271
GLboolean GLboolean GLboolean GLboolean a
Definition glcorearb.h:1233
GLdouble GLdouble GLdouble z
Definition glcorearb.h:843
const float k2
Definition MathUtils.h:72
Global TPC definitions and constants.
Definition SimTraits.h:167
GPUd() void PIDResponse
Definition PIDResponse.h:71
std::vector< o2::mch::ChannelCode > cc
std::vector< int > row