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