]> git.uio.no Git - u/mrichter/AliRoot.git/blobdiff - TPC/AliTPCtracker.cxx
Vertexing with tracks included (Andrea)
[u/mrichter/AliRoot.git] / TPC / AliTPCtracker.cxx
index 4f49515ec2c05d195fca812297d60de6a854deb2..e5891f865a08adee2f143d23695b3df34068aa02 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/*
-$Log$
-Revision 1.27  2003/02/25 16:47:58  hristov
-allow back propagation for more than 1 event (J.Chudoba)
-
-Revision 1.26  2003/02/19 08:49:46  hristov
-Track time measurement (S.Radomski)
-
-Revision 1.25  2003/01/28 16:43:35  hristov
-Additional protection: to be discussed with the Root team (M.Ivanov)
-
-Revision 1.24  2002/11/19 16:13:24  hristov
-stdlib.h included to declare exit() on HP
-
-Revision 1.23  2002/11/19 11:50:08  hristov
-Removing CONTAINERS (Yu.Belikov)
-
-Revision 1.19  2002/07/19 07:31:40  kowal2
-Improvement in tracking by J. Belikov
-
-Revision 1.18  2002/05/13 07:33:52  kowal2
-Added protection in Int_t AliTPCtracker::AliTPCRow::Find(Double_t y) const
-in the case of defined region of interests.
-
-Revision 1.17  2002/03/18 17:59:13  kowal2
-Chnges in the pad geometry - 3 pad lengths introduced.
-
-Revision 1.16  2001/11/08 16:39:03  hristov
-Additional protection (M.Masera)
-
-Revision 1.15  2001/11/08 16:36:33  hristov
-Updated V2 stream of tracking (Yu.Belikov). The new long waited features are: 1) Possibility to pass the primary vertex position to the trackers (both for the TPC and the ITS) 2) Possibility to specify the number of tracking passes together with applying (or not applying) the vertex constraint (ITS only) 3) Possibility to make some use of partial PID provided by the TPC when doing tracking in the ITS (ITS only) 4) V0 reconstruction with a helix minimisation of the DCA. (new macros: AliV0FindVertices.C and AliV0Comparison.C) 4a) ( Consequence of the 4) )  All the efficiencies and resolutions are from now on calculated including *secondary*particles* too. (Don't be surprised by the drop in efficiency etc)
-
-Revision 1.14  2001/10/21 19:04:55  hristov
-Several patches were done to adapt the barel reconstruction to the multi-event case. Some memory leaks were corrected. (Yu.Belikov)
-
-Revision 1.13  2001/07/23 12:01:30  hristov
-Initialisation added
-
-Revision 1.12  2001/07/20 14:32:44  kowal2
-Processing of many events possible now
-
-Revision 1.11  2001/05/23 08:50:10  hristov
-Weird inline removed
-
-Revision 1.10  2001/05/16 14:57:25  alibrary
-New files for folders and Stack
-
-Revision 1.9  2001/05/11 07:16:56  hristov
-Fix needed on Sun and Alpha
-
-Revision 1.8  2001/05/08 15:00:15  hristov
-Corrections for tracking in arbitrary magnenetic field. Changes towards a concept of global Alice track. Back propagation of reconstructed tracks (Yu.Belikov)
-
-Revision 1.5  2000/12/20 07:51:59  kowal2
-Changes suggested by Alessandra and Paolo to avoid overlapped
-data fields in encapsulated classes.
-
-Revision 1.4  2000/11/02 07:27:16  kowal2
-code corrections
-
-Revision 1.2  2000/06/30 12:07:50  kowal2
-Updated from the TPC-PreRelease branch
-
-Revision 1.1.2.1  2000/06/25 08:53:55  kowal2
-Splitted from AliTPCtracking
-
-*/
+/* $Id$ */
 
 //-------------------------------------------------------
 //          Implementation of the TPC tracker
@@ -88,18 +21,19 @@ Splitted from AliTPCtracking
 //   Origin: Iouri Belikov, CERN, Jouri.Belikov@cern.ch 
 //-------------------------------------------------------
 #include <TObjArray.h>
+#include <TError.h>
 #include <TFile.h>
 #include <TTree.h>
-#include <Riostream.h>
+
+#include "AliESD.h"
 
 #include "AliTPCtracker.h"
 #include "AliTPCcluster.h"
 #include "AliTPCParam.h"
 #include "AliClusters.h"
 
-#include "TVector2.h"
+ClassImp(AliTPCtracker)
 
-#include <stdlib.h>
 //_____________________________________________________________________________
 AliTPCtracker::AliTPCtracker(const AliTPCParam *par): 
 AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
@@ -116,6 +50,7 @@ AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
 
   fParam = (AliTPCParam*) par;
   fSeeds=0;
+
 }
 
 //_____________________________________________________________________________
@@ -131,48 +66,6 @@ AliTPCtracker::~AliTPCtracker() {
   }
 }
 
