]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPCtrackerMI.cxx
Gas gain set to 6600 - before it was 20000 (Marian)
[u/mrichter/AliRoot.git] / TPC / AliTPCtrackerMI.cxx
index 44448d118348a91a195216798869a2751c9fb4ca..0b0f2adcc401d8377eb24b47c83aab609fcabf77 100644 (file)
 // The debug level -  different procedure produce tree for numerical debugging
 //                    To enable them set AliTPCReconstructor::SetStreamLevel(n); where nis bigger 1
 //
+
+//
+// Adding systematic errors to the covariance:
+// 
+// The systematic errors due to the misalignment and miscalibration are added to the covariance matrix
+// of the tracks (not to the clusters as they are dependent):
+// The parameters form AliTPCRecoParam are used AliTPCRecoParam::GetSystematicError
+// The systematic errors are expressed there in RMS - position (cm), angle (rad), curvature (1/cm)
+// The default values are 0. 
+//
+// The sytematic errors are added to the covariance matrix in following places:
+//
+// 1. During fisrt itteration - AliTPCtrackerMI::FillESD
+// 2. Second iteration - 
+//      2.a ITS->TPC   - AliTPCtrackerMI::ReadSeeds 
+//      2.b TPC->TRD   - AliTPCtrackerMI::PropagateBack
+// 3. Third iteration  -
+//      3.a TRD->TPC   - AliTPCtrackerMI::ReadSeeds
+//      3.b TPC->ITS   - AliTPCtrackerMI::RefitInward
+//
 // There are several places in the code which can be numerically debuged
 // This code is keeped in order to enable code development and to check the calibration implementtion
 //
 #include "AliTPCReconstructor.h"
 #include "AliTPCpolyTrack.h"
 #include "AliTPCreco.h"
-#include "AliTPCseed.h" 
+#include "AliTPCseed.h"
+
+#include "AliTPCtrackerSector.h" 
 #include "AliTPCtrackerMI.h"
 #include "TStopwatch.h"
 #include "AliTPCReconstructor.h"
 ClassImp(AliTPCtrackerMI)
 
 
+
 class AliTPCFastMath {
 public:
   AliTPCFastMath();  
@@ -280,7 +303,7 @@ Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluste
   Double_t rdistance2  = rdistancey2+rdistancez2;
   //Int_t  accept =0;
   
-  if (AliTPCReconstructor::StreamLevel()>5) {
+  if (AliTPCReconstructor::StreamLevel()>5 && seed->GetNumberOfClusters()>20) {
     Float_t rmsy2 = seed->GetCurrentSigmaY2();
     Float_t rmsz2 = seed->GetCurrentSigmaZ2();
     Float_t rmsy2p30 = seed->GetCMeanSigmaY2p30();
@@ -289,9 +312,18 @@ Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluste
     Float_t rmsz2p30R  = seed->GetCMeanSigmaZ2p30R();
 
     AliExternalTrackParam param(*seed); 
+    static TVectorD gcl(3),gtr(3);
+    Float_t gclf[3];
+    param.GetXYZ(gcl.GetMatrixArray());
+    cluster->GetGlobalXYZ(gclf);
+    gcl[0]=gclf[0];    gcl[1]=gclf[1];    gcl[2]=gclf[2];
+    
+    if (AliTPCReconstructor::StreamLevel()>0) {
     (*fDebugStreamer)<<"ErrParam"<<
       "Cl.="<<cluster<<
       "T.="<<&param<<
+      "gcl.="<<&gcl<<
+      "gtr.="<<&gtr<<
       "erry2="<<sy2<<
       "errz2="<<sz2<<
       "rmsy2="<<rmsy2<<
@@ -301,6 +333,7 @@ Int_t AliTPCtrackerMI::AcceptCluster(AliTPCseed * seed, AliTPCclusterMI * cluste
       "rmsy2p30R="<<rmsy2p30R<<
       "rmsz2p30R="<<rmsz2p30R<<        
       "\n";
+    }
   }
   
   if (rdistance2>16) return 3;
@@ -349,8 +382,8 @@ AliTracker(),
   //---------------------------------------------------------------------
   // The main TPC tracker constructor
   //---------------------------------------------------------------------
-  fInnerSec=new AliTPCSector[fkNIS];         
-  fOuterSec=new AliTPCSector[fkNOS];
+  fInnerSec=new AliTPCtrackerSector[fkNIS];         
+  fOuterSec=new AliTPCtrackerSector[fkNOS];
  
   Int_t i;
   for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
@@ -374,7 +407,9 @@ AliTracker(),
     fYMax[i+nrowlow]      = fXRow[i+nrowlow]*TMath::Tan(0.5*par->GetOuterAngle());
   }
 
-  fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
+  if (AliTPCReconstructor::StreamLevel()>0) {
+    fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
+  }
 }
 //________________________________________________________________________
 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
@@ -438,6 +473,12 @@ void AliTPCtrackerMI::FillESD(TObjArray* arr)
       AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
       if (!pt) continue; 
       pt->UpdatePoints();
+      AddCovariance(pt);
+      if (AliTPCReconstructor::StreamLevel()>1) {
+       (*fDebugStreamer)<<"Track0"<<
+         "Tr0.="<<pt<<
+         "\n";       
+      }
       //      pt->PropagateTo(fParam->GetInnerRadiusLow());
       if (pt->GetKinkIndex(0)<=0){  //don't propagate daughter tracks 
        pt->PropagateTo(fParam->GetInnerRadiusLow());
@@ -1074,7 +1115,7 @@ Int_t  AliTPCtrackerMI::LoadClusters(TObjArray *arr)
     }
     //
     if (clrow->GetArray()->GetEntriesFast()<=0) continue;
-    AliTPCRow * tpcrow=0;
+    AliTPCtrackerRow * tpcrow=0;
     Int_t left=0;
     if (sec<fkNIS*2){
       tpcrow = &(fInnerSec[sec%fkNIS][row]);    
@@ -1104,6 +1145,63 @@ Int_t  AliTPCtrackerMI::LoadClusters(TObjArray *arr)
   return 0;
 }
 
+Int_t  AliTPCtrackerMI::LoadClusters(TClonesArray *arr)
+{
+  //
+  // load clusters to the memory from one 
+  // TClonesArray
+  //
+  AliTPCclusterMI *clust=0;
+  Int_t count[72][96] = { {0} , {0} }; 
+
+  // loop over clusters
+  for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
+    clust = (AliTPCclusterMI*)arr->At(icl);
+    if(!clust) continue;
+    //printf("cluster: det %d, row %d \n", clust->GetDetector(),clust->GetRow());
+
+    // transform clusters
+    Transform(clust);
+
+    // count clusters per pad row
+    count[clust->GetDetector()][clust->GetRow()]++;
+  }
+
+  // insert clusters to sectors
+  for (Int_t icl=0; icl<arr->GetEntriesFast(); icl++) {
+    clust = (AliTPCclusterMI*)arr->At(icl);
+    if(!clust) continue;
+
+    Int_t sec = clust->GetDetector();
+    Int_t row = clust->GetRow();
+
+    // filter overlapping pad rows needed by HLT
+    if(sec<fkNIS*2) { //IROCs
+     if(row == 30) continue;
+    }
+    else { // OROCs
+      if(row == 27 || row == 76) continue;
+    }
+
+    Int_t left=0;
+    if (sec<fkNIS*2){
+      left = sec/fkNIS;
+      fInnerSec[sec%fkNIS].InsertCluster(clust, count[sec][row], fParam);    
+    }
+    else{
+      left = (sec-fkNIS*2)/fkNOS;
+      fOuterSec[(sec-fkNIS*2)%fkNOS].InsertCluster(clust, count[sec][row], fParam);
+    }
+  }
+
+  // Load functions must be called behind LoadCluster(TClonesArray*)
+  // needed by HLT
+  //LoadOuterSectors();
+  //LoadInnerSectors();
+
+  return 0;
+}
+
 
 Int_t  AliTPCtrackerMI::LoadClusters()
 {
@@ -1130,7 +1228,7 @@ Int_t  AliTPCtrackerMI::LoadClusters()
       Transform((AliTPCclusterMI*)(clrow->GetArray()->At(icl)));
     }
     //
