36 return dx * dx + dy * dy + dz *
dz;
45 return dx * dx + dz *
dz;
57 double dS = x * ex + y * ey;
58 if (std::abs(k) > 1.
e-4
f) {
59 dS = atan2(k * dS, 1 + k * (x * ey - y * ex)) /
k;
75 double ax = dx * k + ey;
76 double ay = dy * k - ex;
77 double a = std::sqrt(ax * ax + ay * ay);
78 xp = x0 + (dx - ey * ((dx * dx + dy *
dy) * k - 2 * (-dx * ey + dy * ex)) / (a + 1)) / a;
79 yp = y0 + (dy + ex * ((dx * dx + dy *
dy) * k - 2 * (-dx * ey + dy * ex)) / (a + 1)) / a;
80 double s =
GetS(x, y, Bz);
82 if (std::abs(k) > 1.
e-2
f) {
83 double dZ = std::abs(
GetDzDs() * 2*M_PI / k);
85 zp +=
round((z - zp) / dZ) * dZ;
104 double k = -t0.
QPt() *
Bz;
107 double ey1 = k * dx + ey;
112 if (std::abs(ey1) > maxSinPhi) {
116 ex1 = std::sqrt(1 - ey1 * ey1);
121 double dx2 = dx * dx;
122 double ss = ey + ey1;
123 double cc = ex + ex1;
125 if (std::abs(cc) < 1.
e-4
f || std::abs(ex) < 1.
e-4
f || std::abs(ex1) < 1.
e-4
f) {
132 double dl = dx * std::sqrt(1 + tg * tg);
137 double dSin = dl * k / 2;
144 double dS = (std::abs(k) > 1.e-4
f) ? (2 * asin(dSin) /
k) : dl;
145 double dz = dS * t0.
DzDs();
148 *DL = -dS * std::sqrt(1 + t0.
DzDs() * t0.
DzDs());
151 double cci = 1.f / cc;
152 double exi = 1.f / ex;
153 double ex1i = 1.f / ex1;
163 double h2 = dx * (1 + ey * ey1 + ex * ex1) * exi * ex1i * cci;
164 double h4 = dx2 * (cc + ss * ey1 * ex1i) * cci * cci * (-Bz);
165 double dxBz = dx * (-
Bz);
171 SetPar(0,
Y() + dy + h2 * d[2] + h4 * d[4]);
172 SetPar(1,
Z() + dz + dS * d[3]);
191 mC[0] = c00 + h2 * h2 * c22 + h4 * h4 * c44 + 2 * (h2 * c20 + h4 * c40 + h2 * h4 * c42);
193 mC[1] = c10 + h2 * c21 + h4 * c41 + dS * (c30 + h2 * c32 + h4 * c43);
194 mC[2] = c11 + 2 * dS * c31 + dS * dS * c33;
196 mC[3] = c20 + h2 * c22 + h4 * c42 + dxBz * (c40 + h2 * c42 + h4 * c44);
197 mC[4] = c21 + dS * c32 + dxBz * (c41 + dS * c43);
198 mC[5] = c22 + 2 * dxBz * c42 + dxBz * dxBz * c44;
200 mC[6] = c30 + h2 * c32 + h4 * c43;
201 mC[7] = c31 + dS * c33;
202 mC[8] = c32 + dxBz * c43;
205 mC[10] = c40 + h2 * c42 + h4 * c44;
206 mC[11] = c41 + dS * c43;
207 mC[12] = c42 + dxBz * c44;
227 if (std::abs(ex) < 1.
e-4
f) {
230 double exi = 1.f / ex;
232 double dxBz = dx * (-
Bz);
233 double dS = dx * exi;
234 double h2 = dS * exi * exi;
235 double h4 = .5f * h2 * dxBz;
244 if (maxSinPhi > 0 && std::abs(sinPhi) > maxSinPhi) {
271 mC[0] = std::max(c00, c00 + h2 * h2 * c22 + h4 * h4 * c44 + 2 * (h2 * c20 + h4 * c40 + h2 * h4 * c42));
273 mC[1] = c10 + h2 * c21 + h4 * c41 + dS * (c30 + h2 * c32 + h4 * c43);
274 mC[2] = std::max(c11, c11 + 2 * dS * c31 + dS * dS * c33);
276 mC[3] = c20 + h2 * c22 + h4 * c42 + dxBz * (c40 + h2 * c42 + h4 * c44);
277 mC[4] = c21 + dS * c32 + dxBz * (c41 + dS * c43);
278 mC[5] = c22 + 2 * dxBz * c42 + dxBz * dxBz * c44;
280 mC[6] = c30 + h2 * c32 + h4 * c43;
281 mC[7] = c31 + dS * c33;
282 mC[8] = c32 + dxBz * c43;
285 mC[10] = c40 + h2 * c42 + h4 * c44;
286 mC[11] = c41 + dS * c43;
287 mC[12] = c42 + dxBz * c44;
307 const double kRho = 2.27925e-3
f;
308 const double kRadLen = 14.403f;
309 const double kRhoOverRadLen = kRho / kRadLen;
356 const double mK = 0.307075e-3
f;
357 const double me = 0.511e-3
f;
358 const double rho = kp0;
359 const double x0 = kp1 * 2.303f;
360 const double x1 = kp2 * 2.303f;
361 const double mI = kp3;
362 const double mZA = kp4;
363 const double maxT = 2 * me * bg2;
367 const double x = 0.5f * log(bg2);
368 const double lhwI = log(28.816
f * 1
e-9
f * std::sqrt(rho * mZA) / mI);
370 d2 = lhwI + x - 0.5f;
372 const double r = (x1 -
x) / (x1 - x0);
373 d2 = lhwI + x - 0.5f + (0.5f - lhwI -
x0) * r * r * r;
376 return mK * mZA * (1 + bg2) / bg2 * (0.5
f * log(2 * me * bg2 * maxT / (mI * mI)) - bg2 / (1 + bg2) - d2);
407 const double rho = 2.27925e-3
f;
408 const double x0 = 2.f;
409 const double x1 = 4.f;
410 const double mI = 14.e-9
f;
411 const double mZA = 0.47999f;
426 if (beta2 / (1 - beta2) > 3.5
f * 3.5
f) {
427 return 0.153e-3f / beta2 * (log(3.5f * 5940) + 0.5f * log(beta2 / (1 - beta2)) - beta2);
429 return 0.153e-3f / beta2 * (log(5940 * beta2 / (1 - beta2)) - beta2);
442 double k2 = qpt * qpt;
443 double mass2 = mass *
mass;
444 double beta2 = p2 / (p2 + mass2 * k2);
446 double pp2 = (k2 > 1.e-8
f) ? p2 / k2 : 10000;
450 par.
e = std::sqrt(pp2 + mass2);
451 par.
theta2 = 14.1f * 14.1f / (beta2 * pp2 * 1e6f);
452 par.
EP2 = par.
e / pp2;
456 const double knst = 0.07f;
474 double& mC22 =
mC[5];
475 double& mC33 =
mC[9];
476 double& mC40 =
mC[10];
477 double& mC41 =
mC[11];
478 double& mC42 =
mC[12];
479 double& mC43 =
mC[13];
480 double& mC44 =
mC[14];
484 double dE = par.
bethe * xTimesRho;
485 if (std::abs(dE) > 0.3
f * par.
e) {
488 double corr = (1.f - par.
EP2 * dE);
489 if (corr < 0.3f || corr > 1.3
f) {
499 mC44 += par.
sigmadE2 * std::abs(dE);
503 double theta2 = par.
theta2 * std::abs(xOverX0);
505 mC33 += theta2 * par.
k33;
506 mC43 += theta2 * par.
k43;
507 mC44 += theta2 * par.
k44;
519 double cA = cos(alpha);
520 double sA = sin(alpha);
522 double cosPhi = cP * cA + sP * sA;
523 double sinPhi = -cP * sA + sP * cA;
525 if (std::abs(sinPhi) > maxSinPhi || std::abs(cosPhi) < 1.
e-2
f || std::abs(cP) < 1.
e-2
f) {
529 double j0 = cP / cosPhi;
530 double j2 = cosPhi / cP;
532 SetX(x * cA +
y * sA);
533 SetY(-x * sA +
y * cA);
576 double cA = cos(alpha);
577 double sA = sin(alpha);
579 double cosPhi = cP * cA + sP * sA;
580 double sinPhi = -cP * sA + sP * cA;
582 if (std::abs(sinPhi) > maxSinPhi || std::abs(cosPhi) < 1.
e-2
f || std::abs(cP) < 1.
e-2
f) {
592 double j0 = cP / cosPhi;
593 double j2 = cosPhi / cP;
594 double d[2] = {
Y() - y0,
SinPhi() - sP};
596 SetX(x0 * cA + y0 * sA);
597 SetY(-x0 * sA + y0 * cA + j0 * d[0]);
622 double c00 =
mC[0], c11 =
mC[2], c20 =
mC[3], c31 =
mC[7], c40 =
mC[10];
629 if (err2Y < 1.
e-8
f || err2Z < 1.
e-8
f) {
633 double mS0 = 1.f / err2Y;
634 double mS2 = 1.f / err2Z;
638 double k00, k11, k20, k31, k40;
647 double sinPhi =
GetPar(2) + k20 * z0;
649 if (maxSinPhi > 0 && std::abs(sinPhi) >= maxSinPhi) {
663 mChi2 += mS0 * z0 * z0 + mS2 * z1 * z1;
685 const double*
c =
Cov();
686 for (
int i = 0;
i < 15;
i++) {
687 ok = ok && std::isfinite(c[
i]);
689 for (
int i = 0;
i < 5;
i++) {
690 ok = ok && std::isfinite(
Par()[
i]);
693 if (c[0] <= 0 || c[2] <= 0 || c[5] <= 0 || c[9] <= 0 || c[14] <= 0) {
696 if (c[0] > 5.
f || c[2] > 5.
f || c[5] > 2.
f || c[9] > 2
705 if (std::abs(
QPt()) > 1.
f / 0.05
f) {
709 ok = ok && (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]) && (c[10] * c[10] <= c[14] * c[0]) && (c[11] * c[11] <= c[14] * c[2]) &&
710 (c[12] * c[12] <= c[14] * c[5]) && (c[13] * c[13] <= c[14] * c[9]);
715 #if !defined(GPUCA_GPUCODE)
723 #if !defined(GPUCA_GPUCODE)