New version of the V0 finder (M.Ivanov)
authorhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 29 Apr 2005 03:23:08 +0000 (03:23 +0000)
committerhristov <hristov@f7af4fe6-9843-0410-8265-dc069ae4e863>
Fri, 29 Apr 2005 03:23:08 +0000 (03:23 +0000)
23 files changed:
ITS/AliITStrackV2.cxx
ITS/AliITStrackerMI.cxx
ITS/AliITStrackerMI.h
ITS/AliITStrackerV2.cxx
STEER/AliESDComparisonMI.C
STEER/AliESDV0MI.cxx
STEER/AliESDV0MI.h
STEER/AliESDtrack.cxx
STEER/AliESDtrack.h
STEER/AliESDv0.h
STEER/AliGenInfo.C
STEER/AliGenInfo.h
STEER/AliHelix.cxx
STEER/AliHelix.h
STEER/AliMagFMaps.cxx
STEER/TTreeStream.cxx
TPC/AliTPCseed.cxx [new file with mode: 0644]
TPC/AliTPCseed.h [new file with mode: 0644]
TPC/AliTPCtrack.cxx
TPC/AliTPCtrack.h
TPC/AliTPCtrackerMI.cxx
TPC/AliTPCtrackerMI.h
TPC/libTPCrec.pkg

index 43a394e43ca25123abd0a7bd7d58252e5350a325..af3a54b691aba9a8bea93f5f9b744194ead02245 100644 (file)
@@ -107,9 +107,9 @@ AliKalmanTrack() {
   }
   fESDtrack=&t;
 
+  //  if (!Invariant()) throw "AliITStrackV2: conversion failed !\n";
   for(Int_t i=0; i<4; i++) fdEdxSample[i]=0;
 
