BmnRoot 26.10.0-18-g4e3414f9
Software framework of the BM@N experiment
Loading...
Searching...
No Matches
CbmKFFieldMath.cxx
Go to the documentation of this file.
1#include "CbmKFFieldMath.h"
2
3#include "CbmKFMath.h"
4#include "FairField.h"
5#include "math.h"
6
7// using std::sqrt;
8// using std::finite;
9
10void CbmKFFieldMath::ExtrapolateLine(const Double_t T_in[],
11 const Double_t C_in[],
12 Double_t T_out[],
13 Double_t C_out[],
14 Double_t z_out)
15{
16 if (!T_in)
17 return;
18
19 Double_t dz = z_out - T_in[5];
20
21 if (T_out) {
22 T_out[0] = T_in[0] + dz * T_in[2];
23 T_out[1] = T_in[1] + dz * T_in[3];
24 T_out[2] = T_in[2];
25 T_out[3] = T_in[3];
26 T_out[4] = T_in[4];
27 T_out[5] = z_out;
28 }
29
30 if (C_in && C_out) {
31 const Double_t dzC_in8 = dz * C_in[8];
32
33 C_out[4] = C_in[4] + dzC_in8;
34 C_out[1] = C_in[1] + dz * (C_out[4] + C_in[6]);
35
36 const Double_t C_in3 = C_in[3];
37
38 C_out[3] = C_in3 + dz * C_in[5];
39 C_out[0] = C_in[0] + dz * (C_out[3] + C_in3);
40
41 const Double_t C_in7 = C_in[7];
42
43 C_out[7] = C_in7 + dz * C_in[9];
44 C_out[2] = C_in[2] + dz * (C_out[7] + C_in7);
45 C_out[5] = C_in[5];
46 C_out[6] = C_in[6] + dzC_in8;
47 C_out[8] = C_in[8];
48 C_out[9] = C_in[9];
49
50 C_out[10] = C_in[10];
51 C_out[11] = C_in[11];
52 C_out[12] = C_in[12];
53 C_out[13] = C_in[13];
54 C_out[14] = C_in[14];
55 }
56}
57
58Int_t CbmKFFieldMath::ExtrapolateRK4(const Double_t T_in[], // input track parameters (x,y,tx,ty,Q/p)
59 const Double_t C_in[], // input covariance matrix
60 Double_t T_out[], // output track parameters
61 Double_t C_out[], // output covariance matrix
62 Double_t z_out, // extrapolate to this z position
63 Double_t qp0, // use Q/p linearisation at this value
64 FairField* MF // magnetic field
65)
66{
67 //
68 // Forth-order Runge-Kutta method for solution of the equation
69 // of motion of a particle with parameter qp = Q /P
70 // in the magnetic field B()
71 //
72 // | x | tx
73 // | y | ty
74 // d/dz = | tx| = ft * ( ty * ( B(3)+tx*b(1) ) - (1+tx**2)*B(2) )
75 // | ty| ft * (-tx * ( B(3)+ty+b(2) - (1+ty**2)*B(1) ) ,
76 //
77 // where ft = c_light*qp*sqrt ( 1 + tx**2 + ty**2 ) .
78 //
79 // In the following for RK step :
80 //
81 // x=x[0], y=x[1], tx=x[3], ty=x[4].
82 //
83 //========================================================================
84 //
85 // NIM A395 (1997) 169-184; NIM A426 (1999) 268-282.
86 //
87 // the routine is based on LHC(b) utility code
88 //
89 //========================================================================
90
91 const Double_t c_light = 0.000299792458;
92
93 static Double_t a[4] = {0.0, 0.5, 0.5, 1.0};
94 static Double_t c[4] = {1.0 / 6.0, 1.0 / 3.0, 1.0 / 3.0, 1.0 / 6.0};
95 static Double_t b[4] = {0.0, 0.5, 0.5, 1.0};
96
97 Int_t step4;
98 Double_t k[16], x0[4], x[4], k1[16];
99 Double_t Ax[4], Ay[4], Ax_tx[4], Ay_tx[4], Ax_ty[4], Ay_ty[4];
100
101 //----------------------------------------------------------------
102
103 Double_t qp_in = T_in[4];
104 Double_t z_in = T_in[5];
105 Double_t h = z_out - z_in;
106 Double_t hC = h * c_light;
107 x0[0] = T_in[0];
108 x0[1] = T_in[1];
109 x0[2] = T_in[2];
110 x0[3] = T_in[3];
111 //
112 // Runge-Kutta step
113 //
114
115 Int_t step;
116 Int_t i;
117
118 for (step = 0; step < 4; ++step) {
119 for (i = 0; i < 4; ++i) {
120 if (step == 0) {
121 x[i] = x0[i];
122 } else {
123 x[i] = x0[i] + b[step] * k[step * 4 - 4 + i];
124 }
125 }
126
127 Double_t point[3] = {x[0], x[1], z_in + a[step] * h};
128 Double_t B[3];
129 if (MF)
130 MF->GetFieldValue(point, B);
131 else {
132 B[0] = B[1] = B[2] = 0.;
133 }
134
135 Double_t tx = x[2];
136 Double_t ty = x[3];
137 Double_t tx2 = tx * tx;
138 Double_t ty2 = ty * ty;
139 Double_t txty = tx * ty;
140 Double_t tx2ty21 = 1.0 + tx2 + ty2;
141 if (tx2ty21 > 1.e4)
142 return 1;
143 Double_t I_tx2ty21 = 1.0 / tx2ty21 * qp0;
144 Double_t tx2ty2 = sqrt(tx2ty21);
145 // Double_t I_tx2ty2 = qp0 * hC / tx2ty2 ; unsused ???
146 tx2ty2 *= hC;
147 Double_t tx2ty2qp = tx2ty2 * qp0;
148 Ax[step] = (txty * B[0] + ty * B[2] - (1.0 + tx2) * B[1]) * tx2ty2;
149 Ay[step] = (-txty * B[1] - tx * B[2] + (1.0 + ty2) * B[0]) * tx2ty2;
150
151 Ax_tx[step] = Ax[step] * tx * I_tx2ty21 + (ty * B[0] - 2.0 * tx * B[1]) * tx2ty2qp;
152 Ax_ty[step] = Ax[step] * ty * I_tx2ty21 + (tx * B[0] + B[2]) * tx2ty2qp;
153 Ay_tx[step] = Ay[step] * tx * I_tx2ty21 + (-ty * B[1] - B[2]) * tx2ty2qp;
154 Ay_ty[step] = Ay[step] * ty * I_tx2ty21 + (-tx * B[1] + 2.0 * ty * B[0]) * tx2ty2qp;
155
156 step4 = step * 4;
157 k[step4] = tx * h;
158 k[step4 + 1] = ty * h;
159 k[step4 + 2] = Ax[step] * qp0;
160 k[step4 + 3] = Ay[step] * qp0;
161
162 } // end of Runge-Kutta steps
163
164 for (i = 0; i < 4; ++i) {
165 T_out[i] = x0[i] + c[0] * k[i] + c[1] * k[4 + i] + c[2] * k[8 + i] + c[3] * k[12 + i];
166 }
167 T_out[4] = T_in[4];
168 T_out[5] = z_out;
169 //
170 // Derivatives dx/dqp
171 //
172
173 x0[0] = 0.0;
174 x0[1] = 0.0;
175 x0[2] = 0.0;
176 x0[3] = 0.0;
177
178 //
179 // Runge-Kutta step for derivatives dx/dqp
180
181 for (step = 0; step < 4; ++step) {
182 for (i = 0; i < 4; ++i) {
183 if (step == 0) {
184 x[i] = x0[i];
185 } else {
186 x[i] = x0[i] + b[step] * k1[step * 4 - 4 + i];
187 }
188 }
189 step4 = step * 4;
190 k1[step4] = x[2] * h;
191 k1[step4 + 1] = x[3] * h;
192 k1[step4 + 2] = Ax[step] + Ax_tx[step] * x[2] + Ax_ty[step] * x[3];
193 k1[step4 + 3] = Ay[step] + Ay_tx[step] * x[2] + Ay_ty[step] * x[3];
194
195 } // end of Runge-Kutta steps for derivatives dx/dqp
196
197 Double_t J[25];
198
199 for (i = 0; i < 4; ++i) {
200 J[20 + i] = x0[i] + c[0] * k1[i] + c[1] * k1[4 + i] + c[2] * k1[8 + i] + c[3] * k1[12 + i];
201 }
202 J[24] = 1.;
203 //
204 // end of derivatives dx/dqp
205 //
206
207 // Derivatives dx/tx
208 //
209
210 x0[0] = 0.0;
211 x0[1] = 0.0;
212 x0[2] = 1.0;
213 x0[3] = 0.0;
214
215 //
216 // Runge-Kutta step for derivatives dx/dtx
217 //
218
219 for (step = 0; step < 4; ++step) {
220 for (i = 0; i < 4; ++i) {
221 if (step == 0) {
222 x[i] = x0[i];
223 } else if (i != 2) {
224 x[i] = x0[i] + b[step] * k1[step * 4 - 4 + i];
225 }
226 }
227 step4 = step * 4;
228 k1[step4] = x[2] * h;
229 k1[step4 + 1] = x[3] * h;
230 // k1[step4+2] = Ax_tx[step] * x[2] + Ax_ty[step] * x[3];
231 k1[step4 + 3] = Ay_tx[step] * x[2] + Ay_ty[step] * x[3];
232
233 } // end of Runge-Kutta steps for derivatives dx/dtx
234
235 for (i = 0; i < 4; ++i) {
236 if (i != 2) {
237 J[10 + i] = x0[i] + c[0] * k1[i] + c[1] * k1[4 + i] + c[2] * k1[8 + i] + c[3] * k1[12 + i];
238 }
239 }
240 // end of derivatives dx/dtx
241 J[12] = 1.0;
242 J[14] = 0.0;
243
244 // Derivatives dx/ty
245 //
246
247 x0[0] = 0.0;
248 x0[1] = 0.0;
249 x0[2] = 0.0;
250 x0[3] = 1.0;
251
252 //
253 // Runge-Kutta step for derivatives dx/dty
254 //
255
256 for (step = 0; step < 4; ++step) {
257 for (i = 0; i < 4; ++i) {
258 if (step == 0) {
259 x[i] = x0[i]; // ty fixed
260 } else if (i != 3) {
261 x[i] = x0[i] + b[step] * k1[step * 4 - 4 + i];
262 }
263 }
264 step4 = step * 4;
265 k1[step4] = x[2] * h;
266 k1[step4 + 1] = x[3] * h;
267 k1[step4 + 2] = Ax_tx[step] * x[2] + Ax_ty[step] * x[3];
268 // k1[step4+3] = Ay_tx[step] * x[2] + Ay_ty[step] * x[3];
269
270 } // end of Runge-Kutta steps for derivatives dx/dty
271
272 for (i = 0; i < 3; ++i) {
273 J[15 + i] = x0[i] + c[0] * k1[i] + c[1] * k1[4 + i] + c[2] * k1[8 + i] + c[3] * k1[12 + i];
274 }
275 // end of derivatives dx/dty
276 J[18] = 1.;
277 J[19] = 0.;
278
279 //
280 // derivatives dx/dx and dx/dy
281
282 for (i = 0; i < 10; ++i) {
283 J[i] = 0.;
284 }
285 J[0] = 1.;
286 J[6] = 1.;
287
288 // extrapolate inverse momentum
289
290 T_out[4] = qp_in;
291
292 Double_t dqp = qp_in - qp0;
293
294 {
295 for (Int_t ip = 0; ip < 4; ip++) {
296 T_out[ip] += J[5 * 4 + ip] * dqp;
297 }
298 }
299
300 // covariance matrix transport
301
302 if (C_in && C_out)
303 CbmKFMath::multQtSQ(5, J, C_in, C_out);
304 return 0;
305} // end of RK4
306
307#ifdef XXX
308/****************************************************************
309 *
310 * Field Integration
311 *
312 ****************************************************************/
313
314void CbmKFFieldMath::IntegrateField(CbmMagField* MF,
315 Double_t r0[],
316 Double_t r1[],
317 Double_t r2[],
318 Double_t si[3],
319 Double_t Si[3],
320 Double_t sii[3][3],
321 Double_t Sii[3][3],
322 Double_t siii[3][3][3],
323 Double_t Siii[3][3][3])
324{
325 Double_t dz = r2[2] - r0[2];
326
327 Double_t B[3][3];
328
329 if (MF) {
330 MF->GetFieldValue(r0, B[0]);
331 MF->GetFieldValue(r1, B[1]);
332 MF->GetFieldValue(r2, B[2]);
333 } else
334 B[0][0] = B[0][1] = B[0][2] = B[1][0] = B[1][1] = B[1][2] = B[2][0] = B[2][1] = B[2][2] = 0;
335
336 // coefficients
337
338 Double_t c1[3] = {1, 4, 1} // /=6.
339 ,
340 c2[3][3] = {{5, -4, -1}, {44, 80, -4}, {11, 44, 5}} // /=360.
341 ,
342 c3[3][3][3] = {{{35, 20, -1}, {-124, -256, 20}, {-19, -52, -1}},
343 {{524, 176, -52}, {1760, 2240, -256}, {-52, 176, 20}},
344 {{125, -52, -19}, {1028, 1760, -124}, {125, 524, 35}}}; // /=45360.
345
346 Double_t C1[3] = {1, 2, 0} // /=6.
347 ,
348 C2[3][3] = {{38, 8, -4}, {148, 208, -20}, {3, 36, 3}} // /=2520.
349 ,
350 C3[3][3][3] = {{{85, 28, -5}, {4, -80, 4}, {-17, -20, 1}},
351 {{494, 200, -46}, {1256, 1376, -184}, {-94, 8, 14}},
352 {{15, -12, -3}, {252, 432, -36}, {21, 84, 3}}}; // /=90720.
353
354 // integrate field
355
356 for (Int_t x = 0; x < 3; x++) {
357 if (si) {
358 si[x] = 0;
359 for (Int_t n = 0; n < 3; n++)
360 si[x] += c1[n] * B[n][x];
361 si[x] *= dz / 6.;
362 }
363
364 if (Si) {
365 Si[x] = 0;
366 for (Int_t n = 0; n < 3; n++)
367 Si[x] += C1[n] * B[n][x];
368 Si[x] *= dz * dz / 6.;
369 }
370
371 for (Int_t y = 0; y < 3; y++) {
372 if (sii) {
373 sii[x][y] = 0;
374 for (Int_t n = 0; n < 3; n++)
375 for (Int_t m = 0; m < 3; m++)
376 sii[x][y] += c2[n][m] * B[n][x] * B[m][y];
377 sii[x][y] *= dz * dz / 360.;
378 }
379
380 if (Sii) {
381 Sii[x][y] = 0;
382 for (Int_t n = 0; n < 3; n++)
383 for (Int_t m = 0; m < 3; m++)
384 Sii[x][y] += C2[n][m] * B[n][x] * B[m][y];
385 Sii[x][y] *= dz * dz * dz / 2520.;
386 }
387
388 for (Int_t z = 0; z < 3; z++) {
389 if (siii) {
390 siii[x][y][z] = 0;
391 for (Int_t n = 0; n < 3; n++)
392 for (Int_t m = 0; m < 3; m++)
393 for (Int_t k = 0; k < 3; k++)
394 siii[x][y][z] += c3[n][m][k] * B[n][x] * B[m][y] * B[k][z];
395 siii[x][y][z] *= dz * dz * dz / 45360.;
396 }
397
398 if (Siii) {
399 Siii[x][y][z] = 0;
400 for (Int_t n = 0; n < 3; n++)
401 for (Int_t m = 0; m < 3; m++)
402 for (Int_t k = 0; k < 3; k++)
403 Siii[x][y][z] += C3[n][m][k] * B[n][x] * B[m][y] * B[k][z];
404 Siii[x][y][z] *= dz * dz * dz * dz / 90720.;
405 }
406 }
407 }
408 }
409}
410
411/****************************************************************
412 *
413 * Calculate extrapolation coefficients
414 *
415 ****************************************************************/
416
417void CbmKFFieldMath::GetCoefficients(const Double_t x,
418 const Double_t y,
419 Double_t Xi[3][3],
420 Double_t Yi[3][3],
421 Double_t Xii[3][3][3],
422 Double_t Yii[3][3][3],
423 Double_t Xiii[3][3][3][3],
424 Double_t Yiii[3][3][3][3])
425{
426 Double_t
427
428 xx = x * x,
429 xy = x * y, yy = y * y,
430
431 x2 = x * 2, x4 = x * 4, xx31 = xx * 3 + 1, xx33 = xx * 3 + 3, xx82 = xx * 8 + 2, xx86 = xx * 8 + 6,
432 xx153 = xx * 15 + 3, xx159 = xx * 15 + 9,
433
434 y2 = y * 2, y4 = y * 4, yy31 = yy * 3 + 1, yy33 = yy * 3 + 3, yy82 = yy * 8 + 2, yy86 = yy * 8 + 6,
435 yy153 = yy * 15 + 3, yy159 = yy * 15 + 9, xxyy2 = 2 * (xx - yy),
436
437 xx1053 = y * (30 * xx + xx159), yy1053 = x * (30 * yy + yy159);
438
439 Double_t X1[3] = {xy, -xx - 1, y}
440
441 ,
442 X2[3][3] = {{x * yy31, -y * xx31, 2 * yy + 1}, {-y * xx31, x * xx33, -2 * xy}, {yy - xx, -2 * xy, -x}}
443
444 ,
445 X3[3][3][3] = {{{xy * yy159, -xx * yy153 - yy31, y * yy86},
446 {-xx * yy153 - yy31, xy * xx159, -x * yy82},
447 {y2 * (-xxyy2 + 1), -x * yy82, -3 * xy}},
448 {{-xx * yy153 - yy31, xy * xx159, -x * yy82},
449 {xy * xx159, -3 * (5 * xx * xx + 6 * xx + 1), y * xx82},
450 {x2 * (xxyy2 + 1), y * xx82, xx31}},
451 {{y * (-6 * xx + yy31), x * (xx31 - 6 * yy), -4 * xy},
452 {x * (3 * xx - 6 * yy + 1), 3 * y * xx31, xxyy2},
453 {-4 * xy, xxyy2, -y}}};
454 Double_t X1x[3] = {y, -x2, 0}
455
456 ,
457 X2x[3][3] = {{yy31, -6 * xy, 0}, {-6 * xy, 3 * xx31, -y2}, {-x2, -y2, -1}}
458
459 ,
460 X3x[3][3][3] = {{{y * yy159, -x2 * yy153, 0}, {-x2 * yy153, xx1053, -yy82}, {-8 * xy, -yy82, -3 * y}},
461 {{-x2 * yy153, xx1053, -yy82},
462 {xx1053, -x4 * xx159, 16 * xy},
463 {2 * (6 * xx - 2 * yy + 1), 16 * xy, 6 * x}},
464 {{-12 * xy, 9 * xx - 6 * yy + 1, -y4}, {9 * xx - 6 * yy + 1, 18 * xy, x4}, {-y4, x4, 0}}};
465
466 Double_t X1y[3] = {x, 0, 1}
467
468 ,
469 X2y[3][3] = {{6 * xy, -xx31, y4}, {-xx31, 0, -x2}, {y2, -x2, 0}}
470
471 ,
472 X3y[3][3][3] = {{{yy1053, -y2 * xx153, 16 * yy + yy86},
473 {-y2 * xx153, x * xx159, -16 * xy},
474 {yy82 - 2 * xxyy2, -16 * xy, -3 * x}},
475 {{-y2 * xx153, x * xx159, -16 * xy}, {x * xx159, 0, xx82}, {-8 * xy, xx82, 0}},
476 {{-6 * xx + 9 * yy + 1, -12 * xy, -x4}, {-12 * xy, 3 * xx31, -y4}, {-x4, -y4, -1}}};
477
478 Double_t Y1[3] = {yy + 1, -xy, -x}
479
480 ,
481 Y2[3][3] = {{y * yy33, -x * yy31, -2 * xy}, {-x * yy31, y * xx31, 2 * xx + 1}, {-2 * xy, xx - yy, -y}}
482
483 ,
484 Y3[3][3][3] = {{{3 * (5 * yy * yy + 6 * yy + 1), -xy * yy159, -x * yy82},
485 {-xy * yy159, xx * yy153 + yy31, y * xx82},
486 {-x * yy82, -y2 * (-xxyy2 + 1), -yy31}},
487 {{-xy * yy159, xx * yy153 + yy31, y * xx82},
488 {xx * yy153 + yy31, -xy * xx159, -x * xx86},
489 {y * xx82, -x2 * (xxyy2 + 1), 3 * xy}},
490 {{-3 * x * yy31, y * (6 * xx - yy31), xxyy2},
491 {y * (6 * xx - yy31), x * (6 * yy - xx31), 4 * xy},
492 {xxyy2, 4 * xy, x}}};
493 Double_t Y1x[3] = {0, -y, -1}
494
495 ,
496 Y2x[3][3] = {{0, -yy31, -y2}, {-yy31, 6 * xy, x4}, {-y2, x2, 0}}
497
498 ,
499 Y3x[3][3][3] = {{{0, -y * yy159, -yy82}, {-y * yy159, x2 * yy153, 16 * xy}, {-yy82, 8 * xy, 0}},
500 {{-y * yy159, x2 * yy153, 16 * xy},
501 {x2 * yy153, -xx1053, -16 * xx - xx86},
502 {16 * xy, -2 * (6 * xx - 2 * yy + 1), 3 * y}},
503 {{-3 * yy31, 12 * xy, x4}, {12 * xy, -9 * xx + 6 * yy - 1, y4}, {x4, y4, 1}}};
504 Double_t Y1y[3] = {y2, -x, 0}
505
506 ,
507 Y2y[3][3] = {{3 * yy31, -6 * xy, -x2}, {-6 * xy, xx31, 0}, {-x2, -y2, -1}}
508
509 ,
510 Y3y[3][3][3] = {{{y4 * yy159, -yy1053, -16 * xy},
511 {-yy1053, y2 * xx153, xx82},
512 {-16 * xy, 4 * xx - 12 * yy - 2, -6 * y}},
513 {{-yy1053, y2 * xx153, xx82}, {y2 * xx153, -x * xx159, 0}, {xx82, 8 * xy, 3 * x}},
514 {{-18 * xy, 6 * xx - 9 * yy - 1, -y4}, {6 * xx - 9 * yy - 1, 12 * xy, x4}, {-y4, x4, 0}}};
515
516 if (Xi) {
517 for (Int_t i = 0; i < 3; i++) {
518 Xi[0][i] = X1[i];
519 Xi[1][i] = X1x[i];
520 Xi[2][i] = X1y[i];
521 }
522 }
523 if (Yi) {
524 for (Int_t i = 0; i < 3; i++) {
525 Yi[0][i] = Y1[i];
526 Yi[1][i] = Y1x[i];
527 Yi[2][i] = Y1y[i];
528 }
529 }
530 if (Xii) {
531 for (Int_t i = 0; i < 3; i++)
532 for (Int_t j = 0; j < 3; j++) {
533 Xii[0][i][j] = X2[i][j];
534 Xii[1][i][j] = X2x[i][j];
535 Xii[2][i][j] = X2y[i][j];
536 }
537 }
538 if (Yii) {
539 for (Int_t i = 0; i < 3; i++)
540 for (Int_t j = 0; j < 3; j++) {
541 Yii[0][i][j] = Y2[i][j];
542 Yii[1][i][j] = Y2x[i][j];
543 Yii[2][i][j] = Y2y[i][j];
544 }
545 }
546 if (Xiii) {
547 for (Int_t i = 0; i < 3; i++)
548 for (Int_t j = 0; j < 3; j++)
549 for (Int_t k = 0; k < 3; k++) {
550 Xiii[0][i][j][k] = X3[i][j][k];
551 Xiii[1][i][j][k] = X3x[i][j][k];
552 Xiii[2][i][j][k] = X3y[i][j][k];
553 }
554 }
555 if (Yiii) {
556 for (Int_t i = 0; i < 3; i++)
557 for (Int_t j = 0; j < 3; j++)
558 for (Int_t k = 0; k < 3; k++) {
559 Yiii[0][i][j][k] = Y3[i][j][k];
560 Yiii[1][i][j][k] = Y3x[i][j][k];
561 Yiii[2][i][j][k] = Y3y[i][j][k];
562 }
563 }
564}
565
566/****************************************************************
567 *
568 * ANALYTIC 1-3
569 *
570 ****************************************************************/
571
572void CbmKFFieldMath::ExtrapolateAnalytic(const Double_t T_in[], // input track parameters (x,y,tx,ty,Q/p)
573 const Double_t C_in[], // input covariance matrix
574 Double_t T_out[], // output track parameters
575 Double_t C_out[], // output covariance matrix
576 Double_t z_out, // extrapolate to this z position
577 Double_t qp0, // use Q/p linearisation at this value
578 FairField* MF, // magnetic field
579 Int_t order)
580{
581 //
582 // Part of the exact extrapolation formula with error (c_light*B*dz)^4/4!
583 //
584 // Under investigation, finally supposed to be fast.
585 //
586
587 const Double_t c_light = 0.000299792458;
588
589 Double_t qp_in = T_in[4];
590 Double_t z_in = T_in[5];
591 Double_t dz = z_out - z_in;
592
593 // construct coefficients
594
595 Double_t tx = T_in[2], ty = T_in[3],
596
597 A1[3][3], B1[3][3], A2[3][3][3], B2[3][3][3], A3[3][3][3][3], B3[3][3][3][3];
598
599 GetCoefficients(tx, ty, A1, B1, A2, B2, A3, B3);
600
601 Double_t t2 = 1. + tx * tx + ty * ty, t = sqrt(t2), h = qp0 * c_light, ht = h * t;
602
603 Double_t s1[3], s2[3][3], s3[3][3][3], S1[3], S2[3][3], S3[3][3][3];
604
605 {
606 Double_t r0[3], r1[3], r2[3];
607
608 r0[0] = T_in[0];
609 r0[1] = T_in[1];
610 r0[2] = T_in[5];
611
612 r2[0] = T_in[0] + T_in[2] * dz;
613 r2[1] = T_in[1] + T_in[3] * dz;
614 r2[2] = z_out;
615
616 r1[0] = 0.5 * (r0[0] + r2[0]);
617 r1[1] = 0.5 * (r0[1] + r2[1]);
618 r1[2] = 0.5 * (r0[2] + r2[2]);
619
620 Double_t B[3][3];
621
622 if (MF) {
623 MF->GetFieldValue(r0, B[0]);
624 MF->GetFieldValue(r1, B[1]);
625 MF->GetFieldValue(r2, B[2]);
626 } else
627 B[0][0] = B[0][1] = B[0][2] = B[1][0] = B[1][1] = B[1][2] = B[2][0] = B[2][1] = B[2][2] = 0;
628
629 Double_t Sy = (7 * B[0][1] + 6 * B[1][1] - B[2][1]) * dz * dz / 96.;
630 r1[0] = T_in[0] + tx * dz / 2 + ht * Sy * A1[0][1];
631 r1[1] = T_in[1] + ty * dz / 2 + ht * Sy * B1[0][1];
632
633 Sy = (B[0][1] + 2 * B[1][1]) * dz * dz / 6.;
634 r2[0] = T_in[0] + tx * dz + ht * Sy * A1[0][1];
635 r2[1] = T_in[1] + ty * dz + ht * Sy * B1[0][1];
636
637 IntegrateField(MF, r0, r1, r2, s1, S1, s2, S2, s3, S3);
638 }
639
640 Double_t sA1 = 0, sA2 = 0, sA3 = 0, sB1 = 0, sB2 = 0, sB3 = 0;
641 Double_t sA1x = 0, sA2x = 0, sA3x = 0, sB1x = 0, sB2x = 0, sB3x = 0;
642 Double_t sA1y = 0, sA2y = 0, sA3y = 0, sB1y = 0, sB2y = 0, sB3y = 0;
643
644 Double_t SA1 = 0, SA2 = 0, SA3 = 0, SB1 = 0, SB2 = 0, SB3 = 0;
645 Double_t SA1x = 0, SA2x = 0, SA3x = 0, SB1x = 0, SB2x = 0, SB3x = 0;
646 Double_t SA1y = 0, SA2y = 0, SA3y = 0, SB1y = 0, SB2y = 0, SB3y = 0;
647
648 {
649 for (Int_t n = 0; n < 3; n++) {
650 sA1 += s1[n] * A1[0][n];
651 sA1x += s1[n] * A1[1][n];
652 sA1y += s1[n] * A1[2][n];
653 sB1 += s1[n] * B1[0][n];
654 sB1x += s1[n] * B1[1][n];
655 sB1y += s1[n] * B1[2][n];
656
657 SA1 += S1[n] * A1[0][n];
658 SA1x += S1[n] * A1[1][n];
659 SA1y += S1[n] * A1[2][n];
660 SB1 += S1[n] * B1[0][n];
661 SB1x += S1[n] * B1[1][n];
662 SB1y += S1[n] * B1[2][n];
663
664 if (order > 1)
665 for (Int_t m = 0; m < 3; m++) {
666 sA2 += s2[n][m] * A2[0][n][m];
667 sA2x += s2[n][m] * A2[1][n][m];
668 sA2y += s2[n][m] * A2[2][n][m];
669 sB2 += s2[n][m] * B2[0][n][m];
670 sB2x += s2[n][m] * B2[1][n][m];
671 sB2y += s2[n][m] * B2[2][n][m];
672
673 SA2 += S2[n][m] * A2[0][n][m];
674 SA2x += S2[n][m] * A2[1][n][m];
675 SA2y += S2[n][m] * A2[2][n][m];
676 SB2 += S2[n][m] * B2[0][n][m];
677 SB2x += S2[n][m] * B2[1][n][m];
678 SB2y += S2[n][m] * B2[2][n][m];
679
680 if (order > 2)
681 for (Int_t k = 0; k < 3; k++) {
682 sA3 += s3[n][m][k] * A3[0][n][m][k];
683 sA3x += s3[n][m][k] * A3[1][n][m][k];
684 sA3y += s3[n][m][k] * A3[2][n][m][k];
685 sB3 += s3[n][m][k] * B3[0][n][m][k];
686 sB3x += s3[n][m][k] * B3[1][n][m][k];
687 sB3y += s3[n][m][k] * B3[2][n][m][k];
688
689 SA3 += S3[n][m][k] * A3[0][n][m][k];
690 SA3x += S3[n][m][k] * A3[1][n][m][k];
691 SA3y += S3[n][m][k] * A3[2][n][m][k];
692 SB3 += S3[n][m][k] * B3[0][n][m][k];
693 SB3x += S3[n][m][k] * B3[1][n][m][k];
694 SB3y += S3[n][m][k] * B3[2][n][m][k];
695 }
696 }
697 }
698 }
699
700 T_out[0] = T_in[0] + tx * dz + ht * (SA1 + ht * (SA2 + ht * SA3));
701 T_out[1] = T_in[1] + ty * dz + ht * (SB1 + ht * (SB2 + ht * SB3));
702 T_out[2] = T_in[2] + ht * (sA1 + ht * (sA2 + ht * sA3));
703 T_out[3] = T_in[3] + ht * (sB1 + ht * (sB2 + ht * sB3));
704 T_out[4] = T_in[4];
705 T_out[5] = z_out;
706
707 Double_t J[25];
708
709 // derivatives '_x
710
711 J[0] = 1;
712 J[1] = 0;
713 J[2] = 0;
714 J[3] = 0;
715 J[4] = 0;
716
717 // derivatives '_y
718
719 J[5] = 0;
720 J[6] = 1;
721 J[7] = 0;
722 J[8] = 0;
723 J[9] = 0;
724
725 // derivatives '_tx
726
727 J[10] = dz + h * tx * (1. / t * SA1 + h * (2 * SA2 + 3 * ht * SA3)) + ht * (SA1x + ht * (SA2x + ht * SA3x));
728 J[11] = h * tx * (1. / t * SB1 + h * (2 * SB2 + 3 * ht * SB3)) + ht * (SB1x + ht * (SB2x + ht * SB3x));
729 J[12] = 1 + h * tx * (1. / t * sA1 + h * (2 * sA2 + 3 * ht * sA3)) + ht * (sA1x + ht * (sA2x + ht * sA3x));
730 J[13] = h * tx * (1. / t * sB1 + h * (2 * sB2 + 3 * ht * sB3)) + ht * (sB1x + ht * (sB2x + ht * sB3x));
731 J[14] = 0;
732
733 // derivatives '_ty
734
735 J[15] = h * ty * (1. / t * SA1 + h * (2 * SA2 + 3 * ht * SA3)) + ht * (SA1y + ht * (SA2y + ht * SA3y));
736 J[16] = dz + h * ty * (1. / t * SB1 + h * (2 * SB2 + 3 * ht * SB3)) + ht * (SB1y + ht * (SB2y + ht * SB3y));
737 J[17] = h * ty * (1. / t * sA1 + h * (2 * sA2 + 3 * ht * sA3)) + ht * (sA1y + ht * (sA2y + ht * sA3y));
738 J[18] = 1 + h * ty * (1. / t * sB1 + h * (2 * sB2 + 3 * ht * sB3)) + ht * (sB1y + ht * (sB2y + ht * sB3y));
739 J[19] = 0;
740
741 // derivatives '_qp
742
743 J[20] = c_light * t * (SA1 + ht * (2 * SA2 + 3 * ht * SA3));
744 J[21] = c_light * t * (SB1 + ht * (2 * SB2 + 3 * ht * SB3));
745 J[22] = c_light * t * (sA1 + ht * (2 * sA2 + 3 * ht * sA3));
746 J[23] = c_light * t * (sB1 + ht * (2 * sB2 + 3 * ht * sB3));
747 J[24] = 1;
748
749 // extrapolate inverse momentum
750
751 T_out[4] = qp_in;
752
753 Double_t dqp = qp_in - qp0;
754
755 {
756 for (Int_t i = 0; i < 4; i++) {
757 T_out[i] += J[5 * 4 + i] * dqp;
758 }
759 }
760
761 // covariance matrix transport
762
763 if (C_in && C_out)
764 CbmKFMath::multQtSQ(5, J, C_in, C_out);
765}
766
767/****************************************************************
768 *
769 * Analytic Central
770 *
771 ****************************************************************/
772
773void CbmKFFieldMath::ExtrapolateACentral(const Double_t T_in[], // input track parameters (x,y,tx,ty,Q/p,z)
774 const Double_t C_in[], // input covariance matrix
775 Double_t T_out[], // output track parameters
776 Double_t C_out[], // output covariance matrix
777 Double_t z_out, // extrapolate to this z position
778 Double_t qp0, // use Q/p linearisation at this value
779 CbmMagField* MF // magnetic field
780)
781{
782 //
783 // Part of the exact extrapolation formula with error (c_light*B*dz)^?/?!
784 //
785 // Under investigation, finally supposed to be fast.
786 //
787
788 const Double_t c_light = 0.000299792458;
789
790 Double_t qp_in = T_in[4];
791 Double_t z_in = T_in[5];
792 Double_t dz = z_out - z_in;
793
794 Double_t i1[3], I1[3];
795
796 { // get integrated field
797
798 // get field value at 3 points
799
800 Double_t r[3][3];
801
802 r[0][0] = 0;
803 r[0][1] = 0;
804 r[0][2] = T_in[5];
805
806 r[2][0] = 0;
807 r[2][1] = 0;
808 r[2][2] = z_out;
809
810 r[1][0] = 0;
811 r[1][1] = 0;
812 r[1][2] = 0.5 * (r[0][2] + r[2][2]);
813
814 IntegrateField(MF, r[0], r[1], r[2], i1, I1, 0, 0, 0, 0);
815
816 } // get integrated field
817
818 Double_t x = T_in[2], // tx !!
819 y = T_in[3], // ty !!
820
821 xx = x * x, xy = x * y, yy = y * y,
822
823 x2 = x * 2, y2 = y * 2;
824
825 Double_t X1[3] = {xy, -xx - 1, y};
826 Double_t X1x[3] = {y, -x2, 0};
827 Double_t X1y[3] = {x, 0, 1};
828
829 Double_t Y1[3] = {yy + 1, -xy, -x};
830 Double_t Y1x[3] = {0, -y, -1};
831 Double_t Y1y[3] = {y2, -x, 0};
832
833 Double_t iX1 = 0, iY1 = 0;
834 Double_t iX1x = 0, iY1x = 0;
835 Double_t iX1y = 0, iY1y = 0;
836
837 Double_t IX1 = 0, IY1 = 0;
838 Double_t IX1x = 0, IY1x = 0;
839 Double_t IX1y = 0, IY1y = 0;
840
841 {
842 for (Int_t n = 0; n < 3; n++) {
843 if (n != 1)
844 continue;
845 iX1 += i1[n] * X1[n];
846 iX1x += i1[n] * X1x[n];
847 iX1y += i1[n] * X1y[n];
848 iY1 += i1[n] * Y1[n];
849 iY1x += i1[n] * Y1x[n];
850 iY1y += i1[n] * Y1y[n];
851
852 IX1 += I1[n] * X1[n];
853 IX1x += I1[n] * X1x[n];
854 IX1y += I1[n] * X1y[n];
855 IY1 += I1[n] * Y1[n];
856 IY1x += I1[n] * Y1x[n];
857 IY1y += I1[n] * Y1y[n];
858 }
859 }
860
861 Double_t t2 = 1. + xx + yy, t = sqrt(t2), h = qp0 * c_light, ht = h * t;
862
863 T_out[0] = T_in[0] + x * dz + ht * IX1;
864 T_out[1] = T_in[1] + y * dz + ht * IY1;
865 T_out[2] = T_in[2] + ht * iX1;
866 T_out[3] = T_in[3] + ht * iY1;
867 T_out[4] = T_in[4];
868 T_out[5] = z_out;
869
870 Double_t J[25];
871
872 // derivatives '_x
873
874 J[0] = 1;
875 J[1] = 0;
876 J[2] = 0;
877 J[3] = 0;
878 J[4] = 0;
879
880 // derivatives '_y
881
882 J[5] = 0;
883 J[6] = 1;
884 J[7] = 0;
885 J[8] = 0;
886 J[9] = 0;
887
888 // derivatives '_tx
889
890 J[10] = dz + h * x / t * IX1 + ht * IX1x;
891 J[11] = h * x / t * IY1 + ht * IY1x;
892 J[12] = 1 + h * x / t * iX1 + ht * iX1x;
893 J[13] = h * x / t * iY1 + ht * iY1x;
894 J[14] = 0;
895
896 // derivatives '_ty
897
898 J[15] = h * y / t * IX1 + ht * IX1y;
899 J[16] = dz + h * y / t * IY1 + ht * IY1y;
900 J[17] = h * y / t * iX1 + ht * iX1y;
901 J[18] = 1 + h * y / t * iY1 + ht * iY1y;
902 J[19] = 0;
903
904 // derivatives '_qp
905
906 J[20] = c_light * t * IX1;
907 J[21] = c_light * t * IY1;
908 J[22] = c_light * t * iX1;
909 J[23] = c_light * t * iY1;
910 J[24] = 1;
911
912 // extrapolate inverse momentum
913
914 T_out[4] = qp_in;
915
916 Double_t dqp = qp_in - qp0;
917
918 {
919 for (Int_t i = 0; i < 4; i++) {
920 T_out[i] += J[5 * 4 + i] * dqp;
921 }
922 }
923
924 // covariance matrix transport
925
926 if (C_in && C_out)
927 CbmKFMath::multQtSQ(5, J, C_in, C_out);
928}
929
930#endif
931
932/****************************************************************
933 *
934 * Analytic Light
935 *
936 ****************************************************************/
937
938Int_t CbmKFFieldMath::ExtrapolateALight(const Double_t T_in[], // input track parameters (x,y,tx,ty,Q/p)
939 const Double_t C_in[], // input covariance matrix
940 Double_t T_out[], // output track parameters
941 Double_t C_out[], // output covariance matrix
942 Double_t z_out, // extrapolate to this z position
943 Double_t qp0, // use Q/p linearisation at this value
944 FairField* MF // magnetic field
945)
946{
947 //
948 // Part of the analytic extrapolation formula with error (c_light*B*dz)^4/4!
949 //
950 {
951 bool ok = 1;
952 for (int i = 0; i < 6; i++)
953 ok = ok && finite(T_in[i]) && (T_in[i] < 1.e5);
954 if (C_in)
955 for (int i = 0; i < 15; i++)
956 ok = ok && finite(C_in[i]);
957 if (!ok) {
958 for (int i = 0; i < 6; i++)
959 T_out[i] = 0;
960 if (C_out) {
961 for (int i = 0; i < 15; i++)
962 C_out[i] = 0;
963 C_out[0] = C_out[2] = C_out[5] = C_out[9] = C_out[14] = 100.;
964 }
965 return 1;
966 }
967 }
968
969 const Double_t c_light = 0.000299792458;
970
971 Double_t qp_in = T_in[4];
972 Double_t z_in = T_in[5];
973 Double_t dz = z_out - z_in;
974
975 // construct coefficients
976
977 Double_t x = T_in[2], // tx !!
978 y = T_in[3], // ty !!
979
980 xx = x * x, xy = x * y, yy = y * y, x2 = x * 2, x4 = x * 4, xx31 = xx * 3 + 1, xx159 = xx * 15 + 9, y2 = y * 2;
981
982 Double_t Ax = xy, Ay = -xx - 1, Az = y, Ayy = x * (xx * 3 + 3), Ayz = -2 * xy, Ayyy = -(15 * xx * xx + 18 * xx + 3),
983
984 Ax_x = y, Ay_x = -x2, Az_x = 0, Ayy_x = 3 * xx31, Ayz_x = -y2, Ayyy_x = -x4 * xx159,
985
986 Ax_y = x, Ay_y = 0, Az_y = 1, Ayy_y = 0, Ayz_y = -x2, Ayyy_y = 0,
987
988 Bx = yy + 1, By = -xy, Bz = -x, Byy = y * xx31, Byz = 2 * xx + 1, Byyy = -xy * xx159,
989
990 Bx_x = 0, By_x = -y, Bz_x = -1, Byy_x = 6 * xy, Byz_x = x4, Byyy_x = -y * (45 * xx + 9),
991
992 Bx_y = y2, By_y = -x, Bz_y = 0, Byy_y = xx31, Byz_y = 0, Byyy_y = -x * xx159;
993
994 // end of coefficients calculation
995
996 Double_t t2 = 1. + xx + yy;
997 if (t2 > 1.e4)
998 return 1;
999 Double_t t = sqrt(t2), h = qp0 * c_light, ht = h * t;
1000
1001 Double_t sx = 0, sy = 0, sz = 0, syy = 0, syz = 0, syyy = 0, Sx = 0, Sy = 0, Sz = 0, Syy = 0, Syz = 0, Syyy = 0;
1002
1003 { // get field integrals
1004
1005 Double_t B[3][3];
1006 Double_t r0[3], r1[3], r2[3];
1007
1008 // first order track approximation
1009
1010 r0[0] = T_in[0];
1011 r0[1] = T_in[1];
1012 r0[2] = T_in[5];
1013
1014 r2[0] = T_in[0] + T_in[2] * dz;
1015 r2[1] = T_in[1] + T_in[3] * dz;
1016 r2[2] = z_out;
1017
1018 r1[0] = 0.5 * (r0[0] + r2[0]);
1019 r1[1] = 0.5 * (r0[1] + r2[1]);
1020 r1[2] = 0.5 * (r0[2] + r2[2]);
1021
1022 MF->GetFieldValue(r0, B[0]);
1023 MF->GetFieldValue(r1, B[1]);
1024 MF->GetFieldValue(r2, B[2]);
1025
1026 Sy = (7 * B[0][1] + 6 * B[1][1] - B[2][1]) * dz * dz / 96.;
1027 r1[0] = T_in[0] + x * dz / 2 + ht * Sy * Ay;
1028 r1[1] = T_in[1] + y * dz / 2 + ht * Sy * By;
1029
1030 Sy = (B[0][1] + 2 * B[1][1]) * dz * dz / 6.;
1031 r2[0] = T_in[0] + x * dz + ht * Sy * Ay;
1032 r2[1] = T_in[1] + y * dz + ht * Sy * By;
1033
1034 Sy = 0;
1035
1036 // integrals
1037
1038 MF->GetFieldValue(r0, B[0]);
1039 MF->GetFieldValue(r1, B[1]);
1040 MF->GetFieldValue(r2, B[2]);
1041
1042 sx = (B[0][0] + 4 * B[1][0] + B[2][0]) * dz / 6.;
1043 sy = (B[0][1] + 4 * B[1][1] + B[2][1]) * dz / 6.;
1044 sz = (B[0][2] + 4 * B[1][2] + B[2][2]) * dz / 6.;
1045
1046 Sx = (B[0][0] + 2 * B[1][0]) * dz * dz / 6.;
1047 Sy = (B[0][1] + 2 * B[1][1]) * dz * dz / 6.;
1048 Sz = (B[0][2] + 2 * B[1][2]) * dz * dz / 6.;
1049
1050 Double_t c2[3][3] = {{5, -4, -1}, {44, 80, -4}, {11, 44, 5}}; // /=360.
1051 Double_t C2[3][3] = {{38, 8, -4}, {148, 208, -20}, {3, 36, 3}}; // /=2520.
1052 for (Int_t n = 0; n < 3; n++)
1053 for (Int_t m = 0; m < 3; m++) {
1054 syz += c2[n][m] * B[n][1] * B[m][2];
1055 Syz += C2[n][m] * B[n][1] * B[m][2];
1056 }
1057
1058 syz *= dz * dz / 360.;
1059 Syz *= dz * dz * dz / 2520.;
1060
1061 syy = (B[0][1] + 4 * B[1][1] + B[2][1]) * dz;
1062 syyy = syy * syy * syy / 1296;
1063 syy = syy * syy / 72;
1064
1065 Syy = (B[0][1] * (38 * B[0][1] + 156 * B[1][1] - B[2][1]) + B[1][1] * (208 * B[1][1] + 16 * B[2][1])
1066 + B[2][1] * (3 * B[2][1]))
1067 * dz * dz * dz / 2520.;
1068 Syyy = (B[0][1]
1069 * (B[0][1] * (85 * B[0][1] + 526 * B[1][1] - 7 * B[2][1])
1070 + B[1][1] * (1376 * B[1][1] + 84 * B[2][1]) + B[2][1] * (19 * B[2][1]))
1071 + B[1][1] * (B[1][1] * (1376 * B[1][1] + 256 * B[2][1]) + B[2][1] * (62 * B[2][1]))
1072 + B[2][1] * B[2][1] * (3 * B[2][1]))
1073 * dz * dz * dz * dz / 90720.;
1074 }
1075
1076 Double_t
1077
1078 sA1 = sx * Ax + sy * Ay + sz * Az,
1079 sA1_x = sx * Ax_x + sy * Ay_x + sz * Az_x, sA1_y = sx * Ax_y + sy * Ay_y + sz * Az_y,
1080
1081 sB1 = sx * Bx + sy * By + sz * Bz, sB1_x = sx * Bx_x + sy * By_x + sz * Bz_x,
1082 sB1_y = sx * Bx_y + sy * By_y + sz * Bz_y,
1083
1084 SA1 = Sx * Ax + Sy * Ay + Sz * Az, SA1_x = Sx * Ax_x + Sy * Ay_x + Sz * Az_x,
1085 SA1_y = Sx * Ax_y + Sy * Ay_y + Sz * Az_y,
1086
1087 SB1 = Sx * Bx + Sy * By + Sz * Bz, SB1_x = Sx * Bx_x + Sy * By_x + Sz * Bz_x,
1088 SB1_y = Sx * Bx_y + Sy * By_y + Sz * Bz_y,
1089
1090 sA2 = syy * Ayy + syz * Ayz, sA2_x = syy * Ayy_x + syz * Ayz_x, sA2_y = syy * Ayy_y + syz * Ayz_y,
1091 sB2 = syy * Byy + syz * Byz, sB2_x = syy * Byy_x + syz * Byz_x, sB2_y = syy * Byy_y + syz * Byz_y,
1092
1093 SA2 = Syy * Ayy + Syz * Ayz, SA2_x = Syy * Ayy_x + Syz * Ayz_x, SA2_y = Syy * Ayy_y + Syz * Ayz_y,
1094 SB2 = Syy * Byy + Syz * Byz, SB2_x = Syy * Byy_x + Syz * Byz_x, SB2_y = Syy * Byy_y + Syz * Byz_y,
1095
1096 sA3 = syyy * Ayyy, sA3_x = syyy * Ayyy_x, sA3_y = syyy * Ayyy_y, sB3 = syyy * Byyy, sB3_x = syyy * Byyy_x,
1097 sB3_y = syyy * Byyy_y,
1098
1099 SA3 = Syyy * Ayyy, SA3_x = Syyy * Ayyy_x, SA3_y = Syyy * Ayyy_y, SB3 = Syyy * Byyy, SB3_x = Syyy * Byyy_x,
1100 SB3_y = Syyy * Byyy_y;
1101
1102 T_out[0] = T_in[0] + x * dz + ht * (SA1 + ht * (SA2 + ht * SA3));
1103 T_out[1] = T_in[1] + y * dz + ht * (SB1 + ht * (SB2 + ht * SB3));
1104 T_out[2] = T_in[2] + ht * (sA1 + ht * (sA2 + ht * sA3));
1105 T_out[3] = T_in[3] + ht * (sB1 + ht * (sB2 + ht * sB3));
1106 T_out[4] = T_in[4];
1107 T_out[5] = z_out;
1108
1109 Double_t J[25];
1110
1111 // derivatives '_x
1112
1113 J[0] = 1;
1114 J[1] = 0;
1115 J[2] = 0;
1116 J[3] = 0;
1117 J[4] = 0;
1118
1119 // derivatives '_y
1120
1121 J[5] = 0;
1122 J[6] = 1;
1123 J[7] = 0;
1124 J[8] = 0;
1125 J[9] = 0;
1126
1127 // derivatives '_tx
1128
1129 J[10] = dz + h * x * (1. / t * SA1 + h * (2 * SA2 + 3 * ht * SA3)) + ht * (SA1_x + ht * (SA2_x + ht * SA3_x));
1130 J[11] = h * x * (1. / t * SB1 + h * (2 * SB2 + 3 * ht * SB3)) + ht * (SB1_x + ht * (SB2_x + ht * SB3_x));
1131 J[12] = 1 + h * x * (1. / t * sA1 + h * (2 * sA2 + 3 * ht * sA3)) + ht * (sA1_x + ht * (sA2_x + ht * sA3_x));
1132 J[13] = h * x * (1. / t * sB1 + h * (2 * sB2 + 3 * ht * sB3)) + ht * (sB1_x + ht * (sB2_x + ht * sB3_x));
1133 J[14] = 0;
1134
1135 // derivatives '_ty
1136
1137 J[15] = h * y * (1. / t * SA1 + h * (2 * SA2 + 3 * ht * SA3)) + ht * (SA1_y + ht * (SA2_y + ht * SA3_y));
1138 J[16] = dz + h * y * (1. / t * SB1 + h * (2 * SB2 + 3 * ht * SB3)) + ht * (SB1_y + ht * (SB2_y + ht * SB3_y));
1139 J[17] = h * y * (1. / t * sA1 + h * (2 * sA2 + 3 * ht * sA3)) + ht * (sA1_y + ht * (sA2_y + ht * sA3_y));
1140 J[18] = 1 + h * y * (1. / t * sB1 + h * (2 * sB2 + 3 * ht * sB3)) + ht * (sB1_y + ht * (sB2_y + ht * sB3_y));
1141 J[19] = 0;
1142
1143 // derivatives '_qp
1144
1145 J[20] = c_light * t * (SA1 + ht * (2 * SA2 + 3 * ht * SA3));
1146 J[21] = c_light * t * (SB1 + ht * (2 * SB2 + 3 * ht * SB3));
1147 J[22] = c_light * t * (sA1 + ht * (2 * sA2 + 3 * ht * sA3));
1148 J[23] = c_light * t * (sB1 + ht * (2 * sB2 + 3 * ht * sB3));
1149 J[24] = 1;
1150
1151 // extrapolate inverse momentum
1152
1153 T_out[4] = qp_in;
1154
1155 Double_t dqp = qp_in - qp0;
1156
1157 {
1158 for (Int_t i = 0; i < 4; i++) {
1159 T_out[i] += J[5 * 4 + i] * dqp;
1160 }
1161 }
1162
1163 // covariance matrix transport
1164
1165 if (C_in && C_out)
1166 CbmKFMath::multQtSQ(5, J, C_in, C_out);
1167 return 0;
1168}
i
int i
Definition P4_F32vec4.h:22
m
__m128 m
Definition P4_F32vec4.h:27
Int_t
int Int_t
Double_t
double Double_t
iY1
@ iY1
Definition BmnAlignDefines.h:95
iX1
@ iX1
Definition BmnAlignDefines.h:93
CbmKFFieldMath::ExtrapolateRK4
static Int_t ExtrapolateRK4(const Double_t T_in[], const Double_t C_in[], Double_t T_out[], Double_t C_out[], Double_t z_out, Double_t qp0, FairField *MF)
Definition CbmKFFieldMath.cxx:58
CbmKFFieldMath::ExtrapolateLine
static void ExtrapolateLine(const Double_t T_in[], const Double_t C_in[], Double_t T_out[], Double_t C_out[], Double_t z_out)
Definition CbmKFFieldMath.cxx:10
CbmKFFieldMath::ExtrapolateALight
static Int_t ExtrapolateALight(const Double_t T_in[], const Double_t C_in[], Double_t T_out[], Double_t C_out[], Double_t z_out, Double_t qp0, FairField *MF)
Definition CbmKFFieldMath.cxx:938
CbmKFMath::multQtSQ
static void multQtSQ(Int_t N, const Double_t Q[], const Double_t S[], Double_t S_out[])
Definition CbmKFMath.cxx:68
sqrt
RVec< PromoteType< T > > sqrt(const RVec< T > &v)
n
*brief user defined string literal for JSON values *sa std::size_t n
Definition json.hpp:20362
NS_L1TrackFitter::c_light
const fvec c_light
Definition L1TrackFitter.cxx:29