Added support for NEWIO, merged cern-hlt tree, updated to latest track candidate...
[u/mrichter/AliRoot.git] / HLT / kalman / AliL3KalmanTrack.cxx
1 #include "TMath.h"
2 #include "AliL3KalmanTrack.h"
3 #include "AliL3SpacePointData.h"
4 #include "AliL3Logging.h"
5 #include "AliL3Transform.h"
6 ClassImp(AliL3KalmanTrack)
7
8 // Class for kalman tracks
9
10 AliL3KalmanTrack::AliL3KalmanTrack()
11 {
12   fMaxChi2 = 12;
13   // Constructor
14 }
15
16 AliL3KalmanTrack::~AliL3KalmanTrack()
17 {
18   // Destructor
19 }
20
21 Int_t AliL3KalmanTrack::Init(AliL3Track *track, AliL3SpacePointData *points, UInt_t pos)
22 {
23   fX = points[pos].fX; 
24
25   fP0 = points[pos].fY;
26   fP1 = points[pos].fZ; 
27   fP3 = track->GetTgl(); //cout << fP3; 
28   if (TMath::Abs(fP3) > 1.2) return 0; //From AliTPCtrackerMI
29   fP4 = 0.5*(-track->GetCharge()*1./track->GetPt())*0.0029980*AliL3Transform::GetBField(); // Calculation of the curvature 
30   
31   if (TMath::Abs(fP4) >= 0.0066) return 0; // From AliTPCtrackerMI
32
33   Float_t centerX = track->GetFirstPointX() - ((track->GetPt()/(0.0029980*AliL3Transform::GetBField())) * TMath::Cos(track->GetPsi() + track->GetCharge() * 0.5 * 3.14159265358979323846));
34   Float_t centerY = track->GetFirstPointY() - ((track->GetPt()/(0.0029980*AliL3Transform::GetBField())) * TMath::Sin(track->GetPsi() + track->GetCharge() * 0.5 * 3.14159265358979323846));
35   fP2 = fP4*centerX; // Curvature times center of curvature
36   
37   if (TMath::Abs(fP4*fX - fP2) >= 0.9) // What's this??
38     {
39       return 0;
40     }
41   //cout << track->GetCharge() << endl; 
42   //cout << "AliL3KalmanTrack::Init, " << fP0 << " " << fP1 << " " << fP2 << " " << fP3 << " " << fP4 << endl;
43   //cout << fP2 << " " << fP4 << endl;
44   Float_t num = 12;
45
46   fC00 = points[pos].fSigmaY2; 
47   fC10 = 0; fC11 = points[pos].fSigmaZ2;
48   fC20 = 0; fC21 = 0; fC22 = 5e-05;
49   fC30 = 0; fC31 = 0; fC32 = 0; fC33 = 5e-05;
50   fC40 = 0; fC41 = 0; fC42 = 0; fC43 = 0; fC44 = 5e-09;
51   /*fC00 = points[pos].fSigmaY2;
52   fC10 = 0; fC11 = points[pos].fSigmaZ2;
53   fC20 = 0; fC21 = 0; fC22 = (0.1/TMath::Sqrt(num))*(0.1/TMath::Sqrt(num));
54   fC30 = 0; fC31 = -points[pos].fSigmaZ2; fC32 = 0; fC33 = 2*(points[pos].fSigmaZ2);
55   fC40 = -points[pos].fSigmaY2; fC41 = 0; fC42 = 0; fC43 = 0; fC44 = 2*(points[pos].fSigmaY2);*/
56   //cout << "Init: errors " << fC00 << " " << fC11 << " " << fC22 << " " << fC33 << " " << fC44 << endl;
57
58   fChisq = 0;
59
60   return 1;
61 }
62
63 Int_t AliL3KalmanTrack::MakeTrackSeed(AliL3SpacePointData *points1, UInt_t pos1, AliL3SpacePointData *points2, UInt_t pos2, AliL3SpacePointData *points3, UInt_t pos3)
64 {
65   // Make track seed based on three clusters 
66   fChisq = 0;
67   fX = points1[pos1].fX;
68
69   fP0 = points1[pos1].fY; 
70   fP1 = points1[pos1].fZ; 
71
72   Float_t X2 = points2[pos2].fX;
73   Float_t Y2 = points2[pos2].fY;
74   Float_t Z2 = points2[pos2].fZ;
75
76   Float_t X3 = points3[pos3].fX;
77   Float_t Y3 = points3[pos3].fY;
78   Float_t Z3 = points3[pos3].fZ;
79
80   Float_t ZZ = fP1 - ((fP1 - Z3)/(fX-X3))*(fX-X2);
81   if (TMath::Abs(ZZ - Z2) > 10) return 0; //What's this?? (fP1 - Z3)/(fX-X3)*(fX-X2) is an angle
82
83   // It may make no difference. Check on a big event??.
84   if ((X2-fX)*(0-Y2)-(0-X2)*(Y2-fP0) == 0) return 0; //Straight seed
85
86   // Initial approximation of the state vector
87   fP2 = f2(fX,fP0,X2,Y2,X3,Y3);
88   fP3 = f3(fX,fP0,X2,Y2,fP1,Z2);
89   fP4 = f4(fX,fP0,X2,Y2,X3,Y3);
90   if (TMath::Abs(fP3) > 1.2) return 0; //From AliTPCtrackerMI
91   if (TMath::Abs(fP4) >= 0.0066) return 0; // From AliTPCtrackerMI
92
93   //cout << "AliL3KalmanTrack::MakeTrackSeed, " << fP0 << " " << fP1 << " " << fP2 << " " << fP3 << " " << fP4 << endl;
94   //cout << fP2 << " " << fP4 << endl;
95
96   // Initial approximation of the covariance matrix
97   Float_t sy1=points1[pos1].fSigmaY2;
98   Float_t sz1=points1[pos1].fSigmaZ2;
99   Float_t sy2=points2[pos2].fSigmaY2;
100   Float_t sz2=points2[pos2].fSigmaZ2;
101   Float_t sy3=25000*fP4*fP4+0.1;
102   Float_t sy=0.1;
103   Float_t sz=0.1;
104   
105   Float_t f40=(f4(fX,fP0+sy,X2,Y2,X3,Y3)-fP4)/sy;
106   Float_t f42=(f4(fX,fP0,X2,Y2+sy,X3,Y3)-fP4)/sy;
107   Float_t f43=(f4(fX,fP0,X2,Y2,X3,Y3+sy)-fP4)/sy;
108   Float_t f20=(f2(fX,fP0+sy,X2,Y2,X3,Y3)-fP2)/sy;
109   Float_t f22=(f2(fX,fP0,X2,Y2+sy,X3,Y3)-fP2)/sy;
110   Float_t f23=(f2(fX,fP0,X2,Y2,X3,Y3+sy)-fP2)/sy;
111   Float_t f30=(f3(fX,fP0+sy,X2,Y2,fP1,Z2)-fP3)/sy;
112   Float_t f31=(f3(fX,fP0,X2,Y2,fP1+sz,Z2)-fP3)/sz;
113   Float_t f32=(f3(fX,fP0,X2,Y2+sy,fP1,Z2)-fP3)/sy;
114   Float_t f34=(f3(fX,fP0,X2,Y2,fP1,Z2+sz)-fP3)/sz;
115   
116   fC00 = sy1;
117   fC10 = 0;  
118   fC11 = sz1;
119   fC20 = f20*sy1;  
120   fC21 = 0;  
121   fC22 = f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
122   fC30 = f30*sy1;  
123   fC31 = f31*sz1;  
124   fC32 = f30*sy1*f20+f32*sy2*f22;
125   fC33 = f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
126   fC40 = f40*sy1; 
127   fC41 = 0; 
128   fC42 = f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
129   fC43 = f30*sy1*f40+f32*sy2*f42; 
130   fC44 = f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
131   //cout << "Errors: " << fC00 << " " << fC11 << " " << fC20 << " " << fC22 << " " << " " << fC30 << " " << fC31 << " " << fC32 << " " << fC33 << " " << fC40 << " " << fC42 << " " << fC43 << " " << fC44 << endl;
132   
133   return 1;
134 }
135
136 Int_t AliL3KalmanTrack::Propagate(AliL3SpacePointData *points, UInt_t pos)
137 {
138   // Propagates track to the plane of the next found cluster
139   Float_t Xold = fX; // X position for previous space point
140   Float_t Xnew = points[pos].fX; // X position of current space point
141   Float_t dx = Xnew - Xold;
142   Float_t Yold = fP0; // Y position of old point
143   Float_t Zold = fP1; // Z position of old point
144   Float_t par2old = fP2;  
145   Float_t par3old = fP3;
146   Float_t par4old = fP4;
147   Float_t oldfC00 = fC00;
148   Float_t oldfC10 = fC10;
149   Float_t oldfC11 = fC11;
150   Float_t oldfC20 = fC20;
151   Float_t oldfC21 = fC21;
152   Float_t oldfC22 = fC22;
153   Float_t oldfC30 = fC30;
154   Float_t oldfC31 = fC31;
155   Float_t oldfC32 = fC32;
156   Float_t oldfC33 = fC33;
157   Float_t oldfC40 = fC40;
158   Float_t oldfC41 = fC41;
159   Float_t oldfC42 = fC42;
160   Float_t oldfC43 = fC43;
161   Float_t oldfC44 = fC44;
162
163   if (TMath::Abs(fP4*Xnew - fP2) >= 0.9) // What's this??
164     {
165       //cout << "Propagation failed! Stiff track. " << pos << endl;
166       return 0;
167     }
168   
169   // R must be approx. of the radius of the circle based on
170   // the old and new spacepoint. What then is Cod and Cnew??
171   // C has something to do with curvature at least.
172
173   Float_t Cold = fP4*Xold - fP2; 
174   //if (TMath::Abs(Cold) >= 1.0) Cold = 0.9;
175     /*{
176     cout << "Cold " << endl << fP4*Xnew - fP2 << endl;
177     return 0;
178     }*/
179   Float_t Rold = TMath::Sqrt(1 - Cold*Cold);
180   Float_t Cnew = fP4*Xnew - fP2; 
181   //if (TMath::Abs(Cnew) >= 1.0) Cnew = 0.9;
182     /*{
183     cout << "Cnew " << endl;
184     return 0;
185     }*/
186   Float_t Rnew = TMath::Sqrt(1 - Cnew*Cnew);
187   //if (Rold < 0.9) cout << "Cold = " << Cold << " , Rold = " << Rold << " , Cnew = " << Cnew << " , Rnew = " << Rnew << endl;
188
189   // Prediction of the y- and z- coordinate in the next plane
190   fP0 += dx*(Cold+Cnew)/(Rold+Rnew);
191   fP1 += dx*(Cold+Cnew)/(Cold*Rnew + Cnew*Rold)*fP3; 
192   //cout << "Old point " << Yold << " " << Zold << endl;
193   //cout << "Propagate " << fP0 << " " << fP1 << endl;
194   //cout << "Measured  " << points[pos].fY << " " << points[pos].fZ << endl;
195
196   // f = F - 1 //What is this??
197   // Must be the f-matrix for the prediction, as in eq 1 in ALICE Kalman paper
198   Float_t RR = Rold + Rnew;
199   Float_t CC = Cold + Cnew;
200   Float_t XX = Xold + Xnew;
201
202   Float_t f02 = -dx*(2*RR + CC*(Cold/Rold + Cnew/Rnew))/(RR*RR);
203   Float_t f04 = dx*(RR*XX + CC*(Cold*Xold/Rold + Cnew*Xnew/Rnew))/(RR*RR);
204   Float_t CR = Cold*Rnew + Cnew*Rold;
205   Float_t f12 = -dx*fP3*(2*CR + CC*(Cnew*(Cold/Rold)-Rold + Cold*(Cnew/Rnew)-Rnew))/(CR*CR);
206   Float_t f13 = dx*CC/CR;
207   Float_t f14 = dx*fP3*(CR*XX-CC*(Rold*Xnew-Cnew*Cold*Xold/Rold + Rnew*Xold-Cold*Cnew*Xnew/Rnew))/(CR*CR);
208
209   // b = C*ft // This? 
210   Float_t b00=f02*fC20 + f04*fC40;
211   Float_t b01=f12*fC20 + f14*fC40 + f13*fC30;
212   Float_t b10=f02*fC21 + f04*fC41;
213   Float_t b11=f12*fC21 + f14*fC41 + f13*fC31;
214   Float_t b20=f02*fC22 + f04*fC42;
215   Float_t b21=f12*fC22 + f14*fC42 + f13*fC32;
216   Float_t b30=f02*fC32 + f04*fC43;
217   Float_t b31=f12*fC32 + f14*fC43 + f13*fC33;
218   Float_t b40=f02*fC42 + f04*fC44;
219   Float_t b41=f12*fC42 + f14*fC44 + f13*fC43;
220
221   //a = f*b = f*C*ft
222   Float_t a00 = f02*b20 + f04*b40;
223   Float_t a01 = f02*b21 + f04*b41;
224   Float_t a11 = f12*b21 + f14*b41+f13*b31;
225
226   //F*C*Ft = C + (a + b + bt) /This is the covariance matrix, the samll t 
227   // means transform. Then F must be df/dx 
228   fC00 += a00 + 2*b00;
229   fC10 += a01 + b01 + b10;
230   fC20 += b20;
231   fC30 += b30;
232   fC40 += b40;
233   fC11 += a11 + 2*b11;
234   fC21 += b21;
235   fC31 += b31;
236   fC41 += b41;
237
238   // Multiple scattering (from AliTPCtrack::PropagateTo)
239   // Shold this be included??
240   /*
241   Float_t d = TMath::Sqrt((Xold-Xnew)*(Xold-Xnew)+(Yold-fP0)*(Yold-fP0)+(Zold-fP1)*(Zold-fP1));
242   Float_t Pt = (1e-9*TMath::Abs(fP4)/fP4 + fP4 * (1000/0.299792458/4));
243   if (TMath::Abs(Pt) > 100) {
244     return 0;
245   }
246   if (TMath::Abs(Pt) < 0.01) return 0;
247   
248   Float_t p2 = (1+fP3*fP3)/(Pt*Pt);
249   Float_t beta2 = p2/(p2 +0.14*0.14);
250   Float_t theta2 = 1.0259e-6*10*10/20/(beta2*p2)*d*0.9e-3;
251   
252   Float_t ey=fP4*Xnew - fP2;
253   Float_t ez=fP3;
254   Float_t xz=fP4*ez;
255   Float_t zz1=ez*ez+1;
256   Float_t xy=fP2+ey;
257   
258   fC22 += (2*ey*ez*ez*fP2+1-ey*ey+ez*ez+fP2*fP2*ez*ez)*theta2;
259   fC32 += ez*zz1*xy*theta2;
260   fC33 += zz1*zz1*theta2;
261   fC42 += xz*ez*xy*theta2;
262   fC43 += xz*zz1*theta2;
263   fC44 += xz*xz*theta2;
264   if (TMath::Abs(beta2) >= 1) cout << "here? " << beta2 << ", " << Pt << endl;
265   
266   // Energy loss assuming track is pion
267   Float_t dE = 0.153e-3/beta2*(log(5940*beta2/(1-beta2))-beta2)*d*0.9e-3;
268   if (Xold < Xnew) dE = -dE;
269   CC = fP4;
270   fP4 *= (1 - TMath::Sqrt(p2+0.14*0.14)/p2*dE);
271   fP2 += Xnew*(fP4-CC);
272   */
273   // Update the track parameters with the measured values of the new point
274   Int_t value = UpdateTrack(points, pos);
275
276   /*if (value == 0)
277     {
278       fP0 = Yold;
279       fP1 = Zold;
280       fP2 = par2old;
281       fP3 = par3old;
282       fP4 = par4old;
283       fC00 = oldfC00;
284       fC10 = oldfC10;
285       fC11 = oldfC11;
286       fC20 = oldfC20;
287       fC21 = oldfC21;
288       fC22 = oldfC22;
289       fC30 = oldfC30;
290       fC31 = oldfC31;
291       fC32 = oldfC32;
292       fC33 = oldfC33;
293       fC40 = oldfC40;
294       fC41 = oldfC41;
295       fC42 = oldfC42;
296       fC43 = oldfC43;
297       fC44 = oldfC44;
298
299       return value;
300       }*/
301   
302   //else 
303 return value;
304   //return UpdateTrack(points, pos);
305 }
306
307 Int_t AliL3KalmanTrack::UpdateTrack(AliL3SpacePointData *points, UInt_t pos)
308 {
309   // Update the track parameters with the measured values
310   fX = points[pos].fX; 
311
312   // The errors from the measurement of the spacepoint
313   Float_t sigmaY2 = points[pos].fSigmaY2;
314   Float_t sigmaZ2 = points[pos].fSigmaZ2;
315   Float_t sigmaYZ = 0;
316
317   sigmaY2 += fC00;
318   sigmaZ2 += fC11;
319   sigmaYZ += fC10;
320
321   Float_t det = sigmaY2*sigmaZ2 - sigmaYZ*sigmaYZ;
322   Float_t tmp = sigmaY2;
323   sigmaY2 = sigmaZ2/det;
324   sigmaZ2 = tmp/det;
325   sigmaYZ = -sigmaYZ/det;
326
327   // Kalman gain matrix
328   Float_t k00 = fC00*sigmaY2 + fC10*sigmaYZ;
329   Float_t k01 = fC00*sigmaYZ + fC10*sigmaZ2;
330   Float_t k10 = fC10*sigmaY2 + fC11*sigmaYZ;
331   Float_t k11 = fC10*sigmaYZ + fC11*sigmaZ2;
332   Float_t k20 = fC20*sigmaY2 + fC21*sigmaYZ;
333   Float_t k21 = fC20*sigmaYZ + fC21*sigmaZ2;
334   Float_t k30 = fC30*sigmaY2 + fC31*sigmaYZ;
335   Float_t k31 = fC30*sigmaYZ + fC31*sigmaZ2;
336   Float_t k40 = fC40*sigmaY2 + fC41*sigmaYZ;
337   Float_t k41 = fC40*sigmaYZ + fC41*sigmaZ2;
338
339   // Deviation between the predicted and measured values of y and z 
340   Float_t dy = points[pos].fY-fP0; //cout << "dy = " << dy;
341   Float_t dz = points[pos].fZ-fP1; //cout << ", dz = " << dz << endl; 
342   //cout << "Measured " << points[pos].fY << " " << points[pos].fZ << endl;
343
344   // Prediction of fP2 and fP4
345   Float_t cur = fP4 + k40*dy + k41*dz; 
346   Float_t eta = fP2 + k20*dy + k21*dz;
347
348   if (TMath::Abs(cur*fX - eta) >= 0.9)
349     {
350       //cout << "Update failed! Stiff track. " << pos << endl;
351       return 0;
352     }
353
354   // Filtered state vector
355   fP0 += k00*dy + k01*dz;
356   fP1 += k10*dy + k11*dz;
357   fP2 = eta;
358   fP3 += k30*dy + k31*dz; //cout << "update " << fP3 << endl;
359   fP4 = cur;
360   //cout << "AliL3KalmanTrack::Update, " << fP0 << " " << fP1 << " " << fP2 << " " << fP3 << " " << fP4 << endl;
361   //cout << "Measured, " << points[pos].fY << " " << points[pos].fZ << endl;
362
363   Float_t c10 = fC10;
364   Float_t c20 = fC20;
365   Float_t c30 = fC30;
366   Float_t c40 = fC40;
367   Float_t c21 = fC21;
368   Float_t c31 = fC31;
369   Float_t c41 = fC41;
370
371   // Filtered covariance matrix
372   fC00 -= k00*fC00 + k01*fC10;
373   fC10 -= k00*c10 + k01*fC11;
374   fC11 -= k10*c10 + k11*fC11;
375   fC20 -= k00*c20 + k01*c21;
376   fC21 -= k10*c20 + k11*c21;
377   fC22 -= k20*c20 + k21*c21;
378   fC30 -= k00*c30 + k01*c31;
379   fC31 -= k10*c30 + k11*c31;
380   fC32 -= k20*c30 + k21*c31;
381   fC33 -= k30*c30 + k31*c31;
382   fC40 -= k00*c40 + k01*c41;
383   fC41 -= k10*c40 + k11*c41;
384   fC42 -= k20*c40 + k21*c41;
385   fC43 -= k40*c30 + k41*c31;
386   fC44 -= k40*c40 + k41*c41;
387   //cout << "AliL3KalmanTrack::Update, error " << fC00 << " " << fC11 << " " << fC22 << " " << fC33 << " " << fC44 << endl;
388
389   sigmaY2 = sigmaY2*det;
390   sigmaZ2 = sigmaZ2*det;
391   sigmaYZ = sigmaYZ*det;
392   //cout << "AliL3KalmanTrack::Update, Chi2 = " << GetChisq() << endl;
393   //cout << "AliKalmanTrack::Update, sigmaY2 = " << sigmaY2 << " sigmaZ2 = " << sigmaZ2 << " sigmaYZ = " << sigmaYZ << " dy = " << dy << " dz = " << dz << endl;
394   // Calculate increase of chisquare
395   fChisq = GetChisq() + GetChisqIncrement(points,pos);
396   //cout << "fChisq = " << fChisq << endl;
397 //(dy*sigmaY2*dy + 2*sigmaYZ*dy*dz + dz*sigmaZ2*dz) / (sigmaY2*sigmaZ2 - sigmaYZ*sigmaYZ);
398   // Must at some point make an cut on chisq. Here?
399   //  if (fChisq > fMaxChi2) return 0;
400
401   return 1;
402
403
404 Float_t AliL3KalmanTrack::GetChisqIncrement(AliL3SpacePointData *points, UInt_t pos)
405 {
406   // This function calculates a predicted chi2 increment.
407   //-----------------------------------------------------------------
408   Float_t r00=points[pos].fSigmaY2, r01=0., r11=points[pos].fSigmaZ2;
409   r00+=fC00; r01+=fC10; r11+=fC11;
410
411   Double_t det=r00*r11 - r01*r01;
412   /*if (TMath::Abs(det) < 1.e-10) {
413     Int_t n=GetNumberOfClusters();
414     if (n>4) cerr<<n<<" AliKalmanTrack warning: Singular matrix !\n";
415     return 1e10;
416     }*/
417   Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
418   
419   Double_t dy=points[pos].fY - fP0, dz=points[pos].fZ - fP1;
420   //cout << "AliTPCtrack::GetPredictedChi2, r00 = " << r00 << " r11 = " << r11 << " ro1 = " << r01 << " dy = " << dy << " dz = " << dz << endl;
421   return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
422 }
423
424 Float_t AliL3KalmanTrack::f2(Float_t x1,Float_t y1,
425                              Float_t x2,Float_t y2,
426                              Float_t x3,Float_t y3)
427 {
428   //-----------------------------------------------------------------
429   // Initial approximation of the track curvature times center of curvature
430   // Will be removed ??
431   //-----------------------------------------------------------------
432   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
433   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
434                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
435   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
436                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
437
438   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
439
440   return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
441 }
442
443 Float_t AliL3KalmanTrack::f3(Float_t x1,Float_t y1,
444                              Float_t x2,Float_t y2,
445                              Float_t z1,Float_t z2)
446 {
447   //-----------------------------------------------------------------
448   // Initial approximation of the tangent of the track dip angle
449   // Will be removed ??
450   //-----------------------------------------------------------------
451   return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
452 }
453
454 Float_t AliL3KalmanTrack::f4(Float_t x1,Float_t y1,
455                              Float_t x2,Float_t y2,
456                              Float_t x3,Float_t y3)
457 {
458   //-----------------------------------------------------------------
459   // Initial approximation of the track curvature
460   // Will be removed??
461   //-----------------------------------------------------------------
462   Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
463   Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
464                   (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
465   Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
466                   (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
467
468   Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
469
470   return -xr*yr/sqrt(xr*xr+yr*yr);
471 }
472
473 void AliL3KalmanTrack::Set(AliL3KalmanTrack *track)
474 {
475
476   AliL3KalmanTrack *tpt = (AliL3KalmanTrack*)track;
477   SetX0(tpt->GetX0());
478   SetX1(tpt->GetX1());
479   SetX2(tpt->GetX2());
480   SetX3(tpt->GetX3());
481   SetX4(tpt->GetX4());
482
483   SetC0(tpt->GetC0());
484   SetC1(tpt->GetC1());
485   SetC2(tpt->GetC2());
486   SetC3(tpt->GetC3());
487   SetC4(tpt->GetC4());
488   SetC5(tpt->GetC5());
489   SetC6(tpt->GetC6());
490   SetC7(tpt->GetC7());
491   SetC8(tpt->GetC8());
492   SetC9(tpt->GetC9());
493   SetC10(tpt->GetC10());
494   SetC11(tpt->GetC11());
495   SetC12(tpt->GetC12());
496   SetC13(tpt->GetC13());
497   SetC14(tpt->GetC14());
498
499   SetNHits(tpt->GetNHits());
500 }