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