When calculating a*a-b*b the form (a-b)*(a+b) is usually more numerically stable.
[u/mrichter/AliRoot.git] / STEER / AliVertexerTracks.cxx
index de798a1..97a7626 100644 (file)
 #include <TTree.h>
 #include <TMatrixD.h>
 //---- AliRoot headers -----
+#include "AliLog.h"
 #include "AliStrLine.h"
 #include "AliExternalTrackParam.h"
+#include "AliNeutralTrackParam.h"
+#include "AliVEvent.h"
+#include "AliVTrack.h"
 #include "AliESDEvent.h"
 #include "AliESDtrack.h"
 #include "AliVertexerTracks.h"
@@ -66,8 +70,8 @@ fITSrefit(kTRUE),
 fFiducialR(3.),
 fFiducialZ(30.),
 fnSigmaForUi00(1.5),
-fDebug(0),
-fAlgo(1)
+fAlgo(1),
+fAlgoIter0(4)
 {
 //
 // Default constructor
@@ -100,15 +104,14 @@ fITSrefit(kTRUE),
 fFiducialR(3.),
 fFiducialZ(30.),
 fnSigmaForUi00(1.5),
-fDebug(0),
-fAlgo(1)
+fAlgo(1),
+fAlgoIter0(4)
 {
 //
 // Standard constructor
 //
   SetVtxStart();
   SetVtxStartSigma();
-  SetTPCMode();
 }
 //-----------------------------------------------------------------------------
 AliVertexerTracks::~AliVertexerTracks() 
@@ -121,10 +124,10 @@ AliVertexerTracks::~AliVertexerTracks()
   if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
 }
 //----------------------------------------------------------------------------
-AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliESDEvent *esdEvent)
+AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliVEvent *vEvent)
 {
 //
-// Primary vertex for current ESD event
+// Primary vertex for current ESD or AOD event
 // (Two iterations: 
 //  1st with 5*fNSigma*sigma cut w.r.t. to initial vertex
 //      + cut on sqrt(d0d0+z0z0) if fConstraint=kFALSE  
@@ -132,11 +135,20 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliESDEvent *esdEvent)
 //
   fCurrentVertex = 0;
 
+  TString evtype = vEvent->IsA()->GetName();
+  Bool_t inputAOD = ((evtype=="AliAODEvent") ? kTRUE : kFALSE);
+
+  if(inputAOD && fMode==1) {
+    printf("Error : AliVertexerTracks: no TPC-only vertex from AOD\n"); 
+    TooFewTracks(); 
+    return fCurrentVertex;
+  }
+
   // accept 1-track case only if constraint is available
   if(!fConstraint && fMinTracks==1) fMinTracks=2;
 
-  // read tracks from ESD
-  Int_t nTrks = (Int_t)esdEvent->GetNumberOfTracks();
+  // read tracks from AlivEvent
+  Int_t nTrks = (Int_t)vEvent->GetNumberOfTracks();
   if(nTrks<fMinTracks) {
     TooFewTracks();
     return fCurrentVertex;
@@ -151,38 +163,44 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliESDEvent *esdEvent)
 
   Int_t nTrksOrig=0;
   AliExternalTrackParam *t=0;
