Project
Loading...
Searching...
No Matches
TrackParametrizationWithError.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
12#include "ReconstructionDataFormats/TrackParametrizationWithError.h"
13#include "ReconstructionDataFormats/Vertex.h"
14#include "ReconstructionDataFormats/DCA.h"
15#include "CommonConstants/MathConstants.h"
16#include <GPUCommonLogger.h>
17
18#ifndef GPUCA_GPUCODE_DEVICE
19#include <iostream>
20#endif
21
22#ifndef GPUCA_ALIGPUCODE
23#include <fmt/printf.h>
24#endif
25
26using namespace o2::track;
27using namespace o2::gpu;
28
29//______________________________________________________________
30template <typename value_T>
31GPUd() void TrackParametrizationWithError<value_T>::invert()
32{
33 // Transform this track to the local coord. system rotated by 180 deg.
34 this->invertParam();
35 // since the fP1 and fP2 are not inverted, their covariances with others change sign
36 mC[kSigZY] = -mC[kSigZY];
37 mC[kSigSnpY] = -mC[kSigSnpY];
38 mC[kSigTglZ] = -mC[kSigTglZ];
39 mC[kSigTglSnp] = -mC[kSigTglSnp];
40 mC[kSigQ2PtZ] = -mC[kSigQ2PtZ];
41 mC[kSigQ2PtSnp] = -mC[kSigQ2PtSnp];
42}
43
44//______________________________________________________________
45template <typename value_T>
46GPUd() bool TrackParametrizationWithError<value_T>::propagateTo(value_t xk, value_t bz)
47{
48 //----------------------------------------------------------------
49 // propagate this track to the plane X=xk (cm) in the field "b" (kG)
50 //----------------------------------------------------------------
51 value_t dx = xk - this->getX();
52 if (gpu::CAMath::Abs(dx) < constants::math::Almost0) {
53 return true;
54 }
55 value_t crv = this->getCurvature(bz);
56 value_t x2r = crv * dx;
57 value_t f1 = this->getSnp(), f2 = f1 + x2r;
58 if ((gpu::CAMath::Abs(f1) > constants::math::Almost1) || (gpu::CAMath::Abs(f2) > constants::math::Almost1)) {
59 return false;
60 }
61 value_t r1 = gpu::CAMath::Sqrt((1.f - f1) * (1.f + f1));
62 if (gpu::CAMath::Abs(r1) < constants::math::Almost0) {
63 return false;
64 }
65 value_t r2 = gpu::CAMath::Sqrt((1.f - f2) * (1.f + f2));
66 if (gpu::CAMath::Abs(r2) < constants::math::Almost0) {
67 return false;
68 }
69 double r1pr2Inv = 1. / (r1 + r2);
70 double dy2dx = (f1 + f2) * r1pr2Inv;
71 bool arcz = gpu::CAMath::Abs(x2r) > 0.05f;
72 params_t dP{0.f};
73 if (arcz) {
74 // for small dx/R the linear apporximation of the arc by the segment is OK,
75 // but at large dx/R the error is very large and leads to incorrect Z propagation
76 // angle traversed delta = 2*asin(dist_start_end / R / 2), hence the arc is: R*deltaPhi
77 // The dist_start_end is obtained from sqrt(dx^2+dy^2) = x/(r1+r2)*sqrt(2+f1*f2+r1*r2)
78 // double chord = dx*TMath::Sqrt(1+dy2dx*dy2dx); // distance from old position to new one
79 // double rot = 2*TMath::ASin(0.5*chord*crv); // angular difference seen from the circle center
80 // track1 += rot/crv*track3;
81 //
82 auto arg = r1 * f2 - r2 * f1;
83 if (gpu::CAMath::Abs(arg) > constants::math::Almost1) {
84 return false;
85 }
86 value_t rot = gpu::CAMath::ASin(arg); // more economic version from Yura.
87 if (f1 * f1 + f2 * f2 > 1.f && f1 * f2 < 0.f) { // special cases of large rotations or large abs angles
88 if (f2 > 0.f) {
89 rot = constants::math::PI - rot; //
90 } else {
91 rot = -constants::math::PI - rot;
92 }
93 }
94 dP[kZ] = this->getTgl() / crv * rot;
95 } else {
96 dP[kZ] = dx * (r2 + f2 * dy2dx) * this->getTgl();
97 }
98 this->setX(xk);
99 dP[kY] = dx * dy2dx;
100 dP[kSnp] = x2r;
101
102 this->updateParams(dP); // apply corrections
103
104 value_t &c00 = mC[kSigY2], &c10 = mC[kSigZY], &c11 = mC[kSigZ2], &c20 = mC[kSigSnpY], &c21 = mC[kSigSnpZ],
105 &c22 = mC[kSigSnp2], &c30 = mC[kSigTglY], &c31 = mC[kSigTglZ], &c32 = mC[kSigTglSnp], &c33 = mC[kSigTgl2],
106 &c40 = mC[kSigQ2PtY], &c41 = mC[kSigQ2PtZ], &c42 = mC[kSigQ2PtSnp], &c43 = mC[kSigQ2PtTgl],
107 &c44 = mC[kSigQ2Pt2];
108
109 // evaluate matrix in double prec.
110 value_t kb = bz * constants::math::B2C;
111 double r2inv = 1. / r2, r1inv = 1. / r1;
112 double dx2r1pr2 = dx * r1pr2Inv;
113
114 double hh = dx2r1pr2 * r2inv * (1. + r1 * r2 + f1 * f2), jj = dx * (dy2dx - f2 * r2inv);
115 double f02 = hh * r1inv;
116 double f04 = hh * dx2r1pr2 * kb;
117 double f24 = dx * kb; // x2r/mP[kQ2Pt];
118 double f12 = this->getTgl() * (f02 * f2 + jj);
119 double f13 = dx * (r2 + f2 * dy2dx);
120 double f14 = this->getTgl() * (f04 * f2 + jj * f24);
121
122 // b = C*ft
123 double b00 = f02 * c20 + f04 * c40, b01 = f12 * c20 + f14 * c40 + f13 * c30;
124 double b02 = f24 * c40;
125 double b10 = f02 * c21 + f04 * c41, b11 = f12 * c21 + f14 * c41 + f13 * c31;
126 double b12 = f24 * c41;
127 double b20 = f02 * c22 + f04 * c42, b21 = f12 * c22 + f14 * c42 + f13 * c32;
128 double b22 = f24 * c42;
129 double b40 = f02 * c42 + f04 * c44, b41 = f12 * c42 + f14 * c44 + f13 * c43;
130 double b42 = f24 * c44;
131 double b30 = f02 * c32 + f04 * c43, b31 = f12 * c32 + f14 * c43 + f13 * c33;
132 double b32 = f24 * c43;
133
134 // a = f*b = f*C*ft
135 double a00 = f02 * b20 + f04 * b40, a01 = f02 * b21 + f04 * b41, a02 = f02 * b22 + f04 * b42;
136 double a11 = f12 * b21 + f14 * b41 + f13 * b31, a12 = f12 * b22 + f14 * b42 + f13 * b32;
137 double a22 = f24 * b42;
138
139 // F*C*Ft = C + (b + bt + a)
140 c00 += b00 + b00 + a00;
141 c10 += b10 + b01 + a01;
142 c20 += b20 + b02 + a02;
143 c30 += b30;
144 c40 += b40;
145 c11 += b11 + b11 + a11;
146 c21 += b21 + b12 + a12;
147 c31 += b31;
148 c41 += b41;
149 c22 += b22 + b22 + a22;
150 c32 += b32;
151 c42 += b42;
152
153 checkCovariance();
154
155 return true;
156}
157
158//______________________________________________________________
159template <typename value_T>
160GPUd() bool TrackParametrizationWithError<value_T>::propagateTo(value_t xk, TrackParametrization<value_T>& linRef0, value_t bz)
161{
162 //----------------------------------------------------------------
163 // propagate this track to the plane X=xk (cm) in the field "b" (kG), using linRef as linearization point
164 //----------------------------------------------------------------
165 if (this->getAbsCharge() == 0) {
166 bz = 0;
167 }
168 value_t dx = xk - this->getX();
169 if (gpu::CAMath::Abs(dx) < constants::math::Almost0) {
170 this->setX(xk);
171 linRef0.setX(xk);
172 return true;
173 }
174 // propagate reference track
175 TrackParametrization<value_T> linRef1 = linRef0;
176 if (!linRef1.propagateTo(xk, bz)) {
177 return false;
178 }
179 value_t kb = bz * constants::math::B2C;
180 // evaluate in double prec.
181 double snpRef0 = linRef0.getSnp(), cspRef0 = gpu::CAMath::Sqrt((1 - snpRef0) * (1 + snpRef0));
182 double snpRef1 = linRef1.getSnp(), cspRef1 = gpu::CAMath::Sqrt((1 - snpRef1) * (1 + snpRef1));
183 double cspRef0Inv = 1 / cspRef0, cspRef1Inv = 1 / cspRef1, cc = cspRef0 + cspRef1, ccInv = 1 / cc, dy2dx = (snpRef0 + snpRef1) * ccInv;
184 double dxccInv = dx * ccInv, hh = dxccInv * cspRef1Inv * (1 + cspRef0 * cspRef1 + snpRef0 * snpRef1), jj = dx * (dy2dx - snpRef1 * cspRef1Inv);
185
186 double f02 = hh * cspRef0Inv;
187 double f04 = hh * dxccInv * kb;
188 double f24 = dx * kb;
189 double f12 = linRef0.getTgl() * (f02 * snpRef1 + jj);
190 double f13 = dx * (cspRef1 + snpRef1 * dy2dx); // dS
191 double f14 = linRef0.getTgl() * (f04 * snpRef1 + jj * f24);
192
193 // difference between the current and reference state
194 value_t diff[5];
195 for (int i = 0; i < 5; i++) {
196 diff[i] = this->getParam(i) - linRef0.getParam(i);
197 }
198 value_t snpUpd = snpRef1 + diff[kSnp] + f24 * diff[kQ2Pt];
199 if (gpu::CAMath::Abs(snpUpd) > constants::math::Almost1) {
200 return false;
201 }
202 linRef0 = linRef1; // update reference track
203 this->setX(xk);
204 this->setY(linRef1.getY() + diff[kY] + f02 * diff[kSnp] + f04 * diff[kQ2Pt]);
205 this->setZ(linRef1.getZ() + diff[kZ] + f13 * diff[kTgl] + f14 * diff[kQ2Pt]);
206 this->setSnp(snpUpd);
207 this->setTgl(linRef1.getTgl() + diff[kTgl]);
208 this->setQ2Pt(linRef1.getQ2Pt() + diff[kQ2Pt]);
209
210 value_t &c00 = mC[kSigY2], &c10 = mC[kSigZY], &c11 = mC[kSigZ2], &c20 = mC[kSigSnpY], &c21 = mC[kSigSnpZ],
211 &c22 = mC[kSigSnp2], &c30 = mC[kSigTglY], &c31 = mC[kSigTglZ], &c32 = mC[kSigTglSnp], &c33 = mC[kSigTgl2],
212 &c40 = mC[kSigQ2PtY], &c41 = mC[kSigQ2PtZ], &c42 = mC[kSigQ2PtSnp], &c43 = mC[kSigQ2PtTgl],
213 &c44 = mC[kSigQ2Pt2];
214
215 // b = C*ft
216 double b00 = f02 * c20 + f04 * c40, b01 = f12 * c20 + f14 * c40 + f13 * c30;
217 double b02 = f24 * c40;
218 double b10 = f02 * c21 + f04 * c41, b11 = f12 * c21 + f14 * c41 + f13 * c31;
219 double b12 = f24 * c41;
220 double b20 = f02 * c22 + f04 * c42, b21 = f12 * c22 + f14 * c42 + f13 * c32;
221 double b22 = f24 * c42;
222 double b40 = f02 * c42 + f04 * c44, b41 = f12 * c42 + f14 * c44 + f13 * c43;
223 double b42 = f24 * c44;
224 double b30 = f02 * c32 + f04 * c43, b31 = f12 * c32 + f14 * c43 + f13 * c33;
225 double b32 = f24 * c43;
226
227 // a = f*b = f*C*ft
228 double a00 = f02 * b20 + f04 * b40, a01 = f02 * b21 + f04 * b41, a02 = f02 * b22 + f04 * b42;
229 double a11 = f12 * b21 + f14 * b41 + f13 * b31, a12 = f12 * b22 + f14 * b42 + f13 * b32;
230 double a22 = f24 * b42;
231
232 // F*C*Ft = C + (b + bt + a)
233 c00 += b00 + b00 + a00;
234 c10 += b10 + b01 + a01;
235 c20 += b20 + b02 + a02;
236 c30 += b30;
237 c40 += b40;
238 c11 += b11 + b11 + a11;
239 c21 += b21 + b12 + a12;
240 c31 += b31;
241 c41 += b41;
242 c22 += b22 + b22 + a22;
243 c32 += b32;
244 c42 += b42;
245
246 checkCovariance();
247
248 return true;
249}
250
251//______________________________________________________________
252template <typename value_T>
253GPUd() bool TrackParametrizationWithError<value_T>::testRotate(value_t) const
254{
255 // no ops
256 return true;
257}
258
259//______________________________________________________________
260template <typename value_T>
261GPUd() bool TrackParametrizationWithError<value_T>::rotate(value_t alpha)
262{
263 // rotate to alpha frame
264 if (gpu::CAMath::Abs(this->getSnp()) > constants::math::Almost1) {
265 LOGP(debug, "Precondition is not satisfied: |sin(phi)|>1 ! {:f}", this->getSnp());
266 return false;
267 }
268 //
269 math_utils::detail::bringToPMPi<value_t>(alpha);
270 //
271 value_t ca = 0, sa = 0;
272 math_utils::detail::sincos(alpha - this->getAlpha(), sa, ca);
273 value_t snp = this->getSnp(), csp = gpu::CAMath::Sqrt((1.f - snp) * (1.f + snp)); // Improve precision
274 // RS: check if rotation does no invalidate track model (cos(local_phi)>=0, i.e. particle
275 // direction in local frame is along the X axis
276 if ((csp * ca + snp * sa) < 0) {
277 // LOGP(warning,"Rotation failed: local cos(phi) would become {:.2f}", csp * ca + snp * sa);
278 return false;
279 }
280 //
281
282 value_t updSnp = snp * ca - csp * sa;
283 if (gpu::CAMath::Abs(updSnp) > constants::math::Almost1) {
284 LOGP(debug, "Rotation failed: new snp {:.2f}", updSnp);
285 return false;
286 }
287 value_t xold = this->getX(), yold = this->getY();
288 this->setAlpha(alpha);
289 this->setX(xold * ca + yold * sa);
290 this->setY(-xold * sa + yold * ca);
291 this->setSnp(updSnp);
292
293 if (gpu::CAMath::Abs(csp) < constants::math::Almost0) {
294 LOGP(debug, "Too small cosine value {:f}", csp);
295 csp = constants::math::Almost0;
296 }
297
298 value_t rr = (ca + snp / csp * sa);
299
300 mC[kSigY2] *= (ca * ca);
301 mC[kSigZY] *= ca;
302 mC[kSigSnpY] *= ca * rr;
303 mC[kSigSnpZ] *= rr;
304 mC[kSigSnp2] *= rr * rr;
305 mC[kSigTglY] *= ca;
306 mC[kSigTglSnp] *= rr;
307 mC[kSigQ2PtY] *= ca;
308 mC[kSigQ2PtSnp] *= rr;
309
310 checkCovariance();
311 return true;
312}
313
314//______________________________________________________________
315template <typename value_T>
316GPUd() bool TrackParametrizationWithError<value_T>::rotate(value_t alpha, TrackParametrization<value_T>& linRef0, value_t bz)
317{
318 // RS: similar to int32_t GPUTPCGMPropagator::RotateToAlpha(float newAlpha), i.e. rotate the track to new frame alpha, using linRef as linearization point
319 // rotate to alpha frame the reference (linearization point) trackParam, then align the current track to it
320 if (gpu::CAMath::Abs(this->getSnp()) > constants::math::Almost1) {
321 LOGP(debug, "Precondition is not satisfied: |sin(phi)|>1 ! {:f}", this->getSnp());
322 return false;
323 }
324 //
325 math_utils::detail::bringToPMPi<value_t>(alpha);
326 //
327 value_t ca = 0, sa = 0;
328 TrackParametrization<value_T> linRef1 = linRef0;
329 // rotate the reference, adjusting alpha to +-pi, return precalculated cos and sin of alpha - alphaOld
330 if (!linRef1.rotateParam(alpha, ca, sa)) {
331 return false;
332 }
333
334 value_t trackX = this->getX() * ca + this->getY() * sa; // X of the rotated current track
335 if (!linRef1.propagateParamTo(trackX, bz)) {
336 return false;
337 }
338
339 // now rotate the current track
340 value_t snp = this->getSnp(), csp = gpu::CAMath::Sqrt((1.f - snp) * (1.f + snp)), updSnp = snp * ca - csp * sa;
341 if ((csp * ca + snp * sa) < 0 || gpu::CAMath::Abs(updSnp) > constants::math::Almost1) {
342 // LOGP(warning,"Rotation failed: local cos(phi) would become {:.2f}", csp * ca + snp * sa);
343 return false;
344 }
345 this->setY(-sa * this->getX() + ca * this->getY());
346 this->setX(trackX);
347 this->setSnp(updSnp);
348 this->setAlpha(alpha);
349
350 // rotate covariance, accounting for the extra error from the rotated X
351 value_t snpRef0 = linRef0.getSnp(), cspRef0 = gpu::CAMath::Sqrt((value_t(1) - snpRef0) * (value_t(1) + snpRef0)); // original reference
352 value_t snpRef1 = linRef1.getSnp(), cspRef1 = ca * cspRef0 + sa * snpRef0; // rotated reference
353 value_t rr = cspRef1 / cspRef0; // cos1_ref / cos0_ref
354
355 // "extra row" of the lower triangle of cov. matrix
356 value_t cXSigY = mC[kSigY2] * ca * sa;
357 value_t cXSigZ = mC[kSigZY] * sa;
358 value_t cXSigSnp = mC[kSigSnpY] * rr * sa;
359 value_t cXSigTgl = mC[kSigTglY] * sa;
360 value_t cXSigQ2Pt = mC[kSigQ2PtY] * sa;
361 value_t cSigX2 = mC[kSigY2] * sa * sa;
362
363 // plane rotation of existing cov matrix
364 mC[kSigY2] *= ca * ca;
365 mC[kSigZY] *= ca;
366 mC[kSigSnpY] *= ca * rr;
367 mC[kSigSnpZ] *= rr;
368 mC[kSigSnp2] *= rr * rr;
369 mC[kSigTglY] *= ca;
370 mC[kSigTglSnp] *= rr;
371 mC[kSigQ2PtY] *= ca;
372 mC[kSigQ2PtSnp] *= rr;
373
374 // transport covariance from pseudo 6x6 matrix to usual 5x5, Jacobian (trust to Sergey):
375 auto cspRef1Inv = value_t(1) / cspRef1;
376 auto j3 = -snpRef1 * cspRef1Inv; // -pYmod/pXmod = -tg_pho = -sin_phi_mod / cos_phi_mod
377 auto j4 = -linRef1.getTgl() * cspRef1Inv; // -pZmod/pXmod = -tgl_mod / cos_phi_mod
378 auto j5 = linRef1.getCurvature(bz);
379 // Y Z Sin DzDs q/p X
380 // { { 1, 0, 0, 0, 0, j3 }, // Y
381 // { 0, 1, 0, 0, 0, j4 }, // Z
382 // { 0, 0, 1, 0, 0, j5 }, // snp
383 // { 0, 0, 0, 1, 0, 0 }, // tgl
384 // { 0, 0, 0, 0, 1, 0 } }; // q/pt
385 auto hXSigY = cXSigY + cSigX2 * j3;
386 auto hXSigZ = cXSigZ + cSigX2 * j4;
387 auto hXSigSnp = cXSigSnp + cSigX2 * j5;
388
389 mC[kSigY2] += j3 * (cXSigY + hXSigY);
390 mC[kSigZ2] += j4 * (cXSigZ + hXSigZ);
391 mC[kSigSnpY] += cXSigSnp * j3 + hXSigY * j5;
392 mC[kSigSnp2] += j5 * (cXSigSnp + hXSigSnp);
393 mC[kSigTglZ] += cXSigTgl * j4;
394 mC[kSigQ2PtY] += cXSigQ2Pt * j3;
395 mC[kSigQ2PtSnp] += cXSigQ2Pt * j5;
396
397 mC[kSigZY] += cXSigZ * j3 + hXSigY * j4;
398 mC[kSigSnpZ] += cXSigSnp * j4 + hXSigZ * j5;
399 mC[kSigTglY] += cXSigTgl * j3;
400 mC[kSigTglSnp] += cXSigTgl * j5;
401 mC[kSigQ2PtZ] += cXSigQ2Pt * j4;
402
403 checkCovariance();
404 linRef0 = linRef1;
405
406 return true;
407}
408
409//_______________________________________________________________________
410template <typename value_T>
411GPUd() bool TrackParametrizationWithError<value_T>::propagateToDCA(const o2::dataformats::VertexBase& vtx, value_t b, o2::dataformats::DCA* dca, value_t maxD)
412{
413 // propagate track to DCA to the vertex
414 value_t sn, cs, alp = this->getAlpha();
415 o2::math_utils::detail::sincos(alp, sn, cs);
416 value_t x = this->getX(), y = this->getY(), snp = this->getSnp(), csp = gpu::CAMath::Sqrt((1.f - snp) * (1.f + snp));
417 value_t xv = vtx.getX() * cs + vtx.getY() * sn, yv = -vtx.getX() * sn + vtx.getY() * cs, zv = vtx.getZ();
418 x -= xv;
419 y -= yv;
420 // Estimate the impact parameter neglecting the track curvature
421 value_t d = gpu::CAMath::Abs(x * snp - y * csp);
422 if (d > maxD) {
423 if (dca) { // provide default DCA for failed propag
424 dca->set(o2::track::DefaultDCA, o2::track::DefaultDCA,
425 o2::track::DefaultDCACov, o2::track::DefaultDCACov, o2::track::DefaultDCACov);
426 }
427 return false;
428 }
429 value_t crv = this->getCurvature(b);
430 value_t tgfv = -(crv * x - snp) / (crv * y + csp);
431 sn = tgfv / gpu::CAMath::Sqrt(1.f + tgfv * tgfv);
432 cs = gpu::CAMath::Sqrt((1.f - sn) * (1.f + sn));
433 cs = (gpu::CAMath::Abs(tgfv) > constants::math::Almost0) ? sn / tgfv : constants::math::Almost1;
434
435 x = xv * cs + yv * sn;
436 yv = -xv * sn + yv * cs;
437 xv = x;
438
439 auto tmpT(*this); // operate on the copy to recover after the failure
440 alp += gpu::CAMath::ASin(sn);
441 if (!tmpT.rotate(alp) || !tmpT.propagateTo(xv, b)) {
442#if !defined(GPUCA_ALIGPUCODE)
443 LOG(debug) << "failed to propagate to alpha=" << alp << " X=" << xv << vtx << " | Track is: " << tmpT.asString();
444#endif
445 if (dca) { // provide default DCA for failed propag
446 dca->set(o2::track::DefaultDCA, o2::track::DefaultDCA,
447 o2::track::DefaultDCACov, o2::track::DefaultDCACov, o2::track::DefaultDCACov);
448 }
449 return false;
450 }
451 *this = tmpT;
452 if (dca) {
453 o2::math_utils::detail::sincos(alp, sn, cs);
454 auto s2ylocvtx = vtx.getSigmaX2() * sn * sn + vtx.getSigmaY2() * cs * cs - 2. * vtx.getSigmaXY() * cs * sn;
455 dca->set(this->getY() - yv, this->getZ() - zv, getSigmaY2() + s2ylocvtx, getSigmaZY(), getSigmaZ2() + vtx.getSigmaZ2());
456 }
457 return true;
458}
459
460//______________________________________________________________
461template <typename value_T>
462GPUd() TrackParametrizationWithError<value_T>::TrackParametrizationWithError(const dim3_t& xyz, const dim3_t& pxpypz,
463 const std::array<value_t, kLabCovMatSize>& cv, int charge, bool sectorAlpha, const PID pid)
464{
465 // construct track param and covariance from kinematics and lab errors
466 set(xyz, pxpypz, cv, charge, sectorAlpha, pid);
467}
468
469//______________________________________________________________
470template <typename value_T>
471GPUd() void TrackParametrizationWithError<value_T>::set(const dim3_t& xyz, const dim3_t& pxpypz,
472 const std::array<value_t, kLabCovMatSize>& cv, int charge, bool sectorAlpha, const PID pid)
473{
474 // set track param and covariance from kinematics and lab errors
475
476 // Alpha of the frame is defined as:
477 // sectorAlpha == false : -> angle of pt direction
478 // sectorAlpha == true : -> angle of the sector from X,Y coordinate for r>1
479 // angle of pt direction for r==0
480 //
481 //
482 constexpr value_t kSafe = 1e-5f;
483 value_t radPos2 = xyz[0] * xyz[0] + xyz[1] * xyz[1];
484 value_t alp = 0;
485 if (sectorAlpha || radPos2 < 1) {
486 alp = gpu::CAMath::ATan2(pxpypz[1], pxpypz[0]);
487 } else {
488 alp = gpu::CAMath::ATan2(xyz[1], xyz[0]);
489 }
490 if (sectorAlpha) {
491 alp = math_utils::detail::angle2Alpha<value_t>(alp);
492 }
493 //
494 value_t sn, cs;
495 math_utils::detail::sincos(alp, sn, cs);
496 // protection against cosp<0
497 if (cs * pxpypz[0] + sn * pxpypz[1] < 0) {
498 LOG(debug) << "alpha from phiPos() will invalidate this track parameters, overriding to alpha from phi()";
499 alp = gpu::CAMath::ATan2(pxpypz[1], pxpypz[0]);
500 if (sectorAlpha) {
501 alp = math_utils::detail::angle2Alpha<value_t>(alp);
502 }
503 math_utils::detail::sincos(alp, sn, cs);
504 }
505 // protection: avoid alpha being too close to 0 or +-pi/2
506 if (gpu::CAMath::Abs(sn) < 2.f * kSafe) {
507 if (alp > 0) {
508 alp += alp < constants::math::PIHalf ? 2.f * kSafe : -2.f * kSafe;
509 } else {
510 alp += alp > -constants::math::PIHalf ? -2.f * kSafe : 2.f * kSafe;
511 }
512 math_utils::detail::sincos(alp, sn, cs);
513 } else if (gpu::CAMath::Abs(cs) < 2.f * kSafe) {
514 if (alp > 0) {
515 alp += alp > constants::math::PIHalf ? 2.f * kSafe : -2.f * kSafe;
516 } else {
517 alp += alp > -constants::math::PIHalf ? 2.f * kSafe : -2.f * kSafe;
518 }
519 math_utils::detail::sincos(alp, sn, cs);
520 }
521 // get the vertex of origin and the momentum
522 dim3_t ver{xyz[0], xyz[1], xyz[2]};
523 dim3_t mom{pxpypz[0], pxpypz[1], pxpypz[2]};
524 //
525 // Rotate to the local coordinate system
526 math_utils::detail::rotateZ<value_t>(ver, -alp);
527 math_utils::detail::rotateZ<value_t>(mom, -alp);
528 //
529 value_t pt = gpu::CAMath::Sqrt(mom[0] * mom[0] + mom[1] * mom[1]);
530 value_t ptI = 1.f / pt;
531 this->setX(ver[0]);
532 this->setAlpha(alp);
533 this->setY(ver[1]);
534 this->setZ(ver[2]);
535 this->setSnp(mom[1] * ptI); // cos(phi)
536 this->setTgl(mom[2] * ptI); // tg(lambda)
537 this->setAbsCharge(gpu::CAMath::Abs(charge));
538 this->setQ2Pt(charge ? ptI * charge : ptI);
539 this->setPID(pid);
540 //
541 if (gpu::CAMath::Abs(1.f - this->getSnp()) < kSafe) {
542 this->setSnp(1.f - kSafe); // Protection
543 } else if (gpu::CAMath::Abs(-1.f - this->getSnp()) < kSafe) {
544 this->setSnp(-1.f + kSafe); // Protection
545 }
546 //
547 // Covariance matrix (formulas to be simplified)
548 value_t r = mom[0] * ptI; // cos(phi)
549 value_t cv34 = gpu::CAMath::Sqrt(cv[3] * cv[3] + cv[4] * cv[4]);
550 //
551 int special = 0;
552 value_t sgcheck = r * sn + this->getSnp() * cs;
553 if (gpu::CAMath::Abs(sgcheck) > 1 - kSafe) { // special case: lab phi is +-pi/2
554 special = 1;
555 sgcheck = sgcheck < 0 ? -1.f : 1.f;
556 } else if (gpu::CAMath::Abs(sgcheck) < kSafe) {
557 sgcheck = cs < 0 ? -1.0f : 1.0f;
558 special = 2; // special case: lab phi is 0
559 }
560 //
561 mC[kSigY2] = cv[0] + cv[2];
562 mC[kSigZY] = (-cv[3] * sn) < 0 ? -cv34 : cv34;
563 mC[kSigZ2] = cv[5];
564 //
565 value_t ptI2 = ptI * ptI;
566 value_t tgl2 = this->getTgl() * this->getTgl();
567 if (special == 1) {
568 mC[kSigSnpY] = cv[6] * ptI;
569 mC[kSigSnpZ] = -sgcheck * cv[8] * r * ptI;
570 mC[kSigSnp2] = gpu::CAMath::Abs(cv[9] * r * r * ptI2);
571 mC[kSigTglY] = (cv[10] * this->getTgl() - sgcheck * cv[15]) * ptI / r;
572 mC[kSigTglZ] = (cv[17] - sgcheck * cv[12] * this->getTgl()) * ptI;
573 mC[kSigTglSnp] = (-sgcheck * cv[18] + cv[13] * this->getTgl()) * r * ptI2;
574 mC[kSigTgl2] = gpu::CAMath::Abs(cv[20] - 2 * sgcheck * cv[19] * mC[4] + cv[14] * tgl2) * ptI2;
575 mC[kSigQ2PtY] = cv[10] * ptI2 / r * charge;
576 mC[kSigQ2PtZ] = -sgcheck * cv[12] * ptI2 * charge;
577 mC[kSigQ2PtSnp] = cv[13] * r * ptI * ptI2 * charge;
578 mC[kSigQ2PtTgl] = (-sgcheck * cv[19] + cv[14] * this->getTgl()) * r * ptI2 * ptI;
579 mC[kSigQ2Pt2] = gpu::CAMath::Abs(cv[14] * ptI2 * ptI2);
580 } else if (special == 2) {
581 mC[kSigSnpY] = -cv[10] * ptI * cs / sn;
582 mC[kSigSnpZ] = cv[12] * cs * ptI;
583 mC[kSigSnp2] = gpu::CAMath::Abs(cv[14] * cs * cs * ptI2);
584 mC[kSigTglY] = (sgcheck * cv[6] * this->getTgl() - cv[15]) * ptI / sn;
585 mC[kSigTglZ] = (cv[17] - sgcheck * cv[8] * this->getTgl()) * ptI;
586 mC[kSigTglSnp] = (cv[19] - sgcheck * cv[13] * this->getTgl()) * cs * ptI2;
587 mC[kSigTgl2] = gpu::CAMath::Abs(cv[20] - 2 * sgcheck * cv[18] * this->getTgl() + cv[9] * tgl2) * ptI2;
588 mC[kSigQ2PtY] = sgcheck * cv[6] * ptI2 / sn * charge;
589 mC[kSigQ2PtZ] = -sgcheck * cv[8] * ptI2 * charge;
590 mC[kSigQ2PtSnp] = -sgcheck * cv[13] * cs * ptI * ptI2 * charge;
591 mC[kSigQ2PtTgl] = (-sgcheck * cv[18] + cv[9] * this->getTgl()) * ptI2 * ptI * charge;
592 mC[kSigQ2Pt2] = gpu::CAMath::Abs(cv[9] * ptI2 * ptI2);
593 } else {
594 double m00 = -sn; // m10=cs;
595 double m23 = -pt * (sn + this->getSnp() * cs / r), m43 = -pt * pt * (r * cs - this->getSnp() * sn);
596 double m24 = pt * (cs - this->getSnp() * sn / r), m44 = -pt * pt * (r * sn + this->getSnp() * cs);
597 double m35 = pt, m45 = -pt * pt * this->getTgl();
598 //
599 if (charge) { // RS: this is a hack, proper treatment to be implemented
600 m43 *= charge;
601 m44 *= charge;
602 m45 *= charge;
603 }
604 //
605 double a1 = cv[13] - cv[9] * (m23 * m44 + m43 * m24) / m23 / m43;
606 double a2 = m23 * m24 - m23 * (m23 * m44 + m43 * m24) / m43;
607 double a3 = m43 * m44 - m43 * (m23 * m44 + m43 * m24) / m23;
608 double a4 = cv[14] + 2. * cv[9];
609 double a5 = m24 * m24 - 2. * m24 * m44 * m23 / m43;
610 double a6 = m44 * m44 - 2. * m24 * m44 * m43 / m23;
611 //
612 mC[kSigSnpY] = (cv[10] * m43 - cv[6] * m44) / (m24 * m43 - m23 * m44) / m00;
613 mC[kSigQ2PtY] = (cv[6] / m00 - mC[kSigSnpY] * m23) / m43;
614 mC[kSigTglY] = (cv[15] / m00 - mC[kSigQ2PtY] * m45) / m35;
615 mC[kSigSnpZ] = (cv[12] * m43 - cv[8] * m44) / (m24 * m43 - m23 * m44);
616 mC[kSigQ2PtZ] = (cv[8] - mC[kSigSnpZ] * m23) / m43;
617 mC[kSigTglZ] = cv[17] / m35 - mC[kSigQ2PtZ] * m45 / m35;
618 mC[kSigSnp2] = gpu::CAMath::Abs((a4 * a3 - a6 * a1) / (a5 * a3 - a6 * a2));
619 mC[kSigQ2Pt2] = gpu::CAMath::Abs((a1 - a2 * mC[kSigSnp2]) / a3);
620 mC[kSigQ2PtSnp] = (cv[9] - mC[kSigSnp2] * m23 * m23 - mC[kSigQ2Pt2] * m43 * m43) / m23 / m43;
621 double b1 = cv[18] - mC[kSigQ2PtSnp] * m23 * m45 - mC[kSigQ2Pt2] * m43 * m45;
622 double b2 = m23 * m35;
623 double b3 = m43 * m35;
624 double b4 = cv[19] - mC[kSigQ2PtSnp] * m24 * m45 - mC[kSigQ2Pt2] * m44 * m45;
625 double b5 = m24 * m35;
626 double b6 = m44 * m35;
627 mC[kSigTglSnp] = (b4 - b6 * b1 / b3) / (b5 - b6 * b2 / b3);
628 mC[kSigQ2PtTgl] = b1 / b3 - b2 * mC[kSigTglSnp] / b3;
629 mC[kSigTgl2] = gpu::CAMath::Abs((cv[20] - mC[kSigQ2Pt2] * (m45 * m45) - mC[kSigQ2PtTgl] * 2.f * m35 * m45) / (m35 * m35));
630 }
631 checkCovariance();
632}
633
634//____________________________________________________________
635template <typename value_T>
636GPUd() bool TrackParametrizationWithError<value_T>::propagateTo(value_t xk, const dim3_t& b)
637{
638 //----------------------------------------------------------------
639 // Extrapolate this track to the plane X=xk in the field b[].
640 //
641 // X [cm] is in the "tracking coordinate system" of this track.
642 // b[]={Bx,By,Bz} [kG] is in the Global coordidate system.
643 //----------------------------------------------------------------
644
645 value_t dx = xk - this->getX();
646 if (gpu::CAMath::Abs(dx) < constants::math::Almost0) {
647 return true;
648 }
649 // Do not propagate tracks outside the ALICE detector
650 if (gpu::CAMath::Abs(dx) > 1e5 || gpu::CAMath::Abs(this->getY()) > 1e5 || gpu::CAMath::Abs(this->getZ()) > 1e5) {
651 LOG(warning) << "Anomalous track, target X:" << xk;
652 // print();
653 return false;
654 }
655 value_t crv = (gpu::CAMath::Abs(b[2]) < constants::math::Almost0) ? 0.f : this->getCurvature(b[2]);
656 if (gpu::CAMath::Abs(crv) < constants::math::Almost0) {
657 return propagateTo(xk, 0.);
658 }
659 value_t x2r = crv * dx;
660 value_t f1 = this->getSnp(), f2 = f1 + x2r;
661 if ((gpu::CAMath::Abs(f1) > constants::math::Almost1) || (gpu::CAMath::Abs(f2) > constants::math::Almost1)) {
662 return false;
663 }
664 value_t r1 = gpu::CAMath::Sqrt((1.f - f1) * (1.f + f1));
665 if (gpu::CAMath::Abs(r1) < constants::math::Almost0) {
666 return false;
667 }
668 value_t r2 = gpu::CAMath::Sqrt((1.f - f2) * (1.f + f2));
669 if (gpu::CAMath::Abs(r2) < constants::math::Almost0) {
670 return false;
671 }
672 double r1pr2Inv = 1. / (r1 + r2), r2inv = 1. / r2, r1inv = 1. / r1;
673 double dy2dx = (f1 + f2) * r1pr2Inv, dx2r1pr2 = dx * r1pr2Inv;
674 value_t step = (gpu::CAMath::Abs(x2r) < 0.05f) ? dx * gpu::CAMath::Abs(r2 + f2 * dy2dx) // chord
675 : 2.f * gpu::CAMath::ASin(0.5f * dx * gpu::CAMath::Sqrt(1.f + dy2dx * dy2dx) * crv) / crv; // arc
676 step *= gpu::CAMath::Sqrt(1.f + this->getTgl() * this->getTgl());
677 //
678 // get the track x,y,z,px/p,py/p,pz/p,p,sinAlpha,cosAlpha in the Global System
679 std::array<value_t, 9> vecLab{0.f};
680 if (!this->getPosDirGlo(vecLab)) {
681 return false;
682 }
683 //
684 // matrix transformed with Bz component only
685 value_t &c00 = mC[kSigY2], &c10 = mC[kSigZY], &c11 = mC[kSigZ2], &c20 = mC[kSigSnpY], &c21 = mC[kSigSnpZ],
686 &c22 = mC[kSigSnp2], &c30 = mC[kSigTglY], &c31 = mC[kSigTglZ], &c32 = mC[kSigTglSnp], &c33 = mC[kSigTgl2],
687 &c40 = mC[kSigQ2PtY], &c41 = mC[kSigQ2PtZ], &c42 = mC[kSigQ2PtSnp], &c43 = mC[kSigQ2PtTgl],
688 &c44 = mC[kSigQ2Pt2];
689
690 // evaluate matrix in double prec.
691 value_t kb = b[2] * constants::math::B2C;
692 double hh = dx2r1pr2 * r2inv * (1. + r1 * r2 + f1 * f2), jj = dx * (dy2dx - f2 * r2inv);
693 double f02 = hh * r1inv;
694 double f04 = hh * dx2r1pr2 * kb;
695 double f24 = dx * kb; // x2r/mP[kQ2Pt];
696 double f12 = this->getTgl() * (f02 * f2 + jj);
697 double f13 = dx * (r2 + f2 * dy2dx);
698 double f14 = this->getTgl() * (f04 * f2 + jj * f24);
699
700 // b = C*ft
701 double b00 = f02 * c20 + f04 * c40, b01 = f12 * c20 + f14 * c40 + f13 * c30;
702 double b02 = f24 * c40;
703 double b10 = f02 * c21 + f04 * c41, b11 = f12 * c21 + f14 * c41 + f13 * c31;
704 double b12 = f24 * c41;
705 double b20 = f02 * c22 + f04 * c42, b21 = f12 * c22 + f14 * c42 + f13 * c32;
706 double b22 = f24 * c42;
707 double b40 = f02 * c42 + f04 * c44, b41 = f12 * c42 + f14 * c44 + f13 * c43;
708 double b42 = f24 * c44;
709 double b30 = f02 * c32 + f04 * c43, b31 = f12 * c32 + f14 * c43 + f13 * c33;
710 double b32 = f24 * c43;
711
712 // a = f*b = f*C*ft
713 double a00 = f02 * b20 + f04 * b40, a01 = f02 * b21 + f04 * b41, a02 = f02 * b22 + f04 * b42;
714 double a11 = f12 * b21 + f14 * b41 + f13 * b31, a12 = f12 * b22 + f14 * b42 + f13 * b32;
715 double a22 = f24 * b42;
716
717 // F*C*Ft = C + (b + bt + a)
718 c00 += b00 + b00 + a00;
719 c10 += b10 + b01 + a01;
720 c20 += b20 + b02 + a02;
721 c30 += b30;
722 c40 += b40;
723 c11 += b11 + b11 + a11;
724 c21 += b21 + b12 + a12;
725 c31 += b31;
726 c41 += b41;
727 c22 += b22 + b22 + a22;
728 c32 += b32;
729 c42 += b42;
730
731 checkCovariance();
732
733 // Rotate to the system where Bx=By=0.
734 value_t bxy2 = b[0] * b[0] + b[1] * b[1];
735 value_t bt = gpu::CAMath::Sqrt(bxy2);
736 value_t cosphi = 1.f, sinphi = 0.f;
737 if (bt > constants::math::Almost0) {
738 cosphi = b[0] / bt;
739 sinphi = b[1] / bt;
740 }
741 value_t bb = gpu::CAMath::Sqrt(bxy2 + b[2] * b[2]);
742 value_t costet = 1., sintet = 0.;
743 if (bb > constants::math::Almost0) {
744 costet = b[2] / bb;
745 sintet = bt / bb;
746 }
747 std::array<value_t, 7> vect{costet * cosphi * vecLab[0] + costet * sinphi * vecLab[1] - sintet * vecLab[2],
748 -sinphi * vecLab[0] + cosphi * vecLab[1],
749 sintet * cosphi * vecLab[0] + sintet * sinphi * vecLab[1] + costet * vecLab[2],
750 costet * cosphi * vecLab[3] + costet * sinphi * vecLab[4] - sintet * vecLab[5],
751 -sinphi * vecLab[3] + cosphi * vecLab[4],
752 sintet * cosphi * vecLab[3] + sintet * sinphi * vecLab[4] + costet * vecLab[5],
753 vecLab[6]};
754
755 // Do the helix step
756 value_t q = this->getCharge();
757 g3helx3(q * bb, step, vect);
758
759 // Rotate back to the Global System
760 vecLab[0] = cosphi * costet * vect[0] - sinphi * vect[1] + cosphi * sintet * vect[2];
761 vecLab[1] = sinphi * costet * vect[0] + cosphi * vect[1] + sinphi * sintet * vect[2];
762 vecLab[2] = -sintet * vect[0] + costet * vect[2];
763
764 vecLab[3] = cosphi * costet * vect[3] - sinphi * vect[4] + cosphi * sintet * vect[5];
765 vecLab[4] = sinphi * costet * vect[3] + cosphi * vect[4] + sinphi * sintet * vect[5];
766 vecLab[5] = -sintet * vect[3] + costet * vect[5];
767
768 // Rotate back to the Tracking System
769 value_t sinalp = -vecLab[7], cosalp = vecLab[8];
770 value_t t = cosalp * vecLab[0] - sinalp * vecLab[1];
771 vecLab[1] = sinalp * vecLab[0] + cosalp * vecLab[1];
772 vecLab[0] = t;
773 t = cosalp * vecLab[3] - sinalp * vecLab[4];
774 vecLab[4] = sinalp * vecLab[3] + cosalp * vecLab[4];
775 vecLab[3] = t;
776
777 // Do the final correcting step to the target plane (linear approximation)
778 value_t x = vecLab[0], y = vecLab[1], z = vecLab[2];
779 if (gpu::CAMath::Abs(dx) > constants::math::Almost0) {
780 if (gpu::CAMath::Abs(vecLab[3]) < constants::math::Almost0) {
781 return false;
782 }
783 dx = xk - vecLab[0];
784 x += dx;
785 y += vecLab[4] / vecLab[3] * dx;
786 z += vecLab[5] / vecLab[3] * dx;
787 }
788
789 // Calculate the track parameters
790 t = 1.f / gpu::CAMath::Sqrt(vecLab[3] * vecLab[3] + vecLab[4] * vecLab[4]);
791 this->setX(xk);
792 this->setY(y);
793 this->setZ(z);
794 this->setSnp(vecLab[4] * t);
795 this->setTgl(vecLab[5] * t);
796 this->setQ2Pt(q * t / vecLab[6]);
797
798 return true;
799}
800
801//____________________________________________________________
802template <typename value_T>
803GPUd() bool TrackParametrizationWithError<value_T>::propagateTo(value_t xk, TrackParametrization<value_T>& linRef0, const dim3_t& b)
804{
805 //----------------------------------------------------------------
806 // Extrapolate this track to the plane X=xk in the field b[].
807 //
808 // X [cm] is in the "tracking coordinate system" of this track.
809 // b[]={Bx,By,Bz} [kG] is in the Global coordidate system.
810 //----------------------------------------------------------------
811
812 value_t dx = xk - this->getX();
813 if (gpu::CAMath::Abs(dx) < constants::math::Almost0) {
814 return true;
815 }
816 // Do not propagate tracks outside the ALICE detector
817 if (gpu::CAMath::Abs(dx) > 1e5 || gpu::CAMath::Abs(this->getY()) > 1e5 || gpu::CAMath::Abs(this->getZ()) > 1e5) {
818 LOG(warning) << "Anomalous track, target X:" << xk;
819 // print();
820 return false;
821 }
822 if (gpu::CAMath::Abs(dx) < constants::math::Almost0) {
823 this->setX(xk);
824 linRef0.setX(xk);
825 return true;
826 }
827 // preliminary calculations to find the step size
828 value_t crv = (gpu::CAMath::Abs(b[2]) < constants::math::Almost0) ? 0.f : linRef0.getCurvature(b[2]);
829 if (gpu::CAMath::Abs(crv) < constants::math::Almost0) {
830 return propagateTo(xk, linRef0, 0.);
831 }
832 value_t kb = b[2] * constants::math::B2C, x2r = crv * dx;
833 // evaluate in double prec.
834 value_t snpRef0 = linRef0.getSnp(), snpRef1 = snpRef0 + x2r;
835 if ((gpu::CAMath::Abs(snpRef0) > constants::math::Almost1) || (gpu::CAMath::Abs(snpRef1) > constants::math::Almost1)) {
836 return false;
837 }
838 value_t cspRef0 = gpu::CAMath::Sqrt((1 - snpRef0) * (1 + snpRef0)), cspRef1 = gpu::CAMath::Sqrt((1 - snpRef1) * (1 + snpRef1));
839 if (gpu::CAMath::Abs(cspRef0) < constants::math::Almost0 || gpu::CAMath::Abs(cspRef1) < constants::math::Almost0) {
840 return false;
841 }
842 value_t cspRef0Inv = value_t(1) / cspRef0, cspRef1Inv = value_t(1) / cspRef1, cc = cspRef0 + cspRef1, ccInv = value_t(1) / cc, dy2dx = (snpRef0 + snpRef1) * ccInv;
843 value_t step = (gpu::CAMath::Abs(crv * dx) < 0.05f) ? dx * (cspRef1 + snpRef1 * dy2dx) : 2. * gpu::CAMath::ASin(0.5 * dx * gpu::CAMath::Sqrt(1.f + dy2dx * dy2dx) * crv) / crv; // arc
844 step *= gpu::CAMath::Sqrt(1.f + linRef0.getTgl() * linRef0.getTgl());
845
846 //
847 // get the track x,y,z,px/p,py/p,pz/p,p,sinAlpha,cosAlpha in the Global System
848 std::array<value_t, 9> vecLab{0.f};
849 if (!linRef0.getPosDirGlo(vecLab)) {
850 return false;
851 }
852 //
853 // Rotate to the system where Bx=By=0.
854 value_t bxy2 = b[0] * b[0] + b[1] * b[1];
855 value_t bt = gpu::CAMath::Sqrt(bxy2);
856 value_t cosphi = 1.f, sinphi = 0.f;
857 if (bt > constants::math::Almost0) {
858 cosphi = b[0] / bt;
859 sinphi = b[1] / bt;
860 }
861 value_t bb = gpu::CAMath::Sqrt(bxy2 + b[2] * b[2]);
862 value_t costet = 1., sintet = 0.;
863 if (bb > constants::math::Almost0) {
864 costet = b[2] / bb;
865 sintet = bt / bb;
866 }
867 std::array<value_t, 7> vect{costet * cosphi * vecLab[0] + costet * sinphi * vecLab[1] - sintet * vecLab[2],
868 -sinphi * vecLab[0] + cosphi * vecLab[1],
869 sintet * cosphi * vecLab[0] + sintet * sinphi * vecLab[1] + costet * vecLab[2],
870 costet * cosphi * vecLab[3] + costet * sinphi * vecLab[4] - sintet * vecLab[5],
871 -sinphi * vecLab[3] + cosphi * vecLab[4],
872 sintet * cosphi * vecLab[3] + sintet * sinphi * vecLab[4] + costet * vecLab[5],
873 vecLab[6]};
874
875 // Do the helix step
876 value_t q = this->getCharge();
877 g3helx3(q * bb, step, vect);
878
879 // Rotate back to the Global System
880 vecLab[0] = cosphi * costet * vect[0] - sinphi * vect[1] + cosphi * sintet * vect[2];
881 vecLab[1] = sinphi * costet * vect[0] + cosphi * vect[1] + sinphi * sintet * vect[2];
882 vecLab[2] = -sintet * vect[0] + costet * vect[2];
883
884 vecLab[3] = cosphi * costet * vect[3] - sinphi * vect[4] + cosphi * sintet * vect[5];
885 vecLab[4] = sinphi * costet * vect[3] + cosphi * vect[4] + sinphi * sintet * vect[5];
886 vecLab[5] = -sintet * vect[3] + costet * vect[5];
887
888 // Rotate back to the Tracking System
889 value_t sinalp = -vecLab[7], cosalp = vecLab[8];
890 value_t t = cosalp * vecLab[0] - sinalp * vecLab[1];
891 vecLab[1] = sinalp * vecLab[0] + cosalp * vecLab[1];
892 vecLab[0] = t;
893 t = cosalp * vecLab[3] - sinalp * vecLab[4];
894 vecLab[4] = sinalp * vecLab[3] + cosalp * vecLab[4];
895 vecLab[3] = t;
896
897 // Do the final correcting step to the target plane (linear approximation)
898 value_t x = vecLab[0], y = vecLab[1], z = vecLab[2];
899 if (gpu::CAMath::Abs(dx) > constants::math::Almost0) {
900 if (gpu::CAMath::Abs(vecLab[3]) < constants::math::Almost0) {
901 return false;
902 }
903 auto dxFin = xk - vecLab[0];
904 x += dxFin;
905 y += vecLab[4] / vecLab[3] * dxFin;
906 z += vecLab[5] / vecLab[3] * dxFin;
907 }
908
909 // Calculate the track parameters
910 auto linRef1 = linRef0;
911 t = 1.f / gpu::CAMath::Sqrt(vecLab[3] * vecLab[3] + vecLab[4] * vecLab[4]);
912 linRef1.setX(xk);
913 linRef1.setY(y);
914 linRef1.setZ(z);
915 linRef1.setSnp(snpRef1 = vecLab[4] * t); // reassign snpRef1
916 linRef1.setTgl(vecLab[5] * t);
917 linRef1.setQ2Pt(q * t / vecLab[6]);
918
919 // recalculate parameters of the transported ref track needed for transport of this:
920 cspRef1 = gpu::CAMath::Sqrt((1 - snpRef1) * (1 + snpRef1));
921 cspRef1Inv = value_t(1) / cspRef1;
922 cc = cspRef0 + cspRef1;
923 ccInv = value_t(1) / cc;
924 dy2dx = (snpRef0 + snpRef1) * ccInv;
925 double dxccInv = dx * ccInv, hh = dxccInv * cspRef1Inv * (1 + cspRef0 * cspRef1 + snpRef0 * snpRef1), jj = dx * (dy2dx - snpRef1 * cspRef1Inv);
926 double f02 = hh * cspRef0Inv;
927 double f04 = hh * dxccInv * kb;
928 double f24 = dx * kb;
929 double f12 = linRef0.getTgl() * (f02 * snpRef1 + jj);
930 double f13 = dx * (cspRef1 + snpRef1 * dy2dx); // dS
931 double f14 = linRef0.getTgl() * (f04 * snpRef1 + jj * f24);
932
933 // difference between the current and reference state
934 value_t diff[5];
935 for (int i = 0; i < 5; i++) {
936 diff[i] = this->getParam(i) - linRef0.getParam(i);
937 }
938 value_t snpUpd = snpRef1 + diff[kSnp] + f24 * diff[kQ2Pt];
939 if (gpu::CAMath::Abs(snpUpd) > constants::math::Almost1) {
940 return false;
941 }
942 this->setX(xk);
943 this->setY(linRef1.getY() + diff[kY] + f02 * diff[kSnp] + f04 * diff[kQ2Pt]);
944 this->setZ(linRef1.getZ() + diff[kZ] + f13 * diff[kTgl] + f14 * diff[kQ2Pt]);
945 this->setSnp(snpUpd);
946 this->setTgl(linRef1.getTgl() + diff[kTgl]);
947 this->setQ2Pt(linRef1.getQ2Pt() + diff[kQ2Pt]);
948
949 linRef0 = linRef1; // update reference track
950
951 // matrix transformed with Bz component only
952 value_t &c00 = mC[kSigY2], &c10 = mC[kSigZY], &c11 = mC[kSigZ2], &c20 = mC[kSigSnpY], &c21 = mC[kSigSnpZ],
953 &c22 = mC[kSigSnp2], &c30 = mC[kSigTglY], &c31 = mC[kSigTglZ], &c32 = mC[kSigTglSnp], &c33 = mC[kSigTgl2],
954 &c40 = mC[kSigQ2PtY], &c41 = mC[kSigQ2PtZ], &c42 = mC[kSigQ2PtSnp], &c43 = mC[kSigQ2PtTgl],
955 &c44 = mC[kSigQ2Pt2];
956
957 // b = C*ft
958 double b00 = f02 * c20 + f04 * c40, b01 = f12 * c20 + f14 * c40 + f13 * c30;
959 double b02 = f24 * c40;
960 double b10 = f02 * c21 + f04 * c41, b11 = f12 * c21 + f14 * c41 + f13 * c31;
961 double b12 = f24 * c41;
962 double b20 = f02 * c22 + f04 * c42, b21 = f12 * c22 + f14 * c42 + f13 * c32;
963 double b22 = f24 * c42;
964 double b40 = f02 * c42 + f04 * c44, b41 = f12 * c42 + f14 * c44 + f13 * c43;
965 double b42 = f24 * c44;
966 double b30 = f02 * c32 + f04 * c43, b31 = f12 * c32 + f14 * c43 + f13 * c33;
967 double b32 = f24 * c43;
968
969 // a = f*b = f*C*ft
970 double a00 = f02 * b20 + f04 * b40, a01 = f02 * b21 + f04 * b41, a02 = f02 * b22 + f04 * b42;
971 double a11 = f12 * b21 + f14 * b41 + f13 * b31, a12 = f12 * b22 + f14 * b42 + f13 * b32;
972 double a22 = f24 * b42;
973
974 // F*C*Ft = C + (b + bt + a)
975 c00 += b00 + b00 + a00;
976 c10 += b10 + b01 + a01;
977 c20 += b20 + b02 + a02;
978 c30 += b30;
979 c40 += b40;
980 c11 += b11 + b11 + a11;
981 c21 += b21 + b12 + a12;
982 c31 += b31;
983 c41 += b41;
984 c22 += b22 + b22 + a22;
985 c32 += b32;
986 c42 += b42;
987
988 checkCovariance();
989
990 return true;
991}
992
993//______________________________________________
994template <typename value_T>
995GPUd() void TrackParametrizationWithError<value_T>::checkCorrelations()
996{
997 // This function forces the abs of correlation coefficients to be <1.
998 constexpr float MaxCorr = 0.99;
999 for (int i = kNParams; i--;) {
1000 for (int j = i; j--;) {
1001 auto sig2 = mC[DiagMap[i]] * mC[DiagMap[j]];
1002 auto& cov = mC[CovarMap[i][j]];
1003 if (cov * cov >= MaxCorr * sig2) { // constrain correlation
1004 cov = gpu::CAMath::Sqrt(sig2) * (cov > 0. ? MaxCorr : -MaxCorr);
1005 }
1006 }
1007 }
1008}
1009
1010//______________________________________________
1011template <typename value_T>
1012GPUd() void TrackParametrizationWithError<value_T>::checkCovariance()
1013{
1014 // This function forces the diagonal elements of the covariance matrix to be positive and abs of correlation coefficients to be <1.
1015 // In case the diagonal element is bigger than the maximal allowed value, it is set to
1016 // the limit and the off-diagonal elements that correspond to it are set to zero.
1017
1018 mC[kSigY2] = gpu::CAMath::Abs(mC[kSigY2]);
1019 if (mC[kSigY2] > kCY2max) {
1020 value_t scl = gpu::CAMath::Sqrt(kCY2max / mC[kSigY2]);
1021 mC[kSigY2] = kCY2max;
1022 mC[kSigZY] *= scl;
1023 mC[kSigSnpY] *= scl;
1024 mC[kSigTglY] *= scl;
1025 mC[kSigQ2PtY] *= scl;
1026 }
1027 mC[kSigZ2] = gpu::CAMath::Abs(mC[kSigZ2]);
1028 if (mC[kSigZ2] > kCZ2max) {
1029 value_t scl = gpu::CAMath::Sqrt(kCZ2max / mC[kSigZ2]);
1030 mC[kSigZ2] = kCZ2max;
1031 mC[kSigZY] *= scl;
1032 mC[kSigSnpZ] *= scl;
1033 mC[kSigTglZ] *= scl;
1034 mC[kSigQ2PtZ] *= scl;
1035 }
1036 mC[kSigSnp2] = gpu::CAMath::Abs(mC[kSigSnp2]);
1037 if (mC[kSigSnp2] > kCSnp2max) {
1038 value_t scl = gpu::CAMath::Sqrt(kCSnp2max / mC[kSigSnp2]);
1039 mC[kSigSnp2] = kCSnp2max;
1040 mC[kSigSnpY] *= scl;
1041 mC[kSigSnpZ] *= scl;
1042 mC[kSigTglSnp] *= scl;
1043 mC[kSigQ2PtSnp] *= scl;
1044 }
1045 mC[kSigTgl2] = gpu::CAMath::Abs(mC[kSigTgl2]);
1046 if (mC[kSigTgl2] > kCTgl2max) {
1047 value_t scl = gpu::CAMath::Sqrt(kCTgl2max / mC[kSigTgl2]);
1048 mC[kSigTgl2] = kCTgl2max;
1049 mC[kSigTglY] *= scl;
1050 mC[kSigTglZ] *= scl;
1051 mC[kSigTglSnp] *= scl;
1052 mC[kSigQ2PtTgl] *= scl;
1053 }
1054 mC[kSigQ2Pt2] = gpu::CAMath::Abs(mC[kSigQ2Pt2]);
1055 if (mC[kSigQ2Pt2] > kC1Pt2max) {
1056 value_t scl = gpu::CAMath::Sqrt(kC1Pt2max / mC[kSigQ2Pt2]);
1057 mC[kSigQ2Pt2] = kC1Pt2max;
1058 mC[kSigQ2PtY] *= scl;
1059 mC[kSigQ2PtZ] *= scl;
1060 mC[kSigQ2PtSnp] *= scl;
1061 mC[kSigQ2PtTgl] *= scl;
1062 }
1063}
1064
1065//______________________________________________
1066template <typename value_T>
1067GPUd() void TrackParametrizationWithError<value_T>::resetCovariance(value_t s2)
1068{
1069 // Reset the covarince matrix to "something big"
1070 double d0(kCY2max), d1(kCZ2max), d2(kCSnp2max), d3(kCTgl2max), d4(kC1Pt2max);
1071 if (s2 > constants::math::Almost0) {
1072 d0 = getSigmaY2() * s2;
1073 d1 = getSigmaZ2() * s2;
1074 d2 = getSigmaSnp2() * s2;
1075 d3 = getSigmaTgl2() * s2;
1076 d4 = getSigma1Pt2() * s2;
1077 if (d0 > kCY2max) {
1078 d0 = kCY2max;
1079 }
1080 if (d1 > kCZ2max) {
1081 d1 = kCZ2max;
1082 }
1083 if (d2 > kCSnp2max) {
1084 d2 = kCSnp2max;
1085 }
1086 if (d3 > kCTgl2max) {
1087 d3 = kCTgl2max;
1088 }
1089 if (d4 > kC1Pt2max) {
1090 d4 = kC1Pt2max;
1091 }
1092 }
1093 for (int i = 0; i < kCovMatSize; i++) {
1094 mC[i] = 0;
1095 }
1096 mC[kSigY2] = d0;
1097 mC[kSigZ2] = d1;
1098 mC[kSigSnp2] = d2;
1099 mC[kSigTgl2] = d3;
1100 mC[kSigQ2Pt2] = d4;
1101}
1102
1103//______________________________________________
1104template <typename value_T>
1105GPUd() auto TrackParametrizationWithError<value_T>::getPredictedChi2(const value_t* p, const value_t* cov) const -> value_t
1106{
1107 // Estimate the chi2 of the space point "p" with the cov. matrix "cov"
1108 auto sdd = static_cast<double>(getSigmaY2()) + static_cast<double>(cov[0]);
1109 auto sdz = static_cast<double>(getSigmaZY()) + static_cast<double>(cov[1]);
1110 auto szz = static_cast<double>(getSigmaZ2()) + static_cast<double>(cov[2]);
1111 auto det = sdd * szz - sdz * sdz;
1112
1113 if (gpu::CAMath::Abs(det) < constants::math::Almost0) {
1114 return constants::math::VeryBig;
1115 }
1116
1117 value_t d = this->getY() - p[0];
1118 value_t z = this->getZ() - p[1];
1119 auto chi2 = (d * (szz * d - sdz * z) + z * (sdd * z - d * sdz)) / det;
1120 if (chi2 < 0.) {
1121#ifndef GPUCA_ALIGPUCODE
1122 LOGP(warning, "Negative chi2={}, Cluster: {} {} {} Dy:{} Dz:{} | sdd:{} sdz:{} szz:{} det:{}", chi2, cov[0], cov[1], cov[2], d, z, sdd, sdz, szz, det);
1123 LOGP(warning, "Track: {}", asString());
1124#endif
1125 }
1126 return chi2;
1127}
1128
1129//______________________________________________
1130template <typename value_T>
1131GPUd() auto TrackParametrizationWithError<value_T>::getPredictedChi2Quiet(const value_t* p, const value_t* cov) const -> value_t
1132{
1133 // Estimate the chi2 of the space point "p" with the cov. matrix "cov"
1134 auto sdd = static_cast<double>(getSigmaY2()) + static_cast<double>(cov[0]);
1135 auto sdz = static_cast<double>(getSigmaZY()) + static_cast<double>(cov[1]);
1136 auto szz = static_cast<double>(getSigmaZ2()) + static_cast<double>(cov[2]);
1137 auto det = sdd * szz - sdz * sdz;
1138
1139 if (gpu::CAMath::Abs(det) < constants::math::Almost0) {
1140 return constants::math::VeryBig;
1141 }
1142
1143 value_t d = this->getY() - p[0];
1144 value_t z = this->getZ() - p[1];
1145
1146 return (d * (szz * d - sdz * z) + z * (sdd * z - d * sdz)) / det;
1147}
1148
1149//______________________________________________
1150template <typename value_T>
1151GPUd() auto TrackParametrizationWithError<value_T>::getPredictedChi2(const TrackParametrizationWithError<value_T>& rhs) const -> value_t
1152{
1153 MatrixDSym5 cov; // perform matrix operations in double!
1154 return getPredictedChi2(rhs, cov);
1155}
1156
1157//______________________________________________
1158template <typename value_T>
1159GPUd() void TrackParametrizationWithError<value_T>::buildCombinedCovMatrix(const TrackParametrizationWithError<value_T>& rhs, MatrixDSym5& cov) const
1160{
1161 // fill combined cov.matrix (NOT inverted)
1162 cov(kY, kY) = static_cast<double>(getSigmaY2()) + static_cast<double>(rhs.getSigmaY2());
1163 cov(kZ, kY) = static_cast<double>(getSigmaZY()) + static_cast<double>(rhs.getSigmaZY());
1164 cov(kZ, kZ) = static_cast<double>(getSigmaZ2()) + static_cast<double>(rhs.getSigmaZ2());
1165 cov(kSnp, kY) = static_cast<double>(getSigmaSnpY()) + static_cast<double>(rhs.getSigmaSnpY());
1166 cov(kSnp, kZ) = static_cast<double>(getSigmaSnpZ()) + static_cast<double>(rhs.getSigmaSnpZ());
1167 cov(kSnp, kSnp) = static_cast<double>(getSigmaSnp2()) + static_cast<double>(rhs.getSigmaSnp2());
1168 cov(kTgl, kY) = static_cast<double>(getSigmaTglY()) + static_cast<double>(rhs.getSigmaTglY());
1169 cov(kTgl, kZ) = static_cast<double>(getSigmaTglZ()) + static_cast<double>(rhs.getSigmaTglZ());
1170 cov(kTgl, kSnp) = static_cast<double>(getSigmaTglSnp()) + static_cast<double>(rhs.getSigmaTglSnp());
1171 cov(kTgl, kTgl) = static_cast<double>(getSigmaTgl2()) + static_cast<double>(rhs.getSigmaTgl2());
1172 cov(kQ2Pt, kY) = static_cast<double>(getSigma1PtY()) + static_cast<double>(rhs.getSigma1PtY());
1173 cov(kQ2Pt, kZ) = static_cast<double>(getSigma1PtZ()) + static_cast<double>(rhs.getSigma1PtZ());
1174 cov(kQ2Pt, kSnp) = static_cast<double>(getSigma1PtSnp()) + static_cast<double>(rhs.getSigma1PtSnp());
1175 cov(kQ2Pt, kTgl) = static_cast<double>(getSigma1PtTgl()) + static_cast<double>(rhs.getSigma1PtTgl());
1176 cov(kQ2Pt, kQ2Pt) = static_cast<double>(getSigma1Pt2()) + static_cast<double>(rhs.getSigma1Pt2());
1177}
1178
1179//______________________________________________
1180template <typename value_T>
1181GPUd() auto TrackParametrizationWithError<value_T>::getPredictedChi2(const TrackParametrizationWithError<value_T>& rhs, MatrixDSym5& covToSet) const -> value_t
1182{
1183 // get chi2 wrt other track, which must be defined at the same parameters X,alpha
1184 // Supplied non-initialized covToSet matrix is filled by inverse combined matrix for further use
1185
1186 if (gpu::CAMath::Abs(this->getAlpha() - rhs.getAlpha()) > o2::constants::math::Epsilon) {
1187 LOG(error) << "The reference Alpha of the tracks differ: " << this->getAlpha() << " : " << rhs.getAlpha();
1188 return 2.f * HugeF;
1189 }
1190 if (gpu::CAMath::Abs(this->getX() - rhs.getX()) > o2::constants::math::Epsilon) {
1191 LOG(error) << "The reference X of the tracks differ: " << this->getX() << " : " << rhs.getX();
1192 return 2.f * HugeF;
1193 }
1194 buildCombinedCovMatrix(rhs, covToSet);
1195 if (!covToSet.Invert()) {
1196 LOG(warning) << "Cov.matrix inversion failed: " << covToSet;
1197 return 2.f * HugeF;
1198 }
1199 double chi2diag = 0., chi2ndiag = 0., diff[kNParams];
1200 for (int i = kNParams; i--;) {
1201 diff[i] = this->getParam(i) - rhs.getParam(i);
1202 chi2diag += diff[i] * diff[i] * covToSet(i, i);
1203 }
1204 for (int i = kNParams; i--;) {
1205 for (int j = i; j--;) {
1206 chi2ndiag += diff[i] * diff[j] * covToSet(i, j);
1207 }
1208 }
1209 return chi2diag + 2. * chi2ndiag;
1210}
1211
1212//______________________________________________
1213template <typename value_T>
1214GPUd() bool TrackParametrizationWithError<value_T>::update(const TrackParametrizationWithError<value_T>& rhs, const MatrixDSym5& covInv)
1215{
1216 // update track with other track, the inverted combined cov matrix should be supplied
1217
1218 // consider skipping this check, since it is usually already done upstream
1219 if (gpu::CAMath::Abs(this->getAlpha() - rhs.getAlpha()) > o2::constants::math::Epsilon) {
1220 LOG(error) << "The reference Alpha of the tracks differ: " << this->getAlpha() << " : " << rhs.getAlpha();
1221 return false;
1222 }
1223 if (gpu::CAMath::Abs(this->getX() - rhs.getX()) > o2::constants::math::Epsilon) {
1224 LOG(error) << "The reference X of the tracks differ: " << this->getX() << " : " << rhs.getX();
1225 return false;
1226 }
1227
1228 // gain matrix K = Cov0*H*(Cov0+Cov0)^-1 (for measurement matrix H=I)
1229 MatrixDSym5 matC0;
1230 matC0(kY, kY) = getSigmaY2();
1231 matC0(kZ, kY) = getSigmaZY();
1232 matC0(kZ, kZ) = getSigmaZ2();
1233 matC0(kSnp, kY) = getSigmaSnpY();
1234 matC0(kSnp, kZ) = getSigmaSnpZ();
1235 matC0(kSnp, kSnp) = getSigmaSnp2();
1236 matC0(kTgl, kY) = getSigmaTglY();
1237 matC0(kTgl, kZ) = getSigmaTglZ();
1238 matC0(kTgl, kSnp) = getSigmaTglSnp();
1239 matC0(kTgl, kTgl) = getSigmaTgl2();
1240 matC0(kQ2Pt, kY) = getSigma1PtY();
1241 matC0(kQ2Pt, kZ) = getSigma1PtZ();
1242 matC0(kQ2Pt, kSnp) = getSigma1PtSnp();
1243 matC0(kQ2Pt, kTgl) = getSigma1PtTgl();
1244 matC0(kQ2Pt, kQ2Pt) = getSigma1Pt2();
1245 MatrixD5 matK = matC0 * covInv;
1246
1247 // updated state vector: x = K*(x1-x0)
1248 // RS: why SMatix, SVector does not provide multiplication operators ???
1249 double diff[kNParams];
1250 for (int i = kNParams; i--;) {
1251 diff[i] = rhs.getParam(i) - this->getParam(i);
1252 }
1253 for (int i = kNParams; i--;) {
1254 for (int j = kNParams; j--;) {
1255 this->updateParam(matK(i, j) * diff[j], i);
1256 }
1257 }
1258
1259 // updated covariance: Cov0 = Cov0 - K*Cov0
1260 matK *= o2::math_utils::SMatrix<double, kNParams, kNParams, o2::math_utils::MatRepStd<double, kNParams>>(matC0);
1261 mC[kSigY2] -= matK(kY, kY);
1262 mC[kSigZY] -= matK(kZ, kY);
1263 mC[kSigZ2] -= matK(kZ, kZ);
1264 mC[kSigSnpY] -= matK(kSnp, kY);
1265 mC[kSigSnpZ] -= matK(kSnp, kZ);
1266 mC[kSigSnp2] -= matK(kSnp, kSnp);
1267 mC[kSigTglY] -= matK(kTgl, kY);
1268 mC[kSigTglZ] -= matK(kTgl, kZ);
1269 mC[kSigTglSnp] -= matK(kTgl, kSnp);
1270 mC[kSigTgl2] -= matK(kTgl, kTgl);
1271 mC[kSigQ2PtY] -= matK(kQ2Pt, kY);
1272 mC[kSigQ2PtZ] -= matK(kQ2Pt, kZ);
1273 mC[kSigQ2PtSnp] -= matK(kQ2Pt, kSnp);
1274 mC[kSigQ2PtTgl] -= matK(kQ2Pt, kTgl);
1275 mC[kSigQ2Pt2] -= matK(kQ2Pt, kQ2Pt);
1276
1277 return true;
1278}
1279
1280//______________________________________________
1281template <typename value_T>
1282GPUd() bool TrackParametrizationWithError<value_T>::update(const TrackParametrizationWithError<value_T>& rhs)
1283{
1284 // update track with other track
1285 MatrixDSym5 covI; // perform matrix operations in double!
1286 buildCombinedCovMatrix(rhs, covI);
1287 if (!covI.Invert()) {
1288 LOG(warning) << "Cov.matrix inversion failed: " << covI;
1289 return false;
1290 }
1291 return update(rhs, covI);
1292}
1293
1294//______________________________________________
1295template <typename value_T>
1296GPUd() bool TrackParametrizationWithError<value_T>::update(const value_t* p, const value_t* cov)
1297{
1298 // Update the track parameters with the space point "p" having
1299 // the covariance matrix "cov"
1300
1301 value_t &cm00 = mC[kSigY2], &cm10 = mC[kSigZY], &cm11 = mC[kSigZ2], &cm20 = mC[kSigSnpY], &cm21 = mC[kSigSnpZ],
1302 &cm22 = mC[kSigSnp2], &cm30 = mC[kSigTglY], &cm31 = mC[kSigTglZ], &cm32 = mC[kSigTglSnp], &cm33 = mC[kSigTgl2],
1303 &cm40 = mC[kSigQ2PtY], &cm41 = mC[kSigQ2PtZ], &cm42 = mC[kSigQ2PtSnp], &cm43 = mC[kSigQ2PtTgl],
1304 &cm44 = mC[kSigQ2Pt2];
1305
1306 // use double precision?
1307 double r00 = static_cast<double>(cov[0]) + static_cast<double>(cm00);
1308 double r01 = static_cast<double>(cov[1]) + static_cast<double>(cm10);
1309 double r11 = static_cast<double>(cov[2]) + static_cast<double>(cm11);
1310 double det = r00 * r11 - r01 * r01;
1311
1312 if (gpu::CAMath::Abs(det) < constants::math::Almost0) {
1313 return false;
1314 }
1315 double detI = 1. / det;
1316 double tmp = r00;
1317 r00 = r11 * detI;
1318 r11 = tmp * detI;
1319 r01 = -r01 * detI;
1320
1321 double k00 = cm00 * r00 + cm10 * r01, k01 = cm00 * r01 + cm10 * r11;
1322 double k10 = cm10 * r00 + cm11 * r01, k11 = cm10 * r01 + cm11 * r11;
1323 double k20 = cm20 * r00 + cm21 * r01, k21 = cm20 * r01 + cm21 * r11;
1324 double k30 = cm30 * r00 + cm31 * r01, k31 = cm30 * r01 + cm31 * r11;
1325 double k40 = cm40 * r00 + cm41 * r01, k41 = cm40 * r01 + cm41 * r11;
1326
1327 value_t dy = p[kY] - this->getY(), dz = p[kZ] - this->getZ();
1328 value_t dsnp = k20 * dy + k21 * dz;
1329 if (gpu::CAMath::Abs(this->getSnp() + dsnp) > constants::math::Almost1) {
1330 return false;
1331 }
1332
1333 const params_t dP{value_t(k00 * dy + k01 * dz), value_t(k10 * dy + k11 * dz), dsnp, value_t(k30 * dy + k31 * dz),
1334 value_t(k40 * dy + k41 * dz)};
1335 this->updateParams(dP);
1336
1337 double c01 = cm10, c02 = cm20, c03 = cm30, c04 = cm40;
1338 double c12 = cm21, c13 = cm31, c14 = cm41;
1339
1340 cm00 -= k00 * cm00 + k01 * cm10;
1341 cm10 -= k00 * c01 + k01 * cm11;
1342 cm20 -= k00 * c02 + k01 * c12;
1343 cm30 -= k00 * c03 + k01 * c13;
1344 cm40 -= k00 * c04 + k01 * c14;
1345
1346 cm11 -= k10 * c01 + k11 * cm11;
1347 cm21 -= k10 * c02 + k11 * c12;
1348 cm31 -= k10 * c03 + k11 * c13;
1349 cm41 -= k10 * c04 + k11 * c14;
1350
1351 cm22 -= k20 * c02 + k21 * c12;
1352 cm32 -= k20 * c03 + k21 * c13;
1353 cm42 -= k20 * c04 + k21 * c14;
1354
1355 cm33 -= k30 * c03 + k31 * c13;
1356 cm43 -= k30 * c04 + k31 * c14;
1357
1358 cm44 -= k40 * c04 + k41 * c14;
1359
1360 checkCovariance();
1361
1362 return true;
1363}
1364
1365//______________________________________________
1366template <typename value_T>
1367GPUd() value_T TrackParametrizationWithError<value_T>::update(const o2::dataformats::VertexBase& vtx, value_T maxChi2)
1368{
1369 // Update track with vertex if the track-vertex chi2 does not exceed maxChi2. Track must be already propagated to the DCA to vertex
1370 // return update chi2 or -chi2 if chi2 if chi2 exceeds maxChi2
1371 auto vtLoc = this->getVertexInTrackFrame(vtx);
1372 value_T chi2 = getPredictedChi2(vtLoc.yz, vtLoc.yzerr);
1373 return chi2 < maxChi2 && update(vtLoc.yz, vtLoc.yzerr) ? chi2 : -chi2;
1374}
1375
1376//______________________________________________
1377template <typename value_T>
1378GPUd() bool TrackParametrizationWithError<value_T>::correctForMaterial(value_t x2x0, value_t xrho, bool anglecorr)
1379{
1380 //------------------------------------------------------------------
1381 // This function corrects the track parameters for the crossed material.
1382 // "x2x0" - X/X0, the thickness in units of the radiation length.
1383 // "xrho" - is the product length*density (g/cm^2).
1384 // It should be passed as negative when propagating tracks
1385 // from the intreaction point to the outside of the central barrel.
1386 // "dedx" - mean enery loss (GeV/(g/cm^2), if <=kCalcdEdxAuto : calculate on the fly
1387 // "anglecorr" - switch for the angular correction
1388 //------------------------------------------------------------------
1389 constexpr value_t kMSConst2 = 0.0136f * 0.0136f;
1390 constexpr value_t kMinP = 0.01f; // kill below this momentum
1391
1392 value_t csp2 = (1.f - this->getSnp()) * (1.f + this->getSnp()); // cos(phi)^2
1393 value_t cst2I = (1.f + this->getTgl() * this->getTgl()); // 1/cos(lambda)^2
1394 if (anglecorr) { // Apply angle correction, if requested
1395 value_t angle = gpu::CAMath::Sqrt(cst2I / csp2);
1396 x2x0 *= angle;
1397 xrho *= angle;
1398 }
1399 auto m = this->getPID().getMass();
1400 int charge2 = this->getAbsCharge() * this->getAbsCharge();
1401 value_t p = this->getP(), p0 = p, p02 = p * p, e2 = p02 + this->getPID().getMass2(), massInv = 1. / m, bg = p * massInv, dETot = 0.;
1402 value_t e = gpu::CAMath::Sqrt(e2), e0 = e;
1403 if (m > 0 && xrho != 0.f) {
1404 value_t ekin = e - m, dedx = this->getdEdxBBOpt(bg);
1405#ifdef _BB_NONCONST_CORR_
1406 value_t dedxDer = 0., dedx1 = dedx;
1407#endif
1408 if (charge2 != 1) {
1409 dedx *= charge2;
1410 }
1411 value_t dE = dedx * xrho;
1412 int na = 1 + int(gpu::CAMath::Abs(dE) / ekin * ELoss2EKinThreshInv);
1413 if (na > MaxELossIter) {
1414 na = MaxELossIter;
1415 }
1416 if (na > 1) {
1417 dE /= na;
1418 xrho /= na;
1419#ifdef _BB_NONCONST_CORR_
1420 dedxDer = this->getBetheBlochSolidDerivativeApprox(dedx1, bg); // require correction for non-constantness of dedx vs betagamma
1421 if (charge2 != 1) {
1422 dedxDer *= charge2;
1423 }
1424#endif
1425 }
1426 while (na--) {
1427#ifdef _BB_NONCONST_CORR_
1428 if (dedxDer != 0.) { // correction for non-constantness of dedx vs beta*gamma (in linear approximation): for a single step dE -> dE * [(exp(dedxDer) - 1)/dedxDer]
1429 if (xrho < 0) {
1430 dedxDer = -dedxDer; // E.loss ( -> positive derivative)
1431 }
1432 auto corrC = (gpu::CAMath::Exp(dedxDer) - 1.) / dedxDer;
1433 dE *= corrC;
1434 }
1435#endif
1436 e += dE;
1437 if (e > m) { // stopped
1438 p = gpu::CAMath::Sqrt(e * e - this->getPID().getMass2());
1439 } else {
1440 return false;
1441 }
1442 if (na) {
1443 bg = p * massInv;
1444 dedx = this->getdEdxBBOpt(bg);
1445#ifdef _BB_NONCONST_CORR_
1446 dedxDer = this->getBetheBlochSolidDerivativeApprox(dedx, bg);
1447#endif
1448 if (charge2 != 1) {
1449 dedx *= charge2;
1450#ifdef _BB_NONCONST_CORR_
1451 dedxDer *= charge2;
1452#endif
1453 }
1454 dE = dedx * xrho;
1455 }
1456 }
1457
1458 if (p < kMinP) {
1459 return false;
1460 }
1461 dETot = e - e0;
1462 } // end of e.loss correction
1463
1464 // Calculating the multiple scattering corrections******************
1465 value_t& fC22 = mC[kSigSnp2];
1466 value_t& fC33 = mC[kSigTgl2];
1467 value_t& fC43 = mC[kSigQ2PtTgl];
1468 value_t& fC44 = mC[kSigQ2Pt2];
1469 //
1470 value_t cC22(0.f), cC33(0.f), cC43(0.f), cC44(0.f);
1471 if (x2x0 != 0.f) {
1472 value_t beta2 = p02 / e2, theta2 = kMSConst2 / (beta2 * p02) * gpu::CAMath::Abs(x2x0);
1473 value_t fp34 = this->getTgl();
1474 if (charge2 != 1) {
1475 theta2 *= charge2;
1476 fp34 *= this->getCharge2Pt();
1477 }
1478 if (theta2 > constants::math::PI * constants::math::PI) {
1479 return false;
1480 }
1481 value_t t2c2I = theta2 * cst2I;
1482 cC22 = t2c2I * csp2;
1483 cC33 = t2c2I * cst2I;
1484 cC43 = t2c2I * fp34;
1485 cC44 = theta2 * fp34 * fp34;
1486 // optimize this
1487 // cC22 = theta2*((1.-getSnp())*(1.+getSnp()))*(1. + this->getTgl()*getTgl());
1488 // cC33 = theta2*(1. + this->getTgl()*getTgl())*(1. + this->getTgl()*getTgl());
1489 // cC43 = theta2*getTgl()*this->getQ2Pt()*(1. + this->getTgl()*getTgl());
1490 // cC44 = theta2*getTgl()*this->getQ2Pt()*getTgl()*this->getQ2Pt();
1491 }
1492
1493 // the energy loss correction contribution to cov.matrix: approximate energy loss fluctuation (M.Ivanov)
1494 constexpr value_t knst = 0.0007f; // To be tuned.
1495 value_t sigmadE = knst * gpu::CAMath::Sqrt(gpu::CAMath::Abs(dETot)) * e0 / p02 * this->getCharge2Pt();
1496 cC44 += sigmadE * sigmadE;
1497
1498 // Applying the corrections*****************************
1499 fC22 += cC22;
1500 fC33 += cC33;
1501 fC43 += cC43;
1502 fC44 += cC44;
1503 this->setQ2Pt(this->getQ2Pt() * p0 / p);
1504
1505 checkCovariance();
1506
1507 return true;
1508}
1509
1510//______________________________________________
1511template <typename value_T>
1512GPUd() bool TrackParametrizationWithError<value_T>::correctForMaterial(TrackParametrization<value_T>& linRef, value_t x2x0, value_t xrho, bool anglecorr)
1513{
1514 //------------------------------------------------------------------
1515 // This function corrects the reference and current track parameters for the crossed material
1516 // "x2x0" - X/X0, the thickness in units of the radiation length.
1517 // "xrho" - is the product length*density (g/cm^2).
1518 // It should be passed as negative when propagating tracks
1519 // from the intreaction point to the outside of the central barrel.
1520 // "dedx" - mean enery loss (GeV/(g/cm^2), if <=kCalcdEdxAuto : calculate on the fly
1521 // "anglecorr" - switch for the angular correction
1522 //------------------------------------------------------------------
1523 constexpr value_t kMSConst2 = 0.0136f * 0.0136f;
1524 constexpr value_t kMinP = 0.01f; // kill below this momentum
1525
1526 value_t csp2 = (1.f - linRef.getSnp()) * (1.f + linRef.getSnp()); // cos(phi)^2
1527 value_t cst2I = (1.f + linRef.getTgl() * linRef.getTgl()); // 1/cos(lambda)^2
1528 if (anglecorr) { // Apply angle correction, if requested
1529 value_t angle = gpu::CAMath::Sqrt(cst2I / csp2);
1530 x2x0 *= angle;
1531 xrho *= angle;
1532 }
1533 auto pid = linRef.getPID();
1534 auto m = pid.getMass();
1535 int charge2 = linRef.getAbsCharge() * linRef.getAbsCharge();
1536 value_t p = linRef.getP(), p0 = p, p02 = p * p, e2 = p02 + pid.getMass2(), massInv = 1. / m, bg = p * massInv, dETot = 0.;
1537 value_t e = gpu::CAMath::Sqrt(e2), e0 = e;
1538 if (m > 0 && xrho != 0.f) {
1539 value_t ekin = e - m, dedx = this->getdEdxBBOpt(bg);
1540#ifdef _BB_NONCONST_CORR_
1541 value_t dedxDer = 0., dedx1 = dedx;
1542#endif
1543 if (charge2 != 1) {
1544 dedx *= charge2;
1545 }
1546 value_t dE = dedx * xrho;
1547 int na = 1 + int(gpu::CAMath::Abs(dE) / ekin * ELoss2EKinThreshInv);
1548 if (na > MaxELossIter) {
1549 na = MaxELossIter;
1550 }
1551 if (na > 1) {
1552 dE /= na;
1553 xrho /= na;
1554#ifdef _BB_NONCONST_CORR_
1555 dedxDer = this->getBetheBlochSolidDerivativeApprox(dedx1, bg); // require correction for non-constantness of dedx vs betagamma
1556 if (charge2 != 1) {
1557 dedxDer *= charge2;
1558 }
1559#endif
1560 }
1561 while (na--) {
1562#ifdef _BB_NONCONST_CORR_
1563 if (dedxDer != 0.) { // correction for non-constantness of dedx vs beta*gamma (in linear approximation): for a single step dE -> dE * [(exp(dedxDer) - 1)/dedxDer]
1564 if (xrho < 0) {
1565 dedxDer = -dedxDer; // E.loss ( -> positive derivative)
1566 }
1567 auto corrC = (gpu::CAMath::Exp(dedxDer) - 1.) / dedxDer;
1568 dE *= corrC;
1569 }
1570#endif
1571 e += dE;
1572 if (e > m) { // stopped
1573 p = gpu::CAMath::Sqrt(e * e - pid.getMass2());
1574 } else {
1575 return false;
1576 }
1577 if (na) {
1578 bg = p * massInv;
1579 dedx = this->getdEdxBBOpt(bg);
1580#ifdef _BB_NONCONST_CORR_
1581 dedxDer = this->getBetheBlochSolidDerivativeApprox(dedx, bg);
1582#endif
1583 if (charge2 != 1) {
1584 dedx *= charge2;
1585#ifdef _BB_NONCONST_CORR_
1586 dedxDer *= charge2;
1587#endif
1588 }
1589 dE = dedx * xrho;
1590 }
1591 }
1592
1593 if (p < kMinP) {
1594 return false;
1595 }
1596 dETot = e - e0;
1597 } // end of e.loss correction
1598
1599 // Calculating the multiple scattering corrections******************
1600 value_t& fC22 = mC[kSigSnp2];
1601 value_t& fC33 = mC[kSigTgl2];
1602 value_t& fC43 = mC[kSigQ2PtTgl];
1603 value_t& fC44 = mC[kSigQ2Pt2];
1604 //
1605 value_t cC22(0.f), cC33(0.f), cC43(0.f), cC44(0.f);
1606 if (x2x0 != 0.f) {
1607 value_t beta2 = p02 / e2, theta2 = kMSConst2 / (beta2 * p02) * gpu::CAMath::Abs(x2x0);
1608 value_t fp34 = linRef.getTgl();
1609 if (charge2 != 1) {
1610 theta2 *= charge2;
1611 fp34 *= linRef.getCharge2Pt();
1612 }
1613 if (theta2 > constants::math::PI * constants::math::PI) {
1614 return false;
1615 }
1616 value_t t2c2I = theta2 * cst2I;
1617 cC22 = t2c2I * csp2;
1618 cC33 = t2c2I * cst2I;
1619 cC43 = t2c2I * fp34;
1620 cC44 = theta2 * fp34 * fp34;
1621 // optimize this
1622 // cC22 = theta2*((1.-getSnp())*(1.+getSnp()))*(1. + this->getTgl()*getTgl());
1623 // cC33 = theta2*(1. + this->getTgl()*getTgl())*(1. + this->getTgl()*getTgl());
1624 // cC43 = theta2*getTgl()*this->getQ2Pt()*(1. + this->getTgl()*getTgl());
1625 // cC44 = theta2*getTgl()*this->getQ2Pt()*getTgl()*this->getQ2Pt();
1626 }
1627
1628 // the energy loss correction contribution to cov.matrix: approximate energy loss fluctuation (M.Ivanov)
1629 constexpr value_t knst = 0.0007f; // To be tuned.
1630 value_t sigmadE = knst * gpu::CAMath::Sqrt(gpu::CAMath::Abs(dETot)) * e0 / p02 * linRef.getCharge2Pt();
1631 cC44 += sigmadE * sigmadE;
1632
1633 // Applying the corrections*****************************
1634 fC22 += cC22;
1635 fC33 += cC33;
1636 fC43 += cC43;
1637 fC44 += cC44;
1638 auto pscale = p0 / p;
1639 linRef.setQ2Pt(linRef.getQ2Pt() * pscale);
1640 this->setQ2Pt(this->getQ2Pt() * pscale);
1641
1642 checkCovariance();
1643
1644 return true;
1645}
1646
1647//______________________________________________________________
1648template <typename value_T>
1649GPUd() bool TrackParametrizationWithError<value_T>::getCovXYZPxPyPzGlo(std::array<value_t, kLabCovMatSize>& cv) const
1650{
1651 //---------------------------------------------------------------------
1652 // This function returns the global covariance matrix of the track params
1653 //
1654 // Cov(x,x) ... : cv[0]
1655 // Cov(y,x) ... : cv[1] cv[2]
1656 // Cov(z,x) ... : cv[3] cv[4] cv[5]
1657 // Cov(px,x)... : cv[6] cv[7] cv[8] cv[9]
1658 // Cov(py,x)... : cv[10] cv[11] cv[12] cv[13] cv[14]
1659 // Cov(pz,x)... : cv[15] cv[16] cv[17] cv[18] cv[19] cv[20]
1660 //
1661 // Results for (nearly) straight tracks are meaningless !
1662 //---------------------------------------------------------------------
1663 if (gpu::CAMath::Abs(this->getQ2Pt()) <= constants::math::Almost0 || gpu::CAMath::Abs(this->getSnp()) > constants::math::Almost1) {
1664 for (int i = 0; i < 21; i++) {
1665 cv[i] = 0.;
1666 }
1667 return false;
1668 }
1669
1670 auto pt = this->getPt();
1671 value_t sn, cs;
1672 o2::math_utils::detail::sincos(this->getAlpha(), sn, cs);
1673 auto r = gpu::CAMath::Sqrt((1. - this->getSnp()) * (1. + this->getSnp()));
1674 auto m00 = -sn, m10 = cs;
1675 auto m23 = -pt * (sn + this->getSnp() * cs / r), m43 = -pt * pt * (r * cs - this->getSnp() * sn);
1676 auto m24 = pt * (cs - this->getSnp() * sn / r), m44 = -pt * pt * (r * sn + this->getSnp() * cs);
1677 auto m35 = pt, m45 = -pt * pt * this->getTgl();
1678
1679 if (this->getSign() < 0) {
1680 m43 = -m43;
1681 m44 = -m44;
1682 m45 = -m45;
1683 }
1684
1685 cv[0] = mC[0] * m00 * m00;
1686 cv[1] = mC[0] * m00 * m10;
1687 cv[2] = mC[0] * m10 * m10;
1688 cv[3] = mC[1] * m00;
1689 cv[4] = mC[1] * m10;
1690 cv[5] = mC[2];
1691 cv[6] = m00 * (mC[3] * m23 + mC[10] * m43);
1692 cv[7] = m10 * (mC[3] * m23 + mC[10] * m43);
1693 cv[8] = mC[4] * m23 + mC[11] * m43;
1694 cv[9] = m23 * (mC[5] * m23 + mC[12] * m43) + m43 * (mC[12] * m23 + mC[14] * m43);
1695 cv[10] = m00 * (mC[3] * m24 + mC[10] * m44);
1696 cv[11] = m10 * (mC[3] * m24 + mC[10] * m44);
1697 cv[12] = mC[4] * m24 + mC[11] * m44;
1698 cv[13] = m23 * (mC[5] * m24 + mC[12] * m44) + m43 * (mC[12] * m24 + mC[14] * m44);
1699 cv[14] = m24 * (mC[5] * m24 + mC[12] * m44) + m44 * (mC[12] * m24 + mC[14] * m44);
1700 cv[15] = m00 * (mC[6] * m35 + mC[10] * m45);
1701 cv[16] = m10 * (mC[6] * m35 + mC[10] * m45);
1702 cv[17] = mC[7] * m35 + mC[11] * m45;
1703 cv[18] = m23 * (mC[8] * m35 + mC[12] * m45) + m43 * (mC[13] * m35 + mC[14] * m45);
1704 cv[19] = m24 * (mC[8] * m35 + mC[12] * m45) + m44 * (mC[13] * m35 + mC[14] * m45);
1705 cv[20] = m35 * (mC[9] * m35 + mC[13] * m45) + m45 * (mC[13] * m35 + mC[14] * m45);
1706
1707 return true;
1708}
1709
1710#ifndef GPUCA_ALIGPUCODE
1711//______________________________________________________________
1712template <typename value_T>
1713std::string TrackParametrizationWithError<value_T>::asString() const
1714{
1715 return TrackParametrization<value_t>::asString() +
1716 fmt::format(" Cov: [{:+.3e}] [{:+.3e} {:+.3e}] [{:+.3e} {:+.3e} {:+.3e}] [{:+.3e} {:+.3e} {:+.3e} {:+.3e}] [{:+.3e} {:+.3e} {:+.3e} {:+.3e} {:+.3e}]",
1717 mC[kSigY2], mC[kSigZY], mC[kSigZ2], mC[kSigSnpY], mC[kSigSnpZ], mC[kSigSnp2], mC[kSigTglY],
1718 mC[kSigTglZ], mC[kSigTglSnp], mC[kSigTgl2], mC[kSigQ2PtY], mC[kSigQ2PtZ], mC[kSigQ2PtSnp], mC[kSigQ2PtTgl],
1719 mC[kSigQ2Pt2]);
1720}
1721
1722template <typename value_T>
1723std::string TrackParametrizationWithError<value_T>::asStringHexadecimal()
1724{
1725 return TrackParametrization<value_t>::asStringHexadecimal() +
1726 fmt::format(" Cov: [{:x}] [{:x} {:x}] [{:x} {:x} {:x}] [{:x} {:x} {:x} {:x}] [{:x} {:x} {:x} {:x} {:x}]",
1727 reinterpret_cast<const unsigned int&>(mC[kSigY2]), reinterpret_cast<const unsigned int&>(mC[kSigZY]), reinterpret_cast<const unsigned int&>(mC[kSigZ2]),
1728 reinterpret_cast<const unsigned int&>(mC[kSigSnpY]), reinterpret_cast<const unsigned int&>(mC[kSigSnpZ]), reinterpret_cast<const unsigned int&>(mC[kSigSnp2]),
1729 reinterpret_cast<const unsigned int&>(mC[kSigTglY]), reinterpret_cast<const unsigned int&>(mC[kSigTglZ]), reinterpret_cast<const unsigned int&>(mC[kSigTglSnp]),
1730 reinterpret_cast<const unsigned int&>(mC[kSigTgl2]), reinterpret_cast<const unsigned int&>(mC[kSigQ2PtY]), reinterpret_cast<const unsigned int&>(mC[kSigQ2PtZ]),
1731 reinterpret_cast<const unsigned int&>(mC[kSigQ2PtSnp]), reinterpret_cast<const unsigned int&>(mC[kSigQ2PtTgl]), reinterpret_cast<const unsigned int&>(mC[kSigQ2Pt2]));
1732}
1733#endif
1734
1735//______________________________________________________________
1736template <typename value_T>
1737GPUd() void TrackParametrizationWithError<value_T>::print() const
1738{
1739 // print parameters
1740#ifndef GPUCA_ALIGPUCODE
1741 printf("%s\n", asString().c_str());
1742#elif !defined(GPUCA_GPUCODE_DEVICE) || (!defined(__OPENCL__) && defined(GPUCA_GPU_DEBUG_PRINT))
1743 TrackParametrization<value_T>::printParam();
1744 printf(
1745 " Cov: [%+.3e] [%+.3e %+.3e] [%+.3e %+.3e %+.3e] [%+.3e %+.3e %+.3e %+.3e] [%+.3e %+.3e %+.3e %+.3e %+.3e]\n",
1746 mC[kSigY2], mC[kSigZY], mC[kSigZ2], mC[kSigSnpY], mC[kSigSnpZ], mC[kSigSnp2], mC[kSigTglY],
1747 mC[kSigTglZ], mC[kSigTglSnp], mC[kSigTgl2], mC[kSigQ2PtY], mC[kSigQ2PtZ], mC[kSigQ2PtSnp], mC[kSigQ2PtTgl],
1748 mC[kSigQ2Pt2]);
1749#endif
1750}
1751
1752//______________________________________________________________
1753template <typename value_T>
1754GPUd() void TrackParametrizationWithError<value_T>::printHexadecimal()
1755{
1756 // print parameters
1757#ifndef GPUCA_ALIGPUCODE
1758 printf("%s\n", asStringHexadecimal().c_str());
1759#elif !defined(GPUCA_GPUCODE_DEVICE) || (!defined(__OPENCL__) && defined(GPUCA_GPU_DEBUG_PRINT))
1760 TrackParametrization<value_T>::printParamHexadecimal();
1761 printf(
1762 " Cov: [%x] [%x %x] [%x %x %x] [%x %x %x %x] [%x %x %x %x %x]\n",
1763 gpu::CAMath::Float2UIntReint(mC[kSigY2]),
1764 gpu::CAMath::Float2UIntReint(mC[kSigZY]), gpu::CAMath::Float2UIntReint(mC[kSigZ2]),
1765 gpu::CAMath::Float2UIntReint(mC[kSigSnpY]), gpu::CAMath::Float2UIntReint(mC[kSigSnpZ]), gpu::CAMath::Float2UIntReint(mC[kSigSnp2]),
1766 gpu::CAMath::Float2UIntReint(mC[kSigTglY]), gpu::CAMath::Float2UIntReint(mC[kSigTglZ]), gpu::CAMath::Float2UIntReint(mC[kSigTglSnp]), gpu::CAMath::Float2UIntReint(mC[kSigTgl2]),
1767 gpu::CAMath::Float2UIntReint(mC[kSigQ2PtY]), gpu::CAMath::Float2UIntReint(mC[kSigQ2PtZ]), gpu::CAMath::Float2UIntReint(mC[kSigQ2PtSnp]), gpu::CAMath::Float2UIntReint(mC[kSigQ2PtTgl]), gpu::CAMath::Float2UIntReint(mC[kSigQ2Pt2]));
1768#endif
1769}
1770
1771namespace o2::track
1772{
1773#if !defined(GPUCA_GPUCODE) || defined(GPUCA_GPUCODE_DEVICE) // FIXME: DR: WORKAROUND to avoid CUDA bug creating host symbols for device code.
1774template class TrackParametrizationWithError<float>;
1775#endif
1776#ifndef GPUCA_GPUCODE
1777template class TrackParametrizationWithError<double>;
1778#endif
1779} // namespace o2::track
asString
std::string asString(TDataMember const &dm, char *pointer)
Definition ConfigurableParamHelper.cxx:151
debug
std::ostringstream debug
Definition ExpressionJSONHelpers.cxx:83
charge
int16_t charge
Definition RawEventData.h:5
print
void print() const
i
int32_t i
Definition GPUCommonAlgorithm.h:436
bb
const int16_t bb
Definition GPUTPCGMMerger.cxx:1480
MathConstants.h
useful math constants
j
uint32_t j
Definition RawData.h:0
pid
uint16_t pid
Definition RawData.h:2
o2::track::TrackParametrization::asString
std::string asString() const
alpha
GLfloat GLfloat GLfloat alpha
Definition glcorearb.h:279
x
GLint GLenum GLint x
Definition glcorearb.h:403
m
const GLfloat * m
Definition glcorearb.h:4066
array
GLenum array
Definition glcorearb.h:4274
b
GLboolean GLboolean GLboolean b
Definition glcorearb.h:1233
void
typedef void(APIENTRYP PFNGLCULLFACEPROC)(GLenum mode)
angle
GLfloat angle
Definition glcorearb.h:4071
r
GLboolean r
Definition glcorearb.h:1233
z
GLdouble GLdouble GLdouble z
Definition glcorearb.h:843
invert
GLboolean invert
Definition glcorearb.h:543
o2::algorithm::pageparser::set
void set(T *h, size_t v)
Definition PageParser.h:107
o2::align::params_t
typename trackParam_t::params_t params_t
Definition utils.h:33
o2::align::dim3_t
typename trackParam_t::dim3_t dim3_t
Definition utils.h:32
o2::align::value_t
typename trackParam_t::value_t value_t
Definition utils.h:30
o2::constants::math::Epsilon
constexpr float Epsilon
Definition MathConstants.h:26
o2::constants::math::Almost1
constexpr float Almost1
Definition MathConstants.h:27
o2::its::math_utils::p
float p
Definition MathUtils.h:99
o2::math_utils::detail::rhs
D const SVectorGPU< T, D > & rhs
Definition SMatrixGPU.h:193
o2::mch::raw::p0
std::array< int, 24 > p0
Definition SampaHeader.cxx:491
o2::mch::getX
double * getX(double *xyDxy, int N)
Definition mathieson.cxx:1006
o2::mch::getY
double * getY(double *xyDxy, int N)
Definition mathieson.cxx:1007
o2::rans::this
return * this
Definition DenseHistogram.h:573
o2::track::bg
value_T bg
Definition TrackUtils.h:194
o2::track::kCTgl2max
constexpr float kCTgl2max
Definition TrackParametrization.h:102
o2::track::HugeF
constexpr float HugeF
Definition TrackParametrization.h:117
o2::track::kCovMatSize
constexpr int kCovMatSize
Definition TrackParametrization.h:97
o2::track::kCSnp2max
constexpr float kCSnp2max
Definition TrackParametrization.h:101
o2::track::kNParams
constexpr int kNParams
Definition TrackParametrization.h:97
o2::track::kLabCovMatSize
constexpr int kLabCovMatSize
Definition TrackParametrization.h:97
o2::track::step
value_T step
Definition TrackUtils.h:42
o2::track::f1
value_T f1
Definition TrackUtils.h:91
o2::track::x
const value_T x
Definition TrackUtils.h:136
o2::track::kCY2max
constexpr float kCY2max
Definition TrackParametrization.h:99
o2::track::DefaultDCACov
constexpr float DefaultDCACov
Definition TrackParametrization.h:123
o2::track::d2
value_T d2
Definition TrackUtils.h:135
o2::track::vect
value_T std::array< value_T, 7 > & vect
Definition TrackUtils.h:42
o2::track::beta2
kp1 *kp2 *value_T beta2
Definition TrackUtils.h:131
o2::track::GPUd
GPUd() value_T BetheBlochSolid(value_T bg
Definition TrackUtils.h:150
o2::track::kC1Pt2max
constexpr float kC1Pt2max
Definition TrackParametrization.h:103
o2::track::f2
value_T f2
Definition TrackUtils.h:92
o2::track::DefaultDCA
constexpr float DefaultDCA
Definition TrackParametrization.h:122
o2::track::kCZ2max
constexpr float kCZ2max
Definition TrackParametrization.h:100
o2::track::MaxELossIter
constexpr int MaxELossIter
Definition TrackParametrization.h:121
o2::track::ELoss2EKinThreshInv
constexpr float ELoss2EKinThreshInv
Definition TrackParametrization.h:120
o2
a couple of static helper functions to create timestamp values for CCDB queries or override obsolete ...
Definition BitstreamReader.h:24
std
Defining DataPointCompositeObject explicitly as copiable.
Definition CcdbObjectInfo.h:121
cc
std::vector< o2::mch::ChannelCode > cc
Definition testStatusMap.cxx:37
LOG
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"