RemovetracksFromVertex() method and other new possibilities (Andrea, Francesco)

diff --git a/STEER/AliVertexerTracks.cxx b/STEER/AliVertexerTracks.cxx
index ebde7f38c2887008ee10317ca5238c55f20b3054..41ecac10d35465300b6088adbe9094862ad075e5 100644 (file)
--- a/STEER/AliVertexerTracks.cxx
+++ b/STEER/AliVertexerTracks.cxx
@@ -42,6 +42,8 @@ AliVertexerTracks::AliVertexerTracks():
 TObject(),
 fVert(),
 fCurrentVertex(0),
 TObject(),
 fVert(),
 fCurrentVertex(0),

+fConstraint(kFALSE),
+fOnlyFitter(kFALSE),

 fMinTracks(1),
 fMinITSClusters(5),
 fTrkArray(),
 fMinTracks(1),
 fMinITSClusters(5),
 fTrkArray(),


@@ -50,6 +52,7 @@ fNTrksToSkip(0),
 fDCAcut(0),
 fAlgo(1),
 fNSigma(3),
 fDCAcut(0),
 fAlgo(1),
 fNSigma(3),

+fMaxd0z0(0.5),

 fITSin(kTRUE),
 fITSrefit(kTRUE),
 fDebug(0)
 fITSin(kTRUE),
 fITSrefit(kTRUE),
 fDebug(0)


@@ -62,12 +65,15 @@ fDebug(0)
   SetMinTracks();
   SetMinITSClusters();
   SetNSigmad0(); 
   SetMinTracks();
   SetMinITSClusters();
   SetNSigmad0(); 

+  SetMaxd0z0(); 

 }
 //-----------------------------------------------------------------------------
 AliVertexerTracks::AliVertexerTracks(Double_t xStart, Double_t yStart):
 TObject(),
 fVert(),
 fCurrentVertex(0),
 }
 //-----------------------------------------------------------------------------
 AliVertexerTracks::AliVertexerTracks(Double_t xStart, Double_t yStart):
 TObject(),
 fVert(),
 fCurrentVertex(0),

+fConstraint(kFALSE),
+fOnlyFitter(kFALSE),

 fMinTracks(1),
 fMinITSClusters(5),
 fTrkArray(),
 fMinTracks(1),
 fMinITSClusters(5),
 fTrkArray(),


@@ -76,6 +82,7 @@ fNTrksToSkip(0),
 fDCAcut(0),
 fAlgo(1),
 fNSigma(3),
 fDCAcut(0),
 fAlgo(1),
 fNSigma(3),

+fMaxd0z0(0.5),

 fITSin(kTRUE),
 fITSrefit(kTRUE),
 fDebug(0)
 fITSin(kTRUE),
 fITSrefit(kTRUE),
 fDebug(0)


@@ -88,9 +95,11 @@ fDebug(0)
   SetMinTracks();
   SetMinITSClusters();
   SetNSigmad0(); 
   SetMinTracks();
   SetMinITSClusters();
   SetNSigmad0(); 

+  SetMaxd0z0(); 

 }
 //-----------------------------------------------------------------------------
 }
 //-----------------------------------------------------------------------------

-AliVertexerTracks::~AliVertexerTracks() {
+AliVertexerTracks::~AliVertexerTracks() 
+{

   // Default Destructor
   // The objects pointed by the following pointer are not owned
   // by this class and are not deleted
   // Default Destructor
   // The objects pointed by the following pointer are not owned
   // by this class and are not deleted


@@ -103,12 +112,14 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliESD *esdEvent)
 //
 // Primary vertex for current ESD event
 // (Two iterations: 
 //
 // Primary vertex for current ESD event
 // (Two iterations: 

-//  1st with 5*fNSigma*sigma(pt) cut w.r.t. to initial vertex; 
-//  2nd with fNSigma*sigma(pt) cut w.r.t. to vertex found in 1st iteration) 
+//  1st with 5*fNSigma*sigma cut w.r.t. to initial vertex
+//      + cut on sqrt(d0d0+z0z0) if fConstraint=kFALSE  
+//  2nd with fNSigma*sigma cut w.r.t. to vertex found in 1st iteration) 

 // All ESD tracks with inside the beam pipe are then propagated to found vertex
 //
   fCurrentVertex = 0;
 
 // All ESD tracks with inside the beam pipe are then propagated to found vertex
 //
   fCurrentVertex = 0;
 

+  // read tracks from ESD

   Int_t nTrksTot = (Int_t)esdEvent->GetNumberOfTracks();
   TTree *trkTree = new TTree("TreeT","tracks");
   AliESDtrack *esdTrack = 0;
   Int_t nTrksTot = (Int_t)esdEvent->GetNumberOfTracks();
   TTree *trkTree = new TTree("TreeT","tracks");
   AliESDtrack *esdTrack = 0;


@@ -116,16 +127,16 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliESD *esdEvent)
 
   Bool_t   skipThis;
   for(Int_t i=0; i<nTrksTot; i++) {
 
   Bool_t   skipThis;
   for(Int_t i=0; i<nTrksTot; i++) {

+    AliESDtrack *et = esdEvent->GetTrack(i);
+    esdTrack = new AliESDtrack(*et);

     // check tracks to skip
     skipThis = kFALSE;
     for(Int_t j=0; j<fNTrksToSkip; j++) { 
     // check tracks to skip
     skipThis = kFALSE;
     for(Int_t j=0; j<fNTrksToSkip; j++) { 

-      if(i==fTrksToSkip[j]) {
+      if(et->GetID()==fTrksToSkip[j]) {

        if(fDebug) printf("skipping track: %d\n",i);
        skipThis = kTRUE;
       }
     }
        if(fDebug) printf("skipping track: %d\n",i);
        skipThis = kTRUE;
       }
     }

-    AliESDtrack *et = esdEvent->GetTrack(i);
-    esdTrack = new AliESDtrack(*et);

     if(skipThis) {delete esdTrack;continue;}
     if(fITSin) {
       if(!(esdTrack->GetStatus()&AliESDtrack::kITSin)) {delete esdTrack;continue;}
     if(skipThis) {delete esdTrack;continue;}
     if(fITSin) {
       if(!(esdTrack->GetStatus()&AliESDtrack::kITSin)) {delete esdTrack;continue;}


@@ -136,102 +147,96 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliESD *esdEvent)
     trkTree->Fill();
     delete esdTrack;
   }
     trkTree->Fill();
     delete esdTrack;
   }

- 
-
-  // ITERATION 1
-  // propagate tracks to initVertex
-  // preselect them  (reject for |d0|>5*fNSigma*sigma w.r.t. initVertex)
-  Int_t nTrksPrep;
-  nTrksPrep = PrepareTracks(*trkTree,1);
-  if(fDebug) printf(" tracks prepared - step 1: %d\n",nTrksPrep);
-  if(nTrksPrep < fMinTracks) {
-    if(fDebug) printf("TooFewTracks\n");
-    TooFewTracks(esdEvent);
-    if(fDebug) fCurrentVertex->PrintStatus();
+  
+  // If fConstraint=kFALSE
+  // run VertexFinder(1) to get rough estimate of initVertex (x,y)
+  if(!fConstraint) {
+    // fill fTrkArray, for VertexFinder()
+    if(!fTrkArray.IsEmpty()) fTrkArray.Delete();
+    PrepareTracks(*trkTree,0);
+    Double_t cutsave = fDCAcut;  fDCAcut = 0.1; // 1 mm
+    VertexFinder(1); // using weights, cutting dca < fDCAcut
+    fDCAcut = cutsave;

     fTrkArray.Delete();
     fTrkArray.Delete();

-    delete trkTree;
-    return fCurrentVertex; 
-  }
-
-  if(nTrksPrep==1){
-    if(fDebug) printf("Just one track\n");
-    OneTrackVertFinder();
-  }else{
-  // vertex finder
-    switch (fAlgo) {
-    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;  
+    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]);
+    } else {
+      fNominalPos[0] = 0.;
+      fNominalPos[1] = 0.;
+      fNominalPos[2] = 0.;
+      if(fDebug) printf("No mean vertex and VertexFinder failed\n");

     }
   }
     }
   }

-  if(fDebug) printf(" vertex finding completed\n");
-
-  // vertex fitter
-  VertexFitter(kTRUE);
-  if(fDebug) printf(" vertex fit completed\n");
-  fTrkArray.Delete();
+  
+  // TWO ITERATIONS:
+  //
+  // ITERATION 1
+  // propagate tracks to fNominalPos vertex
+  // preselect them:
+  // if(constraint) reject for |d0|>5*fNSigma*sigma w.r.t. fNominal... vertex
+  // else  reject for |d0|\oplus|z0| > 5 mm w.r.t. fNominal... vertex

   // ITERATION 2
   // propagate tracks to best between initVertex and fCurrentVertex
   // preselect tracks (reject for |d0|>fNSigma*sigma w.r.t. best 
   //                   between initVertex and fCurrentVertex) 
   // ITERATION 2
   // propagate tracks to best between initVertex and fCurrentVertex
   // preselect tracks (reject for |d0|>fNSigma*sigma w.r.t. best 
   //                   between initVertex and fCurrentVertex) 

-  nTrksPrep = PrepareTracks(*trkTree,2);
-  delete trkTree;
-  if(fDebug) printf(" tracks prepared - step 2: %d\n",nTrksPrep);
-  if(nTrksPrep < fMinTracks) {
-    if(fDebug) printf("TooFewTracks\n");
-    TooFewTracks(esdEvent);
-    if(fDebug) fCurrentVertex->PrintStatus();
-    fTrkArray.Delete();
-    return fCurrentVertex; 
-  }
+  for(Int_t iter=0; iter<2; iter++) {
+    if(fOnlyFitter && iter==0) continue; 
+    Int_t nTrksPrep = PrepareTracks(*trkTree,iter+1);
+    if(fDebug) printf(" tracks prepared - iteration %d: %d\n",iter+1,nTrksPrep);
+    if(nTrksPrep < fMinTracks) {
+      if(fDebug) printf("TooFewTracks\n");
+      TooFewTracks(esdEvent);
+      if(fDebug) fCurrentVertex->PrintStatus();
+      fTrkArray.Delete();
+      delete trkTree;
+      return fCurrentVertex; 
+    }

 
 

-  // vertex finder
-  if(nTrksPrep==1){
-    if(fDebug) printf("Just one track\n");
-    OneTrackVertFinder();
-  }else{
-    switch (fAlgo) {
-    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;  
+    // vertex finder
+    if(!fOnlyFitter) {
+      if(nTrksPrep==1){
+       if(fDebug) printf("Just one track\n");
+       OneTrackVertFinder();
+      }else{
+       switch (fAlgo) {
+        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;  
+       }
+      }
+      if(fDebug) printf(" Vertex finding completed\n");

     }
     }

-  }
-  if(fDebug) printf(" vertex finding completed\n");

 
 

-  // fitter
-  VertexFitter(kTRUE);
-  if(fDebug) printf(" vertex fit completed\n");
+    // vertex fitter
+    VertexFitter(fConstraint);
+    if(fDebug) printf(" Vertex fit completed\n");
+    if(iter==0) fTrkArray.Delete();
+  } // end loop on the two iterations

 
 
 
 

-  // take true pos from ESD
+  // take true pos from SPD vertex in ESD and write it in tracks' vertex

   Double_t tp[3];
   esdEvent->GetVertex()->GetTruePos(tp);
   fCurrentVertex->SetTruePos(tp);
   Double_t tp[3];
   esdEvent->GetVertex()->GetTruePos(tp);
   fCurrentVertex->SetTruePos(tp);

-  if(fNominalCov[0]>1.) {
-    fCurrentVertex->SetTitle("VertexerTracksNoConstraint");
-  } else {
-    fCurrentVertex->SetTitle("VertexerTracksWithConstraint");
-  }

 
 

-  
-  // propagate tracks to found vertex
-  if(TMath::Sqrt(fCurrentVertex->GetXv()*fCurrentVertex->GetXv()+fCurrentVertex->GetYv()*fCurrentVertex->GetYv())<3.) {
-    for(Int_t ii=0; ii<nTrksTot; ii++) {
-      AliESDtrack *et = esdEvent->GetTrack(ii);
-      if(!(et->GetStatus()&AliESDtrack::kITSin)) continue;
-      if(et->GetX()>3.) continue;
-      et->RelateToVertex(fCurrentVertex,GetField(),100.);
+  if(fConstraint) {
+    if(fOnlyFitter) {
+      fCurrentVertex->SetTitle("VertexerTracksWithConstraintOnlyFitter");
+    } else {
+      fCurrentVertex->SetTitle("VertexerTracksWithConstraint");

     }
   } else {
     }
   } else {

-    AliWarning("Found vertex outside beam pipe!");
+    fCurrentVertex->SetTitle("VertexerTracksNoConstraint");

   }
 
   }
 

+    

   // set indices of used tracks
   UShort_t *indices = 0;
   AliESDtrack *ett = 0;
   // set indices of used tracks
   UShort_t *indices = 0;
   AliESDtrack *ett = 0;


@@ -244,6 +249,9 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliESD *esdEvent)
     fCurrentVertex->SetIndices(fCurrentVertex->GetNContributors(),indices);
   }
   delete [] indices;
     fCurrentVertex->SetIndices(fCurrentVertex->GetNContributors(),indices);
   }
   delete [] indices;

+
+  delete trkTree;
+  

   fTrkArray.Delete();
 
   if(fTrksToSkip) delete [] fTrksToSkip;
   fTrkArray.Delete();
 
   if(fTrksToSkip) delete [] fTrksToSkip;


@@ -251,18 +259,20 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliESD *esdEvent)
 
   if(fDebug) fCurrentVertex->PrintStatus();
   if(fDebug) fCurrentVertex->PrintIndices();
 
   if(fDebug) fCurrentVertex->PrintStatus();
   if(fDebug) fCurrentVertex->PrintIndices();

- 
+
+  

   return fCurrentVertex;
 }
 //------------------------------------------------------------------------
   return fCurrentVertex;
 }
 //------------------------------------------------------------------------

-Double_t AliVertexerTracks::GetDeterminant3X3(Double_t matr[][3]){
+Double_t AliVertexerTracks::GetDeterminant3X3(Double_t matr[][3])
+{

   //
   Double_t det=matr[0][0]*matr[1][1]*matr[2][2]-matr[0][0]*matr[1][2]*matr[2][1]-matr[0][1]*matr[1][0]*matr[2][2]+matr[0][1]*matr[1][2]*matr[2][0]+matr[0][2]*matr[1][0]*matr[2][1]-matr[0][2]*matr[1][1]*matr[2][0];
  return det;
 }
 //-------------------------------------------------------------------------
   //
   Double_t det=matr[0][0]*matr[1][1]*matr[2][2]-matr[0][0]*matr[1][2]*matr[2][1]-matr[0][1]*matr[1][0]*matr[2][2]+matr[0][1]*matr[1][2]*matr[2][0]+matr[0][2]*matr[1][0]*matr[2][1]-matr[0][2]*matr[1][1]*matr[2][0];
  return det;
 }
 //-------------------------------------------------------------------------

