20#include "GPUParam.inc"
25#if defined(GPUCA_GM_USE_FULL_FIELD)
26#include "AliTracker.h"
36 float gx = getGlobalX(cosAlpha, sinAlpha, X, Y);
37 float gy = getGlobalY(cosAlpha, sinAlpha, X, Y);
39#if defined(GPUCA_GM_USE_FULL_FIELD)
40 const float kCLight = gpu_common_constants::kCLight;
41 double r[3] = {gx, gy, Z};
43 AliTracker::GetBxByBz(
r,
bb);
58 switch (mFieldRegion) {
60 mField->GetFieldIts(gx, gy, Z,
bb);
63 mField->GetFieldTrd(gx, gy, Z,
bb);
67 mField->GetField(gx, gy, Z,
bb);
74 B[0] =
bb[0] * cosAlpha +
bb[1] * sinAlpha;
75 B[1] = -
bb[0] * sinAlpha +
bb[1] * cosAlpha;
81 float gx = getGlobalX(cosAlpha, sinAlpha, X, Y);
82 float gy = getGlobalY(cosAlpha, sinAlpha, X, Y);
84#if defined(GPUCA_GM_USE_FULL_FIELD)
85 const float kCLight = gpu_common_constants::kCLight;
86 double r[3] = {gx, gy, Z};
88 AliTracker::GetBxByBz(
r,
bb);
89 return bb[2] * kCLight;
91 switch (mFieldRegion) {
93 return mField->GetFieldItsBz(gx, gy, Z);
95 return mField->GetFieldTrdBz(gx, gy, Z);
98 return mField->GetFieldBz(gx, gy, Z);
110 float newCosAlpha, newSinAlpha;
111 CAMath::SinCos(newAlpha, newSinAlpha, newCosAlpha);
113 float cc = newCosAlpha * mCosAlpha + newSinAlpha * mSinAlpha;
114 float ss = newSinAlpha * mCosAlpha - newCosAlpha * mSinAlpha;
120 float px0 = mT0.Px();
121 float py0 = mT0.Py();
124 if (CAMath::Abs(mT->GetSinPhi()) >= mMaxSinPhi || CAMath::Abs(px0) < (1 - mMaxSinPhi)) {
129 float px1 = px0 *
cc + py0 * ss;
130 float py1 = -px0 * ss + py0 *
cc;
140 if (CAMath::Abs(py1) > mMaxSinPhi * mT0.GetPt() || CAMath::Abs(px1) < (1 - mMaxSinPhi)) {
145 float trackX =
x0 *
cc + ss * mT->Y();
150 if (mPropagateBzOnly) {
151 Bz = GetBzBase(newCosAlpha, newSinAlpha,
t0.X(),
t0.Y(),
t0.Z());
152 if (
t0.PropagateToXBzLight(trackX, Bz, dLp)) {
157 GetBxByBz(newAlpha,
t0.X(),
t0.Y(),
t0.Z(),
B);
159 if (
t0.PropagateToXBxByBz(trackX,
B[0],
B[1],
B[2], dLp)) {
164 if (CAMath::Abs(
t0.SinPhi()) >= mMaxSinPhi) {
181 float j1 = px1 / px0;
187 mT->Y() = -
x0 * ss +
cc * mT->Y();
189 mT->SinPhi() = -CAMath::Sqrt(1.f - mT->SinPhi() * mT->SinPhi()) * ss + mT->SinPhi() *
cc;
192 float*
c = mT->Cov();
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;
209 if (!mFitInProjections && mT->NDF() > 0) {
216 if (
t0.SetDirectionAlongX()) {
217 mT->SinPhi() = -mT->SinPhi();
218 mT->DzDs() = -mT->DzDs();
219 mT->QPt() = -mT->QPt();
233 float j3 = -
t0.Py() /
t0.Px();
234 float j4 = -
t0.Pz() /
t0.Px();
235 float j5 =
t0.QPt() * Bz;
244 float h15 = c15 + c20 * j3;
245 float h16 = c16 + c20 * j4;
246 float h17 = c17 + c20 * j5;
248 c[0] += j3 * (c15 + h15);
250 c[2] += j4 * (c16 + h16);
252 c[3] += c17 * j3 + h15 * j5;
253 c[5] += j5 * (c17 + h17);
262 if (!mFitInProjections && mT->NDF() > 0) {
263 c[1] += c16 * j3 + h15 * j4;
264 c[4] += c17 * j4 + h16 * j5;
272 mCosAlpha = newCosAlpha;
273 mSinAlpha = newSinAlpha;
281 if (mPropagateBzOnly) {
282 return PropagateToXAlphaBz(posX, posAlpha, inFlyDirection);
284 if (CAMath::Abs(posAlpha - mAlpha) > 1.e-4f) {
285 if (RotateToAlpha(posAlpha) != 0) {
289 if (CAMath::Abs(posX - mT->X()) < 1.e-7f) {
295 GetBxByBz(mT0.X(), mT0.Y(), mT0.Z(),
B);
301 if (t0e.PropagateToXBxByBz(posX,
B[0],
B[1],
B[2], dLp) && t0e.PropagateToXBzLight(posX,
B[2], dLp)) {
305 if (CAMath::Abs(t0e.SinPhi()) >= mMaxSinPhi) {
309 return FollowLinearization(t0e,
B[2], dLp, inFlyDirection);
314 if (CAMath::Abs(posAlpha - mAlpha) > 1.e-4f) {
315 if (RotateToAlpha(posAlpha) != 0) {
319 if (CAMath::Abs(posX - mT->X()) < 1.e-7f) {
324 float Bz = GetBz(mT0.X(), mT0.Y(), mT0.Z());
330 if (t0e.PropagateToXBzLight(posX, Bz, dLp)) {
334 if (CAMath::Abs(t0e.SinPhi()) >= mMaxSinPhi) {
338 return FollowLinearization(t0e, Bz, dLp, inFlyDirection);
347 float dS = dLp * t0e.GetPt();
348 float dL = CAMath::Abs(dLp * t0e.GetP());
350 if (inFlyDirection) {
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();
361 float k = -mT0.QPt() * Bz;
362 float dx = t0e.GetX() - mT0.X();
365 float cci = 1.f /
cc;
367 float dxcci = dx * cci;
368 float hh = dxcci * ex1i * (1.f + ex * ex1 + ey * ey1);
370 float j02 = exi * hh;
371 float j04 = -Bz * dxcci * hh;
373 float j24 = -dx * Bz;
375 float* p = mT->Par();
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();
383 float newSinPhi = ey1 +
d2 + j24 * d4;
393 mT->X() = t0e.GetX();
394 p[0] = t0e.GetY() + d0 + j02 *
d2 + j04 * d4;
395 p[1] = t0e.GetZ() + d1 + j13 * d3;
397 p[3] = t0e.GetDzDs() + d3;
398 p[4] = t0e.GetQPt() + d4;
400 float*
c = mT->Cov();
417 if (mFitInProjections || mT->NDF() <= 0) {
418 float c20ph04c42 = c20 + j04 * c42;
419 float j02c22 = j02 * c22;
420 float j04c44 = j04 * c44;
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;
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);
436 c[8] = c32 + c43 * j24;
447 float xx = 1.f - 0.25f * kdx * kdx * (1.f + tg * tg);
451 xx = CAMath::Sqrt(xx);
452 float yy = CAMath::Sqrt(ss * ss +
cc *
cc);
454 float j12 = dx * mT0.DzDs() * tg * (2.f + tg * (ey * exi + ey1 * ex1i)) / (xx * yy);
456 if (CAMath::Abs(mT0.QPt()) > 1.e-6f) {
457 j14 = (2.f * xx * ex1i * dx / yy - dS) * mT0.DzDs() / mT0.QPt();
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);
462 p[1] += j12 *
d2 + j14 * d4;
464 float h00 = c00 + c20 * j02 + c40 * j04;
466 float h02 = c20 + c22 * j02 + c42 * j04;
468 float h04 = c40 + c42 * j02 + c44 * j04;
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;
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;
482 c[0] = h00 + h02 * j02 + h04 * j04;
484 c[1] = h10 + h12 * j02 + h14 * j04;
485 c[2] = h11 + h12 * j12 + h13 * j13 + h14 * j14;
487 c[3] = h20 + h22 * j02 + h24 * j04;
488 c[4] = h21 + h22 * j12 + h23 * j13 + h24 * j14;
489 c[5] = h22 + h24 * j24;
491 c[6] = c30 + c32 * j02 + c43 * j04;
492 c[7] = c31 + c32 * j12 + c33 * j13 + c43 * j14;
493 c[8] = c32 + c43 * j24;
496 c[10] = c40 + c42 * j02 + c44 * j04;
497 c[11] = c41 + c42 * j12 + c43 * j13 + c44 * j14;
498 c[12] = c42 + c44 * j24;
512 bool maskMS = (CAMath::Abs(dL) < mMaterial.DLMax);
516 float dLabs = CAMath::Abs(dLmask);
521 float corr = 1.f - mMaterial.EP2 * dLmask;
522 float corrInv = 1.f / corr;
536 mC44 = mC44 * corr * corr + dLabs * mMaterial.sigmadE2;
543 mC22 += dLabs * mMaterial.k22 * mT0.CosPhi() * mT0.CosPhi();
544 mC33 += dLabs * mMaterial.k33;
545 mC43 += dLabs * mMaterial.k43;
546 mC44 += dLabs * mMaterial.k44;
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;
563 float ex1 = CAMath::Sqrt(1.f - ey1 * ey1);
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);
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));
580 float dz = dS * mT0.DzDs();
581 projY = mT->Y() + dy;
582 projZ = mT->Z() + dz;
588 GetErr2(err2Y, err2Z,
param, mT0.GetSinPhi(), mT0.DzDs(), posZ, mT->GetX(), mT->GetY(), iRow, clusterState, sector,
time, avgCharge,
charge, mSeedingErrors);
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)
593#ifndef GPUCA_TPC_GEOMETRY_O2
595 param.GetClusterErrorsSeeding2(sector, iRow, posZ, snp, tgl,
time, err2Y, err2Z);
599 param.GetClusterErrors2(sector, iRow, posZ, snp, tgl,
time, avgCharge,
charge, err2Y, err2Z);
601 param.UpdateClusterError2ByState(clusterState, err2Y, err2Z);
602 float statErr2 =
param.GetSystematicClusterErrorIFC2(trackX, trackY, posZ, sector >= (
GPUCA_NSECTORS / 2));
604 statErr2 +=
param.GetSystematicClusterErrorC122(trackX, trackY, sector);
613 GetErr2(err2Y, err2Z,
param, posZ, iRow, clusterState, sector,
time, avgCharge,
charge);
614 return PredictChi2(posY, posZ, err2Y, err2Z);
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];
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;
629 float w0 = mC[2] + err2Z, w1 = mC[1], w2 = mC[0] + err2Y;
631 float det = w0 * w2 - w1 * w1;
632 if (CAMath::Abs(det) < 1.e-10f) {
640 return CAMath::Abs((w0 * z0 + w1 * z1) * z0) + CAMath::Abs((w1 * z0 + w2 * z1) * z1);
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))
647 GetErr2(err2Y, err2Z,
param, posZ, iRow, clusterState, sector,
time, avgInvCharge, invCharge);
650 if (rejectChi2 >= rejectInterFill) {
651 if (rejectChi2 == rejectInterReject && inter->errorY < (GPUCA_PAR_MERGER_INTERPOLATION_ERROR_TYPE_A)0) {
652 rejectChi2 = rejectDirect;
654 int32_t
retVal = InterpolateReject(
param, posY, posZ, clusterState, rejectChi2, inter, err2Y, err2Z);
662 if (mT->NDF() == -5) {
663 mT->ResetCovariance();
664 float* mC = mT->Cov();
665 float* mP = mT->Par();
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);
679 return Update(posY, posZ, clusterState, rejectChi2 == rejectDirect || (
param.rec.tpc.mergerInterpolateRejectAlsoOnCurrentPosition && rejectChi2 == rejectInterReject), err2Y, err2Z, &
param);
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)
686 if (rejectChi2 == rejectInterFill) {
689 inter->errorY = mC[0];
690 inter->errorZ = mC[2];
691 }
else if (rejectChi2 == rejectInterReject) {
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];
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;
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]);
719 const float Iw1 = mC[1] * 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);
737 const float Jw1 = ImC1 * Jdet;
739 chi2Y = CAMath::Abs((Jw0 * Jz0 + Jw1 * Jz1) * Jz0);
740 chi2Z = CAMath::Abs((Jw1 * Jz0 + Jw2 * Jz1) * Jz1);
742 if (RejectCluster(chi2Y *
param.rec.tpc.clusterRejectChi2TolleranceY, chi2Z *
param.rec.tpc.clusterRejectChi2TolleranceZ, clusterState)) {
743 return updateErrorClusterRejectedInInterpolation;
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];
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);
763 w2 = 1.f / (err2Z + d11);
764 chi2Y = w0 * z0 * z0;
765 chi2Z = w2 *
z1 *
z1;
767 w0 = d11 + err2Z, w1 = d10, w2 = d00 + err2Y;
769 float det = w0 * w2 - w1 * w1;
770 if (CAMath::Abs(det) < 1.e-10f) {
771 return updateErrorFitFailed;
778 chi2Y = CAMath::Abs((w0 * z0 + w1 * z1) * z0);
779 chi2Z = CAMath::Abs((w1 * z0 + w2 * z1) * z1);
781 float dChi2 = chi2Y + chi2Z;
783 if (rejectChi2 && RejectCluster(chi2Y *
param->rec.tpc.clusterRejectChi2TolleranceY, chi2Z *
param->rec.tpc.clusterRejectChi2TolleranceZ, clusterState)) {
784 return updateErrorClusterRejectedInUpdate;
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;
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;
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;
828 mC[0] -= k00 * d00 + k01 * d10;
830 mC[2] -= k10 * d01 + k11 * d11;
832 mC[3] -= k20 * d00 + k21 * d10;
833 mC[5] -= k20 * d02 + k21 * d12;
835 mC[7] -= k30 * d01 + k31 * d11;
836 mC[9] -= k30 * d03 + k31 * d13;
838 mC[10] -= k40 * d00 + k41 * d10;
839 mC[12] -= k40 * d02 + k41 * d12;
840 mC[14] -= k40 * d04 + k41 * d14;
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;
864 const float log0 = CAMath::Log(5940.f);
865 const float log1 = CAMath::Log(3.5f * 5940.f);
867 bool bad = (
beta2 >= .999f) || (beta2 < 1.e-8f);
874 float b = 0.5f * CAMath::Log(
a);
875 float d = 0.153e-3f /
beta2;
878 float ret = d * (log0 +
b +
c);
879 float case1 = d * (log1 +
c);
881 if (
a > 3.5f * 3.5f) {
895 const float mass = 0.13957f;
897 float qpt = mT0.GetQPt();
898 if (CAMath::Abs(qpt) > 20) {
902 float w2 = (1.f + mT0.GetDzDs() * mT0.GetDzDs());
903 float pti2 = qpt * qpt;
908 float mass2 = mass * mass;
909 float beta2 = w2 / (w2 + mass2 * pti2);
911 float p2 = w2 / pti2;
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;
916 mMaterial.EP2 = E /
p2;
920 const float knst = 0.0007f;
921 mMaterial.sigmadE2 = knst * mMaterial.EP2 * qpt;
922 mMaterial.sigmadE2 = mMaterial.sigmadE2 * mMaterial.sigmadE2;
924 mMaterial.k22 = theta2 * w2;
925 mMaterial.k33 = mMaterial.k22 * w2;
927 mMaterial.k44 = theta2 * mT0.GetDzDs() * mT0.GetDzDs() * pti2;
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;
940 mT0.Pz() = -mT0.Pz();
945 mT->QPt() = -mT->QPt();
946 mT->DzDs() = -mT->DzDs();
948 mAlpha = mAlpha + CAMath::Pi();
949 while (mAlpha >= CAMath::Pi()) {
950 mAlpha -= CAMath::TwoPi();
952 while (mAlpha < -CAMath::Pi()) {
953 mAlpha += CAMath::TwoPi();
955 mCosAlpha = -mCosAlpha;
956 mSinAlpha = -mSinAlpha;
958 float*
c = mT->Cov();
969 mT0.Py() = -mT0.Py();
970 mT0.Pz() = -mT0.Pz();
972 mT->SinPhi() = -mT->SinPhi();
973 mT->DzDs() = -mT->DzDs();
974 mT->QPt() = -mT->QPt();
977 float*
c = mT->Cov();
990 GetBxByBz(mT0.X(), mT0.Y(), mT0.Z(),
B);
992 if (CAMath::Abs(Bz) < 1.e-8f) {
996 float dy = -2.f * mT0.Q() * mT0.Px() / Bz;
1000 float sa = -mT0.CosPhi();
1006 float sa2 = sa * sa;
1007 const float k2 = 1.f / 6.f;
1008 const float k4 = 3.f / 40.f;
1010 dS = chord + chord * sa2 * (k2 + k4 * sa2);
1014 if (mT0.SinPhi() < 0.f) {
1018 mT0.Y() = mT0.Y() + dy;
1019 mT0.Z() = mT0.Z() + mT0.DzDs() * dS;
1020 mT0.Px() = mT0.Px();
1021 mT->Y() = mT->Y() + dy;
1022 mT->Z() = mT->Z() + mT0.DzDs() * dS;
1029 float dL = CAMath::Abs(dS * mT0.GetDlDs());
1031 if (inFlyDirection) {
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];
1043 bool maskMS = (CAMath::Abs(dL) < mMaterial.DLMax);
1047 float dLabs = CAMath::Abs(dLmask);
1048 float corr = 1.f - mMaterial.EP2 * dLmask;
1050 float corrInv = 1.f / corr;
1051 mT0.Px() *= corrInv;
1052 mT0.Py() *= corrInv;
1053 mT0.Pz() *= corrInv;
1054 mT0.Pt() *= corrInv;
1064 mC44 = mC44 * corr * corr + dLabs * mMaterial.sigmadE2;
1072 return mMatLUT->getMatBudget(
p1[0],
p1[1],
p1[2],
p2[0],
p2[1],
p2[2]);
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()};
#define GPUCA_DEBUG_STREAMER_CHECK(...)
#define GPUCA_MAX_SIN_PHI
GPUd() void GPUTPCGMPropagator
constexpr int p1()
constexpr to accelerate the coordinates changing
GLboolean GLboolean GLboolean b
typedef void(APIENTRYP PFNGLCULLFACEPROC)(GLenum mode)
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat t0
GLboolean GLboolean GLboolean GLboolean a
GLuint GLfloat GLfloat y0
float float float float z1
a couple of static helper functions to create timestamp values for CCDB queries or override obsolete ...
float length
length in material
float meanRho
mean density, g/cm^3
float meanX2X0
fraction of radiaton lenght
std::vector< o2::mch::ChannelCode > cc