]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TRD/AliTRDtrack.cxx
Additional protection and change in the interface (Yu.Belikov)
[u/mrichter/AliRoot.git] / TRD / AliTRDtrack.cxx
index 133f1382ef164db55c37db11663b8679ca9a2428..868da8bf6cfee2f2ebe79196df27afe802cd097d 100644 (file)
 #include "AliTRDgeometry.h" 
 #include "AliTRDcluster.h" 
 #include "AliTRDtrack.h"
-#include "AliESDtrack.h" 
 #include "AliTRDclusterCorrection.h"
+#include "AliTrackReference.h"
 
+ClassImp(AliTRDtracklet)
 ClassImp(AliTRDtrack)
 
+
+  AliTRDtracklet::AliTRDtracklet():fY(0),fZ(0),fX(0),fAlpha(0),fSigma2(0),fP0(0),fP1(0),fNFound(0),fNCross(0),fPlane(0),fExpectedSigma2(0),fChi2(0),fTilt(0),fMaxPos(0),fMaxPos4(0),fMaxPos5(0){
+}
+
 //_____________________________________________________________________________
 
 AliTRDtrack::AliTRDtrack(const AliTRDcluster *c, UInt_t index, 
@@ -44,6 +49,8 @@ AliTRDtrack::AliTRDtrack(const AliTRDcluster *c, UInt_t index,
   fX=xref;
 
   fY=xx[0]; fZ=xx[1]; fE=xx[2]; fT=xx[3]; fC=xx[4];
+  SaveLocalConvConst();
 
   fCyy=cc[0];
   fCzy=cc[1];  fCzz=cc[2];
@@ -55,11 +62,23 @@ AliTRDtrack::AliTRDtrack(const AliTRDcluster *c, UInt_t index,
   SetNumberOfClusters(1);
 
   fdEdx=0.;
+  fdEdxT=0.;
+  fDE=0.;
+  for (Int_t i=0;i<kNPlane;i++){
+      fdEdxPlane[i] = 0.;
+      fTimBinPlane[i] = -1;
+  }
 
   fLhElectron = 0.0;
   fNWrong = 0;
   fNRotate = 0;
   fStopped = 0;
+  fNCross =0;
+  fNLast  =0;
+  fChi2Last=0;
+  fNExpected=0;
+  fNExpectedLast=0;
+  fNdedx=0;
   Double_t q = TMath::Abs(c->GetQ());
   Double_t s = fX*fC - fE, t=fT;
   if(s*s < 1) q *= TMath::Sqrt((1-s*s)/(1+t*t));
@@ -72,6 +91,8 @@ AliTRDtrack::AliTRDtrack(const AliTRDcluster *c, UInt_t index,
     fIndex[i] = 0;
     fIndexBackup[i] = 0;  //bacup indexes MI    
   }
+  for (Int_t i=0;i<3;i++) { fBudget[i]=0;};
+  fBackupTrack =0;  
 }                              
            
 //_____________________________________________________________________________
@@ -85,14 +106,29 @@ AliTRDtrack::AliTRDtrack(const AliTRDtrack& t) : AliKalmanTrack(t) {
 
   SetChi2(t.GetChi2());
   fdEdx=t.fdEdx;
+  fdEdxT=t.fdEdxT;
+  fDE=t.fDE;
+  for (Int_t i=0;i<kNPlane;i++){
+      fdEdxPlane[i] = t.fdEdxPlane[i];
+      fTimBinPlane[i] = t.fTimBinPlane[i];
+      fTracklets[i]   = t.fTracklets[i];
+  }
 
   fLhElectron = 0.0;
   fNWrong = t.fNWrong;
   fNRotate = t.fNRotate;
   fStopped = t.fStopped;
+  fNCross  = t.fNCross;
+  fNExpected = t.fNExpected;
+  fNExpectedLast = t.fNExpectedLast;
+  fNdedx         = t.fNdedx;
+  fNLast     = t.fNLast;
+  fChi2Last  = t.fChi2Last;
+  fBackupTrack =0;
   fAlpha=t.fAlpha;
   fX=t.fX;
 
+
   fY=t.fY; fZ=t.fZ; fE=t.fE; fT=t.fT; fC=t.fC;
 
   fCyy=t.fCyy;
@@ -115,6 +151,10 @@ AliTRDtrack::AliTRDtrack(const AliTRDtrack& t) : AliKalmanTrack(t) {
     fIndex[i] = 0;
     fIndexBackup[i] = 0;  //MI backup indexes
   }
+  for (Int_t i=0;i<6;i++){
+    fTracklets[i] = t.fTracklets[i];
+  }
+  for (Int_t i=0;i<3;i++) { fBudget[i]=t.fBudget[i];};
 }                                
 
 //_____________________________________________________________________________
@@ -129,12 +169,23 @@ AliTRDtrack::AliTRDtrack(const AliKalmanTrack& t, Double_t alpha)
   SetMass(t.GetMass());
   SetNumberOfClusters(0);
 
-  fdEdx=t.GetdEdx();
+  fdEdx=t.GetPIDsignal();
+  for (Int_t i=0;i<kNPlane;i++){
+    fdEdxPlane[i] = 0.0;
+    fTimBinPlane[i] = -1;
+  }
 
   fLhElectron = 0.0;
   fNWrong = 0;
   fNRotate = 0;
   fStopped = 0;
+  fNExpected=0;
+  fNExpectedLast=0;  
+  fNdedx        =0;
+  fNCross =0;
+  fNLast  =0;
+  fChi2Last =0;
+  fBackupTrack =0;
 
   fAlpha = alpha;
   if      (fAlpha < -TMath::Pi()) fAlpha += 2*TMath::Pi();
@@ -144,11 +195,9 @@ AliTRDtrack::AliTRDtrack(const AliKalmanTrack& t, Double_t alpha)
 
   fX=x;
 
-  x = GetConvConst();  
-
   fY=p[0];
   fZ=p[1];
-  fT=p[3];
+  fT=p[3]; x=GetLocalConvConst();
   fC=p[4]/x;
   fE=fC*fX - p[2];   
 
@@ -173,6 +222,8 @@ AliTRDtrack::AliTRDtrack(const AliKalmanTrack& t, Double_t alpha)
     fIndex[i] = 0;
     fIndexBackup[i] = 0;  // MI backup indexes    
   }
+  
+  for (Int_t i=0;i<3;i++) { fBudget[i]=0;};
 }              
 //_____________________________________________________________________________
 AliTRDtrack::AliTRDtrack(const AliESDtrack& t) 
