First version of the parallel TPC tracking (M.Ivanov)
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 15 Nov 2002 14:27:45 +0000 (14:27 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 15 Nov 2002 14:27:45 +0000 (14:27 +0000)
23 files changed:
TPC/AliComplexCluster.cxx
TPC/AliComplexCluster.h
TPC/AliH2F.cxx
TPC/AliTPCCompareClusters.C
TPC/AliTPCCompareTracks.C [new file with mode: 0644]
TPC/AliTPCComparison.C
TPC/AliTPCComparison2.C
TPC/AliTPCFindClusters.C
TPC/AliTPCFindClustersMI.C [new file with mode: 0644]
TPC/AliTPCFindTracks.C
TPC/AliTPCFindTracksMI.C [new file with mode: 0644]
TPC/AliTPCHits2ExactClusters.C
TPC/AliTPCclusterMI.cxx [new file with mode: 0644]
TPC/AliTPCclusterMI.h [new file with mode: 0644]
TPC/AliTPCclustererMI.cxx [new file with mode: 0644]
TPC/AliTPCclustererMI.h [new file with mode: 0644]
TPC/AliTPCtrack.h
TPC/AliTPCtrackerMI.cxx [new file with mode: 0644]
TPC/AliTPCtrackerMI.h [new file with mode: 0644]
TPC/LandauTest.C [new file with mode: 0644]
TPC/TPCLinkDef.h
TPC/alles.h
TPC/libTPC.pkg

index 29d13e5..332a0c0 100644 (file)
@@ -15,6 +15,9 @@
 
 /*
 $Log$
+Revision 1.2  2001/02/05 14:43:13  hristov
+Compare() declared const
+
 Revision 1.1  2000/10/05 16:17:27  kowal2
 New class replacing AliCluster
 
@@ -58,4 +61,10 @@ Bool_t AliComplexCluster::IsSortable() const
 }
 
 ClassImp(AliDigitCluster)
-ClassImp(AliDifCluster)
+ClassImp(AliTPCExactPoint)
+ClassImp(AliTPCClusterPoint)
+ClassImp(AliTPCTrackerPoint)
+ClassImp(AliTPCTrackPoint)
+ClassImp(AliTPCTrackPointRef)
+
+
index 20d5f57..a6bda58 100644 (file)
@@ -6,6 +6,7 @@
 /* $Id$ */
 
 #include "TObject.h"
+#include "TMath.h"
 
 class AliComplexCluster : public TObject {
 public:
@@ -21,7 +22,7 @@ public:
 public:
   AliComplexCluster() {
     fTracks[0]=fTracks[1]=fTracks[2]=0; 
-    fX=fY=fQ=fSigmaX2=fSigmaY2=0.;
+    fX=fY=fQ=fSigmaX2=fSigmaY2=fSigmaXY=fArea=fMax=0.;
   }
   virtual ~AliComplexCluster() {;}
   Bool_t    IsSortable() const;
@@ -39,15 +40,109 @@ public:
   ClassDef(AliDigitCluster,1)  // Tclusters
 };
 
-class AliDifCluster : public AliDigitCluster {
-public:
-  Float_t fDx; //delta x 
-  Float_t fDy; //delta y
-  Float_t fAngleX;//x angle
-  Float_t fAngleY;//y angle
-  Float_t fOrigQ; //original charge
-  Int_t fGenerTrack;  //track number of nearest track
-  ClassDef(AliDifCluster,1)  // // Cluster manager
+
+
+class AliTPCTrackerPoint  {  
+ private:
+  Short_t   fTX;        // x position of the cluster  in cm - 10 mum prec
+  Short_t   fTZ;        // current prolongation in Z  in cm - 10 mum prec.
+  Short_t   fTY;        // current prolongation in Y  in cm - 10 mum prec.
+  Char_t   fTAngleZ;   // angle 
+  Char_t   fTAngleY;   // angle
+ public:
+  AliTPCTrackerPoint(){fTX=0; fTY=0; fTZ=0; fTAngleZ=0; fTAngleY=0;}
+  inline Float_t  GetX() {return (fTX*0.01);}
+  inline Float_t  GetZ() {return (fTZ*0.01);}
+  inline Float_t  GetY() {return (fTY*0.01);}
+  inline Float_t  GetAngleZ() {return (Float_t(fTAngleZ)*0.02);}
+  inline Float_t  GetAngleY() {return (Float_t(fTAngleY)*0.02);}
+  //
+  void     SetX(Float_t x){ fTX = Short_t(TMath::Nint(x*100.));} 
+  void     SetY(Float_t y){ fTY = Short_t(TMath::Nint(y*100.));} 
+  void     SetZ(Float_t z){ fTZ = Short_t(TMath::Nint(z*100.));} 
+  void     SetAngleZ(Float_t anglez) {fTAngleZ = Char_t(TMath::Nint(anglez*50.));}
+  void     SetAngleY(Float_t angley) {fTAngleY = Char_t(TMath::Nint(angley*50.));}
+  //
+ public:
+  ClassDef(AliTPCTrackerPoint,1)  
 };
 
+class AliTPCClusterPoint  {
+ private:
+  Short_t  fCZ;       // current cluster position Z in cm - 100 mum precision
+  Short_t  fCY;       // current cluster position Y in cm - 100 mum precision
+  UChar_t  fErrZ;     // z error estimate - in  mm - 50 mum precision 
+  UChar_t  fErrY;     // y error estimate - in  mm - 50 mum precision 
+  UChar_t  fSigmaZ;   // shape  Z - normalised shape - normaliziation 1 - precision 2 percent
+  UChar_t  fSigmaY;   // shape  Y - normalised shape - normaliziation 1 - precision 2 percent
+  UShort_t fQ;        // total charge in cluster 
+  UShort_t fMax;      // charge at maximum  
+  Char_t   fCType;    // type of the cluster
+ public:
+  AliTPCClusterPoint(){fCZ=fCY=fSigmaZ=fSigmaY=fErrZ=fErrY=fQ=fMax=fCType=0;}
+  inline Float_t  GetZ()    {return (fCZ*0.01);}
+  inline Float_t  GetY()    {return (fCY*0.01);}
+  inline Float_t  GetErrZ()   {return (fErrZ*0.005);}
+  inline Float_t  GetErrY()   {return (fErrY*0.005);}
+  inline Float_t  GetSigmaZ() {return (fSigmaZ*0.02);}
+  inline Float_t  GetSigmaY() {return (fSigmaY*0.02);}  
+  inline Int_t  GetType()   {return fCType;}
+  inline Int_t  GetMax()   {return fMax;}
+  inline Float_t  GetQ()   {return fQ;}
+
+  //
+  void     SetY(Float_t y){ fCY = Short_t(TMath::Nint(y*100.));} 
+  void     SetZ(Float_t z){ fCZ = Short_t(TMath::Nint(z*100.));} 
+  void     SetErrZ(Float_t errz) {fErrZ = UChar_t(TMath::Nint(errz*200.));}
+  void     SetErrY(Float_t erry) {fErrY = UChar_t(TMath::Nint(erry*200.));}
+  void     SetSigmaZ(Float_t sigmaz) {fSigmaZ = UChar_t(TMath::Nint(sigmaz*50.));}
+  void     SetSigmaY(Float_t sigmay) {fSigmaY = UChar_t(TMath::Nint(sigmay*50.));}
+  void     SetQ(Float_t q) {fQ = UShort_t(q);}
+  void     SetMax(Float_t max) {fMax = UShort_t(max);}
+  void     SetType(Char_t type) {fCType = type;}
+
+  //
+ public:
+  ClassDef(AliTPCClusterPoint,1)  
+};
+
+
+class AliTPCExactPoint : public TObject{
+ public:
+  AliTPCExactPoint(){fEZ=fEY=fEAngleZ=fEAngleY=fEAmp=fEPrim=fTrackID=0;}
+  Float_t fEZ;       // current "exact" position according simulation
+  Float_t fEY;       // current "exact" position according simulation
+  Float_t fEAngleZ;  // angle Z
+  Float_t fEAngleY;  // angle Y
+  Float_t fEAmp;     // total charge deposited in row
+  Float_t fEPrim;    // primary charge deposited in row
+  Int_t   fTrackID;  // id of the track
+  ClassDef(AliTPCExactPoint,1)  
+};
+
+
+class AliTPCTrackPoint: public TObject{
+ public:
+  AliTPCTrackPoint(){ fIsShared = kFALSE;}
+  AliTPCClusterPoint & GetCPoint(){return fCPoint;}
+  AliTPCTrackerPoint & GetTPoint(){return fTPoint;}
+ public:
+  AliTPCClusterPoint fCPoint;
+  AliTPCTrackerPoint fTPoint;
+  Bool_t fIsShared;  // indicate sharing of the point between several tracks
+  ClassDef(AliTPCTrackPoint,1)  
+};
+
+
+class AliTPCTrackPointRef: public AliTPCTrackPoint{
+ public:
+  AliTPCExactPoint & GetExactPoint(){return fEPoint;}
+  AliTPCExactPoint & GetNearestPoint(){return fNPoint;}  
+ public:
+  AliTPCExactPoint fEPoint; //exact point belonging to track
+  AliTPCExactPoint fNPoint; //nearest point  
+  ClassDef(AliTPCTrackPointRef,1)  
+};
+
+
 #endif //ALICOMPLEXCLUSTER_H
index b59a819..831f46f 100644 (file)
@@ -15,6 +15,9 @@
 
 /*
 $Log$
+Revision 1.4  2000/10/05 16:04:21  kowal2
+Forward declarations
+
 Revision 1.3  2000/06/30 12:07:49  kowal2
 Updated from the TPC-PreRelease branch
 
@@ -45,7 +48,6 @@ Revision 1.1.4.2  2000/04/10 11:32:37  kowal2
 #include "TClonesArray.h"
 #include "AliTPC.h"
 #include "TRandom.h"
-#include "AliTPCClusterFinder.h"
 
 
 ClassImp(AliH2F)
@@ -87,15 +89,11 @@ AliH2F & AliH2F::operator = (const AliH2F & his)
 
 TClonesArray * AliH2F::FindPeaks(Float_t threshold, Float_t noise)
 {
-  //find peaks and write it in form of AliTPCcluster to array
-    
-  //firstly we need to create object for cluster finding
-  //and fill it with contents of histogram
-  AliTPCClusterFinder cfinder;
-  cfinder.SetThreshold(threshold);
-  cfinder.SetNoise(noise);
-  cfinder.GetHisto(this);
-  return cfinder.FindPeaks3();
+
+  //
+  // not implemented
+  //
+  return 0;
 }
 
 void AliH2F::ClearSpectrum()
index dedf107..b6ccd68 100644 (file)
@@ -1,3 +1,9 @@
+/****************************************************************************
+ *           Origin: M.Ivanov, CERN, Marian.Ivanov@cern.ch                 *
+ ****************************************************************************/
+
+
+
 #ifndef __CINT__
 #include "alles.h"
 #include "AliComplexCluster.h"
diff --git a/TPC/AliTPCCompareTracks.C b/TPC/AliTPCCompareTracks.C
new file mode 100644 (file)
index 0000000..4adac14
--- /dev/null
@@ -0,0 +1,172 @@
+#ifndef __CINT__
+#include "alles.h"
+#include "AliComplexCluster.h"
+//#include "AliTPCclusterM.h"
+#include "AliTPCtrackerMI.h"
+#include "AliTPCclusterMI.h"
+#endif
+
+
+Int_t AliTPCCompareTracks(Int_t eventn, Bool_t all = kFALSE) {
+   cerr<<"Comparing tracks...\n";
+   //CONNECT FILES
+   TFile *file=TFile::Open("galice.root");
+   if (!file->IsOpen()) {cerr<<"Can't open galice.root !\n"; return 1;}
+   //
+   TFile *ftracks=TFile::Open("AliTPCtracks.root","update");
+   if (!ftracks->IsOpen()){cerr<<"Can't open AliTPCtracks.root !\n"; return 3;}
+   //
+   AliTPCParam *param=(AliTPCParam *)file->Get("75x40_100x60_150x60");
+   if (!param) {cerr<<"TPC parameters have not been found !\n"; return 2;}   
+   //
+   TFile *fout=TFile::Open("AliTPCTracksDif.root","new"); 
+   if (!fout->IsOpen()) fout=TFile::Open("AliTPCClustersDif.root","recreate");  
+   //
+   //connect exact clusters 
+   file->cd();
+   AliTPCClustersArray *ce=new AliTPCClustersArray;
+   ce->Setup(param);
+   ce->SetClusterType("AliComplexCluster");
+   char  cname[100];
+   sprintf(cname,"TreeCExact_%d",eventn);
+   ce->ConnectTree(cname);
+   //
+   //connect reconstructed tracks
+   ftracks->cd();
+   sprintf(cname,"Seeds");
+   TTree * treetracks = (TTree*)ftracks->Get(cname);
+   TBranch * branchtracks = treetracks->GetBranch("seeds");   
+   //
+   //load seeds to the memory
+   Int_t trackmap[500000][4];  // map of tracks corresponding to given track number
+   memset(trackmap,0,sizeof(Int_t)*4*500000);
+   Int_t ntracks = treetracks->GetEntries();
+   TObjArray * trackarr= new TObjArray(ntracks);
+   Int_t nproces = TMath::Min(ntracks,4000);
+
+   //   for (Int_t itrack =0; itrack<ntracks; itrack++){
+   for (Int_t itrack =0; itrack<nproces; itrack++){
+     AliTPCseed * seed = new AliTPCseed;  
+     //
+     seed->fPoints = new TClonesArray("AliTPCTrackPoint",200);
+     branchtracks->SetAddress(&seed); 
+     branchtracks->GetEntry(itrack);
+
+     //if (seed->fRemoval>0 && (itrack%4) ) continue;
+     trackarr->AddLast(seed);    
+
+     //crete array with exact position information
+     seed->fEPoints = new TClonesArray("AliTPCTrackPointRef",1);
+     seed->fEPoints->ExpandCreateFast(200);
+
+     Int_t label = TMath::Abs(seed->GetLabel());
+     Int_t i;
+     if (label>500000) {
+       //too big track label
+     }else{
+       for (i=0;i<4;i++) {
+        if ( trackmap[label][i]==0) 
+          break;
+       }
+       if(i<4) trackmap[label][i]=itrack;
+     }
+   }
+
+   //add information about exact positions   
+   Int_t nrows=Int_t(ce->GetTree()->GetEntries());
+
+   for (Int_t n=0; n<nrows; n++) {
+       AliSegmentID *se=ce->LoadEntry(n);
+       Int_t sec,row;
+       param->AdjustSectorRow(se->GetID(),sec,row);
+       //
+       AliTPCClustersRow* clrowe = ce->GetRow(sec,row);
+       //
+       Int_t ncl=clrowe->GetArray()->GetEntriesFast();
+       const Int_t kNIS=param->GetNInnerSector(), kNOS=param->GetNOuterSector();
+       Int_t npads, sign;
+       if (sec < kNIS) {
+          npads = param->GetNPadsLow(row);
+          sign = (sec < kNIS/2) ? 1 : -1;
+       } else {
+          npads = param->GetNPadsUp(row);
+          sign = ((sec-kNIS) < kNOS/2) ? 1 : -1;
+       }
+       
+       Int_t trrow = row;
+       if (sec>=param->GetNInnerSector()) trrow += param->GetNRowLow(); 
+       
+       while (ncl--) {
+          AliComplexCluster *cl=(AliComplexCluster*)(clrowe->GetArray()->At(ncl));
+           Double_t x=param->GetPadRowRadii(sec,row);
+          Double_t y=cl->fY;
+          y = ( y + 0.5 -  0.5*npads) *param->GetPadPitchWidth(sec);
+          Double_t z=cl->fX*param->GetZWidth();
+          z = sign*(param->GetZLength() - z);
+           Float_t cs, sn, tmp;
+           param->AdjustCosSin(sec,cs,sn);
+          for (Int_t i=0;i<3;i++){
+            if (cl->fTracks[0]<500000) if (trackmap[cl->fTracks[0]][i]) {
+              AliTPCseed * seed = (AliTPCseed*)trackarr->At(trackmap[cl->fTracks[0]][i]);
+              TClonesArray * clarr = seed->fPoints;
+              if (!clarr){
+                //printf("Seed %d without info",trackmap[cl->fTracks[0]][i]);
+                continue;
+              }
+              AliTPCTrackPoint    * trcluster = (AliTPCTrackPoint*)(seed->fPoints->At(trrow));
+              AliTPCTrackPointRef * ecluster  = (AliTPCTrackPointRef*)(seed->fEPoints->At(trrow));
+              //
+              ecluster->GetCPoint() = trcluster->GetCPoint();
+              ecluster->GetTPoint() = trcluster->GetTPoint();
+              //
+              AliTPCExactPoint & epoint =  ecluster->GetExactPoint();
+              /*
+                trcluster->fEZ = z;
+                trcluster->fEY = y;
+                trcluster->fAngleY = cl->fSigmaY2*param->GetPadPitchLength(sec);
+                trcluster->fAngleZ = cl->fSigmaX2*param->GetPadPitchLength(sec);
+                trcluster->fEAmp = cl->fQ;
+                trcluster->fEPrim = cl->fMax;
+              */
+              epoint.fEZ = z;
+              epoint.fEY = y;
+              epoint.fEAngleY = cl->fSigmaY2*param->GetPadPitchLength(sec);
+              epoint.fEAngleZ = cl->fSigmaX2*param->GetPadPitchLength(sec);
+              epoint.fEAmp = cl->fQ;
+              epoint.fEPrim = cl->fMax;
+            }
+          } 
+       }
+       //       cr->ClearRow(sec,row);
+       ce->ClearRow(sec,row);
+   }
+
+   // make new tree - with tracks - exact position
+   fout->cd();
+   TTree * treenew = new TTree("Seedref","Seedref");
+   AliTPCseed *  ioseed  =   (AliTPCseed*)trackarr->At(0);   
+   TBranch * br = treenew->Branch("seeds","AliTPCseed",&ioseed,32000,99);
+
+   for (Int_t itrack =0; itrack<ntracks; itrack++){
+     ioseed  =  (AliTPCseed*) trackarr->At(itrack); 
+     br->SetAddress(&ioseed);
+     treenew->Fill();
+   }
+
+   treenew->Write();
+   printf("1\n");
+   delete ce;
+   printf("2\n");
+   //delete cr;
+   //   carray->GetTree()->Write();
+   printf("3\n");
+   //   delete carray;
+   printf("4\n");
+   //   cf->Close(); 
+   printf("5\n");
+   fout->Close(); 
+   printf("6\n");
+   ftracks->Close();
+   printf("7\n");
+   return 0;
+}
index e938a88..2587a1a 100644 (file)
@@ -11,6 +11,7 @@
 #ifndef __CINT__
   #include "alles.h"
   #include "AliTPCtracker.h"
+  #include <AliMagF.h>
 #endif
 
 struct GoodTrackTPC {
@@ -23,7 +24,7 @@ struct GoodTrackTPC {
 Int_t AliTPCComparison(Int_t event=0) {
 
    cerr<<"Doing comparison...\n";
-
+   AliKalmanTrack::SetConvConst(1000/0.299792458/gAlice->Field()->SolenoidField());
    const Int_t MAX=20000;
    Int_t good_tracks_tpc(GoodTrackTPC *gt, const Int_t max, const Int_t event);
 
@@ -352,6 +353,8 @@ Int_t good_tracks_tpc(GoodTrackTPC *gt, const Int_t max, const Int_t event) {
      break;
    case 2:
      {
+       TFile fff("digits.root");
+        
       char dname[100]; sprintf(dname,"TreeD_75x40_100x60_150x60_%d",event);
       TTree *TD=(TTree*)gDirectory->Get(dname);
       AliSimDigits da, *digits=&da;
index c9dc7e6..8486a9b 100644 (file)
@@ -396,7 +396,9 @@ Int_t good_tracks(GoodTrackTPC *gt, Int_t max, Int_t firstev, Int_t eventn) {
     cerr<<"gAlice have not been found on galice.root !\n";
     exit(5);
   }
-   
+
+  TFile *fdigit = TFile::Open("digits.root");
+  file->cd();
 
   AliTPC *TPC=(AliTPC*)gAlice->GetDetector("TPC");
   Int_t ver = TPC->IsVersion(); 
@@ -474,7 +476,8 @@ Int_t good_tracks(GoodTrackTPC *gt, Int_t max, Int_t firstev, Int_t eventn) {
     case 2:
       {
         sprintf(treeName,"TreeD_75x40_100x60_150x60_%d",event);  
-        TD=(TTree*)gDirectory->Get(treeName);
+       //        TD=(TTree*)gDirectory->Get(treeName);
+        TD=(TTree*)fdigit->Get(treeName);
         TD->GetBranch("Segment")->SetAddress(&digits);
         count = new Int_t[np];
         Int_t i;
@@ -485,12 +488,15 @@ Int_t good_tracks(GoodTrackTPC *gt, Int_t max, Int_t firstev, Int_t eventn) {
           Int_t sec,row;
           digp->AdjustSectorRow(digits->GetID(),sec,row);
           digits->First();
+         digits->ExpandTrackBuffer();
           do { //Many thanks to J.Chudoba who noticed this
             Int_t it=digits->CurrentRow(), ip=digits->CurrentColumn();
-            Short_t dig = digits->GetDigit(it,ip);
-            Int_t idx0=digits->GetTrackID(it,ip,0); 
-            Int_t idx1=digits->GetTrackID(it,ip,1);
-            Int_t idx2=digits->GetTrackID(it,ip,2);
+           //            Short_t dig = digits->GetDigit(it,ip);
+            Short_t dig = digits->CurrentDigit();
+
+            Int_t idx0=digits->GetTrackIDFast(it,ip,0)-2; 
+            Int_t idx1=digits->GetTrackIDFast(it,ip,1)-2;
+            Int_t idx2=digits->GetTrackIDFast(it,ip,2)-2;
             if (idx0>=0 && dig>=zero && idx0<np ) count[idx0]+=1;   //background events - higher track ID's
             if (idx1>=0 && dig>=zero && idx1<np ) count[idx1]+=1;
             if (idx2>=0 && dig>=zero && idx2<np ) count[idx2]+=1;
@@ -611,6 +617,8 @@ Int_t FindPrimaries(Int_t nparticles){
   for(Int_t iprim = 0; iprim<nparticles; iprim++){   //loop on  tracks
 
     part = gAlice->Particle(iprim);
+    Double_t ptot=TMath::Sqrt(part->Px()*part->Px()+part->Py()*part->Py()+part->Pz()*part->Pz());
+    if (ptot<0.01) continue;
     char *xxx=strstr(part->GetName(),"XXX");
     if(xxx)continue;
 
@@ -622,7 +630,7 @@ Int_t FindPrimaries(Int_t nparticles){
 
     if(part->T()>timecut)continue;
 
-    Double_t ptot=TMath::Sqrt(part->Px()*part->Px()+part->Py()*part->Py()+part->Pz()*part->Pz());
+    //Double_t ptot=TMath::Sqrt(part->Px()*part->Px()+part->Py()*part->Py()+part->Pz()*part->Pz());
     if(ptot==(TMath::Abs(part->Pz())))continue; // no beam particles
 
     Bool_t prmch = kTRUE;   // candidate primary track
index 40cf79d..bf33e30 100644 (file)
@@ -19,6 +19,9 @@ Int_t AliTPCFindClusters(Int_t n=1) {
    TFile *in=TFile::Open("rfio:galice.root");
    if (!in->IsOpen()) {cerr<<"Can't open galice.root !\n"; return 2;}
 
+   TFile *ind=TFile::Open("digits.root");
+   if (!ind->IsOpen()) {cerr<<"Can't open galice.root !\n"; return 2;}
+
    if (!(gAlice=(AliRun*)in->Get("gAlice"))) {
      cerr<<"gAlice have not been found on galice.root !\n";
      return 3;
@@ -53,7 +56,7 @@ Int_t AliTPCFindClusters(Int_t n=1) {
       {
        // delete gAlice; gAlice=0;
        AliTPCv2 tpc; 
-       tpc.SetParam(dig); timer.Start(); cwd->cd();  
+       tpc.SetParam(dig); timer.Start(); ind->cd();  
        for (Int_t i=0;i<n;i++){
         printf("Processing event %d\n",i);
          tpc.Digits2Clusters(out,i);
diff --git a/TPC/AliTPCFindClustersMI.C b/TPC/AliTPCFindClustersMI.C
new file mode 100644 (file)
index 0000000..aaeb046
--- /dev/null
@@ -0,0 +1,113 @@
+/****************************************************************************
+ *           Origin: M.Ivanov marian.ivanov@cern.ch                         *
+ ****************************************************************************/
+
+/*
+
+  macro to create array of clusters from TPC digits
+  input files - galice.root 
+                digits.root - file with digits - usualy use link to galice.root
+                           - in splitted mode - neccesary to create link to proper file
+                           
+   output file - AliTPCclusters.root
+               - to be used by AliTPCTrackFinderMI.C
+
+  Warning - if cluster file AliTPCclusters.root already exist - macro exit and don't produce anything
+              
+*/
+
+
+#ifndef __CINT__
+#include <iostream.h>
+#include "AliRun.h"
+#include "AliTPCv1.h"
+#include "AliTPCv2.h"
+#include "AliTPCParam.h"
+#include "AliTPCclustererMI.h"
+#include "TFile.h"
+#include "TStopwatch.h"
+#include "TTree.h"
+#endif
+
+Int_t AliTPCFindClustersMI(Int_t n=1) {
+   TFile *out=TFile::Open("AliTPCclusters.root","new");
+   if (!out->IsOpen()) {cerr<<"Delete old AliTPCclusters.root !\n"; return 1;}
+   TFile *in=TFile::Open("rfio:galice.root");
+   if (!in->IsOpen()) {cerr<<"Can't open galice.root !\n"; return 2;}
+
+   TFile *ind=TFile::Open("rfio:digits.root");
+   if (!ind->IsOpen()) {cerr<<"Can't open digits file !\n"; return 2;}
+
+
+   if (!(gAlice=(AliRun*)in->Get("gAlice"))) {
+     cerr<<"gAlice have not been found on galice.root !\n";
+     return 3;
+   }
+
+   TDirectory *cwd = gDirectory;
+
+   AliTPC *TPC = (AliTPC*)gAlice->GetDetector("TPC"); 
+   Int_t ver = TPC->IsVersion(); 
+   cerr<<"TPC version "<<ver<<" has been found !\n";
+
+   AliTPCParam *dig=(AliTPCParam *)in->Get("75x40_100x60_150x60");
+   if (!dig) {cerr<<"TPC parameters have not been found !\n"; return 4;}
+
+   TStopwatch timer;
+
+   switch (ver) {
+   case 1:
+      cerr<<"Making clusters...\n";
+      {
+       AliTPCv1 &tpc=*((AliTPCv1*)TPC);
+       tpc.SetParam(dig); timer.Start(); cwd->cd(); 
+       for(Int_t i=0;i<n;i++){
+         printf("Processing event %d\n",i);
+         gAlice->GetEvent(i);
+         tpc.Hits2Clusters(out,i);
+       } 
+      }
+      break;
+   case 2:
+     cerr<<"Looking for clusters...\n";
+     {
+       // delete gAlice; gAlice=0;
+       AliTPCv2 tpc; 
+       tpc.SetParam(dig); timer.Start(); cwd->cd();  
+       
+       
+       for (Int_t i=0;i<n;i++){ 
+        AliTPCclustererMI clusterer;
+        char dname[100];
+        char cname[100];
+        sprintf(dname,"TreeD_75x40_100x60_150x60_%d",i);
+        sprintf(cname,"TreeC_TPC_%d",i);
+        TTree * input = (TTree*)ind->Get(dname);
+        out->cd();
+        TTree * output = new TTree(cname,cname); 
+        
+        printf("Processing event %d\n",i); 
+        clusterer.SetInput(input);
+        clusterer.SetOutput(output);
+        clusterer.Digits2Clusters(dig, i);
+         //tpc.Digits2Clusters(out,i);
+        //      AliTPCclusterer::Digits2Clusters(dig, out, i);
+       }
+     }
+     break;
+   default:
+     cerr<<"Invalid TPC version !\n";
+     return 5;
+   }
+   
+   timer.Stop(); timer.Print();
+   
+   delete gAlice; gAlice=0;
+
+   out->Close();
+
+   in->Close();
+
+   return 0;
+}
index 445d2d7..745de21 100644 (file)
@@ -14,6 +14,8 @@
 Int_t AliTPCFindTracks(Int_t eventn=1) {
 
    cerr<<"Looking for tracks...\n";
+   TFile f("galice.root");
+   gAlice = (AliRun*)f.Get("gAlice");
 
    TFile *out=TFile::Open("AliTPCtracks.root","new");
    if (!out->IsOpen()) {cerr<<"Delete old AliTPCtracks.root !\n"; return 1;}
diff --git a/TPC/AliTPCFindTracksMI.C b/TPC/AliTPCFindTracksMI.C
new file mode 100644 (file)
index 0000000..7e3acc0
--- /dev/null
@@ -0,0 +1,49 @@
+/****************************************************************************
+ *           Origin: I.Belikov, CERN, Jouri.Belikov@cern.ch                 *
+ ****************************************************************************/
+
+#ifndef __CINT__
+#include <iostream.h>
+#include "AliTPCParam.h"
+#include "AliTPCtrackerMI.h"
+#include "TFile.h"
+#include "TStopwatch.h"
+#include "AliRun.h"
+#include "AliMagF.h"
+#endif
+
+Int_t AliTPCFindTracks(Int_t eventn=1) { 
+   cerr<<"Looking for tracks...\n";
+   TFile f("galice.root");
+   gAlice = (AliRun*)f.Get("gAlice");
+   AliKalmanTrack::SetConvConst(1000/0.299792458/gAlice->Field()->SolenoidField());
+   TFile *out=TFile::Open("AliTPCtracks.root","new");
+   if (!out->IsOpen()) {cerr<<"Delete old AliTPCtracks.root !\n"; return 1;}
+
+   TFile *in=TFile::Open("AliTPCclusters.root");
+   if (!in->IsOpen()) {cerr<<"Can't open AliTPCclusters.root !\n"; return 2;}
+
+   AliTPCParam *par=(AliTPCParam*)in->Get("75x40_100x60_150x60");
+   if (!par) {cerr<<"Can't get TPC parameters !\n"; return 3;}
+   TStopwatch timer;
+  
+   Int_t rc=0;
+   for (Int_t i=0;i<eventn;i++){
+     printf("Processing event %d\n",i);
+     AliTPCtrackerMI *tracker = new AliTPCtrackerMI(par,i);
+     //delete tracker; 
+     //tracker = new AliTPCtrackerMI(par,i);
+     //Double_t xyz[]={0.,0.,0.}; tracker->SetVertex(xyz); //primary vertex
+     rc=tracker->Clusters2Tracks(0,out);
+     delete tracker;
+   }
+   timer.Stop(); timer.Print();
+   delete par; //Thanks to Mariana Bondila
+
+   in->Close();
+   out->Close();
+
+   return rc;
+}
index 34de70c..9575f02 100644 (file)
@@ -1,7 +1,15 @@
+/****************************************************************************
+ *           Origin: M.Ivanov, CERN, Marian.Ivanov@cern.ch                 *
+ ****************************************************************************/
+
 Int_t AliTPCHits2ExactClusters(Int_t nevent=1)
 {
+  //
+  // loop over hits and find itersection of the trak with pad rows
+  // so called exact clusters stored in the structure similar to hits
+  // it's used for comparison purpose - to have benchmark of cluster
+  // finder and tracker
 
-  // new version by J.Belikov
 
   // Connect the Root Galice file containing Geometry, Kine and Hits
 
diff --git a/TPC/AliTPCclusterMI.cxx b/TPC/AliTPCclusterMI.cxx
new file mode 100644 (file)
index 0000000..469ca53
--- /dev/null
@@ -0,0 +1,36 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+
+#include "AliTPCclusterMI.h"
+
+ClassImp(AliTPCclusterMI)
+ClassImp(AliTPCclusterLMI)
+
+Bool_t AliTPCclusterMI::IsSortable() const
+{
+  //
+  //
+  return kTRUE;
+
+}
+
+Int_t AliTPCclusterMI::Compare(const TObject* obj) const
+{
+  //
+  // compare according y
+  AliTPCclusterMI * o2 = (AliTPCclusterMI*)obj;
+  return (o2->GetY()>fY)? -1:1; 
+}
diff --git a/TPC/AliTPCclusterMI.h b/TPC/AliTPCclusterMI.h
new file mode 100644 (file)
index 0000000..f3b5311
--- /dev/null
@@ -0,0 +1,75 @@
+#ifndef ALITPCCLUSTERMI_H
+#define ALITPCCLUSTERMI_H
+
+//-------------------------------------------------------
+//                    TPC Cluster Class
+//
+//   Origin: Marian Ivanov
+//-------------------------------------------------------
+
+#include "AliCluster.h"
+#include "TMath.h"
+
+//_____________________________________________________________________________
+class AliTPCclusterMI : public AliCluster {
+public:
+  AliTPCclusterMI():AliCluster(){fQ=0;}
+  AliTPCclusterMI(Int_t *lab, Float_t *hit) : AliCluster(lab,hit) {fQ = (Short_t)hit[4];}
+  virtual Bool_t IsSortable() const; 
+  virtual Int_t Compare(const TObject* obj) const;
+  void Use() {fQ=-fQ;}
+  void SetQ(Float_t q) {fQ=(Short_t)q;}
+  void SetType(Char_t type) {fType=type;}
+  void SetMax(Short_t max) {fMax=max;}
+  Int_t IsUsed() const {return (fQ<0) ? 1 : 0;}
+
+  Float_t GetQ() const {return TMath::Abs(fQ);}
+  Float_t GetMax() const {return fMax;} 
+  Char_t  GetType()const {return fType;}
+private:
+  Short_t   fQ ;       //Q of cluster (in ADC counts)  
+  Char_t    fType;     //type of the cluster 0 means golden 
+  Short_t   fMax;      //maximal amplitude in cluster
+
+  ClassDef(AliTPCclusterMI,1)  // Time Projection Chamber clusters
+};
+
+class AliTPCclusterLMI  {
+ private:
+  Short_t  fCZ;       // current cluster position Z in cm - 100 mum precision
+  Short_t  fCY;       // current cluster position Y in cm - 100 mum precision
+  UChar_t  fSigmaZ;   // shape  Z - normalised shape - normaliziation 1 - precision 2 percent
+  UChar_t  fSigmaY;   // shape  Y - normalised shape - normaliziation 1 - precision 2 percent
+  UShort_t fQ;        // total charge in cluster 
+  UShort_t fMax;      // charge at maximum  
+  Char_t   fCType;    // type of the cluster
+  Int_t    fLabel[3];    // track indexes 
+ public:
+  AliTPCclusterLMI(){fCZ=fCY=fSigmaZ=fSigmaY=fQ=fMax=fCType=0;}
+  inline Float_t  GetZ()    {return (fCZ*0.01);}
+  inline Float_t  GetY()    {return (fCY*0.01);}
+  inline Float_t  GetSigmaZ() {return (fSigmaZ*0.02);}
+  inline Float_t  GetSigmaY() {return (fSigmaY*0.02);}  
+  inline Int_t  GetType()   {return fCType;}
+  inline Int_t  GetMax()   {return fMax;}
+  inline Float_t  GetQ()   {return fQ;}
+  inline Int_t    GelLabel(Int_t i){return fLabel[i];}
+  //
+  void     SetY(Float_t y){ fCY = Short_t(TMath::Nint(y*100.));} 
+  void     SetZ(Float_t z){ fCZ = Short_t(TMath::Nint(z*100.));} 
+  void     SetSigmaZ(Float_t sigmaz) {fSigmaZ = UChar_t(TMath::Nint(sigmaz*50.));}
+  void     SetSigmaY(Float_t sigmay) {fSigmaY = UChar_t(TMath::Nint(sigmay*50.));}
+  void     SetQ(Float_t q) {fQ = UShort_t(q);}
+  void     SetMax(Float_t max) {fMax = UShort_t(max);}
+  void     SetType(Char_t type) {fCType = type;}
+  void     SetLabels(Int_t labels[3]){fLabel[0] = labels[0];fLabel[1] = labels[1];fLabel[2] = labels[2];}
+  //
+ public:
+  ClassDef(AliTPCclusterLMI,1)  
+};
+
+
+
+#endif
+
+
diff --git a/TPC/AliTPCclustererMI.cxx b/TPC/AliTPCclustererMI.cxx
new file mode 100644 (file)
index 0000000..6b7d814
--- /dev/null
@@ -0,0 +1,588 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+
+
+//-------------------------------------------------------
+//          Implementation of the TPC clusterer
+//
+//   Origin: Marian Ivanov 
+//-------------------------------------------------------
+
+#include "AliTPCclustererMI.h"
+#include "AliTPCclusterMI.h"
+#include <TObjArray.h>
+#include <TFile.h>
+#include "AliTPCClustersArray.h"
+#include "AliTPCClustersRow.h"
+#include "AliDigits.h"
+#include "AliSimDigits.h"
+#include "AliTPCParam.h"
+#include <iostream.h>
+#include <TTree.h>
+
+ClassImp(AliTPCclustererMI)
+
+
+
+AliTPCclustererMI::AliTPCclustererMI()
+{
+  fInput =0;
+  fOutput=0;
+}
+void AliTPCclustererMI::SetInput(TTree * tree)
+{
+  //
+  // set input tree with digits
+  //
+  fInput = tree;  
+  if  (!fInput->GetBranch("Segment")){
+    cerr<<"AliTPC::Digits2Clusters(): no porper input tree !\n";
+    fInput=0;
+    return;
+  }
+}
+
+void AliTPCclustererMI::SetOutput(TTree * tree) 
+{
+  //
+  //
+  fOutput= tree;  
+  AliTPCClustersRow clrow;
+  AliTPCClustersRow *pclrow=&clrow;  
+  clrow.SetClass("AliTPCclusterMI");
+  clrow.SetArray(1); // to make Clones array
+  fOutput->Branch("Segment","AliTPCClustersRow",&pclrow,32000,200);    
+}
+
+
+Float_t  AliTPCclustererMI::GetSigmaY2(Int_t iz){
+  // sigma y2 = in digits  - we don't know the angle
+  Float_t z = iz*fParam->GetZWidth();
+  Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/
+    (fPadWidth*fPadWidth);
+  Float_t sres = 0.25;
+  Float_t res = sd2+sres;
+  return res;
+}
+
+
+Float_t  AliTPCclustererMI::GetSigmaZ2(Int_t iz){
+  //sigma z2 = in digits - angle estimated supposing vertex constraint
+  Float_t z = iz*fZWidth;
+  Float_t sd2 = (z*fParam->GetDiffL()*fParam->GetDiffL())/(fZWidth*fZWidth);
+  Float_t angular = fPadLength*(fParam->GetZLength()-z)/(fRx*fZWidth);
+  angular*=angular;
+  angular/=12.;
+  Float_t sres = fParam->GetZSigma()/fZWidth;
+  sres *=sres;
+  Float_t res = angular +sd2+sres;
+  return res;
+}
+
+void AliTPCclustererMI::MakeCluster(Int_t k,Int_t max,Int_t *bins, UInt_t m,
+AliTPCclusterMI &c) 
+{
+  Int_t i0=k/max;  //central pad
+  Int_t j0=k%max;  //central time bin
+
+  // set pointers to data
+  //Int_t dummy[5] ={0,0,0,0,0};
+  Int_t * matrix[5]; //5x5 matrix with digits  - indexing i = 0 ..4  j = -2..2
+  Int_t * resmatrix[5];
+  for (Int_t di=-2;di<=2;di++){
+    matrix[di+2]  =  &bins[k+di*max];
+    resmatrix[di+2]  =  &fResBins[k+di*max];
+  }
+  //build matrix with virtual charge
+  Float_t sigmay2= GetSigmaY2(j0);
+  Float_t sigmaz2= GetSigmaZ2(j0);
+
+  Float_t vmatrix[5][5];
+  vmatrix[2][2] = matrix[2][0];
+  c.SetType(0);
+  c.SetMax(Short_t(vmatrix[2][2])); // write maximal amplitude
+  for (Int_t di =-1;di <=1;di++)
+    for (Int_t dj =-1;dj <=1;dj++){
+      Float_t amp = matrix[di+2][dj];
+      if ( (amp<2) && (fLoop<2)){
+       // if under threshold  - calculate virtual charge
+       Float_t ratio = TMath::Exp(-1.2*TMath::Abs(di)/sigmay2)*TMath::Exp(-1.2*TMath::Abs(dj)/sigmaz2);
+       amp = ((matrix[2][0]-2)*(matrix[2][0]-2)/(matrix[-di+2][-dj]+2))*ratio;
+       if (amp>2) amp = 2;
+       vmatrix[2+di][2+dj]=amp;
+       vmatrix[2+2*di][2+2*dj]=0;
+       if ( (di*dj)!=0){       
+         //DIAGONAL ELEMENTS
+         vmatrix[2+2*di][2+dj] =0;
+         vmatrix[2+di][2+2*dj] =0;
+       }
+       continue;
+      }
+      if (amp<4){
+       //if small  amplitude - below  2 x threshold  - don't consider other one        
+       vmatrix[2+di][2+dj]=amp;
+       vmatrix[2+2*di][2+2*dj]=0;  // don't take to the account next bin
+       if ( (di*dj)!=0){       
+         //DIAGONAL ELEMENTS
+         vmatrix[2+2*di][2+dj] =0;
+         vmatrix[2+di][2+2*dj] =0;
+       }
+       continue;
+      }
+      //if bigger then take everything
+      vmatrix[2+di][2+dj]=amp;
+      vmatrix[2+2*di][2+2*dj]= matrix[2*di+2][2*dj] ;      
+      if ( (di*dj)!=0){       
+         //DIAGONAL ELEMENTS
+         vmatrix[2+2*di][2+dj] = matrix[2*di+2][dj];
+         vmatrix[2+di][2+2*dj] = matrix[2+di][dj*2];
+       }      
+    }
+
+
+  
+  Float_t sumw=0;
+  Float_t sumiw=0;
+  Float_t sumi2w=0;
+  Float_t sumjw=0;
+  Float_t sumj2w=0;
+  //
+  for (Int_t i=-2;i<=2;i++)
+    for (Int_t j=-2;j<=2;j++){
+      Float_t amp = vmatrix[i+2][j+2];
+
+      sumw    += amp;
+      sumiw   += i*amp;
+      sumi2w  += i*i*amp;
+      sumjw   += j*amp;
+      sumj2w  += j*j*amp;
+    }    
+  //   
+  Float_t meani = sumiw/sumw;
+  Float_t mi2   = sumi2w/sumw-meani*meani;
+  Float_t meanj = sumjw/sumw;
+  Float_t mj2   = sumj2w/sumw-meanj*meanj;
+  //
+  Float_t ry = mi2/sigmay2;
+  Float_t rz = mj2/sigmaz2;
+  
+  //
+  if ( ( (ry<0.6) || (rz<0.6) ) && fLoop==2) return;
+  if ( (ry <1.2) && (rz<1.2) ) {
+    //if cluster looks like expected 
+    //+1.2 deviation from expected sigma accepted
+    //    c.fMax = FitMax(vmatrix,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
+
+    meani +=i0;
+    meanj +=j0;
+    //set cluster parameters
+    c.SetQ(sumw);
+    c.SetY(meani*fPadWidth); 
+    c.SetZ(meanj*fZWidth); 
+    c.SetSigmaY2(mi2);
+    c.SetSigmaZ2(mj2);
+    AddCluster(c);
+    //remove cluster data from data
+    for (Int_t di=-2;di<=2;di++)
+      for (Int_t dj=-2;dj<=2;dj++){
+       resmatrix[di+2][dj] -= Int_t(vmatrix[di+2][dj+2]);
+       if (resmatrix[di+2][dj]<0) resmatrix[di+2][dj]=0;
+      }
+    resmatrix[2][0] =0;
+    return;     
+  }
+  //
+  //unfolding when neccessary  
+  //
+  
+  Int_t * matrix2[7]; //7x7 matrix with digits  - indexing i = 0 ..6  j = -3..3
+  Int_t dummy[7]={0,0,0,0,0,0};
+  for (Int_t di=-3;di<=3;di++){
+    matrix2[di+3] =  &bins[k+di*max];
+    if ((k+di*max)<3)  matrix2[di+3] = &dummy[3];
+    if ((k+di*max)>fMaxBin-3)  matrix2[di+3] = &dummy[3];
+  }
+  Float_t vmatrix2[5][5];
+  Float_t sumu;
+  Float_t overlap;
+  UnfoldCluster(matrix2,vmatrix2,meani,meanj,sumu,overlap);
+  //
+  //  c.fMax = FitMax(vmatrix2,meani,meanj,TMath::Sqrt(sigmay2),TMath::Sqrt(sigmaz2));
+  meani +=i0;
+  meanj +=j0;
+  //set cluster parameters
+  c.SetQ(sumu);
+  c.SetY(meani*fPadWidth); 
+  c.SetZ(meanj*fZWidth); 
+  c.SetSigmaY2(mi2);
+  c.SetSigmaZ2(mj2);
+  c.SetType(Char_t(overlap)+1);
+  AddCluster(c);
+
+  //unfolding 2
+  meani-=i0;
+  meanj-=j0;
+  if (gDebug>4)
+    printf("%f\t%f\n", vmatrix2[2][2], vmatrix[2][2]);
+}
+
+
+
+void AliTPCclustererMI::UnfoldCluster(Int_t * matrix2[7], Float_t recmatrix[5][5], Float_t & meani, Float_t & meanj, 
+                                     Float_t & sumu, Float_t & overlap )
+{
+  //
+  //unfold cluster from input matrix
+  //data corresponding to cluster writen in recmatrix
+  //output meani and meanj
+
+  //take separatelly y and z
+
+  Float_t sum3i[7] = {0,0,0,0,0,0,0};
+  Float_t sum3j[7] = {0,0,0,0,0,0,0};
+
+  for (Int_t k =0;k<7;k++)
+    for (Int_t l = -1; l<=1;l++){
+      sum3i[k]+=matrix2[k][l];
+      sum3j[k]+=matrix2[l+3][k-3];
+    }
+  Float_t mratio[3][3]={{1,1,1},{1,1,1},{1,1,1}};
+  //
+  //unfold  y 
+  Float_t sum3wi    = 0;  //charge minus overlap
+  Float_t sum3wio   = 0;  //full charge
+  Float_t sum3iw    = 0;  //sum for mean value
+  for (Int_t dk=-1;dk<=1;dk++){
+    sum3wio+=sum3i[dk+3];
+    if (dk==0){
+      sum3wi+=sum3i[dk+3];     
+    }
+    else{
+      Float_t ratio =1;
+      if (  ( ((sum3i[dk+3]+3)/(sum3i[3]-3))+1 < (sum3i[2*dk+3]-3)/(sum3i[dk+3]+3))||
+          sum3i[dk+3]<=sum3i[2*dk+3] && sum3i[dk+3]>2 ){
+       Float_t xm2 = sum3i[-dk+3];
+       Float_t xm1 = sum3i[+3];
+       Float_t x1  = sum3i[2*dk+3];
+       Float_t x2  = sum3i[3*dk+3];    
+       Float_t w11   = TMath::Max(4.*xm1-xm2,0.000001);          
+       Float_t w12   = TMath::Max(4 *x1 -x2,0.);
+       ratio = w11/(w11+w12);   
+       for (Int_t dl=-1;dl<=1;dl++)
+         mratio[dk+1][dl+1] *= ratio;
+      }
+      Float_t amp = sum3i[dk+3]*ratio;
+      sum3wi+=amp;
+      sum3iw+= dk*amp;      
+    }
+  }
+  meani = sum3iw/sum3wi;
+  Float_t overlapi = (sum3wio-sum3wi)/sum3wio;
+
+
+
+  //unfold  z 
+  Float_t sum3wj    = 0;  //charge minus overlap
+  Float_t sum3wjo   = 0;  //full charge
+  Float_t sum3jw    = 0;  //sum for mean value
+  for (Int_t dk=-1;dk<=1;dk++){
+    sum3wjo+=sum3j[dk+3];
+    if (dk==0){
+      sum3wj+=sum3j[dk+3];     
+    }
+    else{
+      Float_t ratio =1;
+      if ( ( ((sum3j[dk+3]+3)/(sum3j[3]-3))+1 < (sum3j[2*dk+3]-3)/(sum3j[dk+3]+3)) ||
+          (sum3j[dk+3]<=sum3j[2*dk+3] && sum3j[dk+3]>2)){
+       Float_t xm2 = sum3j[-dk+3];
+       Float_t xm1 = sum3j[+3];
+       Float_t x1  = sum3j[2*dk+3];
+       Float_t x2  = sum3j[3*dk+3];    
+       Float_t w11   = TMath::Max(4.*xm1-xm2,0.000001);          
+       Float_t w12   = TMath::Max(4 *x1 -x2,0.);
+       ratio = w11/(w11+w12);   
+       for (Int_t dl=-1;dl<=1;dl++)
+         mratio[dl+1][dk+1] *= ratio;
+      }
+      Float_t amp = sum3j[dk+3]*ratio;
+      sum3wj+=amp;
+      sum3jw+= dk*amp;      
+    }
+  }
+  meanj = sum3jw/sum3wj;
+  Float_t overlapj = (sum3wjo-sum3wj)/sum3wjo;  
+  overlap = Int_t(100*TMath::Max(overlapi,overlapj)+3);  
+  sumu = (sum3wj+sum3wi)/2.;
+  
+  if (overlap ==3) {
+    //if not overlap detected remove everything
+    for (Int_t di =-2; di<=2;di++)
+      for (Int_t dj =-2; dj<=2;dj++){
+       recmatrix[di+2][dj+2] = matrix2[3+di][dj];
+      }
+  }
+  else{
+    for (Int_t di =-1; di<=1;di++)
+      for (Int_t dj =-1; dj<=1;dj++){
+       Float_t ratio =1;
+       if (mratio[di+1][dj+1]==1){
+         recmatrix[di+2][dj+2]     = matrix2[3+di][dj];
+         if (TMath::Abs(di)+TMath::Abs(dj)>1){
+           recmatrix[2*di+2][dj+2] = matrix2[3+2*di][dj];
+           recmatrix[di+2][2*dj+2] = matrix2[3+di][2*dj];
+         }       
+         recmatrix[2*di+2][2*dj+2] = matrix2[3+2*di][2*dj];
+       }
+       else
+         {
+           //if we have overlap in direction
+           recmatrix[di+2][dj+2] = mratio[di+1][dj+1]* matrix2[3+di][dj];    
+           if (TMath::Abs(di)+TMath::Abs(dj)>1){
+             ratio =  TMath::Min(recmatrix[di+2][dj+2]/(matrix2[3+0*di][1*dj]+1),1.);
+             recmatrix[2*di+2][dj+2] = ratio*recmatrix[di+2][dj+2];
+             //
+             ratio =  TMath::Min(recmatrix[di+2][dj+2]/(matrix2[3+1*di][0*dj]+1),1.);
+             recmatrix[di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
+           }
+           else{
+             ratio =  recmatrix[di+2][dj+2]/matrix2[3][0];
+             recmatrix[2*di+2][2*dj+2] = ratio*recmatrix[di+2][dj+2];
+           }
+         }
+      }
+  }
+  if (gDebug>4) 
+    printf("%f\n", recmatrix[2][2]);
+  
+}
+
+Float_t AliTPCclustererMI::FitMax(Float_t vmatrix[5][5], Float_t y, Float_t z, Float_t sigmay, Float_t sigmaz)
+{
+  //
+  // estimate max
+  Float_t sumteor= 0;
+  Float_t sumamp = 0;
+
+  for (Int_t di = -1;di<=1;di++)
+    for (Int_t dj = -1;dj<=1;dj++){
+      if (vmatrix[2+di][2+dj]>2){
+       Float_t teor = TMath::Gaus(di,y,sigmay*1.2)*TMath::Gaus(dj,z,sigmaz*1.2);
+       sumteor += teor*vmatrix[2+di][2+dj];
+       sumamp  += vmatrix[2+di][2+dj]*vmatrix[2+di][2+dj];
+      }
+    }    
+  Float_t max = sumamp/sumteor;
+  return max;
+}
+
+void AliTPCclustererMI::AddCluster(AliTPCclusterMI &c){
+  //
+  // transform cluster to the global coordinata
+  // add the cluster to the array
+  //
+  Float_t meani = c.GetY()/fPadWidth;
+  Float_t meanj = c.GetZ()/fZWidth;
+
+  Int_t ki = TMath::Nint(meani-3);
+  if (ki<0) ki=0;
+  if (ki>=fMaxPad) ki = fMaxPad-1;
+  Int_t kj = TMath::Nint(meanj-3);
+  if (kj<0) kj=0;
+  if (kj>=fMaxTime-3) kj=fMaxTime-4;
+  // ki and kj shifted to "real" coordinata
+  c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,0)-2,0);
+  c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,1)-2,1);
+  c.SetLabel(fRowDig->GetTrackIDFast(kj,ki,2)-2,2);
+  
+  
+  Float_t s2 = c.GetSigmaY2();
+  Float_t w=fParam->GetPadPitchWidth(fSector);
+  
+  c.SetSigmaY2(s2*w*w);
+  s2 = c.GetSigmaZ2(); 
+  w=fZWidth;
+  c.SetSigmaZ2(s2*w*w);
+  c.SetY((meani - 2.5 - 0.5*fMaxPad)*fParam->GetPadPitchWidth(fSector));
+  c.SetZ(fZWidth*(meanj-3)); 
+  c.SetZ(c.GetZ()   - 3.*fParam->GetZSigma()); // PASA delay 
+  c.SetZ(fSign*(fParam->GetZLength() - c.GetZ()));
+  
+  if (ki<=1 || ki>=fMaxPad-1 || kj==1 || kj==fMaxTime-2) {
+    //c.SetSigmaY2(c.GetSigmaY2()*25.);
+    //c.SetSigmaZ2(c.GetSigmaZ2()*4.);
+    c.SetType(-(c.GetType()+3));  //edge clusters
+  }
+  if (fLoop==2) c.SetType(100);
+
+  TClonesArray * arr = fRowCl->GetArray();
+  AliTPCclusterMI * cl = new ((*arr)[fNcluster]) AliTPCclusterMI(c);
+
+  fNcluster++;
+}
+
+
+//_____________________________________________________________________________
+void AliTPCclustererMI::Digits2Clusters(const AliTPCParam *par, Int_t eventn)
+{
+  //-----------------------------------------------------------------
+  // This is a simple cluster finder.
+  //-----------------------------------------------------------------
+  TDirectory *savedir=gDirectory; 
+
+ if (fInput==0){ 
+    cerr<<"AliTPC::Digits2Clusters(): input tree not initialised !\n";
+    return;
+  }
+  if (fOutput==0) {
+     cerr<<"AliTPC::Digits2Clusters(): output tree not initialised !\n";
+     return;
+  }
+
+  AliSimDigits digarr, *dummy=&digarr;
+  fRowDig = dummy;
+  fInput->GetBranch("Segment")->SetAddress(&dummy);
+  Stat_t nentries = fInput->GetEntries();
+  
+  fMaxTime=par->GetMaxTBin()+6; // add 3 virtual time bins before and 3   after
+    
+  fParam = par;
+  ((AliTPCParam*)par)->Write(par->GetTitle());
+
+  Int_t nclusters  = 0;
+  
+  for (Int_t n=0; n<nentries; n++) {
+    fInput->GetEvent(n);
+    Int_t row;
+    if (!par->AdjustSectorRow(digarr.GetID(),fSector,row)) {
+      cerr<<"AliTPC warning: invalid segment ID ! "<<digarr.GetID()<<endl;
+      continue;
+    }
+        
+    AliTPCClustersRow *clrow= new AliTPCClustersRow();
+    fRowCl = clrow;
+    clrow->SetClass("AliTPCclusterMI");
+    clrow->SetArray(1);
+
+    clrow->SetID(digarr.GetID());
+    fOutput->GetBranch("Segment")->SetAddress(&clrow);
+    fRx=par->GetPadRowRadii(fSector,row);
+    
+    
+    const Int_t kNIS=par->GetNInnerSector(), kNOS=par->GetNOuterSector();
+    fZWidth = fParam->GetZWidth();
+    if (fSector < kNIS) {
+      fMaxPad = par->GetNPadsLow(row);
+      fSign = (fSector < kNIS/2) ? 1 : -1;
+      fPadLength = par->GetPadPitchLength(fSector,row);
+      fPadWidth = par->GetPadPitchWidth();
+    } else {
+      fMaxPad = par->GetNPadsUp(row);
+      fSign = ((fSector-kNIS) < kNOS/2) ? 1 : -1;
+      fPadLength = par->GetPadPitchLength(fSector,row);
+      fPadWidth  = par->GetPadPitchWidth();
+    }
+    
+    
+    fMaxBin=fMaxTime*(fMaxPad+6);  // add 3 virtual pads  before and 3 after
+    fBins    =new Int_t[fMaxBin];
+    fResBins =new Int_t[fMaxBin];  //fBins with residuals after 1 finder loop 
+    memset(fBins,0,sizeof(Int_t)*fMaxBin);
+    
+    if (digarr.First()) //MI change
+      do {
+       Short_t dig=digarr.CurrentDigit();
+       if (dig<=par->GetZeroSup()) continue;
+       Int_t j=digarr.CurrentRow()+3, i=digarr.CurrentColumn()+3;
+       fBins[i*fMaxTime+j]=dig;
+      } while (digarr.Next());
+    digarr.ExpandTrackBuffer();
+
+    //add virtual charge at the edge   
+    for (Int_t i=0; i<fMaxTime; i++){
+      Float_t amp1 = fBins[i+3*fMaxTime]; 
+      Float_t amp0 =0;
+      if (amp1>0){
+       Float_t amp2 = fBins[i+4*fMaxTime];
+       if (amp2==0) amp2=0.5;
+       Float_t sigma2 = GetSigmaY2(i);         
+       amp0 = (amp1*amp1/amp2)*TMath::Exp(-1./sigma2); 
+       if (gDebug>4) printf("\n%f\n",amp0);
+      }  
+      fBins[i+2*fMaxTime] = Int_t(amp0);    
+      amp0 = 0;
+      amp1 = fBins[(fMaxPad+2)*fMaxTime+i];
+      if (amp1>0){
+       Float_t amp2 = fBins[i+(fMaxPad+1)*fMaxTime];
+       if (amp2==0) amp2=0.5;
+       Float_t sigma2 = GetSigmaY2(i);         
+       amp0 = (amp1*amp1/amp2)*TMath::Exp(-1./sigma2);         
+       if (gDebug>4) printf("\n%f\n",amp0);
+      }        
+      fBins[(fMaxPad+3)*fMaxTime+i] = Int_t(amp0);           
+    }
+
+    memcpy(fResBins,fBins, fMaxBin*2);
+    //
+    fNcluster=0;
+    //first loop - for "gold cluster" 
+    fLoop=1;
+    Int_t *b=&fBins[-1]+2*fMaxTime;
+    Int_t crtime = (fParam->GetZLength()-1.05*fRx)/fZWidth-5;
+
+    for (Int_t i=2*fMaxTime; i<fMaxBin-2*fMaxTime; i++) {
+      b++;
+      if (*b<8) continue;   //threshold form maxima
+      if (i%fMaxTime<crtime) {
+       Int_t delta = -(i%fMaxTime)+crtime;
+       b+=delta;
+       i+=delta;
+       continue; 
+      }
+     
+      if (!IsMaximum(*b,fMaxTime,b)) continue;
+      AliTPCclusterMI c;
+      Int_t dummy;
+      MakeCluster(i, fMaxTime, fBins, dummy,c);
+      //}
+    }
+    //memcpy(fBins,fResBins, fMaxBin*2);
+    //second  loop - for rest cluster 
+    /*        
+    fLoop=2;
+    b=&fResBins[-1]+2*fMaxTime;
+    for (Int_t i=2*fMaxTime; i<fMaxBin-2*fMaxTime; i++) {
+      b++;
+      if (*b<25) continue;   // bigger threshold for maxima
+      if (!IsMaximum(*b,fMaxTime,b)) continue;
+      AliTPCclusterMI c;
+      Int_t dummy;
+      MakeCluster(i, fMaxTime, fResBins, dummy,c);
+      //}
+    }
+    */
+
+    fOutput->Fill();
+    delete clrow;    
+    nclusters+=fNcluster;    
+    delete[] fBins;      
+    delete[] fResBins;  
+  }  
+  cerr<<"Number of found clusters : "<<nclusters<<"                        \n";
+  
+  fOutput->Write();
+  savedir->cd();
+}
diff --git a/TPC/AliTPCclustererMI.h b/TPC/AliTPCclustererMI.h
new file mode 100644 (file)
index 0000000..43c750d
--- /dev/null
@@ -0,0 +1,82 @@
+#ifndef ALITPCCLUSTERERMI_H
+#define ALITPCCLUSTERERMI_H
+
+//-------------------------------------------------------
+//                       TPC clusterer
+//
+//   Origin: Marian Ivanov  
+//-------------------------------------------------------
+#include <Rtypes.h>
+#include <TObject.h>
+#define kMAXCLUSTER 2500
+
+class TFile;
+class AliTPCParam;
+class AliTPCclusterMI;
+class AliTPCClustersRow;
+class AliSimDigits;
+class TTree;
+
+class AliTPCclustererMI : public TObject{
+public:
+  AliTPCclustererMI();
+  virtual void Digits2Clusters(const AliTPCParam *par, Int_t eventn=1);
+  virtual void SetInput(TTree * tree);  // set input tree with digits    
+  virtual void SetOutput(TTree * tree); //set output tree with 
+private:
+  Bool_t IsMaximum(Int_t k, Int_t max, const Int_t *bins); 
+  void MakeCluster2(Int_t k,Int_t max,Int_t *bins,UInt_t m,
+   AliTPCclusterMI &c);  
+  void MakeCluster(Int_t k,Int_t max,Int_t *bins,UInt_t m,
+   AliTPCclusterMI &c); 
+  Float_t  GetSigmaY2(Int_t iz);
+  Float_t  GetSigmaZ2(Int_t iz);
+  Float_t  FitMax(Float_t vmatrix[5][5], Float_t y, Float_t z, Float_t sigmay, Float_t sigmaz);
+  void AddCluster(AliTPCclusterMI &c);  // add the cluster to the array
+  void UnfoldCluster(Int_t * matrix[7], Float_t recmatrix[5][5], 
+                    Float_t & meani, Float_t & meanj, Float_t & sum, Float_t &overlap );
+
+
+
+  Int_t * fBins;       //!digits array
+  Int_t * fResBins;    //!digits array with res. after 1 finder
+  Int_t fLoop;         //loop - cf in 2 loops
+  Int_t fMaxBin;
+  Int_t fMaxTime;
+  Int_t fMaxPad;
+  Int_t fSector;      //!current sector
+  Float_t fSign;      //!current sign 
+  Float_t fRx;        // current radius
+  Float_t fPadWidth;  // the width of the pad
+  Float_t fPadLength;  // the width of the pad
+  Float_t fZWidth;     //the z bin width
+
+  TTree * fInput;   //!input  tree with digits - object not owner
+  TTree * fOutput;   //!output tree with digits - object not owner
+  AliTPCClustersRow * fRowCl;  //! current cluster row
+  AliSimDigits * fRowDig;      //! current digits row
+  const AliTPCParam * fParam;        //! tpc parameters
+  Int_t fNcluster;             // number of clusters - for given row
+  ClassDef(AliTPCclustererMI,1)  // Time Projection Chamber digits
+};
+
+inline Bool_t AliTPCclustererMI::IsMaximum(Int_t q,Int_t max,const Int_t *bins){
+  //is this a local maximum ?
+  if (bins[-max] >= q) return kFALSE;
+  if (bins[-1  ] >= q) return kFALSE; 
+  if (bins[+max] > q) return kFALSE; 
+  if (bins[+1  ] > q) return kFALSE; 
+  if (bins[-max-1] >= q) return kFALSE;
+  if (bins[+max-1] >= q) return kFALSE; 
+  if (bins[+max+1] > q) return kFALSE; 
+  if (bins[-max+1] >= q) return kFALSE;
+  return kTRUE; 
+}
+
+
+
+//-----------------------------------------------------------------
+
+#endif
+
+
index 0b8063d..64c1fa8 100644 (file)
@@ -85,7 +85,7 @@ public:
   Int_t Update(const AliCluster* c, Double_t chi2, UInt_t i);
   void ResetCovariance();
 
-private: 
+protected: 
   Double_t fX;              // X-coordinate of this track (reference plane)
   Double_t fAlpha;          // Rotation angle the local (TPC sector)
                             // coordinate system and the global ALICE one.
diff --git a/TPC/AliTPCtrackerMI.cxx b/TPC/AliTPCtrackerMI.cxx
new file mode 100644 (file)
index 0000000..16ac730
--- /dev/null
@@ -0,0 +1,1992 @@
+/**************************************************************************
+ * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ *                                                                        *
+ * Author: The ALICE Off-line Project.                                    *
+ * Contributors are mentioned in the code where appropriate.              *
+ *                                                                        *
+ * Permission to use, copy, modify and distribute this software and its   *
+ * documentation strictly for non-commercial purposes is hereby granted   *
+ * without fee, provided that the above copyright notice appears in all   *
+ * copies and that both the copyright notice and this permission notice   *
+ * appear in the supporting documentation. The authors make no claims     *
+ * about the suitability of this software for any purpose. It is          *
+ * provided "as is" without express or implied warranty.                  *
+ **************************************************************************/
+
+
+
+
+
+
+/*
+  AliTPC parallel tracker - 
+  How to use?  - 
+  run AliTPCFindClusters.C macro - clusters neccessary for tracker are founded
+  run AliTPCFindTracksMI.C macro - to find tracks
+  tracks are written to AliTPCtracks.root file
+  for comparison also seeds are written to the same file - to special branch
+*/
+
+//-------------------------------------------------------
+//          Implementation of the TPC tracker
+//
+//   Origin: Marian Ivanov   Marian.Ivanov@cern.ch
+// 
+//-------------------------------------------------------
+
+#include <TObjArray.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <iostream.h>
+
+#include "AliTPCtrackerMI.h"
+#include "AliTPCclusterMI.h"
+#include "AliTPCParam.h"
+#include "AliTPCClustersRow.h"
+#include "AliComplexCluster.h"
+#include "TStopwatch.h"
+
+
+ClassImp(AliTPCseed)
+
+AliTPCclusterTracks::AliTPCclusterTracks(){
+  // class for storing overlaping info
+  fTrackIndex[0]=-1;
+  fTrackIndex[1]=-1;
+  fTrackIndex[2]=-1;
+  fDistance[0]=1000;
+  fDistance[1]=1000;
+  fDistance[2]=1000;
+}
+
+
+
+
+
+Int_t AliTPCtrackerMI::UpdateTrack(AliTPCseed * track, AliTPCclusterMI* c, Double_t chi2, UInt_t i){
+
+  Int_t sec=(i&0xff000000)>>24; 
+  Int_t row = (i&0x00ff0000)>>16;
+  track->fRow=(i&0x00ff0000)>>16;
+  track->fSector = sec;
+  //  Int_t index = i&0xFFFF;
+  if (sec>=fParam->GetNInnerSector()) track->fRow += fParam->GetNRowLow(); 
+  track->fClusterIndex[track->fRow] = i;
+  track->fFirstPoint = row;
+  //
+
+  AliTPCTrackPoint   *trpoint =track->GetTrackPoint(track->fRow);
+  Float_t angle2 = track->GetSnp()*track->GetSnp();
+  angle2 = TMath::Sqrt(angle2/(1-angle2)); 
+  //
+  //SET NEW Track Point
+  //
+  if (c!=0){
+    //if we have a cluster
+    trpoint->GetCPoint().SetY(c->GetY());
+    trpoint->GetCPoint().SetZ(c->GetZ());    
+    //
+    trpoint->GetCPoint().SetSigmaY(c->GetSigmaY2()/(track->fCurrentSigmaY*track->fCurrentSigmaY));
+    trpoint->GetCPoint().SetSigmaZ(c->GetSigmaZ2()/(track->fCurrentSigmaZ*track->fCurrentSigmaZ));
+    //
+    trpoint->GetCPoint().SetType(c->GetType());
+    trpoint->GetCPoint().SetQ(c->GetQ());
+    trpoint->GetCPoint().SetMax(c->GetMax());
+    //  
+    trpoint->GetCPoint().SetErrY(TMath::Sqrt(track->fErrorY2));
+    trpoint->GetCPoint().SetErrZ(TMath::Sqrt(track->fErrorZ2));
+    //
+  }
+  trpoint->GetTPoint().SetX(track->GetX());
+  trpoint->GetTPoint().SetY(track->GetY());
+  trpoint->GetTPoint().SetZ(track->GetZ());
+  //
+  trpoint->GetTPoint().SetAngleY(angle2);
+  trpoint->GetTPoint().SetAngleZ(track->GetTgl());
+  
+
+  
+  if (chi2>10){
+    //    printf("suspicious chi2 %f\n",chi2);
+  }
+  //  if (track->fIsSeeding){ 
+    track->fErrorY2 *= 1.2;
+    track->fErrorY2 += 0.0064;    
+    track->fErrorZ2 *= 1.2;   
+    track->fErrorY2 += 0.005;    
+    //}
+
+  return track->Update(c,chi2,i);
+
+}
+//_____________________________________________________________________________
+AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCParam *par, Int_t eventn): 
+AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
+{
+  //---------------------------------------------------------------------
+  // The main TPC tracker constructor
+  //---------------------------------------------------------------------
+  fInnerSec=new AliTPCSector[fkNIS];         
+  fOuterSec=new AliTPCSector[fkNOS];
+
+  Int_t i;
+  for (i=0; i<fkNIS; i++) fInnerSec[i].Setup(par,0);
+  for (i=0; i<fkNOS; i++) fOuterSec[i].Setup(par,1);
+
+  fN=0;  fSectors=0;
+
+  fClustersArray.Setup(par);
+  fClustersArray.SetClusterType("AliTPCclusterMI");
+
+  char   cname[100];
+  if (eventn==-1) {
+    sprintf(cname,"TreeC_TPC");
+  }
+  else {
+    sprintf(cname,"TreeC_TPC_%d",eventn);
+  }
+
+  fClustersArray.ConnectTree(cname);
+
+  fEventN = eventn;
+  fSeeds=0;
+  fNtracks = 0;
+  fParam = par;
+}
+
+//_____________________________________________________________________________
+AliTPCtrackerMI::~AliTPCtrackerMI() {
+  //------------------------------------------------------------------
+  // TPC tracker destructor
+  //------------------------------------------------------------------
+  delete[] fInnerSec;
+  delete[] fOuterSec;
+  if (fSeeds) {
+    fSeeds->Delete(); 
+    delete fSeeds;
+  }
+}
+
+
+Double_t AliTPCtrackerMI::ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl){
+  //
+  //
+  Float_t snoise2;
+  Float_t z = fParam->GetZLength()-TMath::Abs(seed->GetZ());
+
+  //cluster "quality"
+  Float_t rsigmay = 1;
+  Int_t ctype = 0;
+
+  //standard if we don't have cluster - take MIP
+  const Float_t chmip = 50.; 
+  Float_t amp = chmip/0.3;  
+  Float_t nel;
+  Float_t nprim;
+  if (cl){
+    amp = cl->GetQ();
+    rsigmay = cl->GetSigmaY2()/(seed->fCurrentSigmaY*seed->fCurrentSigmaY);
+    ctype = cl->GetType();
+  }
+  
+
+  Float_t landau=2 ;    //landau fluctuation part
+  Float_t gg=2;         // gg fluctuation part
+  Float_t padlength= fSectors->GetPadPitchLength(seed->GetX());
+
+
+  if (fSectors==fInnerSec){
+    snoise2 = 0.0004;
+    nel     = 0.268*amp;
+    nprim   = 0.155*amp;
+    gg      = (2+0.0002*amp)/nel;
+    landau  = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim;
+    if (landau>1) landau=1;
+  }
+  else {
+    snoise2 = 0.0004;
+    nel     = 0.3*amp;
+    nprim   = 0.133*amp;
+    gg      = (2+0.0002*amp)/nel;
+    landau  = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim;
+    if (landau>1) landau=1;
+  }
+
+
+  Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z;
+  Float_t angle2 = seed->GetSnp()*seed->GetSnp();
+  angle2 = angle2/(1-angle2); 
+  Float_t angular = landau*angle2*padlength*padlength/12.;
+  Float_t res = sdiff + angular;
+  
+  
+  if ((ctype==0) && (fSectors ==fOuterSec))
+    res *= 0.78 +TMath::Exp(7.4*(rsigmay-1.2));
+
+  if ((ctype==0) && (fSectors ==fInnerSec))
+    res *= 0.72 +TMath::Exp(3.36*(rsigmay-1.2));
+
+
+  if ((ctype>0))
+    res*= TMath::Power((rsigmay+0.5),1.5)+0.0064;
+  
+  if (ctype<0)
+    res*=2.4;  // overestimate error 2 times
+  
+  res+= snoise2;
+  if (res<2*snoise2)
+    res = 2*snoise2;
+
+  seed->SetErrorY2(res);
+  return res;
+}
+
+
+
+
+Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
+  //
+  //
+  Float_t snoise2;
+  Float_t z = fParam->GetZLength()-TMath::Abs(seed->GetZ());
+  //signal quality
+  Float_t rsigmaz=1;
+  Int_t ctype =0;
+
+  const Float_t chmip = 50.;
+  Float_t amp = chmip/0.3;  
+  Float_t nel;
+  Float_t nprim;
+  if (cl){
+    amp = cl->GetQ();
+    rsigmaz = cl->GetSigmaZ2()/(seed->fCurrentSigmaZ*seed->fCurrentSigmaZ);
+    ctype = cl->GetType();
+  }
+
+  //
+  Float_t landau=2 ;    //landau fluctuation part
+  Float_t gg=2;         // gg fluctuation part
+  Float_t padlength= fSectors->GetPadPitchLength(seed->GetX());
+
+  if (fSectors==fInnerSec){
+    snoise2 = 0.0004;
+    nel     = 0.268*amp;
+    nprim   = 0.155*amp;
+    gg      = (2+0.0002*amp)/nel;
+    landau  = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim;
+    if (landau>1) landau=1;
+  }
+  else {
+    snoise2 = 0.0004;
+    nel     = 0.3*amp;
+    nprim   = 0.133*amp;
+    gg      = (2+0.0002*amp)/nel;
+    landau  = (2.+0.12*nprim)*0.5*(amp*amp/40000.+2)/nprim;
+    if (landau>1) landau=1;
+  }
+  Float_t sdiff = gg*fParam->GetDiffT()*fParam->GetDiffT()*z;
+  
+  Float_t angle = seed->GetTgl();
+  Float_t angular = landau*angle*angle*padlength*padlength/12.;
+  Float_t res = sdiff + angular;
+
+  if ((ctype==0) && (fSectors ==fOuterSec))
+    res *= 0.81 +TMath::Exp(6.8*(rsigmaz-1.2));
+
+  if ((ctype==0) && (fSectors ==fInnerSec))
+    res *= 0.72 +TMath::Exp(2.04*(rsigmaz-1.2));
+  if ((ctype>0))
+    res*= TMath::Power(rsigmaz+0.5,1.5)+0.0064;  //0.31+0.147*ctype;
+  if (ctype<0)
+    res*=1.3;
+  if ((ctype<0) &&amp<70)
+    res*=1.3;  
+
+  res += snoise2;
+  if (res<2*snoise2)
+     res = 2*snoise2;
+
+  seed->SetErrorZ2(res);
+  return res;
+}
+
+
+
+Int_t  AliTPCseed::GetProlongation(Double_t xk, Double_t &y, Double_t & z) const
+{
+  //-----------------------------------------------------------------
+  // This function find proloncation of a track to a reference plane x=xk.
+  // doesn't change internal state of the track
+  //-----------------------------------------------------------------
+  
+  Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1;
+  //  Double_t y1=fP0, z1=fP1;
+  Double_t c1=fP4*x1 - fP2, r1=sqrt(1.- c1*c1);
+  Double_t c2=fP4*x2 - fP2, r2=sqrt(1.- c2*c2);
+  
+  y = fP0;
+  z = fP1;
+  y += dx*(c1+c2)/(r1+r2);
+  z += dx*(c1+c2)/(c1*r2 + c2*r1)*fP3;
+  return 0;  
+}
+
+
+//_____________________________________________________________________________
+Double_t AliTPCseed::GetPredictedChi2(const AliTPCclusterMI *c) const 
+{
+  //-----------------------------------------------------------------
+  // This function calculates a predicted chi2 increment.
+  //-----------------------------------------------------------------
+  //Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
+  Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
+  r00+=fC00; r01+=fC10; r11+=fC11;
+
+  Double_t det=r00*r11 - r01*r01;
+  if (TMath::Abs(det) < 1.e-10) {
+    Int_t n=GetNumberOfClusters();
+    if (n>4) cerr<<n<<" AliKalmanTrack warning: Singular matrix !\n";
+    return 1e10;
+  }
+  Double_t tmp=r00; r00=r11; r11=tmp; r01=-r01;
+  
+  Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
+  
+  return (dy*r00*dy + 2*r01*dy*dz + dz*r11*dz)/det;
+}
+
+
+//_________________________________________________________________________________________
+
+
+Int_t AliTPCseed::Compare(const TObject *o) const {
+  //-----------------------------------------------------------------
+  // This function compares tracks according to the sector - for given sector according z
+  //-----------------------------------------------------------------
+  AliTPCseed *t=(AliTPCseed*)o;
+  if (t->fSector>fSector) return -1;
+  if (t->fSector<fSector) return 1;
+
+  Double_t z2 = t->GetZ();
+  Double_t z1 = GetZ();
+  if (z2>z1) return 1;
+  if (z2<z1) return -1;
+  return 0;
+}
+
+
+
+
+
+
+//_____________________________________________________________________________
+Int_t AliTPCseed::Update(const AliTPCclusterMI *c, Double_t chisq, UInt_t index) {
+  //-----------------------------------------------------------------
+  // This function associates a cluster with this track.
+  //-----------------------------------------------------------------
+  //  Double_t r00=c->GetSigmaY2(), r01=0., r11=c->GetSigmaZ2();
+  //Double_t r00=sigmay2, r01=0., r11=sigmaz2;
+  Double_t r00=fErrorY2, r01=0., r11=fErrorZ2;
+
+  r00+=fC00; r01+=fC10; r11+=fC11;
+  Double_t det=r00*r11 - r01*r01;
+  Double_t tmp=r00; r00=r11/det; r11=tmp/det; r01=-r01/det;
+
+  Double_t k00=fC00*r00+fC10*r01, k01=fC00*r01+fC10*r11;
+  Double_t k10=fC10*r00+fC11*r01, k11=fC10*r01+fC11*r11;
+  Double_t k20=fC20*r00+fC21*r01, k21=fC20*r01+fC21*r11;
+  Double_t k30=fC30*r00+fC31*r01, k31=fC30*r01+fC31*r11;
+  Double_t k40=fC40*r00+fC41*r01, k41=fC40*r01+fC41*r11;
+
+  Double_t dy=c->GetY() - fP0, dz=c->GetZ() - fP1;
+  Double_t cur=fP4 + k40*dy + k41*dz, eta=fP2 + k20*dy + k21*dz;
+  if (TMath::Abs(cur*fX-eta) >= 0.99999) {
+    Int_t n=GetNumberOfClusters();
+    if (n>4) cerr<<n<<" AliTPCtrack warning: Filtering failed !\n";
+    return 0;
+  }
+
+  fP0 += k00*dy + k01*dz;
+  fP1 += k10*dy + k11*dz;
+  fP2  = eta;
+  fP3 += k30*dy + k31*dz;
+  fP4  = cur;
+
+  Double_t c01=fC10, c02=fC20, c03=fC30, c04=fC40;
+  Double_t c12=fC21, c13=fC31, c14=fC41;
+
+  fC00-=k00*fC00+k01*fC10; fC10-=k00*c01+k01*fC11;
+  fC20-=k00*c02+k01*c12;   fC30-=k00*c03+k01*c13;
+  fC40-=k00*c04+k01*c14; 
+
+  fC11-=k10*c01+k11*fC11;
+  fC21-=k10*c02+k11*c12;   fC31-=k10*c03+k11*c13;
+  fC41-=k10*c04+k11*c14; 
+
+  fC22-=k20*c02+k21*c12;   fC32-=k20*c03+k21*c13;
+  fC42-=k20*c04+k21*c14; 
+
+  fC33-=k30*c03+k31*c13;
+  fC43-=k40*c03+k41*c13; 
+
+  fC44-=k40*c04+k41*c14; 
+
+  Int_t n=GetNumberOfClusters();
+  fIndex[n]=index;
+  SetNumberOfClusters(n+1);
+  SetChi2(GetChi2()+chisq);
+
+  return 1;
+}
+
+
+
+//_____________________________________________________________________________
+Double_t AliTPCtrackerMI::f1(Double_t x1,Double_t y1,
+                   Double_t x2,Double_t y2,
+                   Double_t x3,Double_t y3) 
+{
+  //-----------------------------------------------------------------
+  // Initial approximation of the track curvature
+  //-----------------------------------------------------------------
+  Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
+  Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
+                  (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
+  Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
+                  (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
+
+  Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
+
+  return -xr*yr/sqrt(xr*xr+yr*yr); 
+}
+
+
+//_____________________________________________________________________________
+Double_t AliTPCtrackerMI::f2(Double_t x1,Double_t y1,
+                   Double_t x2,Double_t y2,
+                   Double_t x3,Double_t y3) 
+{
+  //-----------------------------------------------------------------
+  // Initial approximation of the track curvature times center of curvature
+  //-----------------------------------------------------------------
+  Double_t d=(x2-x1)*(y3-y2)-(x3-x2)*(y2-y1);
+  Double_t a=0.5*((y3-y2)*(y2*y2-y1*y1+x2*x2-x1*x1)-
+                  (y2-y1)*(y3*y3-y2*y2+x3*x3-x2*x2));
+  Double_t b=0.5*((x2-x1)*(y3*y3-y2*y2+x3*x3-x2*x2)-
+                  (x3-x2)*(y2*y2-y1*y1+x2*x2-x1*x1));
+
+  Double_t xr=TMath::Abs(d/(d*x1-a)), yr=d/(d*y1-b);
+  
+  return -a/(d*y1-b)*xr/sqrt(xr*xr+yr*yr);
+}
+
+//_____________________________________________________________________________
+Double_t AliTPCtrackerMI::f3(Double_t x1,Double_t y1, 
+                   Double_t x2,Double_t y2,
+                   Double_t z1,Double_t z2) 
+{
+  //-----------------------------------------------------------------
+  // Initial approximation of the tangent of the track dip angle
+  //-----------------------------------------------------------------
+  return (z1 - z2)/sqrt((x1-x2)*(x1-x2)+(y1-y2)*(y1-y2));
+}
+
+
+void AliTPCtrackerMI::LoadClusters()
+{
+  //
+  // load clusters to the memory
+  Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
+  for (Int_t i=0; i<j; i++) {
+    fClustersArray.LoadEntry(i);
+  }
+}
+
+void AliTPCtrackerMI::UnloadClusters()
+{
+  //
+  // load clusters to the memory
+  Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
+  for (Int_t i=0; i<j; i++) {
+    fClustersArray.ClearSegment(i);
+  }
+}
+
+
+
+//_____________________________________________________________________________
+void AliTPCtrackerMI::LoadOuterSectors() {
+  //-----------------------------------------------------------------
+  // This function fills outer TPC sectors with clusters.
+  //-----------------------------------------------------------------
+  UInt_t index;
+  Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
+  for (Int_t i=0; i<j; i++) {
+    //  AliSegmentID *s=fClustersArray.LoadEntry(i);
+    AliSegmentID *s= fClustersArray.At(i);
+    Int_t sec,row;
+    AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
+    par->AdjustSectorRow(s->GetID(),sec,row);
+    if (sec<fkNIS*2) continue;
+    AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
+    Int_t ncl=clrow->GetArray()->GetEntriesFast();
+    while (ncl--) {
+      AliTPCclusterMI *c=(AliTPCclusterMI*)(*clrow)[ncl];
+      index=(((sec<<8)+row)<<16)+ncl;
+      fOuterSec[(sec-fkNIS*2)%fkNOS][row].InsertCluster(c,index);
+    }
+  }  
+  fN=fkNOS;
+  fSectors=fOuterSec;
+}
+
+
+//_____________________________________________________________________________
+void AliTPCtrackerMI::LoadInnerSectors() {
+  //-----------------------------------------------------------------
+  // This function fills inner TPC sectors with clusters.
+  //-----------------------------------------------------------------
+  UInt_t index;
+  Int_t j=Int_t(fClustersArray.GetTree()->GetEntries());
+  for (Int_t i=0; i<j; i++) {
+    //   AliSegmentID *s=fClustersArray.LoadEntry(i);
+    AliSegmentID *s=fClustersArray.At(i);
+    Int_t sec,row;
+    AliTPCParam *par=(AliTPCParam*)fClustersArray.GetParam();
+    par->AdjustSectorRow(s->GetID(),sec,row);
+    if (sec>=fkNIS*2) continue;
+    AliTPCClustersRow *clrow=fClustersArray.GetRow(sec,row);
+    Int_t ncl=clrow->GetArray()->GetEntriesFast();
+    while (ncl--) {
+      AliTPCclusterMI *c=(AliTPCclusterMI*)(*clrow)[ncl];
+      index=(((sec<<8)+row)<<16)+ncl;
+      fInnerSec[sec%fkNIS][row].InsertCluster(c,index);
+    }
+  }
+
+  fN=fkNIS;
+  fSectors=fInnerSec;
+}
+
+Int_t AliTPCtrackerMI::FollowToNext(AliTPCseed& t, Int_t nr) {
+  //-----------------------------------------------------------------
+  // This function tries to find a track prolongation to next pad row
+  //-----------------------------------------------------------------
+  //  Double_t xt=t.GetX();
+  //  Int_t row = fSectors->GetRowNumber(xt)-1;
+  //  if (row < nr) return 1; // don't prolongate if not information until now -
+  //
+  Double_t  x=fSectors->GetX(nr), ymax=fSectors->GetMaxY(nr);
+  if (!t.PropagateTo(x)) return 0;
+  // update current
+  t.fCurrentSigmaY = GetSigmaY(&t);
+  t.fCurrentSigmaZ = GetSigmaZ(&t);
+  //  
+  AliTPCclusterMI *cl=0;
+  UInt_t index=0;
+  const AliTPCRow &krow=fSectors[t.fRelativeSector][nr];
+  Double_t sy2=ErrY2(&t)*2;
+  Double_t sz2=ErrZ2(&t)*2;
+
+
+  Double_t  roady  =3.*sqrt(t.GetSigmaY2() + sy2);
+  Double_t  roadz = 3 *sqrt(t.GetSigmaZ2() + sz2);
+  Double_t  y=t.GetY(), z=t.GetZ();
+
+  if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
+    t.fInDead = kTRUE;
+    return 0;
+  } 
+  else
+    {
+      t.fNFoundable++;
+    }   
+  //calculate 
+  Float_t maxdistance = roady*roady + roadz*roadz;
+  if (krow) {
+    for (Int_t i=krow.Find(z-roadz); i<krow; i++) {
+      AliTPCclusterMI *c=(AliTPCclusterMI*)(krow[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;
+       index=krow.GetIndex(i);       
+      }
+    }
+  }      
+  if (cl) {
+    //Double_t sy2=
+    ErrY2(&t,cl);
+    //Double_t sz2=
+    ErrZ2(&t,cl);
+    Double_t chi2= t.GetPredictedChi2(cl);    
+    UpdateTrack(&t,cl,chi2,index);   
+  } else {    
+    if (y > ymax) {
+      t.fRelativeSector= (t.fRelativeSector+1) % fN;
+      if (!t.Rotate(fSectors->GetAlpha())) 
+       return 0;
+    } else if (y <-ymax) {
+      t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
+      if (!t.Rotate(-fSectors->GetAlpha())) 
+       return 0;
+    }   
+  }
+  return 1;
+}
+
+
+Int_t AliTPCtrackerMI::UpdateClusters(AliTPCseed& t,Int_t trindex,  Int_t nr) {
+  //-----------------------------------------------------------------
+  // This function tries to find a track prolongation to next pad row
+  //-----------------------------------------------------------------
+  t.fCurrentCluster  = 0;
+  t.fCurrentClusterIndex1 = 0;   
+  t.fCurrentClusterIndex2 = 0;
+   
+  Double_t xt=t.GetX();
+  Int_t row = fSectors->GetRowNumber(xt)-1;
+  if (row < nr) return 1; // don't prolongate if not information until now -
+  Double_t x=fSectors->GetX(nr);
+  if (!t.PropagateTo(x)) return 0;
+  // update current
+  t.fCurrentSigmaY = GetSigmaY(&t);
+  t.fCurrentSigmaZ = GetSigmaZ(&t);
+    
+  AliTPCclusterMI *cl=0;
+  UInt_t index=0;
+  AliTPCRow &krow=fSectors[t.fRelativeSector][nr];
+  //
+  Double_t  y=t.GetY(), z=t.GetZ();
+  Double_t roady = 3.* TMath::Sqrt(t.GetSigmaY2() + t.fCurrentSigmaY*t.fCurrentSigmaY);
+  Double_t roadz = 3.* TMath::Sqrt(t.GetSigmaZ2() + t.fCurrentSigmaZ*t.fCurrentSigmaZ);
+  //
+
+  Float_t maxdistance = 1000000;
+  if (krow) {    
+    for (Int_t i=krow.Find(z-roadz); i<krow; i++) {
+      AliTPCclusterMI *c=(AliTPCclusterMI*)(krow[i]);
+      if (c->GetZ() > z+roadz) break;
+      if (TMath::Abs(c->GetY()-y)>roady) continue;
+            
+      //krow.UpdateClusterTrack(i,trindex,&t);
+
+      Float_t dy2 = (c->GetY()- t.GetY());
+      dy2*=dy2;
+      Float_t dz2 = (c->GetZ()- t.GetZ());
+      dz2*=dz2;
+      //
+      Float_t distance = dy2+dz2;
+      //
+      if (distance > maxdistance) continue;
+      maxdistance = distance;
+      cl=c;
+      index=i;       
+    }
+  }
+  t.fCurrentCluster  = cl;
+  t.fCurrentClusterIndex1 = krow.GetIndex(index);   
+  t.fCurrentClusterIndex2 = index;   
+  return 1;
+}
+
+
+Int_t AliTPCtrackerMI::FollowToNextCluster(Int_t trindex, Int_t nr) {
+  //-----------------------------------------------------------------
+  // This function tries to find a track prolongation to next pad row
+  //-----------------------------------------------------------------
+  AliTPCseed & t  = *((AliTPCseed*)(fSeeds->At(trindex)));
+  AliTPCRow &krow=fSectors[t.fRelativeSector][nr];
+  //  Double_t pt=t.GetConvConst()/(100/0.299792458/0.2)/t.Get1Pt();
+  Double_t y=t.GetY();
+  Double_t ymax=fSectors->GetMaxY(nr);
+
+  if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
+    t.fInDead = kTRUE;
+    return 0;
+  } 
+  else
+    {
+      t.fNFoundable++;
+    }
+  
+  if (t.fCurrentCluster) {
+    //    Float_t l=fSectors->GetPadPitchLength();
+    //    AliTPCclusterTracks * cltrack = krow.GetClusterTracks(t.fCurrentClusterIndex1);
+
+    Double_t sy2=ErrY2(&t,t.fCurrentCluster);
+    Double_t sz2=ErrZ2(&t,t.fCurrentCluster);
+
+
+    Double_t sdistancey = TMath::Sqrt(sy2+t.GetSigmaY2());
+    Double_t sdistancez = TMath::Sqrt(sz2+t.GetSigmaZ2());
+
+    Double_t rdistancey = TMath::Abs(t.fCurrentCluster->GetY()-t.GetY());
+    Double_t rdistancez = TMath::Abs(t.fCurrentCluster->GetZ()-t.GetZ());
+    
+    Double_t rdistance  = TMath::Sqrt(TMath::Power(rdistancey/sdistancey,2)+TMath::Power(rdistancez/sdistancez,2));
+
+
+    //    printf("\t%f\t%f\t%f\n",rdistancey/sdistancey,rdistancez/sdistancez,rdistance);
+    if (rdistance>4) return 0;
+
+    if ((rdistancey/sdistancey>2.5 || rdistancez/sdistancez>2.5) && t.fCurrentCluster->GetType()==0)  
+       return 0;  //suspisiouce - will be changed
+
+    if ((rdistancey/sdistancey>2. || rdistancez/sdistancez>2.0) && t.fCurrentCluster->GetType()>0)  
+       // strict cut on overlaped cluster
+       return 0;  //suspisiouce - will be changed
+
+    if ( (rdistancey/sdistancey>1. || rdistancez/sdistancez>2.5 ||t.fCurrentCluster->GetQ()<70 ) 
+        && t.fCurrentCluster->GetType()<0)
+      return 0;
+
+    //    t.SetSampledEdx(0.3*t.fCurrentCluster->GetQ()/l,t.GetNumberOfClusters(), GetSigmaY(&t), GetSigmaZ(&t));
+    UpdateTrack(&t,t.fCurrentCluster,t.GetPredictedChi2(t.fCurrentCluster),t.fCurrentClusterIndex1);
+   
+  } else {
+    if (y > ymax) {
+      t.fRelativeSector= (t.fRelativeSector+1) % fN;
+      if (!t.Rotate(fSectors->GetAlpha())) 
+       return 0;
+    } else if (y <-ymax) {
+      t.fRelativeSector= (t.fRelativeSector-1+fN) % fN;
+      if (!t.Rotate(-fSectors->GetAlpha())) 
+       return 0;
+    }
+  }
+  return 1;
+}
+
+
+
+
+/*
+Int_t AliTPCtrackerMI::FollowProlongationFast(AliTPCseed& t, Int_t step)
+{
+  //-----------------------------------------------------------------
+  // fast prolongation mathod -
+  // don't update track only after step clusters
+  //-----------------------------------------------------------------
+  Double_t xt=t.GetX();
+  //
+  Double_t alpha=t.GetAlpha(); 
+  alpha =- fSectors->GetAlphaShift();
+  if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();  
+  if (alpha < 0.            ) alpha += 2.*TMath::Pi();  
+  t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha())%fN;
+  Int_t row0 = fSectors->GetRowNumber(xt); 
+  Double_t x    = fSectors->GetX(row0);
+  Double_t ymax = fSectors->GetMaxY(row0);
+  //
+  Double_t sy2=ErrY2(&t)*2;
+  Double_t sz2=ErrZ2(&t)*2;
+  Double_t  roady  =3.*sqrt(t.GetSigmaY2() + sy2);
+  Double_t  roadz = 3 *sqrt(t.GetSigmaZ2() + sz2);
+  Float_t maxdistance = roady*roady + roadz*roadz; 
+  t.fCurrentSigmaY = GetSigmaY(&t);
+  t.fCurrentSigmaZ = GetSigmaZ(&t);
+  //
+  Int_t nclusters = 0;
+  Double_t y;
+  Double_t z;
+  Double_t yy[200];   //track prolongation
+  Double_t zz[200];
+  Double_t cy[200];  // founded cluster position
+  Double_t cz[200];
+  Double_t sy[200];  // founded cluster error
+  Double_t sz[200];
+  Bool_t   hitted[200]; // indication of cluster presence
+  //
+  
+  //
+  for (Int_t drow = step; drow>=0; drow--) {
+    Int_t row = row0-drow;
+    if (row<0) break;
+    Double_t x    = fSectors->GetX(row);
+    Double_t ymax = fSectors->GetMaxY(row);
+    t.GetProlongation(x,y,z);
+    yy[drow] =y;
+    zz[drow] =z;    
+    const AliTPCRow &krow=fSectors[t.fRelativeSector][row];
+    if (TMath::Abs(TMath::Abs(y)-ymax)<krow.fDeadZone){
+      t.fInDead = kTRUE;
+      break;
+    } 
+    else
+      {
+       t.fNFoundable++;
+      }  
+    
+    //find nearest  cluster 
+    AliTPCclusterMI *cl= 0;
+    if (krow) {
+      for (Int_t i=krow.Find(z-roadz); i<krow; i++) {
+       AliTPCclusterMI *c=(AliTPCclusterMI*)(krow[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;
+         //      index=krow.GetIndex(i);       
+       }       
+      }              
+    }  // end of seearch
+    //update cluster information
+    if (cl){ 
+      cy[drow] = cl->GetY();
+      cz[drow] = cl->GetZ();
+      sy[drow] = ErrY2(&t,cl);
+      sz[drow] = ErrZ2(&t,cl);
+      hitted[drow] = kTRUE;
+      nclusters++;
+    }
+    else
+      hitted[drow] = kFALSE;
+  }
+  //if we have information - update track
+  if (nclusters>0){
+    Float_t sumyw0   = 0;
+    Float_t sumzw0   = 0;
+    Float_t sumyw   = 0;
+    Float_t sumzw   = 0;
+    Float_t sumrow  = 0;
+    for (Int_t i=0;i<step;i++)
+      if (hitted[i]){
+       sumyw0+= 1/sy[i];
+       sumzw0+= 1/sz[i];
+       //
+       sumyw+= (cy[i]-yy[i])/sy[i];
+       sumzw+= (cz[i]-zz[i])/sz[i];
+       sumrow+= i;
+      }
+    Float_t dy = sumyw/sumyw0;
+    Float_t dz = sumzw/sumzw0;
+    Float_t mrow = sumrow/nclusters+row0;
+    Float_t x = fSectors->GetX(mrow);
+    t.PropagateTo(x);
+    AliTPCclusterMI cvirtual;
+    cvirtual.SetZ(dz+t.GetZ());
+    cvirtual.SetY(dy+t.GetY());
+    t.SetErrorY2(1.2*t.fErrorY2/TMath::Sqrt(Float_t(nclusters)));
+    t.SetErrorZ2(1.2*t.fErrorZ2/TMath::Sqrt(Float_t(nclusters)));
+    Float_t chi2 = t.GetPredictedChi2(&cvirtual);
+    t.Update(&cvirtual,chi2,0);
+    Int_t ncl = t.GetNumberOfClusters();
+    ncl = ncl-1+nclusters;
+    t.SetN(ncl);
+  }     
+  return  1;
+}   
+*/
+
+//_____________________________________________________________________________
+Int_t AliTPCtrackerMI::FollowProlongation(AliTPCseed& t, Int_t rf) {
+  //-----------------------------------------------------------------
+  // This function tries to find a track prolongation.
+  //-----------------------------------------------------------------
+  Double_t xt=t.GetX();
+  //
+  Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
+  if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();  
+  if (alpha < 0.            ) alpha += 2.*TMath::Pi();  
+  t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha())%fN;
+    
+  for (Int_t nr=fSectors->GetRowNumber(xt)-1; nr>=rf; nr--) {
+   
+    if (FollowToNext(t,nr)==0) {
+    }
+  }   
+  return 1;
+}
+
+
+//_____________________________________________________________________________
+Int_t AliTPCtrackerMI::FollowBackProlongation(AliTPCseed& t, Int_t rf) {
+  //-----------------------------------------------------------------
+  // This function tries to find a track prolongation.
+  //-----------------------------------------------------------------
+  Double_t xt=t.GetX();
+  //
+  Double_t alpha=t.GetAlpha() - fSectors->GetAlphaShift();
+  if (alpha > 2.*TMath::Pi()) alpha -= 2.*TMath::Pi();  
+  if (alpha < 0.            ) alpha += 2.*TMath::Pi();  
+  t.fRelativeSector = Int_t(alpha/fSectors->GetAlpha())%fN;
+    
+  for (Int_t nr=fSectors->GetRowNumber(xt)+1; nr<=rf; nr++) {
+    if (t.GetSnp()>0.2)
+      FollowToNext(t,nr);                                                            
+  }   
+  return 1;
+}
+
+
+
+
+   
+Float_t AliTPCtrackerMI::OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t & sum2)
+{
+  //
+  //
+  sum1=0;
+  sum2=0;
+  Int_t sum=0;
+  if (s1->fSector!=s2->fSector) return 0;
+  //
+  Float_t dz2 =(s1->GetZ() - s2->GetZ());
+  dz2*=dz2;
+  Float_t dy2 =(s1->GetY() - s2->GetY());
+  dy2*=dy2;
+  Float_t distance = TMath::Sqrt(dz2+dy2);
+  if (distance>5.) return 0; // if there are far away  - not overlap - to reduce combinatorics
+  for (Int_t i=0;i<160;i++){
+    if (s1->fClusterIndex[i]>0) sum1++;
+    if (s2->fClusterIndex[i]>0) sum2++;
+    if (s1->fClusterIndex[i]==s2->fClusterIndex[i] && s1->fClusterIndex[i]>0) {
+      sum++;
+    }
+  }
+  Float_t summin = TMath::Min(sum1,sum2);
+  Float_t ratio = sum/Float_t(summin);
+  return ratio;
+}
+
+void  AliTPCtrackerMI::SignShared(AliTPCseed * s1, AliTPCseed * s2)
+{
+  //
+  //
+  if (s1->fSector!=s2->fSector) return;
+  //
+  Float_t dz2 =(s1->GetZ() - s2->GetZ());
+  dz2*=dz2;
+  Float_t dy2 =(s1->GetY() - s2->GetY());
+  dy2*=dy2;
+  Float_t distance = TMath::Sqrt(dz2+dy2);
+  if (distance>15.) return ; // if there are far away  - not overlap - to reduce combinatorics
+  //trpoint = new (pointarray[track->fRow]) AliTPCTrackPoint;
+  //  TClonesArray &pointarray1 = *(s1->fPoints);
+  //TClonesArray &pointarray2 = *(s2->fPoints);
+  //
+  for (Int_t i=0;i<160;i++){
+    if (s1->fClusterIndex[i]==s2->fClusterIndex[i] && s1->fClusterIndex[i]>0) {
+      //  AliTPCTrackPoint *p1  = (AliTPCTrackPoint *)(pointarray1.UncheckedAt(i));
+      //AliTPCTrackPoint *p2  = (AliTPCTrackPoint *)(pointarray2.UncheckedAt(i)); 
+      AliTPCTrackPoint *p1  = s1->GetTrackPoint(i);
+      AliTPCTrackPoint *p2  = s2->GetTrackPoint(i);; 
+      p1->fIsShared = kTRUE;
+      p2->fIsShared = kTRUE;
+    }
+  } 
+}
+
+
+
+
+void  AliTPCtrackerMI::RemoveOverlap(TObjArray * arr, Float_t factor, Int_t removalindex , Bool_t shared=kFALSE){
+
+  
+
+  //
+  // remove overlap - used removal factor - removal index stored in the track
+  arr->Sort();  // sorting according z
+  arr->Expand(arr->GetEntries());
+  Int_t nseed=arr->GetEntriesFast();
+  //  printf("seeds \t%p \t%d\n",arr, nseed);
+  //  arr->Expand(arr->GetEntries());  //remove 0 pointers
+  nseed = arr->GetEntriesFast();
+  Int_t removed = 0;
+  for (Int_t i=0; i<nseed; i++) {
+    AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
+    if (!pt) {
+      continue;
+    }
+    if (!(pt->IsActive())) continue;
+    for (Int_t j=i+1; j<nseed; j++){
+      AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
+      if ((pt2) && pt2->IsActive())
+       if (pt->fSector == pt2->fSector)
+         if (TMath::Abs(pt2->GetZ()-pt->GetZ())<2){
+           Int_t sum1,sum2;
+           Float_t ratio = OverlapFactor(pt,pt2,sum1,sum2);
+           if (ratio>factor){
+             //          if (pt->GetChi2()<pt2->GetChi2()) pt2->Desactivate(removalindex);  // arr->RemoveAt(j);             
+             Float_t ratio2 = (pt->GetChi2()*sum2)/(pt2->GetChi2()*sum1);
+             Float_t ratio3 = Float_t(sum1-sum2)/Float_t(sum1+sum2);
+             removed++;
+             if (TMath::Abs(ratio3)>0.025){  // if much more points  
+               if (sum1>sum2) pt2->Desactivate(removalindex);
+               else {
+                 pt->Desactivate(removalindex); // arr->RemoveAt(i); 
+                 break;
+               }
+             }
+             else{  //decide on mean chi2
+               if (ratio2<1)  
+                 pt2->Desactivate(removalindex);
+               else {
+                 pt->Desactivate(removalindex); // arr->RemoveAt(i); 
+                 break;
+               }           
+             }  
+             
+           }  // if suspicious ratio
+         }
+         else
+           break;
+    }
+  }
+  //  printf("removed\t%d\n",removed);
+  Int_t good =0; 
+  for (Int_t i=0; i<nseed; i++) {
+    AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
+    if (!pt) break;
+    if (pt->GetNumberOfClusters() < pt->fNFoundable*0.5) {
+      //desactivate tracks with small number of points
+      //      printf("%d\t%d\t%f\n", pt->GetNumberOfClusters(), pt->fNFoundable,pt->GetNumberOfClusters()/Float_t(pt->fNFoundable));
+      pt->Desactivate(10);  //desactivate  - small muber of points
+    }
+    if (!(pt->IsActive())) continue;
+    good++;
+  }
+  
+  
+  if (shared)
+    for (Int_t i=0; i<nseed; i++) {
+      AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);    
+      if (!pt) continue;
+      if (!(pt->IsActive())) continue;
+      for (Int_t j=i+1; j<nseed; j++){
+       AliTPCseed *pt2=(AliTPCseed*)arr->UncheckedAt(j);
+       if ((pt2) && pt2->IsActive()) {
+         if ( TMath::Abs(pt->fSector-pt2->fSector)>1) break;
+         SignShared(pt,pt2);
+       }
+      }
+    }
+  fNtracks = good;
+  printf("\n*****\nNumber of good tracks after overlap removal\t%d\n",fNtracks);
+
+
+}
+void AliTPCtrackerMI::MakeSeedsAll()
+{
+  if (fSeeds == 0) fSeeds = new TObjArray;
+  TObjArray * arr;
+  for (Int_t sec=0;sec<fkNOS;sec+=3){
+     arr = MakeSeedsSectors(sec,sec+3);
+     Int_t nseed = arr->GetEntriesFast();
+     for (Int_t i=0;i<nseed;i++)
+       fSeeds->AddLast(arr->RemoveAt(i));
+  }
+   //  fSeeds = MakeSeedsSectors(0,fkNOS);
+}
+
+TObjArray *  AliTPCtrackerMI::MakeSeedsSectors(Int_t sec1, Int_t sec2)
+{
+  //
+  // loop over all  sectors and make seed
+  //find track seeds
+  TStopwatch timer;
+  Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
+  Int_t nrows=nlow+nup;  
+  Int_t gap=Int_t(0.125*nrows), shift=Int_t(0.5*gap);
+  //  if (fSeeds==0) fSeeds = new TObjArray;
+  TObjArray * arr = new TObjArray;
+
+  for (Int_t sec=sec1; sec<sec2;sec++){
+    MakeSeeds(arr, sec, nup-1, nup-1-gap);
+    MakeSeeds(arr, sec, nup-1-shift, nup-1-shift-gap);
+  }
+  Int_t nseed=arr->GetEntriesFast();
+  gap=Int_t(0.3*nrows);
+  // continue seeds 
+  Int_t i;
+  for (i=0; i<nseed; i++) {
+    AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i), &t=*pt; 
+    if (!pt) continue;
+    if (FollowProlongation(t,nup-gap)) {
+      pt->fIsSeeding =kFALSE;
+      continue;
+    }
+    delete arr->RemoveAt(i);
+  }   
+  //
+  //remove seeds which overlaps  
+  RemoveOverlap(arr,0.6,1);      
+  //delete seeds - which were sign  
+  nseed=arr->GetEntriesFast();
+  for (i=0; i<nseed; i++) {
+    AliTPCseed *pt=(AliTPCseed*)arr->UncheckedAt(i);
+    //, &t=*pt;
+    if (!pt) continue;
+    if ((pt->IsActive())  &&  pt->GetNumberOfClusters() > pt->fNFoundable*0.5 ) {
+      //FollowBackProlongation(t,nup);
+      continue;
+    }
+    delete arr->RemoveAt(i);    
+  }  
+  return arr;
+}
+
+
+
+//_____________________________________________________________________________
+void AliTPCtrackerMI::MakeSeeds(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2) {
+  //-----------------------------------------------------------------
+  // This function creates track seeds.
+  //-----------------------------------------------------------------
+  //  if (fSeeds==0) fSeeds=new TObjArray(15000);
+
+  Double_t x[5], c[15];
+
+  Double_t alpha=fOuterSec->GetAlpha(), shift=fOuterSec->GetAlphaShift();
+  Double_t cs=cos(alpha), sn=sin(alpha);
+
+  Double_t x1 =fOuterSec->GetX(i1);
+  Double_t xx2=fOuterSec->GetX(i2);
+  //
+  //  for (Int_t ns=0; ns<fkNOS; ns++) 
+  Int_t ns =sec;
+    {
+    Int_t nl=fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
+    Int_t nm=fOuterSec[ns][i2];
+    Int_t nu=fOuterSec[(ns+1)%fkNOS][i2];
+    const AliTPCRow& kr1=fOuterSec[ns][i1];
+    AliTPCRow&  kr21 = fOuterSec[(ns-1+fkNOS)%fkNOS][i2];
+    AliTPCRow&  kr22 = fOuterSec[(ns)%fkNOS][i2];
+    AliTPCRow&  kr23 = fOuterSec[(ns+1)%fkNOS][i2];
+
+    for (Int_t is=0; is < kr1; is++) {
+      Double_t y1=kr1[is]->GetY(), z1=kr1[is]->GetZ();
+      Double_t x3=GetX(), y3=GetY(), z3=GetZ();
+
+      Float_t anglez = (z1-z3)/(x1-x3); 
+      Float_t extraz = z1 - anglez*(x1-xx2);  // extrapolated z
+
+      for (Int_t js=0; js < nl+nm+nu; js++) {
+       const AliTPCclusterMI *kcl;
+        Double_t x2,   y2,   z2;
+       if (js<nl) {     
+         if (js==0) {
+           js = kr21.Find(extraz-15.);
+           if (js>=nl) continue;
+         }       
+         kcl=kr21[js];
+         z2=kcl->GetZ();
+         if ((extraz-z2)>10) continue;   
+         if ((extraz-z2)<-10) {
+           js = nl;
+           continue;
+         }
+          y2=kcl->GetY(); 
+          x2= xx2*cs+y2*sn;
+          y2=-xx2*sn+y2*cs;
+       } else 
+         if (js<nl+nm) {
+           if (js==nl) {
+             js = nl+kr22.Find(extraz-15.);
+             if (js>=nl+nm) continue;
+           }             
+           kcl=kr22[js-nl];
+           z2=kcl->GetZ();
+           if ((extraz-z2)>10) continue;
+           if ((extraz-z2)<-10) {
+             js = nl+nm;
+             continue;
+           }
+            x2=xx2; y2=kcl->GetY(); 
+         } else {
+           //const AliTPCRow& kr2=fOuterSec[(ns+1)%fkNOS][i2];   
+           if (js==nl+nm) {
+             js = nl+nm+kr23.Find(extraz-15.);
+             if (js>=nl+nm+nu) break;
+           }     
+           kcl=kr23[js-nl-nm];
+           z2=kcl->GetZ(); 
+           if ((extraz-z2)>10) continue;
+           if ((extraz-z2)<-10) {
+             break;        
+           }
+            y2=kcl->GetY();
+            x2=xx2*cs-y2*sn;
+            y2=xx2*sn+y2*cs;
+         }
+
+        Double_t zz=z1 - anglez*(x1-x2); 
+        if (TMath::Abs(zz-z2)>10.) continue;
+
+        Double_t d=(x2-x1)*(0.-y2)-(0.-x2)*(y2-y1);
+        if (d==0.) {cerr<<"MakeSeeds warning: Straight seed !\n"; continue;}
+
+       x[0]=y1;
+       x[1]=z1;
+       x[4]=f1(x1,y1,x2,y2,x3,y3);
+       if (TMath::Abs(x[4]) >= 0.0066) continue;
+       x[2]=f2(x1,y1,x2,y2,x3,y3);
+       //if (TMath::Abs(x[4]*x1-x[2]) >= 0.99999) continue;
+       x[3]=f3(x1,y1,x2,y2,z1,z2);
+       if (TMath::Abs(x[3]) > 1.2) continue;
+       Double_t a=asin(x[2]);
+       Double_t zv=z1 - x[3]/x[4]*(a+asin(x[4]*x1-x[2]));
+       if (TMath::Abs(zv-z3)>10.) continue; 
+
+        Double_t sy1=kr1[is]->GetSigmaY2()*2, sz1=kr1[is]->GetSigmaZ2()*4;
+        Double_t sy2=kcl->GetSigmaY2()*2,     sz2=kcl->GetSigmaZ2()*4;
+       //Double_t sy3=400*3./12., sy=0.1, sz=0.1;
+       Double_t sy3=25000*x[4]*x[4]+0.1, sy=0.1, sz=0.1;
+       //Double_t sy3=25000*x[4]*x[4]*60+0.5, sy=0.1, sz=0.1;
+
+       Double_t f40=(f1(x1,y1+sy,x2,y2,x3,y3)-x[4])/sy;
+       Double_t f42=(f1(x1,y1,x2,y2+sy,x3,y3)-x[4])/sy;
+       Double_t f43=(f1(x1,y1,x2,y2,x3,y3+sy)-x[4])/sy;
+       Double_t f20=(f2(x1,y1+sy,x2,y2,x3,y3)-x[2])/sy;
+       Double_t f22=(f2(x1,y1,x2,y2+sy,x3,y3)-x[2])/sy;
+       Double_t f23=(f2(x1,y1,x2,y2,x3,y3+sy)-x[2])/sy;
+       Double_t f30=(f3(x1,y1+sy,x2,y2,z1,z2)-x[3])/sy;
+       Double_t f31=(f3(x1,y1,x2,y2,z1+sz,z2)-x[3])/sz;
+       Double_t f32=(f3(x1,y1,x2,y2+sy,z1,z2)-x[3])/sy;
+       Double_t f34=(f3(x1,y1,x2,y2,z1,z2+sz)-x[3])/sz;
+
+        c[0]=sy1;
+        c[1]=0.;       c[2]=sz1;
+        c[3]=f20*sy1;  c[4]=0.;       c[5]=f20*sy1*f20+f22*sy2*f22+f23*sy3*f23;
+        c[6]=f30*sy1;  c[7]=f31*sz1;  c[8]=f30*sy1*f20+f32*sy2*f22;
+                       c[9]=f30*sy1*f30+f31*sz1*f31+f32*sy2*f32+f34*sz2*f34;
+        c[10]=f40*sy1; c[11]=0.; c[12]=f40*sy1*f20+f42*sy2*f22+f43*sy3*f23;
+        c[13]=f30*sy1*f40+f32*sy2*f42;
+        c[14]=f40*sy1*f40+f42*sy2*f42+f43*sy3*f43;
+
+        UInt_t index=kr1.GetIndex(is);
+       AliTPCseed *track=new AliTPCseed(index, x, c, x1, ns*alpha+shift);
+       track->fIsSeeding = kTRUE;
+       Int_t rc=FollowProlongation(*track, i2);
+       //FollowProlongationFast(*track, 5);
+       //FollowProlongationFast(*track, 5);
+       //FollowProlongationFast(*track, 5);
+       //FollowProlongationFast(*track, 5);
+       //Int_t rc = 1;
+       
+       track->fLastPoint = i1;  // first cluster in track position
+       if (rc==0 || track->GetNumberOfClusters()<(i1-i2)/4 || track->GetNumberOfClusters() < track->fNFoundable/2. ) delete track;
+        else arr->AddLast(track); 
+      }
+    }
+  }
+}
+
+//_____________________________________________________________________________
+Int_t AliTPCtrackerMI::ReadSeeds(const TFile *inp) {
+  //-----------------------------------------------------------------
+  // This function reades track seeds.
+  //-----------------------------------------------------------------
+  TDirectory *savedir=gDirectory; 
+
+  TFile *in=(TFile*)inp;
+  if (!in->IsOpen()) {
+     cerr<<"AliTPCtrackerMI::ReadSeeds(): input file is not open !\n";
+     return 1;
+  }
+
+  in->cd();
+  TTree *seedTree=(TTree*)in->Get("Seeds");
+  if (!seedTree) {
+     cerr<<"AliTPCtrackerMI::ReadSeeds(): ";
+     cerr<<"can't get a tree with track seeds !\n";
+     return 2;
+  }
+  AliTPCtrack *seed=new AliTPCtrack; 
+  seedTree->SetBranchAddress("tracks",&seed);
+  
+  if (fSeeds==0) fSeeds=new TObjArray(15000);
+
+  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()));
+  }
+  
+  delete seed;
+  delete seedTree; 
+  savedir->cd();
+  return 0;
+}
+
+//_____________________________________________________________________________
+Int_t AliTPCtrackerMI::Clusters2Tracks(const TFile *inp, TFile *out) {
+  //-----------------------------------------------------------------
+  // This is a track finder.
+  //-----------------------------------------------------------------
+  TDirectory *savedir=gDirectory; 
+
+  if (inp) {
+     TFile *in=(TFile*)inp;
+     if (!in->IsOpen()) {
+        cerr<<"AliTPCtrackerMI::Clusters2Tracks(): input file is not open !\n";
+        return 1;
+     }
+  }
+
+  if (!out->IsOpen()) {
+     cerr<<"AliTPCtrackerMI::Clusters2Tracks(): output file is not open !\n";
+     return 2;
+  }
+
+  out->cd();
+
+  char   tname[100];
+  sprintf(tname,"TreeT_TPC_%d",fEventN);
+  TTree tracktree(tname,"Tree with TPC tracks");
+  TTree seedtree("Seeds","Seeds");
+  AliTPCtrack *iotrack=0;
+  AliTPCseed  *ioseed=0;
+  tracktree.Branch("tracks","AliTPCtrack",&iotrack,32000,0);
+  TStopwatch timer;
+  
+  printf("Loading clusters \n");
+  LoadClusters();
+  printf("Time for loading clusters: \t");timer.Print();timer.Start();
+
+  printf("Loading outer sectors\n");
+  LoadOuterSectors();
+  printf("Time for loading outer sectors: \t");timer.Print();timer.Start();
+
+  printf("Loading inner sectors\n");
+  LoadInnerSectors();
+  printf("Time for loading inner sectors: \t");timer.Print();timer.Start();
+  fSectors = fOuterSec;
+  fN=fkNOS;
+  
+  
+  //find track seeds
+  MakeSeedsAll(); 
+  printf("Time for seeding: \t"); timer.Print();timer.Start();
+  Int_t nup=fOuterSec->GetNRows(), nlow=fInnerSec->GetNRows();
+  Int_t nrows=nlow+nup;
+  
+  Int_t gap=Int_t(0.3*nrows);
+  Int_t i;
+  Int_t nseed=fSeeds->GetEntriesFast();
+  // outer sectors parallel tracking
+  ParallelTracking(fSectors->GetNRows()-gap-1,0); 
+  printf("Time for parralel tracking outer sectors: \t"); timer.Print();timer.Start();
+
+  RemoveOverlap(fSeeds, 0.6,3);         
+  printf("Time for removal overlap- outer sectors: \t");timer.Print();timer.Start();
+  //parallel tracking 
+  fSectors = fInnerSec;
+  fN=fkNIS;
+  ParallelTracking(fSectors->GetNRows()-1,0);
+  printf("Number of tracks after  inner tracking  %d\n",fNtracks); 
+  printf("Time for parralel tracking inner sectors: \t"); timer.Print();timer.Start();
+  //
+  RemoveOverlap(fSeeds,0.6,5,kTRUE);  // remove overlap -  shared points signed 
+  printf("Time for removal overlap- inner sectors: \t"); timer.Print();timer.Start();
+  //
+  // 
+  
+  ioseed  = (AliTPCseed*)(fSeeds->UncheckedAt(0));
+  AliTPCseed * vseed = new AliTPCseed;
+  vseed->fPoints = new TClonesArray("AliTPCTrackPoint",1);
+  vseed->fEPoints = new TClonesArray("AliTPCExactPoint",1);
+  vseed->fPoints->ExpandCreateFast(2);
+  
+  TBranch * seedbranch =   seedtree.Branch("seeds","AliTPCseed",&vseed,32000,99);
+  //delete vseed;
+  nseed=fSeeds->GetEntriesFast();
+
+  Int_t found = 0;
+  for (i=0; i<nseed; i++) {
+    AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;    
+    if (!pt) continue;    
+    Int_t nc=t.GetNumberOfClusters();
+    t.CookdEdx(0.02,0.6);
+    CookLabel(pt,0.1); //For comparison only
+    //   if ((pt->IsActive()) && (nc>Int_t(0.4*nrows))){
+    if ((pt->IsActive()) && (nc>Int_t(0.5*t.fNFoundable) && (nc>Int_t(0.3*nrows)))){
+      iotrack=pt;
+      tracktree.Fill();
+     cerr<<found++<<'\r';      
+    }   
+    else 
+      if ( (pt->IsActive())) fNtracks--;
+    pt->RebuildSeed();
+    seedbranch->SetAddress(&pt);
+
+    seedtree.Fill();        
+    for (Int_t j=0;j<160;j++){
+      delete pt->fPoints->RemoveAt(j);
+    }
+    delete pt->fPoints;
+    pt->fPoints =0;
+    delete fSeeds->RemoveAt(i);
+  }
+  printf("Time for track writing and dedx cooking: \t"); timer.Print();timer.Start();
+
+  UnloadClusters();
+  printf("Time for unloading cluster: \t"); timer.Print();timer.Start();
+
+  tracktree.Write();
+  seedtree.Write();
+  cerr<<"Number of found tracks : "<<fNtracks<<"\t"<<found<<endl;
+  
+  savedir->cd();
+  
+  return 0;
+}
+
+
+void  AliTPCtrackerMI::ParallelTracking(Int_t rfirst, Int_t rlast)
+{
+  //
+  // try to track in parralel
+
+  Int_t nseed=fSeeds->GetEntriesFast();
+  //prepare seeds for tracking
+  for (Int_t i=0; i<nseed; i++) {
+    AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt; 
+    if (!pt) continue;
+    if (!t.IsActive()) continue;
+    // follow prolongation to the first layer
+    FollowProlongation(t, rfirst+1);
+  }
+
+
+  //
+  for (Int_t nr=rfirst; nr>=rlast; nr--){      
+    // make indexes with the cluster tracks for given       
+    //    for (Int_t i = 0;i<fN;i++)
+    //  fSectors[i][nr].MakeClusterTracks();
+
+    // find nearest cluster
+    for (Int_t i=0; i<nseed; i++) {
+      AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt; 
+      if (!pt) continue;
+      if (!pt->IsActive()) continue;
+      UpdateClusters(t,i,nr);
+    }
+    // prolonagate to the nearest cluster - if founded
+    for (Int_t i=0; i<nseed; i++) {
+      AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i); 
+      if (!pt) continue;
+      if (!pt->IsActive()) continue;
+      FollowToNextCluster(i,nr);
+    }
+    //    for (Int_t i= 0;i<fN;i++)
+    //  fSectors[i][nr].ClearClusterTracks();
+  }  
+  
+}
+
+Float_t  AliTPCtrackerMI::GetSigmaY(AliTPCseed * seed)
+{
+  //
+  //  
+  Float_t sd2 = (fParam->GetZLength()-TMath::Abs(seed->GetZ()))*fParam->GetDiffL()*fParam->GetDiffL();
+  Float_t padlength =  fParam->GetPadPitchLength(seed->fSector);
+  Float_t sres = (seed->fSector < fParam->GetNSector()/2) ? 0.2 :0.3;
+  Float_t angular  = seed->GetSnp();
+  angular = angular*angular/(1-angular*angular);
+  //  angular*=angular;
+  //angular  = TMath::Sqrt(angular/(1-angular));
+  Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular/12.+sres*sres);
+  return res;
+}
+Float_t  AliTPCtrackerMI::GetSigmaZ(AliTPCseed * seed)
+{
+  //
+  //
+  Float_t sd2 = (fParam->GetZLength()-TMath::Abs(seed->GetZ()))*fParam->GetDiffL()*fParam->GetDiffL();
+  Float_t padlength =  fParam->GetPadPitchLength(seed->fSector);
+  Float_t sres = fParam->GetZSigma();
+  Float_t angular  = seed->GetTgl();
+  Float_t res = TMath::Sqrt(sd2+padlength*padlength*angular*angular/12.+sres*sres);
+  return res;
+}
+
+
+
+//_________________________________________________________________________
+AliTPCclusterMI *AliTPCtrackerMI::GetClusterMI(Int_t index) const {
+  //--------------------------------------------------------------------
+  //       Return pointer to a given cluster
+  //--------------------------------------------------------------------
+  Int_t sec=(index&0xff000000)>>24; 
+  Int_t row=(index&0x00ff0000)>>16; 
+  Int_t ncl=(index&0x0000ffff)>>00;
+
+  AliTPCClustersRow *clrow=((AliTPCtrackerMI *) this)->fClustersArray.GetRow(sec,row);
+  return (AliTPCclusterMI*)(*clrow)[ncl];      
+}
+
+//__________________________________________________________________________
+void AliTPCtrackerMI::CookLabel(AliKalmanTrack *t, Float_t wrong) const {
+  //--------------------------------------------------------------------
+  //This function "cooks" a track label. If label<0, this track is fake.
+  //--------------------------------------------------------------------
+  Int_t noc=t->GetNumberOfClusters();
+  Int_t *lb=new Int_t[noc];
+  Int_t *mx=new Int_t[noc];
+  AliTPCclusterMI **clusters=new AliTPCclusterMI*[noc];
+
+  Int_t i;
+  for (i=0; i<noc; i++) {
+     lb[i]=mx[i]=0;
+     Int_t index=t->GetClusterIndex(i);
+     clusters[i]=GetClusterMI(index);
+  }
+
+  Int_t lab=123456789;
+  for (i=0; i<noc; i++) {
+    AliTPCclusterMI *c=clusters[i];
+    lab=TMath::Abs(c->GetLabel(0));
+    Int_t j;
+    for (j=0; j<noc; j++) if (lb[j]==lab || mx[j]==0) break;
+    lb[j]=lab;
+    (mx[j])++;
+  }
+
+  Int_t max=0;
+  for (i=0; i<noc; i++) if (mx[i]>max) {max=mx[i]; lab=lb[i];}
+    
+  for (i=0; i<noc; i++) {
+    AliTPCclusterMI *c=clusters[i];
+    if (TMath::Abs(c->GetLabel(1)) == lab ||
+        TMath::Abs(c->GetLabel(2)) == lab ) max++;
+  }
+
+  if ((1.- Float_t(max)/noc) > wrong) lab=-lab;
+
+  else {
+     Int_t tail=Int_t(0.10*noc);
+     max=0;
+     for (i=1; i<=tail; i++) {
+       AliTPCclusterMI *c=clusters[noc-i];
+       if (lab == TMath::Abs(c->GetLabel(0)) ||
+           lab == TMath::Abs(c->GetLabel(1)) ||
+           lab == TMath::Abs(c->GetLabel(2))) max++;
+     }
+     if (max < Int_t(0.5*tail)) lab=-lab;
+  }
+
+  t->SetLabel(lab);
+
+  delete[] lb;
+  delete[] mx;
+  delete[] clusters;
+}
+
+//_________________________________________________________________________
+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();
+     fRow=new AliTPCRow[fN];
+     for (Int_t i=0; i<fN; i++) {
+       fRow[i].SetX(par->GetPadRowRadiiLow(i));
+       fRow[i].fDeadZone =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();
+     fRow=new AliTPCRow[fN];
+     for (Int_t i=0; i<fN; i++) {
+       fRow[i].SetX(par->GetPadRowRadiiUp(i)); 
+       fRow[i].fDeadZone =1.5;  // 1.5 cm of dead zone
+     }
+  } 
+}
+
+
+AliTPCtrackerMI::AliTPCRow::~AliTPCRow(){
+  //
+  if (fClusterTracks) delete [] fClusterTracks;
+  fClusterTracks = 0;
+}
+
+void AliTPCtrackerMI::AliTPCRow::MakeClusterTracks(){
+  //create cluster tracks
+  if (fN>0) 
+    fClusterTracks = new AliTPCclusterTracks[fN];
+}
+
+void AliTPCtrackerMI::AliTPCRow::ClearClusterTracks(){
+  if (fClusterTracks) delete[] fClusterTracks;
+  fClusterTracks =0;
+}
+
+
+
+void AliTPCtrackerMI::AliTPCRow::UpdateClusterTrack(Int_t clindex, Int_t trindex, AliTPCseed * seed){
+  //
+  //
+  // update information of the cluster tracks - if track is nearer then other tracks to the 
+  // given track
+  const AliTPCclusterMI * cl = (*this)[clindex];
+  AliTPCclusterTracks * cltracks = GetClusterTracks(clindex);
+  // find the distance of the cluster to the track
+  Float_t dy2 = (cl->GetY()- seed->GetY());
+  dy2*=dy2;
+  Float_t dz2 = (cl->GetZ()- seed->GetZ());
+  dz2*=dz2;
+  //
+  Float_t distance = TMath::Sqrt(dy2+dz2);
+  if (distance> 3) 
+    return;  // MI - to be changed - AliTPCtrackerParam
+  
+  if ( distance < cltracks->fDistance[0]){
+    cltracks->fDistance[2] =cltracks->fDistance[1];
+    cltracks->fDistance[1] =cltracks->fDistance[0];
+    cltracks->fDistance[0] =distance;
+    cltracks->fTrackIndex[2] =cltracks->fTrackIndex[1];
+    cltracks->fTrackIndex[1] =cltracks->fTrackIndex[0];
+    cltracks->fTrackIndex[0] =trindex; 
+  }
+  else
+    if ( distance < cltracks->fDistance[1]){
+      cltracks->fDistance[2] =cltracks->fDistance[1];  
+      cltracks->fDistance[1] =distance;
+      cltracks->fTrackIndex[2] =cltracks->fTrackIndex[1];
+      cltracks->fTrackIndex[1] =trindex; 
+    } else
+      if (distance < cltracks->fDistance[2]){
+       cltracks->fDistance[2] =distance;
+       cltracks->fTrackIndex[2] =trindex;
+      }  
+}
+
+
+//_________________________________________________________________________
+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) {
+    cerr<<"AliTPCRow::InsertCluster(): Too many clusters !\n"; return;
+  }
+  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++;
+}
+
+//___________________________________________________________________
+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;
+}
+
+
+
+AliTPCseed::AliTPCseed():AliTPCtrack(){
+  //
+  fRow=0; 
+  fRemoval =0; 
+  memset(fClusterIndex,0,sizeof(Int_t)*200);
+  fPoints = 0;
+  fEPoints = 0;
+  fNFoundable =0;
+  fNShared  =0;
+  fTrackPoints =0;
+  fRemoval = 0;
+}
+
+AliTPCseed::AliTPCseed(const AliTPCtrack &t):AliTPCtrack(t){
+  fPoints = 0;
+  fEPoints = 0;
+  fNShared  =0; 
+  fTrackPoints =0;
+  fRemoval =0;
+}
+
+AliTPCseed::AliTPCseed(const AliKalmanTrack &t, Double_t a):AliTPCtrack(t,a){
+  fRow=0;
+  memset(fClusterIndex,0,sizeof(Int_t)*200); 
+  fPoints = 0;
+  fEPoints = 0;
+  fNFoundable =0; 
+  fNShared  =0; 
+  fTrackPoints =0;
+  fRemoval =0;
+}
+
+AliTPCseed::AliTPCseed(UInt_t index, const Double_t xx[5], const Double_t cc[15], 
+                                       Double_t xr, Double_t alpha):      
+  AliTPCtrack(index, xx, cc, xr, alpha) {
+  //
+  //
+  fRow =0;
+  memset(fClusterIndex,0,sizeof(Int_t)*200); 
+  fPoints = 0;
+  fEPoints = 0;
+  fNFoundable =0;
+  fNShared  = 0;
+  fTrackPoints =0;
+  fRemoval =0;
+}
+
+AliTPCseed::~AliTPCseed(){
+  if (fPoints) delete fPoints;
+  fPoints =0;
+  fEPoints = 0;
+  if (fTrackPoints){
+    for (Int_t i=0;i<8;i++){
+      delete [] fTrackPoints[i];
+    }
+    delete fTrackPoints;
+    fTrackPoints =0;
+  }
+
+}
+
+AliTPCTrackPoint * AliTPCseed::GetTrackPoint(Int_t i)
+{
+  //
+  // 
+  if (!fTrackPoints) {
+    fTrackPoints = new AliTPCTrackPoint*[8];
+    for ( Int_t i=0;i<8;i++)
+      fTrackPoints[i]=0;
+  }
+  Int_t index1 = i/20;
+  if (!fTrackPoints[index1]) fTrackPoints[index1] = new AliTPCTrackPoint[20];
+  return &(fTrackPoints[index1][i%20]);
+}
+
+void AliTPCseed::RebuildSeed()
+{
+  //
+  // rebuild seed to be ready for storing
+  fPoints = new TClonesArray("AliTPCTrackPoint",160);
+  fPoints->ExpandCreateFast(160);
+  fEPoints = new TClonesArray("AliTPCExactPoint",1);
+  for (Int_t i=0;i<160;i++){
+    AliTPCTrackPoint *trpoint = (AliTPCTrackPoint*)fPoints->UncheckedAt(i);
+    *trpoint = *(GetTrackPoint(i));
+  }
+
+}
+
+//_____________________________________________________________________________
+void AliTPCseed::CookdEdx(Double_t low, Double_t up) {
+  //-----------------------------------------------------------------
+  // This funtion calculates dE/dX within the "low" and "up" cuts.
+  //-----------------------------------------------------------------
+
+  Float_t amp[200];
+  Float_t angular[200];
+  Float_t weight[200];
+  Int_t index[200];
+  //Int_t nc = 0;
+  //  TClonesArray & arr = *fPoints; 
+  Float_t meanlog = 100.;
+  
+  Float_t mean[4]  = {0,0,0,0};
+  Float_t sigma[4] = {1000,1000,1000,1000};
+  Int_t nc[4]      = {0,0,0,0};
+  Float_t norm[4]    = {1000,1000,1000,1000};
+  //
+  //
+  fNShared =0;
+
+  for (Int_t of =0; of<4; of++){    
+    for (Int_t i=of;i<160;i+=4)
+      {
+       //AliTPCTrackPoint * point = (AliTPCTrackPoint *) arr.At(i);
+       AliTPCTrackPoint * point = GetTrackPoint(i);
+       if (point==0) continue;
+       if (point->fIsShared){
+         fNShared++;
+         continue;
+       }
+       if (point->GetCPoint().GetMax()<5) continue;
+       Float_t angley = point->GetTPoint().GetAngleY();
+       Float_t anglez = point->GetTPoint().GetAngleZ();
+       Int_t   type   = point->GetCPoint().GetType();
+       Float_t rsigmay =  point->GetCPoint().GetSigmaY();
+       Float_t rsigmaz =  point->GetCPoint().GetSigmaZ();
+       Float_t rsigma = TMath::Sqrt(rsigmay*rsigmaz);
+
+       Float_t ampc   = 0;     // normalization to the number of electrons
+       if (i>64){
+         ampc = 1.*point->GetCPoint().GetMax();
+         //ampc = 1.*point->GetCPoint().GetQ();          
+         //      AliTPCClusterPoint & p = point->GetCPoint();
+         //      Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.6)) - TMath::Abs(p.GetY()/0.6)+0.5);
+         // Float_t iz =  (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
+         //Float_t dz = 
+         //  TMath::Abs( Int_t(iz) - iz + 0.5);
+         //ampc *= 1.15*(1-0.3*dy);
+         //ampc *= 1.15*(1-0.3*dz);
+
+       }
+       else{ 
+         ampc = 1.0*point->GetCPoint().GetMax(); 
+         //ampc = 1.0*point->GetCPoint().GetQ(); 
+         //AliTPCClusterPoint & p = point->GetCPoint();
+         // Float_t dy = TMath::Abs(Int_t( TMath::Abs(p.GetY()/0.4)) - TMath::Abs(p.GetY()/0.4)+0.5);
+         //Float_t iz =  (250.0-TMath::Abs(p.GetZ())+0.11)/0.566;
+         //Float_t dz = 
+         //  TMath::Abs( Int_t(iz) - iz + 0.5);
+
+         //ampc *= 1.15*(1-0.3*dy);
+         //ampc *= 1.15*(1-0.3*dz);
+
+       }
+       ampc *= 2.0;     // put mean value to channel 50
+       //ampc *= 0.58;     // put mean value to channel 50
+       Float_t w      =  1.;
+       //      if (type>0)  w =  1./(type/2.-0.5); 
+       if (i<64) {
+         ampc /= 0.6;
+       }
+       if (i>128){
+         ampc /=1.5;
+       }
+       if (type<0) {  //amp at the border - lower weight
+         // w*= 2.;
+         continue;
+       }
+       if (rsigma>1.5) ampc/=1.3;  // if big backround
+       amp[nc[of]]        = ampc;
+       angular[nc[of]]    = TMath::Sqrt(1.+angley*angley+anglez*anglez);
+       weight[nc[of]]     = w;
+       nc[of]++;
+      }
+    
+    TMath::Sort(nc[of],amp,index,kFALSE);
+    Float_t sumamp=0;
+    Float_t sumamp2=0;
+    Float_t sumw=0;
+    //meanlog = amp[index[Int_t(nc[of]*0.33)]];
+    meanlog = 200;
+    for (Int_t i=int(nc[of]*low+0.5);i<int(nc[of]*up+0.5);i++){
+      Float_t ampl      = amp[index[i]]/angular[index[i]];
+      ampl              = meanlog*TMath::Log(1.+ampl/meanlog);
+      //
+      sumw    += weight[index[i]]; 
+      sumamp  += weight[index[i]]*ampl;
+      sumamp2 += weight[index[i]]*ampl*ampl;
+      norm[of]    += angular[index[i]]*weight[index[i]];
+    }
+    if (sumw<1){ 
+      SetdEdx(0);  
+    }
+    else {
+      norm[of] /= sumw;
+      mean[of]  = sumamp/sumw;
+      sigma[of] = sumamp2/sumw-mean[of]*mean[of];
+      if (sigma[of]>0.1) 
+       sigma[of] = TMath::Sqrt(sigma[of]);
+      else
+       sigma[of] = 1000;
+      
+    mean[of] = (TMath::Exp(mean[of]/meanlog)-1)*meanlog;
+    //mean  *=(1-0.02*(sigma/(mean*0.17)-1.));
+    //mean *=(1-0.1*(norm-1.));
+    }
+  }
+
+  Float_t dedx =0;
+  fSdEdx =0;
+  fMAngular =0;
+  //  mean[0]*= (1-0.05*(sigma[0]/(0.01+mean[1]*0.18)-1));
+  //  mean[1]*= (1-0.05*(sigma[1]/(0.01+mean[0]*0.18)-1));
+
+  
+  //  dedx = (mean[0]* TMath::Sqrt((1.+nc[0]))+ mean[1]* TMath::Sqrt((1.+nc[1])) )/ 
+  //  (  TMath::Sqrt((1.+nc[0]))+TMath::Sqrt((1.+nc[1])));
+
+  Int_t norm2 = 0;
+  Int_t norm3 = 0;
+  for (Int_t i =0;i<4;i++){
+    if (nc[i]>2&&nc[i]<1000){
+      dedx      += mean[i] *nc[i];
+      fSdEdx    += sigma[i]*(nc[i]-2);
+      fMAngular += norm[i] *nc[i];    
+      norm2     += nc[i];
+      norm3     += nc[i]-2;
+    }
+    fDEDX[i]  = mean[i];             
+    fSDEDX[i] = sigma[i];            
+    fNCDEDX[i]= nc[i]; 
+  }
+
+  if (norm3>0){
+    dedx   /=norm2;
+    fSdEdx /=norm3;
+    fMAngular/=norm2;
+  }
+  else{
+    SetdEdx(0);
+    return;
+  }
+  //  Float_t dedx1 =dedx;
+  
+  dedx =0;
+  for (Int_t i =0;i<4;i++){
+    if (nc[i]>2&&nc[i]<1000){
+      mean[i]   = mean[i]*(1-0.12*(sigma[i]/(fSdEdx)-1.));
+      dedx      += mean[i] *nc[i];
+    }
+    fDEDX[i]  = mean[i];                
+  }
+  dedx /= norm2;
+  
+
+  
+  SetdEdx(dedx);
+    
+  //mi deDX
+
+
+
+  //Very rough PID
+  Double_t p=TMath::Sqrt((1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt()));
+
+  if (p<0.6) {
+    if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
+    if (dedx < 39.+ 12./p/p) { SetMass(0.49368); return;}
+    SetMass(0.93827); return;
+  }
+
+  if (p<1.2) {
+    if (dedx < 39.+ 12./(p+0.25)/(p+0.25)) { SetMass(0.13957); return;}
+    SetMass(0.93827); return;
+  }
+
+  SetMass(0.13957); return;
+
+}
+
+
diff --git a/TPC/AliTPCtrackerMI.h b/TPC/AliTPCtrackerMI.h
new file mode 100644 (file)
index 0000000..eba829b
--- /dev/null
@@ -0,0 +1,251 @@
+#ifndef ALITPCTRACKERMI_H
+#define ALITPCTRACKERMI_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+/* $Id$ */
+
+//-------------------------------------------------------
+//                       TPC trackerMI
+//
+//   Origin: 
+//-------------------------------------------------------
+#include "AliTracker.h"
+#include "AliTPCtrack.h"
+#include "AliTPCClustersArray.h"
+
+#include "AliTPCreco.h" 
+#include "Rtypes.h"
+
+class TFile;
+class AliTPCParam;
+class AliTPCseed;
+class AliTPCclusterMI;
+class AliTPCTrackPoint;
+
+
+
+class AliTPCclusterTracks {
+ public: 
+  AliTPCclusterTracks();
+  Float_t fDistance[3];   // distance to the 3 nerest track if there overlap with cluster
+  Short_t fTrackIndex[3]; // indexes of the  tracks overlapped with clusters 
+};
+
+
+class AliTPCseed : public AliTPCtrack {
+   public:
+     AliTPCseed();
+     virtual ~AliTPCseed();
+     AliTPCseed(const AliTPCtrack &t);
+     AliTPCseed(const AliKalmanTrack &t, Double_t a);
+     Int_t Compare(const TObject *o) const;
+
+     Int_t GetProlongation(Double_t xr, Double_t &y, Double_t & z) const;
+     virtual Double_t GetPredictedChi2(const AliTPCclusterMI *cluster) const;
+     virtual Int_t Update(const AliTPCclusterMI* c, Double_t chi2, UInt_t i);
+     AliTPCTrackPoint * GetTrackPoint(Int_t i);
+     void RebuildSeed(); // rebuild seed to be ready for storing
+     AliTPCseed(UInt_t index, const Double_t xx[5], 
+                const Double_t cc[15], Double_t xr, Double_t alpha);
+     void SetClusterIndex(Int_t index){
+       fClusterIndex[fRow] = index;
+     }
+     void SetErrorY2(Float_t sy2){fErrorY2=sy2;}
+     void SetErrorZ2(Float_t sz2){fErrorZ2=sz2;}
+     void CookdEdx(Double_t low=0.05, Double_t up=0.70);
+     Bool_t IsActive(){ return !(fRemoval);}
+     void Desactivate(Int_t reason){ fRemoval = reason;} 
+     Int_t fRelativeSector;  // ! index of current relative sector
+     Int_t   fClusterIndex[200];  //array of cluster indexes
+    
+     Int_t fRemoval;               //reason - why was track removed - 0 - means still active     
+     TClonesArray * fPoints;              // array with points along the track   
+     TClonesArray * fEPoints;             // array with exact points - calculated in special macro not used in tracking
+     Int_t fRow;                 //! current row number  
+     Int_t fSector;              //!current sector number
+     Float_t fCurrentSigmaY;     //!expected current cluster sigma Y
+     Float_t fCurrentSigmaZ;     //!expected current cluster sigma Z
+     AliTPCclusterMI * fCurrentCluster; //!pointer to the current cluster for prolongation
+     Int_t   fCurrentClusterIndex1; //! index of the current cluster
+     Int_t   fCurrentClusterIndex2; //! index of the current cluster
+    
+     Float_t fErrorY2;   //!sigma of current cluster 
+     Float_t fErrorZ2;   //!sigma of current cluster    
+     Int_t   fNFoundable;      //number of foundable clusters - dead zone taken to the account
+     Bool_t  fInDead;         // indicate if the track is in dead zone
+     Int_t   fFirstPoint;    // first cluster position
+     Int_t   fLastPoint;     // last  cluster position     
+     Int_t   fNShared;       // number of shared points
+     Bool_t  fIsSeeding;     //indicates if it is proces of seeading
+   private:
+     Float_t fSdEdx;           // sigma of dedx
+     Float_t fMAngular;        // mean angular factor
+     AliTPCTrackPoint   ** fTrackPoints;  //!track points - array track points
+     Float_t fDEDX[4];         // dedx according padrows
+     Float_t fSDEDX[4];        // sdedx according padrows
+     Int_t   fNCDEDX[4];       // number of clusters for dedx measurment
+   
+     ClassDef(AliTPCseed,1)  
+};
+
+
+
+
+class AliTPCtrackerMI : public AliTracker {
+public:
+   AliTPCtrackerMI():AliTracker(),fkNIS(0),fkNOS(0) {
+      fInnerSec=fOuterSec=0; fSeeds=0; 
+   }
+   AliTPCtrackerMI(const AliTPCParam *par, Int_t eventn=0);
+  ~AliTPCtrackerMI();
+
+   Int_t ReadSeeds(const TFile *in);
+   void LoadClusters();
+   void UnloadClusters();
+
+   void LoadInnerSectors();
+   void LoadOuterSectors();
+   AliCluster * GetCluster (int) const {return 0;}
+   AliTPCclusterMI *GetClusterMI(Int_t index) const;
+   Int_t Clusters2Tracks(const TFile *in, TFile *out);
+   //   Int_t PropagateBack(const TFile *in, TFile *out);
+
+   virtual void  CookLabel(AliKalmanTrack *t,Float_t wrong) const; 
+
+   virtual Double_t ErrY2(AliTPCseed* seed, AliTPCclusterMI * cl = 0);
+   virtual Double_t ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl = 0);   
+
+   Double_t f1(Double_t x1,Double_t y1, Double_t x2,Double_t y2, Double_t x3,Double_t y3); 
+   Double_t f2(Double_t x1,Double_t y1, Double_t x2,Double_t y2, Double_t x3,Double_t y3); 
+   Double_t f3(Double_t x1,Double_t y1, Double_t x2,Double_t y2, Double_t z1,Double_t z2); 
+public:
+//**************** Internal tracker class ********************** 
+   class AliTPCRow {
+   public:
+     AliTPCRow() {fN=0; fClusterTracks=0;}
+     ~AliTPCRow();
+     void InsertCluster(const AliTPCclusterMI *c, UInt_t index);
+     operator int() const {return fN;}
+     const AliTPCclusterMI* operator[](Int_t i) const {return fClusters[i];}
+     UInt_t GetIndex(Int_t i) const {return fIndex[i];}
+     Int_t Find(Double_t z) const; 
+     void SetX(Double_t x) {fX=x;}
+     Double_t GetX() const {return fX;}
+     AliTPCclusterTracks *  GetClusterTracks(Int_t index){ return ( (index<fN) && fClusterTracks!=0)? &(fClusterTracks[index]):0;}
+     void UpdateClusterTrack(Int_t clindex, Int_t trindex,AliTPCseed * seed); 
+     void MakeClusterTracks();
+     void ClearClusterTracks();
+     Float_t fDeadZone;  // the width of the dead zone
+   private:
+     Int_t fN;                                          //number of clusters 
+     const AliTPCclusterMI *fClusters[kMaxClusterPerRow]; //pointers to clusters
+     UInt_t fIndex[kMaxClusterPerRow];                  //indeces of clusters
+     Double_t fX;                                 //X-coordinate of this row
+     AliTPCclusterTracks * fClusterTracks;        // array of cluster tracks - for overlap calculation
+   private:
+     AliTPCRow(const AliTPCRow& r);            //dummy copy constructor
+     AliTPCRow &operator=(const AliTPCRow& r); //dummy assignment operator
+   };
+
+//**************** Internal tracker class ********************** 
+   class AliTPCSector {
+   public:
+     AliTPCSector() { fN=0; fRow = 0; }
+    ~AliTPCSector() { delete[] fRow; }
+     AliTPCRow& operator[](Int_t i) const { return *(fRow+i); }
+     Int_t GetNRows() const { return fN; }
+     void Setup(const AliTPCParam *par, Int_t flag);
+     Double_t GetX(Int_t l) const {return fRow[l].GetX();}
+     Double_t GetMaxY(Int_t l) const {
+         return GetX(l)*TMath::Tan(0.5*GetAlpha());
+     } 
+     Double_t GetAlpha() const {return fAlpha;}
+     Double_t GetAlphaShift() const {return fAlphaShift;}     
+     Int_t 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);} 
+       }
+     }
+     Double_t GetPadPitchWidth()  const {return fPadPitchWidth;}
+     Double_t GetPadPitchLength() const {return fPadPitchLength;}
+     Double_t GetPadPitchLength(Float_t x) const {return (x<200) ? fPadPitchLength:f2PadPitchLength ;}
+
+   private:
+     Int_t fN;                        //number of pad rows 
+     AliTPCRow *fRow;                    //array of pad rows
+     Double_t fAlpha;                    //opening angle
+     Double_t fAlphaShift;               //shift angle;
+     Double_t fPadPitchWidth;            //pad pitch width
+     Double_t fPadPitchLength;           //pad pitch length
+     Double_t f1PadPitchLength;           //pad pitch length
+     Double_t f2PadPitchLength;           //pad pitch length
+    
+   private:
+     AliTPCSector(const AliTPCSector &s);           //dummy copy contructor
+     AliTPCSector& operator=(const AliTPCSector &s);//dummy assignment operator
+   };
+
+   Float_t OverlapFactor(AliTPCseed * s1, AliTPCseed * s2, Int_t &sum1, Int_t &sum2);
+   void  SignShared(AliTPCseed * s1, AliTPCseed * s2);
+   void  RemoveOverlap(TObjArray * arr, Float_t factor, Int_t removalindex, Bool_t shared=kFALSE);
+private:
+   Float_t  GetSigmaY(AliTPCseed * seed);
+   Float_t  GetSigmaZ(AliTPCseed * seed);
+
+   void MakeSeeds(TObjArray * arr, Int_t sec, Int_t i1, Int_t i2);
+   TObjArray * MakeSeedsSectors(Int_t sec1, Int_t sec2);   // make seeds from all sectors
+   void MakeSeedsAll();
+   Int_t FollowProlongation(AliTPCseed& t, Int_t rf=0);
+   //Int_t FollowProlongationFast(AliTPCseed& t, Int_t step);   
+   Int_t FollowBackProlongation(AliTPCseed& t, Int_t rf=0);
+
+   Int_t FollowToNext(AliTPCseed& t, Int_t nr);
+   Int_t UpdateClusters(AliTPCseed& t, Int_t trindex,  Int_t nr);
+   Int_t FollowToNextCluster( Int_t trindex, Int_t nr);
+
+   virtual Int_t PropagateBack (const TFile *, TFile *){return 0;}
+   void ParallelTracking(Int_t rfirst, Int_t rlast);
+   void SetSampledEdx(AliTPCseed *t, Float_t q, Int_t i) {;}
+   Int_t UpdateTrack(AliTPCseed *t,AliTPCclusterMI* c, Double_t chi2, UInt_t i); //update trackinfo
+
+   //   Int_t FollowBackProlongation(AliTPCseed &s, const AliTPCtrack &t);
+
+   AliTPCtrackerMI(const AliTPCtrackerMI& r);           //dummy copy constructor
+   AliTPCtrackerMI &operator=(const AliTPCtrackerMI& r);//dummy assignment operator
+
+   const Int_t fkNIS;        //number of inner sectors
+   AliTPCSector *fInnerSec;  //array of inner sectors;
+   const Int_t fkNOS;        //number of outer sectors
+   AliTPCSector *fOuterSec;  //array of outer sectors;
+
+   Int_t fN;               //number of loaded sectors
+   AliTPCSector *fSectors; //pointer to loaded sectors;
+
+   Int_t fEventN;                      //event number
+   AliTPCClustersArray fClustersArray; //array of TPC clusters
+   Int_t fNtracks;                     //current number of tracks
+   TObjArray *fSeeds;                  //array of track seeds
+   //   TObjArray * fTrackPointPool;        // ! pool with track points
+   const AliTPCParam *fParam;          //pointer to the parameters
+};
+
+#endif
+
+
diff --git a/TPC/LandauTest.C b/TPC/LandauTest.C
new file mode 100644 (file)
index 0000000..4b9c638
--- /dev/null
@@ -0,0 +1,275 @@
+#include "TH1.h"
+#include "TH2.h"
+#include "TFile.h"
+#include "TTree.h"
+
+#include "TRandom.h"
+#include "TPad.h"
+#include "TCanvas.h"
+
+
+class TLandauMean: public TObject {
+public:
+  void Init(Int_t n, Float_t mean, Float_t sigma);  // initial parameters
+  void Gener();          // gener sample 
+  //  void Anal();
+  //
+  Int_t fNSample;      // number of samples
+  Float_t fLMean;        // landau mean
+  Float_t fLSigma;       // landau sigma
+  //
+  Float_t fTM_0_6[3];    // truncated method  - first 3 momenta
+  Float_t fTM_0_7[3];    // truncated method  - first 3 momenta
+  Float_t fTM_0_8[3];    // truncated method  - first 3 momenta
+  Float_t fTM_0_10[3];   // truncated method  - first 3 momenta
+  //
+  Float_t fLM_0_6[3];    // truncated log.  method  - first 3 momenta
+  Float_t fLM_0_7[3];    // truncated log.  method  - first 3 momenta
+  Float_t fLM_0_8[3];    // truncated log.  method  - first 3 momenta
+  Float_t fLM_0_10[3];   // truncated log.  method  - first 3 momenta
+
+  Float_t fMedian3;      // median 3 value
+private:
+  Float_t Moment3(Float_t sum1, Float_t sum2, Float_t sum3, Int_t n, Float_t m[3]);
+  ClassDef(TLandauMean,1)
+};
+
+ClassImp(TLandauMean)
+
+void TLandauMean::Init(Int_t n, Float_t mean, Float_t sigma)
+{
+  //
+  //init parameters
+  fNSample = n;
+  fLMean   = mean;
+  fLSigma  = sigma;
+}
+
+Float_t TLandauMean::Moment3(Float_t sumi1, Float_t sumi2, Float_t sumi3, Int_t sum, Float_t m[3])
+{
+  Float_t m3=0;
+
+  //  m3 = (sumi3-3*pos*sumi2+3*pos*pos*sumi-pos*pos*pos*sum)/sum; 
+  Float_t pos = sumi1/sum;
+  m[0] = pos;
+  m[1] = sumi2/sum-pos*pos;
+  if (m[1]<=0){
+    printf("pici pici\n");
+  }
+  else
+    m[1] = TMath::Sqrt(m[1]);
+  m3 = (sumi3-3*pos*sumi2+3*pos*pos*sumi1-pos*pos*pos*sum)/sum; 
+  Float_t sign = m3/TMath::Abs(m3);
+  m3 = TMath::Power(sign*m3,1/3.);
+  m3*=sign;
+
+  m[2] = m3;  
+  return m3;
+}
+
+void TLandauMean::Gener()
+{
+  // 
+  // generate sample
+  Float_t * buffer = new Float_t[fNSample];
+  
+  for (Int_t i=0;i<fNSample;i++) {
+    buffer[i] = gRandom->Landau(fLMean,fLSigma); 
+    if (buffer[i]>1000) buffer[i]=1000;
+  }
+
+  Int_t *index = new Int_t[fNSample];
+  TMath::Sort(fNSample,buffer,index,kFALSE);
+
+  //
+  Float_t median = buffer[index[fNSample/3]];
+  //
+  Float_t sum06[4]  = {0.,0.,0.,0.};
+  Float_t sum07[4]  = {0.,0.,0.,0.};
+  Float_t sum08[4]  = {0.,0.,0.,0.};
+  Float_t sum010[4]  = {0.,0.,0.,0.};
+  //
+  Float_t suml06[4] = {0.,0.,0.,0.};
+  Float_t suml07[4] = {0.,0.,0.,0.};
+  Float_t suml08[4] = {0.,0.,0.,0.};
+  Float_t suml010[4] = {0.,0.,0.,0.};
+  //
+
+  for (Int_t i =0; i<fNSample; i++){
+    Float_t amp  = buffer[index[i]];
+    Float_t lamp = median*TMath::Log(1.+amp/median);
+    if (i<0.6*fNSample){
+      sum06[0]+= amp;
+      sum06[1]+= amp*amp;
+      sum06[2]+= amp*amp*amp;
+      sum06[3]++;
+      suml06[0]+= lamp;
+      suml06[1]+= lamp*lamp;
+      suml06[2]+= lamp*lamp*lamp;
+      suml06[3]++;
+    }
+
+    if (i<0.7*fNSample){
+      sum07[0]+= amp;
+      sum07[1]+= amp*amp;
+      sum07[2]+= amp*amp*amp;
+      sum07[3]++;
+      suml07[0]+= lamp;
+      suml07[1]+= lamp*lamp;
+      suml07[2]+= lamp*lamp*lamp;
+      suml07[3]++;
+    }
+    if (i<0.8*fNSample){
+      sum08[0]+= amp;
+      sum08[1]+= amp*amp;
+      sum08[2]+= amp*amp*amp;
+      sum08[3]++;
+      suml08[0]+= lamp;
+      suml08[1]+= lamp*lamp;
+      suml08[2]+= lamp*lamp*lamp;
+      suml08[3]++;
+    }
+    if (i<1*fNSample){
+      sum010[0]+= amp;
+      sum010[1]+= amp*amp;
+      sum010[2]+= amp*amp*amp;
+      sum010[3]++;
+      suml010[0]+= lamp;
+      suml010[1]+= lamp*lamp;
+      suml010[2]+= lamp*lamp*lamp;
+      suml010[3]++;
+
+    }
+  }
+  //  
+  fMedian3 = median;
+  //
+  Moment3(sum06[0],sum06[1],sum06[2],sum06[3],fTM_0_6);  
+  Moment3(sum07[0],sum07[1],sum07[2],sum07[3],fTM_0_7);  
+  Moment3(sum08[0],sum08[1],sum08[2],sum08[3],fTM_0_8);  
+  Moment3(sum010[0],sum010[1],sum010[2],sum010[3],fTM_0_10);  
+  //
+
+  Moment3(suml06[0],suml06[1],suml06[2],suml06[3],fLM_0_6);  
+  Moment3(suml07[0],suml07[1],suml07[2],suml07[3],fLM_0_7);  
+  Moment3(suml08[0],suml08[1],suml08[2],suml08[3],fLM_0_8);  
+  Moment3(suml010[0],suml010[1],suml010[2],suml010[3],fLM_0_10);  
+  //
+  fLM_0_6[0] = (TMath::Exp(fLM_0_6[0]/median)-1.)*median;
+  fLM_0_7[0] = (TMath::Exp(fLM_0_7[0]/median)-1.)*median;
+  fLM_0_8[0] = (TMath::Exp(fLM_0_8[0]/median)-1.)*median;
+  fLM_0_10[0] = (TMath::Exp(fLM_0_10[0]/median)-1.)*median;
+  //
+  delete [] buffer;
+}   
+
+
+void GenerLandau(Int_t nsamples)
+{
+  TLandauMean * landau = new TLandauMean;
+  TFile f("Landau.root","recreate");
+  TTree * tree = new TTree("Landau","Landau");
+  tree->Branch("Landau","TLandauMean",&landau); 
+  
+  for (Int_t i=0;i<nsamples;i++){
+    Int_t   n     = 20 + Int_t(gRandom->Rndm()*150);
+    Float_t mean  = 40. +gRandom->Rndm()*50.;
+    Float_t sigma = 5.  +gRandom->Rndm()*15.;
+    landau->Init(n, mean, sigma);
+    landau->Gener();
+    tree->Fill();
+  }
+  tree->Write();
+  f.Close();
+
+}
+
+
+
+
+
+TH1F *  LandauTest(Float_t meano,  Float_t sigma, Float_t meanlog0, Int_t n,Float_t ratio)
+{ 
+  //
+  // test for different approach of de dx resolution
+  // meano, sigma  - mean value of Landau distribution and sigma
+  // meanlog0      - scaling factor for logarithmic mean value
+  // n             - number of used layers
+  // ratio         - ratio of used amplitudes for truncated mean
+  //
+  
+
+  TCanvas * pad = new TCanvas("Landau test");
+  pad->Divide(2,2);
+  TH1F * h1 = new TH1F("h1","Logarithmic mean",300,0,4*meano);
+  TH1F * h2 = new TH1F("h2","Logarithmic amplitudes",300,0,8*meano);
+  TH1F * h3 = new TH1F("h3","Mean",300,0,4*meano);
+  TH1F * h4 = new TH1F("h4","Amplitudes",300,0,8*meano);
+
+  for(Int_t j=0;j<10000;j++){
+    //generate sample and sort it
+    Float_t * buffer = new Float_t[n];
+    Float_t * buffer2= new Float_t[n];
+    
+    for (Int_t i=0;i<n;i++) {
+      buffer[i] = gRandom->Landau(meano,sigma); 
+      buffer2[i] = buffer[i];    
+    }
+    //add crosstalk
+    for (Int_t i=1;i<n-1;i++) {
+      buffer[i] =    buffer2[i]*1.0+ buffer2[i-1]*0.0+ buffer2[i+1]*0.0;
+      buffer[i] = TMath::Min(buffer[i],1000.); 
+    }
+    Int_t *index = new Int_t[n];
+    TMath::Sort(n,buffer,index,kFALSE);
+
+    //calculate mean
+    Float_t sum;
+    sum=0;
+    Float_t mean;
+    Float_t used = 0;
+    for (Int_t i=0;i<n*ratio;i++) {
+      if (buffer[index[i]]<1000.){
+       Float_t amp = meanlog0*TMath::Log(1+buffer[index[i]]/meanlog0);
+       sum += amp;            
+       used++;
+      }
+    }
+    mean = sum/used;
+    //
+    sum=0;
+    used=0;
+    Float_t sum2=0;
+    Float_t meanlog =meanlog0;
+    for (Int_t i=0;i<n*ratio;i++) {
+      if (buffer[index[i]]<1000.){
+       Float_t amp = meanlog*TMath::Log(1.+buffer[index[i]]/(meanlog));
+       sum +=amp;
+       sum2+=buffer[index[i]];
+       used++;
+       h2->Fill(amp);
+       h4->Fill(buffer[index[i]]);
+      }
+    }
+    mean = sum/used;
+    mean = (TMath::Exp(mean/meanlog)-1)*meanlog;
+    Float_t mean2 = sum2/used;
+    //mean2 = (mean+mean2)/2.;
+    h1->Fill(mean);    
+    h3->Fill(mean2);
+  }
+
+  pad->cd(1);
+  h1->Draw();
+  pad->cd(2);
+  h2->Draw();
+  pad->cd(3);
+  h3->Draw();
+  pad->cd(4);
+  h4->Draw();
+
+
+  return h1;
+
+}
+
index b137239..9317179 100644 (file)
@@ -20,6 +20,7 @@
 #pragma link C++ class  AliTPCtrack+;
 #pragma link C++ class  AliTPCtracker+;
 
