When calculating a*a-b*b the form (a-b)*(a+b) is usually more numerically stable.
[u/mrichter/AliRoot.git] / STEER / AliVertexerTracks.cxx
index 511ce4b..97a7626 100644 (file)
 #include <TTree.h>
 #include <TMatrixD.h>
 //---- AliRoot headers -----
+#include "AliLog.h"
 #include "AliStrLine.h"
 #include "AliExternalTrackParam.h"
+#include "AliNeutralTrackParam.h"
+#include "AliVEvent.h"
+#include "AliVTrack.h"
 #include "AliESDEvent.h"
 #include "AliESDtrack.h"
 #include "AliVertexerTracks.h"
@@ -46,65 +50,68 @@ AliVertexerTracks::AliVertexerTracks():
 TObject(),
 fVert(),
 fCurrentVertex(0),
+fMode(0),
 fFieldkG(-999.),
-fConstraint(kFALSE),
-fOnlyFitter(kFALSE),
-fMinTracks(1),
-fMinITSClusters(5),
 fTrkArraySel(),
 fIdSel(0),
 fTrksToSkip(0),
 fNTrksToSkip(0),
-fDCAcut(0),
-fAlgo(1),
-fNSigma(3),
+fConstraint(kFALSE),
+fOnlyFitter(kFALSE),
+fMinTracks(1),
+fMinClusters(5),
+fDCAcut(0.1),
+fDCAcutIter0(0.1),
+fNSigma(3.),
 fMaxd0z0(0.5),
+fMinDetFitter(100.),
+fMaxTgl(1000.),
 fITSrefit(kTRUE),
+fFiducialR(3.),
+fFiducialZ(30.),
 fnSigmaForUi00(1.5),
-fDebug(0)
+fAlgo(1),
+fAlgoIter0(4)
 {
 //
 // Default constructor
 //
   SetVtxStart();
   SetVtxStartSigma();
-  SetMinTracks();
-  SetMinITSClusters();
-  SetNSigmad0(); 
-  SetMaxd0z0(); 
 }
 //----------------------------------------------------------------------------
 AliVertexerTracks::AliVertexerTracks(Double_t fieldkG):
 TObject(),
 fVert(),
 fCurrentVertex(0),
-fFieldkG(-999.),
-fConstraint(kFALSE),
-fOnlyFitter(kFALSE),
-fMinTracks(1),
-fMinITSClusters(5),
+fMode(0),
+fFieldkG(fieldkG),
 fTrkArraySel(),
 fIdSel(0),
 fTrksToSkip(0),
 fNTrksToSkip(0),
-fDCAcut(0),
-fAlgo(1),
-fNSigma(3),
+fConstraint(kFALSE),
+fOnlyFitter(kFALSE),
+fMinTracks(1),
+fMinClusters(5),
+fDCAcut(0.1),
+fDCAcutIter0(0.1),
+fNSigma(3.),
 fMaxd0z0(0.5),
+fMinDetFitter(100.),
+fMaxTgl(1000.),
 fITSrefit(kTRUE),
+fFiducialR(3.),
+fFiducialZ(30.),
 fnSigmaForUi00(1.5),
-fDebug(0)
+fAlgo(1),
+fAlgoIter0(4)
 {
 //
 // Standard constructor
 //
   SetVtxStart();
   SetVtxStartSigma();
-  SetMinTracks();
-  SetMinITSClusters();
-  SetNSigmad0(); 
-  SetMaxd0z0();
-  SetFieldkG(fieldkG);
 }
 //-----------------------------------------------------------------------------
 AliVertexerTracks::~AliVertexerTracks() 
@@ -117,10 +124,10 @@ AliVertexerTracks::~AliVertexerTracks()
   if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
 }
 //----------------------------------------------------------------------------
-AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliESDEvent *esdEvent)
+AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliVEvent *vEvent)
 {
 //
-// Primary vertex for current ESD event
+// Primary vertex for current ESD or AOD event
 // (Two iterations: 
 //  1st with 5*fNSigma*sigma cut w.r.t. to initial vertex
 //      + cut on sqrt(d0d0+z0z0) if fConstraint=kFALSE  
@@ -128,55 +135,84 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliESDEvent *esdEvent)
 //
   fCurrentVertex = 0;
 
+  TString evtype = vEvent->IsA()->GetName();
+  Bool_t inputAOD = ((evtype=="AliAODEvent") ? kTRUE : kFALSE);
+
+  if(inputAOD && fMode==1) {
+    printf("Error : AliVertexerTracks: no TPC-only vertex from AOD\n"); 
+    TooFewTracks(); 
+    return fCurrentVertex;
+  }
+
   // accept 1-track case only if constraint is available
   if(!fConstraint && fMinTracks==1) fMinTracks=2;
 