@@ -180,7 +231,7 @@ AliTRDtrack::AliTRDtrack(const AliESDtrack& t)
   //
   // Constructor from AliESDtrack
   //
-
+  fDE =0;     
   SetLabel(t.GetLabel());
   SetChi2(0.);
   SetMass(t.GetMass());
@@ -191,11 +242,22 @@ AliTRDtrack::AliTRDtrack(const AliESDtrack& t)
     fIndex[i] = 0; //MI store indexes
   }
   fdEdx=t.GetTRDsignal();
+  for (Int_t i=0;i<kNPlane;i++){
+    fdEdxPlane[i] = t.GetTRDsignals(i);
+    fTimBinPlane[i] = t.GetTRDTimBin(i);
+  }
 
   fLhElectron = 0.0;
   fNWrong = 0;
   fStopped = 0;
   fNRotate = 0;
+  fNExpected =0;
+  fNExpectedLast=0;
+  fNdedx = 0;
+  fNCross =0;
+  fNLast  =0;
+  fChi2Last =0;
+  fBackupTrack =0;
 
   fAlpha = t.GetAlpha();
   if      (fAlpha < -TMath::Pi()) fAlpha += 2*TMath::Pi();
@@ -204,17 +266,18 @@ AliTRDtrack::AliTRDtrack(const AliESDtrack& t)
   Double_t x, p[5]; t.GetExternalParameters(x,p);
   //Conversion of the covariance matrix
   Double_t c[15]; t.GetExternalCovariance(c);
-  if (t.GetStatus()&AliESDtrack::AliESDtrack::kTRDbackup){
-    t.GetTRDExternalParameters(x,p,c);
+  if (t.GetStatus()&AliESDtrack::kTRDbackup){
+    t.GetOuterExternalParameters(fAlpha,x,p);
+    t.GetOuterExternalCovariance(c);
+    if      (fAlpha < -TMath::Pi()) fAlpha += 2*TMath::Pi();
+    else if (fAlpha >= TMath::Pi()) fAlpha -= 2*TMath::Pi();
   }
 
   fX=x;
 
-  x = GetConvConst();  
-
   fY=p[0];
-  fZ=p[1];
-  fT=p[3];
+  fZ=p[1]; SaveLocalConvConst();
+  fT=p[3]; x=GetLocalConvConst();
   fC=p[4]/x;
   fE=fC*fX - p[2];   
 
@@ -237,40 +300,80 @@ AliTRDtrack::AliTRDtrack(const AliESDtrack& t)
     //    fIndex[i] = 0; //MI store indexes
   }
 
