EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PndPidTrackInfo.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PndPidTrackInfo.cxx
1 #include "PndPidCorrelator.h"
2 #include "PndPidCandidate.h"
3 #include "PndTrack.h"
4 #include "FairTrackParH.h"
5 #include <FairGeanePro.h>
6 
7 #if _TODAY_
8 #include "Fitter/PndVtxFitterParticle.h"
9 #endif
10 #include "FairRunAna.h"
11 #include "FairField.h"
12 
13 #include "TVector3.h"
14 #include "TDatabasePDG.h"
15 #include <cmath>
16 
17 //_________________________________________________________________
18 Bool_t PndPidCorrelator::GetTrackInfo(PndTrack* track, PndPidCandidate* pidCand)
19 {
20 #if _TODAY_
21  static PndVtxFitterParticle covTool;
22 #endif
23 
24  Int_t charge = TMath::Sign(1, track->GetParamFirst().GetQ());
25  pidCand->SetCharge(charge);
26 
27  TVector3 first(track->GetParamFirst().GetX(),
28  track->GetParamFirst().GetY(),
29  track->GetParamFirst().GetZ());
30  TVector3 last(track->GetParamLast().GetX(),
31  track->GetParamLast().GetY(),
32  track->GetParamLast().GetZ());
33  FairTrackParP par = track->GetParamFirst();
34  Int_t ierr = 0;
35  FairTrackParH *helix = new FairTrackParH(&par, ierr);
36 
37  if (fGeanePro) // Overwrites vertex if Geane is used
38  {
39  TVector3 momentum, pocaz;
40  FairGeanePro *fPro0 = new FairGeanePro();
41  FairTrackParH *fRes= new FairTrackParH();
42 
43  // track back to Origin
44  // [ralfk: changed to propagate to z axis - 03/2011]
45 #if 1
46  fPro0->SetPoint(TVector3(0,0,0));
47  fPro0->PropagateToPCA(1, -1);
48  //assert(0);
49 #else
50  // Propagatetrack back to z Axis
51  fPro0->PropagateToPCA(2, -1);// track back to z axis
52  TVector3 ex1(0.,0.,-50.); // virtual wire, dimensions chosen arbitrarily
53  TVector3 ex2(0.,0.,100.); // we expect fast decaying tracks to be close to that
54  fPro0->SetWire(ex1,ex2);
55 #endif
56 
57  Bool_t rc = fPro0->Propagate(helix, fRes, fPidHyp*charge);
58  if (!rc)
59  {
60  std::cout << "-W- PndPidCorrelator::GetTrackInfo :: Failed backward propagation" << std::endl;
61  if (fVerbose>0) helix->Print();
62  return kFALSE;
63  }
64 
65  pocaz.SetXYZ(fRes->GetX(), fRes->GetY(), fRes->GetZ()); //cm ??
66  momentum = fRes->GetMomentum();
67  {
68  // A hack to extract just 1/p diagonal component;
69  double buffer[15];
70  fRes->GetCov(buffer);
71  //printf("%f\n", sqrt(buffer[0])*pow(momentum.Mag(), 2));
72  for(unsigned iq=1; iq<15; iq++)
73  buffer[iq] = 0.0;
74  pidCand->SetHelixCov(buffer);
75  }
76 
77  // ********
78  // ********
79 
80  // ******** this block replaces
81 #if _TODAY_
82  TVector3 di = momentum;
83  di.SetMag(1.);
84  TVector3 dj = di.Orthogonal();
85  TVector3 dk = di.Cross(dj);
86  FairTrackParP *fParab = new FairTrackParP(fRes, dj, dk, ierr);
87  // ******* that line (recipe by Lia)
88  //FairTrackParP *fParab = new FairTrackParP(fRes, TVector3(1.,0.,0.), TVector3(0.,1.,0.), ierr);
89  // ******* that line (recipe by Lia)
90 
91  // ********
92  // ********
93 
94  Double_t globalCov[6][6];
95  fParab->GetMARSCov(globalCov);
96 
97  Int_t ii,jj;
98  TMatrixD err(6,6);
99  for (ii=0;ii<6;ii++) for(jj=0;jj<6;jj++) err[ii][jj]=globalCov[ii][jj];
100 #endif
101 
102  //err.Print();
103 
104  TLorentzVector lv;
105  lv.SetVectM(momentum, TDatabasePDG::Instance()->GetParticle(fPidHyp)->Mass()); // set mass hypothesis
106  Float_t energy = lv.E();
107 #if _TODAY_
108  TMatrixD mat = covTool.GetConverted7(covTool.GetFitError(lv, err));
109 #endif
110 
111  pidCand->SetPosition(pocaz);
112  pidCand->SetMomentum(momentum);
113  pidCand->SetEnergy(energy);
114 #if _TODAY_
115  pidCand->SetCov7(mat);
116 
117  // Adding the helix parameters to the candidate (TFitParams)
118  // It is nice to know at analysis stage
119  //rho helix: (D0,Phi0,rho(omega),Z0,tan(dip))
120  //fair helix: (q/p,lambda, phi, y_perp, z_perp)
121  Double_t Q=fRes->GetQ();
122  //if(0==Q) ??? break/return?;
123  Double_t pnt[3], Bf[3];
124  pnt[0]=pocaz.X();
125  pnt[1]=pocaz.Y();
126  pnt[2]=pocaz.Z();
127  FairRunAna::Instance()->GetField()->GetFieldValue(pnt, Bf); //[kGs]
128  //Double_t B = sqrt(Bf[0]*Bf[0]+Bf[1]*Bf[1]+Bf[2]*Bf[2]);
129  Double_t B = Bf[2];
130  Double_t qBc = -0.000299792458*B*Q;
131  Double_t icL = 1. / cos(fRes->GetLambda()); // inverted for practical reasons (better to multiply than to divide)
132  Double_t icLs = icL*icL;
133  Double_t helixparams[5];
134  //helixparams[0]=pocaz.Perp(); //D0
135  helixparams[0]=fRes->GetY_sc() ; //D0
136  //helixparams[1]=fRes->GetPhi(); //phi0
137  helixparams[1]=fRes->GetMomentum().Phi(); //phi0
138  helixparams[2]=qBc/(fRes->GetMomentum().Perp()); //omega=rho=1/R[cm]=-2.998e-4*B[kGs]*Q[e]/p_perp[GeV/c]
139  //helixparams[3]=pocaz.Z(); //z0
140  helixparams[3]=fRes->GetZ_sc()*icL; //z0
141  helixparams[4]=tan(fRes->GetLambda()); //lambda(averey)=cot(theta)=tan(lambda(geane))
142  pidCand->SetHelixParams(helixparams);
143  Double_t fairhelixcov[15];
144  fRes->GetCov(fairhelixcov);
145  Double_t rhohelixcov[15];
146  // in the poca to z axis yperp=D0, x_perp^2+z_perp^2 = z_perp/cos(Lambda)= Z0
147  rhohelixcov[0] = fairhelixcov[12]; // sigma^2 D0
148  rhohelixcov[1] = fairhelixcov[10]; // cov D0 - Phi0
149  rhohelixcov[2] = fairhelixcov[3] * qBc * icL; // cov D0 - rho
150  rhohelixcov[3] = fairhelixcov[13] * icL; // cov D0 - Z0
151  rhohelixcov[4] = fairhelixcov[7] * icLs; // cov D0 - tan(dip)
152  rhohelixcov[5] = fairhelixcov[9]; // sigma^2 Phi0
153  rhohelixcov[6] = fairhelixcov[2] * qBc * icL; // cov Phi0 - rho
154  rhohelixcov[7] = fairhelixcov[11] * icL; // cov Phi0 - Z0
155  rhohelixcov[8] = fairhelixcov[6] * icLs; // cov Phi0 - tan(dip)
156  rhohelixcov[9] = fairhelixcov[0] * qBc * qBc * icLs; // sigma^2 rho
157  rhohelixcov[10] = fairhelixcov[4] * qBc * icLs; // cov rho - Z0
158  rhohelixcov[11] = fairhelixcov[1] * qBc * icL * icLs; // cov rho - tan(dip)
159  rhohelixcov[12] = fairhelixcov[14] * icLs; // sigma^2 Z0
160  rhohelixcov[13] = fairhelixcov[8] * icL * icLs; //cov Z0 - tan(dip)
161  rhohelixcov[14] = fairhelixcov[5] * icLs * icLs; // sigma^2 tan(dip) - from
162  pidCand->SetHelixCov(rhohelixcov);
163 #endif
164  }
165  else
166  {
167  std::cout << std::endl << "-E- PndPidCorrelator::GetTrackInfo :: NO GEANE - no propagation available" << std::endl;
168  }
169  pidCand->SetFirstHit(first);
170  pidCand->SetLastHit(last);
171  pidCand->SetDegreesOfFreedom(track->GetNDF());
172  pidCand->SetFitStatus(track->GetFlag());
173  pidCand->SetChiSquared(track->GetChi2());
174 
175  return kTRUE;
176 }
177 
178 ClassImp(PndPidCorrelator)