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 const value_t pt2 = mom[0] * mom[0] + mom[1] * mom[1];
531 const value_t pt = gpu::CAMath::Sqrt(pt2);
532 const value_t ptI = 1.f / pt;
533 this->setX(ver[0]);
534 this->setAlpha(alp);
535 this->setY(ver[1]);
536 this->setZ(ver[2]);
537 this->setSnp(mom[1] * ptI); // cos(phi)
538 this->setTgl(mom[2] * ptI); // tg(lambda)
539 this->setAbsCharge(gpu::CAMath::Abs(charge));
540 this->setQ2Pt(charge ? ptI * charge : ptI);
541 this->setPID(pid);
542 //
543 if (gpu::CAMath::Abs(1.f - this->getSnp()) < kSafe) {
544 this->setSnp(1.f - kSafe); // Protection
545 } else if (gpu::CAMath::Abs(-1.f - this->getSnp()) < kSafe) {
546 this->setSnp(-1.f + kSafe); // Protection
547 }
548 //
549 // Covariance matrix from the fixed-alpha Jacobian
550 // d(Y,Z,snp,tgl,q/pt) / d(X,Y,Z,Px,Py,Pz).
551 const value_t pt3I = ptI / pt2;
552 const value_t qeff = charge ? static_cast<value_t>(charge) : 1.f;
553
554 value_t cLab[6][6] = {};
555 int idx = 0;
556 for (int i = 0; i < 6; ++i) {
557 for (int j = 0; j <= i; ++j) {
558 cLab[i][j] = cLab[j][i] = cv[idx++];
559 }
560 }
561
562 value_t jac[5][6] = {};
563 jac[kY][0] = -sn;
564 jac[kY][1] = cs;
565 jac[kZ][2] = 1.;
566
567 const value_t u = mom[0];
568 const value_t v = mom[1];
569 const value_t w = mom[2];
570 const value_t dSnpDu = -u * v * pt3I;
571 const value_t dSnpDv = u * u * pt3I;
572 const value_t dTglDu = -w * u * pt3I;
573 const value_t dTglDv = -w * v * pt3I;
574 const value_t dTglDw = ptI;
575 const value_t dQ2PtDu = -qeff * u * pt3I;
576 const value_t dQ2PtDv = -qeff * v * pt3I;
577
578 jac[kSnp][3] = dSnpDu * cs - dSnpDv * sn;
579 jac[kSnp][4] = dSnpDu * sn + dSnpDv * cs;
580 jac[kTgl][3] = dTglDu * cs - dTglDv * sn;
581 jac[kTgl][4] = dTglDu * sn + dTglDv * cs;
582 jac[kTgl][5] = dTglDw;
583 jac[kQ2Pt][3] = dQ2PtDu * cs - dQ2PtDv * sn;
584 jac[kQ2Pt][4] = dQ2PtDu * sn + dQ2PtDv * cs;
585
586 for (int i = 0; i < kNParams; ++i) {
587 for (int j = 0; j <= i; ++j) {
588 value_t cij = 0.;
589 for (int k = 0; k < 6; ++k) {
590 for (int l = 0; l < 6; ++l) {
591 cij += jac[i][k] * cLab[k][l] * jac[j][l];
592 }
593 }
594 mC[CovarMap[i][j]] = cij;
595 }
596 }
597 checkCovariance();
598}
599
600//____________________________________________________________
601template <typename value_T>
602GPUd() bool TrackParametrizationWithError<value_T>::propagateTo(value_t xk, const dim3_t& b)
603{
604 //----------------------------------------------------------------
605 // Extrapolate this track to the plane X=xk in the field b[].
606 //
607 // X [cm] is in the "tracking coordinate system" of this track.
608 // b[]={Bx,By,Bz} [kG] is in the Global coordidate system.
609 //----------------------------------------------------------------
610
611 value_t dx = xk - this->getX();
612 if (gpu::CAMath::Abs(dx) < constants::math::Almost0) {
613 return true;
614 }
615 // Do not propagate tracks outside the ALICE detector
616 if (gpu::CAMath::Abs(dx) > 1e5 || gpu::CAMath::Abs(this->getY()) > 1e5 || gpu::CAMath::Abs(this->getZ()) > 1e5) {
617 LOG(warning) << "Anomalous track, target X:" << xk;
618 // print();
619 return false;
620 }
621 value_t crv = (gpu::CAMath::Abs(b[2]) < constants::math::Almost0) ? 0.f : this->getCurvature(b[2]);
622 if (gpu::CAMath::Abs(crv) < constants::math::Almost0) {
623 return propagateTo(xk, 0.);
624 }
625 value_t x2r = crv * dx;
626 value_t f1 = this->getSnp(), f2 = f1 + x2r;
627 if ((gpu::CAMath::Abs(f1) > constants::math::Almost1) || (gpu::CAMath::Abs(f2) > constants::math::Almost1)) {
628 return false;
629 }
630 value_t r1 = gpu::CAMath::Sqrt((1.f - f1) * (1.f + f1));
631 if (gpu::CAMath::Abs(r1) < constants::math::Almost0) {
632 return false;
633 }
634 value_t r2 = gpu::CAMath::Sqrt((1.f - f2) * (1.f + f2));
635 if (gpu::CAMath::Abs(r2) < constants::math::Almost0) {
636 return false;
637 }
638 double r1pr2Inv = 1. / (r1 + r2), r2inv = 1. / r2, r1inv = 1. / r1;
639 double dy2dx = (f1 + f2) * r1pr2Inv, dx2r1pr2 = dx * r1pr2Inv;
640 value_t step = (gpu::CAMath::Abs(x2r) < 0.05f) ? dx * gpu::CAMath::Abs(r2 + f2 * dy2dx) // chord
641 : 2.f * gpu::CAMath::ASin(0.5f * dx * gpu::CAMath::Sqrt(1.f + dy2dx * dy2dx) * crv) / crv; // arc
642 step *= gpu::CAMath::Sqrt(1.f + this->getTgl() * this->getTgl());
643 //
644 // get the track x,y,z,px/p,py/p,pz/p,p,sinAlpha,cosAlpha in the Global System
645 std::array<value_t, 9> vecLab{0.f};
646 if (!this->getPosDirGlo(vecLab)) {
647 return false;
648 }
649 //
650 // matrix transformed with Bz component only
651 value_t &c00 = mC[kSigY2], &c10 = mC[kSigZY], &c11 = mC[kSigZ2], &c20 = mC[kSigSnpY], &c21 = mC[kSigSnpZ],
652 &c22 = mC[kSigSnp2], &c30 = mC[kSigTglY], &c31 = mC[kSigTglZ], &c32 = mC[kSigTglSnp], &c33 = mC[kSigTgl2],
653 &c40 = mC[kSigQ2PtY], &c41 = mC[kSigQ2PtZ], &c42 = mC[kSigQ2PtSnp], &c43 = mC[kSigQ2PtTgl],
654 &c44 = mC[kSigQ2Pt2];
655
656 // evaluate matrix in double prec.
657 value_t kb = b[2] * constants::math::B2C;
658 double hh = dx2r1pr2 * r2inv * (1. + r1 * r2 + f1 * f2), jj = dx * (dy2dx - f2 * r2inv);
659 double f02 = hh * r1inv;
660 double f04 = hh * dx2r1pr2 * kb;
661 double f24 = dx * kb; // x2r/mP[kQ2Pt];
662 double f12 = this->getTgl() * (f02 * f2 + jj);
663 double f13 = dx * (r2 + f2 * dy2dx);
664 double f14 = this->getTgl() * (f04 * f2 + jj * f24);
665
666 // b = C*ft
667 double b00 = f02 * c20 + f04 * c40, b01 = f12 * c20 + f14 * c40 + f13 * c30;
668 double b02 = f24 * c40;
669 double b10 = f02 * c21 + f04 * c41, b11 = f12 * c21 + f14 * c41 + f13 * c31;
670 double b12 = f24 * c41;
671 double b20 = f02 * c22 + f04 * c42, b21 = f12 * c22 + f14 * c42 + f13 * c32;
672 double b22 = f24 * c42;
673 double b40 = f02 * c42 + f04 * c44, b41 = f12 * c42 + f14 * c44 + f13 * c43;
674 double b42 = f24 * c44;
675 double b30 = f02 * c32 + f04 * c43, b31 = f12 * c32 + f14 * c43 + f13 * c33;
676 double b32 = f24 * c43;
677
678 // a = f*b = f*C*ft
679 double a00 = f02 * b20 + f04 * b40, a01 = f02 * b21 + f04 * b41, a02 = f02 * b22 + f04 * b42;
680 double a11 = f12 * b21 + f14 * b41 + f13 * b31, a12 = f12 * b22 + f14 * b42 + f13 * b32;
681 double a22 = f24 * b42;
682
683 // F*C*Ft = C + (b + bt + a)
684 c00 += b00 + b00 + a00;
685 c10 += b10 + b01 + a01;
686 c20 += b20 + b02 + a02;
687 c30 += b30;
688 c40 += b40;
689 c11 += b11 + b11 + a11;
690 c21 += b21 + b12 + a12;
691 c31 += b31;
692 c41 += b41;
693 c22 += b22 + b22 + a22;
694 c32 += b32;
695 c42 += b42;
696
697 checkCovariance();
698
699 // Rotate to the system where Bx=By=0.
700 value_t bxy2 = b[0] * b[0] + b[1] * b[1];
701 value_t bt = gpu::CAMath::Sqrt(bxy2);
702 value_t cosphi = 1.f, sinphi = 0.f;
703 if (bt > constants::math::Almost0) {
704 cosphi = b[0] / bt;
705 sinphi = b[1] / bt;
706 }
707 value_t bb = gpu::CAMath::Sqrt(bxy2 + b[2] * b[2]);
708 value_t costet = 1., sintet = 0.;
709 if (bb > constants::math::Almost0) {
710 costet = b[2] / bb;
711 sintet = bt / bb;
712 }
713 std::array<value_t, 7> vect{costet * cosphi * vecLab[0] + costet * sinphi * vecLab[1] - sintet * vecLab[2],
714 -sinphi * vecLab[0] + cosphi * vecLab[1],
715 sintet * cosphi * vecLab[0] + sintet * sinphi * vecLab[1] + costet * vecLab[2],
716 costet * cosphi * vecLab[3] + costet * sinphi * vecLab[4] - sintet * vecLab[5],
717 -sinphi * vecLab[3] + cosphi * vecLab[4],
718 sintet * cosphi * vecLab[3] + sintet * sinphi * vecLab[4] + costet * vecLab[5],
719 vecLab[6]};
720
721 // Do the helix step
722 value_t q = this->getCharge();
723 g3helx3(q * bb, step, vect);
724
725 // Rotate back to the Global System
726 vecLab[0] = cosphi * costet * vect[0] - sinphi * vect[1] + cosphi * sintet * vect[2];
727 vecLab[1] = sinphi * costet * vect[0] + cosphi * vect[1] + sinphi * sintet * vect[2];
728 vecLab[2] = -sintet * vect[0] + costet * vect[2];
729
730 vecLab[3] = cosphi * costet * vect[3] - sinphi * vect[4] + cosphi * sintet * vect[5];
731 vecLab[4] = sinphi * costet * vect[3] + cosphi * vect[4] + sinphi * sintet * vect[5];
732 vecLab[5] = -sintet * vect[3] + costet * vect[5];
733
734 // Rotate back to the Tracking System
735 value_t sinalp = -vecLab[7], cosalp = vecLab[8];
736 value_t t = cosalp * vecLab[0] - sinalp * vecLab[1];
737 vecLab[1] = sinalp * vecLab[0] + cosalp * vecLab[1];
738 vecLab[0] = t;
739 t = cosalp * vecLab[3] - sinalp * vecLab[4];
740 vecLab[4] = sinalp * vecLab[3] + cosalp * vecLab[4];
741 vecLab[3] = t;
742
743 // Do the final correcting step to the target plane (linear approximation)
744 value_t x = vecLab[0], y = vecLab[1], z = vecLab[2];
745 if (gpu::CAMath::Abs(x - xk) > constants::math::Almost0) {
746 if (gpu::CAMath::Abs(vecLab[3]) < constants::math::Almost0) {
747 return false;
748 }
749 auto dxFin = xk - vecLab[0];
750 x += dxFin;
751 y += vecLab[4] / vecLab[3] * dxFin;
752 z += vecLab[5] / vecLab[3] * dxFin;
753 }
754
755 // Calculate the track parameters
756 t = 1.f / gpu::CAMath::Sqrt(vecLab[3] * vecLab[3] + vecLab[4] * vecLab[4]);
757 this->setX(xk);
758 this->setY(y);
759 this->setZ(z);
760 this->setSnp(vecLab[4] * t);
761 this->setTgl(vecLab[5] * t);
762 this->setQ2Pt(q * t / vecLab[6]);
763
764 return true;
765}
766
767//____________________________________________________________
768template <typename value_T>
769GPUd() bool TrackParametrizationWithError<value_T>::propagateTo(value_t xk, TrackParametrization<value_T>& linRef0, const dim3_t& b)
770{
771 //----------------------------------------------------------------
772 // Extrapolate this track to the plane X=xk in the field b[].
773 //
774 // X [cm] is in the "tracking coordinate system" of this track.
775 // b[]={Bx,By,Bz} [kG] is in the Global coordidate system.
776 //----------------------------------------------------------------
777
778 value_t dx = xk - this->getX();
779 if (gpu::CAMath::Abs(dx) < constants::math::Almost0) {
780 return true;
781 }
782 // Do not propagate tracks outside the ALICE detector
783 if (gpu::CAMath::Abs(dx) > 1e5 || gpu::CAMath::Abs(this->getY()) > 1e5 || gpu::CAMath::Abs(this->getZ()) > 1e5) {
784 LOG(warning) << "Anomalous track, target X:" << xk;
785 // print();
786 return false;
787 }
788 if (gpu::CAMath::Abs(dx) < constants::math::Almost0) {
789 this->setX(xk);
790 linRef0.setX(xk);
791 return true;
792 }
793 // preliminary calculations to find the step size
794 value_t crv = (gpu::CAMath::Abs(b[2]) < constants::math::Almost0) ? 0.f : linRef0.getCurvature(b[2]);
795 if (gpu::CAMath::Abs(crv) < constants::math::Almost0) {
796 return propagateTo(xk, linRef0, 0.);
797 }
798 value_t kb = b[2] * constants::math::B2C, x2r = crv * dx;
799 // evaluate in double prec.
800 value_t snpRef0 = linRef0.getSnp(), snpRef1 = snpRef0 + x2r;
801 if ((gpu::CAMath::Abs(snpRef0) > constants::math::Almost1) || (gpu::CAMath::Abs(snpRef1) > constants::math::Almost1)) {
802 return false;
803 }
804 value_t cspRef0 = gpu::CAMath::Sqrt((1 - snpRef0) * (1 + snpRef0)), cspRef1 = gpu::CAMath::Sqrt((1 - snpRef1) * (1 + snpRef1));
805 if (gpu::CAMath::Abs(cspRef0) < constants::math::Almost0 || gpu::CAMath::Abs(cspRef1) < constants::math::Almost0) {
806 return false;
807 }
808 value_t cspRef0Inv = value_t(1) / cspRef0, cspRef1Inv = value_t(1) / cspRef1, cc = cspRef0 + cspRef1, ccInv = value_t(1) / cc, dy2dx = (snpRef0 + snpRef1) * ccInv;
809 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
810 step *= gpu::CAMath::Sqrt(1.f + linRef0.getTgl() * linRef0.getTgl());
811
812 //
813 // get the track x,y,z,px/p,py/p,pz/p,p,sinAlpha,cosAlpha in the Global System
814 std::array<value_t, 9> vecLab{0.f};
815 if (!linRef0.getPosDirGlo(vecLab)) {
816 return false;
817 }
818 //
819 // Rotate to the system where Bx=By=0.
820 value_t bxy2 = b[0] * b[0] + b[1] * b[1];
821 value_t bt = gpu::CAMath::Sqrt(bxy2);
822 value_t cosphi = 1.f, sinphi = 0.f;
823 if (bt > constants::math::Almost0) {
824 cosphi = b[0] / bt;
825 sinphi = b[1] / bt;
826 }
827 value_t bb = gpu::CAMath::Sqrt(bxy2 + b[2] * b[2]);
828 value_t costet = 1., sintet = 0.;
829 if (bb > constants::math::Almost0) {
830 costet = b[2] / bb;
831 sintet = bt / bb;
832 }
833 std::array<value_t, 7> vect{costet * cosphi * vecLab[0] + costet * sinphi * vecLab[1] - sintet * vecLab[2],
834 -sinphi * vecLab[0] + cosphi * vecLab[1],
835 sintet * cosphi * vecLab[0] + sintet * sinphi * vecLab[1] + costet * vecLab[2],
836 costet * cosphi * vecLab[3] + costet * sinphi * vecLab[4] - sintet * vecLab[5],
837 -sinphi * vecLab[3] + cosphi * vecLab[4],
838 sintet * cosphi * vecLab[3] + sintet * sinphi * vecLab[4] + costet * vecLab[5],
839 vecLab[6]};
840
841 // Do the helix step
842 value_t q = this->getCharge();
843 g3helx3(q * bb, step, vect);
844
845 // Rotate back to the Global System
846 vecLab[0] = cosphi * costet * vect[0] - sinphi * vect[1] + cosphi * sintet * vect[2];
847 vecLab[1] = sinphi * costet * vect[0] + cosphi * vect[1] + sinphi * sintet * vect[2];
848 vecLab[2] = -sintet * vect[0] + costet * vect[2];
849
850 vecLab[3] = cosphi * costet * vect[3] - sinphi * vect[4] + cosphi * sintet * vect[5];
851 vecLab[4] = sinphi * costet * vect[3] + cosphi * vect[4] + sinphi * sintet * vect[5];
852 vecLab[5] = -sintet * vect[3] + costet * vect[5];
853
854 // Rotate back to the Tracking System
855 value_t sinalp = -vecLab[7], cosalp = vecLab[8];
856 value_t t = cosalp * vecLab[0] - sinalp * vecLab[1];
857 vecLab[1] = sinalp * vecLab[0] + cosalp * vecLab[1];
858 vecLab[0] = t;
859 t = cosalp * vecLab[3] - sinalp * vecLab[4];
860 vecLab[4] = sinalp * vecLab[3] + cosalp * vecLab[4];
861 vecLab[3] = t;
862
863 // Do the final correcting step to the target plane (linear approximation)
864 value_t x = vecLab[0], y = vecLab[1], z = vecLab[2];
865 if (gpu::CAMath::Abs(x - xk) > constants::math::Almost0) {
866 if (gpu::CAMath::Abs(vecLab[3]) < constants::math::Almost0) {
867 return false;
868 }
869 auto dxFin = xk - vecLab[0];
870 x += dxFin;
871 y += vecLab[4] / vecLab[3] * dxFin;
872 z += vecLab[5] / vecLab[3] * dxFin;
873 }
874
875 // Calculate the track parameters
876 auto linRef1 = linRef0;
877 t = 1.f / gpu::CAMath::Sqrt(vecLab[3] * vecLab[3] + vecLab[4] * vecLab[4]);
878 linRef1.setX(xk);
879 linRef1.setY(y);
880 linRef1.setZ(z);
881 linRef1.setSnp(snpRef1 = vecLab[4] * t); // reassign snpRef1
882 linRef1.setTgl(vecLab[5] * t);
883 linRef1.setQ2Pt(q * t / vecLab[6]);
884
885 // recalculate parameters of the transported ref track needed for transport of this:
886 cspRef1 = gpu::CAMath::Sqrt((1 - snpRef1) * (1 + snpRef1));
887 cspRef1Inv = value_t(1) / cspRef1;
888 cc = cspRef0 + cspRef1;
889 ccInv = value_t(1) / cc;
890 dy2dx = (snpRef0 + snpRef1) * ccInv;
891 double dxccInv = dx * ccInv, hh = dxccInv * cspRef1Inv * (1 + cspRef0 * cspRef1 + snpRef0 * snpRef1), jj = dx * (dy2dx - snpRef1 * cspRef1Inv);
892 double f02 = hh * cspRef0Inv;
893 double f04 = hh * dxccInv * kb;
894 double f24 = dx * kb;
895 double f12 = linRef0.getTgl() * (f02 * snpRef1 + jj);
896 double f13 = dx * (cspRef1 + snpRef1 * dy2dx); // dS
897 double f14 = linRef0.getTgl() * (f04 * snpRef1 + jj * f24);
898
899 // difference between the current and reference state
900 value_t diff[5];
901 for (int i = 0; i < 5; i++) {
902 diff[i] = this->getParam(i) - linRef0.getParam(i);
903 }
904 value_t snpUpd = snpRef1 + diff[kSnp] + f24 * diff[kQ2Pt];
905 if (gpu::CAMath::Abs(snpUpd) > constants::math::Almost1) {
906 return false;
907 }
908 this->setX(xk);
909 this->setY(linRef1.getY() + diff[kY] + f02 * diff[kSnp] + f04 * diff[kQ2Pt]);
910 this->setZ(linRef1.getZ() + diff[kZ] + f13 * diff[kTgl] + f14 * diff[kQ2Pt]);
911 this->setSnp(snpUpd);
912 this->setTgl(linRef1.getTgl() + diff[kTgl]);
913 this->setQ2Pt(linRef1.getQ2Pt() + diff[kQ2Pt]);
914
915 linRef0 = linRef1; // update reference track
916
917 // matrix transformed with Bz component only
918 value_t &c00 = mC[kSigY2], &c10 = mC[kSigZY], &c11 = mC[kSigZ2], &c20 = mC[kSigSnpY], &c21 = mC[kSigSnpZ],
919 &c22 = mC[kSigSnp2], &c30 = mC[kSigTglY], &c31 = mC[kSigTglZ], &c32 = mC[kSigTglSnp], &c33 = mC[kSigTgl2],
920 &c40 = mC[kSigQ2PtY], &c41 = mC[kSigQ2PtZ], &c42 = mC[kSigQ2PtSnp], &c43 = mC[kSigQ2PtTgl],
921 &c44 = mC[kSigQ2Pt2];
922
923 // b = C*ft
924 double b00 = f02 * c20 + f04 * c40, b01 = f12 * c20 + f14 * c40 + f13 * c30;
925 double b02 = f24 * c40;
926 double b10 = f02 * c21 + f04 * c41, b11 = f12 * c21 + f14 * c41 + f13 * c31;
927 double b12 = f24 * c41;
928 double b20 = f02 * c22 + f04 * c42, b21 = f12 * c22 + f14 * c42 + f13 * c32;
929 double b22 = f24 * c42;
930 double b40 = f02 * c42 + f04 * c44, b41 = f12 * c42 + f14 * c44 + f13 * c43;
931 double b42 = f24 * c44;
932 double b30 = f02 * c32 + f04 * c43, b31 = f12 * c32 + f14 * c43 + f13 * c33;
933 double b32 = f24 * c43;
934
935 // a = f*b = f*C*ft
936 double a00 = f02 * b20 + f04 * b40, a01 = f02 * b21 + f04 * b41, a02 = f02 * b22 + f04 * b42;
937 double a11 = f12 * b21 + f14 * b41 + f13 * b31, a12 = f12 * b22 + f14 * b42 + f13 * b32;
938 double a22 = f24 * b42;
939
940 // F*C*Ft = C + (b + bt + a)
941 c00 += b00 + b00 + a00;
942 c10 += b10 + b01 + a01;
943 c20 += b20 + b02 + a02;
944 c30 += b30;
945 c40 += b40;
946 c11 += b11 + b11 + a11;
947 c21 += b21 + b12 + a12;
948 c31 += b31;
949 c41 += b41;
950 c22 += b22 + b22 + a22;
951 c32 += b32;
952 c42 += b42;
953
954 checkCovariance();
955
956 return true;
957}
958
959//______________________________________________
960template <typename value_T>
961GPUd() void TrackParametrizationWithError<value_T>::checkCorrelations()
962{
963 // This function forces the abs of correlation coefficients to be <1.
964 constexpr float MaxCorr = 0.99;
965 for (int i = kNParams; i--;) {
966 for (int j = i; j--;) {
967 auto sig2 = mC[DiagMap[i]] * mC[DiagMap[j]];
968 auto& cov = mC[CovarMap[i][j]];
969 if (cov * cov >= MaxCorr * sig2) { // constrain correlation
970 cov = gpu::CAMath::Sqrt(sig2) * (cov > 0. ? MaxCorr : -MaxCorr);
971 }
972 }
973 }
974}
975
976//______________________________________________
977template <typename value_T>
978GPUd() void TrackParametrizationWithError<value_T>::checkCovariance()
979{
980 // This function forces the diagonal elements of the covariance matrix to be positive and abs of correlation coefficients to be <1.
981 // In case the diagonal element is bigger than the maximal allowed value, it is set to
982 // the limit and the off-diagonal elements that correspond to it are set to zero.
983
984 mC[kSigY2] = gpu::CAMath::Abs(mC[kSigY2]);
985 if (mC[kSigY2] > kCY2max) {
986 value_t scl = gpu::CAMath::Sqrt(kCY2max / mC[kSigY2]);
987 mC[kSigY2] = kCY2max;
988 mC[kSigZY] *= scl;
989 mC[kSigSnpY] *= scl;
990 mC[kSigTglY] *= scl;
991 mC[kSigQ2PtY] *= scl;
992 }
993 mC[kSigZ2] = gpu::CAMath::Abs(mC[kSigZ2]);
994 if (mC[kSigZ2] > kCZ2max) {
995 value_t scl = gpu::CAMath::Sqrt(kCZ2max / mC[kSigZ2]);
996 mC[kSigZ2] = kCZ2max;
997 mC[kSigZY] *= scl;
998 mC[kSigSnpZ] *= scl;
999 mC[kSigTglZ] *= scl;
1000 mC[kSigQ2PtZ] *= scl;
1001 }
1002 mC[kSigSnp2] = gpu::CAMath::Abs(mC[kSigSnp2]);
1003 if (mC[kSigSnp2] > kCSnp2max) {
1004 value_t scl = gpu::CAMath::Sqrt(kCSnp2max / mC[kSigSnp2]);
1005 mC[kSigSnp2] = kCSnp2max;
1006 mC[kSigSnpY] *= scl;
1007 mC[kSigSnpZ] *= scl;
1008 mC[kSigTglSnp] *= scl;
1009 mC[kSigQ2PtSnp] *= scl;
1010 }
1011 mC[kSigTgl2] = gpu::CAMath::Abs(mC[kSigTgl2]);
1012 if (mC[kSigTgl2] > kCTgl2max) {
1013 value_t scl = gpu::CAMath::Sqrt(kCTgl2max / mC[kSigTgl2]);
1014 mC[kSigTgl2] = kCTgl2max;
1015 mC[kSigTglY] *= scl;
1016 mC[kSigTglZ] *= scl;
1017 mC[kSigTglSnp] *= scl;
1018 mC[kSigQ2PtTgl] *= scl;
1019 }
1020 mC[kSigQ2Pt2] = gpu::CAMath::Abs(mC[kSigQ2Pt2]);
1021 if (mC[kSigQ2Pt2] > kC1Pt2max) {
1022 value_t scl = gpu::CAMath::Sqrt(kC1Pt2max / mC[kSigQ2Pt2]);
1023 mC[kSigQ2Pt2] = kC1Pt2max;
1024 mC[kSigQ2PtY] *= scl;
1025 mC[kSigQ2PtZ] *= scl;
1026 mC[kSigQ2PtSnp] *= scl;
1027 mC[kSigQ2PtTgl] *= scl;
1028 }
1029}
1030
1031//______________________________________________
1032template <typename value_T>
1033GPUd() void TrackParametrizationWithError<value_T>::resetCovariance(value_t s2)
1034{
1035 // Reset the covarince matrix to "something big"
1036 double d0(kCY2max), d1(kCZ2max), d2(kCSnp2max), d3(kCTgl2max), d4(kC1Pt2max);
1037 if (s2 > constants::math::Almost0) {
1038 d0 = getSigmaY2() * s2;
1039 d1 = getSigmaZ2() * s2;
1040 d2 = getSigmaSnp2() * s2;
1041 d3 = getSigmaTgl2() * s2;
1042 d4 = getSigma1Pt2() * s2;
1043 if (d0 > kCY2max) {
1044 d0 = kCY2max;
1045 }
1046 if (d1 > kCZ2max) {
1047 d1 = kCZ2max;
1048 }
1049 if (d2 > kCSnp2max) {
1050 d2 = kCSnp2max;
1051 }
1052 if (d3 > kCTgl2max) {
1053 d3 = kCTgl2max;
1054 }
1055 if (d4 > kC1Pt2max) {
1056 d4 = kC1Pt2max;
1057 }
1058 }
1059 for (int i = 0; i < kCovMatSize; i++) {
1060 mC[i] = 0;
1061 }
1062 mC[kSigY2] = d0;
1063 mC[kSigZ2] = d1;
1064 mC[kSigSnp2] = d2;
1065 mC[kSigTgl2] = d3;
1066 mC[kSigQ2Pt2] = d4;
1067}
1068
1069//______________________________________________
1070template <typename value_T>
1071GPUd() auto TrackParametrizationWithError<value_T>::getPredictedChi2(const value_t* p, const value_t* cov) const -> value_t
1072{
1073 // Estimate the chi2 of the space point "p" with the cov. matrix "cov"
1074 auto sdd = static_cast<double>(getSigmaY2()) + static_cast<double>(cov[0]);
1075 auto sdz = static_cast<double>(getSigmaZY()) + static_cast<double>(cov[1]);
1076 auto szz = static_cast<double>(getSigmaZ2()) + static_cast<double>(cov[2]);
1077 auto det = sdd * szz - sdz * sdz;
1078
1079 if (gpu::CAMath::Abs(det) < constants::math::Almost0) {
1080 return constants::math::VeryBig;
1081 }
1082
1083 value_t d = this->getY() - p[0];
1084 value_t z = this->getZ() - p[1];
1085 auto chi2 = (d * (szz * d - sdz * z) + z * (sdd * z - d * sdz)) / det;
1086 if (chi2 < 0.) {
1087#ifndef GPUCA_ALIGPUCODE
1088 LOGP(warning, "Negative chi2={}, Cluster: {} {} {} Dy:{} Dz:{} | sdd:{} sdz:{} szz:{} det:{}", chi2, cov[0], cov[1], cov[2], d, z, sdd, sdz, szz, det);
1089 LOGP(warning, "Track: {}", asString());
1090#endif
1091 }
1092 return chi2;
1093}
1094
1095//______________________________________________
1096template <typename value_T>
1097GPUd() auto TrackParametrizationWithError<value_T>::getPredictedChi2Quiet(const value_t* p, const value_t* cov) const -> value_t
1098{
1099 // Estimate the chi2 of the space point "p" with the cov. matrix "cov"
1100 auto sdd = static_cast<double>(getSigmaY2()) + static_cast<double>(cov[0]);
1101 auto sdz = static_cast<double>(getSigmaZY()) + static_cast<double>(cov[1]);
1102 auto szz = static_cast<double>(getSigmaZ2()) + static_cast<double>(cov[2]);
1103 auto det = sdd * szz - sdz * sdz;
1104
1105 if (gpu::CAMath::Abs(det) < constants::math::Almost0) {
1106 return constants::math::VeryBig;
1107 }
1108
1109 value_t d = this->getY() - p[0];
1110 value_t z = this->getZ() - p[1];
1111
1112 return (d * (szz * d - sdz * z) + z * (sdd * z - d * sdz)) / det;
1113}
1114
1115//______________________________________________
1116template <typename value_T>
1117GPUd() auto TrackParametrizationWithError<value_T>::getPredictedChi2(const TrackParametrizationWithError<value_T>& rhs) const -> value_t
1118{
1119 MatrixDSym5 cov; // perform matrix operations in double!
1120 return getPredictedChi2(rhs, cov);
1121}
1122
1123//______________________________________________
1124template <typename value_T>
1125GPUd() void TrackParametrizationWithError<value_T>::buildCombinedCovMatrix(const TrackParametrizationWithError<value_T>& rhs, MatrixDSym5& cov) const
1126{
1127 // fill combined cov.matrix (NOT inverted)
1128 cov(kY, kY) = static_cast<double>(getSigmaY2()) + static_cast<double>(rhs.getSigmaY2());
1129 cov(kZ, kY) = static_cast<double>(getSigmaZY()) + static_cast<double>(rhs.getSigmaZY());
1130 cov(kZ, kZ) = static_cast<double>(getSigmaZ2()) + static_cast<double>(rhs.getSigmaZ2());
1131 cov(kSnp, kY) = static_cast<double>(getSigmaSnpY()) + static_cast<double>(rhs.getSigmaSnpY());
1132 cov(kSnp, kZ) = static_cast<double>(getSigmaSnpZ()) + static_cast<double>(rhs.getSigmaSnpZ());
1133 cov(kSnp, kSnp) = static_cast<double>(getSigmaSnp2()) + static_cast<double>(rhs.getSigmaSnp2());
1134 cov(kTgl, kY) = static_cast<double>(getSigmaTglY()) + static_cast<double>(rhs.getSigmaTglY());
1135 cov(kTgl, kZ) = static_cast<double>(getSigmaTglZ()) + static_cast<double>(rhs.getSigmaTglZ());
1136 cov(kTgl, kSnp) = static_cast<double>(getSigmaTglSnp()) + static_cast<double>(rhs.getSigmaTglSnp());
1137 cov(kTgl, kTgl) = static_cast<double>(getSigmaTgl2()) + static_cast<double>(rhs.getSigmaTgl2());
1138 cov(kQ2Pt, kY) = static_cast<double>(getSigma1PtY()) + static_cast<double>(rhs.getSigma1PtY());
1139 cov(kQ2Pt, kZ) = static_cast<double>(getSigma1PtZ()) + static_cast<double>(rhs.getSigma1PtZ());
1140 cov(kQ2Pt, kSnp) = static_cast<double>(getSigma1PtSnp()) + static_cast<double>(rhs.getSigma1PtSnp());
1141 cov(kQ2Pt, kTgl) = static_cast<double>(getSigma1PtTgl()) + static_cast<double>(rhs.getSigma1PtTgl());
1142 cov(kQ2Pt, kQ2Pt) = static_cast<double>(getSigma1Pt2()) + static_cast<double>(rhs.getSigma1Pt2());
1143}
1144
1145//______________________________________________
1146template <typename value_T>
1147GPUd() auto TrackParametrizationWithError<value_T>::getPredictedChi2(const TrackParametrizationWithError<value_T>& rhs, MatrixDSym5& covToSet) const -> value_t
1148{
1149 // get chi2 wrt other track, which must be defined at the same parameters X,alpha
1150 // Supplied non-initialized covToSet matrix is filled by inverse combined matrix for further use
1151
1152 if (gpu::CAMath::Abs(this->getAlpha() - rhs.getAlpha()) > o2::constants::math::Epsilon) {
1153 LOG(error) << "The reference Alpha of the tracks differ: " << this->getAlpha() << " : " << rhs.getAlpha();
1154 return 2.f * HugeF;
1155 }
1156 if (gpu::CAMath::Abs(this->getX() - rhs.getX()) > o2::constants::math::Epsilon) {
1157 LOG(error) << "The reference X of the tracks differ: " << this->getX() << " : " << rhs.getX();
1158 return 2.f * HugeF;
1159 }
1160 buildCombinedCovMatrix(rhs, covToSet);
1161 if (!covToSet.Invert()) {
1162 LOG(warning) << "Cov.matrix inversion failed: " << covToSet;
1163 return 2.f * HugeF;
1164 }
1165 double chi2diag = 0., chi2ndiag = 0., diff[kNParams];
1166 for (int i = kNParams; i--;) {
1167 diff[i] = this->getParam(i) - rhs.getParam(i);
1168 chi2diag += diff[i] * diff[i] * covToSet(i, i);
1169 }
1170 for (int i = kNParams; i--;) {
1171 for (int j = i; j--;) {
1172 chi2ndiag += diff[i] * diff[j] * covToSet(i, j);
1173 }
1174 }
1175 return chi2diag + 2. * chi2ndiag;
1176}
1177
1178//______________________________________________
1179template <typename value_T>
1180GPUd() bool TrackParametrizationWithError<value_T>::update(const TrackParametrizationWithError<value_T>& rhs, const MatrixDSym5& covInv)
1181{
1182 // update track with other track, the inverted combined cov matrix should be supplied
1183
1184 // consider skipping this check, since it is usually already done upstream
1185 if (gpu::CAMath::Abs(this->getAlpha() - rhs.getAlpha()) > o2::constants::math::Epsilon) {
1186 LOG(error) << "The reference Alpha of the tracks differ: " << this->getAlpha() << " : " << rhs.getAlpha();
1187 return false;
1188 }
1189 if (gpu::CAMath::Abs(this->getX() - rhs.getX()) > o2::constants::math::Epsilon) {
1190 LOG(error) << "The reference X of the tracks differ: " << this->getX() << " : " << rhs.getX();
1191 return false;
1192 }
1193
1194 // gain matrix K = Cov0*H*(Cov0+Cov0)^-1 (for measurement matrix H=I)
1195 MatrixDSym5 matC0;
1196 matC0(kY, kY) = getSigmaY2();
1197 matC0(kZ, kY) = getSigmaZY();
1198 matC0(kZ, kZ) = getSigmaZ2();
1199 matC0(kSnp, kY) = getSigmaSnpY();
1200 matC0(kSnp, kZ) = getSigmaSnpZ();
1201 matC0(kSnp, kSnp) = getSigmaSnp2();
1202 matC0(kTgl, kY) = getSigmaTglY();
1203 matC0(kTgl, kZ) = getSigmaTglZ();
1204 matC0(kTgl, kSnp) = getSigmaTglSnp();
1205 matC0(kTgl, kTgl) = getSigmaTgl2();
1206 matC0(kQ2Pt, kY) = getSigma1PtY();
1207 matC0(kQ2Pt, kZ) = getSigma1PtZ();
1208 matC0(kQ2Pt, kSnp) = getSigma1PtSnp();
1209 matC0(kQ2Pt, kTgl) = getSigma1PtTgl();
1210 matC0(kQ2Pt, kQ2Pt) = getSigma1Pt2();
1211 MatrixD5 matK = matC0 * covInv;
1212
1213 // updated state vector: x = K*(x1-x0)
1214 // RS: why SMatix, SVector does not provide multiplication operators ???
1215 double diff[kNParams];
1216 for (int i = kNParams; i--;) {
1217 diff[i] = rhs.getParam(i) - this->getParam(i);
1218 }
1219 for (int i = kNParams; i--;) {
1220 for (int j = kNParams; j--;) {
1221 this->updateParam(matK(i, j) * diff[j], i);
1222 }
1223 }
1224
1225 // updated covariance: Cov0 = Cov0 - K*Cov0
1226 matK *= o2::math_utils::SMatrix<double, kNParams, kNParams, o2::math_utils::MatRepStd<double, kNParams>>(matC0);
1227 mC[kSigY2] -= matK(kY, kY);
1228 mC[kSigZY] -= matK(kZ, kY);
1229 mC[kSigZ2] -= matK(kZ, kZ);
1230 mC[kSigSnpY] -= matK(kSnp, kY);
1231 mC[kSigSnpZ] -= matK(kSnp, kZ);
1232 mC[kSigSnp2] -= matK(kSnp, kSnp);
1233 mC[kSigTglY] -= matK(kTgl, kY);
1234 mC[kSigTglZ] -= matK(kTgl, kZ);
1235 mC[kSigTglSnp] -= matK(kTgl, kSnp);
1236 mC[kSigTgl2] -= matK(kTgl, kTgl);
1237 mC[kSigQ2PtY] -= matK(kQ2Pt, kY);
1238 mC[kSigQ2PtZ] -= matK(kQ2Pt, kZ);
1239 mC[kSigQ2PtSnp] -= matK(kQ2Pt, kSnp);
1240 mC[kSigQ2PtTgl] -= matK(kQ2Pt, kTgl);
1241 mC[kSigQ2Pt2] -= matK(kQ2Pt, kQ2Pt);
1242
1243 return true;
1244}
1245
1246//______________________________________________
1247template <typename value_T>
1248GPUd() bool TrackParametrizationWithError<value_T>::update(const TrackParametrizationWithError<value_T>& rhs)
1249{
1250 // update track with other track
1251 MatrixDSym5 covI; // perform matrix operations in double!
1252 buildCombinedCovMatrix(rhs, covI);
1253 if (!covI.Invert()) {
1254 LOG(warning) << "Cov.matrix inversion failed: " << covI;
1255 return false;
1256 }
1257 return update(rhs, covI);
1258}
1259
1260//______________________________________________
1261template <typename value_T>
1262GPUd() bool TrackParametrizationWithError<value_T>::update(const value_t* p, const value_t* cov)
1263{
1264 // Update the track parameters with the space point "p" having
1265 // the covariance matrix "cov"
1266
1267 value_t &cm00 = mC[kSigY2], &cm10 = mC[kSigZY], &cm11 = mC[kSigZ2], &cm20 = mC[kSigSnpY], &cm21 = mC[kSigSnpZ],
1268 &cm22 = mC[kSigSnp2], &cm30 = mC[kSigTglY], &cm31 = mC[kSigTglZ], &cm32 = mC[kSigTglSnp], &cm33 = mC[kSigTgl2],
1269 &cm40 = mC[kSigQ2PtY], &cm41 = mC[kSigQ2PtZ], &cm42 = mC[kSigQ2PtSnp], &cm43 = mC[kSigQ2PtTgl],
1270 &cm44 = mC[kSigQ2Pt2];
1271
1272 // use double precision?
1273 double r00 = static_cast<double>(cov[0]) + static_cast<double>(cm00);
1274 double r01 = static_cast<double>(cov[1]) + static_cast<double>(cm10);
1275 double r11 = static_cast<double>(cov[2]) + static_cast<double>(cm11);
1276 double det = r00 * r11 - r01 * r01;
1277
1278 if (gpu::CAMath::Abs(det) < constants::math::Almost0) {
1279 return false;
1280 }
1281 double detI = 1. / det;
1282 double tmp = r00;
1283 r00 = r11 * detI;
1284 r11 = tmp * detI;
1285 r01 = -r01 * detI;
1286
1287 double k00 = cm00 * r00 + cm10 * r01, k01 = cm00 * r01 + cm10 * r11;
1288 double k10 = cm10 * r00 + cm11 * r01, k11 = cm10 * r01 + cm11 * r11;
1289 double k20 = cm20 * r00 + cm21 * r01, k21 = cm20 * r01 + cm21 * r11;
1290 double k30 = cm30 * r00 + cm31 * r01, k31 = cm30 * r01 + cm31 * r11;
1291 double k40 = cm40 * r00 + cm41 * r01, k41 = cm40 * r01 + cm41 * r11;
1292
1293 value_t dy = p[kY] - this->getY(), dz = p[kZ] - this->getZ();
1294 value_t dsnp = k20 * dy + k21 * dz;
1295 if (gpu::CAMath::Abs(this->getSnp() + dsnp) > constants::math::Almost1) {
1296 return false;
1297 }
1298
1299 const params_t dP{value_t(k00 * dy + k01 * dz), value_t(k10 * dy + k11 * dz), dsnp, value_t(k30 * dy + k31 * dz),
1300 value_t(k40 * dy + k41 * dz)};
1301 this->updateParams(dP);
1302
1303 double c01 = cm10, c02 = cm20, c03 = cm30, c04 = cm40;
1304 double c12 = cm21, c13 = cm31, c14 = cm41;
1305
1306 cm00 -= k00 * cm00 + k01 * cm10;
1307 cm10 -= k00 * c01 + k01 * cm11;
1308 cm20 -= k00 * c02 + k01 * c12;
1309 cm30 -= k00 * c03 + k01 * c13;
1310 cm40 -= k00 * c04 + k01 * c14;
1311
1312 cm11 -= k10 * c01 + k11 * cm11;
1313 cm21 -= k10 * c02 + k11 * c12;
1314 cm31 -= k10 * c03 + k11 * c13;
1315 cm41 -= k10 * c04 + k11 * c14;
1316
1317 cm22 -= k20 * c02 + k21 * c12;
1318 cm32 -= k20 * c03 + k21 * c13;
1319 cm42 -= k20 * c04 + k21 * c14;
1320
1321 cm33 -= k30 * c03 + k31 * c13;
1322 cm43 -= k30 * c04 + k31 * c14;
1323
1324 cm44 -= k40 * c04 + k41 * c14;
1325
1326 checkCovariance();
1327
1328 return true;
1329}
1330
1331//______________________________________________
1332template <typename value_T>
1333GPUd() value_T TrackParametrizationWithError<value_T>::update(const o2::dataformats::VertexBase& vtx, value_T maxChi2)
1334{
1335 // Update track with vertex if the track-vertex chi2 does not exceed maxChi2. Track must be already propagated to the DCA to vertex
1336 // return update chi2 or -chi2 if chi2 if chi2 exceeds maxChi2
1337 auto vtLoc = this->getVertexInTrackFrame(vtx);
1338 value_T chi2 = getPredictedChi2(vtLoc.yz, vtLoc.yzerr);
1339 return chi2 < maxChi2 && update(vtLoc.yz, vtLoc.yzerr) ? chi2 : -chi2;
1340}
1341
1342//______________________________________________
1343template <typename value_T>
1344GPUd() bool TrackParametrizationWithError<value_T>::correctForMaterial(value_t x2x0, value_t xrho, bool anglecorr)
1345{
1346 //------------------------------------------------------------------
1347 // This function corrects the track parameters for the crossed material.
1348 // "x2x0" - X/X0, the thickness in units of the radiation length.
1349 // "xrho" - is the product length*density (g/cm^2).
1350 // It should be passed as negative when propagating tracks
1351 // from the intreaction point to the outside of the central barrel.
1352 // "dedx" - mean enery loss (GeV/(g/cm^2), if <=kCalcdEdxAuto : calculate on the fly
1353 // "anglecorr" - switch for the angular correction
1354 //------------------------------------------------------------------
1355 constexpr value_t kMSConst2 = 0.0136f * 0.0136f;
1356 constexpr value_t kMinP = 0.01f; // kill below this momentum
1357
1358 value_t csp2 = (1.f - this->getSnp()) * (1.f + this->getSnp()); // cos(phi)^2
1359 value_t cst2I = (1.f + this->getTgl() * this->getTgl()); // 1/cos(lambda)^2
1360 if (anglecorr) { // Apply angle correction, if requested
1361 value_t angle = gpu::CAMath::Sqrt(cst2I / csp2);
1362 x2x0 *= angle;
1363 xrho *= angle;
1364 }
1365 auto m = this->getPID().getMass();
1366 int charge2 = this->getAbsCharge() * this->getAbsCharge();
1367 value_t p = this->getP(), p0 = p, p02 = p * p, e2 = p02 + this->getPID().getMass2(), massInv = 1. / m, bg = p * massInv, dETot = 0.;
1368 value_t e = gpu::CAMath::Sqrt(e2), e0 = e;
1369 if (m > 0 && xrho != 0.f) {
1370 value_t ekin = e - m, dedx = this->getdEdxBBOpt(bg);
1371#ifdef _BB_NONCONST_CORR_
1372 value_t dedxDer = 0., dedx1 = dedx;
1373#endif
1374 if (charge2 != 1) {
1375 dedx *= charge2;
1376 }
1377 value_t dE = dedx * xrho;
1378 int na = 1 + int(gpu::CAMath::Abs(dE) / ekin * ELoss2EKinThreshInv);
1379 if (na > MaxELossIter) {
1380 na = MaxELossIter;
1381 }
1382 if (na > 1) {
1383 dE /= na;
1384 xrho /= na;
1385#ifdef _BB_NONCONST_CORR_
1386 dedxDer = this->getBetheBlochSolidDerivativeApprox(dedx1, bg); // require correction for non-constantness of dedx vs betagamma
1387 if (charge2 != 1) {
1388 dedxDer *= charge2;
1389 }
1390#endif
1391 }
1392 while (na--) {
1393#ifdef _BB_NONCONST_CORR_
1394 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]
1395 if (xrho < 0) {
1396 dedxDer = -dedxDer; // E.loss ( -> positive derivative)
1397 }
1398 auto corrC = (gpu::CAMath::Exp(dedxDer) - 1.) / dedxDer;
1399 dE *= corrC;
1400 }
1401#endif
1402 e += dE;
1403 if (e > m) { // stopped
1404 p = gpu::CAMath::Sqrt(e * e - this->getPID().getMass2());
1405 } else {
1406 return false;
1407 }
1408 if (na) {
1409 bg = p * massInv;
1410 dedx = this->getdEdxBBOpt(bg);
1411#ifdef _BB_NONCONST_CORR_
1412 dedxDer = this->getBetheBlochSolidDerivativeApprox(dedx, bg);
1413#endif
1414 if (charge2 != 1) {
1415 dedx *= charge2;
1416#ifdef _BB_NONCONST_CORR_
1417 dedxDer *= charge2;
1418#endif
1419 }
1420 dE = dedx * xrho;
1421 }
1422 }
1423
1424 if (p < kMinP) {
1425 return false;
1426 }
1427 dETot = e - e0;
1428 } // end of e.loss correction
1429
1430 // Calculating the multiple scattering corrections******************
1431 value_t& fC22 = mC[kSigSnp2];
1432 value_t& fC33 = mC[kSigTgl2];
1433 value_t& fC43 = mC[kSigQ2PtTgl];
1434 value_t& fC44 = mC[kSigQ2Pt2];
1435 //
1436 value_t cC22(0.f), cC33(0.f), cC43(0.f), cC44(0.f);
1437 if (x2x0 != 0.f) {
1438 value_t beta2 = p02 / e2, theta2 = kMSConst2 / (beta2 * p02) * gpu::CAMath::Abs(x2x0);
1439 value_t fp34 = this->getTgl();
1440 if (charge2 != 1) {
1441 theta2 *= charge2;
1442 fp34 *= this->getCharge2Pt();
1443 }
1444 if (theta2 > constants::math::PI * constants::math::PI) {
1445 return false;
1446 }
1447 value_t t2c2I = theta2 * cst2I;
1448 cC22 = t2c2I * csp2;
1449 cC33 = t2c2I * cst2I;
1450 cC43 = t2c2I * fp34;
1451 cC44 = theta2 * fp34 * fp34;
1452 // optimize this
1453 // cC22 = theta2*((1.-getSnp())*(1.+getSnp()))*(1. + this->getTgl()*getTgl());
1454 // cC33 = theta2*(1. + this->getTgl()*getTgl())*(1. + this->getTgl()*getTgl());
1455 // cC43 = theta2*getTgl()*this->getQ2Pt()*(1. + this->getTgl()*getTgl());
1456 // cC44 = theta2*getTgl()*this->getQ2Pt()*getTgl()*this->getQ2Pt();
1457 }
1458
1459 // the energy loss correction contribution to cov.matrix: approximate energy loss fluctuation (M.Ivanov)
1460 constexpr value_t knst = 0.0007f; // To be tuned.
1461 value_t sigmadE = knst * gpu::CAMath::Sqrt(gpu::CAMath::Abs(dETot)) * e0 / p02 * this->getCharge2Pt();
1462 cC44 += sigmadE * sigmadE;
1463
1464 // Applying the corrections*****************************
1465 fC22 += cC22;
1466 fC33 += cC33;
1467 fC43 += cC43;
1468 fC44 += cC44;
1469 this->setQ2Pt(this->getQ2Pt() * p0 / p);
1470
1471 checkCovariance();
1472
1473 return true;
1474}
1475
1476//______________________________________________
1477template <typename value_T>
1478GPUd() bool TrackParametrizationWithError<value_T>::correctForMaterial(TrackParametrization<value_T>& linRef, value_t x2x0, value_t xrho, bool anglecorr)
1479{
1480 //------------------------------------------------------------------
1481 // This function corrects the reference and current track parameters for the crossed material
1482 // "x2x0" - X/X0, the thickness in units of the radiation length.
1483 // "xrho" - is the product length*density (g/cm^2).
1484 // It should be passed as negative when propagating tracks
1485 // from the intreaction point to the outside of the central barrel.
1486 // "dedx" - mean enery loss (GeV/(g/cm^2), if <=kCalcdEdxAuto : calculate on the fly
1487 // "anglecorr" - switch for the angular correction
1488 //------------------------------------------------------------------
1489 constexpr value_t kMSConst2 = 0.0136f * 0.0136f;
1490 constexpr value_t kMinP = 0.01f; // kill below this momentum
1491
1492 value_t csp2 = (1.f - linRef.getSnp()) * (1.f + linRef.getSnp()); // cos(phi)^2
1493 value_t cst2I = (1.f + linRef.getTgl() * linRef.getTgl()); // 1/cos(lambda)^2
1494 if (anglecorr) { // Apply angle correction, if requested
1495 value_t angle = gpu::CAMath::Sqrt(cst2I / csp2);
1496 x2x0 *= angle;
1497 xrho *= angle;
1498 }
1499 auto pid = linRef.getPID();
1500 auto m = pid.getMass();
1501 int charge2 = linRef.getAbsCharge() * linRef.getAbsCharge();
1502 value_t p = linRef.getP(), p0 = p, p02 = p * p, e2 = p02 + pid.getMass2(), massInv = 1. / m, bg = p * massInv, dETot = 0.;
1503 value_t e = gpu::CAMath::Sqrt(e2), e0 = e;
1504 if (m > 0 && xrho != 0.f) {
1505 value_t ekin = e - m, dedx = this->getdEdxBBOpt(bg);
1506#ifdef _BB_NONCONST_CORR_
1507 value_t dedxDer = 0., dedx1 = dedx;
1508#endif
1509 if (charge2 != 1) {
1510 dedx *= charge2;
1511 }
1512 value_t dE = dedx * xrho;
1513 int na = 1 + int(gpu::CAMath::Abs(dE) / ekin * ELoss2EKinThreshInv);
1514 if (na > MaxELossIter) {
1515 na = MaxELossIter;
1516 }
1517 if (na > 1) {
1518 dE /= na;
1519 xrho /= na;
1520#ifdef _BB_NONCONST_CORR_
1521 dedxDer = this->getBetheBlochSolidDerivativeApprox(dedx1, bg); // require correction for non-constantness of dedx vs betagamma
1522 if (charge2 != 1) {
1523 dedxDer *= charge2;
1524 }
1525#endif
1526 }
1527 while (na--) {
1528#ifdef _BB_NONCONST_CORR_
1529 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]
1530 if (xrho < 0) {
1531 dedxDer = -dedxDer; // E.loss ( -> positive derivative)
1532 }
1533 auto corrC = (gpu::CAMath::Exp(dedxDer) - 1.) / dedxDer;
1534 dE *= corrC;
1535 }
1536#endif
1537 e += dE;
1538 if (e > m) { // stopped
1539 p = gpu::CAMath::Sqrt(e * e - pid.getMass2());
1540 } else {
1541 return false;
1542 }
1543 if (na) {
1544 bg = p * massInv;
1545 dedx = this->getdEdxBBOpt(bg);
1546#ifdef _BB_NONCONST_CORR_
1547 dedxDer = this->getBetheBlochSolidDerivativeApprox(dedx, bg);
1548#endif
1549 if (charge2 != 1) {
1550 dedx *= charge2;
1551#ifdef _BB_NONCONST_CORR_
1552 dedxDer *= charge2;
1553#endif
1554 }
1555 dE = dedx * xrho;
1556 }
1557 }
1558
1559 if (p < kMinP) {
1560 return false;
1561 }
1562 dETot = e - e0;
1563 } // end of e.loss correction
1564
1565 // Calculating the multiple scattering corrections******************
1566 value_t& fC22 = mC[kSigSnp2];
1567 value_t& fC33 = mC[kSigTgl2];
1568 value_t& fC43 = mC[kSigQ2PtTgl];
1569 value_t& fC44 = mC[kSigQ2Pt2];
1570 //
1571 value_t cC22(0.f), cC33(0.f), cC43(0.f), cC44(0.f);
1572 if (x2x0 != 0.f) {
1573 value_t beta2 = p02 / e2, theta2 = kMSConst2 / (beta2 * p02) * gpu::CAMath::Abs(x2x0);
1574 value_t fp34 = linRef.getTgl();
1575 if (charge2 != 1) {
1576 theta2 *= charge2;
1577 fp34 *= linRef.getCharge2Pt();
1578 }
1579 if (theta2 > constants::math::PI * constants::math::PI) {
1580 return false;
1581 }
1582 value_t t2c2I = theta2 * cst2I;
1583 cC22 = t2c2I * csp2;
1584 cC33 = t2c2I * cst2I;
1585 cC43 = t2c2I * fp34;
1586 cC44 = theta2 * fp34 * fp34;
1587 // optimize this
1588 // cC22 = theta2*((1.-getSnp())*(1.+getSnp()))*(1. + this->getTgl()*getTgl());
1589 // cC33 = theta2*(1. + this->getTgl()*getTgl())*(1. + this->getTgl()*getTgl());
1590 // cC43 = theta2*getTgl()*this->getQ2Pt()*(1. + this->getTgl()*getTgl());
1591 // cC44 = theta2*getTgl()*this->getQ2Pt()*getTgl()*this->getQ2Pt();
1592 }
1593
1594 // the energy loss correction contribution to cov.matrix: approximate energy loss fluctuation (M.Ivanov)
1595 constexpr value_t knst = 0.0007f; // To be tuned.
1596 value_t sigmadE = knst * gpu::CAMath::Sqrt(gpu::CAMath::Abs(dETot)) * e0 / p02 * linRef.getCharge2Pt();
1597 cC44 += sigmadE * sigmadE;
1598
1599 // Applying the corrections*****************************
1600 fC22 += cC22;
1601 fC33 += cC33;
1602 fC43 += cC43;
1603 fC44 += cC44;
1604 auto pscale = p0 / p;
1605 linRef.setQ2Pt(linRef.getQ2Pt() * pscale);
1606 this->setQ2Pt(this->getQ2Pt() * pscale);
1607
1608 checkCovariance();
1609
1610 return true;
1611}
1612
1613//______________________________________________________________
1614template <typename value_T>
1615GPUd() bool TrackParametrizationWithError<value_T>::getCovXYZPxPyPzGlo(std::array<value_t, kLabCovMatSize>& cv) const
1616{
1617 //---------------------------------------------------------------------
1618 // This function returns the global covariance matrix of the track params
1619 //
1620 // Cov(x,x) ... : cv[0]
1621 // Cov(y,x) ... : cv[1] cv[2]
1622 // Cov(z,x) ... : cv[3] cv[4] cv[5]
1623 // Cov(px,x)... : cv[6] cv[7] cv[8] cv[9]
1624 // Cov(py,x)... : cv[10] cv[11] cv[12] cv[13] cv[14]
1625 // Cov(pz,x)... : cv[15] cv[16] cv[17] cv[18] cv[19] cv[20]
1626 //
1627 // Results for (nearly) straight tracks are meaningless !
1628 //---------------------------------------------------------------------
1629 if (gpu::CAMath::Abs(this->getQ2Pt()) <= constants::math::Almost0 || gpu::CAMath::Abs(this->getSnp()) > constants::math::Almost1) {
1630 for (int i = 0; i < 21; i++) {
1631 cv[i] = 0.;
1632 }
1633 return false;
1634 }
1635
1636 const value_t pt = this->getPt();
1637 const value_t q2pt = this->getQ2Pt();
1638 value_t sn = 0.f, cs = 0.f;
1639 o2::math_utils::detail::sincos(this->getAlpha(), sn, cs);
1640 const value_t snp = this->getSnp();
1641 const value_t csp = gpu::CAMath::Sqrt((1.f - snp) * (1.f + snp));
1642 const value_t pXLoc = pt * csp;
1643 const value_t pYLoc = pt * snp;
1644 const value_t pZ = pt * this->getTgl();
1645 const value_t pX = cs * pXLoc - sn * pYLoc;
1646 const value_t pY = sn * pXLoc + cs * pYLoc;
1647
1648 value_t cTr[5][5] = {};
1649 for (int i = 0; i < kNParams; ++i) {
1650 for (int j = 0; j <= i; ++j) {
1651 cTr[i][j] = cTr[j][i] = mC[CovarMap[i][j]];
1652 }
1653 }
1654
1655 double jac[6][5] = {};
1656 jac[0][kY] = -sn;
1657 jac[1][kY] = cs;
1658 jac[2][kZ] = 1.f;
1659
1660 const value_t dPxDSnp = -pt * (cs * snp / csp + sn);
1661 const value_t dPyDSnp = pt * (cs - sn * snp / csp);
1662 jac[3][kSnp] = dPxDSnp;
1663 jac[4][kSnp] = dPyDSnp;
1664 jac[5][kTgl] = pt;
1665
1666 jac[3][kQ2Pt] = -pX / q2pt;
1667 jac[4][kQ2Pt] = -pY / q2pt;
1668 jac[5][kQ2Pt] = -pZ / q2pt;
1669
1670 int idx = 0;
1671 for (int i = 0; i < 6; ++i) {
1672 for (int j = 0; j <= i; ++j) {
1673 double cij = 0.f;
1674 for (int k = 0; k < kNParams; ++k) {
1675 for (int l = 0; l < kNParams; ++l) {
1676 cij += jac[i][k] * cTr[k][l] * jac[j][l];
1677 }
1678 }
1679 cv[idx++] = cij;
1680 }
1681 }
1682
1683 return true;
1684}
1685
1686#ifndef GPUCA_ALIGPUCODE
1687//______________________________________________________________
1688template <typename value_T>
1689std::string TrackParametrizationWithError<value_T>::asString() const
1690{
1691 return TrackParametrization<value_t>::asString() +
1692 fmt::format(" Cov: [{:+.3e}] [{:+.3e} {:+.3e}] [{:+.3e} {:+.3e} {:+.3e}] [{:+.3e} {:+.3e} {:+.3e} {:+.3e}] [{:+.3e} {:+.3e} {:+.3e} {:+.3e} {:+.3e}]",
1693 mC[kSigY2], mC[kSigZY], mC[kSigZ2], mC[kSigSnpY], mC[kSigSnpZ], mC[kSigSnp2], mC[kSigTglY],
1694 mC[kSigTglZ], mC[kSigTglSnp], mC[kSigTgl2], mC[kSigQ2PtY], mC[kSigQ2PtZ], mC[kSigQ2PtSnp], mC[kSigQ2PtTgl],
1695 mC[kSigQ2Pt2]);
1696}
1697
1698template <typename value_T>
1699std::string TrackParametrizationWithError<value_T>::asStringHexadecimal()
1700{
1701 return TrackParametrization<value_t>::asStringHexadecimal() +
1702 fmt::format(" Cov: [{:x}] [{:x} {:x}] [{:x} {:x} {:x}] [{:x} {:x} {:x} {:x}] [{:x} {:x} {:x} {:x} {:x}]",
1703 reinterpret_cast<const unsigned int&>(mC[kSigY2]), reinterpret_cast<const unsigned int&>(mC[kSigZY]), reinterpret_cast<const unsigned int&>(mC[kSigZ2]),
1704 reinterpret_cast<const unsigned int&>(mC[kSigSnpY]), reinterpret_cast<const unsigned int&>(mC[kSigSnpZ]), reinterpret_cast<const unsigned int&>(mC[kSigSnp2]),
1705 reinterpret_cast<const unsigned int&>(mC[kSigTglY]), reinterpret_cast<const unsigned int&>(mC[kSigTglZ]), reinterpret_cast<const unsigned int&>(mC[kSigTglSnp]),
1706 reinterpret_cast<const unsigned int&>(mC[kSigTgl2]), reinterpret_cast<const unsigned int&>(mC[kSigQ2PtY]), reinterpret_cast<const unsigned int&>(mC[kSigQ2PtZ]),
1707 reinterpret_cast<const unsigned int&>(mC[kSigQ2PtSnp]), reinterpret_cast<const unsigned int&>(mC[kSigQ2PtTgl]), reinterpret_cast<const unsigned int&>(mC[kSigQ2Pt2]));
1708}
1709#endif
1710
1711//______________________________________________________________
1712template <typename value_T>
1713GPUd() void TrackParametrizationWithError<value_T>::print() const
1714{
1715 // print parameters
1716#ifndef GPUCA_ALIGPUCODE
1717 printf("%s\n", asString().c_str());
1718#elif !defined(GPUCA_GPUCODE_DEVICE) || (!defined(__OPENCL__) && defined(GPUCA_GPU_DEBUG_PRINT))
1719 TrackParametrization<value_T>::printParam();
1720 printf(
1721 " Cov: [%+.3e] [%+.3e %+.3e] [%+.3e %+.3e %+.3e] [%+.3e %+.3e %+.3e %+.3e] [%+.3e %+.3e %+.3e %+.3e %+.3e]\n",
1722 mC[kSigY2], mC[kSigZY], mC[kSigZ2], mC[kSigSnpY], mC[kSigSnpZ], mC[kSigSnp2], mC[kSigTglY],
1723 mC[kSigTglZ], mC[kSigTglSnp], mC[kSigTgl2], mC[kSigQ2PtY], mC[kSigQ2PtZ], mC[kSigQ2PtSnp], mC[kSigQ2PtTgl],
1724 mC[kSigQ2Pt2]);
1725#endif
1726}
1727
1728//______________________________________________________________
1729template <typename value_T>
1730GPUd() void TrackParametrizationWithError<value_T>::printHexadecimal()
1731{
1732 // print parameters
1733#ifndef GPUCA_ALIGPUCODE
1734 printf("%s\n", asStringHexadecimal().c_str());
1735#elif !defined(GPUCA_GPUCODE_DEVICE) || (!defined(__OPENCL__) && defined(GPUCA_GPU_DEBUG_PRINT))
1736 TrackParametrization<value_T>::printParamHexadecimal();
1737 printf(
1738 " Cov: [%x] [%x %x] [%x %x %x] [%x %x %x %x] [%x %x %x %x %x]\n",
1739 gpu::CAMath::Float2UIntReint(mC[kSigY2]),
1740 gpu::CAMath::Float2UIntReint(mC[kSigZY]), gpu::CAMath::Float2UIntReint(mC[kSigZ2]),
1741 gpu::CAMath::Float2UIntReint(mC[kSigSnpY]), gpu::CAMath::Float2UIntReint(mC[kSigSnpZ]), gpu::CAMath::Float2UIntReint(mC[kSigSnp2]),
1742 gpu::CAMath::Float2UIntReint(mC[kSigTglY]), gpu::CAMath::Float2UIntReint(mC[kSigTglZ]), gpu::CAMath::Float2UIntReint(mC[kSigTglSnp]), gpu::CAMath::Float2UIntReint(mC[kSigTgl2]),
1743 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]));
1744#endif
1745}
1746
1747#ifndef GPUCA_ALIGPUCODE
1748//______________________________________________________________
1749template <typename value_T>
1750bool TrackParametrizationWithError<value_T>::toFwdTrackParCov(TrackParCovFwd& t) const
1751{
1752 auto p = this->getXYZGlo();
1753 t.setZ(p.Z());
1754 t.setX(p.X());
1755 t.setY(p.Y());
1756 t.setPhi(this->getPhi());
1757 t.setTanl(this->getTgl());
1758 t.setInvQPt(this->getQ2Pt());
1759 //
1760 if (gpu::CAMath::Abs(this->getSnp()) >= o2::constants::math::Almost1 ||
1761 gpu::CAMath::Abs(this->getTgl()) <= o2::constants::math::Almost0) {
1762 return false;
1763 }
1764 value_T csa, sna, csP, snP, csp = gpu::CAMath::Sqrt((1. - this->getSnp()) * (1. + this->getSnp()));
1765 math_utils::detail::sincos(value_T(this->getAlpha()), sna, csa);
1766 math_utils::detail::sincos(value_T(t.getPhi()), snP, csP);
1767 /*
1768 Jacobian is
1769 /-sna -csP/tgL 0 0 0 \
1770 | csa -snP/tgL 0 0 0 |
1771 | 0 0 1/csp 0 0 |
1772 | 0 0 0 1 0 |
1773 \ 0 0 0 0 1 /
1774 */
1775 auto tgLI = 1 / this->getTgl();
1776 const value_T d1 = -sna;
1777 const value_T d2 = -csP * tgLI;
1778 const value_T e1 = csa;
1779 const value_T e2 = -snP * tgLI;
1780 const value_T f1 = 1 / csp;
1781 SMatrix55Sym C;
1782 C(0, 0) = d1 * d1 * getSigmaY2() + 2 * d1 * d2 * getSigmaZY() + d2 * d2 * getSigmaZ2();
1783 C(0, 1) = d1 * e1 * getSigmaY2() + (d1 * e2 + d2 * e1) * getSigmaZY() + d2 * e2 * getSigmaZ2();
1784 C(1, 1) = e1 * e1 * getSigmaY2() + 2 * e1 * e2 * getSigmaZY() + e2 * e2 * getSigmaZ2();
1785
1786 C(0, 2) = f1 * (d1 * getSigmaSnpY() + d2 * getSigmaSnpZ());
1787 C(1, 2) = f1 * (e1 * getSigmaSnpY() + e2 * getSigmaSnpZ());
1788 C(2, 2) = f1 * f1 * getSigmaSnp2();
1789
1790 C(0, 3) = d1 * getSigmaTglY() + d2 * getSigmaTglZ();
1791 C(1, 3) = e1 * getSigmaTglY() + e2 * getSigmaTglZ();
1792 C(2, 3) = f1 * getSigmaTglSnp();
1793 C(3, 3) = getSigmaTgl2();
1794
1795 C(0, 4) = d1 * getSigma1PtY() + d2 * getSigma1PtZ();
1796 C(1, 4) = e1 * getSigma1PtY() + e2 * getSigma1PtZ();
1797 C(2, 4) = f1 * getSigma1PtSnp();
1798 C(3, 4) = getSigma1PtTgl();
1799 C(4, 4) = getSigma1Pt2();
1800 t.setCovariances(C);
1801 return true;
1802}
1803#endif
1804
1805namespace o2::track
1806{
1807#if !defined(GPUCA_GPUCODE) || defined(GPUCA_GPUCODE_DEVICE) // FIXME: DR: WORKAROUND to avoid CUDA bug creating host symbols for device code.
1808template class TrackParametrizationWithError<float>;
1809#endif
1810#ifndef GPUCA_GPUCODE
1811template class TrackParametrizationWithError<double>;
1812#endif
1813} // 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::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
cc
std::vector< o2::mch::ChannelCode > cc
Definition testStatusMap.cxx:37
LOG
LOG(info)<< "Compressed in "<< sw.CpuTime()<< " s"