+  for (Int_t i=0;i<3;i++) { fBudget[i]=0;};
   if ((t.GetStatus()&AliESDtrack::kTIME) == 0) return;
   StartTimeIntegral();
   Double_t times[10]; t.GetIntegratedTimes(times); SetIntegratedTimes(times);
   SetIntegratedLength(t.GetIntegratedLength());
 
-}              
-//_____________________________________________________________________________
+}  
+
 
-void  AliTRDtrack::GetBarrelTrack(AliBarrelTrack *track) {
+AliTRDtrack * AliTRDtrack::MakeTrack(const AliTrackReference *ref, Double_t mass)
+{
   //
+  // Make dummy track from the track reference 
+  // negative mass means opposite charge 
   //
+  Double_t xx[5];
+  Double_t cc[15];
+  for (Int_t i=0;i<15;i++) cc[i]=0;
+  Double_t x = ref->X(), y = ref->Y(), z = ref->Z();
+  Double_t alpha = TMath::ATan2(y,x);
+  Double_t xr = TMath::Sqrt(x*x+y*y);
+  xx[0] = 0;
+  xx[1] = z;
+  xx[3] = ref->Pz()/ref->Pt();
+  Float_t b[3];
+  Float_t xyz[3]={x,y,z};
+  Float_t convConst = 0;
+  (AliKalmanTrack::GetFieldMap())->Field(xyz,b);
+  convConst=1000/0.299792458/(1e-13 - b[2]);
+  xx[4] = 1./(convConst*ref->Pt());
+  if (mass<0) xx[4]*=-1.;  // negative mass - negative direction
+  Double_t lcos = (x*ref->Px()+y*ref->Py())/(xr*ref->Pt());
+  Double_t lsin = TMath::Sin(TMath::ACos(lcos));
+  if (mass<0) lsin*=-1.;
+  xx[2]   = xr*xx[4]-lsin;
+  AliTRDcluster cl;
+  AliTRDtrack * track = new  AliTRDtrack(&cl,100,xx,cc,xr,alpha);
+  track->SetMass(TMath::Abs(mass));
+  track->StartTimeIntegral();  
+  return track;
+}
+
+
+AliTRDtrack::~AliTRDtrack()
+{
   //
-  
-  if (!track) return;
-  Double_t xr, vec[5], cov[15];
-
-  track->SetLabel(GetLabel());
-  track->SetX(fX, fAlpha);
-  track->SetNClusters(GetNumberOfClusters(), GetChi2());
-  track->SetNWrongClusters(fNWrong);
-  track->SetNRotate(fNRotate);
-  Double_t times[10];
-  GetIntegratedTimes(times);
-  track->SetTime(times, GetIntegratedLength());
-
-  track->SetMass(GetMass());
-  track->SetdEdX(GetdEdx());
-
-  GetExternalParameters(xr, vec);
-  track->SetStateVector(vec);
-
-  GetExternalCovariance(cov);
-  track->SetCovarianceMatrix(cov);
+  //
+
+  if (fBackupTrack) delete fBackupTrack;
+  fBackupTrack=0;
+
+}
+
+
+Float_t    AliTRDtrack::StatusForTOF()
+{
+
+  Float_t res = (0.2 + 0.8*(fN/(fNExpected+5.)))*(0.4+0.6*fTracklets[5].GetN()/20.);
+  res *= (0.25+0.8*40./(40.+fBudget[2]));
+  return res;
+
+  Int_t status=0;
+  if (GetNumberOfClusters()<20) return 0;   //
+  if (fN>110&&fChi2/(Float_t(fN))<3) return 3;            //gold
+  if (fNLast>30&&fChi2Last/(Float_t(fNLast))<3) return 3; //gold
+  if (fNLast>20&&fChi2Last/(Float_t(fNLast))<2) return 3; //gold
+  if (fNLast/(fNExpectedLast+3.)>0.8 && fChi2Last/Float_t(fNLast)<5&&fNLast>20) return 2; //silber
+  if (fNLast>5 &&((fNLast+1.)/(fNExpectedLast+1.))>0.8&&fChi2Last/(fNLast-5.)<6)   return 1; 
+  //
+
+  return status;
 }
+            
+
 //____________________________________________________________________________
 void AliTRDtrack::GetExternalParameters(Double_t& xr, Double_t x[5]) const {
   //
@@ -289,7 +392,7 @@ void AliTRDtrack::GetExternalCovariance(Double_t cc[15]) const {
   //
   // This function returns external representation of the covriance matrix.
   //
-  Double_t a=GetConvConst();
+  Double_t a=GetLocalConvConst();
 
   Double_t c22=fX*fX*fCcc-2*fX*fCce+fCee;
   Double_t c32=fX*fCct-fCte;
@@ -339,31 +442,39 @@ void AliTRDtrack::CookdEdx(Double_t low, Double_t up) {
   //-----------------------------------------------------------------
 
   Int_t i;
-  Int_t nc=GetNumberOfClusters(); 
+  //Int_t nc=GetNumberOfClusters(); 
+  Int_t nc=fNdedx; 
+  if (nc<10)  {
+    SetdEdx(0);
+    return;
+  }
 
   Float_t sorted[kMAX_CLUSTERS_PER_TRACK];
   for (i=0; i < nc; i++) {
     sorted[i]=fdQdl[i];
   }
-
-  Int_t swap; 
-
-  do {
-    swap=0;
-    for (i=0; i<nc-1; i++) {
-      if (sorted[i]<=sorted[i+1]) continue;
-      Float_t tmp=sorted[i];
-      sorted[i]=sorted[i+1]; sorted[i+1]=tmp;
-      swap++;
-    }
-  } while (swap);
-
   Int_t nl=Int_t(low*nc), nu=Int_t(up*nc);
   Float_t dedx=0;
-  for (i=nl; i<=nu; i++) dedx += sorted[i];
-  dedx /= (nu-nl+1);
+  //for (i=nl; i<=nu; i++) dedx += sorted[i];
+  //dedx /= (nu-nl+1);
+  for (i=0; i<nc; i++) dedx += sorted[i];       // ADDED by PS
+  if((nu-nl)) dedx /= (nu-nl);                  // ADDED by PS
 
+  //SetdEdx(dedx);
+  //
+  // now real truncated mean
+  for (i=0; i < nc; i++) {
+    sorted[i]=TMath::Abs(fdQdl[i]);
+  }
+  Int_t * index = new Int_t[nc];
+  TMath::Sort(nc, sorted, index,kFALSE);
+  dedx=0;
+  for (i=nl; i<=nu; i++) dedx += sorted[index[i]];
+  dedx /= (nu-nl+1);
+  fdEdxT = dedx;
+  delete [] index;
   SetdEdx(dedx);
+
 }                     
 
 
@@ -381,6 +492,7 @@ Int_t AliTRDtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho)
     //              << GetPt() << "\t" << GetLabel() << "\t" << GetMass() << endl;
     return 0;
   }
+  Double_t lcc=GetLocalConvConst();
 
   // track Length measurement [SR, GSI, 17.02.2003]
   Double_t oldX = fX, oldY = fY, oldZ = fZ;  
@@ -428,6 +540,12 @@ Int_t AliTRDtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho)
 
   fX=x2;                                                     
 
