EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
FairTrackParP.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file FairTrackParP.cxx
1 // Class for the representation of a track as parabola (SD system)
2 //
3 // Authors: M. Al-Turany, A. Fontana, L. Lavezzi and A. Rotondi
4 //
5 //
6 // GEANE parameters (q/p, v', w', v, w) of Helix track
7 // The Helix can be constructed using the Helix parameter (1/p, v', w', v, w) in SC reference
8 // and the covariance matrix. Or using position and momentum in LAB referance.
9 
10 #include "FairTrackParP.h"
11 #include "FairGeaneUtil.h"
12 #include "FairRunAna.h"
13 #include "FairField.h"
14 #include <cmath>
15 #include "TMath.h"
16 
17 #include <iostream>
18 
19 using std::cout;
20 using std::endl;
21 
22 ClassImp(FairTrackParP)
23 
24 // ----- Default constructor -------------------------------------------
25 FairTrackParP::FairTrackParP()
26  : FairTrackPar(),
27  fU (0.),
28  fV (0.),
29  fW (0.),
30  fTV(0.),
31  fTW(0.),
32  fPx_sd(0.),
33  fPy_sd(0.),
34  fPz_sd(0.),
35  fDU(0.),
36  fDV(0.),
37  fDW(0.),
38  fDTV(0.),
39  fDTW(0.),
40  forigin(TVector3(0,0,0)),
41  fiver(TVector3(0,0,0)),
42  fjver(TVector3(0,0,0)),
43  fkver(TVector3(0,0,0)),
44  fSPU(0.)
45 {
46 // Reset();
47  for(Int_t i=0; i<15; i++) {
48  fCovMatrix[i]=0;
49  }
50  for(int i = 0; i < 6; i++) for(int j = 0; j < 6; j++) { fCovMatrix66[i][j] = 0; }
51 
52 
53 
54 }
55 
56 // ----- Constructor with track parameters in SD -----------------------------------
57 FairTrackParP::FairTrackParP(Double_t v, Double_t w, Double_t Tv,
58  Double_t Tw, Double_t qp,
59  Double_t CovMatrix[15],
60  TVector3 o, TVector3 dj, TVector3 dk)
61  : FairTrackPar(),
62  fU (0.),
63  fV (v),
64  fW (w),
65  fTV(Tv),
66  fTW(Tw),
67  fPx_sd(0.),
68  fPy_sd(0.),
69  fPz_sd(0.),
70  fDU(0.),
71  fDV(0.),
72  fDW(0.),
73  fDTV(0.),
74  fDTW(0.),
75  forigin(TVector3(0,0,0)),
76  fiver(TVector3(0,0,0)),
77  fjver(TVector3(0,0,0)),
78  fkver(TVector3(0,0,0)),
79  fSPU(0.)
80 {
81 // Reset();
82  SetPlane(o, dj, dk);
83 // fV = v;
84 // fW = w;
85 // fTV = Tv;
86 // fTW = Tw;
87  fQp = qp;
88 
89  Double_t P = TMath::Abs(1/fQp);
90  //fq= int (P * fQp);
91  fq = (int)TMath::Sign(1.0, fQp);
92  for(Int_t i=0; i<15; i++) {
93  fCovMatrix[i]=CovMatrix[i];
94  }
95 
96 
97  fPx_sd = TMath::Sqrt( P*P / ( fTV*fTV + fTW*fTW + 1 ) );
98  fPy_sd = fTV * fPx_sd;
99  fPz_sd = fTW * fPx_sd;
100 
101 
102 // fU = 0.;
103 
104  FairGeaneUtil util;
105  TVector3 position = util.FromSDToMARSCoord(TVector3(fU, fV, fW), forigin, fiver, fjver, fkver);
106  fX = position.X();
107  fY = position.Y();
108  fZ = position.Z();
109 
110  fDQp = TMath::Sqrt(fabs(fCovMatrix[0]));
111  fDTV = TMath::Sqrt(fabs(fCovMatrix[5]));
112  fDTW = TMath::Sqrt(fabs(fCovMatrix[9]));
113  fDV = TMath::Sqrt(fabs(fCovMatrix[12]));
114  fDW = TMath::Sqrt(fabs(fCovMatrix[14]));
115 
116 
117  Double_t PD[3],RD[6][6],H[3],PC[3],RC[15], SP1;
118  Int_t CH;
119  PC[0] = fQp;
120  PC[1] = fTV;
121  PC[2] = fTW;
122 
123  for(Int_t i=0; i<15; i++) {
124  RC[i]=fCovMatrix[i];
125  }
126 
127  // retrieve field
128  Double_t pnt[3];
129  pnt[0] = fX;
130  pnt[1] = fY;
131  pnt[2] = fZ;
132  FairRunAna* fRun = FairRunAna::Instance();
133  fRun->GetField()->GetFieldValue(pnt, H);
134 
135  CH=fq;
136 
137  TVector3 momsd = TVector3(fPx_sd,fPy_sd,fPz_sd);
138  SP1 = (momsd.Dot(fjver.Cross(fkver)))/TMath::Abs(momsd.Dot(fjver.Cross(fkver)));
139  fSPU=SP1;
140 
141  util.FromSDToMars(PC, RC, H, CH, SP1, fDJ, fDK, PD, RD);
142 
143  fPx = PD[0];
144  fPy = PD[1];
145  fPz = PD[2];
146 
147  for(int i = 0; i < 6; i++) for(int j = 0; j < 6; j++) { fCovMatrix66[i][j] = RD[i][j]; }
148 
149  fDPx = TMath::Sqrt(fabs(RD[0][0]));
150  fDPy = TMath::Sqrt(fabs(RD[1][1]));
151  fDPz = TMath::Sqrt(fabs(RD[2][2]));
152  fDX = TMath::Sqrt(fabs(RD[3][3]));
153  fDY = TMath::Sqrt(fabs(RD[4][4]));
154  fDZ = TMath::Sqrt(fabs(RD[5][5]));
155 
156 }
157 
158 // ----- Constructor with track parameters in SD -----------------------------------
159 FairTrackParP::FairTrackParP(Double_t v, Double_t w, Double_t Tv,
160  Double_t Tw, Double_t qp,
161  Double_t CovMatrix[15],
162  TVector3 o, TVector3 dj, TVector3 dk, Double_t spu)
163  : FairTrackPar(),
164  fU (0.),
165  fV (v),
166  fW (w),
167  fTV(Tv),
168  fTW(Tw),
169  fPx_sd(0.),
170  fPy_sd(0.),
171  fPz_sd(0.),
172  fDU(0.),
173  fDV(0.),
174  fDW(0.),
175  fDTV(0.),
176  fDTW(0.),
177  forigin(TVector3(0,0,0)),
178  fiver(TVector3(0,0,0)),
179  fjver(TVector3(0,0,0)),
180  fkver(TVector3(0,0,0)),
181  fSPU(spu)
182 
183 {
184 // Reset();
185  SetPlane(o, dj, dk);
186 // fV = v;
187 // fW = w;
188 // fTV = Tv;
189 // fTW = Tw;
190  fQp = qp;
191 // fSPU = spu;
192  Double_t P=0;
193  if(0!=fQp) {
194  P = TMath::Abs(1/fQp);
195  fq = (int)TMath::Sign(1.0, fQp);
196  //fq= int (P * fQp);
197  } else { fq=1; }
198  for(Int_t i=0; i<15; i++) {
199  fCovMatrix[i]=CovMatrix[i];
200  }
201 
202 
203  fPx_sd = fSPU*TMath::Sqrt( P*P / ( fTV*fTV + fTW*fTW + 1 ) );
204  fPy_sd = fTV * fPx_sd;
205  fPz_sd = fTW * fPx_sd;
206 
207 // fU = 0.;
208 
209  FairGeaneUtil util;
210  TVector3 position = util.FromSDToMARSCoord(TVector3(fU, fV, fW), forigin, fiver, fjver, fkver);
211  fX = position.X();
212  fY = position.Y();
213  fZ = position.Z();
214 
215  fDQp = TMath::Sqrt(fabs(fCovMatrix[0]));
216  fDTV = TMath::Sqrt(fabs(fCovMatrix[5]));
217  fDTW = TMath::Sqrt(fabs(fCovMatrix[9]));
218  fDV = TMath::Sqrt(fabs(fCovMatrix[12]));
219  fDW = TMath::Sqrt(fabs(fCovMatrix[14]));
220 
221 
222  Double_t PD[3],RD[6][6],H[3],PC[3],RC[15], SP1;
223  Int_t CH;
224  PC[0] = fQp;
225  PC[1] = fTV;
226  PC[2] = fTW;
227 
228  for(Int_t i=0; i<15; i++) {
229  RC[i]=fCovMatrix[i];
230  }
231 
232  // retrieve field
233  Double_t pnt[3];
234  pnt[0] = fX;
235  pnt[1] = fY;
236  pnt[2] = fZ;
237  FairRunAna* fRun = FairRunAna::Instance();
238  fRun->GetField()->GetFieldValue(pnt, H);
239 
240  CH=fq;
241 
242  TVector3 momsd = TVector3(fPx_sd,fPy_sd,fPz_sd);
243  // SP1 = (momsd.Dot(fjver.Cross(fkver)))/TMath::Abs(momsd.Dot(fjver.Cross(fkver)));
244  SP1 = fSPU; // correct calculation of SP1
245 
246  util.FromSDToMars(PC, RC, H, CH, SP1, fDJ, fDK, PD, RD);
247 
248  fPx = PD[0];
249  fPy = PD[1];
250  fPz = PD[2];
251 
252  for(int i = 0; i < 6; i++) for(int j = 0; j < 6; j++) { fCovMatrix66[i][j] = RD[i][j]; }
253 
254  fDPx = TMath::Sqrt(fabs(RD[0][0]));
255  fDPy = TMath::Sqrt(fabs(RD[1][1]));
256  fDPz = TMath::Sqrt(fabs(RD[2][2]));
257  fDX = TMath::Sqrt(fabs(RD[3][3]));
258  fDY = TMath::Sqrt(fabs(RD[4][4]));
259  fDZ = TMath::Sqrt(fabs(RD[5][5]));
260 }
261 
262 // ----- Constructor with track parameters in LAB -----------------------------------
263 FairTrackParP::FairTrackParP(TVector3 pos, TVector3 Mom, TVector3 posErr, TVector3 MomErr, Int_t Q, TVector3 o, TVector3 dj, TVector3 dk)
264  : FairTrackPar(pos.x(),pos.y(),pos.z(),Mom.x(),Mom.y(),Mom.z(),Q),
265  fU (0.),
266  fV (0),
267  fW (0),
268  fTV(0),
269  fTW(0),
270  fPx_sd(0.),
271  fPy_sd(0.),
272  fPz_sd(0.),
273  fDU(0.),
274  fDV(0.),
275  fDW(0.),
276  fDTV(0.),
277  fDTW(0.),
278  forigin(TVector3(0,0,0)),
279  fiver(TVector3(0,0,0)),
280  fjver(TVector3(0,0,0)),
281  fkver(TVector3(0,0,0)),
282  fSPU(0.)
283 {
284 // Reset();
285  SetPlane(o, dj, dk);
286 
287  SetPx(Mom.x());
288  SetPy(Mom.y());
289  SetPz(Mom.z());
290 
291  SetX(pos.x()); //x (lab)
292  SetY(pos.y()); //y (lab)
293  SetZ(pos.z()); //z (lab)
294  Double_t P = TMath::Sqrt(fPx*fPx +fPy*fPy +fPz*fPz );
295  if(Q!=0) { fq= Q/TMath::Abs(Q); }
296 
297 
298 
299 
300  FairGeaneUtil util;
301  TVector3 positionsd = util.FromMARSToSDCoord(TVector3(fX, fY, fZ), forigin, fiver, fjver, fkver);
302  fU = positionsd.X();
303  fV = positionsd.Y();
304  fW = positionsd.Z();
305  fQp = fq/P;
306 
307 
308 
309 
310  fDPx = MomErr.x(); //dpx
311  fDPy = MomErr.y(); //dpy
312  fDPz = MomErr.z(); //dpz
313  fDX = posErr.x(); //dpx
314  fDY = posErr.y(); //dpy
315  fDZ = posErr.z(); //dpz
316 
317  Double_t PD[3], RD[6][6], H[3], SP1, PC[3], RC[15];
318 
319  Int_t CH, IERR;
320 
321  PD[0] = fPx;
322  PD[1] = fPy;
323  PD[2] = fPz;
324 
325  for(int i = 0; i < 6; i++) for(int j = 0; j < 6; j++) { RD[i][j] = 0.; }
326  RD[0][0] = fDPx * fDPx;
327  RD[1][1] = fDPy * fDPy;
328  RD[2][2] = fDPz * fDPz;
329  RD[3][3] = fDX * fDX;
330  RD[4][4] = fDY * fDY;
331  RD[5][5] = fDZ * fDZ;
332 
333  for(int i = 0; i < 6; i++) for(int j = 0; j < 6; j++) { fCovMatrix66[i][j] = RD[i][j]; }
334 
335  // retrieve field
336  Double_t pnt[3];
337  pnt[0] = fX;
338  pnt[1] = fY;
339  pnt[2] = fZ;
340  FairRunAna* fRun = FairRunAna::Instance();
341  fRun->GetField()->GetFieldValue(pnt, H);
342 
343  CH=fq;
344 
345  util.FromMarsToSD(PD, RD, H, CH, fDJ, fDK, IERR, SP1, PC, RC);
346 
347  for(Int_t i=0; i<15; i++) { fCovMatrix[i] = RC[i]; }
348 
349  fTV = PC[1];
350  fTW = PC[2];
351 
352  fDQp = TMath::Sqrt(fabs(fCovMatrix[0]));
353  fDTV = TMath::Sqrt(fabs(fCovMatrix[5]));
354  fDTW = TMath::Sqrt(fabs(fCovMatrix[9]));
355  fDV = TMath::Sqrt(fabs(fCovMatrix[12]));
356  fDW = TMath::Sqrt(fabs(fCovMatrix[14]));
357 
358  fSPU = SP1;
359 
360  fPx_sd = fSPU*TMath::Sqrt( P*P / ( fTV*fTV + fTW*fTW + 1 ) );
361  fPy_sd = fTV * fPx_sd;
362  fPz_sd = fTW * fPx_sd;
363 
364 }
365 
366 // ----- Constructor with track parameters in LAB (with complete covariance matrix in MARS) -----------------------------------
367 FairTrackParP::FairTrackParP(TVector3 pos, TVector3 Mom,
368  Double_t covMARS[6][6], Int_t Q,
369  TVector3 o, TVector3 dj, TVector3 dk)
370  : FairTrackPar(pos.x(),pos.y(),pos.z(),Mom.x(),Mom.y(),Mom.z(),Q),
371  fU (0.),
372  fV (0),
373  fW (0),
374  fTV(0),
375  fTW(0),
376  fPx_sd(0.),
377  fPy_sd(0.),
378  fPz_sd(0.),
379  fDU(0.),
380  fDV(0.),
381  fDW(0.),
382  fDTV(0.),
383  fDTW(0.),
384  forigin(TVector3(0,0,0)),
385  fiver(TVector3(0,0,0)),
386  fjver(TVector3(0,0,0)),
387  fkver(TVector3(0,0,0)),
388  fSPU(0.)
389 {
390 // Reset();
391  SetPlane(o, dj, dk);
392 
393  SetPx(Mom.x());
394  SetPy(Mom.y());
395  SetPz(Mom.z());
396 
397  SetX(pos.x()); //x (lab)
398  SetY(pos.y()); //y (lab)
399  SetZ(pos.z()); //z (lab)
400  Double_t P = TMath::Sqrt(fPx*fPx +fPy*fPy +fPz*fPz );
401  if(Q!=0) { fq= Q/TMath::Abs(Q); }
402 
403  FairGeaneUtil util;
404  TVector3 positionsd = util.FromMARSToSDCoord(TVector3(fX, fY, fZ), forigin, fiver, fjver, fkver);
405  fU = positionsd.X();
406  fV = positionsd.Y();
407  fW = positionsd.Z();
408  fQp = fq/P;
409 
410  fDPx = TMath::Sqrt(fabs(covMARS[0][0])); //dpx
411  fDPy = TMath::Sqrt(fabs(covMARS[1][1])); //dpy
412  fDPz = TMath::Sqrt(fabs(covMARS[2][2])); //dpz
413  fDX = TMath::Sqrt(fabs(covMARS[3][3])); //dpx
414  fDY = TMath::Sqrt(fabs(covMARS[4][4])); //dpy
415  fDZ = TMath::Sqrt(fabs(covMARS[5][5])); //dpz
416 
417  Double_t PD[3], RD[6][6], H[3], SP1, PC[3], RC[15];
418  Int_t IERR, CH;
419 
420  PD[0] = fPx;
421  PD[1] = fPy;
422  PD[2] = fPz;
423 
424  for(int i = 0; i < 6; i++) for(int j = 0; j < 6; j++) { RD[i][j] = covMARS[i][j]; }
425  for(int i = 0; i < 6; i++) for(int j = 0; j < 6; j++) { fCovMatrix66[i][j] = RD[i][j]; }
426 
427 
428  // retrieve field
429  Double_t pnt[3];
430  pnt[0] = fX;
431  pnt[1] = fY;
432  pnt[2] = fZ;
433  FairRunAna* fRun = FairRunAna::Instance();
434  fRun->GetField()->GetFieldValue(pnt, H);
435 
436  CH=fq;
437 
438  util.FromMarsToSD(PD, RD, H, CH, fDJ, fDK, IERR, SP1, PC, RC);
439 
440  for(Int_t i=0; i<15; i++) { fCovMatrix[i] = RC[i]; }
441 
442  fTV = PC[1];
443  fTW = PC[2];
444 
445  fDQp = TMath::Sqrt(fabs(fCovMatrix[0]));
446  fDTV = TMath::Sqrt(fabs(fCovMatrix[5]));
447  fDTW = TMath::Sqrt(fabs(fCovMatrix[9]));
448  fDV = TMath::Sqrt(fabs(fCovMatrix[12]));
449  fDW = TMath::Sqrt(fabs(fCovMatrix[14]));
450 
451  fSPU = SP1;
452 
453  fPx_sd = fSPU*TMath::Sqrt( P*P / ( fTV*fTV + fTW*fTW + 1 ) );
454  fPy_sd = fTV * fPx_sd;
455  fPz_sd = fTW * fPx_sd;
456 }
457 
458 
459 // ctor from helix to parabola ---------------------------------------------------------------
460 // -- WARNING -- please pay attention to choose a plane which DOES NOT CONTAIN the TRACK, i.e.
461 // the plane normal axis u MUST NOT BE ORTHOGONAL to the particle momentum. If this happens
462 // you will get an error flag:
463 // IERR = 1 [when the particle moves perpendicular to u-axis (i.e. ON the CHOSEN PLANE)
464 // ==> v', w' are not definded.]
465 FairTrackParP::FairTrackParP(FairTrackParH* helix, TVector3 dj, TVector3 dk, Int_t& ierr)
466  : FairTrackPar(),
467  fU (0.),
468  fV (0),
469  fW (0),
470  fTV(0),
471  fTW(0),
472  fPx_sd(0.),
473  fPy_sd(0.),
474  fPz_sd(0.),
475  fDU(0.),
476  fDV(0.),
477  fDW(0.),
478  fDTV(0.),
479  fDTW(0.),
480  forigin(TVector3(0,0,0)),
481  fiver(TVector3(0,0,0)),
482  fjver(TVector3(0,0,0)),
483  fkver(TVector3(0,0,0)),
484  fSPU(0.)
485 {
486 
487  // q/p, lambda, phi --> q/p, v', w'
488  Double_t PC[3] = {helix->GetQp(), helix->GetLambda(), helix->GetPhi()};
489  Double_t RC[15];
490  helix->GetCov(RC);
491  // retrieve field
492  TVector3 xyz(helix->GetX(), helix->GetY(), helix->GetZ());
493  Double_t H[3], pnt[3];
494  pnt[0] = xyz.X();
495  pnt[1] = xyz.Y();
496  pnt[2] = xyz.Z();
497  FairRunAna* fRun = FairRunAna::Instance();
498  fRun->GetField()->GetFieldValue(pnt, H);
499  Int_t CH = helix->GetQ();
500 
501  Double_t DJ[3] = {dj.X(), dj.Y(), dj.Z()};
502  Double_t DK[3] = {dk.X(), dk.Y(), dk.Z()};
503 
504  Int_t IERR = 0;
505  Double_t SPU = 0;
506  Double_t PD[3], RD[15];
507 
508  FairGeaneUtil util;
509  util.FromSCToSD(PC, RC, H, CH, DJ, DK,
510  IERR, SPU, PD, RD);
511 
512  ierr = IERR;
513 
514  // MARS coordinates
515  TVector3 o(xyz.X(), xyz.Y(), xyz.Z());
516  TVector3 di = dj.Cross(dk);
517  TVector3 uvm = util.FromMARSToSDCoord(xyz, o, di, dj, dk);
518 
519  if(ierr == 0) SetTrackPar(uvm.Y(), uvm.Z(),
520  PD[1], PD[2], PD[0],
521  RD,
522  o, di, dj, dk, SPU);
523  else { cout << "FairTrackParP(FairTrackParH *) contructor ERROR: CANNOT convert helix to parabola" << endl; }
524 
525 }
526 
527 //define track in LAB
528 
529 void FairTrackParP::SetTrackPar(Double_t X, Double_t Y, Double_t Z,
530  Double_t Px, Double_t Py, Double_t Pz, Int_t Q,
531  Double_t CovMatrix[15],
532  TVector3 o, TVector3 di, TVector3 dj, TVector3 dk)
533 {
534 
535 
536 
537 
538  Reset();
539  SetPlane(o, dj, dk);
540 
541  Double_t P =TMath::Sqrt(Px*Px+Py*Py+Pz*Pz);
542 
543  if (Q!=0) { fq = TMath::Abs(Q)/Q; }
544  fQp = fq/P;
545 
546  SetX(X);
547  SetY(Y);
548  SetZ(Z);
549 
550  SetPx(Px);
551  SetPy(Py);
552  SetPz(Pz);
553 
554  for(Int_t i=0; i<15; i++) {
555  fCovMatrix[i]=CovMatrix[i];
556  }
557 
558  FairGeaneUtil util;
559 
560  TVector3 positionsd = util.FromMARSToSDCoord(TVector3(fX, fY, fZ), forigin, fiver, fjver, fkver);
561  fU = positionsd.X();
562  fV = positionsd.Y();
563  fW = positionsd.Z();
564  fQp = fq/P;
565 
566 
567  Double_t PD[3], RD[6][6], H[3], SP1, PC[3], RC[15];
568  Int_t IERR, CH;
569  PD[0] = Px;
570  PD[1] = Py;
571  PD[2] = Pz;
572 
573  for(int i = 0; i < 6; i++) for(int j = 0; j < 6; j++) { RD[i][j] = 0.; }
574  // retrieve field
575  Double_t pnt[3];
576  pnt[0] = fX;
577  pnt[1] = fY;
578  pnt[2] = fZ;
579  FairRunAna* fRun = FairRunAna::Instance();
580  fRun->GetField()->GetFieldValue(pnt, H);
581 
582  CH=fq;
583  util.FromMarsToSD(PD, RD, H, CH, fDJ, fDK, IERR, SP1, PC, RC);
584  fTV = PC[1];
585  fTW = PC[2];
586 
587  fDQp = TMath::Sqrt(fabs(fCovMatrix[0]));
588  fDTV = TMath::Sqrt(fabs(fCovMatrix[5]));
589  fDTW = TMath::Sqrt(fabs(fCovMatrix[9]));
590  fDV = TMath::Sqrt(fabs(fCovMatrix[12]));
591  fDW = TMath::Sqrt(fabs(fCovMatrix[14]));
592 
593 
594 
595  for(Int_t i=0; i<15; i++) { RC[i] = fCovMatrix[i]; }
596 
597  util.FromSDToMars(PC, RC, H, CH, SP1, fDJ, fDK, PD, RD);
598 
599  for(int i = 0; i < 6; i++) for(int j = 0; j < 6; j++) { fCovMatrix66[i][j] = RD[i][j]; }
600 
601 
602  fDPx = TMath::Sqrt(fabs(RD[0][0]));
603  fDPy = TMath::Sqrt(fabs(RD[1][1]));
604  fDPz = TMath::Sqrt(fabs(RD[2][2]));
605  fDX = TMath::Sqrt(fabs(RD[3][3]));
606  fDY = TMath::Sqrt(fabs(RD[4][4]));
607  fDZ = TMath::Sqrt(fabs(RD[5][5]));
608 
609  fSPU = SP1;
610 
611 
612 }
613 
614 //define track in SD
615 void FairTrackParP::SetTrackPar(Double_t v, Double_t w, Double_t Tv,
616  Double_t Tw, Double_t qp,
617  Double_t CovMatrix[15],
618  TVector3 o, TVector3 di, TVector3 dj, TVector3 dk, Double_t spu)
619 {
620  Reset();
621  SetPlane(o, dj, dk);
622  fU = 0;
623  fV = v;
624  fW = w;
625  fTV = Tv;
626  fTW = Tw;
627  fQp = qp;
628  fSPU = spu;
629 
630  Double_t P = TMath::Abs(1/fQp);
631  //fq= int (P * fQp);
632  fq = (int)TMath::Sign(1.0, fQp);
633  for(Int_t i=0; i<15; i++) {
634  fCovMatrix[i]=CovMatrix[i];
635  }
636 
637 
638  fPx_sd = fSPU*TMath::Sqrt( P*P / ( fTV*fTV + fTW*fTW + 1 ) );
639  fPy_sd = fTV * fPx_sd;
640  fPz_sd = fTW * fPx_sd;
641 
642  FairGeaneUtil util;
643  TVector3 position = util.FromSDToMARSCoord(TVector3(fU, fV, fW), forigin, fiver, fjver, fkver);
644  fX = position.X();
645  fY = position.Y();
646  fZ = position.Z();
647 
648  fDQp = TMath::Sqrt(fabs(fCovMatrix[0]));
649  fDTV = TMath::Sqrt(fabs(fCovMatrix[5]));
650  fDTW = TMath::Sqrt(fabs(fCovMatrix[9]));
651  fDV = TMath::Sqrt(fabs(fCovMatrix[12]));
652  fDW = TMath::Sqrt(fabs(fCovMatrix[14]));
653 
654 
655  Double_t PD[3],RD[6][6],H[3],PC[3],RC[15], SP1;
656  Int_t CH;
657  PC[0] = fQp;
658  PC[1] = fTV;
659  PC[2] = fTW;
660 
661  for(Int_t i=0; i<15; i++) {
662  RC[i]=fCovMatrix[i];
663  }
664 
665  // retrieve field
666  Double_t pnt[3];
667  pnt[0] = fX;
668  pnt[1] = fY;
669  pnt[2] = fZ;
670  FairRunAna* fRun = FairRunAna::Instance();
671  fRun->GetField()->GetFieldValue(pnt, H);
672 
673  CH=fq;
674 
675 
676  TVector3 momsd = TVector3(fPx_sd,fPy_sd,fPz_sd);
677  // SP1 = (momsd.Dot(fjver.Cross(fkver)))/TMath::Abs(momsd.Dot(fjver.Cross(fkver)));
678  SP1=fSPU;
679 
680  util.FromSDToMars(PC, RC, H, CH, SP1, fDJ, fDK, PD, RD);
681 
682  fPx = PD[0];
683  fPy = PD[1];
684  fPz = PD[2];
685 
686  for(int i = 0; i < 6; i++) for(int j = 0; j < 6; j++) { fCovMatrix66[i][j] = RD[i][j]; }
687 
688  fDPx = TMath::Sqrt(fabs(RD[0][0]));
689  fDPy = TMath::Sqrt(fabs(RD[1][1]));
690  fDPz = TMath::Sqrt(fabs(RD[2][2]));
691  fDX = TMath::Sqrt(fabs(RD[3][3]));
692  fDY = TMath::Sqrt(fabs(RD[4][4]));
693  fDZ = TMath::Sqrt(fabs(RD[5][5]));
694 
695 }
696 
697 void FairTrackParP::CalCov()
698 {
699  //not needed
700 }
701 
702 
703 // ----- Destructor ----------------------------------------------------
704 FairTrackParP::~FairTrackParP() {}
705 // -------------------------------------------------------------------------
706 
707 // ----- Public method Print -------------------------------------------
708 void FairTrackParP::Print(Option_t* option) const
709 {
710  cout << "Position : (";
711  cout.precision(2);
712  cout << fX << ", " << fY << ", " << fZ << ")" << endl;
713  cout << "Slopes : dx/dz = " << fTV << ", dy/dz = " << fTW << endl;
714  cout << "q/p = " << fQp << endl;
715 }
716 // -------------------------------------------------------------------------
717 /*Double_t FairTrackParP::GetDX()
718 { return fDX; }
719 Double_t FairTrackParP::GetDY()
720 { return fDY; }
721 Double_t FairTrackParP::GetDZ()
722 { return fDZ; }
723 */
724 Double_t FairTrackParP::GetV()
725 {
726  return fV;
727 }
728 
729 Double_t FairTrackParP::GetW()
730 {
731  return fW;
732 }
733 Double_t FairTrackParP::GetTV()
734 {
735  return fTV;
736 }
737 Double_t FairTrackParP::GetTW()
738 {
739  return fTW;
740 }
741 Double_t FairTrackParP::GetDV()
742 {
743  return fDV;
744 }
745 
746 Double_t FairTrackParP::GetDW()
747 {
748  return fDW;
749 }
750 Double_t FairTrackParP::GetDTV()
751 {
752  return fDTV;
753 }
754 Double_t FairTrackParP::GetDTW()
755 {
756  return fDTW;
757 }
758 /*
759 Double_t FairTrackParP::GetX()
760 {
761  return fX;
762 }
763 Double_t FairTrackParP::GetY()
764 {
765  return fY;
766 }
767 
768 Double_t FairTrackParP::GetZ()
769 {
770  return fZ;
771 }
772 
773 Double_t FairTrackParP::GetDPx()
774 {
775  return fDPx;
776 }
777 
778 Double_t FairTrackParP::GetDPy()
779 {
780  return fDPy;
781 }
782 
783 Double_t FairTrackParP::GetDPz()
784 {
785  return fDPz;
786 }
787 
788 Double_t FairTrackParP::GetDQp()
789 {
790  return fDQp;
791 }
792 */
793 TVector3 FairTrackParP::GetOrigin()
794 {
795  return forigin;
796 }
797 TVector3 FairTrackParP::GetIVer()
798 {
799  return fiver;
800 }
801 TVector3 FairTrackParP::GetJVer()
802 {
803  return fjver;
804 }
805 
806 TVector3 FairTrackParP::GetKVer()
807 {
808  return fkver;
809 }
810 
811 
812 
813 void FairTrackParP::Reset()
814 {
815  fU = 0.;
816  fV = 0.;
817  fW = 0.;
818  fTV = 0.;
819  fTW = 0.;
820  fDV = 0.;
821  fDW = 0.;
822  fDTV = 0.;
823  fDTW = 0.;
824  for(Int_t i= 0; i<15; i++) { fCovMatrix[i] = 0.; }
825  fPx_sd = 0.;
826  fPy_sd = 0.;
827  fPz_sd = 0.;
828 
829  //base class members
830  fX = fY = fZ = 0.;
831  fDX = fDY = fDZ = 0.;
832  fPx = fPy = fPz = 0.;
833  fDPx = fDPy = fDPz = 0.;
834  fQp = fDQp = fq = 0;
835 
836 }
837 
838 
839 void FairTrackParP::SetPlane(TVector3 o, TVector3 j, TVector3 k)
840 {
841  // origin
842  forigin = TVector3(o.X(), o.Y(), o.Z());
843 
844  // check unity
845  j.SetMag(1.);
846  k.SetMag(1.);
847  // check orthogonality
848  if(j* k != 0) {
849 
850  k = (j.Cross(k)).Cross(j);
851  }
852 
853  // plane
854  // i
855  TVector3 i = j.Cross(k);
856  i.SetMag(1.);
857  fDI[0] = i.X();
858  fDI[1] = i.Y();
859  fDI[2] = i.Z();
860  fiver = TVector3(fDI[0],fDI[1],fDI[2]);
861  // j
862  fDJ[0] = j.X();
863  fDJ[1] = j.Y();
864  fDJ[2] = j.Z();
865  fjver = TVector3(fDJ[0],fDJ[1],fDJ[2]);
866  // k
867  fDK[0] = k.X();
868  fDK[1] = k.Y();
869  fDK[2] = k.Z();
870  fkver = TVector3(fDK[0],fDK[1],fDK[2]);
871 
872 
873 
874 }
875 
876 
877 void FairTrackParP::SetTransportMatrix(Double_t mat[][5])
878 {
879  for(int i = 0; i < 5; i++) for(int j = 0; j < 5; j++) { ftrmat[i][j] = mat[i][j]; }
880 
881 }
882 
883 void FairTrackParP::GetTransportMatrix(Double_t mat[][5])
884 {
885  for(int i = 0; i < 5; i++) for(int j = 0; j < 5; j++) { mat[i][j] = ftrmat[i][j]; }
886 }
887 
888 void FairTrackParP::GetCovQ(Double_t* CovQ)
889 {
890  // return error matrix in 1/p instead of q/p
891 
892  for(int i = 0; i < 15; i++) {
893  CovQ[i] = fCovMatrix[i];
894  if(fq!=0) {
895  if(i == 0) { CovQ[i] = CovQ[i] / (fq * fq); }
896  if(i > 0 && i < 5) { CovQ[i] = CovQ[i] / fq; }
897  }
898  }
899 }