-void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t (*m)[3],Double_t *d){
-
+void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t (*m)[3],Double_t *d)
+{

   //
   Double_t x12=p0[0]-p1[0];
   Double_t y12=p0[1]-p1[1];
   //
   Double_t x12=p0[0]-p1[0];
   Double_t y12=p0[1]-p1[1];


@@ -283,7 +293,8 @@ void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t
 
 }
 //--------------------------------------------------------------------------  
 
 }
 //--------------------------------------------------------------------------  

-void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t *sigmasq,Double_t (*m)[3],Double_t *d){
+void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t *sigmasq,Double_t (*m)[3],Double_t *d)
+{

   //
   Double_t x12=p1[0]-p0[0];
   Double_t y12=p1[1]-p0[1];
   //
   Double_t x12=p1[0]-p0[0];
   Double_t y12=p1[1]-p0[1];


@@ -325,7 +336,8 @@ void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t
 
   }
 //--------------------------------------------------------------------------   
 
   }
 //--------------------------------------------------------------------------   

-Double_t AliVertexerTracks::GetStrLinMinDist(Double_t *p0,Double_t *p1,Double_t *x0){
+Double_t AliVertexerTracks::GetStrLinMinDist(Double_t *p0,Double_t *p1,Double_t *x0)
+{

   //
   Double_t x12=p0[0]-p1[0];
   Double_t y12=p0[1]-p1[1];
   //
   Double_t x12=p0[0]-p1[0];
   Double_t y12=p0[1]-p1[1];


@@ -336,7 +348,8 @@ Double_t AliVertexerTracks::GetStrLinMinDist(Double_t *p0,Double_t *p1,Double_t
   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);
 }
 //---------------------------------------------------------------------------

-void AliVertexerTracks::OneTrackVertFinder() {
+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",fTrkArray.GetEntries());
   AliESDtrack *track1;
   // find vertex for events with 1 track, using DCA to nominal beam axis
   if(fDebug) printf("Number of prepared tracks =%d - Call OneTrackVertFinder",fTrkArray.GetEntries());
   AliESDtrack *track1;


@@ -366,8 +379,8 @@ void AliVertexerTracks::OneTrackVertFinder() {
   fVert.SetNContributors(nContrib);  
 }
 //---------------------------------------------------------------------------
   fVert.SetNContributors(nContrib);  
 }
 //---------------------------------------------------------------------------

-void AliVertexerTracks::HelixVertexFinder() {
-
+void AliVertexerTracks::HelixVertexFinder() 
+{

   // Get estimate of vertex position in (x,y) from tracks DCA
 
 
   // Get estimate of vertex position in (x,y) from tracks DCA
 
 


@@ -397,7 +410,6 @@ void AliVertexerTracks::HelixVertexFinder() {
       track2 = (AliESDtrack*)fTrkArray.At(j);
 
       distCA=track2->PropagateToDCA(track1,field);
       track2 = (AliESDtrack*)fTrkArray.At(j);
 
       distCA=track2->PropagateToDCA(track1,field);

-

       if(fDCAcut<=0 ||(fDCAcut>0&&distCA<fDCAcut)){
        track1->GetExternalParameters(x,par);
        alpha=track1->GetAlpha();
       if(fDCAcut<=0 ||(fDCAcut>0&&distCA<fDCAcut)){
        track1->GetExternalParameters(x,par);
        alpha=track1->GetAlpha();


@@ -446,13 +458,16 @@ void AliVertexerTracks::HelixVertexFinder() {
   fVert.SetNContributors(ncombi);
 }
 //----------------------------------------------------------------------------
   fVert.SetNContributors(ncombi);
 }
 //----------------------------------------------------------------------------

-Int_t AliVertexerTracks::PrepareTracks(TTree &trkTree,Int_t OptImpParCut) {
+Int_t AliVertexerTracks::PrepareTracks(TTree &trkTree,Int_t optImpParCut) 
+{

 //
 // Propagate tracks to initial vertex position and store them in a TObjArray
 //
 //
 // Propagate tracks to initial vertex position and store them in a TObjArray
 //

-  Double_t maxd0rphi = 3.;  
+  Double_t maxd0rphi; 
+  Double_t maxd0z0 = fMaxd0z0; // default is 5 mm  

   Int_t    nTrks    = 0;
   Double_t sigmaCurr[3];
   Int_t    nTrks    = 0;
   Double_t sigmaCurr[3];

+  Double_t normdistx,normdisty;

   Float_t  d0z0[2],covd0z0[3]; 
   Double_t sigma;
   Double_t field=GetField();
   Float_t  d0z0[2],covd0z0[3]; 
   Double_t sigma;
   Double_t field=GetField();


@@ -472,25 +487,41 @@ Int_t AliVertexerTracks::PrepareTracks(TTree &trkTree,Int_t OptImpParCut) {
     trkTree.GetEvent(i);
 
     // propagate track to vertex
     trkTree.GetEvent(i);
 
     // propagate track to vertex

-    if(OptImpParCut==1) { // OptImpParCut==1
+    if(optImpParCut<=1 || fOnlyFitter) { // optImpParCut==1 or 0

       track->RelateToVertex(initVertex,field,100.);
       track->RelateToVertex(initVertex,field,100.);

-    } else {              // OptImpParCut==2
+    } else {              // optImpParCut==2

       fCurrentVertex->GetSigmaXYZ(sigmaCurr);
       fCurrentVertex->GetSigmaXYZ(sigmaCurr);

-      if((sigmaCurr[0]+sigmaCurr[1])<(TMath::Sqrt(fNominalCov[0])+TMath::Sqrt(fNominalCov[2]))) {
+      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[1]);
+      if(normdistx < 3. && normdisty < 3. &&
+        (sigmaCurr[0]+sigmaCurr[1])<(TMath::Sqrt(fNominalCov[0])+TMath::Sqrt(fNominalCov[1]))) {

        track->RelateToVertex(fCurrentVertex,field,100.);
       } else {
        track->RelateToVertex(initVertex,field,100.);
       }
     }
 
        track->RelateToVertex(fCurrentVertex,field,100.);
       } else {
        track->RelateToVertex(initVertex,field,100.);
       }
     }
 

-    // select tracks with d0rphi < maxd0rphi

     track->GetImpactParameters(d0z0,covd0z0);
     sigma = TMath::Sqrt(covd0z0[0]);
     maxd0rphi = fNSigma*sigma;
     track->GetImpactParameters(d0z0,covd0z0);
     sigma = TMath::Sqrt(covd0z0[0]);
     maxd0rphi = fNSigma*sigma;

-    if(OptImpParCut==1) maxd0rphi *= 5.;
+    if(optImpParCut==1) maxd0rphi *= 5.;
+
+

 
 

-    if(fDebug) printf("trk %d; lab %d; |d0| = %f;  cut = %f\n",i,track->GetLabel(),TMath::Abs(d0z0[0]),maxd0rphi);
-    if(TMath::Abs(d0z0[0]) > maxd0rphi) { 
+    if(fDebug) printf("trk %d; lab %d; |d0| = %f;  d0 cut = %f; |z0| = %f; |d0|oplus|z0| = %f; d0z0 cut = %f\n",i,track->GetLabel(),TMath::Abs(d0z0[0]),maxd0rphi,TMath::Abs(d0z0[1]),TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1]),maxd0z0);
+
+    // during iteration 1, if fConstraint=kFALSE,
+    // select tracks with d0oplusz0 < maxd0z0
+    if(optImpParCut==1 && !fConstraint && nEntries>=3 && 
+       fVert.GetNContributors()>0) {
+      if(TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1]) > maxd0z0) { 
+       if(fDebug) printf("     rejected\n");
+       delete track; continue; 
+      }
+    }
+
+    // select tracks with d0rphi < maxd0rphi
+    if(optImpParCut>0 && TMath::Abs(d0z0[0]) > maxd0rphi) { 

       if(fDebug) printf("     rejected\n");
       delete track; continue; 
     }
       if(fDebug) printf("     rejected\n");
       delete track; continue; 
     }


@@ -504,27 +535,157 @@ Int_t AliVertexerTracks::PrepareTracks(TTree &trkTree,Int_t OptImpParCut) {
   return nTrks;
 } 
 //---------------------------------------------------------------------------
   return nTrks;
 } 
 //---------------------------------------------------------------------------

-void AliVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped) {
+AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
+                                                       TTree *trksTree,
+                                                       Float_t *diamondxy) 
+{

 //
 //

-// Mark the tracks not ot be used in the vertex finding
+// Removes tracks in trksTree from fit of inVtx

 //
 //

-  fNTrksToSkip = n;
-  fTrksToSkip = new Int_t[n]; 
+
+  if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
+    printf("WARNING: result of tracks' removal will be only approximately correct\n");
+
+  TMatrixD rv(3,1);
+  rv(0,0) = inVtx->GetXv();
+  rv(1,0) = inVtx->GetYv();
+  rv(2,0) = inVtx->GetZv();
+  TMatrixD vV(3,3);
+  Double_t cov[6];
+  inVtx->GetCovMatrix(cov);
+  vV(0,0) = cov[0];
+  vV(0,1) = cov[1]; vV(1,0) = cov[1];
+  vV(1,1) = cov[2];
+  vV(0,2) = cov[3]; vV(2,0) = cov[3];
+  vV(1,2) = cov[4]; vV(2,1) = cov[4]; 
+  vV(2,2) = cov[5];
+
+  TMatrixD sumWi(TMatrixD::kInverted,vV);
+  TMatrixD sumWiri(sumWi,TMatrixD::kMult,rv);
+
+  Int_t nUsedTrks = inVtx->GetNContributors();
+  Double_t chi2 = inVtx->GetChi2();
+
+  AliESDtrack *track = 0;
+  trksTree->SetBranchAddress("tracks",&track);
+  Int_t ntrks = trksTree->GetEntries();
+  for(Int_t i=0;i<ntrks;i++) {
+    trksTree->GetEvent(i);
+    if(!inVtx->UsesTrack(track->GetID())) {
+      printf("track %d was not used in vertex fit\n",track->GetID());
+      continue;
+    }
+    Double_t alpha = track->GetAlpha();
+    Double_t xl = diamondxy[0]*TMath::Cos(alpha)+diamondxy[1]*TMath::Sin(alpha);
+    track->AliExternalTrackParam::PropagateTo(xl,AliTracker::GetBz()); 
+    // vector of track global coordinates
+    TMatrixD ri(3,1);
+    // covariance matrix of ri
+    TMatrixD wWi(3,3);
+    
+    // get space point from track
+    if(!TrackToPoint(track,ri,wWi)) continue;
+
+    TMatrixD wWiri(wWi,TMatrixD::kMult,ri); 
+
+    sumWi -= wWi;
+    sumWiri -= wWiri;
+
+    // track chi2
+    TMatrixD deltar = rv; deltar -= ri;
+    TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
+    Double_t chi2i = deltar(0,0)*wWideltar(0,0)+
+                     deltar(1,0)*wWideltar(1,0)+
+                    deltar(2,0)*wWideltar(2,0);
+    // remove from total chi2
+    chi2 -= chi2i;
+
+    nUsedTrks--;
+    if(nUsedTrks<2) {
+      printf("Trying to remove too many tracks!\n");
+      return 0x0;
+    }
+  }
+
+  TMatrixD rvnew(3,1);
+  TMatrixD vVnew(3,3);
+
+  // new inverted of weights matrix
+  TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
+  vVnew = invsumWi;
+  // new position of primary vertex
+  rvnew.Mult(vVnew,sumWiri);
+
+  Double_t position[3];
+  position[0] = rvnew(0,0);
+  position[1] = rvnew(1,0);
+  position[2] = rvnew(2,0);
+  cov[0] = vVnew(0,0);
+  cov[1] = vVnew(0,1);
+  cov[2] = vVnew(1,1);
+  cov[3] = vVnew(0,2);
+  cov[4] = vVnew(1,2);
+  cov[5] = vVnew(2,2);
+  
+  // store data in the vertex object
+  AliESDVertex *outVtx = new AliESDVertex(position,cov,chi2,nUsedTrks);
+  outVtx->SetTitle(inVtx->GetTitle());
+  Double_t tp[3];
+  inVtx->GetTruePos(tp);
+  outVtx->SetTruePos(tp);
+  UShort_t *inindices = inVtx->GetIndices();
+  UShort_t *outindices = new UShort_t[outVtx->GetNContributors()];
+  Int_t j=0;
+  Bool_t copyindex;
+  for(Int_t k=0; k<inVtx->GetNIndices(); k++) {
+    copyindex=kTRUE;
+    for(Int_t l=0; l<ntrks; l++) {
+      trksTree->GetEvent(l);
+      if(inindices[k]==track->GetID()) copyindex=kFALSE;
+    }
+    if(copyindex) {
+      outindices[j] = inindices[k]; j++;
+    }
+  }
+  outVtx->SetIndices(outVtx->GetNContributors(),outindices);
+  delete [] outindices;
+
+  if(fDebug) {
+    printf("Vertex before removing tracks:\n");
+    inVtx->PrintStatus();
+    inVtx->PrintIndices();
+    printf("Vertex after removing tracks:\n");
+    outVtx->PrintStatus();
+    outVtx->PrintIndices();
+  }
+
+  return outVtx;
+}
+//---------------------------------------------------------------------------
+void AliVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped) 
+{
+//
+// Mark the tracks not to be used in the vertex reconstruction.
+// Tracks are identified by AliESDtrack::GetID()
+//
+  fNTrksToSkip = n;  fTrksToSkip = new Int_t[n]; 

   for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i]; 
   return; 
 }
 //---------------------------------------------------------------------------
   for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i]; 
   return; 
 }
 //---------------------------------------------------------------------------

-void  AliVertexerTracks::SetVtxStart(AliESDVertex *vtx) { 
+void  AliVertexerTracks::SetVtxStart(AliESDVertex *vtx) 
+{ 

 //
 // Set initial vertex knowledge
 //
   vtx->GetXYZ(fNominalPos);
   vtx->GetCovMatrix(fNominalCov);
 //
 // Set initial vertex knowledge
 //
   vtx->GetXYZ(fNominalPos);
   vtx->GetCovMatrix(fNominalCov);

+  SetConstraintOn();

   return; 
 }
 //---------------------------------------------------------------------------
   return; 
 }
 //---------------------------------------------------------------------------