+  //Change of the magnetic field *************
+  SaveLocalConvConst();
+  cc=fC;
+  fC*=lcc/GetLocalConvConst();
+  fE+=fX*(fC-cc);
+
   //Multiple scattering  ******************
   Double_t d=sqrt((x1-fX)*(x1-fX)+(y1-fY)*(y1-fY)+(z1-fZ)*(z1-fZ));
   Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
@@ -461,16 +579,35 @@ Int_t AliTRDtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho)
   if((5940*beta2/(1-beta2+1e-10) - beta2) < 0) return 0;
 
   Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2+1e-10)) - beta2)*d*rho;
+  //
+  // suspicious part - think about it ?
+  Double_t kinE =  TMath::Sqrt(p2);
+  if (dE>0.8*kinE) dE = 0.8*kinE;  //      
+  if (dE<0)        dE = 0.0;       // not valid region for Bethe bloch 
+  //
+  //
+  fDE+=dE;
   if (x1 < x2) dE=-dE;
   cc=fC;
   fC*=(1.- sqrt(p2+GetMass()*GetMass())/p2*dE);
   fE+=fX*(fC-cc);    
+  //  Double_t sigmade = 0.1*dE*TMath::Sqrt(TMath::Sqrt(1+fT*fT)*90./(d+0.0001));   // 20 percent fluctuation - normalized to some length 
+  Double_t sigmade = 0.07*TMath::Sqrt(TMath::Abs(dE));   // energy loss fluctuation 
+  Double_t sigmac2 = sigmade*sigmade*fC*fC*(p2+GetMass()*GetMass())/(p2*p2);
+  fCcc += sigmac2;
+  fCee += fX*fX*sigmac2;  
 
   // track time measurement [SR, GSI 17.02.2002]
   if (x1 < x2)
   if (IsStartedTimeIntegral()) {
-    Double_t l2 = (fX-oldX)*(fX-oldX) + (fY-oldY)*(fY-oldY) + (fZ-oldZ)*(fZ-oldZ);
-    AddTimeStep(TMath::Sqrt(l2));
+    Double_t l2 = TMath::Sqrt((fX-oldX)*(fX-oldX) + (fY-oldY)*(fY-oldY) + (fZ-oldZ)*(fZ-oldZ));
+    if (TMath::Abs(l2*fC)>0.0001){
+      // make correction for curvature if neccesary
+      l2 = 0.5*TMath::Sqrt((fX-oldX)*(fX-oldX) + (fY-oldY)*(fY-oldY));
+      l2 = 2*TMath::ASin(l2*fC)/fC;
+      l2 = TMath::Sqrt(l2*l2+(fZ-oldZ)*(fZ-oldZ));
+    }
+    AddTimeStep(l2);
   }
 
   return 1;            