-  // read tracks from ESD
-  Int_t nTrks = (Int_t)esdEvent->GetNumberOfTracks();
-  if(nTrks<1) {
+  // read tracks from AlivEvent
+  Int_t nTrks = (Int_t)vEvent->GetNumberOfTracks();
+  if(nTrks<fMinTracks) {
     TooFewTracks();
     return fCurrentVertex;
   } 
   //
 
   TDirectory * olddir = gDirectory;
-  TFile f("VertexerTracks.root","recreate");
+  TFile *f = 0;
+  if(nTrks>500) f = new TFile("VertexerTracks.root","recreate");
   TObjArray trkArrayOrig(nTrks);
   UShort_t *idOrig = new UShort_t[nTrks];
 
   Int_t nTrksOrig=0;
+  AliExternalTrackParam *t=0;
+  // loop on tracks
   for(Int_t i=0; i<nTrks; i++) {
-    AliESDtrack *esdt = esdEvent->GetTrack(i);
+    AliVTrack *track = (AliVTrack*)vEvent->GetTrack(i);
     // check tracks to skip
     Bool_t skipThis = kFALSE;
     for(Int_t j=0; j<fNTrksToSkip; j++) { 
-      if(esdt->GetID()==fTrksToSkip[j]) {
-       if(fDebug) printf("skipping track: %d\n",i);
+      if(track->GetID()==fTrksToSkip[j]) {
+       AliDebug(1,Form("skipping track: %d",i));
        skipThis = kTRUE;
       }
     }
     if(skipThis) continue;
-    if(fITSrefit) {
-      if(!(esdt->GetStatus()&AliESDtrack::kITSrefit)) continue;
-      // check number of clusters in ITS
-      if(esdt->GetNcls(0)<fMinITSClusters) continue;
+
+    if(!inputAOD) {  // ESD
+      AliESDtrack* esdt = (AliESDtrack*)track;
+      if(esdt->GetNcls(fMode) < fMinClusters) continue;
+      if(fMode==0) {        // ITS mode
+       if(fITSrefit && !(esdt->GetStatus()&AliESDtrack::kITSrefit)) continue;
+       Double_t x,p[5],cov[15];
+       esdt->GetExternalParameters(x,p);
+       esdt->GetExternalCovariance(cov);
+       t = new AliExternalTrackParam(x,esdt->GetAlpha(),p,cov);
+      } else if(fMode==1) { // TPC mode
+       t = (AliExternalTrackParam*)esdt->GetTPCInnerParam();
+       if(!t) continue;
+       Double_t radius = 2.8; //something less than the beam pipe radius
+       if(!PropagateTrackTo(t,radius)) continue;
+      }
+    } else {          // AOD (only ITS mode)
+      Int_t ncls0=0;
+      for(Int_t l=0;l<6;l++) if(TESTBIT(track->GetITSClusterMap(),l)) ncls0++;
+      if(ncls0 < fMinClusters) continue;
+      t = new AliExternalTrackParam(track);
     }
-    Double_t x,p[5],cov[15];
-    esdt->GetExternalParameters(x,p);
-    esdt->GetExternalCovariance(cov);
-    AliExternalTrackParam *t = new AliExternalTrackParam(x,esdt->GetAlpha(),p,cov);
     trkArrayOrig.AddLast(t);
-    idOrig[nTrksOrig]=(UShort_t)esdt->GetID();
+    idOrig[nTrksOrig]=(UShort_t)track->GetID();
     nTrksOrig++;
-  }
+  } // end loop on tracks
   
+  // call method that will reconstruct the vertex
   FindPrimaryVertex(&trkArrayOrig,idOrig);
 
-  trkArrayOrig.Delete();
+  if(fMode==0) trkArrayOrig.Delete();
   delete [] idOrig; idOrig=NULL;
-  f.Close();
-  gSystem->Unlink("VertexerTracks.root");
-  olddir->cd();
+
+  if(f) {
+    f->Close(); delete f; f = NULL;
+    gSystem->Unlink("VertexerTracks.root");
+    olddir->cd();
+  }
 
   return fCurrentVertex;
 }
@@ -197,11 +233,11 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(TObjArray *trkArrayOrig,
   // accept 1-track case only if constraint is available
   if(!fConstraint && fMinTracks==1) fMinTracks=2;
 
-  // read tracks from ESD
+  // read tracks from array
   Int_t nTrksOrig = (Int_t)trkArrayOrig->GetEntriesFast();
-  if(fDebug) printf("Initial number of tracks: %d\n",nTrksOrig);
-  if(nTrksOrig<=0) {
-    if(fDebug) printf("TooFewTracks\n");
+  AliDebug(1,Form("Initial number of tracks: %d",nTrksOrig));
+  if(nTrksOrig<fMinTracks) {
+    AliDebug(1,"TooFewTracks");
     TooFewTracks();
     return fCurrentVertex;
   } 
@@ -214,20 +250,28 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(TObjArray *trkArrayOrig,
     PrepareTracks(*trkArrayOrig,idOrig,0);
     if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
     Double_t cutsave = fDCAcut;  
-    fDCAcut = (fITSrefit ? 0.1 : 0.5);
-    VertexFinder(1); // using weights, cutting dca < fDCAcut
+    fDCAcut = fDCAcutIter0;
+    // vertex finder
+    switch (fAlgoIter0) {
+    case 1: StrLinVertexFinderMinDist(1); break;
+    case 2: StrLinVertexFinderMinDist(0); break;
+    case 3: HelixVertexFinder();          break;
+    case 4: VertexFinder(1);              break;
+    case 5: VertexFinder(0);              break;
+    default: printf("Wrong algorithm\n"); break;  
+    }
     fDCAcut = cutsave;
     if(fVert.GetNContributors()>0) {
       fVert.GetXYZ(fNominalPos);
       fNominalPos[0] = fVert.GetXv();
       fNominalPos[1] = fVert.GetYv();
       fNominalPos[2] = fVert.GetZv();
-      if(fDebug) printf("No mean vertex: VertexFinder gives (%f, %f, %f)\n",fNominalPos[0],fNominalPos[1],fNominalPos[2]);
+      AliDebug(1,Form("No mean vertex: VertexFinder gives (%f, %f, %f)",fNominalPos[0],fNominalPos[1],fNominalPos[2]));
     } else {
       fNominalPos[0] = 0.;
       fNominalPos[1] = 0.;
       fNominalPos[2] = 0.;
-      if(fDebug) printf("No mean vertex and VertexFinder failed\n");
+      AliDebug(1,"No mean vertex and VertexFinder failed");
     }
   }
   
@@ -247,7 +291,7 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(TObjArray *trkArrayOrig,
     if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
     fIdSel = new UShort_t[nTrksOrig];
     Int_t nTrksSel = PrepareTracks(*trkArrayOrig,idOrig,iter);
-    if(fDebug) printf("N tracks selected in iteration %d: %d\n",iter,nTrksSel);
+    AliDebug(1,Form("N tracks selected in iteration %d: %d",iter,nTrksSel));
     if(nTrksSel < fMinTracks) {
       TooFewTracks();
       return fCurrentVertex; 
@@ -256,7 +300,7 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(TObjArray *trkArrayOrig,
     // vertex finder
     if(!fOnlyFitter) {
       if(nTrksSel==1) {
-       if(fDebug) printf("Just one track\n");
+       AliDebug(1,"Just one track");
        OneTrackVertFinder();
       } else {
        switch (fAlgo) {
@@ -265,10 +309,10 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(TObjArray *trkArrayOrig,
         case 3: HelixVertexFinder();          break;
         case 4: VertexFinder(1);              break;
         case 5: VertexFinder(0);              break;
-        default: printf("Wrong algorithm\n"); break;  
+        default: printf("Wrong algorithm"); break;  
        }
       }
-      if(fDebug) printf(" Vertex finding completed\n");
+      AliDebug(1," Vertex finding completed");
     }
 
     // vertex fitter
@@ -297,10 +341,7 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(TObjArray *trkArrayOrig,
   fCurrentVertex->SetTitle(title.Data());
   //
 
-  if(fDebug) {
-    fCurrentVertex->PrintStatus();
-    fCurrentVertex->PrintIndices();
-  }
+  AliDebug(1,Form("xyz: %f %f %f; nc %d",fCurrentVertex->GetXv(),fCurrentVertex->GetYv(),fCurrentVertex->GetZv(),fCurrentVertex->GetNContributors()));
 
   // clean up
   delete [] fIdSel; fIdSel=NULL;
@@ -392,21 +433,24 @@ Double_t AliVertexerTracks::GetStrLinMinDist(Double_t *p0,Double_t *p1,Double_t
   Double_t x10=p0[0]-x0[0];
   Double_t y10=p0[1]-x0[1];
   Double_t z10=p0[2]-x0[2];
-  return ((x10*x10+y10*y10+z10*z10)*(x12*x12+y12*y12+z12*z12)-(x10*x12+y10*y12+z10*z12)*(x10*x12+y10*y12+z10*z12))/(x12*x12+y12*y12+z12*z12);
+  //  return ((x10*x10+y10*y10+z10*z10)*(x12*x12+y12*y12+z12*z12)-(x10*x12+y10*y12+z10*z12)*(x10*x12+y10*y12+z10*z12))/(x12*x12+y12*y12+z12*z12);
+  return ((y10*z12-z10*y12)*(y10*z12-z10*y12)+
+         (z10*x12-x10*z12)*(z10*x12-x10*z12)+
+         (x10*y12-y10*x12)*(x10*y12-y10*x12))
+    /(x12*x12+y12*y12+z12*z12);
 }
 //---------------------------------------------------------------------------
 void AliVertexerTracks::OneTrackVertFinder() 
 {
   // find vertex for events with 1 track, using DCA to nominal beam axis
-  if(fDebug) printf("Number of prepared tracks =%d - Call OneTrackVertFinder",fTrkArraySel.GetEntries());
+  AliDebug(1,Form("Number of prepared tracks =%d - Call OneTrackVertFinder",fTrkArraySel.GetEntries()));
   AliExternalTrackParam *track1;
   track1 = (AliExternalTrackParam*)fTrkArraySel.At(0);
-  Double_t field=GetFieldkG();
   Double_t alpha=track1->GetAlpha();
   Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
   Double_t pos[3],dir[3]; 
-  track1->GetXYZAt(mindist,field,pos);
-  track1->GetPxPyPzAt(mindist,field,dir);
+  track1->GetXYZAt(mindist,GetFieldkG(),pos);
+  track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
   AliStrLine *line1 = new AliStrLine(pos,dir);
   Double_t p1[3]={fNominalPos[0],fNominalPos[1],0.}; 
   Double_t p2[3]={fNominalPos[0],fNominalPos[1],10.}; 
@@ -434,7 +478,6 @@ void AliVertexerTracks::HelixVertexFinder()
   Double_t initPos[3];
   initPos[2] = 0.;
   for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
-  Double_t field=GetFieldkG();
 
   Int_t nacc = (Int_t)fTrkArraySel.GetEntriesFast();
 
@@ -456,7 +499,7 @@ void AliVertexerTracks::HelixVertexFinder()
     for(Int_t j=i+1; j<nacc; j++){
       track2 = (AliExternalTrackParam*)fTrkArraySel.At(j);
 
-      distCA=track2->PropagateToDCA(track1,field);
+      distCA=track2->PropagateToDCA(track1,GetFieldkG());
       if(fDCAcut<=0 ||(fDCAcut>0&&distCA<fDCAcut)){
        x=track1->GetX();
        alpha=track1->GetAlpha();
@@ -513,17 +556,15 @@ Int_t AliVertexerTracks::PrepareTracks(TObjArray &trkArrayOrig,
 //
 // Propagate tracks to initial vertex position and store them in a TObjArray
 //
-  if(fDebug) printf(" PrepareTracks()\n");
+  AliDebug(1," PrepareTracks()");
 
   Int_t nTrksOrig = (Int_t)trkArrayOrig.GetEntriesFast();
   Int_t nTrksSel = 0;
   Double_t maxd0rphi; 
-  Double_t maxd0z0 = (fITSrefit ? fMaxd0z0 : 10.*fMaxd0z0);
   Double_t sigmaCurr[3];
   Double_t normdistx,normdisty;
   Double_t d0z0[2],covd0z0[3]; 
   Double_t sigmad0;
-  Double_t field = GetFieldkG();
 
   AliESDVertex *initVertex = new AliESDVertex(fNominalPos,fNominalCov,1,1);
 
@@ -531,66 +572,131 @@ Int_t AliVertexerTracks::PrepareTracks(TObjArray &trkArrayOrig,
 
   AliExternalTrackParam *track=0;
 
+  // loop on tracks
   for(Int_t i=0; i<nTrksOrig; i++) {
-    track = new AliExternalTrackParam(*(AliExternalTrackParam*)trkArrayOrig.At(i));
-    // only TPC tracks in |eta|<0.9
-    if(!fITSrefit && TMath::Abs(track->GetTgl())>1.) {
-      if(fDebug) printf(" rejecting track with tgl = %f\n",track->GetTgl());
-      delete track;
-      continue;
+    AliExternalTrackParam *trackOrig=(AliExternalTrackParam*)trkArrayOrig.At(i);
+    if(trackOrig->Charge()!=0) { // normal tracks
+      track = new AliExternalTrackParam(*(AliExternalTrackParam*)trkArrayOrig.At(i));
+    } else { // neutral tracks (from a V0)
+      track = new AliNeutralTrackParam(*(AliNeutralTrackParam*)trkArrayOrig.At(i));
+    }
+    // tgl cut
+    if(TMath::Abs(track->GetTgl())>fMaxTgl) {
+      AliDebug(1,Form(" rejecting track with tgl = %f",track->GetTgl()));
+      delete track; continue;
     }
 
+    Bool_t propagateOK = kFALSE;
     // propagate track to vertex
     if(optImpParCut<2 || fOnlyFitter) { // optImpParCut==1 or 0
-      track->PropagateToDCA(initVertex,field,100.,d0z0,covd0z0);
+      propagateOK = track->PropagateToDCA(initVertex,GetFieldkG(),100.,d0z0,covd0z0);
     } else {              // optImpParCut==2
       fCurrentVertex->GetSigmaXYZ(sigmaCurr);
       normdistx = TMath::Abs(fCurrentVertex->GetXv()-fNominalPos[0])/TMath::Sqrt(sigmaCurr[0]*sigmaCurr[0]+fNominalCov[0]);
       normdisty = TMath::Abs(fCurrentVertex->GetYv()-fNominalPos[1])/TMath::Sqrt(sigmaCurr[1]*sigmaCurr[1]+fNominalCov[2]);
       if(normdistx < 3. && normdisty < 3. &&
         (sigmaCurr[0]+sigmaCurr[1])<(TMath::Sqrt(fNominalCov[0])+TMath::Sqrt(fNominalCov[2]))) {
-       track->PropagateToDCA(fCurrentVertex,field,100.,d0z0,covd0z0);
+       propagateOK = track->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
       } else {
-       track->PropagateToDCA(initVertex,field,100.,d0z0,covd0z0);
+       propagateOK = track->PropagateToDCA(initVertex,GetFieldkG(),100.,d0z0,covd0z0);
       }
     }
 
+    if(!propagateOK) { 
+      AliDebug(1,"     rejected");
+      delete track; continue; 
+    }
+
     sigmad0 = TMath::Sqrt(covd0z0[0]);
     maxd0rphi = fNSigma*sigmad0;
     if(optImpParCut==1) maxd0rphi *= 5.;
+    maxd0rphi = TMath::Min(maxd0rphi,fFiducialR); 
     //sigmad0z0 = TMath::Sqrt(covd0z0[0]+covd0z0[2]); // for future improvement
     //maxd0z0 = 10.*fNSigma*sigmad0z0;
 
-    if(fDebug) printf("trk %d; id %d; |d0| = %f;  d0 cut = %f; |z0| = %f; |d0|oplus|z0| = %f; d0z0 cut = %f\n",i,(Int_t)idOrig[i],TMath::Abs(d0z0[0]),maxd0rphi,TMath::Abs(d0z0[1]),TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1]),maxd0z0);
+    AliDebug(1,Form("trk %d; id %d; |d0| = %f;  d0 cut = %f; |z0| = %f; |d0|oplus|z0| = %f; d0z0 cut = %f",i,(Int_t)idOrig[i],TMath::Abs(d0z0[0]),maxd0rphi,TMath::Abs(d0z0[1]),TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1]),fMaxd0z0));
+
+
+    //---- track selection based on impact parameters ----//
+
+    // always reject tracks outside fiducial volume
+    if(TMath::Abs(d0z0[0])>fFiducialR || TMath::Abs(d0z0[1])>fFiducialZ) { 
+      AliDebug(1,"     rejected");
+      delete track; continue; 
+    }
+
+    // during iterations 1 and 2 , reject tracks with d0rphi > maxd0rphi
+    if(optImpParCut>0 && TMath::Abs(d0z0[0])>maxd0rphi) { 
+      AliDebug(1,"     rejected");
+      delete track; continue; 
+    }
 
     // if fConstraint=kFALSE, during iterations 1 and 2 ||
     // if fConstraint=kTRUE, during iteration 2,
-    // select tracks with d0oplusz0 < maxd0z0
+    // select tracks with d0oplusz0 < fMaxd0z0
     if((!fConstraint && optImpParCut>0 && fVert.GetNContributors()>0) ||
        ( fConstraint && optImpParCut==2)) {
       if(nTrksOrig>=3 && 
-        TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1])>maxd0z0) { 
-       if(fDebug) printf("     rejected\n");
+        TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1])>fMaxd0z0) { 
+       AliDebug(1,"     rejected");
        delete track; continue; 
       }
     }
     
-
-    // select tracks with d0rphi < maxd0rphi
-    if(optImpParCut>0 && TMath::Abs(d0z0[0])>maxd0rphi) { 
-      if(fDebug) printf("     rejected\n");
-      delete track; continue; 
-    }
-
+    // track passed all selections
     fTrkArraySel.AddLast(track);
     fIdSel[nTrksSel] = idOrig[i];
     nTrksSel++; 
-  }
+  } // end loop on tracks
 
   delete initVertex;
 
   return nTrksSel;
 } 
+//----------------------------------------------------------------------------
+Bool_t AliVertexerTracks::PropagateTrackTo(AliExternalTrackParam *track,
+                                          Double_t xToGo) {
+  //----------------------------------------------------------------
+  // COPIED from AliTracker
+  //
+  // Propagates the track to the plane X=xk (cm).
+  //
+  //  Origin: Marian Ivanov,  Marian.Ivanov@cern.ch
+  //----------------------------------------------------------------
+
+  const Double_t kEpsilon = 0.00001;
+  Double_t xpos = track->GetX();
+  Double_t dir = (xpos<xToGo) ? 1. : -1.;
+  Double_t maxStep = 5;
+  Double_t maxSnp = 0.8;
+  //
+  while ( (xToGo-xpos)*dir > kEpsilon){
+    Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
+    Double_t x    = xpos+step;
+    Double_t xyz0[3],xyz1[3];
+    track->GetXYZ(xyz0);   //starting global position
+
+    if(!track->GetXYZAt(x,GetFieldkG(),xyz1)) return kFALSE;   // no prolongation
+    xyz1[2]+=kEpsilon; // waiting for bug correction in geo
+
+    if(TMath::Abs(track->GetSnpAt(x,GetFieldkG())) >= maxSnp) return kFALSE;
+    if(!track->PropagateTo(x,GetFieldkG()))  return kFALSE;
+
+    if(TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
+    track->GetXYZ(xyz0);   // global position
+    Double_t alphan = TMath::ATan2(xyz0[1], xyz0[0]); 
+    //
+    Double_t ca=TMath::Cos(alphan-track->GetAlpha()), 
+      sa=TMath::Sin(alphan-track->GetAlpha());
+    Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf));
+    Double_t sinNew =  sf*ca - cf*sa;
+    if(TMath::Abs(sinNew) >= maxSnp) return kFALSE;
+    if(!track->Rotate(alphan)) return kFALSE;
+    xpos = track->GetX();
+  }
+  return kTRUE;
+}
 //---------------------------------------------------------------------------
 AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
                                                        TObjArray *trkArray,
@@ -602,11 +708,11 @@ AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
 //
 
   if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraint")) {
-    printf("ERROR: primary vertex has no constraint: cannot remove tracks\n");
+    printf("ERROR: primary vertex has no constraint: cannot remove tracks");
     return 0x0;
   }
   if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
-    printf("WARNING: result of tracks' removal will be only approximately correct\n");
+    printf("WARNING: result of tracks' removal will be only approximately correct");
 
   TMatrixD rv(3,1);
   rv(0,0) = inVtx->GetXv();
@@ -633,7 +739,7 @@ AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
   for(Int_t i=0;i<ntrks;i++) {
     track = (AliExternalTrackParam*)trkArray->At(i);
     if(!inVtx->UsesTrack(id[i])) {
-      printf("track %d was not used in vertex fit\n",id[i]);
+      printf("track %d was not used in vertex fit",id[i]);
       continue;
     }
     Double_t alpha = track->GetAlpha();
@@ -663,7 +769,7 @@ AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
 
     nUsedTrks--;
     if(nUsedTrks<2) {
-      printf("Trying to remove too many tracks!\n");
+      printf("Trying to remove too many tracks!");
       return 0x0;
     }
   }
@@ -707,18 +813,109 @@ AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
   outVtx->SetIndices(nIndices,outindices);
   delete [] outindices;
 
-  if(fDebug) {
-    printf("Vertex before removing tracks:\n");
+  /*
+    printf("Vertex before removing tracks:");
     inVtx->PrintStatus();
     inVtx->PrintIndices();
-    printf("Vertex after removing tracks:\n");
+    printf("Vertex after removing tracks:");
     outVtx->PrintStatus();
     outVtx->PrintIndices();
-  }
+  */
 
   return outVtx;
 }
 //---------------------------------------------------------------------------
+void AliVertexerTracks::SetCuts(Double_t *cuts) 
+{
+//
+//  Cut values
+//
+  SetDCAcut(cuts[0]);
+  SetDCAcutIter0(cuts[1]);
+  SetMaxd0z0(cuts[2]);
+  if(fMode==0 && cuts[3]<0) SetITSrefitNotRequired();
+  SetMinClusters((Int_t)(TMath::Abs(cuts[3])));
+  SetMinTracks((Int_t)(cuts[4]));
+  SetNSigmad0(cuts[5]);
+  SetMinDetFitter(cuts[6]);
+  SetMaxTgl(cuts[7]);
+  SetFiducialRZ(cuts[8],cuts[9]);
+  fAlgo=(Int_t)(cuts[10]);
+  fAlgoIter0=(Int_t)(cuts[11]);
+
+  return;
+}
+//---------------------------------------------------------------------------
+void AliVertexerTracks::SetITSMode(Double_t dcacut,
+                                  Double_t dcacutIter0,
+                                  Double_t maxd0z0,
+                                  Int_t minCls,
+                                  Int_t mintrks,
+                                  Double_t nsigma,
+                                  Double_t mindetfitter,
+                                  Double_t maxtgl,
+                                  Double_t fidR,
+                                  Double_t fidZ,
+                                  Int_t finderAlgo,
+                                  Int_t finderAlgoIter0)
+{
+//
+//  Cut values for ITS mode
+//
+  fMode = 0;
+  if(minCls>0) {
+    SetITSrefitRequired();
+  } else {
+    SetITSrefitNotRequired();
+  }
+  SetDCAcut(dcacut);
+  SetDCAcutIter0(dcacutIter0);
+  SetMaxd0z0(maxd0z0);
+  SetMinClusters(TMath::Abs(minCls));
+  SetMinTracks(mintrks);
+  SetNSigmad0(nsigma);
+  SetMinDetFitter(mindetfitter);
+  SetMaxTgl(maxtgl);
+  SetFiducialRZ(fidR,fidZ);
+  fAlgo=finderAlgo;
+  fAlgoIter0=finderAlgoIter0;
+
+  return; 
+}
+//---------------------------------------------------------------------------
+void AliVertexerTracks::SetTPCMode(Double_t dcacut,
+                                  Double_t dcacutIter0,
+                                  Double_t maxd0z0,
+                                  Int_t minCls,
+                                  Int_t mintrks,
+                                  Double_t nsigma,
+                                  Double_t mindetfitter,
+                                  Double_t maxtgl,
+                                  Double_t fidR,
+                                  Double_t fidZ,
+                                  Int_t finderAlgo,
+                                  Int_t finderAlgoIter0) 
+{
+//
+//  Cut values for TPC mode
+//
+  fMode = 1;
+  SetITSrefitNotRequired();
+  SetDCAcut(dcacut);
+  SetDCAcutIter0(dcacutIter0);
+  SetMaxd0z0(maxd0z0);
+  SetMinClusters(minCls);
+  SetMinTracks(mintrks);
+  SetNSigmad0(nsigma);
+  SetMinDetFitter(mindetfitter);
+  SetMaxTgl(maxtgl);
+  SetFiducialRZ(fidR,fidZ);
+  fAlgo=finderAlgo;
+  fAlgoIter0=finderAlgoIter0;
+
+  return; 
+}
+//---------------------------------------------------------------------------
 void AliVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped) 
 {
 //
@@ -744,23 +941,22 @@ void  AliVertexerTracks::SetVtxStart(AliESDVertex *vtx)
 void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
 {
   AliExternalTrackParam *track1;
-  Double_t field=GetFieldkG();
   const Int_t knacc = (Int_t)fTrkArraySel.GetEntriesFast();
-  TClonesArray *linarray = new TClonesArray("AliStrLine",1000);
-  TClonesArray &lines = *linarray;
+  static TClonesArray linarray("AliStrLine",knacc);
   for(Int_t i=0; i<knacc; i++){
     track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
     Double_t alpha=track1->GetAlpha();
     Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
     Double_t pos[3],dir[3],sigmasq[3]; 
-    track1->GetXYZAt(mindist,field,pos);
-    track1->GetPxPyPzAt(mindist,field,dir);
+    track1->GetXYZAt(mindist,GetFieldkG(),pos);
+    track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
     sigmasq[0]=TMath::Sin(alpha)*TMath::Sin(alpha)*track1->GetSigmaY2();
     sigmasq[1]=TMath::Cos(alpha)*TMath::Cos(alpha)*track1->GetSigmaY2();
     sigmasq[2]=track1->GetSigmaZ2();
     TMatrixD ri(3,1);
     TMatrixD wWi(3,3);
-    if(!TrackToPoint(track1,ri,wWi)) continue;
+    //if(!TrackToPoint(track1,ri,wWi)) {printf("WARNING\n");continue;}
+    if(!TrackToPoint(track1,ri,wWi)) {optUseWeights=kFALSE;printf("WARNING\n");}
     Double_t wmat[9];
     Int_t iel=0;
     for(Int_t ia=0;ia<3;ia++){
@@ -769,17 +965,15 @@ void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
        iel++;
       }    
     }
-    new(lines[i]) AliStrLine(pos,sigmasq,wmat,dir);     
+    new(linarray[i]) AliStrLine(pos,sigmasq,wmat,dir);
   }
-  fVert=TrackletVertexFinder(linarray,optUseWeights);
-  linarray->Delete();
-  delete linarray;
+  fVert=TrackletVertexFinder(&linarray,optUseWeights);
+  linarray.Clear("C");
 }
 //---------------------------------------------------------------------------
 AliESDVertex AliVertexerTracks::TrackletVertexFinder(TClonesArray *lines, Int_t optUseWeights)
 {
   // Calculate the point at minimum distance to prepared tracks 
-  
   const Int_t knacc = (Int_t)lines->GetEntriesFast();
   Double_t initPos[3]={0.,0.,0.};
 
@@ -797,10 +991,10 @@ AliESDVertex AliVertexerTracks::TrackletVertexFinder(TClonesArray *lines, Int_t
   }
 
   for(Int_t i=0; i<knacc; i++){
-    AliStrLine* line1 = (AliStrLine*)lines->At(i);
-    if (!line1) continue;
+    AliStrLine* line1 = (AliStrLine*)lines->At(i); 
     Double_t p0[3],cd[3],sigmasq[3];
     Double_t wmat[9];
+    if(!line1) printf("ERROR %d %d\n",i,knacc);
     line1->GetP0(p0);
     line1->GetCd(cd);
     line1->GetSigma2P0(sigmasq);
@@ -868,7 +1062,7 @@ AliESDVertex AliVertexerTracks::TrackletVertexFinder(TClonesArray *lines, Int_t
       sigma+=GetStrLinMinDist(p0,p1,initPos);
     }
 
-    sigma=TMath::Sqrt(sigma);
+    if(sigma>0.) {sigma=TMath::Sqrt(sigma);}else{sigma=999;}
   }else{
     sigma=999;
   }
@@ -916,9 +1110,9 @@ Bool_t AliVertexerTracks::TrackToPoint(AliExternalTrackParam *t,
     uUi(0,1) = t->GetSigmaZY();
     uUi(1,0) = t->GetSigmaZY();
     uUi(1,1) = t->GetSigmaZ2();
-    //printf(" Ui :\n");
-    //printf(" %f   %f\n",uUi(0,0),uUi(0,1));
-    //printf(" %f   %f\n",uUi(1,0),uUi(1,1));
+    //printf(" Ui :");
+    //printf(" %f   %f",uUi(0,0),uUi(0,1));
+    //printf(" %f   %f",uUi(1,0),uUi(1,1));
 
     if(uUi.Determinant() <= 0.) return kFALSE;
     TMatrixD uUiInv(TMatrixD::kInverted,uUi);
@@ -962,7 +1156,7 @@ Bool_t AliVertexerTracks::TrackToPoint(AliExternalTrackParam *t,
     // covariance matrix of local (x,y,z) - inverted
     TMatrixD uUi(3,3);
     uUi(0,0) = covfVertalongt * fnSigmaForUi00 * fnSigmaForUi00;
-    if(fDebug) printf("=====> sqrtUi00 cm  %f\n",TMath::Sqrt(uUi(0,0)));
+    AliDebug(1,Form("=====> sqrtUi00 cm  %f",TMath::Sqrt(uUi(0,0))));
     uUi(0,1) = 0.;
     uUi(0,2) = 0.;
     uUi(1,0) = 0.;
@@ -994,7 +1188,7 @@ void AliVertexerTracks::TooFewTracks()
 // When the number of tracks is < fMinTracks,
 // deal with vertices not found and prepare to exit
 //
-  if(fDebug) printf("TooFewTracks\n");
+  AliDebug(1,"TooFewTracks");
 
   Double_t pos[3],err[3];
   pos[0] = fNominalPos[0];
@@ -1004,7 +1198,7 @@ void AliVertexerTracks::TooFewTracks()
   pos[2] = fNominalPos[2];
   err[2] = TMath::Sqrt(fNominalCov[5]);
   Int_t    ncontr = (err[0]>1. ? -1 : -3);
-  fCurrentVertex = 0;
+  if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
   fCurrentVertex = new AliESDVertex(pos,err);
   fCurrentVertex->SetNContributors(ncontr);
 
@@ -1039,28 +1233,27 @@ void AliVertexerTracks::VertexFinder(Int_t optUseWeights)
   AliExternalTrackParam *track2;
   Double_t pos[3],dir[3]; 
   Double_t alpha,mindist;
-  Double_t field=GetFieldkG();
 
   for(Int_t i=0; i<nacc; i++){
     track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
     alpha=track1->GetAlpha();
     mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
-    track1->GetXYZAt(mindist,field,pos);
-    track1->GetPxPyPzAt(mindist,field,dir);
+    track1->GetXYZAt(mindist,GetFieldkG(),pos);
+    track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
     AliStrLine *line1 = new AliStrLine(pos,dir); 
 
    //    AliStrLine *line1 = new AliStrLine();
-   //    track1->ApproximateHelixWithLine(mindist,field,line1);
+   //    track1->ApproximateHelixWithLine(mindist,GetFieldkG(),line1);
    
     for(Int_t j=i+1; j<nacc; j++){
       track2 = (AliExternalTrackParam*)fTrkArraySel.At(j);
       alpha=track2->GetAlpha();
       mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
-      track2->GetXYZAt(mindist,field,pos);
-      track2->GetPxPyPzAt(mindist,field,dir);
+      track2->GetXYZAt(mindist,GetFieldkG(),pos);
+      track2->GetPxPyPzAt(mindist,GetFieldkG(),dir);
       AliStrLine *line2 = new AliStrLine(pos,dir); 
     //      AliStrLine *line2 = new AliStrLine();
-    //  track2->ApproximateHelixWithLine(mindist,field,line2);
+    //  track2->ApproximateHelixWithLine(mindist,GetFieldkG(),line2);
       Double_t distCA=line2->GetDCA(line1);
       //printf("%d   %d   %f\n",i,j,distCA);
        if(fDCAcut<=0 || (fDCAcut>0&&distCA<fDCAcut)){
@@ -1140,13 +1333,11 @@ void AliVertexerTracks::VertexFitter()
 
   Int_t nTrksSel = (Int_t)fTrkArraySel.GetEntries();
   if(nTrksSel==1) useConstraint=kTRUE;
-  if(fDebug) { 
-    printf("--- VertexFitter(): start\n");
-    printf(" Number of tracks in array: %d\n",nTrksSel);
-    printf(" Minimum # tracks required in fit: %d\n",fMinTracks);
-    printf(" Vertex position after finder: %f,%f,%f\n",initPos[0],initPos[1],initPos[2]);
-    if(useConstraint) printf(" This vertex will be used in fit: (%f+-%f,%f+-%f)\n",fNominalPos[0],TMath::Sqrt(fNominalCov[0]),fNominalPos[1],TMath::Sqrt(fNominalCov[2])); 
-  }
+  AliDebug(1,"--- VertexFitter(): start");
+  AliDebug(1,Form(" Number of tracks in array: %d\n",nTrksSel));
+  AliDebug(1,Form(" Minimum # tracks required in fit: %d\n",fMinTracks));
+  AliDebug(1,Form(" Vertex position after finder: %f,%f,%f\n",initPos[0],initPos[1],initPos[2]));
+  if(useConstraint) AliDebug(1,Form(" This vertex will be used in fit: (%f+-%f,%f+-%f)\n",fNominalPos[0],TMath::Sqrt(fNominalCov[0]),fNominalPos[1],TMath::Sqrt(fNominalCov[2]))); 
 
   // special treatment for few-tracks fits (e.g. secondary vertices)
   Bool_t uUi3by3 = kFALSE; if(nTrksSel<5 && !useConstraint) uUi3by3 = kTRUE;
@@ -1186,7 +1377,7 @@ void AliVertexerTracks::VertexFitter()
   // 1st - estimate of vtx using all tracks
   // 2nd - estimate of global chi2
   for(step=0; step<2; step++) {
-    if(fDebug) printf(" step = %d\n",step);
+    AliDebug(1,Form(" step = %d\n",step));
     chi2 = 0.;
     nTrksUsed = 0;
 
@@ -1254,9 +1445,8 @@ void AliVertexerTracks::VertexFitter()
     }
 
     Double_t determinant = sumWi.Determinant();
-    //cerr<<" determinant: "<<determinant<<endl;
-    if(determinant < 100.)  { 
-      printf("det(V) = 0\n");       
+    if(determinant < fMinDetFitter)  { 
+      AliDebug(1,Form("det(V) = %f (<%f)\n",determinant,fMinDetFitter));       
       failed=1;
       continue;
     }
@@ -1290,13 +1480,13 @@ void AliVertexerTracks::VertexFitter()
   // for correct chi2/ndf, count constraint as additional "track"
   if(fConstraint) nTrksUsed++;
   // store data in the vertex object
+  if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
   fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);
 
-  if(fDebug) {
-    printf(" Vertex after fit:\n");
-    fCurrentVertex->PrintStatus();
-    printf("--- VertexFitter(): finish\n");
-  }
+  AliDebug(1," Vertex after fit:");
+  AliDebug(1,Form("xyz: %f %f %f; nc %d",fCurrentVertex->GetXv(),fCurrentVertex->GetYv(),fCurrentVertex->GetZv(),fCurrentVertex->GetNContributors()));
+  AliDebug(1,"--- VertexFitter(): finish\n");
 
   return;
 }
@@ -1309,7 +1499,7 @@ AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray,
 //
 // Return vertex from tracks (AliExternalTrackParam) in array
 //
-
+  fCurrentVertex = 0;
   SetConstraintOff();
 
   // get tracks and propagate them to initial vertex position
@@ -1329,7 +1519,7 @@ AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray,
     case 5: VertexFinder(0);              break;
     default: printf("Wrong algorithm\n"); break;  
   }
-  if(fDebug) printf(" Vertex finding completed\n");
+  AliDebug(1," Vertex finding completed\n");
 
   // vertex fitter
   if(optUseFitter) {
@@ -1357,7 +1547,7 @@ AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray,
       if(optPropagate && optUseFitter) {
        if(TMath::Sqrt(fCurrentVertex->GetXv()*fCurrentVertex->GetXv()+fCurrentVertex->GetYv()*fCurrentVertex->GetYv())<3.) {
          t->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
-         if(fDebug) printf("Track %d propagated to found vertex\n",jj);
+         AliDebug(1,Form("Track %d propagated to found vertex",jj));
        } else {
          AliWarning("Found vertex outside beam pipe!");
        }
@@ -1398,4 +1588,3 @@ AliESDVertex* AliVertexerTracks::VertexForSelectedESDTracks(TObjArray *trkArray,
   return fCurrentVertex;
 }
 //--------------------------------------------------------------------------
-