-void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights){
-
+void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
+{

   // Calculate the point at minimum distance to prepared tracks 
   
   Double_t initPos[3];
   // Calculate the point at minimum distance to prepared tracks 
   
   Double_t initPos[3];


@@ -585,27 +746,27 @@ void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights){
   Double_t vett[3][3];
   Double_t det=GetDeterminant3X3(sum);
   
   Double_t vett[3][3];
   Double_t det=GetDeterminant3X3(sum);
   

-   if(det!=0){
-     for(Int_t zz=0;zz<3;zz++){
-       for(Int_t ww=0;ww<3;ww++){
-        for(Int_t kk=0;kk<3;kk++) vett[ww][kk]=sum[ww][kk];
-       }
-       for(Int_t kk=0;kk<3;kk++) vett[kk][zz]=dsum[kk];
-       initPos[zz]=GetDeterminant3X3(vett)/det;
-     }
-
-
-     for(Int_t i=0; i<knacc; i++){
-       Double_t p0[3]={0,0,0},p1[3]={0,0,0};
-       for(Int_t ii=0;ii<3;ii++){
-        p0[ii]=vectP0[i][ii];
-        p1[ii]=vectP1[i][ii];
-       }
-       sigma+=GetStrLinMinDist(p0,p1,initPos);
-     }
-
-     sigma=TMath::Sqrt(sigma);
-   }else{
+  if(det!=0){
+    for(Int_t zz=0;zz<3;zz++){
+      for(Int_t ww=0;ww<3;ww++){
+       for(Int_t kk=0;kk<3;kk++) vett[ww][kk]=sum[ww][kk];
+      }
+      for(Int_t kk=0;kk<3;kk++) vett[kk][zz]=dsum[kk];
+      initPos[zz]=GetDeterminant3X3(vett)/det;
+    }
+
+
+    for(Int_t i=0; i<knacc; i++){
+      Double_t p0[3]={0,0,0},p1[3]={0,0,0};
+      for(Int_t ii=0;ii<3;ii++){
+       p0[ii]=vectP0[i][ii];
+       p1[ii]=vectP1[i][ii];
+      }
+      sigma+=GetStrLinMinDist(p0,p1,initPos);
+    }
+
+    sigma=TMath::Sqrt(sigma);
+  }else{

     Warning("StrLinVertexFinderMinDist","Finder did not succed");
     sigma=999;
   }
     Warning("StrLinVertexFinderMinDist","Finder did not succed");
     sigma=999;
   }


@@ -616,7 +777,54 @@ void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights){
   fVert.SetNContributors(knacc);
 }
 //---------------------------------------------------------------------------
   fVert.SetNContributors(knacc);
 }
 //---------------------------------------------------------------------------

-void AliVertexerTracks::TooFewTracks(const AliESD* esdEvent) {
+Bool_t AliVertexerTracks::TrackToPoint(AliESDtrack *t,
+                                      TMatrixD &ri,TMatrixD &wWi) const 
+{
+//
+// Extract from the AliESDtrack the global coordinates ri and covariance matrix
+// wWi of the space point that it represents (to be used in VertexFitter())
+//
+
+  
+  Double_t rotAngle = t->GetAlpha();
+  if(rotAngle<0.) rotAngle += 2.*TMath::Pi();
+  Double_t cosRot = TMath::Cos(rotAngle);
+  Double_t sinRot = TMath::Sin(rotAngle);
+
+  ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot;
+  ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot;
+  ri(2,0) = t->GetZ();
+
+  // matrix to go from global (x,y,z) to local (y,z);
+  TMatrixD qQi(2,3);
+  qQi(0,0) = -sinRot;
+  qQi(0,1) = cosRot;
+  qQi(0,2) = 0.;
+  qQi(1,0) = 0.;
+  qQi(1,1) = 0.;
+  qQi(1,2) = 1.;
+
+  // covariance matrix of local (y,z) - inverted
+  TMatrixD uUi(2,2);
+  Double_t cc[15];
+  t->GetExternalCovariance(cc);
+  uUi(0,0) = cc[0];
+  uUi(0,1) = cc[1];
+  uUi(1,0) = cc[1];
+  uUi(1,1) = cc[2];
+  if(uUi.Determinant() <= 0.) return kFALSE;
+  TMatrixD uUiInv(TMatrixD::kInverted,uUi);
+  
+  // weights matrix: wWi = qQiT * uUiInv * qQi
+  TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
+  TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
+  wWi = m;
+
+  return kTRUE;
+} 
+//---------------------------------------------------------------------------
+void AliVertexerTracks::TooFewTracks(const AliESD* esdEvent) 
+{

 //
 // When the number of tracks is < fMinTracks
 //
 //
 // When the number of tracks is < fMinTracks
 //


@@ -645,14 +853,17 @@ void AliVertexerTracks::TooFewTracks(const AliESD* esdEvent) {
   Double_t tp[3];
   esdEvent->GetVertex()->GetTruePos(tp);
   fCurrentVertex->SetTruePos(tp);
   Double_t tp[3];
   esdEvent->GetVertex()->GetTruePos(tp);
   fCurrentVertex->SetTruePos(tp);

-  fCurrentVertex->SetTitle("VertexerTracksWithConstraint");
-  if(ncontr==-1||ncontr==-2) 
+  if(fConstraint) {
+    fCurrentVertex->SetTitle("VertexerTracksWithConstraint");
+  } else {

     fCurrentVertex->SetTitle("VertexerTracksNoConstraint");
     fCurrentVertex->SetTitle("VertexerTracksNoConstraint");

+  }

 
   return;
 }
 //---------------------------------------------------------------------------
 
   return;
 }
 //---------------------------------------------------------------------------

-void AliVertexerTracks::VertexFinder(Int_t optUseWeights) {
+void AliVertexerTracks::VertexFinder(Int_t optUseWeights) 
+{

 
   // Get estimate of vertex position in (x,y) from tracks DCA
  
 
   // Get estimate of vertex position in (x,y) from tracks DCA
  


@@ -692,7 +903,8 @@ void AliVertexerTracks::VertexFinder(Int_t optUseWeights) {
     //      AliStrLine *line2 = new AliStrLine();
     //  track2->ApproximateHelixWithLine(mindist,field,line2);
       Double_t distCA=line2->GetDCA(line1);
     //      AliStrLine *line2 = new AliStrLine();
     //  track2->ApproximateHelixWithLine(mindist,field,line2);
       Double_t distCA=line2->GetDCA(line1);

-      if(fDCAcut<=0 || (fDCAcut>0&&distCA<fDCAcut)){
+      //printf("%d   %d   %f\n",i,j,distCA);
+       if(fDCAcut<=0 || (fDCAcut>0&&distCA<fDCAcut)){

        Double_t pnt1[3],pnt2[3],crosspoint[3];
 
        if(optUseWeights<=0){
        Double_t pnt1[3],pnt2[3],crosspoint[3];
 
        if(optUseWeights<=0){


@@ -734,6 +946,7 @@ void AliVertexerTracks::VertexFinder(Int_t optUseWeights) {
   if(ncombi>0){
     for(Int_t jj=0;jj<3;jj++){
       initPos[jj] = aver[jj]/ncombi;
   if(ncombi>0){
     for(Int_t jj=0;jj<3;jj++){
       initPos[jj] = aver[jj]/ncombi;

+      //printf("%f\n",initPos[jj]);

       aversq[jj]/=ncombi;
       sigmasq[jj]=aversq[jj]-initPos[jj]*initPos[jj];
       sigma+=sigmasq[jj];
       aversq[jj]/=ncombi;
       sigmasq[jj]=aversq[jj]-initPos[jj]*initPos[jj];
       sigma+=sigmasq[jj];


@@ -749,20 +962,23 @@ void AliVertexerTracks::VertexFinder(Int_t optUseWeights) {
   fVert.SetNContributors(ncombi);
 }
 //---------------------------------------------------------------------------
   fVert.SetNContributors(ncombi);
 }
 //---------------------------------------------------------------------------

-void AliVertexerTracks::VertexFitter(Bool_t useNominalVtx) {
+void AliVertexerTracks::VertexFitter(Bool_t useConstraint) 
+{

 //
 // The optimal estimate of the vertex position is given by a "weighted 
 // average of tracks positions"
 //
 // The optimal estimate of the vertex position is given by a "weighted 
 // average of tracks positions"

-// Original method: CMS Note 97/0051
+// Original method: V. Karimaki, CMS Note 97/0051

 //
   Double_t initPos[3];
   fVert.GetXYZ(initPos);
 //
   Double_t initPos[3];
   fVert.GetXYZ(initPos);

+  Int_t arrEntries = (Int_t)fTrkArray.GetEntries();
+  if(arrEntries==1) useConstraint=kTRUE;

   if(fDebug) { 
     printf(" VertexFitter(): start\n");
     printf(" Number of tracks in array: %d\n",(Int_t)fTrkArray.GetEntriesFast());
     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(fDebug) { 
     printf(" VertexFitter(): start\n");
     printf(" Number of tracks in array: %d\n",(Int_t)fTrkArray.GetEntriesFast());
     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(useNominalVtx) 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])); 
+    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])); 

   }
 
   Int_t i,j,k,step=0;
   }
 
   Int_t i,j,k,step=0;


@@ -772,12 +988,8 @@ void AliVertexerTracks::VertexFitter(Bool_t useNominalVtx) {
   rv(1,0) = initPos[1];
   rv(2,0) = 0.;
   Double_t xlStart,alpha;
   rv(1,0) = initPos[1];
   rv(2,0) = 0.;
   Double_t xlStart,alpha;

-  Double_t rotAngle;
-  Double_t cosRot,sinRot;
-  Double_t cc[15];

   Int_t nUsedTrks;
   Int_t nUsedTrks;

-  Double_t chi2,chi2i;
-  Int_t arrEntries = (Int_t)fTrkArray.GetEntries();
+  Double_t chi2,chi2i,chi2b;

   AliESDtrack *t = 0;
   Int_t failed = 0;
 
   AliESDtrack *t = 0;
   Int_t failed = 0;
 


@@ -815,12 +1027,19 @@ void AliVertexerTracks::VertexFitter(Bool_t useNominalVtx) {
       for(j=0; j<3; j++) sumWi(j,i) = 0.;
     }
 
       for(j=0; j<3; j++) sumWi(j,i) = 0.;
     }
 

-
-    if(useNominalVtx) {
+    // mean vertex constraint
+    if(useConstraint) {

       for(i=0;i<3;i++) {
        sumWiri(i,0) += vVbInvrb(i,0);
        for(k=0;k<3;k++) sumWi(i,k) += vVbInv(i,k);
       }
       for(i=0;i<3;i++) {
        sumWiri(i,0) += vVbInvrb(i,0);
        for(k=0;k<3;k++) sumWi(i,k) += vVbInv(i,k);
       }

+      // chi2
+      TMatrixD deltar = rv; deltar -= rb;
+      TMatrixD vVbInvdeltar(vVbInv,TMatrixD::kMult,deltar);
+      chi2b = deltar(0,0)*vVbInvdeltar(0,0)+
+              deltar(1,0)*vVbInvdeltar(1,0)+
+             deltar(2,0)*vVbInvdeltar(2,0);
+      chi2 += chi2b;

     }
 
 
     }
 
 


@@ -830,40 +1049,19 @@ void AliVertexerTracks::VertexFitter(Bool_t useNominalVtx) {
       t = (AliESDtrack*)fTrkArray.At(k);
       alpha = t->GetAlpha();
       xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha);
       t = (AliESDtrack*)fTrkArray.At(k);
       alpha = t->GetAlpha();
       xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha);

-      t->AliExternalTrackParam::PropagateTo(xlStart,AliTracker::GetBz());   // to vtxSeed
-      rotAngle = alpha;
-      if(alpha<0.) rotAngle += 2.*TMath::Pi();
-      cosRot = TMath::Cos(rotAngle);
-      sinRot = TMath::Sin(rotAngle);
-      
+      // to vtxSeed (from finder)
+      t->AliExternalTrackParam::PropagateTo(xlStart,AliTracker::GetBz());   
+ 
+

       // vector of track global coordinates
       TMatrixD ri(3,1);
       // vector of track global coordinates
       TMatrixD ri(3,1);

-      ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot;
-      ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot;
-      ri(2,0) = t->GetZ();
-
-      // matrix to go from global (x,y,z) to local (y,z);
-      TMatrixD qQi(2,3);
-      qQi(0,0) = -sinRot;
-      qQi(0,1) = cosRot;
-      qQi(0,2) = 0.;
-      qQi(1,0) = 0.;
-      qQi(1,1) = 0.;
-      qQi(1,2) = 1.;
-
-      // covariance matrix of local (y,z) - inverted
-      TMatrixD uUi(2,2);
-      t->GetExternalCovariance(cc);
-      uUi(0,0) = cc[0];
-      uUi(0,1) = cc[1];
-      uUi(1,0) = cc[1];
-      uUi(1,1) = cc[2];
-
-      // weights matrix: wWi = qQiT * uUiInv * qQi
-      if(uUi.Determinant() <= 0.) continue;
-      TMatrixD uUiInv(TMatrixD::kInverted,uUi);
-      TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
-      TMatrixD wWi(qQi,TMatrixD::kTransposeMult,uUiInvQi);
+      // covariance matrix of ri
+      TMatrixD wWi(3,3);
+
+      // get space point from track
+      if(!TrackToPoint(t,ri,wWi)) continue;
+
+      TMatrixD wWiri(wWi,TMatrixD::kMult,ri); 

 
       // track chi2
       TMatrixD deltar = rv; deltar -= ri;
 
       // track chi2
       TMatrixD deltar = rv; deltar -= ri;


@@ -872,12 +1070,9 @@ void AliVertexerTracks::VertexFitter(Bool_t useNominalVtx) {
               deltar(1,0)*wWideltar(1,0)+
              deltar(2,0)*wWideltar(2,0);
 
               deltar(1,0)*wWideltar(1,0)+
              deltar(2,0)*wWideltar(2,0);
 

-

       // add to total chi2
       chi2 += chi2i;
 
       // add to total chi2
       chi2 += chi2i;
 

-      TMatrixD wWiri(wWi,TMatrixD::kMult,ri); 
-

       sumWiri += wWiri;
       sumWi   += wWi;
 
       sumWiri += wWiri;
       sumWi   += wWi;
 


@@ -938,49 +1133,74 @@ void AliVertexerTracks::VertexFitter(Bool_t useNominalVtx) {
   return;
 }
 //----------------------------------------------------------------------------
   return;
 }
 //----------------------------------------------------------------------------

-AliVertex* AliVertexerTracks::VertexForSelectedTracks(TTree *trkTree) {
+AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(TTree *trkTree,Bool_t optUseFitter,Bool_t optPropagate) 
+{

 //
 // Return vertex from tracks in trkTree
 //
 
 //
 // Return vertex from tracks in trkTree
 //
 

+  SetConstraintOff();
+

   // get tracks and propagate them to initial vertex position
   // get tracks and propagate them to initial vertex position

-  Int_t nTrksPrep = PrepareTracks(*trkTree,1);
+  Int_t nTrksPrep = PrepareTracks(*trkTree,0);

   if(nTrksPrep <  TMath::Max(2,fMinTracks) ) {
     if(fDebug) printf("TooFewTracks\n");
   if(nTrksPrep <  TMath::Max(2,fMinTracks) ) {
     if(fDebug) printf("TooFewTracks\n");

-    Double_t vtx[3]={0,0,0};
-    fVert.SetXYZ(vtx);
-    fVert.SetDispersion(999);
-    fVert.SetNContributors(-5);
-    fTrkArray.Delete();
-    return &fVert;
+    fCurrentVertex = new AliESDVertex(0.,0.,-1);
+    return fCurrentVertex;

   }
  
   }
  

-  // Set initial vertex position from ESD
-  if(fAlgo==1)  StrLinVertexFinderMinDist(1);
-  if(fAlgo==2)  StrLinVertexFinderMinDist(0);
-  if(fAlgo==3)  HelixVertexFinder();
-  if(fAlgo==4)  VertexFinder(1);
-  if(fAlgo==5)  VertexFinder(0);
-
+  switch (fAlgo) {
+    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;  
+  }

   
   

-  // set indices of used tracks
+  if(fDebug) printf(" Vertex finding completed\n");
+
+  // vertex fitter
+  if(optUseFitter){
+    VertexFitter(fConstraint);
+    if(fDebug) printf(" Vertex fit completed\n");
+  }else{
+    Double_t position[3]={fVert.GetXv(),fVert.GetYv(),fVert.GetZv()};
+    Double_t covmatrix[6]={0.,0.,0.,0.,0.,0.};
+    Double_t chi2=99999.;
+    Int_t    nUsedTrks=fVert.GetNContributors();
+    fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nUsedTrks);    
+  }
+  fCurrentVertex->SetDispersion(fVert.GetDispersion());
+
+
+  // set indices of used tracks and propagate track to found vertex

   UShort_t *indices = 0;
   AliESDtrack *eta = 0;
   UShort_t *indices = 0;
   AliESDtrack *eta = 0;

-  if(fVert.GetNContributors()>0) {
-    indices = new UShort_t[fVert.GetNContributors()];
+  if(fCurrentVertex->GetNContributors()>0) {
+    indices = new UShort_t[fCurrentVertex->GetNContributors()];

     for(Int_t jj=0;jj<(Int_t)fTrkArray.GetEntriesFast();jj++) {
       eta = (AliESDtrack*)fTrkArray.At(jj);
       indices[jj] = (UShort_t)eta->GetID();
     for(Int_t jj=0;jj<(Int_t)fTrkArray.GetEntriesFast();jj++) {
       eta = (AliESDtrack*)fTrkArray.At(jj);
       indices[jj] = (UShort_t)eta->GetID();

+      if(optPropagate&&optUseFitter){
+       if(TMath::Sqrt(fCurrentVertex->GetXv()*fCurrentVertex->GetXv()+fCurrentVertex->GetYv()*fCurrentVertex->GetYv())<3.) {
+         eta->RelateToVertex(fCurrentVertex,GetField(),100.);
+         if(fDebug) printf("Track %d propagated to found vertex\n",jj);
+       }else{
+         AliWarning("Found vertex outside beam pipe!");
+       }
+      }

     }
     }

-    fVert.SetIndices(fVert.GetNContributors(),indices);
+    fCurrentVertex->SetIndices(fCurrentVertex->GetNContributors(),indices);

   }
   delete [] indices;
   fTrkArray.Delete();
   
   }
   delete [] indices;
   fTrkArray.Delete();
   

-  return &fVert;
+  return fCurrentVertex;

 }
 //----------------------------------------------------------------------------
 }
 //----------------------------------------------------------------------------

-AliVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray) {
+AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray,Bool_t optUseFitter, Bool_t optPropagate) 
+{

 //
 // Return vertex from array of tracks
 //
 //
 // Return vertex from array of tracks
 //


@@ -989,12 +1209,8 @@ AliVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray) {
   Int_t nTrks = trkArray->GetEntriesFast();
   if(nTrks < TMath::Max(2,fMinTracks) ) {
     if(fDebug) printf("TooFewTracks\n");
   Int_t nTrks = trkArray->GetEntriesFast();
   if(nTrks < TMath::Max(2,fMinTracks) ) {
     if(fDebug) printf("TooFewTracks\n");

-    Double_t vtx[3]={0,0,0};
-    fVert.SetXYZ(vtx);
-    fVert.SetDispersion(999);
-    fVert.SetNContributors(-5);
-    fTrkArray.Delete();
-    return &fVert;
+    fCurrentVertex = new AliESDVertex(0.,0.,-1);
+    return fCurrentVertex;

   }
   TTree *trkTree = new TTree("TreeT","tracks");
   AliESDtrack *esdTrack = 0;
   }
   TTree *trkTree = new TTree("TreeT","tracks");
   AliESDtrack *esdTrack = 0;


@@ -1004,7 +1220,7 @@ AliVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray) {
     trkTree->Fill();
   }
     
     trkTree->Fill();
   }
     

-  AliVertex *vtx =  VertexForSelectedTracks(trkTree);
+  AliESDVertex *vtx =  VertexForSelectedTracks(trkTree,optUseFitter,optPropagate);

   delete trkTree;
   return vtx;
 }
   delete trkTree;
   return vtx;
 }

diff --git a/STEER/AliVertexerTracks.h b/STEER/AliVertexerTracks.h
index e541ea0259bf6cf38ac5ff902b37e91996a9b584..2ae81423b3548a7a1e0c409d8a8f2b178db4a1ca 100644 (file)
--- a/STEER/AliVertexerTracks.h
+++ b/STEER/AliVertexerTracks.h
@@ -22,10 +22,12 @@
  *****************************************************************************/
 
 #include "AliESDVertex.h"
  *****************************************************************************/
 
 #include "AliESDVertex.h"

+#include "AliESDtrack.h"

 #include "AliTracker.h"
 #include "AliLog.h"
 
 #include <TObjArray.h>
 #include "AliTracker.h"
 #include "AliLog.h"
 
 #include <TObjArray.h>

+#include <TMatrixD.h>

 
 class TTree; 
 class AliESD;
 
 class TTree; 
 class AliESD;


@@ -37,27 +39,33 @@ class AliVertexerTracks : public TObject {
   AliVertexerTracks(Double_t xStart, Double_t yStart); 
   virtual ~AliVertexerTracks();
 
   AliVertexerTracks(Double_t xStart, Double_t yStart); 
   virtual ~AliVertexerTracks();
 

-
-  // computes the vertex from the set of tracks in the tree
-  AliVertex* VertexForSelectedTracks(TTree *trkTree);
-  AliVertex* VertexForSelectedTracks(TObjArray *trkArray);

   AliESDVertex* FindPrimaryVertex(const AliESD *esdEvent);
   AliESDVertex* FindPrimaryVertex(const AliESD *esdEvent);

-  void  SetMinTracks(Int_t n=1) { fMinTracks = n; return; }
+  AliESDVertex* VertexForSelectedTracks(TTree *trkTree,Bool_t optUseFitter=kTRUE, Bool_t optPropagate=kTRUE);
+  AliESDVertex* VertexForSelectedTracks(TObjArray *trkArray,Bool_t optUseFitter=kTRUE, Bool_t optPropagate=kTRUE);
+  AliESDVertex* RemoveTracksFromVertex(AliESDVertex *inVtx,TTree *trksTree,Float_t *diamondxy); 
+  void  SetConstraintOff() { fConstraint=kFALSE; return; }
+  void  SetConstraintOn() { fConstraint=kTRUE; return; }
+  void  SetDebug(Int_t optdebug=0) { fDebug=optdebug; return; }
+  void  SetDCAcut(Double_t maxdca) { fDCAcut=maxdca; return; }
+  void  SetFinderAlgorithm(Int_t opt=1) { fAlgo=opt; return; }
+  void  SetITSRequired() { fITSin=kTRUE; return; }
+  void  SetITSrefitRequired() { fITSin=kTRUE;fITSrefit=kTRUE; return; }

   void  SetITSNotRequired() { fITSrefit=kFALSE;fITSin=kFALSE; return; }
   void  SetITSrefitNotRequired() { fITSrefit=kFALSE; return; }
   void  SetITSNotRequired() { fITSrefit=kFALSE;fITSin=kFALSE; return; }
   void  SetITSrefitNotRequired() { fITSrefit=kFALSE; return; }

+  void  SetMaxd0z0(Double_t maxd0z0=0.5) { fMaxd0z0=maxd0z0; return; }

   void  SetMinITSClusters(Int_t n=5) { fMinITSClusters = n; return; }
   void  SetMinITSClusters(Int_t n=5) { fMinITSClusters = n; return; }

+  void  SetMinTracks(Int_t n=1) { fMinTracks = n; return; }
+  void  SetNSigmad0(Double_t n=3) { fNSigma=n; return; }
+  Double_t GetNSigmad0() const { return fNSigma; }
+  void  SetOnlyFitter() { if(!fConstraint) AliFatal("Set constraint first!"); 
+     fOnlyFitter=kTRUE; return; }

   void  SetSkipTracks(Int_t n,Int_t *skipped);
   void  SetSkipTracks(Int_t n,Int_t *skipped);

-  void  SetDebug(Int_t optdebug=0) {fDebug=optdebug;}

   void  SetVtxStart(Double_t x=0,Double_t y=0,Double_t z=0) 
     { fNominalPos[0]=x; fNominalPos[1]=y; fNominalPos[2]=z; return; }
   void  SetVtxStartSigma(Double_t sx=3,Double_t sy=3,Double_t sz=6) 
     { fNominalCov[0]=sx*sx; fNominalCov[2]=sy*sy; fNominalCov[5]=sz*sz;
       fNominalCov[1]=0.; fNominalCov[3]=0.; fNominalCov[4]=0.; return; }
   void  SetVtxStart(AliESDVertex *vtx);
   void  SetVtxStart(Double_t x=0,Double_t y=0,Double_t z=0) 
     { fNominalPos[0]=x; fNominalPos[1]=y; fNominalPos[2]=z; return; }
   void  SetVtxStartSigma(Double_t sx=3,Double_t sy=3,Double_t sz=6) 
     { fNominalCov[0]=sx*sx; fNominalCov[2]=sy*sy; fNominalCov[5]=sz*sz;
       fNominalCov[1]=0.; fNominalCov[3]=0.; fNominalCov[4]=0.; return; }
   void  SetVtxStart(AliESDVertex *vtx);

-  void  SetDCAcut(Double_t maxdca) { fDCAcut=maxdca; return; }
-  void  SetFinderAlgorithm(Int_t opt=1) { fAlgo=opt; return; }
-  void  SetNSigmad0(Double_t n=3) { fNSigma=n; return; }
-  Double_t GetNSigmad0() const { return fNSigma; }

   static Double_t GetStrLinMinDist(Double_t *p0,Double_t *p1,Double_t *x0);
   static Double_t GetDeterminant3X3(Double_t matr[][3]);
   static void GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t (*m)[3],Double_t *d);
   static Double_t GetStrLinMinDist(Double_t *p0,Double_t *p1,Double_t *x0);
   static Double_t GetDeterminant3X3(Double_t matr[][3]);
   static void GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t (*m)[3],Double_t *d);


@@ -68,12 +76,14 @@ class AliVertexerTracks : public TObject {
     if(!AliTracker::GetFieldMap())
       AliFatal("Field map not set; use AliTracker::SetFieldMap()!");
     return AliTracker::GetBz(); } 
     if(!AliTracker::GetFieldMap())
       AliFatal("Field map not set; use AliTracker::SetFieldMap()!");
     return AliTracker::GetBz(); } 

