CbmRoot
CbmKFParticle.cxx
Go to the documentation of this file.
1 //3456789012345678901234567890123456789012345678901234567890123456789012
10 #include "CbmKFParticle.h"
11 #include "CbmKF.h"
12 #include "CbmKFMath.h"
13 #include "CbmKFTrack.h"
14 #include "CbmStsKFTrackFitter.h"
15 #include "TMath.h"
16 //#include "TDatabasePDG.h"
17 
18 #include "CbmKFParticleDatabase.h"
19 
20 #include <cmath>
21 #include <vector>
22 
23 using namespace std;
24 
25 ClassImp(CbmKFParticle)
26 
27  CbmKFParticle::CbmKFParticle(CbmKFTrackInterface* Track,
28  Double_t* z0,
29  Int_t* qHypo,
30  Int_t* pdg)
31  : fId(-1)
32  , fDaughtersIds()
33  , fPDG(-1)
34  , NDF(0)
35  , Chi2(0)
36  , Q(0)
37  , AtProductionVertex(0) {
38 
39  fDaughtersIds.push_back(Track->Id());
40 
41  CbmKFTrack Tr(*Track);
42 
43  Double_t* m = Tr.GetTrack();
44  Double_t* V = Tr.GetCovMatrix();
45 
46  //* Fit of vertex track slopes and momentum (a,b,qp) to xyz0 vertex
47 
48  // if( z0 ) Tr.Extrapolate( *z0 );
49 
50  Double_t a = m[2], b = m[3], qp = m[4];
51 
52  //* convert the track to (x,y,px,py,pz,e,t/p) parameterization
53 
54  double Mass = Tr.GetMass();
55  if (pdg) {
56  Mass = CbmKFParticleDatabase::Instance()->GetMass(*pdg);
57  // Mass=TDatabasePDG::Instance()->GetParticle(*pdg)->Mass();
58  }
59 
60 
61  double c2 = 1. / (1. + a * a + b * b);
62  double pq = 1. / qp;
63  if (qHypo) pq *= *qHypo;
64  double p2 = pq * pq;
65  double pz = sqrt(p2 * c2);
66  double px = a * pz;
67  double py = b * pz;
68  double E = sqrt(Mass * Mass + p2);
69 
70  double H[3] = {-px * c2, -py * c2, -pz * pq};
71  double HE = -pq * p2 / E;
72 
73  r[0] = m[0];
74  r[1] = m[1];
75  r[2] = m[5];
76  r[3] = px;
77  r[4] = py;
78  r[5] = pz;
79  r[6] = E;
80  r[7] = 0;
81 
82  double cxpz = H[0] * V[3] + H[1] * V[6] + H[2] * V[10];
83  double cypz = H[0] * V[4] + H[1] * V[7] + H[2] * V[11];
84  double capz = H[0] * V[5] + H[1] * V[8] + H[2] * V[12];
85  double cbpz = H[0] * V[8] + H[1] * V[9] + H[2] * V[13];
86  double cqpz = H[0] * V[12] + H[1] * V[13] + H[2] * V[14];
87  double cpzpz =
88  H[0] * H[0] * V[5] + H[1] * H[1] * V[9] + H[2] * H[2] * V[14]
89  + 2 * (H[0] * H[1] * V[8] + H[0] * H[2] * V[12] + H[1] * H[2] * V[13]);
90 
91  C[0] = V[0];
92  C[1] = V[1];
93  C[2] = V[2];
94  C[3] = 0;
95  C[4] = 0;
96  C[5] = 0;
97  C[6] = V[3] * pz + a * cxpz;
98  C[7] = V[4] * pz + a * cypz;
99  C[8] = 0;
100  C[9] = V[5] * pz * pz + 2 * a * pz * capz + a * a * cpzpz;
101  C[10] = V[6] * pz + b * cxpz;
102  C[11] = V[7] * pz + b * cypz;
103  C[12] = 0;
104  C[13] = V[8] * pz * pz + a * pz * cbpz + b * pz * capz + a * b * cpzpz;
105  C[14] = V[9] * pz * pz + 2 * b * pz * cbpz + b * b * cpzpz;
106  C[15] = cxpz;
107  C[16] = cypz;
108  C[17] = 0;
109  C[18] = capz * pz + a * cpzpz;
110  C[19] = cbpz * pz + b * cpzpz;
111  C[20] = cpzpz;
112  C[21] = HE * V[10];
113  C[22] = HE * V[11];
114  C[23] = 0;
115  C[24] = HE * (V[12] * pz + a * cqpz);
116  C[25] = HE * (V[13] * pz + b * cqpz);
117  C[26] = HE * cqpz;
118  C[27] = HE * HE * V[14];
119  C[28] = C[29] = C[30] = C[31] = C[32] = C[33] = C[34] = 0;
120  C[35] = 1.;
121 
122  NDF = Tr.GetRefNDF();
123  Chi2 = Tr.GetRefChi2();
124  Q = (qp > 0.) ? 1 : ((qp < 0) ? -1 : 0);
125  if (qHypo) Q *= *qHypo;
126  AtProductionVertex = 1;
127 }
128 
129 void drawcov(double C[]) {
130  double s0 = sqrt(C[0]);
131  double s1 = sqrt(C[2]);
132  double s2 = sqrt(C[5]);
133  double s3 = sqrt(C[9]);
134  double s4 = sqrt(C[14]);
135  cout << s0 << " ";
136  cout << C[1] / s0 / s1 << " " << s1 << " ";
137  cout << C[3] / s0 / s2 << " " << C[4] / s1 / s2 << " " << s2 << " ";
138  cout << C[6] / s0 / s3 << " " << C[7] / s1 / s3 << " " << C[8] / s2 / s3
139  << " " << s3 << " ";
140  cout << C[10] / s0 / s4 << " " << C[11] / s1 / s4 << " " << C[12] / s2 / s4
141  << " " << C[13] / s3 / s4 << " " << s4 << endl;
142 }
143 
144 void CbmKFParticle::Construct(vector<CbmKFTrackInterface*>& vDaughters,
145  CbmKFVertexInterface* Parent,
146  Double_t Mass,
147  Double_t CutChi2) {
148  const Int_t MaxIter = 3;
149 
150  if (CbmKF::Instance()->vTargets.empty()) {
151  r[0] = r[1] = r[2] = 0.;
152  C[0] = C[2] = C[5] = 1.;
153  C[1] = C[3] = C[4] = 0;
154  } else {
155  CbmKFTube& t = CbmKF::Instance()->vTargets[0];
156  r[0] = r[1] = 0.;
157  r[2] = t.z;
158  C[0] = C[2] = t.R * t.R / 9.;
159  C[5] = t.dz / 6.;
160  C[1] = C[3] = C[4] = 0;
161  }
162 
163  r[3] = 0;
164  r[4] = 0;
165  r[5] = 0;
166  r[6] = 0;
167  r[7] = 0;
168 
169  for (Int_t iteration = 0; iteration < MaxIter; ++iteration) {
170 
171  Double_t r0[8], C0[6];
172 
173  for (int i = 0; i < 8; i++)
174  r0[i] = r[i];
175  for (int i = 0; i < 6; i++)
176  C0[i] = C[i];
177 
178  //simple approximation for the vertex position. TODO find better solution
179  if (iteration == 0) {
180  Double_t VertexGuess[3];
181 
182  Double_t fTDaughter[(int) vDaughters.size()][6];
183  Double_t fCDaughter[(int) vDaughters.size()][15];
184 
185  vector<CbmKFTrack*> TrV;
186 
187  Int_t nvect = 0;
188  for (vector<CbmKFTrackInterface*>::iterator tr = vDaughters.begin();
189  tr != vDaughters.end();
190  ++tr) {
191  CbmKFTrack* Tr = new CbmKFTrack(**tr);
192  fTDaughter[nvect][0] = Tr->GetTrack()[0];
193  fTDaughter[nvect][1] = Tr->GetTrack()[1];
194  fTDaughter[nvect][2] = Tr->GetTrack()[2];
195  fTDaughter[nvect][3] = Tr->GetTrack()[3];
196  fTDaughter[nvect][4] = Tr->GetTrack()[4];
197  fTDaughter[nvect][5] = Tr->GetTrack()[5];
198 
199  fCDaughter[nvect][0] = Tr->GetCovMatrix()[0];
200  fCDaughter[nvect][2] = Tr->GetCovMatrix()[2];
201  fCDaughter[nvect][5] = Tr->GetCovMatrix()[5];
202  fCDaughter[nvect][9] = Tr->GetCovMatrix()[9];
203  fCDaughter[nvect][14] = Tr->GetCovMatrix()[14];
204 
205  delete Tr;
206  nvect++;
207  }
208 
209  Double_t z = (fTDaughter[0][3] * fTDaughter[0][5]
210  - fTDaughter[1][3] * fTDaughter[1][5] + fTDaughter[1][1]
211  - fTDaughter[0][1])
212  / (fTDaughter[0][3] - fTDaughter[1][3]);
213 
214  CbmKFTrack* Tr = new CbmKFTrack(*vDaughters[1]);
215  Tr->Extrapolate(z);
216 
217  VertexGuess[0] = Tr->GetTrack()[0];
218  VertexGuess[1] = Tr->GetTrack()[1];
219  VertexGuess[2] = Tr->GetTrack()[5];
220 
221  Double_t dr0[2];
222  dr0[0] = sqrt(Tr->GetCovMatrix()[0]);
223  dr0[1] = sqrt(Tr->GetCovMatrix()[2]);
224  delete Tr;
225  r0[0] = VertexGuess[0];
226  r0[1] = VertexGuess[1];
227  r0[2] = VertexGuess[2];
228  }
229 
230  r[3] = 0;
231  r[4] = 0;
232  r[5] = 0;
233  r[6] = 0;
234 
235  for (Int_t i = 0; i < 36; ++i)
236  C[i] = 0.;
237 
238  C[0] = C[2] = C[5] = 100.;
239  C[35] = 1.;
240 
241  Double_t B[3];
242  {
243  FairField* MF = CbmKF::Instance()->GetMagneticField();
244  MF->GetFieldValue(r0, B);
245  const Double_t c_light = 0.000299792458;
246  B[0] *= c_light;
247  B[1] *= c_light;
248  B[2] *= c_light;
249  }
250 
251  NDF = -3;
252  Chi2 = 0.;
253  Q = 0;
254  bool first = 1;
255  for (vector<CbmKFTrackInterface*>::iterator tr = vDaughters.begin();
256  tr != vDaughters.end();
257  ++tr) {
258  CbmKFParticle Daughter(*tr, &(r0[2]));
259 
260  Daughter.Extrapolate(Daughter.r, Daughter.GetDStoPoint(r0));
261 
262  Double_t* m = Daughter.r;
263  Double_t* Cd = Daughter.C;
264 
265  Double_t d[3] = {r0[0] - m[0], r0[1] - m[1], r0[2] - m[2]};
266  Double_t SigmaS = .1
267  + 10.
268  * sqrt((d[0] * d[0] + d[1] * d[1] + d[2] * d[2])
269  / (m[3] * m[3] + m[4] * m[4] + m[5] * m[5]));
270 
271  Double_t h[6];
272 
273  h[0] = m[3] * SigmaS;
274  h[1] = m[4] * SigmaS;
275  h[2] = m[5] * SigmaS;
276  h[3] = (h[1] * B[2] - h[2] * B[1]) * Daughter.Q;
277  h[4] = (h[2] * B[0] - h[0] * B[2]) * Daughter.Q;
278  h[5] = (h[0] * B[1] - h[1] * B[0]) * Daughter.Q;
279 
280  //* Fit of daughter momentum (x,y,z) to r0[0,1,2] vertex
281  {
282  Double_t zeta[3] = {r0[0] - m[0], r0[1] - m[1], r0[2] - m[2]};
283 
284  Double_t Vv[6] = {Cd[0] + h[0] * h[0],
285  Cd[1] + h[1] * h[0],
286  Cd[2] + h[1] * h[1],
287  Cd[3] + h[2] * h[0],
288  Cd[4] + h[2] * h[1],
289  Cd[5] + h[2] * h[2]};
290 
291  Double_t Vvp[9] = {Cd[6] + h[0] * h[3],
292  Cd[7] + h[1] * h[3],
293  Cd[8] + h[2] * h[3],
294  Cd[10] + h[0] * h[4],
295  Cd[11] + h[1] * h[4],
296  Cd[12] + h[2] * h[4],
297  Cd[15] + h[0] * h[5],
298  Cd[16] + h[1] * h[5],
299  Cd[17] + h[2] * h[5]};
300 
301  if (CutChi2 > 0.) {
302 
303  Double_t Si[6] = {Vv[0] + C0[0],
304  Vv[1] + C0[1],
305  Vv[2] + C0[2],
306  Vv[3] + C0[3],
307  Vv[4] + C0[4],
308  Vv[5] + C0[5]};
309  Double_t S[6] = {Si[2] * Si[5] - Si[4] * Si[4],
310  Si[3] * Si[4] - Si[1] * Si[5],
311  Si[0] * Si[5] - Si[3] * Si[3],
312  Si[1] * Si[4] - Si[2] * Si[3],
313  Si[1] * Si[3] - Si[0] * Si[4],
314  Si[0] * Si[2] - Si[1] * Si[1]};
315  Double_t det = (Si[0] * S[0] + Si[1] * S[1] + Si[3] * S[3]);
316  Double_t chi2 =
317  (+(S[0] * zeta[0] + S[1] * zeta[1] + S[3] * zeta[2]) * zeta[0]
318  + (S[1] * zeta[0] + S[2] * zeta[1] + S[4] * zeta[2]) * zeta[1]
319  + (S[3] * zeta[0] + S[4] * zeta[1] + S[5] * zeta[2]) * zeta[2]);
320  if (chi2 > 2 * CutChi2 * det) continue;
321  }
322 
323  double S[6] = {Vv[2] * Vv[5] - Vv[4] * Vv[4],
324  Vv[3] * Vv[4] - Vv[1] * Vv[5],
325  Vv[0] * Vv[5] - Vv[3] * Vv[3],
326  Vv[1] * Vv[4] - Vv[2] * Vv[3],
327  Vv[1] * Vv[3] - Vv[0] * Vv[4],
328  Vv[0] * Vv[2] - Vv[1] * Vv[1]};
329 
330  double s = (Vv[0] * S[0] + Vv[1] * S[1] + Vv[3] * S[3]);
331  s = (s > 1.E-20) ? 1. / s : 0;
332 
333  S[0] *= s;
334  S[1] *= s;
335  S[2] *= s;
336  S[3] *= s;
337  S[4] *= s;
338  S[5] *= s;
339 
340  Double_t Sz[3] = {(S[0] * zeta[0] + S[1] * zeta[1] + S[3] * zeta[2]),
341  (S[1] * zeta[0] + S[2] * zeta[1] + S[4] * zeta[2]),
342  (S[3] * zeta[0] + S[4] * zeta[1] + S[5] * zeta[2])};
343 
344  Double_t x = m[3] + Vvp[0] * Sz[0] + Vvp[1] * Sz[1] + Vvp[2] * Sz[2];
345  Double_t y = m[4] + Vvp[3] * Sz[0] + Vvp[4] * Sz[1] + Vvp[5] * Sz[2];
346  Double_t z = m[5] + Vvp[6] * Sz[0] + Vvp[7] * Sz[1] + Vvp[8] * Sz[2];
347 
348  h[0] = x * SigmaS;
349  h[1] = y * SigmaS;
350  h[2] = z * SigmaS;
351  h[3] = (h[1] * B[2] - h[2] * B[1]) * Daughter.Q;
352  h[4] = (h[2] * B[0] - h[0] * B[2]) * Daughter.Q;
353  h[5] = (h[0] * B[1] - h[1] * B[0]) * Daughter.Q;
354  }
355 
356  Double_t V[28];
357 
358  V[0] = Cd[0] + h[0] * h[0];
359  V[1] = Cd[1] + h[1] * h[0];
360  V[2] = Cd[2] + h[1] * h[1];
361  V[3] = Cd[3] + h[2] * h[0];
362  V[4] = Cd[4] + h[2] * h[1];
363  V[5] = Cd[5] + h[2] * h[2];
364 
365  V[6] = Cd[6] + h[3] * h[0];
366  V[7] = Cd[7] + h[3] * h[1];
367  V[8] = Cd[8] + h[3] * h[2];
368  V[9] = Cd[9] + h[3] * h[3];
369 
370  V[10] = Cd[10] + h[4] * h[0];
371  V[11] = Cd[11] + h[4] * h[1];
372  V[12] = Cd[12] + h[4] * h[2];
373  V[13] = Cd[13] + h[4] * h[3];
374  V[14] = Cd[14] + h[4] * h[4];
375 
376  V[15] = Cd[15] + h[5] * h[0];
377  V[16] = Cd[16] + h[5] * h[1];
378  V[17] = Cd[17] + h[5] * h[2];
379  V[18] = Cd[18] + h[5] * h[3];
380  V[19] = Cd[19] + h[5] * h[4];
381  V[20] = Cd[20] + h[5] * h[5];
382 
383  V[21] = Cd[21];
384  V[22] = Cd[22];
385  V[23] = Cd[23];
386  V[24] = Cd[24];
387  V[25] = Cd[25];
388  V[26] = Cd[26];
389  V[27] = Cd[27];
390 
391  //*
392 
393  NDF += 2;
394  Q += Daughter.Q;
395 
396  if (0 && first) {
397  first = 0;
398  for (int i = 0; i < 7; i++)
399  r[i] = m[i];
400  for (int i = 0; i < 28; i++)
401  C[i] = V[i];
402  continue;
403  }
404 
405  //*
406 
407  double S[6];
408  {
409  double Si[6] = {C[0] + V[0],
410  C[1] + V[1],
411  C[2] + V[2],
412  C[3] + V[3],
413  C[4] + V[4],
414  C[5] + V[5]};
415 
416  S[0] = Si[2] * Si[5] - Si[4] * Si[4];
417  S[1] = Si[3] * Si[4] - Si[1] * Si[5];
418  S[2] = Si[0] * Si[5] - Si[3] * Si[3];
419  S[3] = Si[1] * Si[4] - Si[2] * Si[3];
420  S[4] = Si[1] * Si[3] - Si[0] * Si[4];
421  S[5] = Si[0] * Si[2] - Si[1] * Si[1];
422 
423  double s = (Si[0] * S[0] + Si[1] * S[1] + Si[3] * S[3]);
424  s = (s > 1.E-20) ? 1. / s : 0;
425  S[0] *= s;
426  S[1] *= s;
427  S[2] *= s;
428  S[3] *= s;
429  S[4] *= s;
430  S[5] *= s;
431  }
432 
433  //* Residual (measured - estimated)
434 
435  Double_t zeta[3] = {m[0] - r[0], m[1] - r[1], m[2] - r[2]};
436 
437  //* Add the daughter momentum to the particle momentum
438 
439  r[3] += m[3];
440  r[4] += m[4];
441  r[5] += m[5];
442  r[6] += m[6];
443 
444  C[9] += V[9];
445  C[13] += V[13];
446  C[14] += V[14];
447  C[18] += V[18];
448  C[19] += V[19];
449  C[20] += V[20];
450  C[24] += V[24];
451  C[25] += V[25];
452  C[26] += V[26];
453  C[27] += V[27];
454 
455  //* CHt = CH' - D'
456 
457  Double_t CHt0[7], CHt1[7], CHt2[7];
458 
459  CHt0[0] = C[0];
460  CHt1[0] = C[1];
461  CHt2[0] = C[3];
462  CHt0[1] = C[1];
463  CHt1[1] = C[2];
464  CHt2[1] = C[4];
465  CHt0[2] = C[3];
466  CHt1[2] = C[4];
467  CHt2[2] = C[5];
468  CHt0[3] = C[6] - V[6];
469  CHt1[3] = C[7] - V[7];
470  CHt2[3] = C[8] - V[8];
471  CHt0[4] = C[10] - V[10];
472  CHt1[4] = C[11] - V[11];
473  CHt2[4] = C[12] - V[12];
474  CHt0[5] = C[15] - V[15];
475  CHt1[5] = C[16] - V[16];
476  CHt2[5] = C[17] - V[17];
477  CHt0[6] = C[21] - V[21];
478  CHt1[6] = C[22] - V[22];
479  CHt2[6] = C[23] - V[23];
480 
481  //* Kalman gain K = CH'*S
482 
483  Double_t K0[7], K1[7], K2[7];
484 
485  for (Int_t i = 0; i < 7; ++i) {
486  K0[i] = CHt0[i] * S[0] + CHt1[i] * S[1] + CHt2[i] * S[3];
487  K1[i] = CHt0[i] * S[1] + CHt1[i] * S[2] + CHt2[i] * S[4];
488  K2[i] = CHt0[i] * S[3] + CHt1[i] * S[4] + CHt2[i] * S[5];
489  }
490 
491  //* New estimation of the vertex position r += K*zeta
492 
493  for (Int_t i = 0; i < 7; ++i)
494  r[i] += K0[i] * zeta[0] + K1[i] * zeta[1] + K2[i] * zeta[2];
495 
496  //* New covariance matrix C -= K*(CH')'
497 
498  for (Int_t i = 0, k = 0; i < 7; ++i) {
499  for (Int_t j = 0; j <= i; ++j, ++k)
500  C[k] -= K0[i] * CHt0[j] + K1[i] * CHt1[j] + K2[i] * CHt2[j];
501  }
502 
503  //* Calculate Chi^2 & NDF
504 
505  Chi2 += (S[0] * zeta[0] + S[1] * zeta[1] + S[3] * zeta[2]) * zeta[0]
506  + (S[1] * zeta[0] + S[2] * zeta[1] + S[4] * zeta[2]) * zeta[1]
507  + (S[3] * zeta[0] + S[4] * zeta[1] + S[5] * zeta[2]) * zeta[2];
508 
509  } // add tracks
510 
511  if (iteration == 0) {
512  r0[3] = r[3];
513  r0[4] = r[4];
514  r0[5] = r[5];
515  r0[6] = r[6];
516  if (Parent) {
517  double dx = Parent->GetRefX() - r[0];
518  double dy = Parent->GetRefY() - r[1];
519  double dz = Parent->GetRefZ() - r[2];
520  r0[7] = r[7] = sqrt((dx * dx + dy * dy + dz * dz)
521  / (r[3] * r[3] + r[4] * r[4] + r[5] * r[5]));
522  }
523  }
524 
525  if (Mass >= 0) MeasureMass(r0, Mass);
526  if (Parent) MeasureProductionVertex(r0, Parent);
527 
528  } // iterations
529 
530  AtProductionVertex = 0;
531 }
532 
533 void CbmKFParticle::ConstructFromKFParticle(vector<CbmKFParticle*>& vDaughters,
534  CbmKFVertexInterface* Parent,
535  Double_t Mass,
536  Double_t CutChi2) {
537  const Int_t MaxIter = 3;
538 
539  if (CbmKF::Instance()->vTargets.empty()) {
540  r[0] = r[1] = r[2] = 0.;
541  C[0] = C[2] = C[5] = 1.;
542  C[1] = C[3] = C[4] = 0;
543  } else {
544  CbmKFTube& t = CbmKF::Instance()->vTargets[0];
545  r[0] = r[1] = 0.;
546  r[2] = t.z;
547  C[0] = C[2] = t.R * t.R / 9.;
548  C[5] = t.dz / 6.;
549  C[1] = C[3] = C[4] = 0;
550  }
551 
552  r[3] = 0;
553  r[4] = 0;
554  r[5] = 0;
555  r[6] = 0;
556  r[7] = 0;
557 
558  for (Int_t iteration = 0; iteration < MaxIter; ++iteration) {
559 
560  Double_t r0[8], C0[6];
561 
562  for (int i = 0; i < 8; i++)
563  r0[i] = r[i];
564  for (int i = 0; i < 6; i++)
565  C0[i] = C[i];
566 
567  if (iteration == 0) {
568  Double_t VertexGuess[3];
569 
570  Double_t fTDaughter[(int) vDaughters.size()][6];
571  Double_t fCDaughter[(int) vDaughters.size()][15];
572 
573  vector<CbmKFTrack*> TrV;
574 
575  Int_t nvect = 0;
576  for (vector<CbmKFParticle*>::iterator tr = vDaughters.begin();
577  tr != vDaughters.end();
578  ++tr) {
579  CbmKFParticle* Part = *tr;
580  CbmKFTrackInterface* TrInt = new CbmKFTrackInterface();
581  Part->GetKFTrack(TrInt);
582  CbmKFTrack* Tr = new CbmKFTrack(*TrInt);
583 
584  fTDaughter[nvect][0] = Tr->GetTrack()[0];
585  fTDaughter[nvect][1] = Tr->GetTrack()[1];
586  fTDaughter[nvect][2] = Tr->GetTrack()[2];
587  fTDaughter[nvect][3] = Tr->GetTrack()[3];
588  fTDaughter[nvect][4] = Tr->GetTrack()[4];
589  fTDaughter[nvect][5] = Tr->GetTrack()[5];
590 
591  fCDaughter[nvect][0] = Tr->GetCovMatrix()[0];
592  fCDaughter[nvect][2] = Tr->GetCovMatrix()[2];
593  fCDaughter[nvect][5] = Tr->GetCovMatrix()[5];
594  fCDaughter[nvect][9] = Tr->GetCovMatrix()[9];
595  fCDaughter[nvect][14] = Tr->GetCovMatrix()[14];
596 
597  delete Tr;
598  delete TrInt;
599  nvect++;
600  }
601 
602  Double_t z = (fTDaughter[0][3] * fTDaughter[0][5]
603  - fTDaughter[1][3] * fTDaughter[1][5] + fTDaughter[1][1]
604  - fTDaughter[0][1])
605  / (fTDaughter[0][3] - fTDaughter[1][3]);
606 
607  CbmKFParticle* Part = vDaughters[0];
608  CbmKFTrackInterface* TrInt = new CbmKFTrackInterface();
609  Part->GetKFTrack(TrInt);
610  CbmKFTrack* Tr = new CbmKFTrack(*TrInt);
611 
612  Tr->Extrapolate(z);
613 
614  VertexGuess[0] = Tr->GetTrack()[0];
615  VertexGuess[1] = Tr->GetTrack()[1];
616  VertexGuess[2] = Tr->GetTrack()[5];
617  delete Tr;
618  r0[0] = VertexGuess[0];
619  r0[1] = VertexGuess[1];
620  r0[2] = VertexGuess[2];
621  }
622 
623  r[3] = 0;
624  r[4] = 0;
625  r[5] = 0;
626  r[6] = 0;
627 
628  for (Int_t i = 0; i < 36; ++i)
629  C[i] = 0.;
630 
631  C[0] = C[2] = C[5] = 100.;
632  C[35] = 1.;
633 
634  Double_t B[3];
635  {
636  FairField* MF = CbmKF::Instance()->GetMagneticField();
637  MF->GetFieldValue(r0, B);
638  const Double_t c_light = 0.000299792458;
639  B[0] *= c_light;
640  B[1] *= c_light;
641  B[2] *= c_light;
642  }
643 
644  NDF = -3;
645  Chi2 = 0.;
646  Q = 0;
647  bool first = 1;
648  for (vector<CbmKFParticle*>::iterator tr = vDaughters.begin();
649  tr != vDaughters.end();
650  ++tr) {
651  CbmKFParticle Daughter = *(*tr);
652 
653  Double_t dx_aprox = Daughter.r[0] - r0[0];
654  Double_t dy_aprox = Daughter.r[1] - r0[1];
655  Double_t dz_aprox = Daughter.r[2] - r0[2];
656  Double_t dS_aprox =
657  sqrt((dx_aprox * dx_aprox + dy_aprox * dy_aprox + dz_aprox * dz_aprox)
658  / (Daughter.r[3] * Daughter.r[3] + Daughter.r[4] * Daughter.r[4]
659  + Daughter.r[5] * Daughter.r[5]));
660 
661  Daughter.Extrapolate(Daughter.r, -dS_aprox);
662 
663  Double_t* m = Daughter.r;
664  Double_t* Cd = Daughter.C;
665 
666  Double_t d[3] = {r0[0] - m[0], r0[1] - m[1], r0[2] - m[2]};
667  Double_t SigmaS = .1
668  + 10.
669  * sqrt((d[0] * d[0] + d[1] * d[1] + d[2] * d[2])
670  / (m[3] * m[3] + m[4] * m[4] + m[5] * m[5]));
671 
672  Double_t h[6];
673 
674  h[0] = m[3] * SigmaS;
675  h[1] = m[4] * SigmaS;
676  h[2] = m[5] * SigmaS;
677  h[3] = (h[1] * B[2] - h[2] * B[1]) * Daughter.Q;
678  h[4] = (h[2] * B[0] - h[0] * B[2]) * Daughter.Q;
679  h[5] = (h[0] * B[1] - h[1] * B[0]) * Daughter.Q;
680 
681  //* Fit of daughter momentum (x,y,z) to r0[0,1,2] vertex
682  {
683  Double_t zeta[3] = {r0[0] - m[0], r0[1] - m[1], r0[2] - m[2]};
684 
685  Double_t Vv[6] = {Cd[0] + h[0] * h[0],
686  Cd[1] + h[1] * h[0],
687  Cd[2] + h[1] * h[1],
688  Cd[3] + h[2] * h[0],
689  Cd[4] + h[2] * h[1],
690  Cd[5] + h[2] * h[2]};
691 
692  Double_t Vvp[9] = {Cd[6] + h[0] * h[3],
693  Cd[7] + h[1] * h[3],
694  Cd[8] + h[2] * h[3],
695  Cd[10] + h[0] * h[4],
696  Cd[11] + h[1] * h[4],
697  Cd[12] + h[2] * h[4],
698  Cd[15] + h[0] * h[5],
699  Cd[16] + h[1] * h[5],
700  Cd[17] + h[2] * h[5]};
701 
702  if (CutChi2 > 0.) {
703 
704  Double_t Si[6] = {Vv[0] + C0[0],
705  Vv[1] + C0[1],
706  Vv[2] + C0[2],
707  Vv[3] + C0[3],
708  Vv[4] + C0[4],
709  Vv[5] + C0[5]};
710  Double_t S[6] = {Si[2] * Si[5] - Si[4] * Si[4],
711  Si[3] * Si[4] - Si[1] * Si[5],
712  Si[0] * Si[5] - Si[3] * Si[3],
713  Si[1] * Si[4] - Si[2] * Si[3],
714  Si[1] * Si[3] - Si[0] * Si[4],
715  Si[0] * Si[2] - Si[1] * Si[1]};
716  Double_t det = (Si[0] * S[0] + Si[1] * S[1] + Si[3] * S[3]);
717  Double_t chi2 =
718  (+(S[0] * zeta[0] + S[1] * zeta[1] + S[3] * zeta[2]) * zeta[0]
719  + (S[1] * zeta[0] + S[2] * zeta[1] + S[4] * zeta[2]) * zeta[1]
720  + (S[3] * zeta[0] + S[4] * zeta[1] + S[5] * zeta[2]) * zeta[2]);
721  if (chi2 > 2 * CutChi2 * det) continue;
722  }
723 
724  double S[6] = {Vv[2] * Vv[5] - Vv[4] * Vv[4],
725  Vv[3] * Vv[4] - Vv[1] * Vv[5],
726  Vv[0] * Vv[5] - Vv[3] * Vv[3],
727  Vv[1] * Vv[4] - Vv[2] * Vv[3],
728  Vv[1] * Vv[3] - Vv[0] * Vv[4],
729  Vv[0] * Vv[2] - Vv[1] * Vv[1]};
730 
731  double s = (Vv[0] * S[0] + Vv[1] * S[1] + Vv[3] * S[3]);
732  s = (s > 1.E-20) ? 1. / s : 0;
733 
734  S[0] *= s;
735  S[1] *= s;
736  S[2] *= s;
737  S[3] *= s;
738  S[4] *= s;
739  S[5] *= s;
740 
741  Double_t Sz[3] = {(S[0] * zeta[0] + S[1] * zeta[1] + S[3] * zeta[2]),
742  (S[1] * zeta[0] + S[2] * zeta[1] + S[4] * zeta[2]),
743  (S[3] * zeta[0] + S[4] * zeta[1] + S[5] * zeta[2])};
744 
745  Double_t x = m[3] + Vvp[0] * Sz[0] + Vvp[1] * Sz[1] + Vvp[2] * Sz[2];
746  Double_t y = m[4] + Vvp[3] * Sz[0] + Vvp[4] * Sz[1] + Vvp[5] * Sz[2];
747  Double_t z = m[5] + Vvp[6] * Sz[0] + Vvp[7] * Sz[1] + Vvp[8] * Sz[2];
748 
749  h[0] = x * SigmaS;
750  h[1] = y * SigmaS;
751  h[2] = z * SigmaS;
752  h[3] = (h[1] * B[2] - h[2] * B[1]) * Daughter.Q;
753  h[4] = (h[2] * B[0] - h[0] * B[2]) * Daughter.Q;
754  h[5] = (h[0] * B[1] - h[1] * B[0]) * Daughter.Q;
755  }
756 
757  Double_t V[28];
758 
759  V[0] = Cd[0] + h[0] * h[0];
760  V[1] = Cd[1] + h[1] * h[0];
761  V[2] = Cd[2] + h[1] * h[1];
762  V[3] = Cd[3] + h[2] * h[0];
763  V[4] = Cd[4] + h[2] * h[1];
764  V[5] = Cd[5] + h[2] * h[2];
765 
766  V[6] = Cd[6] + h[3] * h[0];
767  V[7] = Cd[7] + h[3] * h[1];
768  V[8] = Cd[8] + h[3] * h[2];
769  V[9] = Cd[9] + h[3] * h[3];
770 
771  V[10] = Cd[10] + h[4] * h[0];
772  V[11] = Cd[11] + h[4] * h[1];
773  V[12] = Cd[12] + h[4] * h[2];
774  V[13] = Cd[13] + h[4] * h[3];
775  V[14] = Cd[14] + h[4] * h[4];
776 
777  V[15] = Cd[15] + h[5] * h[0];
778  V[16] = Cd[16] + h[5] * h[1];
779  V[17] = Cd[17] + h[5] * h[2];
780  V[18] = Cd[18] + h[5] * h[3];
781  V[19] = Cd[19] + h[5] * h[4];
782  V[20] = Cd[20] + h[5] * h[5];
783 
784  V[21] = Cd[21];
785  V[22] = Cd[22];
786  V[23] = Cd[23];
787  V[24] = Cd[24];
788  V[25] = Cd[25];
789  V[26] = Cd[26];
790  V[27] = Cd[27];
791 
792  //*
793 
794  NDF += 2;
795  Q += Daughter.Q;
796 
797  if (0 && first) {
798  first = 0;
799  for (int i = 0; i < 7; i++)
800  r[i] = m[i];
801  for (int i = 0; i < 28; i++)
802  C[i] = V[i];
803  continue;
804  }
805 
806  //*
807 
808  double S[6];
809  {
810  double Si[6] = {C[0] + V[0],
811  C[1] + V[1],
812  C[2] + V[2],
813  C[3] + V[3],
814  C[4] + V[4],
815  C[5] + V[5]};
816 
817  S[0] = Si[2] * Si[5] - Si[4] * Si[4];
818  S[1] = Si[3] * Si[4] - Si[1] * Si[5];
819  S[2] = Si[0] * Si[5] - Si[3] * Si[3];
820  S[3] = Si[1] * Si[4] - Si[2] * Si[3];
821  S[4] = Si[1] * Si[3] - Si[0] * Si[4];
822  S[5] = Si[0] * Si[2] - Si[1] * Si[1];
823 
824  double s = (Si[0] * S[0] + Si[1] * S[1] + Si[3] * S[3]);
825  s = (s > 1.E-20) ? 1. / s : 0;
826  S[0] *= s;
827  S[1] *= s;
828  S[2] *= s;
829  S[3] *= s;
830  S[4] *= s;
831  S[5] *= s;
832  }
833 
834  //* Residual (measured - estimated)
835 
836  Double_t zeta[3] = {m[0] - r[0], m[1] - r[1], m[2] - r[2]};
837 
838  //* Add the daughter momentum to the particle momentum
839 
840  r[3] += m[3];
841  r[4] += m[4];
842  r[5] += m[5];
843  r[6] += m[6];
844 
845  C[9] += V[9];
846  C[13] += V[13];
847  C[14] += V[14];
848  C[18] += V[18];
849  C[19] += V[19];
850  C[20] += V[20];
851  C[24] += V[24];
852  C[25] += V[25];
853  C[26] += V[26];
854  C[27] += V[27];
855 
856  //* CHt = CH' - D'
857 
858  Double_t CHt0[7], CHt1[7], CHt2[7];
859 
860  CHt0[0] = C[0];
861  CHt1[0] = C[1];
862  CHt2[0] = C[3];
863  CHt0[1] = C[1];
864  CHt1[1] = C[2];
865  CHt2[1] = C[4];
866  CHt0[2] = C[3];
867  CHt1[2] = C[4];
868  CHt2[2] = C[5];
869  CHt0[3] = C[6] - V[6];
870  CHt1[3] = C[7] - V[7];
871  CHt2[3] = C[8] - V[8];
872  CHt0[4] = C[10] - V[10];
873  CHt1[4] = C[11] - V[11];
874  CHt2[4] = C[12] - V[12];
875  CHt0[5] = C[15] - V[15];
876  CHt1[5] = C[16] - V[16];
877  CHt2[5] = C[17] - V[17];
878  CHt0[6] = C[21] - V[21];
879  CHt1[6] = C[22] - V[22];
880  CHt2[6] = C[23] - V[23];
881 
882  //* Kalman gain K = CH'*S
883 
884  Double_t K0[7], K1[7], K2[7];
885 
886  for (Int_t i = 0; i < 7; ++i) {
887  K0[i] = CHt0[i] * S[0] + CHt1[i] * S[1] + CHt2[i] * S[3];
888  K1[i] = CHt0[i] * S[1] + CHt1[i] * S[2] + CHt2[i] * S[4];
889  K2[i] = CHt0[i] * S[3] + CHt1[i] * S[4] + CHt2[i] * S[5];
890  }
891 
892  //* New estimation of the vertex position r += K*zeta
893 
894  for (Int_t i = 0; i < 7; ++i)
895  r[i] += K0[i] * zeta[0] + K1[i] * zeta[1] + K2[i] * zeta[2];
896 
897  //* New covariance matrix C -= K*(CH')'
898 
899  for (Int_t i = 0, k = 0; i < 7; ++i) {
900  for (Int_t j = 0; j <= i; ++j, ++k)
901  C[k] -= K0[i] * CHt0[j] + K1[i] * CHt1[j] + K2[i] * CHt2[j];
902  }
903 
904  //* Calculate Chi^2 & NDF
905 
906  Chi2 += (S[0] * zeta[0] + S[1] * zeta[1] + S[3] * zeta[2]) * zeta[0]
907  + (S[1] * zeta[0] + S[2] * zeta[1] + S[4] * zeta[2]) * zeta[1]
908  + (S[3] * zeta[0] + S[4] * zeta[1] + S[5] * zeta[2]) * zeta[2];
909 
910  } // add tracks
911 
912  if (iteration == 0) {
913  r0[3] = r[3];
914  r0[4] = r[4];
915  r0[5] = r[5];
916  r0[6] = r[6];
917  if (Parent) {
918  double dx = Parent->GetRefX() - r[0];
919  double dy = Parent->GetRefY() - r[1];
920  double dz = Parent->GetRefZ() - r[2];
921  r0[7] = r[7] = sqrt((dx * dx + dy * dy + dz * dz)
922  / (r[3] * r[3] + r[4] * r[4] + r[5] * r[5]));
923  }
924  }
925 
926  if (Mass >= 0) MeasureMass(r0, Mass);
927  if (Parent) MeasureProductionVertex(r0, Parent);
928 
929  } // iterations
930 
931  AtProductionVertex = 0;
932 }
933 
934 void CbmKFParticle::MeasureMass(Double_t r0[], Double_t Mass) {
935  double H[8];
936  H[0] = H[1] = H[2] = 0.;
937  H[3] = -2 * r0[3];
938  H[4] = -2 * r0[4];
939  H[5] = -2 * r0[5];
940  H[6] = 2 * r0[6];
941  H[7] = 0;
942  double m2 = r0[6] * r0[6] - r0[3] * r0[3] - r0[4] * r0[4] - r0[5] * r0[5];
943 
944  double zeta = Mass * Mass - m2;
945  for (Int_t i = 0; i < 8; ++i)
946  zeta -= H[i] * (r[i] - r0[i]);
947 
948  Double_t S = 0.;
949  Double_t CHt[8];
950  for (Int_t i = 0; i < 8; ++i) {
951  CHt[i] = 0.0;
952  for (Int_t j = 0; j < 8; ++j)
953  CHt[i] += Cij(i, j) * H[j];
954  S += H[i] * CHt[i];
955  }
956 
957  if (S < 1.e-20) return;
958  S = 1. / S;
959  Chi2 += zeta * zeta * S;
960  NDF += 1;
961  for (Int_t i = 0, ii = 0; i < 8; ++i) {
962  Double_t Ki = CHt[i] * S;
963  r[i] += Ki * zeta;
964  for (Int_t j = 0; j <= i; ++j)
965  C[ii++] -= Ki * CHt[j];
966  }
967 }
968 
969 
970 void CbmKFParticle::MeasureProductionVertex(Double_t r0[],
971  CbmKFVertexInterface* Parent) {
972 
973  double m[3], *V;
974  {
975  m[0] = Parent->GetRefX();
976  m[1] = Parent->GetRefY();
977  m[2] = Parent->GetRefZ();
978  V = Parent->GetCovMatrix();
979  }
980  r[7] = r0[7];
981  Extrapolate(r0, -r0[7]);
982  Convert(r0, 1);
983 
984  double Ai[6];
985  Ai[0] = C[4] * C[4] - C[2] * C[5];
986  Ai[1] = C[1] * C[5] - C[3] * C[4];
987  Ai[3] = C[2] * C[3] - C[1] * C[4];
988  double det = 1. / (C[0] * Ai[0] + C[1] * Ai[1] + C[3] * Ai[3]);
989  Ai[0] *= det;
990  Ai[1] *= det;
991  Ai[3] *= det;
992  Ai[2] = (C[3] * C[3] - C[0] * C[5]) * det;
993  Ai[4] = (C[0] * C[4] - C[1] * C[3]) * det;
994  Ai[5] = (C[1] * C[1] - C[0] * C[2]) * det;
995 
996  double B[5][3];
997 
998  B[0][0] = C[6] * Ai[0] + C[7] * Ai[1] + C[8] * Ai[3];
999  B[0][1] = C[6] * Ai[1] + C[7] * Ai[2] + C[8] * Ai[4];
1000  B[0][2] = C[6] * Ai[3] + C[7] * Ai[4] + C[8] * Ai[5];
1001 
1002  B[1][0] = C[10] * Ai[0] + C[11] * Ai[1] + C[12] * Ai[3];
1003  B[1][1] = C[10] * Ai[1] + C[11] * Ai[2] + C[12] * Ai[4];
1004  B[1][2] = C[10] * Ai[3] + C[11] * Ai[4] + C[12] * Ai[5];
1005 
1006  B[2][0] = C[15] * Ai[0] + C[16] * Ai[1] + C[17] * Ai[3];
1007  B[2][1] = C[15] * Ai[1] + C[16] * Ai[2] + C[17] * Ai[4];
1008  B[2][2] = C[15] * Ai[3] + C[16] * Ai[4] + C[17] * Ai[5];
1009 
1010  B[3][0] = C[21] * Ai[0] + C[22] * Ai[1] + C[23] * Ai[3];
1011  B[3][1] = C[21] * Ai[1] + C[22] * Ai[2] + C[23] * Ai[4];
1012  B[3][2] = C[21] * Ai[3] + C[22] * Ai[4] + C[23] * Ai[5];
1013 
1014  B[4][0] = C[28] * Ai[0] + C[29] * Ai[1] + C[30] * Ai[3];
1015  B[4][1] = C[28] * Ai[1] + C[29] * Ai[2] + C[30] * Ai[4];
1016  B[4][2] = C[28] * Ai[3] + C[29] * Ai[4] + C[30] * Ai[5];
1017 
1018  double z[3] = {m[0] - r[0], m[1] - r[1], m[2] - r[2]};
1019 
1020  {
1021  double AV[6] = {C[0] - V[0],
1022  C[1] - V[1],
1023  C[2] - V[2],
1024  C[3] - V[3],
1025  C[4] - V[4],
1026  C[5] - V[5]};
1027  double AVi[6];
1028  AVi[0] = AV[4] * AV[4] - AV[2] * AV[5];
1029  AVi[1] = AV[1] * AV[5] - AV[3] * AV[4];
1030  AVi[2] = AV[3] * AV[3] - AV[0] * AV[5];
1031  AVi[3] = AV[2] * AV[3] - AV[1] * AV[4];
1032  AVi[4] = AV[0] * AV[4] - AV[1] * AV[3];
1033  AVi[5] = AV[1] * AV[1] - AV[0] * AV[2];
1034 
1035  det = 1. / (AV[0] * AVi[0] + AV[1] * AVi[1] + AV[3] * AVi[3]);
1036 
1037  NDF += 2;
1038  Chi2 += (+(AVi[0] * z[0] + AVi[1] * z[1] + AVi[3] * z[2]) * z[0]
1039  + (AVi[1] * z[0] + AVi[2] * z[1] + AVi[4] * z[2]) * z[1]
1040  + (AVi[3] * z[0] + AVi[4] * z[1] + AVi[5] * z[2]) * z[2])
1041  * det;
1042  }
1043 
1044  r[0] = m[0];
1045  r[1] = m[1];
1046  r[2] = m[2];
1047  r[3] += B[0][0] * z[0] + B[0][1] * z[1] + B[0][2] * z[2];
1048  r[4] += B[1][0] * z[0] + B[1][1] * z[1] + B[1][2] * z[2];
1049  r[5] += B[2][0] * z[0] + B[2][1] * z[1] + B[2][2] * z[2];
1050  r[6] += B[3][0] * z[0] + B[3][1] * z[1] + B[3][2] * z[2];
1051  r[7] += B[4][0] * z[0] + B[4][1] * z[1] + B[4][2] * z[2];
1052 
1053  double d0, d1, d2;
1054 
1055  C[0] = V[0];
1056  C[1] = V[1];
1057  C[2] = V[2];
1058  C[3] = V[3];
1059  C[4] = V[4];
1060  C[5] = V[5];
1061 
1062  d0 = B[0][0] * V[0] + B[0][1] * V[1] + B[0][2] * V[3] - C[6];
1063  d1 = B[0][0] * V[1] + B[0][1] * V[2] + B[0][2] * V[4] - C[7];
1064  d2 = B[0][0] * V[3] + B[0][1] * V[4] + B[0][2] * V[5] - C[8];
1065 
1066  C[6] += d0;
1067  C[7] += d1;
1068  C[8] += d2;
1069  C[9] += d0 * B[0][0] + d1 * B[0][1] + d2 * B[0][2];
1070 
1071  d0 = B[1][0] * V[0] + B[1][1] * V[1] + B[1][2] * V[3] - C[10];
1072  d1 = B[1][0] * V[1] + B[1][1] * V[2] + B[1][2] * V[4] - C[11];
1073  d2 = B[1][0] * V[3] + B[1][1] * V[4] + B[1][2] * V[5] - C[12];
1074 
1075  C[10] += d0;
1076  C[11] += d1;
1077  C[12] += d2;
1078  C[13] += d0 * B[0][0] + d1 * B[0][1] + d2 * B[0][2];
1079  C[14] += d0 * B[1][0] + d1 * B[1][1] + d2 * B[1][2];
1080 
1081  d0 = B[2][0] * V[0] + B[2][1] * V[1] + B[2][2] * V[3] - C[15];
1082  d1 = B[2][0] * V[1] + B[2][1] * V[2] + B[2][2] * V[4] - C[16];
1083  d2 = B[2][0] * V[3] + B[2][1] * V[4] + B[2][2] * V[5] - C[17];
1084 
1085  C[15] += d0;
1086  C[16] += d1;
1087  C[17] += d2;
1088  C[18] += d0 * B[0][0] + d1 * B[0][1] + d2 * B[0][2];
1089  C[19] += d0 * B[1][0] + d1 * B[1][1] + d2 * B[1][2];
1090  C[20] += d0 * B[2][0] + d1 * B[2][1] + d2 * B[2][2];
1091 
1092  d0 = B[3][0] * V[0] + B[3][1] * V[1] + B[3][2] * V[3] - C[21];
1093  d1 = B[3][0] * V[1] + B[3][1] * V[2] + B[3][2] * V[4] - C[22];
1094  d2 = B[3][0] * V[3] + B[3][1] * V[4] + B[3][2] * V[5] - C[23];
1095 
1096  C[21] += d0;
1097  C[22] += d1;
1098  C[23] += d2;
1099  C[24] += d0 * B[0][0] + d1 * B[0][1] + d2 * B[0][2];
1100  C[25] += d0 * B[1][0] + d1 * B[1][1] + d2 * B[1][2];
1101  C[26] += d0 * B[2][0] + d1 * B[2][1] + d2 * B[2][2];
1102  C[27] += d0 * B[3][0] + d1 * B[3][1] + d2 * B[3][2];
1103 
1104  d0 = B[4][0] * V[0] + B[4][1] * V[1] + B[4][2] * V[3] - C[28];
1105  d1 = B[4][0] * V[1] + B[4][1] * V[2] + B[4][2] * V[4] - C[29];
1106  d2 = B[4][0] * V[3] + B[4][1] * V[4] + B[4][2] * V[5] - C[30];
1107 
1108  C[28] += d0;
1109  C[29] += d1;
1110  C[30] += d2;
1111  C[31] += d0 * B[0][0] + d1 * B[0][1] + d2 * B[0][2];
1112  C[32] += d0 * B[1][0] + d1 * B[1][1] + d2 * B[1][2];
1113  C[33] += d0 * B[2][0] + d1 * B[2][1] + d2 * B[2][2];
1114  C[34] += d0 * B[3][0] + d1 * B[3][1] + d2 * B[3][2];
1115  C[35] += d0 * B[4][0] + d1 * B[4][1] + d2 * B[4][2];
1116 
1117  Extrapolate(r0, r[7]);
1118  Convert(r0, 0);
1119 }
1120 
1121 
1122 void CbmKFParticle::Convert(Double_t r0[], bool ToProduction) {
1123 
1124  Double_t B[3];
1125  {
1126  FairField* MF = CbmKF::Instance()->GetMagneticField();
1127  MF->GetFieldValue(r0, B);
1128  const Double_t c_light = Q * 0.000299792458;
1129  B[0] *= c_light;
1130  B[1] *= c_light;
1131  B[2] *= c_light;
1132  }
1133 
1134  Double_t h[6];
1135 
1136  h[0] = r0[3];
1137  h[1] = r0[4];
1138  h[2] = r0[5];
1139  if (ToProduction) {
1140  h[0] = -h[0];
1141  h[1] = -h[1];
1142  h[2] = -h[2];
1143  }
1144  h[3] = h[1] * B[2] - h[2] * B[1];
1145  h[4] = h[2] * B[0] - h[0] * B[2];
1146  h[5] = h[0] * B[1] - h[1] * B[0];
1147 
1148  Double_t c;
1149 
1150  c = C[28] + h[0] * C[35];
1151  C[0] += h[0] * (c + C[28]);
1152  C[28] = c;
1153 
1154  C[1] += h[1] * C[28] + h[0] * C[29];
1155  c = C[29] + h[1] * C[35];
1156  C[2] += h[1] * (c + C[29]);
1157  C[29] = c;
1158 
1159  C[3] += h[2] * C[28] + h[0] * C[30];
1160  C[4] += h[2] * C[29] + h[1] * C[30];
1161  c = C[30] + h[2] * C[35];
1162  C[5] += h[2] * (c + C[30]);
1163  C[30] = c;
1164 
1165  C[6] += h[3] * C[28] + h[0] * C[31];
1166  C[7] += h[3] * C[29] + h[1] * C[31];
1167  C[8] += h[3] * C[30] + h[2] * C[31];
1168  c = C[31] + h[3] * C[35];
1169  C[9] += h[3] * (c + C[31]);
1170  C[31] = c;
1171 
1172  C[10] += h[4] * C[28] + h[0] * C[32];
1173  C[11] += h[4] * C[29] + h[1] * C[32];
1174  C[12] += h[4] * C[30] + h[2] * C[32];
1175  C[13] += h[4] * C[31] + h[3] * C[32];
1176  c = C[32] + h[4] * C[35];
1177  C[14] += h[4] * (c + C[32]);
1178  C[32] = c;
1179 
1180  C[15] += h[5] * C[28] + h[0] * C[33];
1181  C[16] += h[5] * C[29] + h[1] * C[33];
1182  C[17] += h[5] * C[30] + h[2] * C[33];
1183  C[18] += h[5] * C[31] + h[3] * C[33];
1184  C[19] += h[5] * C[32] + h[4] * C[33];
1185  c = C[33] + h[5] * C[35];
1186  C[20] += h[5] * (c + C[33]);
1187  C[33] = c;
1188 
1189  C[21] += h[0] * C[34];
1190  C[22] += h[1] * C[34];
1191  C[23] += h[2] * C[34];
1192  C[24] += h[3] * C[34];
1193  C[25] += h[4] * C[34];
1194  C[26] += h[5] * C[34];
1195 }
1196 
1197 void CbmKFParticle::ExtrapolateLine(Double_t r0[], Double_t dS) {
1198  if (r0 && r0 != r) {
1199  r0[0] += dS * r0[3];
1200  r0[1] += dS * r0[4];
1201  r0[2] += dS * r0[5];
1202  }
1203 
1204  r[0] += dS * r[3];
1205  r[1] += dS * r[4];
1206  r[2] += dS * r[5];
1207 
1208  double C6 = C[6] + dS * C[9];
1209  double C11 = C[11] + dS * C[14];
1210  double C17 = C[17] + dS * C[20];
1211  double SC13 = dS * C[13];
1212  double SC18 = dS * C[18];
1213  double SC19 = dS * C[19];
1214 
1215  C[0] += dS * (C[6] + C6);
1216  C[2] += dS * (C[11] + C11);
1217  C[5] += dS * (C[17] + C17);
1218 
1219  C[7] += SC13;
1220  C[8] += SC18;
1221  C[12] += SC19;
1222 
1223  C[1] += dS * (C[10] + C[7]);
1224  C[3] += dS * (C[15] + C[8]);
1225  C[4] += dS * (C[16] + C[12]);
1226 
1227  C[6] = C6;
1228  C[10] += SC13;
1229  C[11] = C11;
1230  C[15] += SC18;
1231  C[16] += SC19;
1232  C[17] = C17;
1233 
1234  C[21] += dS * C[24];
1235  C[22] += dS * C[25];
1236  C[23] += dS * C[26];
1237  C[28] += dS * C[31];
1238  C[29] += dS * C[32];
1239  C[30] += dS * C[33];
1240 }
1241 
1242 void CbmKFParticle::Extrapolate(Double_t r0[], Double_t dS) {
1243 
1244  if (Q == 0) {
1245  ExtrapolateLine(r0, dS);
1246  return;
1247  }
1248 
1249  const Double_t c_light = 0.000299792458;
1250  FairField* MF = CbmKF::Instance()->GetMagneticField();
1251 
1252  Double_t c = Q * c_light;
1253 
1254  // construct coefficients
1255 
1256  Double_t px = r[3], py = r[4], pz = r[5];
1257 
1258  Double_t sx = 0, sy = 0, sz = 0, syy = 0, syz = 0, syyy = 0, Sx = 0, Sy = 0,
1259  Sz = 0, Syy = 0, Syz = 0, Syyy = 0;
1260 
1261  { // get field integrals
1262 
1263  Double_t B[3][3];
1264  Double_t p0[3], p1[3], p2[3];
1265 
1266  // line track approximation
1267 
1268  p0[0] = r[0];
1269  p0[1] = r[1];
1270  p0[2] = r[2];
1271 
1272  p2[0] = r[0] + px * dS;
1273  p2[1] = r[1] + py * dS;
1274  p2[2] = r[2] + pz * dS;
1275 
1276  p1[0] = 0.5 * (p0[0] + p2[0]);
1277  p1[1] = 0.5 * (p0[1] + p2[1]);
1278  p1[2] = 0.5 * (p0[2] + p2[2]);
1279 
1280  // first order track approximation
1281  {
1282  MF->GetFieldValue(p0, B[0]);
1283  MF->GetFieldValue(p1, B[1]);
1284  MF->GetFieldValue(p2, B[2]);
1285 
1286  Double_t Sy1 = (7 * B[0][1] + 6 * B[1][1] - B[2][1]) * c * dS * dS / 96.;
1287  Double_t Sy2 = (B[0][1] + 2 * B[1][1]) * c * dS * dS / 6.;
1288 
1289  p1[0] -= Sy1 * pz;
1290  p1[2] += Sy1 * px;
1291  p2[0] -= Sy2 * pz;
1292  p2[2] += Sy2 * px;
1293  }
1294 
1295  MF->GetFieldValue(p0, B[0]);
1296  MF->GetFieldValue(p1, B[1]);
1297  MF->GetFieldValue(p2, B[2]);
1298 
1299  sx = c * (B[0][0] + 4 * B[1][0] + B[2][0]) * dS / 6.;
1300  sy = c * (B[0][1] + 4 * B[1][1] + B[2][1]) * dS / 6.;
1301  sz = c * (B[0][2] + 4 * B[1][2] + B[2][2]) * dS / 6.;
1302 
1303  Sx = c * (B[0][0] + 2 * B[1][0]) * dS * dS / 6.;
1304  Sy = c * (B[0][1] + 2 * B[1][1]) * dS * dS / 6.;
1305  Sz = c * (B[0][2] + 2 * B[1][2]) * dS * dS / 6.;
1306 
1307  Double_t c2[3][3] = {{5, -4, -1}, {44, 80, -4}, {11, 44, 5}}; // /=360.
1308  Double_t C2[3][3] = {{38, 8, -4}, {148, 208, -20}, {3, 36, 3}}; // /=2520.
1309  for (Int_t n = 0; n < 3; n++)
1310  for (Int_t m = 0; m < 3; m++) {
1311  syz += c2[n][m] * B[n][1] * B[m][2];
1312  Syz += C2[n][m] * B[n][1] * B[m][2];
1313  }
1314 
1315  syz *= c * c * dS * dS / 360.;
1316  Syz *= c * c * dS * dS * dS / 2520.;
1317 
1318  syy = c * (B[0][1] + 4 * B[1][1] + B[2][1]) * dS;
1319  syyy = syy * syy * syy / 1296;
1320  syy = syy * syy / 72;
1321 
1322  Syy = (B[0][1] * (38 * B[0][1] + 156 * B[1][1] - B[2][1])
1323  + B[1][1] * (208 * B[1][1] + 16 * B[2][1]) + B[2][1] * (3 * B[2][1]))
1324  * dS * dS * dS * c * c / 2520.;
1325  Syyy = (B[0][1]
1326  * (B[0][1] * (85 * B[0][1] + 526 * B[1][1] - 7 * B[2][1])
1327  + B[1][1] * (1376 * B[1][1] + 84 * B[2][1])
1328  + B[2][1] * (19 * B[2][1]))
1329  + B[1][1]
1330  * (B[1][1] * (1376 * B[1][1] + 256 * B[2][1])
1331  + B[2][1] * (62 * B[2][1]))
1332  + B[2][1] * B[2][1] * (3 * B[2][1]))
1333  * dS * dS * dS * dS * c * c * c / 90720.;
1334  }
1335 
1336  Double_t J[8][8];
1337  for (int i = 0; i < 8; i++)
1338  for (int j = 0; j < 8; j++)
1339  J[i][j] = 0;
1340 
1341  J[0][0] = 1;
1342  J[0][1] = 0;
1343  J[0][2] = 0;
1344  J[0][3] = dS - Syy;
1345  J[0][4] = Sx;
1346  J[0][5] = Syyy - Sy;
1347  J[1][0] = 0;
1348  J[1][1] = 1;
1349  J[1][2] = 0;
1350  J[1][3] = -Sz;
1351  J[1][4] = dS;
1352  J[1][5] = Sx + Syz;
1353  J[2][0] = 0;
1354  J[2][1] = 0;
1355  J[2][2] = 1;
1356  J[2][3] = Sy - Syyy;
1357  J[2][4] = -Sx;
1358  J[2][5] = dS - Syy;
1359 
1360  J[3][0] = 0;
1361  J[3][1] = 0;
1362  J[3][2] = 0;
1363  J[3][3] = 1 - syy;
1364  J[3][4] = sx;
1365  J[3][5] = syyy - sy;
1366  J[4][0] = 0;
1367  J[4][1] = 0;
1368  J[4][2] = 0;
1369  J[4][3] = -sz;
1370  J[4][4] = 1;
1371  J[4][5] = sx + syz;
1372  J[5][0] = 0;
1373  J[5][1] = 0;
1374  J[5][2] = 0;
1375  J[5][3] = sy - syyy;
1376  J[5][4] = -sx;
1377  J[5][5] = 1 - syy;
1378  J[6][6] = J[7][7] = 1;
1379 
1380  r[0] += J[0][3] * px + J[0][4] * py + J[0][5] * pz;
1381  r[1] += J[1][3] * px + J[1][4] * py + J[1][5] * pz;
1382  r[2] += J[2][3] * px + J[2][4] * py + J[2][5] * pz;
1383  r[3] = J[3][3] * px + J[3][4] * py + J[3][5] * pz;
1384  r[4] = J[4][3] * px + J[4][4] * py + J[4][5] * pz;
1385  r[5] = J[5][3] * px + J[5][4] * py + J[5][5] * pz;
1386 
1387  if (r0 && r0 != r) {
1388  px = r0[3];
1389  py = r0[4];
1390  pz = r0[5];
1391 
1392  r0[0] += J[0][3] * px + J[0][4] * py + J[0][5] * pz;
1393  r0[1] += J[1][3] * px + J[1][4] * py + J[1][5] * pz;
1394  r0[2] += J[2][3] * px + J[2][4] * py + J[2][5] * pz;
1395  r0[3] = J[3][3] * px + J[3][4] * py + J[3][5] * pz;
1396  r0[4] = J[4][3] * px + J[4][4] * py + J[4][5] * pz;
1397  r0[5] = J[5][3] * px + J[5][4] * py + J[5][5] * pz;
1398  }
1399 
1400  CbmKFMath::multQSQt(8, J[0], C, C);
1401 }
1402 
1403 
1404 void CbmKFParticle::TransportToProductionVertex() {
1405  if (AtProductionVertex) return;
1406  Double_t r0[8];
1407  for (int i = 0; i < 8; i++)
1408  r0[i] = r[i];
1409  Extrapolate(r0, -r[7]);
1410  Convert(r0, 1);
1411  AtProductionVertex = 1;
1412 }
1413 
1414 void CbmKFParticle::TransportToDecayVertex() {
1415  if (!AtProductionVertex) return;
1416  Double_t r0[8];
1417  for (int i = 0; i < 8; i++)
1418  r0[i] = r[i];
1419  Extrapolate(r0, r[7]);
1420  Convert(r0, 0);
1421  AtProductionVertex = 0;
1422 }
1423 
1424 
1425 void CbmKFParticle::GetKFTrack(CbmKFTrackInterface* Track) {
1426 
1427  Double_t* T = Track->GetTrack();
1428  Double_t* Cov = Track->GetCovMatrix();
1429  Track->GetRefNDF() = NDF;
1430  Track->GetRefChi2() = Chi2;
1431 
1432  Double_t px = r[3], py = r[4], pz = r[5];
1433  Double_t p = sqrt(px * px + py * py + pz * pz);
1434  Double_t pzi = 1 / pz;
1435  Double_t qp = 1 / p;
1436 
1437  if (Q != 0) qp = Q * qp;
1438 
1439  T[0] = r[0]; // dX = dx -T[2]*dz
1440  T[1] = r[1]; // dY = dy -T[3]*dz
1441  T[2] = px * pzi; // dtx= (dpx - T[2]*dpz)/pz
1442  T[3] = py * pzi; // dty= (dpy - T[3]*dpz)/pz
1443  T[4] = qp; // dq = (px*dpx + py*dpy + pz*dpz)*(-q/p^3)
1444  T[5] = r[2];
1445 
1446  Double_t qp3 = qp * qp * qp;
1447  Double_t pzi2 = pzi * pzi;
1448 
1449  Double_t cXpz = C[15] - T[2] * C[17];
1450  Double_t cYpz = C[16] - T[3] * C[17];
1451  Double_t cqpx = -qp3 * (px * C[9] + py * C[13] + pz * C[18]);
1452  Double_t cqpy = -qp3 * (px * C[13] + py * C[14] + pz * C[19]);
1453  Double_t cqpz = -qp3 * (px * C[18] + py * C[19] + pz * C[20]);
1454 
1455  Cov[0] = C[0] + T[2] * (-2 * C[3] + T[2] * C[5]);
1456  Cov[1] = C[1] - T[2] * C[4] + T[3] * (-C[3] + T[2] * C[5]);
1457  Cov[2] = C[2] + T[3] * (-2 * C[4] + T[3] * C[5]);
1458  Cov[3] = pzi * (C[6] - T[2] * (C[8] + cXpz));
1459  Cov[4] = pzi * (C[7] - T[3] * C[8] - T[2] * cYpz);
1460  Cov[5] = pzi2 * (C[9] + T[2] * (-2 * C[18] + T[2] * C[20]));
1461  Cov[6] = pzi * (C[10] - T[2] * C[12] - T[3] * cXpz);
1462  Cov[7] = pzi * (C[11] - T[3] * (C[12] + cYpz));
1463  Cov[8] = pzi2 * (C[13] - T[2] * C[19] - T[3] * C[18] + T[2] * T[3] * C[20]);
1464  Cov[9] = pzi2 * (C[14] + T[3] * (-2 * C[19] + T[3] * C[20]));
1465  Cov[10] = -qp3 * (px * C[6] + py * C[10] + pz * C[15]) - T[2] * cqpz;
1466  Cov[11] = -qp3 * (px * C[7] + py * C[11] + pz * C[16]) - T[3] * cqpz;
1467  Cov[12] = pzi * (cqpx - T[2] * cqpz);
1468  Cov[13] = pzi * (cqpy - T[3] * cqpz);
1469  Cov[14] = -qp3 * (px * cqpx + py * cqpy + pz * cqpz);
1470 }
1471 
1472 void CbmKFParticle::GetKFVertex(CbmKFVertexInterface& vtx) {
1473  vtx.GetRefX() = r[0];
1474  vtx.GetRefY() = r[1];
1475  vtx.GetRefZ() = r[2];
1476  for (int i = 0; i < 6; i++)
1477  vtx.GetCovMatrix()[i] = C[i];
1478  vtx.GetRefChi2() = Chi2;
1479  vtx.GetRefNDF() = NDF;
1480 }
1481 
1482 void CbmKFParticle::GetMomentum(Double_t& P, Double_t& Error_) {
1483  Double_t x = r[3];
1484  Double_t y = r[4];
1485  Double_t z = r[5];
1486  Double_t x2 = x * x;
1487  Double_t y2 = y * y;
1488  Double_t z2 = z * z;
1489  Double_t p2 = x2 + y2 + z2;
1490  P = sqrt(p2);
1491  Error_ = sqrt((x2 * C[9] + y2 * C[14] + z2 * C[20]
1492  + 2 * (x * y * C[13] + x * z * C[18] + y * z * C[19]))
1493  / p2);
1494 }
1495 
1496 void CbmKFParticle::GetMass(Double_t& M, Double_t& Error_) {
1497  // S = sigma^2 of m2/2
1498 
1499  Double_t S =
1500  (r[3] * r[3] * C[9] + r[4] * r[4] * C[14] + r[5] * r[5] * C[20]
1501  + r[6] * r[6] * C[27]
1502  + 2
1503  * (+r[3] * r[4] * C[13] + r[5] * (r[3] * C[18] + r[4] * C[19])
1504  - r[6] * (r[3] * C[24] + r[4] * C[25] + r[5] * C[26])));
1505  Double_t m2 = r[6] * r[6] - r[3] * r[3] - r[4] * r[4] - r[5] * r[5];
1506  M = (m2 > 1.e-20) ? sqrt(m2) : 0.;
1507  Error_ = (S >= 0 && m2 > 1.e-20) ? sqrt(S / m2) : 1.e4;
1508 }
1509 
1510 
1511 void CbmKFParticle::GetDecayLength(Double_t& L, Double_t& Error_) {
1512  Double_t x = r[3];
1513  Double_t y = r[4];
1514  Double_t z = r[5];
1515  Double_t t = r[7];
1516  Double_t x2 = x * x;
1517  Double_t y2 = y * y;
1518  Double_t z2 = z * z;
1519  Double_t p2 = x2 + y2 + z2;
1520  L = t * sqrt(p2);
1521  Error_ = sqrt(p2 * C[35]
1522  + t * t / p2
1523  * (x2 * C[9] + y2 * C[14] + z2 * C[20]
1524  + 2 * (x * y * C[13] + x * z * C[18] + y * z * C[19]))
1525  + 2 * t * (x * C[31] + y * C[32] + z * C[33]));
1526 }
1527 
1528 void CbmKFParticle::GetLifeTime(Double_t& TauC, Double_t& Error_) {
1529  Double_t m, dm;
1530  GetMass(m, dm);
1531  Double_t cTM = (-r[3] * C[31] - r[4] * C[32] - r[5] * C[33] + r[6] * C[34]);
1532  TauC = r[7] * m;
1533  Error_ = sqrt(m * m * C[35] + 2 * r[7] * cTM + r[7] * r[7] * dm * dm);
1534 }
1535 
1536 Double_t CbmKFParticle::GetRapidity() const {
1537  return 0.5 * TMath::Log((r[6] + r[5]) / (r[6] - r[5]));
1538 }
1539 Double_t CbmKFParticle::GetPt() const {
1540  return TMath::Sqrt(r[3] * r[3] + r[4] * r[4]);
1541 }
1542 Double_t CbmKFParticle::GetTheta() const { return TMath::ATan2(GetPt(), r[5]); }
1543 Double_t CbmKFParticle::GetPhi() const { return TMath::ATan2(r[4], r[5]); }
1544 
1545 
1546 Double_t CbmKFParticle::GetDStoPoint(const Double_t xyz[]) const {
1547  //TODO check the method!
1548  Double_t dx = xyz[0] - r[0];
1549  Double_t dy = xyz[1] - r[1];
1550  Double_t dz = xyz[2] - r[2];
1551  Double_t p2 = r[3] * r[3] + r[4] * r[4] + r[5] * r[5];
1552 
1553  if (Q == 0 && r[2] < xyz[2]) return sqrt((dx * dx + dy * dy + dz * dz) / p2);
1554  if (Q == 0 && r[2] >= xyz[2])
1555  return -sqrt((dx * dx + dy * dy + dz * dz) / p2);
1556 
1557  const Double_t kCLight = 0.000299792458;
1558  Double_t B[3];
1559  FairField* MF = CbmKF::Instance()->GetMagneticField();
1560  MF->GetFieldValue(r, B);
1561 
1562  Double_t bq1 = B[0] * Q * kCLight;
1563  Double_t bq2 = B[1] * Q * kCLight;
1564  Double_t bq3 = B[2] * Q * kCLight;
1565  Double_t a = dx * r[3] + dy * r[4] + dz * r[5];
1566  Double_t B2 = B[0] * B[0] + B[1] * B[1] + B[2] * B[2];
1567  Double_t bq = Q * kCLight * sqrt(B2);
1568  Double_t pB = r[3] * B[0] + r[4] * B[1] + r[5] * B[2];
1569  Double_t rB = dx * B[0] + dy * B[1] + dz * B[2];
1570  Double_t pt2 = p2 - pB * pB / B2;
1571  a = a - pB * rB / B2;
1572 
1573  Double_t dS;
1574 
1575  if (pt2 < 1.e-4) return 0;
1576 
1577  if (TMath::Abs(bq) < 1.e-8)
1578  dS = a / pt2;
1579  else
1580  dS = TMath::ATan2(bq * a,
1581  pt2 + bq1 * (dz * r[4] - dy * r[5])
1582  - bq2 * (dz * r[3] - dx * r[5])
1583  + bq3 * (dy * r[3] - dx * r[4]))
1584  / bq;
1585 
1586  // Double_t dSm = rB/pB;
1587  // cout <<"normalnyj S " << dS << endl;
1588  // cout <<"moj S " << dSm << endl;
1589 
1590  return dS;
1591  //return sqrt(dx*dx+dy*dy+dz*dz);
1592 }
CbmKFParticle::NDF
Int_t NDF
Definition: CbmKFParticle.h:39
CbmKFParticle::GetLifeTime
void GetLifeTime(Double_t &T, Double_t &Error)
Definition: CbmKFParticle.cxx:1528
CbmKFParticle::CbmKFParticle
CbmKFParticle()
Definition: CbmKFParticle.h:45
CbmKFParticle::GetDecayLength
void GetDecayLength(Double_t &L, Double_t &Error)
Definition: CbmKFParticle.cxx:1511
CbmKFParticle::GetKFTrack
void GetKFTrack(CbmKFTrackInterface *Track)
Definition: CbmKFParticle.cxx:1425
CbmKFTube::z
Double_t z
Definition: CbmKFMaterial.h:92
h
Generates beam ions for transport simulation.
Definition: CbmBeamGenerator.h:17
CbmKFParticle::Convert
void Convert(Double_t r0[], bool ToProduction)
Definition: CbmKFParticle.cxx:1122
CbmKFParticle::Extrapolate
void Extrapolate(Double_t r0[], double T)
Definition: CbmKFParticle.cxx:1242
CbmKF.h
CbmKFParticle::TransportToDecayVertex
void TransportToDecayVertex()
Definition: CbmKFParticle.cxx:1414
CbmKFParticle::MeasureProductionVertex
void MeasureProductionVertex(Double_t r0[], CbmKFVertexInterface *Parent)
Definition: CbmKFParticle.cxx:970
CbmKFParticle::r
Double_t r[8]
Definition: CbmKFParticle.h:38
CbmKFTrackInterface::GetTrack
virtual Double_t * GetTrack()
Is it electron.
Definition: CbmKFTrackInterface.cxx:33
CbmKF::GetMagneticField
FairField * GetMagneticField()
Definition: CbmKF.h:88
CbmKFParticle::GetPhi
Double_t GetPhi() const
Definition: CbmKFParticle.cxx:1543
sqrt
friend F32vec4 sqrt(const F32vec4 &a)
Definition: L1/vectors/P4_F32vec4.h:41
CbmKFParticle::Construct
void Construct(std::vector< CbmKFTrackInterface * > &vDaughters, CbmKFVertexInterface *Parent=0, Double_t Mass=-1, Double_t CutChi2=-1)
Definition: CbmKFParticle.cxx:144
CbmKFTrack::GetTrack
Double_t * GetTrack()
Is it electron.
Definition: CbmKFTrack.h:58
i
int i
Definition: L1/vectors/P4_F32vec4.h:25
CbmKFParticle::ConstructFromKFParticle
void ConstructFromKFParticle(std::vector< CbmKFParticle * > &vDaughters, CbmKFVertexInterface *Parent=0, Double_t Mass=-1, Double_t CutChi2=-1)
Definition: CbmKFParticle.cxx:533
CbmKF::vTargets
std::vector< CbmKFTube > vTargets
Definition: CbmKF.h:79
CbmKFParticle::GetKFVertex
void GetKFVertex(CbmKFVertexInterface &vtx)
Definition: CbmKFParticle.cxx:1472
CbmKFTrack::GetMass
Double_t GetMass()
Definition: CbmKFTrack.h:62
CbmKFParticle::GetMass
void GetMass(Double_t &M, Double_t &Error)
Definition: CbmKFParticle.cxx:1496
CbmKFTube::R
Double_t R
Definition: CbmKFMaterial.h:93
CbmKFParticle::GetPt
Double_t GetPt() const
Definition: CbmKFParticle.cxx:1539
CbmKFVertexInterface::GetRefChi2
virtual Double_t & GetRefChi2()
Array[6] of covariance matrix.
Definition: CbmKFVertexInterface.cxx:32
CbmKFParticle::GetDStoPoint
Double_t GetDStoPoint(const Double_t xyz[]) const
Definition: CbmKFParticle.cxx:1546
CbmKFParticle::TransportToProductionVertex
void TransportToProductionVertex()
Definition: CbmKFParticle.cxx:1404
CbmKFMath.h
CbmKFTrack::GetRefChi2
Double_t & GetRefChi2()
array[15] of covariance matrix
Definition: CbmKFTrack.h:60
CbmKF::Instance
static CbmKF * Instance()
Definition: CbmKF.h:39
CbmKFTrack::GetRefNDF
Int_t & GetRefNDF()
Chi^2 after fit.
Definition: CbmKFTrack.h:61
z2
Double_t z2[nSects2]
Definition: pipe_v16a_mvdsts100.h:11
CbmKFMath::multQSQt
static void multQSQt(Int_t N, const Double_t Q[], const Double_t S[], Double_t S_out[])
Definition: CbmKFMath.cxx:42
CbmKFVertexInterface::GetRefZ
virtual Double_t & GetRefZ()
Definition: CbmKFVertexInterface.cxx:30
CbmKFParticle::ExtrapolateLine
void ExtrapolateLine(Double_t r0[], double T)
Definition: CbmKFParticle.cxx:1197
d
double d
Definition: P4_F64vec2.h:24
CbmKFTrackInterface::Id
int Id() const
Definition: CbmKFTrackInterface.h:67
drawcov
void drawcov(double C[])
Definition: CbmKFParticle.cxx:129
CbmKFParticle::Cij
Double_t & Cij(Int_t i, Int_t j)
Definition: CbmKFParticle.h:145
CbmKFTrack::GetCovMatrix
Double_t * GetCovMatrix()
array[6] of track parameters(x,y,tx,ty,qp,z)
Definition: CbmKFTrack.h:59
CbmKFVertexInterface
Definition: CbmKFVertexInterface.h:24
CbmKFTube::dz
Double_t dz
Definition: CbmKFMaterial.h:92
ClassImp
ClassImp(CbmKFParticle) CbmKFParticle
Definition: CbmKFParticle.cxx:25
CbmKFTrack.h
CbmKFTrackInterface::GetRefNDF
virtual Int_t & GetRefNDF()
Chi^2 after fit.
Definition: CbmKFTrackInterface.cxx:36
CbmKFParticle::AtProductionVertex
Bool_t AtProductionVertex
Definition: CbmKFParticle.h:41
CbmKFVertexInterface::GetRefX
virtual Double_t & GetRefX()
Definition: CbmKFVertexInterface.cxx:28
first
bool first
Definition: LKFMinuit.cxx:143
CbmKFParticle::Chi2
Double_t Chi2
Definition: CbmKFParticle.h:40
CbmKFParticle.h
CbmKFTrackInterface
Definition: CbmKFTrackInterface.h:26
CbmKFParticle::Q
Double_t Q
Definition: CbmKFParticle.h:40
CbmKFParticle::C
Double_t C[36]
Definition: CbmKFParticle.h:38
CbmKFVertexInterface::GetRefY
virtual Double_t & GetRefY()
Definition: CbmKFVertexInterface.cxx:29
CbmKFTrackInterface::GetCovMatrix
virtual Double_t * GetCovMatrix()
array[6] of track parameters(x,y,tx,ty,qp,z)
Definition: CbmKFTrackInterface.cxx:34
CbmKFParticle
Definition: CbmKFParticle.h:18
CbmKFParticle::GetRapidity
Double_t GetRapidity() const
Definition: CbmKFParticle.cxx:1536
CbmKFTube
Definition: CbmKFMaterial.h:77
x
Double_t x
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmKFParticle::GetMomentum
void GetMomentum(Double_t &P, Double_t &Error)
Definition: CbmKFParticle.cxx:1482
m
__m128 m
Definition: L1/vectors/P4_F32vec4.h:26
y
Double_t y
Definition: CbmMvdSensorDigiToHitTask.cxx:68
CbmKFParticle::MeasureMass
void MeasureMass(Double_t r0[], Double_t Mass)
Definition: CbmKFParticle.cxx:934
lit::parallel::NDF
unsigned short NDF(const LitScalTrack &track)
Returns number of degrees of freedom for the track.
Definition: LitMath.h:50
CbmKFParticle::GetTheta
Double_t GetTheta() const
Definition: CbmKFParticle.cxx:1542
CbmKFTrack
Definition: CbmKFTrack.h:21
CbmKFVertexInterface::GetRefNDF
virtual Int_t & GetRefNDF()
Chi^2 after fit.
Definition: CbmKFVertexInterface.cxx:33
CbmKFTrackInterface::GetRefChi2
virtual Double_t & GetRefChi2()
array[15] of covariance matrix
Definition: CbmKFTrackInterface.cxx:35
CbmKFVertexInterface::GetCovMatrix
virtual Double_t * GetCovMatrix()
Definition: CbmKFVertexInterface.cxx:31
CbmStsKFTrackFitter.h
CbmKFTrackInterface::Extrapolate
Int_t Extrapolate(Double_t z, Double_t *QP0=0)
Access to i-th hit.
Definition: CbmKFTrackInterface.cxx:39
NS_L1TrackFitter::c_light
const fvec c_light
Definition: L1TrackFitter.cxx:31