-  if (!Invariant()) throw "AliITStrackV2: conversion failed !\n";
 }
 
 void AliITStrackV2::UpdateESDtrack(ULong_t flags) const {
@@ -278,10 +278,12 @@ Int_t AliITStrackV2::PropagateTo(Double_t xk, Double_t d, Double_t x0) {
   //------------------------------------------------------------------
   Double_t x1=fX, x2=xk, dx=x2-x1;
   Double_t f1=fP2, f2=f1 + fP4*dx;
-  if (TMath::Abs(f2) >= 0.9999) {
-    Int_t n=GetNumberOfClusters();
-    if (n>kWARN) 
-       Warning("PropagateTo","Propagation failed !\n",n);
+  if (TMath::Abs(f2) >= 0.98) {   
+    // MI change  - don't propagate highly inclined tracks
+    //              covariance matrix distorted
+    // Int_t n=GetNumberOfClusters();
+    //     if (n>kWARN) 
+    //        Warning("PropagateTo","Propagation failed !\n",n);
     return 0;
   }
 
@@ -504,12 +506,16 @@ Int_t AliITStrackV2::Propagate(Double_t alp,Double_t xk) {
   {
   Double_t dx=xk-fX;
   Double_t f1=fP2, f2=f1 + fP4*dx;
-  if (TMath::Abs(f2) >= 0.9999) {
-    Int_t n=GetNumberOfClusters();
-    if (n>kWARN) 
-       Warning("Propagate","Propagation failed (%d) !\n",n);
+  if (TMath::Abs(f2) >= 0.98) {  
+    // don't propagate highly inclined tracks MI
     return 0;
   }
+ //  if (TMath::Abs(f2) >= 0.9999) {
+//     Int_t n=GetNumberOfClusters();
+//     if (n>kWARN) 
+//        Warning("Propagate","Propagation failed (%d) !\n",n);
+//     return 0;
+//   }
   Double_t r1=TMath::Sqrt(1.- f1*f1), r2=TMath::Sqrt(1.- f2*f2);
   
   fX=xk;
index 998b8c6f25ae9f6a4db91ea86e45c242d10c950d..647ef97b81248989d0141e4b1923070c5d024857 100644 (file)
@@ -100,6 +100,9 @@ AliITStrackerMI::AliITStrackerMI(const AliITSgeom *geom) : AliTracker() {
   for (Int_t i=0;i<100000;i++){
     fBestTrackIndex[i]=0;
   }
+  //
+  fDebugStreamer = new TTreeSRedirector("ITSdebug.root");
+
 }
 
 AliITStrackerMI::~AliITStrackerMI()
@@ -108,6 +111,10 @@ AliITStrackerMI::~AliITStrackerMI()
   //destructor
   //
   if (fCoeficients) delete []fCoeficients;
+  if (fDebugStreamer) {
+    //fDebugStreamer->Close();
+    delete fDebugStreamer;
+  }
 }
 
 void AliITStrackerMI::SetLayersNotToSkip(Int_t *l) {
@@ -231,6 +238,7 @@ Int_t AliITStrackerMI::Clusters2Tracks(AliESD *event) {
   //--------------------------------------------------------------------
   TObjArray itsTracks(15000);
   fOriginal.Clear();
+  fEsd = event;         // store pointer to the esd 
   {/* Read ESD tracks */
     Int_t nentr=event->GetNumberOfTracks();
     Info("Clusters2Tracks", "Number of ESD tracks: %d\n", nentr);
@@ -256,25 +264,30 @@ Int_t AliITStrackerMI::Clusters2Tracks(AliESD *event) {
       // write expected q
       t->fExpQ = TMath::Max(0.8*t->fESDtrack->GetTPCsignal(),30.);
 
-      if (TMath::Abs(t->fD[0])>10) {
-       delete t;
-       continue;
-      }
-
-      if (TMath::Abs(vdist)>20) {
-       delete t;
-       continue;
-      }
-      if (TMath::Abs(1/t->Get1Pt())<0.120) {
-       delete t;
-       continue;
+      if (esd->GetV0Index(0)>0 && t->fD[0]<30){
+       //track - can be  V0 according to TPC
       }
-
-      if (CorrectForDeadZoneMaterial(t)!=0) {
-       //Warning("Clusters2Tracks",
-       //        "failed to correct for the material in the dead zone !\n");
-         delete t;
-         continue;
+      else{    
+       if (TMath::Abs(t->fD[0])>10) {
+         delete t;
+         continue;
+       }
+       
+       if (TMath::Abs(vdist)>20) {
+         delete t;
+         continue;
+       }
+       if (TMath::Abs(1/t->Get1Pt())<0.120) {
+         delete t;
+         continue;
+       }
+       
+       if (CorrectForDeadZoneMaterial(t)!=0) {
+         //Warning("Clusters2Tracks",
+         //        "failed to correct for the material in the dead zone !\n");
+         delete t;
+         continue;
+       }
       }
       t->fReconstructed = kFALSE;
       itsTracks.AddLast(t);
@@ -286,6 +299,7 @@ Int_t AliITStrackerMI::Clusters2Tracks(AliESD *event) {
   fOriginal.Sort();
   Int_t nentr=itsTracks.GetEntriesFast();
   fTrackHypothesys.Expand(nentr);
+  fBestHypothesys.Expand(nentr);
   MakeCoeficients(nentr);
   Int_t ntrk=0;
   for (fPass=0; fPass<2; fPass++) {
@@ -303,9 +317,7 @@ Int_t AliITStrackerMI::Clusters2Tracks(AliESD *event) {
        fI = 6;
        ResetTrackToFollow(*t);
        ResetBestTrack();
-      
-       FollowProlongationTree(t,i);
-
+       FollowProlongationTree(t,i,fConstraint[fPass]);
 
        SortTrackHypothesys(fCurrentEsdTrack,20,0);  //MI change
        //
@@ -336,7 +348,8 @@ Int_t AliITStrackerMI::Clusters2Tracks(AliESD *event) {
   }
 
   //GetBestHypothesysMIP(itsTracks);
-  FindV0(event);
+  UpdateTPCV0(event);
+  FindV02(event);
   fAfterV0 = kTRUE;
   //GetBestHypothesysMIP(itsTracks);
   //
@@ -350,6 +363,7 @@ Int_t AliITStrackerMI::Clusters2Tracks(AliESD *event) {
   }
 
   fTrackHypothesys.Delete();
+  fBestHypothesys.Delete();
   fOriginal.Clear();
   delete []fCoeficients;
   fCoeficients=0;
@@ -483,9 +497,6 @@ Int_t AliITStrackerMI::RefitInward(AliESD *event) {
 
         if (fTrackToFollow.Propagate(fv+a,xv)) {
             fTrackToFollow.UpdateESDtrack(AliESDtrack::kITSrefit);
-           Float_t d=fTrackToFollow.GetD(GetX(),GetY());        
-           Float_t z=fTrackToFollow.GetZ()-GetZ();      
-           fTrackToFollow.GetESDtrack()->SetImpactParameters(d,z);      
             //UseClusters(&fTrackToFollow);
             {
             AliITSclusterV2 c; c.SetY(yv); c.SetZ(GetZ());
@@ -520,12 +531,40 @@ AliCluster *AliITStrackerMI::GetCluster(Int_t index) const {
 }
 
 
-void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex) 
+void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain
 {
   //--------------------------------------------------------------------
   // Follow prolongation tree
   //--------------------------------------------------------------------
+  //
+  AliESDtrack * esd = otrack->fESDtrack;
+  if (esd->GetV0Index(0)>0){
+    //
+    // TEMPORARY SOLLUTION: map V0 indexes to point to proper track
+    //                      mapping of esd track is different as its track in Containers
+    //                      Need something more stable
+    //                      Indexes are set back againg to the ESD track indexes in UpdateTPCV0
+    for (Int_t i=0;i<3;i++){
+      Int_t  index = esd->GetV0Index(i);
+      if (index==0) break;
+      AliESDV0MI * vertex = fEsd->GetV0MI(index);
+      if (vertex->GetStatus()<0) continue;     // rejected V0
+      //
+      if (esd->GetSign()>0) {
+        vertex->SetIndex(0,esdindex);
+      }
+      else{
+        vertex->SetIndex(1,esdindex);
+      }
+    }
+  }
+  TObjArray *bestarray = (TObjArray*)fBestHypothesys.At(esdindex);
+  if (!bestarray){
+    bestarray = new TObjArray(5);
+    fBestHypothesys.AddAt(bestarray,esdindex);
+  }
 
+  //
   //setup tree of the prolongations
   //
   static AliITStrackMI tracks[7][100];
@@ -540,7 +579,7 @@ void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdin
   otrack->fNSkipped=0;
   new (&(tracks[6][0])) AliITStrackMI(*otrack);
   ntracks[6]=1;
-  nindexes[6][0]=0;
+  for (Int_t i=0;i<7;i++) nindexes[i][0]=0;
   // 
   //
   // follow prolongations
@@ -626,7 +665,7 @@ void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdin
       //
       Double_t msz=1./((currenttrack1.GetSigmaZ2() + 16.*kSigmaZ2[ilayer]));
       Double_t msy=1./((currenttrack1.GetSigmaY2() + 16.*kSigmaY2[ilayer]));
-      if (fConstraint[fPass]){
+      if (constrain){
        msy/=60; msz/=60.;
       }
       else{
@@ -687,8 +726,8 @@ void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdin
          Double_t x0;
          Double_t d=layer.GetThickness(updatetrack->GetY(),updatetrack->GetZ(),x0);
          updatetrack->CorrectForMaterial(d,x0);          
-         if (fConstraint[fPass]) {
-           updatetrack->fConstrain = fConstraint[fPass];
+         if (constrain) {
+           updatetrack->fConstrain = constrain;
            fI = ilayer;
            Double_t d=GetEffectiveThickness(0,0); //Think of this !!!!
            Double_t xyz[]={GetX(),GetY(),GetZ()};
@@ -705,8 +744,8 @@ void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdin
          ntracks[ilayer]++;
        }  // create new hypothesy 
       } // loop over possible cluster prolongation      
-      //      if (fConstraint[fPass]&&itrack<2&&currenttrack1.fNSkipped==0 && deadzone==0){    
-      if (fConstraint[fPass]&&itrack<2&&currenttrack1.fNSkipped==0 && deadzone==0&&ntracks[ilayer]<100){       
+      //      if (constrain&&itrack<2&&currenttrack1.fNSkipped==0 && deadzone==0){     
+      if (constrain&&itrack<2&&currenttrack1.fNSkipped==0 && deadzone==0&&ntracks[ilayer]<100){        
        AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
        vtrack->fClIndex[ilayer]=0;
        fI = ilayer;
@@ -718,7 +757,7 @@ void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdin
        ntracks[ilayer]++;
       }
 
-      if (fConstraint[fPass]&&itrack<1&&TMath::Abs(currenttrack1.fP3)>1.1){  //big theta -- for low mult. runs
+      if (constrain&&itrack<1&&TMath::Abs(currenttrack1.fP3)>1.1){  //big theta -- for low mult. runs
        AliITStrackMI* vtrack = new (&tracks[ilayer][ntracks[ilayer]]) AliITStrackMI(currenttrack1);
        vtrack->fClIndex[ilayer]=0;
        fI = ilayer;
@@ -742,8 +781,8 @@ void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdin
       if ( normalizedchi2[itrack]<3+0.5*ilayer) golds++;
       if (ilayer>4) accepted++;
       else{
-       if ( fConstraint[fPass] && normalizedchi2[itrack]<kMaxNormChi2C[ilayer]+1) accepted++;
-       if (!fConstraint[fPass] && normalizedchi2[itrack]<kMaxNormChi2NonC[ilayer]+1) accepted++;
+       if ( constrain && normalizedchi2[itrack]<kMaxNormChi2C[ilayer]+1) accepted++;
+       if (!constrain && normalizedchi2[itrack]<kMaxNormChi2NonC[ilayer]+1) accepted++;
       }
     }
     TMath::Sort(ntracks[ilayer],normalizedchi2,nindexes[ilayer],kFALSE);
@@ -752,20 +791,20 @@ void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdin
     if (ntracks[ilayer]>90) ntracks[ilayer]=90; 
   } //loop over layers
   //printf("%d\t%d\t%d\t%d\t%d\t%d\n",ntracks[0],ntracks[1],ntracks[2],ntracks[3],ntracks[4],ntracks[5]);
-  Int_t max = fConstraint[fPass]? 20: 5;
+  Int_t max = constrain? 20: 5;
 
   for (Int_t i=0;i<TMath::Min(max,ntracks[0]);i++) {
     AliITStrackMI & track= tracks[0][nindexes[0][i]];
     if (track.GetNumberOfClusters()<2) continue;
-    if (!fConstraint[fPass]&&track.fNormChi2[0]>7.)continue;
+    if (!constrain&&track.fNormChi2[0]>7.)continue;
     AddTrackHypothesys(new AliITStrackMI(track), esdindex);
   }
   for (Int_t i=0;i<TMath::Min(2,ntracks[1]);i++) {
     AliITStrackMI & track= tracks[1][nindexes[1][i]];
     if (track.GetNumberOfClusters()<4) continue;
-    if (!fConstraint[fPass]&&track.fNormChi2[1]>7.)continue;
-    if (fConstraint[fPass]) track.fNSkipped+=1;
-    if (!fConstraint[fPass]) {
+    if (!constrain&&track.fNormChi2[1]>7.)continue;
+    if (constrain) track.fNSkipped+=1;
+    if (!constrain) {
       track.fD[0] = track.GetD(GetX(),GetY());   
       track.fNSkipped+=4./(4.+8.*TMath::Abs(track.fD[0]));
       if (track.fN+track.fNDeadZone+track.fNSkipped>6) {
@@ -776,13 +815,13 @@ void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdin
   }
   //}
   
-  if (!fConstraint[fPass]){  
+  if (!constrain){  
     for (Int_t i=0;i<TMath::Min(2,ntracks[2]);i++) {
       AliITStrackMI & track= tracks[2][nindexes[2][i]];
       if (track.GetNumberOfClusters()<3) continue;
-      if (!fConstraint[fPass]&&track.fNormChi2[2]>7.)continue;
-      if (fConstraint[fPass]) track.fNSkipped+=2;      
-      if (!fConstraint[fPass]){
+      if (!constrain&&track.fNormChi2[2]>7.)continue;
+      if (constrain) track.fNSkipped+=2;      
+      if (!constrain){
        track.fD[0] = track.GetD(GetX(),GetY());
        track.fNSkipped+= 7./(7.+8.*TMath::Abs(track.fD[0]));
        if (track.fN+track.fNDeadZone+track.fNSkipped>6) {
@@ -792,6 +831,76 @@ void AliITStrackerMI::FollowProlongationTree(AliITStrackMI * otrack, Int_t esdin
       AddTrackHypothesys(new AliITStrackMI(track), esdindex);
     }
   }
+  
+  if (!constrain){
+    //
+    // register best tracks - important for V0 finder
+    //
+    for (Int_t ilayer=0;ilayer<5;ilayer++){
+      if (ntracks[ilayer]==0) continue;
+      AliITStrackMI & track= tracks[ilayer][nindexes[ilayer][0]];
+      if (track.GetNumberOfClusters()<1) continue;
+      CookLabel(&track,0);
+      bestarray->AddAt(new AliITStrackMI(track),ilayer);
+    }
+  }
+  //
+  // update TPC V0 information
+  //
+  if (otrack->fESDtrack->GetV0Index(0)>0){    
+    Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
+    for (Int_t i=0;i<3;i++){
+      Int_t  index = otrack->fESDtrack->GetV0Index(i); 
+      if (index==0) break;
+      AliESDV0MI * vertex = fEsd->GetV0MI(index);
+      if (vertex->GetStatus()<0) continue;     // rejected V0
+      //
+      if (otrack->fP4>0) {
+       vertex->SetIndex(0,esdindex);
+      }
+      else{
+       vertex->SetIndex(1,esdindex);
+      }
+      //find nearest layer with track info
+      Int_t nearestold  = GetNearestLayer(vertex->GetXrp());
+      Int_t nearest     = nearestold; 
+      for (Int_t ilayer =nearest;ilayer<8;ilayer++){
+       if (ntracks[nearest]==0){
+         nearest = ilayer;
+       }
+      }
+      //
+      AliITStrackMI & track= tracks[nearest][nindexes[nearest][0]];
+      if (nearestold<5&&nearest<5){
+       Bool_t accept = track.fNormChi2[nearest]<10; 
+       if (accept){
+         if (track.fP4>0) {
+           vertex->SetP(track);
+           vertex->Update(fprimvertex);
+           //      vertex->SetIndex(0,track.fESDtrack->GetID()); 
+           if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
+         }else{
+           vertex->SetM(track);
+           vertex->Update(fprimvertex);
+           //vertex->SetIndex(1,track.fESDtrack->GetID());
+           if (track.GetNumberOfClusters()>2) AddTrackHypothesys(new AliITStrackMI(track), esdindex);
+         }
+         vertex->SetStatus(vertex->GetStatus()+1);
+       }else{
+         //  vertex->SetStatus(-2);  // reject V0  - not enough clusters
+       }
+      }
+      // if (nearestold>3){
+//     Int_t indexlayer = (ntracks[0]>0)? 0:1;
+//     if (ntracks[indexlayer]>0){
+//       AliITStrackMI & track= tracks[indexlayer][nindexes[indexlayer][0]];
+//       if (track.GetNumberOfClusters()>4&&track.fNormChi2[indexlayer]<4){
+//         vertex->SetStatus(-1);  // reject V0 - clusters before
+//       }
+//     }
+//      }
+    }
+  }  
 }
 
 
@@ -2389,7 +2498,7 @@ AliITStrackMI * AliITStrackerMI::GetBestHypothesys(Int_t esdindex, AliITStrackMI
       if (!backtrack->Improve(0,xyzv,ersv))         continue;                    
       if (!backtrack->PropagateToVertex())          continue;
       backtrack->ResetCovariance();      
-      if (!backtrack->Improve(0,xyzv,ersv))         continue;                    
+      if (!backtrack->Improve(0,xyzv,ersv))         continue;                          
     }else{
       backtrack->ResetCovariance();
     }
@@ -2975,7 +3084,6 @@ void   AliITStrackerMI::GetDCASigma(AliITStrackMI* track, Float_t & sigmarfi, Fl
   // 
   sigmarfi = 0.004+1.4 *TMath::Abs(track->fP4)+332.*track->fP4*track->fP4;
   sigmaz   = 0.011+4.37*TMath::Abs(track->fP4);
-
 }
 
 
@@ -3064,607 +3172,685 @@ void AliITStrackerMI::UpdateESDtrack(AliITStrackMI* track, ULong_t flags) const
 
 
 
+Int_t AliITStrackerMI::GetNearestLayer(const Double_t *xr) const{
+  //
+  //Get nearest upper layer close to the point xr.
+  // rough approximation 
+  //
+  const Float_t radiuses[6]={4,6.5,15.03,24.,38.5,43.7};
+  Float_t radius = TMath::Sqrt(xr[0]*xr[0]+xr[1]*xr[1]);
+  Int_t res =6;
+  for (Int_t i=0;i<6;i++){
+    if (radius<radiuses[i]){
+      res =i;
+      break;
+    }
+  }
+  return res;
+}
 
 
+void AliITStrackerMI::UpdateTPCV0(AliESD *event){
+  //
+  //try to update, or reject TPC  V0s
+  //
+  Int_t nv0s = event->GetNumberOfV0MIs();
+  Int_t nitstracks = fTrackHypothesys.GetEntriesFast();
 
+  for (Int_t i=0;i<nv0s;i++){
+    AliESDV0MI * vertex = event->GetV0MI(i);
+    Int_t ip = vertex->GetIndex(0);
+    Int_t im = vertex->GetIndex(1);
+    //
+    TObjArray * arrayp = (ip<nitstracks) ? (TObjArray*)fTrackHypothesys.At(ip):0;
+    TObjArray * arraym = (im<nitstracks) ? (TObjArray*)fTrackHypothesys.At(im):0;
+    AliITStrackMI * trackp = (arrayp!=0) ? (AliITStrackMI*)arrayp->At(0):0;
+    AliITStrackMI * trackm = (arraym!=0) ? (AliITStrackMI*)arraym->At(0):0;
+    //
+    //
+    if (trackp){
+      if (trackp->fN+trackp->fNDeadZone>5.5){
+       if (trackp->fConstrain&&trackp->fChi2MIP[0]<3) vertex->SetStatus(-100);
+       if (!trackp->fConstrain&&trackp->fChi2MIP[0]<2) vertex->SetStatus(-100); 
+      }
+    }
 
+    if (trackm){
+      if (trackm->fN+trackm->fNDeadZone>5.5){
+       if (trackm->fConstrain&&trackm->fChi2MIP[0]<3) vertex->SetStatus(-100);
+       if (!trackm->fConstrain&&trackm->fChi2MIP[0]<2) vertex->SetStatus(-100); 
+      }
+    }
+    if (vertex->GetStatus()==-100) continue;
+    //
+    Int_t clayer = GetNearestLayer(vertex->GetXrp());
+    vertex->SetNBefore(clayer);        //
+    vertex->SetChi2Before(9*clayer);   //
+    vertex->SetNAfter(6-clayer);       //
+    vertex->SetChi2After(0);           //
+    //
+    if (clayer >1 ){ // calculate chi2 before vertex
+      Float_t chi2p = 0, chi2m=0;  
+      //
+      if (trackp){
+       for (Int_t ilayer=0;ilayer<clayer;ilayer++){
+         if (trackp->fClIndex[ilayer]>0){
+           chi2p+=trackp->fDy[ilayer]*trackp->fDy[ilayer]/(trackp->fSigmaY[ilayer]*trackp->fSigmaY[ilayer])+
+             trackp->fDz[ilayer]*trackp->fDz[ilayer]/(trackp->fSigmaZ[ilayer]*trackp->fSigmaZ[ilayer]);
+         }
+         else{
+           chi2p+=9;
+         }
+       }
+      }else{
+       chi2p = 9*clayer;
+      }
+      //
+      if (trackm){
+       for (Int_t ilayer=0;ilayer<clayer;ilayer++){
+         if (trackm->fClIndex[ilayer]>0){
+           chi2m+=trackm->fDy[ilayer]*trackm->fDy[ilayer]/(trackm->fSigmaY[ilayer]*trackm->fSigmaY[ilayer])+
+             trackm->fDz[ilayer]*trackm->fDz[ilayer]/(trackm->fSigmaZ[ilayer]*trackm->fSigmaZ[ilayer]);
+         }
+         else{
+           chi2m+=9;
+         }
+       }
+      }else{
+       chi2m = 9*clayer;
+      }
+      vertex->SetChi2Before(TMath::Min(chi2p,chi2m));
+      if (TMath::Min(chi2p,chi2m)/Float_t(clayer)<4) vertex->SetStatus(-10);  // track exist before vertex
+    }
+    
+    if (clayer < 5 ){ // calculate chi2 after vertex
+      Float_t chi2p = 0, chi2m=0;  
+      //
+      if (trackp&&TMath::Abs(trackp->fP3)<1.){
+       for (Int_t ilayer=clayer;ilayer<6;ilayer++){
+         if (trackp->fClIndex[ilayer]>0){
+           chi2p+=trackp->fDy[ilayer]*trackp->fDy[ilayer]/(trackp->fSigmaY[ilayer]*trackp->fSigmaY[ilayer])+
+             trackp->fDz[ilayer]*trackp->fDz[ilayer]/(trackp->fSigmaZ[ilayer]*trackp->fSigmaZ[ilayer]);
+         }
+         else{
+           chi2p+=9;
+         }
+       }
+      }else{
+       chi2p = 0;
+      }
+      //
+      if (trackm&&TMath::Abs(trackm->fP3)<1.){
+       for (Int_t ilayer=clayer;ilayer<6;ilayer++){
+         if (trackm->fClIndex[ilayer]>0){
+           chi2m+=trackm->fDy[ilayer]*trackm->fDy[ilayer]/(trackm->fSigmaY[ilayer]*trackm->fSigmaY[ilayer])+
+             trackm->fDz[ilayer]*trackm->fDz[ilayer]/(trackm->fSigmaZ[ilayer]*trackm->fSigmaZ[ilayer]);
+         }
+         else{
+           chi2m+=9;
+         }
+       }
+      }else{
+       chi2m = 0;
+      }
+      vertex->SetChi2After(TMath::Max(chi2p,chi2m));
+      if (TMath::Max(chi2m,chi2p)/Float_t(6-clayer)>9) vertex->SetStatus(-20);  // track not found in ITS
+    }
+  }
+  //
+}
 
-void  AliITStrackerMI::FindV0(AliESD *event)
+void  AliITStrackerMI::FindV02(AliESD *event)
 {
   //
-  // fast V0 finder
+  // V0 finder
+  //
+  //  Cuts on DCA -  R dependent
+  //          max distance DCA between 2 tracks cut 
+  //          maxDist = TMath::Min(kMaxDist,kMaxDist0+pvertex->GetRr()*kMaxDist);
+  //
+  const Float_t kMaxDist0      = 0.1;    
+  const Float_t kMaxDist1      = 0.1;     
+  const Float_t kMaxDist       = 1;
+  const Float_t kMinPointAngle  = 0.85;
+  const Float_t kMinPointAngle2 = 0.99;
+  const Float_t kMinR           = 0.5;
+  const Float_t kMaxR           = 220;
+  //const Float_t kCausality0Cut   = 0.19;
+  //const Float_t kLikelihood01Cut = 0.25;
+  //const Float_t kPointAngleCut   = 0.9996;
+  const Float_t kCausality0Cut   = 0.19;
+  const Float_t kLikelihood01Cut = 0.45;
+  const Float_t kLikelihood1Cut  = 0.5;
+  const Float_t kCombinedCut     = 0.55;
+
+  //
+  //
+  TTreeSRedirector &cstream = *fDebugStreamer;
+  Int_t ntracks    = event->GetNumberOfTracks(); 
+  Int_t nitstracks = fTrackHypothesys.GetEntriesFast();
+  fOriginal.Expand(ntracks);
+  fTrackHypothesys.Expand(ntracks);
+  fBestHypothesys.Expand(ntracks);
+  //
+  AliHelix * helixes   = new AliHelix[ntracks+2];
+  TObjArray trackarray(ntracks+2);     //array with tracks - with vertex constrain
+  TObjArray trackarrayc(ntracks+2);    //array of "best    tracks" - without vertex constrain
+  TObjArray trackarrayl(ntracks+2);    //array of "longest tracks" - without vertex constrain
+  Bool_t * forbidden   = new Bool_t [ntracks+2];
+  Int_t   *itsmap      = new Int_t  [ntracks+2];
+  Float_t *dist        = new Float_t[ntracks+2];
+  Float_t *normdist0   = new Float_t[ntracks+2];
+  Float_t *normdist1   = new Float_t[ntracks+2];
+  Float_t *normdist    = new Float_t[ntracks+2];
+  Float_t *norm        = new Float_t[ntracks+2];
+  Float_t *maxr        = new Float_t[ntracks+2];
+  Float_t *minr        = new Float_t[ntracks+2];
+  Float_t *minPointAngle= new Float_t[ntracks+2];
+  //
+  AliESDV0MI *pvertex      = new AliESDV0MI;
+  AliITStrackMI * dummy= new AliITStrackMI;
+  dummy->SetLabel(0);
+  AliITStrackMI  trackat0;    //temporary track for DCA calculation
   //
-  //TTreeSRedirector cstream("itsv0.root");
-  Int_t centries=0;
-  AliHelix * helixes = new AliHelix[30000];
-  TObjArray trackarray(30000);
-  TObjArray trackarrayc(30000);
-  Float_t * dist = new Float_t[30000];
-  Float_t * normdist0 = new Float_t[30000];
-  Float_t * normdist1 = new Float_t[30000];
-  Float_t * normdist = new Float_t[30000];
-  Float_t * norm = new Float_t[30000];
-  AliESDV0MI  *vertexarray    = new AliESDV0MI[100000];
-  AliESDV0MI *pvertex     = &vertexarray[0];
-  AliITStrackMI * dummy=0;
+  Float_t primvertex[3]={GetX(),GetY(),GetZ()};
   //
+  // make its -  esd map
   //
-  Int_t entries = fTrackHypothesys.GetEntriesFast();
-  for (Int_t i=0;i<entries;i++){
-    TObjArray * array = (TObjArray*)fTrackHypothesys.At(i);
-    if (!array) continue;
-    // get best track without vertex constrain
-    Int_t hentries = array->GetEntriesFast();
+  for (Int_t itrack=0;itrack<ntracks+2;itrack++) {
+    itsmap[itrack]        = -1;
+    forbidden[itrack]     = kFALSE;
+    maxr[itrack]          = kMaxR;
+    minr[itrack]          = kMinR;
+    minPointAngle[itrack] = kMinPointAngle;
+  }
+  for (Int_t itrack=0;itrack<nitstracks;itrack++){
+    AliITStrackMI * original =   (AliITStrackMI*)(fOriginal.At(itrack));
+    Int_t           esdindex =   original->fESDtrack->GetID();
+    itsmap[esdindex]         =   itrack;
+  }
+  //
+  // create its tracks from esd tracks if not done before
+  //
+  for (Int_t itrack=0;itrack<ntracks;itrack++){
+    if (itsmap[itrack]>=0) continue;
+    AliITStrackMI * tpctrack = new AliITStrackMI(*(event->GetTrack(itrack)));
+    tpctrack->fD[0] = tpctrack->GetD(GetX(),GetY());
+    tpctrack->fD[1] = tpctrack->GetZat(GetX())-GetZ(); 
+    if (tpctrack->fD[0]<20 && tpctrack->fD[1]<20){
+      // tracks which can reach inner part of ITS
+      // propagate track to outer its volume - with correction for material
+      CorrectForDeadZoneMaterial(tpctrack);  
+    }
+    itsmap[itrack] = nitstracks;
+    fOriginal.AddAt(tpctrack,nitstracks);
+    nitstracks++;
+  }
+  //
+  // fill temporary arrays
+  //
+  for (Int_t itrack=0;itrack<ntracks;itrack++){
+    AliESDtrack *  esdtrack = event->GetTrack(itrack);
+    Int_t          itsindex = itsmap[itrack];
+    AliITStrackMI *original = (AliITStrackMI*)fOriginal.At(itsmap[itrack]);
+    if (!original) continue;
+    AliITStrackMI *bestConst  = 0;
+    AliITStrackMI *bestLong   = 0;
+    AliITStrackMI *best       = 0;    
+    //
     //
-    // best with vertex constrain
-    AliITStrackMI * trackc = (AliITStrackMI*)array->At(0);
-    if (trackc&&trackc->fConstrain&&trackc->fN==6&&trackc->fNormChi2[0]<2.) continue;
-    trackc=0;
+    TObjArray * array    = (TObjArray*)  fTrackHypothesys.At(itsindex);
+    Int_t       hentries = (array==0) ?  0 : array->GetEntriesFast();
+    // Get best track with vertex constrain
     for (Int_t ih=0;ih<hentries;ih++){
       AliITStrackMI * trackh = (AliITStrackMI*)array->At(ih);
       if (!trackh->fConstrain) continue;
-      if (trackh->fN<6) continue;
-      trackc = trackh;
-      if (!dummy) dummy = trackc;
-      dummy->SetLabel(0);
+      if (!bestConst) bestConst = trackh;
+      if (trackh->fN>5.0){
+       bestConst  = trackh;                         // full track -  with minimal chi2
+       break;
+      }
+      if (trackh->fN+trackh->fNDeadZone<=bestConst->fN+bestConst->fNDeadZone)  continue;      
+      bestConst = trackh;
       break;
-    }    
-    //
-    // best without vertex
-    AliITStrackMI * track = 0;
+    }
+    // Get best long track without vertex constrain and best track without vertex constrain
     for (Int_t ih=0;ih<hentries;ih++){
       AliITStrackMI * trackh = (AliITStrackMI*)array->At(ih);
       if (trackh->fConstrain) continue;
-      track = trackh;
-      break;
+      if (!best)     best     = trackh;
+      if (!bestLong) bestLong = trackh;
+      if (trackh->fN>5.0){
+       bestLong  = trackh;                         // full track -  with minimal chi2
+       break;
+      }
+      if (trackh->fN+trackh->fNDeadZone<=bestLong->fN+bestLong->fNDeadZone)  continue;      
+      bestLong = trackh;       
     }
-    if (trackc&&track){
-      if (trackc->fChi2MIP[1]<2.) continue;
-      if (trackc->fChi2MIP[0]<2. && trackc->fChi2MIP[1]<2.) continue;
-      trackarrayc.AddAt(trackc,i);
-      if (trackc->fN==6&&track->fN&&trackc->fNormChi2[0] < track->fNormChi2[0]-2) continue;
+    if (!best) {
+      best     = original;
+      bestLong = original;
     }
+    trackat0 = *bestLong;
+    Double_t xx,yy,zz,alpha; 
+    bestLong->GetGlobalXYZat(bestLong->GetX(),xx,yy,zz);
+    alpha = TMath::ATan2(yy,xx);    
+    trackat0.Propagate(alpha,0);      
+    // calculate normalized distances to the vertex 
     //
+    if ( bestLong->fN>3 ){
+      dist[itsindex]      = trackat0.fP0;
+      norm[itsindex]      = TMath::Sqrt(trackat0.fC00);
+      normdist0[itsindex] = TMath::Abs(trackat0.fP0/norm[itsindex]);
+      normdist1[itsindex] = TMath::Abs((trackat0.fP1-primvertex[2])/TMath::Sqrt(trackat0.fC11));
+      normdist[itsindex]  = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
+    }
+    else{      
+      if (bestConst&&bestConst->fN+bestConst->fNDeadZone>4.5){
+       dist[itsindex] = bestConst->fD[0];
+       norm[itsindex] = bestConst->fDnorm[0];
+       normdist0[itsindex] = TMath::Abs(bestConst->fD[0]/norm[itsindex]);
+       normdist1[itsindex] = TMath::Abs(bestConst->fD[0]/norm[itsindex]);
+       normdist[itsindex]  = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
+      }else{
+       dist[itsindex]      = trackat0.fP0;
+       norm[itsindex]      = TMath::Sqrt(trackat0.fC00);
+       normdist0[itsindex] = TMath::Abs(trackat0.fP0/norm[itsindex]);
+       normdist1[itsindex] = TMath::Abs((trackat0.fP1-primvertex[2])/TMath::Sqrt(trackat0.fC11));
+       normdist[itsindex]  = TMath::Sqrt(normdist0[itsindex]*normdist0[itsindex]+normdist1[itsindex]*normdist1[itsindex]);
+       if (TMath::Abs(trackat0.fP3)>1.1){
+         if (normdist[itsindex]<10) forbidden[itsindex]=kTRUE;
+       }
+      }
+    }
     //
+    //-----------------------------------------------------------
+    //Forbid primary track candidates - 
     //
-    if (track){
-      dist[i] = TMath::Sqrt(track->fD[0]*track->fD[0]+track->fD[1]*track->fD[1]);
-      norm[i] = track->fDnorm[0];
-      normdist0[i] = TMath::Abs(track->fD[0]/track->fDnorm[0]);
-      normdist1[i] = TMath::Abs(track->fD[1]/track->fDnorm[1]);
-      normdist[i]  = TMath::Sqrt(normdist0[i]*normdist0[i]+normdist1[i]*normdist1[i]);
-      if (track->IsGoldPrimary()) continue; //primary track
-      if (track->fD[0]<0.02 && (track->fN+track->fNDeadZone>5.8)){
-       if (normdist[i]<3.) continue;  // primary track  - cutoff 3 sigma
-       if (normdist0[i]<2.) continue; //DCA normalized cut 2 sigma
+    //treetr->SetAlias("forbidden0","Tr0.fN<4&&Tr1.fN+Tr1.fNDeadZone>4.5");
+    //treetr->SetAlias("forbidden1","ND<3&&Tr1.fN+Tr1.fNDeadZone>5.5");
+    //treetr->SetAlias("forbidden2","ND<2&&Tr1.fClIndex[0]>0&&Tr1.fClIndex[0]>0");
+    //treetr->SetAlias("forbidden3","ND<1&&Tr1.fClIndex[0]>0");
+    //treetr->SetAlias("forbidden4","ND<4&&Tr1.fNormChi2[0]<2");
+    //treetr->SetAlias("forbidden5","ND<5&&Tr1.fNormChi2[0]<1");
+    //-----------------------------------------------------------
+    if (bestConst){
+      if (bestLong->fN<4       && bestConst->fN+bestConst->fNDeadZone>4.5)               forbidden[itsindex]=kTRUE;
+      if (normdist[itsindex]<3 && bestConst->fN+bestConst->fNDeadZone>5.5)               forbidden[itsindex]=kTRUE;
+      if (normdist[itsindex]<2 && bestConst->fClIndex[0]>0 && bestConst->fClIndex[1]>0 ) forbidden[itsindex]=kTRUE;
+      if (normdist[itsindex]<1 && bestConst->fClIndex[0]>0)                              forbidden[itsindex]=kTRUE;
+      if (normdist[itsindex]<4 && bestConst->fNormChi2[0]<2)                             forbidden[itsindex]=kTRUE;
+      if (normdist[itsindex]<5 && bestConst->fNormChi2[0]<1)                             forbidden[itsindex]=kTRUE;      
+      if (bestConst->fNormChi2[0]<2.5) {
+       minPointAngle[itsindex]= 0.9999;
+       maxr[itsindex]         = 10;
       }
-      trackarray.AddAt(track,i);
-      new (&helixes[i]) AliHelix(*track);
     }
+    //
+    //forbid daughter kink candidates
+    //
+    if (esdtrack->GetKinkIndex(0)>0) forbidden[itsindex] = kTRUE;
+    Bool_t isElectron = kTRUE;
+    Double_t pid[5];
+    esdtrack->GetESDpid(pid);
+    for (Int_t i=1;i<5;i++){
+      if (pid[0]<pid[i]) isElectron= kFALSE;
+    }
+    if (isElectron){
+      forbidden[itsindex]=kFALSE;      
+      normdist[itsindex]+=4;
+    }
+    //
+    // Causality cuts in TPC volume
+    //
+    if (esdtrack->GetTPCdensity(0,10) >0.6)  maxr[itsindex] = TMath::Min(Float_t(110),maxr[itsindex]);
+    if (esdtrack->GetTPCdensity(10,30)>0.6)  maxr[itsindex] = TMath::Min(Float_t(120),maxr[itsindex]);
+    if (esdtrack->GetTPCdensity(20,40)>0.6)  maxr[itsindex] = TMath::Min(Float_t(130),maxr[itsindex]);
+    if (esdtrack->GetTPCdensity(30,50)>0.6)  maxr[itsindex] = TMath::Min(Float_t(140),maxr[itsindex]);
+    //
+    if (esdtrack->GetTPCdensity(0,60)<0.4&&bestLong->fN<3) minr[itsindex]=100;    
+    //
+    //
+    if (0){
+      cstream<<"Track"<<
+       "Tr0.="<<best<<
+       "Tr1.="<<((bestConst)? bestConst:dummy)<<
+       "Tr2.="<<bestLong<<
+       "Tr3.="<<&trackat0<<
+       "Esd.="<<esdtrack<<
+       "Dist="<<dist[itsindex]<<
+       "ND0="<<normdist0[itsindex]<<
+       "ND1="<<normdist1[itsindex]<<
+       "ND="<<normdist[itsindex]<<
+       "Pz="<<primvertex[2]<<
+       "Forbid="<<forbidden[itsindex]<<
+       "\n";
+      //
+    }
+    trackarray.AddAt(best,itsindex);
+    trackarrayc.AddAt(bestConst,itsindex);
+    trackarrayl.AddAt(bestLong,itsindex);
+    new (&helixes[itsindex]) AliHelix(*best);
   }
   //
   //
-  //  Int_t multifound=0;
-  Int_t vertexall =0;
-  AliESDV0MI tempvertex;
-  Float_t primvertex[3]={GetX(),GetY(),GetZ()};
-  
-  
-  for (Int_t itrack0=0;itrack0<entries;itrack0++){
-    //
-    AliITStrackMI * track0 = (AliITStrackMI*)trackarray.At(itrack0);
-    if (!track0) continue;
-    if (track0->fP4>0) continue;
+  //
+  // first iterration of V0 finder
+  //
+  for (Int_t iesd0=0;iesd0<ntracks;iesd0++){
+    Int_t itrack0 = itsmap[iesd0];
+    if (forbidden[itrack0]) continue;
+    AliITStrackMI * btrack0 = (AliITStrackMI*)trackarray.At(itrack0);
+    if (!btrack0) continue;    
+    if (btrack0->fP4>0) continue;
     AliITStrackMI *trackc0 = (AliITStrackMI*)trackarrayc.At(itrack0);
     //
-    TObjArray * array0     = (TObjArray*)fTrackHypothesys.At(itrack0);
-    //
-    Int_t vertexes =0;
-    for (Int_t itrack1=0;itrack1<entries;itrack1++){
-      AliITStrackMI * track1 = (AliITStrackMI*)trackarray.At(itrack1); 
-      if (!track1) continue;
-      if (track1->fP4<0) continue;
-      AliITStrackMI *trackc1 = (AliITStrackMI*)trackarrayc.At(itrack1);
-      if (trackc0&&trackc1){
-       if (TMath::Min(trackc0->fChi2MIP[1],trackc1->fChi2MIP[1])<2.) continue;
+    for (Int_t iesd1=0;iesd1<ntracks;iesd1++){
+      Int_t itrack1 = itsmap[iesd1];
+      if (forbidden[itrack1]) continue;
+
+      AliITStrackMI * btrack1 = (AliITStrackMI*)trackarray.At(itrack1); 
+      if (!btrack1) continue;
+      if (btrack1->fP4<0) continue;
+      Bool_t isGold = kFALSE;
+      if (TMath::Abs(TMath::Abs(btrack0->GetLabel())-TMath::Abs(btrack1->GetLabel()))==1){
+       isGold = kTRUE;
       }
-      if (track1->fNDeadZone+track0->fNDeadZone>1.1) continue;
-      TObjArray * array1     = (TObjArray*)fTrackHypothesys.At(itrack1);
+      AliITStrackMI *trackc1 = (AliITStrackMI*)trackarrayc.At(itrack1);
+      AliHelix &h1 = helixes[itrack0];
+      AliHelix &h2 = helixes[itrack1];
       //
-      //if (normdist0[itrack0]+normdist0[itrack1]<3) continue;
-      //if (normdist[itrack0]+normdist[itrack1]<4) continue;
+      // find linear distance
+      Double_t rmin =0;
       //
       //
-      AliHelix *h1 = &helixes[itrack0];
-      AliHelix *h2 = &helixes[itrack1];
-      Double_t rmin =0;
-      Double_t distance = TestV0(h1,h2,pvertex,rmin);
       //
-      if (distance>0.4) continue;
-      if (pvertex->GetRr()<0.3) continue;
-      if (pvertex->GetRr()>20.) continue;
-      pvertex->SetM(*track0);
-      pvertex->SetP(*track1);
-      pvertex->Update(primvertex);
-      if (pvertex->GetRr()<0.3) continue;
-      if (pvertex->GetRr()>20.) continue;
-      if (track1->fNDeadZone+track0->fNDeadZone>0.5 &&distance>0.12) continue;
-
+      Double_t phase[2][2],radius[2];
+      Int_t  points = h1.GetRPHIintersections(h2, phase, radius);
+      if    (points==0)  continue;
+      Double_t delta[2]={1000,1000};        
+      rmin = radius[0];
+      h1.ParabolicDCA(h2,phase[0][0],phase[0][1],radius[0],delta[0]);
+      if (points==2){    
+       if (radius[1]<rmin) rmin = radius[1];
+       h1.ParabolicDCA(h2,phase[1][0],phase[1][1],radius[1],delta[1]);
+      }
+      rmin = TMath::Sqrt(rmin);
+      Double_t distance = 0;
+      Double_t radiusC  = 0;
+      Int_t    iphase   = 0;
+      if (delta[0]<delta[1]){
+       distance = TMath::Sqrt(delta[0]);
+       radiusC  = TMath::Sqrt(radius[0]);
+      }else{
+       distance = TMath::Sqrt(delta[1]);
+       radiusC  = TMath::Sqrt(radius[1]);
+       iphase=1;
+      }
+      if (radiusC<TMath::Max(minr[itrack0],minr[itrack1]))    continue;
+      if (radiusC>TMath::Min(maxr[itrack0],maxr[itrack1]))     continue; 
+      Float_t maxDist  = TMath::Min(kMaxDist,Float_t(kMaxDist0+radiusC*kMaxDist1));      
+      if (distance>maxDist) continue;
+      Float_t pointAngle = h1.GetPointAngle(h2,phase[iphase],primvertex);
+      if (pointAngle<TMath::Max(minPointAngle[itrack0],minPointAngle[itrack1])) continue;
       //
-
-      if ( TMath::Abs((TMath::Abs(track0->GetLabel())-TMath::Abs(track1->GetLabel())))<2 
-          ||(centries<5000&&(pvertex->GetPointAngle()>0.95))){ 
-       //cstream<<"Iter0"<<track0<<track1<<pvertex<<normdist[itrack0]<<normdist[itrack1]<<"\n";
-       centries++;
-      }
       //
+      //      Double_t distance = TestV0(h1,h2,pvertex,rmin);
       //
-      if (pvertex->GetPointAngle()<0.85) continue;
-      //      if (normdist[itrack0]+normdist[itrack1]<6&&pvertex->GetPointAngle()<0.99) continue;
+      //       if (distance>maxDist)           continue;
+      //       if (pvertex->GetRr()<kMinR)     continue;
+      //       if (pvertex->GetRr()>kMaxR)     continue;
+      AliITStrackMI * track0=btrack0;
+      AliITStrackMI * track1=btrack1;
+      //      if (pvertex->GetRr()<3.5){  
+      if (radiusC<3.5){  
+       //use longest tracks inside the pipe
+       track0 = (AliITStrackMI*)trackarrayl.At(itrack0);
+       track1 = (AliITStrackMI*)trackarrayl.At(itrack1);
+      }      
       //
       //
+      pvertex->SetM(*track0);
+      pvertex->SetP(*track1);
+      pvertex->Update(primvertex);
+      if (pvertex->GetRr()<kMinR) continue;
+      if (pvertex->GetRr()>kMaxR) continue;
+      if (pvertex->GetPointAngle()<kMinPointAngle) continue;
+      if (pvertex->GetDist2()>maxDist) continue;
       pvertex->SetLab(0,track0->GetLabel());
       pvertex->SetLab(1,track1->GetLabel());
-      pvertex->SetIndex(0,track0->GetESDtrack()->GetID());
-      pvertex->SetIndex(1,track1->GetESDtrack()->GetID());
-      // calculate chi2s
-      //
-      pvertex->SetChi2After(0);
-      pvertex->SetChi2Before(0);
-      pvertex->SetNBefore(0);
-      pvertex->SetNAfter(0);
       
-      for (Int_t i=0;i<6;i++){
-       Float_t radius = fgLayers[i].GetR();
-       if (pvertex->GetRr()>radius+0.5){
-         pvertex->SetNBefore(pvertex->GetNBefore()+2.);
-         //
-         if (track0->fClIndex[i]<=0) {
-           pvertex->SetChi2Before(pvertex->GetChi2Before()+9);
-         }else{
-           Float_t chi2 = track0->fDy[i]*track0->fDy[i]/(track0->fSigmaY[i]*track0->fSigmaY[i])+
-             track0->fDz[i]*track0->fDz[i]/(track0->fSigmaZ[i]*track0->fSigmaZ[i]);
-           pvertex->SetChi2Before(pvertex->GetChi2Before()+chi2);
-         }
+      //      
+      AliITStrackMI * htrackc0 = trackc0 ? trackc0:dummy;      
+      AliITStrackMI * htrackc1 = trackc1 ? trackc1:dummy;
 
-         if (track1->fClIndex[i]<=0) {
-           pvertex->SetChi2Before(pvertex->GetChi2Before()+9);
-
-         }else{
-           Float_t chi2 = track1->fDy[i]*track1->fDy[i]/(track1->fSigmaY[i]*track1->fSigmaY[i])+
-             track1->fDz[i]*track1->fDz[i]/(track1->fSigmaZ[i]*track1->fSigmaZ[i]);
-           //        pvertex->fChi2Before+=chi2;
-           pvertex->SetChi2Before(pvertex->GetChi2Before()+chi2);
+      //
+      //
+      TObjArray * array0b     = (TObjArray*)fBestHypothesys.At(itrack0);
+      if (!array0b&&pvertex->GetRr()<40 && TMath::Abs(track0->fP3)<1.1) 
+       FollowProlongationTree((AliITStrackMI*)fOriginal.At(itrack0),itrack0, kFALSE);
+      TObjArray * array1b    = (TObjArray*)fBestHypothesys.At(itrack1);
+      if (!array1b&&pvertex->GetRr()<40 && TMath::Abs(track1->fP3)<1.1) 
+       FollowProlongationTree((AliITStrackMI*)fOriginal.At(itrack1),itrack1, kFALSE);
+      //
+      AliITStrackMI * track0b = (AliITStrackMI*)fOriginal.At(itrack0);       
+      AliITStrackMI * track1b = (AliITStrackMI*)fOriginal.At(itrack1);
+      AliITStrackMI * track0l = (AliITStrackMI*)fOriginal.At(itrack0);       
+      AliITStrackMI * track1l = (AliITStrackMI*)fOriginal.At(itrack1);
+      
+      Float_t minchi2before0=16;
+      Float_t minchi2before1=16;
+      Float_t minchi2after0 =16;
+      Float_t minchi2after1 =16;
+      Int_t maxLayer = GetNearestLayer(pvertex->GetXrp());
+      
+      if (array0b) for (Int_t i=0;i<5;i++){
+       // best track after vertex
+       AliITStrackMI * btrack = (AliITStrackMI*)array0b->At(i);
+       if (!btrack) continue;
+       if (btrack->fN>track0l->fN) track0l = btrack;     
+       if (btrack->fX<pvertex->GetRr()-2) {
+         if (maxLayer>i+2 && btrack->fN==(6-i)&&i<2){
+           Float_t sumchi2= 0;
+           Float_t sumn   = 0;
+           for (Int_t ilayer=i;ilayer<maxLayer;ilayer++){
+             sumn+=1.;
+             if (!btrack->fClIndex[ilayer]){
+               sumchi2+=25;
+               continue;
+             }else{
+               sumchi2+=btrack->fDy[ilayer]*btrack->fDy[ilayer]/(btrack->fSigmaY[ilayer]*btrack->fSigmaY[ilayer]);
+               sumchi2+=btrack->fDz[ilayer]*btrack->fDz[ilayer]/(btrack->fSigmaZ[ilayer]*btrack->fSigmaZ[ilayer]);            
+             }
+           }
+           sumchi2/=sumn;
+           if (sumchi2<minchi2before0) minchi2before0=sumchi2; 
          }
+         continue;   //safety space - Geo manager will give exact layer
        }
-
-       if (pvertex->GetRr()<radius-0.5){
-         pvertex->SetNAfter(pvertex->GetNAfter()+2.);
-         //
-         if (track0->fClIndex[i]<=0) {
-           pvertex->SetChi2After(pvertex->GetChi2After()+9);
-         }else{
-           Float_t chi2 = track0->fDy[i]*track0->fDy[i]/(track0->fSigmaY[i]*track0->fSigmaY[i])+
-             track0->fDz[i]*track0->fDz[i]/(track0->fSigmaZ[i]*track0->fSigmaZ[i]);
-           pvertex->SetChi2After(pvertex->GetChi2After()+chi2);
-         }
-
-         if (track1->fClIndex[i]<=0) {
-           pvertex->SetChi2After(pvertex->GetChi2After()+9.);
-         }else{
-           Float_t chi2 = track1->fDy[i]*track1->fDy[i]/(track1->fSigmaY[i]*track1->fSigmaY[i])+
-             track1->fDz[i]*track1->fDz[i]/(track1->fSigmaZ[i]*track1->fSigmaZ[i]);
-           pvertex->SetChi2After(pvertex->GetChi2After()+chi2);
+       track0b       = btrack;
+       minchi2after0 = btrack->fNormChi2[i];
+       break;
+      }
+      if (array1b) for (Int_t i=0;i<5;i++){
+       // best track after vertex
+       AliITStrackMI * btrack = (AliITStrackMI*)array1b->At(i);
+       if (!btrack) continue;
+       if (btrack->fN>track1l->fN) track1l = btrack;     
+       if (btrack->fX<pvertex->GetRr()-2){
+         if (maxLayer>i+2&&btrack->fN==(6-i)&&(i<2)){
+           Float_t sumchi2= 0;
+           Float_t sumn   = 0;
+           for (Int_t ilayer=i;ilayer<maxLayer;ilayer++){
+             sumn+=1.;
+             if (!btrack->fClIndex[ilayer]){
+               sumchi2+=30;
+               continue;
+             }else{
+               sumchi2+=btrack->fDy[ilayer]*btrack->fDy[ilayer]/(btrack->fSigmaY[ilayer]*btrack->fSigmaY[ilayer]);
+               sumchi2+=btrack->fDz[ilayer]*btrack->fDz[ilayer]/(btrack->fSigmaZ[ilayer]*btrack->fSigmaZ[ilayer]);            
+             }
+           }
+           sumchi2/=sumn;
+           if (sumchi2<minchi2before1) minchi2before1=sumchi2; 
          }
+         continue;   //safety space - Geo manager will give exact layer           
        }
+       track1b = btrack;
+       minchi2after1 = btrack->fNormChi2[i];
+       break;
       }
-      if (pvertex->GetNBefore()>2){
-       if (pvertex->GetChi2Before()/pvertex->GetNBefore()<4.) continue; //have clusters before vetex
-      }
-      Int_t ibest0=0,ibest1=0;
-      AliITStrackMI * ntrack0 = track0;
-      AliITStrackMI * ntrack1 = track1; 
-      //
-      //
-      //PH      Float_t oldistance = pvertex->GetDist2();
-      Bool_t improve  = FindBestPair(itrack0,itrack1,pvertex,ibest0,ibest1);   // try to improve vertex
-      if (pvertex->GetPointAngle()<0.5) continue;
-      distance = pvertex->GetDist2(); 
-      if (improve){
-       ntrack0 = (AliITStrackMI*)array0->At(ibest0);
-       ntrack1 = (AliITStrackMI*)array1->At(ibest1); 
-      }
-      Bool_t accept0 = kFALSE;      // accept ==>  because of pointing angle 
-      if (pvertex->GetPointAngle()>0.999){
-       if (pvertex->GetRr()<3.5 && (ntrack0->fN+ntrack0->fNDeadZone+ntrack1->fN+ntrack1->fNDeadZone)<11.5) continue;
-       if (pvertex->GetRr()>3.5&& pvertex->GetDistNorm()<12) accept0 = kTRUE;
-       if (pvertex->GetRr()>1 && pvertex->GetDist2()<0.1 && pvertex->GetDistNorm()<12) accept0 = kTRUE;
-       if (pvertex->GetPointAngle()>0.9995&&pvertex->GetRr()>5.) accept0 = kTRUE;
-      }
-      Bool_t reject1= kFALSE;      // reject ==> bad kinematic
-      //
-      reject1 |=  TMath::Abs(ntrack0->fN+ntrack0->fNDeadZone-ntrack1->fN-ntrack1->fNDeadZone)>1.02 || 
-                  TMath::Abs(ntrack0->fN-ntrack1->fN)>1.02;  // cut1
-      reject1 |=  ntrack0->fNUsed+ntrack1->fNUsed>1.01;                                 // cut2
-      reject1 |=  pvertex->GetDistNorm()>12;                                                // cut3
-      reject1 |=  pvertex->GetDist2()>0.1 && improve;                                       // cut4
-      reject1 |=  (TMath::Abs(ntrack0->fD[0])+TMath::Abs(ntrack1->fD[0]))/pvertex->GetDist2()<5; //cut5
-      reject1 |=  TMath::Abs(ntrack0->fD[0]/pvertex->GetDist2())<2 || TMath::Abs(ntrack1->fD[0]/pvertex->GetDist2())<2;  //cut 6
       //
-      // small radii cuts  
-      Bool_t reject2 = kFALSE;
-      if (pvertex->GetRr()<3.6){
-       reject2 |=  (TMath::Abs(ntrack0->fN+ntrack0->fNDeadZone-ntrack1->fN-ntrack1->fNDeadZone)>0.01);  // cut7
-       reject2 |=  ntrack0->fNUsed+ntrack1->fNUsed>0.01;                                  // cut8
-       reject2 |=  ntrack0->fN+ntrack0->fNDeadZone+ntrack1->fN+ntrack1->fNDeadZone<11.5;  // cut9
-       reject2 |=  (ntrack0->fN+ntrack1->fN<11.5)&&pvertex->GetRr()<2;                        // cut10
-       reject2 |=  pvertex->GetDist2()>0.1;                                                   // cut11   
-      }        
-      //PH      AliITStrackMI * htrackc0 = trackc0 ? trackc0:dummy;
-      //PH      AliITStrackMI * htrackc1 = trackc1 ? trackc1:dummy;
+      // position resolution - used for DCA cut
+      Float_t sigmad = track0b->fC00+track0b->fC11+track1b->fC00+track1b->fC11+
+       (track0b->fX-pvertex->GetRr())*(track0b->fX-pvertex->GetRr())*(track0b->fC22+track0b->fC33)+
+       (track1b->fX-pvertex->GetRr())*(track1b->fX-pvertex->GetRr())*(track1b->fC22+track1b->fC33);
+      sigmad =TMath::Sqrt(sigmad)+0.04;
+      if (pvertex->GetRr()>50){
+       Double_t cov0[15],cov1[15];
+       track0b->fESDtrack->GetInnerExternalCovariance(cov0);
+       track1b->fESDtrack->GetInnerExternalCovariance(cov1);
+       sigmad = cov0[0]+cov0[2]+cov1[0]+cov1[2]+
+         (80.-pvertex->GetRr())*(80.-pvertex->GetRr())*(cov0[5]+cov0[9])+
+         (80.-pvertex->GetRr())*(80.-pvertex->GetRr())*(cov1[5]+cov1[9]);
+       sigmad =TMath::Sqrt(sigmad)+0.05;
+      }
+      //       
+      AliESDV0MI vertex2;
+      vertex2.SetM(*track0b);
+      vertex2.SetP(*track1b);
+      vertex2.Update(primvertex);
+      if (vertex2.GetDist2()<=pvertex->GetDist2()&&(vertex2.GetPointAngle()>=pvertex->GetPointAngle())){
+       pvertex->SetM(*track0b);
+       pvertex->SetP(*track1b);
+       pvertex->Update(primvertex);
+      }
+      pvertex->SetDistSigma(sigmad);
+      pvertex->SetDistNorm(pvertex->GetDist2()/sigmad);       
       //
+      // define likelihhod and causalities
       //
-      //
-      if ( TMath::Abs((TMath::Abs(track0->GetLabel())-TMath::Abs(track1->GetLabel())))<2 
-          ||(centries<500000)){        
-       /*
-       cstream<<"It1"<<"Tr0.="<<ntrack0<<"TR1.="<<ntrack1<<"V0.="<<pvertex<<"ND0.="<<normdist[itrack0]<<"ND1.="<<
-         normdist[itrack1]<<"D.="<<distance<<"DistOld="<<oldistance<<"Imp.="<<improve<<
-         "A0="<<accept0<<"R1="<<reject1<<"R2="<<reject2<<"Rmin.="<<rmin<<
-         "TrC0.="<<htrackc0<<"TRC1.="<<htrackc1<<"\n";
-       */
-       centries++;
+      Float_t pa0=1, pa1=1, pb0=0.26, pb1=0.26;      
+      if (maxLayer<2){
+       if (pvertex->GetAnglep()[2]>0.2){
+         pb0    =  TMath::Exp(-TMath::Min(normdist[itrack0],Float_t(16.))/12.);
+         pb1    =  TMath::Exp(-TMath::Min(normdist[itrack1],Float_t(16.))/12.);
+       }
+       pvertex->SetChi2Before(normdist[itrack0]);
+       pvertex->SetChi2After(normdist[itrack1]);       
+       pvertex->SetNAfter(0);
+       pvertex->SetNBefore(0);
+      }else{
+       pvertex->SetChi2Before(minchi2before0);
+       pvertex->SetChi2After(minchi2before1);
+       if (pvertex->GetAnglep()[2]>0.2){
+         pb0    =  TMath::Exp(-TMath::Min(minchi2before0,Float_t(16))/12.);
+         pb1    =  TMath::Exp(-TMath::Min(minchi2before1,float_t(16))/12.);
+       }
+       pvertex->SetNAfter(maxLayer);
+       pvertex->SetNBefore(maxLayer);
       }
-  
-      if (!accept0 && (reject1 || reject2))  continue;
-
-//       if (distance>0.5) continue; 
-//       distance = pvertex->GetDist2();
-//       if (pvertex->GetRr()>25 || pvertex->GetRr()<0.2) continue;
-//       if (pvertex->GetRr()/pvertex->fDistSigma<1) continue;
-//       if (pvertex->GetDistNorm()>10) continue;
-//       if (pvertex->GetPointAngle()<0.85) continue;  
-//       if ((normdist[itrack0]<3||normdist[itrack1]<3)){
-//     if (pvertex->GetPointAngle()<0.99||pvertex->GetDist2()>0.15) continue;
-//       }
-//       if (distance>0.05*(0.8+0.2*(0.5+pvertex->GetRr()))) continue;
-//       if (pvertex->GetRr()<0.3) continue;
-//       if (pvertex->GetRr()>27.) continue; 
-
-
-
-      //
-      if (distance<0.3 &&pvertex->GetPointAngle()>0.998){
-       track0->fGoldV0 = kTRUE;
-       track1->fGoldV0 = kTRUE;
+      if (pvertex->GetRr()<90){
+       pa0  *= TMath::Min(track0->fESDtrack->GetTPCdensity(0,60),Float_t(1.));
+       pa1  *= TMath::Min(track1->fESDtrack->GetTPCdensity(0,60),Float_t(1.));
       }
-      vertexes++;
-      vertexall++;
-      if (vertexall>=100000) break;
-      pvertex = &vertexarray[vertexall];
-    }
-  }
-  //  printf("\n\n\nMultifound\t%d\n\n\n",multifound);
-  //
-  // sort vertices according quality
-  Float_t quality[40000];
-  Int_t   indexes[40000];
-  Int_t   trackvertices[40000];
-  for (Int_t i=0;i<entries;i++) trackvertices[i]=0;
-  for (Int_t i=0;i<vertexall;i++) {
-    Float_t norm     = 1.-0.999*TMath::Abs(vertexarray[i].GetPointAngle());
-    Float_t fnormdist = 1./(1+vertexarray[i].GetRr());
-    quality[i] = norm*fnormdist;
-  }
-  //
-  TMath::Sort(vertexall,quality,indexes,kFALSE);
-  
-  for (Int_t i=0;i<vertexall;i++){
-    pvertex= &vertexarray[indexes[i]];
-    Int_t index0 = vertexarray[indexes[i]].GetIndex(0);
-    Int_t index1 = vertexarray[indexes[i]].GetIndex(1);
-    vertexarray[indexes[i]].SetOrder(2,i);
-    vertexarray[indexes[i]].SetOrder(1,trackvertices[index1]);
-    vertexarray[indexes[i]].SetOrder(0,trackvertices[index0]);
-    Int_t v0index = event->AddV0MI(&vertexarray[indexes[i]]);
-    //
-    if (trackvertices[index1]+trackvertices[index0]>5) continue;
-    if (trackvertices[index0]>2) continue;
-    if (trackvertices[index1]>2) continue;
-
-    if (trackvertices[index1]+trackvertices[index0]>0) {
-      //      if (pvertex->GetPointAngle()<0.995)  continue;
-    }
-    trackvertices[index0]++;
-    trackvertices[index1]++;
-    
-    AliESDtrack * ptrack0 = event->GetTrack(vertexarray[indexes[i]].GetIndex(0));
-    AliESDtrack * ptrack1 = event->GetTrack(vertexarray[indexes[i]].GetIndex(1));
-    if (!ptrack0 || !ptrack1){
-      printf("BBBBBBBUUUUUUUUUUGGGGGGGGGG\n");
-    }
-    Int_t v0index0[3]={ptrack0->GetV0Index(0),ptrack0->GetV0Index(1),ptrack0->GetV0Index(2)};
-    Int_t v0index1[3]={ptrack1->GetV0Index(0),ptrack1->GetV0Index(1),ptrack1->GetV0Index(2)};
-    for (Int_t i=0;i<3;i++){
-      if (v0index0[i]<0) {
-       v0index0[i]=v0index;
-       ptrack0->SetV0Indexes(v0index0);
-       break;
-      }
-    }
-    for (Int_t i=0;i<3;i++){
-      if (v0index1[i]<0) {
-       v0index1[i]=v0index;
-       ptrack1->SetV0Indexes(v0index1);
-       break;
+      if (pvertex->GetRr()<20){
+       pa0  *= (0.2+TMath::Exp(-TMath::Min(minchi2after0,Float_t(16))/8.))/1.2;
+       pa1  *= (0.2+TMath::Exp(-TMath::Min(minchi2after1,Float_t(16))/8.))/1.2;
       }
-    }
-  }
-
-  delete [] helixes;
-  delete [] dist;
-  delete [] normdist0;
-  delete [] normdist1;
-  delete [] normdist;
-  delete [] norm;
-
-  delete[] vertexarray;
-}
-
-
-
-Bool_t  AliITStrackerMI::FindBestPair(Int_t esdtrack0, Int_t esdtrack1, AliESDV0MI *vertex, Int_t &i0, Int_t &i1)
-{
-  //
-  // try to find best pair from the tree of track hyp.
-  //
-  TObjArray * array0 = (TObjArray*)fTrackHypothesys.At(esdtrack0);
-  Int_t entries0 = array0->GetEntriesFast();
-  TObjArray * array1 = (TObjArray*)fTrackHypothesys.At(esdtrack1);
-  Int_t entries1 = array1->GetEntriesFast();  
-  //  AliITStrackMI *orig0 = (AliITStrackMI*)fOriginal.At(esdtrack0);
-  //AliITStrackMI *orig1 = (AliITStrackMI*)fOriginal.At(esdtrack1);
-  Double_t criticalradius = vertex->GetRr();
-  AliITStrackMI * track0= 0;
-  AliITStrackMI * track1= 0;
-  i0 = -1;
-  i1 = -1;
-  //
-  //
-  Float_t rfirst0[2000];  //radius position of  the first cluster - track0
-  Float_t rfirst1[2000];  //                                      - track1
-  Float_t maxlocalx0=0;      //local x  for first  track        
-  Float_t maxlocalx1=0;      //local x  for second track       
-  Float_t cs0=1, sn0=0;            //rotations
-  Float_t cs1=1, sn1=0;            //rotations
-
-  //
-  for (Int_t itrack0=0;itrack0<entries0;itrack0++){
-    rfirst0[itrack0]=-1.;
-    AliITStrackMI * htrack0 = (AliITStrackMI*)array0->At(itrack0);
-    if (!htrack0) continue;
-    if (htrack0->fConstrain) continue;    
-    if (i0<0){
-      i0     = itrack0;
-      track0 = htrack0; 
-    }
-    Double_t cs = TMath::Cos(htrack0->fAlpha);
-    Double_t sn = TMath::Sin(htrack0->fAlpha);
-    Double_t x  = htrack0->fX*cs - htrack0->fP0*sn;
-    Double_t y  = htrack0->fX*sn + htrack0->fP0*cs;
-    Double_t radius = TMath::Sqrt(x*x+y*y);
-    if (criticalradius<3  && radius>6&&htrack0->fNDeadZone<0.2) continue;  // all cluster required
-    if (criticalradius>10 && radius<6) continue;  // causality
-    Double_t localx =  TMath::Abs(vertex->GetXr(0)*cs + vertex->GetXr(1)*sn);
-    if (localx>maxlocalx0) {
-      maxlocalx0=localx;
-      cs0 = cs; sn0=sn;
-    }
-    rfirst0[itrack0] = radius;
-  }
-  for (Int_t itrack1=0;itrack1<entries1;itrack1++){
-    rfirst1[itrack1]=-1.;
-    AliITStrackMI * htrack1 = (AliITStrackMI*)array1->At(itrack1);
-    if (!htrack1) continue;
-    if (htrack1->fConstrain) continue;    
-    if (i1<0){
-      i1     = itrack1;
-      track1 = htrack1; 
-    }
-    Double_t cs = TMath::Cos(htrack1->fAlpha);
-    Double_t sn = TMath::Sin(htrack1->fAlpha);
-    //
-    //
-    Double_t x  = htrack1->fX*cs - htrack1->fP0*sn;
-    Double_t y  = htrack1->fX*sn + htrack1->fP0*cs;
-    Double_t radius = TMath::Sqrt(x*x+y*y);
-    if (criticalradius<3  && radius>6&&htrack1->fNDeadZone<0.2) continue; //all clusters required
-    if (criticalradius>10 && radius<6) continue; //causality
-    Double_t localx =  TMath::Abs(vertex->GetXr(0)*cs + vertex->GetXr(1)*sn);
-    if (localx>maxlocalx1) {
-      maxlocalx1=localx;
-      cs1 = cs; sn1=sn;
-    }
-    rfirst1[itrack1] = radius;
-  }
-  //
-  //
-  //
-  const Float_t radiuses[4]={4,6.5,15.03,24.};
-  //
-  //
-  // find the best tracks after decay point
-  Float_t bestquality =100000;
-  Float_t bestradius=0;
-  AliESDV0MI v0;
-  Double_t vpos[3];
-  Float_t v[3]={GetX(),GetY(),GetZ()};  
-  //
-  //
-  for (Int_t itrack0=0;itrack0<entries0;itrack0++){
-    if (rfirst0[itrack0]<0) continue;
-    AliITStrackMI * htrack0 = (AliITStrackMI*)array0->At(itrack0);
-    if (!htrack0) continue;
-    //
-    for (Int_t itrack1=0;itrack1<entries1;itrack1++){
-      if (rfirst1[itrack1]<0) continue;
-      AliITStrackMI * htrack1 = (AliITStrackMI*)array1->At(itrack1);
-      if (!htrack1) continue;
-      if (TMath::Abs(rfirst0[itrack0]-rfirst1[itrack1])>1.) continue;
-      if (htrack0->fClIndex[6-htrack0->fN]==htrack1->fClIndex[6-htrack1->fN]) continue; //sharing of last cluster not allowe
       //
-      if (htrack0->fNUsed+htrack1->fNUsed>0) continue; //sharing of clusters not allowed
+      pvertex->SetCausality(pb0,pb1,pa0,pa1);
       //
-      //    
-      v0.SetM(*htrack0);
-      v0.SetP(*htrack1);
-      //      if (v0.Update(v)==0) continue;
-      v0.Update(v);
-      if (TMath::Min(rfirst0[itrack0],rfirst1[itrack1]) <v0.GetRr()-0.3) continue;
+      //  Likelihood calculations  - apply cuts
+      //         
+      Bool_t v0OK = kTRUE;
+      Float_t p12 = pvertex->GetParamP()->GetParameter()[4]*pvertex->GetParamP()->GetParameter()[4];
+      p12        += pvertex->GetParamM()->GetParameter()[4]*pvertex->GetParamM()->GetParameter()[4];
+      p12         = TMath::Sqrt(p12);                             // "mean" momenta
+      Float_t    sigmap0   = 0.0001+0.001/(0.1+pvertex->GetRr()); 
+      Float_t    sigmap    = 0.5*sigmap0*(0.6+0.4*p12);           // "resolution: of point angle - as a function of radius and momenta
+
+      Float_t CausalityA  = (1.0-pvertex->GetCausalityP()[0])*(1.0-pvertex->GetCausalityP()[1]);
+      Float_t CausalityB  = TMath::Sqrt(TMath::Min(pvertex->GetCausalityP()[2],Float_t(0.7))*
+                                       TMath::Min(pvertex->GetCausalityP()[3],Float_t(0.7)));
       //
+      Float_t Likelihood0 = (TMath::Exp(-pvertex->GetDistNorm())+0.1) *(pvertex->GetDist2()<0.5)*(pvertex->GetDistNorm()<5);
+
+      Float_t Likelihood1 = TMath::Exp(-(1.0001-pvertex->GetPointAngle())/sigmap)+
+       0.4*TMath::Exp(-(1.0001-pvertex->GetPointAngle())/(4.*sigmap))+
+       0.4*TMath::Exp(-(1.0001-pvertex->GetPointAngle())/(8.*sigmap))+
+       0.1*TMath::Exp(-(1.0001-pvertex->GetPointAngle())/0.01);
       //
-      if (v0.GetDist2()>0.3) continue;
-      if (v0.GetRr()<radiuses[1] && ( TMath::Abs(htrack0->fN+htrack0->fNDeadZone-htrack1->fN-htrack1->fNDeadZone)>0.5)) continue;
+      if (CausalityA<kCausality0Cut)                                          v0OK = kFALSE;
+      if (TMath::Sqrt(Likelihood0*Likelihood1)<kLikelihood01Cut)              v0OK = kFALSE;
+      if (Likelihood1<kLikelihood1Cut)                                        v0OK = kFALSE;
+      if (TMath::Power(Likelihood0*Likelihood1*CausalityB,0.33)<kCombinedCut) v0OK = kFALSE;
+      
       //
-      //if (v0.GetRr()<3. && (htrack0->fN+htrack0->fNDeadZone+htrack1->fN+htrack1->fNDeadZone)<11.7) continue;
-      //if (v0.GetRr()<6. && (htrack0->fN+htrack0->fNDeadZone+htrack1->fN+htrack1->fNDeadZone)<9.7) continue;
-      if (v0.GetRr()<3. && (htrack0->fN+htrack1->fN)<11.7) continue;
-      if (v0.GetRr()<6. && (htrack0->fN+htrack1->fN)<9.7) continue;
-      Double_t localx0=v0.GetXr(0)*cs0+v0.GetXr(1)*sn0;
-      Double_t localx1=v0.GetXr(0)*cs1+v0.GetXr(1)*sn1;
-      Double_t maxlocalx = TMath::Max(localx0,localx1);
-      if (maxlocalx<3.4  && (htrack0->fN+htrack1->fN)<11.7) continue;
-      if (maxlocalx<6.1  && (htrack0->fN+htrack1->fN)<9.7) continue;
       //
-      Float_t fnormdist = v0.GetDist2()/0.05;      
-      fnormdist +=(htrack0->fNormChi2[6-htrack0->fN]+htrack1->fNormChi2[6-htrack1->fN]);
-      fnormdist +=3*(htrack0->fNUsed+htrack1->fNUsed);
-      if (TMath::Min(rfirst0[itrack0],rfirst1[itrack1]) <v0.GetRr()+0.2){
-       fnormdist +=(v0.GetRr()+0.2-TMath::Min(rfirst0[itrack0],rfirst1[itrack1]))/0.1;
-      }
+      if (kTRUE){      
+       Bool_t gold = TMath::Abs(TMath::Abs(track0->GetLabel())-TMath::Abs(track1->GetLabel()))==1;
+       cstream<<"It0"<<
+         "Tr0.="<<track0<<                       //best without constrain
+         "Tr1.="<<track1<<                       //best without constrain  
+         "Tr0B.="<<track0b<<                     //best without constrain  after vertex
+         "Tr1B.="<<track1b<<                     //best without constrain  after vertex 
+         "Tr0C.="<<htrackc0<<                    //best with constrain     if exist
+         "Tr1C.="<<htrackc1<<                    //best with constrain     if exist
+         "Tr0L.="<<track0l<<                     //longest best           
+         "Tr1L.="<<track1l<<                     //longest best
+         "Esd0.="<<track0->fESDtrack<<           // esd track0 params
+         "Esd1.="<<track1->fESDtrack<<           // esd track1 params
+         "V0.="<<pvertex<<                       //vertex properties
+         "V0b.="<<&vertex2<<                       //vertex properties at "best" track
+         "ND0="<<normdist[itrack0]<<             //normalize distance for track0
+         "ND1="<<normdist[itrack1]<<             //normalize distance for track1
+         "Gold="<<gold<<                         //
+         //      "RejectBase="<<rejectBase<<             //rejection in First itteration
+         "OK="<<v0OK<<
+         "rmin="<<rmin<<
+         "sigmad="<<sigmad<<
+         "\n";
+      }      
+      //if (rejectBase) continue;
       //
-      Float_t quality = fnormdist;
-      if (quality<bestquality && v0.GetDist2()<vertex->GetDist2()){
-       i0=itrack0;
-       i1=itrack1;
-       track0 =htrack0;
-       track1 =htrack1;
-       bestquality = quality;
-       bestradius  = v0.GetRr();
-       vpos[0] = v0.GetXr(0);
-       vpos[1] = v0.GetXr(1);
-       vpos[2] = v0.GetXr(2);
+      pvertex->SetStatus(0);
+      //      if (rejectBase) {
+      //       pvertex->SetStatus(-100);
+      //}
+      if (pvertex->GetPointAngle()>kMinPointAngle2) {
+       if (v0OK){
+         AliV0vertex vertexjuri(*track0,*track1);
+         vertexjuri.SetESDindexes(track0->fESDtrack->GetID(),track1->fESDtrack->GetID());
+         pvertex->SetESDindexes(track0->fESDtrack->GetID(),track1->fESDtrack->GetID());
+         event->AddV0(&vertexjuri);
+         pvertex->SetStatus(100);
+       }
+       event->AddV0MI(pvertex);
       }
     }
   }
   //
   //
-  if (!track0||!track1) return kFALSE;
-  if (track0->fNUsed+track1->fNUsed>0) return kFALSE;  // sharing of clusters not allowed
-  //
-  //propagate to vertex
-  
-  Double_t alpha  = TMath::ATan2(vpos[1],vpos[0]);
-  //
-  AliITStrackMI track0p = *track0;
-  AliITStrackMI track1p = *track1;
-  //
-  //  RefitAt(bestradius+0.5,&track0p,track0);
-  //RefitAt(bestradius+0.5,&track1p,track1);
-  
-  if ( track0p.Propagate(alpha,bestradius+0.2)) {
-    track0= &track0p;
-  }
-  if (track1p.Propagate(alpha,bestradius+0.2)){
-    track1 = &track1p;
-  }
-  //
-  v0.SetM(*track0);
-  v0.SetP(*track1);
-  v0.Update(v);
-  if (v0.GetDist2()<vertex->GetDist2() && v0.GetRr()<20){
-    vertex->SetM(*track0);
-    vertex->SetP(*track1);
-    vertex->Update(v);
-    return kTRUE;
-  }
-  return kFALSE;
+  // delete temporary arrays
+  //  
+  delete[] minPointAngle;
+  delete[] maxr;
+  delete[] minr;
+  delete[] norm;
+  delete[] normdist;
+  delete[] normdist1;
+  delete[] normdist0;
+  delete[] dist;
+  delete[] itsmap;
+  delete[] helixes;
+  delete   pvertex;
 }
 
 
 
 
-Double_t   AliITStrackerMI::TestV0(AliHelix *helix1, AliHelix *helix2, AliESDV0MI *vertex, Double_t &rmin)
-{
-  //
-  // test the helixes for the distnce calculate vertex characteristic
-  //
-  rmin =0;
-  Float_t distance1,distance2;
-  AliHelix & dhelix1 = *helix1;
-  Double_t pp[3],xx[3];
-  dhelix1.GetMomentum(0,pp,0);
-  dhelix1.Evaluate(0,xx);      
-  AliHelix &mhelix = *helix2;    
-  //
-  //find intersection linear
-  //
-  Double_t phase[2][2],radius[2];
-  Int_t  points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
-  Double_t delta1=10000,delta2=10000;  
-  
-  if (points>0){
-    rmin = radius[0];
-    dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
-    dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
-    dhelix1.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
-  }
-  if (points==2){    
-    if (radius[1]<rmin) rmin = radius[1];
-    dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
-    dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
-    dhelix1.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
-  }
-  rmin = TMath::Sqrt(rmin);
-  distance1 = TMath::Min(delta1,delta2);
-  vertex->SetDist1(TMath::Sqrt(distance1));
-  
-  //
-  //find intersection parabolic
-  //
-  points = dhelix1.GetRPHIintersections(mhelix, phase, radius);
-  delta1=10000,delta2=10000;  
-  
-  if (points>0){
-    dhelix1.ParabolicDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
-  }
-  if (points==2){    
-    dhelix1.ParabolicDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
-  }
-  
-  distance2 = TMath::Min(delta1,delta2);
-  vertex->SetDist2(TMath::Sqrt(distance2));      
-  if (distance2<100){
-    if (delta1<delta2){
-      //get V0 info
-      dhelix1.Evaluate(phase[0][0],vertex->GetXrp());
-      dhelix1.GetMomentum(phase[0][0],vertex->GetPPp());
-      mhelix.GetMomentum(phase[0][1],vertex->GetPMp());
-      dhelix1.GetAngle(phase[0][0],mhelix,phase[0][1],vertex->GetAnglep());
-      vertex->SetRr(TMath::Sqrt(radius[0]));
-    }
-    else{
-      dhelix1.Evaluate(phase[1][0],vertex->GetXrp());
-      dhelix1.GetMomentum(phase[1][0],vertex->GetPPp());
-      mhelix.GetMomentum(phase[1][1],vertex->GetPMp());
-      dhelix1.GetAngle(phase[1][0],mhelix,phase[1][1],vertex->GetAnglep());
-      vertex->SetRr(TMath::Sqrt(radius[1]));
-    }
-  }
-  //            
-  //
-  return  vertex->GetDist2();
-}
+
 
 
 
index 79369a33faff27d4b0897082f94d4a635907cdd1..b6f3d9ad79034cdd53883419a59cbf3ef82f208a 100644 (file)
@@ -23,6 +23,8 @@ class AliHelix;
 class AliV0vertex;
 class AliESDV0MI;
 
+class TTreeSRedirector;
+
 
 
 
@@ -110,7 +112,7 @@ public:
     Int_t GetSkip() const {return fSkip;}
     void  SetSkip(Int_t skip){fSkip=skip;}
     void IncAccepted(){fAccepted++;}
-    Int_t GetAccepted() const {return fAccepted;}
+    Int_t GetAccepted() const {return fAccepted;}    
   protected:
     AliITSlayer(const AliITSlayer& /*layer*/){;}
     Double_t fR;                // mean radius of this layer
@@ -172,13 +174,13 @@ public:
   AliITStrackerMI::AliITSdetector & GetDetector(Int_t layer, Int_t n) const {return GetLayer(layer).GetDetector(n); }
 
 protected:
-  void FindV0(AliESD *event);  //try to find V0
-  Double_t  TestV0(AliHelix *h1, AliHelix *h2, AliESDV0MI *vertex, Double_t &rmin);  //try to find V0 - return DCA
-  Bool_t  FindBestPair(Int_t esdtrack0, Int_t esdtrack1,AliESDV0MI *vertex, Int_t &ibest0, Int_t &ibest1);  // try to find best pair from the tree of track hyp.
+  Int_t GetNearestLayer(const Double_t *xr) const;  //get nearest upper layer close to the point xr
+  void FindV02(AliESD *event);  //try to find V0
+  void UpdateTPCV0(AliESD *event);  //try to update, or reject TPC  V0s
   void CookLabel(AliKalmanTrack *t,Float_t wrong) const;
   void CookLabel(AliITStrackMI *t,Float_t wrong) const;
   Double_t GetEffectiveThickness(Double_t y, Double_t z) const;
-  void FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex);
+  void FollowProlongationTree(AliITStrackMI * otrack, Int_t esdindex, Bool_t constrain);
   void ResetBestTrack() {
      fBestTrack.~AliITStrackMI();
      new(&fBestTrack) AliITStrackMI(fTrackToFollow);
@@ -218,6 +220,7 @@ protected:
   AliITStrackMI fBestTrack;              // "best" track 
   AliITStrackMI fTrackToFollow;          // followed track
   TObjArray     fTrackHypothesys;        // ! array with track hypothesys- ARRAY is the owner of tracks- MI
+  TObjArray     fBestHypothesys;         // ! array with track hypothesys- ARRAY is the owner of tracks- MI
   TObjArray     fOriginal;               // ! array with seeds from the TPC
   Int_t         fBestTrackIndex[100000]; // ! index of the best track
   Int_t         fCurrentEsdTrack;        // ! current esd track           - MI
@@ -227,6 +230,8 @@ protected:
   Int_t fLayersNotToSkip[kMaxLayer];     // layer masks
   Int_t fLastLayerToTrackTo;             // the innermost layer to track to
   Float_t * fCoeficients;                //! working array with errors and mean cluser shape
+  AliESD  * fEsd;                        //! pointer to the ESD event
+  TTreeSRedirector *fDebugStreamer;     //!debug streamer
  private:
   AliITStrackerMI(const AliITStrackerMI * /*tracker*/){;}
   ClassDef(AliITStrackerMI,2)   //ITS tracker MI
index efced05e745921724aeb56ee1269835ab9eb720e..be13f576463971c910c0c4de2f189b5980eb9984 100644 (file)
@@ -670,6 +670,10 @@ Int_t AliITStrackerV2::AliITSlayer::InsertCluster(AliITSclusterV2 *c) {
   //--------------------------------------------------------------------
   Float_t circ=TMath::TwoPi()*fR;
   Int_t sec=Int_t(kNsector*c->GetPhiR()/circ);
+  if (sec>=kNsector) {
+     ::Error("InsertCluster","Wrong sector !\n");
+     return 1;
+  }
   Int_t &n=fN[sec];
   if (n>=kMaxClusterPerSector) {
      ::Error("InsertCluster","Too many clusters !\n");
index 74dc37d633d2aea6a59cfe8e957f54d4f268cfa5..ffaa283387827bc9678c42fb12c6968623c0cb68 100644 (file)
@@ -189,6 +189,49 @@ TProfile prof("prof","prof",10,0.5,5);
 
 
 
+TCut cprim("cprim","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<0.01&&abs(MC.fVDist[2])<0.01");
+TCut citsin("citsin","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)<5");
+TCut csec("csec","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)>0.5");
+TCut crec("crec","fReconstructed==1");
+TCut cteta1("cteta1","abs(MC.fParticle.Theta()/3.1415-0.5)<0.25");
+TCut cteta05("cteta05","abs(MC.fParticle.Theta()/3.1415-0.5)<0.1");
+
+TCut cpos1("cpos1","abs(MC.fParticle.fVz/sqrt(MC.fParticle.fVx*MC.fParticle.fVx+MC.fParticle.fVy*MC.fParticle.fVy))<1");
+TCut csens("csens","abs(sqrt(fVDist[0]**2+fVDist[1]**2)-170)<50");
+TCut cmuon("cmuon","abs(MC.fParticle.fPdgCode==-13)");
+TCut cchi2("cchi2","fESDTrack.fITSchi2MIP[0]<7.&&fESDTrack.fITSchi2MIP[1]<5.&&fESDTrack.fITSchi2MIP[2]<7.&&fESDTrack.fITSchi2MIP[3]<7.5&&fESDTrack.fITSchi2MIP[4]<6.");
+  
+
+
+void MakeAliases(AliESDComparisonDraw&comp)
+{
+  //
+  // aliases definition
+  //
+  comp.fTree->SetAlias("radius","TMath::Sqrt(MC.fVDist[0]**2+MC.fVDist[1]**2)");
+  comp.fTree->SetAlias("direction","MC.fParticle.fVx*MC.fParticle.fPx+MC.fParticle.fVy*MC.fParticle.fPy");
+  comp.fTree->SetAlias("decaydir","MC.fTRdecay.fX*MC.fTRdecay.fPx+MC.fTRdecay.fY*MC.fTRdecay.fPy");
+  comp.fTree->SetAlias("theta","MC.fTrackRef.Theta()");
+  comp.fTree->SetAlias("primdca","sqrt(RC.fITStrack.fD[0]**2+RC.fITStrack.fD[1]**2)");
+  comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
+  comp.fTree->SetAlias("trdchi2","fTRDtrack.fChi2/fTRDtrack.fN");
+  
+  comp.fTree->SetAlias("trddedx","(RC.fESDTrack.fTRDsignals[0]+RC.fESDTrack.fTRDsignals[1]+RC.fESDTrack.fTRDsignals[2]+RC.fESDTrack.fTRDsignals[3]+RC.fESDTrack.fTRDsignals[4]+RC.fESDTrack.fTRDsignals[5])/6.");
+  
+  comp.fTree->SetAlias("dtofmc2","fESDTrack.fTrackTime[2]-(10^12*MC.fTOFReferences[0].fTime)");
+  comp.fTree->SetAlias("dtofrc2","(fESDTrack.fTrackTime[2]-fESDTrack.fTOFsignal)");
+
+  comp.fTree->SetAlias("psum","fESDTrack.fTOFr[4]+fESDTrack.fTOFr[3]+fESDTrack.fTOFr[2]+fESDTrack.fTOFr[1]+fESDTrack.fTOFr[0]");
+  comp.fTree->SetAlias("P0","fESDTrack.fTOFr[0]/psum");
+  comp.fTree->SetAlias("P1","fESDTrack.fTOFr[1]/psum");
+  comp.fTree->SetAlias("P2","fESDTrack.fTOFr[2]/psum");
+  comp.fTree->SetAlias("P3","fESDTrack.fTOFr[3]/psum");
+  comp.fTree->SetAlias("P4","fESDTrack.fTOFr[4]/psum");
+  comp.fTree->SetAlias("MaxP","max(max(max(P0,P1),max(P2,P3)),P4)");
+}
+
+
+
 
 void  AliESDRecInfo::UpdatePoints(AliESDtrack*track)
 {
@@ -410,7 +453,7 @@ void AliESDRecInfo::Update(AliMCInfo* info,AliTPCParam * /*par*/, Bool_t reconst
   }
 
 
-  if (fStatus[1]>0 &&info->fNTPCRef>0){
+  if (fStatus[1]>0 &&info->fNTPCRef>0&&TMath::Abs(fTPCinP0[3])>0.0001){
     //TPC
     fESDTrack.GetInnerXYZ(fTPCinR1);
     fTPCinR1[3] = TMath::Sqrt(fTPCinR1[0]*fTPCinR1[0]+fTPCinR1[1]*fTPCinR1[1]);
@@ -930,6 +973,7 @@ Int_t ESDCmpTr::Exec()
    
   fNextTreeGenEntryToRead = 0;
   fNextKinkToRead = 0;
+  fNextV0ToRead   =0;
   cerr<<"fFirstEventNr, fNEvents: "<<fFirstEventNr<<" "<<fNEvents<<endl;
   for (Int_t eventNr = fFirstEventNr; eventNr < fFirstEventNr+fNEvents;
        eventNr++) {
@@ -1298,6 +1342,9 @@ Int_t ESDCmpTr::TreeGenLoop(Int_t eventNr)
       fRecInfo->fESDTrack =*track; 
       if (track->GetITStrack())
        fRecInfo->fITStrack = *((AliITStrackMI*)track->GetITStrack());
+      else{
+       fRecInfo->fITStrack = *track;
+      }
       if (track->GetTRDtrack()){
        fRecInfo->fTRDtrack = *((AliTRDtrack*)track->GetTRDtrack());
       }
index 0fad89192bd7c6869f8d42614888ef47e28a12d6..13d839543cb61225bf1fa1ba2efde837a6a32992 100644 (file)
@@ -51,23 +51,50 @@ AliESDV0MI::AliESDV0MI() :
   //
   //Dafault constructor
   //
+  for (Int_t i=0;i<4;i++){fCausality[i]=0;}
+}
+
+
+void AliESDV0MI::SetCausality(Float_t pb0, Float_t pb1, Float_t pa0, Float_t pa1)
+{
+  //
+  // set probabilities
+  //
+  fCausality[0] = pb0;     // probability - track 0 exist before vertex
+  fCausality[1] = pb1;     // probability - track 1 exist before vertex
+  fCausality[2] = pa0;     // probability - track 0 exist close after vertex
+  fCausality[3] = pa1;     // probability - track 1 exist close after vertex
 }
 
 void AliESDV0MI::SetP(const AliExternalTrackParam & paramp)  {
   //
-  // set mother
+  // set track +
   //
   fParamP   = paramp;
 }
 
 void AliESDV0MI::SetM(const AliExternalTrackParam & paramm){
   //
-  //set daughter
+  //set track -
   //
   fParamM = paramm;
-
 }
   
+void AliESDV0MI::SetRp(const Double_t *rp){
+  //
+  // set pid +
+  //
+  for (Int_t i=0;i<5;i++) fRP[i]=rp[i];
+}
+
+void AliESDV0MI::SetRm(const Double_t *rm){
+  //
+  // set pid -
+  //
+  for (Int_t i=0;i<5;i++) fRM[i]=rm[i];
+}
+
+
 void  AliESDV0MI::UpdatePID(Double_t pidp[5], Double_t pidm[5])
 {
   //
@@ -136,7 +163,8 @@ void  AliESDV0MI::Update(Float_t vertex[3])
   //
   // updates Kink Info
   //
-  Float_t distance1,distance2;
+  //  Float_t distance1,distance2;
+  Float_t distance2;
   //
   AliHelix phelix(fParamP);
   AliHelix mhelix(fParamM);    
@@ -146,7 +174,7 @@ void  AliESDV0MI::Update(Float_t vertex[3])
   Double_t phase[2][2],radius[2];
   Int_t  points = phelix.GetRPHIintersections(mhelix, phase, radius,200);
   Double_t delta1=10000,delta2=10000;  
-
+  /*
   if (points<=0) return;
   if (points>0){
     phelix.LinearDCA(mhelix,phase[0][0],phase[0][1],radius[0],delta1);
@@ -159,6 +187,7 @@ void  AliESDV0MI::Update(Float_t vertex[3])
     phelix.LinearDCA(mhelix,phase[1][0],phase[1][1],radius[1],delta2);
   }
   distance1 = TMath::Min(delta1,delta2);
+  */
   //
   //find intersection parabolic
   //
@@ -215,13 +244,13 @@ void  AliESDV0MI::Update(Float_t vertex[3])
   fDist2 = TMath::Sqrt(distance2);      
   //            
   //
-  Float_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]};
-  Float_t p[3] = {fPP[0]+fPM[0], fPP[1]+fPM[1],fPP[2]+fPM[2]};
-  Float_t vnorm2 = v[0]*v[0]+v[1]*v[1];
-  Float_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
+  Double_t v[3] = {fXr[0]-vertex[0],fXr[1]-vertex[1],fXr[2]-vertex[2]};
+  Double_t p[3] = {fPP[0]+fPM[0], fPP[1]+fPM[1],fPP[2]+fPM[2]};
+  Double_t vnorm2 = v[0]*v[0]+v[1]*v[1];
+  Double_t vnorm3 = TMath::Sqrt(v[2]*v[2]+vnorm2);
   vnorm2 = TMath::Sqrt(vnorm2);
-  Float_t pnorm2 = p[0]*p[0]+p[1]*p[1];
-  Float_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
+  Double_t pnorm2 = p[0]*p[0]+p[1]*p[1];
+  Double_t pnorm3 = TMath::Sqrt(p[2]*p[2]+pnorm2);
   pnorm2 = TMath::Sqrt(pnorm2);  
   fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
   fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);  
index a0c9878c4583b050540d7ced487e7169c7208d4a..4828d08af182a96055d39cf44c27a53cf0843581 100644 (file)
@@ -22,9 +22,15 @@ public:
   //  friend class AliITStrackerMI;
   AliESDV0MI();             //constructor
   //
+  const AliExternalTrackParam *GetParamP() const {return &fParamP;}
+  const AliExternalTrackParam *GetParamM() const {return &fParamM;}
   void SetP(const AliExternalTrackParam & paramp); 
   void SetM(const AliExternalTrackParam & paramd);
+  void SetRp(const Double_t *rp);
+  void SetRm(const Double_t *rm);
   void UpdatePID(Double_t pidp[5], Double_t pidm[5]);
+  void SetStatus(Int_t status){fStatus=status;}
+  Int_t  GetStatus() const {return fStatus;}
   Float_t GetEffMass(UInt_t p1, UInt_t p2);
   Float_t GetProb(UInt_t p1, UInt_t p2);
   void Update(Float_t vertex[3]);            //update
@@ -58,7 +64,8 @@ public:
   Float_t GetNBefore() const {return fNBefore;}
   void SetNBefore(Float_t nb) {fNBefore=nb;}  
   void SetLab(Int_t i, Int_t lab) {fLab[i]=lab;}
-
+  void SetCausality(Float_t pb0, Float_t pb1, Float_t pa0, Float_t pa1);
+  const Float_t * GetCausalityP() const {return fCausality;}
 private:
   AliExternalTrackParam fParamP;
   AliExternalTrackParam fParamM;
@@ -78,12 +85,13 @@ private:
   Double_t       fXr[3];      //rec. position according helix
   Double_t       fAngle[3];   //three angles
   Double_t       fRr;         //rec position of the vertex 
-  Int_t          fStatus;       //status 
+  Int_t          fStatus;       //status  - 1 - TPC V0  2- ITS V0  4- accepted - 0 -rejected
   Int_t          fRow0;         // critical layer
   Int_t          fOrder[3]; //order of the vertex 
   //  quality information
   Double_t       fDistNorm; //normalized  DCA
   Double_t       fDistSigma; //sigma of distance
+  Float_t        fCausality[4];  // causality information - see comments in SetCausality
   Float_t        fChi2Before;   //chi2 of the tracks before V0
   Float_t        fNBefore;      // number of possible points before V0
   Float_t        fChi2After;   // chi2 of the tracks after V0
@@ -92,9 +100,7 @@ private:
   Float_t        fPointAngleTh; //point angle theta
   Float_t        fPointAngle;   //point angle full
 
-
-
-  ClassDef(AliESDV0MI,1)      // ESD V0 vertex
+  ClassDef(AliESDV0MI,2)      // ESD V0 vertex
 };
 
 
index 86262cbf460b11b07456d530cc4f05c57e8a6a4d..de33b519e3ab9443192e8848ca112699e6c03693 100644 (file)
@@ -833,10 +833,42 @@ Int_t AliESDtrack::GetTPCclusters(Int_t *idx) const {
   return fTPCncls;
 }
 
+Float_t AliESDtrack::GetTPCdensity(Int_t row0, Int_t row1) const{
+  //
+  // GetDensity of the clusters on given region between row0 and row1
+  // Dead zone effect takin into acoount
+  //
+  Int_t good  = 0;
+  Int_t found = 0;
+  //  
+  for (Int_t i=row0;i<=row1;i++){     
+    Int_t index = fTPCindex[i];
+    if (index!=-1)  good++;             // track outside of dead zone
+    if (index>0)    found++;
+  }
+  Float_t density=0.5;
+  if (good>(row1-row0)*0.5) density = Float_t(found)/Float_t(good);
+  return density;
+}
 //_______________________________________________________________________
 void AliESDtrack::SetTPCpid(const Double_t *p) {  
   // Sets values for the probability of each particle type (in TPC)
-  for (Int_t i=0; i<AliPID::kSPECIES; i++) fTPCr[i]=p[i];
+  // normalize probabiluty to 1
+  // 
+  //
+//   Double_t sump=0;
+//   for (Int_t i=0; i<AliPID::kSPECIES; i++) {
+//     sump+=p[i];
+//   }
+//   for (Int_t i=0; i<AliPID::kSPECIES; i++) {
+//     if (sump>0){
+//       fTPCr[i]=p[i]/sump;
+//     }
+//     else{
+//       fTPCr[i]=p[i];
+//     }
+//   }
+  for (Int_t i=0; i<AliPID::kSPECIES; i++) fTPCr[i] = p[i];
   SetStatus(AliESDtrack::kTPCpid);
 }
 
index 596f2c424cf6e016ccd1e195d809cb4f1c6236aa..09c69814293fbae89fb0a019c508595fafdb773e 100644 (file)
@@ -106,6 +106,7 @@ public:
   Float_t GetTPCsignal() const {return fTPCsignal;}
   Float_t GetTPCchi2() const {return fTPCchi2;}
   Int_t GetTPCclusters(Int_t *idx) const;
+  Float_t GetTPCdensity(Int_t row0, Int_t row1) const;
   Int_t GetTPCLabel() const {return fTPCLabel;}
   Int_t GetKinkIndex(Int_t i) const { return fKinkIndexes[i];}
   Int_t GetV0Index(Int_t i) const { return fV0Indexes[i];}
@@ -207,7 +208,7 @@ public:
   }; 
 protected:
   
-  AliESDtrack & operator=(const AliESDtrack & );
+  //AliESDtrack & operator=(const AliESDtrack & );
 
   ULong_t   fFlags;        // Reconstruction status flags 
   Int_t     fLabel;        // Track label
@@ -259,7 +260,7 @@ protected:
   // TPC related track information
   Float_t fTPCchi2;        // chi2 in the TPC
   Int_t   fTPCncls;        // number of clusters assigned in the TPC
-  Int_t  fTPCindex[180];  //! indices of the assigned TPC clusters
+  Int_t  fTPCindex[180];  // indices of the assigned TPC clusters
   TBits   fTPCClusterMap;  // Map of clusters, one bit per padrow; 1 if has a cluster on given padrow
   Float_t fTPCsignal;      // detector's PID signal
   Float_t fTPCr[AliPID::kSPECIES]; // "detector response probabilities" (for the PID)
@@ -287,7 +288,7 @@ protected:
   Float_t fTOFsignal;      // detector's PID signal
   Float_t fTOFr[AliPID::kSPECIES]; // "detector response probabilities" (for the PID)
   Int_t   fTOFLabel[3];       // TOF label 
-  Float_t fTOFInfo[10];       //! TOF informations
+  Float_t fTOFInfo[10];       // TOF informations
 
   // PHOS related track information 
   Float_t fPHOSpos[3]; // position localised by PHOS in global coordinate system
@@ -310,7 +311,7 @@ protected:
   Float_t fRICHdx;         // x of the track impact minus x of the MIP
   Float_t fRICHdy;         // y of the track impact minus y of the MIP
        
-  ClassDef(AliESDtrack,13)  //ESDtrack 
+  ClassDef(AliESDtrack,14)  //ESDtrack 
 };
 
 #endif 
index 20ca2fab15c1f209ae83a90cdfe8d7aea254882a..d9381c9295e72c4a12561bf2eca2e2bef5d4d55c 100644 (file)
@@ -33,9 +33,9 @@ public:
   Double_t GetD(Double_t x0=0.,Double_t y0=0.,Double_t z0=0.) const;
   Int_t GetNindex() const {return fNidx;}
   Int_t GetPindex() const {return fPidx;}
+  void SetESDindexes(Int_t ip, Int_t im){fNidx=ip;fPidx=im;}
   void SetDcaDaughters(Double_t rDcaDaughters=0.);
   Double_t GetDcaDaughters() {return fDcaDaughters;}
-
 protected: 
   Int_t fPdgCode;           // reconstructed V0's type (PDG code)
   Double_t fEffMass;        // reconstructed V0's effective mass
index b560cd5599481d4af7742239c8c72a8eec6707d5..e6e40ffaa346072670c874bdd9403a5d93bda0fa 100644 (file)
@@ -202,6 +202,7 @@ AliMCInfo::AliMCInfo()
   fTRDReferences  = new TClonesArray("AliTrackReference",10);
   fTOFReferences  = new TClonesArray("AliTrackReference",10);
   fTRdecay.SetTrack(-1);
+  fCharge = 0;
 }
 
 AliMCInfo::~AliMCInfo()
@@ -440,23 +441,23 @@ void AliGenV0Info::Update(Float_t vertex[3])
   fPointAngleFi = (v[0]*p[0]+v[1]*p[1])/(vnorm2*pnorm2);
   fPointAngleTh = (v[2]*p[2]+vnorm2*pnorm2)/(vnorm3*pnorm3);  
   fPointAngle   = (v[0]*p[0]+v[1]*p[1]+v[2]*p[2])/(vnorm3*pnorm3);
-  Float_t mass1 = fMCd.fMass;
-  Float_t mass2 = fMCm.fMass;
+  Double_t mass1 = fMCd.fMass;
+  Double_t mass2 = fMCm.fMass;
 
   
-  Float_t e1    = TMath::Sqrt(mass1*mass1+
-                             fMCPdr[0]*fMCPdr[0]+
-                             fMCPdr[1]*fMCPdr[1]+
-                             fMCPdr[2]*fMCPdr[2]);
-  Float_t e2    = TMath::Sqrt(mass2*mass2+
+  Double_t e1    = TMath::Sqrt(mass1*mass1+
+                             fMCPd[0]*fMCPd[0]+
+                             fMCPd[1]*fMCPd[1]+
+                             fMCPd[2]*fMCPd[2]);
+  Double_t e2    = TMath::Sqrt(mass2*mass2+
                              fMCPm[0]*fMCPm[0]+
                              fMCPm[1]*fMCPm[1]+
                              fMCPm[2]*fMCPm[2]);
   
   fInvMass =  
-    (fMCPm[0]+fMCPdr[0])*(fMCPm[0]+fMCPdr[0])+
-    (fMCPm[1]+fMCPdr[1])*(fMCPm[1]+fMCPdr[1])+
-    (fMCPm[2]+fMCPdr[2])*(fMCPm[2]+fMCPdr[2]);
+    (fMCPm[0]+fMCPd[0])*(fMCPm[0]+fMCPd[0])+
+    (fMCPm[1]+fMCPd[1])*(fMCPm[1]+fMCPd[1])+
+    (fMCPm[2]+fMCPd[2])*(fMCPm[2]+fMCPd[2]);
   
   //  fInvMass = TMath::Sqrt((e1+e2)*(e1+e2)-fInvMass);
   fInvMass = (e1+e2)*(e1+e2)-fInvMass;
@@ -467,8 +468,6 @@ void AliGenV0Info::Update(Float_t vertex[3])
 
 
 
-
-
 /////////////////////////////////////////////////////////////////////////////////
 void AliGenKinkInfo::Update()
 {
@@ -1130,6 +1129,7 @@ Int_t  AliGenInfoMaker::BuildV0Info()
     if (dinfo){
       v0info->fMCm = (*info);
       v0info->fMCd = (*dinfo);
+      v0info->fMotherP = (*motherparticle);
       v0info->Update(fVPrim);
       br->SetAddress(&v0info);    
       fTreeV0->Fill();
index 1341b6057f81b669bc6f8bdafebdf3029199b698..bc4acb1e7848e5fe92bf5e7f8394507f16f4f16b 100644 (file)
@@ -80,6 +80,7 @@ class AliGenV0Info: public TObject {
 public:
   AliMCInfo fMCd;      //info about daughter particle - second particle for V0
   AliMCInfo fMCm;      //info about mother particle   - first particle for V0
+  TParticle fMotherP;   //particle info about mother particle
   void Update(Float_t vertex[3]);        // put some derived info to special field 
   Double_t    fMCDist1;    //info about closest distance according closest MC - linear DCA
   Double_t    fMCDist2;    //info about closest distance parabolic DCA
index 670d0c74666b0752ba8d3030518f04e2b269c5cc..1810a0f68f3b412210a284bc4d7e5df39fee44c5 100644 (file)
@@ -156,7 +156,7 @@ AliHelix::AliHelix(Double_t x[3], Double_t p[3], Double_t charge, Double_t conve
 
 }
 
-void  AliHelix::GetMomentum(Double_t phase, Double_t p[4],Double_t conversion)
+void  AliHelix::GetMomentum(Double_t phase, Double_t p[4],Double_t conversion, Double_t *xr)
 {
   // return  momentum at given phase
   Double_t x[3],g[3],gg[3];
@@ -166,6 +166,9 @@ void  AliHelix::GetMomentum(Double_t phase, Double_t p[4],Double_t conversion)
   p[0] = fHelix[8]*g[0]/(mt*conversion);
   p[1] = fHelix[8]*g[1]/(mt*conversion);
   p[2] = fHelix[8]*g[2]/(mt*conversion);
+  if (xr){
+    xr[0] = x[0]; xr[1] = x[1]; xr[2] = x[2];
+  }
 }
 
 void   AliHelix::GetAngle(Double_t t1, AliHelix &h, Double_t t2, Double_t angle[3])
@@ -267,6 +270,26 @@ Int_t     AliHelix::GetClosestPhases(AliHelix &h, Double_t phase[2][2])
   return 1;
 }
 
+Double_t  AliHelix::GetPointAngle(AliHelix &h, Double_t phase[2], const Float_t * vertex)
+{
+  //
+  // get point angle bettwen two helixes
+  // 
+  Double_t r0[3],p0[4];
+  Double_t r1[3],p1[4];
+  GetMomentum(phase[0],p0,1,r0);
+  h.GetMomentum(phase[1],p1,1,r1);
+  //
+  Double_t r[3] = {(r0[0]+r1[0])*0.5-vertex[0],(r0[1]+r1[1])*0.5-vertex[1],(r0[2]+r1[2])*0.5-vertex[2]};
+  //intersection point - relative to the prim vertex
+  Double_t p[3] = { p0[0]+p1[0], p0[1]+p1[1],p0[2]+p1[2]};
+  // derivation vector
+  Double_t normr = TMath::Sqrt(r[0]*r[0]+r[1]*r[1]+r[2]*r[2]);
+  Double_t normp = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);   
+  Double_t pointAngle = (r[0]*p[0]+r[1]*p[1]+r[2]*p[2])/(normr*normp);
+  return pointAngle;
+}
+
 Double_t  AliHelix::GetPhase(Double_t x, Double_t y )
                        
 {
index 804caaeb49001c5512d300ded02fa0e15d937ba2..2a21b04c13a4a5e314359495700be6fd953e32dd 100644 (file)
@@ -31,7 +31,7 @@ public:
   void Evaluate(Double_t t, Double_t r[3],  //radius vector
                      Double_t g[3],  //first defivatives
                Double_t gg[3]);     //second derivatives
-  void GetMomentum(Double_t phase, Double_t p[4], Double_t conversion=0.);  // return  momentum  
+  void GetMomentum(Double_t phase, Double_t p[4], Double_t conversion=0., Double_t *xr=0);  // return  momentum  
   void  GetAngle(Double_t t1, AliHelix &h, Double_t t2, Double_t angle[3]);
   inline Double_t GetHelixR(Double_t phase=0);
   inline Double_t GetHelixZ(Double_t phase=0);
@@ -41,6 +41,7 @@ public:
   Int_t     GetPhase(Double_t r0, Double_t t[2]);               //return phase for the nearest point
   Int_t     GetRPHIintersections(AliHelix &h, Double_t phase[2][2], Double_t ri[2], Double_t cut=3.);
   Int_t     GetClosestPhases(AliHelix &h, Double_t phase[2][2]);
+  Double_t  GetPointAngle(AliHelix &h, Double_t phase[2],const Float_t *vertex);
   Int_t     LinearDCA(AliHelix &h, Double_t &t1, Double_t &t2, 
                      Double_t &R, Double_t &dist);
   //
index 1c78a65184dc4e136d9d02b5134d7d8dcb78d4ff..c4b2d6dd375fb1b94fce06129984c5d3a7d41d43 100644 (file)
@@ -207,7 +207,7 @@ void AliMagFMaps::Field(Float_t *x, Float_t *b) const
       //
       //     Constant L3 field , if this option was selected
       //
-         b[2] = - fSolenoid;
+       b[2] = (- fSolenoid)*fFactor;
          return;
       } 
   } else if (fFieldMap[1]->Inside(xm[0], xm[1], xm[2])) {
index a6cdeba07ec502c1958e087a549931fc49000b62..f25118102651a084958852676de75c688c9b46cd 100644 (file)
@@ -232,7 +232,7 @@ TTreeDataElement:: TTreeDataElement(Char_t type) :
 
 TTreeDataElement:: TTreeDataElement(TDataType* type) :
   fName(),
-  fType(' '),
+  fType(0),
   fDType(type),
   fClass(0),
   fPointer(0)
@@ -244,7 +244,7 @@ TTreeDataElement:: TTreeDataElement(TDataType* type) :
 
 TTreeDataElement:: TTreeDataElement(TClass* cl) :
   fName(),
-  fType(' '),
+  fType(0),
   fDType(0),
   fClass(cl),
   fPointer(0)
diff --git a/TPC/AliTPCseed.cxx b/TPC/AliTPCseed.cxx
new file mode 100644 (file)
index 0000000..3d80259
--- /dev/null
@@ -0,0 +1,922 @@
+/**************************************************************************
+ * 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 seed class
+//        This class is used by the AliTPCtrackerMI class
+//      Origin: Marian Ivanov, CERN, Marian.Ivanov@cern.ch
+//-----------------------------------------------------------------
+#include "TClonesArray.h"
+#include "AliTPCseed.h"
+
+ClassImp(AliTPCseed)
+
+
+
+AliTPCseed::AliTPCseed():AliTPCtrack(){
+  //
+  fRow=0; 
+  fRemoval =0; 
+  for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
+  for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
+  for (Int_t i=0;i<3;i++)   fKinkIndexes[i]=0;
+  for (Int_t i=0;i<5;i++)   fTPCr[i]=0.2;
+  fPoints = 0;
+  fEPoints = 0;
+  fNFoundable =0;
+  fNShared  =0;
+  fRemoval = 0;
+  fSort =0;
+  fFirstPoint =0;
+  fNoCluster =0;
+  fBSigned = kFALSE;
+  fSeed1 =-1;
+  fSeed2 =-1;
+  fCurrentCluster =0;
+  fCurrentSigmaY2=0;
+  fCurrentSigmaZ2=0;
+  fCircular = 0;  // not curling track
+}
+AliTPCseed::AliTPCseed(const AliTPCseed &s):AliTPCtrack(s){
+  //---------------------
+  // dummy copy constructor
+  //-------------------------
+  for (Int_t i=0;i<160;i++) fClusterPointer[i] = s.fClusterPointer[i];
+  for (Int_t i=0;i<160;i++) fIndex[i] = s.fIndex[i];
+
+  fPoints  = 0;
+  fEPoints = 0;
+  fCircular =0;
+}
+AliTPCseed::AliTPCseed(const AliTPCtrack &t):AliTPCtrack(t){
+  //
+  //copy constructor
+  fPoints = 0;
+  fEPoints = 0;
+  fNShared  =0; 
+  //  fTrackPoints =0;
+  fRemoval =0;
+  fSort =0;
+  for (Int_t i=0;i<3;i++)   fKinkIndexes[i]=t.GetKinkIndex(i);
+
+  for (Int_t i=0;i<160;i++) {
+    fClusterPointer[i] = 0;
+    Int_t index = t.GetClusterIndex(i);
+    if (index>=-1){ 
+      SetClusterIndex2(i,index);
+    }
+    else{
+      SetClusterIndex2(i,-3); 
+    }    
+  }
+  fFirstPoint =0;
+  fNoCluster =0;
+  fBSigned = kFALSE;
+  fSeed1 =-1;
+  fSeed2 =-1;
+  fCurrentCluster =0;
+  fCurrentSigmaY2=0;
+  fCurrentSigmaZ2=0;
+  fCircular =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) {
+   //
+  //
+  //constructor
+  fRow =0;
+  for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
+  for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
+  for (Int_t i=0;i<3;i++)   fKinkIndexes[i]=0;
+  for (Int_t i=0;i<5;i++)   fTPCr[i]=0.2;
+
+  fPoints = 0;
+  fEPoints = 0;
+  fNFoundable =0;
+  fNShared  = 0;
+  //  fTrackPoints =0;
+  fRemoval =0;
+  fSort =0;
+  fFirstPoint =0;
+  //  fHelixIn = new TClonesArray("AliHelix",0);
+  //fHelixOut = new TClonesArray("AliHelix",0);
+  fNoCluster =0;
+  fBSigned = kFALSE;
+  fSeed1 =-1;
+  fSeed2 =-1;
+  fCurrentCluster =0;
+  fCurrentSigmaY2=0;
+  fCurrentSigmaZ2=0;
+}
+
+AliTPCseed::~AliTPCseed(){
+  //
+  // destructor
+  if (fPoints) delete fPoints;
+  fPoints =0;
+  if (fEPoints) delete fEPoints;
+  fEPoints = 0;
+  fNoCluster =0;
+}
+
+AliTPCTrackerPoint * AliTPCseed::GetTrackPoint(Int_t i)
+{
+  //
+  // 
+  return &fTrackPoints[i];
+}
+
+void AliTPCseed::RebuildSeed()
+{
+  //
+  // rebuild seed to be ready for storing
+  AliTPCclusterMI cldummy;
+  cldummy.SetQ(0);
+  AliTPCTrackPoint pdummy;
+  pdummy.GetTPoint().fIsShared = 10;
+  for (Int_t i=0;i<160;i++){
+    AliTPCclusterMI * cl0 = fClusterPointer[i];
+    AliTPCTrackPoint *trpoint = (AliTPCTrackPoint*)fPoints->UncheckedAt(i);     
+    if (cl0){
+      trpoint->GetTPoint() = *(GetTrackPoint(i));
+      trpoint->GetCPoint() = *cl0;
+      trpoint->GetCPoint().SetQ(TMath::Abs(cl0->GetQ()));
+    }
+    else{
+      *trpoint = pdummy;
+      trpoint->GetCPoint()= cldummy;
+    }
+    
+  }
+
+}
+
+
+Double_t AliTPCseed::GetDensityFirst(Int_t n)
+{
+  //
+  //
+  // return cluster for n rows bellow first point
+  Int_t nfoundable = 1;
+  Int_t nfound      = 1;
+  for (Int_t i=fLastPoint-1;i>0&&nfoundable<n; i--){
+    Int_t index = GetClusterIndex2(i);
+    if (index!=-1) nfoundable++;
+    if (index>0) nfound++;
+  }
+  if (nfoundable<n) return 0;
+  return Double_t(nfound)/Double_t(nfoundable);
+
+}
+
+
+void AliTPCseed::GetClusterStatistic(Int_t first, Int_t last, Int_t &found, Int_t &foundable, Int_t &shared, Bool_t plus2)
+{
+  // get cluster stat.  on given region
+  //
+  found       = 0;
+  foundable   = 0;
+  shared      =0;
+  for (Int_t i=first;i<last; i++){
+    Int_t index = GetClusterIndex2(i);
+    if (index!=-1) foundable++;
+    if (fClusterPointer[i]) {
+      found++;
+    }
+    else 
+      continue;
+
+    if (fClusterPointer[i]->IsUsed(10)) {
+      shared++;
+      continue;
+    }
+    if (!plus2) continue; //take also neighborhoud
+    //
+    if ( (i>0) && fClusterPointer[i-1]){
+      if (fClusterPointer[i-1]->IsUsed(10)) {
+       shared++;
+       continue;
+      }
+    }
+    if ( fClusterPointer[i+1]){
+      if (fClusterPointer[i+1]->IsUsed(10)) {
+       shared++;
+       continue;
+      }
+    }
+    
+  }
+  //if (shared>found){
+    //Error("AliTPCseed::GetClusterStatistic","problem\n");
+  //}
+}
+
+
+
+
+
+void AliTPCseed::Reset(Bool_t all)
+{
+  //
+  //
+  SetNumberOfClusters(0);
+  fNFoundable = 0;
+  SetChi2(0);
+  ResetCovariance();
+  /*
+  if (fTrackPoints){
+    for (Int_t i=0;i<8;i++){
+      delete [] fTrackPoints[i];
+    }
+    delete fTrackPoints;
+    fTrackPoints =0;
+  }
+  */
+
+  if (all){   
+    for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
+    for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
+  }
+
+}
+
+
+void AliTPCseed::Modify(Double_t factor)
+{
+
+  //------------------------------------------------------------------
+  //This function makes a track forget its history :)  
+  //------------------------------------------------------------------
+  if (factor<=0) {
+    ResetCovariance();
+    return;
+  }
+  fC00*=factor;
+  fC10*=0;  fC11*=factor;
+  fC20*=0;  fC21*=0;  fC22*=factor;
+  fC30*=0;  fC31*=0;  fC32*=0;  fC33*=factor;
+  fC40*=0;  fC41*=0;  fC42*=0;  fC43*=0;  fC44*=factor;
+  SetNumberOfClusters(0);
+  fNFoundable =0;
+  SetChi2(0);
+  fRemoval = 0;
+  fCurrentSigmaY2 = 0.000005;
+  fCurrentSigmaZ2 = 0.000005;
+  fNoCluster     = 0;
+  //fFirstPoint = 160;
+  //fLastPoint  = 0;
+}
+
+
+
+
+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;
+
+  if (TMath::Abs(fP4*xk - fP2) >= 0.999) {   
+    return 0;
+  }
+
+  //  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;
+  
+  Double_t dy = dx*(c1+c2)/(r1+r2);
+  Double_t dz = 0;
+  //
+  Double_t delta = fP4*dx*(c1+c2)/(c1*r2 + c2*r1);
+  /*
+  if (TMath::Abs(delta)>0.0001){
+    dz = fP3*TMath::ASin(delta)/fP4;
+  }else{
+    dz = dx*fP3*(c1+c2)/(c1*r2 + c2*r1);
+  }
+  */
+  //  dz =  fP3*AliTPCFastMath::FastAsin(delta)/fP4;
+  dz =  fP3*TMath::ASin(delta)/fP4;
+  //
+  y+=dy;
+  z+=dz;
+  
+
+  return 1;  
+}
+
+
+//_____________________________________________________________________________
+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 (fSort == 0){
+    if (t->fRelativeSector>fRelativeSector) return -1;
+    if (t->fRelativeSector<fRelativeSector) return 1;
+    Double_t z2 = t->GetZ();
+    Double_t z1 = GetZ();
+    if (z2>z1) return 1;
+    if (z2<z1) return -1;
+    return 0;
+  }
+  else {
+    Float_t f2 =1;
+    f2 = 1-20*TMath::Sqrt(t->fC44)/(TMath::Abs(t->GetC())+0.0066);
+    if (t->fBConstrain) f2=1.2;
+
+    Float_t f1 =1;
+    f1 = 1-20*TMath::Sqrt(fC44)/(TMath::Abs(GetC())+0.0066);
+
+    if (fBConstrain)   f1=1.2;
+    if (t->GetNumberOfClusters()*f2 <GetNumberOfClusters()*f1) return -1;
+    else return +1;
+  }
+}
+
+
+
+
+//_____________________________________________________________________________
+Int_t AliTPCseed::Update(const AliTPCclusterMI *c, Double_t chisq, UInt_t /*index*/) {
+  //-----------------------------------------------------------------
+  // This function associates a cluster with this track.
+  //-----------------------------------------------------------------
+  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.9) {
+    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;
+}
+
+
+
+//_____________________________________________________________________________
+void AliTPCseed::CookdEdx(Double_t low, Double_t up,Int_t i1, Int_t i2, Bool_t onlyused) {
+  //-----------------------------------------------------------------
+  // 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+i1;i<i2;i+=4)
+      {
+       Int_t index = fIndex[i];
+       if (index<0||index&0x8000) continue;
+
+       //AliTPCTrackPoint * point = (AliTPCTrackPoint *) arr.At(i);
+       AliTPCTrackerPoint * point = GetTrackPoint(i);
+       //AliTPCTrackerPoint * pointm = GetTrackPoint(i-1);
+       //AliTPCTrackerPoint * pointp = 0;
+       //if (i<159) pointp = GetTrackPoint(i+1);
+
+       if (point==0) continue;
+       AliTPCclusterMI * cl = fClusterPointer[i];
+       if (cl==0) continue;    
+       if (onlyused && (!cl->IsUsed(10))) continue;
+       if (cl->IsUsed(11)) {
+         fNShared++;
+         continue;
+       }
+       Int_t   type   = cl->GetType();
+       //if (point->fIsShared){
+       //  fNShared++;
+       //  continue;
+       //}
+       //if (pointm) 
+       //  if (pointm->fIsShared) continue;
+       //if (pointp) 
+       //  if (pointp->fIsShared) continue;
+
+       if (type<0) continue;
+       //if (type>10) continue;       
+       //if (point->GetErrY()==0) continue;
+       //if (point->GetErrZ()==0) continue;
+
+       //Float_t ddy = (point->GetY()-cl->GetY())/point->GetErrY();
+       //Float_t ddz = (point->GetZ()-cl->GetZ())/point->GetErrZ();
+       //if ((ddy*ddy+ddz*ddz)>10) continue; 
+
+
+       //      if (point->GetCPoint().GetMax()<5) continue;
+       if (cl->GetMax()<5) continue;
+       Float_t angley = point->GetAngleY();
+       Float_t anglez = point->GetAngleZ();
+
+       Float_t rsigmay2 =  point->GetSigmaY();
+       Float_t rsigmaz2 =  point->GetSigmaZ();
+       /*
+       Float_t ns = 1.;
+       if (pointm){
+         rsigmay +=  pointm->GetTPoint().GetSigmaY();
+         rsigmaz +=  pointm->GetTPoint().GetSigmaZ();
+         ns+=1.;
+       }
+       if (pointp){
+         rsigmay +=  pointp->GetTPoint().GetSigmaY();
+         rsigmaz +=  pointp->GetTPoint().GetSigmaZ();
+         ns+=1.;
+       }
+       rsigmay/=ns;
+       rsigmaz/=ns;
+       */
+
+       Float_t rsigma = TMath::Sqrt(rsigmay2*rsigmaz2);
+
+       Float_t ampc   = 0;     // normalization to the number of electrons
+       if (i>64){
+         //      ampc = 1.*point->GetCPoint().GetMax();
+         ampc = 1.*cl->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);
+         //      Float_t zfactor = (AliTPCReconstructor::GetCtgRange()-0.0004*TMath::Abs(point->GetCPoint().GetZ()));
+         //ampc               *=zfactor; 
+       }
+       else{ 
+         //ampc = 1.0*point->GetCPoint().GetMax(); 
+         ampc = 1.0*cl->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);
+         //    Float_t zfactor = (1.02-0.000*TMath::Abs(point->GetCPoint().GetZ()));
+         //ampc               *=zfactor; 
+
+       }
+       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); 
+       //      Float_t z = TMath::Abs(cl->GetZ());
+       if (i<64) {
+         ampc /= 0.6;
+         //ampc /= (1+0.0008*z);
+       } else
+         if (i>128){
+           ampc /=1.5;
+           //ampc /= (1+0.0008*z);
+         }else{
+           //ampc /= (1+0.0008*z);
+         }
+       
+       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 = 50;
+    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;
+
+}
+Double_t AliTPCseed::Bethe(Double_t bg){
+  //
+  // This is the Bethe-Bloch function normalised to 1 at the minimum
+  //
+  Double_t bg2=bg*bg;
+  Double_t bethe;
+  if (bg<3.5e1) 
+    bethe=(1.+ bg2)/bg2*(log(5940*bg2) - bg2/(1.+ bg2));
+  else // Density effect ( approximately :) 
+    bethe=1.15*(1.+ bg2)/bg2*(log(3.5*5940*bg) - bg2/(1.+ bg2));
+  return bethe/11.091;
+}
+
+void AliTPCseed::CookPID()
+{
+  //
+  // cook PID information according dEdx
+  //
+  Double_t fRange = 10.;
+  Double_t fRes   = 0.1;
+  Double_t fMIP   = 47.;
+  //
+  Int_t ns=AliPID::kSPECIES;
+  Double_t sumr =0;
+  for (Int_t j=0; j<ns; j++) {
+    Double_t mass=AliPID::ParticleMass(j);
+    Double_t mom=P();
+    Double_t dedx=fdEdx/fMIP;
+    Double_t bethe=Bethe(mom/mass); 
+    Double_t sigma=fRes*bethe;
+    if (sigma>0.001){
+      if (TMath::Abs(dedx-bethe) > fRange*sigma) {
+       fTPCr[j]=TMath::Exp(-0.5*fRange*fRange)/sigma;
+       sumr+=fTPCr[j];
+       continue;
+      }
+      fTPCr[j]=TMath::Exp(-0.5*(dedx-bethe)*(dedx-bethe)/(sigma*sigma))/sigma;
+      sumr+=fTPCr[j];
+    }
+    else{
+      fTPCr[j]=1.;
+      sumr+=fTPCr[j];
+    }
+  }
+  for (Int_t j=0; j<ns; j++) {
+    fTPCr[j]/=sumr;           //normalize
+  }
+}
+
+/*
+void AliTPCseed::CookdEdx2(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];
+  Bool_t inlimit[200];
+  for (Int_t i=0;i<200;i++) inlimit[i]=kFALSE;
+  for (Int_t i=0;i<200;i++) amp[i]=10000;
+  for (Int_t i=0;i<200;i++) angular[i]= 1;;
+  
+
+  //
+  Float_t meanlog = 100.;
+  Int_t indexde[4]={0,64,128,160};
+
+  Float_t amean     =0;
+  Float_t asigma    =0;
+  Float_t anc       =0;
+  Float_t anorm     =0;
+
+  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<3; of++){    
+  //  for (Int_t i=indexde[of];i<indexde[of+1];i++)
+  for (Int_t i =0; i<160;i++)
+    {
+       AliTPCTrackPoint * point = GetTrackPoint(i);
+       if (point==0) continue;
+       if (point->fIsShared){
+         fNShared++;     
+         continue;
+       }
+       Int_t   type   = point->GetCPoint().GetType();
+       if (type<0) continue;
+       if (point->GetCPoint().GetMax()<5) continue;
+       Float_t angley = point->GetTPoint().GetAngleY();
+       Float_t anglez = point->GetTPoint().GetAngleZ();
+       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 =  point->GetCPoint().GetMax();
+       }
+       else{ 
+         ampc = point->GetCPoint().GetMax(); 
+       }
+       ampc *= 2.0;     // put mean value to channel 50
+       //      ampc *= 0.565;     // put mean value to channel 50
+
+       Float_t w      =  1.;
+       Float_t z = TMath::Abs(point->GetCPoint().GetZ());
+       if (i<64) {
+         ampc /= 0.63;
+       } else
+         if (i>128){
+           ampc /=1.51;
+         }             
+       if (type<0) {  //amp at the border - lower weight                 
+         continue;
+       }
+       if (rsigma>1.5) ampc/=1.3;  // if big backround
+       angular[i]    = TMath::Sqrt(1.+angley*angley+anglez*anglez);
+       amp[i]        = ampc/angular[i];
+       weight[i]     = w;
+       anc++;
+    }
+
+  TMath::Sort(159,amp,index,kFALSE);
+  for (Int_t i=int(anc*low+0.5);i<int(anc*up+0.5);i++){      
+    inlimit[index[i]] = kTRUE;  // take all clusters
+  }
+  
+  //  meanlog = amp[index[Int_t(anc*0.3)]];
+  meanlog =10000.;
+  for (Int_t of =0; of<3; of++){    
+    Float_t sumamp=0;
+    Float_t sumamp2=0;
+    Float_t sumw=0;    
+   for (Int_t i=indexde[of];i<indexde[of+1];i++)
+      {
+       if (inlimit[i]==kFALSE) continue;
+       Float_t ampl      = amp[i];
+       ///angular[i];
+       ampl              = meanlog*TMath::Log(1.+ampl/meanlog);
+       //
+       sumw    += weight[i]; 
+       sumamp  += weight[i]*ampl;
+       sumamp2 += weight[i]*ampl*ampl;
+       norm[of]    += angular[i]*weight[i];
+       nc[of]++;
+      }
+   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;
+   }
+  }
+    
+  Float_t dedx =0;
+  fSdEdx =0;
+  fMAngular =0;
+  //
+  Int_t norm2 = 0;
+  Int_t norm3 = 0;
+  Float_t www[3] = {12.,14.,17.};
+  //Float_t www[3] = {1.,1.,1.};
+
+  for (Int_t i =0;i<3;i++){
+    if (nc[i]>2&&nc[i]<1000){
+      dedx      += mean[i] *nc[i]*www[i]/sigma[i];
+      fSdEdx    += sigma[i]*(nc[i]-2)*www[i]/sigma[i];
+      fMAngular += norm[i] *nc[i];    
+      norm2     += nc[i]*www[i]/sigma[i];
+      norm3     += (nc[i]-2)*www[i]/sigma[i];
+    }
+    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;
+  Float_t norm4 = 0;
+  for (Int_t i =0;i<3;i++){
+    if (nc[i]>2&&nc[i]<1000&&sigma[i]>3){
+      //mean[i]   = mean[i]*(1+0.08*(sigma[i]/(fSdEdx)-1.));
+      dedx      += mean[i] *(nc[i])/(sigma[i]);
+      norm4     += (nc[i])/(sigma[i]);
+    }
+    fDEDX[i]  = mean[i];                
+  }
+  if (norm4>0) dedx /= norm4;
+  
+
+  
+  SetdEdx(dedx);
+    
+  //mi deDX
+
+}
+*/
diff --git a/TPC/AliTPCseed.h b/TPC/AliTPCseed.h
new file mode 100644 (file)
index 0000000..1b01fbb
--- /dev/null
@@ -0,0 +1,120 @@
+#ifndef ALITPCSEED_H
+#define ALITPCSEED_H
+/* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
+ * See cxx source for full Copyright notice                               */
+
+
+/* $Id$ */
+
+//-------------------------------------------------------
+//   TPC seed
+//   Class  needed for TPC parallel tracking 
+//
+//   Origin: 
+//-------------------------------------------------------
+
+#include <TError.h>
+
+#include "AliTPCtrack.h"
+#include "AliComplexCluster.h"
+#include "AliPID.h"
+
+class TFile;
+class AliTPCParam;
+class AliTPCseed;
+class AliTPCclusterMI;
+class AliTPCTrackerPoint;
+class AliESD;   
+class TClonesArray;
+
+class AliTPCseed : public AliTPCtrack {
+  friend class AliTPCtrackerMI;
+  public:  
+     AliTPCseed();
+     virtual ~AliTPCseed();
+     AliTPCseed(const AliTPCtrack &t);
+     AliTPCseed(const AliTPCseed &s);
+     //AliTPCseed(const AliTPCseed &t, Double_t a);
+     AliTPCseed(UInt_t index, const Double_t xx[5], 
+                const Double_t cc[15], Double_t xr, Double_t alpha);     
+     Int_t Compare(const TObject *o) const;
+     void Reset(Bool_t all = kTRUE);
+     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);
+     AliTPCTrackerPoint * GetTrackPoint(Int_t i);
+     void RebuildSeed(); // rebuild seed to be ready for storing
+     Double_t GetDensityFirst(Int_t n);
+     Double_t GetSigma2C() const {return fC44;};
+     void GetClusterStatistic(Int_t first, Int_t last, Int_t &found, Int_t &foundable, Int_t &shared, Bool_t plus2);
+     
+     void Modify(Double_t factor);
+     void SetClusterIndex2(Int_t row, Int_t index) {
+       fIndex[row] = index;
+     }
+     Int_t  GetClusterIndex2(Int_t row) const {
+       return fIndex[row];
+     }
+     Int_t GetClusterSector(Int_t row) const {
+       Int_t pica = -1;
+       if (fIndex[row]>=0) pica =  ((fIndex[row]&0xff000000)>>24);
+       return pica;
+     }
+    
+     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, Int_t i1=0, Int_t i2=159, Bool_t onlyused = kFALSE);
+     void CookPID();
+     Double_t Bethe(Double_t bg);     // return bethe-bloch
+     //     void CookdEdx2(Double_t low=0.05, Double_t up=0.70);
+     Bool_t IsActive() const { return !(fRemoval);}
+     void Desactivate(Int_t reason){ fRemoval = reason;} 
+     //
+     //
+ private:
+     //     AliTPCseed & operator = (const AliTPCseed &)
+     //  {::Fatal("= operator","Not Implemented\n");return *this;}
+     AliESDtrack * fEsd; //!
+     AliTPCclusterMI*   fClusterPointer[160];  //! array of cluster pointers  - 
+     TClonesArray * fPoints;              // array with points along the track
+     TClonesArray * fEPoints;             // array with exact points - calculated in special macro not used in tracking
+     //---CURRENT VALUES
+     Int_t fRow;                 //!current row number  
+     Int_t fSector;              //!current sector number
+     Int_t fRelativeSector;      //! index of current relative sector
+     Float_t fCurrentSigmaY2;    //!expected current cluster sigma Y
+     Float_t fCurrentSigmaZ2;    //!expected current cluster sigma Z
+     Float_t fErrorY2;           //!sigma of current cluster 
+     Float_t fErrorZ2;           //!sigma of current cluster    
+     AliTPCclusterMI * fCurrentCluster; //!pointer to the current cluster for prolongation
+     Int_t   fCurrentClusterIndex1; //! index of the current cluster
+     Bool_t  fInDead;            //! indicate if the track is in dead zone
+     Bool_t  fIsSeeding;         //!indicates if it is proces of seeading
+     Int_t   fNoCluster;         //!indicates number of rows without clusters
+     Int_t   fSort;              //!indicate criteria for sorting
+     Bool_t  fBSigned;        //indicates that clusters of this trackes are signed to be used
+     //
+     //
+     Float_t fDEDX[4];         // dedx according padrows
+     Float_t fSDEDX[4];        // sdedx according padrows
+     Int_t   fNCDEDX[4];       // number of clusters for dedx measurment
+     Double_t fTPCr[AliPID::kSPECIES];   // rough PID according TPC   
+     //
+     Int_t   fSeedType;         //seeding type
+     Int_t   fSeed1;            //first row for seeding
+     Int_t   fSeed2;            //last row for seeding
+     Int_t   fOverlapLabels[12];  //track labels and the length of the  overlap     
+     Float_t fMAngular;           // mean angular factor
+     Char_t   fCircular;           // indicates curlin track
+     AliTPCTrackerPoint  fTrackPoints[160];  //!track points - array track points
+   
+     ClassDef(AliTPCseed,1)  
+};
+
+
+
+
+
+#endif
+
+
index f514d185c27a9f0ca733b7dcce44bbb82224403a..32c780dff81544183eef4fcf43797ed8e0dbc24c 100644 (file)
@@ -38,7 +38,7 @@ AliTPCtrack::AliTPCtrack(): AliKalmanTrack()
   fX = fP0 = fP1 = fP2 = fP3 = fP3 = fP4 = 0.0;
   fAlpha = fdEdx = 0.0;
   fNumber = 0;  // [SR, 01.04.2003]
-  for (Int_t i=0; i<3;i++) fKinkIndexes[i]=0;
+  for (Int_t i=0; i<3;i++) {fKinkIndexes[i]=0; fV0Indexes[i]=0;}
 }
 
 //_________________________________________________________________________
@@ -77,6 +77,7 @@ const Double_t cc[15], Double_t xref, Double_t alpha) : AliKalmanTrack() {
   fTrackType  = 0;
   fLab2       = 0;
   for (Int_t i=0; i<3;i++) fKinkIndexes[i]=0;
+  for (Int_t i=0; i<3;i++) fV0Indexes[i]=0;
 }
 
 //_____________________________________________________________________________
@@ -88,6 +89,7 @@ AliTPCtrack::AliTPCtrack(const AliESDtrack& t) : AliKalmanTrack() {
   SetLabel(t.GetLabel());
   SetMass(t.GetMass());
   for (Int_t i=0; i<3;i++) fKinkIndexes[i]=t.GetKinkIndex(i);
+  for (Int_t i=0; i<3;i++) fV0Indexes[i]=t.GetV0Index(i);
 
   fdEdx  = t.GetTPCsignal();
   fAlpha = t.GetAlpha();
@@ -165,6 +167,7 @@ AliTPCtrack::AliTPCtrack(const AliTPCtrack& t) : AliKalmanTrack(t) {
   fTrackType  = t.fTrackType;
   fLab2       = t.fLab2;
   for (Int_t i=0; i<3;i++) fKinkIndexes[i]=t.fKinkIndexes[i];
+  for (Int_t i=0; i<3;i++) fV0Indexes[i]=t.fV0Indexes[i];
 }
 
 //_____________________________________________________________________________
index c339ce07d00e33707e4f7f1c6e00d6ac23464337..a99cdaddc97e470b3cdb3154129e174a7aca17f6 100644 (file)
@@ -108,8 +108,10 @@ public:
   Double_t GetD(Double_t x=0, Double_t y=0) const;
   AliExternalTrackParam & GetReference(){ return fReference;}
   void UpdateReference(){ new (&fReference) AliExternalTrackParam(*this);}
-  Int_t  GetKinkIndex(Int_t i) const{ return fKinkIndexes[i];}
+  Int_t   GetKinkIndex(Int_t i) const{ return fKinkIndexes[i];}
   Int_t*  GetKinkIndexes() { return fKinkIndexes;}
+  Int_t   GetV0Index(Int_t i) const{ return fV0Indexes[i];}
+  Int_t*  GetV0Indexes() { return fV0Indexes;}
 protected: 
   Double_t fX;              // X-coordinate of this track (reference plane)
   Double_t fAlpha;          // Rotation angle the local (TPC sector)
@@ -147,6 +149,7 @@ protected:
   AliExternalTrackParam   fReference; // track parameters at the middle of the chamber
   Float_t  fKinkPoint[12];      //radius, of kink,  dfi and dtheta
   Int_t    fKinkIndexes[3];     // kink indexes - minus = mother + daughter
+  Int_t    fV0Indexes[3];     // kink indexes - minus = mother + daughter
 
   ClassDef(AliTPCtrack,2)   // Time Projection Chamber reconstructed tracks
 };
index f35ec1290eb5c9cd1f96e24c3d37c3230a2c14f1..5ef0349b5383c9411b60c68cc0c95f1017f5ef50 100644 (file)
 #include "AliTPCReconstructor.h"
 #include "AliTPCclusterMI.h"
 #include "AliTPCpolyTrack.h"
-#include "AliTPCreco.h" 
+#include "AliTPCreco.h"
+#include "AliTPCseed.h" 
 #include "AliTPCtrackerMI.h"
 #include "TStopwatch.h"
+#include "AliTPCReconstructor.h"
+#include "AliESDkink.h"
+#include "AliPID.h"
 #include "TTreeStream.h"
 //
 
-ClassImp(AliTPCseed)
 ClassImp(AliTPCtrackerMI)
 
 
@@ -255,6 +258,7 @@ AliTracker(), fkNIS(par->GetNInnerSector()/2), fkNOS(par->GetNOuterSector()/2)
   fNewIO     =0;
   fDebug     =0;
   fEvent     =0;
+  fDebugStreamer = new TTreeSRedirector("TPCdebug.root");
 }
 //________________________________________________________________________
 AliTPCtrackerMI::AliTPCtrackerMI(const AliTPCtrackerMI &t):
@@ -283,6 +287,7 @@ AliTPCtrackerMI::~AliTPCtrackerMI() {
     fSeeds->Delete(); 
     delete fSeeds;
   }
+  if (fDebugStreamer) delete fDebugStreamer;
 }
 
 void AliTPCtrackerMI::SetIO()
@@ -372,6 +377,8 @@ void AliTPCtrackerMI::FillESD(TObjArray* arr)
        iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);
        iotrack.SetTPCPoints(pt->GetPoints());
        iotrack.SetKinkIndexes(pt->GetKinkIndexes());
+       iotrack.SetV0Indexes(pt->GetV0Indexes());
+       //      iotrack.SetTPCpid(pt->fTPCr);
        //iotrack.SetTPCindex(i);
        fEvent->AddTrack(&iotrack);
        continue;
@@ -383,6 +390,8 @@ void AliTPCtrackerMI::FillESD(TObjArray* arr)
        iotrack.SetTPCPoints(pt->GetPoints());
        //iotrack.SetTPCindex(i);
        iotrack.SetKinkIndexes(pt->GetKinkIndexes());
+       iotrack.SetV0Indexes(pt->GetV0Indexes());
+       //      iotrack.SetTPCpid(pt->fTPCr);
        fEvent->AddTrack(&iotrack);
        continue;
       } 
@@ -398,6 +407,8 @@ void AliTPCtrackerMI::FillESD(TObjArray* arr)
          //iotrack.SetTPCindex(i);
          iotrack.SetTPCPoints(pt->GetPoints());
          iotrack.SetKinkIndexes(pt->GetKinkIndexes());
+         iotrack.SetV0Indexes(pt->GetV0Indexes());
+         //iotrack.SetTPCpid(pt->fTPCr);
          fEvent->AddTrack(&iotrack);
          continue;
        }
@@ -413,6 +424,8 @@ void AliTPCtrackerMI::FillESD(TObjArray* arr)
        iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);      
        iotrack.SetTPCPoints(pt->GetPoints());
        iotrack.SetKinkIndexes(pt->GetKinkIndexes());
+       iotrack.SetV0Indexes(pt->GetV0Indexes());
+       //iotrack.SetTPCpid(pt->fTPCr);
        //iotrack.SetTPCindex(i);
        fEvent->AddTrack(&iotrack);
        continue;
@@ -427,6 +440,8 @@ void AliTPCtrackerMI::FillESD(TObjArray* arr)
          iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);    
          iotrack.SetTPCPoints(pt->GetPoints());
          iotrack.SetKinkIndexes(pt->GetKinkIndexes());
+         iotrack.SetV0Indexes(pt->GetV0Indexes());
+         //iotrack.SetTPCpid(pt->fTPCr);       
          //iotrack.SetTPCindex(i);
          fEvent->AddTrack(&iotrack);
          continue;
@@ -444,6 +459,8 @@ void AliTPCtrackerMI::FillESD(TObjArray* arr)
        iotrack.UpdateTrackParams(pt,AliESDtrack::kTPCin);      
        iotrack.SetTPCPoints(pt->GetPoints());
        iotrack.SetKinkIndexes(pt->GetKinkIndexes());
+       iotrack.SetV0Indexes(pt->GetV0Indexes());
+       //      iotrack.SetTPCpid(pt->fTPCr);
        //iotrack.SetTPCindex(i);
        fEvent->AddTrack(&iotrack);
        continue;
@@ -895,160 +912,6 @@ Double_t AliTPCtrackerMI::ErrZ2(AliTPCseed* seed, AliTPCclusterMI * cl){
 
 
 
-void AliTPCseed::Reset(Bool_t all)
-{
-  //
-  //
-  SetNumberOfClusters(0);
-  fNFoundable = 0;
-  SetChi2(0);
-  ResetCovariance();
-  /*
-  if (fTrackPoints){
-    for (Int_t i=0;i<8;i++){
-      delete [] fTrackPoints[i];
-    }
-    delete fTrackPoints;
-    fTrackPoints =0;
-  }
-  */
-
-  if (all){   
-    for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
-    for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
-  }
-
-}
-
-
-void AliTPCseed::Modify(Double_t factor)
-{
-
-  //------------------------------------------------------------------
-  //This function makes a track forget its history :)  
-  //------------------------------------------------------------------
-  if (factor<=0) {
-    ResetCovariance();
-    return;
-  }
-  fC00*=factor;
-  fC10*=0;  fC11*=factor;
-  fC20*=0;  fC21*=0;  fC22*=factor;
-  fC30*=0;  fC31*=0;  fC32*=0;  fC33*=factor;
-  fC40*=0;  fC41*=0;  fC42*=0;  fC43*=0;  fC44*=factor;
-  SetNumberOfClusters(0);
-  fNFoundable =0;
-  SetChi2(0);
-  fRemoval = 0;
-  fCurrentSigmaY2 = 0.000005;
-  fCurrentSigmaZ2 = 0.000005;
-  fNoCluster     = 0;
-  //fFirstPoint = 160;
-  //fLastPoint  = 0;
-}
-
-
-
-
-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;
-
-  if (TMath::Abs(fP4*xk - fP2) >= 0.999) {   
-    return 0;
-  }
-
-  //  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;
-  
-  Double_t dy = dx*(c1+c2)/(r1+r2);
-  Double_t dz = 0;
-  //
-  Double_t delta = fP4*dx*(c1+c2)/(c1*r2 + c2*r1);
-  /*
-  if (TMath::Abs(delta)>0.0001){
-    dz = fP3*TMath::ASin(delta)/fP4;
-  }else{
-    dz = dx*fP3*(c1+c2)/(c1*r2 + c2*r1);
-  }
-  */
-  dz =  fP3*AliTPCFastMath::FastAsin(delta)/fP4;
-  //
-  y+=dy;
-  z+=dz;
-  
-
-  return 1;  
-}
-
-
-//_____________________________________________________________________________
-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 (fSort == 0){
-    if (t->fRelativeSector>fRelativeSector) return -1;
-    if (t->fRelativeSector<fRelativeSector) return 1;
-    Double_t z2 = t->GetZ();
-    Double_t z1 = GetZ();
-    if (z2>z1) return 1;
-    if (z2<z1) return -1;
-    return 0;
-  }
-  else {
-    Float_t f2 =1;
-    f2 = 1-20*TMath::Sqrt(t->fC44)/(TMath::Abs(t->GetC())+0.0066);
-    if (t->fBConstrain) f2=1.2;
-
-    Float_t f1 =1;
-    f1 = 1-20*TMath::Sqrt(fC44)/(TMath::Abs(GetC())+0.0066);
-
-    if (fBConstrain)   f1=1.2;
-    if (t->GetNumberOfClusters()*f2 <GetNumberOfClusters()*f1) return -1;
-    else return +1;
-  }
-}
 
 void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
 {
@@ -1071,65 +934,6 @@ void AliTPCtrackerMI::RotateToLocal(AliTPCseed *seed)
 
 
 
-
-//_____________________________________________________________________________
-Int_t AliTPCseed::Update(const AliTPCclusterMI *c, Double_t chisq, UInt_t /*index*/) {
-  //-----------------------------------------------------------------
-  // This function associates a cluster with this track.
-  //-----------------------------------------------------------------
-  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.9) {
-    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::F1old(Double_t x1,Double_t y1,
                    Double_t x2,Double_t y2,
@@ -4083,6 +3887,7 @@ void  AliTPCtrackerMI::FindKinks(TObjArray * array, AliESD *esd)
   //
 
   TObjArray *kinks= new TObjArray(10000);
+  //  TObjArray *v0s= new TObjArray(10000);
   Int_t nentries = array->GetEntriesFast();
   AliHelix *helixes      = new AliHelix[nentries];
   Int_t    *sign         = new Int_t[nentries];
@@ -4095,6 +3900,8 @@ void  AliTPCtrackerMI::FindKinks(TObjArray * array, AliESD *esd)
   Float_t  *fim          = new Float_t[nentries];
   Float_t  *shared       = new Float_t[nentries];
   Bool_t   *circular     = new Bool_t[nentries];
+  Float_t *dca          = new Float_t[nentries];
+  //const AliESDVertex * primvertex = esd->GetVertex();
   //
   //  nentries = array->GetEntriesFast();
   //
@@ -4106,6 +3913,7 @@ void  AliTPCtrackerMI::FindKinks(TObjArray * array, AliESD *esd)
     usage[i]=0;
     AliTPCseed* track = (AliTPCseed*)array->At(i);    
     if (!track) continue;
+    track->fCircular =0;
     shared[i] = kFALSE;
     track->UpdatePoints();
     if (( track->GetPoints()[2]- track->GetPoints()[0])>5 && track->GetPoints()[3]>0.8){
@@ -4128,7 +3936,8 @@ void  AliTPCtrackerMI::FindKinks(TObjArray * array, AliESD *esd)
     }   
     z0[i]=1000;
     circular[i]= kFALSE;
-    if (track->GetProlongation(0,y,z))  z0[i] = TMath::Abs(z);
+    if (track->GetProlongation(0,y,z))  z0[i] = z;
+    dca[i] = track->GetD(0,0);    
   }
   //
   //
@@ -4141,7 +3950,7 @@ void  AliTPCtrackerMI::FindKinks(TObjArray * array, AliESD *esd)
 
   //
   // Find circling track
-  //  TTreeSRedirector cstream("circling.root");
+  TTreeSRedirector &cstream = *fDebugStreamer;
   //
   for (Int_t i0=0;i0<nentries;i0++){
     AliTPCseed * track0 = (AliTPCseed*)array->At(i0);
@@ -4151,17 +3960,23 @@ void  AliTPCtrackerMI::FindKinks(TObjArray * array, AliESD *esd)
     for (Int_t i1=i0+1;i1<nentries;i1++){
       AliTPCseed * track1 = (AliTPCseed*)array->At(i1);
       if (!track1) continue;
-      if (TMath::Abs(1./track1->fP4)>200) continue;
-      if (track1->fN<40) continue;
-      if (z0[i0]<20&&z0[i1]<20) continue;
-      if (track1->fP4*track0->fP4>0) continue;
-      if (track1->fP3*track0->fP3>0) continue;
+      if (track1->fN<40)                  continue;
+      if ( TMath::Abs(track1->fP3+track0->fP3)>0.1) continue;
       if (track0->fBConstrain&&track1->fBConstrain) continue;
+      if (TMath::Abs(1./track1->fP4)>200) continue;
+      if (track1->fP4*track0->fP4>0)      continue;
+      if (track1->fP3*track0->fP3>0)      continue;
       if (max(TMath::Abs(1./track0->fP4),TMath::Abs(1./track1->fP4))>190) continue;
       if (track0->fBConstrain&&TMath::Abs(track1->fP4)<TMath::Abs(track0->fP4)) continue; //returning - lower momenta
       if (track1->fBConstrain&&TMath::Abs(track0->fP4)<TMath::Abs(track1->fP4)) continue; //returning - lower momenta
       //
-      if ( TMath::Abs(track1->fP3+track0->fP3)>0.1) continue;
+      Float_t mindcar = TMath::Min(TMath::Abs(dca[i0]),TMath::Abs(dca[i1]));
+      if (mindcar<5)   continue;
+      Float_t mindcaz = TMath::Min(TMath::Abs(z0[i0]-GetZ()),TMath::Abs(z0[i1]-GetZ()));
+      if (mindcaz<5) continue;
+      if (mindcar+mindcaz<20) continue;
+      //
+      //
       Float_t xc0 = helixes[i0].GetHelix(6);
       Float_t yc0 = helixes[i0].GetHelix(7);
       Float_t r0  = helixes[i0].GetHelix(8);
@@ -4171,7 +3986,7 @@ void  AliTPCtrackerMI::FindKinks(TObjArray * array, AliESD *esd)
        
       Float_t rmean = (r0+r1)*0.5;
       Float_t delta =TMath::Sqrt((xc1-xc0)*(xc1-xc0)+(yc1-yc0)*(yc1-yc0));
-      if (delta>30) continue;
+      //if (delta>30) continue;
       if (delta>rmean*0.25) continue;
       if (TMath::Abs(r0-r1)/rmean>0.3) continue; 
       //
@@ -4198,44 +4013,66 @@ void  AliTPCtrackerMI::FindKinks(TObjArray * array, AliESD *esd)
       */
       if (npoints>0){
        Int_t ibest=0;
-       helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],5);
+       helixes[i0].ParabolicDCA(helixes[i1],phase[0][0],phase[0][1],radius[0],deltah[0],2);
        if (npoints==2){
-         helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],5);
+         helixes[i0].ParabolicDCA(helixes[i1],phase[1][0],phase[1][1],radius[1],deltah[1],2);
          if (deltah[1]<deltah[0]) ibest=1;
        }
        deltabest  = TMath::Sqrt(deltah[ibest]);
        helixes[i0].Evaluate(phase[ibest][0],xyz0);
        helixes[i1].Evaluate(phase[ibest][1],xyz1);
        helixes[i0].GetAngle(phase[ibest][0],helixes[i1],phase[ibest][1],hangles);
-       //Double_t radiusbest = TMath::Sqrt(radius[ibest]);
+       Double_t radiusbest = TMath::Sqrt(radius[ibest]);
        //
-       if (deltabest<25){        
-         /*
-         cstream<<"CC1"<<track0->fLab<<track1->fLab<<   //0-1
-           track0->fP3<<track1->fP3<<                   //2-3 
-           track0->fP4<<track1->fP4<<                   //4-5
-           delta<<rmean<<npoints<<                      //6-8
-           hangles[0]<<hangles[2]<<                     //9-10
-           xyz0[2]<<xyz1[2]<<radiusbest<<deltabest<< //11--14
-           track0->fFirstPoint<<track1->fFirstPoint<<   //15-16 
-           track0->fLastPoint<<track1->fLastPoint<<     //17-18
-           track0<<track1<<                             //19-20
-           phase[ibest][0]<<phase[ibest][1]<<"\n"; //21-22       
-         */
-       }
+       if (deltabest>6) continue;
+       if (mindcar+mindcaz<40 && (hangles[2]<3.12||deltabest>3)) continue;
        Bool_t sign =kFALSE;
-       if (deltabest<3&&hangles[2]>3.06) sign =kTRUE;
-       if (TMath::Min(TMath::Abs(track0->GetD()),TMath::Abs(track1->GetD()))>10&&deltabest<6&&hangles[2]>3.) sign= kTRUE;
+       if (hangles[2]>3.06) sign =kTRUE;
+       //
        if (sign){
          circular[i0] = kTRUE;
-         circular[i1] = kTRUE;           
+         circular[i1] = kTRUE;
+         if (TMath::Abs(track0->fP4)<TMath::Abs(track1->fP4)){
+           track0->fCircular += 1;
+           track1->fCircular += 2;
+         }
+         else{
+           track1->fCircular += 1;
+           track0->fCircular += 2;
+         }
+       }               
+       if (sign){        
+         //debug stream
+         cstream<<"Curling"<<
+           "lab0="<<track0->fLab<<
+           "lab1="<<track1->fLab<<   
+           "Tr0.="<<track0<<
+           "Tr1.="<<track1<<      
+           "dca0="<<dca[i0]<<
+           "dca1="<<dca[i1]<<
+           "mindcar="<<mindcar<<
+           "mindcaz="<<mindcaz<<
+           "delta="<<delta<<
+           "rmean="<<rmean<<
+           "npoints="<<npoints<<                      
+           "hangles0="<<hangles[0]<<
+           "hangles2="<<hangles[2]<<                    
+           "xyz0="<<xyz0[2]<<
+           "xyzz1="<<xyz1[2]<<
+           "z0="<<z0[i0]<<
+           "z1="<<z0[i1]<<
+           "radius="<<radiusbest<<
+           "deltabest="<<deltabest<< 
+           "phase0="<<phase[ibest][0]<<
+           "phase1="<<phase[ibest][1]<<
+           "\n";                 
        }
       }
     }
   }
   //
-  //
-  //
+  //  Finf kinks loop
+  // 
   //
   for (Int_t i =0;i<nentries;i++){
     if (sign[i]==0) continue;
@@ -4624,6 +4461,12 @@ void  AliTPCtrackerMI::FindKinks(TObjArray * array, AliESD *esd)
   //
   RemoveUsed2(array,0.5,0.4,30);
   UnsignClusters();
+  for (Int_t i=0;i<nentries;i++){
+    AliTPCseed * track0 = (AliTPCseed*)array->At(i);
+    if (!track0) continue;
+    track0->CookdEdx(0.02,0.6);
+    track0->CookPID();
+  }
   //
   for (Int_t i=0;i<nentries;i++){
     AliTPCseed * track0 = (AliTPCseed*)array->At(i);
@@ -4681,6 +4524,7 @@ void  AliTPCtrackerMI::FindKinks(TObjArray * array, AliESD *esd)
   delete [] mothers;
   //
   //
+  delete [] dca;
   delete []circular;
   delete []shared;
   delete []quality;
@@ -4702,6 +4546,342 @@ void  AliTPCtrackerMI::FindKinks(TObjArray * array, AliESD *esd)
   timer.Print();
 }
 
+void  AliTPCtrackerMI::FindV0s(TObjArray * array, AliESD *esd)
+{
+  //
+  //  find V0s
+  //
+  //
+  TObjArray *tpcv0s      = new TObjArray(100000);
+  Int_t     nentries     = array->GetEntriesFast();
+  AliHelix *helixes      = new AliHelix[nentries];
+  Int_t    *sign         = new Int_t[nentries];
+  Float_t  *alpha        = new Float_t[nentries];
+  Float_t  *z0           = new Float_t[nentries]; 
+  Float_t  *dca          = new Float_t[nentries];
+  Float_t  *sdcar        = new Float_t[nentries];
+  Float_t  *cdcar        = new Float_t[nentries];
+  Float_t  *pulldcar     = new Float_t[nentries];
+  Float_t  *pulldcaz     = new Float_t[nentries];
+  Float_t  *pulldca      = new Float_t[nentries];
+  Bool_t   *isPrim       = new Bool_t[nentries];  
+  const AliESDVertex * primvertex = esd->GetVertex();
+  Double_t             zvertex = primvertex->GetZv(); 
+  //
+  //  nentries = array->GetEntriesFast();
+  //
+  for (Int_t i=0;i<nentries;i++){
+    sign[i]=0;
+    isPrim[i]=0;
+    AliTPCseed* track = (AliTPCseed*)array->At(i);    
+    if (!track) continue;
+    track->GetV0Indexes()[0] = 0;  //rest v0 indexes
+    track->GetV0Indexes()[1] = 0;  //rest v0 indexes
+    track->GetV0Indexes()[2] = 0;  //rest v0 indexes
+    //
+    alpha[i] = track->GetAlpha();
+    new (&helixes[i]) AliHelix(*track);
+    Double_t xyz[3];
+    helixes[i].Evaluate(0,xyz);
+    sign[i] = (track->GetC()>0) ? -1:1;
+    Double_t x,y,z;
+    x=160;
+    z0[i]=1000;
+    if (track->GetProlongation(0,y,z))  z0[i] = z;
+    dca[i] = track->GetD(0,0);
+    // 
+    // dca error parrameterezation + pulls
+    //
+    sdcar[i]      = TMath::Sqrt(0.150*0.150+(100*track->fP4)*(100*track->fP4));
+    if (TMath::Abs(track->fP3)>1) sdcar[i]*=2.5;
+    cdcar[i]      = TMath::Exp((TMath::Abs(track->fP4)-0.0106)*525.3);
+    pulldcar[i]   = (dca[i]-cdcar[i])/sdcar[i];
+    pulldcaz[i]   = (z0[i]-zvertex)/sdcar[i];
+    pulldca[i]    = TMath::Sqrt(pulldcar[i]*pulldcar[i]+pulldcaz[i]*pulldcaz[i]);
+    if (track->fTPCr[1]+track->fTPCr[2]+track->fTPCr[3]>0.5) {
+      if (pulldca[i]<3.) isPrim[i]=kTRUE;  //pion, muon and Kaon  3 sigma cut
+    }
+    if (track->fTPCr[4]>0.5) {
+      if (pulldca[i]<0.5) isPrim[i]=kTRUE;  //proton 0.5 sigma cut
+    }
+    if (track->fTPCr[0]>0.4) {
+      isPrim[i]=kFALSE;  //electron no  sigma cut
+    }
+  }
+  //
+  //
+  TStopwatch timer;
+  timer.Start();
+  Int_t ncandidates =0;
+  Int_t nall =0;
+  Int_t ntracks=0; 
+  Double_t phase[2][2],radius[2];
+  //
+  //  Finf V0s loop
+  // 
+  //
+  // //  
+  TTreeSRedirector &cstream = *fDebugStreamer; 
+  Float_t fprimvertex[3]={GetX(),GetY(),GetZ()};
+  AliESDV0MI vertex; 
+  Double_t cradius0 = 10*10;
+  Double_t cradius1 = 200*200;
+  Double_t cdist1=3.;
+  Double_t cdist2=4.;
+  Double_t cpointAngle = 0.95;
+  //
+  Double_t delta[2]={10000,10000};
+  for (Int_t i =0;i<nentries;i++){
+    if (sign[i]==0) continue;
+    AliTPCseed * track0 = (AliTPCseed*)array->At(i);
+    if (!track0) continue;
+    cstream<<"Tracks"<<
+      "Tr0.="<<track0<<
+      "dca="<<dca[i]<<
+      "z0="<<z0[i]<<
+      "zvertex="<<zvertex<<
+      "sdcar0="<<sdcar[i]<<
+      "cdcar0="<<cdcar[i]<<
+      "pulldcar0="<<pulldcar[i]<<
+      "pulldcaz0="<<pulldcaz[i]<<
+      "pulldca0="<<pulldca[i]<<
+      "isPrim="<<isPrim[i]<<
+      "\n";
+    //
+    if (track0->fP4<0) continue;
+    if (track0->GetKinkIndex(0)>0||isPrim[i]) continue;   //daughter kink
+    //
+    if (TMath::Abs(helixes[i].GetHelix(4))<0.000000001) continue;
+    ntracks++;
+    // debug output
+    
+    
+    for (Int_t j =0;j<nentries;j++){
+      AliTPCseed * track1 = (AliTPCseed*)array->At(j);
+      if (!track1) continue;
+      if (track1->GetKinkIndex(0)>0 || isPrim[j]) continue; //daughter kink
+      if (sign[j]*sign[i]>0) continue; 
+      if (TMath::Abs(helixes[j].GetHelix(4))<0.000001) continue;
+      if (track0->fCircular+track1->fCircular>1) continue;    //circling -returning track
+      nall++;
+      //
+      // DCA to prim vertex cut
+      //
+      //
+      delta[0]=10000;
+      delta[1]=10000;
+
+      Int_t npoints = helixes[i].GetRPHIintersections(helixes[j], phase, radius,cdist2);
+      if (npoints<1) continue;
+      Int_t iclosest=0;
+      // cuts on radius      
+      if (npoints==1){
+       if (radius[0]<cradius0||radius[0]>cradius1) continue;
+       helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);
+       if (delta[0]>cdist1) continue;
+      }
+      else{
+       if (TMath::Max(radius[0],radius[1])<cradius0|| TMath::Min(radius[0],radius[1])>cradius1) continue;      
+       helixes[i].LinearDCA(helixes[j],phase[0][0],phase[0][1],radius[0],delta[0]);    
+       helixes[i].LinearDCA(helixes[j],phase[1][0],phase[1][1],radius[1],delta[1]);
+       if (delta[1]<delta[0]) iclosest=1;
+       if (delta[iclosest]>cdist1) continue;
+      }
+      helixes[i].ParabolicDCA(helixes[j],phase[iclosest][0],phase[iclosest][1],radius[iclosest],delta[iclosest]);
+      if (radius[iclosest]<cradius0 || radius[iclosest]>cradius1 || delta[iclosest]>cdist1) continue;
+      //
+      Double_t pointAngle = helixes[i].GetPointAngle(helixes[j],phase[iclosest],fprimvertex);
+      if (pointAngle<cpointAngle) continue;
+      //
+      Bool_t isGamma = kFALSE;
+      vertex.SetP(*track0); //track0 - plus
+      vertex.SetM(*track1); //track1 - minus
+      vertex.Update(fprimvertex);
+      if (track0->fTPCr[0]>0.3&&track1->fTPCr[0]>0.3&&vertex.GetAnglep()[2]<0.15) isGamma=kTRUE;              // gamma conversion candidate
+      Double_t pointAngle2 = vertex.GetPointAngle();
+      //continue;
+      if (vertex.GetPointAngle()<cpointAngle && (!isGamma)) continue; // point angle cut
+      if (vertex.GetDist2()>2&&(!isGamma)) continue;         // point angle cut
+      Float_t sigmae     = 0.15*0.15;
+      if (vertex.GetRr()<80) 
+       sigmae += (sdcar[i]*sdcar[i]+sdcar[j]*sdcar[j])*(1.-vertex.GetRr()/80.)*(1.-vertex.GetRr()/80.);
+      sigmae+= TMath::Sqrt(sigmae);
+      if (vertex.GetDist2()/sigmae>3.&&(!isGamma)) continue; 
+      Float_t densb0=0,densb1=0,densa0=0,densa1=0;
+      Int_t row0 = GetRowNumber(vertex.GetRr());
+      if (row0>15){
+       if (vertex.GetDist2()>0.2) continue;             
+       densb0     = track0->Density2(0,row0-5);          
+       densb1     = track1->Density2(0,row0-5);         
+       if (densb0>0.3|| densb1>0.3) continue;            //clusters before vertex
+       densa0     = track0->Density2(row0+5,row0+40);    
+       densa1     = track1->Density2(row0+5,row0+40);    
+       if ((densa0<0.4|| densa1<0.4)&&(!isGamma)) continue;            //missing clusters after vertex
+      }
+      else{
+       densa0     = track0->Density2(0,40);  //cluster density
+       densa1     = track1->Density2(0,40);  //cluster density
+       if ((vertex.GetRr()<80&&densa0+densa1<1.)&&(!isGamma)) continue;
+      }
+      vertex.SetLab(0,track0->GetLabel());
+      vertex.SetLab(1,track1->GetLabel());
+      vertex.SetChi2After((densa0+densa1)*0.5);
+      vertex.SetChi2Before((densb0+densb1)*0.5);
+      vertex.SetIndex(0,i);
+      vertex.SetIndex(1,j);
+      vertex.SetStatus(1); // TPC v0 candidate
+      vertex.SetRp(track0->fTPCr);
+      vertex.SetRm(track1->fTPCr);
+      tpcv0s->AddLast(new AliESDV0MI(vertex));      
+      ncandidates++;
+      {
+       Int_t eventNr = esd->GetEventNumber();
+       Double_t radiusm= (delta[0]<delta[1])? TMath::Sqrt(radius[0]):TMath::Sqrt(radius[1]);  
+       Double_t deltam= (delta[0]<delta[1])? TMath::Sqrt(delta[0]):TMath::Sqrt(delta[1]);  
+       cstream<<"V0"<<
+         "Event="<<eventNr<<
+         "vertex.="<<&vertex<<
+         "Tr0.="<<track0<<
+         "lab0="<<track0->fLab<<
+         "Helix0.="<<&helixes[i]<<     
+         "Tr1.="<<track1<<
+         "lab1="<<track1->fLab<<
+         "Helix1.="<<&helixes[j]<<
+         "pointAngle="<<pointAngle<<
+         "pointAngle2="<<pointAngle2<<
+         "dca0="<<dca[i]<<
+         "dca1="<<dca[j]<<
+         "z0="<<z0[i]<<
+         "z1="<<z0[j]<<
+         "zvertex="<<zvertex<<
+         "circular0="<<track0->fCircular<<
+         "circular1="<<track1->fCircular<<
+         "npoints="<<npoints<<
+         "radius0="<<radius[0]<<
+         "delta0="<<delta[0]<<
+         "radius1="<<radius[1]<<
+         "delta1="<<delta[1]<<
+         "radiusm="<<radiusm<<
+         "deltam="<<deltam<<
+         "sdcar0="<<sdcar[i]<<
+         "sdcar1="<<sdcar[j]<<
+         "cdcar0="<<cdcar[i]<<
+         "cdcar1="<<cdcar[j]<<
+         "pulldcar0="<<pulldcar[i]<<
+         "pulldcar1="<<pulldcar[j]<<
+         "pulldcaz0="<<pulldcaz[i]<<
+         "pulldcaz1="<<pulldcaz[j]<<
+         "pulldca0="<<pulldca[i]<<
+         "pulldca1="<<pulldca[j]<<
+         "densb0="<<densb0<<
+         "densb1="<<densb1<<
+         "densa0="<<densa0<<
+         "densa1="<<densa1<<
+         "sigmae="<<sigmae<<
+         "\n";
+      }
+    }
+  }    
+  Float_t *quality = new Float_t[ncandidates];
+  Int_t *indexes = new Int_t[ncandidates];
+  Int_t naccepted =0;
+  for (Int_t i=0;i<ncandidates;i++){
+    quality[i]     = 0; 
+    AliESDV0MI *v0 = (AliESDV0MI*)tpcv0s->At(i);
+    quality[i]     = 1./(1.00001-v0->GetPointAngle());   //base point angle
+    // quality[i]    /= (0.5+v0->GetDist2());  
+    // quality[i]    *= v0->GetChi2After();               //density factor
+    Double_t minpulldca = TMath::Min(2.+pulldca[v0->GetIndex(0)],(2.+pulldca[v0->GetIndex(1)]) );     //pull
+    Int_t index0 = v0->GetIndex(0);
+    Int_t index1 = v0->GetIndex(1);
+    AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
+    AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
+    if (track0->fTPCr[0]>0.3&&track1->fTPCr[0]>0.3&&v0->GetAnglep()[2]<0.15) quality[i]+=1000000;              // gamma conversion candidate
+    if (track0->fTPCr[4]>0.9||track1->fTPCr[4]>0.9&&minpulldca>4) quality[i]*=10;    // lambda candidate candidate
+  }
+
+  TMath::Sort(ncandidates,quality,indexes,kTRUE);
+  //
+  //
+  for (Int_t i=0;i<ncandidates;i++){
+    AliESDV0MI * v0 = (AliESDV0MI*)tpcv0s->At(indexes[i]);
+    if (!v0) continue;
+    Int_t index0 = v0->GetIndex(0);
+    Int_t index1 = v0->GetIndex(1);
+    AliTPCseed * track0 = (AliTPCseed*)array->At(index0);
+    AliTPCseed * track1 = (AliTPCseed*)array->At(index1);
+    if (!track0||!track1) {
+      printf("Bug\n");
+      continue;
+    }
+    Bool_t accept =kTRUE;  //default accept
+    Int_t *v0indexes0 = track0->GetV0Indexes();
+    Int_t *v0indexes1 = track1->GetV0Indexes();
+    //
+    Int_t order0 = (v0indexes0[0]!=0) ? 1:0;
+    Int_t order1 = (v0indexes1[0]!=0) ? 1:0;    
+    if (v0indexes0[1]!=0) order0 =2;
+    if (v0indexes1[1]!=0) order1 =2;      
+    //
+    if (v0indexes0[2]!=0) {order0=3; accept=kFALSE;}
+    if (v0indexes0[2]!=0) {order1=3; accept=kFALSE;}
+    //
+    AliESDV0MI * v02 = v0;
+    if (accept){
+      v0->SetOrder(0,order0);
+      v0->SetOrder(1,order1);
+      v0->SetOrder(1,order0+order1);     
+      Int_t index = esd->AddV0MI(v0);
+      v02 = esd->GetV0MI(index);
+      v0indexes0[order0]=index;
+      v0indexes1[order1]=index;
+      naccepted++;
+    }
+    {
+      Int_t eventNr = esd->GetEventNumber();
+      cstream<<"V02"<<
+       "Event="<<eventNr<<
+       "vertex.="<<v0<<        
+       "vertex2.="<<v02<<
+       "Tr0.="<<track0<<
+       "lab0="<<track0->fLab<<
+       "Tr1.="<<track1<<
+       "lab1="<<track1->fLab<<
+       "dca0="<<dca[index0]<<
+       "dca1="<<dca[index1]<<
+       "order0="<<order0<<
+       "order1="<<order1<<
+       "accept="<<accept<<
+       "quality="<<quality[i]<<
+       "pulldca0="<<pulldca[index0]<<
+       "pulldca1="<<pulldca[index1]<<
+       "index="<<i<<
+       "\n";
+    }
+  }    
+
+
+  //
+  //
+  delete []quality;
+  delete []indexes;
+//
+  delete [] isPrim;
+  delete [] pulldca;
+  delete [] pulldcaz;
+  delete [] pulldcar;
+  delete [] cdcar;
+  delete [] sdcar;
+  delete [] dca;
+  //
+  delete[] z0;
+  delete[] alpha;
+  delete[] sign;
+  delete[] helixes;
+  printf("TPC V0 finder : naccepted\t%d\tncandidates\t%d\tntracks\t%d\tnall\t%d\n",naccepted,ncandidates,ntracks,nall);
+  timer.Print();
+}
+
 Int_t AliTPCtrackerMI::RefitKink(AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &kink)
 {
   //
@@ -5115,12 +5295,47 @@ Int_t AliTPCtrackerMI::Clusters2Tracks() {
   //
   RemoveUsed2(fSeeds,0.85,0.85,0);
   FindKinks(fSeeds,fEvent);
-  RemoveUsed2(fSeeds,0.5,0.4,50);
+  RemoveUsed2(fSeeds,0.5,0.4,20);
+ //  //
+//   // refit short tracks
+//   //
+  Int_t nseed=fSeeds->GetEntriesFast();
+//   for (Int_t i=0; i<nseed; i++) {
+//     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;    
+//     if (!pt) continue;    
+//     Int_t nc=t.GetNumberOfClusters();
+//     if (nc<15) {
+//       delete fSeeds->RemoveAt(i);
+//       continue;
+//     }
+//     if (pt->GetKinkIndexes()[0]!=0) continue; // ignore kink candidates
+//     if (nc>100) continue;                     // hopefully, enough for ITS
+//     AliTPCseed *seed2 = new AliTPCseed(*pt);
+//     //seed2->Reset(kFALSE);
+//     //pt->ResetCovariance();
+//     seed2->Modify(1);
+//     FollowBackProlongation(*seed2,158);
+//     //seed2->Reset(kFALSE);
+//     seed2->Modify(10);
+//     FollowProlongation(*seed2,0);
+//     TTreeSRedirector &cstream = *fDebugStreamer;
+//     cstream<<"Crefit"<<
+//       "Tr0.="<<pt<<
+//       "Tr1.="<<seed2<<
+//       "\n";     
+//     if (seed2->fN>pt->fN){
+//       delete fSeeds->RemoveAt(i);
+//       fSeeds->AddAt(seed2,i);
+//     }else{
+//       delete seed2;
+//     }
+//   }
+//   RemoveUsed2(fSeeds,0.6,0.6,50);
+
+//  FindV0s(fSeeds,fEvent);  
   //RemoveDouble(fSeeds,0.2,0.6,11);
-  //  RemoveUsed(fSeeds,0.5,0.5,6);
 
   //
-  Int_t nseed=fSeeds->GetEntriesFast();
   Int_t found = 0;
   for (Int_t i=0; i<nseed; i++) {
     AliTPCseed *pt=(AliTPCseed*)fSeeds->UncheckedAt(i), &t=*pt;    
@@ -6116,654 +6331,12 @@ AliTPCclusterMI * AliTPCtrackerMI::AliTPCRow::FindNearest3(Double_t y, Double_t
 
 
 
-AliTPCseed::AliTPCseed():AliTPCtrack(){
-  //
-  fRow=0; 
-  fRemoval =0; 
-  for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
-  for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
-  for (Int_t i=0;i<3;i++)   fKinkIndexes[i]=0;
-  fPoints = 0;
-  fEPoints = 0;
-  fNFoundable =0;
-  fNShared  =0;
-  fRemoval = 0;
-  fSort =0;
-  fFirstPoint =0;
-  fNoCluster =0;
-  fBSigned = kFALSE;
-  fSeed1 =-1;
-  fSeed2 =-1;
-  fCurrentCluster =0;
-  fCurrentSigmaY2=0;
-  fCurrentSigmaZ2=0;
-}
-AliTPCseed::AliTPCseed(const AliTPCseed &s):AliTPCtrack(s){
-  //---------------------
-  // dummy copy constructor
-  //-------------------------
-  for (Int_t i=0;i<160;i++) fClusterPointer[i] = s.fClusterPointer[i];
-  for (Int_t i=0;i<160;i++) fIndex[i] = s.fIndex[i];
-  fPoints  = 0;
-  fEPoints = 0;
-}
-AliTPCseed::AliTPCseed(const AliTPCtrack &t):AliTPCtrack(t){
-  //
-  //copy constructor
-  fPoints = 0;
-  fEPoints = 0;
-  fNShared  =0; 
-  //  fTrackPoints =0;
-  fRemoval =0;
-  fSort =0;
-  for (Int_t i=0;i<3;i++)   fKinkIndexes[i]=t.GetKinkIndex(i);
-
-  for (Int_t i=0;i<160;i++) {
-    fClusterPointer[i] = 0;
-    Int_t index = t.GetClusterIndex(i);
-    if (index>=-1){ 
-      SetClusterIndex2(i,index);
-    }
-    else{
-      SetClusterIndex2(i,-3); 
-    }    
-  }
-  fFirstPoint =0;
-  fNoCluster =0;
-  fBSigned = kFALSE;
-  fSeed1 =-1;
-  fSeed2 =-1;
-  fCurrentCluster =0;
-  fCurrentSigmaY2=0;
-  fCurrentSigmaZ2=0;
-}
-/*
-AliTPCseed::AliTPCseed(const AliTPCtrack &t, Double_t a):AliTPCtrack(t,a){
-  //
-  //copy constructor
-  fRow=0;
-  for (Int_t i=0;i<160;i++) {
-    fClusterPointer[i] = 0;
-    Int_t index = t.GetClusterIndex(i);
-    SetClusterIndex2(i,index);
-  }
-  for (Int_t i=0;i<3;i++)   fKinkIndexes[i]=0;
-  fPoints = 0;
-  fEPoints = 0;
-  fNFoundable =0; 
-  fNShared  =0; 
-  //  fTrackPoints =0;
-  fRemoval =0;
-  fSort = 0;
-  fFirstPoint =0;
-  fNoCluster =0;
-  fBSigned = kFALSE;
-  fSeed1 =-1;
-  fSeed2 =-1;
-  fCurrentCluster =0;
-  fCurrentSigmaY2=0;
-  fCurrentSigmaZ2=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) {
-   //
-  //
-  //constructor
-  fRow =0;
-  for (Int_t i=0;i<200;i++) SetClusterIndex2(i,-3);
-  for (Int_t i=0;i<160;i++) fClusterPointer[i]=0;
-  for (Int_t i=0;i<3;i++)   fKinkIndexes[i]=0;
-  fPoints = 0;
-  fEPoints = 0;
-  fNFoundable =0;
-  fNShared  = 0;
-  //  fTrackPoints =0;
-  fRemoval =0;
-  fSort =0;
-  fFirstPoint =0;
-  //  fHelixIn = new TClonesArray("AliHelix",0);
-  //fHelixOut = new TClonesArray("AliHelix",0);
-  fNoCluster =0;
-  fBSigned = kFALSE;
-  fSeed1 =-1;
-  fSeed2 =-1;
-  fCurrentCluster =0;
-  fCurrentSigmaY2=0;
-  fCurrentSigmaZ2=0;
-}
-
-AliTPCseed::~AliTPCseed(){
-  //
-  // destructor
-  if (fPoints) delete fPoints;
-  fPoints =0;
-  if (fEPoints) delete fEPoints;
-  fEPoints = 0;
-  fNoCluster =0;
-}
-
-AliTPCTrackerPoint * AliTPCseed::GetTrackPoint(Int_t i)
-{
-  //
-  // 
-  return &fTrackPoints[i];
-}
-
-void AliTPCseed::RebuildSeed()
-{
-  //
-  // rebuild seed to be ready for storing
-  AliTPCclusterMI cldummy;
-  cldummy.SetQ(0);
-  AliTPCTrackPoint pdummy;
-  pdummy.GetTPoint().fIsShared = 10;
-  for (Int_t i=0;i<160;i++){
-    AliTPCclusterMI * cl0 = fClusterPointer[i];
-    AliTPCTrackPoint *trpoint = (AliTPCTrackPoint*)fPoints->UncheckedAt(i);     
-    if (cl0){
-      trpoint->GetTPoint() = *(GetTrackPoint(i));
-      trpoint->GetCPoint() = *cl0;
-      trpoint->GetCPoint().SetQ(TMath::Abs(cl0->GetQ()));
-    }
-    else{
-      *trpoint = pdummy;
-      trpoint->GetCPoint()= cldummy;
-    }
-    
-  }
-
-}
-
-
-Double_t AliTPCseed::GetDensityFirst(Int_t n)
-{
-  //
-  //
-  // return cluster for n rows bellow first point
-  Int_t nfoundable = 1;
-  Int_t nfound      = 1;
-  for (Int_t i=fLastPoint-1;i>0&&nfoundable<n; i--){
-    Int_t index = GetClusterIndex2(i);
-    if (index!=-1) nfoundable++;
-    if (index>0) nfound++;
-  }
-  if (nfoundable<n) return 0;
-  return Double_t(nfound)/Double_t(nfoundable);
-
-}
-
-
-void AliTPCseed::GetClusterStatistic(Int_t first, Int_t last, Int_t &found, Int_t &foundable, Int_t &shared, Bool_t plus2)
-{
-  // get cluster stat.  on given region
-  //
-  found       = 0;
-  foundable   = 0;
-  shared      =0;
-  for (Int_t i=first;i<last; i++){
-    Int_t index = GetClusterIndex2(i);
-    if (index!=-1) foundable++;
-    if (fClusterPointer[i]) {
-      found++;
-    }
-    else 
-      continue;
-
-    if (fClusterPointer[i]->IsUsed(10)) {
-      shared++;
-      continue;
-    }
-    if (!plus2) continue; //take also neighborhoud
-    //
-    if ( (i>0) && fClusterPointer[i-1]){
-      if (fClusterPointer[i-1]->IsUsed(10)) {
-       shared++;
-       continue;
-      }
-    }
-    if ( fClusterPointer[i+1]){
-      if (fClusterPointer[i+1]->IsUsed(10)) {
-       shared++;
-       continue;
-      }
-    }
-    
-  }
-  //if (shared>found){
-    //Error("AliTPCseed::GetClusterStatistic","problem\n");
-  //}
-}
-
-//_____________________________________________________________________________
-void AliTPCseed::CookdEdx(Double_t low, Double_t up,Int_t i1, Int_t i2, Bool_t onlyused) {
-  //-----------------------------------------------------------------
-  // 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+i1;i<i2;i+=4)
-      {
-       Int_t index = fIndex[i];
-       if (index<0||index&0x8000) continue;
-
-       //AliTPCTrackPoint * point = (AliTPCTrackPoint *) arr.At(i);
-       AliTPCTrackerPoint * point = GetTrackPoint(i);
-       //AliTPCTrackerPoint * pointm = GetTrackPoint(i-1);
-       //AliTPCTrackerPoint * pointp = 0;
-       //if (i<159) pointp = GetTrackPoint(i+1);
-
-       if (point==0) continue;
-       AliTPCclusterMI * cl = fClusterPointer[i];
-       if (cl==0) continue;    
-       if (onlyused && (!cl->IsUsed(10))) continue;
-       if (cl->IsUsed(11)) {
-         fNShared++;
-         continue;
-       }
-       Int_t   type   = cl->GetType();
-       //if (point->fIsShared){
-       //  fNShared++;
-       //  continue;
-       //}
-       //if (pointm) 
-       //  if (pointm->fIsShared) continue;
-       //if (pointp) 
-       //  if (pointp->fIsShared) continue;
-
-       if (type<0) continue;
-       //if (type>10) continue;       
-       //if (point->GetErrY()==0) continue;
-       //if (point->GetErrZ()==0) continue;
-
-       //Float_t ddy = (point->GetY()-cl->GetY())/point->GetErrY();
-       //Float_t ddz = (point->GetZ()-cl->GetZ())/point->GetErrZ();
-       //if ((ddy*ddy+ddz*ddz)>10) continue; 
-
-
-       //      if (point->GetCPoint().GetMax()<5) continue;
-       if (cl->GetMax()<5) continue;
-       Float_t angley = point->GetAngleY();
-       Float_t anglez = point->GetAngleZ();
-
-       Float_t rsigmay2 =  point->GetSigmaY();
-       Float_t rsigmaz2 =  point->GetSigmaZ();
-       /*
-       Float_t ns = 1.;
-       if (pointm){
-         rsigmay +=  pointm->GetTPoint().GetSigmaY();
-         rsigmaz +=  pointm->GetTPoint().GetSigmaZ();
-         ns+=1.;
-       }
-       if (pointp){
-         rsigmay +=  pointp->GetTPoint().GetSigmaY();
-         rsigmaz +=  pointp->GetTPoint().GetSigmaZ();
-         ns+=1.;
-       }
-       rsigmay/=ns;
-       rsigmaz/=ns;
-       */
-
-       Float_t rsigma = TMath::Sqrt(rsigmay2*rsigmaz2);
-
-       Float_t ampc   = 0;     // normalization to the number of electrons
-       if (i>64){
-         //      ampc = 1.*point->GetCPoint().GetMax();
-         ampc = 1.*cl->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);
-         //      Float_t zfactor = (AliTPCReconstructor::GetCtgRange()-0.0004*TMath::Abs(point->GetCPoint().GetZ()));
-         //ampc               *=zfactor; 
-       }
-       else{ 
-         //ampc = 1.0*point->GetCPoint().GetMax(); 
-         ampc = 1.0*cl->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);
-         //    Float_t zfactor = (1.02-0.000*TMath::Abs(point->GetCPoint().GetZ()));
-         //ampc               *=zfactor; 
-
-       }
-       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); 
-       //      Float_t z = TMath::Abs(cl->GetZ());
-       if (i<64) {
-         ampc /= 0.6;
-         //ampc /= (1+0.0008*z);
-       } else
-         if (i>128){
-           ampc /=1.5;
-           //ampc /= (1+0.0008*z);
-         }else{
-           //ampc /= (1+0.0008*z);
-         }
-       
-       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 = 50;
-    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;
-
-}
-
-
-
-/*
-
-
 
-void AliTPCseed::CookdEdx2(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];
-  Bool_t inlimit[200];
-  for (Int_t i=0;i<200;i++) inlimit[i]=kFALSE;
-  for (Int_t i=0;i<200;i++) amp[i]=10000;
-  for (Int_t i=0;i<200;i++) angular[i]= 1;;
-  
-
-  //
-  Float_t meanlog = 100.;
-  Int_t indexde[4]={0,64,128,160};
+// AliTPCTrackerPoint * AliTPCseed::GetTrackPoint(Int_t i)
+// {
+//   //
+//   // 
+//   return &fTrackPoints[i];
+// }
 
-  Float_t amean     =0;
-  Float_t asigma    =0;
-  Float_t anc       =0;
-  Float_t anorm     =0;
-
-  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<3; of++){    
-  //  for (Int_t i=indexde[of];i<indexde[of+1];i++)
-  for (Int_t i =0; i<160;i++)
-    {
-       AliTPCTrackPoint * point = GetTrackPoint(i);
-       if (point==0) continue;
-       if (point->fIsShared){
-         fNShared++;     
-         continue;
-       }
-       Int_t   type   = point->GetCPoint().GetType();
-       if (type<0) continue;
-       if (point->GetCPoint().GetMax()<5) continue;
-       Float_t angley = point->GetTPoint().GetAngleY();
-       Float_t anglez = point->GetTPoint().GetAngleZ();
-       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 =  point->GetCPoint().GetMax();
-       }
-       else{ 
-         ampc = point->GetCPoint().GetMax(); 
-       }
-       ampc *= 2.0;     // put mean value to channel 50
-       //      ampc *= 0.565;     // put mean value to channel 50
-
-       Float_t w      =  1.;
-       Float_t z = TMath::Abs(point->GetCPoint().GetZ());
-       if (i<64) {
-         ampc /= 0.63;
-       } else
-         if (i>128){
-           ampc /=1.51;
-         }             
-       if (type<0) {  //amp at the border - lower weight                 
-         continue;
-       }
-       if (rsigma>1.5) ampc/=1.3;  // if big backround
-       angular[i]    = TMath::Sqrt(1.+angley*angley+anglez*anglez);
-       amp[i]        = ampc/angular[i];
-       weight[i]     = w;
-       anc++;
-    }
 
-  TMath::Sort(159,amp,index,kFALSE);
-  for (Int_t i=int(anc*low+0.5);i<int(anc*up+0.5);i++){      
-    inlimit[index[i]] = kTRUE;  // take all clusters
-  }
-  
-  //  meanlog = amp[index[Int_t(anc*0.3)]];
-  meanlog =10000.;
-  for (Int_t of =0; of<3; of++){    
-    Float_t sumamp=0;
-    Float_t sumamp2=0;
-    Float_t sumw=0;    
-   for (Int_t i=indexde[of];i<indexde[of+1];i++)
-      {
-       if (inlimit[i]==kFALSE) continue;
-       Float_t ampl      = amp[i];
-       ///angular[i];
-       ampl              = meanlog*TMath::Log(1.+ampl/meanlog);
-       //
-       sumw    += weight[i]; 
-       sumamp  += weight[i]*ampl;
-       sumamp2 += weight[i]*ampl*ampl;
-       norm[of]    += angular[i]*weight[i];
-       nc[of]++;
-      }
-   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;
-   }
-  }
-    
-  Float_t dedx =0;
-  fSdEdx =0;
-  fMAngular =0;
-  //
-  Int_t norm2 = 0;
-  Int_t norm3 = 0;
-  Float_t www[3] = {12.,14.,17.};
-  //Float_t www[3] = {1.,1.,1.};
-
-  for (Int_t i =0;i<3;i++){
-    if (nc[i]>2&&nc[i]<1000){
-      dedx      += mean[i] *nc[i]*www[i]/sigma[i];
-      fSdEdx    += sigma[i]*(nc[i]-2)*www[i]/sigma[i];
-      fMAngular += norm[i] *nc[i];    
-      norm2     += nc[i]*www[i]/sigma[i];
-      norm3     += (nc[i]-2)*www[i]/sigma[i];
-    }
-    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;
-  Float_t norm4 = 0;
-  for (Int_t i =0;i<3;i++){
-    if (nc[i]>2&&nc[i]<1000&&sigma[i]>3){
-      //mean[i]   = mean[i]*(1+0.08*(sigma[i]/(fSdEdx)-1.));
-      dedx      += mean[i] *(nc[i])/(sigma[i]);
-      norm4     += (nc[i])/(sigma[i]);
-    }
-    fDEDX[i]  = mean[i];                
-  }
-  if (norm4>0) dedx /= norm4;
-  
-
-  
-  SetdEdx(dedx);
-    
-  //mi deDX
-
-}
-
-*/
index 211e941497e902570cdde04df6620cf2b5fa2ae2..53836a25209bb915bc2197073be93ecbbf5dde55 100644 (file)
 //   Origin: 
 //-------------------------------------------------------
 
+#include <TError.h>
 #include "AliTracker.h"
-#include "AliTPCtrack.h"
-#include "AliComplexCluster.h"
+#include "AliTPCreco.h"
+#include "AliPID.h"
+
 
 class TFile;
 class AliTPCParam;
@@ -25,87 +27,7 @@ class AliTPCTrackerPoint;
 class AliESD;   
 class TTree;
 class AliESDkink;
-
-class AliTPCseed : public AliTPCtrack {
-  friend class AliTPCtrackerMI;
-  public:  
-     AliTPCseed();
-     virtual ~AliTPCseed();
-     AliTPCseed(const AliTPCtrack &t);
-     AliTPCseed(const AliTPCseed &s);
-     //AliTPCseed(const AliTPCseed &t, Double_t a);
-     AliTPCseed(UInt_t index, const Double_t xx[5], 
-                const Double_t cc[15], Double_t xr, Double_t alpha);     
-     Int_t Compare(const TObject *o) const;
-     void Reset(Bool_t all = kTRUE);
-     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);
-     AliTPCTrackerPoint * GetTrackPoint(Int_t i);
-     void RebuildSeed(); // rebuild seed to be ready for storing
-     Double_t GetDensityFirst(Int_t n);
-     Double_t GetSigma2C() const {return fC44;};
-     void GetClusterStatistic(Int_t first, Int_t last, Int_t &found, Int_t &foundable, Int_t &shared, Bool_t plus2);
-     
-     void Modify(Double_t factor);
-     void SetClusterIndex2(Int_t row, Int_t index) {
-       fIndex[row] = index;
-     }
-     Int_t  GetClusterIndex2(Int_t row) const {
-       return fIndex[row];
-     }
-     Int_t GetClusterSector(Int_t row) const {
-       Int_t pica = -1;
-       if (fIndex[row]>=0) pica =  ((fIndex[row]&0xff000000)>>24);
-       return pica;
-     }
-    
-     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, Int_t i1=0, Int_t i2=159, Bool_t onlyused = kFALSE);
-     //     void CookdEdx2(Double_t low=0.05, Double_t up=0.70);
-     Bool_t IsActive() const { return !(fRemoval);}
-     void Desactivate(Int_t reason){ fRemoval = reason;} 
-     //
-     //
- private:
-     AliTPCseed & operator = (const AliTPCseed &);
-     AliESDtrack * fEsd; //!
-     AliTPCclusterMI*   fClusterPointer[160];  //! array of cluster pointers  - 
-     TClonesArray * fPoints;              // array with points along the track
-     TClonesArray * fEPoints;             // array with exact points - calculated in special macro not used in tracking
-     //---CURRENT VALUES
-     Int_t fRow;                 //!current row number  
-     Int_t fSector;              //!current sector number
-     Int_t fRelativeSector;      //! index of current relative sector
-     Float_t fCurrentSigmaY2;    //!expected current cluster sigma Y
-     Float_t fCurrentSigmaZ2;    //!expected current cluster sigma Z
-     Float_t fErrorY2;           //!sigma of current cluster 
-     Float_t fErrorZ2;           //!sigma of current cluster    
-     AliTPCclusterMI * fCurrentCluster; //!pointer to the current cluster for prolongation
-     Int_t   fCurrentClusterIndex1; //! index of the current cluster
-     Bool_t  fInDead;            //! indicate if the track is in dead zone
-     Bool_t  fIsSeeding;         //!indicates if it is proces of seeading
-     Int_t   fNoCluster;         //!indicates number of rows without clusters
-     Int_t   fSort;              //!indicate criteria for sorting
-     Bool_t  fBSigned;        //indicates that clusters of this trackes are signed to be used
-     //
-     //
-     Float_t fDEDX[4];         // dedx according padrows
-     Float_t fSDEDX[4];        // sdedx according padrows
-     Int_t   fNCDEDX[4];       // number of clusters for dedx measurment
-     //
-     Int_t   fSeedType;         //seeding type
-     Int_t   fSeed1;            //first row for seeding
-     Int_t   fSeed2;            //last row for seeding
-     Int_t   fOverlapLabels[12];  //track labels and the length of the  overlap     
-     Float_t fMAngular;        // mean angular factor
-     AliTPCTrackerPoint  fTrackPoints[160];  //!track points - array track points
-   
-     ClassDef(AliTPCseed,1)  
-};
-
-
+class TTreeSRedirector;
 
 
 class AliTPCtrackerMI : public AliTracker {
@@ -134,6 +56,7 @@ public:
   void DeleteSeeds();
   void SetDebug(Int_t debug){ fDebug = debug;}
   void FindKinks(TObjArray * array, AliESD * esd);
+  void FindV0s(TObjArray * array, AliESD * esd);
   void UpdateKinkQualityM(AliTPCseed * seed);
   void UpdateKinkQualityD(AliTPCseed * seed);
   Int_t CheckKinkPoint(AliTPCseed*seed, AliTPCseed &mother, AliTPCseed &daughter, AliESDkink &kink);
@@ -338,6 +261,7 @@ private:
    Double_t fYMax[200];                // max y for given pad row
    Double_t fPadLength[200];                // max y for given pad row
    const AliTPCParam *fParam;          //pointer to the parameters
+   TTreeSRedirector *fDebugStreamer;     //!debug streamer
    ClassDef(AliTPCtrackerMI,1) 
 };
 
index a61fe3a360634b9fc1ef00108bc7731658471bc3..2bf4acb6217bbd27da5a7db917325e1a9ea37c3c 100644 (file)
@@ -4,7 +4,7 @@ SRCS:=  AliTPCcluster.cxx \
         AliClustersArray.cxx AliTPCClustersArray.cxx \
        AliTPCclusterer.cxx AliTPCclustererMI.cxx \
        AliTPCtrack.cxx AliTPCtracker.cxx \
-       AliTPCpolyTrack.cxx AliTPCtrackerMI.cxx \
+       AliTPCpolyTrack.cxx  AliTPCseed.cxx AliTPCtrackerMI.cxx \
        AliTPCPid.cxx AliTPCtrackPid.cxx AliTPCpidESD.cxx \
        AliTPCReconstructor.cxx