ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
GPUTPCTrackParam.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file GPUTPCTrackParam.cxx
1 // Copyright CERN and copyright holders of ALICE O2. This software is
2 // distributed under the terms of the GNU General Public License v3 (GPL
3 // Version 3), copied verbatim in the file "COPYING".
4 //
5 // See http://alice-o2.web.cern.ch/license for full licensing information.
6 //
7 // In applying this license CERN does not waive the privileges and immunities
8 // granted to it by virtue of its status as an Intergovernmental Organization
9 // or submit itself to any jurisdiction.
10 
13 
15 #include "GPUTPCTrackParam.h"
16 #include <cmath>
17 #include <algorithm>
18 
19 //
20 // Circle in XY:
21 //
22 // kCLight = 0.000299792458;
23 // Kappa = -Bz*kCLight*QPt;
24 // R = 1/fstd::absf(Kappa);
25 // Xc = X - sin(Phi)/Kappa;
26 // Yc = Y + cos(Phi)/Kappa;
27 //
28 
30 {
31  // get squared distance between tracks
32 
33  double dx = GetX() - t.GetX();
34  double dy = GetY() - t.GetY();
35  double dz = GetZ() - t.GetZ();
36  return dx * dx + dy * dy + dz * dz;
37 }
38 
40 {
41  // get squared distance between tracks in X&Z
42 
43  double dx = GetX() - t.GetX();
44  double dz = GetZ() - t.GetZ();
45  return dx * dx + dz * dz;
46 }
47 
48 double GPUTPCTrackParam::GetS(double x, double y, double Bz) const
49 {
50  //* Get XY path length to the given point
51 
52  double k = GetKappa(Bz);
53  double ex = GetCosPhi();
54  double ey = GetSinPhi();
55  x -= GetX();
56  y -= GetY();
57  double dS = x * ex + y * ey;
58  if (std::abs(k) > 1.e-4f) {
59  dS = atan2(k * dS, 1 + k * (x * ey - y * ex)) / k;
60  }
61  return dS;
62 }
63 
64 void GPUTPCTrackParam::GetDCAPoint(double x, double y, double z, double& xp, double& yp, double& zp, double Bz) const
65 {
66  //* Get the track point closest to the (x,y,z)
67 
68  double x0 = GetX();
69  double y0 = GetY();
70  double k = GetKappa(Bz);
71  double ex = GetCosPhi();
72  double ey = GetSinPhi();
73  double dx = x - x0;
74  double dy = y - y0;
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);
81  zp = GetZ() + GetDzDs() * s;
82  if (std::abs(k) > 1.e-2f) {
83  double dZ = std::abs(GetDzDs() * 2*M_PI / k);
84  if (dZ > .1f) {
85  zp += round((z - zp) / dZ) * dZ;
86  }
87  }
88 }
89 
90 //*
91 //* Transport routines
92 //*
93 
94 bool GPUTPCTrackParam::TransportToX(double x, GPUTPCTrackLinearisation& t0, double Bz, double maxSinPhi, double* DL)
95 {
96  //* Transport the track parameters to X=x, using linearization at t0, and the field value Bz
97  //* maxSinPhi is the max. allowed value for |t0.SinPhi()|
98  //* linearisation of trajectory t0 is also transported to X=x,
99  //* returns 1 if OK
100  //*
101 
102  double ex = t0.CosPhi();
103  double ey = t0.SinPhi();
104  double k = -t0.QPt() * Bz;
105  double dx = x - X();
106 
107  double ey1 = k * dx + ey;
108  double ex1;
109 
110  // check for intersection with X=x
111 
112  if (std::abs(ey1) > maxSinPhi) {
113  return 0;
114  }
115 
116  ex1 = std::sqrt(1 - ey1 * ey1);
117  if (ex < 0) {
118  ex1 = -ex1;
119  }
120 
121  double dx2 = dx * dx;
122  double ss = ey + ey1;
123  double cc = ex + ex1;
124 
125  if (std::abs(cc) < 1.e-4f || std::abs(ex) < 1.e-4f || std::abs(ex1) < 1.e-4f) {
126  return 0;
127  }
128 
129  double tg = ss / cc; // tanf((phi1+phi)/2)
130 
131  double dy = dx * tg;
132  double dl = dx * std::sqrt(1 + tg * tg);
133 
134  if (cc < 0) {
135  dl = -dl;
136  }
137  double dSin = dl * k / 2;
138  if (dSin > 1) {
139  dSin = 1;
140  }
141  if (dSin < -1) {
142  dSin = -1;
143  }
144  double dS = (std::abs(k) > 1.e-4f) ? (2 * asin(dSin) / k) : dl;
145  double dz = dS * t0.DzDs();
146 
147  if (DL) {
148  *DL = -dS * std::sqrt(1 + t0.DzDs() * t0.DzDs());
149  }
150 
151  double cci = 1.f / cc;
152  double exi = 1.f / ex;
153  double ex1i = 1.f / ex1;
154 
155  double d[5] = {0, 0, GetPar(2) - t0.SinPhi(), GetPar(3) - t0.DzDs(), GetPar(4) - t0.QPt()};
156 
157  // double H0[5] = { 1,0, h2, 0, h4 };
158  // double H1[5] = { 0, 1, 0, dS, 0 };
159  // double H2[5] = { 0, 0, 1, 0, dxBz };
160  // double H3[5] = { 0, 0, 0, 1, 0 };
161  // double H4[5] = { 0, 0, 0, 0, 1 };
162 
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);
166 
167  t0.SetCosPhi(ex1);
168  t0.SetSinPhi(ey1);
169 
170  SetX(X() + dx);
171  SetPar(0, Y() + dy + h2 * d[2] + h4 * d[4]);
172  SetPar(1, Z() + dz + dS * d[3]);
173  SetPar(2, t0.SinPhi() + d[2] + dxBz * d[4]);
174 
175  double c00 = mC[0];
176  double c10 = mC[1];
177  double c11 = mC[2];
178  double c20 = mC[3];
179  double c21 = mC[4];
180  double c22 = mC[5];
181  double c30 = mC[6];
182  double c31 = mC[7];
183  double c32 = mC[8];
184  double c33 = mC[9];
185  double c40 = mC[10];
186  double c41 = mC[11];
187  double c42 = mC[12];
188  double c43 = mC[13];
189  double c44 = mC[14];
190 
191  mC[0] = c00 + h2 * h2 * c22 + h4 * h4 * c44 + 2 * (h2 * c20 + h4 * c40 + h2 * h4 * c42);
192 
193  mC[1] = c10 + h2 * c21 + h4 * c41 + dS * (c30 + h2 * c32 + h4 * c43);
194  mC[2] = c11 + 2 * dS * c31 + dS * dS * c33;
195 
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;
199 
200  mC[6] = c30 + h2 * c32 + h4 * c43;
201  mC[7] = c31 + dS * c33;
202  mC[8] = c32 + dxBz * c43;
203  mC[9] = c33;
204 
205  mC[10] = c40 + h2 * c42 + h4 * c44;
206  mC[11] = c41 + dS * c43;
207  mC[12] = c42 + dxBz * c44;
208  mC[13] = c43;
209  mC[14] = c44;
210 
211  return 1;
212 }
213 
214 bool GPUTPCTrackParam::TransportToX(double x, double sinPhi0, double cosPhi0, double Bz, double maxSinPhi)
215 {
216  //* Transport the track parameters to X=x, using linearization at phi0 with 0 curvature,
217  //* and the field value Bz
218  //* maxSinPhi is the max. allowed value for |t0.SinPhi()|
219  //* linearisation of trajectory t0 is also transported to X=x,
220  //* returns 1 if OK
221  //*
222 
223  double ex = cosPhi0;
224  double ey = sinPhi0;
225  double dx = x - X();
226 
227  if (std::abs(ex) < 1.e-4f) {
228  return 0;
229  }
230  double exi = 1.f / ex;
231 
232  double dxBz = dx * (-Bz);
233  double dS = dx * exi;
234  double h2 = dS * exi * exi;
235  double h4 = .5f * h2 * dxBz;
236 
237  // double H0[5] = { 1,0, h2, 0, h4 };
238  // double H1[5] = { 0, 1, 0, dS, 0 };
239  // double H2[5] = { 0, 0, 1, 0, dxBz };
240  // double H3[5] = { 0, 0, 0, 1, 0 };
241  // double H4[5] = { 0, 0, 0, 0, 1 };
242 
243  double sinPhi = SinPhi() + dxBz * QPt();
244  if (maxSinPhi > 0 && std::abs(sinPhi) > maxSinPhi) {
245  return 0;
246  }
247 
248  SetX(X() + dx);
249  SetPar(0, GetPar(0) + dS * ey + h2 * (SinPhi() - ey) + h4 * QPt());
250  SetPar(1, GetPar(1) + dS * DzDs());
251  SetPar(2, sinPhi);
252 
253  double c00 = mC[0];
254  double c10 = mC[1];
255  double c11 = mC[2];
256  double c20 = mC[3];
257  double c21 = mC[4];
258  double c22 = mC[5];
259  double c30 = mC[6];
260  double c31 = mC[7];
261  double c32 = mC[8];
262  double c33 = mC[9];
263  double c40 = mC[10];
264  double c41 = mC[11];
265  double c42 = mC[12];
266  double c43 = mC[13];
267  double c44 = mC[14];
268 
269  // This is not the correct propagation!!! The max ensures the positional error does not decrease during track finding.
270  // A different propagation method is used for the fit.
271  mC[0] = std::max(c00, c00 + h2 * h2 * c22 + h4 * h4 * c44 + 2 * (h2 * c20 + h4 * c40 + h2 * h4 * c42));
272 
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);
275 
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;
279 
280  mC[6] = c30 + h2 * c32 + h4 * c43;
281  mC[7] = c31 + dS * c33;
282  mC[8] = c32 + dxBz * c43;
283  mC[9] = c33;
284 
285  mC[10] = c40 + h2 * c42 + h4 * c44;
286  mC[11] = c41 + dS * c43;
287  mC[12] = c42 + dxBz * c44;
288  mC[13] = c43;
289  mC[14] = c44;
290 
291  return 1;
292 }
293 
294 bool GPUTPCTrackParam::TransportToX(double x, double Bz, double maxSinPhi)
295 {
296  //* Transport the track parameters to X=x
297  GPUTPCTrackLinearisation t0(*this);
298  return TransportToX(x, t0, Bz, maxSinPhi);
299 }
300 
302 {
303  //* Transport the track parameters to X=x taking into account material budget
304 
305 // const double kRho = 1.025e-3f; // 0.9e-3;
306 // const double kRadLen = 29.532f; // 28.94;
307  const double kRho = 2.27925e-3f;
308  const double kRadLen = 14.403f;
309  const double kRhoOverRadLen = kRho / kRadLen;
310  double dl;
311 
312  if (!TransportToX(x, t0, Bz, maxSinPhi, &dl)) {
313  return 0;
314  }
315 
316  CorrectForMeanMaterial(dl * kRhoOverRadLen, dl * kRho, par);
317  return 1;
318 }
319 
320 bool GPUTPCTrackParam::TransportToXWithMaterial(double x, GPUTPCTrackFitParam& par, double Bz, double maxSinPhi)
321 {
322  //* Transport the track parameters to X=x taking into account material budget
323 
324  GPUTPCTrackLinearisation t0(*this);
325  return TransportToXWithMaterial(x, t0, par, Bz, maxSinPhi);
326 }
327 
328 bool GPUTPCTrackParam::TransportToXWithMaterial(double x, double Bz, double maxSinPhi)
329 {
330  //* Transport the track parameters to X=x taking into account material budget
331 
334  return TransportToXWithMaterial(x, par, Bz, maxSinPhi);
335 }
336 
337 //*
338 //* Multiple scattering and energy losses
339 //*
340 double GPUTPCTrackParam::BetheBlochGeant(double bg2, double kp0, double kp1, double kp2, double kp3, double kp4)
341 {
342  //
343  // This is the parameterization of the Bethe-Bloch formula inspired by Geant.
344  //
345  // bg2 - (beta*gamma)^2
346  // kp0 - density [g/cm^3]
347  // kp1 - density effect first junction point
348  // kp2 - density effect second junction point
349  // kp3 - mean excitation energy [GeV]
350  // kp4 - mean Z/A
351  //
352  // The default values for the kp* parameters are for silicon.
353  // The returned value is in [GeV/(g/cm^2)].
354  //
355 
356  const double mK = 0.307075e-3f; // [GeV*cm^2/g]
357  const double me = 0.511e-3f; // [GeV/c^2]
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; // neglecting the electron mass
364 
365  //*** Density effect
366  double d2 = 0.f;
367  const double x = 0.5f * log(bg2);
368  const double lhwI = log(28.816f * 1e-9f * std::sqrt(rho * mZA) / mI);
369  if (x > x1) {
370  d2 = lhwI + x - 0.5f;
371  } else if (x > x0) {
372  const double r = (x1 - x) / (x1 - x0);
373  d2 = lhwI + x - 0.5f + (0.5f - lhwI - x0) * r * r * r;
374  }
375 
376  return mK * mZA * (1 + bg2) / bg2 * (0.5f * log(2 * me * bg2 * maxT / (mI * mI)) - bg2 / (1 + bg2) - d2);
377 }
378 
380 {
381  //------------------------------------------------------------------
382  // This is an approximation of the Bethe-Bloch formula,
383  // reasonable for solid materials.
384  // All the parameters are, in fact, for Si.
385  // The returned value is in [GeV]
386  //------------------------------------------------------------------
387 
388  return BetheBlochGeant(bg);
389 }
390 
392 {
393  //------------------------------------------------------------------
394  // This is an approximation of the Bethe-Bloch formula,
395  // reasonable for gas materials.
396  // All the parameters are, in fact, for Ne.
397  // The returned value is in [GeV]
398  //------------------------------------------------------------------
399 
400 // const double rho = 0.9e-3f;
401 // const double x0 = 2.f;
402 // const double x1 = 4.f;
403 // const double mI = 140.e-9f;
404 // const double mZA = 0.49555f;
405 
406  // these are for the sPHENIX TPC gas mixture
407  const double rho = 2.27925e-3f;
408  const double x0 = 2.f;
409  const double x1 = 4.f;
410  const double mI = 14.e-9f;
411  const double mZA = 0.47999f;
412  return BetheBlochGeant(bg, rho, x0, x1, mI, mZA);
413 }
414 
416 {
417  //------------------------------------------------------------------
418  // This is an approximation of the Bethe-Bloch formula with
419  // the density effect taken into account at beta*gamma > 3.5
420  // (the approximation is reasonable only for solid materials)
421  //------------------------------------------------------------------
422  if (beta2 >= 1) {
423  return 0;
424  }
425 
426  if (beta2 / (1 - beta2) > 3.5f * 3.5f) {
427  return 0.153e-3f / beta2 * (log(3.5f * 5940) + 0.5f * log(beta2 / (1 - beta2)) - beta2);
428  }
429  return 0.153e-3f / beta2 * (log(5940 * beta2 / (1 - beta2)) - beta2);
430 }
431 
433 {
434  //*!
435 
436  double qpt = GetPar(4);
437  if (mC[14] >= 1.f) {
438  qpt = 1.f / 0.35f;
439  }
440 
441  double p2 = (1.f + GetPar(3) * GetPar(3));
442  double k2 = qpt * qpt;
443  double mass2 = mass * mass;
444  double beta2 = p2 / (p2 + mass2 * k2);
445 
446  double pp2 = (k2 > 1.e-8f) ? p2 / k2 : 10000; // impuls 2
447 
448  par.bethe = BetheBlochGas( pp2/mass2);
449  //par.bethe = ApproximateBetheBloch(pp2 / mass2);
450  par.e = std::sqrt(pp2 + mass2);
451  par.theta2 = 14.1f * 14.1f / (beta2 * pp2 * 1e6f);
452  par.EP2 = par.e / pp2;
453 
454  // Approximate energy loss fluctuation (M.Ivanov)
455 
456  const double knst = 0.07f; // To be tuned.
457  par.sigmadE2 = knst * par.EP2 * qpt;
458  par.sigmadE2 = par.sigmadE2 * par.sigmadE2;
459 
460  par.k22 = (1.f + GetPar(3) * GetPar(3));
461  par.k33 = par.k22 * par.k22;
462  par.k43 = 0;
463  par.k44 = GetPar(3) * GetPar(3) * k2;
464 }
465 
466 bool GPUTPCTrackParam::CorrectForMeanMaterial(double xOverX0, double xTimesRho, const GPUTPCTrackFitParam& par)
467 {
468  //------------------------------------------------------------------
469  // This function corrects the track parameters for the crossed material.
470  // "xOverX0" - X/X0, the thickness in units of the radiation length.
471  // "xTimesRho" - is the product length*density (g/cm^2).
472  //------------------------------------------------------------------
473 
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];
481 
482  // Energy losses************************
483 
484  double dE = par.bethe * xTimesRho;
485  if (std::abs(dE) > 0.3f * par.e) {
486  return 0; // 30% energy loss is too much!
487  }
488  double corr = (1.f - par.EP2 * dE);
489  if (corr < 0.3f || corr > 1.3f) {
490  return 0;
491  }
492 
493  SetPar(4, GetPar(4) * corr);
494  mC40 *= corr;
495  mC41 *= corr;
496  mC42 *= corr;
497  mC43 *= corr;
498  mC44 *= corr * corr;
499  mC44 += par.sigmadE2 * std::abs(dE);
500 
501  // Multiple scattering******************
502 
503  double theta2 = par.theta2 * std::abs(xOverX0);
504  mC22 += theta2 * par.k22 * (1.f - GetPar(2)) * (1.f + GetPar(2));
505  mC33 += theta2 * par.k33;
506  mC43 += theta2 * par.k43;
507  mC44 += theta2 * par.k44;
508 
509  return 1;
510 }
511 
512 //*
513 //* Rotation
514 //*
515 bool GPUTPCTrackParam::Rotate(double alpha, double maxSinPhi)
516 {
517  //* Rotate the coordinate system in XY on the angle alpha
518 
519  double cA = cos(alpha);
520  double sA = sin(alpha);
521  double x = X(), y = Y(), sP = SinPhi(), cP = GetCosPhi();
522  double cosPhi = cP * cA + sP * sA;
523  double sinPhi = -cP * sA + sP * cA;
524 
525  if (std::abs(sinPhi) > maxSinPhi || std::abs(cosPhi) < 1.e-2f || std::abs(cP) < 1.e-2f) {
526  return 0;
527  }
528 
529  double j0 = cP / cosPhi;
530  double j2 = cosPhi / cP;
531 
532  SetX(x * cA + y * sA);
533  SetY(-x * sA + y * cA);
534  SetSignCosPhi(cosPhi);
535  SetSinPhi(sinPhi);
536 
537  // double J[5][5] = { { j0, 0, 0, 0, 0 }, // Y
538  // { 0, 1, 0, 0, 0 }, // Z
539  // { 0, 0, j2, 0, 0 }, // SinPhi
540  // { 0, 0, 0, 1, 0 }, // DzDs
541  // { 0, 0, 0, 0, 1 } }; // Kappa
542  // cout<<"alpha="<<alpha<<" "<<x<<" "<<y<<" "<<sP<<" "<<cP<<" "<<j0<<" "<<j2<<endl;
543  // cout<<" "<<mC[0]<<" "<<mC[1]<<" "<<mC[6]<<" "<<mC[10]<<" "<<mC[4]<<" "<<mC[5]<<" "<<mC[8]<<" "<<mC[12]<<endl;
544  mC[0] *= j0 * j0;
545  mC[1] *= j0;
546  mC[3] *= j0;
547  mC[6] *= j0;
548  mC[10] *= j0;
549 
550  mC[3] *= j2;
551  mC[4] *= j2;
552  mC[5] *= j2 * j2;
553  mC[8] *= j2;
554  mC[12] *= j2;
555 
556  if (cosPhi < 0) {
557  SetSinPhi(-SinPhi());
558  SetDzDs(-DzDs());
559  SetQPt(-QPt());
560  mC[3] = -mC[3];
561  mC[4] = -mC[4];
562  mC[6] = -mC[6];
563  mC[7] = -mC[7];
564  mC[10] = -mC[10];
565  mC[11] = -mC[11];
566  }
567 
568  // cout<<" "<<mC[0]<<" "<<mC[1]<<" "<<mC[6]<<" "<<mC[10]<<" "<<mC[4]<<" "<<mC[5]<<" "<<mC[8]<<" "<<mC[12]<<endl;
569  return 1;
570 }
571 
572 bool GPUTPCTrackParam::Rotate(double alpha, GPUTPCTrackLinearisation& t0, double maxSinPhi)
573 {
574  //* Rotate the coordinate system in XY on the angle alpha
575 
576  double cA = cos(alpha);
577  double sA = sin(alpha);
578  double x0 = X(), y0 = Y(), sP = t0.SinPhi(), cP = t0.CosPhi();
579  double cosPhi = cP * cA + sP * sA;
580  double sinPhi = -cP * sA + sP * cA;
581 
582  if (std::abs(sinPhi) > maxSinPhi || std::abs(cosPhi) < 1.e-2f || std::abs(cP) < 1.e-2f) {
583  return 0;
584  }
585 
586  // double J[5][5] = { { j0, 0, 0, 0, 0 }, // Y
587  // { 0, 1, 0, 0, 0 }, // Z
588  // { 0, 0, j2, 0, 0 }, // SinPhi
589  // { 0, 0, 0, 1, 0 }, // DzDs
590  // { 0, 0, 0, 0, 1 } }; // Kappa
591 
592  double j0 = cP / cosPhi;
593  double j2 = cosPhi / cP;
594  double d[2] = {Y() - y0, SinPhi() - sP};
595 
596  SetX(x0 * cA + y0 * sA);
597  SetY(-x0 * sA + y0 * cA + j0 * d[0]);
598  t0.SetCosPhi(cosPhi);
599  t0.SetSinPhi(sinPhi);
600 
601  SetSinPhi(sinPhi + j2 * d[1]);
602 
603  mC[0] *= j0 * j0;
604  mC[1] *= j0;
605  mC[3] *= j0;
606  mC[6] *= j0;
607  mC[10] *= j0;
608 
609  mC[3] *= j2;
610  mC[4] *= j2;
611  mC[5] *= j2 * j2;
612  mC[8] *= j2;
613  mC[12] *= j2;
614 
615  return 1;
616 }
617 
618 bool GPUTPCTrackParam::Filter(double y, double z, double err2Y, double err2Z, double maxSinPhi, bool paramOnly)
619 {
620  //* Add the y,z measurement with the Kalman filter
621 
622  double c00 = mC[0], c11 = mC[2], c20 = mC[3], c31 = mC[7], c40 = mC[10];
623 
624  err2Y += c00;
625  err2Z += c11;
626 
627  double z0 = y - GetPar(0), z1 = z - GetPar(1);
628 
629  if (err2Y < 1.e-8f || err2Z < 1.e-8f) {
630  return 0;
631  }
632 
633  double mS0 = 1.f / err2Y;
634  double mS2 = 1.f / err2Z;
635 
636  // K = CHtS
637 
638  double k00, k11, k20, k31, k40;
639 
640  k00 = c00 * mS0;
641  k20 = c20 * mS0;
642  k40 = c40 * mS0;
643 
644  k11 = c11 * mS2;
645  k31 = c31 * mS2;
646 
647  double sinPhi = GetPar(2) + k20 * z0;
648 
649  if (maxSinPhi > 0 && std::abs(sinPhi) >= maxSinPhi) {
650  return 0;
651  }
652 
653  SetPar(0, GetPar(0) + k00 * z0);
654  SetPar(1, GetPar(1) + k11 * z1);
655  SetPar(2, sinPhi);
656  SetPar(3, GetPar(3) + k31 * z1);
657  SetPar(4, GetPar(4) + k40 * z0);
658  if (paramOnly) {
659  return true;
660  }
661 
662  mNDF += 2;
663  mChi2 += mS0 * z0 * z0 + mS2 * z1 * z1;
664 
665  mC[0] -= k00 * c00;
666  mC[3] -= k20 * c00;
667  mC[5] -= k20 * c20;
668  mC[10] -= k40 * c00;
669  mC[12] -= k40 * c20;
670  mC[14] -= k40 * c40;
671 
672  mC[2] -= k11 * c11;
673  mC[7] -= k31 * c11;
674  mC[9] -= k31 * c31;
675 
676  return 1;
677 }
678 
680 {
681  //* Check that the track parameters and covariance matrix are reasonable
682 
684 
685  const double* c = Cov();
686  for (int i = 0; i < 15; i++) {
687  ok = ok && std::isfinite(c[i]);
688  }
689  for (int i = 0; i < 5; i++) {
690  ok = ok && std::isfinite(Par()[i]);
691  }
692 
693  if (c[0] <= 0 || c[2] <= 0 || c[5] <= 0 || c[9] <= 0 || c[14] <= 0) {
694  ok = 0;
695  }
696  if (c[0] > 5.f || c[2] > 5.f || c[5] > 2.f || c[9] > 2
697  //|| ( std::abs( QPt() ) > 1.e-2 && c[14] > 2. )
698  ) {
699  ok = 0;
700  }
701 
702  if (std::abs(SinPhi()) > GPUCA_MAX_SIN_PHI) {
703  ok = 0;
704  }
705  if (std::abs(QPt()) > 1.f / 0.05f) {
706  ok = 0;
707  }
708  if (ok) {
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]);
711  }
712  return ok;
713 }
714 
715 #if !defined(GPUCA_GPUCODE)
716 #include <iostream>
717 #endif
718 
720 {
721  //* print parameters
722 
723 #if !defined(GPUCA_GPUCODE)
724  std::cout << "track: x=" << GetX() << " c=" << GetSignCosPhi() << ", P= " << GetY() << " " << GetZ() << " " << GetSinPhi() << " " << GetDzDs() << " " << GetQPt() << std::endl;
725  std::cout << "errs2: " << GetErr2Y() << " " << GetErr2Z() << " " << GetErr2SinPhi() << " " << GetErr2DzDs() << " " << GetErr2QPt() << std::endl;
726 #endif
727 }