ATLAS Offline Software
Loading...
Searching...
No Matches
RungeKuttaPropagator.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3*/
4
6// RungeKuttaPropagator.cxx, (c) ATLAS Detector software
8
9
10#include "TrkExRungeKuttaPropagator/RungeKuttaPropagator.h"
11//
12#include "TrkExUtils/JacobianHelper.h"
13#include "TrkExUtils/RungeKuttaUtils.h"
14//
15#include "TrkSurfaces/ConeSurface.h"
16#include "TrkSurfaces/CylinderSurface.h"
17#include "TrkSurfaces/DiscSurface.h"
18#include "TrkSurfaces/PerigeeSurface.h"
19#include "TrkSurfaces/PlaneSurface.h"
20#include "TrkSurfaces/StraightLineSurface.h"
21//
22#include "TrkGeometry/MagneticFieldProperties.h"
23//
24#include "EventPrimitives/EventPrimitivesCovarianceHelpers.h"
25//
26#include "TrkEventPrimitives/TransportJacobian.h"
27//
28#include "TrkPatternParameters/PatternTrackParameters.h"
29
30#include "CxxUtils/restrict.h"
31
32namespace {
33/*
34 * All internal implementation methods
35 * are in the anonymous namespace
36 */
37
39// get magnetic field methods
41
42using Cache = Trk::RungeKuttaPropagator::Cache;
43
44void
45getField(Cache& cache, const double* ATH_RESTRICT R, double* ATH_RESTRICT H)
46{
47
48 if (cache.m_solenoid) {
49 cache.m_fieldCache.getFieldZR(R, H);
50 } else {
51 cache.m_fieldCache.getField(R, H);
52 }
53}
54
55void
56getFieldGradient(Cache& cache,
57 const double* ATH_RESTRICT R,
58 double* ATH_RESTRICT H,
59 double* ATH_RESTRICT dH)
60{
61 if (cache.m_solenoid) {
62 cache.m_fieldCache.getFieldZR(R, H, dH);
63 } else {
64 cache.m_fieldCache.getField(R, H, dH);
65 }
66}
67
69// Runge Kutta trajectory model (units->mm,MeV,kGauss)
70// Uses Nystroem algorithm (See Handbook Net. Bur. ofStandards, procedure 25.5.20)
72
73double
74rungeKuttaStep(Cache& cache, bool Jac, double S, double* ATH_RESTRICT P, bool& InS)
75{
76 double* R = &P[0]; // Coordinates
77 double* A = &P[3]; // Directions
78 double* sA = &P[42];
79 const double Pi = 149.89626 * P[6]; // Invert mometum/2.
80 double dltm = cache.m_dlt * .03;
81
82 double f0[3];
83 double f[3];
84 if (cache.m_newfield)
85 getField(cache, R, f0);
86 else {
87 f0[0] = cache.m_field[0];
88 f0[1] = cache.m_field[1];
89 f0[2] = cache.m_field[2];
90 }
91
92 bool Helix = false;
93 if (std::abs(S) < cache.m_helixStep)
94 Helix = true;
95 while (S != 0.) {
96
97 const double S3 = (1. / 3.) * S;
98 const double S4 = .25 * S;
99 const double PS2 = Pi * S;
100
101 // First point
102 //
103 const double H0[3] = { f0[0] * PS2, f0[1] * PS2, f0[2] * PS2 };
104 const double A0 = A[1] * H0[2] - A[2] * H0[1];
105 const double B0 = A[2] * H0[0] - A[0] * H0[2];
106 const double C0 = A[0] * H0[1] - A[1] * H0[0];
107 const double A2 = A0 + A[0];
108 const double B2 = B0 + A[1];
109 const double C2 = C0 + A[2];
110 const double A1 = A2 + A[0];
111 const double B1 = B2 + A[1];
112 const double C1 = C2 + A[2];
113
114 // Second point
115 //
116 if (!Helix) {
117 const double gP[3] = { R[0] + A1 * S4, R[1] + B1 * S4, R[2] + C1 * S4 };
118 getField(cache, gP, f);
119 } else {
120 f[0] = f0[0];
121 f[1] = f0[1];
122 f[2] = f0[2];
123 }
124
125 const double H1[3] = { f[0] * PS2, f[1] * PS2, f[2] * PS2 };
126 const double A3 = (A[0] + B2 * H1[2]) - C2 * H1[1];
127 const double B3 = (A[1] + C2 * H1[0]) - A2 * H1[2];
128 const double C3 = (A[2] + A2 * H1[1]) - B2 * H1[0];
129 const double A4 = (A[0] + B3 * H1[2]) - C3 * H1[1];
130 const double B4 = (A[1] + C3 * H1[0]) - A3 * H1[2];
131 const double C4 = (A[2] + A3 * H1[1]) - B3 * H1[0];
132 const double A5 = 2. * A4 - A[0];
133 const double B5 = 2. * B4 - A[1];
134 const double C5 = 2. * C4 - A[2];
135
136 // Last point
137 //
138 if (!Helix) {
139 const double gP[3] = { R[0] + S * A4, R[1] + S * B4, R[2] + S * C4 };
140 getField(cache, gP, f);
141 } else {
142 f[0] = f0[0];
143 f[1] = f0[1];
144 f[2] = f0[2];
145 }
146
147 const double H2[3] = { f[0] * PS2, f[1] * PS2, f[2] * PS2 };
148 const double A6 = B5 * H2[2] - C5 * H2[1];
149 const double B6 = C5 * H2[0] - A5 * H2[2];
150 const double C6 = A5 * H2[1] - B5 * H2[0];
151
152 // Test approximation quality on give step and possible step reduction
153 //
154 const double EST = std::abs((A1 + A6) - (A3 + A4)) +
155 std::abs((B1 + B6) - (B3 + B4)) +
156 std::abs((C1 + C6) - (C3 + C4));
157 if (EST > cache.m_dlt) {
158 S *= .5;
159 dltm = 0.;
160 continue;
161 }
162 EST < dltm ? InS = true : InS = false;
163
164 // Parameters calculation
165 //
166 const double A00 = A[0];
167 const double A11 = A[1];
168 const double A22 = A[2];
169
170 const double Aarr[3]{ A00, A11, A22 };
171 const double A0arr[3]{ A0, B0, C0 };
172 const double A3arr[3]{ A3, B3, C3 };
173 const double A4arr[3]{ A4, B4, C4 };
174 const double A6arr[3]{ A6, B6, C6 };
175
176 A[0] = 2. * A3 + (A0 + A5 + A6);
177 A[1] = 2. * B3 + (B0 + B5 + B6);
178 A[2] = 2. * C3 + (C0 + C5 + C6);
179
180 double D = (A[0] * A[0] + A[1] * A[1]) + (A[2] * A[2] - 9.);
181 const double Sl = 2. / S;
182 D = (1. / 3.) - ((1. / 648.) * D) * (12. - D);
183
184 R[0] += (A2 + A3 + A4) * S3;
185 R[1] += (B2 + B3 + B4) * S3;
186 R[2] += (C2 + C3 + C4) * S3;
187 A[0] *= D;
188 A[1] *= D;
189 A[2] *= D;
190 sA[0] = A6 * Sl;
191 sA[1] = B6 * Sl;
192 sA[2] = C6 * Sl;
193
194 cache.m_field[0] = f[0];
195 cache.m_field[1] = f[1];
196 cache.m_field[2] = f[2];
197 cache.m_newfield = false;
198
199 if (!Jac)
200 return S;
201
202 // Jacobian calculation - outsourced into a helper also used by SiTrajectoryElement_xk
203 Trk::propJacobian(P, H0, H1, H2, Aarr, A0arr, A3arr, A4arr, A6arr, S3);
204 return S;
205 }
206 return S;
207}
208
210// Straight line trajectory model
212
213double
214straightLineStep(bool Jac, double S, double* ATH_RESTRICT P)
215{
216 double* R = &P[0]; // Start coordinates
217 const double* A = &P[3]; // Start directions
218 double* sA = &P[42];
219
220 // Track parameters in last point
221 //
222 R[0] += (A[0] * S);
223 R[1] += (A[1] * S);
224 R[2] += (A[2] * S);
225 if (!Jac)
226 return S;
227
228 // Derivatives of track parameters in last point
229 //
230 for (int i = 7; i < 42; i += 7) {
231
232 double* dR = &P[i];
233 const double* dA = &P[i + 3];
234 dR[0] += (dA[0] * S);
235 dR[1] += (dA[1] * S);
236 dR[2] += (dA[2] * S);
237 }
238 sA[0] = sA[1] = sA[2] = 0.;
239 return S;
240}
241
243// Runge Kutta trajectory model (units->mm,MeV,kGauss)
244// Uses Nystroem algorithm (See Handbook Net. Bur. ofStandards, procedure 25.5.20)
245// Where magnetic field information iS
246// f[ 0],f[ 1],f[ 2] - Hx , Hy and Hz of the magnetic field
247// f[ 3],f[ 4],f[ 5] - dHx/dx, dHx/dy and dHx/dz
248// f[ 6],f[ 7],f[ 8] - dHy/dx, dHy/dy and dHy/dz
249// f[ 9],f[10],f[11] - dHz/dx, dHz/dy and dHz/dz
250//
252double
253rungeKuttaStepWithGradient(Cache& cache, double S, double* ATH_RESTRICT P, bool& InS)
254{
255 const double C33 = 1. / 3.;
256 double* R = &P[0]; // Coordinates
257 double* A = &P[3]; // Directions
258 double* sA = &P[42];
259 const double Pi = 149.89626 * P[6]; // Invert mometum/2.
260 double dltm = cache.m_dlt * .03;
261
262 double f0[3];
263 double f1[3];
264 double f2[3];
265 double g0[9];
266 double g1[9];
267 double g2[9];
268 double H0[12];
269 double H1[12];
270 double H2[12];
271 getFieldGradient(cache, R, f0, g0);
272
273 while (S != 0.) {
274
275 const double S3 = C33 * S;
276 const double S4 = .25 * S;
277 const double PS2 = Pi * S;
278
279 // First point
280 //
281 H0[0] = f0[0] * PS2;
282 H0[1] = f0[1] * PS2;
283 H0[2] = f0[2] * PS2;
284 const double A0 = A[1] * H0[2] - A[2] * H0[1];
285 const double B0 = A[2] * H0[0] - A[0] * H0[2];
286 const double C0 = A[0] * H0[1] - A[1] * H0[0];
287 const double A2 = A[0] + A0;
288 const double B2 = A[1] + B0;
289 const double C2 = A[2] + C0;
290 const double A1 = A2 + A[0];
291 const double B1 = B2 + A[1];
292 const double C1 = C2 + A[2];
293
294 // Second point
295 //
296 const double gP1[3] = { R[0] + A1 * S4, R[1] + B1 * S4, R[2] + C1 * S4 };
297 getFieldGradient(cache, gP1, f1, g1);
298
299 H1[0] = f1[0] * PS2;
300 H1[1] = f1[1] * PS2;
301 H1[2] = f1[2] * PS2;
302 const double A3 = B2 * H1[2] - C2 * H1[1] + A[0];
303 const double B3 = C2 * H1[0] - A2 * H1[2] + A[1];
304 const double C3 = A2 * H1[1] - B2 * H1[0] + A[2];
305 const double A4 = B3 * H1[2] - C3 * H1[1] + A[0];
306 const double B4 = C3 * H1[0] - A3 * H1[2] + A[1];
307 const double C4 = A3 * H1[1] - B3 * H1[0] + A[2];
308 const double A5 = A4 - A[0] + A4;
309 const double B5 = B4 - A[1] + B4;
310 const double C5 = C4 - A[2] + C4;
311
312 // Last point
313 //
314 const double gP2[3] = { R[0] + S * A4, R[1] + S * B4, R[2] + S * C4 };
315 getFieldGradient(cache, gP2, f2, g2);
316
317 H2[0] = f2[0] * PS2;
318 H2[1] = f2[1] * PS2;
319 H2[2] = f2[2] * PS2;
320 const double A6 = B5 * H2[2] - C5 * H2[1];
321 const double B6 = C5 * H2[0] - A5 * H2[2];
322 const double C6 = A5 * H2[1] - B5 * H2[0];
323
324 // Test approximation quality on give step and possible step reduction
325 //
326 const double EST = std::abs((A1 + A6) - (A3 + A4)) + std::abs((B1 + B6) - (B3 + B4)) +
327 std::abs((C1 + C6) - (C3 + C4));
328 if (EST > cache.m_dlt) {
329 S *= .5;
330 dltm = 0.;
331 continue;
332 }
333 EST < dltm ? InS = true : InS = false;
334
335 // Parameters calculation
336 //
337 const double A00 = A[0];
338 const double A11 = A[1];
339 const double A22 = A[2];
340 R[0] += (A2 + A3 + A4) * S3;
341 A[0] = ((A0 + 2. * A3) + (A5 + A6)) * C33;
342 R[1] += (B2 + B3 + B4) * S3;
343 A[1] = ((B0 + 2. * B3) + (B5 + B6)) * C33;
344 R[2] += (C2 + C3 + C4) * S3;
345 A[2] = ((C0 + 2. * C3) + (C5 + C6)) * C33;
346 const double CBA = 1. / std::sqrt(A[0] * A[0] + A[1] * A[1] + A[2] * A[2]);
347 A[0] *= CBA;
348 A[1] *= CBA;
349 A[2] *= CBA;
350
351 const double Sl = 2. / S;
352 sA[0] = A6 * Sl;
353 sA[1] = B6 * Sl;
354 sA[2] = C6 * Sl;
355
356 // Jacobian calculation
357 //
358 for (int i = 3; i != 12; ++i) {
359 H0[i] = g0[i - 3] * PS2;
360 H1[i] = g1[i - 3] * PS2;
361 H2[i] = g2[i - 3] * PS2;
362 }
363
364 for (int i = 7; i < 35; i += 7) {
365
366 double* dR = &P[i];
367 double* dA = &P[i + 3];
368 double dH0 = H0[3] * dR[0] + H0[4] * dR[1] + H0[5] * dR[2]; // dHx/dp
369 double dH1 = H0[6] * dR[0] + H0[7] * dR[1] + H0[8] * dR[2]; // dHy/dp
370 double dH2 = H0[9] * dR[0] + H0[10] * dR[1] + H0[11] * dR[2]; // dHz/dp
371
372 const double dA0 = (H0[2] * dA[1] - H0[1] * dA[2]) + (A[1] * dH2 - A[2] * dH1); // dA0/dp
373 const double dB0 = (H0[0] * dA[2] - H0[2] * dA[0]) + (A[2] * dH0 - A[0] * dH2); // dB0/dp
374 const double dC0 = (H0[1] * dA[0] - H0[0] * dA[1]) + (A[0] * dH1 - A[1] * dH0); // dC0/dp
375 const double dA2 = dA0 + dA[0];
376 double dX = dR[0] + (dA2 + dA[0]) * S4; // dX /dp
377 const double dB2 = dB0 + dA[1];
378 double dY = dR[1] + (dB2 + dA[1]) * S4; // dY /dp
379 const double dC2 = dC0 + dA[2];
380 double dZ = dR[2] + (dC2 + dA[2]) * S4; // dZ /dp
381 dH0 = H1[3] * dX + H1[4] * dY + H1[5] * dZ; // dHx/dp
382 dH1 = H1[6] * dX + H1[7] * dY + H1[8] * dZ; // dHy/dp
383 dH2 = H1[9] * dX + H1[10] * dY + H1[11] * dZ; // dHz/dp
384
385 const double dA3 = (dA[0] + dB2 * H1[2] - dC2 * H1[1]) + (B2 * dH2 - C2 * dH1); // dA3/dp
386 const double dB3 = (dA[1] + dC2 * H1[0] - dA2 * H1[2]) + (C2 * dH0 - A2 * dH2); // dB3/dp
387 const double dC3 = (dA[2] + dA2 * H1[1] - dB2 * H1[0]) + (A2 * dH1 - B2 * dH0); // dC3/dp
388 const double dA4 = (dA[0] + dB3 * H1[2] - dC3 * H1[1]) + (B3 * dH2 - C3 * dH1); // dA4/dp
389 const double dB4 = (dA[1] + dC3 * H1[0] - dA3 * H1[2]) + (C3 * dH0 - A3 * dH2); // dB4/dp
390 const double dC4 = (dA[2] + dA3 * H1[1] - dB3 * H1[0]) + (A3 * dH1 - B3 * dH0); // dC4/dp
391 const double dA5 = dA4 + dA4 - dA[0];
392 dX = dR[0] + dA4 * S; // dX /dp
393 const double dB5 = dB4 + dB4 - dA[1];
394 dY = dR[1] + dB4 * S; // dY /dp
395 const double dC5 = dC4 + dC4 - dA[2];
396 dZ = dR[2] + dC4 * S; // dZ /dp
397 dH0 = H2[3] * dX + H2[4] * dY + H2[5] * dZ; // dHx/dp
398 dH1 = H2[6] * dX + H2[7] * dY + H2[8] * dZ; // dHy/dp
399 dH2 = H2[9] * dX + H2[10] * dY + H2[11] * dZ; // dHz/dp
400
401 const double dA6 = (dB5 * H2[2] - dC5 * H2[1]) + (B5 * dH2 - C5 * dH1); // dA6/dp
402 const double dB6 = (dC5 * H2[0] - dA5 * H2[2]) + (C5 * dH0 - A5 * dH2); // dB6/dp
403 const double dC6 = (dA5 * H2[1] - dB5 * H2[0]) + (A5 * dH1 - B5 * dH0); // dC6/dp
404 dR[0] += (dA2 + dA3 + dA4) * S3;
405 dA[0] = ((dA0 + 2. * dA3) + (dA5 + dA6)) * C33;
406 dR[1] += (dB2 + dB3 + dB4) * S3;
407 dA[1] = ((dB0 + 2. * dB3) + (dB5 + dB6)) * C33;
408 dR[2] += (dC2 + dC3 + dC4) * S3;
409 dA[2] = ((dC0 + 2. * dC3) + (dC5 + dC6)) * C33;
410 }
411
412 double* dR = &P[35];
413 double* dA = &P[38];
414
415 double dH0 = H0[3] * dR[0] + H0[4] * dR[1] + H0[5] * dR[2]; // dHx/dp
416 double dH1 = H0[6] * dR[0] + H0[7] * dR[1] + H0[8] * dR[2]; // dHy/dp
417 double dH2 = H0[9] * dR[0] + H0[10] * dR[1] + H0[11] * dR[2]; // dHz/dp
418 const double dA0 = (H0[2] * dA[1] - H0[1] * dA[2]) + (A[1] * dH2 - A[2] * dH1 + A0); // dA0/dp
419 const double dB0 = (H0[0] * dA[2] - H0[2] * dA[0]) + (A[2] * dH0 - A[0] * dH2 + B0); // dB0/dp
420 const double dC0 = (H0[1] * dA[0] - H0[0] * dA[1]) + (A[0] * dH1 - A[1] * dH0 + C0); // dC0/dp
421 const double dA2 = dA0 + dA[0];
422 double dX = dR[0] + (dA2 + dA[0]) * S4; // dX /dp
423 const double dB2 = dB0 + dA[1];
424 double dY = dR[1] + (dB2 + dA[1]) * S4; // dY /dp
425 const double dC2 = dC0 + dA[2];
426 double dZ = dR[2] + (dC2 + dA[2]) * S4; // dZ /dp
427 dH0 = H1[3] * dX + H1[4] * dY + H1[5] * dZ; // dHx/dp
428 dH1 = H1[6] * dX + H1[7] * dY + H1[8] * dZ; // dHy/dp
429 dH2 = H1[9] * dX + H1[10] * dY + H1[11] * dZ; // dHz/dp
430
431 const double dA3 =
432 (dA[0] + dB2 * H1[2] - dC2 * H1[1]) + ((B2 * dH2 - C2 * dH1) + (A3 - A00)); // dA3/dp
433 const double dB3 =
434 (dA[1] + dC2 * H1[0] - dA2 * H1[2]) + ((C2 * dH0 - A2 * dH2) + (B3 - A11)); // dB3/dp
435 const double dC3 =
436 (dA[2] + dA2 * H1[1] - dB2 * H1[0]) + ((A2 * dH1 - B2 * dH0) + (C3 - A22)); // dC3/dp
437 const double dA4 =
438 (dA[0] + dB3 * H1[2] - dC3 * H1[1]) + ((B3 * dH2 - C3 * dH1) + (A4 - A00)); // dA4/dp
439 const double dB4 =
440 (dA[1] + dC3 * H1[0] - dA3 * H1[2]) + ((C3 * dH0 - A3 * dH2) + (B4 - A11)); // dB4/dp
441 const double dC4 =
442 (dA[2] + dA3 * H1[1] - dB3 * H1[0]) + ((A3 * dH1 - B3 * dH0) + (C4 - A22)); // dC4/dp
443 const double dA5 = dA4 + dA4 - dA[0];
444 dX = dR[0] + dA4 * S; // dX /dp
445 const double dB5 = dB4 + dB4 - dA[1];
446 dY = dR[1] + dB4 * S; // dY /dp
447 const double dC5 = dC4 + dC4 - dA[2];
448 dZ = dR[2] + dC4 * S; // dZ /dp
449 dH0 = H2[3] * dX + H2[4] * dY + H2[5] * dZ; // dHx/dp
450 dH1 = H2[6] * dX + H2[7] * dY + H2[8] * dZ; // dHy/dp
451 dH2 = H2[9] * dX + H2[10] * dY + H2[11] * dZ; // dHz/dp
452
453 const double dA6 = (dB5 * H2[2] - dC5 * H2[1]) + (B5 * dH2 - C5 * dH1 + A6); // dA6/dp
454 const double dB6 = (dC5 * H2[0] - dA5 * H2[2]) + (C5 * dH0 - A5 * dH2 + B6); // dB6/dp
455 const double dC6 = (dA5 * H2[1] - dB5 * H2[0]) + (A5 * dH1 - B5 * dH0 + C6); // dC6/dp
456
457 dR[0] += (dA2 + dA3 + dA4) * S3;
458 dA[0] = ((dA0 + 2. * dA3) + (dA5 + dA6)) * C33;
459 dR[1] += (dB2 + dB3 + dB4) * S3;
460 dA[1] = ((dB0 + 2. * dB3) + (dB5 + dB6)) * C33;
461 dR[2] += (dC2 + dC3 + dC4) * S3;
462 dA[2] = ((dC0 + 2. * dC3) + (dC5 + dC6)) * C33;
463 return S;
464 }
465 return S;
466}
467
469// Test new cross point
471bool
472newCrossPoint(const Trk::CylinderSurface& Su,
473 const double* ATH_RESTRICT Ro,
474 const double* ATH_RESTRICT P)
475{
476 const double pi = 3.1415927;
477 const double pi2 = 2. * pi;
478 const Amg::Transform3D& T = Su.transform();
479 const double Ax[3] = { T(0, 0), T(1, 0), T(2, 0) };
480 const double Ay[3] = { T(0, 1), T(1, 1), T(2, 1) };
481
482 const double R = Su.bounds().r();
483 double x = Ro[0] - T(0, 3);
484 double y = Ro[1] - T(1, 3);
485 double z = Ro[2] - T(2, 3);
486
487 double RC = x * Ax[0] + y * Ax[1] + z * Ax[2];
488 double RS = x * Ay[0] + y * Ay[1] + z * Ay[2];
489
490 if ((RC * RC + RS * RS) <= (R * R))
491 return false;
492
493 x = P[0] - T(0, 3);
494 y = P[1] - T(1, 3);
495 z = P[2] - T(2, 3);
496 RC = x * Ax[0] + y * Ax[1] + z * Ax[2];
497 RS = x * Ay[0] + y * Ay[1] + z * Ay[2];
498 double dF = std::abs(std::atan2(RS, RC) - Su.bounds().averagePhi());
499 if (dF > pi)
500 dF = pi2 - pi;
501 return dF > Su.bounds().halfPhiSector();
502}
503
505// Build new track parameters without propagation
507std::unique_ptr<Trk::TrackParameters>
508buildTrackParametersWithoutPropagation(const Trk::TrackParameters& Tp,
509 double* ATH_RESTRICT Jac)
510{
511 Jac[0] = Jac[6] = Jac[12] = Jac[18] = Jac[20] = 1.;
512 Jac[1] = Jac[2] = Jac[3] = Jac[4] = Jac[5] = Jac[7] = Jac[8] = Jac[9] =
513 Jac[10] = Jac[11] = Jac[13] = Jac[14] = Jac[15] = Jac[16] = Jac[17] = Jac[19] = 0.;
514 return std::unique_ptr<Trk::TrackParameters>(Tp.clone());
515}
516
518// Build new neutral track parameters without propagation
520std::unique_ptr<Trk::NeutralParameters>
521buildTrackParametersWithoutPropagation(const Trk::NeutralParameters& Tp,
522 double* ATH_RESTRICT Jac)
523{
524 Jac[0] = Jac[6] = Jac[12] = Jac[18] = Jac[20] = 1.;
525 Jac[1] = Jac[2] = Jac[3] = Jac[4] = Jac[5] = Jac[7] = Jac[8] = Jac[9] =
526 Jac[10] = Jac[11] = Jac[13] = Jac[14] = Jac[15] = Jac[16] = Jac[17] = Jac[19] = 0.;
527 return std::unique_ptr<Trk::NeutralParameters>(Tp.clone());
528}
529
531// Track parameters in cross point preparation
533std::unique_ptr<Trk::TrackParameters>
534crossPoint(const Trk::TrackParameters& Tp,
535 std::vector<Trk::DestSurf>& SU,
536 std::vector<unsigned int>& So,
537 double* ATH_RESTRICT P,
538 std::pair<double, int>& SN)
539{
540 double* R = &P[0]; // Start coordinates
541 double* A = &P[3]; // Start directions
542 const double* SA = &P[42]; // d(directions)/dStep
543 const double Step = SN.first;
544 int N = SN.second;
545
546 double As[3];
547 double Rs[3];
548
549 As[0] = A[0] + SA[0] * Step;
550 As[1] = A[1] + SA[1] * Step;
551 As[2] = A[2] + SA[2] * Step;
552 const double CBA = 1. / std::sqrt(As[0] * As[0] + As[1] * As[1] + As[2] * As[2]);
553
554 Rs[0] = R[0] + Step * (As[0] - .5 * Step * SA[0]);
555 As[0] *= CBA;
556 Rs[1] = R[1] + Step * (As[1] - .5 * Step * SA[1]);
557 As[1] *= CBA;
558 Rs[2] = R[2] + Step * (As[2] - .5 * Step * SA[2]);
559 As[2] *= CBA;
560
561 const Amg::Vector3D pos(Rs[0], Rs[1], Rs[2]);
562 const Amg::Vector3D dir(As[0], As[1], As[2]);
563
564 const Trk::DistanceSolution ds = SU[N].first->straightLineDistanceEstimate(pos, dir, SU[N].second);
565 if (ds.currentDistance(false) > .010)
566 return nullptr;
567
568 P[0] = Rs[0];
569 A[0] = As[0];
570 P[1] = Rs[1];
571 A[1] = As[1];
572 P[2] = Rs[2];
573 A[2] = As[2];
574
575 So.push_back(N);
576
577 // Transformation track parameters
578 //
579 bool useJac = 0;
580 Tp.covariance() ? useJac = true : useJac = false;
581
582 if (useJac) {
583 const double d = 1. / P[6];
584 P[35] *= d;
585 P[36] *= d;
586 P[37] *= d;
587 P[38] *= d;
588 P[39] *= d;
589 P[40] *= d;
590 }
591 double p[5];
592 double Jac[25];
593 Trk::RungeKuttaUtils::transformGlobalToLocal(SU[N].first, useJac, P, p, Jac);
594
595 if (!useJac){
596 return SU[N].first->createUniqueTrackParameters(p[0], p[1], p[2], p[3], p[4], std::nullopt);
597 }
598
599 AmgSymMatrix(5) e = Trk::RungeKuttaUtils::newCovarianceMatrix(Jac, *Tp.covariance());
600 if (!Amg::hasPositiveOrZeroDiagElems(e)) {
601 return nullptr;
602 }
603 return SU[N].first->createUniqueTrackParameters(p[0], p[1], p[2], p[3], p[4],
604 std::move(e));
605}
606
608// Step estimator take into accout curvature
610double
611stepEstimatorWithCurvature(Cache& cache,
612 int kind,
613 double* ATH_RESTRICT Su,
614 const double* ATH_RESTRICT P,
615 bool& Q)
616{
617 // Straight step estimation
618 const double Step = Trk::RungeKuttaUtils::stepEstimator(kind, Su, P, Q);
619 if (!Q)
620 return 0.;
621 const double AStep = std::abs(Step);
622 if (kind || AStep < cache.m_straightStep || !cache.m_mcondition)
623 return Step;
624
625 const double* SA = &P[42]; // Start direction
626 const double S = .5 * Step;
627
628 const double Ax = P[3] + S * SA[0];
629 const double Ay = P[4] + S * SA[1];
630 const double Az = P[5] + S * SA[2];
631 const double As = 1. / std::sqrt(Ax * Ax + Ay * Ay + Az * Az);
632
633 const double PN[6] = { P[0], P[1], P[2], Ax * As, Ay * As, Az * As };
634 const double StepN = Trk::RungeKuttaUtils::stepEstimator(kind, Su, PN, Q);
635 if (!Q) {
636 Q = true;
637 return Step;
638 }
639 if (std::abs(StepN) < AStep)
640 return StepN;
641 return Step;
642}
643
645// Runge Kutta main program for propagation with or without Jacobian
647bool
648propagateWithJacobian(Cache& cache,
649 bool Jac,
650 int kind,
651 double* ATH_RESTRICT Su,
652 double* ATH_RESTRICT P,
653 double& ATH_RESTRICT W)
654{
655 const double Smax = 1000.; // max. step allowed
656 const double Wmax = cache.m_maxPath; // Max way allowed
657 const double Wwrong = 500.; // Max way with wrong direction
658 double* R = &P[0]; // Start coordinates
659 double* A = &P[3]; // Start directions
660 double* SA = &P[42];
661 SA[0] = SA[1] = SA[2] = 0.;
662 cache.m_maxPathLimit = false;
663
664 if (cache.m_mcondition && std::abs(P[6]) > .1)
665 return false;
666
667 // Step estimation until surface
668 //
669 bool Q = 0;
670 double S = 0;
671 double Step = Trk::RungeKuttaUtils::stepEstimator(kind, Su, P, Q);
672 if (!Q)
673 return false;
674
675 bool dir = true;
676 if (cache.m_mcondition && cache.m_direction && cache.m_direction * Step < 0.) {
677 Step = -Step;
678 dir = false;
679 }
680
681 Step > Smax ? S = Smax : Step < -Smax ? S = -Smax : S = Step;
682 double So = std::abs(S);
683 int iS = 0;
684
685 bool InS = false;
686
687 // Rkuta extrapolation
688 //
689 int niter = 0;
690 cache.m_newfield = true;
691 while (std::abs(Step) > cache.m_straightStep) {
692
693 if (++niter > 10000)
694 return false;
695
696 if (cache.m_mcondition) {
697
698 if (!cache.m_needgradient)
699 W += (S = rungeKuttaStep(cache, Jac, S, P, InS));
700 else
701 W += (S = rungeKuttaStepWithGradient(cache, S, P, InS));
702 } else {
703 W += (S = straightLineStep(Jac, S, P));
704 }
705
706 Step = stepEstimatorWithCurvature(cache, kind, Su, P, Q);
707 if (!Q)
708 return false;
709
710 if (!dir) {
711 if (cache.m_direction && cache.m_direction * Step < 0.)
712 Step = -Step;
713 else
714 dir = true;
715 }
716
717 if (S * Step < 0.) {
718 S = -S;
719 ++iS;
720 }
721
722 const double aS = std::abs(S);
723 const double aStep = std::abs(Step);
724 if (aS > aStep)
725 S = Step;
726 else if (!iS && InS && aS * 2. < aStep)
727 S *= 2.;
728 if (!dir && std::abs(W) > Wwrong)
729 return false;
730
731 if (iS > 10 || (iS > 3 && std::abs(S) >= So)) {
732 if (!kind)
733 break;
734 return false;
735 }
736 const double dW = Wmax - std::abs(W);
737 if (std::abs(S) > dW) {
738 S > 0. ? S = dW : S = -dW;
739 Step = S;
740 cache.m_maxPathLimit = true;
741 }
742 So = std::abs(S);
743 }
744
745 // Output track parameteres
746 //
747 W += Step;
748
749 if (std::abs(Step) < .001)
750 return true;
751
752 A[0] += (SA[0] * Step);
753 A[1] += (SA[1] * Step);
754 A[2] += (SA[2] * Step);
755 const double CBA = 1. / std::sqrt(A[0] * A[0] + A[1] * A[1] + A[2] * A[2]);
756
757 R[0] += Step * (A[0] - .5 * Step * SA[0]);
758 A[0] *= CBA;
759 R[1] += Step * (A[1] - .5 * Step * SA[1]);
760 A[1] *= CBA;
761 R[2] += Step * (A[2] - .5 * Step * SA[2]);
762 A[2] *= CBA;
763 return true;
764}
765
766bool propagateWithJacobianSwitch(Cache& cache,
767 const Trk::Surface& Su,
768 bool useJac,
769 double* ATH_RESTRICT P,
770 double& ATH_RESTRICT Step) {
771
772 const Amg::Transform3D& T = Su.transform();
773 const Trk::SurfaceType ty = Su.type();
774
775 switch (ty) {
776 case Trk::SurfaceType::Line:
777 case Trk::SurfaceType::Perigee: {
778 double s[6] = {T(0, 3), T(1, 3), T(2, 3), T(0, 2), T(1, 2), T(2, 2)};
779 return propagateWithJacobian(cache, useJac, 0, s, P, Step);
780 }
781 case Trk::SurfaceType::Plane:
782 case Trk::SurfaceType::Disc: {
783 double s[4];
784 const double d =
785 T(0, 3) * T(0, 2) + T(1, 3) * T(1, 2) + T(2, 3) * T(2, 2);
786
787 if (d >= 0.) {
788 s[0] = T(0, 2);
789 s[1] = T(1, 2);
790 s[2] = T(2, 2);
791 s[3] = d;
792 } else {
793 s[0] = -T(0, 2);
794 s[1] = -T(1, 2);
795 s[2] = -T(2, 2);
796 s[3] = -d;
797 }
798 return propagateWithJacobian(cache, useJac, 1, s, P, Step);
799 }
800 case Trk::SurfaceType::Cylinder: {
801 const Trk::CylinderSurface* cyl =
802 static_cast<const Trk::CylinderSurface*>(&Su);
803
804 const double r0[3] = {P[0], P[1], P[2]};
805 double s[9] = {T(0, 3), T(1, 3), T(2, 3),
806 T(0, 2), T(1, 2), T(2, 2),
807 cyl->bounds().r(), cache.m_direction, 0.};
808
809 const bool status = propagateWithJacobian(cache, useJac, 2, s, P, Step);
810 // For cylinder we do test for next cross point
811 if (status && cyl->bounds().halfPhiSector() < 3.1 &&
812 newCrossPoint(*cyl, r0, P)) {
813 s[8] = 0.;
814 return propagateWithJacobian(cache, useJac, 2, s, P, Step);
815 }
816 return status;
817 }
818 case Trk::SurfaceType::Cone: {
819 double k = static_cast<const Trk::ConeSurface*>(&Su)->bounds().tanAlpha();
820 k = k * k + 1.;
821 double s[9] = {T(0, 3), T(1, 3), T(2, 3), T(0, 2), T(1, 2),
822 T(2, 2), k, cache.m_direction, 0.};
823 return propagateWithJacobian(cache, useJac, 3, s, P, Step);
824 }
825 default: {
826 return false;
827 }
828 }
829}
830
832// Main function for neutral track parameters propagation with or without jacobian
834std::unique_ptr<Trk::NeutralParameters>
835propagateStraightLine(Cache& cache,
836 bool useJac,
837 const Trk::NeutralParameters& Tp,
838 const Trk::Surface& Su,
839 Trk::PropDirection D,
840 const Trk::BoundaryCheck& B,
841 double* ATH_RESTRICT Jac,
842 bool returnCurv)
843{
844 const Trk::Surface* su = &Su;
845 if (su == &Tp.associatedSurface()){
846 return buildTrackParametersWithoutPropagation(Tp, Jac);
847 }
848
849 cache.m_direction = D;
850 cache.m_mcondition = false;
851
852 double P[64];
853 double Step = 0;
854 if (!Trk::RungeKuttaUtils::transformLocalToGlobal(useJac, Tp, P)){
855 return nullptr;
856 }
857
858 if (!propagateWithJacobianSwitch(cache,Su,useJac,P,Step)){
859 return nullptr;
860 }
861
862 if (cache.m_direction != 0. && (cache.m_direction * Step) < 0.) {
863 return nullptr;
864 }
865
866 // Common transformation for all surfaces (angles and momentum)
867 //
868 if (useJac) {
869 const double p = 1. / P[6];
870 P[35] *= p;
871 P[36] *= p;
872 P[37] *= p;
873 P[38] *= p;
874 P[39] *= p;
875 P[40] *= p;
876 }
877
878 if (cache.m_maxPathLimit)
879 returnCurv = true;
880
881 bool uJ = useJac;
882 if (returnCurv)
883 uJ = false;
884
885 double p[5];
886 Trk::RungeKuttaUtils::transformGlobalToLocal(su, uJ, P, p, Jac);
887
888 if (B) {
889 const Amg::Vector2D L(p[0], p[1]);
890 if (!Su.insideBounds(L, 0.))
891 return nullptr;
892 }
893
894 // Transformation to curvilinear presentation
895 //
896 if (returnCurv) {
897 Trk::RungeKuttaUtils::transformGlobalToCurvilinear(useJac, P, p, Jac);
898 }
899
900 if (!useJac || !Tp.covariance()) {
901
902 if (!returnCurv) {
903 return Su.createUniqueNeutralParameters(p[0], p[1], p[2], p[3], p[4], std::nullopt);
904 } else {
905 const Amg::Vector3D gp(P[0], P[1], P[2]);
906 return std::make_unique<Trk::NeutralCurvilinearParameters>(gp, p[2], p[3], p[4]);
907 }
908 }
909
910 AmgSymMatrix(5) e = Trk::RungeKuttaUtils::newCovarianceMatrix(Jac, *Tp.covariance());
911
912 if (!Amg::hasPositiveOrZeroDiagElems(e)) {
913 return nullptr;
914 }
915
916 if (!returnCurv) {
917 return Su.createUniqueNeutralParameters(p[0], p[1], p[2], p[3], p[4], std::move(e));
918 } else {
919 const Amg::Vector3D gp(P[0], P[1], P[2]);
920 return std::make_unique<Trk::NeutralCurvilinearParameters>(gp, p[2], p[3], p[4], std::move(e));
921 }
922}
923
925// Global positions calculation inside CylinderBounds (one side)
926// where mS - max step allowed if mS > 0 propagate along momentum
927// mS < 0 propogate opposite momentum
929void
930globalOneSidePositions(Cache& cache,
931 std::deque<Amg::Vector3D>& GP,
932 const double* ATH_RESTRICT P,
933 const Trk::MagneticFieldProperties& M,
934 const Trk::CylinderBounds& CB,
935 double mS)
936{
937 M.magneticFieldMode() == Trk::FastField ? cache.m_solenoid = true : cache.m_solenoid = false;
938 M.magneticFieldMode() != Trk::NoField ? cache.m_mcondition = true : cache.m_mcondition = false;
939
940 double Pm[45];
941 for (int i = 0; i != 7; ++i)
942 Pm[i] = P[i];
943
944 double W = 0.; // way
945 const double R2max = CB.r() * CB.r(); // max. radius**2 of region
946 const double Zmax = CB.halflengthZ(); // max. Z of region
947 double R2 = P[0] * P[0] + P[1] * P[1]; // Start radius**2
948 const double Dir = P[0] * P[3] + P[1] * P[4]; // Direction
949 double S = mS; // max step allowed
950 double R2m = R2;
951
952 if (cache.m_mcondition && std::abs(P[6]) > .1)
953 return;
954
955 // Test position of the track
956 //
957 if (std::abs(P[2]) > Zmax || R2 > R2max)
958 return;
959
960 const Amg::Vector3D g0(P[0], P[1], P[2]);
961 GP.push_back(g0);
962
963 bool per = false;
964 if (std::abs(Dir) < .00001)
965 per = true;
966 bool InS = false;
967
968 int niter = 0;
969 int sm = 1;
970 for (int s = 0; s != 2; ++s) {
971
972 cache.m_newfield = true;
973
974 if (s) {
975 if (per)
976 break;
977 S = -S;
978 }
979 double p[45];
980 for (int i = 0; i != 7; ++i)
981 p[i] = P[i];
982 p[42] = p[43] = p[44] = 0.;
983
984 while (W < 100000. && ++niter < 1000) {
985
986 if (cache.m_mcondition) {
987 W += (S = rungeKuttaStep(cache, 0, S, p, InS));
988 } else {
989 W += (S = straightLineStep(0, S, p));
990 }
991 if (InS && std::abs(2. * S) < mS)
992 S *= 2.;
993
994 const Amg::Vector3D g(p[0], p[1], p[2]);
995 if (!s)
996 GP.push_back(g);
997 else
998 GP.push_front(g);
999
1000 // Test position of the track
1001 //
1002 R2 = p[0] * p[0] + p[1] * p[1];
1003
1004 if (R2 < R2m) {
1005 Pm[0] = p[0];
1006 Pm[1] = p[1];
1007 Pm[2] = p[2];
1008 Pm[3] = p[3];
1009 Pm[4] = p[4];
1010 Pm[5] = p[5];
1011 R2m = R2;
1012 sm = s;
1013 }
1014 if (std::abs(p[2]) > Zmax || R2 > R2max)
1015 break;
1016 if (!s && P[3] * p[3] + P[4] * p[4] < 0.)
1017 break;
1018
1019 // Test perigee
1020 //
1021 if ((p[0] * p[3] + p[1] * p[4]) * Dir < 0.) {
1022 if (s)
1023 break;
1024 per = true;
1025 }
1026 }
1027 }
1028
1029 if (R2m < 400.)
1030 return;
1031
1032 // Propagate to perigee
1033 //
1034 Pm[42] = Pm[43] = Pm[44] = 0.;
1035 per = false;
1036
1037 for (int s = 0; s != 3; ++s) {
1038
1039 cache.m_newfield = true;
1040
1041 const double A = (1. - Pm[5]) * (1. + Pm[5]);
1042 if (A == 0.)
1043 break;
1044 S = -(Pm[0] * Pm[3] + Pm[1] * Pm[4]) / A;
1045 if (std::abs(S) < 1. || ++niter > 1000)
1046 break;
1047
1048 if (cache.m_mcondition) {
1049 W += rungeKuttaStep(cache, 0, S, Pm, InS);
1050 } else {
1051 W += straightLineStep(0, S, Pm);
1052 }
1053 per = true;
1054 }
1055
1056 if (per) {
1057 if (sm) {
1058 const Amg::Vector3D gf(Pm[0], Pm[1], Pm[2]);
1059 GP.front() = gf;
1060 } else {
1061 const Amg::Vector3D gf(Pm[0], Pm[1], Pm[2]);
1062 GP.back() = gf;
1063 }
1064 } else {
1065 const double x = GP.front().x();
1066 const double y = GP.front().y();
1067 if ((x * x + y * y) > (Pm[0] * Pm[0] + Pm[1] * Pm[1])) {
1068 if (sm)
1069 GP.pop_front();
1070 else
1071 GP.pop_back();
1072 }
1073 }
1074}
1075
1077// Global positions calculation inside CylinderBounds (one side)
1078// where mS - max step allowed if mS > 0 propagate along momentum
1079// mS < 0 propogate opposite momentum
1081void
1082globalTwoSidePositions(Cache& cache,
1083 std::deque<Amg::Vector3D>& GP,
1084 const double* ATH_RESTRICT P,
1085 const Trk::MagneticFieldProperties& M,
1086 const Trk::CylinderBounds& CB,
1087 double mS)
1088{
1089 M.magneticFieldMode() == Trk::FastField ? cache.m_solenoid = true : cache.m_solenoid = false;
1090 M.magneticFieldMode() != Trk::NoField ? cache.m_mcondition = true : cache.m_mcondition = false;
1091
1092 double W = 0.; // way
1093 const double R2max = CB.r() * CB.r(); // max. radius**2 of region
1094 const double Zmax = CB.halflengthZ(); // max. Z of region
1095 double R2 = P[0] * P[0] + P[1] * P[1]; // Start radius**2
1096 double S = mS; // max step allowed
1097
1098 if (cache.m_mcondition && std::abs(P[6]) > .1)
1099 return;
1100
1101 // Test position of the track
1102 //
1103 if (std::abs(P[2]) > Zmax || R2 > R2max)
1104 return;
1105
1106 const Amg::Vector3D g0(P[0], P[1], P[2]);
1107 GP.push_back(g0);
1108
1109 bool InS = false;
1110
1111 int niter = 0;
1112 for (int s = 0; s != 2; ++s) {
1113
1114 cache.m_newfield = true;
1115
1116 if (s) {
1117 S = -S;
1118 }
1119 double p[45];
1120 for (int i = 0; i != 7; ++i)
1121 p[i] = P[i];
1122 p[42] = p[43] = p[44] = 0.;
1123
1124 while (W < 100000. && ++niter < 1000) {
1125
1126 if (cache.m_mcondition) {
1127 W += (S = rungeKuttaStep(cache, 0, S, p, InS));
1128 } else {
1129 W += (S = straightLineStep(0, S, p));
1130 }
1131 if (InS && std::abs(2. * S) < mS)
1132 S *= 2.;
1133
1134 const Amg::Vector3D g(p[0], p[1], p[2]);
1135 if (!s)
1136 GP.push_back(g);
1137 else
1138 GP.push_front(g);
1139
1140 // Test position of the track
1141 //
1142 R2 = p[0] * p[0] + p[1] * p[1];
1143 if (std::abs(p[2]) > Zmax || R2 > R2max)
1144 break;
1145 }
1146 }
1147}
1149// Main function for charged track parameters propagation with or without
1150// jacobian
1152std::unique_ptr<Trk::TrackParameters>
1153propagateRungeKutta(Cache& cache,
1154 bool useJac,
1155 const Trk::TrackParameters& Tp,
1156 const Trk::Surface& Su,
1157 Trk::PropDirection D,
1158 const Trk::BoundaryCheck& B,
1159 const Trk::MagneticFieldProperties& M,
1160 double* Jac,
1161 bool returnCurv)
1162{
1163 const Trk::Surface* su = &Su;
1164
1165 cache.m_direction = D;
1166
1167 M.magneticFieldMode() == Trk::FastField ? cache.m_solenoid = true
1168 : cache.m_solenoid = false;
1169 (useJac && cache.m_usegradient) ? cache.m_needgradient = true
1170 : cache.m_needgradient = false;
1171 M.magneticFieldMode() != Trk::NoField ? cache.m_mcondition = true
1172 : cache.m_mcondition = false;
1173
1174 if (su == &Tp.associatedSurface())
1175 return buildTrackParametersWithoutPropagation(Tp, Jac);
1176
1177 double P[64];
1178 double Step = 0.;
1179 if (!Trk::RungeKuttaUtils::transformLocalToGlobal(useJac, Tp, P)){
1180 return nullptr;
1181 }
1182
1183 if (!propagateWithJacobianSwitch(cache,Su,useJac,P,Step)){
1184 return nullptr;
1185 }
1186
1187 if (cache.m_direction && (cache.m_direction * Step) < 0.) {
1188 return nullptr;
1189 }
1190 cache.m_step = Step;
1191
1192 // Common transformation for all surfaces (angles and momentum)
1193 //
1194 if (useJac) {
1195 const double p = 1. / P[6];
1196 P[35] *= p;
1197 P[36] *= p;
1198 P[37] *= p;
1199 P[38] *= p;
1200 P[39] *= p;
1201 P[40] *= p;
1202 }
1203
1204 if (cache.m_maxPathLimit)
1205 returnCurv = true;
1206
1207 bool uJ = useJac;
1208 if (returnCurv)
1209 uJ = false;
1210 double p[5];
1211 Trk::RungeKuttaUtils::transformGlobalToLocal(su, uJ, P, p, Jac);
1212
1213 if (B) {
1214 Amg::Vector2D const L(p[0], p[1]);
1215 if (!Su.insideBounds(L, 0.))
1216 return nullptr;
1217 }
1218
1219 // Transformation to curvilinear presentation
1220 //
1221 if (returnCurv)
1222 Trk::RungeKuttaUtils::transformGlobalToCurvilinear(useJac, P, p, Jac);
1223
1224 if (!useJac || !Tp.covariance()) {
1225
1226 if (!returnCurv) {
1227 return Su.createUniqueTrackParameters(
1228 p[0], p[1], p[2], p[3], p[4], std::nullopt);
1229 } else {
1230 const Amg::Vector3D gp(P[0], P[1], P[2]);
1231 return std::make_unique<Trk::CurvilinearParameters>(gp, p[2], p[3], p[4]);
1232 }
1233 }
1234
1235 AmgSymMatrix(5) e =
1236 Trk::RungeKuttaUtils::newCovarianceMatrix(Jac, *Tp.covariance());
1237
1238 if (!Amg::hasPositiveOrZeroDiagElems(e)) {
1239 return nullptr;
1240 }
1241
1242 if (!returnCurv) {
1243 return Su.createUniqueTrackParameters(
1244 p[0], p[1], p[2], p[3], p[4], std::move(e));
1245 } else {
1246 const Amg::Vector3D gp(P[0], P[1], P[2]);
1247 return std::make_unique<Trk::CurvilinearParameters>(
1248 gp, p[2], p[3], p[4], std::move(e));
1249 }
1250}
1251
1253// Main function for simple track propagation with or without jacobian
1254// Ta->Su = Tb for pattern track parameters
1256bool
1257propagateRungeKutta(Cache& cache,
1258 bool useJac,
1259 Trk::PatternTrackParameters& Ta,
1260 const Trk::Surface& Su,
1261 Trk::PatternTrackParameters& Tb,
1262 Trk::PropDirection D,
1263 const Trk::MagneticFieldProperties& M,
1264 double& Step)
1265{
1266 const Trk::Surface* su = &Su;
1267 if (!su)
1268 return false;
1269 if (su == &Ta.associatedSurface()) {
1270 Tb = Ta;
1271 return true;
1272 }
1273
1274 cache.m_direction = D;
1275
1276 if (useJac && !Ta.iscovariance())
1277 useJac = false;
1278
1279 M.magneticFieldMode() == Trk::FastField ? cache.m_solenoid = true
1280 : cache.m_solenoid = false;
1281 (useJac && cache.m_usegradient) ? cache.m_needgradient = true
1282 : cache.m_needgradient = false;
1283 M.magneticFieldMode() != Trk::NoField ? cache.m_mcondition = true
1284 : cache.m_mcondition = false;
1285
1286 double P[45];
1287 if (!Trk::RungeKuttaUtils::transformLocalToGlobal(useJac, Ta, P)){
1288 return false;
1289 }
1290 Step = 0.;
1291
1292 if (!propagateWithJacobianSwitch(cache,Su,useJac,P,Step)){
1293 return false;
1294 }
1295
1296 if (cache.m_maxPathLimit ||
1297 (cache.m_direction && (cache.m_direction * Step) < 0.)){
1298 return false;
1299 }
1300
1301 // Common transformation for all surfaces (angles and momentum)
1302 //
1303 if (useJac) {
1304 const double p = 1. / P[6];
1305 P[35] *= p;
1306 P[36] *= p;
1307 P[37] *= p;
1308 P[38] *= p;
1309 P[39] *= p;
1310 P[40] *= p;
1311 }
1312
1313 double p[5];
1314 double Jac[21];
1315 Trk::RungeKuttaUtils::transformGlobalToLocal(su, useJac, P, p, Jac);
1316
1317 // New simple track parameters production
1318 //
1319 if (useJac) {
1320 const AmgSymMatrix(5) newCov =
1321 Trk::RungeKuttaUtils::newCovarianceMatrix(Jac, *Ta.covariance());
1322 Tb.setParametersWithCovariance(&Su, p, newCov);
1323 const AmgSymMatrix(5)& cv = *Tb.covariance();
1324 if (!Amg::hasPositiveOrZeroDiagElems(cv))
1325 return false;
1326 } else {
1327 Tb.setParameters(&Su, p);
1328 }
1329 return true;
1330}
1331
1333// GlobalPostions and steps for set surfaces
1335void
1336globalPositionsImpl(
1337 Cache& cache,
1338 const Trk::PatternTrackParameters& Tp,
1339 std::vector<const Trk::Surface*>& SU,
1340 std::vector<std::pair<Amg::Vector3D, double>>& GP,
1341 const Trk::MagneticFieldProperties& M)
1342{
1343 cache.m_direction = 0.;
1344 cache.m_mcondition = false;
1345 cache.m_maxPath = 10000.;
1346 cache.m_needgradient = false;
1347 M.magneticFieldMode() == Trk::FastField ? cache.m_solenoid = true : cache.m_solenoid = false;
1348 M.magneticFieldMode() != Trk::NoField ? cache.m_mcondition = true : cache.m_mcondition = false;
1349
1350 double Step = 0.;
1351 double P[64];
1352 if (!Trk::RungeKuttaUtils::transformLocalToGlobal(false, Tp, P)){
1353 return;
1354 }
1355
1356 auto su = SU.begin();
1357 auto sue = SU.end();
1358 // Loop trough all input surfaces
1359 for (; su != sue; ++su) {
1360
1361 if (!propagateWithJacobianSwitch(cache, **su, false, P, Step)) {
1362 return;
1363 }
1364
1365 if (cache.m_maxPathLimit) {
1366 return;
1367 }
1368
1369 const Amg::Vector3D gp(P[0], P[1], P[2]);
1370 GP.emplace_back(gp, Step);
1371 }
1372}
1373
1374/*
1375 * end of anonymous namespace
1376 * with internal implementation methods
1377 */
1378}
1379
1380Trk::RungeKuttaPropagator::RungeKuttaPropagator(const std::string& p,
1381 const std::string& n,
1382 const IInterface* t)
1383 : AthAlgTool(p, n, t)
1384{
1385 declareInterface<Trk::IPropagator>(this);
1386 declareInterface<Trk::IPatternParametersPropagator>(this);
1387}
1388
1389StatusCode
1390Trk::RungeKuttaPropagator::initialize()
1391{
1392 // Read handle for AtlasFieldCacheCondObj
1393 ATH_CHECK(m_fieldCondObjInputKey.initialize());
1394 ATH_MSG_DEBUG("initialize() init key: " << m_fieldCondObjInputKey.key());
1395 return StatusCode::SUCCESS;
1396}
1397
1399// Main function for NeutralParameters propagation
1401std::unique_ptr<Trk::NeutralParameters>
1402Trk::RungeKuttaPropagator::propagate(const Trk::NeutralParameters& Tp,
1403 const Trk::Surface& Su,
1404 Trk::PropDirection D,
1405 const Trk::BoundaryCheck& B,
1406 bool returnCurv) const
1407{
1408 double J[25];
1409 Cache cache{};
1410 cache.m_dlt = m_dlt;
1411 cache.m_helixStep = m_helixStep;
1412 cache.m_straightStep = m_straightStep;
1413 cache.m_usegradient = m_usegradient;
1414 return propagateStraightLine(cache, true, Tp, Su, D, B, J, returnCurv);
1415}
1416
1418// Main function for track parameters and covariance matrix propagation
1419// without transport Jacobian production
1421std::unique_ptr<Trk::TrackParameters>
1422Trk::RungeKuttaPropagator::propagate(const ::EventContext& ctx,
1423 const Trk::TrackParameters& Tp,
1424 const Trk::Surface& Su,
1425 Trk::PropDirection D,
1426 const Trk::BoundaryCheck& B,
1427 const MagneticFieldProperties& M,
1428 ParticleHypothesis,
1429 bool returnCurv,
1430 const TrackingVolume*) const
1431{
1432 double J[25];
1433 Cache cache = getInitializedCache(ctx);
1434 return propagateRungeKutta(cache, true, Tp, Su, D, B, M, J, returnCurv);
1435}
1436
1438// Main function for MultiComponentState propagation used by the GSF
1440Trk::MultiComponentState
1441Trk::RungeKuttaPropagator::multiStatePropagate(
1442 const ::EventContext& ctx,
1443 const Trk::MultiComponentState& multiComponentState,
1444 const Trk::Surface& surface,
1445 const Trk::MagneticFieldProperties& fieldProperties,
1446 const Trk::PropDirection direction,
1447 const Trk::BoundaryCheck& boundaryCheck,
1448 const Trk::ParticleHypothesis) const
1449{
1450 Cache cache = getInitializedCache(ctx);
1451
1452 Trk::MultiComponentState propagatedState{};
1453 propagatedState.reserve(multiComponentState.size());
1454 double sumw(0); // sum of the weights of the propagated parameters
1455 double J[25];
1456 Trk::MultiComponentState::const_iterator component =
1457 multiComponentState.begin();
1458 for (; component != multiComponentState.end(); ++component) {
1459 const Trk::TrackParameters* currentParameters = component->params.get();
1460 if (!currentParameters) {
1461 continue;
1462 }
1463 auto propagatedParameters =
1464 propagateRungeKutta(cache, true, *currentParameters, surface, direction,
1465 boundaryCheck, fieldProperties, J, false);
1466
1467 if (!propagatedParameters) {
1468 continue;
1469 }
1470 sumw += component->weight;
1471 // Propagation does not affect the weightings of the states
1472 propagatedState.push_back({std::move(propagatedParameters),
1473 component->weight});
1474 }
1475 // Protect low weight propagation
1476 constexpr double minPropWeight = (1./12.);
1477 if (sumw < minPropWeight) {
1478 propagatedState.clear();
1479 }
1480 return propagatedState;
1481}
1482
1484// Main function for track parameters and covariance matrix propagation
1485// with transport Jacobian production
1487std::unique_ptr<Trk::TrackParameters>
1488Trk::RungeKuttaPropagator::propagate(const ::EventContext& ctx,
1489 const Trk::TrackParameters& Tp,
1490 const Trk::Surface& Su,
1491 Trk::PropDirection D,
1492 const Trk::BoundaryCheck& B,
1493 const MagneticFieldProperties& M,
1494 std::optional<TransportJacobian>& Jac,
1495 double& pathLength,
1496 ParticleHypothesis,
1497 bool returnCurv,
1498 const TrackingVolume*) const
1499{
1500 double J[25];
1501 Cache cache = getInitializedCache(ctx);
1502 pathLength < 0. ? cache.m_maxPath = 10000. : cache.m_maxPath = pathLength;
1503 auto Tpn = propagateRungeKutta(cache, true, Tp, Su, D, B, M, J, returnCurv);
1504 pathLength = cache.m_step;
1505
1506 if (Tpn) {
1507 J[24] = J[20];
1508 J[23] = 0.;
1509 J[22] = 0.;
1510 J[21] = 0.;
1511 J[20] = 0.;
1512 Jac = std::make_optional<Trk::TransportJacobian>(J);
1513 } else
1514 Jac.reset();
1515 return Tpn;
1516}
1517
1519// Main function to finds the closest surface
1521std::unique_ptr<Trk::TrackParameters>
1522Trk::RungeKuttaPropagator::propagate(const ::EventContext& ctx,
1523 const TrackParameters& Tp,
1524 std::vector<DestSurf>& DS,
1525 PropDirection D,
1526 const MagneticFieldProperties& M,
1527 ParticleHypothesis,
1528 std::vector<unsigned int>& Sol,
1529 double& Path,
1530 bool usePathLim,
1531 bool,
1532 const TrackingVolume*) const
1533{
1534
1535 Cache cache = getInitializedCache(ctx);
1536 Sol.erase(Sol.begin(), Sol.end());
1537 Path = 0.;
1538 if (DS.empty())
1539 return nullptr;
1540 cache.m_direction = D;
1541
1542 // Test is it measured track parameters
1543 //
1544 bool useJac = 0;
1545 Tp.covariance() ? useJac = true : useJac = false;
1546
1547 // Magnetic field information preparation
1548 //
1549 M.magneticFieldMode() == Trk::FastField ? cache.m_solenoid = true : cache.m_solenoid = false;
1550 (useJac && m_usegradient) ? cache.m_needgradient = true : cache.m_needgradient = false;
1551 M.magneticFieldMode() != Trk::NoField ? cache.m_mcondition = true : cache.m_mcondition = false;
1552
1553 // Transform to global presentation
1554 //
1555
1556 double Po[45];
1557 double Pn[45];
1558 if (!Trk::RungeKuttaUtils::transformLocalToGlobal(useJac, Tp, Po))
1559 return nullptr;
1560 Po[42] = Po[43] = Po[44] = 0.;
1561
1562 // Straight line track propagation for small step
1563 //
1564 if (D != 0) {
1565 double S = cache.m_straightStep;
1566 if (D < 0)
1567 S = -S;
1568 S = straightLineStep(useJac, S, Po);
1569 }
1570 double Wmax = 50000.; // Max pass
1571 double W = 0.; // Current pass
1572 double Smax = 100.; // Max step
1573 if (D < 0)
1574 Smax = -Smax;
1575 if (usePathLim)
1576 Wmax = std::abs(Path);
1577
1578 std::multimap<double, int> DN;
1579 double Scut[3];
1580 int Nveto = Trk::RungeKuttaUtils::fillDistancesMap(DS, DN, Po, W, &Tp.associatedSurface(), Scut);
1581
1582 // Test conditions tor start propagation and chocse direction if D == 0
1583 //
1584 if (DN.empty())
1585 return nullptr;
1586
1587 if (D == 0 && std::abs(Scut[0]) < std::abs(Scut[1]))
1588 Smax = -Smax;
1589
1590 if (Smax < 0. && Scut[0] > Smax)
1591 Smax = Scut[0];
1592 if (Smax > 0. && Scut[1] < Smax)
1593 Smax = Scut[1];
1594 if (Wmax > 3. * Scut[2])
1595 Wmax = 3. * Scut[2];
1596
1597 double Sl = Smax;
1598 double St = Smax;
1599 bool InS = false;
1600
1601 for (int i = 0; i != 45; ++i)
1602 Pn[i] = Po[i];
1603
1604 //----------------------------------Niels van Eldik patch
1605 double last_St = 0.;
1606 bool last_InS = !InS;
1607 bool reverted_P = false;
1608 //----------------------------------
1609
1610 cache.m_newfield = true;
1619 constexpr unsigned int max_back_forth_flips{ 100 };
1620 unsigned int flips{0};
1621 int last_surface{-1};
1622
1623 while (std::abs(W) < Wmax) {
1624
1625 std::pair<double, int> SN;
1626 double S = 0;
1627
1628 if (cache.m_mcondition) {
1629
1630 //----------------------------------Niels van Eldik patch
1631 if (reverted_P && std::abs(St - last_St) <= DBL_EPSILON &&
1632 InS == last_InS) {
1633 // inputs are not changed will get same result.
1634 break;
1635 }
1636 last_St = St;
1637 last_InS = InS;
1638 //----------------------------------
1639
1640 if (!cache.m_needgradient)
1641 S = rungeKuttaStep(cache, useJac, St, Pn, InS);
1642 else
1643 S = rungeKuttaStepWithGradient(cache, St, Pn, InS);
1644 } else {
1645
1646 //----------------------------------Niels van Eldik patch
1647 if (reverted_P && std::abs(St - last_St) <= DBL_EPSILON) {
1648 // inputs are not changed will get same result.
1649 break;
1650 }
1651 last_St = St;
1652 last_InS = InS;
1653 //----------------------------------
1654
1655 S = straightLineStep(useJac, St, Pn);
1656 }
1657 //----------------------------------Niels van Eldik patch
1658 reverted_P = false;
1659 //----------------------------------
1660
1661 bool next{ false };
1662 SN = Trk::RungeKuttaUtils::stepEstimator(DS, DN, Po, Pn, W, cache.m_straightStep, Nveto, next);
1663 if (next) {
1664 for (int i = 0; i != 45; ++i)
1665 Po[i] = Pn[i];
1666 W += S;
1667 Nveto = -1;
1668 } else {
1669 for (int i = 0; i != 45; ++i)
1670 Pn[i] = Po[i];
1671 reverted_P = true;
1672 cache.m_newfield = true;
1673 }
1674
1675 if (std::abs(S) + 1. < std::abs(St))
1676 Sl = S;
1677 InS ? St = 2. * S : St = S;
1678
1679 if (SN.second >= 0) {
1680
1681 const double Sa = std::abs(SN.first);
1682
1686 flips += last_surface == SN.second ? std::abs(last_St + SN.first) < 1.e-6
1687 : -flips;
1688 last_surface = SN.second;
1689
1690
1691 if (Sa > cache.m_straightStep) {
1692 if (std::abs(St) > Sa)
1693 St = SN.first;
1694 } else {
1695 Path = W + SN.first;
1696 if (auto To{ crossPoint(Tp, DS, Sol, Pn, SN) }; To){
1697 return To;
1698 }
1699 Nveto = SN.second;
1700 St = Sl;
1701 }
1702 } else if (std::abs(S) < DBL_EPSILON){
1703 return nullptr;
1704 }
1705
1706 if (flips > max_back_forth_flips) {
1707 return nullptr;
1708 }
1709 }
1710 return nullptr;
1711}
1712
1714// Main function for track parameters propagation without covariance matrix
1715// without transport Jacobian production
1717std::unique_ptr<Trk::TrackParameters>
1718Trk::RungeKuttaPropagator::propagateParameters(const ::EventContext& ctx,
1719 const Trk::TrackParameters& Tp,
1720 const Trk::Surface& Su,
1721 Trk::PropDirection D,
1722 const Trk::BoundaryCheck& B,
1723 const MagneticFieldProperties& M,
1724 ParticleHypothesis,
1725 bool returnCurv,
1726 const TrackingVolume*) const
1727{
1728 double J[25];
1729 Cache cache = getInitializedCache(ctx);
1730 return propagateRungeKutta(cache, false, Tp, Su, D, B, M, J, returnCurv);
1731}
1732
1734// Main function for track parameters propagation without covariance matrix
1735// with transport Jacobian production
1737std::unique_ptr<Trk::TrackParameters>
1738Trk::RungeKuttaPropagator::propagateParameters(const ::EventContext& ctx,
1739 const Trk::TrackParameters& Tp,
1740 const Trk::Surface& Su,
1741 Trk::PropDirection D,
1742 const Trk::BoundaryCheck& B,
1743 const MagneticFieldProperties& M,
1744 std::optional<TransportJacobian>& Jac,
1745 ParticleHypothesis,
1746 bool returnCurv,
1747 const TrackingVolume*) const
1748{
1749 double J[25];
1750 Cache cache = getInitializedCache(ctx);
1751 auto Tpn{ propagateRungeKutta(cache, true, Tp, Su, D, B, M, J, returnCurv) };
1752
1753 if (Tpn) {
1754 J[24] = J[20];
1755 J[23] = 0.;
1756 J[22] = 0.;
1757 J[21] = 0.;
1758 J[20] = 0.;
1759 Jac = std::make_optional<Trk::TransportJacobian>(J);
1760 } else
1761 Jac.reset();
1762 return Tpn;
1763}
1764
1766// Global positions calculation inside CylinderBounds
1767// where mS - max step allowed if mS > 0 propagate along momentum
1768// mS < 0 propogate opposite momentum
1770void
1771Trk::RungeKuttaPropagator::globalPositions(const ::EventContext& ctx,
1772 std::deque<Amg::Vector3D>& GP,
1773 const TrackParameters& Tp,
1774 const MagneticFieldProperties& M,
1775 const CylinderBounds& CB,
1776 double mS,
1777 ParticleHypothesis,
1778 const TrackingVolume*) const
1779{
1780 double P[45];
1781 if (!Trk::RungeKuttaUtils::transformLocalToGlobal(false, Tp, P))
1782 return;
1783 Cache cache = getInitializedCache(ctx);
1784
1785 cache.m_direction = std::abs(mS);
1786 if (mS > 0.)
1787 globalOneSidePositions(cache, GP, P, M, CB, mS);
1788 else
1789 globalTwoSidePositions(cache, GP, P, M, CB, -mS);
1790}
1791
1793// Global position together with direction of the trajectory on the surface
1795std::optional<Trk::TrackSurfaceIntersection>
1796Trk::RungeKuttaPropagator::intersect(const ::EventContext& ctx,
1797 const Trk::TrackParameters& Tp,
1798 const Trk::Surface& Su,
1799 const MagneticFieldProperties& M,
1800 ParticleHypothesis,
1801 const TrackingVolume*) const
1802{
1803 bool const nJ = false;
1804 const Trk::Surface* su = &Su;
1805 Cache cache = getInitializedCache(ctx);
1806 cache.m_direction = 0.;
1807 cache.m_needgradient = false;
1808
1809 M.magneticFieldMode() == Trk::FastField ? cache.m_solenoid = true : cache.m_solenoid = false;
1810 M.magneticFieldMode() != Trk::NoField ? cache.m_mcondition = true : cache.m_mcondition = false;
1811
1812 double P[64];
1813 if (!Trk::RungeKuttaUtils::transformLocalToGlobal(false, Tp, P)){
1814 return std::nullopt;
1815 }
1816 double Step = 0.;
1817 if (!propagateWithJacobianSwitch(cache, (*su), nJ, P, Step)) {
1818 return std::nullopt;
1819 }
1820
1821 const Amg::Vector3D Glo(P[0], P[1], P[2]);
1822 const Amg::Vector3D Dir(P[3], P[4], P[5]);
1823 return std::make_optional<Trk::TrackSurfaceIntersection>(Glo, Dir, Step);
1824}
1825
1827// Main function for simple track parameters and covariance matrix propagation
1828// Ta->Su = Tb
1830bool
1831Trk::RungeKuttaPropagator::propagate(const ::EventContext& ctx,
1832 Trk::PatternTrackParameters& Ta,
1833 const Trk::Surface& Su,
1834 Trk::PatternTrackParameters& Tb,
1835 Trk::PropDirection D,
1836 const MagneticFieldProperties& M,
1837 ParticleHypothesis) const
1838{
1839 double S = 0;
1840 Cache cache = getInitializedCache(ctx);
1841 return propagateRungeKutta(cache, true, Ta, Su, Tb, D, M, S);
1842}
1843
1845// Main function for simple track parameters and covariance matrix propagation
1846// Ta->Su = Tb with step to surface calculation
1848bool
1849Trk::RungeKuttaPropagator::propagate(const ::EventContext& ctx,
1850 Trk::PatternTrackParameters& Ta,
1851 const Trk::Surface& Su,
1852 Trk::PatternTrackParameters& Tb,
1853 Trk::PropDirection D,
1854 const MagneticFieldProperties& M,
1855 double& S,
1856 ParticleHypothesis) const
1857{
1858 Cache cache = getInitializedCache(ctx);
1859 return propagateRungeKutta(cache, true, Ta, Su, Tb, D, M, S);
1860}
1861
1863// Main function for simple track parameters propagation without covariance matrix
1864// Ta->Su = Tb
1866bool
1867Trk::RungeKuttaPropagator::propagateParameters(const ::EventContext& ctx,
1868 Trk::PatternTrackParameters& Ta,
1869 const Trk::Surface& Su,
1870 Trk::PatternTrackParameters& Tb,
1871 Trk::PropDirection D,
1872 const MagneticFieldProperties& M,
1873 ParticleHypothesis) const
1874{
1875 double S = 0;
1876 Cache cache = getInitializedCache(ctx);
1877 return propagateRungeKutta(cache, false, Ta, Su, Tb, D, M, S);
1878}
1879
1881// Main function for simple track parameters propagation without covariance matrix
1882// Ta->Su = Tb with step calculation
1884bool
1885Trk::RungeKuttaPropagator::propagateParameters(const ::EventContext& ctx,
1886 Trk::PatternTrackParameters& Ta,
1887 const Trk::Surface& Su,
1888 Trk::PatternTrackParameters& Tb,
1889 Trk::PropDirection D,
1890 const MagneticFieldProperties& M,
1891 double& S,
1892 ParticleHypothesis) const
1893{
1894 Cache cache = getInitializedCache(ctx);
1895 return propagateRungeKutta(cache, false, Ta, Su, Tb, D, M, S);
1896}
1897
1899// GlobalPostions and steps for set surfaces
1901void
1902Trk::RungeKuttaPropagator::globalPositions(const ::EventContext& ctx,
1903 const PatternTrackParameters& Tp,
1904 std::vector<const Trk::Surface*>& SU,
1905 std::vector<std::pair<Amg::Vector3D, double>>& GP,
1906 const Trk::MagneticFieldProperties& M,
1907 ParticleHypothesis) const
1908{
1909 Cache cache = getInitializedCache(ctx);
1910 globalPositionsImpl(cache, Tp, SU, GP, M);
1911}
1912
1913Cache
1914Trk::RungeKuttaPropagator::getInitializedCache(const EventContext& ctx) const
1915{
1916 SG::ReadCondHandle<AtlasFieldCacheCondObj> readHandle{m_fieldCondObjInputKey,
1917 ctx};
1918 const AtlasFieldCacheCondObj* fieldCondObj{*readHandle};
1919 Cache cache{};
1920 fieldCondObj->getInitializedCache(cache.m_fieldCache);
1921 cache.m_dlt = m_dlt;
1922 cache.m_helixStep = m_helixStep;
1923 cache.m_straightStep = m_straightStep;
1924 cache.m_usegradient = m_usegradient;
1925 return cache;
1926}
ATH_CHECK
#define ATH_CHECK
Evaluate an expression and check for errors.
Definition AthCheckMacros.h:40
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition AthMsgStreamMacros.h:29
B0
static const int B0
Definition AtlasPID.h:122
AmgSymMatrix
#define AmgSymMatrix(dim)
Definition EventPrimitives.h:50
Pi
static const double Pi
Definition LArCellBinning.cxx:17
Tp
static Double_t Tp(Double_t *t, Double_t *par)
Definition LArPhysWaveHECTool.cxx:354
Ro
static Double_t Ro(Double_t t)
Definition LArPhysWaveHECTool.cxx:298
P
static Double_t P(Double_t *tt, Double_t *par)
Definition LArPhysWaveHECTool.cxx:280
if
if(febId1==febId2)
Definition LArRodBlockPhysicsV0.cxx:567
H
#define H(x, y, z)
Definition MD5.cxx:114
H2
#define H2(x, y, z)
Definition MD5.cxx:115
S3
struct TBPatternUnitContext S3
pi
#define pi
Definition TileMuonFitter.cxx:65
y
#define y
x
#define x
z
#define z
AthAlgTool::AthAlgTool
AthAlgTool(const std::string &type, const std::string &name, const IInterface *parent)
Constructor with parameters:
Definition AthAlgTool.cxx:16
AtlasFieldCacheCondObj::getInitializedCache
void getInitializedCache(MagField::AtlasFieldCache &cache) const
get B field cache for evaluation as a function of 2-d or 3-d position.
Definition AtlasFieldCacheCondObj.h:32
Trk::BoundaryCheck
The BoundaryCheck class allows to steer the way surface boundaries are used for inside/outside checks...
Definition BoundaryCheck.h:51
Trk::ConeSurface
Class for a conical surface in the ATLAS detector.
Definition ConeSurface.h:51
Trk::CylinderBounds
Bounds for a cylindrical Surface.
Definition CylinderBounds.h:46
Trk::CylinderBounds::r
virtual double r() const override final
This method returns the radius.
Trk::CylinderBounds::halfPhiSector
double halfPhiSector() const
This method returns the halfPhiSector angle.
Trk::CylinderBounds::averagePhi
double averagePhi() const
This method returns the average phi.
Trk::CylinderSurface
Class for a CylinderSurface in the ATLAS detector.
Definition CylinderSurface.h:55
Trk::CylinderSurface::bounds
virtual const CylinderBounds & bounds() const override final
This method returns the CylinderBounds by reference (NoBounds is not possible for cylinder)
Trk::DistanceSolution
Access to distance solutions.
Definition DistanceSolution.h:25
Trk::MagneticFieldProperties
magnetic field properties to steer the behavior of the extrapolation
Definition MagneticFieldProperties.h:31
Trk::PatternTrackParameters::associatedSurface
virtual const Surface & associatedSurface() const override final
Access to the Surface associated to the Parameters.
Trk::PatternTrackParameters::setParameters
void setParameters(const Surface *, const double *)
Trk::RungeKuttaPropagator::m_fieldCondObjInputKey
SG::ReadCondHandleKey< AtlasFieldCacheCondObj > m_fieldCondObjInputKey
Definition RungeKuttaPropagator.h:387
Trk::RungeKuttaPropagator::initialize
virtual StatusCode initialize() override final
Definition RungeKuttaPropagator.cxx:1390
Trk::RungeKuttaPropagator::RungeKuttaPropagator
RungeKuttaPropagator(const std::string &, const std::string &, const IInterface *)
Definition RungeKuttaPropagator.cxx:1380
Trk::RungeKuttaPropagator::multiStatePropagate
virtual Trk::MultiComponentState multiStatePropagate(const EventContext &ctx, const MultiComponentState &multiComponentState, const Surface &surface, const MagneticFieldProperties &fieldProperties, const PropDirection direction=Trk::anyDirection, const BoundaryCheck &boundaryCheck=true, const ParticleHypothesis particleHypothesis=nonInteracting) const override final
Main propagation method for Multi Component state.
Definition RungeKuttaPropagator.cxx:1441
Trk::RungeKuttaPropagator::propagateParameters
virtual std::unique_ptr< TrackParameters > propagateParameters(const EventContext &ctx, const TrackParameters &, const Surface &, const PropDirection, const BoundaryCheck &, const MagneticFieldProperties &, ParticleHypothesis, bool, const TrackingVolume *) const override final
Main propagation method for parameters only.
Trk::RungeKuttaPropagator::propagate
virtual std::unique_ptr< NeutralParameters > propagate(const NeutralParameters &, const Surface &, PropDirection, const BoundaryCheck &, bool) const override final
Main propagation method for NeutralParameters.
Definition RungeKuttaPropagator.cxx:1402
Trk::RungeKuttaPropagator::getInitializedCache
Cache getInitializedCache(const EventContext &ctx) const
Definition RungeKuttaPropagator.cxx:1914
Trk::RungeKuttaPropagator::intersect
virtual std::optional< Trk::TrackSurfaceIntersection > intersect(const EventContext &ctx, const TrackParameters &, const Surface &, const MagneticFieldProperties &, ParticleHypothesis, const TrackingVolume *tvol=nullptr) const override final
Global position together with direction of the trajectory on the surface.
Definition RungeKuttaPropagator.cxx:1796
Trk::RungeKuttaPropagator::globalPositions
virtual void globalPositions(const EventContext &ctx, std::deque< Amg::Vector3D > &, const TrackParameters &, const MagneticFieldProperties &, const CylinderBounds &, double, ParticleHypothesis, const TrackingVolume *tvol=nullptr) const override final
GlobalPositions list interface:
Trk::Surface
Abstract Base Class for tracking surfaces.
Definition Tracking/TrkDetDescr/TrkSurfaces/TrkSurfaces/Surface.h:79
Trk::Surface::insideBounds
virtual bool insideBounds(const Amg::Vector2D &locpos, double tol1=0., double tol2=0.) const =0
virtual methods to be overwritten by the inherited surfaces
Trk::Surface::createUniqueNeutralParameters
virtual NeutralTrackParametersUniquePtr createUniqueNeutralParameters(double l1, double l2, double phi, double theat, double qop, std::optional< AmgSymMatrix(5)> cov=std::nullopt) const =0
Use the Surface as a ParametersBase constructor, from local parameters - neutral.
Trk::Surface::createUniqueTrackParameters
virtual ChargedTrackParametersUniquePtr createUniqueTrackParameters(double l1, double l2, double phi, double theat, double qop, std::optional< AmgSymMatrix(5)> cov=std::nullopt) const =0
Use the Surface as a ParametersBase constructor, from local parameters - charged.
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.
Trk::TrackingVolume
Full Volume description used in Tracking, it inherits from Volume to get the geometrical structure,...
Definition TrackingVolume.h:119
const
Amg
Definition of ATLAS Math & Geometry primitives (Amg)
Definition AmgStringHelpers.h:19
Amg::Transform3D
Eigen::Affine3d Transform3D
Definition GeoPrimitives.h:46
Amg::Vector2D
Eigen::Matrix< double, 2, 1 > Vector2D
Definition GeoPrimitives.h:48
Amg::Vector3D
Eigen::Matrix< double, 3, 1 > Vector3D
Definition GeoPrimitives.h:47
RungeKuttaUtils
Trk::RungeKuttaUtils is set algorithms for track parameters transformation from local to global and g...
TRT_PAI_physicsConstants::r0
const double r0
electron radius{cm}
Definition TRT_PAI_physicsConstants.h:22
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::stepEstimator
double stepEstimator(int kind, double *ATH_RESTRICT Su, const double *ATH_RESTRICT P, bool &Q)
Definition RungeKuttaUtils.cxx:1212
Trk
Ensure that the ATLAS eigen extensions are properly loaded.
Definition FakeTrackBuilder.h:9
Trk::PropDirection
PropDirection
PropDirection, enum for direction of the propagation.
Definition PropDirection.h:19
Trk::SurfaceType
SurfaceType
This enumerator simplifies the persistency & calculations,.
Definition SurfaceTypes.h:17
Trk::propJacobian
ATH_ALWAYS_INLINE void propJacobian(double *ATH_RESTRICT P, const double *ATH_RESTRICT H0, const double *ATH_RESTRICT H1, const double *ATH_RESTRICT H2, const double *ATH_RESTRICT A, const double *ATH_RESTRICT A0, const double *ATH_RESTRICT A3, const double *ATH_RESTRICT A4, const double *ATH_RESTRICT A6, const double S3)
This provides an inline helper function for updating the jacobian during Runge-Kutta propagation.
Definition JacobianHelper.h:41
Trk::next
@ next
Definition BinningData.h:33
Trk::MultiComponentState
std::vector< ComponentParameters > MultiComponentState
Definition ComponentParameters.h:27
Trk::NeutralParameters
ParametersBase< NeutralParametersDim, Neutral > NeutralParameters
Definition NeutralParameters.h:26
Trk::FastField
@ FastField
call the fast field access method of the FieldSvc
Definition MagneticFieldMode.h:20
Trk::NoField
@ NoField
Field is set to 0., 0., 0.,.
Definition MagneticFieldMode.h:18
Trk::ParticleHypothesis
ParticleHypothesis
Enumeration for Particle hypothesis respecting the interaction with material.
Definition ParticleHypothesis.h:28
Trk::TrackParameters
ParametersBase< TrackParametersDim, Charged > TrackParameters
Definition Tracking/TrkEvent/TrkParameters/TrkParameters/TrackParameters.h:27
VertexShift::Zmax
const float Zmax
Definition VertexShift.h:26
hist_file_dump.d
d
Definition hist_file_dump.py:142
hist_file_dump.f
f
Definition hist_file_dump.py:140
merge.status
status
Definition merge.py:16
read_hist_ntuple.f2
f2
Definition read_hist_ntuple.py:20
read_hist_ntuple.f1
f1
Definition read_hist_ntuple.py:4
restrict.h
Macro wrapping the nonstandard restrict keyword.
ATH_RESTRICT
#define ATH_RESTRICT
Definition restrict.h:31
A
hold the test vectors and ease the comparison
Trk::RungeKuttaPropagator::Cache::m_fieldCache
MagField::AtlasFieldCache m_fieldCache
Definition RungeKuttaPropagator.h:360
Generated on for ATLAS Offline Software by doxygen 1.14.0