CbmRoot
L1TrackParFit.cxx
Go to the documentation of this file.
1 #include "L1TrackParFit.h"
2 
3 
4 #define cnst const fvec
5 
6 void L1TrackParFit::Filter(L1UMeasurementInfo& info, fvec u, fvec w) {
7  fvec wi, zeta, zetawi, HCH;
8  fvec F0, F1, F2, F3, F4, F5;
9  fvec K1, K2, K3, K4, K5;
10 
11  zeta = info.cos_phi * fx + info.sin_phi * fy - u;
12 
13  // F = CH'
14  F0 = info.cos_phi * C00 + info.sin_phi * C10;
15  F1 = info.cos_phi * C10 + info.sin_phi * C11;
16 
17  HCH = (F0 * info.cos_phi + F1 * info.sin_phi);
18 
19  F2 = info.cos_phi * C20 + info.sin_phi * C21;
20  F3 = info.cos_phi * C30 + info.sin_phi * C31;
21  F4 = info.cos_phi * C40 + info.sin_phi * C41;
22  F5 = info.cos_phi * C50 + info.sin_phi * C51;
23 
24 #if 0 // use mask
25  const fvec mask = (HCH < info.sigma2 * 16.);
26  wi = w/( (mask & info.sigma2) +HCH );
27  zetawi = zeta *wi;
28  chi2 += mask & (zeta * zetawi);
29 #else
30  wi = w / (info.sigma2 + HCH);
31  zetawi = zeta * wi;
32  chi2 += zeta * zetawi;
33 #endif // 0
34  NDF += w;
35 
36  K1 = F1 * wi;
37  K2 = F2 * wi;
38  K3 = F3 * wi;
39  K4 = F4 * wi;
40  K5 = F5 * wi;
41 
42  fx -= F0 * zetawi;
43  fy -= F1 * zetawi;
44  ftx -= F2 * zetawi;
45  fty -= F3 * zetawi;
46  fqp -= F4 * zetawi;
47  ft -= F5 * zetawi;
48 
49  C00 -= F0 * F0 * wi;
50  C10 -= K1 * F0;
51  C11 -= K1 * F1;
52  C20 -= K2 * F0;
53  C21 -= K2 * F1;
54  C22 -= K2 * F2;
55  C30 -= K3 * F0;
56  C31 -= K3 * F1;
57  C32 -= K3 * F2;
58  C33 -= K3 * F3;
59  C40 -= K4 * F0;
60  C41 -= K4 * F1;
61  C42 -= K4 * F2;
62  C43 -= K4 * F3;
63  C44 -= K4 * F4;
64  C50 -= K5 * F0;
65  C51 -= K5 * F1;
66  C52 -= K5 * F2;
67  C53 -= K5 * F3;
68  C54 -= K5 * F4;
69  C55 -= K5 * F5;
70 }
71 
72 void L1TrackParFit::FilterNoP(L1UMeasurementInfo& info, fvec u, fvec w) {
73  fvec wi, zeta, zetawi, HCH;
74  fvec F0, F1, F2, F3, F4, F5;
75  fvec K1, K2, K3, K4, K5;
76 
77  zeta = info.cos_phi * fx + info.sin_phi * fy - u;
78 
79  // F = CH'
80  F0 = info.cos_phi * C00 + info.sin_phi * C10;
81  F1 = info.cos_phi * C10 + info.sin_phi * C11;
82 
83  HCH = (F0 * info.cos_phi + F1 * info.sin_phi);
84 
85  F2 = info.cos_phi * C20 + info.sin_phi * C21;
86  F3 = info.cos_phi * C30 + info.sin_phi * C31;
87  F4 = info.cos_phi * C40 + info.sin_phi * C41;
88  F5 = info.cos_phi * C50 + info.sin_phi * C51;
89 
90 #if 0 // use mask
91  const fvec mask = (HCH < info.sigma2 * 16.);
92  wi = w/( (mask & info.sigma2) +HCH );
93  zetawi = zeta *wi;
94  chi2 += mask & (zeta * zetawi);
95 #else
96  wi = w / (info.sigma2 + HCH);
97  zetawi = zeta * wi;
98  chi2 += zeta * zetawi;
99 #endif // 0
100  NDF += w;
101 
102  K1 = F1 * wi;
103  K2 = F2 * wi;
104  K3 = F3 * wi;
105  K4 = F4 * wi;
106  K5 = F5 * wi;
107 
108  fx -= F0 * zetawi;
109  fy -= F1 * zetawi;
110  ftx -= F2 * zetawi;
111  fty -= F3 * zetawi;
112  // fqp -= F4*zetawi;
113  ft -= F5 * zetawi;
114 
115  C00 -= F0 * F0 * wi;
116  C10 -= K1 * F0;
117  C11 -= K1 * F1;
118  C20 -= K2 * F0;
119  C21 -= K2 * F1;
120  C22 -= K2 * F2;
121  C30 -= K3 * F0;
122  C31 -= K3 * F1;
123  C32 -= K3 * F2;
124  C33 -= K3 * F3;
125  // C40-= K4*F0;
126  // C41-= K4*F1;
127  // C42-= K4*F2;
128  // C43-= K4*F3;
129  // C44-= K4*F4;
130  C50 -= K5 * F0;
131  C51 -= K5 * F1;
132  C52 -= K5 * F2;
133  C53 -= K5 * F3;
134  C54 -= K5 * F4;
135  C55 -= K5 * F5;
136 }
137 
138 void L1TrackParFit::Filter(fvec t0, fvec dt0, fvec w) {
139  fvec wi, zeta, zetawi, HCH;
140  fvec F0, F1, F2, F3, F4, F5;
141  fvec K1, K2, K3, K4, K5;
142 
143  zeta = ft - t0;
144 
145  // F = CH'
146  F0 = C50;
147  F1 = C51;
148 
149  HCH = C55;
150 
151  F2 = C52;
152  F3 = C53;
153  F4 = C54;
154  F5 = C55;
155 
156 #if 0 // use mask
157  const fvec mask = (HCH < info.sigma2 * 16.);
158  wi = w/( (mask & info.sigma2) +HCH );
159  zetawi = zeta *wi;
160  chi2 += mask & (zeta * zetawi);
161 #else
162  wi = w / (dt0 * dt0 + HCH);
163  zetawi = zeta * wi;
164  chi2 += zeta * zetawi;
165 #endif // 0
166  NDF += w;
167 
168  K1 = F1 * wi;
169  K2 = F2 * wi;
170  K3 = F3 * wi;
171  K4 = F4 * wi;
172  K5 = F5 * wi;
173 
174  fx -= F0 * zetawi;
175  fy -= F1 * zetawi;
176  ftx -= F2 * zetawi;
177  fty -= F3 * zetawi;
178  fqp -= F4 * zetawi;
179  ft -= F5 * zetawi;
180 
181  C00 -= F0 * F0 * wi;
182  C10 -= K1 * F0;
183  C11 -= K1 * F1;
184  C20 -= K2 * F0;
185  C21 -= K2 * F1;
186  C22 -= K2 * F2;
187  C30 -= K3 * F0;
188  C31 -= K3 * F1;
189  C32 -= K3 * F2;
190  C33 -= K3 * F3;
191  C40 -= K4 * F0;
192  C41 -= K4 * F1;
193  C42 -= K4 * F2;
194  C43 -= K4 * F3;
195  C44 -= K4 * F4;
196  C50 -= K5 * F0;
197  C51 -= K5 * F1;
198  C52 -= K5 * F2;
199  C53 -= K5 * F3;
200  C54 -= K5 * F4;
201  C55 -= K5 * F5;
202 }
203 
204 void L1TrackParFit::ExtrapolateLine(fvec z_out, fvec* w) {
205 
206  cnst ZERO = 0.0, ONE = 1.;
207  cnst c_light = 29.9792458;
208 
209  fvec initialised = ZERO;
210  if (w) //TODO use operator {?:}
211  {
212  const fvec zero = ZERO;
213  initialised = fvec(zero < *w);
214  } else {
215  const fvec one = ONE;
216  const fvec zero = ZERO;
217  initialised = fvec(zero < one);
218  }
219 
220  fvec dz = (z_out - fz);
221 
222  fx += (ftx * dz & initialised);
223  fy += (fty * dz & initialised);
224  fz += (dz & initialised);
225  ft += ((dz * sqrt(1 + ftx * ftx + fty * fty) / c_light) & initialised);
226 
227  const fvec k1 = ftx * dz / (c_light * sqrt((ftx * ftx) + (fty * fty) + 1));
228  const fvec k2 = fty * dz / (c_light * sqrt((ftx * ftx) + (fty * fty) + 1));
229 
230  const fvec dzC32_in = dz * C32;
231 
232  C21 += (dzC32_in & initialised);
233  C10 += (dz * (C21 + C30) & initialised);
234 
235  const fvec C20_in = C20;
236 
237  C20 += ((dz * C22) & initialised);
238  C00 += (dz * (C20 + C20_in) & initialised);
239 
240  const fvec C31_in = C31;
241 
242  C31 += ((dz * C33) & initialised);
243  C11 += ((dz * (C31 + C31_in)) & initialised);
244  C30 += dzC32_in & initialised;
245 
246  C40 += (dz * C42 & initialised);
247  C41 += (dz * C43 & initialised);
248 
249  fvec c52 = C52;
250  fvec c53 = C53;
251 
252  C50 = ((k1 * C20 + k2 * C30) & initialised) + C50;
253  C51 = ((k1 * C21 + k2 * C31) & initialised) + C51;
254  C52 = ((k1 * C22 + k2 * C32) & initialised) + C52;
255  C53 = ((k1 * C32 + k2 * C33) & initialised) + C53;
256  C54 = ((k1 * C42 + k2 * C43) & initialised) + C54;
257  C55 = ((k1 * (C52 + c52) + k2 * (C53 + c53)) & initialised) + C55;
258 }
259 
260 void L1TrackParFit::ExtrapolateLine1(fvec z_out, fvec* w, fvec v) {
261 
262  cnst ZERO = 0.0, ONE = 1.;
263  cnst c_light = 29.9792458;
264 
265  fvec initialised = ZERO;
266  if (w) //TODO use operator {?:}
267  {
268  const fvec zero = ZERO;
269  initialised = fvec(zero < *w);
270  } else {
271  const fvec one = ONE;
272  const fvec zero = ZERO;
273  initialised = fvec(zero < one);
274  }
275 
276  fvec dz = (z_out - fz);
277 
278 
279  fx += (ftx * dz & initialised);
280  fy += (fty * dz & initialised);
281  fz += (dz & initialised);
282 
283 
284  ft += ((dz * sqrt(1 + ftx * ftx + fty * fty) / (v * c_light)) & initialised);
285 
286  const fvec k1 =
287  ftx * dz / ((v * c_light) * sqrt((ftx * ftx) + (fty * fty) + 1));
288  const fvec k2 =
289  fty * dz / ((v * c_light) * sqrt((ftx * ftx) + (fty * fty) + 1));
290 
291  const fvec dzC32_in = dz * C32;
292 
293  C21 += (dzC32_in & initialised);
294  C10 += (dz * (C21 + C30) & initialised);
295 
296  const fvec C20_in = C20;
297 
298  C20 += ((dz * C22) & initialised);
299  C00 += (dz * (C20 + C20_in) & initialised);
300 
301  const fvec C31_in = C31;
302 
303  C31 += ((dz * C33) & initialised);
304  C11 += ((dz * (C31 + C31_in)) & initialised);
305  C30 += dzC32_in & initialised;
306 
307  C40 += (dz * C42 & initialised);
308  C41 += (dz * C43 & initialised);
309 
310  fvec c52 = C52;
311  fvec c53 = C53;
312 
313  C50 = ((k1 * C20 + k2 * C30) & initialised) + C50;
314  C51 = ((k1 * C21 + k2 * C31) & initialised) + C51;
315  C52 = ((k1 * C22 + k2 * C32) & initialised) + C52;
316  C53 = ((k1 * C32 + k2 * C33) & initialised) + C53;
317  C54 = ((k1 * C42 + k2 * C43) & initialised) + C54;
318  C55 = ((k1 * (C52 + c52) + k2 * (C53 + c53)) & initialised) + C55;
319 
320  // fz = ( z_out & initialised) + ( (!initialised) & fz); //TEST
321  // cout << "fz = " << fz << endl;
322 }
323 
324 void L1TrackParFit::
325  Extrapolate // extrapolates track parameters and returns jacobian for extrapolation of CovMatrix
326  (fvec z_out, // extrapolate to this z position
327  fvec qp0, // use Q/p linearisation at this value
328  const L1FieldRegion& F,
329  fvec* w) {
330  //
331  // Forth-order Runge-Kutta method for solution of the equation
332  // of motion of a particle with parameter qp = Q /P
333  // in the magnetic field B()
334  //
335  // | x | tx
336  // | y | ty
337  // d/dz = | tx| = ft * ( ty * ( B(3)+tx*b(1) ) - (1+tx**2)*B(2) )
338  // | ty| ft * (-tx * ( B(3)+ty+b(2) - (1+ty**2)*B(1) ) ,
339  //
340  // where ft = c_light*qp*sqrt ( 1 + tx**2 + ty**2 ) .
341  //
342  // In the following for RK step :
343  //
344  // x=x[0], y=x[1], tx=x[3], ty=x[4].
345  //
346  //========================================================================
347  //
348  // NIM A395 (1997) 169-184; NIM A426 (1999) 268-282.
349  //
350  // the routine is based on LHC(b) utility code
351  //
352  //========================================================================
353 
354 
355  cnst ZERO = 0.0, ONE = 1.;
356  cnst c_light = 0.000299792458;
357 
358  const fvec a[4] = {0.0f, 0.5f, 0.5f, 1.0f};
359  const fvec c[4] = {1.0f / 6.0f, 1.0f / 3.0f, 1.0f / 3.0f, 1.0f / 6.0f};
360  const fvec b[4] = {0.0f, 0.5f, 0.5f, 1.0f};
361 
362  int step4;
363  fvec k[20], x0[5], x[5], k1[20];
364  fvec Ax[4], Ay[4], Ax_tx[4], Ay_tx[4], Ax_ty[4], Ay_ty[4], At[4], At_tx[4],
365  At_ty[4];
366 
367  //----------------------------------------------------------------
368 
369  fvec qp_in = fqp;
370  const fvec z_in = fz;
371  const fvec h = (z_out - fz);
372 
373  // cout<<h<<" h"<<endl;
374  // cout<<ftx<<" ftx"<<endl;
375  // cout<<fty<<" fty"<<endl;
376 
377  fvec hC = h * c_light;
378  x0[0] = fx;
379  x0[1] = fy;
380  x0[2] = ftx;
381  x0[3] = fty;
382  x0[4] = ft;
383  //
384  // Runge-Kutta step
385  //
386 
387  int step;
388  int i;
389 
390  fvec B[4][3];
391  for (step = 0; step < 4; ++step) {
392  F.Get(z_in + a[step] * h, B[step]);
393  }
394 
395  for (step = 0; step < 4; ++step) {
396  for (i = 0; i < 5; ++i) {
397  if (step == 0) {
398  x[i] = x0[i];
399  } else {
400  x[i] = x0[i] + b[step] * k[step * 5 - 5 + i];
401  }
402  }
403 
404  fvec tx = x[2];
405  fvec ty = x[3];
406  fvec tx2 = tx * tx;
407  fvec ty2 = ty * ty;
408  fvec txty = tx * ty;
409  fvec tx2ty21 = 1.f + tx2 + ty2;
410  // if( tx2ty21 > 1.e4 ) return 1;
411  fvec I_tx2ty21 = 1.f / tx2ty21 * qp0;
412  fvec tx2ty2 = sqrt(tx2ty21);
413  // fvec I_tx2ty2 = qp0 * hC / tx2ty2 ; unsused ???
414  tx2ty2 *= hC;
415  fvec tx2ty2qp = tx2ty2 * qp0;
416 
417  // cout<<B[step][0]<<" B["<<step<<"][0] "<<B[step][2]<<" B["<<step<<"][2] "<<B[step][1]<<" B["<<step<<"][1]"<<endl;
418  Ax[step] =
419  (txty * B[step][0] + ty * B[step][2] - (1.f + tx2) * B[step][1]) * tx2ty2;
420  Ay[step] = (-txty * B[step][1] - tx * B[step][2] + (1.f + ty2) * B[step][0])
421  * tx2ty2;
422 
423  Ax_tx[step] = Ax[step] * tx * I_tx2ty21
424  + (ty * B[step][0] - 2.f * tx * B[step][1]) * tx2ty2qp;
425  Ax_ty[step] =
426  Ax[step] * ty * I_tx2ty21 + (tx * B[step][0] + B[step][2]) * tx2ty2qp;
427  Ay_tx[step] =
428  Ay[step] * tx * I_tx2ty21 + (-ty * B[step][1] - B[step][2]) * tx2ty2qp;
429  Ay_ty[step] = Ay[step] * ty * I_tx2ty21
430  + (-tx * B[step][1] + 2.f * ty * B[step][0]) * tx2ty2qp;
431 
432  fvec p2 = 1.f / (qp0 * qp0);
433  fvec m2 = 0.1395679f * 0.1395679f;
434  fvec v = 29.9792458f * sqrt(p2 / (m2 + p2));
435  At[step] = sqrt(1.f + tx * tx + ty * ty) / v;
436  At_tx[step] = tx / sqrt(1.f + tx * tx + ty * ty) / v;
437  At_ty[step] = ty / sqrt(1.f + tx * tx + ty * ty) / v;
438 
439  // cout<<Ax[step]<<" Ax[step] "<<Ay[step]<<" ay "<<At[step]<<" At[step] "<<qp0<<" qp0 "<<h<<" h"<<endl;
440 
441  step4 = step * 5;
442  k[step4] = tx * h;
443  k[step4 + 1] = ty * h;
444  k[step4 + 2] = Ax[step] * qp0;
445  k[step4 + 3] = Ay[step] * qp0;
446  k[step4 + 4] = At[step] * h;
447  } // end of Runge-Kutta steps
448 
449  fvec initialised = ZERO;
450  if (w) //TODO use operator {?:}
451  {
452  const fvec zero = ZERO;
453  initialised = fvec(zero < *w);
454  } else {
455  const fvec one = ONE;
456  const fvec zero = ZERO;
457  initialised = fvec(zero < one);
458  }
459 
460  {
461 
462 
463  // cout<<x0[0]<<" x0[0] "<<c[0]<<" c 0 "<<k[0]<<" k0 "<<c[1]<<" c1 "<<k[5+0]<<" k5 "<<c[2]<<" c2 "<<k[10+0]<<" k10"<<c[3]<<" c3 "<<k[15+0]<<" k15"<<endl;
464 
465  // cout << "w = " << *w << "; ";
466  // cout << "initialised = " << initialised << "; ";
467  // cout << "fx = " << fx;
468 
469  fx = ((x0[0] + c[0] * k[0] + c[1] * k[5 + 0] + c[2] * k[10 + 0]
470  + c[3] * k[15 + 0])
471  & initialised)
472  + ((!initialised) & fx);
473  fy = ((x0[1] + c[0] * k[1] + c[1] * k[5 + 1] + c[2] * k[10 + 1]
474  + c[3] * k[15 + 1])
475  & initialised)
476  + ((!initialised) & fy);
477  ftx = ((x0[2] + c[0] * k[2] + c[1] * k[5 + 2] + c[2] * k[10 + 2]
478  + c[3] * k[15 + 2])
479  & initialised)
480  + ((!initialised) & ftx);
481  fty = ((x0[3] + c[0] * k[3] + c[1] * k[5 + 3] + c[2] * k[10 + 3]
482  + c[3] * k[15 + 3])
483  & initialised)
484  + ((!initialised) & fty);
485  ft = ((x0[4] + c[0] * k[4] + c[1] * k[5 + 4] + c[2] * k[10 + 4]
486  + c[3] * k[15 + 4])
487  & initialised)
488  + ((!initialised) & ft);
489  fz = (z_out & initialised) + ((!initialised) & fz);
490 
491  // cout << "; fx = " << fx << endl;
492  }
493  // cout<<fx<<" fx"<<endl;
494 
495  //
496  // Derivatives dx/dqp
497  //
498 
499  x0[0] = 0.f;
500  x0[1] = 0.f;
501  x0[2] = 0.f;
502  x0[3] = 0.f;
503  x0[4] = 0.f;
504 
505  //
506  // Runge-Kutta step for derivatives dx/dqp
507 
508  for (step = 0; step < 4; ++step) {
509  for (i = 0; i < 5; ++i) {
510  if (step == 0) {
511  x[i] = x0[i];
512  } else {
513  x[i] = x0[i] + b[step] * k1[step * 5 - 5 + i];
514  }
515  }
516  step4 = step * 5;
517  k1[step4] = x[2] * h;
518  k1[step4 + 1] = x[3] * h;
519  k1[step4 + 2] = Ax[step] + Ax_tx[step] * x[2] + Ax_ty[step] * x[3];
520  k1[step4 + 3] = Ay[step] + Ay_tx[step] * x[2] + Ay_ty[step] * x[3];
521  k1[step4 + 4] = At[step] + At_tx[step] * x[2] + At_ty[step] * x[3];
522 
523  } // end of Runge-Kutta steps for derivatives dx/dqp
524 
525  fvec J[36];
526 
527  for (i = 0; i < 4; ++i) {
528  J[24 + i] = x0[i] + c[0] * k1[i] + c[1] * k1[5 + i] + c[2] * k1[10 + i]
529  + c[3] * k1[15 + i];
530  }
531  J[28] = 1.;
532  J[29] = x0[4] + c[0] * k1[4] + c[1] * k1[5 + 4] + c[2] * k1[10 + 4]
533  + c[3] * k1[15 + 4];
534  //
535  // end of derivatives dx/dqp
536  //
537 
538  // Derivatives dx/tx
539  //
540 
541  x0[0] = 0.f;
542  x0[1] = 0.f;
543  x0[2] = 1.f;
544  x0[3] = 0.f;
545  x0[4] = 0.f;
546 
547  //
548  // Runge-Kutta step for derivatives dx/dtx
549  //
550 
551  for (step = 0; step < 4; ++step) {
552  for (i = 0; i < 5; ++i) {
553  if (step == 0) {
554  x[i] = x0[i];
555  } else if (i != 2) {
556  x[i] = x0[i] + b[step] * k1[step * 5 - 5 + i];
557  }
558  }
559  step4 = step * 5;
560  k1[step4] = x[2] * h;
561  k1[step4 + 1] = x[3] * h;
562  // k1[step4+2] = Ax_tx[step] * x[2] + Ax_ty[step] * x[3];
563  k1[step4 + 3] = Ay_tx[step] * x[2] + Ay_ty[step] * x[3];
564  k1[step4 + 4] = At_tx[step] * x[2] + At_ty[step] * x[3];
565  } // end of Runge-Kutta steps for derivatives dx/dtx
566 
567  for (i = 0; i < 4; ++i) {
568  if (i != 2) {
569  J[12 + i] = x0[i] + c[0] * k1[i] + c[1] * k1[5 + i] + c[2] * k1[10 + i]
570  + c[3] * k1[15 + i];
571  }
572  }
573  // end of derivatives dx/dtx
574  J[14] = 1.f;
575  J[16] = 0.f;
576  J[17] = x0[4] + c[0] * k1[4] + c[1] * k1[5 + 4] + c[2] * k1[10 + 4]
577  + c[3] * k1[15 + 4];
578 
579  // Derivatives dx/ty
580  //
581 
582  x0[0] = 0.f;
583  x0[1] = 0.f;
584  x0[2] = 0.f;
585  x0[3] = 1.f;
586  x0[4] = 0.f;
587 
588  //
589  // Runge-Kutta step for derivatives dx/dty
590  //
591 
592  for (step = 0; step < 4; ++step) {
593  for (i = 0; i < 5; ++i) {
594  if (step == 0) {
595  x[i] = x0[i]; // ty fixed
596  } else if (i != 3) {
597  x[i] = x0[i] + b[step] * k1[step * 5 - 5 + i];
598  }
599  }
600  step4 = step * 5;
601  k1[step4] = x[2] * h;
602  k1[step4 + 1] = x[3] * h;
603  k1[step4 + 2] = Ax_tx[step] * x[2] + Ax_ty[step] * x[3];
604  // k1[step4+3] = Ay_tx[step] * x[2] + Ay_ty[step] * x[3];
605  k1[step4 + 4] = At_tx[step] * x[2] + At_ty[step] * x[3];
606  } // end of Runge-Kutta steps for derivatives dx/dty
607 
608  for (i = 0; i < 3; ++i) {
609  J[18 + i] = x0[i] + c[0] * k1[i] + c[1] * k1[5 + i] + c[2] * k1[10 + i]
610  + c[3] * k1[15 + i];
611  }
612  // end of derivatives dx/dty
613  J[21] = 1.;
614  J[22] = 0.;
615  J[23] = x0[4] + c[0] * k1[4] + c[1] * k1[5 + 4] + c[2] * k1[10 + 4]
616  + c[3] * k1[15 + 4];
617  //
618  // derivatives dx/dx and dx/dy
619 
620  for (i = 0; i < 12; ++i)
621  J[i] = 0.;
622  J[0] = 1.;
623  J[7] = 1.;
624  for (i = 30; i < 35; i++)
625  J[i] = 0.f;
626  J[35] = 1.f;
627 
628  fvec dqp = qp_in - qp0;
629 
630  { // update parameters
631  fx += ((J[6 * 4 + 0] * dqp) & initialised);
632  fy += ((J[6 * 4 + 1] * dqp) & initialised);
633  ftx += ((J[6 * 4 + 2] * dqp) & initialised);
634  fty += ((J[6 * 4 + 3] * dqp) & initialised);
635  ft += ((J[6 * 4 + 5] * dqp) & initialised);
636  }
637  // cout<<fx<<" fx"<<endl;
638  // covariance matrix transport
639 
640  //cout<< (ft - ft_old)<<" ft dt "<<endl;
641 
642  // // if ( C_in&&C_out ) CbmKFMath::multQtSQ( 5, J, C_in, C_out); // TODO
643  // j(0,2) = J[5*2 + 0];
644  // j(1,2) = J[5*2 + 1];
645  // j(2,2) = J[5*2 + 2];
646  // j(3,2) = J[5*2 + 3];
647  //
648  // j(0,3) = J[5*3 + 0];
649  // j(1,3) = J[5*3 + 1];
650  // j(2,3) = J[5*3 + 2];
651  // j(3,3) = J[5*3 + 3];
652  //
653  // j(0,4) = J[5*4 + 0];
654  // j(1,4) = J[5*4 + 1];
655  // j(2,4) = J[5*4 + 2];
656  // j(3,4) = J[5*4 + 3];
657 
658  const fvec c42 = C42, c43 = C43;
659 
660  const fvec cj00 =
661  C00 + C20 * J[6 * 2 + 0] + C30 * J[6 * 3 + 0] + C40 * J[6 * 4 + 0];
662  const fvec cj10 =
663  C10 + C21 * J[6 * 2 + 0] + C31 * J[6 * 3 + 0] + C41 * J[6 * 4 + 0];
664  const fvec cj20 =
665  C20 + C22 * J[6 * 2 + 0] + C32 * J[6 * 3 + 0] + c42 * J[6 * 4 + 0];
666  const fvec cj30 =
667  C30 + C32 * J[6 * 2 + 0] + C33 * J[6 * 3 + 0] + c43 * J[6 * 4 + 0];
668  const fvec cj40 =
669  C40 + C42 * J[6 * 2 + 0] + C43 * J[6 * 3 + 0] + C44 * J[6 * 4 + 0];
670  const fvec cj50 =
671  C50 + C52 * J[6 * 2 + 0] + C53 * J[6 * 3 + 0] + C54 * J[6 * 4 + 0];
672 
673  // const fvec cj01 = C10 + C20*J[6*2 + 1] + C30*J[6*3 + 1] + C40*J[6*4 + 1];
674  const fvec cj11 =
675  C11 + C21 * J[6 * 2 + 1] + C31 * J[6 * 3 + 1] + C41 * J[6 * 4 + 1];
676  const fvec cj21 =
677  C21 + C22 * J[6 * 2 + 1] + C32 * J[6 * 3 + 1] + c42 * J[6 * 4 + 1];
678  const fvec cj31 =
679  C31 + C32 * J[6 * 2 + 1] + C33 * J[6 * 3 + 1] + c43 * J[6 * 4 + 1];
680  const fvec cj41 =
681  C41 + C42 * J[6 * 2 + 1] + C43 * J[6 * 3 + 1] + C44 * J[6 * 4 + 1];
682  const fvec cj51 =
683  C51 + C52 * J[6 * 2 + 1] + C53 * J[6 * 3 + 1] + C54 * J[6 * 4 + 1];
684 
685  // const fvec cj02 = C20*J[6*2 + 2] + C30*J[6*3 + 2] + C40*J[6*4 + 2];
686  // const fvec cj12 = C21*J[6*2 + 2] + C31*J[6*3 + 2] + C41*J[6*4 + 2];
687  const fvec cj22 =
688  C22 * J[6 * 2 + 2] + C32 * J[6 * 3 + 2] + c42 * J[6 * 4 + 2];
689  const fvec cj32 =
690  C32 * J[6 * 2 + 2] + C33 * J[6 * 3 + 2] + c43 * J[6 * 4 + 2];
691  const fvec cj42 =
692  C42 * J[6 * 2 + 2] + C43 * J[6 * 3 + 2] + C44 * J[6 * 4 + 2];
693  const fvec cj52 =
694  C52 * J[6 * 2 + 2] + C53 * J[6 * 3 + 2] + C54 * J[6 * 4 + 2];
695 
696  // const fvec cj03 = C20*J[6*2 + 3] + C30*J[6*3 + 3] + C40*J[6*4 + 3];
697  // const fvec cj13 = C21*J[6*2 + 3] + C31*J[6*3 + 3] + C41*J[6*4 + 3];
698  const fvec cj23 =
699  C22 * J[6 * 2 + 3] + C32 * J[6 * 3 + 3] + c42 * J[6 * 4 + 3];
700  const fvec cj33 =
701  C32 * J[6 * 2 + 3] + C33 * J[6 * 3 + 3] + c43 * J[6 * 4 + 3];
702  const fvec cj43 =
703  C42 * J[6 * 2 + 3] + C43 * J[6 * 3 + 3] + C44 * J[6 * 4 + 3];
704  const fvec cj53 =
705  C52 * J[6 * 2 + 3] + C53 * J[6 * 3 + 3] + C54 * J[6 * 4 + 3];
706 
707  const fvec cj24 = C42;
708  const fvec cj34 = C43;
709  const fvec cj44 = C44;
710  const fvec cj54 = C54;
711 
712  // const fvec cj05 = C50 + C20*J[17] + C30*J[23] + C40*J[29];
713  // const fvec cj15 = C51 + C21*J[17] + C31*J[23] + C41*J[29];
714  const fvec cj25 = C52 + C22 * J[17] + C32 * J[23] + C42 * J[29];
715  const fvec cj35 = C53 + C32 * J[17] + C33 * J[23] + C43 * J[29];
716  const fvec cj45 = C54 + C42 * J[17] + C43 * J[23] + C44 * J[29];
717  const fvec cj55 = C55 + C52 * J[17] + C53 * J[23] + C54 * J[29];
718 
719 
720  C00 = ((cj00 + cj20 * J[12] + cj30 * J[18] + cj40 * J[24]) & initialised)
721  + ((!initialised) & C00);
722 
723  C10 = ((cj10 + cj20 * J[13] + cj30 * J[19] + cj40 * J[25]) & initialised)
724  + ((!initialised) & C10);
725  C11 = ((cj11 + cj21 * J[13] + cj31 * J[19] + cj41 * J[25]) & initialised)
726  + ((!initialised) & C11);
727 
728  C20 = ((cj20 + cj30 * J[20] + cj40 * J[26]) & initialised)
729  + ((!initialised) & C20);
730  C21 = ((cj21 + cj31 * J[20] + cj41 * J[26]) & initialised)
731  + ((!initialised) & C21);
732  C22 = ((cj22 + cj32 * J[20] + cj42 * J[26]) & initialised)
733  + ((!initialised) & C22);
734 
735  C30 = ((cj30 + cj20 * J[15] + cj40 * J[27]) & initialised)
736  + ((!initialised) & C30);
737  C31 = ((cj31 + cj21 * J[15] + cj41 * J[27]) & initialised)
738  + ((!initialised) & C31);
739  C32 = ((cj32 + cj22 * J[15] + cj42 * J[27]) & initialised)
740  + ((!initialised) & C32);
741  C33 = ((cj33 + cj23 * J[15] + cj43 * J[27]) & initialised)
742  + ((!initialised) & C33);
743 
744  C40 = ((cj40) &initialised) + ((!initialised) & C40);
745  C41 = ((cj41) &initialised) + ((!initialised) & C41);
746  C42 = ((cj42) &initialised) + ((!initialised) & C42);
747  C43 = ((cj43) &initialised) + ((!initialised) & C43);
748  C44 = ((cj44) &initialised) + ((!initialised) & C44);
749 
750  C50 = ((cj50 + cj20 * J[17] + cj30 * J[23] + cj40 * J[29]) & initialised)
751  + ((!initialised) & C50);
752  C51 = ((cj51 + cj21 * J[17] + cj31 * J[23] + cj41 * J[29]) & initialised)
753  + ((!initialised) & C51);
754  C52 = ((cj52 + cj22 * J[17] + cj32 * J[23] + cj42 * J[29]) & initialised)
755  + ((!initialised) & C52);
756  C53 = ((cj53 + cj23 * J[17] + cj33 * J[23] + cj43 * J[29]) & initialised)
757  + ((!initialised) & C53);
758  C54 = ((cj54 + cj24 * J[17] + cj34 * J[23] + cj44 * J[29]) & initialised)
759  + ((!initialised) & C54);
760  C55 = ((cj55 + cj25 * J[17] + cj35 * J[23] + cj45 * J[29]) & initialised)
761  + ((!initialised) & C55);
762 }
763 
764 void L1TrackParFit::L1AddPipeMaterial(fvec qp0, fvec w, fvec mass2) {
765  cnst ONE = 1.f;
766 
767  // static const fscal RadThick=0.0009f;//0.5/18.76;
768  // static const fscal logRadThick=log(RadThick);
769  //const fscal RadThick=0.0009f;//0.5/18.76;
770 
771  const fscal logRadThick = log(PipeRadThick[0]);
772  fvec tx = ftx;
773  fvec ty = fty;
774  fvec txtx = ftx * ftx;
775  fvec tyty = fty * fty;
776  fvec txtx1 = txtx + ONE;
777  fvec h = txtx + tyty;
778  fvec t = sqrt(txtx1 + tyty);
779  fvec h2 = h * h;
780  fvec qp0t = qp0 * t;
781 
782  cnst c1 = 0.0136f, c2 = c1 * 0.038f, c3 = c2 * 0.5f, c4 = -c3 / 2.0f,
783  c5 = c3 / 3.0f, c6 = -c3 / 4.0f;
784  fvec s0 =
785  (c1 + c2 * fvec(logRadThick) + c3 * h + h2 * (c4 + c5 * h + c6 * h2))
786  * qp0t;
787  //fvec a = ( (ONE+mass2*qp0*qp0t)*RadThick*s0*s0 );
788  fvec a = ((t + mass2 * qp0 * qp0t) * PipeRadThick * s0 * s0);
789 
790  C22 += w * txtx1 * a;
791  C32 += w * tx * ty * a;
792  C33 += w * (ONE + tyty) * a;
793 }
794 
795 
796 void L1TrackParFit::L1AddMaterial(fvec radThick, fvec qp0, fvec w, fvec mass2) {
797  cnst ONE = 1.;
798 
799  fvec tx = ftx;
800  fvec ty = fty;
801  fvec txtx = tx * tx;
802  fvec tyty = ty * ty;
803  fvec txtx1 = txtx + ONE;
804  fvec h = txtx + tyty;
805  fvec t = sqrt(txtx1 + tyty);
806  fvec h2 = h * h;
807  fvec qp0t = qp0 * t;
808 
809  cnst c1 = 0.0136f, c2 = c1 * 0.038f, c3 = c2 * 0.5f, c4 = -c3 / 2.0f,
810  c5 = c3 / 3.0f, c6 = -c3 / 4.0f;
811 
812  fvec s0 =
813  (c1 + c2 * log(radThick) + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
814  //fvec a = ( (ONE+mass2*qp0*qp0t)*radThick*s0*s0 );
815  fvec a = ((t + mass2 * qp0 * qp0t) * radThick * s0 * s0);
816 
817  C22 += w * txtx1 * a;
818  C32 += w * tx * ty * a;
819  C33 += w * (ONE + tyty) * a;
820 }
821 
822 void L1TrackParFit::L1AddThickMaterial(fvec radThick,
823  fvec qp0,
824  fvec w,
825  fvec mass2,
826  fvec thickness,
827  bool fDownstream) {
828  cnst ONE = 1.;
829 
830  fvec tx = ftx;
831  fvec ty = fty;
832  fvec txtx = tx * tx;
833  fvec tyty = ty * ty;
834  fvec txtx1 = txtx + ONE;
835  fvec h = txtx + tyty;
836  fvec t = sqrt(txtx1 + tyty);
837  fvec h2 = h * h;
838  fvec qp0t = qp0 * t;
839 
840  cnst c1 = 0.0136f, c2 = c1 * 0.038f, c3 = c2 * 0.5f, c4 = -c3 / 2.0f,
841  c5 = c3 / 3.0f, c6 = -c3 / 4.0f;
842 
843  fvec s0 =
844  (c1 + c2 * log(radThick) + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
845  //fvec a = ( (ONE+mass2*qp0*qp0t)*radThick*s0*s0 );
846  fvec a = ((t + mass2 * qp0 * qp0t) * radThick * s0 * s0);
847 
848  fvec D = (fDownstream) ? 1. : -1.;
849  fvec T23 = (thickness * thickness) / 3.0;
850  fvec T2 = thickness / 2.0;
851 
852  C00 += w * txtx1 * a * T23;
853  C10 += w * tx * ty * a * T23;
854  C20 += w * txtx1 * a * D * T2;
855  C30 += w * tx * ty * a * D * T2;
856 
857  C11 += w * (ONE + tyty) * a * T23;
858  C21 += w * tx * ty * a * D * T2;
859  C31 += w * (ONE + tyty) * a * D * T2;
860 
861  C22 += w * txtx1 * a;
862  C32 += w * tx * ty * a;
863  C33 += w * (ONE + tyty) * a;
864 }
865 
866 
867 void L1TrackParFit::L1AddMaterial(L1MaterialInfo& info,
868  fvec qp0,
869  fvec w,
870  fvec mass2) {
871  cnst ONE = 1.f;
872 
873  fvec tx = ftx;
874  fvec ty = fty;
875  // fvec time = ft;
876  fvec txtx = tx * tx;
877  fvec tyty = ty * ty;
878  fvec txtx1 = txtx + ONE;
879  fvec h = txtx + tyty;
880  fvec t = sqrt(txtx1 + tyty);
881  fvec h2 = h * h;
882  fvec qp0t = qp0 * t;
883 
884  cnst c1 = 0.0136f, c2 = c1 * 0.038f, c3 = c2 * 0.5f, c4 = -c3 / 2.0f,
885  c5 = c3 / 3.0f, c6 = -c3 / 4.0f;
886 
887  fvec s0 =
888  (c1 + c2 * info.logRadThick + c3 * h + h2 * (c4 + c5 * h + c6 * h2)) * qp0t;
889  //fvec a = ( (ONE+mass2*qp0*qp0t)*info.RadThick*s0*s0 );
890  fvec a = ((t + mass2 * qp0 * qp0t) * info.RadThick * s0 * s0);
891 
892  C22 += w * txtx1 * a;
893  C32 += w * tx * ty * a;
894  C33 += w * (ONE + tyty) * a;
895 }
896 
897 
898 void L1TrackParFit::EnergyLossCorrection(const fvec& mass2,
899  const fvec& radThick,
900  fvec& qp0,
901  fvec direction,
902  fvec w) {
903  const fvec& p2 = 1.f / (fqp * fqp);
904  const fvec& E2 = mass2 + p2;
905 
906  const fvec& bethe = ApproximateBetheBloch(p2 / mass2);
907 
908  fvec tr = sqrt(fvec(1.f) + ftx * ftx + fty * fty);
909 
910  const fvec& dE = bethe * radThick * tr * 2.33f * 9.34961f;
911 
912  const fvec& E2Corrected =
913  (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
914  fvec corr = sqrt(p2 / (E2Corrected - mass2));
915  fvec init = fvec(!(corr == corr)) | fvec(w < 1);
916  corr = fvec(fvec(1.f) & init) + fvec(corr & fvec(!(init)));
917 
918  qp0 *= corr;
919  fqp *= corr;
920  C40 *= corr;
921  C41 *= corr;
922  C42 *= corr;
923  C43 *= corr;
924  // C54 *= corr;
925  C44 *= corr * corr;
926 }
927 
928 void L1TrackParFit::EnergyLossCorrectionIron(const fvec& mass2,
929  const fvec& radThick,
930  fvec& qp0,
931  fvec direction,
932  fvec w) {
933  const fvec& p2 = 1.f / (qp0 * qp0);
934  const fvec& E2 = mass2 + p2;
935 
936  int atomicZ = 26;
937  float atomicA = 55.845f;
938  float rho = 7.87;
939  float radLen = 1.758f;
940 
941  fvec i;
942  if (atomicZ < 13)
943  i = (12. * atomicZ + 7.) * 1.e-9;
944  else
945  i = (9.76 * atomicZ + 58.8 * std::pow(atomicZ, -0.19)) * 1.e-9;
946 
947  const fvec& bethe =
948  ApproximateBetheBloch(p2 / mass2, rho, 0.20, 3.00, i, atomicZ / atomicA);
949 
950  fvec tr = sqrt(fvec(1.f) + ftx * ftx + fty * fty);
951 
952 
953  fvec dE = bethe * radThick * tr * radLen * rho;
954 
955  const fvec& E2Corrected =
956  (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
957  fvec corr = sqrt(p2 / (E2Corrected - mass2));
958  fvec init = fvec(!(corr == corr)) | fvec(w < 1);
959  corr = fvec(fvec(1.f) & init) + fvec(corr & fvec(!(init)));
960 
961  qp0 *= corr;
962  fqp *= corr;
963 
964  float P = fabs(1. / qp0[0]); // GeV
965 
966  float Z = atomicZ;
967  float A = atomicA;
968  float RHO = rho;
969 
970  fvec STEP = radThick * tr * radLen;
971  fvec EMASS = 0.511 * 1e-3; // GeV
972 
973  fvec BETA = P / sqrt(E2Corrected);
974  fvec GAMMA = sqrt(E2Corrected) / sqrt(mass2);
975 
976  // Calculate xi factor (KeV).
977  fvec XI = (153.5 * Z * STEP * RHO) / (A * BETA * BETA);
978 
979  // Maximum energy transfer to atomic electron (KeV).
980  fvec ETA = BETA * GAMMA;
981  fvec ETASQ = ETA * ETA;
982  fvec RATIO = EMASS / sqrt(mass2);
983  fvec F1 = 2. * EMASS * ETASQ;
984  fvec F2 = 1. + 2. * RATIO * GAMMA + RATIO * RATIO;
985  fvec EMAX = 1e6 * F1 / F2;
986 
987  fvec DEDX2 = XI * EMAX * (1. - (BETA * BETA / 2.)) * 1e-12;
988 
989  fvec SDEDX = ((E2) *DEDX2) / std::pow(P, 6);
990 
991  // T.C40 *= corr;
992  // T.C41 *= corr;
993  // T.C42 *= corr;
994  // T.C43 *= corr;
995  // T.C44 *= corr*corr;
996  C44 += fabs(SDEDX);
997 }
998 
999 void L1TrackParFit::EnergyLossCorrectionCarbon(const fvec& mass2,
1000  const fvec& radThick,
1001  fvec& qp0,
1002  fvec direction,
1003  fvec w) {
1004  const fvec& p2 = 1.f / (qp0 * qp0);
1005  const fvec& E2 = mass2 + p2;
1006 
1007  int atomicZ = 6;
1008  float atomicA = 12.011f;
1009  float rho = 2.265;
1010  float radLen = 18.76f;
1011 
1012  fvec i;
1013  if (atomicZ < 13)
1014  i = (12. * atomicZ + 7.) * 1.e-9;
1015  else
1016  i = (9.76 * atomicZ + 58.8 * std::pow(atomicZ, -0.19)) * 1.e-9;
1017 
1018  const fvec& bethe =
1019  ApproximateBetheBloch(p2 / mass2, rho, 0.20, 3.00, i, atomicZ / atomicA);
1020 
1021  fvec tr = sqrt(fvec(1.f) + ftx * ftx + fty * fty);
1022 
1023 
1024  fvec dE = bethe * radThick * tr * radLen * rho;
1025 
1026  const fvec& E2Corrected =
1027  (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
1028  fvec corr = sqrt(p2 / (E2Corrected - mass2));
1029  fvec init = fvec(!(corr == corr)) | fvec(w < 1);
1030  corr = fvec(fvec(1.f) & init) + fvec(corr & fvec(!(init)));
1031 
1032  qp0 *= corr;
1033  fqp *= corr;
1034 
1035  float P = fabs(1. / qp0[0]); // GeV
1036 
1037  float Z = atomicZ;
1038  float A = atomicA;
1039  float RHO = rho;
1040 
1041  fvec STEP = radThick * tr * radLen;
1042  fvec EMASS = 0.511 * 1e-3; // GeV
1043 
1044  fvec BETA = P / sqrt(E2Corrected);
1045  fvec GAMMA = sqrt(E2Corrected) / sqrt(mass2);
1046 
1047  // Calculate xi factor (KeV).
1048  fvec XI = (153.5 * Z * STEP * RHO) / (A * BETA * BETA);
1049 
1050  // Maximum energy transfer to atomic electron (KeV).
1051  fvec ETA = BETA * GAMMA;
1052  fvec ETASQ = ETA * ETA;
1053  fvec RATIO = EMASS / sqrt(mass2);
1054  fvec F1 = 2. * EMASS * ETASQ;
1055  fvec F2 = 1. + 2. * RATIO * GAMMA + RATIO * RATIO;
1056  fvec EMAX = 1e6 * F1 / F2;
1057 
1058  fvec DEDX2 = XI * EMAX * (1. - (BETA * BETA / 2.)) * 1e-12;
1059 
1060  fvec SDEDX = ((E2) *DEDX2) / std::pow(P, 6);
1061 
1062  // T.C40 *= corr;
1063  // T.C41 *= corr;
1064  // T.C42 *= corr;
1065  // T.C43 *= corr;
1066  // T.C44 *= corr*corr;
1067  C44 += fabs(SDEDX);
1068 }
1069 
1070 void L1TrackParFit::EnergyLossCorrectionAl(const fvec& mass2,
1071  const fvec& radThick,
1072  fvec& qp0,
1073  fvec direction,
1074  fvec w) {
1075  const fvec& p2 = 1.f / (qp0 * qp0);
1076  const fvec& E2 = mass2 + p2;
1077 
1078  int atomicZ = 13;
1079  float atomicA = 26.981f;
1080  float rho = 2.70f;
1081  float radLen = 2.265f;
1082 
1083  fvec i;
1084  if (atomicZ < 13)
1085  i = (12. * atomicZ + 7.) * 1.e-9;
1086  else
1087  i = (9.76 * atomicZ + 58.8 * std::pow(atomicZ, -0.19)) * 1.e-9;
1088 
1089  const fvec& bethe =
1090  ApproximateBetheBloch(p2 / mass2, rho, 0.20, 3.00, i, atomicZ / atomicA);
1091 
1092  fvec tr = sqrt(fvec(1.f) + ftx * ftx + fty * fty);
1093 
1094 
1095  fvec dE = bethe * radThick * tr * radLen * rho;
1096 
1097  const fvec& E2Corrected =
1098  (sqrt(E2) + direction * dE) * (sqrt(E2) + direction * dE);
1099  fvec corr = sqrt(p2 / (E2Corrected - mass2));
1100  fvec init = fvec(!(corr == corr)) | fvec(w < 1);
1101  corr = fvec(fvec(1.f) & init) + fvec(corr & fvec(!(init)));
1102 
1103  qp0 *= corr;
1104  fqp *= corr;
1105 
1106  float P = fabs(1. / qp0[0]); // GeV
1107 
1108  float Z = atomicZ;
1109  float A = atomicA;
1110  float RHO = rho;
1111 
1112  fvec STEP = radThick * tr * radLen;
1113  fvec EMASS = 0.511 * 1e-3; // GeV
1114 
1115  fvec BETA = P / sqrt(E2Corrected);
1116  fvec GAMMA = sqrt(E2Corrected) / sqrt(mass2);
1117 
1118  // Calculate xi factor (KeV).
1119  fvec XI = (153.5 * Z * STEP * RHO) / (A * BETA * BETA);
1120 
1121  // Maximum energy transfer to atomic electron (KeV).
1122  fvec ETA = BETA * GAMMA;
1123  fvec ETASQ = ETA * ETA;
1124  fvec RATIO = EMASS / sqrt(mass2);
1125  fvec F1 = 2. * EMASS * ETASQ;
1126  fvec F2 = 1. + 2. * RATIO * GAMMA + RATIO * RATIO;
1127  fvec EMAX = 1e6 * F1 / F2;
1128 
1129  fvec DEDX2 = XI * EMAX * (1. - (BETA * BETA / 2.)) * 1e-12;
1130 
1131  fvec SDEDX = ((E2) *DEDX2) / std::pow(P, 6);
1132 
1133  // T.C40 *= corr;
1134  // T.C41 *= corr;
1135  // T.C42 *= corr;
1136  // T.C43 *= corr;
1137  // T.C44 *= corr*corr;
1138  C44 += fabs(SDEDX);
1139 }
1140 
1141 #undef cnst
L1MaterialInfo::RadThick
fvec RadThick
Definition: L1MaterialInfo.h:11
fscal
float fscal
Definition: L1/vectors/P4_F32vec4.h:250
h
Generates beam ions for transport simulation.
Definition: CbmBeamGenerator.h:17
L1TrackParFit::L1AddPipeMaterial
void L1AddPipeMaterial(fvec qp0, fvec w=1, fvec mass2=0.1395679f *0.1395679f)
Definition: L1TrackParFit.cxx:764
L1TrackParFit::C31
fvec C31
Definition: L1TrackParFit.h:15
f
float f
Definition: L1/vectors/P4_F32vec4.h:24
L1TrackParFit::C43
fvec C43
Definition: L1TrackParFit.h:16
L1TrackParFit::C41
fvec C41
Definition: L1TrackParFit.h:16
F32vec4
Definition: L1/vectors/P4_F32vec4.h:47
L1TrackParFit::ft
fvec ft
Definition: L1TrackParFit.h:15
L1TrackParFit::C50
fvec C50
Definition: L1TrackParFit.h:16
cnst
#define cnst
Definition: L1TrackParFit.cxx:4
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
L1TrackParFit::chi2
fvec chi2
Definition: L1TrackParFit.h:16
fvec
F32vec4 fvec
Definition: L1/vectors/P4_F32vec4.h:249
L1TrackParFit::fy
fvec fy
Definition: L1TrackParFit.h:15
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
L1TrackParFit::L1AddThickMaterial
void L1AddThickMaterial(fvec radThick, fvec qp0, fvec w=1, fvec mass2=0.1395679f *0.1395679f, fvec thickness=0, bool fDownstream=1)
Definition: L1TrackParFit.cxx:822
L1MaterialInfo
Definition: L1MaterialInfo.h:7
L1TrackParFit::C30
fvec C30
Definition: L1TrackParFit.h:15
L1TrackParFit::EnergyLossCorrectionAl
void EnergyLossCorrectionAl(const fvec &mass2, const fvec &radThick, fvec &qp0, fvec direction, fvec w=1)
Definition: L1TrackParFit.cxx:1070
L1TrackParFit::C10
fvec C10
Definition: L1TrackParFit.h:15
L1UMeasurementInfo::sigma2
fvec sigma2
Definition: L1UMeasurementInfo.h:12
L1TrackParFit::L1AddMaterial
void L1AddMaterial(L1MaterialInfo &info, fvec qp0, fvec w=1, fvec mass2=0.1395679f *0.1395679f)
Definition: L1TrackParFit.cxx:867
L1TrackParFit::C33
fvec C33
Definition: L1TrackParFit.h:16
L1MaterialInfo::logRadThick
fvec logRadThick
Definition: L1MaterialInfo.h:11
L1UMeasurementInfo::sin_phi
fvec sin_phi
Definition: L1UMeasurementInfo.h:12
L1TrackParFit::C42
fvec C42
Definition: L1TrackParFit.h:16
L1TrackParFit::C54
fvec C54
Definition: L1TrackParFit.h:16
L1TrackParFit::C00
fvec C00
Definition: L1TrackParFit.h:15
L1TrackParFit::C40
fvec C40
Definition: L1TrackParFit.h:16
L1FieldRegion
Definition: L1Field.h:85
L1TrackParFit::C20
fvec C20
Definition: L1TrackParFit.h:15
ApproximateBetheBloch
fvec ApproximateBetheBloch(const fvec &bg2)
Definition: L1AddMaterial.h:120
L1TrackParFit.h
L1TrackParFit::fx
fvec fx
Definition: L1TrackParFit.h:15
log
friend F32vec4 log(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:135
L1TrackParFit::C52
fvec C52
Definition: L1TrackParFit.h:16
L1TrackParFit::Filter
void Filter(L1UMeasurementInfo &info, fvec u, fvec w=1.)
Definition: L1TrackParFit.cxx:6
L1TrackParFit::C51
fvec C51
Definition: L1TrackParFit.h:16
L1TrackParFit::C11
fvec C11
Definition: L1TrackParFit.h:15
L1TrackParFit::ftx
fvec ftx
Definition: L1TrackParFit.h:15
L1TrackParFit::ExtrapolateLine
void ExtrapolateLine(fvec z_out, fvec *w=0)
Definition: L1TrackParFit.cxx:204
L1TrackParFit::FilterNoP
void FilterNoP(L1UMeasurementInfo &info, fvec u, fvec w=1.)
Definition: L1TrackParFit.cxx:72
L1UMeasurementInfo
Definition: L1UMeasurementInfo.h:7
L1TrackParFit::C21
fvec C21
Definition: L1TrackParFit.h:15
L1TrackParFit::NDF
fvec NDF
Definition: L1TrackParFit.h:16
L1TrackParFit::C32
fvec C32
Definition: L1TrackParFit.h:16
v
__m128 v
Definition: L1/vectors/P4_F32vec4.h:1
L1TrackParFit::C53
fvec C53
Definition: L1TrackParFit.h:16
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:60
L1TrackParFit::fqp
fvec fqp
Definition: L1TrackParFit.h:15
L1TrackParFit::C55
fvec C55
Definition: L1TrackParFit.h:16
L1TrackParFit::C44
fvec C44
Definition: L1TrackParFit.h:16
L1UMeasurementInfo::cos_phi
fvec cos_phi
Definition: L1UMeasurementInfo.h:12
L1TrackParFit::EnergyLossCorrectionIron
void EnergyLossCorrectionIron(const fvec &mass2, const fvec &radThick, fvec &qp0, fvec direction, fvec w=1)
Definition: L1TrackParFit.cxx:928
PipeRadThick
const fvec PipeRadThick
Definition: L1AddMaterial.h:11
NS_L1TrackFitter::ZERO
const fvec ZERO
Definition: L1TrackFitter.cxx:32
L1TrackParFit::ExtrapolateLine1
void ExtrapolateLine1(fvec z_out, fvec *w=0, fvec v=0)
Definition: L1TrackParFit.cxx:260
NS_L1TrackFitter::ONE
const fvec ONE
Definition: L1TrackFitter.cxx:33
L1TrackParFit::fz
fvec fz
Definition: L1TrackParFit.h:15
L1TrackParFit::EnergyLossCorrectionCarbon
void EnergyLossCorrectionCarbon(const fvec &mass2, const fvec &radThick, fvec &qp0, fvec direction, fvec w=1)
Definition: L1TrackParFit.cxx:999
L1TrackParFit::EnergyLossCorrection
void EnergyLossCorrection(const fvec &mass2, const fvec &radThick, fvec &qp0, fvec direction, fvec w=1)
Definition: L1TrackParFit.cxx:898
L1TrackParFit::fty
fvec fty
Definition: L1TrackParFit.h:15
NS_L1TrackFitter::c_light
const fvec c_light
Definition: L1TrackFitter.cxx:31
L1TrackParFit::C22
fvec C22
Definition: L1TrackParFit.h:15