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