Project
Loading...
Searching...
No Matches
GPUTPCGMPropagator.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 "GPUTPCGMPropagator.h"
17#include "GPUParam.h"
19#include "GPUO2DataTypes.h"
20#include "GPUParam.inc"
21#include "GPUTPCGMMergerTypes.h"
22#include "GPUDebugStreamer.h"
23#include "GPUTPCGMMerger.h"
24
25#if defined(GPUCA_GM_USE_FULL_FIELD)
26#include "AliTracker.h"
27#include "AliMagF.h"
28#endif
29
30using namespace o2::gpu;
31
32GPUd() void GPUTPCGMPropagator::GetBxByBzBase(float cosAlpha, float sinAlpha, float X, float Y, float Z, float B[3]) const
33{
34 // get global coordinates
35
36 float gx = getGlobalX(cosAlpha, sinAlpha, X, Y);
37 float gy = getGlobalY(cosAlpha, sinAlpha, X, Y);
38
39#if defined(GPUCA_GM_USE_FULL_FIELD)
40 const float kCLight = gpu_common_constants::kCLight;
41 double r[3] = {gx, gy, Z};
42 double bb[3];
43 AliTracker::GetBxByBz(r, bb);
44 bb[0] *= kCLight;
45 bb[1] *= kCLight;
46 bb[2] *= kCLight;
47/*
48 cout<<"AliTracker::GetBz()= "<<AliTracker::GetBz()<<endl;
49 cout<<"AliTracker::UniformField() "<<AliTracker::UniformField()<<endl;
50 AliMagF* fld = (AliMagF*)TGeoGlobalMagField::Instance()->GetField();
51 cout<<"Fast field = "<<(void*) fld->GetFastField()<<endl;
52 AliMagF::BMap_t type = fld->GetMapType() ;
53 cout<<"Field type: "<<type<<endl;
54 // fMapType==k2BMap_t
55*/
56#else
57 float bb[3];
58 switch (mFieldRegion) {
59 case ITS:
60 mField->GetFieldIts(gx, gy, Z, bb);
61 break;
62 case TRD:
63 mField->GetFieldTrd(gx, gy, Z, bb);
64 break;
65 case TPC:
66 default:
67 mField->GetField(gx, gy, Z, bb);
68 }
69
70#endif
71
72 // rotate field to local coordinates
73
74 B[0] = bb[0] * cosAlpha + bb[1] * sinAlpha;
75 B[1] = -bb[0] * sinAlpha + bb[1] * cosAlpha;
76 B[2] = bb[2];
77}
78
79GPUd() float GPUTPCGMPropagator::GetBzBase(float cosAlpha, float sinAlpha, float X, float Y, float Z) const
80{
81 float gx = getGlobalX(cosAlpha, sinAlpha, X, Y);
82 float gy = getGlobalY(cosAlpha, sinAlpha, X, Y);
83
84#if defined(GPUCA_GM_USE_FULL_FIELD)
85 const float kCLight = gpu_common_constants::kCLight;
86 double r[3] = {gx, gy, Z};
87 double bb[3];
88 AliTracker::GetBxByBz(r, bb);
89 return bb[2] * kCLight;
90#else
91 switch (mFieldRegion) {
92 case ITS:
93 return mField->GetFieldItsBz(gx, gy, Z);
94 case TRD:
95 return mField->GetFieldTrdBz(gx, gy, Z);
96 case TPC:
97 default:
98 return mField->GetFieldBz(gx, gy, Z);
99 }
100#endif
101}
102
103GPUd() int32_t GPUTPCGMPropagator::RotateToAlpha(float newAlpha)
104{
105 //
106 // Rotate the track coordinate system in XY to the angle newAlpha
107 // return value is error code (0==no error)
108 //
109
110 float newCosAlpha, newSinAlpha;
111 CAMath::SinCos(newAlpha, newSinAlpha, newCosAlpha);
112
113 float cc = newCosAlpha * mCosAlpha + newSinAlpha * mSinAlpha; // CAMath::Cos(newAlpha - mAlpha);
114 float ss = newSinAlpha * mCosAlpha - newCosAlpha * mSinAlpha; // CAMath::Sin(newAlpha - mAlpha);
115
117
118 float x0 = mT0.X();
119 float y0 = mT0.Y();
120 float px0 = mT0.Px();
121 float py0 = mT0.Py();
122 // float pt0 = mT0.GetPt();
123
124 if (CAMath::Abs(mT->GetSinPhi()) >= mMaxSinPhi || CAMath::Abs(px0) < (1 - mMaxSinPhi)) {
125 return -1;
126 }
127
128 // rotate t0 track
129 float px1 = px0 * cc + py0 * ss;
130 float py1 = -px0 * ss + py0 * cc;
131
132 {
133 t0.X() = x0 * cc + y0 * ss;
134 t0.Y() = -x0 * ss + y0 * cc;
135 t0.Px() = px1;
136 t0.Py() = py1;
137 t0.UpdateValues();
138 }
139
140 if (CAMath::Abs(py1) > mMaxSinPhi * mT0.GetPt() || CAMath::Abs(px1) < (1 - mMaxSinPhi)) {
141 return -1;
142 }
143
144 // calculate X of rotated track:
145 float trackX = x0 * cc + ss * mT->Y();
146
147 // transport t0 to trackX
148 float dLp = 0;
149 float Bz;
150 if (mPropagateBzOnly) {
151 Bz = GetBzBase(newCosAlpha, newSinAlpha, t0.X(), t0.Y(), t0.Z());
152 if (t0.PropagateToXBzLight(trackX, Bz, dLp)) {
153 return -1;
154 }
155 } else {
156 float B[3];
157 GetBxByBz(newAlpha, t0.X(), t0.Y(), t0.Z(), B);
158 Bz = B[2];
159 if (t0.PropagateToXBxByBz(trackX, B[0], B[1], B[2], dLp)) {
160 return -1;
161 }
162 }
163
164 if (CAMath::Abs(t0.SinPhi()) >= mMaxSinPhi) {
165 return -1;
166 }
167
168 // now t0 is rotated and propagated, all checks are passed
169
170 // Rotate track using mT0 for linearisation. After rotation X is not fixed, but has a covariance
171
172 // Y Z Sin DzDs q/p
173 // Jacobian J0 = { { j0, 0, 0, 0, 0 }, // Y
174 // { 0, 1, 0, 0, 0 }, // Z
175 // { 0, 0, j1, 0, 0 }, // SinPhi
176 // { 0, 0, 0, 1, 0 }, // DzDs
177 // { 0, 0, 0, 0, 1 }, // q/p
178 // { j2, 0, 0, 0, 0 } }// X (rotated )
179
180 float j0 = cc;
181 float j1 = px1 / px0;
182 float j2 = ss;
183 // float dy = mT->Y() - y0;
184 // float ds = mT->SinPhi() - mT0.SinPhi();
185
186 mT->X() = trackX; // == x0*cc + ss*mT->Y() == t0.X() + j0*dy;
187 mT->Y() = -x0 * ss + cc * mT->Y(); //== t0.Y() + j0*dy;
188 // mT->SinPhi() = py1/pt0 + j1*ds; // == t0.SinPhi() + j1*ds; // use py1, since t0.SinPhi can have different sign
189 mT->SinPhi() = -CAMath::Sqrt(1.f - mT->SinPhi() * mT->SinPhi()) * ss + mT->SinPhi() * cc;
190
191 // Rotate cov. matrix Cr = J0 x C x J0T. Cr has one more row+column for X:
192 float* c = mT->Cov();
193
194 float c15 = c[0] * j0 * j2;
195 float c16 = c[1] * j2;
196 float c17 = c[3] * j1 * j2;
197 float c18 = c[6] * j2;
198 float c19 = c[10] * j2;
199 float c20 = c[0] * j2 * j2;
200
201 c[0] *= j0 * j0;
202 c[3] *= j0;
203 c[10] *= j0;
204
205 c[3] *= j1;
206 c[5] *= j1 * j1;
207 c[12] *= j1;
208
209 if (!mFitInProjections && mT->NDF() > 0) {
210 c[1] *= j0;
211 c[6] *= j0;
212 c[4] *= j1;
213 c[8] *= j1;
214 }
215
216 if (t0.SetDirectionAlongX()) { // change direction if Px < 0
217 mT->SinPhi() = -mT->SinPhi();
218 mT->DzDs() = -mT->DzDs();
219 mT->QPt() = -mT->QPt();
220 c[3] = -c[3]; // covariances with SinPhi
221 c[4] = -c[4];
222 c17 = -c17;
223 c[6] = -c[6]; // covariances with DzDs
224 c[7] = -c[7];
225 c18 = -c18;
226 c[10] = -c[10]; // covariances with QPt
227 c[11] = -c[11];
228 c19 = -c19;
229 }
230
231 // Now fix the X coordinate: so to say, transport track T to fixed X = mT->X().
232 // only covariance changes. Use rotated and transported t0 for linearisation
233 float j3 = -t0.Py() / t0.Px();
234 float j4 = -t0.Pz() / t0.Px();
235 float j5 = t0.QPt() * Bz;
236
237 // Y Z Sin DzDs q/p X
238 // Jacobian J1 = { { 1, 0, 0, 0, 0, j3 }, // Y
239 // { 0, 1, 0, 0, 0, j4 }, // Z
240 // { 0, 0, 1, 0, 0, j5 }, // SinPhi
241 // { 0, 0, 0, 1, 0, 0 }, // DzDs
242 // { 0, 0, 0, 0, 1, 0 } }; // q/p
243
244 float h15 = c15 + c20 * j3;
245 float h16 = c16 + c20 * j4;
246 float h17 = c17 + c20 * j5;
247
248 c[0] += j3 * (c15 + h15);
249
250 c[2] += j4 * (c16 + h16);
251
252 c[3] += c17 * j3 + h15 * j5;
253 c[5] += j5 * (c17 + h17);
254
255 c[7] += c18 * j4;
256 // c[ 9] = c[ 9];
257
258 c[10] += c19 * j3;
259 c[12] += c19 * j5;
260 // c[14] = c[14];
261
262 if (!mFitInProjections && mT->NDF() > 0) {
263 c[1] += c16 * j3 + h15 * j4;
264 c[4] += c17 * j4 + h16 * j5;
265 c[6] += c18 * j3;
266 c[8] += c18 * j5;
267 c[11] += c19 * j4;
268 // c[13] = c[13];
269 }
270
271 mAlpha = newAlpha;
272 mCosAlpha = newCosAlpha;
273 mSinAlpha = newSinAlpha;
274 mT0 = t0;
275
276 return 0;
277}
278
279GPUd() int32_t GPUTPCGMPropagator::PropagateToXAlpha(float posX, float posAlpha, bool inFlyDirection)
280{
281 if (mPropagateBzOnly) {
282 return PropagateToXAlphaBz(posX, posAlpha, inFlyDirection);
283 }
284 if (CAMath::Abs(posAlpha - mAlpha) > 1.e-4f) {
285 if (RotateToAlpha(posAlpha) != 0) {
286 return -2;
287 }
288 }
289 if (CAMath::Abs(posX - mT->X()) < 1.e-7f) {
290 mT->SetX(posX);
291 return 0;
292 }
293
294 float B[3];
295 GetBxByBz(mT0.X(), mT0.Y(), mT0.Z(), B);
296
297 // propagate mT0 to t0e
298
300 float dLp = 0;
301 if (t0e.PropagateToXBxByBz(posX, B[0], B[1], B[2], dLp) && t0e.PropagateToXBzLight(posX, B[2], dLp)) {
302 return 1;
303 }
304
305 if (CAMath::Abs(t0e.SinPhi()) >= mMaxSinPhi) {
306 return -3;
307 }
308
309 return FollowLinearization(t0e, B[2], dLp, inFlyDirection);
310}
311
312GPUd() int32_t GPUTPCGMPropagator::PropagateToXAlphaBz(float posX, float posAlpha, bool inFlyDirection)
313{
314 if (CAMath::Abs(posAlpha - mAlpha) > 1.e-4f) {
315 if (RotateToAlpha(posAlpha) != 0) {
316 return -2;
317 }
318 }
319 if (CAMath::Abs(posX - mT->X()) < 1.e-7f) {
320 mT->SetX(posX);
321 return 0;
322 }
323
324 float Bz = GetBz(mT0.X(), mT0.Y(), mT0.Z());
325
326 // propagate mT0 to t0e
327
329 float dLp = 0;
330 if (t0e.PropagateToXBzLight(posX, Bz, dLp)) {
331 return 1;
332 }
333
334 if (CAMath::Abs(t0e.SinPhi()) >= mMaxSinPhi) {
335 return -3;
336 }
337
338 return FollowLinearization(t0e, Bz, dLp, inFlyDirection);
339}
340
341GPUd() int32_t GPUTPCGMPropagator::FollowLinearization(const GPUTPCGMPhysicalTrackModel& GPUrestrict() t0e, float Bz, float dLp, bool inFlyDirection)
342{
343 // t0e is alrerady extrapolated t0
344
345 // propagate track and cov matrix with derivatives for (0,0,Bz) field
346
347 float dS = dLp * t0e.GetPt();
348 float dL = CAMath::Abs(dLp * t0e.GetP());
349
350 if (inFlyDirection) {
351 dL = -dL;
352 }
353
354 float ey = mT0.SinPhi();
355 float ex = mT0.CosPhi();
356 float exi = mT0.SecPhi();
357 float ey1 = t0e.GetSinPhi();
358 float ex1 = t0e.GetCosPhi();
359 float ex1i = t0e.GetSecPhi();
360
361 float k = -mT0.QPt() * Bz;
362 float dx = t0e.GetX() - mT0.X();
363 float kdx = k * dx;
364 float cc = ex + ex1;
365 float cci = 1.f / cc;
366
367 float dxcci = dx * cci;
368 float hh = dxcci * ex1i * (1.f + ex * ex1 + ey * ey1);
369
370 float j02 = exi * hh;
371 float j04 = -Bz * dxcci * hh;
372 float j13 = dS;
373 float j24 = -dx * Bz;
374
375 float* p = mT->Par();
376
377 float d0 = p[0] - mT0.Y();
378 float d1 = p[1] - mT0.Z();
379 float d2 = p[2] - mT0.SinPhi();
380 float d3 = p[3] - mT0.DzDs();
381 float d4 = p[4] - mT0.QPt();
382
383 float newSinPhi = ey1 + d2 + j24 * d4;
384 if (mT->NDF() >= 15 && CAMath::Abs(newSinPhi) > GPUCA_MAX_SIN_PHI) {
385 return -4;
386 }
387
388 if (mMatLUT) {
389 UpdateMaterial(t0e);
390 }
391
392 mT0 = t0e;
393 mT->X() = t0e.GetX();
394 p[0] = t0e.GetY() + d0 + j02 * d2 + j04 * d4;
395 p[1] = t0e.GetZ() + d1 + j13 * d3;
396 p[2] = newSinPhi;
397 p[3] = t0e.GetDzDs() + d3;
398 p[4] = t0e.GetQPt() + d4;
399
400 float* c = mT->Cov();
401
402 float c20 = c[3];
403 float c21 = c[4];
404 float c22 = c[5];
405
406 float c30 = c[6];
407 float c31 = c[7];
408 float c32 = c[8];
409 float c33 = c[9];
410
411 float c40 = c[10];
412 float c41 = c[11];
413 float c42 = c[12];
414 float c43 = c[13];
415 float c44 = c[14];
416
417 if (mFitInProjections || mT->NDF() <= 0) {
418 float c20ph04c42 = c20 + j04 * c42;
419 float j02c22 = j02 * c22;
420 float j04c44 = j04 * c44;
421
422 float n6 = c30 + j02 * c32 + j04 * c43;
423 float n7 = c31 + j13 * c33;
424 float n10 = c40 + j02 * c42 + j04c44;
425 float n11 = c41 + j13 * c43;
426 float n12 = c42 + j24 * c44;
427
428 c[0] += j02 * j02c22 + j04 * j04c44 + 2.f * (j02 * c20ph04c42 + j04 * c40);
429 c[1] += j02 * c21 + j04 * c41 + j13 * n6;
430 c[2] += j13 * (c31 + n7);
431 c[3] = c20ph04c42 + j02c22 + j24 * n10;
432 c[4] = c21 + j13 * c32 + j24 * n11;
433 c[5] = c22 + j24 * (c42 + n12);
434 c[6] = n6;
435 c[7] = n7;
436 c[8] = c32 + c43 * j24;
437 c[10] = n10;
438 c[11] = n11;
439 c[12] = n12;
440 } else {
441 float c00 = c[0];
442 float c10 = c[1];
443 float c11 = c[2];
444
445 float ss = ey + ey1;
446 float tg = ss * cci;
447 float xx = 1.f - 0.25f * kdx * kdx * (1.f + tg * tg);
448 if (xx < 1.e-8f) {
449 return -1;
450 }
451 xx = CAMath::Sqrt(xx);
452 float yy = CAMath::Sqrt(ss * ss + cc * cc);
453
454 float j12 = dx * mT0.DzDs() * tg * (2.f + tg * (ey * exi + ey1 * ex1i)) / (xx * yy);
455 float j14 = 0;
456 if (CAMath::Abs(mT0.QPt()) > 1.e-6f) {
457 j14 = (2.f * xx * ex1i * dx / yy - dS) * mT0.DzDs() / mT0.QPt();
458 } else {
459 j14 = -mT0.DzDs() * Bz * dx * dx * exi * exi * exi * (0.5f * ey + (1.f / 3.f) * kdx * (1 + 2.f * ey * ey) * exi * exi);
460 }
461
462 p[1] += j12 * d2 + j14 * d4;
463
464 float h00 = c00 + c20 * j02 + c40 * j04;
465 // float h01 = c10 + c21*j02 + c41*j04;
466 float h02 = c20 + c22 * j02 + c42 * j04;
467 // float h03 = c30 + c32*j02 + c43*j04;
468 float h04 = c40 + c42 * j02 + c44 * j04;
469
470 float h10 = c10 + c20 * j12 + c30 * j13 + c40 * j14;
471 float h11 = c11 + c21 * j12 + c31 * j13 + c41 * j14;
472 float h12 = c21 + c22 * j12 + c32 * j13 + c42 * j14;
473 float h13 = c31 + c32 * j12 + c33 * j13 + c43 * j14;
474 float h14 = c41 + c42 * j12 + c43 * j13 + c44 * j14;
475
476 float h20 = c20 + c40 * j24;
477 float h21 = c21 + c41 * j24;
478 float h22 = c22 + c42 * j24;
479 float h23 = c32 + c43 * j24;
480 float h24 = c42 + c44 * j24;
481
482 c[0] = h00 + h02 * j02 + h04 * j04;
483
484 c[1] = h10 + h12 * j02 + h14 * j04;
485 c[2] = h11 + h12 * j12 + h13 * j13 + h14 * j14;
486
487 c[3] = h20 + h22 * j02 + h24 * j04;
488 c[4] = h21 + h22 * j12 + h23 * j13 + h24 * j14;
489 c[5] = h22 + h24 * j24;
490
491 c[6] = c30 + c32 * j02 + c43 * j04;
492 c[7] = c31 + c32 * j12 + c33 * j13 + c43 * j14;
493 c[8] = c32 + c43 * j24;
494 // c[ 9] = c33;
495
496 c[10] = c40 + c42 * j02 + c44 * j04;
497 c[11] = c41 + c42 * j12 + c43 * j13 + c44 * j14;
498 c[12] = c42 + c44 * j24;
499 // c[13] = c43;
500 // c[14] = c44;
501 }
502
503 float& mC22 = c[5];
504 float& mC33 = c[9];
505 float& mC40 = c[10];
506 float& mC41 = c[11];
507 float& mC42 = c[12];
508 float& mC43 = c[13];
509 float& mC44 = c[14];
510
511 float dLmask = 0.f;
512 bool maskMS = (CAMath::Abs(dL) < mMaterial.DLMax);
513 if (maskMS) {
514 dLmask = dL;
515 }
516 float dLabs = CAMath::Abs(dLmask);
517
518 // Energy Loss
519 if (true) {
520 // std::cout<<"APPLY ENERGY LOSS!!!"<<std::endl;
521 float corr = 1.f - mMaterial.EP2 * dLmask;
522 float corrInv = 1.f / corr;
523 mT0.Px() *= corrInv;
524 mT0.Py() *= corrInv;
525 mT0.Pz() *= corrInv;
526 mT0.Pt() *= corrInv;
527 mT0.P() *= corrInv;
528 mT0.QPt() *= corr;
529
530 p[4] *= corr;
531
532 mC40 *= corr;
533 mC41 *= corr;
534 mC42 *= corr;
535 mC43 *= corr;
536 mC44 = mC44 * corr * corr + dLabs * mMaterial.sigmadE2;
537 } else {
538 // std::cout<<"DONT APPLY ENERGY LOSS!!!"<<std::endl;
539 }
540
541 // Multiple Scattering
542 if (true) {
543 mC22 += dLabs * mMaterial.k22 * mT0.CosPhi() * mT0.CosPhi();
544 mC33 += dLabs * mMaterial.k33;
545 mC43 += dLabs * mMaterial.k43;
546 mC44 += dLabs * mMaterial.k44;
547 }
548 return 0;
549}
550
551GPUd() int32_t GPUTPCGMPropagator::GetPropagatedYZ(float x, float& GPUrestrict() projY, float& GPUrestrict() projZ)
552{
553 float bz = GetBz(mT->X(), mT->Y(), mT->Z());
554 float k = mT0.QPt() * bz;
555 float dx = x - mT->X();
556 float ex = mT0.CosPhi();
557 float ey = mT0.SinPhi();
558 float ey1 = ey - k * dx;
559 if (CAMath::Abs(ey1) > GPUCA_MAX_SIN_PHI) {
560 return 1;
561 }
562 float ss = ey + ey1;
563 float ex1 = CAMath::Sqrt(1.f - ey1 * ey1);
564 if (ex < 0) {
565 ex1 = -ex1;
566 }
567 float cc = ex + ex1;
568 float dxcci = dx / cc;
569 float dy = dxcci * ss;
570 float norm2 = 1.f + ey * ey1 + ex * ex1;
571 float dl = dxcci * CAMath::Sqrt(norm2 + norm2);
572 float dS;
573 {
574 float dSin = 0.5f * k * dl;
575 float a = dSin * dSin;
576 const float k2 = 1.f / 6.f;
577 const float k4 = 3.f / 40.f;
578 dS = dl + dl * a * (k2 + a * (k4));
579 }
580 float dz = dS * mT0.DzDs();
581 projY = mT->Y() + dy;
582 projZ = mT->Z() + dz;
583 return 0;
584}
585
586GPUd() void GPUTPCGMPropagator::GetErr2(float& GPUrestrict() err2Y, float& GPUrestrict() err2Z, const GPUParam& GPUrestrict() param, float posZ, int32_t iRow, int16_t clusterState, int8_t sector, float time, float avgCharge, float charge) const
587{
588 GetErr2(err2Y, err2Z, param, mT0.GetSinPhi(), mT0.DzDs(), posZ, mT->GetX(), mT->GetY(), iRow, clusterState, sector, time, avgCharge, charge, mSeedingErrors);
589}
590
591GPUd() void GPUTPCGMPropagator::GetErr2(float& GPUrestrict() err2Y, float& GPUrestrict() err2Z, const GPUParam& GPUrestrict() param, float snp, float tgl, float posZ, float trackX, float trackY, int32_t iRow, int16_t clusterState, int8_t sector, float time, float avgCharge, float charge, bool seedingErrors)
592{
593#ifndef GPUCA_TPC_GEOMETRY_O2
594 if (seedingErrors) {
595 param.GetClusterErrorsSeeding2(sector, iRow, posZ, snp, tgl, time, err2Y, err2Z);
596 } else
597#endif
598 {
599 param.GetClusterErrors2(sector, iRow, posZ, snp, tgl, time, avgCharge, charge, err2Y, err2Z);
600 }
601 param.UpdateClusterError2ByState(clusterState, err2Y, err2Z);
602 float statErr2 = param.GetSystematicClusterErrorIFC2(trackX, trackY, posZ, sector >= (GPUCA_NSECTORS / 2));
603 if (sector >= GPUCA_NSECTORS / 2 + 1 && sector <= GPUCA_NSECTORS / 2 + 2) {
604 statErr2 += param.GetSystematicClusterErrorC122(trackX, trackY, sector);
605 }
606 err2Y += statErr2;
607 err2Z += statErr2;
608}
609
610GPUd() float GPUTPCGMPropagator::PredictChi2(float posY, float posZ, int32_t iRow, const GPUParam& GPUrestrict() param, int16_t clusterState, int8_t sector, float time, float avgCharge, float charge) const
611{
612 float err2Y, err2Z;
613 GetErr2(err2Y, err2Z, param, posZ, iRow, clusterState, sector, time, avgCharge, charge);
614 return PredictChi2(posY, posZ, err2Y, err2Z);
615}
616
617GPUd() float GPUTPCGMPropagator::PredictChi2(float posY, float posZ, float err2Y, float err2Z) const
618{
619 const float* mC = mT->Cov();
620 const float* mP = mT->Par();
621 const float z0 = posY - mP[0];
622 const float z1 = posZ - mP[1];
623
624 if (!mFitInProjections || mT->NDF() <= 0) {
625 const float w0 = 1.f / (err2Y + mC[0]);
626 const float w2 = 1.f / (err2Z + mC[2]);
627 return w0 * z0 * z0 + w2 * z1 * z1;
628 } else {
629 float w0 = mC[2] + err2Z, w1 = mC[1], w2 = mC[0] + err2Y;
630 { // Invert symmetric matrix
631 float det = w0 * w2 - w1 * w1;
632 if (CAMath::Abs(det) < 1.e-10f) {
633 det = 1.e-10f;
634 }
635 det = 1.f / det;
636 w0 = w0 * det;
637 w1 = -w1 * det;
638 w2 = w2 * det;
639 }
640 return CAMath::Abs((w0 * z0 + w1 * z1) * z0) + CAMath::Abs((w1 * z0 + w2 * z1) * z1);
641 }
642}
643
644GPUd() int32_t GPUTPCGMPropagator::Update(float posY, float posZ, int32_t iRow, const GPUParam& GPUrestrict() param, int16_t clusterState, int8_t rejectChi2, gputpcgmmergertypes::InterpolationErrorHit* inter, bool refit, int8_t sector, float time, float avgInvCharge, float invCharge GPUCA_DEBUG_STREAMER_CHECK(, DebugStreamerVals* debugVals))
645{
646 float err2Y, err2Z;
647 GetErr2(err2Y, err2Z, param, posZ, iRow, clusterState, sector, time, avgInvCharge, invCharge);
648 GPUCA_DEBUG_STREAMER_CHECK(if (debugVals) { debugVals->err2Y = err2Y; debugVals->err2Z = err2Z; });
649
650 if (rejectChi2 >= rejectInterFill) {
651 if (rejectChi2 == rejectInterReject && inter->errorY < (GPUCA_PAR_MERGER_INTERPOLATION_ERROR_TYPE_A)0) {
652 rejectChi2 = rejectDirect;
653 } else {
654 int32_t retVal = InterpolateReject(param, posY, posZ, clusterState, rejectChi2, inter, err2Y, err2Z);
655 GPUCA_DEBUG_STREAMER_CHECK(if (debugVals) { debugVals->retVal = retVal; });
656 if (retVal) {
657 return retVal;
658 }
659 }
660 }
661
662 if (mT->NDF() == -5) { // first measurement: no need to filter, as the result is known in advance. just set it.
663 mT->ResetCovariance();
664 float* mC = mT->Cov();
665 float* mP = mT->Par();
666 if (refit) {
667 mC[14] = CAMath::Max(0.5f, CAMath::Abs(mP[4]));
668 mC[5] = CAMath::Max(0.2f, CAMath::Abs(mP[2]) / 2);
669 mC[9] = CAMath::Max(0.5f, CAMath::Abs(mP[3]) / 2);
670 }
671 mP[0] = posY;
672 mP[1] = posZ;
673 mC[0] = err2Y;
674 mC[2] = err2Z;
675 mT->NDF() = -3;
676 return 0;
677 }
678
679 return Update(posY, posZ, clusterState, rejectChi2 == rejectDirect || (param.rec.tpc.mergerInterpolateRejectAlsoOnCurrentPosition && rejectChi2 == rejectInterReject), err2Y, err2Z, &param);
680}
681
682GPUd() int32_t GPUTPCGMPropagator::InterpolateReject(const GPUParam& GPUrestrict() param, float posY, float posZ, int16_t clusterState, int8_t rejectChi2, gputpcgmmergertypes::InterpolationErrorHit* inter, float err2Y, float err2Z)
683{
684 float* GPUrestrict() mC = mT->Cov();
685 float* GPUrestrict() mP = mT->Par();
686 if (rejectChi2 == rejectInterFill) {
687 inter->posY = mP[0];
688 inter->posZ = mP[1];
689 inter->errorY = mC[0];
690 inter->errorZ = mC[2];
691 } else if (rejectChi2 == rejectInterReject) {
692 float chi2Y, chi2Z;
693 if (mFitInProjections || mT->NDF() <= 0) {
694 const float Iz0 = inter->posY - mP[0];
695 const float Iz1 = inter->posZ - mP[1];
696 const float Iw0 = 1.f / (mC[0] + (float)inter->errorY);
697 const float Iw2 = 1.f / (mC[2] + (float)inter->errorZ);
698 const float Ik00 = mC[0] * Iw0;
699 const float Ik11 = mC[2] * Iw2;
700 const float ImP0 = mP[0] + Ik00 * Iz0;
701 const float ImP1 = mP[1] + Ik11 * Iz1;
702 const float ImC0 = mC[0] - Ik00 * mC[0];
703 const float ImC2 = mC[2] - Ik11 * mC[2];
704 // printf("\t%21sInterpo ----- abde artaf%16s Y %8.3f, Z %8.3f (Errors %f <-- (%f, %f) %f <-- (%f, %f))\n", "", "", ImP0, ImP1, sqrtf(ImC0), sqrtf(mC[0]), sqrtf(inter->errorY), sqrtf(ImC2), sqrtf(mC[2]), sqrtf(inter->errorZ));
705 const float Jz0 = posY - ImP0;
706 const float Jz1 = posZ - ImP1;
707 const float Jw0 = 1.f / (ImC0 + err2Y);
708 const float Jw2 = 1.f / (ImC2 + err2Z);
709 chi2Y = Jw0 * Jz0 * Jz0;
710 chi2Z = Jw2 * Jz1 * Jz1;
711 } else {
712 const float Iz0 = inter->posY - mP[0];
713 const float Iz1 = inter->posZ - mP[1];
714 float Iw0 = mC[2] + (float)inter->errorZ;
715 float Iw2 = mC[0] + (float)inter->errorY;
716 float Idet = CAMath::Max(1e-10f, Iw0 * Iw2 - mC[1] * mC[1]);
717 Idet = 1.f / Idet;
718 Iw0 *= Idet;
719 const float Iw1 = mC[1] * Idet;
720 Iw2 *= Idet;
721 const float Ik00 = mC[0] * Iw0 + mC[1] * Iw1;
722 const float Ik01 = mC[0] * Iw1 + mC[1] * Iw2;
723 const float Ik10 = mC[1] * Iw0 + mC[2] * Iw1;
724 const float Ik11 = mC[1] * Iw1 + mC[2] * Iw2;
725 const float ImP0 = mP[0] + Ik00 * Iz0 + Ik01 * Iz1;
726 const float ImP1 = mP[1] + Ik10 * Iz0 + Ik11 * Iz1;
727 const float ImC0 = mC[0] - Ik00 * mC[0] + Ik01 * mC[1];
728 const float ImC1 = mC[1] - Ik10 * mC[0] + Ik11 * mC[1];
729 const float ImC2 = mC[2] - Ik10 * mC[1] + Ik11 * mC[2];
730 const float Jz0 = posY - ImP0;
731 const float Jz1 = posZ - ImP1;
732 float Jw0 = ImC2 + err2Z;
733 float Jw2 = ImC0 + err2Y;
734 float Jdet = CAMath::Max(1e-10f, Jw0 * Jw2 - ImC1 * ImC1);
735 Jdet = 1.f / Jdet;
736 Jw0 *= Jdet;
737 const float Jw1 = ImC1 * Jdet;
738 Jw2 *= Jdet;
739 chi2Y = CAMath::Abs((Jw0 * Jz0 + Jw1 * Jz1) * Jz0);
740 chi2Z = CAMath::Abs((Jw1 * Jz0 + Jw2 * Jz1) * Jz1);
741 }
742 if (RejectCluster(chi2Y * param.rec.tpc.clusterRejectChi2TolleranceY, chi2Z * param.rec.tpc.clusterRejectChi2TolleranceZ, clusterState)) { // TODO: Relative Pt resolution decreases slightly, why?
743 return updateErrorClusterRejectedInInterpolation;
744 }
745 }
746 return 0;
747}
748
749GPUd() int32_t GPUTPCGMPropagator::Update(float posY, float posZ, int16_t clusterState, bool rejectChi2, float err2Y, float err2Z, const GPUParam* GPUrestrict() param)
750{
751 float* GPUrestrict() mC = mT->Cov();
752 float* GPUrestrict() mP = mT->Par();
753
754 const float d00 = mC[0], d01 = mC[1], d02 = mC[3], d03 = mC[6], d04 = mC[10];
755 const float d10 = mC[1], d11 = mC[2], d12 = mC[4], d13 = mC[7], d14 = mC[11];
756
757 const float z0 = posY - mP[0];
758 const float z1 = posZ - mP[1];
759 float w0, w1, w2, chi2Y, chi2Z;
760 if (mFitInProjections || mT->NDF() <= 0) {
761 w0 = 1.f / (err2Y + d00);
762 w1 = 0;
763 w2 = 1.f / (err2Z + d11);
764 chi2Y = w0 * z0 * z0;
765 chi2Z = w2 * z1 * z1;
766 } else {
767 w0 = d11 + err2Z, w1 = d10, w2 = d00 + err2Y;
768 { // Invert symmetric matrix
769 float det = w0 * w2 - w1 * w1;
770 if (CAMath::Abs(det) < 1.e-10f) {
771 return updateErrorFitFailed;
772 }
773 det = 1.f / det;
774 w0 = w0 * det;
775 w1 = -w1 * det;
776 w2 = w2 * det;
777 }
778 chi2Y = CAMath::Abs((w0 * z0 + w1 * z1) * z0);
779 chi2Z = CAMath::Abs((w1 * z0 + w2 * z1) * z1);
780 }
781 float dChi2 = chi2Y + chi2Z;
782 // GPUInfo("hits %d chi2 %f, new %f %f (dy %f dz %f)", N, mChi2, chi2Y, chi2Z, z0, z1);
783 if (rejectChi2 && RejectCluster(chi2Y * param->rec.tpc.clusterRejectChi2TolleranceY, chi2Z * param->rec.tpc.clusterRejectChi2TolleranceZ, clusterState)) {
784 return updateErrorClusterRejectedInUpdate;
785 }
786 mT->Chi2() += dChi2;
787 mT->NDF() += 2;
788
789 if (mFitInProjections || mT->NDF() <= 0) {
790 const float k00 = d00 * w0;
791 const float k20 = d02 * w0;
792 const float k40 = d04 * w0;
793 const float k11 = d11 * w2;
794 const float k31 = d13 * w2;
795 mP[0] += k00 * z0;
796 mP[1] += k11 * z1;
797 mP[2] += k20 * z0;
798 mP[3] += k31 * z1;
799 mP[4] += k40 * z0;
800
801 mC[0] -= k00 * d00;
802 mC[2] -= k11 * d11;
803 mC[3] -= k20 * d00;
804 mC[5] -= k20 * d02;
805 mC[7] -= k31 * d11;
806 mC[9] -= k31 * d13;
807 mC[10] -= k00 * d04;
808 mC[12] -= k40 * d02;
809 mC[14] -= k40 * d04;
810 } else {
811 const float k00 = d00 * w0 + d01 * w1;
812 const float k01 = d00 * w1 + d10 * w2;
813 const float k10 = d01 * w0 + d11 * w1;
814 const float k11 = d01 * w1 + d11 * w2;
815 const float k20 = d02 * w0 + d12 * w1;
816 const float k21 = d02 * w1 + d12 * w2;
817 const float k30 = d03 * w0 + d13 * w1;
818 const float k31 = d03 * w1 + d13 * w2;
819 const float k40 = d04 * w0 + d14 * w1;
820 const float k41 = d04 * w1 + d14 * w2;
821
822 mP[0] += k00 * z0 + k01 * z1;
823 mP[1] += k10 * z0 + k11 * z1;
824 mP[2] += k20 * z0 + k21 * z1;
825 mP[3] += k30 * z0 + k31 * z1;
826 mP[4] += k40 * z0 + k41 * z1;
827
828 mC[0] -= k00 * d00 + k01 * d10;
829
830 mC[2] -= k10 * d01 + k11 * d11;
831
832 mC[3] -= k20 * d00 + k21 * d10;
833 mC[5] -= k20 * d02 + k21 * d12;
834
835 mC[7] -= k30 * d01 + k31 * d11;
836 mC[9] -= k30 * d03 + k31 * d13;
837
838 mC[10] -= k40 * d00 + k41 * d10;
839 mC[12] -= k40 * d02 + k41 * d12;
840 mC[14] -= k40 * d04 + k41 * d14;
841
842 mC[1] -= k10 * d00 + k11 * d10;
843 mC[4] -= k20 * d01 + k21 * d11;
844 mC[6] -= k30 * d00 + k31 * d10;
845 mC[8] -= k30 * d02 + k31 * d12;
846 mC[11] -= k40 * d01 + k41 * d11;
847 mC[13] -= k40 * d03 + k41 * d13;
848 }
849 return 0;
850}
851
852//*
853//* Multiple scattering and energy losses
854//*
855
856GPUd() float GPUTPCGMPropagator::ApproximateBetheBloch(float beta2)
857{
858 //------------------------------------------------------------------
859 // This is an approximation of the Bethe-Bloch formula with
860 // the density effect taken into account at beta*gamma > 3.5
861 // (the approximation is reasonable only for solid materials)
862 //------------------------------------------------------------------
863
864 const float log0 = CAMath::Log(5940.f);
865 const float log1 = CAMath::Log(3.5f * 5940.f);
866
867 bool bad = (beta2 >= .999f) || (beta2 < 1.e-8f);
868
869 if (bad) {
870 beta2 = 0.5f;
871 }
872
873 float a = beta2 / (1.f - beta2);
874 float b = 0.5f * CAMath::Log(a);
875 float d = 0.153e-3f / beta2;
876 float c = b - beta2;
877
878 float ret = d * (log0 + b + c);
879 float case1 = d * (log1 + c);
880
881 if (a > 3.5f * 3.5f) {
882 ret = case1;
883 }
884 if (bad) {
885 ret = 0.f;
886 }
887
888 return ret;
889}
890
891GPUd() void GPUTPCGMPropagator::CalculateMaterialCorrection()
892{
893 //*!
894
895 const float mass = 0.13957f;
896
897 float qpt = mT0.GetQPt();
898 if (CAMath::Abs(qpt) > 20) {
899 qpt = 20;
900 }
901
902 float w2 = (1.f + mT0.GetDzDs() * mT0.GetDzDs()); //==(P/pt)2
903 float pti2 = qpt * qpt;
904 if (pti2 < 1.e-4f) {
905 pti2 = 1.e-4f;
906 }
907
908 float mass2 = mass * mass;
909 float beta2 = w2 / (w2 + mass2 * pti2);
910
911 float p2 = w2 / pti2; // impuls 2
912 float betheRho = ApproximateBetheBloch(p2 / mass2) * mMaterial.rho;
913 float E = CAMath::Sqrt(p2 + mass2);
914 float theta2 = (14.1f * 14.1f / 1.e6f) / (beta2 * p2) * mMaterial.radLenInv;
915
916 mMaterial.EP2 = E / p2;
917
918 // Approximate energy loss fluctuation (M.Ivanov)
919
920 const float knst = 0.0007f; // To be tuned.
921 mMaterial.sigmadE2 = knst * mMaterial.EP2 * qpt;
922 mMaterial.sigmadE2 = mMaterial.sigmadE2 * mMaterial.sigmadE2;
923
924 mMaterial.k22 = theta2 * w2;
925 mMaterial.k33 = mMaterial.k22 * w2;
926 mMaterial.k43 = 0.f;
927 mMaterial.k44 = theta2 * mT0.GetDzDs() * mT0.GetDzDs() * pti2;
928
929 float br = (betheRho > 1.e-8f) ? betheRho : 1.e-8f;
930 mMaterial.DLMax = 0.3f * E / br;
931 mMaterial.EP2 *= betheRho;
932 mMaterial.sigmadE2 = mMaterial.sigmadE2 * betheRho; // + mMaterial.fK44;
933}
934
935GPUd() void GPUTPCGMPropagator::Rotate180()
936{
937 mT0.X() = -mT0.X();
938 mT0.Y() = -mT0.Y();
939 mT0.Q() = -mT0.Q();
940 mT0.Pz() = -mT0.Pz();
941 mT0.UpdateValues();
942
943 mT->X() = -mT->X();
944 mT->Y() = -mT->Y();
945 mT->QPt() = -mT->QPt();
946 mT->DzDs() = -mT->DzDs();
947
948 mAlpha = mAlpha + CAMath::Pi();
949 while (mAlpha >= CAMath::Pi()) {
950 mAlpha -= CAMath::TwoPi();
951 }
952 while (mAlpha < -CAMath::Pi()) {
953 mAlpha += CAMath::TwoPi();
954 }
955 mCosAlpha = -mCosAlpha;
956 mSinAlpha = -mSinAlpha;
957
958 float* c = mT->Cov();
959 c[6] = -c[6];
960 c[7] = -c[7];
961 c[8] = -c[8];
962 c[10] = -c[10];
963 c[11] = -c[11];
964 c[12] = -c[12];
965}
966
967GPUd() void GPUTPCGMPropagator::ChangeDirection()
968{
969 mT0.Py() = -mT0.Py();
970 mT0.Pz() = -mT0.Pz();
971 mT0.Q() = -mT0.Q();
972 mT->SinPhi() = -mT->SinPhi();
973 mT->DzDs() = -mT->DzDs();
974 mT->QPt() = -mT->QPt();
975 mT0.UpdateValues();
976
977 float* c = mT->Cov();
978 c[3] = -c[3];
979 c[4] = -c[4];
980 c[6] = -c[6];
981 c[7] = -c[7];
982 c[10] = -c[10];
983 c[11] = -c[11];
984}
985
986GPUd() void GPUTPCGMPropagator::Mirror(bool inFlyDirection)
987{
988 // mirror the track and the track approximation to the point which has the same X, but located on the other side of trajectory
989 float B[3];
990 GetBxByBz(mT0.X(), mT0.Y(), mT0.Z(), B);
991 float Bz = B[2];
992 if (CAMath::Abs(Bz) < 1.e-8f) {
993 Bz = 1.e-8f;
994 }
995
996 float dy = -2.f * mT0.Q() * mT0.Px() / Bz;
997 float dS; // path in XY
998 {
999 float chord = dy; // chord to the extrapolated point == |dy|*sign(x direction)
1000 float sa = -mT0.CosPhi(); // sin( half of the rotation angle ) == (chord/2) / radius
1001
1002 // dS = (Pt/b)*2*arcsin( sa )
1003 // = (Pt/b)*2*sa*(1 + 1/6 sa^2 + 3/40 sa^4 + 5/112 sa^6 +... )
1004 // = chord*(1 + 1/6 sa^2 + 3/40 sa^4 + 5/112 sa^6 +... )
1005
1006 float sa2 = sa * sa;
1007 const float k2 = 1.f / 6.f;
1008 const float k4 = 3.f / 40.f;
1009 // const float k6 = 5.f/112.f;
1010 dS = chord + chord * sa2 * (k2 + k4 * sa2);
1011 // dS = CAMath::Sqrt(pt2)/b*2.f*CAMath::ASin( sa );
1012 }
1013
1014 if (mT0.SinPhi() < 0.f) {
1015 dS = -dS;
1016 }
1017
1018 mT0.Y() = mT0.Y() + dy;
1019 mT0.Z() = mT0.Z() + mT0.DzDs() * dS;
1020 mT0.Px() = mT0.Px(); // should be positive
1021 mT->Y() = mT->Y() + dy;
1022 mT->Z() = mT->Z() + mT0.DzDs() * dS;
1023 ChangeDirection();
1024
1025 // Energy Loss
1026 if (true) {
1027 // std::cout<<"MIRROR: APPLY ENERGY LOSS!!!"<<std::endl;
1028
1029 float dL = CAMath::Abs(dS * mT0.GetDlDs());
1030
1031 if (inFlyDirection) {
1032 dL = -dL;
1033 }
1034
1035 float* c = mT->Cov();
1036 float& mC40 = c[10];
1037 float& mC41 = c[11];
1038 float& mC42 = c[12];
1039 float& mC43 = c[13];
1040 float& mC44 = c[14];
1041
1042 float dLmask = 0.f;
1043 bool maskMS = (CAMath::Abs(dL) < mMaterial.DLMax);
1044 if (maskMS) {
1045 dLmask = dL;
1046 }
1047 float dLabs = CAMath::Abs(dLmask);
1048 float corr = 1.f - mMaterial.EP2 * dLmask;
1049
1050 float corrInv = 1.f / corr;
1051 mT0.Px() *= corrInv;
1052 mT0.Py() *= corrInv;
1053 mT0.Pz() *= corrInv;
1054 mT0.Pt() *= corrInv;
1055 mT0.P() *= corrInv;
1056 mT0.QPt() *= corr;
1057
1058 mT->QPt() *= corr;
1059
1060 mC40 *= corr;
1061 mC41 *= corr;
1062 mC42 *= corr;
1063 mC43 *= corr;
1064 mC44 = mC44 * corr * corr + dLabs * mMaterial.sigmadE2;
1065 } else {
1066 // std::cout<<"MIRROR: DONT APPLY ENERGY LOSS!!!"<<std::endl;
1067 }
1068}
1069
1070GPUd() o2::base::MatBudget GPUTPCGMPropagator::getMatBudget(const float* p1, const float* p2)
1071{
1072 return mMatLUT->getMatBudget(p1[0], p1[1], p1[2], p2[0], p2[1], p2[2]);
1073}
1074
1076{
1077 float xyz1[3] = {getGlobalX(mT0.GetX(), mT0.GetY()), getGlobalY(mT0.GetX(), mT0.GetY()), mT0.GetZ()};
1078 float xyz2[3] = {getGlobalX(t0e.GetX(), t0e.GetY()), getGlobalY(t0e.GetX(), t0e.GetY()), t0e.GetZ()};
1079 o2::base::MatBudget mat = getMatBudget(xyz1, xyz2);
1080 if (mat.meanX2X0 > 1.e-8) {
1081 SetMaterial(mat.length / mat.meanX2X0, mat.meanRho);
1082 } else {
1083 SetMaterialTPC();
1084 }
1085}
int16_t charge
Definition RawEventData.h:5
int16_t time
Definition RawEventData.h:4
#define GPUdic(...)
#define GPUrestrict()
#define GPUCA_DEBUG_STREAMER_CHECK(...)
#define GPUCA_MAX_SIN_PHI
int32_t retVal
const int16_t bb
GPUd() void GPUTPCGMPropagator
#define GPUCA_NSECTORS
constexpr int p2()
constexpr int p1()
constexpr to accelerate the coordinates changing
uint32_t c
Definition RawData.h:2
Definition B.h:16
GLint GLenum GLint x
Definition glcorearb.h:403
GLboolean GLboolean GLboolean b
Definition glcorearb.h:1233
typedef void(APIENTRYP PFNGLCULLFACEPROC)(GLenum mode)
GLuint GLfloat x0
Definition glcorearb.h:5034
GLboolean r
Definition glcorearb.h:1233
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat t0
Definition glcorearb.h:5034
GLenum GLfloat param
Definition glcorearb.h:271
GLboolean GLboolean GLboolean GLboolean a
Definition glcorearb.h:1233
GLuint GLfloat GLfloat y0
Definition glcorearb.h:5034
float float float float z1
Definition MathUtils.h:44
value_T d2
Definition TrackUtils.h:135
kp1 *kp2 *value_T beta2
Definition TrackUtils.h:131
constexpr uint32_t TPC
Definition Triggers.h:39
a couple of static helper functions to create timestamp values for CCDB queries or override obsolete ...
float length
length in material
Definition MatCell.h:55
float meanRho
mean density, g/cm^3
Definition MatCell.h:30
float meanX2X0
fraction of radiaton lenght
Definition MatCell.h:31
std::vector< o2::mch::ChannelCode > cc