BmnRoot
Loading...
Searching...
No Matches
CbmKFParticleInterface.cxx
Go to the documentation of this file.
1#include "CbmKFParticleInterface.h"
2
3#define cnst static const fvec
4
5//for particle finding
6#include "CbmL1Track.h"
7#include "CbmStsTrack.h"
8#include "CbmKFTrack.h"
9
10#include <map>
11using std::map;
12
13const float CbmKFParticleInterface::DefaultCuts[2][3] = {{3.,3.,-100.},{3.,3.,-100.}};
14
18
23
24CbmKFParticleInterface::CbmKFParticleInterface(CbmKFParticle* part[], const int nPart):KFPart()
25{
26 KFPart.Create(part,nPart);
27}
28
32
37
43
44void CbmKFParticleInterface::Construct(CbmKFTrackInterface* vDaughters[][fvecLen],int NDaughters,
45 CbmKFVertexInterface *Parent,float Mass,float CutChi2)
46{
47 const int ND = NDaughters;
48
49 CbmKFParticle_simd* vDaughters_part = new CbmKFParticle_simd[ND];
50
51 for( int tr=0; tr<NDaughters; ++tr )
52 vDaughters_part[tr].Create(vDaughters[tr]);
53
54 Construct( vDaughters_part, NDaughters, Parent, Mass, CutChi2 );
55 delete[] vDaughters_part;
56}
57
58void CbmKFParticleInterface::Construct( CbmKFParticle_simd* vDaughters, int NDaughters,
59 CbmKFVertexInterface *Parent, float Mass, float CutChi2, bool isAtVtxGuess )
60{
61 KFPart.CleanDaughtersId();
62 KFPart.SetNDaughters(NDaughters);
63 for( int tr=0; tr<NDaughters; tr++ )
64 {
65 KFPart.AddDaughterId( vDaughters[tr].Id() );
66 }
67
68 int MaxIter=3;
69 if(!(KFPart.fIsVtxGuess))
70 {
71 fvec Ty[2];
72 fvec Y[2];
73 fvec Z[2];
74
75 for( int tr=0; tr<2; tr++ )
76 {
77 Y[tr] = vDaughters[tr].GetY();
78 Z[tr] = vDaughters[tr].GetZ();
79 Ty[tr] = vDaughters[tr].GetPy()/vDaughters[tr].GetPz();
80 }
81
82 fvec r0[3];
83 r0[0] = 0.;
84 r0[1] = 0.;
85 r0[2] = (Ty[0]*Z[0]-Ty[1]*Z[1] + Y[1] -Y[0])/(Ty[0]-Ty[1]);
86
87 CbmKFParticle_simd Daughter0 = vDaughters[0];
88 Extrapolate( &Daughter0, Daughter0.r , GetDStoPoint(Daughter0, r0) );
89
90 KFPart.fVtxGuess[0] = Daughter0.GetX();
91 KFPart.fVtxGuess[1] = Daughter0.GetY();
92 KFPart.fVtxGuess[2] = Daughter0.GetZ();
93
94 if(!(KFPart.fIsVtxErrGuess))
95 {
96 KFPart.fVtxErrGuess[0] = 10.*sqrt(Daughter0.C[0] );
97 KFPart.fVtxErrGuess[1] = 10.*sqrt(Daughter0.C[2] );
98 KFPart.fVtxErrGuess[2] = 10.*sqrt(Daughter0.C[5] );
99 }
100 // KFPart.fVtxErrGuess[0] = 0.01*KFPart.fVtxGuess[0];
101 // KFPart.fVtxErrGuess[1] = 0.01*KFPart.fVtxGuess[1];
102 // KFPart.fVtxErrGuess[2] = 0.01*KFPart.fVtxGuess[2];
103 }
104 else
105 {
106 MaxIter = 1;
107 }
108
109 cnst ZERO = 0.0, ONE = 1.;
110 cnst TWO = 2.;
124 KFPart.r[3] = ZERO;
125 KFPart.r[4] = ZERO;
126 KFPart.r[5] = ZERO;
127 KFPart.r[6] = ZERO;
128 KFPart.r[7] = ZERO;
129
130 KFPart.r[0] = KFPart.fVtxGuess[0];
131 KFPart.r[1] = KFPart.fVtxGuess[1];
132 KFPart.r[2] = KFPart.fVtxGuess[2];
133
134 KFPart.C[0] = KFPart.fVtxErrGuess[0] * KFPart.fVtxErrGuess[0];
135 KFPart.C[2] = KFPart.fVtxErrGuess[1] * KFPart.fVtxErrGuess[1];
136 KFPart.C[5] = KFPart.fVtxErrGuess[2] * KFPart.fVtxErrGuess[2];
137 KFPart.C[1] = KFPart.C[3] = KFPart.C[4] = ZERO;
138
139
140 for ( Int_t iteration =0; iteration<MaxIter;++iteration )
141 {
142 fvec r0[8], C0[6];
143
144 for( int i=0; i<8; i++) r0[i] = KFPart.r[i];
145 for( int i=0; i<6; i++) C0[i] = KFPart.C[i];
146
147 KFPart.r[3] = ZERO;
148 KFPart.r[4] = ZERO;
149 KFPart.r[5] = ZERO;
150 KFPart.r[6] = ZERO;
151
152 for(Int_t i=0;i<36;++i) KFPart.C[i]=ZERO;
153
154 KFPart.C[0] = KFPart.fVtxErrGuess[0]*KFPart.fVtxErrGuess[0];
155 KFPart.C[2] = KFPart.fVtxErrGuess[1]*KFPart.fVtxErrGuess[1];
156 KFPart.C[5] = KFPart.fVtxErrGuess[2]*KFPart.fVtxErrGuess[2];
157
158 KFPart.C[35] = ONE;
159
160 L1FieldValue B = vDaughters[0].fField.Get(r0[2]);
161 B.x *= 0.000299792458; //multiply by c
162 B.y *= 0.000299792458;
163 B.z *= 0.000299792458;
164
165 KFPart.NDF = -3.;
166 KFPart.Chi2 = ZERO;
167 KFPart.Q = ZERO;
168
169 for( int tr=0; tr<NDaughters; ++tr )
170 {
171 CbmKFParticle_simd Daughter = vDaughters[tr] ;
172 if(!isAtVtxGuess)
173 Extrapolate(&Daughter, Daughter.r , GetDStoPoint(Daughter, r0));
174
175 fvec *m = Daughter.r;
176 fvec *Cd = Daughter.C;
177
178 fvec d[3] = { r0[0]-m[0], r0[1]-m[1], r0[2]-m[2] };
179 fvec SigmaS = .1+10.*sqrt( (d[0]*d[0]+d[1]*d[1]+d[2]*d[2])/
180 (m[3]*m[3]+m[4]*m[4]+m[5]*m[5]) );
181
182 fvec h[6];
183
184 h[0] = m[3]*SigmaS;
185 h[1] = m[4]*SigmaS;
186 h[2] = m[5]*SigmaS;
187 h[3] = ( h[1]*B.z-h[2]*B.y )*Daughter.Q;
188 h[4] = ( h[2]*B.x-h[0]*B.z )*Daughter.Q;
189 h[5] = ( h[0]*B.y-h[1]*B.x )*Daughter.Q;
190
191 //* Fit of daughter momentum (x,y,z) to r0[0,1,2] vertex
192 {
193 fvec zeta[3] = { r0[0]-m[0], r0[1]-m[1], r0[2]-m[2] };
194
195 fvec Vv[6] =
196 { Cd[ 0] + h[0]*h[0],
197 Cd[ 1] + h[1]*h[0], Cd[ 2] + h[1]*h[1],
198 Cd[ 3] + h[2]*h[0], Cd[ 4] + h[2]*h[1], Cd[ 5] + h[2]*h[2] };
199
200 fvec Vvp[9]=
201 { Cd[ 6] + h[0]*h[3], Cd[ 7] + h[1]*h[3], Cd[ 8] + h[2]*h[3],
202 Cd[10] + h[0]*h[4], Cd[11] + h[1]*h[4], Cd[12] + h[2]*h[4],
203 Cd[15] + h[0]*h[5], Cd[16] + h[1]*h[5], Cd[17] + h[2]*h[5] };
204
224 fvec S[6] =
225 { Vv[2]*Vv[5] - Vv[4]*Vv[4],
226 Vv[3]*Vv[4] - Vv[1]*Vv[5], Vv[0]*Vv[5] - Vv[3]*Vv[3],
227 Vv[1]*Vv[4] - Vv[2]*Vv[3], Vv[1]*Vv[3] - Vv[0]*Vv[4],
228 Vv[0]*Vv[2] - Vv[1]*Vv[1] };
229
230 fvec s = ( Vv[0]*S[0] + Vv[1]*S[1] + Vv[3]*S[3] );
231
232 fvec init = fvec(s > 1.E-20);
233 s = rcp(s) & init;
234// for(int j = 0; j<fvecLen; j++)
235// s[j] = ( s[j] > 1.E-20 ) ?1./s[j] :0;
236
237 S[0]*=s; S[1]*=s; S[2]*=s; S[3]*=s; S[4]*=s; S[5]*=s;
238
239 fvec Sz[3] = { (S[0]*zeta[0]+S[1]*zeta[1]+S[3]*zeta[2]),
240 (S[1]*zeta[0]+S[2]*zeta[1]+S[4]*zeta[2]),
241 (S[3]*zeta[0]+S[4]*zeta[1]+S[5]*zeta[2]) };
242
243 fvec x = m[3] + Vvp[0]*Sz[0] + Vvp[1]*Sz[1] + Vvp[2]*Sz[2];
244 fvec y = m[4] + Vvp[3]*Sz[0] + Vvp[4]*Sz[1] + Vvp[5]*Sz[2];
245 fvec z = m[5] + Vvp[6]*Sz[0] + Vvp[7]*Sz[1] + Vvp[8]*Sz[2];
246
247 h[0] = x*SigmaS;
248 h[1] = y*SigmaS;
249 h[2] = z*SigmaS;
250 h[3] = ( h[1]*B.z-h[2]*B.y )*Daughter.Q;
251 h[4] = ( h[2]*B.x-h[0]*B.z )*Daughter.Q;
252 h[5] = ( h[0]*B.y-h[1]*B.x )*Daughter.Q;
253 }
254
255 fvec* V = Cd;
256
257 V[ 0] += h[0]*h[0];
258 V[ 1] += h[1]*h[0];
259 V[ 2] += h[1]*h[1];
260 V[ 3] += h[2]*h[0];
261 V[ 4] += h[2]*h[1];
262 V[ 5] += h[2]*h[2];
263
264 V[ 6] += h[3]*h[0];
265 V[ 7] += h[3]*h[1];
266 V[ 8] += h[3]*h[2];
267 V[ 9] += h[3]*h[3];
268
269 V[10] += h[4]*h[0];
270 V[11] += h[4]*h[1];
271 V[12] += h[4]*h[2];
272 V[13] += h[4]*h[3];
273 V[14] += h[4]*h[4];
274
275 V[15] += h[5]*h[0];
276 V[16] += h[5]*h[1];
277 V[17] += h[5]*h[2];
278 V[18] += h[5]*h[3];
279 V[19] += h[5]*h[4];
280 V[20] += h[5]*h[5];
281
282// V[21] = Cd[21];
283// V[22] = Cd[22];
284// V[23] = Cd[23];
285// V[24] = Cd[24];
286// V[25] = Cd[25];
287// V[26] = Cd[26];
288// V[27] = Cd[27];
289
290 //*
291 KFPart.NDF += TWO;
292 KFPart.Q += Daughter.Q;
293
294 //*
295
296 fvec S[6];
297 {
298 fvec Si[6] = { KFPart.C[0]+V[0],
299 KFPart.C[1]+V[1], KFPart.C[2]+V[2],
300 KFPart.C[3]+V[3], KFPart.C[4]+V[4], KFPart.C[5]+V[5] };
301
302 S[0] = Si[2]*Si[5] - Si[4]*Si[4];
303 S[1] = Si[3]*Si[4] - Si[1]*Si[5];
304 S[2] = Si[0]*Si[5] - Si[3]*Si[3];
305 S[3] = Si[1]*Si[4] - Si[2]*Si[3];
306 S[4] = Si[1]*Si[3] - Si[0]*Si[4];
307 S[5] = Si[0]*Si[2] - Si[1]*Si[1];
308
309 fvec s = ( Si[0]*S[0] + Si[1]*S[1] + Si[3]*S[3] );
310
311 fvec init = fvec(s > 1.E-20);
312 s = rcp(s) & init;
313// for(int j=0; j<fvecLen; j++)
314// s[j] = ( s[j] > 1.E-20 ) ?1./s[j] :0;
315
316 S[0]*=s;
317 S[1]*=s;
318 S[2]*=s;
319 S[3]*=s;
320 S[4]*=s;
321 S[5]*=s;
322 }
323
324 //* Residual (measured - estimated)
325
326 fvec zeta[3] = { m[0]-KFPart.r[0], m[1]-KFPart.r[1], m[2]-KFPart.r[2] };
327
328 //* Add the daughter momentum to the particle momentum
329
330 KFPart.r[ 3] += m[ 3];
331 KFPart.r[ 4] += m[ 4];
332 KFPart.r[ 5] += m[ 5];
333 KFPart.r[ 6] += m[ 6];
334
335 KFPart.C[ 9] += V[ 9];
336 KFPart.C[13] += V[13];
337 KFPart.C[14] += V[14];
338 KFPart.C[18] += V[18];
339 KFPart.C[19] += V[19];
340 KFPart.C[20] += V[20];
341 KFPart.C[24] += V[24];
342 KFPart.C[25] += V[25];
343 KFPart.C[26] += V[26];
344 KFPart.C[27] += V[27];
345
346 //* CHt = CH' - D'
347
348 fvec CHt0[7], CHt1[7], CHt2[7];
349
350 CHt0[0]=KFPart.C[ 0] ; CHt1[0]=KFPart.C[ 1] ; CHt2[0]=KFPart.C[ 3] ;
351 CHt0[1]=KFPart.C[ 1] ; CHt1[1]=KFPart.C[ 2] ; CHt2[1]=KFPart.C[ 4] ;
352 CHt0[2]=KFPart.C[ 3] ; CHt1[2]=KFPart.C[ 4] ; CHt2[2]=KFPart.C[ 5] ;
353 CHt0[3]=KFPart.C[ 6]-V[ 6]; CHt1[3]=KFPart.C[ 7]-V[ 7]; CHt2[3]=KFPart.C[ 8]-V[ 8];
354 CHt0[4]=KFPart.C[10]-V[10]; CHt1[4]=KFPart.C[11]-V[11]; CHt2[4]=KFPart.C[12]-V[12];
355 CHt0[5]=KFPart.C[15]-V[15]; CHt1[5]=KFPart.C[16]-V[16]; CHt2[5]=KFPart.C[17]-V[17];
356 CHt0[6]=KFPart.C[21]-V[21]; CHt1[6]=KFPart.C[22]-V[22]; CHt2[6]=KFPart.C[23]-V[23];
357
358 //* Kalman gain K = CH'*S
359
360 fvec K0[7], K1[7], K2[7];
361
362 for(Int_t i=0;i<7;++i){
363 K0[i] = CHt0[i]*S[0] + CHt1[i]*S[1] + CHt2[i]*S[3];
364 K1[i] = CHt0[i]*S[1] + CHt1[i]*S[2] + CHt2[i]*S[4];
365 K2[i] = CHt0[i]*S[3] + CHt1[i]*S[4] + CHt2[i]*S[5];
366 }
367
368 //* New estimation of the vertex position r += K*zeta
369
370 for(Int_t i=0;i<7;++i)
371 KFPart.r[i] += K0[i]*zeta[0] + K1[i]*zeta[1] + K2[i]*zeta[2];
372
373 //* New covariance matrix C -= K*(CH')'
374
375 for(Int_t i=0, k=0;i<7;++i){
376 for(Int_t j=0;j<=i;++j,++k)
377 KFPart.C[k] -= K0[i]*CHt0[j] + K1[i]*CHt1[j] + K2[i]*CHt2[j];
378 }
379
380 //* Calculate Chi^2 & NDF
381
382 KFPart.Chi2 += (S[0]*zeta[0] + S[1]*zeta[1] + S[3]*zeta[2])*zeta[0]
383 + (S[1]*zeta[0] + S[2]*zeta[1] + S[4]*zeta[2])*zeta[1]
384 + (S[3]*zeta[0] + S[4]*zeta[1] + S[5]*zeta[2])*zeta[2];
385
386 } // add tracks
387
388 if( iteration==0 ){
389 r0[3] = KFPart.r[3];
390 r0[4] = KFPart.r[4];
391 r0[5] = KFPart.r[5];
392 r0[6] = KFPart.r[6];
393 if( Parent ){
394 fvec dx, dy, dz;
395 dx = Parent->GetRefX()-KFPart.r[0];
396 dy = Parent->GetRefY()-KFPart.r[1];
397 dz = Parent->GetRefZ()-KFPart.r[2];
398
399 r0[7] = KFPart.r[7] = sqrt( (dx*dx+dy*dy+dz*dz)
400 /(KFPart.r[3]*KFPart.r[3]+KFPart.r[4]*KFPart.r[4]+KFPart.r[5]*KFPart.r[5]) );
401 }
402 }
403
404 if( Mass>=0 )
405 {
406 fvec vMass = Mass;
407 MeasureMass( &KFPart, r0, vMass );
408 }
409 if( Parent ) MeasureProductionVertex( &KFPart, r0, Parent );
410 }// iterations
411
412 // calculate a field region for the constructed particle
413 L1FieldValue field[3];
414 fvec zField[3] = {0, KFPart.r[2]/2, KFPart.r[2]};
415
416 for(int iPoint=0; iPoint<3; iPoint++)
417 {
418 for(int iD=0; iD<NDaughters; ++iD)
419 {
420 L1FieldValue b = vDaughters[iD].fField.Get(zField[iPoint]);
421 field[iPoint].x += b.x;
422 field[iPoint].y += b.y;
423 field[iPoint].z += b.z;
424 }
425 field[iPoint].x /= NDaughters;
426 field[iPoint].y /= NDaughters;
427 field[iPoint].z /= NDaughters;
428 }
429
430 KFPart.fField.Set( field[2], zField[2], field[1], zField[1], field[0], zField[0] );
431
432 KFPart.AtProductionVertex = 0;
433}
434
435void CbmKFParticleInterface::MeasureMass(CbmKFParticle_simd* Particle, fvec r0[], fvec Mass )
436{
437 fvec *r = Particle->GetParameters();
438 fvec *C = Particle->GetCovMatrix();
439
440 fvec H[8] = {0};
441 H[0] = H[1] = H[2] = 0.;
442 H[3] = -2*r0[3];
443 H[4] = -2*r0[4];
444 H[5] = -2*r0[5];
445 H[6] = 2*r0[6];
446 H[7] = 0;
447 fvec m2 = r0[6]*r0[6] - r0[3]*r0[3] - r0[4]*r0[4] - r0[5]*r0[5];
448
449 fvec zeta = Mass*Mass - m2;
450 for(Int_t i=0;i<8;++i) zeta -= H[i]*(r[i]-r0[i]);
451
452 fvec S = 0.;
453 fvec CHt[8];
454 for (Int_t i=0;i<8;++i ){
455 CHt[i] = 0.0;
456 for (Int_t j=0;j<8;++j) CHt[i] += Cij(i,j)*H[j];
457 S += H[i]*CHt[i];
458 }
459
460 fvec init = fvec(S > 1.E-20);
461 S = rcp(S) & init;
462
463 Particle->Chi2 += zeta*zeta*S;
464 Particle->NDF += (fvec(1.) & init);
465 for( Int_t i=0, ii=0; i<8; ++i ){
466 fvec Ki = CHt[i]*S;
467 r[i] += Ki*zeta;
468 for(Int_t j=0;j<=i;++j) C[ii++] -= Ki*CHt[j];
469 }
470}
471
472
473void CbmKFParticleInterface::MeasureProductionVertex( CbmKFParticle_simd* Particle, fvec r0[],
474 CbmKFVertexInterface *Parent)
475{
476 fvec *r = Particle->GetParameters();
477 fvec *C = Particle->GetCovMatrix();
478
479 fvec m[3], V[6];
480
481 m[0] = Parent->GetRefX();
482 m[1] = Parent->GetRefY();
483 m[2] = Parent->GetRefZ();
484 double *V1 = Parent->GetCovMatrix();
485 for(int i=0; i<6; i++) V[i]=V1[i];
486
487 r[7] = r0[7];
488 Extrapolate(r0, -r0[7]);
489 Convert(Particle, r0, 1);
490
491 fvec Ai[6];
492 Ai[0] = C[4]*C[4] - C[2]*C[5];
493 Ai[1] = C[1]*C[5] - C[3]*C[4];
494 Ai[3] = C[2]*C[3] - C[1]*C[4];
495 fvec det = 1./(C[0]*Ai[0] + C[1]*Ai[1] + C[3]*Ai[3]);
496 Ai[0] *= det;
497 Ai[1] *= det;
498 Ai[3] *= det;
499 Ai[2] = ( C[3]*C[3] - C[0]*C[5] )*det;
500 Ai[4] = ( C[0]*C[4] - C[1]*C[3] )*det;
501 Ai[5] = ( C[1]*C[1] - C[0]*C[2] )*det;
502
503 fvec B[5][3];
504
505 B[0][0] = C[ 6]*Ai[0] + C[ 7]*Ai[1] + C[ 8]*Ai[3];
506 B[0][1] = C[ 6]*Ai[1] + C[ 7]*Ai[2] + C[ 8]*Ai[4];
507 B[0][2] = C[ 6]*Ai[3] + C[ 7]*Ai[4] + C[ 8]*Ai[5];
508
509 B[1][0] = C[10]*Ai[0] + C[11]*Ai[1] + C[12]*Ai[3];
510 B[1][1] = C[10]*Ai[1] + C[11]*Ai[2] + C[12]*Ai[4];
511 B[1][2] = C[10]*Ai[3] + C[11]*Ai[4] + C[12]*Ai[5];
512
513 B[2][0] = C[15]*Ai[0] + C[16]*Ai[1] + C[17]*Ai[3];
514 B[2][1] = C[15]*Ai[1] + C[16]*Ai[2] + C[17]*Ai[4];
515 B[2][2] = C[15]*Ai[3] + C[16]*Ai[4] + C[17]*Ai[5];
516
517 B[3][0] = C[21]*Ai[0] + C[22]*Ai[1] + C[23]*Ai[3];
518 B[3][1] = C[21]*Ai[1] + C[22]*Ai[2] + C[23]*Ai[4];
519 B[3][2] = C[21]*Ai[3] + C[22]*Ai[4] + C[23]*Ai[5];
520
521 B[4][0] = C[28]*Ai[0] + C[29]*Ai[1] + C[30]*Ai[3];
522 B[4][1] = C[28]*Ai[1] + C[29]*Ai[2] + C[30]*Ai[4];
523 B[4][2] = C[28]*Ai[3] + C[29]*Ai[4] + C[30]*Ai[5];
524
525 fvec z[3] = { m[0]-r[0], m[1]-r[1], m[2]-r[2] };
526
527 {
528 fvec AV[6] = { C[0]-V[0], C[1]-V[1], C[2]-V[2],
529 C[3]-V[3], C[4]-V[4], C[5]-V[5] };
530 fvec AVi[6];
531 AVi[0] = AV[4]*AV[4] - AV[2]*AV[5];
532 AVi[1] = AV[1]*AV[5] - AV[3]*AV[4];
533 AVi[2] = AV[3]*AV[3] - AV[0]*AV[5] ;
534 AVi[3] = AV[2]*AV[3] - AV[1]*AV[4];
535 AVi[4] = AV[0]*AV[4] - AV[1]*AV[3] ;
536 AVi[5] = AV[1]*AV[1] - AV[0]*AV[2] ;
537
538 det = 1./( AV[0]*AVi[0] + AV[1]*AVi[1] + AV[3]*AVi[3] );
539
540 Particle->NDF += 2.;
541 Particle->Chi2 +=
542 ( +(AVi[0]*z[0] + AVi[1]*z[1] + AVi[3]*z[2])*z[0]
543 +(AVi[1]*z[0] + AVi[2]*z[1] + AVi[4]*z[2])*z[1]
544 +(AVi[3]*z[0] + AVi[4]*z[1] + AVi[5]*z[2])*z[2] )*det;
545 }
546
547 r[0] = m[0];
548 r[1] = m[1];
549 r[2] = m[2];
550 r[3]+= B[0][0]*z[0] + B[0][1]*z[1] + B[0][2]*z[2];
551 r[4]+= B[1][0]*z[0] + B[1][1]*z[1] + B[1][2]*z[2];
552 r[5]+= B[2][0]*z[0] + B[2][1]*z[1] + B[2][2]*z[2];
553 r[6]+= B[3][0]*z[0] + B[3][1]*z[1] + B[3][2]*z[2];
554 r[7]+= B[4][0]*z[0] + B[4][1]*z[1] + B[4][2]*z[2];
555
556 fvec d0, d1, d2;
557
558 C[0] = V[0];
559 C[1] = V[1];
560 C[2] = V[2];
561 C[3] = V[3];
562 C[4] = V[4];
563 C[5] = V[5];
564
565 d0= B[0][0]*V[0] + B[0][1]*V[1] + B[0][2]*V[3] - C[ 6];
566 d1= B[0][0]*V[1] + B[0][1]*V[2] + B[0][2]*V[4] - C[ 7];
567 d2= B[0][0]*V[3] + B[0][1]*V[4] + B[0][2]*V[5] - C[ 8];
568
569 C[ 6]+= d0;
570 C[ 7]+= d1;
571 C[ 8]+= d2;
572 C[ 9]+= d0*B[0][0] + d1*B[0][1] + d2*B[0][2];
573
574 d0= B[1][0]*V[0] + B[1][1]*V[1] + B[1][2]*V[3] - C[10];
575 d1= B[1][0]*V[1] + B[1][1]*V[2] + B[1][2]*V[4] - C[11];
576 d2= B[1][0]*V[3] + B[1][1]*V[4] + B[1][2]*V[5] - C[12];
577
578 C[10]+= d0;
579 C[11]+= d1;
580 C[12]+= d2;
581 C[13]+= d0*B[0][0] + d1*B[0][1] + d2*B[0][2];
582 C[14]+= d0*B[1][0] + d1*B[1][1] + d2*B[1][2];
583
584 d0= B[2][0]*V[0] + B[2][1]*V[1] + B[2][2]*V[3] - C[15];
585 d1= B[2][0]*V[1] + B[2][1]*V[2] + B[2][2]*V[4] - C[16];
586 d2= B[2][0]*V[3] + B[2][1]*V[4] + B[2][2]*V[5] - C[17];
587
588 C[15]+= d0;
589 C[16]+= d1;
590 C[17]+= d2;
591 C[18]+= d0*B[0][0] + d1*B[0][1] + d2*B[0][2];
592 C[19]+= d0*B[1][0] + d1*B[1][1] + d2*B[1][2];
593 C[20]+= d0*B[2][0] + d1*B[2][1] + d2*B[2][2];
594
595 d0= B[3][0]*V[0] + B[3][1]*V[1] + B[3][2]*V[3] - C[21];
596 d1= B[3][0]*V[1] + B[3][1]*V[2] + B[3][2]*V[4] - C[22];
597 d2= B[3][0]*V[3] + B[3][1]*V[4] + B[3][2]*V[5] - C[23];
598
599 C[21]+= d0;
600 C[22]+= d1;
601 C[23]+= d2;
602 C[24]+= d0*B[0][0] + d1*B[0][1] + d2*B[0][2];
603 C[25]+= d0*B[1][0] + d1*B[1][1] + d2*B[1][2];
604 C[26]+= d0*B[2][0] + d1*B[2][1] + d2*B[2][2];
605 C[27]+= d0*B[3][0] + d1*B[3][1] + d2*B[3][2];
606
607 d0= B[4][0]*V[0] + B[4][1]*V[1] + B[4][2]*V[3] - C[28];
608 d1= B[4][0]*V[1] + B[4][1]*V[2] + B[4][2]*V[4] - C[29];
609 d2= B[4][0]*V[3] + B[4][1]*V[4] + B[4][2]*V[5] - C[30];
610
611 C[28]+= d0;
612 C[29]+= d1;
613 C[30]+= d2;
614 C[31]+= d0*B[0][0] + d1*B[0][1] + d2*B[0][2];
615 C[32]+= d0*B[1][0] + d1*B[1][1] + d2*B[1][2];
616 C[33]+= d0*B[2][0] + d1*B[2][1] + d2*B[2][2];
617 C[34]+= d0*B[3][0] + d1*B[3][1] + d2*B[3][2];
618 C[35]+= d0*B[4][0] + d1*B[4][1] + d2*B[4][2];
619
620 Extrapolate(Particle, r0, r[7]);
621 Convert(Particle, r0, 0);
622}
623
624
625void CbmKFParticleInterface::Convert( CbmKFParticle_simd* Particle, fvec r0[], bool ToProduction )
626{
627 fvec *C = Particle->GetCovMatrix();
628
629 L1FieldValue B = KFPart.fField.Get(r0[2]);
630 B.x *= 0.000299792458; //multiply by c
631 B.y *= 0.000299792458;
632 B.z *= 0.000299792458;
633
634 fvec h[6];
635
636 h[0] = r0[3];
637 h[1] = r0[4];
638 h[2] = r0[5];
639 if( ToProduction ){ h[0]=-h[0]; h[1]=-h[1]; h[2]=-h[2]; }
640 h[3] = h[1]*B.z-h[2]*B.y;
641 h[4] = h[2]*B.x-h[0]*B.z;
642 h[5] = h[0]*B.y-h[1]*B.x;
643
644 fvec c;
645
646 c = C[28]+h[0]*C[35];
647 C[ 0]+= h[0]*(c+C[28]);
648 C[28] = c;
649
650 C[ 1]+= h[1]*C[28] + h[0]*C[29];
651 c = C[29]+h[1]*C[35];
652 C[ 2]+= h[1]*(c+C[29]);
653 C[29] = c;
654
655 C[ 3]+= h[2]*C[28] + h[0]*C[30];
656 C[ 4]+= h[2]*C[29] + h[1]*C[30];
657 c = C[30]+h[2]*C[35];
658 C[ 5]+= h[2]*(c+C[30]);
659 C[30] = c;
660
661 C[ 6]+= h[3]*C[28] + h[0]*C[31];
662 C[ 7]+= h[3]*C[29] + h[1]*C[31];
663 C[ 8]+= h[3]*C[30] + h[2]*C[31];
664 c = C[31]+h[3]*C[35];
665 C[ 9]+= h[3]*(c+C[31]);
666 C[31] = c;
667
668 C[10]+= h[4]*C[28] + h[0]*C[32];
669 C[11]+= h[4]*C[29] + h[1]*C[32];
670 C[12]+= h[4]*C[30] + h[2]*C[32];
671 C[13]+= h[4]*C[31] + h[3]*C[32];
672 c = C[32]+h[4]*C[35];
673 C[14]+= h[4]*(c+C[32]);
674 C[32] = c;
675
676 C[15]+= h[5]*C[28] + h[0]*C[33];
677 C[16]+= h[5]*C[29] + h[1]*C[33];
678 C[17]+= h[5]*C[30] + h[2]*C[33];
679 C[18]+= h[5]*C[31] + h[3]*C[33];
680 C[19]+= h[5]*C[32] + h[4]*C[33];
681 c = C[33]+h[5]*C[35];
682 C[20]+= h[5]*(c+C[33]);
683 C[33] = c;
684
685 C[21]+= h[0]*C[34];
686 C[22]+= h[1]*C[34];
687 C[23]+= h[2]*C[34];
688 C[24]+= h[3]*C[34];
689 C[25]+= h[4]*C[34];
690 C[26]+= h[5]*C[34];
691}
692
693void CbmKFParticleInterface::ExtrapolateLine(CbmKFParticle_simd* Particle, fvec r0[], fvec dS )
694{
695 fvec *r = Particle->GetParameters();
696 fvec *C = Particle->GetCovMatrix();
697
698 if( r0 && r0!=r ){
699 r0[0]+= dS*r0[3];
700 r0[1]+= dS*r0[4];
701 r0[2]+= dS*r0[5];
702 }
703
704 r[0]+= dS*r[3];
705 r[1]+= dS*r[4];
706 r[2]+= dS*r[5];
707
708 fvec C6 = C[ 6] + dS*C[ 9];
709 fvec C11 = C[11] + dS*C[14];
710 fvec C17 = C[17] + dS*C[20];
711 fvec SC13 = dS*C[13];
712 fvec SC18 = dS*C[18];
713 fvec SC19 = dS*C[19];
714
715 C[ 0] += dS*( C[ 6] + C6 );
716 C[ 2] += dS*( C[11] + C11 );
717 C[ 5] += dS*( C[17] + C17 );
718
719 C[ 7] += SC13;
720 C[ 8] += SC18;
721 C[12] += SC19;
722
723 C[ 1] += dS*( C[10] + C[ 7] );
724 C[ 3] += dS*( C[15] + C[ 8] );
725 C[ 4] += dS*( C[16] + C[12] );
726
727 C[ 6] = C6;
728 C[10] += SC13;
729 C[11] = C11;
730 C[15] += SC18;
731 C[16] += SC19;
732 C[17] = C17;
733
734 C[21] += dS*C[24];
735 C[22] += dS*C[25];
736 C[23] += dS*C[26];
737 C[28] += dS*C[31];
738 C[29] += dS*C[32];
739 C[30] += dS*C[33];
740}
741
742void CbmKFParticleInterface::Extrapolate(CbmKFParticle_simd* Particle, fvec r0[], fvec dS )
743{
744 if( Particle->Q[0]==0 ){
745 ExtrapolateLine( Particle, r0, dS );
746 return;
747 }
748
749 fvec *r = Particle->GetParameters();
750 fvec *C = Particle->GetCovMatrix();
751
752 const float c_light = 0.000299792458;
753
754 fvec c = Particle->Q*c_light;
755
756 // construct coefficients
757
758 fvec
759 px = r[3],
760 py = r[4],
761 pz = r[5];
762
763 fvec sx=0, sy=0, sz=0, syy=0, syz=0, syyy=0, Sx=0, Sy=0, Sz=0, Syy=0, Syz=0, Syyy=0;
764
765 { // get field integrals
766
767 L1FieldValue B[3];
768 fvec p0[3], p1[3], p2[3];
769
770 // line track approximation
771
772 p0[0] = r[0];
773 p0[1] = r[1];
774 p0[2] = r[2];
775
776 p2[0] = r[0] + px*dS;
777 p2[1] = r[1] + py*dS;
778 p2[2] = r[2] + pz*dS;
779
780 p1[0] = 0.5*(p0[0]+p2[0]);
781 p1[1] = 0.5*(p0[1]+p2[1]);
782 p1[2] = 0.5*(p0[2]+p2[2]);
783
784 B[0] = Particle->fField.Get(p0[2]);
785 B[1] = Particle->fField.Get(p1[2]);
786 B[2] = Particle->fField.Get(p2[2]);
787
788 // first order track approximation
789
790 fvec Sy1 = ( 7*B[0].y + 6*B[1].y-B[2].y )*c*dS*dS/96.;
791 fvec Sy2 = ( B[0].y + 2*B[1].y )*c*dS*dS/6.;
792
793 p1[0] -= Sy1*pz;
794 p1[2] += Sy1*px;
795 p2[0] -= Sy2*pz;
796 p2[2] += Sy2*px;
797
798 B[0] = Particle->fField.Get(p0[2]);
799 B[1] = Particle->fField.Get(p1[2]);
800 B[2] = Particle->fField.Get(p2[2]);
801
802 sx = c*( B[0].x + 4*B[1].x + B[2].x )*dS/6.;
803 sy = c*( B[0].y + 4*B[1].y + B[2].y )*dS/6.;
804 sz = c*( B[0].z + 4*B[1].z + B[2].z )*dS/6.;
805
806 Sx = c*( B[0].x + 2*B[1].x)*dS*dS/6.;
807 Sy = c*( B[0].y + 2*B[1].y)*dS*dS/6.;
808 Sz = c*( B[0].z + 2*B[1].z)*dS*dS/6.;
809
810 fvec c2[3][3] = { { 5, -4, -1},{ 44, 80, -4},{ 11, 44, 5} }; // /=360.
811 fvec C2[3][3] = { { 38, 8, -4},{ 148, 208, -20},{ 3, 36, 3} }; // /=2520.
812 for(Int_t n=0; n<3; n++)
813 for(Int_t m=0; m<3; m++)
814 {
815 syz += c2[n][m]*B[n].y*B[m].z;
816 Syz += C2[n][m]*B[n].y*B[m].z;
817 }
818
819 syz *= c*c*dS*dS/360.;
820 Syz *= c*c*dS*dS*dS/2520.;
821
822 syy = c*( B[0].y + 4*B[1].y + B[2].y )*dS;
823 syyy = syy*syy*syy / 1296;
824 syy = syy*syy/72;
825
826 Syy = ( B[0].y*( 38*B[0].y + 156*B[1].y - B[2].y )+
827 B[1].y*( 208*B[1].y +16*B[2].y )+
828 B[2].y*( 3*B[2].y )
829 )*dS*dS*dS*c*c/2520.;
830 Syyy =
831 (
832 B[0].y*( B[0].y*( 85*B[0].y + 526*B[1].y - 7*B[2].y )+
833 B[1].y*( 1376*B[1].y +84*B[2].y )+
834 B[2].y*( 19*B[2].y ) )+
835 B[1].y*( B[1].y*( 1376*B[1].y +256*B[2].y )+
836 B[2].y*( 62*B[2].y ) )+
837 B[2].y*B[2].y *( 3*B[2].y )
838 )*dS*dS*dS*dS*c*c*c/90720.;
839 }
840
841 fvec J[11];
842// for( int i=0; i<8; i++ ) for( int j=0; j<8; j++) J[i][j]=0;
843//
844// J[0][0]=1; J[0][1]=0; J[0][2]=0; J[0][3]=dS-Syy; J[0][4]=Sx; J[0][5]=Syyy-Sy;
845// J[1][0]=0; J[1][1]=1; J[1][2]=0; J[1][3]=-Sz; J[1][4]=dS; J[1][5]=Sx+Syz;
846// J[2][0]=0; J[2][1]=0; J[2][2]=1; J[2][3]=Sy-Syyy; J[2][4]=-Sx; J[2][5]=dS-Syy;
847//
848// J[3][0]=0; J[3][1]=0; J[3][2]=0; J[3][3]=1-syy; J[3][4]=sx; J[3][5]=syyy-sy;
849// J[4][0]=0; J[4][1]=0; J[4][2]=0; J[4][3]=-sz; J[4][4]=1; J[4][5]=sx+syz;
850// J[5][0]=0; J[5][1]=0; J[5][2]=0; J[5][3]=sy-syyy; J[5][4]=-sx; J[5][5]=1-syy;
851// J[6][6] = J[7][7] = 1;
852
853
854 J[0]=dS-Syy; J[1]=Sx; J[2]=Syyy-Sy;
855 J[3]=-Sz; J[4]=dS; J[5]=Sx+Syz;
856
857 J[6]=1-syy; J[7]=sx; J[8]=syyy-sy;
858 J[9]=-sz; J[10]=sx+syz;
859
860
861
862 r[0]+= J[0]*px + J[1]*py + J[2]*pz;
863 r[1]+= J[3]*px + J[4]*py + J[5]*pz;
864 r[2]+= -J[2]*px - J[1]*py + J[0]*pz;
865 r[3] = J[6]*px + J[7]*py + J[8]*pz;
866 r[4] = J[9]*px + py + J[10]*pz;
867 r[5] = -J[8]*px - J[7]*py + J[6]*pz;
868
869 if( r0 && r0!=r ){
870 px = r0[3];
871 py = r0[4];
872 pz = r0[5];
873
874 r0[0]+= J[0]*px + J[1]*py + J[2]*pz;
875 r0[1]+= J[3]*px + J[4]*py + J[5]*pz;
876 r0[2]+= -J[2]*px - J[1]*py + J[0]*pz;
877 r0[3] = J[6]*px + J[7]*py + J[8]*pz;
878 r0[4] = J[9]*px + py + J[10]*pz;
879 r0[5] = -J[8]*px - J[7]*py + J[6]*pz;
880 }
881
882 multQSQt1( J, C);
883}
884
885void CbmKFParticleInterface::multQSQt( const fvec Q[], fvec S[] )
886{
887 fvec A[64]={0};
888
889 for( Int_t i=0, n=0; i<8; i++ )
890 for( Int_t j=0; j<8; j++, ++n )
891 for( Int_t k=0; k<8; ++k ) A[n]+= S[IJ(i,k)] * Q[ j*8+k];
892
893 for( Int_t i=0; i<8; i++ )
894 {
895 for( Int_t j=0; j<=i; j++ )
896 {
897 Int_t n = IJ(i,j);
898 S[n] = 0 ;
899 for( Int_t k=0; k<8; k++ ) S[n] += Q[ i*8+k ] * A[ k*8+j ];
900 }
901 }
902}
903
904void CbmKFParticleInterface::multQSQt1( const fvec J[11], fvec S[] )
905{
906 const fvec A00 = S[0 ]+S[6 ]*J[0]+S[10]*J[1]+S[15]*J[2];
907 const fvec A10 = S[1 ]+S[7 ]*J[0]+S[11]*J[1]+S[16]*J[2];
908 const fvec A20 = S[3 ]+S[8 ]*J[0]+S[12]*J[1]+S[17]*J[2];
909 const fvec A30 = S[6 ]+S[9 ]*J[0]+S[13]*J[1]+S[18]*J[2];
910 const fvec A40 = S[10]+S[13]*J[0]+S[14]*J[1]+S[19]*J[2];
911 const fvec A50 = S[15]+S[18]*J[0]+S[19]*J[1]+S[20]*J[2];
912 const fvec A60 = S[21]+S[24]*J[0]+S[25]*J[1]+S[26]*J[2];
913 const fvec A70 = S[28]+S[31]*J[0]+S[32]*J[1]+S[33]*J[2];
914
915 S[0 ] = A00+J[0]*A30+J[1]*A40+J[2]*A50;
916 S[1 ] = A10+J[3]*A30+J[4]*A40+J[5]*A50;
917 S[3 ] = A20-J[2]*A30-J[1]*A40+J[0]*A50;
918 S[6 ] = J[6]*A30+J[7]*A40+J[8]*A50;
919 S[10] = J[9]*A30+ A40+J[10]*A50;
920 S[15] = -J[8]*A30-J[7]*A40+J[6]*A50;
921 S[21] = A60;
922 S[28] = A70;
923
924// const fvec A01 = S[1 ]+S[6 ]*J[3]+S[10]*J[4]+S[15]*J[5];
925 const fvec A11 = S[2 ]+S[7 ]*J[3]+S[11]*J[4]+S[16]*J[5];
926 const fvec A21 = S[4 ]+S[8 ]*J[3]+S[12]*J[4]+S[17]*J[5];
927 const fvec A31 = S[7 ]+S[9 ]*J[3]+S[13]*J[4]+S[18]*J[5];
928 const fvec A41 = S[11]+S[13]*J[3]+S[14]*J[4]+S[19]*J[5];
929 const fvec A51 = S[16]+S[18]*J[3]+S[19]*J[4]+S[20]*J[5];
930 const fvec A61 = S[22]+S[24]*J[3]+S[25]*J[4]+S[26]*J[5];
931 const fvec A71 = S[29]+S[31]*J[3]+S[32]*J[4]+S[33]*J[5];
932
933 S[2 ] = A11+J[3]*A31+J[4]*A41+J[5]*A51;
934 S[4 ] = A21-J[2]*A31-J[1]*A41+J[0]*A51;
935 S[7 ] = J[6]*A31+J[7]*A41+J[8]*A51;
936 S[11] = J[9]*A31+ A41+J[10]*A51;
937 S[16] = -J[8]*A31-J[7]*A41+J[6]*A51;
938 S[22] = A61;
939 S[29] = A71;
940
941// const fvec A02 = S[3 ]+S[6 ]*J[2][3]+S[10]*J[2][4]+S[15]*J[2][5];
942// const fvec A12 = S[4 ]+S[7 ]*J[2][3]+S[11]*J[2][4]+S[16]*J[2][5];
943 const fvec A22 = S[5 ]-S[8 ]*J[2]-S[12]*J[1]+S[17]*J[0];
944 const fvec A32 = S[8 ]-S[9 ]*J[2]-S[13]*J[1]+S[18]*J[0];
945 const fvec A42 = S[12]-S[13]*J[2]-S[14]*J[1]+S[19]*J[0];
946 const fvec A52 = S[17]-S[18]*J[2]-S[19]*J[1]+S[20]*J[0];
947 const fvec A62 = S[23]-S[24]*J[2]-S[25]*J[1]+S[26]*J[0];
948 const fvec A72 = S[30]-S[31]*J[2]-S[32]*J[1]+S[33]*J[0];
949
950 S[5 ] = A22-J[2]*A32-J[1]*A42+J[0]*A52;
951 S[8 ] = J[6]*A32+J[7]*A42+J[8]*A52;
952 S[12] = J[9]*A32+ A42+J[10]*A52;
953 S[17] = -J[8]*A32-J[7]*A42+J[6]*A52;
954 S[23] = A62;
955 S[30] = A72;
956
957// const fvec A03 = S[6] *J[6]+S[10]*J[7]+S[15]*J[8];
958// const fvec A13 = S[7] *J[6]+S[11]*J[7]+S[16]*J[8];
959// const fvec A23 = S[8] *J[6]+S[12]*J[7]+S[17]*J[8];
960 const fvec A33 = S[9] *J[6]+S[13]*J[7]+S[18]*J[8];
961 const fvec A43 = S[13]*J[6]+S[14]*J[7]+S[19]*J[8];
962 const fvec A53 = S[18]*J[6]+S[19]*J[7]+S[20]*J[8];
963 const fvec A63 = S[24]*J[6]+S[25]*J[7]+S[26]*J[8];
964 const fvec A73 = S[31]*J[6]+S[32]*J[7]+S[33]*J[8];
965
966// const fvec A04 = S[6] *J[9]+S[10]*J[4][4]+S[15]*J[10];
967// const fvec A14 = S[7] *J[9]+S[11]*J[4][4]+S[16]*J[10];
968// const fvec A24 = S[8] *J[9]+S[12]*J[4][4]+S[17]*J[10];
969 const fvec A34 = S[9] *J[9]+S[13]+S[18]*J[10];
970 const fvec A44 = S[13]*J[9]+S[14]+S[19]*J[10];
971 const fvec A54 = S[18]*J[9]+S[19]+S[20]*J[10];
972 const fvec A64 = S[24]*J[9]+S[25]+S[26]*J[10];
973 const fvec A74 = S[31]*J[9]+S[32]+S[33]*J[10];
974
975// const fvec A05 = S[6] *J[5][3]+S[10]*J[5][4]+S[15]*J[6];
976// const fvec A15 = S[7] *J[5][3]+S[11]*J[5][4]+S[16]*J[6];
977// const fvec A25 = S[8] *J[5][3]+S[12]*J[5][4]+S[17]*J[6];
978 const fvec A35 = -S[9] *J[8]-S[13]*J[7]+S[18]*J[6];
979 const fvec A45 = -S[13]*J[8]-S[14]*J[7]+S[19]*J[6];
980 const fvec A55 = -S[18]*J[8]-S[19]*J[7]+S[20]*J[6];
981 const fvec A65 = -S[24]*J[8]-S[25]*J[7]+S[26]*J[6];
982 const fvec A75 = -S[31]*J[8]-S[32]*J[7]+S[33]*J[6];
983
984
985 S[9 ] = J[6]*A33+J[7]*A43+J[8]*A53;
986 S[13] = J[9]*A33+ A43+J[10]*A53;
987 S[18] = -J[8]*A33-J[7]*A43+J[6]*A53;
988 S[24] = A63;
989 S[31] = A73;
990
991 S[14] = J[9]*A34+ A44+J[10]*A54;
992 S[19] = -J[8]*A34-J[7]*A44+J[6]*A54;
993 S[25] = A64;
994 S[32] = A74;
995
996 S[20] = -J[8]*A35-J[7]*A45+J[6]*A55;
997 S[26] = A65;
998 S[33] = A75;
999
1000//S[27] = S27;
1001//S[34] = S34;
1002//S[35] = S35;
1003}
1004
1005void CbmKFParticleInterface::Extrapolate( fvec r0[], fvec dS )
1006{
1007 Extrapolate(&KFPart, r0, dS);
1008}
1009
1010void CbmKFParticleInterface::TransportToProductionVertex()
1011{
1012 if( KFPart.AtProductionVertex ) return;
1013 fvec r0[8];
1014 for( int i=0; i<8; i++ ) r0[i] = KFPart.r[i];
1015 Extrapolate(r0, -(KFPart.r[7]));
1016 Convert(&KFPart, r0, 1);
1017 KFPart.AtProductionVertex = 1;
1018}
1019
1020void CbmKFParticleInterface::TransportToDecayVertex()
1021{
1022 if( !KFPart.AtProductionVertex ) return;
1023 fvec r0[8];
1024 for( int i=0; i<8; i++ ) r0[i] = KFPart.r[i];
1025 Extrapolate(r0, KFPart.r[7]);
1026 Convert(&KFPart, r0, 0);
1027 KFPart.AtProductionVertex = 0;
1028}
1029
1030fvec CbmKFParticleInterface::GetDStoPoint( const CbmKFParticle_simd& Particle, const fvec xyz[] ) const
1031{
1032// const fvec *r = Particle.GetParameters();
1033//
1034// const fvec dx = xyz[0] - r[0];
1035// const fvec dy = xyz[1] - r[1];
1036// const fvec dz = xyz[2] - r[2];
1037// const fvec p2 = r[3]*r[3] + r[4]*r[4]+r[5]*r[5];
1038//
1039// L1FieldValue B = const_cast<L1FieldRegion*>(&Particle.fField)->Get(r[2]);
1040// B.x *= 0.000299792458; //multiply by c
1041// B.y *= 0.000299792458;
1042// B.z *= 0.000299792458;
1043//
1044// fvec bq1 = B.x*Particle.Q;
1045// fvec bq2 = B.y*Particle.Q;
1046// fvec bq3 = B.z*Particle.Q;
1047// fvec a = dx*r[3]+dy*r[4]+dz*r[5];
1048// fvec B2 = B.x*B.x+B.y*B.y+B.z*B.z;
1049// fvec bq = Particle.Q*sqrt(B2);
1050// fvec pB = r[3]*B.x+r[4]*B.y+r[5]*B.z;
1051// fvec rB = dx*B.x+dy*B.y+dz*B.z;
1052//
1053// fvec small = 0.0001;
1054// fvec init = fvec(B2 > 0.0001);
1055// B2 = (init & B2) + (fvec(!init) & small);
1056// fvec pt2 = p2 - pB*pB/B2;
1057// a = a - pB*rB/B2;
1058//
1059// fvec dS = atan2(bq*a, pt2 + bq1*(dz*r[4]-dy*r[5]) - bq2*(dz*r[3]-dx*r[5])+bq3*(dy*r[3]-dx*r[4]))/bq;
1060// fvec dS2 = a/pt2;
1061// fvec smallBq = 1.e-8;
1062// fvec initBq = fvec(fabs(bq) > smallBq);
1063// dS = (initBq & dS) + (fvec(!initBq) & dS2);
1064// fvec initPt2 = fvec(pt2 > small);
1065// dS = (initPt2 & dS) + (fvec(!initPt2) & fvec(0));
1066//
1067// //fvec dSm = rB/pB;
1068// return dS;
1069
1070 const fvec *r = Particle.GetParameters();
1071 fvec p2 = r[3]*r[3] + r[4]*r[4] + r[5]*r[5];
1072// if( p2<1.e-4 ) p2 = 1;
1073 return ( r[3]*(xyz[0]-r[0]) + r[4]*(xyz[1]-r[1]) + r[5]*(xyz[2]-r[2]) )/p2;
1074}
1075
1076fvec CbmKFParticleInterface::GetDStoPoint( const fvec xyz[] ) const
1077{
1078 return GetDStoPoint( KFPart, xyz );
1079}
1080
1081void CbmKFParticleInterface::GetKFVertex( CbmKFVertex vtx[] )
1082{
1083 float r0[fvecLen], r1[fvecLen], r2[fvecLen];
1084 float C_temp[6][fvecLen];
1085 float Chi2_temp[fvecLen];
1086 int NDF_temp[fvecLen];
1087
1088 for(int j=0; j<fvecLen; j++)
1089 {
1090 r0[j] = KFPart.r[0][j];
1091 r1[j] = KFPart.r[1][j];
1092 r2[j] = KFPart.r[2][j];
1093 for( int i=0; i<6; i++) C_temp[i][j] = KFPart.C[i][j];
1094
1095 Chi2_temp[j] = KFPart.Chi2[j];
1096 NDF_temp[j] = (int) KFPart.NDF[j];
1097 }
1098
1099 for(int j=0; j<fvecLen; j++)
1100 vtx[j] = GetKFVertexJ(j);
1101}
1102
1103CbmKFVertex CbmKFParticleInterface::GetKFVertexJ(int j)
1104{
1105 CbmKFVertex vtx;
1106 vtx.GetRefX() = KFPart.r[0][j];
1107 vtx.GetRefY() = KFPart.r[1][j];
1108 vtx.GetRefZ() = KFPart.r[2][j];
1109 for( int i=0; i<6; i++) vtx.GetCovMatrix()[i] = KFPart.C[i][j];
1110 vtx.GetRefChi2() = KFPart.Chi2[j];
1111 vtx.GetRefNDF() = (int) KFPart.NDF[j];
1112
1113 return vtx;
1114}
1115
1116void CbmKFParticleInterface::GetKFVertexJ(int j, CbmKFVertex *vtx)
1117{
1118 vtx->GetRefX() = KFPart.r[0][j];
1119 vtx->GetRefY() = KFPart.r[1][j];
1120 vtx->GetRefZ() = KFPart.r[2][j];
1121 for( int i=0; i<6; i++) vtx->GetCovMatrix()[i] = KFPart.C[i][j];
1122 vtx->GetRefChi2() = KFPart.Chi2[j];
1123 vtx->GetRefNDF() = (int) KFPart.NDF[j];
1124}
1125
1126void CbmKFParticleInterface::GetKFParticle(CbmKFParticle &Part, int iPart)
1127{
1128 Part.SetId(static_cast<int>(KFPart.Id()[iPart]));
1129
1130 Part.CleanDaughtersId();
1131 Part.SetNDaughters(KFPart.DaughterIds().size());
1132 for( unsigned int i = 0; i < KFPart.DaughterIds().size(); i++ )
1133 Part.AddDaughter(static_cast<int>(KFPart.DaughterIds()[i][iPart]));
1134
1135 Part.SetPDG( static_cast<int>(KFPart.GetPDG()[iPart]) );
1136
1137 for(int iP=0; iP<8; iP++)
1138 Part.GetParameters()[iP] = KFPart.GetParameters()[iP][iPart];
1139 for(int iC=0; iC<36; iC++)
1140 Part.GetCovMatrix()[iC] = KFPart.GetCovMatrix()[iC][iPart];
1141
1142 Part.NDF = static_cast<int>(KFPart.GetNDF()[iPart]);
1143 Part.Chi2 = KFPart.GetChi2()[iPart];
1144 Part.Q = KFPart.GetQ()[iPart];
1145 Part.AtProductionVertex = KFPart.AtProductionVertex;
1146
1147 Part.fieldRegion[0] = KFPart.fField.cx0[iPart];
1148 Part.fieldRegion[1] = KFPart.fField.cx1[iPart];
1149 Part.fieldRegion[2] = KFPart.fField.cx2[iPart];
1150 Part.fieldRegion[3] = KFPart.fField.cy0[iPart];
1151 Part.fieldRegion[4] = KFPart.fField.cy1[iPart];
1152 Part.fieldRegion[5] = KFPart.fField.cy2[iPart];
1153 Part.fieldRegion[6] = KFPart.fField.cz0[iPart];
1154 Part.fieldRegion[7] = KFPart.fField.cz1[iPart];
1155 Part.fieldRegion[8] = KFPart.fField.cz2[iPart];
1156 Part.fieldRegion[9] = KFPart.fField.z0[iPart];
1157}
1158
1159void CbmKFParticleInterface::GetKFParticle(CbmKFParticle* Part, int nPart)
1160{
1161 for(int i=0; i<nPart; i++)
1162 GetKFParticle(Part[i],i);
1163}
1164
1165void CbmKFParticleInterface::SetVtxGuess( fvec &x, fvec &y, fvec &z )
1166{
1167 KFPart.SetVtxGuess(x,y,z);
1168}
1169
1170void CbmKFParticleInterface::SetVtxErrGuess( fvec &dx, fvec &dy, fvec &dz )
1171{
1172 KFPart.SetVtxErrGuess(dx,dy,dz);
1173}
1174
1175void CbmKFParticleInterface::GetMomentum( fvec &P, fvec &Error )
1176{
1177 KFPart.GetMomentum( P, Error );
1178}
1179
1180void CbmKFParticleInterface::GetMass( fvec &M, fvec &Error )
1181{
1182 KFPart.GetMass ( M, Error );
1183}
1184
1185//template<class T>
1186void CbmKFParticleInterface::FindParticles(vector<CbmKFTrack> &vRTracks, vector<float>& ChiToPrimVtx,
1187 vector<L1FieldRegion>& vField, vector<CbmKFParticle>& Particles,
1188 CbmKFVertex& PrimVtx, const vector<int>& vTrackPDG, const float cuts[2][3])
1189{
1190 //* Finds particles (K0s and Lambda) from a given set of tracks
1191
1192 static const int NTrackTypes = 8;
1193
1194 int pdgPos[NTrackTypes]={-11,-13, 211, 321, 2212, 211, 321, 2212};
1195 int pdgNeg[NTrackTypes]={ 11, 13,-211,-321,-2212, -211,-321,-2212};
1196
1197 vector<short> idPosSec[NTrackTypes];
1198 vector<short> idNegSec[NTrackTypes];
1199
1200 vector<short> idPosPrim[NTrackTypes];
1201 vector<short> idNegPrim[NTrackTypes];
1202
1203 for(unsigned short iTr=0; iTr < vRTracks.size(); iTr++)
1204 {
1205 CbmKFTrack &kfTrack = vRTracks[iTr];
1206 bool ok = 1;
1207 for(unsigned short iT=0; iT<5; iT++)
1208 ok = ok && finite(kfTrack.GetTrack()[iT]);
1209 for(unsigned short iC=0; iC<15; iC++)
1210 ok = ok && finite(kfTrack.GetCovMatrix()[iC]);
1211 ok = ok && (kfTrack.GetCovMatrix()[0] < 1. && kfTrack.GetCovMatrix()[0] > 0.)
1212 && (kfTrack.GetCovMatrix()[2] < 1. && kfTrack.GetCovMatrix()[2] > 0.)
1213 && (kfTrack.GetCovMatrix()[5] < 1. && kfTrack.GetCovMatrix()[5] > 0.)
1214 && (kfTrack.GetCovMatrix()[9] < 1. && kfTrack.GetCovMatrix()[9] > 0.)
1215 && (kfTrack.GetCovMatrix()[14] < 1. && kfTrack.GetCovMatrix()[14] > 0.);
1216 ok = ok && kfTrack.GetRefChi2() < 10*kfTrack.GetRefNDF();
1217 if(!ok) continue;
1218
1219 const int pdg = abs(vTrackPDG[iTr]);
1220
1221 short pdgIndex = -1;
1222 switch (pdg)
1223 {
1224 case 11: pdgIndex = 0; break;
1225 case 13: pdgIndex = 1; break;
1226 case 211: pdgIndex = 2; break;
1227 case 321: pdgIndex = 3; break;
1228 case 2212: pdgIndex = 4; break;
1229 }
1230
1231 short incr = 3;
1232 short pdgIndexMax = pdgIndex+incr;
1233
1234 if(pdgIndex<2)
1235 {
1236 incr = 1;
1237 pdgIndexMax = pdgIndex;
1238 }
1239
1240 if(pdgIndex < 0)
1241 {
1242 pdgIndex = 5;
1243 pdgIndexMax = 7;
1244 incr = 1;
1245 }
1246
1247 if( ChiToPrimVtx[iTr] < cuts[0][0] )
1248 {
1249 if(kfTrack.GetTrack()[4] >= 0.)
1250 {
1251 for(int ipdg = pdgIndex; ipdg<=pdgIndexMax; ipdg+=incr )
1252 idPosPrim[ipdg].push_back(iTr);
1253 }
1254 if(kfTrack.GetTrack()[4] < 0.)
1255 for(int ipdg = pdgIndex; ipdg<=pdgIndexMax; ipdg+=incr )
1256 idNegPrim[ipdg].push_back(iTr);
1257 }
1258 else
1259 {
1260 if(kfTrack.GetTrack()[4] >= 0.)
1261 for(int ipdg = pdgIndex; ipdg<=pdgIndexMax; ipdg+=incr )
1262 idPosSec[ipdg].push_back(iTr);
1263 if(kfTrack.GetTrack()[4] < 0.)
1264 for(int ipdg = pdgIndex; ipdg<=pdgIndexMax; ipdg+=incr )
1265 idNegSec[ipdg].push_back(iTr);
1266 }
1267 }
1268
1269 const int nPart = idPosSec[5].size() * idNegSec[5].size()+
1270 idPosSec[5].size() * idNegSec[7].size()+
1271 idPosSec[7].size() * idNegSec[5].size()+
1272 idPosPrim[2].size() * idNegPrim[3].size() +
1273 idPosPrim[3].size() * idNegPrim[2].size() +
1274 idPosPrim[3].size() * idNegPrim[3].size() +
1275 idPosPrim[4].size() * idNegPrim[3].size() +
1276 idPosPrim[3].size() * idNegPrim[4].size() +
1277 idPosPrim[0].size() * idNegPrim[0].size() +
1278 idPosPrim[1].size() * idNegPrim[1].size();
1279
1280//std::cout << "NPart estim " << nPart << std::endl;
1281 Particles.reserve(vRTracks.size() + nPart);
1282
1283 const float massLambdaPDG = 1.115683;
1284 const float massK0sPDG = 0.497614;
1285 const float massKsiPDG = 1.32171;
1286
1287 for(unsigned short iTr=0; iTr < vRTracks.size(); iTr++) {
1288 CbmKFTrack& kfTrack = vRTracks[iTr] ;
1289 kfTrack.SetId(iTr);
1290 CbmKFParticle tmp(&kfTrack);
1291 tmp.SetPDG(211);//vTrackPDG[iTr]);
1292 tmp.SetId(Particles.size());
1293 Particles.push_back(tmp);
1294 }
1295
1296 vector<float> vLambdaTopoChi2Ndf;
1297 vector<CbmKFParticle> vLambdaSec;
1298 vector<CbmKFParticle> vLambdaPrim;
1299
1300 vector<float> vLambdaBarTopoChi2Ndf;
1301 vector<CbmKFParticle> vLambdaBarSec;
1302 vector<CbmKFParticle> vLambdaBarPrim;
1303
1304 vector<float> vK0sTopoChi2Ndf;
1305 vector<CbmKFParticle> vK0sPrim;
1306
1307 vector<CbmKFParticle> vXiPrim;
1308 vector<CbmKFParticle> vXiSec;
1309 vector<CbmKFParticle> vXiBarPrim;
1310
1311 vector<CbmKFParticle> vXiStarPrim;
1312 vector<CbmKFParticle> vXiStarBarPrim;
1313
1314 const float SecCuts[3] = {3.f,5.f,4.f};
1315 //K0s -> pi+ pi-
1316 Find2DaughterDecay(vRTracks, vField, Particles, pdgNeg[5], pdgPos[5], 310,
1317 idNegSec[5], idPosSec[5], PrimVtx, cuts[0], 0, &vK0sTopoChi2Ndf,
1318 SecCuts, massK0sPDG, 0.0022, &vK0sPrim, 0);
1319 //Lambda -> p pi-
1320 Find2DaughterDecay(vRTracks, vField, Particles, pdgNeg[5], pdgPos[7], 3122,
1321 idNegSec[5], idPosSec[7], PrimVtx, cuts[1], 0,&vLambdaTopoChi2Ndf,
1322 SecCuts, massLambdaPDG, 0.0012, &vLambdaPrim, &vLambdaSec);
1323 //Lambda_bar -> pi+ p-
1324 Find2DaughterDecay(vRTracks, vField, Particles, pdgPos[5], pdgNeg[4], -3122,
1325 idPosSec[5], idNegSec[4], PrimVtx, cuts[1], 0, &vLambdaBarTopoChi2Ndf,
1326 SecCuts, massLambdaPDG, 0.0012, &vLambdaBarPrim, &vLambdaBarSec);
1327 //K*0 -> K+ pi-
1328 Find2DaughterDecay(vRTracks, vField, Particles, pdgPos[3], pdgNeg[2], 313,
1329 idPosPrim[3], idNegPrim[2], PrimVtx, cuts[1], 1);
1330 //K*0_bar -> pi+ K-
1331 Find2DaughterDecay(vRTracks, vField, Particles, pdgNeg[3], pdgPos[2], -313,
1332 idNegPrim[3], idPosPrim[2], PrimVtx, cuts[1], 1);
1333 //Lambda* -> p K-
1334 Find2DaughterDecay(vRTracks, vField, Particles, pdgNeg[3], pdgPos[4], 3124,
1335 idNegPrim[3], idPosPrim[4], PrimVtx, cuts[1], 1);
1336 //Lambda*_bar -> K+ p-
1337 Find2DaughterDecay(vRTracks, vField, Particles, pdgNeg[4], pdgPos[3], -3124,
1338 idNegPrim[4], idPosPrim[3], PrimVtx, cuts[1], 1);
1339 //phi -> K+ K-
1340 Find2DaughterDecay(vRTracks, vField, Particles, pdgNeg[3], pdgPos[3], 333,
1341 idNegPrim[3], idPosPrim[3], PrimVtx, cuts[1], 1);
1342// //rho -> pi+ pi-
1343// Find2DaughterDecay(vRTracks, vField, Particles, pdgNeg[2], pdgPos[2], 113,
1344// idNegPrim[2], idPosPrim[2],
1345// PrimVtx, cuts[1]);
1346 //gamma -> e+ e-
1347 Find2DaughterDecay(vRTracks, vField, Particles, pdgNeg[0], pdgPos[0], 22,
1348 idNegPrim[0], idPosPrim[0], PrimVtx, cuts[1], 1);
1349 Find2DaughterDecay(vRTracks, vField, Particles, pdgNeg[0], pdgPos[0], 22,
1350 idNegSec[0], idPosSec[0], PrimVtx, cuts[1], 0);
1351 Find2DaughterDecay(vRTracks, vField, Particles, pdgNeg[0], pdgPos[0], 22,
1352 idNegSec[0], idPosPrim[0], PrimVtx, cuts[1], 0);
1353 Find2DaughterDecay(vRTracks, vField, Particles, pdgNeg[0], pdgPos[0], 22,
1354 idNegPrim[0], idPosSec[0], PrimVtx, cuts[1], 0);
1355 //JPsi-> e+ e-
1356 Find2DaughterDecay(vRTracks, vField, Particles, pdgNeg[0], pdgPos[0], 443,
1357 idNegPrim[0], idPosPrim[0], PrimVtx, cuts[1], 1, 1.f);
1358 //JPsi-> mu+ mu-
1359 Find2DaughterDecay(vRTracks, vField, Particles, pdgNeg[1], pdgPos[1], 100443,
1360 idNegPrim[1], idPosPrim[1], PrimVtx, cuts[1], 1, 1.f);
1361 //rho, omega, phi -> e+ e-
1362 Find2DaughterDecay(vRTracks, vField, Particles, pdgNeg[0], pdgPos[0], 100113,
1363 idNegPrim[0], idPosPrim[0], PrimVtx, cuts[1], 1, 0.2f);
1364 //rho, omega, phi -> mu+ mu-
1365 const float PCut = 1.f;
1366 Find2DaughterDecay(vRTracks, vField, Particles, pdgNeg[1], pdgPos[1], 200113,
1367 idNegPrim[1], idPosPrim[1], PrimVtx, cuts[1], 1, 0.2f, -100, 0, &PCut);
1368
1369 ExtrapolateToPV(vLambdaPrim,PrimVtx);
1370 ExtrapolateToPV(vLambdaBarPrim,PrimVtx);
1371 ExtrapolateToPV(vK0sPrim,PrimVtx);
1372
1373 // Find Xi-
1374 float cutXi[3] = {3.,5.,6.};
1375 FindTrackV0Decay(3312, Particles, vLambdaSec, vRTracks, vField, pdgNeg[5], idNegSec[5],
1376 PrimVtx, cutXi, 0, 0, &vXiPrim, massKsiPDG, 0.002 );
1377
1378
1379 float cutLL[3] = {6.,10000000.,3.};
1380 float cutLL2[3] = {6.,3.,3.};
1381 vector<CbmKFParticle> vLL;
1382 FindTrackV0Decay(3002, vLL, vLambdaSec, vRTracks, vField, pdgNeg[5], idNegSec[5],
1383 PrimVtx, cutLL, 0, &ChiToPrimVtx);
1384 // Find H0->Lambda p pi-
1385 //Find Omega*-
1386 FindTrackV0Decay(3001, Particles, vLL, vRTracks, vField, pdgPos[4], idPosSec[4],
1387 PrimVtx, cutLL2, 0, &ChiToPrimVtx);
1388 // Find Xi+
1389 float cutXiPlus[3] = {3.,5.,6.};
1390 FindTrackV0Decay(-3312, Particles, vLambdaBarSec, vRTracks, vField, pdgPos[5], idPosSec[5],
1391 PrimVtx, cutXiPlus, 0, 0, &vXiBarPrim, massKsiPDG, 0.002);
1392 //Find Omega-
1393 float cutOmega[3] = {3.,3.,3.};
1394 FindTrackV0Decay(3334, Particles, vLambdaSec, vRTracks, vField, pdgNeg[6], idNegSec[6],
1395 PrimVtx, cutOmega, 0, &ChiToPrimVtx);
1396 //Find Omega+
1397 float cutOmegaPlus[3] = {3.,3.,3.};
1398 FindTrackV0Decay(-3334, Particles, vLambdaBarSec, vRTracks, vField, pdgPos[6], idPosSec[6],
1399 PrimVtx, cutOmegaPlus, 0, &ChiToPrimVtx);
1400 //Find Xi*-
1401 float cutXiStarMinus[3] = {-100.,10000.,3.};
1402 FindTrackV0Decay(1003314, Particles, vLambdaPrim, vRTracks, vField, pdgNeg[3], idNegPrim[3],
1403 PrimVtx, cutXiStarMinus, 1);
1404 //Find Xi*+
1405 float cutXiStarPlus[3] = {-100.,10000.,3.};
1406 FindTrackV0Decay(-1003314, Particles, vLambdaBarPrim, vRTracks, vField, pdgPos[3], idPosPrim[3],
1407 PrimVtx, cutXiStarPlus, 1);
1408
1409 ExtrapolateToPV(vXiPrim,PrimVtx);
1410 ExtrapolateToPV(vXiBarPrim,PrimVtx);
1411
1412 //Find Xi*0
1413 float cutXiStar0[3] = {-100.,10000.,3.};
1414 FindTrackV0Decay(3324, Particles, vXiPrim, vRTracks, vField, pdgPos[5], idPosPrim[5],
1415 PrimVtx, cutXiStar0, 1, 0, &vXiStarPrim);
1416 //Find Xi*0 bar
1417 float cutXiBarStar0[3] = {-100.,10000.,3.};
1418 FindTrackV0Decay(-3324, Particles, vXiBarPrim, vRTracks, vField, pdgNeg[5], idNegPrim[5],
1419 PrimVtx, cutXiBarStar0, 1, 0, &vXiStarBarPrim);
1420 //Find Omega*-
1421 const float cutOmegaStar[2] = {-100., 3.};
1422 for(unsigned int iPart=0; iPart<vXiStarPrim.size(); iPart++)
1423 CombineTrackPart(vRTracks, vField, Particles, vXiStarPrim[iPart], pdgNeg[3],
1424 1003334, idNegPrim[3], cutOmegaStar, 0);
1425 //Find Omega*+
1426 for(unsigned int iPart=0; iPart<vXiStarBarPrim.size(); iPart++)
1427 CombineTrackPart(vRTracks, vField, Particles, vXiStarBarPrim[iPart], pdgPos[3],
1428 -1003334, idPosPrim[3], cutOmegaStar, 0);
1429 // Find K*+
1430 float cutKStarPlus[3] = {-100.,10000.,3.};
1431 FindTrackV0Decay(323, Particles, vK0sPrim, vRTracks, vField, pdgPos[5], idPosPrim[5],
1432 PrimVtx, cutKStarPlus, 1, 0);
1433 // Find K*-
1434 float cutKStarMinus[3] = {-100.,10000.,3.};
1435 FindTrackV0Decay(-323, Particles, vK0sPrim, vRTracks, vField, pdgNeg[5], idNegPrim[5],
1436 PrimVtx, cutKStarMinus, 1, 0);
1437 // Find Sigma*+
1438 float cutSigmaStarPlus[3] = {-100.,10000.,3.};
1439 FindTrackV0Decay(3224, Particles, vLambdaPrim, vRTracks, vField, pdgPos[5], idPosPrim[5],
1440 PrimVtx, cutSigmaStarPlus, 1, 0);
1441 // Find Sigma*+ bar
1442 float cutSigmaStarPlusBar[3] = {-100.,10000.,3.};
1443 FindTrackV0Decay(-3114, Particles, vLambdaBarPrim, vRTracks, vField, pdgPos[5], idPosPrim[5],
1444 PrimVtx, cutSigmaStarPlusBar, 1, 0);
1445 // Find Sigma*-
1446 float cutSigmaStarMinus[3] = {-100.,10000.,3.};
1447 FindTrackV0Decay(3114, Particles, vLambdaPrim, vRTracks, vField, pdgNeg[5], idNegPrim[5],
1448 PrimVtx, cutSigmaStarMinus, 1, 0);
1449 // Find Sigma*- bar
1450 float cutSigmaStarMinusBar[3] = {-100.,10000.,3.};
1451 FindTrackV0Decay(-3224, Particles, vLambdaBarPrim, vRTracks, vField, pdgNeg[5], idNegPrim[5],
1452 PrimVtx, cutSigmaStarMinusBar, 1, 0);
1453 // Find H-dibarion
1454 vector<CbmKFParticle> vHdibarion;
1455 float cutHdb[3] = {3.,3.,3.};
1456 for(unsigned short iL=0; iL < vLambdaSec.size(); iL++)
1457 {
1458 CbmKFParticle_simd vDaughters[2] = {CbmKFParticle_simd(),CbmKFParticle_simd(vLambdaSec[iL])};
1459
1460 vector<int> daughterIds;
1461 for(unsigned int iD=0; iD<vLambdaSec[iL].DaughterIds().size(); iD++)
1462 daughterIds.push_back(vLambdaSec[iL].DaughterIds()[iD]);
1463 FindHyperons(3000, vDaughters, daughterIds, vLambdaSec, vHdibarion, PrimVtx, cutHdb, iL+1);
1464 }
1465
1466 for(unsigned int iH=0; iH<vHdibarion.size(); iH++)
1467 {
1468 vHdibarion[iH].SetId(Particles.size());
1469 Particles.push_back(vHdibarion[iH]);
1470 }
1471 // chi2_prim chi2_geo z pt chi2_topo
1472 const float cutsD[8][8] = {{ 6., 3., 0.04, 0.3, 3.}, //D0 -> pi+ K-
1473 { 6., 3., 0.04, 0.3, 3.}, //D+ -> K- pi+ pi+
1474 { 6., 3., 0.04, 0.3, 3.}, //D0 -> pi+ pi+ pi- K-
1475 { 6., 3., 0.04, 0.3, 3.}, //Ds+ -> K- K+ pi+
1476 { 6., 3., 0.04, 0.3, 3.}, //Lambdac -> pi+ K- p
1477 { 6., 3., -100., 0.3, -100.}, //D*0 -> D+ pi-
1478 { 6., 3., -100., 0.3, -100.}, //D*+ -> D0 pi+
1479 { 6., 3., -100., 0.3, -100.}}; //D*+4 -> D04 pi+
1480
1481 const int DMesLambdcDaughterPDG[5] = { 211, -321, -211, 321, 2212 };
1482 const int DMesLambdcMotherPDG[8] = { 421, 411, 100421, 431, 4122, 10421, 10411, 20411 };
1483 vector<short>* DMesLambdcIdTrack[5] = {&idPosSec[5],
1484 &idNegSec[3],
1485 &idNegSec[5],
1486 &idPosSec[3],
1487 &idPosSec[4]};
1488 FindDMesLambdac(vRTracks, vField, Particles, DMesLambdcDaughterPDG, DMesLambdcMotherPDG,
1489 DMesLambdcIdTrack, PrimVtx, cutsD, ChiToPrimVtx);
1490
1491 const int DMesLambdcBarDaughterPDG[5] = { -211, 321, 211, -321, -2212 };
1492 const int DMesLambdcBarMotherPDG[8] = { -421, -411, -100421, -431, -4122, -10421, -10411, -20411 };
1493 vector<short>* DMesLambdcBarIdTrack[5] = {&idNegSec[5],
1494 &idPosSec[3],
1495 &idPosSec[5],
1496 &idNegSec[3],
1497 &idNegSec[4]};
1498 FindDMesLambdac(vRTracks, vField, Particles, DMesLambdcBarDaughterPDG, DMesLambdcBarMotherPDG,
1499 DMesLambdcBarIdTrack, PrimVtx, cutsD, ChiToPrimVtx);
1500// std::cout << "NPart " << Particles.size() << std::endl;
1501}
1502
1503// void CbmKFParticleInterface::FindParticles(vector<CbmL1Track>& vRTracks, vector<CbmKFParticle>& Particles,
1504// CbmKFVertex& PrimVtx, const vector<int>& vTrackPDG, const float cuts[2][3])
1505// {
1506// FindParticlesT(vRTracks, Particles, PrimVtx, vTrackPDG, cuts);
1507// }
1508//
1509// void CbmKFParticleInterface::FindParticles(vector<CbmStsTrack>& vRTracks, vector<CbmKFParticle>& Particles,
1510// CbmKFVertex& PrimVtx, const vector<int>& vTrackPDG, const float cuts[2][3])
1511// {
1512// FindParticlesT(vRTracks, Particles, PrimVtx, vTrackPDG, cuts);
1513// }
1514
1515void CbmKFParticleInterface::ExtrapolateToPV(vector<CbmKFParticle>& vParticles, CbmKFVertex& PrimVtx)
1516{
1517 fvec pVtx[3] = {PrimVtx.GetRefX(),PrimVtx.GetRefY(),PrimVtx.GetRefZ()};
1518 for(unsigned int iL=0; iL<vParticles.size(); iL += fvecLen)
1519 {
1520 CbmKFParticle* parts[fvecLen];
1521 unsigned int nPart = vParticles.size();
1522
1523 unsigned int nEntries = (iL + fvecLen < nPart) ? fvecLen : (nPart - iL);
1524
1525 for(unsigned short iv=0; iv<nEntries; iv++)
1526 parts[iv] = &vParticles[iL+iv];
1527
1528
1529 CbmKFParticleInterface tmp(parts,nEntries);
1530 tmp.Extrapolate( tmp.KFPart.r, tmp.GetDStoPoint(pVtx) );
1531
1532 for(unsigned int iv=0; iv<nEntries; iv++)
1533 {
1534 tmp.GetKFParticle(vParticles[iL+iv], iv);
1535 }
1536 }
1537}
1538
1539fvec CbmKFParticleInterface::GetChi2BetweenParticles(CbmKFParticle_simd &p1, CbmKFParticle_simd &p2)
1540{
1541 const fvec& x1 = p1.GetX();
1542 const fvec& y1 = p1.GetY();
1543 const fvec& z1 = p1.GetZ();
1544
1545 const fvec& x2 = p2.GetX();
1546 const fvec& y2 = p2.GetY();
1547 const fvec& z2 = p2.GetZ();
1548
1549 const fvec dx = x1 - x2;
1550 const fvec dy = y1 - y2;
1551 const fvec dz = z1 - z2;
1552
1553 const fvec& c0 = p1.GetCovariance(0) + p2.GetCovariance(0);
1554 const fvec& c1 = p1.GetCovariance(1) + p2.GetCovariance(1);
1555 const fvec& c2 = p1.GetCovariance(2) + p2.GetCovariance(2);
1556 const fvec& c3 = p1.GetCovariance(3) + p2.GetCovariance(3);
1557 const fvec& c4 = p1.GetCovariance(4) + p2.GetCovariance(4);
1558 const fvec& c5 = p1.GetCovariance(5) + p2.GetCovariance(5);
1559
1560 const fvec r2 = dx*dx + dy*dy + dz*dz;
1561 const fvec err2 = c0*dx*dx + c2*dy*dy + c5*dz*dz + 2*( c1*dx*dy + c3*dx*dz + c4*dy*dz );
1562
1563 return (r2*r2/err2);
1564}
1565
1566//template<class T>
1567void CbmKFParticleInterface::Find2DaughterDecay(vector<CbmKFTrack>& vTracks,
1568 const vector<L1FieldRegion>& vField,
1569 vector<CbmKFParticle>& Particles,
1570 const int DaughterNegPDG,
1571 const int DaughterPosPDG,
1572 const int MotherPDG,
1573 vector<short>& idNeg,
1574 vector<short>& idPos,
1575 CbmKFVertex& PrimVtx,
1576 const float* cuts,
1577 bool isPrimary,
1578 vector<float>* vMotherTopoChi2Ndf,
1579 const float* secCuts,
1580 const float massMotherPDG,
1581 const float massMotherPDGSigma,
1582 vector<CbmKFParticle>* vMotherPrim,
1583 vector<CbmKFParticle>* vMotherSec )
1584{
1585 const int NPositive = idPos.size();
1586 CbmKFParticle mother_temp;
1587
1588 CbmKFParticleInterface mother;
1589 mother.SetPDG( MotherPDG );
1590 CbmKFParticleInterface motherTopo;
1591
1592 const short NPosVect = NPositive%fvecLen == 0 ? NPositive/fvecLen : NPositive/fvecLen+1;
1593 vector<CbmKFParticle_simd> posPart(NPosVect);
1594
1595 // create particles from tracks
1596 CbmKFTrackInterface* vvPos[fvecLen];
1597
1598 int iPosVect = 0;
1599 for(unsigned short iTrP=0; iTrP < NPositive; iTrP += fvecLen)
1600 {
1601 const unsigned short NTracks = (iTrP + fvecLen < NPositive) ? fvecLen : (NPositive - iTrP);
1602
1603 L1FieldRegion posField;
1604 for(unsigned short iv=0; iv<NTracks; iv++)
1605 {
1606 vvPos[iv] = &vTracks[idPos[iTrP+iv]];
1607
1608 int entrSIMDPos = idPos[iTrP+iv] % fvecLen;
1609 int entrVecPos = idPos[iTrP+iv] / fvecLen;
1610 posField.SetOneEntry(iv,vField[entrVecPos],entrSIMDPos);
1611 }
1612
1613 posPart[iPosVect].Create(vvPos,NTracks,0,&DaughterPosPDG);
1614 posPart[iPosVect].SetField(posField);
1615 fvec posId;
1616 for(int iv=0; iv<NTracks; iv++)
1617 posId[iv] = idPos[iTrP+iv];
1618 posPart[iPosVect].SetId(posId);
1619 iPosVect++;
1620 }
1621
1622 for(unsigned short iTrN=0; iTrN < idNeg.size(); iTrN++)
1623 {
1624 CbmKFTrack &kfTrackNeg = vTracks[idNeg[iTrN]];
1625 CbmKFParticle_simd vDaughters[2] = {CbmKFParticle_simd(kfTrackNeg,0,&DaughterNegPDG),
1626 CbmKFParticle_simd()};
1627 int entrSIMD = idNeg[iTrN] % fvecLen;
1628 int entrVec = idNeg[iTrN] / fvecLen;
1629 vDaughters[0].SetField(vField[entrVec],1,entrSIMD);
1630 vDaughters[0].SetId(idNeg[iTrN]);
1631
1632
1633 for(unsigned short iTrP=0; iTrP < NPositive; iTrP += fvecLen)
1634 {
1635 const unsigned short NTracks = (iTrP + fvecLen < NPositive) ? fvecLen : (NPositive - iTrP);
1636
1637 vDaughters[1] = posPart[iTrP/fvecLen];
1638
1639 if(isPrimary)
1640 {
1641 fvec vtxGuess[3] = {PrimVtx.GetRefX(),
1642 PrimVtx.GetRefY(),
1643 PrimVtx.GetRefZ()};
1644 fvec errGuess[3] = {100*sqrt(PrimVtx.GetCovMatrix()[0]),
1645 100*sqrt(PrimVtx.GetCovMatrix()[2]),
1646 100*sqrt(PrimVtx.GetCovMatrix()[5])};
1647 mother.SetVtxGuess(vtxGuess[0], vtxGuess[1], vtxGuess[2]);
1648 mother.SetVtxErrGuess(errGuess[0], errGuess[1], errGuess[2]);
1649 mother.Construct(vDaughters, 2, 0, -1, -1, 1);
1650 }
1651 else
1652 mother.Construct(vDaughters, 2, 0);
1653
1654 if(vMotherTopoChi2Ndf)
1655 {
1656 motherTopo = mother;
1657 motherTopo.MeasureProductionVertex(&motherTopo.KFPart, motherTopo.GetParameters(), &PrimVtx);
1658 }
1659
1660// if( (fabs(DaughterNegPDG) == 211) && (fabs(DaughterPosPDG) == 211) )
1661// {
1662// fvec pl0[3] = {vDaughters[0].GetPx(),vDaughters[0].GetPy(),vDaughters[0].GetPz()};
1663// fvec pl1[3] = {vDaughters[1].GetPx(),vDaughters[1].GetPy(),vDaughters[1].GetPz()};
1664//
1665// fvec pcm[3] = {pl0[0]+pl1[0], pl0[1]+pl1[1], pl0[2]+pl1[2]};
1666// fvec ecm = vDaughters[0].GetE()+vDaughters[1].GetE();
1667//
1668// fvec gamma = sqrt(fvec(1.) - (pcm[0]*pcm[0]+pcm[1]*pcm[1]+pcm[2]*pcm[2])/(ecm*ecm));
1669//
1670// fvec e0 = gamma*(vDaughters[0].GetE() - (pcm[0]*pl0[0]+pcm[1]*pl0[1]+pcm[2]*pl0[2])/ecm);
1671// fvec e1 = gamma*(vDaughters[1].GetE() - (pcm[0]*pl1[0]+pcm[1]*pl1[1]+pcm[2]*pl1[2])/ecm);
1672// if( fabs(e0[0]-e1[0]) > 1.e-7 )
1673// std::cout << e0 << " " << e1 << std::endl;
1674// }
1675
1676 // CbmKFParticle mother_temp[fvecLen];
1677 // mother.GetKFParticle(mother_temp, NTracks);
1678
1679 for(int iv=0; iv<NTracks; iv++)
1680 {
1681 if(!finite(mother.GetChi2()[iv])) continue;
1682 if(!(mother.GetChi2()[iv] > 0.0f)) continue;
1683 if(!(mother.GetChi2()[iv]==mother.GetChi2()[iv])) continue;
1684 if(mother.GetZ()[iv] < cuts[2]) continue;
1685
1686 if( mother.GetChi2()[iv]/mother.GetNDF()[iv] < cuts[1] )
1687 {
1688 mother.GetKFParticle(mother_temp, iv);
1689 mother_temp.SetId(Particles.size());
1690 Particles.push_back(mother_temp);
1691
1692 float motherTopoChi2Ndf(0);
1693 if(vMotherTopoChi2Ndf)
1694 motherTopoChi2Ndf = motherTopo.GetChi2()[iv]/motherTopo.GetNDF()[iv];
1695// vMotherTopoChi2Ndf->push_back(motherTopo.GetChi2()[iv]/motherTopo.GetNDF()[iv]);
1696
1697 if(vMotherPrim)
1698 {
1699 Double_t mass, errMass;
1700
1701 mother_temp.GetMass(mass, errMass);
1702 if( (fabs(mass - massMotherPDG)/massMotherPDGSigma) > secCuts[0] ) continue;
1703
1704 if( motherTopoChi2Ndf < secCuts[1] )
1705 {
1706 vMotherPrim->push_back(mother_temp);
1707 }
1708 else if(vMotherSec)
1709 {
1710 if( mother_temp.GetZ() < secCuts[2] ) continue;
1711 vMotherSec->push_back(mother_temp);
1712 }
1713 }
1714 }
1715 }
1716 }
1717 }
1718}
1719
1720//template<class T>
1721void CbmKFParticleInterface::Find2DaughterDecay(vector<CbmKFTrack>& vTracks,
1722 const vector<L1FieldRegion>& vField,
1723 vector<CbmKFParticle>& Particles,
1724 const int DaughterNegPDG,
1725 const int DaughterPosPDG,
1726 const int MotherPDG,
1727 vector<short>& idNeg,
1728 vector<short>& idPos,
1729 CbmKFVertex& PrimVtx,
1730 const float* cuts,
1731 bool isPrimary,
1732 const float PtCut,
1733 const float Chi2PrimCut,
1734 vector<float>* ChiToPrimVtx,
1735 const float* PCut)
1736{
1737 vector<short> idPosPt;
1738 vector<short> idNegPt;
1739
1740 for(unsigned int iEl=0; iEl<idPos.size(); iEl++)
1741 {
1742 CbmKFTrack& kfTrack = vTracks[idPos[iEl]];
1743 const float tx = kfTrack.GetTrack()[2];
1744 const float ty = kfTrack.GetTrack()[3];
1745 const float p = fabs(1/kfTrack.GetTrack()[4]);
1746 const float t2 = tx*tx+ty*ty;
1747 float pt = p*(sqrt(t2/(1+t2)));
1748 if(PCut)
1749 if(p < *PCut) continue;
1750 if(pt<PtCut) continue;
1751 if(Chi2PrimCut > -1)
1752 if(ChiToPrimVtx->at(idPos[iEl]) < Chi2PrimCut) continue;
1753 idPosPt.push_back(idPos[iEl]);
1754 }
1755
1756 for(unsigned int iEl=0; iEl<idNeg.size(); iEl++)
1757 {
1758 CbmKFTrack& kfTrack = vTracks[idNeg[iEl]];
1759 const float tx = kfTrack.GetTrack()[2];
1760 const float ty = kfTrack.GetTrack()[3];
1761 const float p = fabs(1/kfTrack.GetTrack()[4]);
1762 const float t2 = tx*tx+ty*ty;
1763 float pt = p*(sqrt(t2/(1+t2)));
1764 if(pt<PtCut) continue;
1765 if(Chi2PrimCut > -1)
1766 if(ChiToPrimVtx->at(idNeg[iEl]) < Chi2PrimCut) continue;
1767 idNegPt.push_back(idNeg[iEl]);
1768 }
1769
1770 Find2DaughterDecay(vTracks, vField,
1771 Particles, DaughterNegPDG, DaughterPosPDG, MotherPDG,
1772 idNegPt, idPosPt,
1773 PrimVtx, cuts, isPrimary);
1774}
1775
1776//template<class T>
1777void CbmKFParticleInterface::FindTrackV0Decay(const int MotherPDG,
1778 vector<CbmKFParticle>& Particles,
1779 vector<CbmKFParticle>& vV0,
1780 vector<CbmKFTrack>& vTracks,
1781 const vector<L1FieldRegion>& vField,
1782 const int DaughterPDG,
1783 vector<short>& idTrack,
1784 CbmKFVertex& PrimVtx,
1785 const float* cuts,
1786 bool isPrimary,
1787 vector<float>* ChiToPrimVtx,
1788 vector<CbmKFParticle>* vHyperonPrim,
1789 float hyperonPrimMass,
1790 float hyperonPrimMassErr,
1791 vector<CbmKFParticle>* vHyperonSec)
1792{
1793 CbmKFParticle hyperon_temp;
1794 CbmKFParticleInterface hyperon;
1795 hyperon.SetPDG( MotherPDG );
1796
1797 for(unsigned short iV0=0; iV0 < vV0.size(); iV0++)
1798 {
1799 unsigned short nElements = 0;
1800 CbmKFParticle_simd vDaughters[2]= {CbmKFParticle_simd(vV0[iV0]),CbmKFParticle_simd()};
1801
1802 CbmKFTrackInterface* vvTr[fvecLen];
1803 L1FieldRegion field;
1804 fvec trId;
1805
1806 for(unsigned short iTr=0; iTr < idTrack.size(); iTr++)
1807 {
1808 bool ok = 1;
1809 if(ChiToPrimVtx)
1810 if( (ChiToPrimVtx->at(idTrack[iTr]) < 7) ) ok=0; //TODO 7 for Omega
1811
1812 if(ok)
1813 {
1814 trId[nElements] = idTrack[iTr];
1815 vvTr[nElements] = &vTracks[idTrack[iTr]];
1816
1817 int entrSIMD = idTrack[iTr] % fvecLen;
1818 int entrVec = idTrack[iTr] / fvecLen;
1819 field.SetOneEntry(nElements,vField[entrVec],entrSIMD);
1820
1821 nElements++;
1822 }
1823 else if( (iTr != idTrack.size()-1) ) continue;
1824
1825 if( (nElements == fvecLen) || ((iTr == idTrack.size()-1)&&(nElements>0)) )
1826 {
1827 vDaughters[1].Create(vvTr,nElements,0,&DaughterPDG);
1828 vDaughters[1].SetField(field);
1829 vDaughters[1].SetId(trId);
1830
1831 if(isPrimary)
1832 {
1833 fvec vtxGuess[3] = {PrimVtx.GetRefX(),
1834 PrimVtx.GetRefY(),
1835 PrimVtx.GetRefZ()};
1836 fvec errGuess[3] = {100*sqrt(PrimVtx.GetCovMatrix()[0]),
1837 100*sqrt(PrimVtx.GetCovMatrix()[2]),
1838 100*sqrt(PrimVtx.GetCovMatrix()[5])};
1839 hyperon.SetVtxGuess(vtxGuess[0], vtxGuess[1], vtxGuess[2]);
1840 hyperon.SetVtxErrGuess(errGuess[0], errGuess[1], errGuess[2]);
1841 hyperon.Construct(vDaughters, 2, 0, -1, -1, 1);
1842 }
1843 else
1844 hyperon.Construct(vDaughters, 2, 0);
1845
1846 CbmKFParticleInterface hyperonTopo(hyperon);
1847 hyperonTopo.MeasureProductionVertex(&hyperonTopo.KFPart, hyperonTopo.GetParameters(), &PrimVtx);
1848
1849 for(unsigned int iv=0; iv<nElements; iv++)
1850 {
1851 bool isSameTrack = 0;
1852 for(unsigned short iD=0; iD<vV0[iV0].DaughterIds().size(); iD++)
1853 if(vV0[iV0].DaughterIds()[iD] == trId[iv]) isSameTrack=1;
1854
1855 if(isSameTrack) continue;
1856 if(!finite(hyperon.GetChi2()[iv])) continue;
1857 if(!(hyperon.GetChi2()[iv] > 0.0f)) continue;
1858 if(!(hyperon.GetChi2()[iv]==hyperon.GetChi2()[iv])) continue;
1859 if((hyperon.GetZ()[iv] < cuts[0]) || ((vV0[iV0].GetZ() - hyperon.GetZ()[iv]) < 0)) continue;
1860
1861 if(hyperonTopo.GetChi2()[iv]/hyperonTopo.GetNDF()[iv] > cuts[1] ) continue;
1862
1863 if( hyperon.GetChi2()[iv]/hyperon.GetNDF()[iv] > cuts[2] ) continue;
1864 hyperon.GetKFParticle(hyperon_temp, iv);
1865
1866 hyperon_temp.SetId(Particles.size());
1867 Particles.push_back(hyperon_temp);
1868
1869 if(vHyperonPrim)
1870 {
1871 Double_t mass, errMass;
1872
1873 hyperon_temp.GetMass(mass, errMass);
1874 if( (fabs(mass - hyperonPrimMass)/hyperonPrimMassErr) <= 3 )
1875 vHyperonPrim->push_back(hyperon_temp);
1876 else
1877 if( vHyperonSec )
1878 vHyperonSec->push_back(hyperon_temp);
1879 }
1880 }
1881 nElements=0;
1882 }
1883 }
1884 }
1885}
1886
1887void CbmKFParticleInterface::FindHyperons(int PDG,
1888 CbmKFParticle_simd vDaughters[2],
1889 vector<int>& daughterIds,
1890 vector<CbmKFParticle>& vLambdaSec,
1891 vector<CbmKFParticle>& vHyperon,
1892 CbmKFVertex& PrimVtx,
1893 const float *cuts,
1894 int startIndex)
1895{
1896 CbmKFParticle* lambdas[fvecLen];
1897 int nLambdasSec = vLambdaSec.size();
1898
1899 for(unsigned short iL=startIndex; iL < vLambdaSec.size(); iL += fvecLen)
1900 {
1901 unsigned int nEntries = (iL + fvecLen < nLambdasSec) ? fvecLen : (nLambdasSec - iL);
1902
1903 for(unsigned short iv=0; iv<nEntries; iv++)
1904 lambdas[iv] = &vLambdaSec[iL+iv];
1905
1906 vDaughters[0].Create(lambdas,nEntries);
1907
1908 CbmKFParticleInterface hyperon;
1909
1910 hyperon.SetPDG( PDG );
1911 hyperon.Construct(vDaughters, 2, 0);
1912
1913 CbmKFParticleInterface hyperonTopo(hyperon);
1914 hyperonTopo.MeasureProductionVertex(&hyperonTopo.KFPart, hyperonTopo.GetParameters(), &PrimVtx);
1915
1916 for(unsigned int iv=0; iv<nEntries; iv++)
1917 {
1918 bool isSameTrack = 0;
1919 for(unsigned short iD=0; iD<lambdas[iv]->DaughterIds().size(); iD++)
1920 for(unsigned short iD0=0; iD0<daughterIds.size(); iD0++)
1921 if(lambdas[iv]->DaughterIds()[iD] == daughterIds[iD0]) isSameTrack=1;
1922
1923 if(isSameTrack) continue;
1924 if(!finite(hyperon.GetChi2()[iv])) continue;
1925 if(!(hyperon.GetChi2()[iv] > 0.0f)) continue;
1926 if(!(hyperon.GetChi2()[iv]==hyperon.GetChi2()[iv])) continue;
1927 if((hyperon.GetZ()[iv] < cuts[0]) || ((lambdas[iv]->GetZ() - hyperon.GetZ()[iv]) < 0)) continue;
1928
1929 if(hyperonTopo.GetChi2()[iv]/hyperonTopo.GetNDF()[iv] > cuts[1] ) continue;
1930
1931 if( hyperon.GetChi2()[iv]/hyperon.GetNDF()[iv] > cuts[2] ) continue;
1932 CbmKFParticle hyperon_temp;
1933 hyperon.GetKFParticle(hyperon_temp, iv);
1934 vHyperon.push_back(hyperon_temp);
1935 }
1936 }
1937}
1938
1939//template<class T>
1940void CbmKFParticleInterface::FindDMesLambdac(vector<CbmKFTrack>& vTracks,
1941 const vector<L1FieldRegion>& vField,
1942 vector<CbmKFParticle>& Particles,
1943 const int DaughterPDG[5], //pi, K_b, pi_b, K, p
1944 const int MotherPDG[8], //D0, D+, D0_4, Ds+, Lc
1945 vector<short>* idTrack[5], //pi, K_b, pi_b, K, p
1946 CbmKFVertex& PrimVtx,
1947 const float cuts[8][8],
1948 vector<float> ChiToPrimVtx)
1949{
1950 vector<short> id[5]; //pi, K_b, pi_b, K, p
1951 for(int iId=0; iId<5; iId++)
1952 {
1953 for(unsigned int iTr=0; iTr<idTrack[iId]->size(); iTr++)
1954 {
1955 CbmKFTrack& kfTrack = vTracks[idTrack[iId]->at(iTr)];
1956 const float tx = kfTrack.GetTrack()[2];
1957 const float ty = kfTrack.GetTrack()[3];
1958 const float p = fabs(1/kfTrack.GetTrack()[4]);
1959 const float t2 = tx*tx+ty*ty;
1960 float pt = p*(sqrt(t2/(1+t2)));
1961 if(pt<cuts[0][3]) continue;
1962 if(ChiToPrimVtx[idTrack[iId]->at(iTr)] < cuts[0][0]) continue;
1963 id[iId].push_back(idTrack[iId]->at(iTr));
1964 }
1965 }
1966
1967 vector<CbmKFParticle> kpi;
1968 vector<CbmKFParticle> kpipi;
1969 vector<CbmKFParticle> kpipipi;
1970 vector<CbmKFParticle> kpik;
1971 vector<CbmKFParticle> kpip;
1972
1973 const float cutskpi[3] = {3., 3., -100.};
1974 Find2DaughterDecay(vTracks, vField, kpi, DaughterPDG[0], DaughterPDG[1], MotherPDG[0],
1975 id[0], id[1], PrimVtx, cutskpi, 0);
1976
1977 for(unsigned int iKPi=0; iKPi<kpi.size(); iKPi++)
1978 {
1979 unsigned short startPiIndex = kpi[iKPi].DaughterIds()[0]+1;
1980
1981 CombineTrackPart(vTracks,vField,kpipi,kpi[iKPi],DaughterPDG[0], MotherPDG[1], id[0], cuts[1], startPiIndex, 1);
1982 CombineTrackPart(vTracks,vField,kpik ,kpi[iKPi],DaughterPDG[3], MotherPDG[3], id[3], cuts[3], 0, 1);
1983 CombineTrackPart(vTracks,vField,kpip ,kpi[iKPi],DaughterPDG[4], MotherPDG[4], id[4], cuts[4], 0, 1);
1984 }
1985 for(unsigned int iKPiPi=0; iKPiPi<kpipi.size(); iKPiPi++)
1986 CombineTrackPart(vTracks,vField,kpipipi,kpipi[iKPiPi],DaughterPDG[2], MotherPDG[2], id[2], cuts[2], 0, 1);
1987
1988 SelectParticleCandidates(Particles, kpi, PrimVtx, cuts[0]);
1989 SelectParticleCandidates(Particles, kpipi, PrimVtx, cuts[1]);
1990 SelectParticleCandidates(Particles, kpipipi, PrimVtx, cuts[2]);
1991 SelectParticleCandidates(Particles, kpik, PrimVtx, cuts[3]);
1992 SelectParticleCandidates(Particles, kpip, PrimVtx, cuts[4]);
1993
1994 vector<CbmKFParticle> d0pi;
1995 vector<CbmKFParticle> d2pi;
1996 vector<CbmKFParticle> d4pi;
1997
1998 for(unsigned int iKPiPi=0; iKPiPi<kpipi.size(); iKPiPi++)
1999 CombineTrackPart(vTracks,vField,d0pi,kpipi[iKPiPi],DaughterPDG[2], MotherPDG[5], id[2], cuts[5], 0, 0);
2000
2001 for(unsigned int iKPi=0; iKPi<kpi.size(); iKPi++)
2002 CombineTrackPart(vTracks,vField,d2pi,kpi[iKPi],DaughterPDG[0], MotherPDG[6], id[0], cuts[6], 0, 0);
2003
2004 for(unsigned int iKPiPiPi=0; iKPiPiPi<kpipipi.size(); iKPiPiPi++)
2005 CombineTrackPart(vTracks,vField,d4pi,kpipipi[iKPiPiPi],DaughterPDG[0], MotherPDG[7], id[0], cuts[7], 0, 0);
2006
2007 SelectParticleCandidates(Particles, d0pi, PrimVtx, cuts[5]);
2008 SelectParticleCandidates(Particles, d2pi, PrimVtx, cuts[6]);
2009 SelectParticleCandidates(Particles, d4pi, PrimVtx, cuts[7]);
2010}
2011
2012//template<class T>
2013void CbmKFParticleInterface::CombineTrackPart(vector<CbmKFTrack>& vTracks,
2014 const vector<L1FieldRegion>& vField,
2015 vector<CbmKFParticle>& Particles,
2016 CbmKFParticle& part,
2017 const int DaughterPDG,
2018 const int MotherPDG,
2019 vector<short>& id,
2020 const float* cuts,
2021 const unsigned short startIndex,
2022 const bool IsSamePart)
2023{
2024 CbmKFParticle_simd vDaughters[2] = {CbmKFParticle_simd(part),CbmKFParticle_simd()};
2025
2026 CbmKFTrackInterface* vTr[fvecLen];
2027
2028 fvec vtxGuess[3] = {part.GetX(), part.GetY(), part.GetZ()};
2029 fvec errGuess[3] = {fvec(10.f*sqrt(part.GetCovMatrix()[0])),
2030 fvec(10.f*sqrt(part.GetCovMatrix()[2])),
2031 fvec(10.f*sqrt(part.GetCovMatrix()[5]))};
2032
2033 const unsigned short NTr = id.size();
2034 for(unsigned short iTr=startIndex; iTr < NTr; iTr += fvecLen)
2035 {
2036 const unsigned short NTracks = (iTr + fvecLen < NTr) ? fvecLen : (NTr - iTr);
2037
2038 L1FieldRegion trField;
2039 for(unsigned short iv=0; iv<NTracks; iv++)
2040 {
2041 vTr[iv] = &vTracks[id[iTr+iv]];
2042
2043 int entrSIMDPos = id[iTr+iv] % fvecLen;
2044 int entrVecPos = id[iTr+iv] / fvecLen;
2045 trField.SetOneEntry(iv,vField[entrVecPos],entrSIMDPos);
2046 }
2047
2048 vDaughters[1].Create(vTr,NTracks,0,&DaughterPDG);
2049 vDaughters[1].SetField(trField);
2050 fvec trId;
2051 for(int iv=0; iv<NTracks; iv++)
2052 trId[iv] = id[iTr+iv];
2053 vDaughters[1].SetId(trId);
2054
2055 CbmKFParticleInterface mother;
2056 mother.SetPDG( MotherPDG );
2057 if(IsSamePart)
2058 {
2059 mother.SetVtxGuess(vtxGuess[0], vtxGuess[1], vtxGuess[2]);
2060 mother.SetVtxErrGuess(errGuess[0], errGuess[1], errGuess[2]);
2061 }
2062 mother.Construct(vDaughters, 2, 0);
2063
2064 for(int iv=0; iv<NTracks; iv++)
2065 {
2066 if(!finite(mother.GetChi2()[iv])) continue;
2067 if(!(mother.GetChi2()[iv] > 0.0f)) continue;
2068 if(!(mother.GetChi2()[iv]==mother.GetChi2()[iv])) continue;
2069
2070 if( mother.GetChi2()[iv]/mother.GetNDF()[iv] > cuts[1] ) continue;
2071
2072 bool isSameTrack = 0;
2073 for(unsigned short iD=0; iD<part.DaughterIds().size(); iD++)
2074 if(part.DaughterIds()[iD] == id[iTr+iv]) isSameTrack=1;
2075 if(isSameTrack) continue;
2076
2077 CbmKFParticle mother_temp;
2078 mother.GetKFParticle(mother_temp, iv);
2079 mother_temp.CleanDaughtersId();
2080 mother_temp.AddDaughter(id[iTr+iv]);
2081 for(unsigned short iD=0; iD<part.DaughterIds().size(); iD++)
2082 mother_temp.AddDaughter(part.DaughterIds()[iD]);
2083 Particles.push_back(mother_temp);
2084 }
2085 }
2086}
2087
2088void CbmKFParticleInterface::SelectParticleCandidates(vector<CbmKFParticle>& Particles,
2089 vector<CbmKFParticle>& vCandidates,
2090 CbmKFVertex& PrimVtx,
2091 const float cuts[5])
2092{
2093 CbmKFParticle* cand[fvecLen];
2094 int nCand = vCandidates.size();
2095
2096 for(unsigned short iC=0; iC < nCand; iC += fvecLen)
2097 {
2098 unsigned int nEntries = (iC + fvecLen < nCand) ? fvecLen : (nCand - iC);
2099
2100 for(unsigned short iv=0; iv<nEntries; iv++)
2101 cand[iv] = &vCandidates[iC+iv];
2102
2103 CbmKFParticleInterface candTopo(cand,nEntries);
2104
2105 candTopo.MeasureProductionVertex(&candTopo.KFPart, candTopo.GetParameters(), &PrimVtx);
2106
2107 for(unsigned int iv=0; iv<nEntries; iv++)
2108 {
2109 if(!finite(candTopo.GetChi2()[iv])) continue;
2110 if(!(candTopo.GetChi2()[iv] > 0.0f)) continue;
2111 if(!(candTopo.GetChi2()[iv]==candTopo.GetChi2()[iv])) continue;
2112 if(candTopo.GetZ()[iv] < cuts[2]) continue;
2113
2114 if(candTopo.GetChi2()[iv]/candTopo.GetNDF()[iv] > cuts[4] ) continue;
2115
2116 Particles.push_back(vCandidates[iC+iv]);
2117 }
2118 }
2119}
2120
2121
2122
2123//template<class T>
2124void CbmKFParticleInterface::ConstructPVT(vector<CbmKFTrack>& vRTracks)
2125{
2126
2127}
2128
2129// void CbmKFParticleInterface::ConstructPV(vector<CbmL1Track>& vRTracks)
2130// {
2131// ConstructPVT(vRTracks);
2132// }
2133//
2134// void CbmKFParticleInterface::ConstructPV(vector<CbmStsTrack>& vRTracks)
2135// {
2136// ConstructPVT(vRTracks);
2137// }
2138
2139#undef cnst
d
const Float_t d
Z-ccordinate of the first GEM-station.
Definition BmnMwpcHit.cxx:7
cnst
#define cnst
Definition L1AddMaterial.h:9
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition P4_F32vec4.h:34
fabs
friend F32vec4 fabs(const F32vec4 &a)
Definition P4_F32vec4.h:52
fvecLen
const int fvecLen
Definition P4_F32vec4.h:233
i
int i
Definition P4_F32vec4.h:22
m
__m128 m
Definition P4_F32vec4.h:27
rcp
friend F32vec4 rcp(const F32vec4 &a)
Definition P4_F32vec4.h:45
fvec
F32vec4 fvec
Definition P4_F32vec4.h:231
CbmKFParticleInterface::multQSQt
void multQSQt(const fvec Q[], fvec S[])
Definition CbmKFParticleInterface.cxx:885
CbmKFParticleInterface::MeasureMass
void MeasureMass(CbmKFParticle_simd *Particle, fvec r0[], fvec Mass)
Definition CbmKFParticleInterface.cxx:435
CbmKFParticleInterface::GetKFVertex
void GetKFVertex(CbmKFVertex *vtx)
Definition CbmKFParticleInterface.cxx:1081
CbmKFParticleInterface::Find2DaughterDecay
static void Find2DaughterDecay(vector< CbmKFTrack > &vTracks, const vector< L1FieldRegion > &vField, vector< CbmKFParticle > &Particles, const int DaughterNegPDG, const int DaughterPosPDG, const int MotherPDG, vector< short > &idNeg, vector< short > &idPos, CbmKFVertex &PrimVtx, const float *cuts=0, bool isPrimary=0, vector< float > *vMotherTopoChi2Ndf=0, const float *secCuts=0, const float massMotherPDG=0, const float massMotherPDGSigma=0, vector< CbmKFParticle > *vMotherPrim=0, vector< CbmKFParticle > *vMotherSec=0)
Definition CbmKFParticleInterface.cxx:1567
CbmKFParticleInterface::FindParticles
static void FindParticles(vector< CbmKFTrack > &vRTracks, vector< float > &ChiToPrimVtx, vector< L1FieldRegion > &vField, vector< CbmKFParticle > &Particles, CbmKFVertex &PrimVtx, const vector< int > &vTrackPDG, const float cuts[2][3]=DefaultCuts)
Definition CbmKFParticleInterface.cxx:1186
CbmKFParticleInterface::Construct
void Construct(CbmKFTrackInterface *vDaughters[][fvecLen], int NDaughters, CbmKFVertexInterface *Parent=0, float Mass=-1, float CutChi2=-1)
Definition CbmKFParticleInterface.cxx:44
CbmKFParticleInterface::GetChi2BetweenParticles
static fvec GetChi2BetweenParticles(CbmKFParticle_simd &p1, CbmKFParticle_simd &p2)
Definition CbmKFParticleInterface.cxx:1539
CbmKFParticleInterface::Cij
fvec & Cij(Int_t i, Int_t j)
Definition CbmKFParticleInterface.h:231
CbmKFParticleInterface::SetVtxErrGuess
void SetVtxErrGuess(fvec &d_x, fvec &d_y, fvec &d_z)
Definition CbmKFParticleInterface.cxx:1170
CbmKFParticleInterface::GetMomentum
void GetMomentum(fvec &P, fvec &Error)
Definition CbmKFParticleInterface.cxx:1175
CbmKFParticleInterface::SelectParticleCandidates
static void SelectParticleCandidates(vector< CbmKFParticle > &Particles, vector< CbmKFParticle > &vCandidates, CbmKFVertex &PrimVtx, const float cuts[5])
Definition CbmKFParticleInterface.cxx:2088
CbmKFParticleInterface::FindTrackV0Decay
static void FindTrackV0Decay(const int MotherPDG, vector< CbmKFParticle > &Particles, vector< CbmKFParticle > &vV0, vector< CbmKFTrack > &vTracks, const vector< L1FieldRegion > &field, const int DaughterPDG, vector< short > &idTrack, CbmKFVertex &PrimVtx, const float *cuts=0, bool isPrimary=0, vector< float > *ChiToPrimVtx=0, vector< CbmKFParticle > *vHyperonPrim=0, float hyperonPrimMass=0, float hyperonPrimMassErr=0, vector< CbmKFParticle > *vHyperonSec=0)
Definition CbmKFParticleInterface.cxx:1777
CbmKFParticleInterface::GetKFParticle
void GetKFParticle(CbmKFParticle &Part, int iPart=0)
Definition CbmKFParticleInterface.cxx:1126
CbmKFParticleInterface::multQSQt1
void multQSQt1(const fvec J[11], fvec S[])
Definition CbmKFParticleInterface.cxx:904
CbmKFParticleInterface::FindHyperons
static void FindHyperons(int PDG, CbmKFParticle_simd vDaughters[2], vector< int > &daughterIds, vector< CbmKFParticle > &vLambdaSec, vector< CbmKFParticle > &vHyperon, CbmKFVertex &PrimVtx, const float *cuts=0, int startIndex=0)
Definition CbmKFParticleInterface.cxx:1887
CbmKFParticleInterface::operator=
CbmKFParticleInterface operator=(const CbmKFParticleInterface &c)
Definition CbmKFParticleInterface.cxx:38
CbmKFParticleInterface::GetMass
void GetMass(fvec &M, fvec &Error)
Definition CbmKFParticleInterface.cxx:1180
CbmKFParticleInterface::ExtrapolateLine
void ExtrapolateLine(CbmKFParticle_simd *Particle, fvec r0[], fvec T)
Definition CbmKFParticleInterface.cxx:693
CbmKFParticleInterface::IJ
static Int_t IJ(Int_t i, Int_t j)
Definition CbmKFParticleInterface.h:235
CbmKFParticleInterface::ExtrapolateToPV
static void ExtrapolateToPV(vector< CbmKFParticle > &vParticles, CbmKFVertex &PrimVtx)
Definition CbmKFParticleInterface.cxx:1515
CbmKFParticleInterface::CombineTrackPart
static void CombineTrackPart(vector< CbmKFTrack > &vTracks, const vector< L1FieldRegion > &vField, vector< CbmKFParticle > &Particles, CbmKFParticle &part, const int DaughterPDG, const int MotherPDG, vector< short > &id, const float *cuts, const unsigned short startIndex=0, const bool IsSamePart=0)
Definition CbmKFParticleInterface.cxx:2013
CbmKFParticleInterface::Convert
void Convert(CbmKFParticle_simd *Particle, fvec r0[], bool ToProduction)
Definition CbmKFParticleInterface.cxx:625
CbmKFParticleInterface::ConstructPVT
void ConstructPVT(vector< CbmKFTrack > &vRTracks)
Definition CbmKFParticleInterface.cxx:2124
CbmKFParticleInterface::MeasureProductionVertex
void MeasureProductionVertex(CbmKFParticle_simd *Particle, fvec r0[], CbmKFVertexInterface *Parent)
Definition CbmKFParticleInterface.cxx:473
CbmKFParticleInterface::SetVtxGuess
void SetVtxGuess(fvec &x, fvec &y, fvec &z)
Definition CbmKFParticleInterface.cxx:1165
CbmKFParticleInterface::GetDStoPoint
fvec GetDStoPoint(const CbmKFParticle_simd &Particle, const fvec xyz[]) const
Definition CbmKFParticleInterface.cxx:1030
CbmKFParticleInterface::FindDMesLambdac
static void FindDMesLambdac(vector< CbmKFTrack > &vTracks, const vector< L1FieldRegion > &vField, vector< CbmKFParticle > &Particles, const int DaughterPDG[5], const int MotherPDG[8], vector< short > *idTrack[5], CbmKFVertex &PrimVtx, const float cuts[8][8], vector< float > ChiToPrimVtx)
Definition CbmKFParticleInterface.cxx:1940
CbmKFParticleInterface::Extrapolate
void Extrapolate(CbmKFParticle_simd *Particle, fvec r0[], fvec T)
Definition CbmKFParticleInterface.cxx:742
CbmKFParticle_simd
Definition CbmKFParticle_simd.h:24
CbmKFParticle_simd::Create
void Create(CbmKFTrackInterface *Track[], int Ntracks=fvecLen, Int_t *qHypo=0, const Int_t *pdg=0)
Definition CbmKFParticle_simd.cxx:214
CbmKFParticle_simd::SetId
void SetId(fvec id)
Definition CbmKFParticle_simd.h:120
CbmKFParticle_simd::SetVtxErrGuess
void SetVtxErrGuess(fvec &d_x, fvec &d_y, fvec &d_z)
Definition CbmKFParticle_simd.cxx:434
CbmKFParticle_simd::GetChi2
fvec GetChi2() const
Definition CbmKFParticle_simd.h:56
CbmKFParticle_simd::fIsVtxErrGuess
bool fIsVtxErrGuess
Definition CbmKFParticle_simd.h:147
CbmKFParticle_simd::GetNDF
fvec GetNDF() const
Definition CbmKFParticle_simd.h:57
CbmKFParticle_simd::GetPz
fvec GetPz() const
Definition CbmKFParticle_simd.h:52
CbmKFParticle_simd::SetField
void SetField()
CbmKFParticle_simd::C
fvec C[36]
Definition CbmKFParticle_simd.h:136
CbmKFParticle_simd::GetMomentum
void GetMomentum(fvec &P, fvec &Error)
Definition CbmKFParticle_simd.cxx:369
CbmKFParticle_simd::SetNDaughters
void SetNDaughters(int n)
Definition CbmKFParticle_simd.h:121
CbmKFParticle_simd::AddDaughterId
void AddDaughterId(fvec id)
Definition CbmKFParticle_simd.h:122
CbmKFParticle_simd::GetCovariance
fvec GetCovariance(Int_t i) const
Definition CbmKFParticle_simd.h:62
CbmKFParticle_simd::GetCovMatrix
fvec * GetCovMatrix()
Definition CbmKFParticle_simd.h:60
CbmKFParticle_simd::SetVtxGuess
void SetVtxGuess(fvec &x, fvec &y, fvec &z)
Definition CbmKFParticle_simd.cxx:426
CbmKFParticle_simd::GetY
fvec GetY() const
Definition CbmKFParticle_simd.h:48
CbmKFParticle_simd::GetPDG
const fvec & GetPDG() const
Definition CbmKFParticle_simd.h:118
CbmKFParticle_simd::fField
L1FieldRegion fField
Definition CbmKFParticle_simd.h:149
CbmKFParticle_simd::fVtxErrGuess
fvec fVtxErrGuess[3]
Definition CbmKFParticle_simd.h:138
CbmKFParticle_simd::r
fvec r[8]
Definition CbmKFParticle_simd.h:135
CbmKFParticle_simd::GetMass
void GetMass(fvec &M, fvec &Error)
Definition CbmKFParticle_simd.cxx:383
CbmKFParticle_simd::fIsVtxGuess
bool fIsVtxGuess
Definition CbmKFParticle_simd.h:146
CbmKFParticle_simd::Id
fvec Id() const
Definition CbmKFParticle_simd.h:125
CbmKFParticle_simd::CleanDaughtersId
void CleanDaughtersId()
Definition CbmKFParticle_simd.h:123
CbmKFParticle_simd::DaughterIds
vector< fvec > & DaughterIds()
Definition CbmKFParticle_simd.h:127
CbmKFParticle_simd::GetParameters
fvec * GetParameters()
Definition CbmKFParticle_simd.h:58
CbmKFParticle_simd::GetQ
fvec GetQ() const
Definition CbmKFParticle_simd.h:55
CbmKFParticle_simd::AtProductionVertex
bool AtProductionVertex
Definition CbmKFParticle_simd.h:145
CbmKFParticle_simd::fVtxGuess
fvec fVtxGuess[3]
Definition CbmKFParticle_simd.h:137
CbmKFParticle_simd::GetX
fvec GetX() const
Definition CbmKFParticle_simd.h:47
CbmKFParticle_simd::GetZ
fvec GetZ() const
Definition CbmKFParticle_simd.h:49
CbmKFParticle_simd::GetPy
fvec GetPy() const
Definition CbmKFParticle_simd.h:51
CbmKFParticle::SetNDaughters
void SetNDaughters(int n)
Definition CbmKFParticle.h:27
CbmKFParticle::GetMass
void GetMass(Double_t &M, Double_t &Error)
Definition CbmKFParticle.cxx:1446
CbmKFParticle::GetParameters
Double_t * GetParameters()
Definition CbmKFParticle.h:81
CbmKFParticle::SetId
void SetId(int id)
Definition CbmKFParticle.h:26
CbmKFParticle::CleanDaughtersId
void CleanDaughtersId()
Definition CbmKFParticle.h:29
CbmKFParticle::GetY
Double_t GetY() const
Definition CbmKFParticle.h:71
CbmKFParticle::GetZ
Double_t GetZ() const
Definition CbmKFParticle.h:72
CbmKFParticle::SetPDG
void SetPDG(int pdg)
Definition CbmKFParticle.h:129
CbmKFParticle::GetX
Double_t GetX() const
Definition CbmKFParticle.h:70
CbmKFParticle::DaughterIds
const vector< int > & DaughterIds() const
Definition CbmKFParticle.h:24
CbmKFParticle::fieldRegion
float fieldRegion[10]
Definition CbmKFParticle.h:45
CbmKFParticle::AtProductionVertex
Bool_t AtProductionVertex
Definition CbmKFParticle.h:43
CbmKFParticle::AddDaughter
void AddDaughter(int id)
Definition CbmKFParticle.h:28
CbmKFParticle::GetCovMatrix
Double_t * GetCovMatrix()
Definition CbmKFParticle.h:82
CbmKFTrack::GetTrack
Double_t * GetTrack()
Is it electron.
Definition CbmKFTrack.h:51
CbmKFTrack::GetRefChi2
Double_t & GetRefChi2()
array[15] of covariance matrix
Definition CbmKFTrack.h:53
CbmKFTrack::GetRefNDF
Int_t & GetRefNDF()
Chi^2 after fit.
Definition CbmKFTrack.h:54
CbmKFTrack::GetCovMatrix
Double_t * GetCovMatrix()
array[6] of track parameters(x,y,tx,ty,qp,z)
Definition CbmKFTrack.h:52
CbmKFVertexInterface::GetRefY
virtual Double_t & GetRefY()
Definition CbmKFVertexInterface.cxx:28
CbmKFVertexInterface::GetRefX
virtual Double_t & GetRefX()
Definition CbmKFVertexInterface.cxx:24
CbmKFVertexInterface::GetRefZ
virtual Double_t & GetRefZ()
Definition CbmKFVertexInterface.cxx:32
CbmKFVertexInterface::GetCovMatrix
virtual Double_t * GetCovMatrix()
Definition CbmKFVertexInterface.cxx:36
CbmKFVertex::GetRefZ
Double_t & GetRefZ()
Definition CbmKFVertex.h:21
CbmKFVertex::GetRefX
Double_t & GetRefX()
Definition CbmKFVertex.h:19
CbmKFVertex::GetCovMatrix
Double_t * GetCovMatrix()
Definition CbmKFVertex.h:22
CbmKFVertex::GetRefChi2
Double_t & GetRefChi2()
Array[6] of covariance matrix.
Definition CbmKFVertex.h:23
CbmKFVertex::GetRefY
Double_t & GetRefY()
Definition CbmKFVertex.h:20
CbmKFVertex::GetRefNDF
Int_t & GetRefNDF()
Chi^2 after fit.
Definition CbmKFVertex.h:24
L1FieldRegion::SetOneEntry
void SetOneEntry(const int i0, const L1FieldRegion &f1, const int i1)
Definition L1Field.h:176
L1FieldRegion::Get
L1FieldValue Get(const fvec z)
Definition L1Field.h:106
L1FieldRegion::Set
void Set(const L1FieldValue &B0, const fvec B0z, const L1FieldValue &B1, const fvec B1z, const L1FieldValue &B2, const fvec B2z)
Definition L1Field.h:116
NS_L1TrackFitter::ZERO
const fvec ZERO
Definition L1TrackFitter.cxx:32
NS_L1TrackFitter::ONE
const fvec ONE
Definition L1TrackFitter.cxx:33
check_success.s
s
Definition check_success.py:47