-    AliTPCRow * tpcrow=0;
+    AliTPCtrackerRow * tpcrow=0;
     Int_t left=0;
     if (sec<fkNIS*2){
       tpcrow = &(fInnerSec[sec%fkNIS][row]);    
@@ -1169,7 +1267,7 @@ void AliTPCtrackerMI::UnloadClusters()
   Int_t nrows = fOuterSec->GetNRows();
   for (Int_t sec = 0;sec<fkNOS;sec++)
     for (Int_t row = 0;row<nrows;row++){
-      AliTPCRow*  tpcrow = &(fOuterSec[sec%fkNOS][row]);
+      AliTPCtrackerRow*  tpcrow = &(fOuterSec[sec%fkNOS][row]);
       //      if (tpcrow){
       //       if (tpcrow->fClusters1) delete []tpcrow->fClusters1; 
       //       if (tpcrow->fClusters2) delete []tpcrow->fClusters2; 
@@ -1180,7 +1278,7 @@ void AliTPCtrackerMI::UnloadClusters()
   nrows = fInnerSec->GetNRows();
   for (Int_t sec = 0;sec<fkNIS;sec++)
     for (Int_t row = 0;row<nrows;row++){
-      AliTPCRow*  tpcrow = &(fInnerSec[sec%fkNIS][row]);
+      AliTPCtrackerRow*  tpcrow = &(fInnerSec[sec%fkNIS][row]);
       //if (tpcrow){
       //       if (tpcrow->fClusters1) delete []tpcrow->fClusters1; 
       //if (tpcrow->fClusters2) delete []tpcrow->fClusters2; 
@@ -1191,6 +1289,32 @@ void AliTPCtrackerMI::UnloadClusters()
   return ;
 }
 
+void AliTPCtrackerMI::FillClusterArray(TObjArray* array) const{
+  //
+  // Filling cluster to the array - For visualization purposes
+  //
+  Int_t nrows=0;
+  nrows = fOuterSec->GetNRows();
+  for (Int_t sec = 0;sec<fkNOS;sec++)
+    for (Int_t row = 0;row<nrows;row++){
+      AliTPCtrackerRow*  tpcrow = &(fOuterSec[sec%fkNOS][row]);
+      if (!tpcrow) continue;
+      for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
+       array->AddLast((TObject*)((*tpcrow)[icl]));
+      }
+    } 
+  nrows = fInnerSec->GetNRows();
+  for (Int_t sec = 0;sec<fkNIS;sec++)
+    for (Int_t row = 0;row<nrows;row++){
+      AliTPCtrackerRow*  tpcrow = &(fInnerSec[sec%fkNIS][row]);
+      if (!tpcrow) continue;
+      for (Int_t icl = 0;icl<tpcrow->GetN();icl++){
+       array->AddLast((TObject*)(*tpcrow)[icl]);
+      }
+    }
+}
+
+
 void   AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
   //
   //
@@ -1202,11 +1326,9 @@ void   AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
   Double_t x[3]={cluster->GetRow(),cluster->GetPad(),cluster->GetTimeBin()};
   Int_t i[1]={cluster->GetDetector()};
   transform->Transform(x,i,0,1);  
-  if (!AliTPCReconstructor::GetRecoParam()->GetBYMirror()){
-    if (cluster->GetDetector()%36>17){
-      x[1]*=-1;
-    }
-  }
+  //  if (cluster->GetDetector()%36>17){
+  //  x[1]*=-1;
+  //}
 
   //
   // in debug mode  check the transformation
@@ -1214,8 +1336,10 @@ void   AliTPCtrackerMI::Transform(AliTPCclusterMI * cluster){
   if (AliTPCReconstructor::StreamLevel()>1) {
     Float_t gx[3];
     cluster->GetGlobalXYZ(gx);
+    Int_t event = (fEvent==NULL)? 0: fEvent->GetEventNumberInFile();
     TTreeSRedirector &cstream = *fDebugStreamer;
     cstream<<"Transform"<<
+      "event="<<event<<
       "x0="<<x[0]<<
       "x1="<<x[1]<<
       "x2="<<x[2]<<
@@ -1256,7 +1380,7 @@ Int_t AliTPCtrackerMI::LoadOuterSectors() {
   UInt_t index=0;
   for (Int_t sec = 0;sec<fkNOS;sec++)
     for (Int_t row = 0;row<nrows;row++){
-      AliTPCRow*  tpcrow = &(fOuterSec[sec%fkNOS][row]);  
+      AliTPCtrackerRow*  tpcrow = &(fOuterSec[sec%fkNOS][row]);  
       Int_t sec2 = sec+2*fkNIS;
       //left
       Int_t ncl = tpcrow->GetN1();
@@ -1304,7 +1428,7 @@ Int_t  AliTPCtrackerMI::LoadInnerSectors() {
   UInt_t index=0;
   for (Int_t sec = 0;sec<fkNIS;sec++)
     for (Int_t row = 0;row<nrows;row++){
-      AliTPCRow*  tpcrow = &(fInnerSec[sec%fkNIS][row]);
+      AliTPCtrackerRow*  tpcrow = &(fInnerSec[sec%fkNIS][row]);
       //
       //left
       Int_t ncl = tpcrow->GetN1();
@@ -1356,7 +1480,7 @@ AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
   Int_t row=(index&0x00ff0000)>>16; 
   Int_t ncl=(index&0x00007fff)>>00;
 
-  const AliTPCRow * tpcrow=0;
+  const AliTPCtrackerRow * tpcrow=0;
   AliTPCclusterMI * clrow =0;
 
   if (sec<0 || sec>=fkNIS*4) {
@@ -1493,7 +1617,7 @@ Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
   
   if (!isActive || !isActive2) return 0;
 
-  const AliTPCRow &krow=GetRow(t.GetRelativeSector(),nr);
+  const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
   if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
   Double_t  roady  =1.;
   Double_t  roadz = 1.;
@@ -1590,7 +1714,7 @@ Int_t AliTPCtrackerMI::FollowToNextFast(AliTPCseed& t, Int_t nr) {
   
   
   //Int_t nr2 = nr;
-  const AliTPCRow &krow=GetRow(t.GetRelativeSector(),nr);
+  const AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
   if ( (t.GetSigmaY2()<0) || t.GetSigmaZ2()<0) return 0;
   Double_t  roady  =1.;
   Double_t  roadz = 1.;
@@ -1732,7 +1856,7 @@ Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t,  Int_t nr) {
   }
   //AliInfo(Form("A - Sector%d phi %f - alpha %f", t.fRelativeSector,y/x, t.GetAlpha()));
 
-  AliTPCRow &krow=GetRow(t.GetRelativeSector(),nr);
+  AliTPCtrackerRow &krow=GetRow(t.GetRelativeSector(),nr);
 
   if (TMath::Abs(TMath::Abs(y)-ymax)<krow.GetDeadZone()){
     t.SetInDead(kTRUE);
@@ -2518,6 +2642,7 @@ Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
 
     seed->PropagateTo(fParam->GetInnerRadiusLow());
     seed->UpdatePoints();
+    AddCovariance(seed);
     MakeBitmaps(seed);
     AliESDtrack *esd=event->GetTrack(i);
     seed->CookdEdx(0.02,0.6);
@@ -2557,7 +2682,6 @@ Int_t AliTPCtrackerMI::RefitInward(AliESDEvent *event)
   }
   //FindKinks(fSeeds,event);
   Info("RefitInward","Number of refitted tracks %d",ntracks);
-  fEvent =0;
   return 0;
 }
 
@@ -2584,6 +2708,7 @@ Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
     if (!seed) continue;
     if (seed->GetKinkIndex(0)<0)  UpdateKinkQualityM(seed);  // update quality informations for kinks
     seed->UpdatePoints();
+    AddCovariance(seed);
     AliESDtrack *esd=event->GetTrack(i);
     seed->CookdEdx(0.02,0.6);
     CookLabel(seed,0.1); //For comparison only
@@ -2598,11 +2723,12 @@ Int_t AliTPCtrackerMI::PropagateBack(AliESDEvent *event)
       ntracks++;
       Int_t eventnumber = event->GetEventNumberInFile();// patch 28 fev 06
       // This is most likely NOT the event number you'd like to use. It has nothing to do with the 'real' event number      
-      if (AliTPCReconstructor::StreamLevel()>1) {
+      if (AliTPCReconstructor::StreamLevel()>1 && esd) {
        (*fDebugStreamer)<<"Cback"<<
          "Tr0.="<<seed<<
+         "esd.="<<esd<<
          "EventNrInFile="<<eventnumber<<
-         "\n"; // patch 28 fev 06         
+         "\n";       
       }
     }
   }
@@ -2654,6 +2780,7 @@ void AliTPCtrackerMI::ReadSeeds(AliESDEvent *event, Int_t direction)
     //    AliTPCseed *seed = new AliTPCseed(t,t.GetAlpha());
     AliTPCseed *seed = new AliTPCseed(t/*,t.GetAlpha()*/);
     seed->SetUniqueID(esd->GetID());
+    AddCovariance(seed);   //add systematic ucertainty
     for (Int_t ikink=0;ikink<3;ikink++) {
       Int_t index = esd->GetKinkIndex(ikink);
       seed->GetKinkIndexes()[ikink] = index;
@@ -2769,11 +2896,11 @@ void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2,
 
   Int_t imiddle = (i2+i1)/2;    //middle pad row index
   Double_t xm = GetXrow(imiddle); // radius of middle pad-row
-  const AliTPCRow& krm=GetRow(sec,imiddle); //middle pad -row
+  const AliTPCtrackerRow& krm=GetRow(sec,imiddle); //middle pad -row
   //
   Int_t ns =sec;   
 
-  const AliTPCRow& kr1=GetRow(ns,i1);
+  const AliTPCtrackerRow& kr1=GetRow(ns,i1);
   Double_t ymax  = GetMaxY(i1)-kr1.GetDeadZone()-1.5;  
   Double_t ymaxm = GetMaxY(imiddle)-kr1.GetDeadZone()-1.5;  
 
@@ -2825,10 +2952,10 @@ void AliTPCtrackerMI::MakeSeeds3(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2,
     for (Int_t dsec = dsec1; dsec<=dsec2;dsec++){
       Int_t sec2 = sec + dsec;
       // 
-      //      AliTPCRow&  kr2  = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
-      //AliTPCRow&  kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
-      AliTPCRow&  kr2  = GetRow((sec2+fkNOS)%fkNOS,i2);
-      AliTPCRow&  kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
+      //      AliTPCtrackerRow&  kr2  = fOuterSec[(sec2+fkNOS)%fkNOS][i2];
+      //AliTPCtrackerRow&  kr2m = fOuterSec[(sec2+fkNOS)%fkNOS][imiddle];
+      AliTPCtrackerRow&  kr2  = GetRow((sec2+fkNOS)%fkNOS,i2);
+      AliTPCtrackerRow&  kr2m = GetRow((sec2+fkNOS)%fkNOS,imiddle);
       Int_t  index1 = TMath::Max(kr2.Find(extraz-0.6-dddz1*TMath::Abs(z1)*0.05)-1,0);
       Int_t  index2 = TMath::Min(kr2.Find(extraz+0.6+dddz2*TMath::Abs(z1)*0.05)+1,kr2);
 
@@ -3081,25 +3208,25 @@ void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2,
 
   // first 3 padrows
   Double_t x1 = GetXrow(i1-1);
-  const    AliTPCRow& kr1=GetRow(sec,i1-1);
+  const    AliTPCtrackerRow& kr1=GetRow(sec,i1-1);
   Double_t y1max  = GetMaxY(i1-1)-kr1.GetDeadZone()-1.5;  
   //
   Double_t x1p = GetXrow(i1);
-  const    AliTPCRow& kr1p=GetRow(sec,i1);
+  const    AliTPCtrackerRow& kr1p=GetRow(sec,i1);
   //
   Double_t x1m = GetXrow(i1-2);
-  const    AliTPCRow& kr1m=GetRow(sec,i1-2);
+  const    AliTPCtrackerRow& kr1m=GetRow(sec,i1-2);
 
   //
   //last 3 padrow for seeding
-  AliTPCRow&  kr3  = GetRow((sec+fkNOS)%fkNOS,i1-7);
+  AliTPCtrackerRow&  kr3  = GetRow((sec+fkNOS)%fkNOS,i1-7);
   Double_t    x3   =  GetXrow(i1-7);
   //  Double_t    y3max= GetMaxY(i1-7)-kr3.fDeadZone-1.5;  
   //
-  AliTPCRow&  kr3p  = GetRow((sec+fkNOS)%fkNOS,i1-6);
+  AliTPCtrackerRow&  kr3p  = GetRow((sec+fkNOS)%fkNOS,i1-6);
   Double_t    x3p   = GetXrow(i1-6);
   //
-  AliTPCRow&  kr3m  = GetRow((sec+fkNOS)%fkNOS,i1-8);
+  AliTPCtrackerRow&  kr3m  = GetRow((sec+fkNOS)%fkNOS,i1-8);
   Double_t    x3m   = GetXrow(i1-8);
 
   //
@@ -3107,7 +3234,7 @@ void AliTPCtrackerMI::MakeSeeds5(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2,
   // middle padrow
   Int_t im = i1-4;                           //middle pad row index
   Double_t xm         = GetXrow(im);         // radius of middle pad-row
-  const AliTPCRow& krm=GetRow(sec,im);   //middle pad -row
+  const AliTPCtrackerRow& krm=GetRow(sec,im);   //middle pad -row
   //  Double_t ymmax = GetMaxY(im)-kr1.fDeadZone-1.5;  
   //
   //
@@ -3331,8 +3458,8 @@ void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2,
   //  Double_t cs=cos(alpha), sn=sin(alpha);
   Int_t row0 = (i1+i2)/2;
   Int_t drow = (i1-i2)/2;
-  const AliTPCRow& kr0=fSectors[sec][row0];
-  AliTPCRow * kr=0;
+  const AliTPCtrackerRow& kr0=fSectors[sec][row0];
+  AliTPCtrackerRow * kr=0;
 
   AliTPCpolyTrack polytrack;
   Int_t nclusters=fSectors[sec][row0];
@@ -3345,8 +3472,8 @@ void AliTPCtrackerMI::MakeSeeds2(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2,
     Int_t nfound =0;
     Int_t nfoundable =0;
     for (Int_t iter =1; iter<2; iter++){   //iterations
-      const AliTPCRow& krm=fSectors[sec][row0-iter];
-      const AliTPCRow& krp=fSectors[sec][row0+iter];      
+      const AliTPCtrackerRow& krm=fSectors[sec][row0-iter];
+      const AliTPCtrackerRow& krp=fSectors[sec][row0+iter];      
       const AliTPCclusterMI * cl= kr0[is];
       
       if (cl->IsUsed(10)) {
@@ -4034,7 +4161,6 @@ void  AliTPCtrackerMI::FindMultiMC(TObjArray * array, AliESDEvent */*esd*/, Int_
     }
     if (ncl>0) xm[i]/=Float_t(ncl);
   }  
-  TTreeSRedirector &cstream = *fDebugStreamer;
   //
   for (Int_t i0=0;i0<nentries;i0++){
     AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
@@ -4095,6 +4221,8 @@ void  AliTPCtrackerMI::FindMultiMC(TObjArray * array, AliESDEvent */*esd*/, Int_
        }
       }
       //
+      if (AliTPCReconstructor::StreamLevel()>0) {
+      TTreeSRedirector &cstream = *fDebugStreamer;
       cstream<<"Multi"<<
        "iter="<<iter<<
        "lab0="<<lab0<<
@@ -4130,6 +4258,7 @@ void  AliTPCtrackerMI::FindMultiMC(TObjArray * array, AliESDEvent */*esd*/, Int_
        "r1="<<r1<<
        "rc1="<<rc1<<
        "\n";
+       }
     }
   }    
   delete [] helixes;
@@ -4216,7 +4345,6 @@ void  AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int
     }
     if (ncl>0) xm[i]/=Float_t(ncl);
   }  
-  TTreeSRedirector &cstream = *fDebugStreamer;
   //
   for (Int_t is0=0;is0<nentries;is0++){
     Int_t i0 = indexes[is0];
@@ -4298,39 +4426,43 @@ void  AliTPCtrackerMI::FindSplitted(TObjArray * array, AliESDEvent */*esd*/, Int
       //
       Int_t lab0=track0->GetLabel();
       Int_t lab1=track1->GetLabel();
-      cstream<<"Splitted"<<
-       "iter="<<iter<<
-       "lab0="<<lab0<<
-       "lab1="<<lab1<<   
-       "Tr0.="<<track0<<       // seed0
-       "Tr1.="<<track1<<       // seed1
-       "h0.="<<&helixes[i0]<<
-       "h1.="<<&helixes[i1]<<
-       //
-       "sum="<<sum<<           //the sum of rows with cl in both
-       "sum0="<<sum0<<           //the sum of rows with cl in 0 track
-       "sum1="<<sum1<<           //the sum of rows with cl in 1 track
-       "sums="<<sums<<         //the sum of shared clusters
-       "xm0="<<xm[i0]<<        // the center of track
-       "xm1="<<xm[i1]<<        // the x center of track
-       // General cut variables                   
-       "dfi="<<dfi<<           // distance in fi angle
-       "dtheta="<<dtheta<<     // distance int theta angle
-       //
-       //
-       "dist0="<<dist[4][0]<<     //distance x
-       "dist1="<<dist[4][1]<<     //distance y
-       "dist2="<<dist[4][2]<<     //distance z
-       "mdist0="<<mdist[0]<<   //distance x
-       "mdist1="<<mdist[1]<<   //distance y
-       "mdist2="<<mdist[2]<<   //distance z
-
-       "\n";
+      if( AliTPCReconstructor::StreamLevel()>5){
+      TTreeSRedirector &cstream = *fDebugStreamer;
+       cstream<<"Splitted"<<
+         "iter="<<iter<<
+         "lab0="<<lab0<<
+         "lab1="<<lab1<<   
+         "Tr0.="<<track0<<       // seed0
+         "Tr1.="<<track1<<       // seed1
+         "h0.="<<&helixes[i0]<<
+         "h1.="<<&helixes[i1]<<
+         //
+         "sum="<<sum<<           //the sum of rows with cl in both
+         "sum0="<<sum0<<           //the sum of rows with cl in 0 track
+         "sum1="<<sum1<<           //the sum of rows with cl in 1 track
+         "sums="<<sums<<         //the sum of shared clusters
+         "xm0="<<xm[i0]<<        // the center of track
+         "xm1="<<xm[i1]<<        // the x center of track
+         // General cut variables                   
+         "dfi="<<dfi<<           // distance in fi angle
+         "dtheta="<<dtheta<<     // distance int theta angle
+         //
+         //
+         "dist0="<<dist[4][0]<<     //distance x
+         "dist1="<<dist[4][1]<<     //distance y
+         "dist2="<<dist[4][2]<<     //distance z
+         "mdist0="<<mdist[0]<<   //distance x
+         "mdist1="<<mdist[1]<<   //distance y
+         "mdist2="<<mdist[2]<<   //distance z    
+         "\n";
+      }
       delete array->RemoveAt(i1);
     }
   }    
   delete [] helixes;
   delete [] xm;
+  delete [] quality;
+  delete [] indexes;
   AliInfo("Time for splitted tracks removal");
   timer.Print();
 }
@@ -4402,7 +4534,6 @@ void  AliTPCtrackerMI::FindCurling(TObjArray * array, AliESDEvent */*esd*/, Int_
   //
   // Find tracks
   //
-  TTreeSRedirector &cstream = *fDebugStreamer;
   //
   for (Int_t i0=0;i0<nentries;i0++){
     AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
@@ -4500,6 +4631,7 @@ void  AliTPCtrackerMI::FindCurling(TObjArray * array, AliESDEvent */*esd*/, Int_
          //debug stream to tune "fine" cuts      
          Int_t lab0=track0->GetLabel();
          Int_t lab1=track1->GetLabel();
+          TTreeSRedirector &cstream = *fDebugStreamer;
          cstream<<"Curling2"<<
            "iter="<<iter<<
            "lab0="<<lab0<<
@@ -4614,7 +4746,6 @@ void  AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
 
   //
   // Find circling track
-  TTreeSRedirector &cstream = *fDebugStreamer;
   //
   for (Int_t i0=0;i0<nentries;i0++){
     AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
@@ -4701,6 +4832,7 @@ void  AliTPCtrackerMI::FindKinks(TObjArray * array, AliESDEvent *esd)
          //debug stream          
          Int_t lab0=track0->GetLabel();
          Int_t lab1=track1->GetLabel();
+          TTreeSRedirector &cstream = *fDebugStreamer;
          cstream<<"Curling"<<
            "lab0="<<lab0<<
            "lab1="<<lab1<<   
@@ -5305,7 +5437,6 @@ void  AliTPCtrackerMI::FindV0s(TObjArray * array, AliESDEvent *esd)
   // 
   //
   // //  
-  TTreeSRedirector &cstream = *fDebugStreamer; 
   Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
   AliV0 vertex; 
   Double_t cradius0 = 10*10;
@@ -5320,6 +5451,7 @@ void  AliTPCtrackerMI::FindV0s(TObjArray * array, AliESDEvent *esd)
     AliTPCseed * track0 = (AliTPCseed*)array->At(i);
     if (!track0) continue;
     if (AliTPCReconstructor::StreamLevel()>1){
+      TTreeSRedirector &cstream = *fDebugStreamer;
       cstream<<"Tracks"<<
        "Tr0.="<<track0<<
        "dca="<<dca[i]<<
@@ -5433,6 +5565,7 @@ void  AliTPCtrackerMI::FindV0s(TObjArray * array, AliESDEvent *esd)
          Int_t lab1=track1->GetLabel();
          Char_t c0=track0->GetCircular();
          Char_t c1=track1->GetCircular();
+          TTreeSRedirector &cstream = *fDebugStreamer;
          cstream<<"V0"<<
          "Event="<<eventNr<<
          "vertex.="<<&vertex<<
@@ -5543,6 +5676,7 @@ void  AliTPCtrackerMI::FindV0s(TObjArray * array, AliESDEvent *esd)
       if (AliTPCReconstructor::StreamLevel()>1) {
        Int_t lab0=track0->GetLabel();
        Int_t lab1=track1->GetLabel();
+        TTreeSRedirector &cstream = *fDebugStreamer;
        cstream<<"V02"<<
        "Event="<<eventNr<<
        "vertex.="<<v0<<        
@@ -6834,29 +6968,6 @@ Int_t AliTPCtrackerMI::CookLabel(AliTPCseed *t, Float_t wrong,Int_t first, Int_t
 }
 
 
-Int_t  AliTPCtrackerMI::AliTPCSector::GetRowNumber(Double_t x) const 
-{
-  //return pad row number for this x
-  Double_t r;
-  if (fN < 64){
-    r=fRow[fN-1].GetX();
-    if (x > r) return fN;
-    r=fRow[0].GetX();
-    if (x < r) return -1;
-    return Int_t((x-r)/fPadPitchLength + 0.5);}
-  else{    
-    r=fRow[fN-1].GetX();
-    if (x > r) return fN;
-    r=fRow[0].GetX();
-    if (x < r) return -1;
-    Double_t r1=fRow[64].GetX();
-    if(x<r1){       
-      return Int_t((x-r)/f1PadPitchLength + 0.5);}
-    else{
-      return (Int_t((x-r1)/f2PadPitchLength + 0.5)+64);} 
-  }
-}
-
 Int_t  AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const 
 {
   //return pad row number for given x vector
@@ -6869,193 +6980,6 @@ Int_t  AliTPCtrackerMI::GetRowNumber(Double_t x[3]) const
   return GetRowNumber(localx);
 }
 
-//_________________________________________________________________________
-void AliTPCtrackerMI::AliTPCSector::Setup(const AliTPCParam *par, Int_t f) {
-  //-----------------------------------------------------------------------
-  // Setup inner sector
-  //-----------------------------------------------------------------------
-  if (f==0) {
-     fAlpha=par->GetInnerAngle();
-     fAlphaShift=par->GetInnerAngleShift();
-     fPadPitchWidth=par->GetInnerPadPitchWidth();
-     fPadPitchLength=par->GetInnerPadPitchLength();
-     fN=par->GetNRowLow();
-     if(fRow)delete [] fRow;fRow = 0;
-     fRow=new AliTPCRow[fN];
-     for (Int_t i=0; i<fN; i++) {
-       fRow[i].SetX(par->GetPadRowRadiiLow(i));
-       fRow[i].SetDeadZone(1.5);  //1.5 cm of dead zone
-     }
-  } else {
-     fAlpha=par->GetOuterAngle();
-     fAlphaShift=par->GetOuterAngleShift();
-     fPadPitchWidth  = par->GetOuterPadPitchWidth();
-     fPadPitchLength = par->GetOuter1PadPitchLength();
-     f1PadPitchLength = par->GetOuter1PadPitchLength();
-     f2PadPitchLength = par->GetOuter2PadPitchLength();
-     fN=par->GetNRowUp();
-     if(fRow)delete [] fRow;fRow = 0;
-     fRow=new AliTPCRow[fN];
-     for (Int_t i=0; i<fN; i++) {
-       fRow[i].SetX(par->GetPadRowRadiiUp(i)); 
-       fRow[i].SetDeadZone(1.5);  // 1.5 cm of dead zone
-     }
-  } 
-}
-
-AliTPCtrackerMI::AliTPCRow::AliTPCRow():
-  fDeadZone(0.),
-  fClusters1(0),
-  fN1(0),
-  fClusters2(0),
-  fN2(0),
-  fN(0),
-  fX(0.)
-{
-  //
-  // default constructor
-  //
-}
-
-AliTPCtrackerMI::AliTPCRow::~AliTPCRow(){
-  //
-  for (Int_t i = 0; i < fN1; i++)
-    fClusters1[i].~AliTPCclusterMI();
-  delete [] fClusters1;
-  for (Int_t i = 0; i < fN2; i++)
-    fClusters2[i].~AliTPCclusterMI();
-  delete [] fClusters2;
-}
-
-
-
-//_________________________________________________________________________
-void 
-AliTPCtrackerMI::AliTPCRow::InsertCluster(const AliTPCclusterMI* c, UInt_t index) {
-  //-----------------------------------------------------------------------
-  // Insert a cluster into this pad row in accordence with its y-coordinate
-  //-----------------------------------------------------------------------
-  if (fN==kMaxClusterPerRow) {
-    //AliInfo("AliTPCRow::InsertCluster(): Too many clusters"); 
-    return;
-  }
-  if (fN>=fN1+fN2) {
-    //AliInfo("AliTPCRow::InsertCluster(): Too many clusters !");
-  }
-
-  if (fN==0) {fIndex[0]=index; fClusters[fN++]=c; return;}
-  Int_t i=Find(c->GetZ());
-  memmove(fClusters+i+1 ,fClusters+i,(fN-i)*sizeof(AliTPCclusterMI*));
-  memmove(fIndex   +i+1 ,fIndex   +i,(fN-i)*sizeof(UInt_t));
-  fIndex[i]=index; fClusters[i]=c; fN++;
-}
-
-void AliTPCtrackerMI::AliTPCRow::ResetClusters() {
-   //
-   // reset clusters
-   // MvL: Need to call destructors for AliTPCclusterMI, to delete fInfo
-   for (Int_t i = 0; i < fN1; i++)
-     fClusters1[i].~AliTPCclusterMI();
-   delete [] fClusters1;  fClusters1=0;
-   for (Int_t i = 0; i < fN2; i++)
-     fClusters2[i].~AliTPCclusterMI();
-   delete [] fClusters2;  fClusters2=0;
-
-   fN  = 0; 
-   fN1 = 0;
-   fN2 = 0;
-   //delete[] fClusterArray; 
-
-   //fClusterArray=0;
-}
-
-
-//___________________________________________________________________
-Int_t AliTPCtrackerMI::AliTPCRow::Find(Double_t z) const {
-  //-----------------------------------------------------------------------
-  // Return the index of the nearest cluster 
-  //-----------------------------------------------------------------------
-  if (fN==0) return 0;
-  if (z <= fClusters[0]->GetZ()) return 0;
-  if (z > fClusters[fN-1]->GetZ()) return fN;
-  Int_t b=0, e=fN-1, m=(b+e)/2;
-  for (; b<e; m=(b+e)/2) {
-    if (z > fClusters[m]->GetZ()) b=m+1;
-    else e=m; 
-  }
-  return m;
-}
-
-
-
-//___________________________________________________________________
-AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest(Double_t y, Double_t z, Double_t roady, Double_t roadz) const {
-  //-----------------------------------------------------------------------
-  // Return the index of the nearest cluster in z y 
-  //-----------------------------------------------------------------------
-  Float_t maxdistance = roady*roady + roadz*roadz;
-
-  AliTPCclusterMI *cl =0;
-  for (Int_t i=Find(z-roadz); i<fN; i++) {
-      AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
-      if (c->GetZ() > z+roadz) break;
-      if ( (c->GetY()-y) >  roady ) continue;
-      Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
-      if (maxdistance>distance) {
-       maxdistance = distance;
-       cl=c;       
-      }
-  }
-  return cl;      
-}
-
-AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest2(Double_t y, Double_t z, Double_t roady, Double_t roadz,UInt_t & index) const 
-{
-  //-----------------------------------------------------------------------
-  // Return the index of the nearest cluster in z y 
-  //-----------------------------------------------------------------------
-  Float_t maxdistance = roady*roady + roadz*roadz;
-  AliTPCclusterMI *cl =0;
-
-  //PH Check boundaries. 510 is the size of fFastCluster
-  Int_t iz1 = Int_t(z-roadz+254.5);
-  if (iz1<0 || iz1>=510) return cl;
-  iz1 = TMath::Max(GetFastCluster(iz1)-1,0);
-  Int_t iz2 = Int_t(z+roadz+255.5);
-  if (iz2<0 || iz2>=510) return cl;
-  iz2 = TMath::Min(GetFastCluster(iz2)+1,fN);
-  Bool_t skipUsed = !(AliTPCReconstructor::GetRecoParam()->GetClusterSharing());
-  //FindNearest3(y,z,roady,roadz,index);
-  //  for (Int_t i=Find(z-roadz); i<fN; i++) {
-  for (Int_t i=iz1; i<iz2; i++) {
-      AliTPCclusterMI *c=(AliTPCclusterMI*)(fClusters[i]);
-      if (c->GetZ() > z+roadz) break;
-      if ( c->GetY()-y >  roady ) continue;
-      if ( y-c->GetY() >  roady ) continue;
-      if (skipUsed && c->IsUsed(11)) continue;
-      Float_t distance = (c->GetZ()-z)*(c->GetZ()-z)+(c->GetY()-y)*(c->GetY()-y);
-      if (maxdistance>distance) {
-       maxdistance = distance;
-       cl=c;       
-       index =i;
-       //roady = TMath::Sqrt(maxdistance);
-      }
-  }
-  return cl;      
-}
-
-
-void AliTPCtrackerMI::AliTPCRow::SetFastCluster(Int_t i, Short_t cl){
-  //
-  // Set cluster info for fast navigation
-  //
-  if (i>510|| i<0){
-  }else{
-    fFastCluster[i]=cl;
-  }
-}
-
-
 
 
 void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
@@ -7083,3 +7007,27 @@ void AliTPCtrackerMI::MakeBitmaps(AliTPCseed *t)
     }
   }
 }
+
+
+
+void AliTPCtrackerMI::AddCovariance(AliTPCseed * seed){
+  //
+  // Adding systematic error
+  // !!!! the systematic error for element 4 is in 1/cm not in pt 
+
+  const Double_t *param = AliTPCReconstructor::GetRecoParam()->GetSystematicError();
+  Double_t covar[15];
+  for (Int_t i=0;i<15;i++) covar[i]=0;
+  // 0
+  // 1    2
+  // 3    4    5
+  // 6    7    8    9 
+  // 10   11   12   13   14
+  covar[0] = param[0]*param[0];
+  covar[2] = param[1]*param[1];
+  covar[5] = param[2]*param[2];
+  covar[9] = param[3]*param[3];
+  Double_t facC =  AliTracker::GetBz()*kB2C;
+  covar[14]= param[4]*param[4]*facC*facC;
+  seed->AddCovariance(covar);
+}