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 }