@@ -537,8 +674,9 @@ Int_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, UInt_t index,
 
     r00=c->GetSigmaY2()+errang+add, r01=0., r11=c->GetSigmaZ2()*xu_factor; 
     r00+=(fCyy+2.0*h01*fCzy+h01*h01*fCzz);
+    r01+=(fCzy+h01*fCzz);
+    r11+=fCzz;
 
-    r01+=(fCzy+h01*fCzz);  
     det=r00*r11 - r01*r01;
     tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
 
@@ -600,7 +738,7 @@ Int_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, UInt_t index,
 }                     
 //_____________________________________________________________________________
 Int_t AliTRDtrack::UpdateMI(const AliTRDcluster *c, Double_t chisq, UInt_t index, Double_t h01, 
-                           Int_t plane)
+                           Int_t /*plane*/)
 {
   // Assignes found cluster to the track and updates track information
 
@@ -615,8 +753,15 @@ Int_t AliTRDtrack::UpdateMI(const AliTRDcluster *c, Double_t chisq, UInt_t index
   }
   Double_t tangent = TMath::Sqrt(tangent2);
   if ((fC*fX-fE)<0) tangent*=-1;
-  Double_t correction = 0*plane;
-  Double_t errang = tangent2*0.04; //
+  //  Double_t correction = 0*plane;
+  Double_t errang = tangent2*0.04;  //
+  Double_t errsys =0.025*0.025*20;  //systematic error part 
+  Float_t extend =1;
+  if (c->GetNPads()==4) extend=2;
+  //if (c->GetNPads()==5)  extend=3;
+  //if (c->GetNPads()==6)  extend=3;
+  //if (c->GetQ()<15) return 1;
+
   /*
   if (corrector!=0){
   //if (0){
@@ -630,9 +775,9 @@ Int_t AliTRDtrack::UpdateMI(const AliTRDcluster *c, Double_t chisq, UInt_t index
   }
   */
   //
-  Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12.);
+  //  Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12.);
 
-  Double_t r00=c->GetSigmaY2() +errang, r01=0., r11=c->GetSigmaZ2()*10000.;
+  Double_t r00=(c->GetSigmaY2() +errang+errsys)*extend, r01=0., r11=c->GetSigmaZ2()*10000.;
   r00+=fCyy; r01+=fCzy; r11+=fCzz;
   Double_t det=r00*r11 - r01*r01;
   Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
@@ -660,19 +805,29 @@ Int_t AliTRDtrack::UpdateMI(const AliTRDcluster *c, Double_t chisq, UInt_t index
     fC  = cur;
   }
   else {
+    Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12);
+  
     Double_t xu_factor = 1000.;  // empirical factor set by C.Xu
                                 // in the first tilt version      
     dy=c->GetY() - fY; dz=c->GetZ() - fZ;     