+
 #pragma link C++ class  AliTPCParam+;
 #pragma link C++ class  AliTPCParamSR-;
 #pragma link C++ class  AliTPCParamCR-;
@@ -31,7 +32,6 @@
 
 #pragma link C++ class AliComplexCluster+;
 #pragma link C++ class AliDigitCluster+;
-#pragma link C++ class AliDifCluster+;
 #pragma link C++ class AliClusters+;
 #pragma link C++ class AliClustersArray+;
 
@@ -44,7 +44,6 @@
 #pragma link C++ class AliTPCDigitsArray+;
 
 
-#pragma link C++ class AliTPCClusterFinder+;
 #pragma link C++ class AliH2F+;
 
 #pragma link C++ class AliTPCTrackHits+;
 #pragma link C++ class AliTPCtrackerParam; 
 #pragma link C++ class AliTPCkineGrid;
 
+// points used in new cluster finder
+#pragma link C++ class  AliTPCExactPoint+;
+#pragma link C++ class  AliTPCTrackerPoint+;
+#pragma link C++ class  AliTPCClusterPoint+;
+#pragma link C++ class  AliTPCTrackPoint+;
+#pragma link C++ class  AliTPCTrackPointRef+;
+
+#pragma link C++ class  AliTPCclustererMI+;
+#pragma link C++ class  AliTPCclusterMI+;
+#pragma link C++ class  AliTPCclusterLMI+;
+#pragma link C++ class  AliTPCtrackerMI+;
+#pragma link C++ class  AliTPCseed+;
+
+
+
 
 #endif
 