-//_____________________________________________________________________________
-Double_t SigmaY2(Double_t r, Double_t tgl, Double_t pt)
-{
-  //
-  // Parametrised error of the cluster reconstruction (pad direction)   
-  //
-  // Sigma rphi
-  const Float_t kArphi=0.41818e-2;
-  const Float_t kBrphi=0.17460e-4;
-  const Float_t kCrphi=0.30993e-2;
-  const Float_t kDrphi=0.41061e-3;
-  
-  pt=TMath::Abs(pt)*1000.;
-  Double_t x=r/pt;
-  tgl=TMath::Abs(tgl);
-  Double_t s=kArphi - kBrphi*r*tgl + kCrphi*x*x + kDrphi*x;
-  if (s<0.4e-3) s=0.4e-3;
-  s*=1.3; //Iouri Belikov
-
-  return s;
-}
-
-//_____________________________________________________________________________
-Double_t SigmaZ2(Double_t r, Double_t tgl) 
-{
-  //
-  // Parametrised error of the cluster reconstruction (drift direction)
-  //
-  // Sigma z
-  const Float_t kAz=0.39614e-2;
-  const Float_t kBz=0.22443e-4;
-  const Float_t kCz=0.51504e-1;
-  
-
-  tgl=TMath::Abs(tgl);
-  Double_t s=kAz - kBz*r*tgl + kCz*tgl*tgl;
-  if (s<0.4e-3) s=0.4e-3;
-  s*=1.3; //Iouri Belikov
-
-  return s;
-}
-
 //_____________________________________________________________________________
 Double_t f1(Double_t x1,Double_t y1,
                    Double_t x2,Double_t y2,
@@ -224,31 +117,16 @@ Double_t f3(Double_t x1,Double_t y1,
 }
 
 //_____________________________________________________________________________
-void AliTPCtracker::LoadClusters() {
+Int_t AliTPCtracker::LoadClusters(TTree *cTree) {
   //-----------------------------------------------------------------
   // This function loads TPC clusters.
   //-----------------------------------------------------------------
-  if (!gFile->IsOpen()) {
-    cerr<<"AliTPCtracker::LoadClusters : "<<
-      "file with clusters has not been open !\n";
-    return;
-  }
-
-  Char_t name[100];
-  sprintf(name,"TreeC_TPC_%d",GetEventNumber());
-  TTree *cTree=(TTree*)gFile->Get(name);
-  if (!cTree) {
-    cerr<<"AliTPCtracker::LoadClusters : "<<
-      "can't get the tree with TPC clusters !\n";
-    return;
-  }
-
   TBranch *branch=cTree->GetBranch("Segment");
   if (!branch) {
-    cerr<<"AliTPCtracker::LoadClusters : "<<
-      "can't get the segment branch !\n";
-    return;
+     Error("LoadClusters","Can't get the branch !");
+     return 1;
   }
+
   AliClusters carray, *addr=&carray;
   carray.SetClass("AliTPCcluster");
   carray.SetArray(0);
@@ -264,10 +142,8 @@ void AliTPCtracker::LoadClusters() {
       Int_t nir=fInnerSec->GetNRows(), nor=fOuterSec->GetNRows();
       Int_t id=carray.GetID();
       if ((id<0) || (id>2*(fkNIS*nir + fkNOS*nor))) {
-         cerr<<"AliTPCtracker::LoadClusters : "<<
-              "wrong index !\n";
-         exit(1);
-      }        
+        Fatal("LoadClusters","Wrong index !");
+      }
       Int_t outindex = 2*fkNIS*nir;
       if (id<outindex) {
          Int_t sec = id/nir;
@@ -277,7 +153,7 @@ void AliTPCtracker::LoadClusters() {
          while (ncl--) {
            AliTPCcluster *c=(AliTPCcluster*)carray[ncl];
            padrow.InsertCluster(c,sec,row);
-         }           
+         }
       } else {
          id -= outindex;
          Int_t sec = id/nor;
@@ -289,10 +165,10 @@ void AliTPCtracker::LoadClusters() {
            padrow.InsertCluster(c,sec+fkNIS,row);
          }
       }
-
       carray.GetArray()->Clear();
   }
-  delete cTree;
+
+  return 0;
 }
 
 //_____________________________________________________________________________
@@ -336,14 +212,13 @@ Int_t AliTPCtracker::FollowProlongation(AliTPCseed& t, Int_t rf) {
     Double_t maxchi2=kMaxCHI2;
     const AliTPCRow &krow=fSectors[s][nr];
     Double_t pt=t.GetConvConst()/(100/0.299792458/0.2)/t.Get1Pt();
-    Double_t sy2=SigmaY2(t.GetX(),t.GetTgl(),pt);
-    Double_t sz2=SigmaZ2(t.GetX(),t.GetTgl());
+    Double_t sy2=AliTPCcluster::SigmaY2(t.GetX(),t.GetTgl(),pt);
+    Double_t sz2=AliTPCcluster::SigmaZ2(t.GetX(),t.GetTgl());
     Double_t road=4.*sqrt(t.GetSigmaY2() + sy2), y=t.GetY(), z=t.GetZ();
 
     if (road>kMaxROAD) {
       if (t.GetNumberOfClusters()>4) 
-        cerr<<t.GetNumberOfClusters()
-        <<"FindProlongation warning: Too broad road !\n"; 
+         Warning("FindProlongation","Too broad road !"); 
       return 0;
     }
 
@@ -387,85 +262,50 @@ Int_t AliTPCtracker::FollowProlongation(AliTPCseed& t, Int_t rf) {
 Int_t AliTPCtracker::FollowRefitInward(AliTPCseed *seed, AliTPCtrack *track) {
   //
   // This function propagates seed inward TPC using old clusters
-  // from track.
+  // from the track.
   // 
   // Sylwester Radomski, GSI
   // 26.02.2003
   //
 
-  Double_t alpha=seed->GetAlpha() - fSectors->GetAlphaShift();
-  TVector2::Phi_0_2pi(alpha);
-  Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
-
-  Int_t idx=-1, sec=-1, row=-1;
-  Int_t nc = seed->GetLabel(); // index to start with
-
-  idx=track->GetClusterIndex(nc);
-  sec=(idx&0xff000000)>>24; 
-  row=(idx&0x00ff0000)>>16;
-  
-  if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; } 
-  else { if (sec <  fkNIS) row=-1; }   
+  // loop over rows
 
-  Int_t nr=fSectors->GetNRows();
-  for (Int_t i=0; i<nr; i++) {
+  Int_t nRows = fSectors->GetNRows();
+  for (Int_t iRow = nRows-1; iRow >= 0; iRow--) {
 
-    Double_t x=fSectors->GetX(i); 
+    Double_t x = fSectors->GetX(iRow);
     if (!seed->PropagateTo(x)) return 0;
 
-    // Change sector and rotate seed if necessary
-
-    Double_t ymax=fSectors->GetMaxY(i);
-    Double_t y=seed->GetY();
-
-    if (y > ymax) {
-       s = (s+1) % fN;
-       if (!seed->Rotate(fSectors->GetAlpha())) return 0;
-    } else if (y <-ymax) {
-       s = (s-1+fN) % fN;
-       if (!seed->Rotate(-fSectors->GetAlpha())) return 0;
+    // try to find an assigned cluster in this row
+
+    AliTPCcluster* cluster = NULL;
+    Int_t idx = -1;
+    Int_t sec = -1;
+    for (Int_t iCluster = 0; iCluster < track->GetNumberOfClusters(); iCluster++) {
+      idx = track->GetClusterIndex(iCluster); 
+      sec = (idx&0xff000000)>>24; 
+      Int_t row = (idx&0x00ff0000)>>16;
+      if (((fSectors == fInnerSec) && (sec >= fkNIS)) ||
+         ((fSectors == fOuterSec) && (sec < fkNIS))) continue;
+      if (row == iRow) {
+       cluster = (AliTPCcluster*) GetCluster(idx);
+       break;
+      }
     }
 
-    // try to find a cluster
-    
-    AliTPCcluster *cl=0;
-    Int_t index = 0;
-    Double_t maxchi2 = kMaxCHI2;
-    
-    //cout << i << " " << row << " " << nc << endl;
-
-    if (row==i) {
+    // update the track seed with the found cluster
 
-      // accept already found cluster
-      AliTPCcluster *c=(AliTPCcluster*)GetCluster(idx);
-      Double_t chi2 = seed->GetPredictedChi2(c);
-      
-      index=idx; 
-      cl=c; 
-      maxchi2=chi2;
-      
-      if (++nc < track->GetNumberOfClusters() ) {
-        idx=track->GetClusterIndex(nc); 
-        sec=(idx&0xff000000)>>24; 
-        row=(idx&0x00ff0000)>>16;
+    if (cluster) {
+      Double_t dAlpha = fParam->GetAngle(sec) - seed->GetAlpha();
+      if (TMath::Abs(dAlpha) > 0.0001) {
+       if (!seed->Rotate(dAlpha)) return 0;
+       if (!seed->PropagateTo(x)) return 0;
       }
-            
-      if (fSectors==fInnerSec) { if (sec >= fkNIS) row=-1; }
-      else { if (sec < fkNIS) row=-1; }   
-      
-    }
-    
-    if (cl) {
-      
-      //cout << "Assigned" << endl;
-      Float_t l=fSectors->GetPadPitchWidth();
-      Float_t corr=1.; if (i>63) corr=0.67; // new (third) pad response !
-      seed->SetSampledEdx(cl->GetQ()/l*corr,seed->GetNumberOfClusters());
-      seed->Update(cl,maxchi2,index);
+
+      seed->Update(cluster, seed->GetPredictedChi2(cluster), idx);
     }
   }
 
-  seed->SetLabel(nc);
   return 1;
 }
 
@@ -481,7 +321,8 @@ Int_t AliTPCtracker::FollowBackProlongation
   Int_t s=Int_t(alpha/fSectors->GetAlpha())%fN;
 
   Int_t idx=-1, sec=-1, row=-1;
-  Int_t nc=seed.GetLabel(); //index of the cluster to start with
+  Int_t nc=seed.GetNumber();
+
   if (nc--) {
      idx=track.GetClusterIndex(nc);
      sec=(idx&0xff000000)>>24; row=(idx&0x00ff0000)>>16;
@@ -492,10 +333,8 @@ Int_t AliTPCtracker::FollowBackProlongation
   Int_t nr=fSectors->GetNRows();
   for (Int_t i=0; i<nr; i++) {
     Double_t x=fSectors->GetX(i), ymax=fSectors->GetMaxY(i);
-
-    if (!seed.PropagateTo(x)) return 0;
-
-    Double_t y=seed.GetY();
+    Double_t y=seed.GetYat(x);
     if (y > ymax) {
        s = (s+1) % fN;
        if (!seed.Rotate(fSectors->GetAlpha())) return 0;
@@ -504,20 +343,20 @@ Int_t AliTPCtracker::FollowBackProlongation
        if (!seed.Rotate(-fSectors->GetAlpha())) return 0;
     }
 
+    if (!seed.PropagateTo(x)) return 0;
+
     AliTPCcluster *cl=0;
     Int_t index=0;
     Double_t maxchi2=kMaxCHI2;
     Double_t pt=seed.GetConvConst()/(100/0.299792458/0.2)/seed.Get1Pt();
-    Double_t sy2=SigmaY2(seed.GetX(),seed.GetTgl(),pt);
-    Double_t sz2=SigmaZ2(seed.GetX(),seed.GetTgl());
+    Double_t sy2=AliTPCcluster::SigmaY2(seed.GetX(),seed.GetTgl(),pt);
+    Double_t sz2=AliTPCcluster::SigmaZ2(seed.GetX(),seed.GetTgl());
     Double_t road=4.*sqrt(seed.GetSigmaY2() + sy2), z=seed.GetZ();
     if (road>kMaxROAD) {
-      cerr<<seed.GetNumberOfClusters()
-          <<"AliTPCtracker::FollowBackProlongation: Too broad road !\n"; 
+      Warning("FollowBackProlongation","Too broad road !"); 
       return 0;
     }
 
-
     Int_t accepted=seed.GetNumberOfClusters();
     if (row==i) {
        //try to accept already found cluster
@@ -563,7 +402,7 @@ Int_t AliTPCtracker::FollowBackProlongation
 
   }
 
-  seed.SetLabel(nc);
+  seed.SetNumber(nc);
 
   return 1;
 }
@@ -616,8 +455,10 @@ void AliTPCtracker::MakeSeeds(Int_t i1, Int_t i2) {
         if (TMath::Abs(zz-z2)>5.) continue;
 
         Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
-        if (d==0.) {cerr<<"MakeSeeds warning: Straight seed !\n"; continue;}
-
+        if (d==0.) {
+           Warning("MakeSeeds","Straight seed !"); 
+           continue;
+        }
        x[0]=y1;
        x[1]=z1;
        x[4]=f1(x1,y1,x2,y2,x3,y3);
@@ -677,15 +518,14 @@ Int_t AliTPCtracker::ReadSeeds(const TFile *inp) {
 
   TFile *in=(TFile*)inp;
   if (!in->IsOpen()) {
-     cerr<<"AliTPCtracker::ReadSeeds(): input file is not open !\n";
+     Error("ReadSeeds","Input file has not been open !");
      return 1;
   }
 
   in->cd();
   TTree *seedTree=(TTree*)in->Get("Seeds");
   if (!seedTree) {
-     cerr<<"AliTPCtracker::ReadSeeds(): ";
-     cerr<<"can't get a tree with track seeds !\n";
+     Error("ReadSeeds","Can't get a tree with track seeds !");
      return 2;
   }
   AliTPCtrack *seed=new AliTPCtrack; 
@@ -696,7 +536,8 @@ Int_t AliTPCtracker::ReadSeeds(const TFile *inp) {
   Int_t n=(Int_t)seedTree->GetEntries();
   for (Int_t i=0; i<n; i++) {
      seedTree->GetEvent(i);
-     fSeeds->AddLast(new AliTPCseed(*seed,seed->GetAlpha()));
+     seed->ResetClusters();
+     fSeeds->AddLast(new AliTPCseed(*seed));
   }
   
   delete seed;
@@ -709,34 +550,11 @@ Int_t AliTPCtracker::ReadSeeds(const TFile *inp) {
 }
 
 //_____________________________________________________________________________
-Int_t AliTPCtracker::Clusters2Tracks(const TFile *inp, TFile *out) {
+Int_t AliTPCtracker::Clusters2Tracks(AliESD *event) {
   //-----------------------------------------------------------------
   // This is a track finder.
+  // The clusters must be already loaded ! 
   //-----------------------------------------------------------------
-  TDirectory *savedir=gDirectory; 
-
-  if (inp) {
-     TFile *in=(TFile*)inp;
-     if (!in->IsOpen()) {
-        cerr<<"AliTPCtracker::Clusters2Tracks(): input file is not open !\n";
-        return 1;
-     }
-  }
-
-  if (!out->IsOpen()) {
-     cerr<<"AliTPCtracker::Clusters2Tracks(): output file is not open !\n";
-     return 2;
-  }
-
-  LoadClusters();
-
-  out->cd();
-
-  char   tname[100];
-  sprintf(tname,"TreeT_TPC_%d",GetEventNumber());
-  TTree tracktree(tname,"Tree with TPC tracks");
-  AliTPCtrack *iotrack=0;
-  tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,3);
 
   //find track seeds
   Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
@@ -744,13 +562,12 @@ Int_t AliTPCtracker::Clusters2Tracks(const TFile *inp, TFile *out) {
   if (fSeeds==0) {
      Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
      fSectors=fOuterSec; fN=fkNOS;
-     fSeeds=new TObjArray(15000);     
+     fSeeds=new TObjArray(15000);
      MakeSeeds(nup-1, nup-1-gap);
      MakeSeeds(nup-1-shift, nup-1-shift-gap);
   }
   fSeeds->Sort();
 
-  Int_t found=0;
   Int_t nseed=fSeeds->GetEntriesFast();
   for (Int_t i=0; i<nseed; i++) {
     //tracking in the outer sectors
@@ -777,322 +594,149 @@ Int_t AliTPCtracker::Clusters2Tracks(const TFile *inp, TFile *out) {
         if (t.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
           t.CookdEdx();
           CookLabel(pt,0.1); //For comparison only
-          iotrack=pt;
-          tracktree.Fill();
+          pt->PropagateTo(fParam->GetInnerRadiusLow());
+          AliESDtrack iotrack;
+          iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
+
+          event->AddTrack(&iotrack);
+
           UseClusters(&t);
-          cerr<<found++<<'\r';
         }
       }
     }
-    delete fSeeds->RemoveAt(i); 
+    delete fSeeds->RemoveAt(i);
   }
-  
-  tracktree.Write();
-
-  cerr<<"Number of found tracks : "<<found<<endl;
 
-  savedir->cd();
+  Info("Clusters2Tracks","Number of found tracks : %d",
+       event->GetNumberOfTracks());
 
-  UnloadClusters();
   fSeeds->Clear(); delete fSeeds; fSeeds=0;
 
   return 0;
 }
-//_____________________________________________________________________________
 
-Int_t AliTPCtracker::RefitInward(TFile *in, TFile *out) {
+//_____________________________________________________________________________
+Int_t AliTPCtracker::RefitInward(AliESD* event) {
   //
   // The function propagates tracks throught TPC inward
   // using already associated clusters.
-  //
-  // in - file with back propagated tracks
-  // outs  - file with inward propagation
-  //
-  // Sylwester Radomski, GSI
-  // 26.02.2003
+  // The clusters must be already loaded !
   //
 
+  Int_t nTracks = event->GetNumberOfTracks();
+  Int_t nRefited = 0;
 
-  if (!in->IsOpen()) {
-    cout << "Input File not open !\n" << endl;
-    return 1;
-  }
-  
-  if (!out->IsOpen()) {
-    cout << "Output File not open !\n" << endl;
-    return 2;
-  }
-
-  TDirectory *savedir = gDirectory; 
-  LoadClusters();
+  for (Int_t i = 0; i < nTracks; i++) {
+    AliESDtrack* track = event->GetTrack(i);
+    ULong_t status = track->GetStatus();
 
-  // Connect to input tree
-  in->cd();
-  TTree *inputTree = (TTree*)in->Get("seedsTPCtoTRD_0");
-  TBranch *inBranch = inputTree->GetBranch("tracks");
-  AliTPCtrack *inTrack = 0;
-  inBranch->SetAddress(&inTrack);
+    if ( (status & AliESDtrack::kTPCrefit) != 0 ) continue;    
+    if ( (status & AliESDtrack::kTPCout ) == 0 ) continue;
 
-  Int_t nTracks = (Int_t)inputTree->GetEntries();
+    if ( (status & AliESDtrack::kTRDout ) != 0 ) 
+      if ( (status & AliESDtrack::kTRDrefit ) == 0 ) continue;
 
-  // Create output tree
-  out->cd(); 
-  AliTPCtrack *outTrack = new AliTPCtrack();
-  TTree *outputTree = new TTree("tracksTPC_0","Refited TPC tracks");
-  outputTree->Branch("tracks", "AliTPCtrack", &outTrack, 32000, 3);
+    AliTPCtrack* tpcTrack = new AliTPCtrack(*track);
+    AliTPCseed* seed=new AliTPCseed(*tpcTrack); seed->ResetClusters(); 
 
-  Int_t nRefited = 0;
+    if ( (status & AliESDtrack::kTRDrefit) == 0 ) seed->ResetCovariance();
 
-  for(Int_t t=0; t < nTracks; t++) {
-    
-    cout << t << "\r";
-
-    inputTree->GetEvent(t);
-    AliTPCseed *seed = new AliTPCseed(*inTrack, inTrack->GetAlpha());
-
-    seed->ResetCovariance();
+    fSectors = fOuterSec;
 
-    seed->SetLabel(0);    
-    fSectors = fInnerSec;
-    Int_t res = FollowRefitInward(seed, inTrack);
+    Int_t res = FollowRefitInward(seed, tpcTrack);
     UseClusters(seed);
     Int_t nc = seed->GetNumberOfClusters();
 
-    fSectors = fOuterSec;
-    res = FollowRefitInward(seed, inTrack);
+    fSectors = fInnerSec;
+
+    res = FollowRefitInward(seed, tpcTrack);
     UseClusters(seed, nc);
 
     if (res) {
       seed->PropagateTo(fParam->GetInnerRadiusLow());
-      outTrack = (AliTPCtrack*)seed;
-      outTrack->SetLabel(inTrack->GetLabel());
-      outputTree->Fill();
+      seed->SetLabel(tpcTrack->GetLabel());
+      seed->SetdEdx(tpcTrack->GetdEdx());
+      track->UpdateTrackParams(seed, AliESDtrack::kTPCrefit);
       nRefited++;
     }
 
     delete seed;
+    delete tpcTrack;
   }
 
-  cout << "Refitted " << nRefited << " tracks." << endl;
-
-  outputTree->Write();
-  
-  if (inputTree) delete inputTree;
-  if (outputTree) delete outputTree;
-
-  savedir->cd();
-  UnloadClusters();
+  Info("RefitInward","Number of refitted tracks : %d",nRefited);
 
   return 0;
 }
 
-//_____________________________________________________________________________
-Int_t AliTPCtracker::PropagateBack(const TFile *inp, TFile *out) {
+Int_t AliTPCtracker::PropagateBack(AliESD *event) {
   //-----------------------------------------------------------------
   // This function propagates tracks back through the TPC.
+  // The clusters must be already loaded !
   //-----------------------------------------------------------------
-  fSeeds=new TObjArray(15000);
+  Int_t nentr=event->GetNumberOfTracks();
+  Info("PropagateBack", "Number of ESD tracks: %d\n", nentr);
 
-  TFile *in=(TFile*)inp;
-  TDirectory *savedir=gDirectory; 
+  Int_t ntrk=0;
+  for (Int_t i=0; i<nentr; i++) {
+    AliESDtrack *esd=event->GetTrack(i);
+    ULong_t status=esd->GetStatus();
 
-  if (!in->IsOpen()) {
-     cerr<<"AliTPCtracker::PropagateBack(): ";
-     cerr<<"file with back propagated ITS tracks is not open !\n";
-     return 1;
-  }
+    if ( (status & AliESDtrack::kTPCin ) == 0 ) continue;
+    if ( (status & AliESDtrack::kTPCout) != 0 ) continue;
+    if ( (status & AliESDtrack::kITSin) != 0 )
+       if ( (status & AliESDtrack::kITSout) == 0 ) continue;
 
-  if (!out->IsOpen()) {
-     cerr<<"AliTPCtracker::PropagateBack(): ";
-     cerr<<"file for back propagated TPC tracks is not open !\n";
-     return 2;
-  }
+    const AliTPCtrack t(*esd);
+    AliTPCseed s(t); s.ResetClusters();
 
-  LoadClusters();
+    if ( (status & AliESDtrack::kITSout) == 0 ) s.ResetCovariance();
 
-  in->cd();
-  char tName[100];
-  sprintf(tName,"TreeT_ITSb_%d",GetEventNumber());
-  TTree *bckTree=(TTree*)in->Get(tName);
-  if (!bckTree) {
-     cerr<<"AliTPCtracker::PropagateBack() ";
-     cerr<<"can't get a tree with back propagated ITS tracks !\n";
-     return 3;
-  }
-  AliTPCtrack *bckTrack=new AliTPCtrack; 
-  bckTree->SetBranchAddress("tracks",&bckTrack);
-
-  sprintf(tName,"TreeT_TPC_%d",GetEventNumber());
-  TTree *tpcTree=(TTree*)in->Get(tName);
-  if (!tpcTree) {
-     cerr<<"AliTPCtracker::PropagateBack() ";
-     cerr<<"can't get a tree with TPC tracks !\n";
-     return 4;
-  }
-  AliTPCtrack *tpcTrack=new AliTPCtrack; 
-  tpcTree->SetBranchAddress("tracks",&tpcTrack);
+    Int_t nc=t.GetNumberOfClusters();
+    s.SetNumber(nc); //set number of the cluster to start with
 
-  //*** Prepare an array of tracks to be back propagated
-  Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
-  Int_t nrows=nlow+nup;
+    //inner sectors
+    fSectors=fInnerSec; fN=fkNIS;
 
-  //
-  // Match ITS tracks with old TPC tracks
-  // the tracks do not have to be sorted [SR, GSI, 18.02.2003]
-  //
-  // the algorithm is linear and uses LUT for sorting
-  // cost of the algorithm: 2 * number of tracks
-  //
+    Double_t alpha=s.GetAlpha() - fSectors->GetAlphaShift();
+    if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
+    if (alpha < 0.            ) alpha += 2.*TMath::Pi();
+    Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
+    alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
+    alpha-=s.GetAlpha();
 
-  TObjArray tracks(15000);
-  Int_t i=0;
-  Int_t tpcN= (Int_t)tpcTree->GetEntries();
-  Int_t bckN= (bckTree)? (Int_t)bckTree->GetEntries() : 0;
+    if (!s.Rotate(alpha)) continue;
+    if (!FollowBackProlongation(s,t)) continue;
 
-  // look up table   
-  const Int_t nLab = 20000; // limit on number of primaries (arbitrary)
-  Int_t lut[nLab];
-  for(Int_t i=0; i<nLab; i++) lut[i] = -1;
-    
-  if (bckTree) {
-    for(Int_t i=0; i<bckN; i++) {
-      bckTree->GetEvent(i);
-      Int_t lab = TMath::Abs(bckTrack->GetLabel());
-      if (lab < nLab) lut[lab] = i;
-    }
-  }
-  
-  for (Int_t i=0; i<tpcN; i++) {
-    tpcTree->GetEvent(i);
-    Double_t alpha=tpcTrack->GetAlpha();
-    
-    if (!bckTree) { 
+    UseClusters(&s);
 
-      // No ITS - use TPC track only
-      
-      fSeeds->AddLast(new AliTPCseed(*tpcTrack, tpcTrack->GetAlpha()));
-      tracks.AddLast(new AliTPCtrack(*tpcTrack));
-    
-    } else {  
-      
-      // with ITS
-      // discard not prolongated tracks (to be discussed)
+    //outer sectors
+    fSectors=fOuterSec; fN=fkNOS;
 
-      Int_t lab = TMath::Abs(tpcTrack->GetLabel());
-      if (lab > nLab || lut[lab] < 0) continue;
+    nc=s.GetNumberOfClusters();
 
-      bckTree->GetEvent(lut[lab]);   
-      bckTrack->Rotate(alpha-bckTrack->GetAlpha());
+    alpha=s.GetAlpha() - fSectors->GetAlphaShift();
+    if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
+    if (alpha < 0.            ) alpha += 2.*TMath::Pi();
+    ns=Int_t(alpha/fSectors->GetAlpha())%fN;
 
-      fSeeds->AddLast(new AliTPCseed(*bckTrack,bckTrack->GetAlpha()));
-      tracks.AddLast(new AliTPCtrack(*tpcTrack));
-    }
-  }
+    alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
+    alpha-=s.GetAlpha();
 
-  // end of matching
-
-  /*
-  TObjArray tracks(15000);
-  Int_t i=0, j=0;
-  Int_t tpcN=(Int_t)tpcTree->GetEntries();
-  Int_t bckN=(Int_t)bckTree->GetEntries();
-  if (j<bckN) bckTree->GetEvent(j++);
-  for (i=0; i<tpcN; i++) {
-     tpcTree->GetEvent(i);
-     Double_t alpha=tpcTrack->GetAlpha();
-     if (j<bckN &&
-     TMath::Abs(tpcTrack->GetLabel())==TMath::Abs(bckTrack->GetLabel())) {
-        if (!bckTrack->Rotate(alpha-bckTrack->GetAlpha())) continue;
-        fSeeds->AddLast(new AliTPCseed(*bckTrack,bckTrack->GetAlpha()));
-        bckTree->GetEvent(j++);
-     } else {
-        tpcTrack->ResetCovariance();
-        fSeeds->AddLast(new AliTPCseed(*tpcTrack,alpha));
-     }
-     tracks.AddLast(new AliTPCtrack(*tpcTrack));
+    if (!s.Rotate(alpha)) continue;
+    if (!FollowBackProlongation(s,t)) continue;
+    {
+    Int_t nrows=fOuterSec->GetNRows()+fInnerSec->GetNRows();
+    if (s.GetNumberOfClusters() < Int_t(0.4*nrows)) continue;
+    }
+    s.PropagateTo(fParam->GetOuterRadiusUp());
+    s.CookdEdx();
+    CookLabel(&s,0.1); //For comparison only
+    UseClusters(&s,nc);
+    esd->UpdateTrackParams(&s,AliESDtrack::kTPCout);
+    ntrk++;
   }
-  */
-  
-  out->cd();
-
-  // tree name seedsTPCtoTRD as expected by TRD and as 
-  // discussed and decided in Strasbourg (May 2002)
-  // [SR, GSI, 18.02.2003]
-  
-  sprintf(tName,"seedsTPCtoTRD_%d",GetEventNumber());
-  TTree backTree(tName,"Tree with back propagated TPC tracks");
-
-  AliTPCtrack *otrack=0;
-  backTree.Branch("tracks","AliTPCtrack",&otrack,32000,0);
-
-  //*** Back propagation through inner sectors
-  fSectors=fInnerSec; fN=fkNIS;
-  Int_t nseed=fSeeds->GetEntriesFast();
-  for (i=0; i<nseed; i++) {
-     AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
-     const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
-
-     Int_t nc=t.GetNumberOfClusters();
-     s.SetLabel(nc-1); //set number of the cluster to start with
-
-     if (FollowBackProlongation(s,t)) {
-       UseClusters(&s);
-        continue;
-     }
-     delete fSeeds->RemoveAt(i);
-  }  
-
-//*** Back propagation through outer sectors
-  fSectors=fOuterSec; fN=fkNOS;
-  Int_t found=0;
-  for (i=0; i<nseed; i++) {
-     AliTPCseed *ps=(AliTPCseed*)fSeeds->UncheckedAt(i), &s=*ps;
-     if (!ps) continue;
-     Int_t nc=s.GetNumberOfClusters();
-     const AliTPCtrack *pt=(AliTPCtrack*)tracks.UncheckedAt(i), &t=*pt;
-
-     Double_t alpha=s.GetAlpha() - fSectors->GetAlphaShift();
-     if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();
-     if (alpha < 0.            ) alpha += 2.*TMath::Pi();
-     Int_t ns=Int_t(alpha/fSectors->GetAlpha())%fN;
-
-     alpha =ns*fSectors->GetAlpha() + fSectors->GetAlphaShift();
-     alpha-=s.GetAlpha();
-
-     if (s.Rotate(alpha)) {
-       if (FollowBackProlongation(s,t)) {
-         if (s.GetNumberOfClusters() >= Int_t(0.4*nrows)) {
-           s.CookdEdx();
-           s.SetLabel(t.GetLabel());
-           UseClusters(&s,nc);
-           
-           // Propagate to outer reference plane for comparison reasons
-           // reason for keeping fParam object [SR, GSI, 18.02.2003] 
-         
-           ps->PropagateTo(fParam->GetOuterRadiusUp()); 
-           otrack=ps;
-           backTree.Fill();
-           cerr<<found++<<'\r';
-           continue;
-         }
-       }
-     }
-     delete fSeeds->RemoveAt(i);
-  }  
-
-  backTree.Write();
-  savedir->cd();
-  cerr<<"Number of seeds: "<<nseed<<endl;
-  cerr<<"Number of back propagated ITS tracks: "<<bckN<<endl;
-  cerr<<"Number of back propagated TPC tracks: "<<found<<endl;
-
-  delete bckTrack;
-  delete tpcTrack;
-
-  delete bckTree; //Thanks to Mariana Bondila
-  delete tpcTree; //Thanks to Mariana Bondila
-
-  UnloadClusters();  
+  Info("PropagateBack","Number of back propagated tracks: %d",ntrk);
 
   return 0;
 }
@@ -1212,12 +856,15 @@ AliTPCRow::InsertCluster(const AliTPCcluster* c, Int_t sec, Int_t row) {
   // Insert a cluster into this pad row in accordence with its y-coordinate
   //-----------------------------------------------------------------------
   if (fN==kMaxClusterPerRow) {
-    cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
+     ::Error("InsertCluster","Too many clusters !"); 
+     return;
   }
 
   Int_t index=(((sec<<8)+row)<<16)+fN;
 
   if (fN==0) {
+     fSize=kMaxClusterPerRow/8;
+     fClusterArray=new AliTPCcluster[fSize];
      fIndex[0]=index;
      fClusterArray[0]=*c; fClusters[fN++]=fClusterArray; 
      return;
@@ -1282,23 +929,7 @@ void AliTPCtracker::AliTPCseed::CookdEdx(Double_t low, Double_t up) {
   dedx /= (nu-nl+1);
   SetdEdx(dedx);
 
-  //Very rough PID
-  Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
-
-  Double_t log1=TMath::Log(p+0.45), log2=TMath::Log(p+0.12);
-  if (p<0.6) {
-    if (dedx < 34 + 30/(p+0.45)/(p+0.45) + 24*log1) {SetMass(0.13957); return;}
-    if (dedx < 34 + 30/(p+0.12)/(p+0.12) + 24*log2) {SetMass(0.49368); return;}
-    SetMass(0.93827); return;
-  }
-
-  if (p<1.2) {
-    if (dedx < 34 + 30/(p+0.12)/(p+0.12) + 24*log2) {SetMass(0.13957); return;}
-    SetMass(0.93827); return;
-  }
-
-  SetMass(0.13957); return;
-
+  return;
 }