-  Int_t    PrepareTracks(TTree &trkTree,Int_t OptImpParCut);
+  void     HelixVertexFinder();

   void     OneTrackVertFinder();
   void     OneTrackVertFinder();

+  Int_t    PrepareTracks(TTree &trkTree,Int_t optImpParCut);
+  Bool_t   TrackToPoint(AliESDtrack *t,
+                       TMatrixD &ri,TMatrixD &wWi) const;     

   void     VertexFinder(Int_t optUseWeights=0);
   void     VertexFinder(Int_t optUseWeights=0);

-  void     HelixVertexFinder();
-  void     StrLinVertexFinderMinDist(Int_t OptUseWeights=0);
-  void     VertexFitter(Bool_t useNominaVtx=kFALSE);
+  void     VertexFitter(Bool_t useConstraint=kFALSE);
+  void     StrLinVertexFinderMinDist(Int_t optUseWeights=0);

   void     TooFewTracks(const AliESD *esdEvent);
 
    
   void     TooFewTracks(const AliESD *esdEvent);
 
    


@@ -81,6 +91,10 @@ class AliVertexerTracks : public TObject {
   AliESDVertex *fCurrentVertex;  // ESD vertex after fitter
   Double_t  fNominalPos[3];   // initial knowledge on vertex position
   Double_t  fNominalCov[6];   // initial knowledge on vertex position
   AliESDVertex *fCurrentVertex;  // ESD vertex after fitter
   Double_t  fNominalPos[3];   // initial knowledge on vertex position
   Double_t  fNominalCov[6];   // initial knowledge on vertex position

+  Bool_t    fConstraint;      // true when "mean vertex" was set in 
+                              // fNominal ... and must be used in the fit
+  Bool_t    fOnlyFitter;      // primary with one fitter shot only
+                              // (use only with beam constraint)

   Int_t     fMinTracks;       // minimum number of tracks
   Int_t     fMinITSClusters;  // minimum number of ITS clusters per track
   TObjArray fTrkArray;        // array with tracks to be processed
   Int_t     fMinTracks;       // minimum number of tracks
   Int_t     fMinITSClusters;  // minimum number of ITS clusters per track
   TObjArray fTrkArray;        // array with tracks to be processed


@@ -89,6 +103,8 @@ class AliVertexerTracks : public TObject {
   Double_t  fDCAcut;          // maximum DCA between 2 tracks used for vertex
   Int_t     fAlgo;            // option for vertex finding algorythm
   Double_t  fNSigma;          // number of sigmas for d0 cut in PrepareTracks()
   Double_t  fDCAcut;          // maximum DCA between 2 tracks used for vertex
   Int_t     fAlgo;            // option for vertex finding algorythm
   Double_t  fNSigma;          // number of sigmas for d0 cut in PrepareTracks()

+  Double_t  fMaxd0z0;         // value [mm] for sqrt(d0d0+z0z0) cut 
+                              // in PrepareTracks(1) if fConstraint=kFALSE

   Bool_t    fITSin;           // if kTRUE (default), use only kITSin tracks
                               // if kFALSE, use all tracks (also TPC only)
   Bool_t    fITSrefit;        // if kTRUE (default), use only kITSrefit tracks
   Bool_t    fITSin;           // if kTRUE (default), use only kITSin tracks
                               // if kFALSE, use all tracks (also TPC only)
   Bool_t    fITSrefit;        // if kTRUE (default), use only kITSrefit tracks


@@ -111,7 +127,7 @@ class AliVertexerTracks : public TObject {
   AliVertexerTracks(const AliVertexerTracks & source);
   AliVertexerTracks & operator=(const AliVertexerTracks & source);
 
   AliVertexerTracks(const AliVertexerTracks & source);
   AliVertexerTracks & operator=(const AliVertexerTracks & source);
 

-  ClassDef(AliVertexerTracks,5) // 3D Vertexing with ESD tracks 
+  ClassDef(AliVertexerTracks,6) // 3D Vertexing with ESD tracks 

 };
 
 #endif
 };
 
 #endif

diff --git a/STEER/AliVertexerTracksTest.C b/STEER/AliVertexerTracksTest.C
index 7e91dc0e32b77891e7187b0615b1210836194143..4d2fdecff3ae6fa9d9cc56ec63b6068ab1f3b6f1 100644 (file)
--- a/STEER/AliVertexerTracksTest.C
+++ b/STEER/AliVertexerTracksTest.C
@@ -1,16 +1,17 @@
-//------------------------------------------------------------------------
-// Test macro for AliVertexerTracks.
-// Contains several functions. 
-//
-// A. Dainese - INFN Legnaro
-//------------------------------------------------------------------------

 void AliVertexerTracksTest(Double_t nSigma=3.,
 void AliVertexerTracksTest(Double_t nSigma=3.,

-                          Bool_t useMeanVtx=kTRUE,
-                          TString outname="AliVertexerTracksTest.root") {
+                          Bool_t useMeanVtx=kFALSE,
+                          Int_t itsin=1,
+                          Int_t itsrefit=1,
+                          Double_t maxd0z0=0.5,
+                          Int_t minitscls=5,
+                          Int_t mintracks=1,
+                          TString outname="AliVertexerTracksTest.root",
+                          Bool_t onlyfit=kFALSE) {

 //------------------------------------------------------------------------
 // Run vertex reconstruction and store results in histos and ntuple
 //------------------------------------------------------------------------
 //------------------------------------------------------------------------
 // Run vertex reconstruction and store results in histos and ntuple
 //------------------------------------------------------------------------

-
+// WITHOUT Kinematics.root
+//

   if (gAlice) {
     delete gAlice->GetRunLoader();
     delete gAlice; 
   if (gAlice) {
     delete gAlice->GetRunLoader();
     delete gAlice; 


@@ -35,7 +36,9 @@ void AliVertexerTracksTest(Double_t nSigma=3.,
     return;
   }
   gAlice=rl->GetAliRun();
     return;
   }
   gAlice=rl->GetAliRun();

-  rl->LoadKinematics();
+
+  //rl->LoadKinematics();
+

   // Get field from galice.root
   AliMagF *fiel = (AliMagF*)gAlice->Field();
   // Set the conversion constant between curvature and Pt
   // Get field from galice.root
   AliMagF *fiel = (AliMagF*)gAlice->Field();
   // Set the conversion constant between curvature and Pt


@@ -51,8 +54,11 @@ void AliVertexerTracksTest(Double_t nSigma=3.,
   Double_t errvtx[3]; 
   Double_t diffX,diffY,diffZ;
   Double_t diffXerr,diffYerr,diffZerr;
   Double_t errvtx[3]; 
   Double_t diffX,diffY,diffZ;
   Double_t diffXerr,diffYerr,diffZerr;

-  Double_t multiplicity;
+  Float_t multiplicity;
+  Float_t chi2red;
+  Int_t ntracklets;

   Int_t nc;
   Int_t nc;

+  Int_t triggered;

 
   // Histograms
   TH1F *hdiffX = new TH1F("hdiffX","x_{found} - x_{true} [#mum]",1000,-5000,5000);
 
   // Histograms
   TH1F *hdiffX = new TH1F("hdiffX","x_{found} - x_{true} [#mum]",1000,-5000,5000);


@@ -64,7 +70,7 @@ void AliVertexerTracksTest(Double_t nSigma=3.,
   TH1F *hdiffZerr = new TH1F("hdiffZerr","#Delta z/#sigma_{z}",100,-5,5);
 
   //TNtple
   TH1F *hdiffZerr = new TH1F("hdiffZerr","#Delta z/#sigma_{z}",100,-5,5);
 
   //TNtple

-  TNtuple *ntVtxRes = new TNtuple("ntVtxRes","Vertex Resolution's Ntupla with multiplicity","ntracks:nitstracks5or6:nitstracksFromStrange:nitstracksFromStrange5or6:diffX:diffY:diffZ:diffXerr:diffYerr:diffZerr:multiplicity:nc");
+  TNtuple *ntVtxRes = new TNtuple("ntVtxRes","Vertex Resolution's Ntupla with multiplicity","triggered:ntracks:nitstracks:nitstracks5or6:diffX:diffY:diffZ:diffXerr:diffYerr:diffZerr:multiplicity:nc:zMC:chi2red");

 
 
 
 
 
 


@@ -82,14 +88,20 @@ void AliVertexerTracksTest(Double_t nSigma=3.,
     initVertex = new AliESDVertex(pos,err);
   }
   AliVertexerTracks *vertexer = new AliVertexerTracks();
     initVertex = new AliESDVertex(pos,err);
   }
   AliVertexerTracks *vertexer = new AliVertexerTracks();

-  vertexer->SetITSrefitNotRequired();
+  vertexer->SetDebug(1); // set to 1 to see what it does

   vertexer->SetVtxStart(initVertex);
   vertexer->SetVtxStart(initVertex);

-  vertexer->SetMinTracks(2);
+  if(!useMeanVtx) vertexer->SetNoConstraint();
+  if(onlyfit) vertexer->SetOnlyFitter();

   vertexer->SetNSigmad0(nSigma);
   vertexer->SetNSigmad0(nSigma);

+  if(!itsrefit) vertexer->SetITSrefitNotRequired();
+  if(!itsin) vertexer->SetITSNotRequired();
+  vertexer->SetMinTracks(mintracks);
+  vertexer->SetMinITSClusters(minitscls);
+  vertexer->SetMaxd0z0(maxd0z0);

   //cout<<" Nsigma:  "<<vertexer->GetNSigmad0()<<endl;
   //cout<<" Nsigma:  "<<vertexer->GetNSigmad0()<<endl;

-  vertexer->SetDebug(1); // set to 1 to see what it does

   //----------------------------------------------------------
  
   //----------------------------------------------------------
  

+

   if(gSystem->AccessPathName("AliESDs.root",kFileExists)) {
     printf("AliESDs.root not found!\n"); 
     return;
   if(gSystem->AccessPathName("AliESDs.root",kFileExists)) {
     printf("AliESDs.root not found!\n"); 
     return;


@@ -131,8 +143,6 @@ void AliVertexerTracksTest(Double_t nSigma=3.,
     Double_t dNchdy = 0.;
     multiplicity=0.;
     Int_t nitstracks = 0;
     Double_t dNchdy = 0.;
     multiplicity=0.;
     Int_t nitstracks = 0;

-    Int_t nitstracksFromStrange = 0;
-    Int_t nitstracksFromStrange5or6 = 0;

     Int_t nitstracks1 = 0;
     Int_t nitstracks2 = 0;
     Int_t nitstracks3 = 0;
     Int_t nitstracks1 = 0;
     Int_t nitstracks2 = 0;
     Int_t nitstracks3 = 0;


@@ -140,42 +150,38 @@ void AliVertexerTracksTest(Double_t nSigma=3.,
     Int_t nitstracks5 = 0;
     Int_t nitstracks6 = 0;
     Int_t nitstracks5or6 = 0;
     Int_t nitstracks5 = 0;
     Int_t nitstracks6 = 0;
     Int_t nitstracks5or6 = 0;

-  
-    // loop on particles
-    for(Int_t pa=0; pa<npart; pa++) {   
-      TParticle *part = (TParticle*)gAlice->GetMCApp()->Particle(pa); 
-      Int_t pdg = part->GetPdgCode();      
-      Int_t apdg = TMath::Abs(pdg);
-      Double_t energy  = part->Energy();
-      if(energy>6900.) continue; // reject incoming protons
-      Double_t pz = part->Pz();
-      Double_t y = 0.5*TMath::Log((energy+pz+1.e-13)/(energy-pz+1.e-13)); 
-                                                 
-      // consider charged pi,K,p,el,mu       
-      if(apdg!=11 && apdg!=13 && apdg!=211 && apdg!=321 && apdg!=2212) continue;
-      // reject secondaries
-      if(TMath::Sqrt((part->Vx()-o[0])*(part->Vx()-o[0])+(part->Vy()-o[1])*(part->Vy()-o[1]))>0.0010) continue;
-      if(TMath::Abs(y)<1.0) dNchdy += 0.5; // count 1/2 of particles in |y|<1
-    } // end loop on particles
-    multiplicity=(Double_t)dNchdy;
-
-    printf(" dNch/dy = %f\n",dNchdy);
-    //===============================================================

 
 
     esdTree->GetEvent(iev);
 
 
     esdTree->GetEvent(iev);

+    triggered=0;
+    chi2red=0.;
+    ULong64_t triggerMask = esd->GetTriggerMask();
+    if(triggerMask&32 && ((triggerMask&1) || (triggerMask&2))) triggered=1;
+
+    // get # of tracklets
+    AliMultiplicity *alimult = esd->GetMultiplicity();
+    if(alimult) {
+      ntracklets = alimult->GetNumberOfTracklets();
+      for(Int_t l=0;l<alimult->GetNumberOfTracklets();l++) 
+       if(alimult->GetDeltaPhi(l)<-9998.) ntracklets--;
+    } else {
+      ntracklets = 0;
+    }
+    printf("Number of SPD tracklets in ESD %d  :  %d\n",iev,ntracklets);
+    multiplicity = (Float_t)ntracklets;
+

     Int_t ntracks = esd->GetNumberOfTracks();
     printf("Number of tracks in ESD %d   :   %d\n",iev,ntracks);
     if(ntracks==0) { 
     Int_t ntracks = esd->GetNumberOfTracks();
     printf("Number of tracks in ESD %d   :   %d\n",iev,ntracks);
     if(ntracks==0) { 

-      ntVtxRes->Fill(ntracks,nitstracks5or6,nitstracksFromStrange,nitstracksFromStrange5or6,diffX,diffY,diffZ,diffXerr,diffYerr,diffZerr,multiplicity,nc); 
+      ntVtxRes->Fill(triggered,ntracks,nitstracks,nitstracks5or6,diffX,diffY,diffZ,diffXerr,diffYerr,diffZerr,multiplicity,nc,truevtx[2],chi2red); 

       continue; 
     }
     
     // VERTEX    
     AliESDVertex *vertex = (AliESDVertex*)vertexer->FindPrimaryVertex(esd);
       continue; 
     }
     
     // VERTEX    
     AliESDVertex *vertex = (AliESDVertex*)vertexer->FindPrimaryVertex(esd);

-    //AliESDVertex *vertex = (AliESDVertex*)vertexer->FindPrimaryVertexOld(esd);

 
     nc = (Int_t)vertex->GetNContributors();
 
     nc = (Int_t)vertex->GetNContributors();

+    if(nc>=2) chi2red = vertex->GetChi2toNDF(); 

     printf("nc = %d\n",nc);
 
 
     printf("nc = %d\n",nc);
 
 


@@ -193,15 +199,6 @@ void AliVertexerTracksTest(Double_t nSigma=3.,
       if(npts==3) nitstracks3++;
       if(npts==2) nitstracks2++;
       if(npts==1) nitstracks1++;
       if(npts==3) nitstracks3++;
       if(npts==2) nitstracks2++;
       if(npts==1) nitstracks1++;

-      TParticle *part = (TParticle*)gAlice->GetMCApp()->Particle(TMath::Abs(esdtrack->GetLabel()));
-      if(part->GetFirstMother()>=0) {
-       Int_t mpdg = TMath::Abs(gAlice->GetMCApp()->Particle(part->GetFirstMother())->GetPdgCode());
-       //cout<<TMath::Abs(esdtrack->GetLabel())<<"   "<<mpdg<<endl;
-       if(mpdg==321||mpdg==311||mpdg==310||mpdg==3122||mpdg==3312||mpdg==3334) {
-         nitstracksFromStrange++;
-         if(npts>=5) nitstracksFromStrange5or6++;
-       }
-      }

     }
     printf("Number of kITSin tracks in ESD %d   :   %d\n",iev,nitstracks);
     printf("           6 pts: %d\n",nitstracks6);
     }
     printf("Number of kITSin tracks in ESD %d   :   %d\n",iev,nitstracks);
     printf("           6 pts: %d\n",nitstracks6);


@@ -210,10 +207,9 @@ void AliVertexerTracksTest(Double_t nSigma=3.,
     printf("           3 pts: %d\n",nitstracks3);
     printf("           2 pts: %d\n",nitstracks2);
     printf("           1 pts: %d\n",nitstracks1);
     printf("           3 pts: %d\n",nitstracks3);
     printf("           2 pts: %d\n",nitstracks2);
     printf("           1 pts: %d\n",nitstracks1);

-    printf("Number of kITSin tracks from s:   %d\n",nitstracksFromStrange);

 
 
 
 

-    if(nc>=2) {
+    if(nc>=1) {

       vertex->GetXYZ(vtx);
       vertex->GetSigmaXYZ(errvtx); 
       diffX = 10000.* (vtx[0] - truevtx[0]);
       vertex->GetXYZ(vtx);
       vertex->GetSigmaXYZ(errvtx); 
       diffX = 10000.* (vtx[0] - truevtx[0]);


@@ -230,7 +226,7 @@ void AliVertexerTracksTest(Double_t nSigma=3.,
       hdiffZerr->Fill(diffZerr);
     } 
 
       hdiffZerr->Fill(diffZerr);
     } 
 

-    ntVtxRes->Fill(nitstracks,nitstracks5or6,nitstracksFromStrange,nitstracksFromStrange5or6,diffX,diffY,diffZ,diffXerr,diffYerr,diffZerr,multiplicity,(Float_t)nc);   
+    ntVtxRes->Fill(triggered,ntracks,nitstracks,nitstracks5or6,diffX,diffY,diffZ,diffXerr,diffYerr,diffZerr,multiplicity,(Float_t)nc,truevtx[2],chi2red);   

 
   } // END LOOP ON EVENTS
 
 
   } // END LOOP ON EVENTS
 


@@ -254,6 +250,154 @@ void AliVertexerTracksTest(Double_t nSigma=3.,
   return;
 } 
 //----------------------------------------------------------------------------
   return;
 } 
 //----------------------------------------------------------------------------

+//------------------------------------------------------------------------
+void VertexForOneEvent(Int_t iev=0,
+                      Double_t nSigma=3.,
+                      Bool_t useMeanVtx=kFALSE) {
+//------------------------------------------------------------------------
+// Run vertex reconstruction for event iev
+//------------------------------------------------------------------------
+
+  if (gAlice) {
+    delete gAlice->GetRunLoader();
+    delete gAlice; 
+    gAlice=0;
+  }
+  
+  AliRunLoader *rl = AliRunLoader::Open("galice.root");
+  if (rl == 0x0) {
+    cerr<<"Can not open session"<<endl;
+    return;
+  }
+  Int_t retval = rl->LoadgAlice();
+  if (retval) {
+    cerr<<"LoadgAlice returned error"<<endl;
+    delete rl;
+    return;
+  }
+  retval = rl->LoadHeader();
+  if (retval) {
+    cerr<<"LoadHeader returned error"<<endl;
+    delete rl;
+    return;
+  }
+  gAlice=rl->GetAliRun();
+  //rl->LoadKinematics();
+  // Get field from galice.root
+  AliMagF *fiel = (AliMagF*)gAlice->Field();
+  // Set the conversion constant between curvature and Pt
+  AliTracker::SetFieldMap(fiel,kTRUE);
+ 
+  Double_t truevtx[3];
+  Double_t vtx[3]; 
+  Double_t errvtx[3]; 
+  Double_t diffX,diffY,diffZ;
+  Double_t diffXerr,diffYerr,diffZerr;
+  Double_t multiplicity;
+  Int_t nc;
+
+  // -----------   Create vertexer --------------------------
+  AliESDVertex *initVertex = 0;
+  if(useMeanVtx) {
+    // open the mean vertex
+    TFile *invtx = new TFile("AliESDVertexMean.root");
+    initVertex = (AliESDVertex*)invtx->Get("vtxmean");
+    invtx->Close();
+    delete invtx;
+  } else {
+    Double_t pos[3]={0.,-0.,0.}; 
+    Double_t err[3]={3.,3.,5.};
+    initVertex = new AliESDVertex(pos,err);
+  }
+  AliVertexerTracks *vertexer = new AliVertexerTracks();
+  //vertexer->SetITSrefitNotRequired();
+  vertexer->SetVtxStart(initVertex);
+  if(!useMeanVtx) vertexer->SetNoConstraint();
+  //vertexer->SetMinTracks(2);
+  //vertexer->SetNSigmad0(nSigma);
+  //cout<<" Nsigma:  "<<vertexer->GetNSigmad0()<<endl;
+  //vertexer->SetFinderAlgorithm(2);
+  //vertexer->SetDCAcut(0.1); // 1 mm
+  vertexer->SetDebug(1); // set to 1 to see what it does
+  //----------------------------------------------------------
+ 
+  if(gSystem->AccessPathName("AliESDs.root",kFileExists)) {
+    printf("AliESDs.root not found!\n"); 
+    return;
+  }
+  TFile *fin = TFile::Open("AliESDs.root");
+  TTree *esdTree = (TTree*)fin->Get("esdTree");
+  if(!esdTree) return;
+  AliESD *esd = 0;
+  esdTree->SetBranchAddress("ESD",&esd);
+
+  TArrayF o;
+
+  nc=0;
+
+  //=================================================
+  //         LOOK IN THE SIMULATION
+  // true vertex position
+  Int_t npart = (Int_t)gAlice->GetEvent(iev);
+  printf("particles  %d\n",npart);
+  AliHeader *header = (AliHeader*)rl->GetHeader();
+  AliGenPythiaEventHeader *pyheader = (AliGenPythiaEventHeader*)header->GenEventHeader();
+  pyheader->PrimaryVertex(o);
+  truevtx[0] = o[0];
+  truevtx[1] = o[1];
+  truevtx[2] = o[2];
+  printf("True Vertex: (%f, %f, %f)\n",o[0],o[1],o[2]);
+
+  esdTree->GetEvent(iev);
+  Int_t ntracks = esd->GetNumberOfTracks();
+  printf("Number of tracks in ESD %d   :   %d\n",iev,ntracks);
+    
+  // VERTEX    
+  AliESDVertex *vertex = (AliESDVertex*)vertexer->FindPrimaryVertex(esd);
+  nc = (Int_t)vertex->GetNContributors();
+
+
+  Int_t nitstracks = 0;
+  Int_t nitstracks1 = 0;
+  Int_t nitstracks2 = 0;
+  Int_t nitstracks3 = 0;
+  Int_t nitstracks4 = 0;
+  Int_t nitstracks5 = 0;
+  Int_t nitstracks6 = 0;
+  Int_t nitstracks5or6 = 0;
+
+  // count ITS tracks
+  for(Int_t itrk=0; itrk<ntracks; itrk++) {
+    AliESDtrack *esdtrack = (AliESDtrack*)esd->GetTrack(itrk);
+    UInt_t status = esdtrack->GetStatus();
+    // only tracks found also in ITS
+    if(! (status&AliESDtrack::kITSrefit) ) continue;      
+    nitstracks++;
+    Int_t npts = (Int_t)esdtrack->GetNcls(0);
+    if(npts==6) {nitstracks6++;nitstracks5or6++;}
+    if(npts==5) {nitstracks5++;nitstracks5or6++;}
+    if(npts==4) nitstracks4++;
+    if(npts==3) nitstracks3++;
+    if(npts==2) nitstracks2++;
+    if(npts==1) nitstracks1++;
+  }
+  printf("Number of kITSrefit tracks in ESD %d   :   %d\n",iev,nitstracks);
+  printf("           6 pts: %d\n",nitstracks6);
+  printf("           5 pts: %d\n",nitstracks5);
+  printf("           4 pts: %d\n",nitstracks4);
+  printf("           3 pts: %d\n",nitstracks3);
+  printf("           2 pts: %d\n",nitstracks2);
+  printf("           1 pts: %d\n",nitstracks1);
+  
+
+  delete esdTree;
+  delete fin;
+  delete vertexer;
+  delete initVertex;
+
+  return;
+} 
+//----------------------------------------------------------------------------

 Double_t FitFunc(Double_t *x,Double_t *par) {
   return par[0]+par[1]/TMath::Sqrt(x[0]);
 }
 Double_t FitFunc(Double_t *x,Double_t *par) {
   return par[0]+par[1]/TMath::Sqrt(x[0]);
 }


@@ -265,20 +409,52 @@ Int_t GetBin(Float_t mult) {
   if(mult<22.) return 4;
   return 5;
 }
   if(mult<22.) return 4;
   return 5;
 }

-void PlotVtxRes(const Char_t *inname="AliVertexerTracksTest.root") {
+Int_t GetBinTracklets(Float_t tracklets) {
+  if(tracklets<2.*5.)  return 0;
+  if(tracklets<2.*7.)  return 1;
+  if(tracklets<2.*12.) return 2;
+  if(tracklets<2.*15.) return 3;
+  if(tracklets<2.*22.) return 4;
+  return 5;
+}
+Int_t GetBinZ(Float_t z) {
+  if(z<-13.) return 0;
+  if(z<-11.) return 1;
+  if(z<-9.) return 2;
+  if(z<-7.) return 3;
+  if(z<-5.) return 4;
+  if(z<-3.) return 5;
+  if(z<-1.) return 6;
+  if(z<1.) return 7;
+  if(z<3.) return 8;
+  if(z<5.) return 9;
+  if(z<7.) return 10;
+  if(z<9.) return 11;
+  if(z<11.) return 12;
+  if(z<13.) return 13;
+  return 14;
+}
+//----------------------------------------------------------------------------
+void PlotVtxRes(TString inname="AliVertexerTracksTest.root") {

   //---------------------------------------------------------------------
   // Plot efficiency, resolutions, pulls 
   // [at least 10k pp events are needed]
   //---------------------------------------------------------------------
 
   //---------------------------------------------------------------------
   // Plot efficiency, resolutions, pulls 
   // [at least 10k pp events are needed]
   //---------------------------------------------------------------------
 

+  Float_t zMCmin=0.;
+  Float_t zMCmax=20.;
+  Float_t ncminforzshift=1.;
+

   gStyle->SetOptStat(0);
   gStyle->SetOptFit(10001);
 
   Int_t   nEvVtx=0;
   Int_t   nEvNoVtx=0;
   Int_t   ev,nUsedTrks;
   gStyle->SetOptStat(0);
   gStyle->SetOptFit(10001);
 
   Int_t   nEvVtx=0;
   Int_t   nEvNoVtx=0;
   Int_t   ev,nUsedTrks;

-  Float_t nTrks,nTrksFromStrange,nTrks5or6,nTrksFromStrange5or6,dNchdy;
-  Float_t diffX,diffY,diffZ,diffXerr,diffYerr,diffZerr;
+  Float_t nESDTrks,nTrks,nTrks5or6,ntracklets;
+  Float_t diffX,diffY,diffZ,diffXerr,diffYerr,diffZerr,zMC,nc;
+  Float_t triggered;
+

 
   //
   // HISTOGRAMS
 
   //
   // HISTOGRAMS


@@ -291,10 +467,28 @@ void PlotVtxRes(const Char_t *inname="AliVertexerTracksTest.root") {
   hdy->SetXTitle("#Delta y [#mu m]");
   hdy->SetYTitle("events");
   //
   hdy->SetXTitle("#Delta y [#mu m]");
   hdy->SetYTitle("events");
   //

-  TH1F *hdz = new TH1F("hdz","",50,-600,600);
+  TH1F *hdz = new TH1F("hdz","",200,-600,600);

   hdz->SetXTitle("#Delta z [#mu m]");
   hdz->SetYTitle("events");
   //
   hdz->SetXTitle("#Delta z [#mu m]");
   hdz->SetYTitle("events");
   //

+  TH1F *hzMCVtx = new TH1F("hzMCVtx","events w/ vertex",30,-15,15);
+  hzMCVtx->SetXTitle("z_{true} [cm]");
+  hzMCVtx->SetYTitle("events");
+  hzMCVtx->Sumw2();
+  //
+  TH1F *hzMCNoVtx = new TH1F("hzMCNoVtx","events w/o vertex",30,-15,15);
+  hzMCNoVtx->SetXTitle("z_{true} [cm]");
+  hzMCNoVtx->SetYTitle("events");
+  hzMCNoVtx->Sumw2();
+  //
+  TH1F *hTrackletsVtx = new TH1F("hTrackletsVtx","events w/ vertex",51,-0.5,50.5);
+  hTrackletsVtx->SetXTitle("SPD tracklets");
+  hTrackletsVtx->SetYTitle("events");
+  //
+  TH1F *hTrackletsNoVtx = new TH1F("hTrackletsNoVtx","events w/o vertex",51,-0.5,50.5);
+  hTrackletsNoVtx->SetXTitle("SPD tracklets");
+  hTrackletsNoVtx->SetYTitle("events");
+  //

   TH1F *hTracksVtx = new TH1F("hTracksVtx","events w/ vertex",51,-0.5,50.5);
   hTracksVtx->SetXTitle("Number of reco tracks in TPC&ITS");
   hTracksVtx->SetYTitle("events");
   TH1F *hTracksVtx = new TH1F("hTracksVtx","events w/ vertex",51,-0.5,50.5);
   hTracksVtx->SetXTitle("Number of reco tracks in TPC&ITS");
   hTracksVtx->SetYTitle("events");


@@ -304,21 +498,13 @@ void PlotVtxRes(const Char_t *inname="AliVertexerTracksTest.root") {
   hTracksNoVtx->SetYTitle("events");
   //
   TH1F *hTracksVtx5or6 = new TH1F("hTracksVtx5or6","events w/ vertex",51,-0.5,50.5);
   hTracksNoVtx->SetYTitle("events");
   //
   TH1F *hTracksVtx5or6 = new TH1F("hTracksVtx5or6","events w/ vertex",51,-0.5,50.5);

-  hTracksVtx5or6->SetXTitle("Number of reco tracks in TPC&ITS(pts>=5)");
+  hTracksVtx5or6->SetXTitle("Number of reco tracks in TPC&ITS(cls>=5)");

   hTracksVtx5or6->SetYTitle("events");
   //
   TH1F *hTracksNoVtx5or6 = new TH1F("hTracksNoVtx5or6","events w/o vertex",51,-0.5,50.5);
   hTracksVtx5or6->SetYTitle("events");
   //
   TH1F *hTracksNoVtx5or6 = new TH1F("hTracksNoVtx5or6","events w/o vertex",51,-0.5,50.5);

-  hTracksNoVtx5or6->SetXTitle("Number of reco tracks in TPC&ITS(pts>=5)");
+  hTracksNoVtx5or6->SetXTitle("Number of reco tracks in TPC&ITS(cls>=5)");

   hTracksNoVtx5or6->SetYTitle("events");
   //
   hTracksNoVtx5or6->SetYTitle("events");
   //

-  TH1F *hTracksVtx5or6nonS = new TH1F("hTracksVtx5or6nonS","events w/ vertex",51,-0.5,50.5);
-  hTracksVtx5or6nonS->SetXTitle("Number of reco tracks in TPC&ITS(pts>=5) NOT from s");
-  hTracksVtx5or6nonS->SetYTitle("events");
-  //
-  TH1F *hTracksNoVtx5or6nonS = new TH1F("hTracksNoVtx5or6nonS","events w/o vertex",51,-0.5,50.5);
-  hTracksNoVtx5or6nonS->SetXTitle("Number of reco tracks in TPC&ITS(pts>=5) NOT from s");
-  hTracksNoVtx5or6nonS->SetYTitle("events");
-  //

   TH1F *hPullsx = new TH1F("hPullsx","",50,-7.,7.);
   hPullsx->SetXTitle("#Delta x/#sqrt{<x,x>}");
   hPullsx->SetYTitle("events");
   TH1F *hPullsx = new TH1F("hPullsx","",50,-7.,7.);
   hPullsx->SetXTitle("#Delta x/#sqrt{<x,x>}");
   hPullsx->SetYTitle("events");


@@ -331,6 +517,17 @@ void PlotVtxRes(const Char_t *inname="AliVertexerTracksTest.root") {
   hPullsz->SetXTitle("#Delta z/#sqrt{<z,z>}");
   hPullsz->SetYTitle("events");
   //
   hPullsz->SetXTitle("#Delta z/#sqrt{<z,z>}");
   hPullsz->SetYTitle("events");
   //

+  TProfile *hntrackletsVSzMC = new TProfile("hntrackletsVSzMC","",30,-15,15,0,30);
+  hntrackletsVSzMC->SetXTitle("z_{true} [cm]");
+  hntrackletsVSzMC->SetYTitle("<SPD tracklets>");
+  //
+  TProfile *hnitstracksVSzMC = new TProfile("hnitstracksVSzMC","",30,-15,15,0,30);
+  hnitstracksVSzMC->SetXTitle("z_{true} [cm]");
+  hnitstracksVSzMC->SetYTitle("<tracks TPC&ITS>");
+  //
+  TProfile *hnitstracks5or6VSzMC = new TProfile("hnitstracks5or6VSzMC","",30,-15,15,0,30);
+  hnitstracks5or6VSzMC->SetXTitle("z_{true} [cm]");
+  hnitstracks5or6VSzMC->SetYTitle("<tracks TPC&ITS(cls>=5)>");

 
   Double_t mult[6]={0,0,0,0,0,0};
   Double_t mult2[6]={0,0,0,0,0,0};
 
   Double_t mult[6]={0,0,0,0,0,0};
   Double_t mult2[6]={0,0,0,0,0,0};


@@ -358,7 +555,7 @@ void PlotVtxRes(const Char_t *inname="AliVertexerTracksTest.root") {
   Int_t    all[6]={0,0,0,0,0,0};
   Double_t probVtx[6]={0,0,0,0,0,0};
   Double_t eprobVtx[6]={0,0,0,0,0,0};
   Int_t    all[6]={0,0,0,0,0,0};
   Double_t probVtx[6]={0,0,0,0,0,0};
   Double_t eprobVtx[6]={0,0,0,0,0,0};

-  Int_t    bin;
+  Int_t    bin,zbin;

   //
   // OTHER HISTOGRAMS
   //
   //
   // OTHER HISTOGRAMS
   //


@@ -380,6 +577,10 @@ void PlotVtxRes(const Char_t *inname="AliVertexerTracksTest.root") {
   TH1F *hdz3 = new TH1F("hdz3","z resolution - bin 3",50,-600,600);
   TH1F *hdz4 = new TH1F("hdz4","z resolution - bin 4",50,-500,500);
   TH1F *hdz5 = new TH1F("hdz5","z resolution - bin 5",50,-500,500);
   TH1F *hdz3 = new TH1F("hdz3","z resolution - bin 3",50,-600,600);
   TH1F *hdz4 = new TH1F("hdz4","z resolution - bin 4",50,-500,500);
   TH1F *hdz5 = new TH1F("hdz5","z resolution - bin 5",50,-500,500);

+  //
+  TH1F *hdz_z = NULL; hdz_z = new TH1F[15];
+  for(Int_t l=0;l<15;l++) hdz_z[l] = *hdz;
+

 
   TH1F *hpx0 = new TH1F("hpx0","x pull - bin 0",50,-7,7);
   TH1F *hpx1 = new TH1F("hpx1","x pull - bin 1",50,-7,7);
 
   TH1F *hpx0 = new TH1F("hpx0","x pull - bin 0",50,-7,7);
   TH1F *hpx1 = new TH1F("hpx1","x pull - bin 1",50,-7,7);


@@ -403,31 +604,42 @@ void PlotVtxRes(const Char_t *inname="AliVertexerTracksTest.root") {
 
   TH1F *hmult = new TH1F("hmult","hmult",50,-0.5,49.5);
 
 
   TH1F *hmult = new TH1F("hmult","hmult",50,-0.5,49.5);
 

-  TFile *in = new TFile(inname);
+  TFile *in = new TFile(inname.Data());

   TNtuple *nt = (TNtuple*)in->Get("ntVtxRes");
   TNtuple *nt = (TNtuple*)in->Get("ntVtxRes");

-  nt->SetBranchAddress("ntracks",&nTrks);
+  nt->SetBranchAddress("triggered",&triggered);
+  nt->SetBranchAddress("ntracks",&nESDTrks);
+  nt->SetBranchAddress("nitstracks",&nTrks);

   nt->SetBranchAddress("nitstracks5or6",&nTrks5or6);
   nt->SetBranchAddress("nitstracks5or6",&nTrks5or6);

-  nt->SetBranchAddress("nitstracksFromStrange",&nTrksFromStrange);
-  nt->SetBranchAddress("nitstracksFromStrange5or6",&nTrksFromStrange5or6);

   nt->SetBranchAddress("diffX",&diffX);
   nt->SetBranchAddress("diffY",&diffY);
   nt->SetBranchAddress("diffZ",&diffZ);
   nt->SetBranchAddress("diffXerr",&diffXerr);
   nt->SetBranchAddress("diffYerr",&diffYerr);
   nt->SetBranchAddress("diffZerr",&diffZerr);
   nt->SetBranchAddress("diffX",&diffX);
   nt->SetBranchAddress("diffY",&diffY);
   nt->SetBranchAddress("diffZ",&diffZ);
   nt->SetBranchAddress("diffXerr",&diffXerr);
   nt->SetBranchAddress("diffYerr",&diffYerr);
   nt->SetBranchAddress("diffZerr",&diffZerr);

-  nt->SetBranchAddress("multiplicity",&dNchdy);
+  nt->SetBranchAddress("multiplicity",&ntracklets);
+  nt->SetBranchAddress("zMC",&zMC);
+  nt->SetBranchAddress("nc",&nc);

   Int_t entries = (Int_t)nt->GetEntries();
   Int_t nbytes=0;
 
   for(Int_t i=0; i<entries; i++) {
     nbytes += nt->GetEvent(i);
   Int_t entries = (Int_t)nt->GetEntries();
   Int_t nbytes=0;
 
   for(Int_t i=0; i<entries; i++) {
     nbytes += nt->GetEvent(i);

-    bin = GetBin(dNchdy);
-    mult[bin] += dNchdy;
-    hmult->Fill(dNchdy);
+    // only triggered events
+    if(triggered<0.5) continue;
+    // consider only events with zMCmin<|zMC|<zMCmax
+    if(TMath::Abs(zMC)<zMCmin || TMath::Abs(zMC)>zMCmax) continue;
+    //
+    bin = GetBinTracklets(ntracklets);
+    zbin = GetBinZ(zMC);
+    mult[bin] += ntracklets;
+    hmult->Fill(ntracklets);

     totevts[bin]++;
     totevts[bin]++;

-    if(diffX < 90000.) { // vtx found
+    hntrackletsVSzMC->Fill(zMC,ntracklets);
+    hnitstracksVSzMC->Fill(zMC,nTrks);
+    hnitstracks5or6VSzMC->Fill(zMC,nTrks5or6);
+    if(TMath::Abs(diffX) < 1000.) { // vtx found - closer than 1 mm

       nEvVtx++;
       nEvVtx++;

-      mult2[bin] += dNchdy;
+      mult2[bin] += ntracklets;

       vtx[bin]++;
       tottracks[bin] += nTrks;
       evts[bin]++;
       vtx[bin]++;
       tottracks[bin] += nTrks;
       evts[bin]++;


@@ -452,9 +664,11 @@ void PlotVtxRes(const Char_t *inname="AliVertexerTracksTest.root") {
       hdx->Fill(diffX);
       hdy->Fill(diffY);
       hdz->Fill(diffZ);
       hdx->Fill(diffX);
       hdy->Fill(diffY);
       hdz->Fill(diffZ);

+      if(nc>=ncminforzshift) hdz_z[zbin].Fill(diffZ);
+      hzMCVtx->Fill(zMC);
+      hTrackletsVtx->Fill(ntracklets);

       hTracksVtx->Fill(nTrks);
       hTracksVtx5or6->Fill(nTrks5or6);
       hTracksVtx->Fill(nTrks);
       hTracksVtx5or6->Fill(nTrks5or6);

-      hTracksVtx5or6nonS->Fill(nTrks5or6-nTrksFromStrange5or6);

       hPullsx->Fill(diffXerr);
       hPullsy->Fill(diffYerr);
       hPullsz->Fill(diffZerr);
       hPullsx->Fill(diffXerr);
       hPullsy->Fill(diffYerr);
       hPullsz->Fill(diffZerr);


@@ -478,9 +692,10 @@ void PlotVtxRes(const Char_t *inname="AliVertexerTracksTest.root") {
       if(bin==5) hpz5->Fill(diffZerr);
     } else {
       nEvNoVtx++;
       if(bin==5) hpz5->Fill(diffZerr);
     } else {
       nEvNoVtx++;

+      hzMCNoVtx->Fill(zMC);
+      hTrackletsNoVtx->Fill(ntracklets);

       hTracksNoVtx->Fill(nTrks);
       hTracksNoVtx5or6->Fill(nTrks5or6);
       hTracksNoVtx->Fill(nTrks);
       hTracksNoVtx5or6->Fill(nTrks5or6);

-      hTracksNoVtx5or6nonS->Fill(nTrks5or6-nTrksFromStrange5or6);

     }
     all[bin]++;
   }
     }
     all[bin]++;
   }


@@ -498,9 +713,9 @@ void PlotVtxRes(const Char_t *inname="AliVertexerTracksTest.root") {
   // calculate spread in multiplicity
   for(Int_t i=0; i<entries; i++) {
     nbytes += nt->GetEvent(i);
   // calculate spread in multiplicity
   for(Int_t i=0; i<entries; i++) {
     nbytes += nt->GetEvent(i);

-    bin = GetBin(dNchdy);
-    sigmamult[bin] += (dNchdy-mult[bin])*(dNchdy-mult[bin])/totevts[bin];
-    if(diffX < 90000.) sigmamult2[bin] += (dNchdy-mult2[bin])*(dNchdy-mult2[bin])/evts[bin];
+    bin = GetBinTracklets(ntracklets);
+    sigmamult[bin] += (ntracklets-mult[bin])*(ntracklets-mult[bin])/totevts[bin];
+    if(diffX < 90000.) sigmamult2[bin] += (ntracklets-mult2[bin])*(ntracklets-mult2[bin])/evts[bin];

   }
 
   for(Int_t k=0;k<6;k++) {
   }
 
   for(Int_t k=0;k<6;k++) {


@@ -518,7 +733,7 @@ void PlotVtxRes(const Char_t *inname="AliVertexerTracksTest.root") {
   g->SetLineWidth(3);
   TF1 *line = new TF1("line","pol1");
   line->SetLineWidth(3);
   g->SetLineWidth(3);
   TF1 *line = new TF1("line","pol1");
   line->SetLineWidth(3);

-  TF1 *func = new TF1("func",FitFunc,2,35,2);
+  TF1 *func = new TF1("func",FitFunc,2,80,2);

   func->SetLineWidth(1);
 
   Double_t x1,y1,sigmafit;
   func->SetLineWidth(1);
 
   Double_t x1,y1,sigmafit;


@@ -529,6 +744,23 @@ void PlotVtxRes(const Char_t *inname="AliVertexerTracksTest.root") {
   //
   gStyle->SetOptFit(1111);
   
   //
   gStyle->SetOptFit(1111);
   

+  // tracks VS zMC
+  TCanvas *c0 = new TCanvas("c0","c0",0,0,1000,400);
+  c0->Divide(3,1);
+  c0->cd(1); 
+  hntrackletsVSzMC->SetMinimum(0);
+  hntrackletsVSzMC->SetMaximum(14);
+  hntrackletsVSzMC->Draw("profs");
+  c0->cd(2); 
+  hnitstracksVSzMC->SetMinimum(0);
+  hnitstracksVSzMC->SetMaximum(14);
+  hnitstracksVSzMC->Draw("profs");
+  c0->cd(3); 
+  hnitstracks5or6VSzMC->SetMinimum(0);
+  hnitstracks5or6VSzMC->SetMaximum(14);
+  hnitstracks5or6VSzMC->Draw("profs");
+
+

   // Tracks in ITS for events w/ and w/o vertex
   TCanvas *c1 = new TCanvas("c1","c1",0,0,1000,500);
   c1->Divide(2,1);
   // Tracks in ITS for events w/ and w/o vertex
   TCanvas *c1 = new TCanvas("c1","c1",0,0,1000,500);
   c1->Divide(2,1);


@@ -537,6 +769,16 @@ void PlotVtxRes(const Char_t *inname="AliVertexerTracksTest.root") {
   c1->cd(2); 
   hTracksNoVtx->Draw();
 
   c1->cd(2); 
   hTracksNoVtx->Draw();
 

+  // probability vs ntracklets
+  TCanvas *c1a = new TCanvas("c1a","c1a",0,0,500,500);
+  TH1F *hTotTracklets = (TH1F*)hTrackletsVtx->Clone("hTotTracklets");
+  hTotTracklets->Add(hTrackletsNoVtx);
+  TH1F *hProbTracklets = (TH1F*)hTrackletsVtx->Clone("hProbTracklets");
+  hProbTracklets->Divide(hTotTracklets);
+  hProbTracklets->SetYTitle("Probability");
+  hProbTracklets->SetTitle("Probability to find the vertex");
+  hProbTracklets->Draw();
+

   // probability vs ntracks
   TCanvas *c1b = new TCanvas("c1b","c1b",0,0,500,500);
   TH1F *hTotTracks = (TH1F*)hTracksVtx->Clone("hTotTracks");
   // probability vs ntracks
   TCanvas *c1b = new TCanvas("c1b","c1b",0,0,500,500);
   TH1F *hTotTracks = (TH1F*)hTracksVtx->Clone("hTotTracks");


@@ -556,14 +798,16 @@ void PlotVtxRes(const Char_t *inname="AliVertexerTracksTest.root") {
   hProbTracks5or6->SetTitle("Probability to find the vertex");
   hProbTracks5or6->Draw();
 
   hProbTracks5or6->SetTitle("Probability to find the vertex");
   hProbTracks5or6->Draw();
 

-  TCanvas *c1d = new TCanvas("c1d","c1d",0,0,500,500);
-  TH1F *hTotTracks5or6nonS = (TH1F*)hTracksVtx5or6nonS->Clone("hTotTracks5or6nonS");
-  hTotTracks5or6nonS->Add(hTracksNoVtx5or6nonS);
-  TH1F *hProbTracks5or6nonS = (TH1F*)hTracksVtx5or6nonS->Clone("hProbTracks5or6nonS");
-  hProbTracks5or6nonS->Divide(hTotTracks5or6nonS);
-  hProbTracks5or6nonS->SetYTitle("Probability");
-  hProbTracks5or6nonS->SetTitle("Probability to find the vertex");
-  hProbTracks5or6nonS->Draw();
+
+  // probability vs zMC
+  TCanvas *c1e = new TCanvas("c1e","c1e",0,0,500,500);
+  TH1F *hTotzMC = (TH1F*)hzMCVtx->Clone("hTotzMC");
+  hTotzMC->Add(hzMCNoVtx);
+  TH1F *hProbzMC = (TH1F*)hzMCVtx->Clone("hProbzMC");
+  hProbzMC->Divide(hTotzMC);
+  hProbzMC->SetYTitle("Probability");
+  hProbzMC->SetTitle("Probability to find the vertex");
+  hProbzMC->Draw();

 
   // Global resolutions
   TCanvas *c2 = new TCanvas("c2","c2",0,0,1000,400);
 
   // Global resolutions
   TCanvas *c2 = new TCanvas("c2","c2",0,0,1000,400);


@@ -599,8 +843,8 @@ void PlotVtxRes(const Char_t *inname="AliVertexerTracksTest.root") {
   
   // probability VS multiplicity
   TCanvas *c5 = new TCanvas("c5","c5",0,0,600,500);
   
   // probability VS multiplicity
   TCanvas *c5 = new TCanvas("c5","c5",0,0,600,500);

-  TH2F *fr5 = new TH2F("fr5","",2,0,35,2,0,1.1); 
-  fr5->SetXTitle("dN_{ch}/dy");
+  TH2F *fr5 = new TH2F("fr5","",2,0,80,2,0,1.1); 
+  fr5->SetXTitle("SPD tracklets");

   fr5->SetYTitle("Probability");
   fr5->Draw();
   TGraphErrors *gr5 = new TGraphErrors(6,mult,probVtx,sigmamult,eprobVtx);
   fr5->SetYTitle("Probability");
   fr5->Draw();
   TGraphErrors *gr5 = new TGraphErrors(6,mult,probVtx,sigmamult,eprobVtx);


@@ -852,21 +1096,21 @@ void PlotVtxRes(const Char_t *inname="AliVertexerTracksTest.root") {
 
   gStyle->SetOptFit(0);
 
 
   gStyle->SetOptFit(0);
 

-  // res VS dNchdy
+  // res VS ntracklets

   TCanvas *c6 = new TCanvas("c6","c6",0,0,600,500);
   TCanvas *c6 = new TCanvas("c6","c6",0,0,600,500);

-  TH2F *fr6 = new TH2F("fr6","",2,0,35,2,0,200); 
-  fr6->SetXTitle("dN_{ch}/dy");
+  TH2F *fr6 = new TH2F("fr6","",2,0,80,2,0,240); 
+  fr6->SetXTitle("SPD tracklets");

   fr6->SetYTitle("#sigma [#mu m]");
   fr6->Draw();
   sigmamult2[0]=sigmamult2[1]=sigmamult2[2]=sigmamult2[3]=sigmamult2[4]=sigmamult2[5]=0.;
   TGraphErrors *gr6x = new TGraphErrors(6,mult2,xres,sigmamult2,exres);
   gr6x->Draw("p");
   gr6x->SetMarkerStyle(22);
   fr6->SetYTitle("#sigma [#mu m]");
   fr6->Draw();
   sigmamult2[0]=sigmamult2[1]=sigmamult2[2]=sigmamult2[3]=sigmamult2[4]=sigmamult2[5]=0.;
   TGraphErrors *gr6x = new TGraphErrors(6,mult2,xres,sigmamult2,exres);
   gr6x->Draw("p");
   gr6x->SetMarkerStyle(22);

-  gr6x->Fit("func","E,R");
+  //gr6x->Fit("func","E,R");

   TGraphErrors *gr6z = new TGraphErrors(6,mult2,zres,sigmamult2,ezres);
   gr6z->Draw("p");
   gr6z->SetMarkerStyle(26);
   TGraphErrors *gr6z = new TGraphErrors(6,mult2,zres,sigmamult2,ezres);
   gr6z->Draw("p");
   gr6z->SetMarkerStyle(26);

-  gr6z->Fit("func","E,R");
+  //gr6z->Fit("func","E,R");

   TLegend *leg6 = new TLegend(0.6,0.75,0.9,0.9);
   leg6->AddEntry(gr6x,"x/y coordinate","p");
   leg6->AddEntry(gr6z,"z coordinate","p");
   TLegend *leg6 = new TLegend(0.6,0.75,0.9,0.9);
   leg6->AddEntry(gr6x,"x/y coordinate","p");
   leg6->AddEntry(gr6z,"z coordinate","p");


@@ -874,10 +1118,10 @@ void PlotVtxRes(const Char_t *inname="AliVertexerTracksTest.root") {
   leg6->SetBorderSize(0);
   leg6->Draw();
 
   leg6->SetBorderSize(0);
   leg6->Draw();
 

-  // pull VS dNchdy
+  // pull VS ntracklets

   TCanvas *c8 = new TCanvas("c8","c8",0,0,600,500);
   TCanvas *c8 = new TCanvas("c8","c8",0,0,600,500);

-  TH2F *fr8 = new TH2F("fr8","",2,0,35,2,0,2.); 
-  fr8->SetXTitle("dN_{ch}/dy");
+  TH2F *fr8 = new TH2F("fr8","",2,0,80,2,0,2.); 
+  fr8->SetXTitle("SPD tracklets");

   fr8->SetYTitle("pull");
   fr8->Draw();
   TGraphErrors *gr8x = new TGraphErrors(6,mult2,xpull,sigmamult2,expull);
   fr8->SetYTitle("pull");
   fr8->Draw();
   TGraphErrors *gr8x = new TGraphErrors(6,mult2,xpull,sigmamult2,expull);


@@ -892,13 +1136,373 @@ void PlotVtxRes(const Char_t *inname="AliVertexerTracksTest.root") {
   leg8->SetFillStyle(0);
   leg8->SetBorderSize(0);
   leg8->Draw();
   leg8->SetFillStyle(0);
   leg8->SetBorderSize(0);
   leg8->Draw();

-  TLine *l8 = new TLine(0,1,35,1);
+  TLine *l8 = new TLine(0,1,80,1);

   l8->SetLineStyle(2);
   l8->Draw();
 
   l8->SetLineStyle(2);
   l8->Draw();
 

+
+  // mean z res VS zMC
+  Float_t zmc[15]={-14,-12,-10,-8,-6,-4,-2,0,2,4,6,8,10,12,14};
+  Float_t ezmc[15]; 
+  Float_t dzmean[15],edzmean[15];
+  Int_t dummy;
+  TCanvas *zz = new TCanvas("zz","zz",0,0,1000,1000);
+  zz->Divide(5,3);
+  for(l=0;l<15;l++) {
+    zz->cd(l+1);
+    hdz_z[l].Draw(); 
+    g->SetRange(-3.*hdz_z[l].GetRMS(),+3.*hdz_z[l].GetRMS());
+    hdz_z[l].Fit("g","Q"); 
+    dzmean[l]=g->GetParameter(1);
+    edzmean[l]=g->GetParError(1);
+    ezmc[l]=1.;
+  }
+  TCanvas *zzz = new TCanvas("zzz","zzz",0,0,600,500);
+  TH2F *fr9 = new TH2F("fr9","",2,-15,15,2,-50,50); 
+  fr9->SetXTitle("z_{true} [cm]");
+  fr9->SetYTitle("<z_{rec} - z_{true}> [#mu m]");
+  fr9->Draw();
+  TGraphErrors *grz = new TGraphErrors(15,zmc,dzmean,ezmc,edzmean);
+  grz->Draw("pz");
+  grz->SetMarkerStyle(20);
+
+

   return;
 }
 //----------------------------------------------------------------------------
   return;
 }
 //----------------------------------------------------------------------------

+void SaveFigures(TString name="dummy") {
+
+  TString namefull;
+
+  namefull = name.Data();
+  namefull.Prepend("plots/residualsx_");
+  namefull.Append(".gif");
+  ax->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/residualsx_");
+  namefull.Append(".ps");
+  ax->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/residualsx_");
+  namefull.Append(".eps");
+  ax->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/residualsx_");
+  namefull.Append(".root");
+  ax->Print(namefull.Data());
+
+  namefull = name.Data();
+  namefull.Prepend("plots/pullsx_");
+  namefull.Append(".gif");
+  bx->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/pullsx_");
+  namefull.Append(".ps");
+  bx->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/pullsx_");
+  namefull.Append(".eps");
+  bx->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/pullsx_");
+  namefull.Append(".root");
+  bx->Print(namefull.Data());
+
+  namefull = name.Data();
+  namefull.Prepend("plots/residualsz_");
+  namefull.Append(".gif");
+  az->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/residualsz_");
+  namefull.Append(".ps");
+  az->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/residualsz_");
+  namefull.Append(".eps");
+  az->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/residualsz_");
+  namefull.Append(".root");
+  az->Print(namefull.Data());
+
+  namefull = name.Data();
+  namefull.Prepend("plots/pullsz_");
+  namefull.Append(".gif");
+  bz->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/pullsz_");
+  namefull.Append(".ps");
+  bz->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/pullsz_");
+  namefull.Append(".eps");
+  bz->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/pullsz_");
+  namefull.Append(".root");
+  bz->Print(namefull.Data());
+
+  namefull = name.Data();
+  namefull.Prepend("plots/effVSntracklets_");
+  namefull.Append(".gif");
+  c1a->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/effVSntracklets_");
+  namefull.Append(".ps");
+  c1a->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/effVSntracklets_");
+  namefull.Append(".eps");
+  c1a->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/effVSntracklets_");
+  namefull.Append(".root");
+  c1a->Print(namefull.Data());
+
+
+  namefull = name.Data();
+  namefull.Prepend("plots/effVSnitstracks_");
+  namefull.Append(".gif");
+  c1b->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/effVSnitstracks_");
+  namefull.Append(".ps");
+  c1b->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/effVSnitstracks_");
+  namefull.Append(".eps");
+  c1b->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/effVSnitstracks_");
+  namefull.Append(".root");
+  c1b->Print(namefull.Data());
+
+
+  namefull = name.Data();
+  namefull.Prepend("plots/effVSnitstracks5or6_");
+  namefull.Append(".gif");
+  c1c->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/effVSnitstracks5or6_");
+  namefull.Append(".ps");
+  c1c->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/effVSnitstracks5or6_");
+  namefull.Append(".eps");
+  c1c->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/effVSnitstracks5or6_");
+  namefull.Append(".root");
+  c1c->Print(namefull.Data());
+
+
+  namefull = name.Data();
+  namefull.Prepend("plots/effVSzMC_");
+  namefull.Append(".gif");
+  c1e->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/effVSzMC_");
+  namefull.Append(".ps");
+  c1e->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/effVSzMC_");
+  namefull.Append(".eps");
+  c1e->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/effVSzMC_");
+  namefull.Append(".root");
+  c1e->Print(namefull.Data());
+
+
+
+  namefull = name.Data();
+  namefull.Prepend("plots/effVSntracklets2_");
+  namefull.Append(".gif");
+  c5->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/effVSntracklets2_");
+  namefull.Append(".ps");
+  c5->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/effVSntracklets2_");
+  namefull.Append(".eps");
+  c5->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/effVSntracklets2_");
+  namefull.Append(".root");
+  c5->Print(namefull.Data());
+
+
+  namefull = name.Data();
+  namefull.Prepend("plots/sigmaVSntracklets_");
+  namefull.Append(".gif");
+  c6->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/sigmaVSntracklets_");
+  namefull.Append(".ps");
+  c6->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/sigmaVSntracklets_");
+  namefull.Append(".eps");
+  c6->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/sigmaVSntracklets_");
+  namefull.Append(".root");
+  c6->Print(namefull.Data());
+
+
+  namefull = name.Data();
+  namefull.Prepend("plots/pullsVSntracklets_");
+  namefull.Append(".gif");
+  c8->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/pullsVSntracklets_");
+  namefull.Append(".ps");
+  c8->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/pullsVSntracklets_");
+  namefull.Append(".eps");
+  c8->Print(namefull.Data());
+  namefull = name.Data();
+  namefull.Prepend("plots/pullsVSntracklets_");
+  namefull.Append(".root");
+  c8->Print(namefull.Data());
+
+
+  return;
+}
+//----------------------------------------------------------------------------
+void TestRmTracks(Int_t iev=0) {
+
+  if (gAlice) {
+    delete gAlice->GetRunLoader();
+    delete gAlice; 
+    gAlice=0;
+  }
+  
+  AliRunLoader *rl = AliRunLoader::Open("galice.root");
+  if (rl == 0x0) {
+    cerr<<"Can not open session"<<endl;
+    return;
+  }
+  Int_t retval = rl->LoadgAlice();
+  if (retval) {
+    cerr<<"LoadgAlice returned error"<<endl;
+    delete rl;
+    return;
+  }
+  retval = rl->LoadHeader();
+  if (retval) {
+    cerr<<"LoadHeader returned error"<<endl;
+    delete rl;
+    return;
+  }
+  gAlice=rl->GetAliRun();
+  // Get field from galice.root
+  AliMagF *fiel = (AliMagF*)gAlice->Field();
+  // Set the conversion constant between curvature and Pt
+  AliTracker::SetFieldMap(fiel,kTRUE);
+ 
+  Double_t truevtx[3];
+  Double_t vtx[3]; 
+  Double_t errvtx[3]; 
+  Double_t diffX,diffY,diffZ;
+  Double_t diffXerr,diffYerr,diffZerr;
+  Double_t multiplicity;
+  Int_t nc;
+
+  // -----------   Create vertexer --------------------------
+  AliESDVertex *initVertex = 0;
+  TFile *invtx = new TFile("AliESDVertexMean.root");
+  initVertex = (AliESDVertex*)invtx->Get("vtxmean");
+  invtx->Close();
+  delete invtx;
+
+  AliVertexerTracks *vertexer = new AliVertexerTracks();
+  vertexer->SetVtxStart(initVertex);
+  vertexer->SetOnlyFitter();
+  vertexer->SetDebug(0); // set to 1 to see what it does
+
+  Float_t diamondxy[2];
+  diamondxy[0]=initVertex->GetXv();
+  diamondxy[1]=initVertex->GetYv();
+  //----------------------------------------------------------
+ 
+  if(gSystem->AccessPathName("AliESDs.root",kFileExists)) {
+    printf("AliESDs.root not found!\n"); 
+    return;
+  }
+  TFile *fin = TFile::Open("AliESDs.root");
+  TTree *esdTree = (TTree*)fin->Get("esdTree");
+  if(!esdTree) return;
+  AliESD *esd = 0;
+  esdTree->SetBranchAddress("ESD",&esd);
+
+  TArrayF o;
+
+  nc=0;
+
+
+  esdTree->GetEvent(iev);
+  Int_t ntracks = esd->GetNumberOfTracks();
+  printf("Number of tracks in ESD %d   :   %d\n",iev,ntracks);
+  
+  // example: do vertex without tracks 1 and 2 (AliESDtrack::GetID())
+  //
+  // 1: tell vertexer to skip those two tracks
+  //
+  Int_t nskip=2;
+  Int_t skipped[2];
+  skipped[0]=1;
+  skipped[1]=2;
+  vertexer->SetSkipTracks(nskip,skipped);  
+  AliESDVertex *vertex1 = (AliESDVertex*)vertexer->FindPrimaryVertex(esd);
+  vertex1->PrintStatus();
+  vertex1->PrintIndices();
+  delete vertex1;
+  vertex1 = 0;
+  //
+  // 2: do vertex with all tracks
+  //
+  skipped[0]=-99;
+  skipped[1]=-99;
+  vertexer->SetSkipTracks(nskip,skipped);  
+  AliESDVertex *vertex2 = (AliESDVertex*)vertexer->FindPrimaryVertex(esd);
+  vertex2->PrintStatus();
+  vertex2->PrintIndices();
+  //
+  // 3: now remove those two tracks from vertex2
+  //
+  TTree *trksTree = new TTree("trksTree","tree with tracks to be removed");
+  AliESDtrack *esdTrack = 0;
+  trksTree->Branch("tracks","AliESDtrack",&esdTrack);
+  AliESDtrack *et = esd->GetTrack(1);
+  esdTrack = new AliESDtrack(*et);
+  trksTree->Fill();
+  delete esdTrack;
+  AliESDtrack *et = esd->GetTrack(2);
+  esdTrack = new AliESDtrack(*et);
+  trksTree->Fill();
+  delete esdTrack;
+  AliVertexerTracks vt;
+  AliESDVertex *vertex3 = vt.RemoveTracksFromVertex(vertex2,trksTree,diamondxy);
+  vertex3->PrintStatus();
+  vertex3->PrintIndices();
+  delete vertex2;
+  vertex2 = 0;
+  delete vertex3;
+  vertex3 = 0;
+
+
+  delete esdTree;
+  delete fin;
+  vertexer =0;
+  delete vertexer;
+  delete initVertex;
+
+  return;
+} 
+//----------------------------------------------------------------------------

 void AliComputeVtxMeanFromESD(TString file="AliESDs.root",  
                              Int_t mincontr=20) {
   //-----------------------------------------------------------------------
 void AliComputeVtxMeanFromESD(TString file="AliESDs.root",  
                              Int_t mincontr=20) {
   //-----------------------------------------------------------------------