-    dy=dy+h01*dz+correction;
+    //dy=dy+h01*dz+correction;
+    
+    Double_t tiltdz = dz;
+    if (TMath::Abs(tiltdz)>padlength/2.) {
+      tiltdz = TMath::Sign(padlength/2,dz);
+    }
+    //    dy=dy+h01*dz;
+    dy=dy+h01*tiltdz;
+
     Double_t add=0;
     if (TMath::Abs(dz)>padlength/2.){
       //Double_t dy2 = c->GetY() - fY;
       //Double_t sign = (dz>0) ? -1.: 1.;
       //dy2-=h01*sign*padlength/2.;    
       //dy = dy2;
-      add =1.;
+      add =1;
     }
-    Double_t s00 = c->GetSigmaY2()+errang+add;  // error pad
+    Double_t s00 = (c->GetSigmaY2()+errang)*extend+errsys+add;  // error pad
     Double_t s11 = c->GetSigmaZ2()*xu_factor;   // error pad-row
     //
     r00 = fCyy + 2*fCzy*h01 + fCzz*h01*h01+s00;
@@ -750,12 +905,127 @@ Int_t AliTRDtrack::UpdateMI(const AliTRDcluster *c, Double_t chisq, UInt_t index
 }                     
 
 
+
+//_____________________________________________________________________________
+Int_t AliTRDtrack::UpdateMI(const AliTRDtracklet &tracklet)
+{
+  //
+  // Assignes found tracklet to the track and updates track information
+  //
+  //
+  Double_t r00=(tracklet.GetTrackletSigma2()), r01=0., r11= 10000.;
+  r00+=fCyy; r01+=fCzy; r11+=fCzz;
+  //
+  Double_t det=r00*r11 - r01*r01;
+  Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
+  //
+
+  Double_t dy=tracklet.GetY() - fY, dz=tracklet.GetZ() - fZ;
+
+  
+  Double_t s00 = tracklet.GetTrackletSigma2();  // error pad
+  Double_t s11 = 100000;   // error pad-row
+  Float_t  h01 = tracklet.GetTilt();
+  //
+  //  r00 = fCyy + 2*fCzy*h01 + fCzz*h01*h01+s00;
+  r00 = fCyy + fCzz*h01*h01+s00;
+  //  r01 = fCzy + fCzz*h01;
+  r01 = fCzy ;
+  r11 = fCzz + s11;
+  det = r00*r11 - r01*r01;
+  // inverse matrix
+  tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
+
+  Double_t k00=fCyy*r00+fCzy*r01, k01=fCyy*r01+fCzy*r11;
+  Double_t k10=fCzy*r00+fCzz*r01, k11=fCzy*r01+fCzz*r11;
+  Double_t k20=fCey*r00+fCez*r01, k21=fCey*r01+fCez*r11;
+  Double_t k30=fCty*r00+fCtz*r01, k31=fCty*r01+fCtz*r11;
+  Double_t k40=fCcy*r00+fCcz*r01, k41=fCcy*r01+fCcz*r11;
+  
+  // K matrix
+//   k00=fCyy*r00+fCzy*(r01+h01*r00),k01=fCyy*r01+fCzy*(r11+h01*r01);
+//   k10=fCzy*r00+fCzz*(r01+h01*r00),k11=fCzy*r01+fCzz*(r11+h01*r01);
+//   k20=fCey*r00+fCez*(r01+h01*r00),k21=fCey*r01+fCez*(r11+h01*r01);
+//   k30=fCty*r00+fCtz*(r01+h01*r00),k31=fCty*r01+fCtz*(r11+h01*r01);
+//   k40=fCcy*r00+fCcz*(r01+h01*r00),k41=fCcy*r01+fCcz*(r11+h01*r01);  
+  //
+  //Update measurement
+  Double_t cur=fC + k40*dy + k41*dz, eta=fE + k20*dy + k21*dz;  
+  //  cur=fC + k40*dy + k41*dz; eta=fE + k20*dy + k21*dz;
+  if (TMath::Abs(cur*fX-eta) >= 0.90000) {
+    //Int_t n=GetNumberOfClusters();
+    //      if (n>4) cerr<<n<<" AliTRDtrack warning: Filtering failed !\n";
+    return 0;
+  }                           
+//   k01+=h01*k00;
+//   k11+=h01*k10;
+//   k21+=h01*k20;
+//   k31+=h01*k30;
+//   k41+=h01*k40;  
+
+
+  fY += k00*dy + k01*dz;
+  fZ += k10*dy + k11*dz;
+  fE  = eta;
+  fT += k30*dy + k31*dz;
+  fC  = cur;
+    
+  
+  //Update covariance
+  //
+  //
+  Double_t oldyy = fCyy, oldzz = fCzz; //, oldee=fCee, oldcc =fCcc;
+  Double_t oldzy = fCzy, oldey = fCey, oldty=fCty, oldcy =fCcy;
+  Double_t oldez = fCez, oldtz = fCtz, oldcz=fCcz;
+  //Double_t oldte = fCte, oldce = fCce;
+  //Double_t oldct = fCct;
+
+  fCyy-=k00*oldyy+k01*oldzy;   
+  fCzy-=k10*oldyy+k11*oldzy;
+  fCey-=k20*oldyy+k21*oldzy;   
+  fCty-=k30*oldyy+k31*oldzy;
+  fCcy-=k40*oldyy+k41*oldzy;  
+  //
+  fCzz-=k10*oldzy+k11*oldzz;
+  fCez-=k20*oldzy+k21*oldzz;   
+  fCtz-=k30*oldzy+k31*oldzz;
+  fCcz-=k40*oldzy+k41*oldzz;
+  //
+  fCee-=k20*oldey+k21*oldez;   
+  fCte-=k30*oldey+k31*oldez;
+  fCce-=k40*oldey+k41*oldez;
+  //
+  fCtt-=k30*oldty+k31*oldtz;
+  fCct-=k40*oldty+k41*oldtz;
+  //
+  fCcc-=k40*oldcy+k41*oldcz;                 
+  //
+  /*
+  Int_t n=GetNumberOfClusters();
+  fIndex[n]=index;
+  SetNumberOfClusters(n+1);
+
+  SetChi2(GetChi2()+chisq);
+  //  cerr<<"in update: fIndex["<<fN<<"] = "<<index<<endl;
+  */
+  return 1;      
+}                     
+
+
+
 //_____________________________________________________________________________
