BmnRoot
Loading...
Searching...
No Matches
CbmKFMath.cxx
Go to the documentation of this file.
1#include "CbmKFMath.h"
2
3#include "FairField.h"
4#include "FairTrackParam.h"
5
6#include <cmath>
7
8using std::sqrt;
9using std::exp;
10using std::log;
11using std::fabs;
12
13Bool_t CbmKFMath::intersectCone( Double_t zCone, Double_t ZCone, Double_t rCone, Double_t RCone ,
14 const Double_t x[],
15 Double_t *z1, Double_t *z2
16 )
17{
18 Double_t t = (RCone-rCone)/(ZCone-zCone);
19 Double_t A = (RCone*zCone-rCone*ZCone)/(ZCone-zCone);
20 Double_t x0 = x[0]-x[5]*x[2];
21 Double_t y0 = x[1]-x[5]*x[3];
22 Double_t tx = x[2];
23 Double_t ty = x[3];
24 Double_t a = tx*tx+ty*ty-t*t;
25 Double_t b = 2*(tx*x0+ty*y0 + t*A);
26 Double_t c = x0*x0+y0*y0-A*A;
27
28 if ( fabs(a)<1.E-4 ) return 1;
29 Double_t D = b*b - 4*a*c;
30 if ( D<0. ) return 1;
31 D = sqrt(D);
32 *z1 = (-b + D )/(2*a);
33 *z2 = (-b - D )/(2*a);
34 return 0;
35}
36
37
38void CbmKFMath::multQSQt( Int_t N, const Double_t Q[], const Double_t S[], Double_t S_out[] )
39 {
40 Double_t A[N*N];
41
42 {for( Int_t i=0, n=0; i<N; i++ )
43 {
44 for( Int_t j=0; j<N; j++, ++n )
45 {
46 A[n] = 0 ;
47 for( Int_t k=0; k<N; ++k ) A[n]+= S[indexS(i,k)] * Q[ j*N+k];
48 }
49 }}
50
51 {for( Int_t i=0; i<N; i++ )
52 {
53 for( Int_t j=0; j<=i; j++ )
54 {
55 Int_t n = indexS(i,j);
56 S_out[n] = 0 ;
57 for( Int_t k=0; k<N; k++ ) S_out[n] += Q[ i*N+k ] * A[ k*N+j ];
58 }
59 }}
60 }
61
62void CbmKFMath::multQtSQ( Int_t N, const Double_t Q[], const Double_t S[], Double_t S_out[] )
63{
64 Double_t A[N*N];
65
66 for( Int_t i=0, n=0; i<N; i++ ){
67 for( Int_t j=0; j<N; j++, ++n ){
68 A[n] = 0 ;
69 for( Int_t k=0; k<N; ++k ) A[n]+= S[indexS(i,k)] * Q[ k*N+j];
70 }
71 }
72
73 for( Int_t i=0; i<N; i++ ){
74 for( Int_t j=0; j<=i; j++ ){
75 Int_t n = indexS(i,j);
76 S_out[n] = 0 ;
77 for( Int_t k=0; k<N; k++ ) S_out[n] += Q[ k*N+i ] * A[ k*N+j ];
78 }
79 }
80}
81
82void CbmKFMath::multSSQ( const Double_t *A, const Double_t *B, Double_t *C, Int_t n )
83 {
84 for( Int_t i=0; i<n; ++i )
85 {
86 for( Int_t j=0; j<n; ++j )
87 {
88 Int_t ind = i*n+j;
89 C[ind] = 0;
90 for( Int_t k=0; k<n; ++k ) C[ind] += A[indexS(i,k)] * B[indexS(k,j)];
91 }
92 }
93 }
94
95void CbmKFMath::four_dim_inv( Double_t a[4][4] )
96 {
97 /**** Gaussian algorithm for 4x4 matrix inversion ****/
98 Int_t i,j,k,l;
99 Double_t factor;
100 Double_t temp[4];
101 Double_t b[4][4];
102
103 // Set b to I
104 for(i=0;i<4;i++) for(j=0;j<4;j++)
105 if(i==j) b[i][j]=1.0; else b[i][j]=0.0;
106
107 for(i=0;i<4;i++){
108 for(j=i+1;j<4;j++)
109 if (fabs(a[i][i])<fabs(a[j][i])){
110 for(l=0;l<4;l++) temp[l]=a[i][l];
111 for(l=0;l<4;l++) a[i][l]=a[j][l];
112 for(l=0;l<4;l++) a[j][l]=temp[l];
113 for(l=0;l<4;l++) temp[l]=b[i][l];
114 for(l=0;l<4;l++) b[i][l]=b[j][l];
115 for(l=0;l<4;l++) b[j][l]=temp[l];
116 }
117 factor=a[i][i];
118 for(j=4-1;j>-1;j--) {
119 b[i][j]/=factor;a[i][j]/=factor;
120 }
121 for(j=i+1;j<4;j++) {
122 factor=-a[j][i];
123 for(k=0;k<4;k++){
124 a[j][k]+=a[i][k]*factor;b[j][k]+=b[i][k]*factor;
125 }
126 }
127 }
128 for(i=3;i>0;i--){
129 for(j=i-1;j>-1;j--){
130 factor=-a[j][i];
131 for(k=0;k<4;k++){
132 a[j][k]+=a[i][k]*factor;b[j][k]+=b[i][k]*factor;
133 }
134 }
135 }
136
137 // copy b to a and return
138 for(i=0;i<4;i++) for(j=0;j<4;j++) a[i][j]=b[i][j];
139 }
140
141void CbmKFMath::five_dim_inv(Double_t a[5][5])
142{
143 /**** Gaussian algorithm for 5x5 matrix inversion ****/
144 Int_t i,j,k,l;
145 Double_t factor;
146 Double_t temp[5];
147 Double_t b[5][5];
148
149 // Set b to I
150 for(i=0;i<5;i++)
151 for(j=0;j<5;j++)
152 if(i==j) b[i][j]=1.0; else b[i][j]=0.0;
153
154 for(i=0;i<5;i++){
155 for(j=i+1;j<5;j++)
156 if (fabs(a[i][i])<fabs(a[j][i])){
157 for(l=0;l<5;l++) temp[l]=a[i][l];
158 for(l=0;l<5;l++) a[i][l]=a[j][l];
159 for(l=0;l<5;l++) a[j][l]=temp[l];
160 for(l=0;l<5;l++) temp[l]=b[i][l];
161 for(l=0;l<5;l++) b[i][l]=b[j][l];
162 for(l=0;l<5;l++) b[j][l]=temp[l];
163 }
164 factor=a[i][i];
165 // cout<<"Highest element "<<a[i][i]<<endl;
166 for(j=5-1;j>-1;j--){
167 b[i][j]/=factor;
168 a[i][j]/=factor;
169 }
170 for(j=i+1;j<5;j++) {
171 factor=-a[j][i];
172 for(k=0;k<5;k++){
173 a[j][k]+=a[i][k]*factor;
174 b[j][k]+=b[i][k]*factor;
175 }
176 }
177 }
178 for(i=4;i>0;i--){
179 for(j=i-1;j>-1;j--){
180 factor=-a[j][i];
181 for(k=0;k<5;k++){
182 a[j][k]+=a[i][k]*factor;
183 b[j][k]+=b[i][k]*factor;
184 }
185 }
186 }
187
188 // copy b to a and return
189 for(i=0;i<5;i++) for(j=0;j<5;j++) a[i][j]=b[i][j];
190 }
191
192Bool_t CbmKFMath::invS( Double_t A[], Int_t N )
193{
194
195 Bool_t ret = 0;
196
197 const Double_t ZERO = 1.E-20;
198
199 // input: simmetric > 0 NxN matrix A = {a11,a21,a22,a31..a33,..}
200 // output: inverse A, in case of problems fill zero and return 1
201
202 // A->low triangular Anew : A = Anew x Anew^T
203 // method:
204 // for(j=1,N) for(i=j,N) Aij=(Aii-sum_{k=1}^{j-1}Aik*Ajk )/Ajj
205 //
206
207 {
208 Double_t *j1 = A, *jj = A;
209 for( Int_t j=1; j<=N; j1+=j++, jj+=j ){
210 Double_t *ik = j1, x = 0;
211 while( ik!=jj ){
212 x -= (*ik) * (*ik);
213 ik++;
214 }
215 x += *ik;
216 if( x > ZERO ){
217 x = sqrt(x);
218 *ik = x;
219 ik++;
220 x = 1 / x;
221 for( Int_t step=1; step<=N-j; ik+=++step ){ // ik==Ai1
222 Double_t sum = 0;
223 for( Double_t *jk=j1; jk!=jj; sum += *(jk++) * *(ik++) );
224 *ik = (*ik - sum) * x; // ik == Aij
225 }
226 }else{
227 Double_t *ji=jj;
228 for( Int_t i=j; i<N; i++ ) *(ji+=i) = 0.;
229 ret = 1;
230 }
231 }
232 }
233
234 // A -> Ainv
235 // method :
236 // for(i=1,N){
237 // Aii = 1/Aii;
238 // for(j=1,i-1) Aij=-(sum_{k=j}^{i-1} Aik * Akj) / Aii ;
239 // }
240
241 {
242 Double_t *ii=A,*ij=A;
243 for( Int_t i = 1; i<=N; ij=ii+1, ii+=++i ){
244 if( *ii > ZERO ){
245 Double_t x = -(*ii = 1./ *ii);
246 {
247 Double_t *jj = A;
248 for( Int_t j=1; j<i; jj+=++j, ij++ ){
249 Double_t *ik = ij, *kj = jj, sum = 0.;
250 for( Int_t k=j; ik!=ii; kj+=k++, ik++ ){
251 sum += *ik * *kj;
252 }
253 *kj = sum * x;
254 }
255 }
256 }else{
257 for( Double_t *ik = ij; ik!=ii+1; ik++ ){
258 *ik = 0.;
259 }
260 ret = 1;
261 }
262 }
263 }
264
265 // A -> A^T x A
266 // method:
267 // Aij = sum_{k=i}^N Aki * Akj
268
269 {
270 Double_t *ii=A, *ij=A;
271 for( Int_t i=1; i<=N; ii+=++i ){
272 do{
273 Double_t *ki = ii, *kj = ij, sum = 0.;
274 for( Int_t k=i; k<=N; ki+=k, kj+=k++ ) sum += (*ki) * (*kj);
275 *ij = sum;
276 }while( (ij++)!=ii );
277 }
278 }
279 return ret;
280 }
281
282Double_t CbmKFMath::getDeviation( Double_t x, Double_t y, Double_t C[], Double_t vx, Double_t vy, Double_t Cv[] )
283 {
284 Double_t dx = x - vx;
285 Double_t dy = y - vy;
286 Double_t c[3] = { 0, 0, 0 };
287 if( C ) { c[0] = C [0]; c[1] = C [1]; c[2] = C [2]; }
288 if( Cv ) { c[0]+= Cv[0]; c[1]+= Cv[1]; c[2]+= Cv[2]; }
289 Double_t d = c[0]*c[2] - c[1]*c[1] ;
290 if( fabs(d)<1.e-20 ) return 0;
291 return sqrt( fabs( 0.5*(dx*dx*c[0]-2*dx*dy*c[1]+dy*dy*c[2])/d ));
292 }
293
294Double_t CbmKFMath::AnalyticQP(
295 const Double_t T[], // track parameters (x,y,tx,ty,Q/p,z)
296 const Double_t V[], // vertex parameters (x,y,z)
297 FairField *MagneticField // magnetic field
298 )
299{
300
301 const Double_t c_light = 0.000299792458;
302
303 Double_t
304 x = T[0],
305 y = T[1],
306 tx = T[2],
307 ty = T[3],
308 qp0 = T[4],
309 z = T[5],
310
311 vx = V[0],
312 //vy = V[1],
313 vz = V[2],
314
315 txtx = tx*tx,
316 tyty = ty*ty,
317 txty = tx*ty;
318
319 Double_t
320 Ax = txty,
321 Ay = -txtx-1,
322 Az = ty,
323 Ayy = tx*(txtx*3+3),
324 Ayz = -2*txty,
325 Azy = Ayz,
326 Ayyy = -(15*txtx*txtx+18*txtx+3),
327
328 Bx = tyty+1,
329 By = -txty,
330 Bz = -tx,
331 Byy = ty*(txtx*3+1),
332 Byz = 2*txtx+1,
333 Bzy = txtx-tyty,
334 Byyy = -txty*(txtx*15+9);
335
336 Double_t
337 t = c_light*sqrt( 1. + txtx + tyty ),
338 h = t*qp0;
339
340 // integrals
341
342 Double_t Sx=0, Sy=0, Sz=0, Syy=0, Syz=0, Szy=0, Syyy=0;
343
344 {
345 Double_t sx=0,sy=0,sz=0,syy=0,syz=0, szy=0, syyy=0;
346
347 Double_t r[3] = { x, y, z };
348 Double_t B[3][3];
349
350 Int_t n = int( fabs(vz-z)/5.);
351 if( n<2 ) n = 2;
352 Double_t dz = (vz-z)/n;
353
354 {
355 MagneticField->GetFieldValue( r, B[0] );
356 r[0] += tx*dz;
357 r[1] += ty*dz;
358 r[2] += dz;
359 MagneticField->GetFieldValue( r, B[1] );
360 r[0] += tx*dz;
361 r[1] += ty*dz;
362 r[2] += dz;
363 MagneticField->GetFieldValue( r, B[2] );
364
365 sx = ( B[0][0] + 4*B[1][0] + B[2][0] )*dz*2/6.;
366 sy = ( B[0][1] + 4*B[1][1] + B[2][1] )*dz*2/6.;
367 sz = ( B[0][2] + 4*B[1][2] + B[2][2] )*dz*2/6.;
368
369 Sx = ( B[0][0] + 2*B[1][0])*dz*dz*4/6.;
370 Sy = ( B[0][1] + 2*B[1][1])*dz*dz*4/6.;
371 Sz = ( B[0][2] + 2*B[1][2])*dz*dz*4/6.;
372
373 Double_t c2[3][3] = { { 5, -4, -1},{ 44, 80, -4},{ 11, 44, 5} }; // /=360.
374 Double_t C2[3][3] = { { 38, 8, -4},{ 148, 208, -20},{ 3, 36, 3} }; // /=2520.
375 for(Int_t k=0; k<3; k++)
376 for(Int_t m=0; m<3; m++)
377 {
378 syz += c2[k][m]*B[k][1]*B[m][2];
379 Syz += C2[k][m]*B[k][1]*B[m][2];
380 szy += c2[k][m]*B[k][2]*B[m][1];
381 Szy += C2[k][m]*B[k][2]*B[m][1];
382 }
383
384 syz *= dz*dz*4/360.;
385 Syz *= dz*dz*dz*8/2520.;
386
387 szy *= dz*dz*4/360.;
388 Szy *= dz*dz*dz*8/2520.;
389
390 syy = ( B[0][1] + 4*B[1][1] + B[2][1] )*dz*2;
391 syyy = syy*syy*syy / 1296;
392 syy = syy*syy/72;
393
394 Syy = ( B[0][1]*( 38*B[0][1] + 156*B[1][1] - B[2][1] )+
395 B[1][1]*( 208*B[1][1] +16*B[2][1] )+
396 B[2][1]*( 3*B[2][1] )
397 )*dz*dz*dz*8/2520.;
398 Syyy =
399 (
400 B[0][1]*( B[0][1]*( 85*B[0][1] + 526*B[1][1] - 7*B[2][1] )+
401 B[1][1]*( 1376*B[1][1] +84*B[2][1] )+
402 B[2][1]*( 19*B[2][1] ) )+
403 B[1][1]*( B[1][1]*( 1376*B[1][1] +256*B[2][1] )+
404 B[2][1]*( 62*B[2][1] ) )+
405 B[2][1]*B[2][1] *( 3*B[2][1] )
406 )*dz*dz*dz*dz*16/90720.;
407
408 x += tx*2*dz
409 + h*( Sx*Ax + Sy*Ay + Sz*Az )
410 + h*h*( Syy*Ayy + Syz*Ayz + Szy*Azy )
411 + h*h*h*Syyy*Ayyy;
412 y += ty*2*dz
413 + h*( Sx*Bx + Sy*By + Sz*Bz )
414 + h*h*( Syy*Byy + Syz*Byz + Szy*Bzy )
415 + h*h*h*Syyy*Byyy;
416 tx += h*( sx*Ax + sy*Ay + sz*Az )
417 + h*h*( syy*Ayy + syz*Ayz + szy*Azy )
418 + h*h*h*syyy*Ayyy;
419 ty += h*( sx*Bx + sy*By + sz*Bz )
420 + h*h*( syy*Byy + syz*Byz + szy*Bzy )
421 + h*h*h*syyy*Byyy;
422 z += 2*dz;
423 txtx = tx*tx;
424 tyty = ty*ty;
425 txty = tx*ty;
426
427 Ax = txty;
428 Ay = -txtx-1;
429 Az = ty;
430 Ayy = tx*(txtx*3+3);
431 Ayz = -2*txty;
432 Azy = Ayz;
433 Ayyy = -(15*txtx*txtx+18*txtx+3);
434
435 Bx = tyty+1;
436 By = -txty;
437 Bz = -tx;
438 Byy = ty*(txtx*3+1);
439 Byz = 2*txtx+1;
440 Bzy = txtx-tyty;
441 Byyy = -txty*(txtx*15+9);
442
443 t = c_light*sqrt( 1. + txtx + tyty );
444 h = t*qp0;
445 }
446
447 for( Int_t i=2; i<n; i++ )
448 {
449
450 Double_t B0[3] = { B[0][0], B[0][1], B[0][2] };
451
452 B[0][0] = B[1][0];
453 B[0][1] = B[1][1];
454 B[0][2] = B[1][2];
455
456 B[1][0] = B[2][0];
457 B[1][1] = B[2][1];
458 B[1][2] = B[2][2];
459
460 // estimate B[2] at +dz
461
462 B[2][0] = B0[0] -3*B[0][0] + 3*B[1][0];
463 B[2][1] = B0[1] -3*B[0][1] + 3*B[1][1];
464 B[2][2] = B0[2] -3*B[0][2] + 3*B[1][2];
465
466 Double_t Sx_, Sy_, Sz_, Syy_, Syz_, Szy_, Syyy_;
467
468 Sx_ = ( -B[0][0] + 10*B[1][0] + 3*B[2][0] )*dz*dz*4/96.;
469 Sy_ = ( -B[0][1] + 10*B[1][1] + 3*B[2][1] )*dz*dz*4/96.;
470 Sz_ = ( -B[0][2] + 10*B[1][2] + 3*B[2][2] )*dz*dz*4/96.;
471
472 Syz_ = Szy_ = 0;
473
474 Double_t c2[3][3] = { { 5, -52, -13},{ -28, 320, 68},{ -37, 332, 125} }; // /=5760
475 Double_t C2[3][3] = { { 13, -152, -29},{ -82, 1088, 170},{ -57, 576, 153} }; // /=80640
476 { for(Int_t k=0; k<3; k++)
477 for(Int_t m=0; m<3; m++)
478 {
479 Syz_ += C2[k][m]*B[k][1]*B[m][2];
480 Szy_ += C2[k][m]*B[k][2]*B[m][1];
481 }
482 }
483
484 Syz_ *= dz*dz*dz*8/80640.;
485 Szy_ *= dz*dz*dz*8/80640.;
486
487 Syy_ = ( B[0][1]*( 13*B[0][1] - 234*B[1][1] - 86*B[2][1] )+
488 B[1][1]*( 1088*B[1][1] +746*B[2][1] )+
489 B[2][1]*( 153*B[2][1] )
490 )*dz*dz*dz*8/80640.;
491
492 Syyy_ =
493 (
494 B[0][1]*( B[0][1]*( -43*B[0][1] + 1118*B[1][1]+451*B[2][1] )+
495 B[1][1]*( -9824*B[1][1] -7644*B[2][1] )+
496 B[2][1]*( -1669*B[2][1] ) )+
497 B[1][1]*( B[1][1]*( 29344*B[1][1] +32672*B[2][1] )+
498 B[2][1]*( 13918*B[2][1] ) )+
499 B[2][1]*B[2][1] *( 2157*B[2][1] )
500 )*dz*dz*dz*dz*16/23224320.;
501
502 r[2] += dz;
503 r[0] =
504 x + tx*dz
505 + h*( Sx_*Ax + Sy_*Ay + Sz_*Az )
506 + h*h*( Syy_*Ayy + Syz_*Ayz + Szy_*Azy )
507 + h*h*h*Syyy_*Ayyy;
508 r[1] = y + ty*dz
509 + h*( Sx_*Bx + Sy_*By + Sz_*Bz )
510 + h*h*( Syy_*Byy + Syz_*Byz + Szy_*Bzy )
511 + h*h*h*Syyy_*Byyy;
512 /*
513 Syyy_+= dz*syyy+Sy_*syy+Syy_*sy+Syyy;
514 Syy_ += dz*syy + sy*Sy_ + Syy;
515 Syz_ += dz*syz + sz*Sy_ + Syz;
516 Szy_ += dz*szy + sy*Sz_ + Szy;
517 Sx_ += dz*sx + Sx;
518 Sy_ += dz*sy + Sy;
519 Sz_ += dz*sz + Sz;
520
521 r[2] += dz;
522 r[0] =
523 x + tx*(r[2]-z)
524 + h*( Sx_*Ax + Sy_*Ay + Sz_*Az )
525 + h*h*( Syy_*Ayy + Syz_*Ayz + Szy_*Azy )
526 + h*h*h*Syyy_*Ayyy;
527 r[1] = y + ty*(r[2]-z)
528 + h*( Sx_*Bx + Sy_*By + Sz_*Bz )
529 + h*h*( Syy_*Byy + Syz_*Byz + Szy_*Bzy )
530 + h*h*h*Syyy_*Byyy;
531 */
532 MagneticField->GetFieldValue( r, B[2] );
533
534 Double_t sx_, sy_, sz_, syy_, syz_, szy_, syyy_;
535
536 sx_ = ( -B[0][0] + 8*B[1][0] + 5*B[2][0] )*dz*2/24.;
537 sy_ = ( -B[0][1] + 8*B[1][1] + 5*B[2][1] )*dz*2/24.;
538 sz_ = ( -B[0][2] + 8*B[1][2] + 5*B[2][2] )*dz*2/24.;
539
540 Sx_ = ( -B[0][0] + 10*B[1][0] + 3*B[2][0] )*dz*dz*4/96.;
541 Sy_ = ( -B[0][1] + 10*B[1][1] + 3*B[2][1] )*dz*dz*4/96.;
542 Sz_ = ( -B[0][2] + 10*B[1][2] + 3*B[2][2] )*dz*dz*4/96.;
543
544 syz_ = Syz_ = szy_ = Szy_ = 0;
545
546 { for(Int_t k=0; k<3; k++)
547 for(Int_t m=0; m<3; m++)
548 {
549 syz_ += c2[k][m]*B[k][1]*B[m][2];
550 Syz_ += C2[k][m]*B[k][1]*B[m][2];
551 szy_ += c2[k][m]*B[k][2]*B[m][1];
552 Szy_ += C2[k][m]*B[k][2]*B[m][1];
553 }
554 }
555 syz_ *= dz*dz*4/5760.;
556 Syz_ *= dz*dz*dz*8/80640.;
557
558 szy_ *= dz*dz*4/5760.;
559 Szy_ *= dz*dz*dz*8/80640.;
560
561 syy_ = ( B[0][1] -8*B[1][1] - 5*B[2][1] )*dz*2;
562 syyy_ = -syy_*syy_*syy_/82944;
563 syy_ = syy_*syy_/1152;
564
565 Syy_ = ( B[0][1]*( 13*B[0][1] - 234*B[1][1] - 86*B[2][1] )+
566 B[1][1]*( 1088*B[1][1] +746*B[2][1] )+
567 B[2][1]*( 153*B[2][1] )
568 )*dz*dz*dz*8/80640.;
569
570 Syyy_ =
571 (
572 B[0][1]*( B[0][1]*( -43*B[0][1] + 1118*B[1][1]+451*B[2][1] )+
573 B[1][1]*( -9824*B[1][1] -7644*B[2][1] )+
574 B[2][1]*( -1669*B[2][1] ) )+
575 B[1][1]*( B[1][1]*( 29344*B[1][1] +32672*B[2][1] )+
576 B[2][1]*( 13918*B[2][1] ) )+
577 B[2][1]*B[2][1] *( 2157*B[2][1] )
578 )*dz*dz*dz*dz*16/23224320.;
579
580 x += tx*dz
581 + h*( Sx_*Ax + Sy_*Ay + Sz_*Az )
582 + h*h*( Syy_*Ayy + Syz_*Ayz + Szy_*Azy )
583 + h*h*h*Syyy_*Ayyy;
584 y += ty*dz
585 + h*( Sx_*Bx + Sy_*By + Sz_*Bz )
586 + h*h*( Syy_*Byy + Syz_*Byz + Szy_*Bzy )
587 + h*h*h*Syyy_*Byyy;
588 tx += h*( sx_*Ax + sy_*Ay + sz_*Az )
589 + h*h*( syy_*Ayy + syz_*Ayz + szy_*Azy )
590 + h*h*h*syyy_*Ayyy;
591 ty += h*( sx_*Bx + sy_*By + sz_*Bz )
592 + h*h*( syy_*Byy + syz_*Byz + szy_*Bzy )
593 + h*h*h*syyy_*Byyy;
594 z += dz;
595 txtx = tx*tx;
596 tyty = ty*ty;
597 txty = tx*ty;
598
599 Ax = txty;
600 Ay = -txtx-1;
601 Az = ty;
602 Ayy = tx*(txtx*3+3);
603 Ayz = -2*txty;
604 Azy = Ayz;
605 Ayyy = -(15*txtx*txtx+18*txtx+3);
606
607 Bx = tyty+1;
608 By = -txty;
609 Bz = -tx;
610 Byy = ty*(txtx*3+1);
611 Byz = 2*txtx+1;
612 Bzy = txtx-tyty;
613 Byyy = -txty*(txtx*15+9);
614
615 t = c_light*sqrt( 1. + txtx + tyty );
616 h = t*qp0;
617
618
619 Syyy+= dz*syyy+Sy_*syy+Syy_*sy+Syyy_;
620 Syy += dz*syy + sy*Sy_ + Syy_;
621 Syz += dz*syz + sz*Sy_ + Syz_;
622 Szy += dz*szy + sy*Sz_ + Szy_;
623 Sx += dz*sx + Sx_;
624 Sy += dz*sy + Sy_;
625 Sz += dz*sz + Sz_;
626
627 syyy+= sy_*syy+syy_*sy+syyy_;
628 syy += sy*sy_ + syy_;
629 syz += sz*sy_ + syz_;
630 szy += sz*sz_ + szy_;
631 sx += sx_;
632 sy += sy_;
633 sz += sz_;
634 }
635 //cout<<x-vx<<" "<<y-vy<<" "<<z-vz<<" "<<1./qp0<<endl;
636 }
637
638 x = T[0];
639 y = T[1];
640 tx = T[2];
641 ty = T[3];
642 z = T[5];
643 txtx = tx*tx;
644 tyty = ty*ty;
645 txty = tx*ty;
646
647 Ax = txty;
648 Ay = -txtx-1;
649 Az = ty;
650 Ayy = tx*(txtx*3+3);
651 Ayz = -2*txty;
652 Azy = Ayz;
653 Ayyy = -(15*txtx*txtx+18*txtx+3);
654
655 Bx = tyty+1;
656 By = -txty;
657 Bz = -tx;
658 Byy = ty*(txtx*3+1);
659 Byz = 2*txtx+1;
660 Bzy = txtx-tyty;
661 Byyy = -txty*(txtx*15+9);
662
663 t = c_light*sqrt( 1. + txtx + tyty );
664 h = t*qp0;
665
666 Double_t
667
668 c = (x-vx) + tx*(vz-z),
669 b = t*( Sx*Ax + Sy*Ay + Sz*Az ) ,
670 a = t*t*( Syy*Ayy + Syz*Ayz + Szy*Azy ),
671 d = t*t*t*( Syyy*Ayyy);
672
673 Double_t
674 C = d*qp0*qp0*qp0 + a*qp0*qp0 + b*qp0 + c,
675 B = 3*d*qp0*qp0+2*a*qp0+b,
676 A = 3*d + a;
677
678 Double_t D = B*B-4*A*C;
679 if( D<0 ) D = 0;
680 D = sqrt( D );
681
682 Double_t dqp1 = (-B+D)/2/A;
683 Double_t dqp2 = (-B-D)/2/A;
684 Double_t dqp = (fabs(dqp1)<fabs(dqp2)) ?dqp1 : dqp2;
685
686 return qp0+dqp;
687}
688
689
690Bool_t CbmKFMath::GetThickness( Double_t z1, Double_t z2, Double_t mz, Double_t mthick,
691 Double_t *mz_out, Double_t *mthick_out )
692{
693 Double_t z, Z;
694
695 if ( z1<z2 )
696 {
697 z = z1;
698 Z = z2;
699 }
700 else
701 {
702 Z = z1;
703 z = z2;
704 }
705
706 Double_t
707 tmin = mz - mthick/2,
708 tmax = mz + mthick/2;
709 if ( tmin >= Z || z >= tmax ) return 1;
710 if ( z<=tmin && tmax<=Z ) // case z |**| Z
711 {
712 *mz_out = mz;
713 *mthick_out = mthick;
714 }
715 else if ( z<=tmin && tmin<Z && Z<=tmax ) // case z |*Z*|
716 {
717 *mz_out = (tmin + Z)/2;
718 *mthick_out = Z - tmin;
719 }
720 else if ( tmax<=Z && tmin<=z && z<tmax ) // case |*z*| Z
721 {
722 *mz_out = (tmax + z)/2;
723 *mthick_out = tmax - z;
724 }
725 else if ( tmin<=z && Z<=tmax ) // case |*zZ*|
726 {
727 *mz_out = (Z + z)/2;
728 *mthick_out = Z - z;
729 }
730 return 0;
731}
732
733
734Int_t CbmKFMath::GetNoise( Double_t Lrl, Double_t F,Double_t Fe,
735 Double_t tx, Double_t ty, Double_t qp,
736 Double_t mass, Bool_t is_electron,
737 Bool_t downstream_direction,
738 Double_t *Q5, Double_t *Q8, Double_t *Q9, Double_t *Ecor )
739{
740 *Q5 = 0;
741 *Q8 = 0;
742 *Q9 = 0;
743 *Ecor = 1.;
744 Double_t t = sqrt(1.+tx*tx+ty*ty);
745 if( !finite(t) || t>1.e4 ) return 1;
746 t = sqrt(t);
747 Double_t l = t*Lrl;
748 if ( l > 1.) l = 1.; // protect against l too big
749 Double_t s0 = (l>exp(-1./0.038) ) ? F*.0136*(1+0.038*log(l))*qp :0.;
750 Double_t a = (1.+mass*mass*qp*qp)*s0*s0*t*t*l;
751
752 *Q5 = a*(1.+tx*tx);
753 *Q8 = a*tx*ty;
754 *Q9 = a*(1.+ty*ty);
755
756 // Apply correction for dE/dx energy loss
757
758 Double_t L = downstream_direction ?l :-l;
759
760 if( 0&&is_electron ) // Bremsstrahlung effect for e+,e-
761 {
762 *Ecor = exp( L );
763 }
764
765 if(1) // ionisation energy loss
766 {
767 Double_t m_energyLoss = Fe;
768// Double_t m_energyLoss = 0.02145;//0.060;
769 Double_t corr = (1.- fabs(qp) * L * m_energyLoss );
770 if( corr>0.001*fabs(qp)) *Ecor *= 1./corr;
771 else *Ecor = 1000./fabs(qp);
772 }
773 return 0;
774}
775
776
777void CbmKFMath::CopyTC2TrackParam(FairTrackParam* par, Double_t T[], Double_t C[] )
778{
779 if( T )
780 {
781 par->SetX (T[0]);
782 par->SetY (T[1]);
783 par->SetZ (T[5]);
784 par->SetTx(T[2]);
785 par->SetTy(T[3]);
786 par->SetQp(T[4]);
787 }
788 if( C )
789 {
790 for (Int_t i=0,iCov=0; i<5; i++)
791 for (Int_t j=0; j<=i; j++,iCov++)
792 par->SetCovariance(i,j,C[iCov]);
793 }
794}
795
796void CbmKFMath::CopyTrackParam2TC(FairTrackParam* par, Double_t T[], Double_t C[] )
797{
798 if( T )
799 {
800 T[0]=par->GetX();
801 T[1]=par->GetY();
802 T[2]=par->GetTx();
803 T[3]=par->GetTy();
804 T[4]=par->GetQp();
805 T[5]=par->GetZ();
806 }
807 if( C )
808 {
809 for (Int_t i=0,iCov=0; i<5; i++)
810 for (Int_t j=0; j<=i; j++,iCov++)
811 C[iCov] = par->GetCovariance(i,j);
812 }
813}
d
const Float_t d
Z-ccordinate of the first GEM-station.
Definition BmnMwpcHit.cxx:7
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition P4_F32vec4.h:34
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition P4_F32vec4.h:52
i
int i
Definition P4_F32vec4.h:22
m
__m128 m
Definition P4_F32vec4.h:27
log
friend F32vec4 log(const F32vec4 &a)
Definition P4_F32vec4.h:123
exp
friend F32vec4 exp(const F32vec4 &a)
Definition P4_F32vec4.h:122
CbmKFMath::multQtSQ
static void multQtSQ(Int_t N, const Double_t Q[], const Double_t S[], Double_t S_out[])
Definition CbmKFMath.cxx:62
CbmKFMath::CopyTC2TrackParam
static void CopyTC2TrackParam(FairTrackParam *par, Double_t T[], Double_t C[])
Definition CbmKFMath.cxx:777
CbmKFMath::getDeviation
static Double_t getDeviation(Double_t x, Double_t y, Double_t C[], Double_t vx, Double_t vy, Double_t Cv[]=0)
Definition CbmKFMath.cxx:282
CbmKFMath::CopyTrackParam2TC
static void CopyTrackParam2TC(FairTrackParam *par, Double_t T[], Double_t C[])
Definition CbmKFMath.cxx:796
CbmKFMath::GetThickness
static Bool_t GetThickness(Double_t z1, Double_t z2, Double_t mz, Double_t mthick, Double_t *mz_out, Double_t *mthick_out)
Definition CbmKFMath.cxx:690
CbmKFMath::intersectCone
static Bool_t intersectCone(Double_t zCone, Double_t ZCone, Double_t rCone, Double_t RCone, const Double_t x[], Double_t *z1, Double_t *z2)
Definition CbmKFMath.cxx:13
CbmKFMath::AnalyticQP
static Double_t AnalyticQP(const Double_t T[], const Double_t V[], FairField *MagneticField)
Definition CbmKFMath.cxx:294
CbmKFMath::multSSQ
static void multSSQ(const Double_t *A, const Double_t *B, Double_t *C, Int_t n)
Definition CbmKFMath.cxx:82
CbmKFMath::five_dim_inv
static void five_dim_inv(Double_t a[5][5])
Definition CbmKFMath.cxx:141
CbmKFMath::invS
static Bool_t invS(Double_t A[], Int_t N)
Definition CbmKFMath.cxx:192
CbmKFMath::GetNoise
static Int_t GetNoise(Double_t Lrl, Double_t F, Double_t Fe, Double_t tx, Double_t ty, Double_t qp, Double_t mass, Bool_t is_electron, Bool_t downstream_direction, Double_t *Q5, Double_t *Q8, Double_t *Q9, Double_t *Ecor)
Definition CbmKFMath.cxx:734
CbmKFMath::four_dim_inv
static void four_dim_inv(Double_t a[4][4])
Definition CbmKFMath.cxx:95
CbmKFMath::multQSQt
static void multQSQt(Int_t N, const Double_t Q[], const Double_t S[], Double_t S_out[])
Definition CbmKFMath.cxx:38
CbmKFMath::indexS
static Int_t indexS(Int_t i, Int_t j)
Definition CbmKFMath.h:33