a19b177aacf85b0091b8d3dcab11148e87e374b2
[u/mrichter/AliRoot.git] / HMPID / AliHMPIDtrack.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 #include "AliESDtrack.h" 
17 #include "AliTracker.h" 
18 #include "AliHMPIDtrack.h" 
19
20 ClassImp(AliHMPIDtrack)
21
22 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
23 AliHMPIDtrack::AliHMPIDtrack():AliKalmanTrack()
24 {
25   //
26   // def. ctor
27   //
28 }                                
29 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
30 AliHMPIDtrack::AliHMPIDtrack(const AliHMPIDtrack& t):AliKalmanTrack(t)
31 {
32   //
33   // cctor.
34   //
35  }                                
36 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
37 AliHMPIDtrack::AliHMPIDtrack(const AliESDtrack& t):AliKalmanTrack()
38 {
39   //
40   // Constructor from AliESDtrack
41   //
42   SetLabel(t.GetLabel());
43   SetChi2(0.);
44   SetMass(t.GetMass());
45
46   Set(t.GetX(),t.GetAlpha(),t.GetParameter(),t.GetCovariance());
47
48   if ((t.GetStatus()&AliESDtrack::kTIME) == 0) return;
49   StartTimeIntegral();
50   Double_t times[10]; t.GetIntegratedTimes(times); SetIntegratedTimes(times);
51   SetIntegratedLength(t.GetIntegratedLength());
52 }              
53 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
54 AliHMPIDtrack& AliHMPIDtrack::operator=(const AliHMPIDtrack &/*source*/)
55 {
56   // ass. op.
57   
58   return *this;
59 }
60 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
61 Bool_t AliHMPIDtrack::PropagateTo(Double_t xk, Double_t xx0, Double_t xrho)
62 {
63   //
64   // Propagates this track to a reference plane defined by "xk" [cm] 
65   // correcting for the mean crossed material.
66   // Arguments:
67   // "xx0"  - thickness/rad.length [units of the radiation length] 
68   // "xrho" - thickness*density    [g/cm^2] 
69   //  Returns: kTRUE if the track propagates to plane, else kFALSE
70  
71   if (xk == GetX()) {
72     return kTRUE;
73   }
74   Double_t bz   = GetBz();
75   if (!AliExternalTrackParam::PropagateTo(xk,bz)) {
76     return kFALSE;
77   }
78   if (!AliExternalTrackParam::CorrectForMeanMaterial(xx0,xrho,GetMass())) { 
79     return kFALSE;
80   }
81   return kTRUE;            
82 }//PropagateTo()  
83 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
84 Bool_t AliHMPIDtrack::Rotate(Double_t alpha, Bool_t absolute)
85 {
86   //
87   // Rotates track parameters in R*phi plane
88   // if absolute rotation alpha is in global system
89   // otherwise alpha rotation is relative to the current rotation angle
90   //  
91
92   if (absolute) {
93     alpha -= GetAlpha();
94   }
95  
96
97   return AliExternalTrackParam::Rotate(GetAlpha()+alpha);
98
99 }   
100 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++  
101 Int_t AliHMPIDtrack::GetProlongation(Double_t xk, Double_t &y, Double_t &z)
102 {
103   //
104   // Find a prolongation at given x
105   // Return 0 if it does not exist
106   //  
107
108   Double_t bz = GetBz();
109
110   if (!AliExternalTrackParam::GetYAt(xk,bz,y)) {
111     return 0;
112   }
113   if (!AliExternalTrackParam::GetZAt(xk,bz,z)) {
114     return 0;
115   }
116
117   return 1;  
118
119 }
120  
121 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
122 Bool_t   AliHMPIDtrack::PropagateToR(Double_t r,Double_t step)
123 {
124   //
125   // Propagate track to the radial position
126   // Rotation always connected to the last track position
127   //
128
129   Double_t xyz0[3];
130   Double_t xyz1[3];
131   Double_t y;
132   Double_t z; 
133
134   Double_t radius = TMath::Sqrt(GetX()*GetX() + GetY()*GetY());
135   // Direction +-
136   Double_t dir    = (radius > r) ? -1.0 : 1.0;   
137
138   for (Double_t x = radius+dir*step; dir*x < dir*r; x += dir*step) {
139
140     GetXYZ(xyz0);       
141     Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
142     Rotate(alpha,kTRUE);
143     GetXYZ(xyz0);       
144     if (!GetProlongation(x,y,z)) return kFALSE;
145     xyz1[0] = x * TMath::Cos(alpha) + y * TMath::Sin(alpha); 
146     xyz1[1] = x * TMath::Sin(alpha) - y * TMath::Cos(alpha);
147     xyz1[2] = z;
148     Double_t param[7];
149     AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
150     if (param[1] <= 0) {
151       param[1] = 100000000;
152     }
153     PropagateTo(x,param[1],param[0]*param[4]);
154
155   } 
156
157   GetXYZ(xyz0); 
158   Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
159   Rotate(alpha,kTRUE);
160   GetXYZ(xyz0); 
161   if (!GetProlongation(r,y,z)) return kFALSE;
162   xyz1[0] = r * TMath::Cos(alpha) + y * TMath::Sin(alpha); 
163   xyz1[1] = r * TMath::Sin(alpha) - y * TMath::Cos(alpha);
164   xyz1[2] = z;
165   Double_t param[7];
166   AliTracker::MeanMaterialBudget(xyz0,xyz1,param);
167
168   if (param[1] <= 0) {
169     param[1] = 100000000;
170   }
171  
172   return PropagateTo(r,param[1],param[0]*param[4]);
173 }
174 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
175 Double_t AliHMPIDtrack::GetPredictedChi2(const AliCluster3D *c) const {
176   //
177   // Arguments: AliCluster3D
178   // Returns:   Chi2 of track for the cluster
179   Double_t      p[3]={c->GetX(),       c->GetY(),       c->GetZ()};
180   Double_t  covyz[3]={c->GetSigmaY2(), c->GetSigmaYZ(), c->GetSigmaZ2()};
181   Double_t covxyz[3]={c->GetSigmaX2(), c->GetSigmaXY(), c->GetSigmaXZ()};
182   return AliExternalTrackParam::GetPredictedChi2(p, covyz, covxyz);
183 }//GetPredictedChi2()
184 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
185 Bool_t AliHMPIDtrack::PropagateTo(const AliCluster3D *c) {
186   //
187   // Arguments: AliCluster3D
188   // Returns: kTRUE if the track propagates to the plane of the cluster
189   Double_t      oldX=GetX(),   oldY=GetY(), oldZ=GetZ();
190   Double_t      p[3]={c->GetX(), c->GetY(), c->GetZ()};
191   Double_t  covyz[3]={c->GetSigmaY2(), c->GetSigmaYZ(), c->GetSigmaZ2()};
192   Double_t covxyz[3]={c->GetSigmaX2(), c->GetSigmaXY(), c->GetSigmaXZ()};
193   Double_t bz=GetBz();
194   
195   if(!AliExternalTrackParam::PropagateTo(p, covyz, covxyz, bz)) return kFALSE;
196   if(IsStartedTimeIntegral()) 
197     {
198       Double_t d = TMath::Sqrt((GetX()-oldX)*(GetX()-oldX) + (GetY()-oldY)*(GetY()-oldY) + (GetZ()-oldZ)*(GetZ()-oldZ));
199       if (GetX()<oldX) d=-d;
200       AddTimeStep(d);
201   }
202   return kTRUE;
203 }//PropagateTo()
204 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
205 Double_t AliHMPIDtrack::GetBz() const {
206   // 
207   // Arguments: none
208   // Returns: Bz component of the magnetic field (kG)
209   //
210   if (AliTracker::UniformField()) return AliTracker::GetBz();
211   Double_t r[3]; GetXYZ(r);
212   return AliTracker::GetBz(r);
213 }//GetBz()
214 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
215 Bool_t AliHMPIDtrack::Intersect(Double_t pnt[3], Double_t norm[3], Double_t bz) const {
216   //+++++++++++++++++++++++++++++++++++++++++    
217   // Origin: K. Shileev (Kirill.Shileev@cern.ch)
218   // Finds point of intersection (if exists) of the helix with the plane. 
219   // Stores result in fX and fP.   
220   // Arguments: planePoint,planeNorm - the plane defined by any plane's point 
221   // and vector, normal to the plane
222   // Returns: kTrue if helix intersects the plane, kFALSE otherwise.
223   //+++++++++++++++++++++++++++++++++++++++++    
224   Double_t x0[3]; GetXYZ(x0); //get track position in MARS
225   
226   //estimates initial helix length up to plane
227   Double_t s=(pnt[0]-x0[0])*norm[0] + (pnt[1]-x0[1])*norm[1] + (pnt[2]-x0[2])*norm[2];
228   Double_t dist=99999,distPrev=dist;
229   Double_t x[3],p[3]; 
230   while(TMath::Abs(dist)>0.00001){
231     //calculates helix at the distance s from x0 ALONG the helix
232     Propagate(s,x,p,bz);
233     //distance between current helix position and plane
234     dist=(x[0]-pnt[0])*norm[0]+(x[1]-pnt[1])*norm[1]+(x[2]-pnt[2])*norm[2];  
235     if(TMath::Abs(dist) >= TMath::Abs(distPrev)) {return kFALSE;}
236     distPrev=dist;
237     s-=dist;
238   }
239   //on exit pnt is intersection point,norm is track vector at that point, 
240   //all in MARS
241   for (Int_t i=0; i<3; i++) {pnt[i]=x[i]; norm[i]=p[i];}
242   return kTRUE;
243 }//Intersect()
244 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
245 Bool_t AliHMPIDtrack::Intersect(AliHMPIDtrack *pTrk,Double_t pnt[3], Double_t norm[3]) {
246   // Finds point of intersection (if exists) of the helix with the plane. 
247   // Stores result in fX and fP.   
248   // Arguments: planePoint,planeNorm - the plane defined by any plane's point 
249   // and vector, normal to the plane
250   // Returns: kTrue if helix intersects the plane, kFALSE otherwise.
251  
252   Double_t x0[3]; pTrk->GetXYZ(x0); //get track position in MARS
253   Double_t dist=99999,distPrev=dist;
254   Double_t x[3],p[3],
255   pntrad= TMath::Sqrt(pnt[0]*pnt[0]+pnt[1]*pnt[1]);
256   while(TMath::Abs(dist)> 0.000001){//0.00001){
257     //calculates helix at the distance s from x0 ALONG the helix
258     pTrk->PropagateTo(pntrad);pTrk->GetXYZ(x);pTrk->GetPxPyPz(p); 
259     //distance between current helix position and plane
260     dist=(x[0]-pnt[0])*norm[0]+(x[1]-pnt[1])*norm[1]+(x[2]-pnt[2])*norm[2];
261     pntrad=pntrad-dist*0.7;
262     //Printf("--- 111 --- dist %lf",dist);
263     if(TMath::Abs(2.0*dist) >= TMath::Abs(distPrev)) {return kFALSE;}
264     distPrev=dist;
265   }
266   //on exit pnt is intersection point,norm is track vector at that point, 
267   //Printf("--- 222 --- dist %lf",dist);
268 //  Printf("");
269   for (Int_t i=0; i<3; i++) {pnt[i]=x[i]; norm[i]=p[i];}
270   return kTRUE;
271 }//Intersect()
272 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
273 void AliHMPIDtrack::Propagate(Double_t len, Double_t x[3],Double_t p[3], Double_t bz) const {
274   //+++++++++++++++++++++++++++++++++++++++++    
275   // Origin: K. Shileev (Kirill.Shileev@cern.ch)
276   // Extrapolate track along simple helix in magnetic field
277   // Arguments: len -distance alogn helix, [cm]
278   //            bz  - mag field, [kGaus]   
279   // Returns: x and p contain extrapolated positon and momentum  
280   // The momentum returned for straight-line tracks is meaningless !
281   //+++++++++++++++++++++++++++++++++++++++++    
282   GetXYZ(x);    
283   if (OneOverPt() < kAlmost0 || TMath::Abs(bz) < kAlmost0Field ){ //straight-line tracks
284      Double_t unit[3]; GetDirection(unit);
285      x[0]+=unit[0]*len;   
286      x[1]+=unit[1]*len;   
287      x[2]+=unit[2]*len;
288
289      p[0]=unit[0]/kAlmost0;   
290      p[1]=unit[1]/kAlmost0;   
291      p[2]=unit[2]/kAlmost0;   
292   } else {
293      GetPxPyPz(p);
294      Double_t pp=GetP();
295      Double_t a = -kB2C*bz*GetSign(); ////////// what is kB2C
296      Double_t rho = a/pp;
297      x[0] += p[0]*TMath::Sin(rho*len)/a - p[1]*(1-TMath::Cos(rho*len))/a;
298      x[1] += p[1]*TMath::Sin(rho*len)/a + p[0]*(1-TMath::Cos(rho*len))/a;
299      x[2] += p[2]*len/pp;
300      Double_t p0=p[0];
301      p[0] = p0  *TMath::Cos(rho*len) - p[1]*TMath::Sin(rho*len);
302      p[1] = p[1]*TMath::Cos(rho*len) + p0  *TMath::Sin(rho*len);
303   }
304 }//Propagate()
305 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
306 Bool_t AliHMPIDtrack::Update(const AliHMPIDCluster *pClu, Double_t /*chisq*/, Int_t /*index*/)
307 {
308   //
309   // Arguments: AliCluster3D, chi sq, and clu index
310   // Returns: kTRUE if the track parameters are successfully updated
311   Double_t p[2]={pClu->GetY(), pClu->GetZ()};
312   Double_t cov[3]={pClu->GetSigmaY2(), 0., pClu->GetSigmaZ2()};
313   if (!AliExternalTrackParam::Update(p,cov)) return kFALSE;
314
315   /*
316   AliTracker::FillResiduals(this,p,cov,pClu->GetVolumeId());
317
318   Int_t n=GetNumberOfClusters();
319   fIndex[n]=index;
320   SetNumberOfClusters(n+1);
321   SetChi2(GetChi2()+chisq);
322 */
323   return kTRUE;
324
325 }//Update()
326 //++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++