GetClusterFast function implemented (No getter before) (Marian)
[u/mrichter/AliRoot.git] / TPC / AliTPCtrack.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 /* $Id$ */
17
18 //-----------------------------------------------------------------
19 //           Implementation of the TPC track class
20 //        This class is used by the AliTPCtracker class
21 //      Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch
22 //-----------------------------------------------------------------
23
24 #include <Riostream.h>
25
26 #include "AliTPCtrack.h"
27 #include "AliCluster.h"
28 #include "AliESDtrack.h"
29 #include "AliTPCReconstructor.h"
30
31 ClassImp(AliTPCtrack)
32
33 //_________________________________________________________________________
34 AliTPCtrack::AliTPCtrack(): 
35   AliKalmanTrack(),
36   fX(0),
37   fAlpha(0),
38   fdEdx(0),
39   fP0(0),
40   fP1(0),
41   fP2(0),
42   fP3(0),
43   fP4(0),
44   fC00(1e10),
45   fC10(0),
46   fC11(1e10),
47   fC20(0),
48   fC21(0),
49   fC22(1e10),
50   fC30(0),
51   fC31(0),
52   fC32(0),
53   fC33(1e10),
54   fC40(0),
55   fC41(0),
56   fC42(0),
57   fC43(0),
58   fC44(0),
59   fNumber(0),
60   fSdEdx(1e10),
61   fNFoundable(0),
62   fBConstrain(kFALSE),
63   fLastPoint(-1),
64   fFirstPoint(-1),
65   fRemoval(0),
66   fTrackType(0),
67   fLab2(-1),
68   fNShared(0),
69   fReference()
70 {
71   //-------------------------------------------------
72   // default constructor
73   //-------------------------------------------------
74   for (Int_t i=0; i<kMaxRow;i++) fIndex[i]=-2;
75   for (Int_t i=0; i<4;i++) fPoints[i]=0.;
76   for (Int_t i=0; i<12;i++) fKinkPoint[i]=0.;
77   for (Int_t i=0; i<3;i++) fKinkIndexes[i]=0;
78   for (Int_t i=0; i<3;i++) fV0Indexes[i]=0;
79 }
80
81 //_________________________________________________________________________
82
83
84
85 AliTPCtrack::AliTPCtrack(UInt_t index, const Double_t xx[5],
86 const Double_t cc[15], Double_t xref, Double_t alpha) :
87   AliKalmanTrack(),
88   fX(xref),
89   fdEdx(0),
90   fP0(xx[0]),
91   fP1(xx[1]),
92   fP2(xx[2]),
93   fP3(xx[3]),
94   fP4(xx[4]),
95   fC00(cc[0]),
96   fC10(cc[1]),
97   fC11(cc[2]),
98   fC20(cc[3]),
99   fC21(cc[4]),
100   fC22(cc[5]),
101   fC30(cc[6]),
102   fC31(cc[7]),
103   fC32(cc[8]),
104   fC33(cc[9]),
105   fC40(cc[10]),
106   fC41(cc[11]),
107   fC42(cc[12]),
108   fC43(cc[13]),
109   fC44(cc[14]),
110   fNumber(1),
111   fSdEdx(1e10),
112   fNFoundable(0),
113   fBConstrain(kFALSE),
114   fLastPoint(0),
115   fFirstPoint(0),
116   fRemoval(0),
117   fTrackType(0),
118   fLab2(0),
119   fNShared(0),
120   fReference()
121 {
122   //-----------------------------------------------------------------
123   // This is the main track constructor.
124   //-----------------------------------------------------------------
125   fAlpha=alpha;
126   if (fAlpha<-TMath::Pi()) fAlpha += 2*TMath::Pi();
127   if (fAlpha>=TMath::Pi()) fAlpha -= 2*TMath::Pi();
128
129   SaveLocalConvConst();
130
131   SetNumberOfClusters(1);
132   
133   fIndex[0]=index;
134   for (Int_t i=1; i<kMaxRow;i++) fIndex[i]=-2;
135   for (Int_t i=0; i<4;i++) fPoints[i]=0.;
136   for (Int_t i=0; i<12;i++) fKinkPoint[i]=0.;
137   for (Int_t i=0; i<3;i++) fKinkIndexes[i]=0;
138   for (Int_t i=0; i<3;i++) fV0Indexes[i]=0;
139 }
140
141 //_____________________________________________________________________________
142 AliTPCtrack::AliTPCtrack(const AliESDtrack& t) :
143   AliKalmanTrack(),
144   fSdEdx(1e10),
145   fNFoundable(0),
146   fBConstrain(kFALSE),
147   fLastPoint(0),
148   fFirstPoint(0),
149   fRemoval(0),
150   fTrackType(0),
151   fLab2(0),
152   fNShared(0),
153   fReference()
154 {
155   //-----------------------------------------------------------------
156   // Conversion AliESDtrack -> AliTPCtrack.
157   //-----------------------------------------------------------------
158   SetNumberOfClusters(t.GetTPCclusters(fIndex));
159   SetLabel(t.GetLabel());
160   SetMass(t.GetMass());
161   for (Int_t i=0; i<4;i++) fPoints[i]=0.;
162   for (Int_t i=0; i<12;i++) fKinkPoint[i]=0.;
163   for (Int_t i=0; i<3;i++) fKinkIndexes[i]=0;
164   for (Int_t i=0; i<3;i++) fV0Indexes[i]=0;
165
166   fdEdx  = t.GetTPCsignal();
167   fAlpha = t.GetAlpha();
168   if      (fAlpha < -TMath::Pi()) fAlpha += 2*TMath::Pi();
169   else if (fAlpha >= TMath::Pi()) fAlpha -= 2*TMath::Pi();
170
171   //Conversion of the track parameters
172   Double_t x,p[5]; t.GetExternalParameters(x,p);
173   Double_t c[15]; t.GetExternalCovariance(c);
174
175   fX=x;    
176   fP0=p[0];
177   fP1=p[1];    SaveLocalConvConst();
178   fP3=p[3];    x=GetLocalConvConst();
179   fP4=p[4]/x;
180   fP2=fP4*fX - p[2];
181
182   //Conversion of the covariance matrix
183   c[10]/=x; c[11]/=x; c[12]/=x; c[13]/=x; c[14]/=x*x;
184
185   Double_t c22=fX*fX*c[14] - 2*fX*c[12] + c[5];
186   Double_t c32=fX*c[13] - c[8];
187   Double_t c20=fX*c[10] - c[3], c21=fX*c[11] - c[4], c42=fX*c[14] - c[12];
188
189   fC00=c[0 ];
190   fC10=c[1 ];   fC11=c[2 ];
191   fC20=c20;     fC21=c21;     fC22=c22;
192   fC30=c[6 ];   fC31=c[7 ];   fC32=c32;   fC33=c[9 ];
193   fC40=c[10];   fC41=c[11];   fC42=c42;   fC43=c[13]; fC44=c[14];
194
195   if ((t.GetStatus()&AliESDtrack::kTIME) == 0) return;
196   StartTimeIntegral();
197   Double_t times[10]; t.GetIntegratedTimes(times); SetIntegratedTimes(times);
198   SetIntegratedLength(t.GetIntegratedLength());
199 }
200
201 //_____________________________________________________________________________
202 AliTPCtrack::AliTPCtrack(const AliTPCtrack& t) :
203   AliKalmanTrack(t),
204   fX(t.fX),
205   fAlpha(t.fAlpha),
206   fdEdx(t.fdEdx),
207   fP0(t.fP0),
208   fP1(t.fP1),
209   fP2(t.fP2),
210   fP3(t.fP3),
211   fP4(t.fP4),
212   fC00(t.fC00),
213   fC10(t.fC10),
214   fC11(t.fC11),
215   fC20(t.fC20),
216   fC21(t.fC21),
217   fC22(t.fC22),
218   fC30(t.fC30),
219   fC31(t.fC31),
220   fC32(t.fC32),
221   fC33(t.fC33),
222   fC40(t.fC40),
223   fC41(t.fC41),
224   fC42(t.fC42),
225   fC43(t.fC43),
226   fC44(t.fC44),
227   fNumber(t.fNumber),
228   fSdEdx(t.fSdEdx),
229   fNFoundable(t.fNFoundable),
230   fBConstrain(t.fBConstrain),
231   fLastPoint(t.fLastPoint),
232   fFirstPoint(t.fFirstPoint),
233   fRemoval(t.fRemoval),
234   fTrackType(t.fTrackType),
235   fLab2(t.fLab2),
236   fNShared(t.fNShared),
237   fReference(t.fReference)
238
239 {
240   //-----------------------------------------------------------------
241   // This is a track copy constructor.
242   //-----------------------------------------------------------------
243   for (Int_t i=0; i<kMaxRow; i++) fIndex[i]=t.fIndex[i];
244   for (Int_t i=0; i<4;i++) fPoints[i]=t.fPoints[i];
245   for (Int_t i=0; i<12;i++) fKinkPoint[i]=t.fKinkPoint[i];
246   for (Int_t i=0; i<3;i++) fKinkIndexes[i]=t.fKinkIndexes[i];
247   for (Int_t i=0; i<3;i++) fV0Indexes[i]=t.fV0Indexes[i];
248 }
249
250 //_____________________________________________________________________________
251 Int_t AliTPCtrack::Compare(const TObject *o) const {
252   //-----------------------------------------------------------------
253   // This function compares tracks according to the their curvature
254   //-----------------------------------------------------------------
255   AliTPCtrack *t=(AliTPCtrack*)o;
256   //Double_t co=TMath::Abs(t->Get1Pt());
257   //Double_t c =TMath::Abs(Get1Pt());
258   Double_t co=t->GetSigmaY2()*t->GetSigmaZ2();
259   Double_t c =GetSigmaY2()*GetSigmaZ2();
260   if (c>co) return 1;
261   else if (c<co) return -1;
262   return 0;
263 }
264
265 //_____________________________________________________________________________
266 void AliTPCtrack::GetExternalCovariance(Double_t cc[15]) const {
267  //-------------------------------------------------------------------------
268   // This function returns an external representation of the covriance matrix.
269   //   (See comments in AliTPCtrack.h about external track representation)
270   //-------------------------------------------------------------------------
271   Double_t a=GetLocalConvConst();
272
273   Double_t c22=fX*fX*fC44-2*fX*fC42+fC22;
274   Double_t c32=fX*fC43-fC32;
275   Double_t c20=fX*fC40-fC20, c21=fX*fC41-fC21, c42=fX*fC44-fC42;
276   
277   cc[0 ]=fC00;
278   cc[1 ]=fC10;   cc[2 ]=fC11;
279   cc[3 ]=c20;    cc[4 ]=c21;    cc[5 ]=c22;
280   cc[6 ]=fC30;   cc[7 ]=fC31;   cc[8 ]=c32;   cc[9 ]=fC33; 
281   cc[10]=fC40*a; cc[11]=fC41*a; cc[12]=c42*a; cc[13]=fC43*a; cc[14]=fC44*a*a;
282
283 }
284
285 //_____________________________________________________________________________
286 Double_t AliTPCtrack::GetPredictedChi2(const AliCluster *c) const 
287 {
288   //-----------------------------------------------------------------
289   // This function calculates a predicted chi2 increment.
290   //-----------------------------------------------------------------
291   Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
292   r00+=fC00; r01+=fC10; r11+=fC11;
293
294   Double_t det=r00*r11 - r01*r01;
295   if (TMath::Abs(det) < 1.e-10) {
296     Int_t n=GetNumberOfClusters();
297     if (n>4) cerr<<n<<" AliKalmanTrack warning: Singular matrix !\n";
298     return 1e10;
299   }
300   Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
301   
302   Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
303   
304   return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
305 }
306
307 Double_t AliTPCtrack::GetYat(Double_t xk) const {
308 //-----------------------------------------------------------------
309 // This function calculates the Y-coordinate of a track at the plane x=xk.
310 //-----------------------------------------------------------------
311   if (TMath::Abs(fP4*fX - fP2)>AliTPCReconstructor::GetMaxSnpTrack()) return 0.; //patch 01 jan 06
312     Double_t c1=fP4*fX - fP2, r1=TMath::Sqrt(1.- c1*c1);
313     Double_t c2=fP4*xk - fP2;
314     if (c2*c2>AliTPCReconstructor::GetMaxSnpTrack()) {
315       //       Int_t n=GetNumberOfClusters();
316       // if (n>4) cerr<<n<<"AliTPCtrack::GetYat: can't evaluate the y-coord !\n";
317        return 0;
318     } 
319     Double_t r2=TMath::Sqrt(1.- c2*c2);
320     return fP0 + (xk-fX)*(c1+c2)/(r1+r2);
321 }
322
323 //_____________________________________________________________________________
324 Int_t AliTPCtrack::PropagateTo(Double_t xk,Double_t /*x0*/,Double_t rho) {
325   //-----------------------------------------------------------------
326   // This function propagates a track to a reference plane x=xk.
327   //-----------------------------------------------------------------
328   if (TMath::Abs(fP4*xk - fP2) >= AliTPCReconstructor::GetMaxSnpTrack()) {
329     //    Int_t n=GetNumberOfClusters();
330     //if (n>4) cerr<<n<<" AliTPCtrack warning: Propagation failed !\n";
331     return 0;
332   }
333   Double_t lcc=GetLocalConvConst();  
334
335   // old position for time [SR, GSI 17.02.2003]
336   Double_t oldX = fX;
337   Double_t oldY = fP0;
338   Double_t oldZ = fP1;
339   //
340
341   Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1, y1=fP0, z1=fP1;
342   Double_t c1=fP4*x1 - fP2, r1=sqrt(1.- c1*c1);
343   Double_t c2=fP4*x2 - fP2, r2=sqrt(1.- c2*c2);
344   
345   fP0 += dx*(c1+c2)/(r1+r2);
346   fP1 += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3;
347
348   //f = F - 1
349   Double_t rr=r1+r2, cc=c1+c2, xx=x1+x2;
350   Double_t f02=-dx*(2*rr + cc*(c1/r1 + c2/r2))/(rr*rr);
351   Double_t f04= dx*(rr*xx + cc*(c1*x1/r1+c2*x2/r2))/(rr*rr);
352   Double_t cr=c1*r2+c2*r1;
353   Double_t f12=-dx*fP3*(2*cr + cc*(c2*c1/r1-r1 + c1*c2/r2-r2))/(cr*cr);
354   Double_t f13= dx*cc/cr; 
355   Double_t f14=dx*fP3*(cr*xx-cc*(r1*x2-c2*c1*x1/r1+r2*x1-c1*c2*x2/r2))/(cr*cr);
356
357   //b = C*ft
358   Double_t b00=f02*fC20 + f04*fC40, b01=f12*fC20 + f14*fC40 + f13*fC30;
359   Double_t b10=f02*fC21 + f04*fC41, b11=f12*fC21 + f14*fC41 + f13*fC31;
360   Double_t b20=f02*fC22 + f04*fC42, b21=f12*fC22 + f14*fC42 + f13*fC32;
361   Double_t b30=f02*fC32 + f04*fC43, b31=f12*fC32 + f14*fC43 + f13*fC33;
362   Double_t b40=f02*fC42 + f04*fC44, b41=f12*fC42 + f14*fC44 + f13*fC43;
363   
364   //a = f*b = f*C*ft
365   Double_t a00=f02*b20+f04*b40,a01=f02*b21+f04*b41,a11=f12*b21+f14*b41+f13*b31;
366
367   //F*C*Ft = C + (a + b + bt)
368   fC00 += a00 + 2*b00;
369   fC10 += a01 + b01 + b10; 
370   fC20 += b20;
371   fC30 += b30;
372   fC40 += b40;
373   fC11 += a11 + 2*b11;
374   fC21 += b21; 
375   fC31 += b31; 
376   fC41 += b41; 
377
378   fX=x2;
379
380   //Change of the magnetic field *************
381   SaveLocalConvConst();
382   cc=fP4;
383   fP4*=lcc/GetLocalConvConst();
384   fP2+=fX*(fP4-cc);
385
386   //Multiple scattering ******************
387   Double_t d=sqrt((x1-fX)*(x1-fX)+(y1-fP0)*(y1-fP0)+(z1-fP1)*(z1-fP1));
388   Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
389   Double_t beta2=p2/(p2 + GetMass()*GetMass());
390   beta2 = TMath::Min(beta2,0.99999999999);
391   //Double_t theta2=14.1*14.1/(beta2*p2*1e6)*d/x0*rho;
392   Double_t theta2=1.0259e-6*10*10/20/(beta2*p2)*d*rho;
393
394   Double_t ey=fP4*fX - fP2, ez=fP3;
395   Double_t xz=fP4*ez, zz1=ez*ez+1, xy=fP2+ey;
396     
397   fC22 += (2*ey*ez*ez*fP2+1-ey*ey+ez*ez+fP2*fP2*ez*ez)*theta2;
398   fC32 += ez*zz1*xy*theta2;
399   fC33 += zz1*zz1*theta2;
400   fC42 += xz*ez*xy*theta2;
401   fC43 += xz*zz1*theta2;
402   fC44 += xz*xz*theta2;
403   /*
404   //
405   //MI coeficients
406   Double_t dc22 = (1-ey*ey+xz*xz*fX*fX)*theta2;
407   Double_t dc32 = (xz*fX*zz1)*theta2;
408   Double_t dc33 = (zz1*zz1)*theta2;
409   Double_t dc42 = (xz*fX*xz)*theta2;
410   Double_t dc43 = (zz1*xz)*theta2;
411   Double_t dc44 = (xz*xz)*theta2; 
412   fC22 += dc22;
413   fC32 += dc32;
414   fC33 += dc33;
415   fC42 += dc42;
416   fC43 += dc43;
417   fC44 += dc44;
418   */
419   //Energy losses ************************
420   Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2)) - beta2)*d*rho;
421   if (x1 < x2) dE=-dE;
422   cc=fP4;
423
424   //Double_t E = sqrt(p2+GetMass()*GetMass());
425   //Double_t mifac  = TMath::Sqrt(1.+dE*dE/p2+2*E*dE/p2)-1;
426   //Double_t belfac = E*dE/p2;
427                                //
428   fP4*=(1.- sqrt(p2+GetMass()*GetMass())/p2*dE);
429   fP2+=fX*(fP4-cc);
430
431   // Integrated Time [SR, GSI, 17.02.2003]
432  if (x1 < x2)
433  if (IsStartedTimeIntegral()) {
434     Double_t l2 = (fX-oldX)*(fX-oldX)+(fP0-oldY)*(fP0-oldY)+(fP1-oldZ)*(fP1-oldZ);
435     AddTimeStep(TMath::Sqrt(l2));
436   }
437   //
438
439   return 1;
440 }
441
442 //_____________________________________________________________________________
443 Int_t AliTPCtrack::PropagateToVertex(Double_t x0,Double_t rho) 
444 {
445   //-----------------------------------------------------------------
446   // This function propagates tracks to the "vertex".
447   //-----------------------------------------------------------------
448   Double_t c=fP4*fX - fP2;
449   Double_t tgf=-fP2/(fP4*fP0 + sqrt(1-c*c));
450   Double_t snf=tgf/sqrt(1.+ tgf*tgf);
451   Double_t xv=(fP2+snf)/fP4;
452   return PropagateTo(xv,x0,rho);
453 }
454
455 //_____________________________________________________________________________
456 Int_t AliTPCtrack::Update(const AliCluster *c, Double_t chisq, UInt_t index) {
457   //-----------------------------------------------------------------
458   // This function associates a cluster with this track.
459   //-----------------------------------------------------------------
460   Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
461   r00+=fC00; r01+=fC10; r11+=fC11;
462   Double_t det=r00*r11 - r01*r01;
463   Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
464
465   Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
466   Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
467   Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
468   Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
469   Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
470
471   Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
472   Double_t cur=fP4 + k40*dy + k41*dz, eta=fP2 + k20*dy + k21*dz;
473   if (TMath::Abs(cur*fX-eta) >= AliTPCReconstructor::GetMaxSnpTrack()) {
474     //    Int_t n=GetNumberOfClusters();
475     //if (n>4) cerr<<n<<" AliTPCtrack warning: Filtering failed !\n";
476     return 0;
477   }
478
479   fP0 += k00*dy + k01*dz;
480   fP1 += k10*dy + k11*dz;
481   fP2  = eta;
482   fP3 += k30*dy + k31*dz;
483   fP4  = cur;
484
485   Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
486   Double_t c12=fC21, c13=fC31, c14=fC41;
487
488   fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
489   fC20-=k00*c02+k01*c12;   fC30-=k00*c03+k01*c13;
490   fC40-=k00*c04+k01*c14; 
491
492   fC11-=k10*c01+k11*fC11;
493   fC21-=k10*c02+k11*c12;   fC31-=k10*c03+k11*c13;
494   fC41-=k10*c04+k11*c14; 
495
496   fC22-=k20*c02+k21*c12;   fC32-=k20*c03+k21*c13;
497   fC42-=k20*c04+k21*c14; 
498
499   fC33-=k30*c03+k31*c13;
500   fC43-=k40*c03+k41*c13; 
501
502   fC44-=k40*c04+k41*c14; 
503
504   Int_t n=GetNumberOfClusters();
505   fIndex[n]=index;
506   SetNumberOfClusters(n+1);
507   SetChi2(GetChi2()+chisq);
508
509   return 1;
510 }
511
512 //_____________________________________________________________________________
513 Int_t AliTPCtrack::Rotate(Double_t alpha)
514 {
515   //-----------------------------------------------------------------
516   // This function rotates this track.
517   //-----------------------------------------------------------------
518   fAlpha += alpha;
519   if (fAlpha<-TMath::Pi()) fAlpha += 2*TMath::Pi();
520   if (fAlpha>=TMath::Pi()) fAlpha -= 2*TMath::Pi();
521   
522   Double_t x1=fX, y1=fP0;
523   Double_t ca=cos(alpha), sa=sin(alpha);
524   Double_t r1=fP4*fX - fP2;
525   
526   if (TMath::Abs(r1)>=AliTPCReconstructor::GetMaxSnpTrack()) return 0; //patch 01 jan 06
527
528   fX = x1*ca + y1*sa;
529   fP0=-x1*sa + y1*ca;
530   fP2=fP2*ca + (fP4*y1 + sqrt(1.- r1*r1))*sa;
531   
532   Double_t r2=fP4*fX - fP2;
533   if (TMath::Abs(r2) >= AliTPCReconstructor::GetMaxSnpTrack()) {
534     //Int_t n=GetNumberOfClusters();
535     //    if (n>4) cerr<<n<<" AliTPCtrack warning: Rotation failed !\n";
536     return 0;
537   }
538   
539   Double_t y0=fP0 + sqrt(1.- r2*r2)/fP4;
540   if ((fP0-y0)*fP4 >= 0.) {
541     //Int_t n=GetNumberOfClusters();
542     //    if (n>4) cerr<<n<<" AliTPCtrack warning: Rotation failed !!!\n";
543     return 0;
544   }
545
546   //f = F - 1
547   Double_t f00=ca-1,    f24=(y1 - r1*x1/sqrt(1.- r1*r1))*sa, 
548            f20=fP4*sa,  f22=(ca + sa*r1/sqrt(1.- r1*r1))-1;
549
550   //b = C*ft
551   Double_t b00=fC00*f00, b02=fC00*f20+fC40*f24+fC20*f22;
552   Double_t b10=fC10*f00, b12=fC10*f20+fC41*f24+fC21*f22;
553   Double_t b20=fC20*f00, b22=fC20*f20+fC42*f24+fC22*f22;
554   Double_t b30=fC30*f00, b32=fC30*f20+fC43*f24+fC32*f22;
555   Double_t b40=fC40*f00, b42=fC40*f20+fC44*f24+fC42*f22;
556
557   //a = f*b = f*C*ft
558   Double_t a00=f00*b00, a02=f00*b02, a22=f20*b02+f24*b42+f22*b22;
559
560   // *** Double_t dy2=fCyy;
561
562   //F*C*Ft = C + (a + b + bt)
563   fC00 += a00 + 2*b00;
564   fC10 += b10;
565   fC20 += a02+b20+b02;
566   fC30 += b30;
567   fC40 += b40;
568   fC21 += b12;
569   fC32 += b32;
570   fC22 += a22 + 2*b22;
571   fC42 += b42; 
572
573   // *** fCyy+=dy2*sa*sa*r1*r1/(1.- r1*r1);
574   // *** fCzz+=d2y*sa*sa*fT*fT/(1.- r1*r1);
575
576   return 1;
577 }
578
579 void AliTPCtrack::ResetCovariance() {
580   //------------------------------------------------------------------
581   //This function makes a track forget its history :)  
582   //------------------------------------------------------------------
583
584   fC00*=10.;
585   fC10=0.;  fC11*=10.;
586   fC20=0.;  fC21=0.;  fC22*=10.;
587   fC30=0.;  fC31=0.;  fC32=0.;  fC33*=10.;
588   fC40=0.;  fC41=0.;  fC42=0.;  fC43=0.;  fC44*=10.;
589
590 }
591
592 ////////////////////////////////////////////////////////////////////////
593 Double_t AliTPCtrack::Phi() const {
594 //
595 //
596 //
597   Double_t phi =  TMath::ASin(GetSnp()) + fAlpha;
598   if (phi<0) phi+=2*TMath::Pi();
599   if (phi>=2*TMath::Pi()) phi-=2*TMath::Pi();
600   return phi;
601 }
602 ////////////////////////////////////////////////////////////////////////
603
604
605
606 ////////////////////////////////////////////////////////////////////////
607 // MI ADDITION
608
609 Float_t AliTPCtrack::Density(Int_t row0, Int_t row1)
610 {
611   //
612   // calculate cluster density
613   Int_t good  = 0;
614   Int_t found = 0;
615   //if (row0<fFirstPoint) row0 = fFirstPoint;
616   if (row1>fLastPoint) row1 = fLastPoint;
617
618   
619   for (Int_t i=row0;i<=row1;i++){ 
620     //    Int_t index = fClusterIndex[i];
621     Int_t index = fIndex[i];
622     if (index!=-1)  good++;
623     if (index>0)    found++;
624   }
625   Float_t density=0;
626   if (good>0) density = Float_t(found)/Float_t(good);
627   return density;
628 }
629
630
631 Float_t AliTPCtrack::Density2(Int_t row0, Int_t row1)
632 {
633   //
634   // calculate cluster density
635   Int_t good  = 0;
636   Int_t found = 0;
637   //  
638   for (Int_t i=row0;i<=row1;i++){     
639     Int_t index = fIndex[i];
640     if (index!=-1)  good++;
641     if (index>0)    found++;
642   }
643   Float_t density=0;
644   if (good>0) density = Float_t(found)/Float_t(good);
645   return density;
646 }
647
648
649 Double_t AliTPCtrack::GetZat0() const
650 {
651   //
652   // return virtual z - supposing that x = 0
653   if (TMath::Abs(fP2)>1) return 0;
654   if (TMath::Abs(fX*fP4-fP2)>1) return 0;
655   Double_t vz = fP1+fP3/fP4*(asin(-fP2)-asin(fX*fP4-fP2));
656   return vz;
657 }
658
659
660 Double_t AliTPCtrack::GetD(Double_t x, Double_t y) const {
661   //------------------------------------------------------------------
662   // This function calculates the transverse impact parameter
663   // with respect to a point with global coordinates (x,y)
664   //------------------------------------------------------------------
665   //Double_t xt=fX, yt=fP0;
666
667   Double_t sn=TMath::Sin(fAlpha), cs=TMath::Cos(fAlpha);
668   Double_t a = x*cs + y*sn;
669   y = -x*sn + y*cs; x=a;
670   //
671   Double_t r  = TMath::Abs(1/fP4);
672   Double_t x0 = TMath::Abs(fP2*r);
673   Double_t y0 = fP0;
674   y0= fP0+TMath::Sqrt(1-(fP4*fX-fP2)*(fP4*fX-fP2))/fP4;
675   
676   Double_t  delta = TMath::Sqrt((x-x0)*(x-x0)+(y-y0)*(y-y0));  
677   //  Double_t  delta = TMath::Sqrt(TMath::Abs(x*x-2*x0*x+x0*x0+ y*y-2*y*y0+y0*y0));
678   delta -= TMath::Abs(r);
679   return delta;  
680 }
681
682 //
683 //
684
685 void  AliTPCtrack::UpdatePoints()
686 {
687   //--------------------------------------------------
688   //calculates first ,amx dens and last points
689   //--------------------------------------------------
690   Float_t density[160];
691   for (Int_t i=0;i<160;i++) density[i]=-1.;
692   fPoints[0]= 160;
693   fPoints[1] = -1;
694   //
695   Int_t ngood=0;
696   Int_t undeff=0;
697   Int_t nall =0;
698   Int_t range=20;
699   for (Int_t i=0;i<160;i++){
700     Int_t last = i-range;
701     if (nall<range) nall++;
702     if (last>=0){
703       if (fIndex[last]>0&& (fIndex[last]&0x8000)==0) ngood--;
704       if (fIndex[last]==-1) undeff--;
705     }
706     if (fIndex[i]>0&& (fIndex[i]&0x8000)==0)   ngood++;
707     if (fIndex[i]==-1) undeff++;
708     if (nall==range &&undeff<range/2) density[i-range/2] = Float_t(ngood)/Float_t(nall-undeff);
709   }
710   Float_t maxdens=0;
711   Int_t indexmax =0;
712   for (Int_t i=0;i<160;i++){
713     if (density[i]<0) continue;
714     if (density[i]>maxdens){
715       maxdens=density[i];
716       indexmax=i;
717     }
718   }
719   //
720   //max dens point
721   fPoints[3] = maxdens;
722   fPoints[1] = indexmax;
723   //
724   // last point
725   for (Int_t i=indexmax;i<160;i++){
726     if (density[i]<0) continue;
727     if (density[i]<maxdens/2.) {
728       break;
729     }
730     fPoints[2]=i;
731   }
732   //
733   // first point
734   for (Int_t i=indexmax;i>0;i--){
735     if (density[i]<0) continue;
736     if (density[i]<maxdens/2.) {
737       break;
738     }
739     fPoints[0]=i;
740   }
741   //
742 }
743
744 Double_t AliTPCtrack::Get1Pt() const {
745   //--------------------------------------------------------------
746   // Returns the inverse Pt (1/GeV/c)
747   // (or 1/"most probable pt", if the field is too weak)
748   //--------------------------------------------------------------
749   if (TMath::Abs(GetLocalConvConst()) > kVeryBigConvConst)
750       return 1./kMostProbableMomentum/TMath::Sqrt(1.+ GetTgl()*GetTgl());
751   return (TMath::Sign(1e-9,fP4) + fP4)*GetLocalConvConst();
752 }