Updated version of the TRD reconstruction (M.Ivanov)
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 21 Jun 2005 15:26:01 +0000 (15:26 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Tue, 21 Jun 2005 15:26:01 +0000 (15:26 +0000)
TRD/AliTRDclusterizerV1.cxx
TRD/AliTRDclusterizerV1.h
TRD/AliTRDtrack.cxx
TRD/AliTRDtrack.h
TRD/AliTRDtracker.cxx
TRD/AliTRDtracker.h
TRD/TRDrecLinkDef.h

index 24991b2..c41b255 100644 (file)
@@ -164,10 +164,10 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
   }
   fPar->Init();
 
-  Float_t timeBinSize = fPar->GetDriftVelocity()
-                      / fPar->GetSamplingFrequency();
+  //Float_t timeBinSize = fPar->GetDriftVelocity()
+  //                    / fPar->GetSamplingFrequency();
   // Half of ampl.region
-  const Float_t kAmWidth = AliTRDgeometry::AmThick()/2.; 
+  //  const Float_t kAmWidth = AliTRDgeometry::AmThick()/2.; 
 
   Float_t omegaTau = fPar->GetOmegaTau();
   if (fVerbose > 0) {
@@ -186,7 +186,6 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
   Int_t maxThresh = fPar->GetClusMaxThresh();   
   // Threshold value for the digit signal
   Int_t sigThresh = fPar->GetClusSigThresh();   
-
   // Iteration limit for unfolding procedure
   const Float_t kEpsilon = 0.01;             
 
@@ -195,11 +194,11 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
   const Int_t   kNtrack  = 3 * kNclus;
 
   Int_t    iType         = 0;
-  Int_t    iUnfold       = 0;
-
+  Int_t    iUnfold       = 0;  
   Double_t ratioLeft     = 1.0;
   Double_t ratioRight    = 1.0;
 
+  //
   Double_t padSignal[kNsig];   
   Double_t clusterSignal[kNclus];
   Double_t clusterPads[kNclus];   
@@ -260,12 +259,19 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
               Int_t signalM = TMath::Abs(digits->GetDataUnchecked(row,col-1,time));
               Int_t signalR = TMath::Abs(digits->GetDataUnchecked(row,col-2,time));
  
+//           // Look for the maximum
+//               if (signalM >= maxThresh) {
+//                 if (((signalL >= sigThresh) &&
+//                      (signalL <  signalM))  ||
+//                     ((signalR >= sigThresh) &&
+//                      (signalR <  signalM))) {
+//                   // Maximum found, mark the position by a negative signal
+//                   digits->SetDataUnchecked(row,col-1,time,-signalM);
+//             }
+//           }
              // Look for the maximum
               if (signalM >= maxThresh) {
-                if (((signalL >= sigThresh) &&
-                     (signalL <  signalM))  ||
-                    ((signalR >= sigThresh) &&
-                     (signalR <  signalM))) {
+                if ( (signalL<=signalM) && (signalR<=signalM) && (signalL+signalR)>sigThresh ) {
                   // Maximum found, mark the position by a negative signal
                   digits->SetDataUnchecked(row,col-1,time,-signalM);
                }
@@ -336,10 +342,7 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
                   break;
                };
 
-               // Don't analyze large clusters
-                //if (iType == 4) continue;
-
-                // Look for 5 pad cluster with minimum in the middle
+                // Look for 5 pad cluster with minimum in the middle
                 Bool_t fivePadCluster = kFALSE;
                 if (col < nColMax-3) {
                   if (digits->GetDataUnchecked(row,col+2,time) < 0) {
@@ -362,7 +365,7 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
                // of the cluster which remains from a previous unfolding
                 if (iUnfold) {
                   clusterSignal[0] *= ratioLeft;
-                  iType   = 3;
+                  iType   = 5;
                   iUnfold = 0;
                }
 
@@ -378,10 +381,11 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
                   ratioRight        = Unfold(kEpsilon,iplan,padSignal);
                   ratioLeft         = 1.0 - ratioRight; 
                   clusterSignal[2] *= ratioRight;
-                  iType   = 3;
+                  iType   = 5;
                   iUnfold = 1;
                 }
 
+
                 Double_t clusterCharge = clusterSignal[0]
                                        + clusterSignal[1]
                                        + clusterSignal[2];
@@ -392,29 +396,31 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
                 clusterPads[2] = time - nTimeBefore + 0.5;
 
                 if (fPar->LUTOn()) {
-
                  // Calculate the position of the cluster by using the
                  // lookup table method
-//                   clusterPads[1] = col + 0.5
-//                                  + fPar->LUTposition(iplan,clusterSignal[0]
-//                                                           ,clusterSignal[1]
-//                                                       ,clusterSignal[2]);
-                  clusterPads[1] = 0.5
-                                 + fPar->LUTposition(iplan,clusterSignal[0]
+                  clusterPads[1] =
+                                  fPar->LUTposition(iplan,clusterSignal[0]
                                                           ,clusterSignal[1]
                                                          ,clusterSignal[2]);
-
                }
                else {
-
                  // Calculate the position of the cluster by using the
                  // center of gravity method
-//                   clusterPads[1] = col + 0.5 
-//                                  + (clusterSignal[2] - clusterSignal[0]) 
-//                              / clusterCharge;
-                  clusterPads[1] = 0.5 
-                                 + (clusterSignal[2] - clusterSignal[0]) 
-                                / clusterCharge;
+                 for (Int_t i=0;i<5;i++) padSignal[i]=0;
+                 padSignal[2] = TMath::Abs(digits->GetDataUnchecked(row,col,time));   // central  pad
+                 padSignal[1] = TMath::Abs(digits->GetDataUnchecked(row,col-1,time)); // left     pad
+                 padSignal[3] = TMath::Abs(digits->GetDataUnchecked(row,col+1,time)); // right    pad
+                 if (col>2 &&TMath::Abs(digits->GetDataUnchecked(row,col-2,time)<padSignal[1])){
+                   padSignal[0] = TMath::Abs(digits->GetDataUnchecked(row,col-2,time));
+                 }
+                 if (col<nColMax-3 &&TMath::Abs(digits->GetDataUnchecked(row,col+2,time)<padSignal[3])){
+                   padSignal[4] = TMath::Abs(digits->GetDataUnchecked(row,col+2,time));
+                 }               
+                 clusterPads[1] =  GetCOG(padSignal);
+                 Double_t check = fPar->LUTposition(iplan,clusterSignal[0]
+                                                          ,clusterSignal[1]
+                                                         ,clusterSignal[2]);
+                 //              clusterPads[1] = check;
 
                }
 
@@ -424,49 +430,16 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
                 Double_t clusterSigmaY2 = (q1*(q0+q2)+4*q0*q2) /
                                           (clusterCharge*clusterCharge);
 
-                if (fVerbose > 1) {
-                  printf("-----------------------------------------------------------\n");
-                  printf("Create cluster no. %d\n",nClusters);
-                  printf("Position: row = %f, col = %f, time = %f\n",clusterPads[0]
-                                                                   ,clusterPads[1]
-                                                                    ,clusterPads[2]);
-                  printf("Indices: %d, %d, %d\n",clusterDigit[0]
-                                               ,clusterDigit[1]
-                                                ,clusterDigit[2]);
-                  printf("Total charge = %f\n",clusterCharge);
-                  printf("Tracks: pad0 %d, %d, %d\n",clusterTracks[0]
-                                                   ,clusterTracks[1]
-                                                    ,clusterTracks[2]);
-                  printf("        pad1 %d, %d, %d\n",clusterTracks[3]
-                                                   ,clusterTracks[4]
-                                                    ,clusterTracks[5]);
-                  printf("        pad2 %d, %d, %d\n",clusterTracks[6]
-                                                   ,clusterTracks[7]
-                                                    ,clusterTracks[8]);
-                  printf("Type = %d, Number of pads = %d\n",iType,nPadCount);
-                }
-
                // Calculate the position and the error
+                Double_t colSize = padPlane->GetColSize(col);
+                Double_t rowSize = padPlane->GetRowSize(row);
                 Double_t clusterPos[3];
-//                 clusterPos[0] = clusterPads[1] * colSize + col0;
-//                 clusterPos[1] = clusterPads[0] * rowSize + row0;
-                clusterPos[0] = padPlane->GetColPos(col) - clusterPads[1];
-                clusterPos[1] = padPlane->GetRowPos(row) - clusterPads[0];
+               clusterPos[0] = padPlane->GetColPos(col) + (clusterPads[1]-0.5)*colSize;  // MI change
+               clusterPos[1] = padPlane->GetRowPos(row) -0.5*rowSize; //MI change
                 clusterPos[2] = clusterPads[2];
                 Double_t clusterSig[2];
-                Double_t colSize = padPlane->GetColSize(col);
-                Double_t rowSize = padPlane->GetRowSize(row);
                 clusterSig[0] = (clusterSigmaY2 + 1./12.) * colSize*colSize;
-                clusterSig[1] = rowSize * rowSize / 12.;
-
-                // Correct for ExB displacement
-                if (fPar->ExBOn()) { 
-                  Int_t    local_time_bin = (Int_t) clusterPads[2];
-                  Double_t driftLength    = local_time_bin * timeBinSize + kAmWidth;
-                  Double_t deltaY         = omegaTau * driftLength;
-                  clusterPos[1]           = clusterPos[1] - deltaY;
-                }
-                                       
+                clusterSig[1] = rowSize * rowSize / 12.;                                       
                 // Add the cluster to the output array 
                 AddCluster(clusterPos
                          ,idet
@@ -483,24 +456,12 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
        // Compress the arrays
         digits->Compress(1,0);
         track0->Compress(1,0);
- track1->Compress(1,0);
+       track1->Compress(1,0);
         track2->Compress(1,0);
 
         // Write the cluster and reset the array
        WriteClusters(idet);
        ResetRecPoints();
-
-        if (fVerbose > 0) {
-          printf("<AliTRDclusterizerV1::MakeCluster> ");
-          printf("Found %d clusters in total.\n"
-                ,nClusters);
-          printf("                                    2pad:  %d\n",nClusters2pad);
-          printf("                                    3pad:  %d\n",nClusters3pad);
-          printf("                                    4pad:  %d\n",nClusters4pad);
-          printf("                                    5pad:  %d\n",nClusters5pad);
-          printf("                                    Large: %d\n",nClustersLarge);
-       }
-
       }    
     }      
   }        
@@ -514,6 +475,18 @@ Bool_t AliTRDclusterizerV1::MakeClusters()
 
 }
 
+Double_t AliTRDclusterizerV1::GetCOG(Double_t signal[5])
+{
+  //
+  // get COG position
+  // used for clusters with more than 3 pads - where LUT not applicable
+  Double_t sum = signal[0]+signal[1]+signal[2]+signal[3]+signal[4];
+  Double_t res = (0.0*(-signal[0]+signal[4])+(-signal[1]+signal[3]))/sum;
+  return res;            
+}
+
+
+
 //_____________________________________________________________________________
 Double_t AliTRDclusterizerV1::Unfold(Double_t eps, Int_t plane, Double_t* padSignal)
 {
index 5d5e461..52c7a14 100644 (file)
@@ -33,9 +33,8 @@ class AliTRDclusterizerV1 : public AliTRDclusterizer {
   AliTRDdigitsManager *fDigitsManager;      //! TRD digits manager
 
  private:
-
   virtual Double_t Unfold(Double_t eps, Int_t plane, Double_t *padSignal);
-
+  Double_t GetCOG(Double_t signal[5]);      // get COG position
   ClassDef(AliTRDclusterizerV1,5)           // TRD-Cluster finder, slow simulator
 
 };
index 214f814..aaba24d 100644 (file)
 #include "AliTRDtrack.h"
 #include "AliTRDclusterCorrection.h"
 
+ClassImp(AliTRDtracklet)
 ClassImp(AliTRDtrack)
 
+
+  AliTRDtracklet::AliTRDtracklet():fY(0),fX(0),fAlpha(0),fSigma2(0),fP0(0),fP1(0),fNFound(0),fNCross(0),fPlane(0),fExpectedSigma2(0),fChi2(0){
+}
+
 //_____________________________________________________________________________
 
 AliTRDtrack::AliTRDtrack(const AliTRDcluster *c, UInt_t index, 
                          const Double_t xx[5], const Double_t cc[15], 
-                         Double_t xref, Double_t alpha) : AliKalmanTrack() 
-{
+                         Double_t xref, Double_t alpha) : AliKalmanTrack() {
   //-----------------------------------------------------------------
   // This is the main track constructor.
   //-----------------------------------------------------------------
@@ -77,7 +81,7 @@ AliTRDtrack::AliTRDtrack(const AliTRDcluster *c, UInt_t index,
   fdQdl[0] = q;
   
   // initialisation [SR, GSI 18.02.2003] (i startd for 1)
-  for(UInt_t i=1; i<kMAXCLUSTERSPERTRACK; i++) {
+  for(UInt_t i=1; i<kMAX_CLUSTERS_PER_TRACK; i++) {
     fdQdl[i] = 0;
     fIndex[i] = 0;
     fIndexBackup[i] = 0;  //bacup indexes MI    
@@ -133,11 +137,14 @@ AliTRDtrack::AliTRDtrack(const AliTRDtrack& t) : AliKalmanTrack(t) {
   }
 
   // initialisation (i starts from n) [SR, GSI, 18.02.2003]
-  for(UInt_t i=n; i<kMAXCLUSTERSPERTRACK; i++) {
+  for(UInt_t i=n; i<kMAX_CLUSTERS_PER_TRACK; i++) {
     fdQdl[i] = 0;
     fIndex[i] = 0;
     fIndexBackup[i] = 0;  //MI backup indexes
   }
+  for (Int_t i=0;i<6;i++){
+    fTracklets[i]=t.fTracklets[i];
+  }
 }                                
 
 //_____________________________________________________________________________
@@ -202,14 +209,13 @@ AliTRDtrack::AliTRDtrack(const AliKalmanTrack& t, Double_t alpha)
   fCcy=c[10];   fCcz=c[11];   fCce=c42;   fCct=c[13]; fCcc=c[14];  
 
   // Initialization [SR, GSI, 18.02.2003]
-  for(UInt_t i=0; i<kMAXCLUSTERSPERTRACK; i++) {
+  for(UInt_t i=0; i<kMAX_CLUSTERS_PER_TRACK; i++) {
     fdQdl[i] = 0;
     fIndex[i] = 0;
     fIndexBackup[i] = 0;  // MI backup indexes    
   }
   
 }              
-
 //_____________________________________________________________________________
 AliTRDtrack::AliTRDtrack(const AliESDtrack& t) 
            :AliKalmanTrack() {
@@ -222,7 +228,7 @@ AliTRDtrack::AliTRDtrack(const AliESDtrack& t)
   SetMass(t.GetMass());
   SetNumberOfClusters(t.GetTRDclusters(fIndex)); 
   Int_t ncl = t.GetTRDclusters(fIndexBackup);
-  for (UInt_t i=ncl;i<kMAXCLUSTERSPERTRACK;i++) {
+  for (UInt_t i=ncl;i<kMAX_CLUSTERS_PER_TRACK;i++) {
     fIndexBackup[i]=0;
     fIndex[i] = 0; //MI store indexes
   }
@@ -281,7 +287,7 @@ AliTRDtrack::AliTRDtrack(const AliESDtrack& t)
   fCcy=c[10];   fCcz=c[11];   fCce=c42;   fCct=c[13]; fCcc=c[14];  
 
   // Initialization [SR, GSI, 18.02.2003]
-  for(UInt_t i=0; i<kMAXCLUSTERSPERTRACK; i++) {
+  for(UInt_t i=0; i<kMAX_CLUSTERS_PER_TRACK; i++) {
     fdQdl[i] = 0;
     //    fIndex[i] = 0; //MI store indexes
   }
@@ -293,11 +299,9 @@ AliTRDtrack::AliTRDtrack(const AliESDtrack& t)
 
 }  
 
-//_____________________________________________________________________________
 AliTRDtrack::~AliTRDtrack()
 {
   //
-  // TRD track destructor
   //
 
   if (fBackupTrack) delete fBackupTrack;
@@ -305,25 +309,24 @@ AliTRDtrack::~AliTRDtrack()
 
 }
 
-//_____________________________________________________________________________
-Float_t AliTRDtrack::StatusForTOF()
-{
-  //
-  // Returns a status flag for TOF
-  //
 
+Float_t    AliTRDtrack::StatusForTOF()
+{
   Int_t status=0;
+  if (fNRotate>2) return -1;   // sure it's stopped
   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; 
-  //
+
+  //comp->fTree->SetAlias("nlast2","track.fTracklets[5].fNFound+track.fTracklets[4].fNFound");
+  //comp->fTree->SetAlias("goldtrack","abs((track.fTracklets[5].fP1+track.fTracklets[4].fP1))<0.5&&nlast2>14");
+  Int_t nlast2 = fTracklets[5].fNFound+fTracklets[4].fNFound;
+  if (abs((fTracklets[5].fP1+fTracklets[4].fP1))<0.3 &&nlast2>20) return 3;
+  if (abs((fTracklets[5].fP1+fTracklets[4].fP1))<0.3 &&nlast2>14) return 2;
+  if (abs((fTracklets[5].fP1+fTracklets[4].fP1))<0.5 &&nlast2>14) return 1;
 
   return status;
 }
             
+
 //____________________________________________________________________________
 void AliTRDtrack::GetExternalParameters(Double_t& xr, Double_t x[5]) const {
   //
@@ -356,12 +359,9 @@ void AliTRDtrack::GetExternalCovariance(Double_t cc[15]) const {
   
 }               
                        
+
 //_____________________________________________________________________________
-void AliTRDtrack::GetCovariance(Double_t cc[15]) const 
-{
-  //
-  // Returns the covariance matrix
-  //
+void AliTRDtrack::GetCovariance(Double_t cc[15]) const {
 
   cc[0]=fCyy;
   cc[1]=fCzy;  cc[2]=fCzz;
@@ -402,33 +402,22 @@ void AliTRDtrack::CookdEdx(Double_t low, Double_t up) {
     return;
   }
 
-  Float_t sorted[kMAXCLUSTERSPERTRACK];
+  Float_t sorted[kMAX_CLUSTERS_PER_TRACK];
   for (i=0; i < nc; i++) {
-    sorted[i]=fdQdl[i];
+    sorted[i]=TMath::Abs(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 * index = new Int_t[nc];
+  TMath::Sort(nc, sorted, index,kFALSE);
+
   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=0; i<nc; i++) dedx += sorted[i];       // ADDED by PS
-  if((nu-nl)) dedx /= (nu-nl);                  // ADDED by PS
+  for (i=nl; i<=nu; i++) dedx += sorted[index[i]];
+  dedx /= (nu-nl+1);
 
   SetdEdx(dedx);
 }                     
 
+
 //_____________________________________________________________________________
 Int_t AliTRDtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho)
 {
@@ -495,7 +484,9 @@ Int_t AliTRDtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho)
   Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
   Double_t beta2=p2/(p2 + GetMass()*GetMass());
   Double_t theta2=14.1*14.1/(beta2*p2*1e6)*d/x0*rho;
-
+  theta2*= 3.;                                      // magic const - to normalize pools - waiting for geo manager
+  if (p2>2.) beta2*= 0.4;                           // magic const - theta2 for relativistic particles 
+                                                    // - not valid for electrons
   Double_t ey=fC*fX - fE, ez=fT;
   Double_t xz=fC*ez, zz1=ez*ez+1, xy=fE+ey;
   
@@ -522,12 +513,18 @@ Int_t AliTRDtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho)
   //Energy losses************************
   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;
+  Double_t dE=0.153e-3/beta2*(log(5940*beta2/(1-beta2+1e-10)) - beta2)*d*rho; 
+  dE *= 1.2;                                      // magic const - to normalize pools - waiting for geo manager
   if (x1 < x2) dE=-dE;
   cc=fC;
   fC*=(1.- sqrt(p2+GetMass()*GetMass())/p2*dE);
   fE+=fX*(fC-cc);    
-
+  //
+  Double_t deltac = fC-cc;
+  fCcc   += 4*deltac*deltac;                    // fluctuation of energy losses
+  fCee   += 4*fX*fX*deltac*deltac;              // local angle unchanged
+  fCce   += 4*fX*deltac*deltac;                 // correlation 1
+  //
   // track time measurement [SR, GSI 17.02.2002]
   if (x1 < x2)
   if (IsStartedTimeIntegral()) {
@@ -538,6 +535,7 @@ Int_t AliTRDtrack::PropagateTo(Double_t xk,Double_t x0,Double_t rho)
   return 1;            
 }     
 
+
 //_____________________________________________________________________________
 Int_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, UInt_t index,                          Double_t h01)
 {
@@ -581,7 +579,7 @@ Int_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, UInt_t index,
     fC  = cur;
   }
   else {
-    Double_t xuFactor = 100.;  // empirical factor set by C.Xu
+    Double_t xu_factor = 100.;  // empirical factor set by C.Xu
                                 // in the first tilt version      
     dy=c->GetY() - fY; dz=c->GetZ() - fZ;     
     dy=dy+h01*dz;
@@ -596,7 +594,7 @@ Int_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, UInt_t index,
    
 
 
-    r00=c->GetSigmaY2()+errang+add, r01=0., r11=c->GetSigmaZ2()*xuFactor; 
+    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);  
@@ -659,7 +657,6 @@ Int_t AliTRDtrack::Update(const AliTRDcluster *c, Double_t chisq, UInt_t index,
 
   return 1;     
 }                     
-
 //_____________________________________________________________________________
 Int_t AliTRDtrack::UpdateMI(const AliTRDcluster *c, Double_t chisq, UInt_t index, Double_t h01, 
                            Int_t /*plane*/)
@@ -677,14 +674,7 @@ 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 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;
+  Double_t errsys =(0.025*0.025*20)*(1+tangent2);  //systematic error part 
 
   /*
   if (corrector!=0){
@@ -699,9 +689,7 @@ Int_t AliTRDtrack::UpdateMI(const AliTRDcluster *c, Double_t chisq, UInt_t index
   }
   */
   //
-  //  Double_t padlength = TMath::Sqrt(c->GetSigmaZ2()*12.);
-
-  Double_t r00=(c->GetSigmaY2() +errang+errsys)*extend, r01=0., r11=c->GetSigmaZ2()*10000.;
+  Double_t r00=(c->GetSigmaY2()+errsys), 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;
@@ -729,30 +717,17 @@ 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 padlength = TMath::Sqrt(c->GetSigmaZ2()*12);
   
-    Double_t xuFactor = 1000.;  // empirical factor set by C.Xu
+    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;
     
     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;
-    }
-    Double_t s00 = (c->GetSigmaY2()+errang)*extend+errsys+add;  // error pad
-    Double_t s11 = c->GetSigmaZ2()*xuFactor;   // error pad-row
+    Double_t s00 = c->GetSigmaY2()+errsys;  // error pad
+    Double_t s11 = c->GetSigmaZ2()*xu_factor;   // error pad-row
     //
     r00 = fCyy + 2*fCzy*h01 + fCzz*h01*h01+s00;
     r01 = fCzy + fCzz*h01;
@@ -828,6 +803,116 @@ Int_t AliTRDtrack::UpdateMI(const AliTRDcluster *c, Double_t chisq, UInt_t index
   return 1;      
 }                     
 
+
+
+
+
+//_____________________________________________________________________________
+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)
 {
@@ -891,13 +976,11 @@ Int_t AliTRDtrack::Rotate(Double_t alpha)
   return 1;                            
 }                         
 
+
 //_____________________________________________________________________________
 Double_t AliTRDtrack::GetPredictedChi2(const AliTRDcluster *c, Double_t h01) const
 {
-  //
-  // Returns the predicted chi^2
-  //  
-
+  
   Bool_t fNoTilt = kTRUE;
   if(TMath::Abs(h01) > 0.003) fNoTilt = kFALSE;
   Double_t chi2, dy, r00, r01, r11;
@@ -931,9 +1014,9 @@ Double_t AliTRDtrack::GetPredictedChi2(const AliTRDcluster *c, Double_t h01) con
     chi2 = (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det; 
   }
   return chi2;
-
 }      
 
+
 //_________________________________________________________________________
 void AliTRDtrack::GetPxPyPz(Double_t& px, Double_t& py, Double_t& pz) const
 {
@@ -982,7 +1065,6 @@ void AliTRDtrack::ResetCovariance() {
   fCcy=0.;  fCcz=0.;  fCce=0.;  fCct=0.;  fCcc*=10.;  
 }                                                         
 
-//_________________________________________________________________________
 void AliTRDtrack::ResetCovariance(Float_t mult) {
   //
   // Resets covariance matrix
@@ -995,31 +1077,29 @@ void AliTRDtrack::ResetCovariance(Float_t mult) {
   fCcy*=0.;  fCcz*=0.;  fCce*=0.;  fCct*=0.;  fCcc*=mult;  
 }                                                         
 
-//_________________________________________________________________________
+
+
+
+
 void AliTRDtrack::MakeBackupTrack()
 {
   //
-  // Creates a backup track
   //
-
   if (fBackupTrack) delete fBackupTrack;
   fBackupTrack = new AliTRDtrack(*this);
   
 }
 
-//_________________________________________________________________________
-Int_t  AliTRDtrack::GetProlongation(Double_t xk, Double_t &y, Double_t &z)
-{
+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;
+  if (TMath::Abs(c1)>0.9) return 0;
   Double_t r1=TMath::Sqrt(1.- c1*c1);
   Double_t c2=fC*xk - fE;
-  if (TMath::Abs(c2)>1.) return 0;  
+  if (TMath::Abs(c2)>0.9) 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;
index 579d4f0..d592e91 100644 (file)
@@ -15,16 +15,59 @@ class AliTRDcluster;
 class AliTPCtrack;
 class AliESDtrack;
 
-const unsigned kMAXCLUSTERSPERTRACK=210; 
+const unsigned kMAX_CLUSTERS_PER_TRACK=210; 
+
+class AliTRDtracklet :public TObject{
+  friend class AliTRDtrack;
+ public:
+  AliTRDtracklet();
+  void Set(Float_t x, Float_t y, Float_t z, Float_t alpha, Float_t error2){fX=x; fY=y; fZ=z; fAlpha=alpha; fSigma2= error2;}
+  void SetP0(Float_t p0){fP0=p0;}
+  void SetP1(Float_t p1){fP1=p1;}
+  void SetN(Int_t n){fNFound=n;}
+  void SetNCross(Int_t nc){fNCross=nc;}
+  void SetPlane(Int_t plane){fPlane=plane;}
+  void SetSigma2(Float_t sigma2){fExpectedSigma2=sigma2;}
+  void SetChi2(Float_t chi2){fChi2=chi2;}
+  void SetTilt(Float_t tilt){fTilt=tilt;}
+  Float_t GetX() const { return fX;}
+  Float_t GetY() const { return fY;}
+  Float_t GetZ() const {return fZ;}
+  Float_t GetAlpha() const { return fAlpha;}
+  Float_t GetTrackletSigma2() const { return fSigma2;}
+  //
+  Float_t GetP0() const {return fP0;}
+  Float_t GetP1() const {return fP1;}
+  Int_t GetN() const {return fNFound;}
+  Int_t GetNCross() const {return fNCross;}  
+  Int_t GetPlane() const {return fPlane;}
+  Float_t GetClusterSigma2() const {return fExpectedSigma2;}
+  Float_t GetChi2() const {return fChi2;}
+  Float_t GetTilt() const {return fTilt;}
+ protected:
+  Float_t fY;                 // y position
+  Float_t fZ;                 // z position
+  Float_t fX;                 // x position
+  Float_t fAlpha;             // rotation angle
+  Float_t fSigma2;            // expected error of tracklet position
+  Float_t fP0;                // offset in y
+  Float_t fP1;                // offset in tangent
+  Int_t   fNFound;            // number of found clusters
+  Int_t   fNCross;            // number of crosses
+  Int_t   fPlane;             // plane number
+  Float_t fExpectedSigma2;    // expected sigma of residual distribution of clusters
+  Float_t fChi2;              // chi2 of the tracklet
+  Float_t fTilt;              // tilt factor 
+  ClassDef(AliTRDtracklet,2)
+};
 
-class AliTRDtrack : public AliKalmanTrack {
 
-  //////////////////////////////////////////////////////////////////
-  // Represents a reconstructed TRD track                         //
-  //////////////////////////////////////////////////////////////////
 
-  friend class AliTRDtracker;
 
+class AliTRDtrack : public AliKalmanTrack {
+
+// Represents reconstructed TRD track
+  friend class AliTRDtracker;
 public:
 
    AliTRDtrack():AliKalmanTrack(){fBackupTrack=0;}
@@ -36,7 +79,7 @@ public:
    ~AliTRDtrack();
    Int_t    Compare(const TObject *o) const;
    void     CookdEdx(Double_t low=0.05, Double_t up=0.55);   
-   Float_t  StatusForTOF();
+   Float_t    StatusForTOF();
    Double_t GetAlpha() const {return fAlpha;}
    Int_t    GetSector() const {
      //if (fabs(fAlpha) < AliTRDgeometry::GetAlpha()/2) return 0;
@@ -49,8 +92,8 @@ public:
    void     GetCovariance(Double_t cov[15]) const;  
    Double_t GetdEdx()  const {return fdEdx;}
    Double_t GetPIDsignal()  const {return GetdEdx();}
-   Float_t  GetPIDsignals(Int_t i) const {return fdEdxPlane[i];}
-   Int_t    GetPIDTimBin(Int_t i) const {return fTimBinPlane[i];}
+   Float_t GetPIDsignals(Int_t i) const {return fdEdxPlane[i];}
+   Int_t  GetPIDTimBin(Int_t i) const {return fTimBinPlane[i];}
    Double_t GetEta()   const {return fE;}
 
    void     GetExternalCovariance(Double_t cov[15]) const ;   
@@ -62,8 +105,7 @@ public:
    Double_t GetP()     const {  
      return TMath::Abs(GetPt())*sqrt(1.+GetTgl()*GetTgl());
    }
-   virtual  Double_t GetPredictedChi2(const AliCluster*) const { return 0.0; };
-   Double_t GetPredictedChi2(const AliTRDcluster *c, Double_t h01) const ;
+   Double_t GetPredictedChi2(const AliTRDcluster*, Double_t h01) const ;
    Double_t GetPt()    const {return 1./Get1Pt();}   
    void     GetPxPyPz(Double_t &px, Double_t &py, Double_t &pz) const ;
    void     GetGlobalXYZ(Double_t &x, Double_t &y, Double_t &z) const ;
@@ -79,29 +121,29 @@ public:
    Double_t GetZ()    const {return fZ;}
    UInt_t * GetBackupIndexes()  {return fIndexBackup;}
    UInt_t * GetIndexes()  {return fIndex;}
-   //-----------------------------------------------------------------
-   // This function calculates the Y-coordinate of a track at the plane x=xk.
-   // Needed for matching with the TOF (I.Belikov)
    Double_t GetYat(Double_t xk) const {     
+//-----------------------------------------------------------------
+// This function calculates the Y-coordinate of a track at the plane x=xk.
+// Needed for matching with the TOF (I.Belikov)
+//-----------------------------------------------------------------
       Double_t c1=fC*fX - fE, r1=TMath::Sqrt(1.- c1*c1);
       Double_t c2=fC*xk - fE, r2=TMath::Sqrt(1.- c2*c2);
       return fY + (xk-fX)*(c1+c2)/(r1+r2);
    }
-   //-----------------------------------------------------------------
-   Int_t    GetProlongation(Double_t xk, Double_t &y, Double_t &z);
+   Int_t GetProlongation(Double_t xk, Double_t &y, Double_t &z);
 
-   void     SetStop(Bool_t stop) {fStopped=stop;}
-   Bool_t   GetStop() const {return fStopped;}
+   void SetStop(Bool_t stop) {fStopped=stop;}
+   Bool_t GetStop() const {return fStopped;}
 
    Int_t    PropagateTo(Double_t xr, Double_t x0=8.72, Double_t rho=5.86e-3);
    void     ResetCovariance();   
    void     ResetCovariance(Float_t mult);   
-   void     ResetClusters() { SetChi2(0.); SetNumberOfClusters(0); }
+   void ResetClusters() { SetChi2(0.); SetNumberOfClusters(0); }
    Int_t    Rotate(Double_t angle);
 
    void     SetdEdx(Float_t dedx) {fdEdx=dedx;}  
-   void     SetPIDsignals(Float_t dedx, Int_t i) {fdEdxPlane[i]=dedx;}
-   void     SetPIDTimBin(Int_t timbin, Int_t i) {fTimBinPlane[i]=timbin;}
+   void SetPIDsignals(Float_t dedx, Int_t i) {fdEdxPlane[i]=dedx;}
+   void  SetPIDTimBin(Int_t timbin, Int_t i) {fTimBinPlane[i]=timbin;}
    void     SetLikelihoodElectron(Float_t l) { fLhElectron = l; };  
 
    void     SetSampledEdx(Float_t q, Int_t i) {
@@ -118,31 +160,34 @@ public:
 
    void     SetSeedLabel(Int_t lab) { fSeedLab=lab; }
 
-   Int_t    Update(const AliCluster*, Double_t, UInt_t) { return 0; };
    Int_t    Update(const AliTRDcluster* c, Double_t chi2, UInt_t i, 
                    Double_t h01);
    Int_t    UpdateMI(const AliTRDcluster* c, Double_t chi2, UInt_t i, 
                    Double_t h01, Int_t plane);
+   Int_t    UpdateMI(const AliTRDtracklet & tracklet);
 
-  void      AddNWrong() {fNWrong++;}
+  //
+  void AddNWrong() {fNWrong++;}
   
-  Int_t     GetNWrong() const {return fNWrong;}
-  Int_t     GetNRotate() const {return fNRotate;}
-  Int_t     GetNCross() const {return fNCross;}
-  void      IncCross() {fNCross++; if (fBackupTrack) fBackupTrack->IncCross();}
+  Int_t GetNWrong() const {return fNWrong;}
+  Int_t GetNRotate() const {return fNRotate;}
+  Int_t GetNCross() const {return fNCross;}
+  void  IncCross() {fNCross++; if (fBackupTrack) fBackupTrack->IncCross();}
   AliTRDtrack *  GetBackupTrack(){return fBackupTrack;}
-  void      MakeBackupTrack();
+  void    MakeBackupTrack();
+  //
+
 
 protected:
 
-   Int_t    fSeedLab;               // track label taken from seeding  
-   Float_t  fdEdx;                  // dE/dx 
-   Float_t  fdEdxPlane[kNPlane];    // dE/dx from all 6 planes
-   Int_t    fTimBinPlane[kNPlane];  // time bin of Max cluster from all 6 planes
+   Int_t    fSeedLab;     // track label taken from seeding  
+   Float_t  fdEdx;        // dE/dx 
+   Float_t  fdEdxPlane[kNPlane];  // dE/dx from all 6 planes
+   Int_t  fTimBinPlane[kNPlane];  // time bin of Max cluster from all 6 planes
 
-   Double_t fAlpha;         // rotation angle
-   Double_t fX;             // running local X-coordinate of the track (time bin)
-   Bool_t   fStopped;       // track stop indication
+   Double_t fAlpha;       // rotation angle
+   Double_t fX;           // running local X-coordinate of the track (time bin)
+   Bool_t   fStopped;     // track stop indication
 
    Double_t fY;             // Y-coordinate of the track
    Double_t fZ;             // Z-coordinate of the track
@@ -156,25 +201,23 @@ protected:
    Double_t fCty, fCtz, fCte, fCtt;       // track
    Double_t fCcy, fCcz, fCce, fCct, fCcc; // parameters   
    
-   UInt_t  fIndex[kMAXCLUSTERSPERTRACK];       // global indexes of clusters  
-   UInt_t  fIndexBackup[kMAXCLUSTERSPERTRACK]; // backup indexes of clusters - used in iterations
-   Float_t fdQdl[kMAXCLUSTERSPERTRACK];        // cluster amplitudes corrected 
-                                               // for track angles    
+   UInt_t  fIndex[kMAX_CLUSTERS_PER_TRACK];  // global indexes of clusters  
+   UInt_t  fIndexBackup[kMAX_CLUSTERS_PER_TRACK]; //backup indexes of clusters - used in iterations
+   Float_t fdQdl[kMAX_CLUSTERS_PER_TRACK];   // cluster amplitudes corrected 
+                                             // for track angles    
                            
    Float_t fLhElectron;    // Likelihood to be an electron    
-   Int_t   fNWrong;        // number of wrong clusters
-   Int_t   fNRotate;       // number of rotation
-   Int_t   fNCross;        // number of the cross materials
-   Int_t   fNExpected;     // expected number of cluster
-   Int_t   fNLast;         // number of clusters in last 2 layers
-   Int_t   fNExpectedLast; // number of expected clusters on last 2 layers
-   Int_t   fNdedx;         // number of clusters for dEdx measurment
-   Float_t fChi2Last;      // chi2 in the  last 2 layers
-
+   Int_t fNWrong;    // number of wrong clusters
+   Int_t fNRotate;   // number of rotation
+   Int_t fNCross;     // number of the cross materials
+   Int_t fNExpected;  //expected number of cluster
+   Int_t fNLast;      //number of clusters in last 2 layers
+   Int_t fNExpectedLast; //number of expected clusters on last 2 layers
+   Int_t      fNdedx;      //number of clusters for dEdx measurment
+   Float_t fChi2Last;      //chi2 in the  last 2 layers
+   AliTRDtracklet fTracklets[6]; //tracklets
    AliTRDtrack * fBackupTrack; //! backup track
-
    ClassDef(AliTRDtrack,2) // TRD reconstructed tracks
-
 };                     
 
 
index f287b70..33e1175 100644 (file)
 #include "AliTRDgeometry.h"
 #include "AliTRDparameter.h"
 #include "AliTRDpadPlane.h"
-#include "AliTRDgeometryFull.h"
+#include "AliTRDgeometryDetail.h"
 #include "AliTRDcluster.h" 
 #include "AliTRDtrack.h"
 #include "AliESD.h"
 
+#include "TTreeStream.h"
+#include "TGraph.h"
 #include "AliTRDtracker.h"
+//
 
 ClassImp(AliTRDtracker) 
 
@@ -64,10 +67,26 @@ ClassImp(AliTRDtracker)
 
   const  Double_t    AliTRDtracker::fgkMaxChi2            = 12.; 
 
+//   const  Double_t    AliTRDtracker::fgkOffset             = -0.012;
+//   const  Double_t    AliTRDtracker::fgkOffsetX            = 0.35;
+//   const  Double_t    AliTRDtracker::fgkCoef               = 0.00;
+//   const  Double_t    AliTRDtracker::fgkMean               = 8.;
+//   const  Double_t    AliTRDtracker::fgkDriftCorrection    = 1.07;
+//   const  Double_t    AliTRDtracker::fgkExB                = 0.072;
+
+  const  Double_t    AliTRDtracker::fgkOffset             = -0.015;
+const  Double_t    AliTRDtracker::fgkOffsetX            = 0.26;       // "time offset"  
+  const  Double_t    AliTRDtracker::fgkCoef               = 0.0096;   // angular shift 
+  const  Double_t    AliTRDtracker::fgkMean               = 0.;
+  const  Double_t    AliTRDtracker::fgkDriftCorrection    = 1.04;   // drift coefficient correction
+  const  Double_t    AliTRDtracker::fgkExB                = 0.072;  // ExB angle - for error parameterization
+
+
+//     poscorrection =  fgkCoef*(GetLocalTimeBin() - fgkMean)+fgkOffset; 
+
 const Int_t AliTRDtracker::fgkFirstPlane = 5;
 const Int_t AliTRDtracker::fgkLastPlane = 17;
 
-
 //____________________________________________________________________
 AliTRDtracker::AliTRDtracker():AliTracker(),
                               fGeom(0),
@@ -91,6 +110,7 @@ AliTRDtracker::AliTRDtracker():AliTracker(),
   for(Int_t i=0;i<kTrackingSectors;i++) fTrSec[i]=0;
   for(Int_t j=0;j<5;j++)
     for(Int_t k=0;k<18;k++) fHoles[j][k]=kFALSE;
+  fDebugStreamer = 0;
 } 
 //____________________________________________________________________
 AliTRDtracker::AliTRDtracker(const TFile *geomfile):AliTracker()
@@ -127,7 +147,7 @@ AliTRDtracker::AliTRDtracker(const TFile *geomfile):AliTracker()
     printf("AliTRDtracker::AliTRDtracker(): can't find TRD geometry!\n");
     //printf("The DETAIL TRD geometry will be used\n");
     //fGeom = new AliTRDgeometryDetail();
-    fGeom = new AliTRDgeometryFull();
+    fGeom = new AliTRDgeometryDetail();
     fGeom->SetPHOShole();
     fGeom->SetRICHhole();    
   } 
@@ -135,12 +155,14 @@ AliTRDtracker::AliTRDtracker(const TFile *geomfile):AliTracker()
   if (!fPar) {  
     printf("AliTRDtracker::AliTRDtracker(): can't find TRD parameter!\n");
     printf("The DEFAULT TRD parameter will be used\n");
-    fPar = new AliTRDparameter();
+    fPar = new AliTRDparameter("Pica","Vyjebana");
   }
+  fPar = new AliTRDparameter("Pica","Vyjebana");
   fPar->Init();
 
   savedir->cd();  
 
+
   //  fGeom->SetT0(fTzero);
 
   fNclusters = 0;
@@ -157,9 +179,9 @@ AliTRDtracker::AliTRDtracker(const TFile *geomfile):AliTracker()
       fHoles[icham][trS]=fGeom->IsHole(0,icham,geomS);
     }
   }
-
   AliTRDpadPlane *padPlane = fPar->GetPadPlane(0,0);
-  Float_t tiltAngle = TMath::Abs(padPlane->GetTiltingAngle()); 
+  Float_t tiltAngle = TMath::Abs(padPlane->GetTiltingAngle());
+  //  Float_t tiltAngle = TMath::Abs(fPar->GetTiltingAngle()); 
   if(tiltAngle < 0.1) {
     fNoTilt = kTRUE;
   }
@@ -175,8 +197,9 @@ AliTRDtracker::AliTRDtracker(const TFile *geomfile):AliTracker()
   Double_t dxAmp = (Double_t) fGeom->CamHght();   // Amplification region
   Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region
 
-  Double_t dx = (Double_t) fPar->GetDriftVelocity()
+  Double_t dx = fgkDriftCorrection*(Double_t) fPar->GetDriftVelocity()
                          / fPar->GetSamplingFrequency();
+
   Int_t tbAmp = fPar->GetTimeBefore();
   Int_t maxAmp = (Int_t) ((dxAmp+0.000001)/dx);
   if(kTRUE) maxAmp = 0;  // intentional until we change the parameter class 
@@ -190,6 +213,8 @@ AliTRDtracker::AliTRDtracker(const TFile *geomfile):AliTracker()
   fMaxGap = (Int_t) (fTimeBinsPerPlane * fGeom->Nplan() * fgkSkipDepth);
 
   fVocal = kFALSE;
+  
+  fDebugStreamer = new TTreeSRedirector("TRDdebug.root");
 
   savedir->cd();
 }   
@@ -219,7 +244,10 @@ AliTRDtracker::~AliTRDtracker()
   for(Int_t geomS = 0; geomS < kTrackingSectors; geomS++) {
     delete fTrSec[geomS];
   }
-
+  if (fDebugStreamer) {    
+    //fDebugStreamer->Close();
+    delete fDebugStreamer;
+  }
 }   
 
 //_____________________________________________________________________
@@ -300,7 +328,7 @@ inline Double_t f3trd(Double_t x1,Double_t y1,
 }            
 
 
-AliTRDcluster * AliTRDtracker::GetCluster(AliTRDtrack * track, Int_t plane, Int_t timebin){
+AliTRDcluster * AliTRDtracker::GetCluster(AliTRDtrack * track, Int_t plane, Int_t timebin, UInt_t &index){
   //
   //try to find cluster in the backup list
   //
@@ -314,6 +342,7 @@ AliTRDcluster * AliTRDtracker::GetCluster(AliTRDtrack * track, Int_t plane, Int_
     Int_t iplane = fGeom->GetPlane(cli->GetDetector());
     if (iplane==plane) {
       cl = cli;
+      index = indexes[i];
       break;
     }
   }
@@ -360,7 +389,7 @@ Int_t AliTRDtracker::Clusters2Tracks(AliESD* event)
     if ( (status & AliESDtrack::kTRDout ) == 0 ) continue;
     if ( (status & AliESDtrack::kTRDin) != 0 ) continue;
     nseed++;
-
+    
     AliTRDtrack* seed2 = new AliTRDtrack(*seed);
     //seed2->ResetCovariance(); 
     AliTRDtrack *pt = new AliTRDtrack(*seed2,seed2->GetAlpha());
@@ -512,11 +541,13 @@ Int_t AliTRDtracker::PropagateBack(AliESD* event) {
        Bool_t isGold = kFALSE;
        
        if (track->GetChi2()/track->GetNumberOfClusters()<5) {  //full gold track
-         seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
+         // seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
+          if (track->GetBackupTrack()) seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
          isGold = kTRUE;
        }
        if (!isGold && track->GetNCross()==0&&track->GetChi2()/track->GetNumberOfClusters()<7){ //almost gold track
-         seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
+         //      seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
+         if (track->GetBackupTrack()) seed->UpdateTrackParams(track->GetBackupTrack(), AliESDtrack::kTRDbackup);
          isGold = kTRUE;
        }
        if (!isGold && track->GetBackupTrack()){
@@ -526,11 +557,14 @@ Int_t AliTRDtracker::PropagateBack(AliESD* event) {
            isGold = kTRUE;
          }
        }
+       if (track->StatusForTOF()>0 &&track->fNCross==0 && Float_t(track->fN)/Float_t(track->fNExpected)>0.4){
+         seed->UpdateTrackParams(track, AliESDtrack::kTRDbackup);
+       }
       }
     }
     //
     //Propagation to the TOF (I.Belikov)
-    CookdEdxTimBin(*track);
+    
     if (track->GetStop()==kFALSE){
       
       Double_t xtof=371.;
@@ -570,7 +604,7 @@ Int_t AliTRDtracker::PropagateBack(AliESD* event) {
            seed->SetTRDsignals(track->GetPIDsignals(i),i);
            seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
         }
-       seed->SetTRDtrack(new AliTRDtrack(*track));
+       //      seed->SetTRDtrack(new AliTRDtrack(*track));
        if (track->GetNumberOfClusters()>foundMin) found++;
       }
     }else{
@@ -581,13 +615,32 @@ Int_t AliTRDtracker::PropagateBack(AliESD* event) {
            seed->SetTRDsignals(track->GetPIDsignals(i),i);
            seed->SetTRDTimBin(track->GetPIDTimBin(i),i);
         }
-       seed->SetTRDtrack(new AliTRDtrack(*track));
+       //seed->SetTRDtrack(new AliTRDtrack(*track));
        found++;
       }
     }
-    seed->SetTRDQuality(track->StatusForTOF());
+    seed->SetTRDQuality(track->StatusForTOF());    
+    //
+    // Debug part of tracking
+    TTreeSRedirector& cstream = *fDebugStreamer;
+    Int_t eventNr = event->GetEventNumber();
+    if (track->GetBackupTrack()){
+      cstream<<"Tracks"<<
+       "EventNr="<<eventNr<<
+       "ESD.="<<seed<<
+       "trd.="<<track<<
+       "trdback.="<<track->GetBackupTrack()<<  
+       "\n";
+    }else{
+      cstream<<"Tracks"<<
+       "EventNr="<<eventNr<<
+       "ESD.="<<seed<<
+       "trd.="<<track<<
+       "trdback.="<<track<<
+       "\n";
+    }
     delete track;
-    
+    //
     //End of propagation to the TOF
     //if (foundClr>foundMin)
     //  seed->UpdateTrackParams(track, AliESDtrack::kTRDout);
@@ -598,6 +651,9 @@ Int_t AliTRDtracker::PropagateBack(AliESD* event) {
   cerr<<"Number of seeds: "<<fNseeds<<endl;  
   cerr<<"Number of back propagated TRD tracks: "<<found<<endl;
 
+  //  MakeSeedsMI(3,5); //new seeding
+
+
   fSeeds->Clear(); fNseeds=0;
   delete [] index;
   delete [] quality;
@@ -641,22 +697,22 @@ Int_t AliTRDtracker::RefitInward(AliESD* event)
       continue;
     }
     nseed++;    
-    if (1/seed2.Get1Pt()>5.&& seed2.GetX()>260.) {
-      Double_t oldx = seed2.GetX();
-      seed2.PropagateTo(500.);
-      seed2.ResetCovariance(1.);
-      seed2.PropagateTo(oldx);
-    }
-    else{
-      seed2.ResetCovariance(5.); 
-    }
+//     if (1/seed2.Get1Pt()>1.5&& seed2.GetX()>260.) {
+//       Double_t oldx = seed2.GetX();
+//       seed2.PropagateTo(500.);
+//       seed2.ResetCovariance(1.);
+//       seed2.PropagateTo(oldx);
+//     }
+//     else{
+//       seed2.ResetCovariance(5.); 
+//     }
 
     AliTRDtrack *pt = new AliTRDtrack(seed2,seed2.GetAlpha());
     UInt_t * indexes2 = seed2.GetIndexes();
-//     for (Int_t i=0;i<kNPlane;i++) {
-//       pt->SetPIDsignals(seed2.GetPIDsignals(i),i);
-//       pt->SetPIDTimBin(seed2.GetPIDTimBin(i),i);
-//     }
+    for (Int_t i=0;i<kNPlane;i++) {
+      pt->SetPIDsignals(seed2.GetPIDsignals(i),i);
+      pt->SetPIDTimBin(seed2.GetPIDTimBin(i),i);
+    }
 
     UInt_t * indexes3 = pt->GetBackupIndexes();
     for (Int_t i=0;i<200;i++) {
@@ -669,17 +725,18 @@ Int_t AliTRDtracker::RefitInward(AliESD* event)
     if (t.GetNumberOfClusters() >= foundMin) {
       //      UseClusters(&t);
       //CookLabel(pt, 1-fgkLabelFraction);
-      //      t.CookdEdx();
+      t.CookdEdx();
+      CookdEdxTimBin(t);
     }
     found++;
 //    cout<<found<<'\r';     
 
     if(PropagateToTPC(t)) {
       seed->UpdateTrackParams(pt, AliESDtrack::kTRDrefit);
-     //  for (Int_t i=0;i<kNPlane;i++) {
-//         seed->SetTRDsignals(pt->GetPIDsignals(i),i);
-//         seed->SetTRDTimBin(pt->GetPIDTimBin(i),i);
-//       }
+      for (Int_t i=0;i<kNPlane;i++) {
+        seed->SetTRDsignals(pt->GetPIDsignals(i),i);
+        seed->SetTRDTimBin(pt->GetPIDTimBin(i),i);
+      }
     }else{
       //if not prolongation to TPC - propagate without update
       AliTRDtrack* seed2 = new AliTRDtrack(*seed);
@@ -690,10 +747,10 @@ Int_t AliTRDtracker::RefitInward(AliESD* event)
         pt2->CookdEdx(0.,1.);
         CookdEdxTimBin(*pt2);
        seed->UpdateTrackParams(pt2, AliESDtrack::kTRDrefit);
-       //  for (Int_t i=0;i<kNPlane;i++) {
-//           seed->SetTRDsignals(pt2->GetPIDsignals(i),i);
-//           seed->SetTRDTimBin(pt2->GetPIDTimBin(i),i);
-//         }
+        for (Int_t i=0;i<kNPlane;i++) {
+          seed->SetTRDsignals(pt2->GetPIDsignals(i),i);
+          seed->SetTRDTimBin(pt2->GetPIDTimBin(i),i);
+        }
       }
       delete pt2;
     }  
@@ -854,7 +911,7 @@ Int_t AliTRDtracker::FollowProlongation(AliTRDtrack& t, Int_t rf)
        Int_t plane = fGeom->GetPlane(cl0->GetDetector());
        if (plane>lastplane) continue;
        Int_t timebin = cl0->GetLocalTimeBin();
-       AliTRDcluster * cl2= GetCluster(&t,plane, timebin);
+       AliTRDcluster * cl2= GetCluster(&t,plane, timebin,index);
        if (cl2) {
          cl =cl2;      
          Double_t h01 = GetTiltFactor(cl);
@@ -867,44 +924,6 @@ Int_t AliTRDtracker::FollowProlongation(AliTRDtrack& t, Int_t rf)
          continue;
        }             
 
-       /*
-       if(!cl){
-         Int_t maxn = timeBin;
-         for (Int_t i=timeBin.Find(y-road); i<maxn; i++) {
-           AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
-           if (c->GetY() > y+road) break;
-           if (c->IsUsed() > 0) continue;
-           if((c->GetZ()-z)*(c->GetZ()-z) > 3 * sz2) continue;
-           
-           Double_t h01 = GetTiltFactor(c);
-           Double_t chi2=t.GetPredictedChi2(c,h01);
-           
-           if (chi2 > maxChi2) continue;
-           maxChi2=chi2;
-           cl=c;
-           index=timeBin.GetIndex(i);
-         }               
-       }
-
-        if(!cl) {
-         Int_t maxn = timeBin;
-          for (Int_t i=timeBin.Find(y-road); i<maxn; i++) {
-            AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
-            
-            if (c->GetY() > y+road) break;
-            if (c->IsUsed() > 0) continue;
-            if((c->GetZ()-z)*(c->GetZ()-z) > 12 * sz2) continue;
-            
-            Double_t h01 = GetTiltFactor(c);
-            Double_t chi2=t.GetPredictedChi2(c, h01);
-            
-            if (chi2 > maxChi2) continue;
-            maxChi2=chi2;
-            cl=c;
-            index=timeBin.GetIndex(i);
-          }
-        } 
-       */       
         if (cl) {
          
           wYclosest = cl->GetY();
@@ -912,9 +931,7 @@ Int_t AliTRDtracker::FollowProlongation(AliTRDtrack& t, Int_t rf)
           Double_t h01 = GetTiltFactor(cl);
 
           if (cl->GetNPads()<5) 
-           t.SetSampledEdx(cl->GetQ()/dx); 
-         //printf("Track   position\t%f\t%f\t%f\n",t.GetX(),t.GetY(),t.GetZ());
-         //printf("Cluster position\t%d\t%f\t%f\n",cl->GetLocalTimeBin(),cl->GetY(),cl->GetZ());
+           t.SetSampledEdx(TMath::Abs(cl->GetQ()/dx)); 
          Int_t det = cl->GetDetector();    
          Int_t plane = fGeom->GetPlane(det);
 
@@ -945,13 +962,6 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t)
   // layers confirms prolongation if a close cluster is found. 
   // Returns the number of clusters expected to be found in sensitive layers
 
-
-  Float_t  wIndex, wTB, wChi2;
-  Float_t  wYrt, wYclosest, wYcorrect, wYwindow;
-  Float_t  wZrt, wZclosest, wZcorrect, wZwindow;
-  Float_t  wSigmaC2, wSigmaTgl2, wSigmaY2, wSigmaZ2;
-
-  Int_t trackIndex = t.GetLabel();  
   Int_t tryAgain=fMaxGap;
 
   Double_t alpha=t.GetAlpha();
@@ -963,7 +973,8 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t)
   for (Int_t i=0;i<1000;i++) clusters[i]=-1;
 
   Int_t outerTB = fTrSec[0]->GetOuterTimeBin();
-  Double_t radLength, rho, x, dx, y, ymax = 0, z;
+  //Double_t radLength, rho, x, dx, y, ymax = 0, z;
+  Double_t radLength, rho, x, dx, y, z;
   Bool_t lookForCluster;
 
   Int_t expectedNumberOfClusters = 0;
@@ -971,18 +982,15 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t)
 
   alpha=AliTRDgeometry::GetAlpha();  // note: change in meaning
 
-  Int_t nRefPlane = fgkFirstPlane;
-  Bool_t isNewLayer = kFALSE; 
-
-  Double_t chi2;
-  Double_t minDY;
-  Int_t zone =-10;
+  //  Int_t zone =0;
   Int_t nr;
+  Float_t ratio0=0;
+  AliTRDtracklet tracklet;
+  //
   for (nr=fTrSec[0]->GetLayerNumber(t.GetX()); nr<outerTB+1; nr++) { 
-    
+   
     y = t.GetY(); 
     z = t.GetZ();
-
     // first propagate to the outer surface of the current time bin 
 
     s = t.GetSector();
@@ -992,210 +1000,83 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t)
     z = t.GetZ();
 
     if(!t.PropagateTo(x,radLength,rho)) break;
-    //    if (!AdjustSector(&t)) break;
     //
     // MI -fix untill correct material desription will be implemented
     //
-    Float_t angle =  t.GetAlpha();  // MI - if rotation - we go through the material 
-    if (!AdjustSector(&t)) break;
-    Int_t cross = kFALSE;
-    Int_t crosz = kFALSE;
-    if (TMath::Abs(angle -  t.GetAlpha())>0.000001) cross = kTRUE; //better to stop track
-    Int_t currentzone = fTrSec[s]->GetLayer(nr)->GetZone(z);
-    if (currentzone==-10) {cross = kTRUE,crosz=kTRUE;}  // we are in the frame
-    if (currentzone>-10){   // layer knows where we are
-      if (zone==-10) zone = currentzone;
-      if (zone!=currentzone) {
-       cross=kTRUE;  
-       crosz=kTRUE;
-      }
-    }
-    if (TMath::Abs(t.GetSnp())>0.8 && t.GetBackupTrack()==0) t.MakeBackupTrack();
-    if (cross) {
-      if (t.GetNCross()==0 && t.GetBackupTrack()==0) t.MakeBackupTrack();
-      t.IncCross();
-      if (t.GetNCross()>4) break;
-    }
+    //Int_t nrotate = t.GetNRotate();
+    if (!AdjustSector(&t)) break;    
     //
     //
-    s = t.GetSector();
-    if (!t.PropagateTo(x,radLength,rho)) break;
-
     y = t.GetY();
     z = t.GetZ();
-
-    // Barrel Tracks [SR, 04.04.2003]
-
     s = t.GetSector();
-    if (fTrSec[s]->GetLayer(nr)->IsSensitive() != 
-        fTrSec[s]->GetLayer(nr+1)->IsSensitive() ) {
-
-//      if (IsStoringBarrel()) StoreBarrelTrack(&t, nRefPlane++, kTrackBack);
-    }
-
-    if (fTrSec[s]->GetLayer(nr-1)->IsSensitive() && 
-          ! fTrSec[s]->GetLayer(nr)->IsSensitive()) {
-      isNewLayer = kTRUE;
-    } else {isNewLayer = kFALSE;}
-
-    y = t.GetY();
-    z = t.GetZ();
 
     // now propagate to the middle plane of the next time bin 
     fTrSec[s]->GetLayer(nr+1)->GetPropagationParameters(y,z,dx,rho,radLength,lookForCluster);
-    if (crosz){
-      rho = 1000*2.7; radLength = 24.01;  //TEMPORARY - aluminium in between z - will be detected using GeoModeler in future versions
-    }
+//     if (nrotate!=t.GetNRotate()){
+//       rho = 1000*2.7; radLength = 24.01;  //TEMPORARY - aluminium in between z - will be detected using GeoModeler in future versions
+//     }
     x = fTrSec[s]->GetLayer(nr+1)->GetX(); 
     if(!t.PropagateTo(x,radLength,rho)) break;
     if (!AdjustSector(&t)) break;
     s = t.GetSector();
-    if(!t.PropagateTo(x,radLength,rho)) break;
+    //    if(!t.PropagateTo(x,radLength,rho)) break;
     
     if (TMath::Abs(t.GetSnp())>0.95) break;
 
     y = t.GetY();
     z = t.GetZ();
 
-    //    if(fVocal) printf("nr+1=%d, x %f, z %f, y %f, ymax %f\n",nr+1,x,z,y,ymax);
-    //    printf("label %d, pl %d, lookForCluster %d \n",
-    //     trackIndex, nr+1, lookForCluster);
-
     if(lookForCluster) {
-//       if (clusters[nr]==-1) {
-//     FindClusters(s,nr,nr+30,&t,clusters);
-//       }
+      if (clusters[nr]==-1) {
+       Float_t  ncl   = FindClusters(s,nr,nr+30,&t,clusters,tracklet);
+       ratio0 = ncl/Float_t(fTimeBinsPerPlane);
+       Float_t  ratio1 = Float_t(t.fN+1)/Float_t(t.fNExpected+1.);     
+       if (tracklet.GetChi2()<18.&&ratio0>0.8 && ratio1>0.6 && ratio0+ratio1>1.5 && t.GetNCross()==0 && TMath::Abs(t.GetSnp())<0.85&&t.fN>20){
+         t.MakeBackupTrack();                            // make backup of the track until is gold
+       }
+//     if (ncl>4){
+//       t.PropagateTo(tracklet.GetX());
+//       t.UpdateMI(tracklet);
+//       nr = fTrSec[0]->GetLayerNumber(t.GetX())+1;
+//       continue;
+//     }
+      }
 
       expectedNumberOfClusters++;       
       t.fNExpected++;
       if (t.fX>345) t.fNExpectedLast++;
-      wIndex = (Float_t) t.GetLabel();
-      wTB = fTrSec[s]->GetLayer(nr+1)->GetTimeBinIndex();
 
       AliTRDpropagationLayer& timeBin=*(fTrSec[s]->GetLayer(nr+1));
       Double_t sy2=ExpectedSigmaY2(t.GetX(),t.GetTgl(),t.GetPt());
-      Double_t sz2=ExpectedSigmaZ2(t.GetX(),t.GetTgl());
-      if((t.GetSigmaY2() + sy2) < 0) break;
+      if((t.GetSigmaY2() + sy2) < 0) {
+       printf("problem\n");
+       break;
+      }
       Double_t road = 10.*sqrt(t.GetSigmaY2() + sy2); 
-      Double_t y=t.GetY(), z=t.GetZ();
-
-      wYrt = (Float_t) y;
-      wZrt = (Float_t) z;
-      wYwindow = (Float_t) road;
-      wSigmaC2 = (Float_t) t.GetSigmaC2();
-      wSigmaTgl2    = (Float_t) t.GetSigmaTgl2();
-      wSigmaY2 = (Float_t) t.GetSigmaY2();
-      wSigmaZ2 = (Float_t) t.GetSigmaZ2();
-      wChi2 = -1;            
       
       if (road>fgkWideRoad) {
-        if (t.GetNumberOfClusters()>4)
-          cerr<<t.GetNumberOfClusters()
-              <<"FindProlongation warning: Too broad road !\n";
         return 0;
       }      
 
       AliTRDcluster *cl=0;
       UInt_t index=0;
-
       Double_t maxChi2=fgkMaxChi2;
-
-      if (isNewLayer) { 
-        road = 3 * road;
-        //sz2 = 3 * sz2;
-        maxChi2 = 10 * fgkMaxChi2;
-      }
-      
-      if (nRefPlane == fgkFirstPlane) maxChi2 = 20 * fgkMaxChi2; 
-      if (nRefPlane == fgkFirstPlane+2) maxChi2 = 15 * fgkMaxChi2;
-      if (t.GetNRotate() > 0) maxChi2 = 3 * maxChi2;
       
-
-      wYclosest = 12345678;
-      wYcorrect = 12345678;
-      wZclosest = 12345678;
-      wZcorrect = 12345678;
-      wZwindow  = TMath::Sqrt(2.25 * 12 * sz2);   
-
-      // Find the closest correct cluster for debugging purposes
-      if (timeBin&&fVocal) {
-        minDY = 1000000;
-        for (Int_t i=0; i<timeBin; i++) {
-          AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
-          if((c->GetLabel(0) != trackIndex) &&
-             (c->GetLabel(1) != trackIndex) &&
-             (c->GetLabel(2) != trackIndex)) continue;
-          if(TMath::Abs(c->GetY() - y) > minDY) continue;
-          //minDY = TMath::Abs(c->GetY() - y);
-          minDY = c->GetY() - y;
-          wYcorrect = c->GetY();
-          wZcorrect = c->GetZ();
-
-          Double_t h01 = GetTiltFactor(c);
-          wChi2 = t.GetPredictedChi2(c, h01);
-        }
-      }                    
-
       // Now go for the real cluster search
 
       if (timeBin) {
-       /*
-         if (clusters[nr+1]>0) {
+       
+       if (clusters[nr+1]>0) {
          index = clusters[nr+1];
          cl    = (AliTRDcluster*)GetCluster(index);
          Double_t h01 = GetTiltFactor(cl);
           maxChi2=t.GetPredictedChi2(cl,h01);          
-         }
-       */
-
-       if (!cl){
-         Int_t maxn = timeBin;
-         for (Int_t i=timeBin.Find(y-road); i<maxn; i++) {
-           AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
-           if (c->GetY() > y+road) break;
-           if (c->IsUsed() > 0) continue;
-           if((c->GetZ()-z)*(c->GetZ()-z) > 3 * sz2) continue;
-
-           Double_t h01 = GetTiltFactor(c);
-           chi2=t.GetPredictedChi2(c,h01);
-           
-           if (chi2 > maxChi2) continue;
-           maxChi2=chi2;
-           cl=c;
-           index=timeBin.GetIndex(i);
-           
-          //check is correct
-           if((c->GetLabel(0) != trackIndex) &&
-              (c->GetLabel(1) != trackIndex) &&
-              (c->GetLabel(2) != trackIndex)) t.AddNWrong();
-         }             
        }
-        if(!cl) {
-         Int_t maxn = timeBin;
-          for (Int_t i=timeBin.Find(y-road); i<maxn; i++) {
-            AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
-            
-            if (c->GetY() > y+road) break;
-            if (c->IsUsed() > 0) continue;
-           //  if((c->GetZ()-z)*(c->GetZ()-z) > 2.25 * 12 * sz2) continue;
-           if((c->GetZ()-z)*(c->GetZ()-z) > 12. * sz2) continue;
-           //
-           //            
-            Double_t h01 = GetTiltFactor(c);
-            chi2=t.GetPredictedChi2(c,h01);
-            
-            if (chi2 > maxChi2) continue;
-            maxChi2=chi2;
-            cl=c;
-            index=timeBin.GetIndex(i);
-          }
-        }        
-        
+       
         if (cl) {
-          wYclosest = cl->GetY();
-          wZclosest = cl->GetZ();
          if (cl->GetNPads()<5) 
-           t.SetSampledEdx(cl->GetQ()/dx); 
+           t.SetSampledEdx(TMath::Abs(cl->GetQ()/dx)); 
           Double_t h01 = GetTiltFactor(cl);
          Int_t det = cl->GetDetector();    
          Int_t plane = fGeom->GetPlane(det);
@@ -1205,18 +1086,26 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t)
          }
          if(!t.UpdateMI(cl,maxChi2,index,h01,plane)) {
            if(!t.Update(cl,maxChi2,index,h01)) {
-             if(!tryAgain--) return 0;
+             //if(!tryAgain--) return 0;
            }
           }  
           else tryAgain=fMaxGap;
-        }
+         //
+         
+         if (cl->GetLocalTimeBin()==1&&t.fN>20 && float(t.fChi2)/float(t.fN)<5){
+           Float_t  ratio1 = Float_t(t.fN)/Float_t(t.fNExpected);      
+           if (tracklet.GetChi2()<18&&ratio0>0.8&&ratio1>0.6 &&ratio0+ratio1>1.5 && t.GetNCross()==0 && TMath::Abs(t.GetSnp())<0.85){
+             t.MakeBackupTrack();                            // make backup of the track until is gold
+           }
+         }
+         
+       }                       
         else {
-          if (tryAgain==0) break; 
-          tryAgain--;                                                                               
+         // if (tryAgain==0) break; 
+          //tryAgain--;                                                                               
         }
-
-        isNewLayer = kFALSE;
-
+       
+       
       }
     }  
   }
