CbmRoot
LitExtrapolation.h
Go to the documentation of this file.
1 
18 #ifndef LITEXTRAPOLATION_H_
19 #define LITEXTRAPOLATION_H_
20 
21 #include "LitFieldGrid.h"
22 #include "LitFieldRegion.h"
23 #include "LitMath.h"
24 #include "LitTrackParam.h"
25 
26 namespace lit {
27  namespace parallel {
28 
36  template<class T>
37  inline void LitLineExtrapolation(LitTrackParam<T>& par, T zOut) {
38  T dz = zOut - par.Z;
39 
40  // transport state vector F*X*F.T()
41  par.X += dz * par.Tx;
42  par.Y += dz * par.Ty;
43 
44  // transport covariance matrix F*C*F.T()
45  T t3 = par.C2 + dz * par.C9;
46  T t7 = dz * par.C10;
47  T t8 = par.C3 + t7;
48  T t19 = par.C7 + dz * par.C12;
49  par.C0 += dz * par.C2 + t3 * dz;
50  par.C1 += dz * par.C6 + t8 * dz;
51  par.C2 = t3;
52  par.C3 = t8;
53  par.C4 += dz * par.C11;
54  par.C5 += dz * par.C7 + t19 * dz;
55  par.C6 += t7;
56  par.C7 = t19;
57  par.C8 += dz * par.C13;
58 
59  par.Z = zOut;
60  }
61 
72  template<class T>
73  inline void LitRK4Extrapolation(LitTrackParam<T>& par,
74  T zOut,
75  const LitFieldGrid& field1,
76  const LitFieldGrid& field2,
77  const LitFieldGrid& field3) {
78  static const T fC = 0.000299792458;
79  static const T ZERO = 0., ONE = 1., TWO = 2., C1_3 = 1. / 3.,
80  C1_6 = 1. / 6.;
81 
82  T coef[4] = {0.0, 0.5, 0.5, 1.0};
83 
84  T Ax[4], Ay[4];
85  T dAx_dtx[4], dAy_dtx[4], dAx_dty[4], dAy_dty[4];
86  T kx[4];
87  T ky[4];
88  T ktx[4];
89  T kty[4];
90 
91  T zIn = par.Z;
92  T h = zOut - zIn;
93  T hC = h * fC;
94  T hCqp = h * fC * par.Qp;
95  T x0[4];
96 
97  T x[4] = {par.X, par.Y, par.Tx, par.Ty};
98 
99  T F[25]; // derivatives, transport matrix
100 
101  // Field values for each step
102  LitFieldValue<T> Bfield[4];
103  // Field grid for each step
104  const LitFieldGrid* Bgrid[4] = {&field1, &field2, &field2, &field3};
105  // field.GetFieldValue(zIn + coef[0] * h, B[0]);
106  // field.GetFieldValue(zIn + coef[1] * h, B[1]);
107  // B[2] = B[1];
108  // field.GetFieldValue(zIn + coef[3] * h, B[3]);
109  // B[0] = field1;
110  // B[1] = field2;
111  // B[2] = B[1];
112  // B[3] = field3;
113 
114  // Calculation for zero step
115  {
116  Bgrid[0]->GetFieldValue(x[0], x[1], Bfield[0]);
117 
118  T Bx = Bfield[0].Bx;
119  T By = Bfield[0].By;
120  T Bz = Bfield[0].Bz;
121 
122  T tx = x[2];
123  T ty = x[3];
124  T txtx = tx * tx;
125  T tyty = ty * ty;
126  T txty = tx * ty;
127  T txtxtyty1 = ONE + txtx + tyty;
128  T t1 = sqrt(txtxtyty1);
129  T t2 = rcp(txtxtyty1); //1.0 / txtxtyty1;
130 
131  Ax[0] = (txty * Bx + ty * Bz - (ONE + txtx) * By) * t1;
132  Ay[0] = (-txty * By - tx * Bz + (ONE + tyty) * Bx) * t1;
133 
134  dAx_dtx[0] = Ax[0] * tx * t2 + (ty * Bx - TWO * tx * By) * t1;
135  dAx_dty[0] = Ax[0] * ty * t2 + (tx * Bx + Bz) * t1;
136  dAy_dtx[0] = Ay[0] * tx * t2 + (-ty * By - Bz) * t1;
137  dAy_dty[0] = Ay[0] * ty * t2 + (-tx * By + TWO * ty * Bx) * t1;
138 
139  kx[0] = tx * h;
140  ky[0] = ty * h;
141  ktx[0] = Ax[0] * hCqp;
142  kty[0] = Ay[0] * hCqp;
143  }
144  // end calculation for zero step
145 
146  // Calculate for steps starting from 1
147  for (unsigned int iStep = 1; iStep < 4; iStep++) { // 1
148  x[0] = par.X + coef[iStep] * kx[iStep - 1];
149  x[1] = par.Y + coef[iStep] * ky[iStep - 1];
150  x[2] = par.Tx + coef[iStep] * ktx[iStep - 1];
151  x[3] = par.Ty + coef[iStep] * kty[iStep - 1];
152 
153  Bgrid[iStep]->GetFieldValue(x[0], x[1], Bfield[iStep]);
154 
155  T Bx = Bfield[iStep].Bx;
156  T By = Bfield[iStep].By;
157  T Bz = Bfield[iStep].Bz;
158 
159  T tx = x[2];
160  T ty = x[3];
161  T txtx = tx * tx;
162  T tyty = ty * ty;
163  T txty = tx * ty;
164  T txtxtyty1 = ONE + txtx + tyty;
165  T t1 = sqrt(txtxtyty1);
166  T t2 = rcp(txtxtyty1); //1.0 / txtxtyty1;
167 
168  Ax[iStep] = (txty * Bx + ty * Bz - (ONE + txtx) * By) * t1;
169  Ay[iStep] = (-txty * By - tx * Bz + (ONE + tyty) * Bx) * t1;
170 
171  dAx_dtx[iStep] = Ax[iStep] * tx * t2 + (ty * Bx - TWO * tx * By) * t1;
172  dAx_dty[iStep] = Ax[iStep] * ty * t2 + (tx * Bx + Bz) * t1;
173  dAy_dtx[iStep] = Ay[iStep] * tx * t2 + (-ty * By - Bz) * t1;
174  dAy_dty[iStep] = Ay[iStep] * ty * t2 + (-tx * By + TWO * ty * Bx) * t1;
175 
176  kx[iStep] = tx * h;
177  ky[iStep] = ty * h;
178  ktx[iStep] = Ax[iStep] * hCqp;
179  kty[iStep] = Ay[iStep] * hCqp;
180 
181  } // 1
182 
183  // Calculate output state vector
184  // for (unsigned int i = 0; i < 4; i++) xOut[i] = CalcOut(xIn[i], k[i]);
185  par.X += kx[0] * C1_6 + kx[1] * C1_3 + kx[2] * C1_3 + kx[3] * C1_6;
186  par.Y += ky[0] * C1_6 + ky[1] * C1_3 + ky[2] * C1_3 + ky[3] * C1_6;
187  par.Tx += ktx[0] * C1_6 + ktx[1] * C1_3 + ktx[2] * C1_3 + ktx[3] * C1_6;
188  par.Ty += kty[0] * C1_6 + kty[1] * C1_3 + kty[2] * C1_3 + kty[3] * C1_6;
189  // xOut[4] = xIn[4];
190 
191 
192  // Calculation of the derivatives
193 
194  // derivatives dx/dx and dx/dy
195  // dx/dx
196  F[0] = ONE;
197  F[5] = ZERO;
198  F[10] = ZERO;
199  F[15] = ZERO;
200  F[20] = ZERO;
201  // dx/dy
202  F[1] = ZERO;
203  F[6] = ONE;
204  F[11] = ZERO;
205  F[16] = ZERO;
206  F[21] = ZERO;
207  // end of derivatives dx/dx and dx/dy
208 
209  // Derivatives dx/tx
210  x[0] = x0[0] = ZERO;
211  x[1] = x0[1] = ZERO;
212  x[2] = x0[2] = ONE;
213  x[3] = x0[3] = ZERO;
214 
215  //Calculate for zero step
216  kx[0] = x[2] * h;
217  ky[0] = x[3] * h;
218  //ktx[0] = (dAx_dtx[iStep] * x[2] + dAx_dty[iStep] * x[3]) * hCqp;
219  kty[0] = (dAy_dtx[0] * x[2] + dAy_dty[0] * x[3]) * hCqp;
220  // Calculate for steps starting from 1
221  for (unsigned int iStep = 1; iStep < 4; iStep++) { // 2
222  x[0] = x0[0] + coef[iStep] * kx[iStep - 1];
223  x[1] = x0[1] + coef[iStep] * ky[iStep - 1];
224  // x[2] = x0[2] + coef[iStep] * ktx[iStep - 1];
225  x[3] = x0[3] + coef[iStep] * kty[iStep - 1];
226 
227  kx[iStep] = x[2] * h;
228  ky[iStep] = x[3] * h;
229  //ktx[iStep] = (dAx_dtx[iStep] * x[2] + dAx_dty[iStep] * x[3]) * hCqp;
230  kty[iStep] = (dAy_dtx[iStep] * x[2] + dAy_dty[iStep] * x[3]) * hCqp;
231  } // 2
232 
233  F[2] = x0[0] + kx[0] * C1_6 + kx[1] * C1_3 + kx[2] * C1_3 + kx[3] * C1_6;
234  F[7] = x0[1] + ky[0] * C1_6 + ky[1] * C1_3 + ky[2] * C1_3 + ky[3] * C1_6;
235  F[12] = ONE;
236  F[17] =
237  x0[3] + kty[0] * C1_6 + kty[1] * C1_3 + kty[2] * C1_3 + kty[3] * C1_6;
238  F[22] = ZERO;
239  // end of derivatives dx/dtx
240 
241  // Derivatives dx/ty
242  x[0] = x0[0] = ZERO;
243  x[1] = x0[1] = ZERO;
244  x[2] = x0[2] = ZERO;
245  x[3] = x0[3] = ONE;
246 
247  //Calculate for zero step
248  kx[0] = x[2] * h;
249  ky[0] = x[3] * h;
250  ktx[0] = (dAx_dtx[0] * x[2] + dAx_dty[0] * x[3]) * hCqp;
251  //kty[0] = (dAy_dtx[iStep] * x[2] + dAy_dty[iStep] * x[3]) * hCqp;
252  //Calculate for steps starting from 1
253  for (unsigned int iStep = 1; iStep < 4; iStep++) { // 4
254  x[0] = x0[0] + coef[iStep] * kx[iStep - 1];
255  x[1] = x0[1] + coef[iStep] * ky[iStep - 1];
256  x[2] = x0[2] + coef[iStep] * ktx[iStep - 1];
257  // x[3] = x0[0] + coef[iStep] * kty[iStep - 1];
258 
259  kx[iStep] = x[2] * h;
260  ky[iStep] = x[3] * h;
261  ktx[iStep] = (dAx_dtx[iStep] * x[2] + dAx_dty[iStep] * x[3]) * hCqp;
262  //kty[iStep] = (dAy_dtx[iStep] * x[2] + dAy_dty[iStep] * x[3]) * hCqp;
263  } // 4
264 
265  F[3] = x0[0] + kx[0] * C1_6 + kx[1] * C1_3 + kx[2] * C1_3
266  + kx[3] * C1_6; // CalcOut(x0[0], k[0]);
267  F[8] = x0[1] + ky[0] * C1_6 + ky[1] * C1_3 + ky[2] * C1_3
268  + ky[3] * C1_6; //CalcOut(x0[1], k[1]);
269  F[13] = x0[2] + ktx[0] * C1_6 + ktx[1] * C1_3 + ktx[2] * C1_3
270  + ktx[3] * C1_6; //CalcOut(x0[2], k[2]);
271  F[18] = ONE;
272  F[23] = ZERO;
273  // end of derivatives dx/dty
274 
275  // Derivatives dx/dqp
276  x[0] = x0[0] = ZERO;
277  x[1] = x0[1] = ZERO;
278  x[2] = x0[2] = ZERO;
279  x[3] = x0[3] = ZERO;
280 
281  //Calculate for zero step
282  kx[0] = x[2] * h;
283  ky[0] = x[3] * h;
284  ktx[0] = Ax[0] * hC + hCqp * (dAx_dtx[0] * x[2] + dAx_dty[0] * x[3]);
285  kty[0] = Ay[0] * hC + hCqp * (dAy_dtx[0] * x[2] + dAy_dty[0] * x[3]);
286  //Calculate for steps starting from 1
287  for (unsigned int iStep = 1; iStep < 4; iStep++) { // 4
288  x[0] = x0[0] + coef[iStep] * kx[iStep - 1];
289  x[1] = x0[1] + coef[iStep] * ky[iStep - 1];
290  x[2] = x0[2] + coef[iStep] * ktx[iStep - 1];
291  x[3] = x0[3] + coef[iStep] * kty[iStep - 1];
292 
293  kx[iStep] = x[2] * h;
294  ky[iStep] = x[3] * h;
295  ktx[iStep] = Ax[iStep] * hC
296  + hCqp * (dAx_dtx[iStep] * x[2] + dAx_dty[iStep] * x[3]);
297  kty[iStep] = Ay[iStep] * hC
298  + hCqp * (dAy_dtx[iStep] * x[2] + dAy_dty[iStep] * x[3]);
299  } // 4
300 
301  F[4] = x0[0] + kx[0] * C1_6 + kx[1] * C1_3 + kx[2] * C1_3 + kx[3] * C1_6;
302  F[9] = x0[1] + ky[0] * C1_6 + ky[1] * C1_3 + ky[2] * C1_3 + ky[3] * C1_6;
303  F[14] =
304  x0[2] + ktx[0] * C1_6 + ktx[1] * C1_3 + ktx[2] * C1_3 + ktx[3] * C1_6;
305  F[19] =
306  x0[3] + kty[0] * C1_6 + kty[1] * C1_3 + kty[2] * C1_3 + kty[3] * C1_6;
307  F[24] = 1.;
308  // end of derivatives dx/dqp
309 
310  // end calculation of the derivatives
311 
312 
313  // Transport C matrix
314  {
315  T cIn[15] = {par.C0,
316  par.C1,
317  par.C2,
318  par.C3,
319  par.C4,
320  par.C5,
321  par.C6,
322  par.C7,
323  par.C8,
324  par.C9,
325  par.C10,
326  par.C11,
327  par.C12,
328  par.C13,
329  par.C14};
330  // F*C*Ft
331  T A = cIn[2] + F[2] * cIn[9] + F[3] * cIn[10] + F[4] * cIn[11];
332  T B = cIn[3] + F[2] * cIn[10] + F[3] * cIn[12] + F[4] * cIn[13];
333  T C = cIn[4] + F[2] * cIn[11] + F[3] * cIn[13] + F[4] * cIn[14];
334 
335  T D = cIn[6] + F[7] * cIn[9] + F[8] * cIn[10] + F[9] * cIn[11];
336  T E = cIn[7] + F[7] * cIn[10] + F[8] * cIn[12] + F[9] * cIn[13];
337  T G = cIn[8] + F[7] * cIn[11] + F[8] * cIn[13] + F[9] * cIn[14];
338 
339  T H = cIn[9] + F[13] * cIn[10] + F[14] * cIn[11];
340  T I = cIn[10] + F[13] * cIn[12] + F[14] * cIn[13];
341  T J = cIn[11] + F[13] * cIn[13] + F[14] * cIn[14];
342 
343  T K = cIn[13] + F[17] * cIn[11] + F[19] * cIn[14];
344 
345  par.C0 = cIn[0] + F[2] * cIn[2] + F[3] * cIn[3] + F[4] * cIn[4]
346  + A * F[2] + B * F[3] + C * F[4];
347  par.C1 = cIn[1] + F[2] * cIn[6] + F[3] * cIn[7] + F[4] * cIn[8]
348  + A * F[7] + B * F[8] + C * F[9];
349  par.C2 = A + B * F[13] + C * F[14];
350  par.C3 = B + A * F[17] + C * F[19];
351  par.C4 = C;
352 
353  par.C5 = cIn[5] + F[7] * cIn[6] + F[8] * cIn[7] + F[9] * cIn[8]
354  + D * F[7] + E * F[8] + G * F[9];
355  par.C6 = D + E * F[13] + G * F[14];
356  par.C7 = E + D * F[17] + G * F[19];
357  par.C8 = G;
358 
359  par.C9 = H + I * F[13] + J * F[14];
360  par.C10 = I + H * F[17] + J * F[19];
361  par.C11 = J;
362 
363  par.C12 = cIn[12] + F[17] * cIn[10] + F[19] * cIn[13]
364  + (F[17] * cIn[9] + cIn[10] + F[19] * cIn[11]) * F[17]
365  + K * F[19];
366  par.C13 = K;
367 
368  par.C14 = cIn[14];
369  }
370  //end transport C matrix
371 
372  par.Z = zOut;
373  }
374 
375 
381  template<class T>
382  inline void LitRK4Extrapolation(LitTrackParam<T>& par,
383  T zOut,
384  const LitFieldValue<T>& field1,
385  const LitFieldValue<T>& field2,
386  const LitFieldValue<T>& field3) {
387  static const T fC = 0.000299792458;
388  static const T ZERO = 0., ONE = 1., TWO = 2., C1_3 = 1. / 3.,
389  C1_6 = 1. / 6.;
390 
391  T coef[4] = {0.0, 0.5, 0.5, 1.0};
392 
393  T Ax[4], Ay[4];
394  T dAx_dtx[4], dAy_dtx[4], dAx_dty[4], dAy_dty[4];
395  T kx[4];
396  T ky[4];
397  T ktx[4];
398  T kty[4];
399 
400  T zIn = par.Z;
401  T h = zOut - zIn;
402  T hC = h * fC;
403  T hCqp = h * fC * par.Qp;
404  T x0[4];
405 
406  T x[4] = {par.X, par.Y, par.Tx, par.Ty};
407 
408  T F[25]; // derivatives, transport matrix
409 
410  // Field values for each step
411  LitFieldValue<T> Bfield[4];
412  // Field grid for each step
413  //const LitFieldGrid* Bgrid[4] = {&field1, &field2, &field2, &field3};
414  // field.GetFieldValue(zIn + coef[0] * h, B[0]);
415  // field.GetFieldValue(zIn + coef[1] * h, B[1]);
416  // B[2] = B[1];
417  // field.GetFieldValue(zIn + coef[3] * h, B[3]);
418  Bfield[0] = field1;
419  Bfield[1] = field2;
420  Bfield[2] = Bfield[1];
421  Bfield[3] = field3;
422 
423  // Calculation for zero step
424  {
425  // Bgrid[0]->GetFieldValue(x[0], x[1], B[0]);
426 
427  T Bx = Bfield[0].Bx;
428  T By = Bfield[0].By;
429  T Bz = Bfield[0].Bz;
430 
431  T tx = x[2];
432  T ty = x[3];
433  T txtx = tx * tx;
434  T tyty = ty * ty;
435  T txty = tx * ty;
436  T txtxtyty1 = ONE + txtx + tyty;
437  T t1 = sqrt(txtxtyty1);
438  T t2 = rcp(txtxtyty1); //1.0 / txtxtyty1;
439 
440  Ax[0] = (txty * Bx + ty * Bz - (ONE + txtx) * By) * t1;
441  Ay[0] = (-txty * By - tx * Bz + (ONE + tyty) * Bx) * t1;
442 
443  dAx_dtx[0] = Ax[0] * tx * t2 + (ty * Bx - TWO * tx * By) * t1;
444  dAx_dty[0] = Ax[0] * ty * t2 + (tx * Bx + Bz) * t1;
445  dAy_dtx[0] = Ay[0] * tx * t2 + (-ty * By - Bz) * t1;
446  dAy_dty[0] = Ay[0] * ty * t2 + (-tx * By + TWO * ty * Bx) * t1;
447 
448  kx[0] = tx * h;
449  ky[0] = ty * h;
450  ktx[0] = Ax[0] * hCqp;
451  kty[0] = Ay[0] * hCqp;
452  }
453  // end calculation for zero step
454 
455  // Calculate for steps starting from 1
456  for (unsigned int iStep = 1; iStep < 4; iStep++) { // 1
457  x[0] = par.X + coef[iStep] * kx[iStep - 1];
458  x[1] = par.Y + coef[iStep] * ky[iStep - 1];
459  x[2] = par.Tx + coef[iStep] * ktx[iStep - 1];
460  x[3] = par.Ty + coef[iStep] * kty[iStep - 1];
461 
462  //Bgrid[iStep]->GetFieldValue(x[0], x[1], B[iStep]);
463 
464  T Bx = Bfield[iStep].Bx;
465  T By = Bfield[iStep].By;
466  T Bz = Bfield[iStep].Bz;
467 
468  T tx = x[2];
469  T ty = x[3];
470  T txtx = tx * tx;
471  T tyty = ty * ty;
472  T txty = tx * ty;
473  T txtxtyty1 = ONE + txtx + tyty;
474  T t1 = sqrt(txtxtyty1);
475  T t2 = rcp(txtxtyty1); //1.0 / txtxtyty1;
476 
477  Ax[iStep] = (txty * Bx + ty * Bz - (ONE + txtx) * By) * t1;
478  Ay[iStep] = (-txty * By - tx * Bz + (ONE + tyty) * Bx) * t1;
479 
480  dAx_dtx[iStep] = Ax[iStep] * tx * t2 + (ty * Bx - TWO * tx * By) * t1;
481  dAx_dty[iStep] = Ax[iStep] * ty * t2 + (tx * Bx + Bz) * t1;
482  dAy_dtx[iStep] = Ay[iStep] * tx * t2 + (-ty * By - Bz) * t1;
483  dAy_dty[iStep] = Ay[iStep] * ty * t2 + (-tx * By + TWO * ty * Bx) * t1;
484 
485  kx[iStep] = tx * h;
486  ky[iStep] = ty * h;
487  ktx[iStep] = Ax[iStep] * hCqp;
488  kty[iStep] = Ay[iStep] * hCqp;
489 
490  } // 1
491 
492  // Calculate output state vector
493  // for (unsigned int i = 0; i < 4; i++) xOut[i] = CalcOut(xIn[i], k[i]);
494  par.X += kx[0] * C1_6 + kx[1] * C1_3 + kx[2] * C1_3 + kx[3] * C1_6;
495  par.Y += ky[0] * C1_6 + ky[1] * C1_3 + ky[2] * C1_3 + ky[3] * C1_6;
496  par.Tx += ktx[0] * C1_6 + ktx[1] * C1_3 + ktx[2] * C1_3 + ktx[3] * C1_6;
497  par.Ty += kty[0] * C1_6 + kty[1] * C1_3 + kty[2] * C1_3 + kty[3] * C1_6;
498  // xOut[4] = xIn[4];
499 
500 
501  // Calculation of the derivatives
502 
503  // derivatives dx/dx and dx/dy
504  // dx/dx
505  F[0] = ONE;
506  F[5] = ZERO;
507  F[10] = ZERO;
508  F[15] = ZERO;
509  F[20] = ZERO;
510  // dx/dy
511  F[1] = ZERO;
512  F[6] = ONE;
513  F[11] = ZERO;
514  F[16] = ZERO;
515  F[21] = ZERO;
516  // end of derivatives dx/dx and dx/dy
517 
518  // Derivatives dx/tx
519  x[0] = x0[0] = ZERO;
520  x[1] = x0[1] = ZERO;
521  x[2] = x0[2] = ONE;
522  x[3] = x0[3] = ZERO;
523 
524  //Calculate for zero step
525  kx[0] = x[2] * h;
526  ky[0] = x[3] * h;
527  //ktx[0] = (dAx_dtx[iStep] * x[2] + dAx_dty[iStep] * x[3]) * hCqp;
528  kty[0] = (dAy_dtx[0] * x[2] + dAy_dty[0] * x[3]) * hCqp;
529  // Calculate for steps starting from 1
530  for (unsigned int iStep = 1; iStep < 4; iStep++) { // 2
531  x[0] = x0[0] + coef[iStep] * kx[iStep - 1];
532  x[1] = x0[1] + coef[iStep] * ky[iStep - 1];
533  // x[2] = x0[2] + coef[iStep] * ktx[iStep - 1];
534  x[3] = x0[3] + coef[iStep] * kty[iStep - 1];
535 
536  kx[iStep] = x[2] * h;
537  ky[iStep] = x[3] * h;
538  //ktx[iStep] = (dAx_dtx[iStep] * x[2] + dAx_dty[iStep] * x[3]) * hCqp;
539  kty[iStep] = (dAy_dtx[iStep] * x[2] + dAy_dty[iStep] * x[3]) * hCqp;
540  } // 2
541 
542  F[2] = x0[0] + kx[0] * C1_6 + kx[1] * C1_3 + kx[2] * C1_3 + kx[3] * C1_6;
543  F[7] = x0[1] + ky[0] * C1_6 + ky[1] * C1_3 + ky[2] * C1_3 + ky[3] * C1_6;
544  F[12] = ONE;
545  F[17] =
546  x0[3] + kty[0] * C1_6 + kty[1] * C1_3 + kty[2] * C1_3 + kty[3] * C1_6;
547  F[22] = ZERO;
548  // end of derivatives dx/dtx
549 
550  // Derivatives dx/ty
551  x[0] = x0[0] = ZERO;
552  x[1] = x0[1] = ZERO;
553  x[2] = x0[2] = ZERO;
554  x[3] = x0[3] = ONE;
555 
556  //Calculate for zero step
557  kx[0] = x[2] * h;
558  ky[0] = x[3] * h;
559  ktx[0] = (dAx_dtx[0] * x[2] + dAx_dty[0] * x[3]) * hCqp;
560  //kty[0] = (dAy_dtx[iStep] * x[2] + dAy_dty[iStep] * x[3]) * hCqp;
561  //Calculate for steps starting from 1
562  for (unsigned int iStep = 1; iStep < 4; iStep++) { // 4
563  x[0] = x0[0] + coef[iStep] * kx[iStep - 1];
564  x[1] = x0[1] + coef[iStep] * ky[iStep - 1];
565  x[2] = x0[2] + coef[iStep] * ktx[iStep - 1];
566  // x[3] = x0[0] + coef[iStep] * kty[iStep - 1];
567 
568  kx[iStep] = x[2] * h;
569  ky[iStep] = x[3] * h;
570  ktx[iStep] = (dAx_dtx[iStep] * x[2] + dAx_dty[iStep] * x[3]) * hCqp;
571  //kty[iStep] = (dAy_dtx[iStep] * x[2] + dAy_dty[iStep] * x[3]) * hCqp;
572  } // 4
573 
574  F[3] = x0[0] + kx[0] * C1_6 + kx[1] * C1_3 + kx[2] * C1_3
575  + kx[3] * C1_6; // CalcOut(x0[0], k[0]);
576  F[8] = x0[1] + ky[0] * C1_6 + ky[1] * C1_3 + ky[2] * C1_3
577  + ky[3] * C1_6; //CalcOut(x0[1], k[1]);
578  F[13] = x0[2] + ktx[0] * C1_6 + ktx[1] * C1_3 + ktx[2] * C1_3
579  + ktx[3] * C1_6; //CalcOut(x0[2], k[2]);
580  F[18] = ONE;
581  F[23] = ZERO;
582  // end of derivatives dx/dty
583 
584  // Derivatives dx/dqp
585  x[0] = x0[0] = ZERO;
586  x[1] = x0[1] = ZERO;
587  x[2] = x0[2] = ZERO;
588  x[3] = x0[3] = ZERO;
589 
590  //Calculate for zero step
591  kx[0] = x[2] * h;
592  ky[0] = x[3] * h;
593  ktx[0] = Ax[0] * hC + hCqp * (dAx_dtx[0] * x[2] + dAx_dty[0] * x[3]);
594  kty[0] = Ay[0] * hC + hCqp * (dAy_dtx[0] * x[2] + dAy_dty[0] * x[3]);
595  //Calculate for steps starting from 1
596  for (unsigned int iStep = 1; iStep < 4; iStep++) { // 4
597  x[0] = x0[0] + coef[iStep] * kx[iStep - 1];
598  x[1] = x0[1] + coef[iStep] * ky[iStep - 1];
599  x[2] = x0[2] + coef[iStep] * ktx[iStep - 1];
600  x[3] = x0[3] + coef[iStep] * kty[iStep - 1];
601 
602  kx[iStep] = x[2] * h;
603  ky[iStep] = x[3] * h;
604  ktx[iStep] = Ax[iStep] * hC
605  + hCqp * (dAx_dtx[iStep] * x[2] + dAx_dty[iStep] * x[3]);
606  kty[iStep] = Ay[iStep] * hC
607  + hCqp * (dAy_dtx[iStep] * x[2] + dAy_dty[iStep] * x[3]);
608  } // 4
609 
610  F[4] = x0[0] + kx[0] * C1_6 + kx[1] * C1_3 + kx[2] * C1_3 + kx[3] * C1_6;
611  F[9] = x0[1] + ky[0] * C1_6 + ky[1] * C1_3 + ky[2] * C1_3 + ky[3] * C1_6;
612  F[14] =
613  x0[2] + ktx[0] * C1_6 + ktx[1] * C1_3 + ktx[2] * C1_3 + ktx[3] * C1_6;
614  F[19] =
615  x0[3] + kty[0] * C1_6 + kty[1] * C1_3 + kty[2] * C1_3 + kty[3] * C1_6;
616  F[24] = 1.;
617  // end of derivatives dx/dqp
618 
619  // end calculation of the derivatives
620 
621 
622  // Transport C matrix
623  {
624  T cIn[15] = {par.C0,
625  par.C1,
626  par.C2,
627  par.C3,
628  par.C4,
629  par.C5,
630  par.C6,
631  par.C7,
632  par.C8,
633  par.C9,
634  par.C10,
635  par.C11,
636  par.C12,
637  par.C13,
638  par.C14};
639  // F*C*Ft
640  T A = cIn[2] + F[2] * cIn[9] + F[3] * cIn[10] + F[4] * cIn[11];
641  T B = cIn[3] + F[2] * cIn[10] + F[3] * cIn[12] + F[4] * cIn[13];
642  T C = cIn[4] + F[2] * cIn[11] + F[3] * cIn[13] + F[4] * cIn[14];
643 
644  T D = cIn[6] + F[7] * cIn[9] + F[8] * cIn[10] + F[9] * cIn[11];
645  T E = cIn[7] + F[7] * cIn[10] + F[8] * cIn[12] + F[9] * cIn[13];
646  T G = cIn[8] + F[7] * cIn[11] + F[8] * cIn[13] + F[9] * cIn[14];
647 
648  T H = cIn[9] + F[13] * cIn[10] + F[14] * cIn[11];
649  T I = cIn[10] + F[13] * cIn[12] + F[14] * cIn[13];
650  T J = cIn[11] + F[13] * cIn[13] + F[14] * cIn[14];
651 
652  T K = cIn[13] + F[17] * cIn[11] + F[19] * cIn[14];
653 
654  par.C0 = cIn[0] + F[2] * cIn[2] + F[3] * cIn[3] + F[4] * cIn[4]
655  + A * F[2] + B * F[3] + C * F[4];
656  par.C1 = cIn[1] + F[2] * cIn[6] + F[3] * cIn[7] + F[4] * cIn[8]
657  + A * F[7] + B * F[8] + C * F[9];
658  par.C2 = A + B * F[13] + C * F[14];
659  par.C3 = B + A * F[17] + C * F[19];
660  par.C4 = C;
661 
662  par.C5 = cIn[5] + F[7] * cIn[6] + F[8] * cIn[7] + F[9] * cIn[8]
663  + D * F[7] + E * F[8] + G * F[9];
664  par.C6 = D + E * F[13] + G * F[14];
665  par.C7 = E + D * F[17] + G * F[19];
666  par.C8 = G;
667 
668  par.C9 = H + I * F[13] + J * F[14];
669  par.C10 = I + H * F[17] + J * F[19];
670  par.C11 = J;
671 
672  par.C12 = cIn[12] + F[17] * cIn[10] + F[19] * cIn[13]
673  + (F[17] * cIn[9] + cIn[10] + F[19] * cIn[11]) * F[17]
674  + K * F[19];
675  par.C13 = K;
676 
677  par.C14 = cIn[14];
678  }
679  //end transport C matrix
680 
681  par.Z = zOut;
682  }
683 
684 
693  template<class T>
694  inline void LitRK4Extrapolation(LitTrackParam<T>& par,
695  T zOut,
696  const LitFieldRegion<T>& field) {
697  static const T C1_2 = 0.5;
698  T zIn = par.Z;
699  T h = zOut - zIn;
700  LitFieldValue<T> v1, v2, v3;
701  field.GetFieldValue(zIn, v1);
702  field.GetFieldValue(zIn + C1_2 * h, v2);
703  field.GetFieldValue(zIn + h, v3);
704 
705  LitRK4Extrapolation(par, zOut, v1, v2, v3);
706  }
707 
708  } // namespace parallel
709 } // namespace lit
710 #endif /* LITEXTRAPOLATION_H_ */
lit::parallel::LitTrackParam::C2
T C2
Definition: LitTrackParam.h:100
LitFieldGrid.h
Class stores a grid of magnetic field values in XY slice at Z position.
LitMath.h
Useful math functions.
lit::parallel::LitTrackParam::Z
T Z
Definition: LitTrackParam.h:94
lit::parallel::LitTrackParam::C10
T C10
Definition: LitTrackParam.h:100
lit::parallel::LitFieldGrid::GetFieldValue
void GetFieldValue(fscal x, fscal y, LitFieldValue< fscal > &B) const
Return field value for (X, Y) position (scalar version).
Definition: LitFieldGrid.h:113
lit::parallel::LitTrackParam::C13
T C13
Definition: LitTrackParam.h:100
lit::parallel::LitTrackParam::X
T X
Definition: LitTrackParam.h:92
lit::parallel::LitTrackParam::C8
T C8
Definition: LitTrackParam.h:100
lit::parallel::LitTrackParam::Y
T Y
Definition: LitTrackParam.h:93
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
lit::parallel::LitFieldValue::Bx
T Bx
Definition: LitFieldValue.h:57
lit::parallel::LitTrackParam::C0
T C0
Definition: LitTrackParam.h:100
LitTrackParam.h
Track parameters data class.
lit::parallel::LitTrackParam::C9
T C9
Definition: LitTrackParam.h:100
lit::parallel::LitTrackParam::C3
T C3
Definition: LitTrackParam.h:100
lit::parallel::LitTrackParam::C14
T C14
Definition: LitTrackParam.h:100
lit::parallel::LitTrackParam::C4
T C4
Definition: LitTrackParam.h:100
lit::parallel::LitTrackParam
Track parameters data class.
Definition: LitTrackParam.h:34
lit::parallel::LitTrackParam::C12
T C12
Definition: LitTrackParam.h:100
lit::parallel::LitFieldGrid
Class stores a grid of magnetic field values in XY slice at Z position.
Definition: LitFieldGrid.h:44
lit::parallel::LitFieldRegion
Storage for field approximation along Z.
Definition: LitFieldRegion.h:26
lit::parallel::LitFieldValue
Magnetic field value at a certain point in the space.
Definition: LitFieldValue.h:29
h
Data class with information on a STS local track.
lit::parallel::LitTrackParam::Ty
T Ty
Definition: LitTrackParam.h:96
lit::parallel::LitTrackParam::Tx
T Tx
Definition: LitTrackParam.h:95
lit::parallel::LitRK4Extrapolation
void LitRK4Extrapolation(LitTrackParam< T > &par, T zOut, const LitFieldGrid &field1, const LitFieldGrid &field2, const LitFieldGrid &field3)
Definition: LitExtrapolation.h:73
lit::parallel::LitTrackParam::C5
T C5
Definition: LitTrackParam.h:100
lit::parallel::LitTrackParam::C1
T C1
Definition: LitTrackParam.h:100
lit::parallel::LitFieldValue::By
T By
Definition: LitFieldValue.h:57
lit::parallel::LitTrackParam::C11
T C11
Definition: LitTrackParam.h:100
lit::parallel::LitTrackParam::C7
T C7
Definition: LitTrackParam.h:100
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
LitFieldRegion.h
lit::parallel::LitTrackParam::C6
T C6
Definition: LitTrackParam.h:100
lit::parallel::LitTrackParam::Qp
T Qp
Definition: LitTrackParam.h:97
lit::parallel::rcp
fscal rcp(const fscal &a)
Returns reciprocal.
Definition: LitMath.h:26
lit::parallel::LitLineExtrapolation
void LitLineExtrapolation(LitTrackParam< T > &par, T zOut)
Line track extrapolation for the field free regions.
Definition: LitExtrapolation.h:37
lit::parallel::LitFieldRegion::GetFieldValue
void GetFieldValue(const T &z, LitFieldValue< T > &B) const
Returns field value at a certain Z position.
Definition: LitFieldRegion.h:121
NS_L1TrackFitter::ZERO
const fvec ZERO
Definition: L1TrackFitter.cxx:32
lit::parallel::LitFieldValue::Bz
T Bz
Definition: LitFieldValue.h:57
NS_L1TrackFitter::ONE
const fvec ONE
Definition: L1TrackFitter.cxx:33
lit
Definition: LitTrackFinderNNVecElectron.h:19