index e96a42f..30dd617 100644 (file)
@@ -39,9 +39,8 @@
 #include "AliClusters.h"
 #include "AliTPCClustersRow.h"
 #include "AliTPCClustersArray.h"
-#include "AliTPCClusterFinder.h"
 #include "AliTPCcluster.h"
-
+#include "AliH2F.h"
 #include "TMinuit.h"
 #include "AliTPC.h"
 #include "AliTPCv1.h"
index 639b0ff..45f6a46 100644 (file)
@@ -1,4 +1,4 @@
-SRCS:=  AliTPCClusterFinder.cxx  AliTPC.cxx  AliTPCv0.cxx    AliTPCv1.cxx      AliTPCv2.cxx  \
+SRCS:=  AliTPC.cxx  AliTPCv0.cxx    AliTPCv1.cxx      AliTPCv2.cxx  \
        AliTPCv3.cxx AliDetectorParam.cxx AliTPCParam.cxx \
        AliTPCParamSR.cxx AliTPCParamCR.cxx  AliTPCPRF2D.cxx \
        AliTPCRF1D.cxx \
@@ -10,8 +10,8 @@ SRCS:=  AliTPCClusterFinder.cxx  AliTPC.cxx  AliTPCv0.cxx    AliTPCv1.cxx      A
         AliTPCtrack.cxx AliTPCtracker.cxx \
        AliTPCTrackHits.cxx\
        AliTPCDigitizer.cxx\
-       AliTPCTrackHitsV2.cxx AliTPCtrackerParam.cxx AliTPCkineGrid.cxx
-
+       AliTPCTrackHitsV2.cxx AliTPCtrackerParam.cxx AliTPCkineGrid.cxx \
+       AliTPCclustererMI.cxx AliTPCclusterMI.cxx  AliTPCtrackerMI.cxx
 HDRS:= $(SRCS:.cxx=.h)