ATLAS Offline Software
Loading...
Searching...
No Matches
RungeKuttaUtils.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
3*/
4
6// RungeKuttaUtils.cxx , (c) ATLAS Detector software
7// author Igor Gavrilenko
9
10#include "TrkExUtils/RungeKuttaUtils.h"
11#include "TrkPatternParameters/PatternTrackParameters.h"
12#include "TrkSurfaces/ConeSurface.h"
13#include "TrkSurfaces/CylinderSurface.h"
14#include "TrkSurfaces/DiscSurface.h"
15#include "TrkSurfaces/DistanceSolution.h"
16#include "TrkSurfaces/PerigeeSurface.h"
17#include "TrkSurfaces/PlaneSurface.h"
18#include "TrkSurfaces/StraightLineSurface.h"
19
20#include "CxxUtils/vec.h"
21#include "CxxUtils/inline_hints.h"
22
23namespace {
24/*
25 * Hide all internal implementation methods
26 * inside an anonymous namespace
27 */
28
29/*The notation of this package
30for array P[42].
31 /dL0 /dL1 /dPhi /dThe /dCM
32X ->P[0] dX / P[ 7] P[14] P[21] P[28] P[35]
33Y ->P[1] dY / P[ 8] P[15] P[22] P[29] P[36]
34Z ->P[2] dZ / P[ 9] P[16] P[23] P[30] P[37]
35Ax ->P[3] dAx/ P[10] P[17] P[24] P[31] P[38]
36Ay ->P[4] dAy/ P[11] P[18] P[25] P[32] P[39]
37Az ->P[5] dAz/ P[12] P[19] P[26] P[33] P[40]
38CM ->P[6] dCM/ P[13] P[20] P[27] P[34] P[41]
39*/
40
41inline void
42globalToLocalVecHelper(double* ATH_RESTRICT P,
43 const double s0,
44 const double s1,
45 const double s2,
46 const double s3,
47 const double s4)
48{
49 using namespace CxxUtils;
50 using vec2 = CxxUtils::vec<double, 2>;
51 using vec4 = CxxUtils::vec<double, 4>;
52 /* The calculation (original form)
53 P[ 7]-=(s0*P[ 3]); P[ 8]-=(s0*P[ 4]); P[ 9]-=(s0*P[ 5]);
54 P[10]-=(s0*P[42]); P[11]-=(s0*P[43]); P[12]-=(s0*P[44]);
55 P[14]-=(s1*P[ 3]); P[15]-=(s1*P[ 4]); P[16]-=(s1*P[ 5]);
56 P[17]-=(s1*P[42]); P[18]-=(s1*P[43]); P[19]-=(s1*P[44]);
57 P[21]-=(s2*P[ 3]); P[22]-=(s2*P[ 4]); P[23]-=(s2*P[ 5]);
58 P[24]-=(s2*P[42]); P[25]-=(s2*P[43]); P[26]-=(s2*P[44]);
59 P[28]-=(s3*P[ 3]); P[29]-=(s3*P[ 4]); P[30]-=(s3*P[ 5]);
60 P[31]-=(s3*P[42]); P[32]-=(s3*P[43]); P[33]-=(s3*P[44]);
61 P[35]-=(s4*P[ 3]); P[36]-=(s4*P[ 4]); P[37]-=(s4*P[ 5]);
62 P[38]-=(s4*P[42]); P[39]-=(s4*P[43]); P[40]-=(s4*P[44]);
63 */
64 /*
65 * The naming convention we follow is
66 * A_B_C_D -->SIMD vector {A,B , C, D}
67 * For example :
68 * dXdTheta_dYdTheta_dZdTheta_dAxdTheta --> {dX/dTheta, dY/dTheta, dZ/dTheta,
69 * dAx/dTheta}
70 * --> {P[28,P[29],P[30],P[31]}
71 */
72 vec4 Pmult1_2 = {P[3], P[4], P[5], P[42]};
73 vec2 Pmult3 = {P[43], P[44]};
74
75 vec4 dXdL0_dYdL0_dZdL0_dAxdL0;
76 vload(dXdL0_dYdL0_dZdL0_dAxdL0, &P[7]);
77 vec2 dAydL0_dAzdL0;
78 vload(dAydL0_dAzdL0, &P[11]);
79 dXdL0_dYdL0_dZdL0_dAxdL0 -= s0 * Pmult1_2;
80 dAydL0_dAzdL0 -= s0 * Pmult3;
81 vstore(&P[7], dXdL0_dYdL0_dZdL0_dAxdL0);
82 vstore(&P[11], dAydL0_dAzdL0);
83
84 vec4 dXdL1_dYdL1_dZdL1_dAxdL1;
85 vload(dXdL1_dYdL1_dZdL1_dAxdL1, &P[14]);
86 vec2 dAydL1_dAzdL1;
87 vload(dAydL1_dAzdL1, &P[18]);
88 dXdL1_dYdL1_dZdL1_dAxdL1 -= s1 * Pmult1_2;
89 dAydL1_dAzdL1 -= s1 * Pmult3;
90 vstore(&P[14], dXdL1_dYdL1_dZdL1_dAxdL1);
91 vstore(&P[18], dAydL1_dAzdL1);
92
93 vec4 dXdPhi_dYdPhi_dZdPhi_dAxdPhi;
94 vload(dXdPhi_dYdPhi_dZdPhi_dAxdPhi, &P[21]);
95 vec2 dAydPhi_dAzdPhi;
96 vload(dAydPhi_dAzdPhi, &P[25]);
97 dXdPhi_dYdPhi_dZdPhi_dAxdPhi -= s2 * Pmult1_2;
98 dAydPhi_dAzdPhi -= s2 * Pmult3;
99 vstore(&P[21], dXdPhi_dYdPhi_dZdPhi_dAxdPhi);
100 vstore(&P[25], dAydPhi_dAzdPhi);
101
102 vec4 dXdTheta_dYdTheta_dZdTheta_dAxdTheta;
103 vload(dXdTheta_dYdTheta_dZdTheta_dAxdTheta, &P[28]);
104 vec2 dAydTheta_dAzdTheta;
105 vload(dAydTheta_dAzdTheta, &P[32]);
106 dXdTheta_dYdTheta_dZdTheta_dAxdTheta -= s3 * Pmult1_2;
107 dAydTheta_dAzdTheta -= s3 * Pmult3;
108 vstore(&P[28], dXdTheta_dYdTheta_dZdTheta_dAxdTheta);
109 vstore(&P[32], dAydTheta_dAzdTheta);
110
111 vec4 dXdCM_dYdCM_dZdCM_dAxdCM;
112 vload(dXdCM_dYdCM_dZdCM_dAxdCM, &P[35]);
113 vec2 AydCM_dAzdCM;
114 vload(AydCM_dAzdCM, &P[39]);
115 dXdCM_dYdCM_dZdCM_dAxdCM -= s4 * Pmult1_2;
116 AydCM_dAzdCM -= s4 * Pmult3;
117 vstore(&P[35], dXdCM_dYdCM_dZdCM_dAxdCM);
118 vstore(&P[39], AydCM_dAzdCM);
119}
120
121/* Helper to replace repeated calculation of
122 * 5x1 = 5x3 * 3X1
123 * for the Jacobian
124 *
125 * E.g a calculation like :
126 * Jac[ 0] = Ax[0]*P[ 7]+Ax[1]*P[ 8]+Ax[2]*P[ 9]; // dL0/dL0
127 * Jac[ 1] = Ax[0]*P[14]+Ax[1]*P[15]+Ax[2]*P[16]; // dL0/dL1
128 * Jac[ 2] = Ax[0]*P[21]+Ax[1]*P[22]+Ax[2]*P[23]; // dL0/dPhi
129 * Jac[ 3] = Ax[0]*P[28]+Ax[1]*P[29]+Ax[2]*P[30]; // dL0/dThe
130 * Jac[ 4] = Ax[0]*P[35]+Ax[1]*P[36]+Ax[2]*P[37]; // dL0/dCM
131 * Jac[ 5] = Ay[0]*P[ 7]+Ay[1]*P[ 8]+Ay[2]*P[ 9]; // dL1/dL0
132 * Jac[ 6] = Ay[0]*P[14]+Ay[1]*P[15]+Ay[2]*P[16]; // dL1/dL1
133 * Jac[ 7] = Ay[0]*P[21]+Ay[1]*P[22]+Ay[2]*P[23]; // dL1/dPhi
134 * Jac[ 8] = Ay[0]*P[28]+Ay[1]*P[29]+Ay[2]*P[30]; // dL1/dThe
135 * Jac[ 9] = Ay[0]*P[35]+Ay[1]*P[36]+Ay[2]*P[37]; // dL1/dCM
136 * is replaces with
137 * mult3x5Helper(&Jac[0],Ax,&P[7]);
138 * mult3x5Helper(&Jac[5],Ay,&P[7]);
139 */
140inline void
141mult3x5Helper(double* ATH_RESTRICT Jac,
142 const double* ATH_RESTRICT V,
143 const double* ATH_RESTRICT P)
144{
145 /*
146 * |Jac[0] |= |V[0]| * |P[0]| + |V[1]| * |P[1] | + |V[2]| * |P[2] |
147 * |Jac[1] |= |V[0]| * |P[7]| + |V[1]| * |P[8] | + |V[2]| * |P[9] |
148 * |Jac[2] |= |V[0]| * |P[14]| + |V[1]| * |P[15]| + |V[2]| * |P[16]|
149 * |Jac[3] |= |V[0]| * |P[21]| + |V[1]| * |P[22]| + |V[2]| * |P[23]|
150 *
151 * Jac[4] = V[0] * P[28] + V[1] * P[29] + V[2] * P[30];
152 *
153 * The first 4 we can do in vertical SIMD fashion
154 * {Jac[0] | Jac[1] Jac[2] | Jac[3] } =
155 * V[0] * {P[0] | P[7] | P[14] | P[21]} +
156 * V[1] * {P[1] | P[8] | P[15] | P[22]} +
157 * V[2] * {P[2] | P[9] | P[16] | P[23]}
158 * Where {} is a SIMD size 4 vector
159 *
160 * The remaining odd (5th) element is done at the end
161 * Jac[4] = V[0] * P[28] + V[1] * P[29] + V[2] * P[30];
162 */
163
164 using vec4 = CxxUtils::vec<double, 4>;
165 // 1st 4 elements
166 vec4 P0 = { P[0], P[7], P[14], P[21] };
167 vec4 res = V[0] * P0;
168
169 vec4 P1 = { P[1], P[8], P[15], P[22] };
170 res += V[1] * P1;
171
172 vec4 P2 = { P[2], P[9], P[16], P[23] };
173 res += V[2] * P2;
174
175 CxxUtils::vstore(&Jac[0], res);
176
177 // The 5th element
178 Jac[4] = V[0] * P[28] + V[1] * P[29] + V[2] * P[30];
179}
180
181void
182transformGlobalToPlane(const Amg::Transform3D& T,
183 bool useJac,
184 double* ATH_RESTRICT P,
185 double* ATH_RESTRICT par,
186 double* ATH_RESTRICT Jac)
187{
188 const double* Ax = T.matrix().col(0).data(); // Oth column
189 const double* Ay = T.matrix().col(1).data(); // 1st column
190
191 const double d[3] = { P[0] - T(0, 3), P[1] - T(1, 3), P[2] - T(2, 3) };
192
193 par[0] = d[0] * Ax[0] + d[1] * Ax[1] + d[2] * Ax[2];
194 par[1] = d[0] * Ay[0] + d[1] * Ay[1] + d[2] * Ay[2];
195
196 if (!useJac)
197 return;
198
199 // Condition trajectory on surface
200 //
201 double S[3] = { T(0, 2), T(1, 2), T(2, 2) };
202
203 double A = P[3] * S[0] + P[4] * S[1] + P[5] * S[2];
204 if (A != 0.)
205 A = 1. / A;
206 S[0] *= A;
207 S[1] *= A;
208 S[2] *= A;
209
210 double s[5] = {};
211 mult3x5Helper(s, S, &P[7]);
212 globalToLocalVecHelper(P, s[0], s[1], s[2], s[3], s[4]);
213 // Jacobian production
214 mult3x5Helper(&Jac[0], Ax, &P[7]);
215 mult3x5Helper(&Jac[5], Ay, &P[7]);
216}
217
219// Global position transformation to local Disc system coordinate
221
222void
223transformGlobalToDisc(const Amg::Transform3D& T,
224 bool useJac,
225 double* ATH_RESTRICT P,
226 double* ATH_RESTRICT par,
227 double* ATH_RESTRICT Jac)
228{
229 const double* Ax = T.matrix().col(0).data(); // Oth column
230 const double* Ay = T.matrix().col(1).data(); // 1st column
231
232 const double d[3] = { P[0] - T(0, 3), P[1] - T(1, 3), P[2] - T(2, 3) };
233
234 const double RC = d[0] * Ax[0] + d[1] * Ax[1] + d[2] * Ax[2];
235 const double RS = d[0] * Ay[0] + d[1] * Ay[1] + d[2] * Ay[2];
236 const double R2 = RC * RC + RS * RS;
237 par[0] = std::sqrt(R2);
238 par[1] = atan2(RS, RC);
239
240 if (!useJac)
241 return;
242
243 // Condition trajectory on surface
244 //
245 double S[3] = { T(0, 2), T(1, 2), T(2, 2) };
246
247 double A = P[3] * S[0] + P[4] * S[1] + P[5] * S[2];
248 if (A != 0.)
249 A = 1. / A;
250 S[0] *= A;
251 S[1] *= A;
252 S[2] *= A;
253
254 double s[5] = {};
255 mult3x5Helper(s, S, &P[7]);
256 globalToLocalVecHelper(P, s[0], s[1], s[2], s[3], s[4]);
257 // Jacobian production
258 //
259 double Ri = 1. / par[0];
260
261 const double Av[3] = { (RC * Ax[0] + RS * Ay[0]) * Ri,
262 (RC * Ax[1] + RS * Ay[1]) * Ri,
263 (RC * Ax[2] + RS * Ay[2]) * Ri };
264 const double Bv[3] = { (RC * Ay[0] - RS * Ax[0]) * (Ri = 1. / R2),
265 (RC * Ay[1] - RS * Ax[1]) * Ri,
266 (RC * Ay[2] - RS * Ax[2]) * Ri };
267
268 mult3x5Helper(&Jac[0], Av, &P[7]);
269 mult3x5Helper(&Jac[5], Bv, &P[7]);
270}
272// Global position transformation to local Cylinder system coordinate
274void
275transformGlobalToCylinder(const Amg::Transform3D& T,
276 double R,
277 bool useJac,
278 double* ATH_RESTRICT P,
279 double* ATH_RESTRICT par,
280 double* ATH_RESTRICT Jac)
281{
282
283 const double* Ax = T.matrix().col(0).data(); // Oth column
284 const double* Ay = T.matrix().col(1).data(); // 1st column
285 const double* Az = T.matrix().col(2).data(); // 2nd column
286
287 double x = P[0] - T(0, 3);
288 double y = P[1] - T(1, 3);
289 double z = P[2] - T(2, 3);
290 double RC = x * Ax[0] + y * Ax[1] + z * Ax[2];
291 double RS = x * Ay[0] + y * Ay[1] + z * Ay[2];
292 par[0] = atan2(RS, RC) * R;
293 par[1] = x * Az[0] + y * Az[1] + z * Az[2];
294
295 if (!useJac)
296 return;
297
298 // Condition trajectory on surface
299 //
300 double B = par[1];
301 double C = P[3] * Az[0] + P[4] * Az[1] + P[5] * Az[2];
302 const double ax = P[3] - Az[0] * C;
303 x -= (B * Az[0]);
304 const double ay = P[4] - Az[1] * C;
305 y -= (B * Az[1]);
306 const double az = P[5] - Az[2] * C;
307 z -= (B * Az[2]);
308 double A = (ax * x + ay * y + az * z);
309 if (A != 0.)
310 A = 1. / A;
311 x *= A;
312 y *= A;
313 z *= A;
314
315 const double S[3] = { x, y, z };
316 double s[5] = {};
317 mult3x5Helper(s, S, &P[7]);
318 globalToLocalVecHelper(P, s[0], s[1], s[2], s[3], s[4]);
319 // Jacobian production
320 //
321 const double Av[3] = { (RC * Ay[0] - RS * Ax[0]) * (R = 1. / R),
322 (RC * Ay[1] - RS * Ax[1]) * R,
323 (RC * Ay[2] - RS * Ax[2]) * R };
324
325 mult3x5Helper(&Jac[0], Av, &P[7]);
326 mult3x5Helper(&Jac[5], Az, &P[7]);
327}
328
330// Global position transformation to local Straight line system coordinate
332
333void
334transformGlobalToLine(const Amg::Transform3D& T,
335 bool useJac,
336 double* ATH_RESTRICT P,
337 double* ATH_RESTRICT par,
338 double* ATH_RESTRICT Jac)
339{
340 const double* Az = T.matrix().col(2).data(); // 2nd column
341 double Bx = Az[1] * P[5] - Az[2] * P[4];
342 double By = Az[2] * P[3] - Az[0] * P[5];
343 double Bz = Az[0] * P[4] - Az[1] * P[3];
344 const double B2 = Bx * Bx + By * By + Bz * Bz;
345 const double x = P[0] - T(0, 3);
346 const double y = P[1] - T(1, 3);
347 const double z = P[2] - T(2, 3);
348 if (B2 > 1e-14) {
349 const double Bn = 1. / std::sqrt(B2);
350 Bx *= Bn;
351 By *= Bn;
352 Bz *= Bn;
353 par[0] = x * Bx + y * By + z * Bz;
354 } else {
355 // surface and trajectory are parallel, calculate the distance
356 // of closest approach differently.
357 par[0] = std::sqrt(
358 ( std::pow(Az[1] * z - Az[2] * y, 2.)
359 + std::pow(Az[2] * x - Az[0] * z, 2.)
360 + std::pow(Az[0] * y - Az[1] * x, 2.))
361 / (Az[0] * Az[0] + Az[1] * Az[1] + Az[2] * Az[2]));
362 }
363 par[1] = x * Az[0] + y * Az[1] + z * Az[2];
364
365 if (!useJac) {
366 return;
367 }
368
369 // Condition trajectory on surface
370 //
371 const double d = P[3] * Az[0] + P[4] * Az[1] + P[5] * Az[2];
372 double a = (1. - d) * (1. + d);
373 if (a != 0.) {
374 a = 1. / a;
375 }
376 const double X = d * Az[0] - P[3];
377 const double Y = d * Az[1] - P[4];
378 const double Z = d * Az[2] - P[5];
379
380 double D[5] = {};
381 mult3x5Helper(D, Az, &P[10]);
382 const double s0 =
383 (((P[7] * X + P[8] * Y + P[9] * Z) + x * (D[0] * Az[0] - P[10])) +
384 (y * (D[0] * Az[1] - P[11]) + z * (D[0] * Az[2] - P[12]))) *
385 a;
386 const double s1 =
387 (((P[14] * X + P[15] * Y + P[16] * Z) + x * (D[1] * Az[0] - P[17])) +
388 (y * (D[1] * Az[1] - P[18]) + z * (D[1] * Az[2] - P[19]))) *
389 a;
390 const double s2 =
391 (((P[21] * X + P[22] * Y + P[23] * Z) + x * (D[2] * Az[0] - P[24])) +
392 (y * (D[2] * Az[1] - P[25]) + z * (D[2] * Az[2] - P[26]))) *
393 a;
394 const double s3 =
395 (((P[28] * X + P[29] * Y + P[30] * Z) + x * (D[3] * Az[0] - P[31])) +
396 (y * (D[3] * Az[1] - P[32]) + z * (D[3] * Az[2] - P[33]))) *
397 a;
398 const double s4 =
399 (((P[35] * X + P[36] * Y + P[37] * Z) + x * (D[4] * Az[0] - P[38])) +
400 (y * (D[4] * Az[1] - P[39]) + z * (D[4] * Az[2] - P[40]))) *
401 a;
402
403 // pass -1 (As we want do add rather subtract in the helper)
404 globalToLocalVecHelper(P, -1. * s0, -1. * s1, -1. * s2, -1. * s3, -1. * s4);
405 // Jacobian production
406 const double B[3] = { Bx, By, Bz };
407 mult3x5Helper(&Jac[0], B, &P[7]);
408 mult3x5Helper(&Jac[5], Az, &P[7]);
409}
410
412// Global position transformation to local Cone system coordinate
414
415void
416transformGlobalToCone(const Amg::Transform3D& T,
417 double tA,
418 bool useJac,
419 const double* ATH_RESTRICT P,
420 double* ATH_RESTRICT par,
421 double* ATH_RESTRICT Jac)
422{
423 const double* Ax = T.matrix().col(0).data(); // Oth column
424 const double* Ay = T.matrix().col(1).data(); // 1st column
425 const double* Az = T.matrix().col(2).data(); // 2nd column
426
427 const double x = P[0] - T(0, 3);
428 const double y = P[1] - T(1, 3);
429 const double z = P[2] - T(2, 3);
430 const double RC = x * Ax[0] + y * Ax[1] + z * Ax[2];
431 const double RS = x * Ay[0] + y * Ay[1] + z * Ay[2];
432 par[1] = x * Az[0] + y * Az[1] + z * Az[2];
433 par[0] = atan2(RS, RC) * (par[1] * tA);
434
435 if (!useJac)
436 return;
437
438 Jac[0] = 0.; // dL0/dL0
439 Jac[1] = 0.; // dL0/dL1
440 Jac[2] = 0.; // dL0/dPhi
441 Jac[3] = 0.; // dL0/dThe
442 Jac[4] = 0.; // dL0/dCM
443 Jac[5] = 0.; // dL1/dL0
444 Jac[6] = 0.; // dL1/dL1
445 Jac[7] = 0.; // dL1/dPhi
446 Jac[8] = 0.; // dL1/dThe
447 Jac[9] = 0.; // dL1/dCM
448}
449
451// Plane local position transformation to global system coordinate
453
454void
455transformPlaneToGlobal(bool useJac,
456 const Amg::Transform3D& T,
457 const AmgVector(5) & ATH_RESTRICT p,
458 double* ATH_RESTRICT P)
459{
460 const double* Ax = T.matrix().col(0).data(); // Oth column
461 const double* Ay = T.matrix().col(1).data(); // 1st column
462
463 P[0] = p[0] * Ax[0] + p[1] * Ay[0] + T(0, 3); // X
464 P[1] = p[0] * Ax[1] + p[1] * Ay[1] + T(1, 3); // Y
465 P[2] = p[0] * Ax[2] + p[1] * Ay[2] + T(2, 3); // Z
466
467 if (!useJac)
468 return;
469
470 // /dL1 | /dL2 | /dPhi | /dThe |
471 P[7] = Ax[0];
472 P[14] = Ay[0];
473 P[21] = 0.;
474 P[28] = 0.; // dX/
475 P[8] = Ax[1];
476 P[15] = Ay[1];
477 P[22] = 0.;
478 P[29] = 0.; // dY/
479 P[9] = Ax[2];
480 P[16] = Ay[2];
481 P[23] = 0.;
482 P[30] = 0.; // dZ/
483}
484
486// Disc local position transformation to global system coordinate
488
489void
490transformDiscToGlobal(bool useJac,
491 const Amg::Transform3D& T,
492 const AmgVector(5) & ATH_RESTRICT p,
493 double* ATH_RESTRICT P)
494{
495 const double* Ax = T.matrix().col(0).data(); // Oth column
496 const double* Ay = T.matrix().col(1).data(); // 1st column
497
498 double Sf, Cf;
499 sincos(p[1], &Sf, &Cf);
500
501 const double d0 = Cf * Ax[0] + Sf * Ay[0];
502 const double d1 = Cf * Ax[1] + Sf * Ay[1];
503 const double d2 = Cf * Ax[2] + Sf * Ay[2];
504 P[0] = p[0] * d0 + T(0, 3); // X
505 P[1] = p[0] * d1 + T(1, 3); // Y
506 P[2] = p[0] * d2 + T(2, 3); // Z
507
508 if (!useJac)
509 return;
510
511 // /dL1 | /dL2 | /dPhi | /dThe |
512 P[7] = d0;
513 P[14] = p[0] * (Cf * Ay[0] - Sf * Ax[0]);
514 P[21] = 0.;
515 P[28] = 0.; // dX/
516 P[8] = d1;
517 P[15] = p[0] * (Cf * Ay[1] - Sf * Ax[1]);
518 P[22] = 0.;
519 P[29] = 0.; // dY/
520 P[9] = d2;
521 P[16] = p[0] * (Cf * Ay[2] - Sf * Ax[2]);
522 P[23] = 0.;
523 P[30] = 0.; // dZ/
524}
525
527// Cylinder local position transformation to global system coordinate
529
530void
531transformCylinderToGlobal(bool useJac,
532 const Amg::Transform3D& T,
533 double R,
534 const AmgVector(5) & ATH_RESTRICT p,
535 double* ATH_RESTRICT P)
536{
537 const double* Ax = T.matrix().col(0).data(); // Oth column
538 const double* Ay = T.matrix().col(1).data(); // 1st column
539 const double* Az = T.matrix().col(2).data(); // 2nd column
540
541 const double fr = p[0] / R;
542 double Sf, Cf;
543 sincos(fr, &Sf, &Cf);
544
545 P[0] = R * (Cf * Ax[0] + Sf * Ay[0]) + p[1] * Az[0] + T(0, 3); // X
546 P[1] = R * (Cf * Ax[1] + Sf * Ay[1]) + p[1] * Az[1] + T(1, 3); // Y
547 P[2] = R * (Cf * Ax[2] + Sf * Ay[2]) + p[1] * Az[2] + T(2, 3); // Z
548
549 if (!useJac)
550 return;
551
552 // /dL1 | /dL2 | /dPhi | /dThe |
553 P[7] = Cf * Ay[0] - Sf * Ax[0];
554 P[14] = Az[0];
555 P[21] = 0.;
556 P[28] = 0.; // dX/
557 P[8] = Cf * Ay[1] - Sf * Ax[1];
558 P[15] = Az[1];
559 P[22] = 0.;
560 P[29] = 0.; // dY/
561 P[9] = Cf * Ay[2] - Sf * Ax[2];
562 P[16] = Az[2];
563 P[23] = 0.;
564 P[30] = 0.; // dZ/
565}
566
568// Straight line local position transformation to global system coordinate
570
571void
572transformLineToGlobal(bool useJac,
573 const Amg::Transform3D& T,
574 const AmgVector(5) & ATH_RESTRICT p,
575 double* ATH_RESTRICT P)
576{
577 const double* Az = T.matrix().col(2).data(); // 2nd column
578
579 double Bx = Az[1] * P[5] - Az[2] * P[4];
580 double By = Az[2] * P[3] - Az[0] * P[5];
581 double Bz = Az[0] * P[4] - Az[1] * P[3];
582 const double tmp_B2 = Bx * Bx + By * By + Bz * Bz;
583 constexpr double epsilon2 = 1e-14;
584 // assume that B := Bx, By, Boa is null vector if its norm is very small wrt. P:=P[3],P[4],P[5]:
585 // || B || < epsilon * || P|| => || B || = 0
586 // P seems to be always normalised, therefor it is assumed its norm is 1.
587 const double Bn = ( tmp_B2 > epsilon2 ? 1. / std::sqrt(tmp_B2) : 0. );
588 Bx *= Bn;
589 By *= Bn;
590 Bz *= Bn;
591 P[0] = p[1] * Az[0] + Bx * p[0] + T(0, 3); // X
592 P[1] = p[1] * Az[1] + By * p[0] + T(1, 3); // Y
593 P[2] = p[1] * Az[2] + Bz * p[0] + T(2, 3); // Z
594
595 if (!useJac)
596 return;
597
598 double Bx2 = -Az[2] * P[25];
599 double Bx3 = Az[1] * P[33] - Az[2] * P[32];
600 double By2 = Az[2] * P[24];
601 double By3 = Az[2] * P[31] - Az[0] * P[33];
602 double Bz2 = Az[0] * P[25] - Az[1] * P[24];
603 double Bz3 = Az[0] * P[32] - Az[1] * P[31];
604 const double B2 = Bx * Bx2 + By * By2 + Bz * Bz2;
605 const double B3 = Bx * Bx3 + By * By3 + Bz * Bz3;
606 Bx2 = (Bx2 - Bx * B2) * Bn;
607 Bx3 = (Bx3 - Bx * B3) * Bn;
608 By2 = (By2 - By * B2) * Bn;
609 By3 = (By3 - By * B3) * Bn;
610 Bz2 = (Bz2 - Bz * B2) * Bn;
611 Bz3 = (Bz3 - Bz * B3) * Bn;
612
613 // /dL1 | /dL2 | /dPhi | /dThe |
614 P[7] = Bx;
615 P[14] = Az[0];
616 P[21] = Bx2 * p[0];
617 P[28] = Bx3 * p[0]; // dX/
618 P[8] = By;
619 P[15] = Az[1];
620 P[22] = By2 * p[0];
621 P[29] = By3 * p[0]; // dY/
622 P[9] = Bz;
623 P[16] = Az[2];
624 P[23] = Bz2 * p[0];
625 P[30] = Bz3 * p[0]; // dZ/
626}
627
629// Jacobian of transformations from curvilinear to Plane system coordinates
631
632void
633jacobianTransformCurvilinearToPlane(double* ATH_RESTRICT P,
634 double* ATH_RESTRICT Jac)
635{
636 const double* At = &P[4];
637 double* Au = &P[7];
638 double* Av = &P[10];
639 const double* Ax = &P[13];
640 const double* Ay = &P[16];
641 double* S = &P[19];
642
643 double A = At[0] * S[0] + At[1] * S[1] + At[2] * S[2];
644 if (A != 0.)
645 A = 1. / A;
646 S[0] *= A;
647 S[1] *= A;
648 S[2] *= A;
649
650 const double s1 = Au[0] * S[0] + Au[1] * S[1];
651 const double s2 = Av[0] * S[0] + Av[1] * S[1] + Av[2] * S[2];
652
653 Au[0] -= (s1 * At[0]);
654 Au[1] -= (s1 * At[1]);
655 Au[2] -= (s1 * At[2]);
656 Av[0] -= (s2 * At[0]);
657 Av[1] -= (s2 * At[1]);
658 Av[2] -= (s2 * At[2]);
659
660 Jac[0] = Ax[0] * Au[0] + Ax[1] * Au[1] + Ax[2] * Au[2]; // dL0/dL0
661 Jac[1] = Ax[0] * Av[0] + Ax[1] * Av[1] + Ax[2] * Av[2]; // dL0/dL1
662 Jac[2] = 0.; // dL0/dPhi
663 Jac[3] = 0.; // dL0/dThe
664 Jac[5] = Ay[0] * Au[0] + Ay[1] * Au[1] + Ay[2] * Au[2]; // dL1/dL0
665 Jac[6] = Ay[0] * Av[0] + Ay[1] * Av[1] + Ay[2] * Av[2]; // dL1/dL1
666 Jac[7] = 0.; // dL1/dPhi
667 Jac[8] = 0.; // dL1/dThe
668}
669
671// Jacobian of transformations from curvilinear to Disc system coordinates
673
674void
675jacobianTransformCurvilinearToDisc(double* ATH_RESTRICT P,
676 double* ATH_RESTRICT Jac)
677{
678 const double* p = &P[0];
679 const double* At = &P[4];
680 double* Au = &P[7];
681 double* Av = &P[10];
682 const double* Ax = &P[13];
683 const double* Ay = &P[16];
684 double* S = &P[19];
685
686 // Condition trajectory on surface
687 //
688 double A = At[0] * S[0] + At[1] * S[1] + At[2] * S[2];
689 if (A != 0.)
690 A = 1. / A;
691 S[0] *= A;
692 S[1] *= A;
693 S[2] *= A;
694
695 const double s1 = Au[0] * S[0] + Au[1] * S[1];
696 const double s2 = Av[0] * S[0] + Av[1] * S[1] + Av[2] * S[2];
697
698 Au[0] -= (s1 * At[0]);
699 Au[1] -= (s1 * At[1]);
700 Au[2] -= (s1 * At[2]);
701 Av[0] -= (s2 * At[0]);
702 Av[1] -= (s2 * At[1]);
703 Av[2] -= (s2 * At[2]);
704
705 // Jacobian production
706 //
707 double Sf, Cf;
708 sincos(p[1], &Sf, &Cf);
709
710 const double Ri = 1. / p[0];
711 const double A0 = (Cf * Ax[0] + Sf * Ay[0]);
712 const double A1 = (Cf * Ax[1] + Sf * Ay[1]);
713 const double A2 = (Cf * Ax[2] + Sf * Ay[2]);
714 const double B0 = (Cf * Ay[0] - Sf * Ax[0]) * Ri;
715 const double B1 = (Cf * Ay[1] - Sf * Ax[1]) * Ri;
716 const double B2 = (Cf * Ay[2] - Sf * Ax[2]) * Ri;
717
718 Jac[0] = A0 * Au[0] + A1 * Au[1] + A2 * Au[2]; // dL0/dL0
719 Jac[1] = A0 * Av[0] + A1 * Av[1] + A2 * Av[2]; // dL0/dL1
720 Jac[2] = 0.; // dL0/dPhi
721 Jac[3] = 0.; // dL0/dThe
722 Jac[5] = B0 * Au[0] + B1 * Au[1] + B2 * Au[2]; // dL1/dL0
723 Jac[6] = B0 * Av[0] + B1 * Av[1] + B2 * Av[2]; // dL1/dL1
724 Jac[7] = 0.; // dL1/dPhi
725 Jac[8] = 0.; // dL1/dThe
726}
727
729// Jacobian of transformations from curvilinear to Cylinder system coordinates
731
732void
733jacobianTransformCurvilinearToCylinder(double* ATH_RESTRICT P,
734 double* ATH_RESTRICT Jac)
735{
736 const double* p = &P[0];
737 const double* At = &P[4];
738 double* Au = &P[7];
739 double* Av = &P[10];
740 const double* Ax = &P[13];
741 const double* Ay = &P[16];
742 const double* Az = &P[19];
743 const double R = P[22];
744
745 const double fr = p[0] / R;
746 double Sf, Cf;
747 sincos(fr, &Sf, &Cf);
748
749 const double x = Cf * Ax[0] + Sf * Ay[0];
750 const double y = Cf * Ax[1] + Sf * Ay[1];
751 const double z = Cf * Ax[2] + Sf * Ay[2];
752
753 // Condition trajectory on surface
754 //
755 const double C = At[0] * Az[0] + At[1] * Az[1] + At[2] * Az[2];
756 const double ax = At[0] - Az[0] * C;
757 const double ay = At[1] - Az[1] * C;
758 const double az = At[2] - Az[2] * C;
759 double A = (ax * x + ay * y + az * z);
760 if (A != 0.)
761 A = 1. / A;
762
763 const double s1 = (Au[0] * x + Au[1] * y) * A;
764 const double s2 = (Av[0] * x + Av[1] * y + Av[2] * z) * A;
765
766 Au[0] -= (s1 * At[0]);
767 Au[1] -= (s1 * At[1]);
768 Au[2] -= (s1 * At[2]);
769 Av[0] -= (s2 * At[0]);
770 Av[1] -= (s2 * At[1]);
771 Av[2] -= (s2 * At[2]);
772
773 // Jacobian production
774 //
775 const double A0 = (Cf * Ay[0] - Sf * Ax[0]);
776 const double A1 = (Cf * Ay[1] - Sf * Ax[1]);
777 const double A2 = (Cf * Ay[2] - Sf * Ax[2]);
778
779 Jac[0] = A0 * Au[0] + A1 * Au[1] + A2 * Au[2]; // dL0/dL0
780 Jac[1] = A0 * Av[0] + A1 * Av[1] + A2 * Av[2]; // dL0/dL1
781 Jac[2] = 0.; // dL0/dPhi
782 Jac[3] = 0.; // dL0/dThe
783 Jac[5] = Az[0] * Au[0] + Az[1] * Au[1] + Az[2] * Au[2]; // dL1/dL0
784 Jac[6] = Az[0] * Av[0] + Az[1] * Av[1] + Az[2] * Av[2]; // dL1/dL1
785 Jac[7] = 0.; // dL1/dPhi
786 Jac[8] = 0.; // dL1/dThe
787}
788
790// Jacobian of transformations from curvilinear to Straight Line system
791// coordinates
793
794void
795jacobianTransformCurvilinearToStraightLine(const double* ATH_RESTRICT P,
796 double* ATH_RESTRICT Jac)
797{
798 const double* p = &P[0];
799 const double* At = &P[4];
800 const double* Au = &P[7];
801 const double* Av = &P[10];
802 const double* A = &P[19];
803
804 double Bx = A[1] * At[2] - A[2] * At[1];
805 double By = A[2] * At[0] - A[0] * At[2];
806 double Bz = A[0] * At[1] - A[1] * At[0];
807 const double Bn = 1. / std::sqrt(Bx * Bx + By * By + Bz * Bz);
808 Bx *= Bn;
809 By *= Bn;
810 Bz *= Bn;
811
812 // Condition trajectory on surface
813 //
814 const double d = At[0] * A[0] + At[1] * A[1] + At[2] * A[2];
815 double a = (1. - d) * (1. + d);
816 if (a != 0.)
817 a = 1. / a;
818 const double X = d * A[0] - At[0], Y = d * A[1] - At[1], Z = d * A[2] - At[2];
819
820 const double s1 = (Au[0] * X + Au[1] * Y) * a;
821 const double s2 = (Av[0] * X + Av[1] * Y + Av[2] * Z) * a;
822 const double s3 = p[0] * (Bx * At[1] - By * At[0]) * a;
823 const double s4 = p[0] * (Bx * Av[0] + By * Av[1] + Bz * Av[2]) * a;
824
825 // Jacobian production
826 //
827 Jac[0] = Bx * Au[0] + By * Au[1]; // dL0/dL0
828 Jac[1] = Bx * Av[0] + By * Av[1] + Bz * Av[2]; // dL0/dL1
829 Jac[2] = 0.; // dL0/dPhi
830 Jac[3] = 0.; // dL0/dThe
831 Jac[5] = A[0] * Au[0] + A[1] * Au[1] + s1 * d; // dL1/dL0
832 Jac[6] = (A[0] * Av[0] + A[1] * Av[1] + A[2] * Av[2]) + s2 * d; // dL1/dL1
833 Jac[7] = s3 * d; // dL1/dPhi
834 Jac[8] = s4 * d; // dL1/dThe
835}
836
838// Step estimation to Plane
840double
841stepEstimatorToPlane(const double* ATH_RESTRICT S,
842 const double* ATH_RESTRICT P,
843 bool& Q)
844{
845 const double* r = &P[0]; // Start coordinate
846 const double* a = &P[3]; // Start direction
847
848 const double A = a[0] * S[0] + a[1] * S[1] + a[2] * S[2];
849 if (A == 0.) {
850 Q = false;
851 return 1000000.;
852 }
853 const double D = (S[3] - r[0] * S[0]) - (r[1] * S[1] + r[2] * S[2]);
854 Q = true;
855 return (D / A);
856}
857
859// Step estimation to Cylinder
861double
862stepEstimatorToCylinder(double* ATH_RESTRICT S,
863 const double* ATH_RESTRICT P,
864 bool& Q)
865{
866 const double* r = &P[0]; // Start coordinate
867 const double* a = &P[3]; // Start direction
868
869 double dx = r[0] - S[0];
870 double dy = r[1] - S[1];
871 double dz = r[2] - S[2];
872 double B = dx * S[3] + dy * S[4] + dz * S[5];
873 double C = a[0] * S[3] + a[1] * S[4] + a[2] * S[5];
874 const double ax = a[0] - S[3] * C;
875 dx -= (B * S[3]);
876 const double ay = a[1] - S[4] * C;
877 dy -= (B * S[4]);
878 const double az = a[2] - S[5] * C;
879 dz -= (B * S[5]);
880 double A = 2. * (ax * ax + ay * ay + az * az);
881 if (A == 0.) {
882 Q = false;
883 return 0.;
884 }
885 B = 2. * (ax * dx + ay * dy + az * dz);
886 double g = dx + dy + dz;
887 C = (g - S[6]) * (g + S[6]) - 2. * (dx * (dy + dz) + dy * dz);
888 double Sq = B * B - 2. * A * C;
889
890 double Smin = -B / A, Smax = Smin;
891
892 if (Sq > 0.) {
893
894 Sq = std::sqrt(Sq) / A;
895 if (B > 0.) {
896 Smin += Sq;
897 Smax -= Sq;
898 } else {
899 Smin -= Sq;
900 Smax += Sq;
901 }
902 } else {
903 if (std::fabs(Smax) < .1) {
904 Q = false;
905 return 0.;
906 }
907 }
908
909 Q = true;
910 const double inside = Smin * Smax;
911
912 if (S[8] != 0.) {
913
914 if (inside > 0. || S[8] > 0.)
915 return Smin;
916 if (S[7] >= 0.) {
917 if (Smin >= 0.)
918 return Smin;
919 return Smax;
920 }
921 if (Smin <= 0.)
922 return Smin;
923 return Smax;
924 }
925
926 if (inside < 0.) {
927
928 S[8] = -1.;
929 if (S[7] >= 0.) {
930 if (Smin >= 0.)
931 return Smin;
932 return Smax;
933 }
934 if (Smin <= 0.)
935 return Smin;
936 return Smax;
937 }
938
939 // if(std::fabs(Smin) < .001) {S[8]=-1.; return Smax;}
940
941 S[8] = 1.;
942 return Smin;
943}
944
946// Step estimation to Straight Line
948double
949stepEstimatorToStraightLine(const double* ATH_RESTRICT S,
950 const double* ATH_RESTRICT P,
951 bool& Q)
952{
953 const double* r = &P[0]; // Start coordinate
954 const double* a = &P[3]; // Start direction
955
956 double D = a[0] * S[3] + a[1] * S[4] + a[2] * S[5];
957 const double A = (1. - D) * (1. + D);
958 if (A == 0.) {
959 Q = true;
960 return 1000000.;
961 }
962 D = (r[0] - S[0]) * (D * S[3] - a[0]) + (r[1] - S[1]) * (D * S[4] - a[1]) +
963 (r[2] - S[2]) * (D * S[5] - a[2]);
964 Q = true;
965 return (D / A);
966}
967
969// Step estimation to Cone
971double
972stepEstimatorToCone(double* ATH_RESTRICT S,
973 const double* ATH_RESTRICT P,
974 bool& Q)
975{
976 const double* r = &P[0]; // Start coordinate
977 const double* a = &P[3]; // Start direction
978
979 const double dx = r[0] - S[0];
980 const double dy = r[1] - S[1];
981 const double dz = r[2] - S[2];
982 const double A = dx * S[3] + dy * S[4] + dz * S[5];
983 const double B = a[0] * S[3] + a[1] * S[4] + a[2] * S[5];
984 const double C = a[0] * dx + a[1] * dy + a[2] * dz;
985 const double k = S[6];
986
987 const double KB = 1. - k * B * B;
988 const double KABC = k * A * B - C;
989 double Smin, Smax;
990
991 if (KB != 0.) {
992
993 Smin = KABC / KB;
994 Smax = Smin;
995 double Sq = KABC * KABC + (k * A * A - dx * dx - dy * dy - dz * dz) * KB;
996 if (Sq >= 0.) {
997 Sq = std::sqrt(Sq) / KB;
998 if (KABC > 0.) {
999 Smin -= Sq;
1000 Smax += Sq;
1001 } else {
1002 Smin += Sq;
1003 Smax -= Sq;
1004 }
1005 int n = 2;
1006 if (A + B * Smin < 0.) {
1007 --n;
1008 Smin = Smax;
1009 }
1010 if (A + B * Smax < 0.) {
1011 --n;
1012 Smax = Smin;
1013 }
1014 if (!n) {
1015 Q = false;
1016 return 0.;
1017 }
1018
1019 } else {
1020 Q = false;
1021 return 0.;
1022 }
1023 } else {
1024 Smin = (dx * dx + dy * dy + dz * dz - k * A * A) / (2. * KABC);
1025 Smax = Smin;
1026 if (A + B * Smin < 0.) {
1027 Q = false;
1028 return 0.;
1029 }
1030 }
1031
1032 Q = true;
1033 const double inside = Smin * Smax;
1034
1035 if (S[8] != 0.) {
1036
1037 if (inside > 0. || S[8] > 0.)
1038 return Smin;
1039 if (S[7] >= 0.) {
1040 if (Smin >= 0.)
1041 return Smin;
1042 return Smax;
1043 }
1044 if (Smin <= 0.)
1045 return Smin;
1046 return Smax;
1047 }
1048
1049 if (inside < 0.) {
1050
1051 S[8] = -1.;
1052 if (S[7] >= 0.) {
1053 if (Smin >= 0.)
1054 return Smin;
1055 return Smax;
1056 }
1057 if (Smin <= 0.)
1058 return Smin;
1059 return Smax;
1060 }
1061
1062 S[8] = 1.;
1063 return Smin;
1064}
1065
1066
1067} //end of anonymous namespace
1068
1069
1071// Common transformation from local to global system coordinates for all
1072// surfaces for charged track parameters
1074
1075bool
1076Trk::RungeKuttaUtils::transformLocalToGlobal(bool useJac,
1077 const Trk::TrackParameters& Tp,
1078 double* ATH_RESTRICT P)
1079{
1080 return transformLocalToGlobal(
1081 useJac, &Tp.associatedSurface(), Tp.parameters(), P);
1082}
1083
1085// Common transformation from local to global system coordinates for all
1086// surfaces for neutral track parameters
1088
1089bool
1090Trk::RungeKuttaUtils::transformLocalToGlobal(bool useJac,
1091 const Trk::NeutralParameters& Tp,
1092 double* ATH_RESTRICT P)
1093{
1094 return transformLocalToGlobal(
1095 useJac, &Tp.associatedSurface(), Tp.parameters(), P);
1096}
1097
1099// Common transformation from local to global system coordinates for all
1100// surfaces for pattern parameters
1102
1103bool
1104Trk::RungeKuttaUtils::transformLocalToGlobal(
1105 bool useJac,
1106 const Trk::PatternTrackParameters& Tp,
1107 double* P)
1108{
1109 return transformLocalToGlobal(
1110 useJac, &Tp.associatedSurface(), Tp.parameters(), P);
1111}
1112
1114// Common transformation from global to local system coordinates for all
1115// surfaces
1117
1118void
1119Trk::RungeKuttaUtils::transformGlobalToLocal(double* ATH_RESTRICT P,
1120 double* ATH_RESTRICT par)
1121{
1122 par[2] = atan2(P[4], P[3]);
1123 par[3] = acos(P[5]);
1124 par[4] = P[6];
1125}
1126
1127void
1128Trk::RungeKuttaUtils::transformGlobalToLocal(const Trk::Surface* su,
1129 bool useJac,
1130 double* ATH_RESTRICT P,
1131 double* ATH_RESTRICT par,
1132 double* ATH_RESTRICT Jac)
1133{
1134 par[2] = atan2(P[4], P[3]);
1135 par[3] = acos(P[5]);
1136 par[4] = P[6];
1137
1138 const Trk::SurfaceType ty = su->type();
1139 const Amg::Transform3D& T = su->transform();
1140 switch (ty) {
1141 case Trk::SurfaceType::Plane: {
1142 transformGlobalToPlane(T, useJac, P, par, Jac);
1143 break;
1144 }
1145 // Perigee call the same as line
1146 case Trk::SurfaceType::Line:
1147 case Trk::SurfaceType::Perigee: {
1148 transformGlobalToLine(T, useJac, P, par, Jac);
1149 break;
1150 }
1151 case Trk::SurfaceType::Cylinder: {
1152 transformGlobalToCylinder(
1153 T,
1154 static_cast<const Trk::CylinderSurface*>(su)->bounds().r(),
1155 useJac,
1156 P,
1157 par,
1158 Jac);
1159 break;
1160 }
1161 case Trk::SurfaceType::Disc: {
1162 transformGlobalToDisc(T, useJac, P, par, Jac);
1163 break;
1164 }
1165 case Trk::SurfaceType::Cone: {
1166 transformGlobalToCone(
1167 T,
1168 static_cast<const Trk::ConeSurface*>(su)->bounds().tanAlpha(),
1169 useJac,
1170 P,
1171 par,
1172 Jac);
1173 break;
1174 }
1175 default: {
1176 break;
1177 }
1178 }
1179
1180 if (!useJac)
1181 return;
1182
1183 double P3, P4, C = P[3] * P[3] + P[4] * P[4];
1184 if (C > 1.e-20) {
1185 C = 1. / C;
1186 P3 = P[3] * C;
1187 P4 = P[4] * C;
1188 C = -std::sqrt(C);
1189 } else {
1190 C = -1.e10;
1191 P3 = 1.;
1192 P4 = 0.;
1193 }
1194
1195 Jac[10] = P3 * P[11] - P4 * P[10]; // dPhi/dL0
1196 Jac[11] = P3 * P[18] - P4 * P[17]; // dPhi/dL1
1197 Jac[12] = P3 * P[25] - P4 * P[24]; // dPhi/dPhi
1198 Jac[13] = P3 * P[32] - P4 * P[31]; // dPhi/dThe
1199 Jac[14] = P3 * P[39] - P4 * P[38]; // dPhi/dCM
1200 Jac[15] = C * P[12]; // dThe/dL0
1201 Jac[16] = C * P[19]; // dThe/dL1
1202 Jac[17] = C * P[26]; // dThe/dPhi
1203 Jac[18] = C * P[33]; // dThe/dThe
1204 Jac[19] = C * P[40]; // dThe/dCM
1205 Jac[20] = P[41]; // dCM /dCM
1206}
1207
1209// step estimation to surfaces
1211double
1212Trk::RungeKuttaUtils::stepEstimator(
1213 int kind,
1214 double* ATH_RESTRICT Su,
1215 const double* ATH_RESTRICT P,
1216 bool& Q)
1217{
1218 switch (kind) {
1219 case 0: {
1220 return stepEstimatorToStraightLine(Su, P, Q);
1221 }
1222 case 1: {
1223 return stepEstimatorToPlane(Su, P, Q);
1224 }
1225 case 2: {
1226 return stepEstimatorToCylinder(Su, P, Q);
1227 }
1228 case 3: {
1229 return stepEstimatorToCone(Su, P, Q);
1230 }
1231 default: {
1232 return 1000000.;
1233 }
1234 }
1235}
1236double
1237Trk::RungeKuttaUtils::stepEstimator(
1238 Trk::SurfaceType surfaceType,
1239 double* ATH_RESTRICT Su,
1240 const double* ATH_RESTRICT P,
1241 bool& Q)
1242{
1243 if (surfaceType == Trk::SurfaceType::Plane ||
1244 surfaceType == Trk::SurfaceType::Disc) {
1245 return stepEstimatorToPlane(Su, P, Q);
1246 }
1247 if (surfaceType == Trk::SurfaceType::Cylinder) {
1248 return stepEstimatorToCylinder(Su, P, Q);
1249 }
1250
1251 if (surfaceType == Trk::SurfaceType::Line ||
1252 surfaceType == Trk::SurfaceType::Perigee) {
1253 return stepEstimatorToStraightLine(Su, P, Q);
1254 }
1255
1256 if (surfaceType == Trk::SurfaceType::Cone) {
1257 return stepEstimatorToCone(Su, P, Q);
1258 }
1259 // presumably curvilinear
1260 return stepEstimatorToPlane(Su, P, Q);
1261}
1262
1264// Step estimation to surfaces
1265//
1266// Input iformation
1267//
1268// SU - input vector surfaces and conditions for boundary check
1269// DN - is work map
1270// Pinp - array of track parameters in global system coordinates before some
1271// step Pout - array of track parameters in global system coordinates before
1272// some step W - trajectory total pass for Pinp position So - max. step
1273// for straight line propagation model Nv - surface number in vector SU which
1274// we don't need to use
1275//
1276// Output information
1277//
1278// pair(step to surface,surface number from vector SU)
1279// next - condition to use Pout for next step
1281
1282std::pair<double, int>
1283Trk::RungeKuttaUtils::stepEstimator(
1284 std::vector<std::pair<const Trk::Surface*, Trk::BoundaryCheck>>& SU,
1285 std::multimap<double, int>& DN,
1286 const double* ATH_RESTRICT Pinp,
1287 const double* ATH_RESTRICT Pout,
1288 double W,
1289 double So,
1290 int Nv,
1291 bool& next)
1292{
1293 W = std::fabs(W);
1294 next = false;
1295 int N = -1;
1296 const double D[3] = { Pout[0] - Pinp[0], Pout[1] - Pinp[1], Pout[2] - Pinp[2] };
1297 double Smax = std::sqrt(D[0] * D[0] + D[1] * D[1] + D[2] * D[2]);
1298 double Sign = D[0] * Pinp[3] + D[1] * Pinp[4] + D[2] * Pinp[5];
1299 // The magnitude of the vector is essentially 0. No
1301 if (Smax < DBL_EPSILON) {
1302 next = true;
1303 return std::make_pair(0., -1);
1304 }
1305 Eigen::Map<const Amg::Vector3D> pos(&Pinp[0],3,1);
1306 const Amg::Vector3D dir(D[0] / Smax, D[1] / Smax, D[2] / Smax);
1307
1308 double Smin = 2. * Smax;
1309
1310 std::vector<std::pair<double, int>> LD;
1311 std::multimap<double, int>::iterator i = DN.begin(), ie = DN.end();
1312
1313 for (; i != ie; ++i) {
1314
1315 if ((*i).first > W + Smin)
1316 break;
1317
1318 int j = (*i).second;
1319 Trk::DistanceSolution ds =
1320 SU[j].first->straightLineDistanceEstimate(pos, dir, SU[j].second);
1321 LD.emplace_back(ds.currentDistance(false) + W, j);
1322
1323 const int n = ds.numberOfSolutions();
1324 if (!n)
1325 continue;
1326
1327 for (int i = 0; i != n; ++i) {
1328 double s;
1329 i == 0 ? s = ds.first() : s = ds.second();
1330 if (s < 0. || s > Smin || (j == Nv && s < So))
1331 continue;
1332 Smin = s;
1333 N = j;
1334 }
1335 }
1336
1337 if (!LD.empty()) {
1338
1339 DN.erase(DN.begin(), i);
1340 std::vector<std::pair<double, int>>::iterator l = LD.begin(), le = LD.end();
1341 for (; l != le; ++l)
1342 DN.insert((*l));
1343 }
1344
1345 if (N < 0) {
1346 next = true;
1347 return std::make_pair(0., -1);
1348 }
1349
1350 double Sm = Smin;
1351
1352 if (Sign < 0.)
1353 Sm = -Sm;
1354 if (Smin < So || (Smax - Smin) > 2. * So)
1355 return std::make_pair(Sm, N);
1356
1357 Eigen::Map<const Amg::Vector3D> posn(&Pout[0], 3, 1);
1358 Eigen::Map<const Amg::Vector3D> dirn(&Pout[3], 3, 1);
1359
1360 Trk::DistanceSolution ds =
1361 SU[N].first->straightLineDistanceEstimate(posn, dirn, SU[N].second);
1362 const int n = ds.numberOfSolutions();
1363 if (!n)
1364 return std::make_pair(Sm, N);
1365 Sm = 10000.;
1366
1367 for (int i = 0; i != n; ++i) {
1368
1369 double sa, s;
1370 i == 0 ? s = ds.first() : s = ds.second();
1371 sa = std::fabs(s);
1372
1373 if (s * Sign < 0.) {
1374 // if(sa < So ) {next = true; return std::make_pair(s,N);}
1375 next = true;
1376 return std::make_pair(s, N);
1377 }
1378 if (sa < Sm) {
1379 Sm = s;
1380 next = true;
1381 }
1382 }
1383 return std::make_pair(Sm, N);
1384}
1385
1387// New covariance matrix calculation from old matrix and jacobian
1388// We use Eigen so we use flatten to avoid out of line calls
1390ATH_FLATTEN
1391AmgSymMatrix(5) Trk::RungeKuttaUtils::newCovarianceMatrix(
1392 const double* ATH_RESTRICT J,
1393 const AmgSymMatrix(5) & ATH_RESTRICT M)
1394{
1395 AmgSymMatrix(5) nM;
1396 AmgSymMatrix(5)& m = nM;
1397
1398 Eigen::Map<const AmgVector(5)> JacMap0(&J[0], 5, 1);
1399 const AmgVector(5) a1 = M * JacMap0; //(5x5 * 5x1)
1400 m(0, 0) = a1.dot(JacMap0); // dot product
1401
1402 Eigen::Map<const AmgVector(5)> JacMap5(&J[5]);
1403 const AmgVector(5) a2 = M * JacMap5;
1404 m(1, 0) = a2.dot(JacMap0);
1405 m(1, 1) = a2.dot(JacMap5);
1406 m(0, 1) = m(1, 0);
1407
1408 Eigen::Map<const AmgVector(5)> JacMap10(&J[10], 5, 1);
1409 const AmgVector(5) a3 = M * JacMap10;
1410 m(2, 0) = a3.dot(JacMap0);
1411 m(2, 1) = a3.dot(JacMap5);
1412 m(2, 2) = a3.dot(JacMap10);
1413 m(0, 2) = m(2, 0);
1414 m(1, 2) = m(2, 1);
1415
1416 Eigen::Map<const AmgVector(5)> JacMap15(&J[15], 5, 1);
1417 const AmgVector(5) a4 = M * JacMap15;
1418 m(3, 0) = a4.dot(JacMap0);
1419 m(3, 1) = a4.dot(JacMap5);
1420 m(3, 2) = a4.dot(JacMap10);
1421 m(3, 3) = a4.dot(JacMap15);
1422 m(0, 3) = m(3, 0);
1423 m(1, 3) = m(3, 1);
1424 m(2, 3) = m(3, 2);
1425
1426 const AmgVector(5) a5 = M.row(4) * J[20];
1427 m(4, 0) = a5.dot(JacMap0);
1428 m(4, 1) = a5.dot(JacMap5);
1429 m(4, 2) = a5.dot(JacMap10);
1430 m(4, 3) = a5.dot(JacMap15);
1431 m(4, 4) = a5[4] * J[20];
1432 m(0, 4) = m(4, 0);
1433 m(1, 4) = m(4, 1);
1434 m(2, 4) = m(4, 2);
1435 m(3, 4) = m(4, 3);
1436
1437 return nM;
1438}
1439
1441// Tramsform from local to global for all track parameters
1443
1444bool
1445Trk::RungeKuttaUtils::transformLocalToGlobal(bool useJac,
1446 const Trk::Surface* Su,
1447 const AmgVector(5) & ATH_RESTRICT p,
1448 double* ATH_RESTRICT P)
1449{
1450 if (!Su)
1451 return false;
1452
1453 double Sf, Cf, Ce, Se;
1454 sincos(p[2], &Sf, &Cf);
1455 sincos(p[3], &Se, &Ce);
1456
1457 P[3] = Cf * Se; // Ax
1458 P[4] = Sf * Se; // Ay
1459 P[5] = Ce; // Az
1460 P[6] = p[4]; // CM
1461 if (std::fabs(P[6]) < .000000000000001) {
1462 P[6] < 0. ? P[6] = -.000000000000001 : P[6] = .000000000000001;
1463 }
1464
1465 if (useJac) {
1466
1467 // /dL1 | /dL2 | /dPhi | /dThe | /dCM |
1468 P[35] = 0.; // dX /
1469 P[36] = 0.; // dY /
1470 P[37] = 0.; // dZ /
1471 P[10] = 0.;
1472 P[17] = 0.;
1473 P[24] = -P[4];
1474 P[31] = Cf * Ce;
1475 P[38] = 0.; // dAx/
1476 P[11] = 0.;
1477 P[18] = 0.;
1478 P[25] = P[3];
1479 P[32] = Sf * Ce;
1480 P[39] = 0.; // dAy/
1481 P[12] = 0.;
1482 P[19] = 0.;
1483 P[26] = 0.;
1484 P[33] = -Se;
1485 P[40] = 0.; // dAz/
1486 P[13] = 0.;
1487 P[20] = 0.;
1488 P[27] = 0.;
1489 P[34] = 0.;
1490 P[41] = 1.; // dCM/
1491 P[42] = 0.;
1492 P[43] = 0.;
1493 P[44] = 0.; // d(Ax,Ay,Az)/ds
1494 }
1495
1496 const Trk::SurfaceType ty = Su->type();
1497 const Amg::Transform3D& T = Su->transform();
1498 switch (ty) {
1499 case Trk::SurfaceType::Plane: {
1500 transformPlaneToGlobal(useJac, T, p, P);
1501 return true;
1502 }
1503 case Trk::SurfaceType::Line:
1504 case Trk::SurfaceType::Perigee: {
1505 transformLineToGlobal(useJac, T, p, P);
1506 return true;
1507 }
1508 case Trk::SurfaceType::Cylinder: {
1509 transformCylinderToGlobal(
1510 useJac,
1511 T,
1512 static_cast<const Trk::CylinderSurface*>(Su)->bounds().r(),
1513 p,
1514 P);
1515 return true;
1516 }
1517 case Trk::SurfaceType::Disc: {
1518 transformDiscToGlobal(useJac, T, p, P);
1519 return true;
1520 }
1521 default: {
1522 return false;
1523 }
1524 }
1525}
1526
1528// Global position transformation to curvilinear system coordinate
1530
1531void
1532Trk::RungeKuttaUtils::transformGlobalToCurvilinear(bool useJac,
1533 double* ATH_RESTRICT P,
1534 double* ATH_RESTRICT par,
1535 double* ATH_RESTRICT Jac)
1536{
1537 par[0] = 0.;
1538 par[1] = 0.;
1539
1540 if (!useJac)
1541 return;
1542
1543 const double An = std::sqrt(P[3] * P[3] + P[4] * P[4]);
1544 double Ax[3];
1545 if (An != 0.) {
1546 Ax[0] = -P[4] / An;
1547 Ax[1] = P[3] / An;
1548 Ax[2] = 0.;
1549 } else {
1550 Ax[0] = 1.;
1551 Ax[1] = 0.;
1552 Ax[2] = 0.;
1553 }
1554
1555 const double Ay[3] = { -Ax[1] * P[5], Ax[0] * P[5], An };
1556 double S[3] = { P[3], P[4], P[5] };
1557
1558 double A = P[3] * S[0] + P[4] * S[1] + P[5] * S[2];
1559 if (A != 0.)
1560 A = 1. / A;
1561 S[0] *= A;
1562 S[1] *= A;
1563 S[2] *= A;
1564
1565 double s[5] = {};
1566 mult3x5Helper(s, S, &P[7]);
1567 globalToLocalVecHelper(P, s[0], s[1], s[2], s[3], s[4]);
1568
1569 double P3, P4, C = P[3] * P[3] + P[4] * P[4];
1570 if (C > 1.e-20) {
1571 C = 1. / C;
1572 P3 = P[3] * C;
1573 P4 = P[4] * C;
1574 C = -std::sqrt(C);
1575 } else {
1576 C = -1.e10;
1577 P3 = 1.;
1578 P4 = 0.;
1579 }
1580
1581 // Jacobian production
1582 //
1583 Jac[0] = Ax[0] * P[7] + Ax[1] * P[8]; // dL0/dL0
1584 Jac[1] = Ax[0] * P[14] + Ax[1] * P[15]; // dL0/dL1
1585 Jac[2] = Ax[0] * P[21] + Ax[1] * P[22]; // dL0/dPhi
1586 Jac[3] = Ax[0] * P[28] + Ax[1] * P[29]; // dL0/dThe
1587 Jac[4] = Ax[0] * P[35] + Ax[1] * P[36]; // dL0/dCM
1588 Jac[5] = Ay[0] * P[7] + Ay[1] * P[8] + Ay[2] * P[9]; // dL1/dL0
1589 Jac[6] = Ay[0] * P[14] + Ay[1] * P[15] + Ay[2] * P[16]; // dL1/dL1
1590 Jac[7] = Ay[0] * P[21] + Ay[1] * P[22] + Ay[2] * P[23]; // dL1/dPhi
1591 Jac[8] = Ay[0] * P[28] + Ay[1] * P[29] + Ay[2] * P[30]; // dL1/dThe
1592 Jac[9] = Ay[0] * P[35] + Ay[1] * P[36] + Ay[2] * P[37]; // dL1/dCM
1593 Jac[10] = P3 * P[11] - P4 * P[10]; // dPhi/dL0
1594 Jac[11] = P3 * P[18] - P4 * P[17]; // dPhi/dL1
1595 Jac[12] = P3 * P[25] - P4 * P[24]; // dPhi/dPhi
1596 Jac[13] = P3 * P[32] - P4 * P[31]; // dPhi/dThe
1597 Jac[14] = P3 * P[39] - P4 * P[38]; // dPhi/dCM
1598 Jac[15] = C * P[12]; // dThe/dL0
1599 Jac[16] = C * P[19]; // dThe/dL1
1600 Jac[17] = C * P[26]; // dThe/dPhi
1601 Jac[18] = C * P[33]; // dThe/dThe
1602 Jac[19] = C * P[40]; // dThe/dCM
1603 Jac[20] = P[41]; // dCM /dCM
1604}
1605
1607// Curvilinear covariance transformation to global system coordinate
1609
1610void
1611Trk::RungeKuttaUtils::transformCurvilinearToGlobal(double* ATH_RESTRICT p,
1612 double* ATH_RESTRICT P)
1613{
1614 double Sf, Cf, Ce, Se;
1615 sincos(p[2], &Sf, &Cf);
1616 sincos(p[3], &Se, &Ce);
1617
1618 const double Ax[3] = { -Sf, Cf, 0. };
1619 const double Ay[3] = { -Cf * Ce, -Sf * Ce, Se };
1620
1621 // /dL1 | /dL2 | /dPhi | /dThe | /dCM |
1622 //
1623 P[7] = Ax[0];
1624 P[14] = Ay[0];
1625 P[21] = 0.;
1626 P[28] = 0.;
1627 P[35] = 0.; // dX /
1628 P[8] = Ax[1];
1629 P[15] = Ay[1];
1630 P[22] = 0.;
1631 P[29] = 0.;
1632 P[36] = 0.; // dY /
1633 P[9] = Ax[2];
1634 P[16] = Ay[2];
1635 P[23] = 0.;
1636 P[30] = 0.;
1637 P[37] = 0.; // dZ /
1638 P[10] = 0.;
1639 P[17] = 0.;
1640 P[24] = -Sf * Se;
1641 P[31] = -Ay[0];
1642 P[38] = 0.; // dAx/
1643 P[11] = 0.;
1644 P[18] = 0.;
1645 P[25] = Cf * Se;
1646 P[32] = -Ay[1];
1647 P[39] = 0.; // dAy/
1648 P[12] = 0.;
1649 P[19] = 0.;
1650 P[26] = 0.;
1651 P[33] = -Ay[2];
1652 P[40] = 0.; // dAz/
1653 P[13] = 0.;
1654 P[20] = 0.;
1655 P[27] = 0.;
1656 P[34] = 0.;
1657 P[41] = 1.; // dCM/
1658 P[42] = 0.;
1659 P[43] = 0.;
1660 P[44] = 0.;
1661}
1662
1664// Jacobian of transformations from curvilinear to local system coordinates
1665// for Trk::TrackParameters
1667
1668void
1669Trk::RungeKuttaUtils::jacobianTransformCurvilinearToLocal(
1670 const Trk::TrackParameters& Tp,
1671 double* ATH_RESTRICT Jac)
1672{
1673 const AmgVector(5)& Vp = Tp.parameters();
1674 double P[23];
1675 P[0] = Vp[0];
1676 P[1] = Vp[1];
1677 P[2] = Vp[2];
1678 P[3] = Vp[3];
1679 jacobianTransformCurvilinearToLocal(P, &Tp.associatedSurface(), Jac);
1680}
1681
1683// Jacobian of transformations from curvilinear to local system coordinates
1684// for Trk::PatternTrackParameters
1686
1687void
1688Trk::RungeKuttaUtils::jacobianTransformCurvilinearToLocal(
1689 const Trk::PatternTrackParameters& Tp,
1690 double* ATH_RESTRICT Jac)
1691{
1692 double P[23];
1693 const AmgVector(5)& p = Tp.parameters();
1694 P[0] = p[0];
1695 P[1] = p[1];
1696 P[2] = p[2];
1697 P[3] = p[3];
1698 jacobianTransformCurvilinearToLocal(P, &Tp.associatedSurface(), Jac);
1699}
1700
1702// Jacobian of transformations from curvilinear to local system coordinates
1704
1705void
1706Trk::RungeKuttaUtils::jacobianTransformCurvilinearToLocal(
1707 double* ATH_RESTRICT P,
1708 const Trk::Surface* su,
1709 double* ATH_RESTRICT Jac)
1710{
1711 // Common for all surfaces terms of jacobian
1712 //
1713 Jac[4] = 0.; // dL0/dCM
1714 Jac[9] = 0.; // dL1/dCM
1715 Jac[10] = 0.; // dPhi/dL0
1716 Jac[11] = 0.; // dPhi/dL1
1717 Jac[12] = 1.; // dPhi/dPhi
1718 Jac[13] = 0.; // dPhi/dThe
1719 Jac[14] = 0.; // dPhi/dCM
1720 Jac[15] = 0.; // dThe/dL0
1721 Jac[16] = 0.; // dThe/dL1
1722 Jac[17] = 0.; // dThe/dPhi
1723 Jac[18] = 1.; // dThe/dThe
1724 Jac[19] = 0.; // dThe/dCM
1725 Jac[20] = 1.; // dCM /dCM
1726
1727 const Amg::Transform3D& T = su->transform();
1728
1729 double Sf, Cf, Ce, Se;
1730 sincos(P[2], &Sf, &Cf);
1731 sincos(P[3], &Se, &Ce);
1732
1733 P[4] = Cf * Se;
1734 P[5] = Sf * Se;
1735 P[6] = Ce; // At
1736 P[7] = -Sf;
1737 P[8] = Cf;
1738 P[9] = 0.; // Au
1739 P[10] = -Cf * Ce;
1740 P[11] = -Sf * Ce;
1741 P[12] = Se; // Av
1742 P[13] = T(0, 0);
1743 P[14] = T(1, 0);
1744 P[15] = T(2, 0); // Ax
1745 P[16] = T(0, 1);
1746 P[17] = T(1, 1);
1747 P[18] = T(2, 1); // Ay
1748 P[19] = T(0, 2);
1749 P[20] = T(1, 2);
1750 P[21] = T(2, 2); // Az
1751
1752 const Trk::SurfaceType ty = su->type();
1753
1754 switch (ty) {
1755 case Trk::SurfaceType::Plane: {
1756 jacobianTransformCurvilinearToPlane(P, Jac);
1757 return;
1758 }
1759 case Trk::SurfaceType::Line:
1760 case Trk::SurfaceType::Perigee: {
1761 jacobianTransformCurvilinearToStraightLine(P, Jac);
1762 return;
1763 }
1764 case Trk::SurfaceType::Cylinder: {
1765 P[22] = static_cast<const Trk::CylinderSurface*>(su)->bounds().r();
1766 jacobianTransformCurvilinearToCylinder(P, Jac);
1767 return;
1768 }
1769 case Trk::SurfaceType::Disc: {
1770 jacobianTransformCurvilinearToDisc(P, Jac);
1771 ;
1772 return;
1773 }
1774 default: {
1775 return;
1776 }
1777 }
1778
1779 Jac[0] = Jac[3] = 1.;
1780 Jac[1] = Jac[2] = 0.;
1781}
1782
1784// Fill distances map (initial map<distance,N of SU> preparation
1785// Output
1786// Step[0] - min negative step
1787// Step[1] - min positive step
1788// Step[2] - max distance
1789// int - vector SU element number with initial surface
1791
1792int
1793Trk::RungeKuttaUtils::fillDistancesMap(
1794 std::vector<std::pair<const Trk::Surface*, Trk::BoundaryCheck>>& SU,
1795 std::multimap<double, int>& DN,
1796 const double* ATH_RESTRICT Pinp,
1797 double W,
1798 const Trk::Surface* So,
1799 double* ATH_RESTRICT Step)
1800{
1801 int Ns = -1;
1802 DN.erase(DN.begin(), DN.end());
1803
1804 Eigen::Map<const Amg::Vector3D> pos(&Pinp[0], 3, 1);
1805 Eigen::Map<const Amg::Vector3D> dir(&Pinp[3], 3, 1);
1806
1807 Step[0] = -1.e+20;
1808 Step[1] = 1.e+20;
1809 Step[2] = 0.;
1810
1811 int N = SU.size();
1812 W = std::fabs(W);
1813
1814 for (int i = 0; i != N; ++i) {
1815
1816 const Trk::Surface* s = SU[i].first;
1817 if (!s)
1818 continue;
1819 if (s == So)
1820 Ns = i;
1821 Trk::DistanceSolution ds =
1822 s->straightLineDistanceEstimate(pos, dir, SU[i].second);
1823 double dist = ds.currentDistance(false);
1824 DN.insert(std::make_pair(dist + W, i));
1825
1826 if (dist > Step[2])
1827 Step[2] = dist;
1828
1829 int n = ds.numberOfSolutions();
1830 if (!n)
1831 continue;
1832 for (int i = 0; i != n; ++i) {
1833
1834 double st;
1835 i == 0 ? st = ds.first() : st = ds.second();
1836
1837 if (s == So && std::fabs(st) <= .001)
1838 continue;
1839
1840 if (st < 0.) {
1841 if (st > Step[0])
1842 Step[0] = st;
1843 } else {
1844 if (st < Step[1])
1845 Step[1] = st;
1846 }
1847 }
1848 }
1849 return Ns;
1850}
B0
static const int B0
Definition AtlasPID.h:122
AmgSymMatrix
#define AmgSymMatrix(dim)
Definition EventPrimitives.h:50
AmgVector
#define AmgVector(rows)
Definition EventPrimitives.h:55
Sign
int Sign(int in)
Definition JSSTaggerUtils.cxx:23
res
std::pair< std::vector< unsigned int >, bool > res
Definition JetGroupProductTest.cxx:11
s0
static Double_t s0
Definition LArPhysWaveHECTool.cxx:37
a
static Double_t a
Definition LArPhysWaveHECTool.cxx:38
Tp
static Double_t Tp(Double_t *t, Double_t *par)
Definition LArPhysWaveHECTool.cxx:354
P
static Double_t P(Double_t *tt, Double_t *par)
Definition LArPhysWaveHECTool.cxx:280
y
#define y
x
#define x
z
#define z
Trk::ConeSurface
Class for a conical surface in the ATLAS detector.
Definition ConeSurface.h:51
Trk::CylinderSurface
Class for a CylinderSurface in the ATLAS detector.
Definition CylinderSurface.h:55
Trk::DistanceSolution
Access to distance solutions.
Definition DistanceSolution.h:25
Trk::Surface
Abstract Base Class for tracking surfaces.
Definition Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/Surface.h:79
Trk::Surface::transform
const Amg::Transform3D & transform() const
Returns HepGeom::Transform3D by reference.
Trk::Surface::type
virtual constexpr SurfaceType type() const =0
Returns the Surface type to avoid dynamic casts.
r
int r
Definition globals.cxx:22
inline_hints.h
ATH_FLATTEN
#define ATH_FLATTEN
Definition inline_hints.h:52
C
struct color C
const
Amg::Transform3D
Eigen::Affine3d Transform3D
Definition GeoPrimitives.h:46
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition GeoPrimitives.h:47
AnalysisUtils::Delta::R
double R(const INavigable4Momentum *p1, const double v_eta, const double v_phi)
Definition AnalysisMisc.h:49
CxxUtils::vec
typename vecDetail::vec_typedef< T, N >::type vec
Define a nice alias for the vectorized type.
Definition vec.h:207
CxxUtils::vstore
ATH_ALWAYS_INLINE void vstore(vec_type_t< VEC > *dst, const VEC &src)
Definition vec.h:290
CxxUtils::vload
ATH_ALWAYS_INLINE void vload(VEC &dst, vec_type_t< VEC > const *src)
Definition vec.h:272
D3PDMakerTestInstan::vec2
std::vector< D3PDTest::MyVec2 > vec2
Definition D3PDMakerTestDict.h:14
PrintTrkAnaSummary.l
l
Printing final latex table to .tex output file.
Definition PrintTrkAnaSummary.py:370
RungeKuttaUtils
Trk::RungeKuttaUtils is set algorithms for track parameters transformation from local to global and g...
Trk::RungeKuttaUtils::transformGlobalToCurvilinear
void transformGlobalToCurvilinear(bool, double *ATH_RESTRICT, double *ATH_RESTRICT, double *ATH_RESTRICT)
Trk::RungeKuttaUtils::transformLocalToGlobal
bool transformLocalToGlobal(bool, const Trk::TrackParameters &, double *ATH_RESTRICT)
Trk::RungeKuttaUtils::transformGlobalToLocal
void transformGlobalToLocal(double *ATH_RESTRICT, double *ATH_RESTRICT)
Trk::RungeKuttaUtils::fillDistancesMap
int fillDistancesMap(std::vector< std::pair< const Trk::Surface *, Trk::BoundaryCheck > > &, std::multimap< double, int > &, const double *ATH_RESTRICT, double, const Trk::Surface *, double *ATH_RESTRICT)
Trk::RungeKuttaUtils::jacobianTransformCurvilinearToLocal
void jacobianTransformCurvilinearToLocal(const Trk::TrackParameters &, double *ATH_RESTRICT)
Trk::RungeKuttaUtils::stepEstimator
double stepEstimator(int kind, double *ATH_RESTRICT Su, const double *ATH_RESTRICT P, bool &Q)
Definition RungeKuttaUtils.cxx:1212
Trk::RungeKuttaUtils::transformCurvilinearToGlobal
void transformCurvilinearToGlobal(double *ATH_RESTRICT, double *ATH_RESTRICT)
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition FakeTrackBuilder.h:9
Trk::SurfaceType
SurfaceType
This enumerator simplifies the persistency & calculations,.
Definition SurfaceTypes.h:17
Trk::next
@ next
Definition BinningData.h:33
Trk::NeutralParameters
ParametersBase< NeutralParametersDim, Neutral > NeutralParameters
Definition NeutralParameters.h:26
Trk::TrackParameters
ParametersBase< TrackParametersDim, Charged > TrackParameters
Definition Tracking/TrkEvent/TrkParameters/TrkParameters/TrackParameters.h:27
dq_defect_virtual_defect_validation.d1
d1
Definition dq_defect_virtual_defect_validation.py:79
dq_defect_virtual_defect_validation.d2
d2
Definition dq_defect_virtual_defect_validation.py:81
dqt_zlumi_alleff_HIST.B
B
Definition dqt_zlumi_alleff_HIST.py:110
hist_file_dump.d
d
Definition hist_file_dump.py:142
ATH_RESTRICT
#define ATH_RESTRICT
Definition restrict.h:31
A
hold the test vectors and ease the comparison
CxxUtils::sincos
Helper to simultaneously calculate sin and cos of the same angle.
Definition sincos.h:39
vec.h
Vectorization helpers.
Generated on for ATLAS Offline Software by doxygen 1.14.0