-Int_t AliTRDtrack::Rotate(Double_t alpha)
+Int_t AliTRDtrack::Rotate(Double_t alpha, Bool_t absolute)
 {
   // Rotates track parameters in R*phi plane
+  // if absolute rotation alpha is in global system
+  // otherwise alpha rotation is relative to the current rotation angle
   
-  fNRotate++;
+  if (absolute) {
+    alpha -= fAlpha;
+  }
+  else{
+    fNRotate++;
+  }
 
   fAlpha += alpha;
   if (fAlpha<-TMath::Pi()) fAlpha += 2*TMath::Pi();
@@ -828,6 +1098,8 @@ Double_t AliTRDtrack::GetPredictedChi2(const AliTRDcluster *c, Double_t h01) con
     chi2 = (dy*dy)/r00;    
   }
   else {
+    Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12);
+    //
     r00=c->GetSigmaY2(); r01=0.; r11=c->GetSigmaZ2();
     r00+=fCyy; r01+=fCzy; r11+=fCzz;
 
@@ -839,7 +1111,12 @@ Double_t AliTRDtrack::GetPredictedChi2(const AliTRDcluster *c, Double_t h01) con
     }
     Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
     Double_t dy=c->GetY() - fY, dz=c->GetZ() - fZ;
-    dy=dy+h01*dz;
+    Double_t tiltdz = dz;
+    if (TMath::Abs(tiltdz)>padlength/2.) {
+      tiltdz = TMath::Sign(padlength/2,dz);
+    }
+    //    dy=dy+h01*dz;
+    dy=dy+h01*tiltdz;
 
     chi2 = (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det; 
   }
@@ -901,8 +1178,121 @@ void AliTRDtrack::ResetCovariance(Float_t mult) {
   //
 
   fCyy*=mult;
-  fCzy*=0.;  fCzz*=mult;
+  fCzy*=0.;  fCzz*=1.;
   fCey*=0.;  fCez*=0.;  fCee*=mult;
-  fCty*=0.;  fCtz*=0.;  fCte*=0.;  fCtt*=mult;
+  fCty*=0.;  fCtz*=0.;  fCte*=0.;  fCtt*=1.;
   fCcy*=0.;  fCcz*=0.;  fCce*=0.;  fCct*=0.;  fCcc*=mult;  
 }                                                         