-  // loop on ESD tracks
+  // loop on tracks
   for(Int_t i=0; i<nTrks; i++) {
-    AliESDtrack *esdt = esdEvent->GetTrack(i);
+    AliVTrack *track = (AliVTrack*)vEvent->GetTrack(i);
     // check tracks to skip
     Bool_t skipThis = kFALSE;
     for(Int_t j=0; j<fNTrksToSkip; j++) { 
-      if(esdt->GetID()==fTrksToSkip[j]) {
-       if(fDebug) printf("skipping track: %d\n",i);
+      if(track->GetID()==fTrksToSkip[j]) {
+       AliDebug(1,Form("skipping track: %d",i));
        skipThis = kTRUE;
       }
     }
     if(skipThis) continue;
 
-    // check number of clusters in ITS or TPC
-    if(esdt->GetNcls(fMode) < fMinClusters) continue;
-
-    if(fMode==0) {        // ITS mode
-      if(fITSrefit && !(esdt->GetStatus()&AliESDtrack::kITSrefit)) continue;
-      Double_t x,p[5],cov[15];
-      esdt->GetExternalParameters(x,p);
-      esdt->GetExternalCovariance(cov);
-      t = new AliExternalTrackParam(x,esdt->GetAlpha(),p,cov);
-    } else if(fMode==1) { // TPC mode
-      t = (AliExternalTrackParam*)esdt->GetTPCInnerParam();
-      if(!t) continue;
-      Double_t radius = 2.8; //something less than the beam pipe radius
-      if(!PropagateTrackTo(t,radius)) continue;
+    if(!inputAOD) {  // ESD
+      AliESDtrack* esdt = (AliESDtrack*)track;
+      if(esdt->GetNcls(fMode) < fMinClusters) continue;
+      if(fMode==0) {        // ITS mode
+       if(fITSrefit && !(esdt->GetStatus()&AliESDtrack::kITSrefit)) continue;
+       Double_t x,p[5],cov[15];
+       esdt->GetExternalParameters(x,p);
+       esdt->GetExternalCovariance(cov);
+       t = new AliExternalTrackParam(x,esdt->GetAlpha(),p,cov);
+      } else if(fMode==1) { // TPC mode
+       t = (AliExternalTrackParam*)esdt->GetTPCInnerParam();
+       if(!t) continue;
+       Double_t radius = 2.8; //something less than the beam pipe radius
+       if(!PropagateTrackTo(t,radius)) continue;
+      }
+    } else {          // AOD (only ITS mode)
+      Int_t ncls0=0;
+      for(Int_t l=0;l<6;l++) if(TESTBIT(track->GetITSClusterMap(),l)) ncls0++;
+      if(ncls0 < fMinClusters) continue;
+      t = new AliExternalTrackParam(track);
     }
     trkArrayOrig.AddLast(t);
-    idOrig[nTrksOrig]=(UShort_t)esdt->GetID();
+    idOrig[nTrksOrig]=(UShort_t)track->GetID();
     nTrksOrig++;
-  } // end loop on ESD tracks
+  } // end loop on tracks
   
   // call method that will reconstruct the vertex
   FindPrimaryVertex(&trkArrayOrig,idOrig);
@@ -217,9 +235,9 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(TObjArray *trkArrayOrig,
 
   // read tracks from array
   Int_t nTrksOrig = (Int_t)trkArrayOrig->GetEntriesFast();
-  if(fDebug) printf("Initial number of tracks: %d\n",nTrksOrig);
+  AliDebug(1,Form("Initial number of tracks: %d",nTrksOrig));
   if(nTrksOrig<fMinTracks) {
-    if(fDebug) printf("TooFewTracks\n");
+    AliDebug(1,"TooFewTracks");
     TooFewTracks();
     return fCurrentVertex;
   } 
@@ -233,19 +251,27 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(TObjArray *trkArrayOrig,
     if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
     Double_t cutsave = fDCAcut;  
     fDCAcut = fDCAcutIter0;
-    VertexFinder(1); // using weights, cutting dca < fDCAcutIter0
+    // vertex finder
+    switch (fAlgoIter0) {
+    case 1: StrLinVertexFinderMinDist(1); break;
+    case 2: StrLinVertexFinderMinDist(0); break;
+    case 3: HelixVertexFinder();          break;
+    case 4: VertexFinder(1);              break;
+    case 5: VertexFinder(0);              break;
+    default: printf("Wrong algorithm\n"); break;  
+    }
     fDCAcut = cutsave;
     if(fVert.GetNContributors()>0) {
       fVert.GetXYZ(fNominalPos);
       fNominalPos[0] = fVert.GetXv();
       fNominalPos[1] = fVert.GetYv();
       fNominalPos[2] = fVert.GetZv();
-      if(fDebug) printf("No mean vertex: VertexFinder gives (%f, %f, %f)\n",fNominalPos[0],fNominalPos[1],fNominalPos[2]);
+      AliDebug(1,Form("No mean vertex: VertexFinder gives (%f, %f, %f)",fNominalPos[0],fNominalPos[1],fNominalPos[2]));
     } else {
       fNominalPos[0] = 0.;
       fNominalPos[1] = 0.;
       fNominalPos[2] = 0.;
-      if(fDebug) printf("No mean vertex and VertexFinder failed\n");
+      AliDebug(1,"No mean vertex and VertexFinder failed");
     }
   }
   
@@ -265,7 +291,7 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(TObjArray *trkArrayOrig,
     if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
     fIdSel = new UShort_t[nTrksOrig];
     Int_t nTrksSel = PrepareTracks(*trkArrayOrig,idOrig,iter);
-    if(fDebug) printf("N tracks selected in iteration %d: %d\n",iter,nTrksSel);
+    AliDebug(1,Form("N tracks selected in iteration %d: %d",iter,nTrksSel));
     if(nTrksSel < fMinTracks) {
       TooFewTracks();
       return fCurrentVertex; 
@@ -274,7 +300,7 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(TObjArray *trkArrayOrig,
     // vertex finder
     if(!fOnlyFitter) {
       if(nTrksSel==1) {
-       if(fDebug) printf("Just one track\n");
+       AliDebug(1,"Just one track");
        OneTrackVertFinder();
       } else {
        switch (fAlgo) {
@@ -283,10 +309,10 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(TObjArray *trkArrayOrig,
         case 3: HelixVertexFinder();          break;
         case 4: VertexFinder(1);              break;
         case 5: VertexFinder(0);              break;
-        default: printf("Wrong algorithm\n"); break;  
+        default: printf("Wrong algorithm"); break;  
        }
       }
-      if(fDebug) printf(" Vertex finding completed\n");
+      AliDebug(1," Vertex finding completed");
     }
 
     // vertex fitter
@@ -315,10 +341,7 @@ AliESDVertex* AliVertexerTracks::FindPrimaryVertex(TObjArray *trkArrayOrig,
   fCurrentVertex->SetTitle(title.Data());
   //
 
-  if(fDebug) {
-    fCurrentVertex->PrintStatus();
-    fCurrentVertex->PrintIndices();
-  }
+  AliDebug(1,Form("xyz: %f %f %f; nc %d",fCurrentVertex->GetXv(),fCurrentVertex->GetYv(),fCurrentVertex->GetZv(),fCurrentVertex->GetNContributors()));
 
   // clean up
   delete [] fIdSel; fIdSel=NULL;
@@ -420,7 +443,7 @@ Double_t AliVertexerTracks::GetStrLinMinDist(Double_t *p0,Double_t *p1,Double_t
 void AliVertexerTracks::OneTrackVertFinder() 
 {
   // find vertex for events with 1 track, using DCA to nominal beam axis
-  if(fDebug) printf("Number of prepared tracks =%d - Call OneTrackVertFinder",fTrkArraySel.GetEntries());
+  AliDebug(1,Form("Number of prepared tracks =%d - Call OneTrackVertFinder",fTrkArraySel.GetEntries()));
   AliExternalTrackParam *track1;
   track1 = (AliExternalTrackParam*)fTrkArraySel.At(0);
   Double_t alpha=track1->GetAlpha();
@@ -533,7 +556,7 @@ Int_t AliVertexerTracks::PrepareTracks(TObjArray &trkArrayOrig,
 //
 // Propagate tracks to initial vertex position and store them in a TObjArray
 //
-  if(fDebug) printf(" PrepareTracks()\n");
+  AliDebug(1," PrepareTracks()");
 
   Int_t nTrksOrig = (Int_t)trkArrayOrig.GetEntriesFast();
   Int_t nTrksSel = 0;
@@ -551,28 +574,39 @@ Int_t AliVertexerTracks::PrepareTracks(TObjArray &trkArrayOrig,
 
   // loop on tracks
   for(Int_t i=0; i<nTrksOrig; i++) {
-    track = new AliExternalTrackParam(*(AliExternalTrackParam*)trkArrayOrig.At(i));
+    AliExternalTrackParam *trackOrig=(AliExternalTrackParam*)trkArrayOrig.At(i);
+    if(trackOrig->Charge()!=0) { // normal tracks
+      track = new AliExternalTrackParam(*(AliExternalTrackParam*)trkArrayOrig.At(i));
+    } else { // neutral tracks (from a V0)
+      track = new AliNeutralTrackParam(*(AliNeutralTrackParam*)trkArrayOrig.At(i));
+    }
     // tgl cut
     if(TMath::Abs(track->GetTgl())>fMaxTgl) {
-      if(fDebug) printf(" rejecting track with tgl = %f\n",track->GetTgl());
+      AliDebug(1,Form(" rejecting track with tgl = %f",track->GetTgl()));
       delete track; continue;
     }
 
+    Bool_t propagateOK = kFALSE;
     // propagate track to vertex
     if(optImpParCut<2 || fOnlyFitter) { // optImpParCut==1 or 0
-      track->PropagateToDCA(initVertex,GetFieldkG(),100.,d0z0,covd0z0);
+      propagateOK = track->PropagateToDCA(initVertex,GetFieldkG(),100.,d0z0,covd0z0);
     } else {              // optImpParCut==2
       fCurrentVertex->GetSigmaXYZ(sigmaCurr);
       normdistx = TMath::Abs(fCurrentVertex->GetXv()-fNominalPos[0])/TMath::Sqrt(sigmaCurr[0]*sigmaCurr[0]+fNominalCov[0]);
       normdisty = TMath::Abs(fCurrentVertex->GetYv()-fNominalPos[1])/TMath::Sqrt(sigmaCurr[1]*sigmaCurr[1]+fNominalCov[2]);
       if(normdistx < 3. && normdisty < 3. &&
         (sigmaCurr[0]+sigmaCurr[1])<(TMath::Sqrt(fNominalCov[0])+TMath::Sqrt(fNominalCov[2]))) {
-       track->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
+       propagateOK = track->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
       } else {
-       track->PropagateToDCA(initVertex,GetFieldkG(),100.,d0z0,covd0z0);
+       propagateOK = track->PropagateToDCA(initVertex,GetFieldkG(),100.,d0z0,covd0z0);
       }
     }
 
+    if(!propagateOK) { 
+      AliDebug(1,"     rejected");
+      delete track; continue; 
+    }
+
     sigmad0 = TMath::Sqrt(covd0z0[0]);
     maxd0rphi = fNSigma*sigmad0;
     if(optImpParCut==1) maxd0rphi *= 5.;
@@ -580,20 +614,20 @@ Int_t AliVertexerTracks::PrepareTracks(TObjArray &trkArrayOrig,
     //sigmad0z0 = TMath::Sqrt(covd0z0[0]+covd0z0[2]); // for future improvement
     //maxd0z0 = 10.*fNSigma*sigmad0z0;
 
-    if(fDebug) printf("trk %d; id %d; |d0| = %f;  d0 cut = %f; |z0| = %f; |d0|oplus|z0| = %f; d0z0 cut = %f\n",i,(Int_t)idOrig[i],TMath::Abs(d0z0[0]),maxd0rphi,TMath::Abs(d0z0[1]),TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1]),fMaxd0z0);
+    AliDebug(1,Form("trk %d; id %d; |d0| = %f;  d0 cut = %f; |z0| = %f; |d0|oplus|z0| = %f; d0z0 cut = %f",i,(Int_t)idOrig[i],TMath::Abs(d0z0[0]),maxd0rphi,TMath::Abs(d0z0[1]),TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1]),fMaxd0z0));
 
 
     //---- track selection based on impact parameters ----//
 
     // always reject tracks outside fiducial volume
     if(TMath::Abs(d0z0[0])>fFiducialR || TMath::Abs(d0z0[1])>fFiducialZ) { 
-      if(fDebug) printf("     rejected\n");
+      AliDebug(1,"     rejected");
       delete track; continue; 
     }
 
     // during iterations 1 and 2 , reject tracks with d0rphi > maxd0rphi
     if(optImpParCut>0 && TMath::Abs(d0z0[0])>maxd0rphi) { 
-      if(fDebug) printf("     rejected\n");
+      AliDebug(1,"     rejected");
       delete track; continue; 
     }
 
@@ -604,7 +638,7 @@ Int_t AliVertexerTracks::PrepareTracks(TObjArray &trkArrayOrig,
        ( fConstraint && optImpParCut==2)) {
       if(nTrksOrig>=3 && 
         TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1])>fMaxd0z0) { 
-       if(fDebug) printf("     rejected\n");
+       AliDebug(1,"     rejected");
        delete track; continue; 
       }
     }
@@ -654,7 +688,7 @@ Bool_t AliVertexerTracks::PropagateTrackTo(AliExternalTrackParam *track,
     //
     Double_t ca=TMath::Cos(alphan-track->GetAlpha()), 
       sa=TMath::Sin(alphan-track->GetAlpha());
-    Double_t sf=track->GetSnp(), cf=TMath::Sqrt(1.- sf*sf);
+    Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf));
     Double_t sinNew =  sf*ca - cf*sa;
     if(TMath::Abs(sinNew) >= maxSnp) return kFALSE;
     if(!track->Rotate(alphan)) return kFALSE;
@@ -674,11 +708,11 @@ AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
 //
 
   if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraint")) {
-    printf("ERROR: primary vertex has no constraint: cannot remove tracks\n");
+    printf("ERROR: primary vertex has no constraint: cannot remove tracks");
     return 0x0;
   }
   if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
-    printf("WARNING: result of tracks' removal will be only approximately correct\n");
+    printf("WARNING: result of tracks' removal will be only approximately correct");
 
   TMatrixD rv(3,1);
   rv(0,0) = inVtx->GetXv();
@@ -705,7 +739,7 @@ AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
   for(Int_t i=0;i<ntrks;i++) {
     track = (AliExternalTrackParam*)trkArray->At(i);
     if(!inVtx->UsesTrack(id[i])) {
-      printf("track %d was not used in vertex fit\n",id[i]);
+      printf("track %d was not used in vertex fit",id[i]);
       continue;
     }
     Double_t alpha = track->GetAlpha();
@@ -735,7 +769,7 @@ AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
 
     nUsedTrks--;
     if(nUsedTrks<2) {
-      printf("Trying to remove too many tracks!\n");
+      printf("Trying to remove too many tracks!");
       return 0x0;
     }
   }
@@ -779,14 +813,14 @@ AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
   outVtx->SetIndices(nIndices,outindices);
   delete [] outindices;
 
-  if(fDebug) {
-    printf("Vertex before removing tracks:\n");
+  /*
+    printf("Vertex before removing tracks:");
     inVtx->PrintStatus();
     inVtx->PrintIndices();
-    printf("Vertex after removing tracks:\n");
+    printf("Vertex after removing tracks:");
     outVtx->PrintStatus();
     outVtx->PrintIndices();
-  }
+  */
 
   return outVtx;
 }
@@ -799,12 +833,15 @@ void AliVertexerTracks::SetCuts(Double_t *cuts)
   SetDCAcut(cuts[0]);
   SetDCAcutIter0(cuts[1]);
   SetMaxd0z0(cuts[2]);
-  SetMinClusters((Int_t)(cuts[3]));
+  if(fMode==0 && cuts[3]<0) SetITSrefitNotRequired();
+  SetMinClusters((Int_t)(TMath::Abs(cuts[3])));
   SetMinTracks((Int_t)(cuts[4]));
   SetNSigmad0(cuts[5]);
   SetMinDetFitter(cuts[6]);
   SetMaxTgl(cuts[7]);
   SetFiducialRZ(cuts[8],cuts[9]);
+  fAlgo=(Int_t)(cuts[10]);
+  fAlgoIter0=(Int_t)(cuts[11]);
 
   return;
 }
@@ -818,22 +855,30 @@ void AliVertexerTracks::SetITSMode(Double_t dcacut,
                                   Double_t mindetfitter,
                                   Double_t maxtgl,
                                   Double_t fidR,
-                                  Double_t fidZ) 
+                                  Double_t fidZ,
+                                  Int_t finderAlgo,
+                                  Int_t finderAlgoIter0)
 {
 //
 //  Cut values for ITS mode
 //
   fMode = 0;
-  SetITSrefitRequired();
+  if(minCls>0) {
+    SetITSrefitRequired();
+  } else {
+    SetITSrefitNotRequired();
+  }
   SetDCAcut(dcacut);
   SetDCAcutIter0(dcacutIter0);
   SetMaxd0z0(maxd0z0);
-  SetMinClusters(minCls);
+  SetMinClusters(TMath::Abs(minCls));
   SetMinTracks(mintrks);
   SetNSigmad0(nsigma);
   SetMinDetFitter(mindetfitter);
   SetMaxTgl(maxtgl);
   SetFiducialRZ(fidR,fidZ);
+  fAlgo=finderAlgo;
+  fAlgoIter0=finderAlgoIter0;
 
   return; 
 }
@@ -847,7 +892,9 @@ void AliVertexerTracks::SetTPCMode(Double_t dcacut,
                                   Double_t mindetfitter,
                                   Double_t maxtgl,
                                   Double_t fidR,
-                                  Double_t fidZ) 
+                                  Double_t fidZ,
+                                  Int_t finderAlgo,
+                                  Int_t finderAlgoIter0) 
 {
 //
 //  Cut values for TPC mode
@@ -863,6 +910,8 @@ void AliVertexerTracks::SetTPCMode(Double_t dcacut,
   SetMinDetFitter(mindetfitter);
   SetMaxTgl(maxtgl);
   SetFiducialRZ(fidR,fidZ);
+  fAlgo=finderAlgo;
+  fAlgoIter0=finderAlgoIter0;
 
   return; 
 }
@@ -906,7 +955,8 @@ void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
     sigmasq[2]=track1->GetSigmaZ2();
     TMatrixD ri(3,1);
     TMatrixD wWi(3,3);
-    if(!TrackToPoint(track1,ri,wWi)) continue;
+    //if(!TrackToPoint(track1,ri,wWi)) {printf("WARNING\n");continue;}
+    if(!TrackToPoint(track1,ri,wWi)) {optUseWeights=kFALSE;printf("WARNING\n");}
     Double_t wmat[9];
     Int_t iel=0;
     for(Int_t ia=0;ia<3;ia++){
@@ -915,7 +965,7 @@ void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
        iel++;
       }    
     }
-    new(linarray[i]) AliStrLine(pos,sigmasq,wmat,dir);     
+    new(linarray[i]) AliStrLine(pos,sigmasq,wmat,dir);
   }
   fVert=TrackletVertexFinder(&linarray,optUseWeights);
   linarray.Clear("C");
@@ -924,7 +974,6 @@ void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
 AliESDVertex AliVertexerTracks::TrackletVertexFinder(TClonesArray *lines, Int_t optUseWeights)
 {
   // Calculate the point at minimum distance to prepared tracks 
-  
   const Int_t knacc = (Int_t)lines->GetEntriesFast();
   Double_t initPos[3]={0.,0.,0.};
 
@@ -945,6 +994,7 @@ AliESDVertex AliVertexerTracks::TrackletVertexFinder(TClonesArray *lines, Int_t
     AliStrLine* line1 = (AliStrLine*)lines->At(i); 
     Double_t p0[3],cd[3],sigmasq[3];
     Double_t wmat[9];
+    if(!line1) printf("ERROR %d %d\n",i,knacc);
     line1->GetP0(p0);
     line1->GetCd(cd);
     line1->GetSigma2P0(sigmasq);
@@ -1060,9 +1110,9 @@ Bool_t AliVertexerTracks::TrackToPoint(AliExternalTrackParam *t,
     uUi(0,1) = t->GetSigmaZY();
     uUi(1,0) = t->GetSigmaZY();
     uUi(1,1) = t->GetSigmaZ2();
-    //printf(" Ui :\n");
-    //printf(" %f   %f\n",uUi(0,0),uUi(0,1));
-    //printf(" %f   %f\n",uUi(1,0),uUi(1,1));
+    //printf(" Ui :");
+    //printf(" %f   %f",uUi(0,0),uUi(0,1));
+    //printf(" %f   %f",uUi(1,0),uUi(1,1));
 
     if(uUi.Determinant() <= 0.) return kFALSE;
     TMatrixD uUiInv(TMatrixD::kInverted,uUi);
@@ -1106,7 +1156,7 @@ Bool_t AliVertexerTracks::TrackToPoint(AliExternalTrackParam *t,
     // covariance matrix of local (x,y,z) - inverted
     TMatrixD uUi(3,3);
     uUi(0,0) = covfVertalongt * fnSigmaForUi00 * fnSigmaForUi00;
-    if(fDebug) printf("=====> sqrtUi00 cm  %f\n",TMath::Sqrt(uUi(0,0)));
+    AliDebug(1,Form("=====> sqrtUi00 cm  %f",TMath::Sqrt(uUi(0,0))));
     uUi(0,1) = 0.;
     uUi(0,2) = 0.;
     uUi(1,0) = 0.;
@@ -1138,7 +1188,7 @@ void AliVertexerTracks::TooFewTracks()
 // When the number of tracks is < fMinTracks,
 // deal with vertices not found and prepare to exit
 //
-  if(fDebug) printf("TooFewTracks\n");
+  AliDebug(1,"TooFewTracks");
 
   Double_t pos[3],err[3];
   pos[0] = fNominalPos[0];
@@ -1148,7 +1198,7 @@ void AliVertexerTracks::TooFewTracks()
   pos[2] = fNominalPos[2];
   err[2] = TMath::Sqrt(fNominalCov[5]);
   Int_t    ncontr = (err[0]>1. ? -1 : -3);
-  fCurrentVertex = 0;
+  if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
   fCurrentVertex = new AliESDVertex(pos,err);
   fCurrentVertex->SetNContributors(ncontr);
 
@@ -1283,13 +1333,11 @@ void AliVertexerTracks::VertexFitter()
 
   Int_t nTrksSel = (Int_t)fTrkArraySel.GetEntries();
   if(nTrksSel==1) useConstraint=kTRUE;
-  if(fDebug) { 
-    printf("--- VertexFitter(): start\n");
-    printf(" Number of tracks in array: %d\n",nTrksSel);
-    printf(" Minimum # tracks required in fit: %d\n",fMinTracks);
-    printf(" Vertex position after finder: %f,%f,%f\n",initPos[0],initPos[1],initPos[2]);
-    if(useConstraint) printf(" This vertex will be used in fit: (%f+-%f,%f+-%f)\n",fNominalPos[0],TMath::Sqrt(fNominalCov[0]),fNominalPos[1],TMath::Sqrt(fNominalCov[2])); 
-  }
+  AliDebug(1,"--- VertexFitter(): start");
+  AliDebug(1,Form(" Number of tracks in array: %d\n",nTrksSel));
+  AliDebug(1,Form(" Minimum # tracks required in fit: %d\n",fMinTracks));
+  AliDebug(1,Form(" Vertex position after finder: %f,%f,%f\n",initPos[0],initPos[1],initPos[2]));
+  if(useConstraint) AliDebug(1,Form(" This vertex will be used in fit: (%f+-%f,%f+-%f)\n",fNominalPos[0],TMath::Sqrt(fNominalCov[0]),fNominalPos[1],TMath::Sqrt(fNominalCov[2]))); 
 
   // special treatment for few-tracks fits (e.g. secondary vertices)
   Bool_t uUi3by3 = kFALSE; if(nTrksSel<5 && !useConstraint) uUi3by3 = kTRUE;
@@ -1329,7 +1377,7 @@ void AliVertexerTracks::VertexFitter()
   // 1st - estimate of vtx using all tracks
   // 2nd - estimate of global chi2
   for(step=0; step<2; step++) {
-    if(fDebug) printf(" step = %d\n",step);
+    AliDebug(1,Form(" step = %d\n",step));
     chi2 = 0.;
     nTrksUsed = 0;
 
@@ -1398,7 +1446,7 @@ void AliVertexerTracks::VertexFitter()
 
     Double_t determinant = sumWi.Determinant();
     if(determinant < fMinDetFitter)  { 
-      printf("det(V) = %f (<%f)\n",determinant,fMinDetFitter);       
+      AliDebug(1,Form("det(V) = %f (<%f)\n",determinant,fMinDetFitter));       
       failed=1;
       continue;
     }
@@ -1432,13 +1480,13 @@ void AliVertexerTracks::VertexFitter()
   // for correct chi2/ndf, count constraint as additional "track"
   if(fConstraint) nTrksUsed++;
   // store data in the vertex object
+  if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
   fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);
 
-  if(fDebug) {
-    printf(" Vertex after fit:\n");
-    fCurrentVertex->PrintStatus();
-    printf("--- VertexFitter(): finish\n");
-  }
+  AliDebug(1," Vertex after fit:");
+  AliDebug(1,Form("xyz: %f %f %f; nc %d",fCurrentVertex->GetXv(),fCurrentVertex->GetYv(),fCurrentVertex->GetZv(),fCurrentVertex->GetNContributors()));
+  AliDebug(1,"--- VertexFitter(): finish\n");
 
   return;
 }
@@ -1451,7 +1499,7 @@ AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray,
 //
 // Return vertex from tracks (AliExternalTrackParam) in array
 //
-
+  fCurrentVertex = 0;
   SetConstraintOff();
 
   // get tracks and propagate them to initial vertex position
@@ -1471,7 +1519,7 @@ AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray,
     case 5: VertexFinder(0);              break;
     default: printf("Wrong algorithm\n"); break;  
   }
-  if(fDebug) printf(" Vertex finding completed\n");
+  AliDebug(1," Vertex finding completed\n");
 
   // vertex fitter
   if(optUseFitter) {
@@ -1499,7 +1547,7 @@ AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray,
       if(optPropagate && optUseFitter) {
        if(TMath::Sqrt(fCurrentVertex->GetXv()*fCurrentVertex->GetXv()+fCurrentVertex->GetYv()*fCurrentVertex->GetYv())<3.) {
          t->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
-         if(fDebug) printf("Track %d propagated to found vertex\n",jj);
+         AliDebug(1,Form("Track %d propagated to found vertex",jj));
        } else {
          AliWarning("Found vertex outside beam pipe!");
        }