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