+
+
+
+
+
+void AliTRDtrack::MakeBackupTrack()
+{
+  //
+  //
+  if (fBackupTrack) delete fBackupTrack;
+  fBackupTrack = new AliTRDtrack(*this);
+  
+}
+
+Int_t  AliTRDtrack::GetProlongation(Double_t xk, Double_t &y, Double_t &z){
+  //
+  // Find prolongation at given x
+  // return 0 if not exist
+  
+  Double_t c1=fC*fX - fE;
+  if (TMath::Abs(c1)>1.) return 0;
+  Double_t r1=TMath::Sqrt(1.- c1*c1);
+  Double_t c2=fC*xk - fE;
+  if (TMath::Abs(c2)>1.) return 0;  
+  Double_t r2=TMath::Sqrt(1.- c2*c2);
+  y =fY + (xk-fX)*(c1+c2)/(r1+r2);
+  z =fZ + (xk-fX)*(c1+c2)/(c1*r2 + c2*r1)*fT;
+
+  return 1;
+  
+}
+
+
+Int_t   AliTRDtrack::PropagateToX(Double_t xr, Double_t step)
+{
+  //
+  // Propagate track to given x  position 
+  // works inside of the 20 degree segmentation (local cooordinate frame for TRD , TPC, TOF)
+  // 
+  // material budget from geo manager
+  // 
+  Double_t  xyz0[3], xyz1[3],y,z;
+  const Double_t alphac = TMath::Pi()/9.;   
+  const Double_t talphac = TMath::Tan(alphac*0.5);
+  // critical alpha  - cross sector indication
+  //
+  Double_t dir = (fX>xr) ? -1.:1.;
+  // direction +-
+  for (Double_t x=fX+dir*step;dir*x<dir*xr;x+=dir*step){
+    //
+    GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]);     
+    GetProlongation(x,y,z);
+    xyz1[0] = x*TMath::Cos(fAlpha)+y*TMath::Sin(fAlpha); 
+    xyz1[1] = x*TMath::Sin(fAlpha)-y*TMath::Cos(fAlpha);
+    xyz1[2] = z;
+    Double_t param[7];
+    AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
+    //
+    if (param[0]>0&&param[1]>0) PropagateTo(x,param[1],param[0]);
+    if (fY>fX*talphac){
+      Rotate(-alphac);
+    }
+    if (fY<-fX*talphac){
+      Rotate(alphac);
+    }
+  }
+  //
+  PropagateTo(xr);
+  return 0;
+}
+
+
+Int_t   AliTRDtrack::PropagateToR(Double_t r,Double_t step)
+{
+  //
+  // propagate track to the radial position
+  // rotation always connected to the last track position
+  //
+  Double_t  xyz0[3], xyz1[3],y,z; 
+  Double_t radius = TMath::Sqrt(fX*fX+fY*fY);
+  Double_t dir = (radius>r) ? -1.:1.;   // direction +-
+  //
+  for (Double_t x=radius+dir*step;dir*x<dir*r;x+=dir*step){
+    GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]);     
+    Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
+    Rotate(alpha,kTRUE);
+    GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]);     
+    GetProlongation(x,y,z);
+    xyz1[0] = x*TMath::Cos(alpha)+y*TMath::Sin(alpha); 
+    xyz1[1] = x*TMath::Sin(alpha)-y*TMath::Cos(alpha);
+    xyz1[2] = z;
+    Double_t param[7];
+    AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
+    if (param[1]<=0) param[1] =100000000;
+    PropagateTo(x,param[1],param[0]);
+  } 
+  GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]);       
+  Double_t alpha = TMath::ATan2(xyz0[1],xyz0[0]);
+  Rotate(alpha,kTRUE);
+  GetGlobalXYZ(xyz0[0],xyz0[1],xyz0[2]);       
+  GetProlongation(r,y,z);
+  xyz1[0] = r*TMath::Cos(alpha)+y*TMath::Sin(alpha); 
+  xyz1[1] = r*TMath::Sin(alpha)-y*TMath::Cos(alpha);
+  xyz1[2] = z;
+  Double_t param[7];
+  AliKalmanTrack::MeanMaterialBudget(xyz0,xyz1,param);
+  //
+  if (param[1]<=0) param[1] =100000000;
+  PropagateTo(r,param[1],param[0]);
+  return 0;
+}
+
+