@@ -1225,8 +1114,8 @@ Int_t AliTRDtracker::FollowBackProlongation(AliTRDtrack& t)
   else
     t.SetStop(kFALSE);
   return expectedNumberOfClusters;
-
-
+  
+  
 }         
 
 //---------------------------------------------------------------------------
@@ -1318,7 +1207,7 @@ Int_t AliTRDtracker::Refit(AliTRDtrack& t, Int_t rf)
     AliTRDcluster *cl=(AliTRDcluster*)GetCluster(iCluster[nr-1]);
     Double_t h01 = GetTiltFactor(cl);
     Double_t chi2=t.GetPredictedChi2(cl, h01);
-    if (cl->GetNPads()<5) t.SetSampledEdx(cl->GetQ()/dx); 
+    if (cl->GetNPads()<5) t.SetSampledEdx(TMath::Abs(cl->GetQ()/dx)); 
 
       //t.SetSampledEdx(cl->GetQ()/dx,t.GetNumberOfClusters()); 
     t.Update(cl,chi2,iCluster[nr-1],h01);
@@ -1699,6 +1588,171 @@ void AliTRDtracker::MakeSeeds(Int_t inner, Int_t outer, Int_t turn)
     }
   }
 }            
+//__________________________________________________________________________
+void AliTRDtracker::MakeSeedsMI(Int_t /*inner*/, Int_t /*outer*/)
+{
+  //
+  // Creates  seeds using clusters between  position inner plane  and outer plane 
+  //
+
+  const Double_t maxtheta = 2;
+  const Double_t maxphi   = 1.5;
+  Int_t maxSec=AliTRDgeometry::kNsect;
+
+  //  
+  // find the maximal and minimal layer for the planes
+  // fucking "object oriented" geometry - find the time bin range for different planes
+  //
+  Int_t layers[6][2];
+  for (Int_t i=0;i<6;i++){layers[i][0]=10000; layers[i][1]=0;}
+
+  for (Int_t ns=0;ns<maxSec;ns++){
+    for (Int_t ilayer=0;ilayer<fTrSec[ns]->GetNumberOfLayers();ilayer++){
+      AliTRDpropagationLayer& layer=*(fTrSec[ns]->GetLayer(ilayer));
+      if (layer==0) continue;
+      Int_t det   = layer[0]->GetDetector();    
+      Int_t plane = fGeom->GetPlane(det);
+      if (ilayer<layers[plane][0]) layers[plane][0] = ilayer;
+      if (ilayer>layers[plane][1]) layers[plane][1] = ilayer;
+    }
+  }
+  //
+      
+  Int_t ilayer1 = layers[5][1];  // time bin in mplification region
+  Int_t ilayer2 = layers[3][1];  //
+  Int_t ilayerM = layers[4][1];  //
+  //    
+  Double_t x1 = fTrSec[0]->GetX(ilayer1);
+  Double_t x2 = fTrSec[0]->GetX(ilayer2);
+  Double_t xm = fTrSec[0]->GetX(ilayerM);
+  Double_t dist = x2-x1;
+  //  Int_t indexes1[20];
+  //Int_t indexes2[20];
+  AliTRDcluster *clusters1[15],*clusters2[15],*clustersM[15];
+  //
+  //      
+  for (Int_t ns=0; ns<maxSec; ns++) {
+    AliTRDpropagationLayer& layer1=*(fTrSec[ns]->GetLayer(ilayer1));  //select propagation layers
+    AliTRDpropagationLayer& layer2=*(fTrSec[ns]->GetLayer(ilayer2));
+    //  
+    for (Int_t icl1=0;icl1<layer1;icl1++){
+      AliTRDcluster *cl1 = layer1[icl1];
+      if (!cl1) continue;
+      Double_t y1 = cl1->GetY();
+      Double_t z1 = cl1->GetZ();
+      //
+      for (Int_t icl2=0;icl2<layer2;icl2++){
+       AliTRDcluster *cl2 = layer2[icl2];
+       if (!cl2) continue;
+       Double_t y2 = cl2->GetY();
+       Double_t z2 = cl2->GetZ();      
+       Double_t tanphi   = (y2-y1)/dist; 
+       Double_t tantheta = (z2-z1)/dist; 
+       if (TMath::Abs(tanphi)>maxphi) continue;
+       if (TMath::Abs(tantheta)>maxtheta) continue;
+       //
+       clusters1[0] = cl1;
+       clusters2[0] = cl2;
+       Double_t road =  0.5+TMath::Abs(tanphi)*1;
+       Int_t ncl=0;
+       Double_t sum1=0, sumx1=0,sum2x1=0,sumxy1=0, sumy1=0;
+       Double_t sum2=0, sumx2=0,sum2x2=0,sumxy2=0, sumy2=0;
+       //
+       for (Int_t dlayer=1;dlayer<15;dlayer++){
+         clusters1[dlayer]=0;
+         clusters2[dlayer]=0;
+         AliTRDpropagationLayer& layer1C=*(fTrSec[ns]->GetLayer(ilayer1-dlayer));  //select propagation layers
+         AliTRDpropagationLayer& layer2C=*(fTrSec[ns]->GetLayer(ilayer2-dlayer));  //
+         Double_t yy1 = y1+(tanphi)  *(layer1C.GetX()-x1);
+         Double_t zz1 = z1+(tantheta)*(layer1C.GetX()-x1);
+         Double_t yy2 = y1+(tanphi)  *(layer2C.GetX()-x1);
+         Double_t zz2 = z1+(tantheta)*(layer2C.GetX()-x1);
+         Int_t index1 = layer1C.FindNearestCluster(yy1,zz1,road);
+         Int_t index2 = layer2C.FindNearestCluster(yy2,zz2,road);
+         if (index1>=0) {
+           clusters1[dlayer]= (AliTRDcluster*)GetCluster(index1);
+           ncl++;
+           sum1++;
+           Double_t dx = layer1C.GetX()-x1;
+           sumx1 +=dx;
+           sum2x1+=dx*dx;
+           sumxy1+=dx*clusters1[dlayer]->GetY();       
+           sumy1 +=clusters1[dlayer]->GetY();
+         }
+         if (index2>=0) {
+           clusters2[dlayer]= (AliTRDcluster*)GetCluster(index2);
+           ncl++;
+           sum2++;
+           Double_t dx = layer2C.GetX()-x2;
+           sumx2 +=dx;
+           sum2x2+=dx*dx;
+           sumxy2+=dx*clusters2[dlayer]->GetY();
+           sumy2 +=clusters2[dlayer]->GetY();
+         }
+       }
+       if (sum1<10) continue;
+       if (sum2<10) continue;
+       //
+       Double_t det1   = sum1*sum2x1-sumx1*sumx1;
+       Double_t angle1 = (sum1*sumxy1-sumx1*sumy1)/det1;
+       Double_t pos1   = (sum2x1*sumy1-sumx1*sumxy1)/det1;  // at x1 
+       //
+       Double_t det2   = sum2*sum2x2-sumx2*sumx2;
+       Double_t angle2 = (sum2*sumxy2-sumx2*sumy2)/det2;
+       Double_t pos2   = (sum2x2*sumy2-sumx2*sumxy2)/det2;  // at x2
+       //
+       //
+
+       Double_t sumM=0, sumxM=0,sum2xM=0,sumxyM=0, sumyM=0;
+       //
+       for (Int_t dlayer=1;dlayer<15;dlayer++){
+         clustersM[dlayer]=0;
+         AliTRDpropagationLayer& layerM=*(fTrSec[ns]->GetLayer(ilayerM-dlayer));  //select propagation layers
+         Double_t yyM = y1+(tanphi)  *(layerM.GetX()-x1);
+         Double_t zzM = z1+(tantheta)*(layerM.GetX()-x1);
+         Int_t indexM = layerM.FindNearestCluster(yyM,zzM,road);
+         if (indexM>=0) {
+           clustersM[dlayer]= (AliTRDcluster*)GetCluster(indexM);
+           ncl++;
+           sumM++;
+           Double_t dx = layerM.GetX()-xm;
+           sumxM +=dx;
+           sum2xM+=dx*dx;
+           sumxyM+=dx*clustersM[dlayer]->GetY();       
+           sumyM +=clustersM[dlayer]->GetY();
+         }
+       }
+       Double_t detM   = sumM*sum2xM-sumxM*sumxM;
+       Double_t posM=0, angleM=0;
+       if (TMath::Abs(detM)>0.0000001){
+         angleM = (sumM*sumxyM-sumxM*sumyM)/detM;
+         posM   = (sum2xM*sumyM-sumxM*sumxyM)/detM;  // at xm
+       }
+       //
+
+       if (ncl>15){
+         TTreeSRedirector& cstream = *fDebugStreamer;
+         cstream<<"Seeds"<<
+           "Ncl="<<ncl<<
+           "SumM="<<sumM<<
+           "x1="<<x1<<
+           "x2="<<x2<<
+           "Cl1.="<<cl1<<
+           "Cl2.="<<cl2<<
+           "Phi="<<tanphi<<
+           "Theta="<<tantheta<<
+           "Pos1="<<pos1<<
+           "Pos2="<<pos2<<
+           "PosM="<<posM<<
+           "Angle1="<<angle1<<
+           "Angle2="<<angle2<<
+           "AngleM="<<angleM<<
+           "\n";
+       }
+      }
+    }
+  }
+}            
 
 //_____________________________________________________________________________
 Int_t AliTRDtracker::ReadClusters(TObjArray *array, TTree *ClusterTree) const
@@ -1737,10 +1791,10 @@ Int_t AliTRDtracker::ReadClusters(TObjArray *array, TTree *ClusterTree) const
     // Loop through all TRD digits
     for (Int_t iCluster = 0; iCluster < nCluster; iCluster++) { 
       c = (AliTRDcluster*)clusterArray->UncheckedAt(iCluster);
-      if (c->GetNPads()>3&&(iCluster%3>0)) {
-       delete clusterArray->RemoveAt(iCluster);
-       continue;
-      }
+//       if (c->GetNPads()>3&&(iCluster%3>0)) {
+//     delete clusterArray->RemoveAt(iCluster);
+//     continue;
+//       }
       //      AliTRDcluster *co = new AliTRDcluster(*c);  //remove unnecesary coping - + clusters are together in memory
       AliTRDcluster *co = c;
       co->SetSigmaY2(c->GetSigmaY2() * fSY2corr);
@@ -1939,7 +1993,6 @@ AliTRDtracker::AliTRDtrackingSector::AliTRDtrackingSector(AliTRDgeometry* geo, I
   //
   // AliTRDtrackingSector Constructor
   //
-
   AliTRDpadPlane *padPlane = 0;
 
   fGeom = geo;
@@ -2097,11 +2150,12 @@ AliTRDtracker::AliTRDtrackingSector::AliTRDtrackingSector(AliTRDgeometry* geo, I
     }
 
     ymax          = fGeom->GetChamberWidth(plane)/2.;
-    //
     // Modidified for new pad plane class, 22.04.05 (C.B.)
     // ymaxsensitive = (fPar->GetColPadSize(plane)*fPar->GetColMax(plane)-4)/2.;
     padPlane = fPar->GetPadPlane(plane,0);
     ymaxsensitive = (padPlane->GetColSize(1)*padPlane->GetNcols()-4)/2.;
+
+    //    ymaxsensitive = (fPar->GetColPadSize(plane)*fPar->GetColMax(plane)-4)/2.;
     
     for(Int_t ch = 0; ch < kNchambers; ch++) {
       zmax[ch] = fGeom->GetChamberLength(plane,ch)/2;
@@ -2109,16 +2163,18 @@ AliTRDtracker::AliTRDtrackingSector::AliTRDtrackingSector(AliTRDgeometry* geo, I
       // Modidified for new pad plane class, 22.04.05 (C.B.)
       //Float_t pad = fPar->GetRowPadSize(plane,ch,0);
       Float_t pad = padPlane->GetRowSize(1);
+      //Float_t pad = fPar->GetRowPadSize(plane,ch,0);
       Float_t row0 = fPar->GetRow0(plane,ch,0);
       Int_t nPads = fPar->GetRowMax(plane,ch,0);
       zmaxsensitive[ch] = Float_t(nPads)*pad/2.;      
       //      zc[ch] = (pad * nPads)/2 + row0 - pad/2;
-      zc[ch] = (pad * nPads)/2 + row0;
+      //      zc[ch] = (pad * nPads)/2 + row0;
+      zc[ch] = -(pad * nPads)/2 + row0;
       //zc[ch] = row0+zmax[ch]-AliTRDgeometry::RpadW();
 
     }
 
-    dx  = fPar->GetDriftVelocity()
+    dx  = fgkDriftCorrection*fPar->GetDriftVelocity()
         / fPar->GetSamplingFrequency();
     rho = 0.00295 * 0.85; radLength = 11.0;  
 
@@ -2130,7 +2186,7 @@ AliTRDtracker::AliTRDtrackingSector::AliTRDtrackingSector(AliTRDgeometry* geo, I
     steps = (Int_t) (dxAmp/dx);
 
     for(tb = 0; tb < steps; tb++) {
-      x = x0 + tb * dx + dx/2;
+      x = x0 + tb * dx + dx/2+ fgkOffsetX;
       tbIndex = CookTimeBinIndex(plane, -tb-1);
       ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex);
       ppl->SetYmax(ymax,ymaxsensitive);
@@ -2148,12 +2204,13 @@ AliTRDtracker::AliTRDtrackingSector::AliTRDtrackingSector(AliTRDgeometry* geo, I
     InsertLayer(ppl);
 
     // Drift region
-    dx = fPar->GetDriftVelocity()
+
+    dx = fgkDriftCorrection*fPar->GetDriftVelocity()
        / fPar->GetSamplingFrequency();
     steps = (Int_t) (dxDrift/dx);
 
     for(tb = 0; tb < steps; tb++) {
-      x = x0 - tb * dx - dx/2;
+      x = x0 - tb * dx - dx/2 + fgkOffsetX;                      //temporary fix - fix it the parameters
       tbIndex = CookTimeBinIndex(plane, tb);
 
       ppl = new AliTRDpropagationLayer(x,dx,rho,radLength,tbIndex);
@@ -2229,7 +2286,8 @@ Int_t  AliTRDtracker::AliTRDtrackingSector::CookTimeBinIndex(Int_t plane, Int_t
 
   Double_t dxAmp = (Double_t) fGeom->CamHght();   // Amplification region
   Double_t dxDrift = (Double_t) fGeom->CdrHght(); // Drift region  
-  Double_t dx = (Double_t) fPar->GetDriftVelocity()
+  
+  Double_t dx = fgkDriftCorrection*(Double_t) fPar->GetDriftVelocity()
                          / fPar->GetSamplingFrequency();
 
   Int_t tbAmp = fPar->GetTimeBefore();
@@ -2397,6 +2455,10 @@ Int_t AliTRDtracker::AliTRDtrackingSector::Find(Double_t x) const
   return m;
 }             
 
+
+
+
+
 //______________________________________________________
 void AliTRDtracker::AliTRDpropagationLayer::SetZ(Double_t* center, Double_t *w, Double_t *wsensitive )
 {
@@ -2424,7 +2486,7 @@ void AliTRDtracker::AliTRDpropagationLayer::SetHoles(Bool_t *holes)
 
 
 
-void AliTRDtracker::AliTRDpropagationLayer::GetPropagationParameters(
+Bool_t AliTRDtracker::AliTRDpropagationLayer::GetPropagationParameters(
         Double_t y, Double_t z, Double_t &dx, Double_t &rho, Double_t &radLength, 
         Bool_t &lookForCluster) const
 {
@@ -2467,15 +2529,15 @@ void AliTRDtracker::AliTRDpropagationLayer::GetPropagationParameters(
        radLength = 36.66;
       }
     }else{
-      rho = 2.7; radLength = 24.01;  //aluminium in between
+      cross = kTRUE; rho = 2.7; radLength = 24.01;  //aluminium in between
     }        
   }
   //
-  if (fTimeBinIndex>=0) return;
+  if (fTimeBinIndex>=0) return cross;
   //
   //
   // check hole
-  if (fHole==kFALSE) return;
+  if (fHole==kFALSE) return cross;
   //
   for(Int_t ch = 0; ch < (Int_t) kZones; ch++) {
     if  (TMath::Abs(z - fZc[ch]) < fZmax[ch]){ 
@@ -2486,7 +2548,7 @@ void AliTRDtracker::AliTRDpropagationLayer::GetPropagationParameters(
       }    
     }
   }
-  return;
+  return cross;
 }
 
 Int_t  AliTRDtracker::AliTRDpropagationLayer::GetZone( Double_t z) const
@@ -2543,15 +2605,42 @@ Int_t AliTRDtracker::AliTRDpropagationLayer::Find(Double_t y) const {
   return m;
 }    
 
+Int_t AliTRDtracker::AliTRDpropagationLayer::FindNearestCluster(Double_t y, Double_t z, Double_t maxroad) const 
+{
+  //
+  // Returns index of the cluster nearest to the given y,z
+  //
+  Int_t index = -1;
+  Int_t maxn = fN;
+  Double_t mindist = maxroad;                  
+  Float_t padlength =-1;
+  //
+  for (Int_t i=Find(y-maxroad); i<maxn; i++) {
+    AliTRDcluster* c=(AliTRDcluster*)(fClusters[i]);
+    if (padlength<0){
+      padlength = TMath::Sqrt(c->GetSigmaZ2()*12); 
+    }
+    //
+    if (c->GetY() > y+maxroad) break;
+    if((c->GetZ()-z)*(c->GetZ()-z) > padlength*0.75) continue;      
+    if (TMath::Abs(c->GetY()-y)<mindist){
+      mindist = TMath::Abs(c->GetY()-y);
+      index = GetIndex(i);
+    }        
+  }                                            
+  return index;
+}             
+
+
 //---------------------------------------------------------
 
 Double_t AliTRDtracker::GetTiltFactor(const AliTRDcluster* c) {
 //
 //  Returns correction factor for tilted pads geometry 
 //
-
   AliTRDpadPlane *padPlane = fPar->GetPadPlane(0,0);
   Double_t h01 = sin(TMath::Pi() / 180.0 * padPlane->GetTiltingAngle());
+  //Double_t h01 = sin(TMath::Pi() / 180.0 * fPar->GetTiltingAngle());
   Int_t det = c->GetDetector();    
   Int_t plane = fGeom->GetPlane(det);
 
@@ -2630,169 +2719,448 @@ void AliTRDtracker::CookdEdxTimBin(AliTRDtrack& TRDtrack)
 } // end of function
 
 
-Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1, AliTRDtrack * track, Int_t *clusters)
+Int_t AliTRDtracker::FindClusters(Int_t sector, Int_t t0, Int_t t1, AliTRDtrack * track, Int_t *clusters,AliTRDtracklet&tracklet)
 {
   //
   //
   //  try to find nearest clusters to the track in timebins from t0 to t1 
   //  
   //
-  Double_t x[1000],yt[1000],zt[10000];
-  Double_t dz[1000],dy[10000];
-  Int_t    indexes[2][10000];
-  AliTRDcluster *cl[2][10000];
-  
-  for (Int_t it=t0;it<t1; it++){
-    clusters[it]=-2;
-    indexes[0][it]=-2;              //reset indexes1
-    indexes[1][it]=-2;              //reset indexes1
+  //  
+  // correction coeficients   - depends on TRD parameters  - to be changed according it
+  //
+
+  Double_t x[100],yt[100],zt[100];
+  Double_t xmean=0;   //reference x
+  Double_t dz[10][100],dy[10][100];
+  Float_t zmean[100], nmean[100];
+  Int_t    clfound=0;
+  Int_t    indexes[10][100];    // indexes of the clusters in the road
+  AliTRDcluster *cl[10][100];   // pointers to the clusters in the road
+  Int_t    best[10][100];       // index of best matching cluster 
+  //
+  //
+  TClonesArray array0("AliTRDcluster",1);
+  TClonesArray array1("AliTRDcluster",1);
+  for (Int_t it=0;it<t1-t0; it++){
     x[it]=0;
     yt[it]=0;
     zt[it]=0;
-    dz[it]=0;
-    dy[it]=0;
-    cl[0][it]=0;
-    cl[1][it]=0;
+    clusters[it+t0]=-2;
+    zmean[it]=0;
+    nmean[it]=0;
+    //
+    for (Int_t ih=0;ih<10;ih++){
+      indexes[ih][it]=-2;              //reset indexes1
+      cl[ih][it]=0;
+      dz[ih][it]=-100;
+      dy[ih][it]=-100;
+      best[ih][it]=0;
+    }
   }  
   //
   Double_t x0 = track->GetX();
-  Double_t sy2=ExpectedSigmaY2(x0,track->GetTgl(),track->GetPt());
-  Double_t sz2=ExpectedSigmaZ2(x0,track->GetTgl());
-  Double_t road = 10.*sqrt(track->GetSigmaY2() + sy2);
+  Double_t sigmaz = TMath::Sqrt(track->GetSigmaZ2());
   Int_t nall=0;
   Int_t nfound=0;
+  Double_t h01 =0;
+  Int_t plane =-1;
+  Float_t padlength=0;
+  AliTRDtrack track2(*track);
+  Float_t snpy = track->GetSnp();
+  Float_t tany = TMath::Sqrt(snpy*snpy/(1.-snpy*snpy)); 
+  if (snpy<0) tany*=-1;
+  //
+  Double_t sy2=ExpectedSigmaY2(x0,track->GetTgl(),track->GetPt());
+  Double_t sz2=ExpectedSigmaZ2(x0,track->GetTgl());
+  Double_t road = 15.*sqrt(track->GetSigmaY2() + sy2);
+  if (road>6.) road=6.;
 
-  for (Int_t it=t0;it<t1;it++){
-    Double_t maxChi2=fgkMaxChi2;      
-    AliTRDpropagationLayer& timeBin=*(fTrSec[sector]->GetLayer(it));
+  //
+  for (Int_t it=0;it<t1-t0;it++){
+    Double_t maxChi2[2]={fgkMaxChi2,fgkMaxChi2};      
+    AliTRDpropagationLayer& timeBin=*(fTrSec[sector]->GetLayer(it+t0));
     if (timeBin==0) continue;  // no indexes1
     Int_t maxn = timeBin;
     x[it] = timeBin.GetX();
-    Double_t  y,z;
-    if (!track->GetProlongation(x[it],y,z)) continue;  
-    yt[it]=y;
-    zt[it]=z;
+    track2.PropagateTo(x[it]);
+    yt[it] = track2.GetY();
+    zt[it] = track2.GetZ();
+    
+    Double_t  y=yt[it],z=zt[it];
     Double_t chi2 =1000000;
     nall++;
     //
-    // find nearest cluster at given pad
+    // find 2 nearest cluster at given time bin
+    // 
+    // 
     for (Int_t i=timeBin.Find(y-road); i<maxn; i++) {
       AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);
-      Double_t h01 = GetTiltFactor(c);
-      if (c->GetY() > y+road) break;
-      if (c->IsUsed() > 0) continue;
-      if((c->GetZ()-z)*(c->GetZ()-z) > 3 * sz2) continue;      
-      chi2=track->GetPredictedChi2(c,h01);
-      if (chi2 > maxChi2) continue;
-      maxChi2=chi2;
-      cl[0][it]=c;
-      indexes[0][it] =timeBin.GetIndex(i);         
-    }         
-    //
-    // find nearest cluster also in adjacent 2 pads
-    //
-    for (Int_t i=timeBin.Find(y-road); i<maxn; i++) {
-      AliTRDcluster* c=(AliTRDcluster*)(timeBin[i]);   
-      Double_t h01 = GetTiltFactor(c);
+      h01 = GetTiltFactor(c);
+      if (plane<0){
+       Int_t det = c->GetDetector();    
+       plane = fGeom->GetPlane(det);
+       padlength = TMath::Sqrt(c->GetSigmaZ2()*12.);
+      }
+      //      if (c->GetLocalTimeBin()==0) continue;
       if (c->GetY() > y+road) break;
-      if (c->IsUsed() > 0) continue;
-      if((c->GetZ()-z)*(c->GetZ()-z) >  12 * sz2) continue;
-      chi2=track->GetPredictedChi2(c,h01);     
-      if (chi2 > maxChi2) continue;
-      maxChi2=chi2;
-      cl[1][it]=c;
-      indexes[1][it]=timeBin.GetIndex(i);
-      if (!cl[0][it]) {
-       cl[0][it]=c;
-       indexes[0][it]=timeBin.GetIndex(i);
+      if((c->GetZ()-z)*(c->GetZ()-z) > 12. * sz2) continue;      
+
+      Double_t dist = TMath::Abs(c->GetZ()-z);
+      if (dist> (0.5*padlength+6.*sigmaz)) continue;   // 6 sigma boundary cut
+      Double_t cost = 0;
+      //
+      if (dist> (0.5*padlength-sigmaz)){   //  sigma boundary cost function
+       cost =  (dist-0.5*padlength)/(2.*sigmaz);
+       if (cost>-1) cost= (cost+1.)*(cost+1.);
+       else cost=0;
+      }      
+      //      Int_t label = TMath::Abs(track->GetLabel());
+      //      if (c->GetLabel(0)!=label && c->GetLabel(1)!=label&&c->GetLabel(2)!=label) continue;
+      chi2=track2.GetPredictedChi2(c,h01)+cost;
+      //
+      clfound++;      
+      if (chi2 > maxChi2[1]) continue;
+      
+      for (Int_t ih=2;ih<9; ih++){  //store the clusters in the road
+       if (cl[ih][it]==0){
+         cl[ih][it] = c;
+         indexes[ih][it] =timeBin.GetIndex(i);   // index - 9 - reserved for outliers
+         break;
+       }
       }
-    }    
-    if (cl[0][it]||cl[1][it]) nfound++;
+      //
+      if (chi2 <maxChi2[0]){
+       maxChi2[1]     = maxChi2[0];
+       maxChi2[0]     = chi2;
+       indexes[1][it] = indexes[0][it];
+       cl[1][it]      = cl[0][it];
+       indexes[0][it] = timeBin.GetIndex(i);
+       cl[0][it]      = c;
+       continue;
+      }
+      maxChi2[1]=chi2;
+      cl[1][it] = c;
+      indexes[1][it] =timeBin.GetIndex(i); 
+    }         
+    if (cl[0][it]){
+      nfound++;
+      xmean += x[it];
+    }
   }
   //
-  if (nfound<5) return 0;  
+  if (nfound<4) return 0;  
+  xmean /=Float_t(nfound);     // middle x
+  track2.PropagateTo(xmean);   // propagate track to the center
   //
   // choose one of the variants
   //
-  Int_t changes[2]={0,0};
-  Float_t sigmay[2]={1000,1000};
-  Float_t meany[2] ={1000,1000};
-  Float_t meanz[2] ={1000,1000};
-  Int_t sumall[2]  ={0,0};
-  Int_t ngood[2] ={0,0};
-  Int_t nbad[2]  ={0,0};
-  //
-  //
-  for (Int_t ih=0; ih<2;ih++){    
-    Float_t lastz    =-10000;
-    Float_t sumz   =0;
-    Float_t sum    =0;
-    Double_t sumdy = 0;
-    Double_t sumdy2= 0;
-    //
-    // how many changes ++ mean z ++mean y ++ sigma y
-    for (Int_t it=t0;it<t1;it++){
+  Int_t changes[10];
+  Float_t sumz      = 0;
+  Float_t sum       = 0;
+  Double_t sumdy    = 0;
+  Double_t sumdy2   = 0;
+  Double_t sumx     = 0;
+  Double_t sumxy    = 0;
+  Double_t sumx2    = 0;
+  Double_t mpads    = 0;
+  //
+  Int_t   ngood[10];
+  Int_t   nbad[10];
+  //
+  Double_t meanz[10];
+  Double_t moffset[10];    // mean offset
+  Double_t mean[10];       // mean value
+  Double_t angle[10];      // angle
+  //
+  Double_t smoffset[10];   // sigma of mean offset
+  Double_t smean[10];      // sigma of mean value
+  Double_t sangle[10];     // sigma of angle
+  Double_t smeanangle[10]; // correlation
+  //
+  Double_t sigmas[10];     
+  Double_t tchi2s[10];      // chi2s for tracklet
+  //
+  // calculate zmean
+  //
+  for (Int_t it=0;it<t1-t0;it++){
+    if (!cl[0][it]) continue;
+    for (Int_t dt=-3;dt<=3;dt++){
+      if (it+dt<0) continue;
+      if (it+dt>t1) continue;
+      if (!cl[0][it+dt]) continue;
+      zmean[it]+=cl[0][it+dt]->GetZ();
+      nmean[it]+=1.;
+    }
+    zmean[it]/=nmean[it]; 
+  }
+  //
+  for (Int_t it=0; it<t1-t0;it++){
+    best[0][it]=0;
+    for (Int_t ih=0;ih<10;ih++){
+      dz[ih][it]=-100;
+      dy[ih][it]=-100;
       if (!cl[ih][it]) continue;
-      sumall[ih]++;
-      if (lastz<-9999) lastz = cl[ih][it]->GetZ();
-      if (TMath::Abs(lastz-cl[ih][it]->GetZ())>1) {
-       lastz = cl[ih][it]->GetZ();
-       changes[ih]++;
+      Float_t poscor =  fgkCoef*(cl[ih][it]->GetLocalTimeBin() - fgkMean)+fgkOffset;      
+      dz[ih][it]  = cl[ih][it]->GetZ()- zt[it];                               // calculate distance from track  in z
+      dy[ih][it]  = cl[ih][it]->GetY()+ dz[ih][it]*h01 - poscor  -yt[it];     //                                in y
+    }
+    // minimize changes
+    if (!cl[0][it]) continue;
+    if (TMath::Abs(cl[0][it]->GetZ()-zmean[it])> padlength*0.8 &&cl[1][it])
+      if (TMath::Abs(cl[1][it]->GetZ()-zmean[it])< padlength*0.5){
+       best[0][it]=1;
       }
-      sumz = cl[ih][it]->GetZ();
-      sum++;
-      Double_t h01 = GetTiltFactor(cl[ih][it]);
-      dz[it]  = cl[ih][it]->GetZ()- zt[it]; 
-      dy[it]  = cl[ih][it]->GetY()+ dz[it]*h01 -yt[it];
+  }
+  //
+  // iterative choosing of "best path"
+  //
+  //
+  Int_t label = TMath::Abs(track->GetLabel());
+  Int_t bestiter=0;
+  //
+  for (Int_t iter=0;iter<9;iter++){
+    //
+    changes[iter]= 0;
+    sumz      = 0; sum=0; sumdy=0;sumdy2=0;sumx=0;sumx2=0;sumxy=0;mpads=0; ngood[iter]=0; nbad[iter]=0; 
+    // linear fit
+    for (Int_t it=0;it<t1-t0;it++){
+      if (!cl[best[iter][it]][it]) continue;
+      //calculates pad-row changes
+      Double_t zbefore= cl[best[iter][it]][it]->GetZ();
+      Double_t zafter = cl[best[iter][it]][it]->GetZ();
+      for (Int_t itd = it-1; itd>=0;itd--) {
+       if (cl[best[iter][itd]][itd]) {
+         zbefore= cl[best[iter][itd]][itd]->GetZ();
+         break;
+       }
+      }
+      for (Int_t itd = it+1; itd<t1-t0;itd++) {
+       if (cl[best[iter][itd]][itd]) {
+         zafter= cl[best[iter][itd]][itd]->GetZ();
+         break;
+       }
+      }
+      if (TMath::Abs(cl[best[iter][it]][it]->GetZ()-zbefore)>0.1&&TMath::Abs(cl[best[iter][it]][it]->GetZ()-zafter)>0.1) changes[iter]++;
+      //
+      Double_t dx = x[it]-xmean;  // distance to reference x
+      sumz += cl[best[iter][it]][it]->GetZ();      
       sum++;
-      sumdy += dy[it];
-      sumdy2+= dy[it]*dy[it];
-      Int_t label = TMath::Abs(track->GetLabel());
-      if (cl[ih][it]->GetLabel(0)==label || cl[ih][it]->GetLabel(1)==label||cl[ih][it]->GetLabel(2)==label){
-       ngood[ih]++;
+      sumdy += dy[best[iter][it]][it];
+      sumdy2+= dy[best[iter][it]][it]*dy[best[iter][it]][it];
+      sumx  += dx;
+      sumx2 += dx*dx;
+      sumxy  += dx*dy[best[iter][it]][it];
+      mpads += cl[best[iter][it]][it]->GetNPads();
+      if (cl[best[iter][it]][it]->GetLabel(0)==label || cl[best[iter][it]][it]->GetLabel(1)==label||cl[best[iter][it]][it]->GetLabel(2)==label){
+       ngood[iter]++;
       }
       else{
-       nbad[ih]++;
+       nbad[iter]++;
       }
     }
-    if (sumall[ih]>4){
-      meanz[ih]  = sumz/sum;    
-      meany[ih]  = sumdy/sum;
-      sigmay[ih] = TMath::Sqrt(sumdy2/sum-meany[ih]*meany[ih]);
+    //
+    // calculates line parameters
+    //
+    Double_t det  = sum*sumx2-sumx*sumx;
+    angle[iter]   = (sum*sumxy-sumx*sumdy)/det;
+    mean[iter]    = (sumx2*sumdy-sumx*sumxy)/det;
+    meanz[iter]   = sumz/sum;    
+    moffset[iter] = sumdy/sum;
+    mpads        /= sum;                         // mean number of pads
+    //
+    //
+    Double_t  sigma2 = 0;   // normalized residuals - for line fit
+    Double_t  sigma1 = 0;   // normalized residuals - constant fit
+    //
+    for (Int_t it=0;it<t1-t0;it++){
+      if (!cl[best[iter][it]][it]) continue;
+      Double_t dx = x[it]-xmean;
+      Double_t ytr = mean[iter]+angle[iter]*dx;
+      sigma2 += (dy[best[iter][it]][it]-ytr)*(dy[best[iter][it]][it]-ytr);
+      sigma1 +=  (dy[best[iter][it]][it]-moffset[iter])*(dy[best[iter][it]][it]-moffset[iter]);
+      sum++;
     }
-  }
-  Int_t best =0;
-  /*
-    if (sumall[0]<sumall[1]-2&&sigmay[1]<0.1){
-    if (sigmay[1]<sigmay[0]) best = 1;               // if sigma is better
+    sigma2      /=(sum-2);                    // normalized residuals
+    sigma1      /=(sum-1);                    // normalized residuals
+    //
+    smean[iter]       = sigma2*(sumx2/det);   // estimated error2 of mean
+    sangle[iter]      = sigma2*(sum/det);     // estimated error2 of angle
+    smeanangle[iter]  = sigma2*(-sumx/det);   // correlation
+    //
+    //
+    sigmas[iter]  = TMath::Sqrt(sigma1);      //
+    smoffset[iter]= (sigma1/sum)+0.01*0.01;             // sigma of mean offset + unisochronity sigma 
+    //
+    // iterative choosing of "better path"
+    //
+    for (Int_t it=0;it<t1-t0;it++){
+      if (!cl[best[iter][it]][it]) continue;
+      //
+      Double_t sigmatr2 = smoffset[iter]+0.5*tany*tany;             //add unisochronity + angular effect contribution
+      Double_t sweight  = 1./sigmatr2+1./track->fCyy;
+      Double_t weighty  = (moffset[iter]/sigmatr2)/sweight;         // weighted mean
+      Double_t sigmacl  = TMath::Sqrt(sigma1*sigma1+track->fCyy);   //
+      Double_t mindist=100000; 
+      Int_t ihbest=0;
+      for (Int_t ih=0;ih<10;ih++){
+       if (!cl[ih][it]) break;
+       Double_t dist2 = (dy[ih][it]-weighty)/sigmacl;
+       dist2*=dist2;    //chi2 distance
+       if (dist2<mindist){
+         mindist = dist2;
+         ihbest =ih;
+       }
+      }
+      best[iter+1][it]=ihbest;
     }
-  */
-  Float_t quality0 = (sigmay[0]+TMath::Abs(meany[0]))*(1.+Float_t(changes[0])/Float_t(sumall[0]));
-  Float_t quality1 = (sigmay[1]+TMath::Abs(meany[1]))*(1.+Float_t(changes[1])/Float_t(sumall[1]));
-
-  if (quality0>quality1){
-    best = 1;
-  }
-
-
-  //
-  for (Int_t it=t0;it<t1;it++){
-    if (!cl[best][it]) continue;
-    Double_t h01 = GetTiltFactor(cl[best][it]);
-    dz[it]  = cl[best][it]->GetZ()- zt[it]; 
-    dy[it]  = cl[best][it]->GetY()+ dz[it]*h01 -yt[it];
     //
-    if (TMath::Abs(dy[it])<2.5*sigmay[best])
-      clusters[it] = indexes[best][it];    
-  }
+    //  update best hypothesy if better chi2 according tracklet position and angle
+    //
+    Double_t sy2 = smean[iter]  + track->fCyy;
+    Double_t sa2 = sangle[iter] + track->fCee;
+    Double_t say = track->fCey;
+    //    Double_t chi20 = mean[bestiter]*mean[bestiter]/sy2+angle[bestiter]*angle[bestiter]/sa2;
+    // Double_t chi21 = mean[iter]*mean[iter]/sy2+angle[iter]*angle[iter]/sa2;
+
+    Double_t detchi    = sy2*sa2-say*say;
+    Double_t invers[3] = {sa2/detchi, sy2/detchi, -say/detchi};   //inverse value of covariance matrix  
     
-  if (nbad[0]>4){
-    nbad[0] = nfound;
+    Double_t chi20 = mean[bestiter]*mean[bestiter]*invers[0]+angle[bestiter]*angle[bestiter]*invers[1]+
+      2.*mean[bestiter]*angle[bestiter]*invers[2];
+    Double_t chi21 = mean[iter]*mean[iter]*invers[0]+angle[iter]*angle[iter]*invers[1]+
+      2*mean[iter]*angle[iter]*invers[2];
+    tchi2s[iter] =chi21;
+    //
+    if (changes[iter]<=changes[bestiter] && chi21<chi20) {
+      bestiter =iter;      
+    }
+  }
+  //
+  //set clusters 
+  //
+  Double_t sigma2 = sigmas[0];   // choose as sigma  from 0 iteration
+  //if (tchi2s[bestiter]>25.) sigma2*=tchi2s[bestiter]/25.;
+  //if (tchi2s[bestiter]>25.) sigma2=1000.;  // dont'accept
+
+  Double_t expectederr = sigma2*sigma2+0.01*0.01;
+  if (mpads>3.5) expectederr  +=   (mpads-3.5)*0.04;
+  if (changes[bestiter]>1) expectederr+=   changes[bestiter]*0.01; 
+  expectederr+=(0.03*(tany-fgkExB)*(tany-fgkExB))*15;
+  //  if (tchi2s[bestiter]>18.) expectederr*= tchi2s[bestiter]/18.;
+  //expectederr+=10000;
+  for (Int_t it=0;it<t1-t0;it++){
+    if (!cl[best[bestiter][it]][it]) continue;
+    Float_t poscor =  fgkCoef*(cl[best[bestiter][it]][it]->GetLocalTimeBin() - fgkMean)+fgkOffset;          
+    cl[best[bestiter][it]][it]->SetSigmaY2(expectederr);  // set cluster error
+    if (!cl[best[bestiter][it]][it]->IsUsed()){
+      cl[best[bestiter][it]][it]->SetY( cl[best[bestiter][it]][it]->GetY()-poscor);  // ExB corrction correction    
+      cl[best[bestiter][it]][it]->Use();
+    }
+    clusters[it+t0] = indexes[best[bestiter][it]][it];    
+  } 
+  //
+  // set tracklet parameters
+  //
+  Double_t trackleterr2 = smoffset[bestiter]+0.01*0.01;
+  if (mpads>3.5) trackleterr2  +=   (mpads-3.5)*0.04;
+  trackleterr2+=   changes[bestiter]*0.01;
+  trackleterr2*=   TMath::Max(14.-nfound,1.);
+  trackleterr2+=   0.2*(tany-fgkExB)*(tany-fgkExB); 
+  //
+  tracklet.Set(xmean, track2.GetY()+moffset[bestiter], meanz[bestiter], track2.GetAlpha(), trackleterr2);  //set tracklet parameters
+  tracklet.SetTilt(h01);
+  tracklet.SetP0(mean[bestiter]);
+  tracklet.SetP1(angle[bestiter]);
+  tracklet.SetN(nfound);
+  tracklet.SetNCross(changes[bestiter]);
+  tracklet.SetPlane(plane);
+  tracklet.SetSigma2(expectederr);
+  tracklet.SetChi2(tchi2s[bestiter]);
+  track->fTracklets[plane] = tracklet;
+  track->fNWrong+=nbad[0];
+  //
+  // Debuging part
+  //
+  TTreeSRedirector& cstream = *fDebugStreamer;
+  AliTRDcluster dummy;
+  Double_t dy0[100];
+  Double_t dyb[100];  
+  for (Int_t it=0;it<t1-t0;it++){
+    dy0[it] = dy[0][it];
+    dyb[it] = dy[best[bestiter][it]][it];
+    if(cl[0][it]) {
+      new(array0[it]) AliTRDcluster(*cl[0][it]);
+    }
+    else{
+      new(array0[it]) AliTRDcluster(dummy);
+    }
+    if(cl[best[bestiter][it]][it]) {
+      new(array1[it]) AliTRDcluster(*cl[best[bestiter][it]][it]);
+    }
+    else{
+      new(array1[it]) AliTRDcluster(dummy);
+    }
   }
+  TGraph graph0(t1-t0,x,dy0);
+  TGraph graph1(t1-t0,x,dyb);
+  TGraph graphy(t1-t0,x,yt);
+  TGraph graphz(t1-t0,x,zt);
+  //
+  //
+  cstream<<"tracklet"<<
+    "track.="<<track<<                                       // track parameters
+    "tany="<<tany<<                                          // tangent of the local track angle 
+    "xmean="<<xmean<<                                        // xmean - reference x of tracklet  
+    "tilt="<<h01<<                                           // tilt angle
+    "nall="<<nall<<                                          // number of foundable clusters 
+    "nfound="<<nfound<<                                      // number of found clusters
+    "clfound="<<clfound<<                                    // total number of found clusters in road 
+    "mpads="<<mpads<<                                        // mean number of pads per cluster
+    "plane="<<plane<<                                        // plane number 
+    "road="<<road<<                                          // the width of the used road
+    "graph0.="<<&graph0<<                                    // x - y = dy for closest cluster
+    "graph1.="<<&graph1<<                                    // x - y = dy for second closest cluster    
+    "graphy.="<<&graphy<<                                    // y position of the track
+    "graphz.="<<&graphz<<                                    // z position of the track
+    "fCl.="<<&array0<<                                       // closest cluster
+    "fCl2.="<<&array1<<                                      // second closest cluster
+    //
+    "bestiter="<<bestiter<<                                  // best iteration number 
+    "tracklet.="<<&tracklet<<                                // corrspond to the best iteration
+    "tchi20="<<tchi2s[0]<<                                   // chi2 of cluster in the 0 iteration
+    "tchi2b="<<tchi2s[bestiter]<<                            // chi2 of cluster in the best  iteration
+    "sigmas0="<<sigmas[0]<<                                  // residuals sigma 
+    "sigmasb="<<sigmas[bestiter]<<                           // residulas sigma
+    //
+    "ngood0="<<ngood[0]<<                                    // number of good clusters in 0 iteration
+    "nbad0="<<nbad[0]<<                                      // number of bad clusters in 0 iteration
+    "ngoodb="<<ngood[bestiter]<<                             //                        in  best iteration    
+    "nbadb="<<nbad[bestiter]<<                               //                        in  best iteration
+    //
+    "changes0="<<changes[0]<<                                // changes of pardrows in iteration number 0 
+    "changesb="<<changes[bestiter]<<                         // changes of pardrows in best iteration
+    //
+    "moffset0="<<moffset[0]<<                                // offset fixing angle in iter=0
+    "smoffset0="<<smoffset[0]<<                              // sigma of offset fixing angle in iter=0
+    "moffsetb="<<moffset[bestiter]<<                         // offset fixing angle in iter=best
+    "smoffsetb="<<smoffset[bestiter]<<                       // sigma of offset fixing angle in iter=best
+    //
+    "mean0="<<mean[0]<<                                      // mean dy in iter=0;
+    "smean0="<<smean[0]<<                                    // sigma of mean dy in iter=0
+    "meanb="<<mean[bestiter]<<                               // mean dy in iter=best
+    "smeanb="<<smean[bestiter]<<                             // sigma of mean dy in iter=best
+    //
+    "angle0="<<angle[0]<<                                    // angle deviation in the iteration number 0 
+    "sangle0="<<sangle[0]<<                                  // sigma of angular deviation in iteration number 0
+    "angleb="<<angle[bestiter]<<                             // angle deviation in the best iteration   
+    "sangleb="<<sangle[bestiter]<<                           // sigma of angle deviation in the best iteration   
+    //
+    "expectederr="<<expectederr<<                            // expected error of cluster position
+    "\n";
+  //
+  //
   return nfound;
 }
 
 
-
-
-
index fe71d27..8e8a494 100644 (file)
@@ -15,8 +15,10 @@ class TParticlePDG;
 class AliTRDgeometry;
 class AliTRDparameter;
 class AliTRDtrack;
+class AliTRDtracklet;
 class AliTRDcluster;
 class AliESD;
+class TTreeSRedirector;
 
 const unsigned kMaxLayersPerSector = 1000;  
 const unsigned kMaxTimeBinIndex = 216;  // (30 drift + 6 ampl) * 6 planes  
@@ -24,10 +26,14 @@ const unsigned kMaxClusterPerTimeBin = 7000;
 const unsigned kZones = 5; 
 const Int_t    kTrackingSectors = 18; 
 
+
+
+
 class AliTRDtracker : public AliTracker { 
 
  public:
 
+
   AliTRDtracker();
   AliTRDtracker(const TFile *in);
   virtual ~AliTRDtracker(); 
@@ -35,7 +41,7 @@ class AliTRDtracker : public AliTracker {
   Int_t         Clusters2Tracks(AliESD* event);
   Int_t         PropagateBack(AliESD* event);
   Int_t         RefitInward(AliESD* event);
-
+  
   Int_t         LoadClusters(TTree *cTree);
   void          UnloadClusters();
   AliCluster   *GetCluster(Int_t index) const { if (index >= fNclusters) return NULL; 
@@ -50,9 +56,9 @@ class AliTRDtracker : public AliTracker {
 
   Int_t         ReadClusters(TObjArray *array, TTree *in) const;
   Int_t         CookSectorIndex(Int_t gs) const { return kTrackingSectors - 1 - gs; }
-  AliTRDcluster * GetCluster(AliTRDtrack * track, Int_t plane, Int_t timebin);
+  AliTRDcluster * GetCluster(AliTRDtrack * track, Int_t plane, Int_t timebin, UInt_t &index);
   Int_t         GetLastPlane(AliTRDtrack * track); //return last updated plane
-  Int_t FindClusters(Int_t sector, Int_t t0, Int_t t1, AliTRDtrack * track, Int_t *clusters);
+  Int_t FindClusters(Int_t sector, Int_t t0, Int_t t1, AliTRDtrack * track, Int_t *clusters, AliTRDtracklet& tracklet);
 
   Float_t  GetSeedGap()       const {return fgkSeedGap;}   
   Int_t    GetMaxGap()        const {return fMaxGap;}   
@@ -101,11 +107,12 @@ class AliTRDtracker : public AliTracker {
      Double_t       GetRho() const { return fRho; }
      Double_t       GetX0() const { return fX0; }
      Int_t          GetTimeBinIndex() const { return fTimeBinIndex; }     
-     void           GetPropagationParameters(Double_t y, Double_t z,
+     Bool_t         GetPropagationParameters(Double_t y, Double_t z,
                                 Double_t &dx, Double_t &rho, Double_t &x0, 
-                                Bool_t &lookForCluster) const;
+                                            Bool_t &lookForCluster) const;
      Int_t          GetZone( Double_t z) const;
      Int_t          Find(Double_t y) const; 
+     Int_t          FindNearestCluster(Double_t y, Double_t z, Double_t maxroad) const;
      void           SetZmax(Int_t cham, Double_t center, Double_t w)
                       { fZc[cham] = center;  fZmax[cham] = w; }
      void           SetZ(Double_t* center, Double_t *w, Double_t *wsensitive);
@@ -181,6 +188,8 @@ class AliTRDtracker : public AliTracker {
 
  protected:
 
+  friend class AliTRDtracker::AliTRDtrackingSector;
+
   AliTRDgeometry     *fGeom;            // Pointer to TRD geometry
   AliTRDparameter    *fPar;             // Pointer to TRD parameter
 
@@ -231,6 +240,13 @@ class AliTRDtracker : public AliTracker {
   static const Double_t fgkSeedErrorSZ;    // sz parameter in MakeSeeds
   static const Float_t  fgkLabelFraction;  // min fraction of same label
   static const Float_t  fgkWideRoad;       // max road width in FindProlongation
+  //
+  static const Double_t  fgkDriftCorrection; // correction coeficients for drift velocity
+  static const Double_t  fgkOffset;         // correction coeficients 
+  static const Double_t  fgkOffsetX;         // correction coeficients offset in X 
+  static const Double_t  fgkCoef;           // correction coeficients
+  static const Double_t  fgkMean;           // correction coeficients
+  static const Double_t  fgkExB;            // correction coeficients
 
   Bool_t                fVocal;            // Whatever...
   Bool_t                fAddTRDseeds;      // Something else
@@ -249,6 +265,7 @@ class AliTRDtracker : public AliTracker {
  private:
 
   virtual void  MakeSeeds(Int_t inner, Int_t outer, Int_t turn);
+  void  MakeSeedsMI(Int_t inner, Int_t outer);
 
   Int_t         FollowProlongation(AliTRDtrack& t, Int_t rf);
   Int_t         FollowBackProlongation(AliTRDtrack& t);
@@ -262,7 +279,7 @@ class AliTRDtracker : public AliTracker {
   void          SetSZ2corr(Float_t w)    {fSZ2corr = w;}
   Double_t      ExpectedSigmaY2(Double_t r, Double_t tgl, Double_t pt) const;
   Double_t      ExpectedSigmaZ2(Double_t r, Double_t tgl) const;
-
+  TTreeSRedirector *fDebugStreamer;     //!debug streamer
   ClassDef(AliTRDtracker,2)           // manager base class  
 
 };
index 6deeeb4..35dac45 100644 (file)
@@ -21,6 +21,7 @@
 #pragma link C++ class  AliTRDclusterCorrection+;
 #pragma link C++ class  AliTRDclusterizerMI+;
 
+#pragma link C++ class  AliTRDtracklet+;
 #pragma link C++ class  AliTRDtrack+;
 #pragma link C++ class  AliTRDtracker+;