STEER/ESD/AliVertexerTracks.cxx

   1 /**************************************************************************
   2  * Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. *
   3  *                                                                        *
   4  * Author: The ALICE Off-line Project.                                    *
   5  * Contributors are mentioned in the code where appropriate.              *
   6  *                                                                        *
   7  * Permission to use, copy, modify and distribute this software and its   *
   8  * documentation strictly for non-commercial purposes is hereby granted   *
   9  * without fee, provided that the above copyright notice appears in all   *
  10  * copies and that both the copyright notice and this permission notice   *
  11  * appear in the supporting documentation. The authors make no claims     *
  12  * about the suitability of this software for any purpose. It is          *
  13  * provided "as is" without express or implied warranty.                  *
  14  **************************************************************************/
  15
  16 //-----------------------------------------------------------------
  17 //    Implementation of the vertexer from tracks
  18 //
  19 // Origin: AliITSVertexerTracks
  20 //         A.Dainese, Padova,
  21 //         andrea.dainese@pd.infn.it
  22 //         M.Masera,  Torino,
  23 //         massimo.masera@to.infn.it
  24 // Moved to STEER and adapted to ESD tracks:
  25 //         F.Prino,  Torino, prino@to.infn.it
  26 //-----------------------------------------------------------------
  27
  28 //---- Root headers --------
  29 #include <TSystem.h>
  30 #include <TClonesArray.h>
  31 #include <TDirectory.h>
  32 #include <TFile.h>
  33 //---- AliRoot headers -----
  34 #include "AliStrLine.h"
  35 #include "AliExternalTrackParam.h"
  36 #include "AliNeutralTrackParam.h"
  37 #include "AliVEvent.h"
  38 #include "AliVTrack.h"
  39 #include "AliESDtrack.h"
  40 #include "AliESDEvent.h"
  41 #include "AliVertexerTracks.h"
  42
  43
  44 ClassImp(AliVertexerTracks)
  45
  46
  47 //----------------------------------------------------------------------------
  48 AliVertexerTracks::AliVertexerTracks():
  49 TObject(),
  50 fVert(),
  51 fCurrentVertex(0),
  52 fMode(0),
  53 fFieldkG(-999.),
  54 fTrkArraySel(),
  55 fIdSel(0),
  56 fTrksToSkip(0),
  57 fNTrksToSkip(0),
  58 fConstraint(kFALSE),
  59 fOnlyFitter(kFALSE),
  60 fMinTracks(1),
  61 fMinClusters(5),
  62 fDCAcut(0.1),
  63 fDCAcutIter0(0.1),
  64 fNSigma(3.),
  65 fMaxd0z0(0.5),
  66 fMinDetFitter(100.),
  67 fMaxTgl(1000.),
  68 fITSrefit(kTRUE),
  69 fITSpureSA(kFALSE),
  70 fFiducialR(3.),
  71 fFiducialZ(30.),
  72 fnSigmaForUi00(1.5),
  73 fAlgo(1),
  74 fAlgoIter0(4),
  75 fSelectOnTOFBunchCrossing(kFALSE),
  76 fKeepAlsoUnflaggedTOFBunchCrossing(kTRUE),
  77 fMVWSum(0),
  78 fMVWE2(0),
  79 fMVTukey2(6.),
  80 fMVSigma2(1.),
  81 fMVSig2Ini(1.e3),
  82 fMVMaxSigma2(3.),
  83 fMVMinSig2Red(0.005),
  84 fMVMinDst(10.e-4),
  85 fMVScanStep(3.),
  86 fMVMaxWghNtr(10.),
  87 fMVFinalWBinary(kTRUE),
  88 fBCSpacing(50),
  89 fMVVertices(0)
  90 {
  91 //
  92 // Default constructor
  93 //
  94   SetVtxStart();
  95   SetVtxStartSigma();
  96 }
  97 //----------------------------------------------------------------------------
  98 AliVertexerTracks::AliVertexerTracks(Double_t fieldkG):
  99 TObject(),
 100 fVert(),
 101 fCurrentVertex(0),
 102 fMode(0),
 103 fFieldkG(fieldkG),
 104 fTrkArraySel(),
 105 fIdSel(0),
 106 fTrksToSkip(0),
 107 fNTrksToSkip(0),
 108 fConstraint(kFALSE),
 109 fOnlyFitter(kFALSE),
 110 fMinTracks(1),
 111 fMinClusters(5),
 112 fDCAcut(0.1),
 113 fDCAcutIter0(0.1),
 114 fNSigma(3.),
 115 fMaxd0z0(0.5),
 116 fMinDetFitter(100.),
 117 fMaxTgl(1000.),
 118 fITSrefit(kTRUE),
 119 fITSpureSA(kFALSE),
 120 fFiducialR(3.),
 121 fFiducialZ(30.),
 122 fnSigmaForUi00(1.5),
 123 fAlgo(1),
 124 fAlgoIter0(4),
 125 fSelectOnTOFBunchCrossing(kFALSE),
 126 fKeepAlsoUnflaggedTOFBunchCrossing(kTRUE),
 127 fMVWSum(0),
 128 fMVWE2(0),
 129 fMVTukey2(6.),
 130 fMVSigma2(1.),
 131 fMVSig2Ini(1.e3),
 132 fMVMaxSigma2(3.),
 133 fMVMinSig2Red(0.005),
 134 fMVMinDst(10.e-4),
 135 fMVScanStep(3.),
 136 fMVMaxWghNtr(10.),
 137 fMVFinalWBinary(kTRUE),
 138 fBCSpacing(50),
 139 fMVVertices(0)
 140 {
 141 //
 142 // Standard constructor
 143 //
 144   SetVtxStart();
 145   SetVtxStartSigma();
 146 }
 147 //-----------------------------------------------------------------------------
 148 AliVertexerTracks::~AliVertexerTracks()
 149 {
 150   // Default Destructor
 151   // The objects pointed by the following pointer are not owned
 152   // by this class and are not deleted
 153   fCurrentVertex = 0;
 154   if(fTrksToSkip) { delete [] fTrksToSkip; fTrksToSkip=NULL; fNTrksToSkip=0;}
 155   if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
 156   if(fMVVertices) delete fMVVertices;
 157 }
 158
 159 //----------------------------------------------------------------------------
 160 AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliVEvent *vEvent)
 161 {
 162 //
 163 // Primary vertex for current ESD or AOD event
 164 // (Two iterations:
 165 //  1st with 5*fNSigma*sigma cut w.r.t. to initial vertex
 166 //      + cut on sqrt(d0d0+z0z0) if fConstraint=kFALSE
 167 //  2nd with fNSigma*sigma cut w.r.t. to vertex found in 1st iteration)
 168 //
 169   fCurrentVertex = 0;
 170   TString evtype = vEvent->IsA()->GetName();
 171   Bool_t inputAOD = ((evtype=="AliAODEvent") ? kTRUE : kFALSE);
 172
 173   if(inputAOD && fMode==1) {
 174     printf("Error : AliVertexerTracks: no TPC-only vertex from AOD\n");
 175     TooFewTracks();
 176     return fCurrentVertex;
 177   }
 178
 179   // accept 1-track case only if constraint is available
 180   if(!fConstraint && fMinTracks==1) fMinTracks=2;
 181
 182   // read tracks from AlivEvent
 183   Int_t nTrks = (Int_t)vEvent->GetNumberOfTracks();
 184   if(nTrks<fMinTracks) {
 185     TooFewTracks();
 186     return fCurrentVertex;
 187   }
 188   //
 189   int bcRound = fBCSpacing/25;   // profit from larger than 25ns spacing and set correct BC
 190   TDirectory * olddir = gDirectory;
 191   TFile *f = 0;
 192   if(nTrks>500) f = new TFile("VertexerTracks.root","recreate");
 193   TObjArray trkArrayOrig(nTrks);
 194   UShort_t *idOrig = new UShort_t[nTrks];
 195
 196   Int_t nTrksOrig=0;
 197   AliExternalTrackParam *t=0;
 198   // loop on tracks
 199   for(Int_t i=0; i<nTrks; i++) {
 200     AliVTrack *track = (AliVTrack*)vEvent->GetTrack(i);
 201     // check tracks to skip
 202     Bool_t skipThis = kFALSE;
 203     for(Int_t j=0; j<fNTrksToSkip; j++) {
 204       if(track->GetID()==fTrksToSkip[j]) {
 205         AliDebug(1,Form("skipping track: %d",i));
 206         skipThis = kTRUE;
 207       }
 208     }
 209     if(skipThis) continue;
 210
 211     // skip pure ITS SA tracks (default)
 212     if(!fITSpureSA && (track->GetStatus()&AliESDtrack::kITSpureSA)) continue;
 213     // or use only pure ITS SA tracks
 214     if(fITSpureSA && !(track->GetStatus()&AliESDtrack::kITSpureSA)) continue;
 215
 216     // kITSrefit
 217     if(fMode==0 && fITSrefit && !(track->GetStatus()&AliESDtrack::kITSrefit)) continue;
 218
 219     if(!inputAOD) {  // ESD
 220       AliESDtrack* esdt = (AliESDtrack*)track;
 221       if(esdt->GetNcls(fMode) < fMinClusters) continue;
 222       if(fMode==0) {        // ITS mode
 223         Double_t x,p[5],cov[15];
 224         esdt->GetExternalParameters(x,p);
 225         esdt->GetExternalCovariance(cov);
 226         t = new AliExternalTrackParam(x,esdt->GetAlpha(),p,cov);
 227       } else if(fMode==1) { // TPC mode
 228         t = (AliExternalTrackParam*)esdt->GetTPCInnerParam();
 229         if(!t) continue;
 230         Double_t radius = 2.8; //something less than the beam pipe radius
 231         if(!PropagateTrackTo(t,radius)) continue;
 232       }
 233     } else {          // AOD (only ITS mode)
 234       Int_t ncls0=0;
 235       for(Int_t l=0;l<6;l++) if(TESTBIT(track->GetITSClusterMap(),l)) ncls0++;
 236       if(ncls0 < fMinClusters) continue;
 237       t = new AliExternalTrackParam(track);
 238     }
 239
 240     // use TOF info about bunch crossing
 241     if(fSelectOnTOFBunchCrossing) {
 242       int bc = track->GetTOFBunchCrossing(fFieldkG);
 243       if (bc>AliVTrack::kTOFBCNA) bc /= bcRound;
 244       t->SetUniqueID(UInt_t(bc + kTOFBCShift));
 245     }
 246     //
 247     trkArrayOrig.AddLast(t);
 248     idOrig[nTrksOrig]=(UShort_t)track->GetID();
 249     nTrksOrig++;
 250   } // end loop on tracks
 251
 252   // call method that will reconstruct the vertex
 253   FindPrimaryVertex(&trkArrayOrig,idOrig);
 254
 255   if(!inputAOD) AnalyzePileUp((AliESDEvent*)vEvent);
 256
 257   if(fMode==0) trkArrayOrig.Delete();
 258   delete [] idOrig; idOrig=NULL;
 259
 260   if(f) {
 261     f->Close(); delete f; f = NULL;
 262     gSystem->Unlink("VertexerTracks.root");
 263     olddir->cd();
 264   }
 265
 266   // set vertex ID for tracks used in the fit
 267   // (only for ESD)
 268   if(!inputAOD && fCurrentVertex) {
 269     Int_t nIndices = fCurrentVertex->GetNIndices();
 270     UShort_t *indices = fCurrentVertex->GetIndices();
 271     for(Int_t ind=0; ind<nIndices; ind++) {
 272       AliESDtrack *esdt = (AliESDtrack*)vEvent->GetTrack(indices[ind]);
 273       esdt->SetVertexID(-1);
 274     }
 275   }
 276
 277   return fCurrentVertex;
 278 }
 279 //----------------------------------------------------------------------------
 280 AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const TObjArray *trkArrayOrig,
 281                                                    UShort_t *idOrig)
 282 {
 283 //
 284 // Primary vertex using the AliExternalTrackParam's in the TObjArray.
 285 // idOrig must contain the track IDs from AliESDtrack::GetID()
 286 // (Two iterations:
 287 //  1st with 5*fNSigma*sigma cut w.r.t. to initial vertex
 288 //      + cut on sqrt(d0d0+z0z0) if fConstraint=kFALSE
 289 //  2nd with fNSigma*sigma cut w.r.t. to vertex found in 1st iteration)
 290 //
 291   fCurrentVertex = 0;
 292   // accept 1-track case only if constraint is available
 293   if(!fConstraint && fMinTracks==1) fMinTracks=2;
 294
 295   // read tracks from array
 296   Int_t nTrksOrig = (Int_t)trkArrayOrig->GetEntriesFast();
 297   AliDebug(1,Form("Initial number of tracks: %d",nTrksOrig));
 298   if(nTrksOrig<fMinTracks) {
 299     AliDebug(1,"TooFewTracks");
 300     TooFewTracks();
 301     return fCurrentVertex;
 302   }
 303
 304   // If fConstraint=kFALSE
 305   // run VertexFinder(1) to get rough estimate of initVertex (x,y)
 306   if(!fConstraint) {
 307     // fill fTrkArraySel, for VertexFinder()
 308     fIdSel = new UShort_t[nTrksOrig];
 309     PrepareTracks(*trkArrayOrig,idOrig,0);
 310     if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
 311     Double_t cutsave = fDCAcut;
 312     fDCAcut = fDCAcutIter0;
 313     // vertex finder
 314     switch (fAlgoIter0) {
 315     case 1: StrLinVertexFinderMinDist(1); break;
 316     case 2: StrLinVertexFinderMinDist(0); break;
 317     case 3: HelixVertexFinder();          break;
 318     case 4: VertexFinder(1);              break;
 319     case 5: VertexFinder(0);              break;
 320     default: printf("Wrong algorithm\n"); break;
 321     }
 322     fDCAcut = cutsave;
 323     if(fVert.GetNContributors()>0) {
 324       fVert.GetXYZ(fNominalPos);
 325       fNominalPos[0] = fVert.GetXv();
 326       fNominalPos[1] = fVert.GetYv();
 327       fNominalPos[2] = fVert.GetZv();
 328       AliDebug(1,Form("No mean vertex: VertexFinder gives (%f, %f, %f)",fNominalPos[0],fNominalPos[1],fNominalPos[2]));
 329     } else {
 330       fNominalPos[0] = 0.;
 331       fNominalPos[1] = 0.;
 332       fNominalPos[2] = 0.;
 333       AliDebug(1,"No mean vertex and VertexFinder failed");
 334     }
 335   }
 336
 337   // TWO ITERATIONS:
 338   //
 339   // ITERATION 1
 340   // propagate tracks to fNominalPos vertex
 341   // preselect them:
 342   // if(constraint) reject for |d0|>5*fNSigma*sigma w.r.t. fNominal... vertex
 343   // else  reject for |d0|\oplus|z0| > 5 mm w.r.t. fNominal... vertex
 344   // ITERATION 2
 345   // propagate tracks to best between initVertex and fCurrentVertex
 346   // preselect tracks (reject for |d0|>fNSigma*sigma w.r.t. best
 347   //                   between initVertex and fCurrentVertex)
 348   Bool_t multiMode = kFALSE;
 349   for(Int_t iter=1; iter<=2; iter++) {
 350     if (fAlgo==kMultiVertexer && iter==2) break; // multivertexer does not need 2 iterations
 351     if(fOnlyFitter && iter==1) continue;
 352     if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
 353     fIdSel = new UShort_t[nTrksOrig];
 354     Int_t nTrksSel = PrepareTracks(*trkArrayOrig,idOrig,iter);
 355     AliDebug(1,Form("N tracks selected in iteration %d: %d",iter,nTrksSel));
 356     if(nTrksSel < fMinTracks) {
 357       TooFewTracks();
 358       return fCurrentVertex;
 359     }
 360
 361     // vertex finder
 362     if(!fOnlyFitter) {
 363       if(nTrksSel==1 && fAlgo!=kMultiVertexer) {
 364         AliDebug(1,"Just one track");
 365         OneTrackVertFinder();
 366       } else {
 367         switch (fAlgo) {
 368         case kStrLinVertexFinderMinDist1: StrLinVertexFinderMinDist(1); break;
 369         case kStrLinVertexFinderMinDist0: StrLinVertexFinderMinDist(0); break;
 370         case kHelixVertexFinder         : HelixVertexFinder();          break;
 371         case kVertexFinder1             : VertexFinder(1);              break;
 372         case kVertexFinder0             : VertexFinder(0);              break;
 373         case kMultiVertexer             : FindVerticesMV(); multiMode = kTRUE; break;
 374         default: printf("Wrong algorithm"); break;
 375         }
 376       }
 377       AliDebug(1," Vertex finding completed");
 378     }
 379     if (multiMode) break; // // multivertexer does not need 2nd iteration
 380     // vertex fitter
 381     VertexFitter();
 382   } // end loop on the two iterations
 383
 384   if (!multiMode || fMVVertices->GetEntries()==0) { // in multi-vertex mode this is already done for found vertices
 385     // set indices of used tracks
 386     UShort_t *indices = 0;
 387     if(fCurrentVertex->GetNContributors()>0) {
 388       Int_t nIndices = (Int_t)fTrkArraySel.GetEntriesFast();
 389       indices = new UShort_t[nIndices];
 390       for(Int_t jj=0; jj<nIndices; jj++)
 391         indices[jj] = fIdSel[jj];
 392       fCurrentVertex->SetIndices(nIndices,indices);
 393     }
 394     if (indices) {delete [] indices; indices=NULL;}
 395     //
 396     // set vertex title
 397     TString title="VertexerTracksNoConstraint";
 398     if(fConstraint) {
 399       title="VertexerTracksWithConstraint";
 400       if(fOnlyFitter) title.Append("OnlyFitter");
 401     }
 402     fCurrentVertex->SetTitle(title.Data());
 403     //
 404     AliDebug(1,Form("xyz: %f %f %f; nc %d",fCurrentVertex->GetXv(),fCurrentVertex->GetYv(),fCurrentVertex->GetZv(),fCurrentVertex->GetNContributors()));
 405   }
 406   // clean up
 407   delete [] fIdSel; fIdSel=NULL;
 408   fTrkArraySel.Delete();
 409   if(fTrksToSkip) { delete [] fTrksToSkip; fTrksToSkip=NULL; }
 410   //
 411
 412   return fCurrentVertex;
 413 }
 414 //------------------------------------------------------------------------
 415 Double_t AliVertexerTracks::GetDeterminant3X3(Double_t matr[][3])
 416 {
 417   //
 418   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];
 419  return det;
 420 }
 421 //-------------------------------------------------------------------------
 422 void AliVertexerTracks::GetStrLinDerivMatrix(const Double_t *p0,const Double_t *p1,Double_t (*m)[3],Double_t *d)
 423 {
 424   //
 425   Double_t x12=p0[0]-p1[0];
 426   Double_t y12=p0[1]-p1[1];
 427   Double_t z12=p0[2]-p1[2];
 428   Double_t kk=x12*x12+y12*y12+z12*z12;
 429   m[0][0]=2-2/kk*x12*x12;
 430   m[0][1]=-2/kk*x12*y12;
 431   m[0][2]=-2/kk*x12*z12;
 432   m[1][0]=-2/kk*x12*y12;
 433   m[1][1]=2-2/kk*y12*y12;
 434   m[1][2]=-2/kk*y12*z12;
 435   m[2][0]=-2/kk*x12*z12;
 436   m[2][1]=-2*y12*z12;
 437   m[2][2]=2-2/kk*z12*z12;
 438   d[0]=2*p0[0]-2/kk*p0[0]*x12*x12-2/kk*p0[2]*x12*z12-2/kk*p0[1]*x12*y12;
 439   d[1]=2*p0[1]-2/kk*p0[1]*y12*y12-2/kk*p0[0]*x12*y12-2/kk*p0[2]*z12*y12;
 440   d[2]=2*p0[2]-2/kk*p0[2]*z12*z12-2/kk*p0[0]*x12*z12-2/kk*p0[1]*z12*y12;
 441
 442 }
 443 //--------------------------------------------------------------------------
 444 void AliVertexerTracks::GetStrLinDerivMatrix(const Double_t *p0,const Double_t *p1,const Double_t *sigmasq,Double_t (*m)[3],Double_t *d)
 445 {
 446   //
 447   Double_t x12=p1[0]-p0[0];
 448   Double_t y12=p1[1]-p0[1];
 449   Double_t z12=p1[2]-p0[2];
 450
 451   Double_t den= x12*x12*sigmasq[1]*sigmasq[2]+y12*y12*sigmasq[0]*sigmasq[2]+z12*z12*sigmasq[0]*sigmasq[1];
 452
 453   Double_t kk= 2*(x12*x12/sigmasq[0]+y12*y12/sigmasq[1]+z12*z12/sigmasq[2]);
 454
 455   Double_t cc[3];
 456   cc[0]=-x12/sigmasq[0];
 457   cc[1]=-y12/sigmasq[1];
 458   cc[2]=-z12/sigmasq[2];
 459
 460   Double_t ww=(-p0[0]*x12*sigmasq[1]*sigmasq[2]-p0[1]*y12*sigmasq[0]*sigmasq[2]-p0[2]*z12*sigmasq[0]*sigmasq[1])/den;
 461
 462   Double_t ss= -p0[0]*cc[0]-p0[1]*cc[1]-p0[2]*cc[2];
 463
 464   Double_t aa[3];
 465   aa[0]=x12*sigmasq[1]*sigmasq[2]/den;
 466   aa[1]=y12*sigmasq[0]*sigmasq[2]/den;
 467   aa[2]=z12*sigmasq[0]*sigmasq[1]/den;
 468
 469   m[0][0]=aa[0]*(aa[0]*kk+2*cc[0])+2*cc[0]*aa[0]+2/sigmasq[0];
 470   m[0][1]=aa[1]*(aa[0]*kk+2*cc[0])+2*cc[1]*aa[0];
 471   m[0][2]=aa[2]*(aa[0]*kk+2*cc[0])+2*cc[2]*aa[0];
 472
 473   m[1][0]=aa[0]*(aa[1]*kk+2*cc[1])+2*cc[0]*aa[1];
 474   m[1][1]=aa[1]*(aa[1]*kk+2*cc[1])+2*cc[1]*aa[1]+2/sigmasq[1];
 475   m[1][2]=aa[2]*(aa[1]*kk+2*cc[1])+2*cc[2]*aa[1];
 476
 477   m[2][0]=aa[0]*(aa[2]*kk+2*cc[2])+2*cc[0]*aa[2];
 478   m[2][1]=aa[1]*(aa[2]*kk+2*cc[2])+2*cc[1]*aa[2];
 479   m[2][2]=aa[2]*(aa[2]*kk+2*cc[2])+2*cc[2]*aa[2]+2/sigmasq[2];
 480
 481   d[0]=-ww*(aa[0]*kk+2*cc[0])-2*ss*aa[0]+2*p0[0]/sigmasq[0];
 482   d[1]=-ww*(aa[1]*kk+2*cc[1])-2*ss*aa[1]+2*p0[1]/sigmasq[1];
 483   d[2]=-ww*(aa[2]*kk+2*cc[2])-2*ss*aa[2]+2*p0[2]/sigmasq[2];
 484
 485   }
 486 //--------------------------------------------------------------------------
 487 Double_t AliVertexerTracks::GetStrLinMinDist(const Double_t *p0,const Double_t *p1,const Double_t *x0)
 488 {
 489   //
 490   Double_t x12=p0[0]-p1[0];
 491   Double_t y12=p0[1]-p1[1];
 492   Double_t z12=p0[2]-p1[2];
 493   Double_t x10=p0[0]-x0[0];
 494   Double_t y10=p0[1]-x0[1];
 495   Double_t z10=p0[2]-x0[2];
 496   //  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);
 497   return ((y10*z12-z10*y12)*(y10*z12-z10*y12)+
 498           (z10*x12-x10*z12)*(z10*x12-x10*z12)+
 499           (x10*y12-y10*x12)*(x10*y12-y10*x12))
 500     /(x12*x12+y12*y12+z12*z12);
 501 }
 502 //---------------------------------------------------------------------------
 503 void AliVertexerTracks::OneTrackVertFinder()
 504 {
 505   // find vertex for events with 1 track, using DCA to nominal beam axis
 506   AliDebug(1,Form("Number of prepared tracks =%d - Call OneTrackVertFinder",fTrkArraySel.GetEntries()));
 507   AliExternalTrackParam *track1;
 508   track1 = (AliExternalTrackParam*)fTrkArraySel.At(0);
 509   Double_t alpha=track1->GetAlpha();
 510   Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
 511   Double_t pos[3],dir[3];
 512   track1->GetXYZAt(mindist,GetFieldkG(),pos);
 513   track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
 514   AliStrLine *line1 = new AliStrLine(pos,dir);
 515   Double_t p1[3]={fNominalPos[0],fNominalPos[1],0.};
 516   Double_t p2[3]={fNominalPos[0],fNominalPos[1],10.};
 517   AliStrLine *zeta=new AliStrLine(p1,p2,kTRUE);
 518   Double_t crosspoint[3]={0.,0.,0.};
 519   Double_t sigma=999.;
 520   Int_t nContrib=-1;
 521   Int_t retcode = zeta->Cross(line1,crosspoint);
 522   if(retcode>=0){
 523     sigma=line1->GetDistFromPoint(crosspoint);
 524     nContrib=1;
 525   }
 526   delete zeta;
 527   delete line1;
 528   fVert.SetXYZ(crosspoint);
 529   fVert.SetDispersion(sigma);
 530   fVert.SetNContributors(nContrib);
 531 }
 532 //---------------------------------------------------------------------------
 533 void AliVertexerTracks::HelixVertexFinder()
 534 {
 535   // Get estimate of vertex position in (x,y) from tracks DCA
 536
 537
 538   Double_t initPos[3];
 539   initPos[2] = 0.;
 540   for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
 541
 542   Int_t nacc = (Int_t)fTrkArraySel.GetEntriesFast();
 543
 544   Double_t aver[3]={0.,0.,0.};
 545   Double_t averquad[3]={0.,0.,0.};
 546   Double_t sigmaquad[3]={0.,0.,0.};
 547   Double_t sigma=0;
 548   Int_t ncombi = 0;
 549   AliExternalTrackParam *track1;
 550   AliExternalTrackParam *track2;
 551   Double_t distCA;
 552   Double_t x;
 553   Double_t alpha, cs, sn;
 554   Double_t crosspoint[3];
 555   for(Int_t i=0; i<nacc; i++){
 556     track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
 557
 558
 559     for(Int_t j=i+1; j<nacc; j++){
 560       track2 = (AliExternalTrackParam*)fTrkArraySel.At(j);
 561
 562       distCA=track2->PropagateToDCA(track1,GetFieldkG());
 563       if(fDCAcut<=0 ||(fDCAcut>0&&distCA<fDCAcut)){
 564         x=track1->GetX();
 565         alpha=track1->GetAlpha();
 566         cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
 567         Double_t x1=x*cs - track1->GetY()*sn;
 568         Double_t y1=x*sn + track1->GetY()*cs;
 569         Double_t z1=track1->GetZ();
 570
 571         Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
 572         x=track2->GetX();
 573         alpha=track2->GetAlpha();
 574         cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
 575         Double_t x2=x*cs - track2->GetY()*sn;
 576         Double_t y2=x*sn + track2->GetY()*cs;
 577         Double_t z2=track2->GetZ();
 578         Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
 579         Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
 580         Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
 581         Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
 582         Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
 583         crosspoint[0]=wx1*x1 + wx2*x2;
 584         crosspoint[1]=wy1*y1 + wy2*y2;
 585         crosspoint[2]=wz1*z1 + wz2*z2;
 586
 587         ncombi++;
 588         for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
 589         for(Int_t jj=0;jj<3;jj++)averquad[jj]+=(crosspoint[jj]*crosspoint[jj]);
 590       }
 591     }
 592
 593   }
 594   if(ncombi>0){
 595     for(Int_t jj=0;jj<3;jj++){
 596       initPos[jj] = aver[jj]/ncombi;
 597       averquad[jj]/=ncombi;
 598       sigmaquad[jj]=averquad[jj]-initPos[jj]*initPos[jj];
 599       sigma+=sigmaquad[jj];
 600     }
 601     sigma=TMath::Sqrt(TMath::Abs(sigma));
 602   }
 603   else {
 604     Warning("HelixVertexFinder","Finder did not succed");
 605     sigma=999;
 606   }
 607   fVert.SetXYZ(initPos);
 608   fVert.SetDispersion(sigma);
 609   fVert.SetNContributors(ncombi);
 610 }
 611 //----------------------------------------------------------------------------
 612 Int_t AliVertexerTracks::PrepareTracks(const TObjArray &trkArrayOrig,
 613                                        const UShort_t *idOrig,
 614                                        Int_t optImpParCut)
 615 {
 616 //
 617 // Propagate tracks to initial vertex position and store them in a TObjArray
 618 //
 619   AliDebug(1," PrepareTracks()");
 620
 621   Int_t nTrksOrig = (Int_t)trkArrayOrig.GetEntriesFast();
 622   Int_t nTrksSel = 0;
 623   Double_t maxd0rphi;
 624   Double_t sigmaCurr[3];
 625   Double_t normdistx,normdisty;
 626   Double_t d0z0[2],covd0z0[3];
 627   Double_t sigmad0;
 628
 629   AliESDVertex *initVertex = new AliESDVertex(fNominalPos,fNominalCov,1,1);
 630
 631   if(!fTrkArraySel.IsEmpty()) fTrkArraySel.Delete();
 632
 633   AliExternalTrackParam *track=0;
 634
 635   // loop on tracks
 636   for(Int_t i=0; i<nTrksOrig; i++) {
 637     AliExternalTrackParam *trackOrig=(AliExternalTrackParam*)trkArrayOrig.At(i);
 638     if(trackOrig->Charge()!=0) { // normal tracks
 639       track = new AliExternalTrackParam(*(AliExternalTrackParam*)trkArrayOrig.At(i));
 640     } else { // neutral tracks (from a V0)
 641       track = new AliNeutralTrackParam(*(AliNeutralTrackParam*)trkArrayOrig.At(i));
 642     }
 643     track->SetUniqueID(trackOrig->GetUniqueID());
 644     // tgl cut
 645     if(TMath::Abs(track->GetTgl())>fMaxTgl) {
 646       AliDebug(1,Form(" rejecting track with tgl = %f",track->GetTgl()));
 647       delete track; continue;
 648     }
 649
 650     Bool_t propagateOK = kFALSE, cutond0z0 = kTRUE;
 651     // propagate track to vertex
 652     if(optImpParCut<2 || fOnlyFitter) { // optImpParCut==1 or 0
 653       propagateOK = track->PropagateToDCA(initVertex,GetFieldkG(),100.,d0z0,covd0z0);
 654     } else {              // optImpParCut==2
 655       fCurrentVertex->GetSigmaXYZ(sigmaCurr);
 656       normdistx = TMath::Abs(fCurrentVertex->GetXv()-fNominalPos[0])/TMath::Sqrt(sigmaCurr[0]*sigmaCurr[0]+fNominalCov[0]);
 657       normdisty = TMath::Abs(fCurrentVertex->GetYv()-fNominalPos[1])/TMath::Sqrt(sigmaCurr[1]*sigmaCurr[1]+fNominalCov[2]);
 658       AliDebug(1,Form("normdistx %f  %f    %f",fCurrentVertex->GetXv(),fNominalPos[0],TMath::Sqrt(sigmaCurr[0]*sigmaCurr[0]+fNominalCov[0])));
 659       AliDebug(1,Form("normdisty %f  %f    %f",fCurrentVertex->GetYv(),fNominalPos[1],TMath::Sqrt(sigmaCurr[1]*sigmaCurr[1]+fNominalCov[2])));
 660       AliDebug(1,Form("sigmaCurr %f %f    %f",sigmaCurr[0],sigmaCurr[1],TMath::Sqrt(fNominalCov[0])+TMath::Sqrt(fNominalCov[2])));
 661       if(normdistx < 3. && normdisty < 3. &&
 662          (sigmaCurr[0]+sigmaCurr[1])<(TMath::Sqrt(fNominalCov[0])+TMath::Sqrt(fNominalCov[2]))) {
 663         propagateOK = track->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
 664       } else {
 665         propagateOK = track->PropagateToDCA(initVertex,GetFieldkG(),100.,d0z0,covd0z0);
 666         if(fConstraint) cutond0z0=kFALSE;
 667       }
 668     }
 669
 670     if(!propagateOK) {
 671       AliDebug(1,"     rejected");
 672       delete track; continue;
 673     }
 674
 675     sigmad0 = TMath::Sqrt(covd0z0[0]);
 676     maxd0rphi = fNSigma*sigmad0;
 677     if(optImpParCut==1) maxd0rphi *= 5.;
 678     maxd0rphi = TMath::Min(maxd0rphi,fFiducialR);
 679     //sigmad0z0 = TMath::Sqrt(covd0z0[0]+covd0z0[2]);
 680
 681     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));
 682
 683
 684     //---- track selection based on impact parameters ----//
 685
 686     // always reject tracks outside fiducial volume
 687     if(TMath::Abs(d0z0[0])>fFiducialR || TMath::Abs(d0z0[1])>fFiducialZ) {
 688       AliDebug(1,"     rejected");
 689       delete track; continue;
 690     }
 691
 692     // during iterations 1 and 2 , reject tracks with d0rphi > maxd0rphi
 693     if(optImpParCut>0 && TMath::Abs(d0z0[0])>maxd0rphi) {
 694       AliDebug(1,"     rejected");
 695       delete track; continue;
 696     }
 697
 698     // if fConstraint=kFALSE, during iterations 1 and 2 ||
 699     // if fConstraint=kTRUE, during iteration 2,
 700     // select tracks with d0oplusz0 < fMaxd0z0
 701     if((!fConstraint && optImpParCut>0 && fVert.GetNContributors()>0) ||
 702        ( fConstraint && optImpParCut==2 && cutond0z0)) {
 703       if(nTrksOrig>=3 &&
 704          TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1])>fMaxd0z0) {
 705         AliDebug(1,"     rejected");
 706         delete track; continue;
 707       }
 708     }
 709
 710     // track passed all selections
 711     fTrkArraySel.AddLast(track);
 712     fIdSel[nTrksSel] = idOrig[i];
 713     nTrksSel++;
 714   } // end loop on tracks
 715
 716   delete initVertex;
 717
 718   return nTrksSel;
 719 }
 720 //----------------------------------------------------------------------------
 721 Bool_t AliVertexerTracks::PropagateTrackTo(AliExternalTrackParam *track,
 722                                            Double_t xToGo) {
 723   //----------------------------------------------------------------
 724   // COPIED from AliTracker
 725   //
 726   // Propagates the track to the plane X=xk (cm).
 727   //
 728   //  Origin: Marian Ivanov,  Marian.Ivanov@cern.ch
 729   //----------------------------------------------------------------
 730
 731   const Double_t kEpsilon = 0.00001;
 732   Double_t xpos = track->GetX();
 733   Double_t dir = (xpos<xToGo) ? 1. : -1.;
 734   Double_t maxStep = 5;
 735   Double_t maxSnp = 0.8;
 736   //
 737   while ( (xToGo-xpos)*dir > kEpsilon){
 738     Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
 739     Double_t x    = xpos+step;
 740     Double_t xyz0[3],xyz1[3];
 741     track->GetXYZ(xyz0);   //starting global position
 742
 743     if(!track->GetXYZAt(x,GetFieldkG(),xyz1)) return kFALSE;   // no prolongation
 744     xyz1[2]+=kEpsilon; // waiting for bug correction in geo
 745
 746     if(TMath::Abs(track->GetSnpAt(x,GetFieldkG())) >= maxSnp) return kFALSE;
 747     if(!track->PropagateTo(x,GetFieldkG()))  return kFALSE;
 748
 749     if(TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
 750     track->GetXYZ(xyz0);   // global position
 751     Double_t alphan = TMath::ATan2(xyz0[1], xyz0[0]);
 752     //
 753     Double_t ca=TMath::Cos(alphan-track->GetAlpha()),
 754       sa=TMath::Sin(alphan-track->GetAlpha());
 755     Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf));
 756     Double_t sinNew =  sf*ca - cf*sa;
 757     if(TMath::Abs(sinNew) >= maxSnp) return kFALSE;
 758     if(!track->Rotate(alphan)) return kFALSE;
 759
 760     xpos = track->GetX();
 761   }
 762   return kTRUE;
 763 }
 764 //---------------------------------------------------------------------------
 765 AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
 766                                                         const TObjArray *trkArray,
 767                                                         UShort_t *id,
 768                                                         const Float_t *diamondxy) const
 769 {
 770 //
 771 // Removes tracks in trksTree from fit of inVtx
 772 //
 773
 774   if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraint")) {
 775     printf("ERROR: primary vertex has no constraint: cannot remove tracks");
 776     return 0x0;
 777   }
 778   if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
 779     printf("WARNING: result of tracks' removal will be only approximately correct");
 780
 781   TMatrixD rv(3,1);
 782   rv(0,0) = inVtx->GetXv();
 783   rv(1,0) = inVtx->GetYv();
 784   rv(2,0) = inVtx->GetZv();
 785   TMatrixD vV(3,3);
 786   Double_t cov[6];
 787   inVtx->GetCovMatrix(cov);
 788   vV(0,0) = cov[0];
 789   vV(0,1) = cov[1]; vV(1,0) = cov[1];
 790   vV(1,1) = cov[2];
 791   vV(0,2) = cov[3]; vV(2,0) = cov[3];
 792   vV(1,2) = cov[4]; vV(2,1) = cov[4];
 793   vV(2,2) = cov[5];
 794
 795   TMatrixD sumWi(TMatrixD::kInverted,vV);
 796   TMatrixD sumWiri(sumWi,TMatrixD::kMult,rv);
 797
 798   Int_t nUsedTrks = inVtx->GetNIndices();
 799   Double_t chi2 = inVtx->GetChi2();
 800
 801   AliExternalTrackParam *track = 0;
 802   Int_t ntrks = (Int_t)trkArray->GetEntriesFast();
 803   for(Int_t i=0;i<ntrks;i++) {
 804     track = (AliExternalTrackParam*)trkArray->At(i);
 805     if(!inVtx->UsesTrack(id[i])) {
 806       printf("track %d was not used in vertex fit",id[i]);
 807       continue;
 808     }
 809     Double_t alpha = track->GetAlpha();
 810     Double_t xl = diamondxy[0]*TMath::Cos(alpha)+diamondxy[1]*TMath::Sin(alpha);
 811     track->PropagateTo(xl,GetFieldkG());
 812     // vector of track global coordinates
 813     TMatrixD ri(3,1);
 814     // covariance matrix of ri
 815     TMatrixD wWi(3,3);
 816
 817     // get space point from track
 818     if(!TrackToPoint(track,ri,wWi)) continue;
 819
 820     TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
 821
 822     sumWi -= wWi;
 823     sumWiri -= wWiri;
 824
 825     // track contribution to chi2
 826     TMatrixD deltar = rv; deltar -= ri;
 827     TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
 828     Double_t chi2i = deltar(0,0)*wWideltar(0,0)+
 829                      deltar(1,0)*wWideltar(1,0)+
 830                      deltar(2,0)*wWideltar(2,0);
 831     // remove from total chi2
 832     chi2 -= chi2i;
 833
 834     nUsedTrks--;
 835     if(nUsedTrks<2) {
 836       printf("Trying to remove too many tracks!");
 837       return 0x0;
 838     }
 839   }
 840
 841   TMatrixD rvnew(3,1);
 842   TMatrixD vVnew(3,3);
 843
 844   // new inverted of weights matrix
 845   TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
 846   vVnew = invsumWi;
 847   // new position of primary vertex
 848   rvnew.Mult(vVnew,sumWiri);
 849
 850   Double_t position[3];
 851   position[0] = rvnew(0,0);
 852   position[1] = rvnew(1,0);
 853   position[2] = rvnew(2,0);
 854   cov[0] = vVnew(0,0);
 855   cov[1] = vVnew(0,1);
 856   cov[2] = vVnew(1,1);
 857   cov[3] = vVnew(0,2);
 858   cov[4] = vVnew(1,2);
 859   cov[5] = vVnew(2,2);
 860
 861   // store data in the vertex object
 862   AliESDVertex *outVtx = new AliESDVertex(position,cov,chi2,nUsedTrks+1); // the +1 is for the constraint
 863   outVtx->SetTitle(inVtx->GetTitle());
 864   UShort_t *inindices = inVtx->GetIndices();
 865   Int_t nIndices = nUsedTrks;
 866   UShort_t *outindices = new UShort_t[nIndices];
 867   Int_t j=0;
 868   for(Int_t k=0; k<inVtx->GetNIndices(); k++) {
 869     Bool_t copyindex=kTRUE;
 870     for(Int_t l=0; l<ntrks; l++) {
 871       if(inindices[k]==id[l]) copyindex=kFALSE;
 872     }
 873     if(copyindex) {
 874       outindices[j] = inindices[k]; j++;
 875     }
 876   }
 877   outVtx->SetIndices(nIndices,outindices);
 878   if (outindices) delete [] outindices;
 879
 880   /*
 881     printf("Vertex before removing tracks:");
 882     inVtx->PrintStatus();
 883     inVtx->PrintIndices();
 884     printf("Vertex after removing tracks:");
 885     outVtx->PrintStatus();
 886     outVtx->PrintIndices();
 887   */
 888
 889   return outVtx;
 890 }
 891 //---------------------------------------------------------------------------
 892 AliESDVertex* AliVertexerTracks::RemoveConstraintFromVertex(AliESDVertex *inVtx,
 893                                                      Float_t *diamondxyz,
 894                                                      Float_t *diamondcov) const
 895 {
 896 //
 897 // Removes diamond constraint from fit of inVtx
 898 //
 899
 900   if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraint")) {
 901     printf("ERROR: primary vertex has no constraint: cannot remove it\n");
 902     return 0x0;
 903   }
 904   if(inVtx->GetNContributors()<3) {
 905     printf("ERROR: primary vertex has less than 2 tracks: cannot remove contraint\n");
 906     return 0x0;
 907   }
 908
 909   // diamond constraint
 910   TMatrixD vVb(3,3);
 911   vVb(0,0) = diamondcov[0];
 912   vVb(0,1) = diamondcov[1];
 913   vVb(0,2) = 0.;
 914   vVb(1,0) = diamondcov[1];
 915   vVb(1,1) = diamondcov[2];
 916   vVb(1,2) = 0.;
 917   vVb(2,0) = 0.;
 918   vVb(2,1) = 0.;
 919   vVb(2,2) = diamondcov[5];
 920   TMatrixD vVbInv(TMatrixD::kInverted,vVb);
 921   TMatrixD rb(3,1);
 922   rb(0,0) = diamondxyz[0];
 923   rb(1,0) = diamondxyz[1];
 924   rb(2,0) = diamondxyz[2];
 925   TMatrixD vVbInvrb(vVbInv,TMatrixD::kMult,rb);
 926
 927   // input vertex
 928   TMatrixD rv(3,1);
 929   rv(0,0) = inVtx->GetXv();
 930   rv(1,0) = inVtx->GetYv();
 931   rv(2,0) = inVtx->GetZv();
 932   TMatrixD vV(3,3);
 933   Double_t cov[6];
 934   inVtx->GetCovMatrix(cov);
 935   vV(0,0) = cov[0];
 936   vV(0,1) = cov[1]; vV(1,0) = cov[1];
 937   vV(1,1) = cov[2];
 938   vV(0,2) = cov[3]; vV(2,0) = cov[3];
 939   vV(1,2) = cov[4]; vV(2,1) = cov[4];
 940   vV(2,2) = cov[5];
 941   TMatrixD vVInv(TMatrixD::kInverted,vV);
 942   TMatrixD vVInvrv(vVInv,TMatrixD::kMult,rv);
 943
 944
 945   TMatrixD sumWi = vVInv - vVbInv;
 946
 947
 948   TMatrixD sumWiri = vVInvrv - vVbInvrb;
 949
 950   TMatrixD rvnew(3,1);
 951   TMatrixD vVnew(3,3);
 952
 953   // new inverted of weights matrix
 954   TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
 955   vVnew = invsumWi;
 956   // new position of primary vertex
 957   rvnew.Mult(vVnew,sumWiri);
 958
 959   Double_t position[3];
 960   position[0] = rvnew(0,0);
 961   position[1] = rvnew(1,0);
 962   position[2] = rvnew(2,0);
 963   cov[0] = vVnew(0,0);
 964   cov[1] = vVnew(0,1);
 965   cov[2] = vVnew(1,1);
 966   cov[3] = vVnew(0,2);
 967   cov[4] = vVnew(1,2);
 968   cov[5] = vVnew(2,2);
 969
 970
 971   Double_t chi2 = inVtx->GetChi2();
 972
 973   // diamond constribution to chi2
 974   TMatrixD deltar = rv; deltar -= rb;
 975   TMatrixD vVbInvdeltar(vVbInv,TMatrixD::kMult,deltar);
 976   Double_t chi2b = deltar(0,0)*vVbInvdeltar(0,0)+
 977                    deltar(1,0)*vVbInvdeltar(1,0)+
 978                    deltar(2,0)*vVbInvdeltar(2,0);
 979   // remove from total chi2
 980   chi2 -= chi2b;
 981
 982   // store data in the vertex object
 983   AliESDVertex *outVtx = new AliESDVertex(position,cov,chi2,inVtx->GetNContributors()-1);
 984   outVtx->SetTitle("VertexerTracksNoConstraint");
 985   UShort_t *inindices = inVtx->GetIndices();
 986   Int_t nIndices = inVtx->GetNIndices();
 987   outVtx->SetIndices(nIndices,inindices);
 988
 989   return outVtx;
 990 }
 991 //---------------------------------------------------------------------------
 992 void AliVertexerTracks::SetCuts(Double_t *cuts)
 993 {
 994 //
 995 //  Cut values
 996 //
 997   SetDCAcut(cuts[0]);
 998   SetDCAcutIter0(cuts[1]);
 999   SetMaxd0z0(cuts[2]);
1000   if(fMode==0 && cuts[3]<0) SetITSrefitNotRequired();
1001   SetMinClusters((Int_t)(TMath::Abs(cuts[3])));
1002   SetMinTracks((Int_t)(cuts[4]));
1003   SetNSigmad0(cuts[5]);
1004   SetMinDetFitter(cuts[6]);
1005   SetMaxTgl(cuts[7]);
1006   SetFiducialRZ(cuts[8],cuts[9]);
1007   fAlgo=(Int_t)(cuts[10]);
1008   fAlgoIter0=(Int_t)(cuts[11]);
1009   //
1010   if (cuts[12]>1.)   SetMVTukey2(cuts[12]);
1011   if (cuts[13]>1.)   SetMVSig2Ini(cuts[13]);
1012   if (cuts[14]>0.1)  SetMVMaxSigma2(cuts[14]);
1013   if (cuts[15]>1e-5) SetMVMinSig2Red(cuts[15]);
1014   if (cuts[16]>1e-5) SetMVMinDst(cuts[16]);
1015   if (cuts[17]>0.5)  SetMVScanStep(cuts[17]);
1016   SetMVMaxWghNtr(cuts[18]);
1017   SetMVFinalWBinary(cuts[19]>0);
1018   if (cuts[20]>20.)  SetBCSpacing(int(cuts[20]));
1019   //
1020   if (fAlgo==kMultiVertexer) SetSelectOnTOFBunchCrossing(kTRUE,kTRUE);
1021   else                       SetSelectOnTOFBunchCrossing(kFALSE,kTRUE);
1022   return;
1023 }
1024 //---------------------------------------------------------------------------
1025 void AliVertexerTracks::SetITSMode(Double_t dcacut,
1026                                    Double_t dcacutIter0,
1027                                    Double_t maxd0z0,
1028                                    Int_t minCls,
1029                                    Int_t mintrks,
1030                                    Double_t nsigma,
1031                                    Double_t mindetfitter,
1032                                    Double_t maxtgl,
1033                                    Double_t fidR,
1034                                    Double_t fidZ,
1035                                    Int_t finderAlgo,
1036                                    Int_t finderAlgoIter0)
1037 {
1038 //
1039 //  Cut values for ITS mode
1040 //
1041   fMode = 0;
1042   if(minCls>0) {
1043     SetITSrefitRequired();
1044   } else {
1045     SetITSrefitNotRequired();
1046   }
1047   SetDCAcut(dcacut);
1048   SetDCAcutIter0(dcacutIter0);
1049   SetMaxd0z0(maxd0z0);
1050   SetMinClusters(TMath::Abs(minCls));
1051   SetMinTracks(mintrks);
1052   SetNSigmad0(nsigma);
1053   SetMinDetFitter(mindetfitter);
1054   SetMaxTgl(maxtgl);
1055   SetFiducialRZ(fidR,fidZ);
1056   fAlgo=finderAlgo;
1057   fAlgoIter0=finderAlgoIter0;
1058
1059   return;
1060 }
1061 //---------------------------------------------------------------------------
1062 void AliVertexerTracks::SetTPCMode(Double_t dcacut,
1063                                    Double_t dcacutIter0,
1064                                    Double_t maxd0z0,
1065                                    Int_t minCls,
1066                                    Int_t mintrks,
1067                                    Double_t nsigma,
1068                                    Double_t mindetfitter,
1069                                    Double_t maxtgl,
1070                                    Double_t fidR,
1071                                    Double_t fidZ,
1072                                    Int_t finderAlgo,
1073                                    Int_t finderAlgoIter0)
1074 {
1075 //
1076 //  Cut values for TPC mode
1077 //
1078   fMode = 1;
1079   SetITSrefitNotRequired();
1080   SetDCAcut(dcacut);
1081   SetDCAcutIter0(dcacutIter0);
1082   SetMaxd0z0(maxd0z0);
1083   SetMinClusters(minCls);
1084   SetMinTracks(mintrks);
1085   SetNSigmad0(nsigma);
1086   SetMinDetFitter(mindetfitter);
1087   SetMaxTgl(maxtgl);
1088   SetFiducialRZ(fidR,fidZ);
1089   fAlgo=finderAlgo;
1090   fAlgoIter0=finderAlgoIter0;
1091
1092   return;
1093 }
1094 //---------------------------------------------------------------------------
1095 void AliVertexerTracks::SetSkipTracks(Int_t n,const Int_t *skipped)
1096 {
1097 //
1098 // Mark the tracks not to be used in the vertex reconstruction.
1099 // Tracks are identified by AliESDtrack::GetID()
1100 //
1101   fNTrksToSkip = n;  fTrksToSkip = new Int_t[n];
1102   for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i];
1103   return;
1104 }
1105 //---------------------------------------------------------------------------
1106 void  AliVertexerTracks::SetVtxStart(AliESDVertex *vtx)
1107 {
1108 //
1109 // Set initial vertex knowledge
1110 //
1111   vtx->GetXYZ(fNominalPos);
1112   vtx->GetCovMatrix(fNominalCov);
1113   SetConstraintOn();
1114   return;
1115 }
1116 //---------------------------------------------------------------------------
1117 void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
1118 {
1119   AliExternalTrackParam *track1;
1120   const Int_t knacc = (Int_t)fTrkArraySel.GetEntriesFast();
1121   AliStrLine **linarray = new AliStrLine* [knacc];
1122   for(Int_t i=0; i<knacc; i++){
1123     track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
1124     Double_t alpha=track1->GetAlpha();
1125     Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
1126     Double_t pos[3],dir[3],sigmasq[3];
1127     track1->GetXYZAt(mindist,GetFieldkG(),pos);
1128     track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
1129     sigmasq[0]=TMath::Sin(alpha)*TMath::Sin(alpha)*track1->GetSigmaY2();
1130     sigmasq[1]=TMath::Cos(alpha)*TMath::Cos(alpha)*track1->GetSigmaY2();
1131     sigmasq[2]=track1->GetSigmaZ2();
1132     TMatrixD ri(3,1);
1133     TMatrixD wWi(3,3);
1134     if(!TrackToPoint(track1,ri,wWi)) {optUseWeights=kFALSE;printf("WARNING\n");}
1135     Double_t wmat[9];
1136     Int_t iel=0;
1137     for(Int_t ia=0;ia<3;ia++){
1138       for(Int_t ib=0;ib<3;ib++){
1139         wmat[iel]=wWi(ia,ib);
1140         iel++;
1141       }
1142     }
1143     linarray[i] = new AliStrLine(pos,sigmasq,wmat,dir);
1144   }
1145   fVert=TrackletVertexFinder(linarray,knacc,optUseWeights);
1146   for(Int_t i=0; i<knacc; i++) delete linarray[i];
1147   delete [] linarray;
1148 }
1149 //---------------------------------------------------------------------------
1150 AliESDVertex AliVertexerTracks::TrackletVertexFinder(const TClonesArray *lines, Int_t optUseWeights)
1151 {
1152   // Calculate the point at minimum distance to prepared tracks (TClonesArray)
1153   const Int_t knacc = (Int_t)lines->GetEntriesFast();
1154   AliStrLine** lines2 = new AliStrLine* [knacc];
1155   for(Int_t i=0; i<knacc; i++){
1156     lines2[i]= (AliStrLine*)lines->At(i);
1157   }
1158   AliESDVertex vert = TrackletVertexFinder(lines2,knacc,optUseWeights);
1159   delete [] lines2;
1160   return vert;
1161 }
1162
1163 //---------------------------------------------------------------------------
1164 AliESDVertex AliVertexerTracks::TrackletVertexFinder(AliStrLine **lines, const Int_t knacc, Int_t optUseWeights)
1165 {
1166   // Calculate the point at minimum distance to prepared tracks (array of AliStrLine)
1167
1168   Double_t initPos[3]={0.,0.,0.};
1169
1170   Double_t (*vectP0)[3]=new Double_t [knacc][3];
1171   Double_t (*vectP1)[3]=new Double_t [knacc][3];
1172
1173   Double_t sum[3][3];
1174   Double_t dsum[3]={0,0,0};
1175   TMatrixD sumWi(3,3);
1176   for(Int_t i=0;i<3;i++){
1177     for(Int_t j=0;j<3;j++){
1178       sum[i][j]=0;
1179       sumWi(i,j)=0.;
1180     }
1181   }
1182
1183   for(Int_t i=0; i<knacc; i++){
1184     AliStrLine *line1 = lines[i];
1185     Double_t p0[3],cd[3],sigmasq[3];
1186     Double_t wmat[9];
1187     if(!line1) printf("ERROR %d %d\n",i,knacc);
1188     line1->GetP0(p0);
1189     line1->GetCd(cd);
1190     line1->GetSigma2P0(sigmasq);
1191     line1->GetWMatrix(wmat);
1192     TMatrixD wWi(3,3);
1193     Int_t iel=0;
1194     for(Int_t ia=0;ia<3;ia++){
1195       for(Int_t ib=0;ib<3;ib++){
1196         wWi(ia,ib)=wmat[iel];
1197         iel++;
1198       }
1199     }
1200
1201     sumWi+=wWi;
1202
1203     Double_t p1[3]={p0[0]+cd[0],p0[1]+cd[1],p0[2]+cd[2]};
1204     vectP0[i][0]=p0[0];
1205     vectP0[i][1]=p0[1];
1206     vectP0[i][2]=p0[2];
1207     vectP1[i][0]=p1[0];
1208     vectP1[i][1]=p1[1];
1209     vectP1[i][2]=p1[2];
1210
1211     Double_t matr[3][3];
1212     Double_t dknow[3];
1213     if(optUseWeights==0)GetStrLinDerivMatrix(p0,p1,matr,dknow);
1214     else GetStrLinDerivMatrix(p0,p1,sigmasq,matr,dknow);
1215
1216
1217     for(Int_t iii=0;iii<3;iii++){
1218       dsum[iii]+=dknow[iii];
1219       for(Int_t lj=0;lj<3;lj++) sum[iii][lj]+=matr[iii][lj];
1220     }
1221   }
1222
1223   TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
1224   Double_t covmatrix[6];
1225   covmatrix[0] = invsumWi(0,0);
1226   covmatrix[1] = invsumWi(0,1);
1227   covmatrix[2] = invsumWi(1,1);
1228   covmatrix[3] = invsumWi(0,2);
1229   covmatrix[4] = invsumWi(1,2);
1230   covmatrix[5] = invsumWi(2,2);
1231
1232   Double_t vett[3][3];
1233   Double_t det=GetDeterminant3X3(sum);
1234   Double_t sigma=0;
1235
1236   if(TMath::Abs(det) > kAlmost0){
1237     for(Int_t zz=0;zz<3;zz++){
1238       for(Int_t ww=0;ww<3;ww++){
1239         for(Int_t kk=0;kk<3;kk++) vett[ww][kk]=sum[ww][kk];
1240       }
1241       for(Int_t kk=0;kk<3;kk++) vett[kk][zz]=dsum[kk];
1242       initPos[zz]=GetDeterminant3X3(vett)/det;
1243     }
1244
1245
1246     for(Int_t i=0; i<knacc; i++){
1247       Double_t p0[3]={0,0,0},p1[3]={0,0,0};
1248       for(Int_t ii=0;ii<3;ii++){
1249         p0[ii]=vectP0[i][ii];
1250         p1[ii]=vectP1[i][ii];
1251       }
1252       sigma+=GetStrLinMinDist(p0,p1,initPos);
1253     }
1254
1255     if(sigma>0.) {sigma=TMath::Sqrt(sigma);}else{sigma=999;}
1256   }else{
1257     sigma=999;
1258   }
1259   AliESDVertex theVert(initPos,covmatrix,99999.,knacc);
1260   theVert.SetDispersion(sigma);
1261   delete [] vectP0;
1262   delete [] vectP1;
1263   return theVert;
1264 }
1265
1266 //---------------------------------------------------------------------------
1267 Bool_t AliVertexerTracks::TrackToPoint(AliExternalTrackParam *t,
1268                                        TMatrixD &ri,TMatrixD &wWi,
1269                                        Bool_t uUi3by3) const
1270 {
1271 //
1272 // Extract from the AliExternalTrackParam the global coordinates ri and covariance matrix
1273 // wWi of the space point that it represents (to be used in VertexFitter())
1274 //
1275
1276
1277   Double_t rotAngle = t->GetAlpha();
1278   if(rotAngle<0.) rotAngle += 2.*TMath::Pi();
1279   Double_t cosRot = TMath::Cos(rotAngle);
1280   Double_t sinRot = TMath::Sin(rotAngle);
1281   /*
1282   // RS >>>
1283   Double_t lambda = TMath::ATan(t->GetTgl());
1284   Double_t cosLam   = TMath::Cos(lambda);
1285   Double_t sinLam   = TMath::Sin(lambda);
1286   // RS <<<
1287   */
1288   ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot;
1289   ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot;
1290   ri(2,0) = t->GetZ();
1291
1292   if(!uUi3by3) {
1293     // matrix to go from global (x,y,z) to local (y,z);
1294     TMatrixD qQi(2,3);
1295     qQi(0,0) = -sinRot;
1296     qQi(0,1) = cosRot;
1297     qQi(0,2) = 0.;
1298     //
1299     qQi(1,0) = 0.;
1300     qQi(1,1) = 0.;
1301     qQi(1,2) = 1.;
1302     //
1303     // RS: Added polar inclination
1304     /*
1305     qQi(1,0) = -sinLam*cosRot;
1306     qQi(1,1) = -sinLam*sinRot;
1307     qQi(1,2) = cosLam;
1308     */
1309     // covariance matrix of local (y,z) - inverted
1310     TMatrixD uUi(2,2);
1311     uUi(0,0) = t->GetSigmaY2();
1312     uUi(0,1) = t->GetSigmaZY();
1313     uUi(1,0) = t->GetSigmaZY();
1314     uUi(1,1) = t->GetSigmaZ2();
1315     //printf(" Ui :");
1316     //printf(" %f   %f",uUi(0,0),uUi(0,1));
1317     //printf(" %f   %f",uUi(1,0),uUi(1,1));
1318
1319     if(uUi.Determinant() <= 0.) return kFALSE;
1320     TMatrixD uUiInv(TMatrixD::kInverted,uUi);
1321
1322     // weights matrix: wWi = qQiT * uUiInv * qQi
1323     TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
1324     TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
1325     wWi = m;
1326   } else {
1327     if(fVert.GetNContributors()<1) AliFatal("Vertex from finder is empty");
1328     // matrix to go from global (x,y,z) to local (x,y,z);
1329     TMatrixD qQi(3,3);
1330     qQi(0,0) = cosRot;
1331     qQi(0,1) = sinRot;
1332     qQi(0,2) = 0.;
1333     qQi(1,0) = -sinRot;
1334     qQi(1,1) = cosRot;
1335     qQi(1,2) = 0.;
1336     qQi(2,0) = 0.;
1337     qQi(2,1) = 0.;
1338     qQi(2,2) = 1.;
1339
1340     // covariance of fVert along the track
1341     Double_t p[3],pt,ptot;
1342     t->GetPxPyPz(p);
1343     pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
1344     ptot = TMath::Sqrt(pt*pt+p[2]*p[2]);
1345     Double_t cphi = p[0]/pt;               //cos(phi)=px/pt
1346     Double_t sphi = p[1]/pt;               //sin(phi)=py/pt
1347     Double_t clambda = pt/ptot;            //cos(lambda)=pt/ptot
1348     Double_t slambda = p[2]/ptot;            //sin(lambda)=pz/ptot
1349     Double_t covfVert[6];
1350     fVert.GetCovMatrix(covfVert);
1351     Double_t covfVertalongt =
1352        covfVert[0]*cphi*cphi*clambda*clambda
1353       +covfVert[1]*2.*cphi*sphi*clambda*clambda
1354       +covfVert[3]*2.*cphi*clambda*slambda
1355       +covfVert[2]*sphi*sphi*clambda*clambda
1356       +covfVert[4]*2.*sphi*clambda*slambda
1357       +covfVert[5]*slambda*slambda;
1358     // covariance matrix of local (x,y,z) - inverted
1359     TMatrixD uUi(3,3);
1360     uUi(0,0) = covfVertalongt * fnSigmaForUi00 * fnSigmaForUi00;
1361     AliDebug(1,Form("=====> sqrtUi00 cm  %f",TMath::Sqrt(uUi(0,0))));
1362     uUi(0,1) = 0.;
1363     uUi(0,2) = 0.;
1364     uUi(1,0) = 0.;
1365     uUi(1,1) = t->GetSigmaY2();
1366     uUi(1,2) = t->GetSigmaZY();
1367     uUi(2,0) = 0.;
1368     uUi(2,1) = t->GetSigmaZY();
1369     uUi(2,2) = t->GetSigmaZ2();
1370     //printf(" Ui :\n");
1371     //printf(" %f   %f\n",uUi(0,0),uUi(0,1));
1372     //printf(" %f   %f\n",uUi(1,0),uUi(1,1));
1373
1374     if(uUi.Determinant() <= 0.) return kFALSE;
1375     TMatrixD uUiInv(TMatrixD::kInverted,uUi);
1376
1377     // weights matrix: wWi = qQiT * uUiInv * qQi
1378     TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
1379     TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
1380     wWi = m;
1381   }
1382
1383
1384   return kTRUE;
1385 }
1386 //---------------------------------------------------------------------------
1387 void AliVertexerTracks::TooFewTracks()
1388 {
1389 //
1390 // When the number of tracks is < fMinTracks,
1391 // deal with vertices not found and prepare to exit
1392 //
1393   AliDebug(1,"TooFewTracks");
1394
1395   Double_t pos[3],err[3];
1396   pos[0] = fNominalPos[0];
1397   err[0] = TMath::Sqrt(fNominalCov[0]);
1398   pos[1] = fNominalPos[1];
1399   err[1] = TMath::Sqrt(fNominalCov[2]);
1400   pos[2] = fNominalPos[2];
1401   err[2] = TMath::Sqrt(fNominalCov[5]);
1402   Int_t    ncontr = (err[0]>1. ? -1 : -3);
1403   if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
1404   fCurrentVertex = new AliESDVertex(pos,err);
1405   fCurrentVertex->SetNContributors(ncontr);
1406
1407   if(fConstraint) {
1408     fCurrentVertex->SetTitle("VertexerTracksWithConstraint");
1409   } else {
1410     fCurrentVertex->SetTitle("VertexerTracksNoConstraint");
1411   }
1412
1413   if(!fTrkArraySel.IsEmpty()) fTrkArraySel.Delete();
1414   if(fIdSel) {delete [] fIdSel; fIdSel=NULL;}
1415   if(fTrksToSkip) {delete [] fTrksToSkip; fTrksToSkip=NULL;}
1416
1417   return;
1418 }
1419 //---------------------------------------------------------------------------
1420 void AliVertexerTracks::VertexFinder(Int_t optUseWeights)
1421 {
1422
1423   // Get estimate of vertex position in (x,y) from tracks DCA
1424
1425   Double_t initPos[3];
1426   initPos[2] = 0.;
1427   for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
1428   Int_t nacc = (Int_t)fTrkArraySel.GetEntriesFast();
1429   Double_t aver[3]={0.,0.,0.};
1430   Double_t aversq[3]={0.,0.,0.};
1431   Double_t sigmasq[3]={0.,0.,0.};
1432   Double_t sigma=0;
1433   Int_t ncombi = 0;
1434   AliExternalTrackParam *track1;
1435   AliExternalTrackParam *track2;
1436   Double_t pos[3],dir[3];
1437   Double_t alpha,mindist;
1438
1439   for(Int_t i=0; i<nacc; i++){
1440     track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
1441     alpha=track1->GetAlpha();
1442     mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
1443     track1->GetXYZAt(mindist,GetFieldkG(),pos);
1444     track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
1445     AliStrLine *line1 = new AliStrLine(pos,dir);
1446
1447    //    AliStrLine *line1 = new AliStrLine();
1448    //    track1->ApproximateHelixWithLine(mindist,GetFieldkG(),line1);
1449
1450     for(Int_t j=i+1; j<nacc; j++){
1451       track2 = (AliExternalTrackParam*)fTrkArraySel.At(j);
1452       alpha=track2->GetAlpha();
1453       mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
1454       track2->GetXYZAt(mindist,GetFieldkG(),pos);
1455       track2->GetPxPyPzAt(mindist,GetFieldkG(),dir);
1456       AliStrLine *line2 = new AliStrLine(pos,dir);
1457     //      AliStrLine *line2 = new AliStrLine();
1458     //  track2->ApproximateHelixWithLine(mindist,GetFieldkG(),line2);
1459       Double_t distCA=line2->GetDCA(line1);
1460       //printf("%d   %d   %f\n",i,j,distCA);
1461        if(fDCAcut<=0 || (fDCAcut>0&&distCA<fDCAcut)){
1462         Double_t pnt1[3],pnt2[3],crosspoint[3];
1463
1464         if(optUseWeights<=0){
1465           Int_t retcode = line2->Cross(line1,crosspoint);
1466           if(retcode>=0){
1467             ncombi++;
1468             for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
1469             for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
1470           }
1471         }
1472         if(optUseWeights>0){
1473           Int_t retcode = line1->CrossPoints(line2,pnt1,pnt2);
1474           if(retcode>=0){
1475             Double_t cs, sn;
1476             alpha=track1->GetAlpha();
1477             cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
1478             Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
1479             alpha=track2->GetAlpha();
1480             cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
1481             Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
1482             Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
1483             Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
1484             Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
1485             Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
1486             crosspoint[0]=wx1*pnt1[0] + wx2*pnt2[0];
1487             crosspoint[1]=wy1*pnt1[1] + wy2*pnt2[1];
1488             crosspoint[2]=wz1*pnt1[2] + wz2*pnt2[2];
1489
1490             ncombi++;
1491             for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
1492             for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
1493           }
1494         }
1495       }
1496       delete line2;
1497     }
1498     delete line1;
1499   }
1500   if(ncombi>0){
1501     for(Int_t jj=0;jj<3;jj++){
1502       initPos[jj] = aver[jj]/ncombi;
1503       //printf("%f\n",initPos[jj]);
1504       aversq[jj]/=ncombi;
1505       sigmasq[jj]=aversq[jj]-initPos[jj]*initPos[jj];
1506       sigma+=sigmasq[jj];
1507     }
1508     sigma=TMath::Sqrt(TMath::Abs(sigma));
1509   }
1510   else {
1511     Warning("VertexFinder","Finder did not succed");
1512     sigma=999;
1513   }
1514   fVert.SetXYZ(initPos);
1515   fVert.SetDispersion(sigma);
1516   fVert.SetNContributors(ncombi);
1517 }
1518 //---------------------------------------------------------------------------
1519 void AliVertexerTracks::VertexFitter(Bool_t vfit, Bool_t chiCalc,Int_t useWeights)
1520 {
1521 //
1522 // The optimal estimate of the vertex position is given by a "weighted
1523 // average of tracks positions".
1524 // Original method: V. Karimaki, CMS Note 97/0051
1525 //
1526   const double kTiny = 1e-9;
1527   Bool_t useConstraint = fConstraint;
1528   Double_t initPos[3];
1529   if(!fOnlyFitter) {
1530     fVert.GetXYZ(initPos);
1531   } else {
1532     initPos[0]=fNominalPos[0];
1533     initPos[1]=fNominalPos[1];
1534     initPos[2]=fNominalPos[2];
1535   }
1536
1537   Int_t nTrksSel = (Int_t)fTrkArraySel.GetEntries();
1538   if(nTrksSel==1) useConstraint=kTRUE;
1539   AliDebug(1,Form("--- VertexFitter(): start (%d,%d,%d)",vfit,chiCalc,useWeights));
1540   AliDebug(1,Form(" Number of tracks in array: %d\n",nTrksSel));
1541   AliDebug(1,Form(" Minimum # tracks required in fit: %d\n",fMinTracks));
1542   AliDebug(1,Form(" Vertex position after finder: %f,%f,%f\n",initPos[0],initPos[1],initPos[2]));
1543   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])));
1544
1545   // special treatment for few-tracks fits (e.g. secondary vertices)
1546   Bool_t uUi3by3 = kFALSE; if(nTrksSel<5 && !useConstraint && !useWeights) uUi3by3 = kTRUE;
1547
1548   Int_t i,j,k,step=0;
1549   static TMatrixD rv(3,1);
1550   static TMatrixD vV(3,3);
1551   rv(0,0) = initPos[0];
1552   rv(1,0) = initPos[1];
1553   rv(2,0) = initPos[2];
1554   Double_t xlStart,alpha;
1555   Int_t nTrksUsed = 0;
1556   Double_t chi2=0,chi2i,chi2b;
1557   AliExternalTrackParam *t = 0;
1558   Int_t failed = 0;
1559
1560   // initial vertex covariance matrix
1561   TMatrixD vVb(3,3);
1562   vVb(0,0) = fNominalCov[0];
1563   vVb(0,1) = fNominalCov[1];
1564   vVb(0,2) = 0.;
1565   vVb(1,0) = fNominalCov[1];
1566   vVb(1,1) = fNominalCov[2];
1567   vVb(1,2) = 0.;
1568   vVb(2,0) = 0.;
1569   vVb(2,1) = 0.;
1570   vVb(2,2) = fNominalCov[5];
1571   TMatrixD vVbInv(TMatrixD::kInverted,vVb);
1572   TMatrixD rb(3,1);
1573   rb(0,0) = fNominalPos[0];
1574   rb(1,0) = fNominalPos[1];
1575   rb(2,0) = fNominalPos[2];
1576   TMatrixD vVbInvrb(vVbInv,TMatrixD::kMult,rb);
1577   //
1578   int currBC = fVert.GetBC();
1579   //
1580   // 2 steps:
1581   // 1st - estimate of vtx using all tracks
1582   // 2nd - estimate of global chi2
1583   //
1584   for(step=0; step<2; step++) {
1585     if      (step==0 && !vfit)    continue;
1586     else if (step==1 && !chiCalc) continue;
1587     chi2 = 0.;
1588     nTrksUsed = 0;
1589     fMVWSum = fMVWE2 = 0;
1590     if(step==1) { initPos[0]=rv(0,0); initPos[1]=rv(1,0); initPos[2]=rv(2,0);}
1591     AliDebug(2,Form("Step%d: inipos: %+f %+f %+f MinTr: %d, Sig2:%.2f)",step,initPos[0],initPos[1],initPos[2],fMinTracks,fMVSigma2));
1592     //
1593     TMatrixD sumWiri(3,1);
1594     TMatrixD sumWi(3,3);
1595     for(i=0; i<3; i++) {
1596       sumWiri(i,0) = 0.;
1597       for(j=0; j<3; j++) sumWi(j,i) = 0.;
1598     }
1599
1600     // mean vertex constraint
1601     if(useConstraint) {
1602       for(i=0;i<3;i++) {
1603         sumWiri(i,0) += vVbInvrb(i,0);
1604         for(k=0;k<3;k++) sumWi(i,k) += vVbInv(i,k);
1605       }
1606       // chi2
1607       TMatrixD deltar = rv; deltar -= rb;
1608       TMatrixD vVbInvdeltar(vVbInv,TMatrixD::kMult,deltar);
1609       chi2b = deltar(0,0)*vVbInvdeltar(0,0)+
1610               deltar(1,0)*vVbInvdeltar(1,0)+
1611               deltar(2,0)*vVbInvdeltar(2,0);
1612       chi2 += chi2b;
1613     }
1614
1615     // loop on tracks
1616     for(k=0; k<nTrksSel; k++) {
1617       //
1618       // get track from track array
1619       t = (AliExternalTrackParam*)fTrkArraySel.At(k);
1620       if (useWeights && t->TestBit(kBitUsed)) continue;
1621       //
1622       int tBC = int(t->GetUniqueID()) - kTOFBCShift;    // BC assigned to this track
1623       if (fSelectOnTOFBunchCrossing) {
1624         if (!fKeepAlsoUnflaggedTOFBunchCrossing) continue;   // don't consider tracks with undefined BC
1625         if (currBC!=AliVTrack::kTOFBCNA && tBC!=AliVTrack::kTOFBCNA && tBC!=currBC) continue;  // track does not match to current BCid
1626       }
1627       alpha = t->GetAlpha();
1628       xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha);
1629       // to vtxSeed (from finder)
1630       t->PropagateTo(xlStart,GetFieldkG());
1631
1632       // vector of track global coordinates
1633       TMatrixD ri(3,1);
1634       // covariance matrix of ri
1635       TMatrixD wWi(3,3);
1636
1637       // get space point from track
1638       if(!TrackToPoint(t,ri,wWi,uUi3by3)) continue;
1639
1640       // track chi2
1641       TMatrixD deltar = rv; deltar -= ri;
1642       TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
1643       chi2i = deltar(0,0)*wWideltar(0,0)+
1644               deltar(1,0)*wWideltar(1,0)+
1645               deltar(2,0)*wWideltar(2,0);
1646       //
1647       if (useWeights) {
1648         //double sg = TMath::Sqrt(fMVSigma2);
1649         //double chi2iw = (deltar(0,0)*wWideltar(0,0)+deltar(1,0)*wWideltar(1,0))/sg + deltar(2,0)*wWideltar(2,0)/fMVSigma2;
1650         //double wgh = (1-chi2iw/fMVTukey2);
1651         double chi2iw = chi2i;
1652         double wgh = (1-chi2iw/fMVTukey2/fMVSigma2);
1653
1654         if (wgh<kTiny) wgh = 0;
1655         else if (useWeights==2) wgh = 1.; // use as binary weight
1656         if (step==1) ((AliESDtrack*)t)->SetBit(kBitUsed, wgh>0);
1657         if (wgh<kTiny) continue; // discard the track
1658         wWi *= wgh;  // RS: use weight?
1659         fMVWSum += wgh;
1660         fMVWE2  += wgh*chi2iw;
1661       }
1662       // add to total chi2
1663       if (fSelectOnTOFBunchCrossing && tBC!=AliVTrack::kTOFBCNA && currBC<0) currBC = tBC;
1664       //
1665       chi2 += chi2i;
1666       TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
1667       sumWiri += wWiri;
1668       sumWi   += wWi;
1669
1670       nTrksUsed++;
1671     } // end loop on tracks
1672
1673     if(nTrksUsed < fMinTracks) {
1674       failed=1;
1675       continue;
1676     }
1677
1678     Double_t determinant = sumWi.Determinant();
1679     if(determinant < fMinDetFitter)  {
1680       AliDebug(1,Form("det(V) = %f (<%f)\n",determinant,fMinDetFitter));
1681       failed=1;
1682       continue;
1683     }
1684
1685     if(step==0) {
1686       // inverted of weights matrix
1687       TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
1688       vV = invsumWi;
1689       // position of primary vertex
1690       rv.Mult(vV,sumWiri);
1691     }
1692   } // end loop on the 2 steps
1693
1694   if(failed) {
1695     fVert.SetNContributors(-1);
1696     if (chiCalc) {
1697       TooFewTracks();
1698       if (fCurrentVertex) fVert = *fCurrentVertex;  // RS
1699     }
1700     return;
1701   }
1702   //
1703   Double_t position[3];
1704   position[0] = rv(0,0);
1705   position[1] = rv(1,0);
1706   position[2] = rv(2,0);
1707   Double_t covmatrix[6];
1708   covmatrix[0] = vV(0,0);
1709   covmatrix[1] = vV(0,1);
1710   covmatrix[2] = vV(1,1);
1711   covmatrix[3] = vV(0,2);
1712   covmatrix[4] = vV(1,2);
1713   covmatrix[5] = vV(2,2);
1714
1715   // for correct chi2/ndf, count constraint as additional "track"
1716   if(fConstraint) nTrksUsed++;
1717   //
1718   if (vfit && !chiCalc) { // RS: special mode for multi-vertex finder
1719     fVert.SetXYZ(position);
1720     fVert.SetCovarianceMatrix(covmatrix);
1721     fVert.SetNContributors(nTrksUsed);
1722     return;
1723   }
1724   // store data in the vertex object
1725   if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
1726   fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);
1727   fCurrentVertex->SetBC(currBC);
1728   fVert = *fCurrentVertex;  // RS
1729   AliDebug(1," Vertex after fit:");
1730   AliDebug(1,Form("xyz: %f %f %f; nc %d",fCurrentVertex->GetXv(),fCurrentVertex->GetYv(),fCurrentVertex->GetZv(),fCurrentVertex->GetNContributors()));
1731   AliDebug(1,"--- VertexFitter(): finish\n");
1732
1733
1734   return;
1735 }
1736 //----------------------------------------------------------------------------
1737 AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(const TObjArray *trkArray,
1738                                                          UShort_t *id,
1739                                                          Bool_t optUseFitter,
1740                                                          Bool_t optPropagate,
1741                                                          Bool_t optUseDiamondConstraint)
1742 {
1743 //
1744 // Return vertex from tracks (AliExternalTrackParam) in array
1745 //
1746   fCurrentVertex = 0;
1747   // set optUseDiamondConstraint=TRUE only if you are reconstructing the
1748   // primary vertex!
1749   if(optUseDiamondConstraint) {
1750     SetConstraintOn();
1751   } else {
1752     SetConstraintOff();
1753   }
1754
1755   // get tracks and propagate them to initial vertex position
1756   fIdSel = new UShort_t[(Int_t)trkArray->GetEntriesFast()];
1757   Int_t nTrksSel = PrepareTracks(*trkArray,id,0);
1758   if((!optUseDiamondConstraint && nTrksSel<TMath::Max(2,fMinTracks)) ||
1759      (optUseDiamondConstraint && nTrksSel<1)) {
1760     TooFewTracks();
1761     return fCurrentVertex;
1762   }
1763
1764   // vertex finder
1765   if(nTrksSel==1) {
1766     AliDebug(1,"Just one track");
1767     OneTrackVertFinder();
1768   } else {
1769     switch (fAlgo) {
1770       case 1: StrLinVertexFinderMinDist(1); break;
1771       case 2: StrLinVertexFinderMinDist(0); break;
1772       case 3: HelixVertexFinder();          break;
1773       case 4: VertexFinder(1);              break;
1774       case 5: VertexFinder(0);              break;
1775       default: printf("Wrong algorithm\n"); break;
1776     }
1777   }
1778   AliDebug(1," Vertex finding completed\n");
1779
1780   // vertex fitter
1781   if(optUseFitter) {
1782     VertexFitter();
1783   } else {
1784     Double_t position[3]={fVert.GetXv(),fVert.GetYv(),fVert.GetZv()};
1785     Double_t covmatrix[6];
1786     fVert.GetCovMatrix(covmatrix);
1787     Double_t chi2=99999.;
1788     Int_t    nTrksUsed=fVert.GetNContributors();
1789     fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);
1790   }
1791   fCurrentVertex->SetDispersion(fVert.GetDispersion());
1792
1793
1794   // set indices of used tracks and propagate track to found vertex
1795   UShort_t *indices = 0;
1796   Double_t d0z0[2],covd0z0[3];
1797   AliExternalTrackParam *t = 0;
1798   if(fCurrentVertex->GetNContributors()>0) {
1799     indices = new UShort_t[fTrkArraySel.GetEntriesFast()];
1800     for(Int_t jj=0; jj<(Int_t)fTrkArraySel.GetEntriesFast(); jj++) {
1801       indices[jj] = fIdSel[jj];
1802       t = (AliExternalTrackParam*)fTrkArraySel.At(jj);
1803       if(optPropagate && optUseFitter) {
1804         if(TMath::Sqrt(fCurrentVertex->GetXv()*fCurrentVertex->GetXv()+fCurrentVertex->GetYv()*fCurrentVertex->GetYv())<3.) {
1805           t->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
1806           AliDebug(1,Form("Track %d propagated to found vertex",jj));
1807         } else {
1808           AliWarning("Found vertex outside beam pipe!");
1809         }
1810       }
1811     }
1812     fCurrentVertex->SetIndices(fCurrentVertex->GetNContributors(),indices);
1813   }
1814
1815   // clean up
1816   if (indices) {delete [] indices; indices=NULL;}
1817   delete [] fIdSel; fIdSel=NULL;
1818   fTrkArraySel.Delete();
1819
1820   return fCurrentVertex;
1821 }
1822
1823 //----------------------------------------------------------------------------
1824 AliESDVertex* AliVertexerTracks::VertexForSelectedESDTracks(TObjArray *trkArray,
1825                                                              Bool_t optUseFitter,
1826                                                             Bool_t optPropagate,
1827                                                             Bool_t optUseDiamondConstraint)
1828
1829 {
1830 //
1831 // Return vertex from array of ESD tracks
1832 //
1833
1834   Int_t nTrks = (Int_t)trkArray->GetEntriesFast();
1835   UShort_t *id = new UShort_t[nTrks];
1836
1837   AliESDtrack *esdt = 0;
1838   for(Int_t i=0; i<nTrks; i++){
1839     esdt = (AliESDtrack*)trkArray->At(i);
1840     id[i] = (UShort_t)esdt->GetID();
1841   }
1842
1843   VertexForSelectedTracks(trkArray,id,optUseFitter,optPropagate,optUseDiamondConstraint);
1844
1845   delete [] id; id=NULL;
1846
1847   return fCurrentVertex;
1848 }
1849
1850 //______________________________________________________
1851 Bool_t AliVertexerTracks::FindNextVertexMV()
1852 {
1853   // try to find a new vertex
1854   fMVSigma2 = fMVSig2Ini;
1855   double prevSig2 = fMVSigma2;
1856   double minDst2 = fMVMinDst*fMVMinDst;
1857   const double kSigLimit = 1.0;
1858   const double kSigLimitE = kSigLimit+1e-6;
1859   const double kPushFactor = 0.5;
1860   const int kMaxIter = 20;
1861   double push = kPushFactor;
1862   //
1863   int iter = 0;
1864   double posP[3]={0,0,0},pos[3]={0,0,0};
1865   fVert.GetXYZ(posP);
1866   //
1867
1868   do {
1869     fVert.SetBC(AliVTrack::kTOFBCNA);
1870     VertexFitter(kTRUE,kFALSE,1);
1871     if (fVert.GetNContributors()<fMinTracks) {
1872       AliDebug(3,Form("Failed in iteration %d: No Contributirs",iter));
1873       break;
1874     } // failed
1875     if (fMVWSum>0) fMVSigma2 = TMath::Max(fMVWE2/fMVWSum,kSigLimit);
1876     else {
1877       AliDebug(3,Form("Failed %d, no weithgs",iter));
1878       iter = kMaxIter+1;
1879       break;
1880     } // failed
1881     //
1882     double sigRed = (prevSig2-fMVSigma2)/prevSig2; // reduction of sigma2
1883     //
1884     fVert.GetXYZ(pos);
1885     double dst2 = (pos[0]-posP[0])*(pos[0]-posP[0])+(pos[1]-posP[1])*(pos[1]-posP[1])+(pos[2]-posP[2])*(pos[2]-posP[2]);
1886     AliDebug(3,Form("It#%2d Vtx: %+f %+f %+f Dst:%f Sig: %f [%.2f/%.2f] SigRed:%f",iter,pos[0],pos[1],pos[2],TMath::Sqrt(dst2),fMVSigma2,fMVWE2,fMVWSum,sigRed));
1887     if ( (++iter<kMaxIter) && (sigRed<0 || sigRed<fMVMinSig2Red) && fMVSigma2>fMVMaxSigma2) {
1888       fMVSigma2 *= push; // stuck, push little bit
1889       push *= kPushFactor;
1890       if (fMVSigma2<1.) fMVSigma2 = 1.;
1891       AliDebug(3,Form("Pushed sigma2 to %f",fMVSigma2));
1892     }
1893     else if (dst2<minDst2 && ((sigRed<0 || sigRed<fMVMinSig2Red) || fMVSigma2<kSigLimitE)) break;
1894     //
1895     fVert.GetXYZ(posP); // fetch previous vertex position
1896     prevSig2 = fMVSigma2;
1897   } while(iter<kMaxIter);
1898   //
1899   if (fVert.GetNContributors()<0 || iter>kMaxIter || fMVSigma2>fMVMaxSigma2) {
1900     return kFALSE;
1901   }
1902   else {
1903     VertexFitter(kFALSE,kTRUE,fMVFinalWBinary ? 2:1); // final chi2 calculation
1904     int nv = fMVVertices->GetEntries();
1905     // create indices
1906     int ntrk = fTrkArraySel.GetEntries();
1907     int nindices = fCurrentVertex->GetNContributors() - (fConstraint ? 1:0);
1908     UShort_t *indices = 0;
1909     if (nindices>0) indices = new UShort_t[nindices];
1910     int nadded = 0;
1911     for (int itr=0;itr<ntrk;itr++) {
1912       AliExternalTrackParam* t = (AliExternalTrackParam*)fTrkArraySel[itr];
1913       if (t->TestBit(kBitAccounted) || !t->TestBit(kBitUsed)) continue;   // already belongs to some vertex
1914       t->SetBit(kBitAccounted);
1915       indices[nadded++] = fIdSel[itr];
1916     }
1917     if (nadded!=nindices) printf("Mismatch : NInd: %d Nadd: %d\n",nindices,nadded);
1918     fCurrentVertex->SetIndices(nindices,indices);
1919     // set vertex title
1920     TString title="VertexerTracksMVNoConstraint";
1921     if(fConstraint) title="VertexerTracksMVWithConstraint";
1922     fCurrentVertex->SetTitle(title.Data());
1923     fMVVertices->AddLast(fCurrentVertex);
1924     AliDebug(3,Form("Added new vertex #%d NCont:%d XYZ: %f %f %f",nindices,nv,fCurrentVertex->GetX(),fCurrentVertex->GetY(),fCurrentVertex->GetZ()));
1925     if (indices) delete[] indices;
1926     fCurrentVertex = 0; // already attached to fMVVertices
1927     return kTRUE;
1928   }
1929 }
1930
1931 //______________________________________________________
1932 void AliVertexerTracks::FindVerticesMV()
1933 {
1934   // find and fit multiple vertices
1935   //
1936   double step = fMVScanStep>1 ?  fMVScanStep : 1.;
1937   double zmx = 3*TMath::Sqrt(fNominalCov[5]);
1938   double zmn = -zmx;
1939   int nz = (zmx-zmn)/step; if (nz<1) nz=1;
1940   double dz = (zmx-zmn)/nz;
1941   int izStart=0;
1942   AliDebug(2,Form("%d seeds between %f and %f",nz,zmn+dz/2,zmx+dz/2));
1943   //
1944   if (!fMVVertices) fMVVertices = new TObjArray(10);
1945   fMVVertices->Clear();
1946   //
1947   int ntrLeft = (Int_t)fTrkArraySel.GetEntries();
1948   //
1949   double sig2Scan = fMVSig2Ini;
1950   Bool_t runMore = kTRUE;
1951   int cntWide = 0;
1952   while (runMore) {
1953     fMVSig2Ini = sig2Scan*1e3;  // try wide search
1954     Bool_t found = kFALSE;
1955     cntWide++;
1956     fVert.SetNContributors(-1);
1957     fVert.SetXYZ(fNominalPos);
1958     AliDebug(3,Form("Wide search #%d Z= %f Sigma2=%f",cntWide,fNominalPos[2],fMVSig2Ini));
1959     if (FindNextVertexMV()) { // are there tracks left to consider?
1960       AliESDVertex* vtLast = (AliESDVertex*)fMVVertices->Last();
1961       if (vtLast && vtLast->GetNContributors()>0) ntrLeft -= vtLast->GetNContributors()-(fConstraint ? 1:0);
1962       if (ntrLeft<1) runMore = kFALSE;
1963       found = kTRUE;
1964       continue;
1965     }
1966     // if nothing is found, do narrow sig2ini scan
1967     fMVSig2Ini = sig2Scan;
1968     for (;izStart<nz;izStart++) {
1969       double zSeed = zmn+dz*(izStart+0.5);
1970       AliDebug(3,Form("Seed %d: Z= %f Sigma2=%f",izStart,zSeed,fMVSig2Ini));
1971       fVert.SetNContributors(-1);
1972       fVert.SetXYZ(fNominalPos);
1973       fVert.SetZv(zSeed);
1974       if (FindNextVertexMV()) { // are there tracks left to consider?
1975         AliESDVertex* vtLast = (AliESDVertex*)fMVVertices->Last();
1976         if (vtLast && vtLast->GetNContributors()>0) ntrLeft -= vtLast->GetNContributors()-(fConstraint ? 1:0);
1977         if (ntrLeft<1) runMore = kFALSE;
1978         found = kTRUE;
1979         break;
1980       }
1981     }
1982     runMore = found; // if nothing was found, no need for new iteration
1983   }
1984   fMVSig2Ini = sig2Scan;
1985   int nvFound = fMVVertices->GetEntriesFast();
1986   AliDebug(2,Form("Number of found vertices: %d",nvFound));
1987   if (nvFound<1) TooFewTracks();
1988   if (AliLog::GetGlobalDebugLevel()>0) fMVVertices->Print();
1989   //
1990 }
1991
1992 //______________________________________________________
1993 void AliVertexerTracks::AnalyzePileUp(AliESDEvent* esdEv)
1994 {
1995   // if multiple vertices are found, try to find the primary one and attach it as fCurrentVertex
1996   // then attach pile-up vertices directly to esdEv
1997   //
1998   int nFND = (fMVVertices && fMVVertices->GetEntriesFast()) ? fMVVertices->GetEntriesFast() : 0;
1999   if (nFND<1) { if (!fCurrentVertex) TooFewTracks(); return;} // no multiple vertices
2000   //
2001   int indCont[nFND];
2002   int nIndx[nFND];
2003   for (int iv=0;iv<nFND;iv++) {
2004     AliESDVertex* fnd = (AliESDVertex*)fMVVertices->At(iv);
2005     indCont[iv] = iv;
2006     nIndx[iv]   = fnd->GetNIndices();
2007   }
2008   TMath::Sort(nFND, nIndx, indCont, kTRUE); // sort in decreasing order of Nindices
2009   double dists[nFND];
2010   int    distOrd[nFND],indx[nFND];
2011   for (int iv=0;iv<nFND;iv++) {
2012     AliESDVertex* fndI = (AliESDVertex*)fMVVertices->At(indCont[iv]);
2013     if (fndI->GetStatus()<1) continue; // discarded
2014     int ncomp = 0;
2015     for (int jv=iv+1;jv<nFND;jv++) {
2016       AliESDVertex* fndJ = (AliESDVertex*)fMVVertices->At(indCont[jv]);
2017       if (fndJ->GetStatus()<1) continue;
2018       dists[ncomp] = fndI->GetWDist(fndJ)*fndJ->GetNIndices();
2019       distOrd[ncomp] = indCont[jv];
2020       indx[ncomp]  = ncomp;
2021       ncomp++;
2022     }
2023     if (ncomp<1) continue;
2024     TMath::Sort(ncomp, dists, indx, kFALSE); // sort in increasing distance order
2025     for (int jv0=0;jv0<ncomp;jv0++) {
2026       int jv = distOrd[indx[jv0]];
2027       AliESDVertex* fndJ = (AliESDVertex*)fMVVertices->At(jv);
2028       if (dists[indx[jv0]]<fMVMaxWghNtr) { // candidate for split vertex
2029         //before eliminating the small close vertex, check if we should transfere its BC to largest one
2030         if (fndJ->GetBC()!=AliVTrack::kTOFBCNA && fndI->GetBC()==AliVTrack::kTOFBCNA) fndI->SetBC(fndJ->GetBC());
2031         //
2032         // leave the vertex only if both vertices have definit but different BC's, then this is not splitting.
2033         if ( (fndJ->GetBC()==fndI->GetBC()) || (fndJ->GetBC()==AliVTrack::kTOFBCNA)) fndJ->SetNContributors(-fndJ->GetNContributors());
2034       }
2035     }
2036   }
2037   //
2038   // select as a primary the largest vertex with BC=0, or the largest with BC non-ID
2039   int primBC0=-1,primNoBC=-1;
2040   for (int iv0=0;iv0<nFND;iv0++) {
2041     int iv = indCont[iv0];
2042     AliESDVertex* fndI = (AliESDVertex*)fMVVertices->At(iv);
2043     if (!fndI) continue;
2044     if (fndI->GetStatus()<1) {fMVVertices->RemoveAt(iv); delete fndI; continue;}
2045     if (primBC0<0  && fndI->GetBC()==0) primBC0 = iv;
2046     if (primNoBC<0 && fndI->GetBC()==AliVTrack::kTOFBCNA)  primNoBC = iv;
2047   }
2048   //
2049   if (primBC0>=0) fCurrentVertex = (AliESDVertex*)fMVVertices->At(primBC0);
2050   if (!fCurrentVertex && primNoBC>=0) fCurrentVertex = (AliESDVertex*)fMVVertices->At(primNoBC);
2051   if (fCurrentVertex) fMVVertices->Remove(fCurrentVertex);
2052   else {  // all vertices have BC>0, no primary vertex
2053     fCurrentVertex = new AliESDVertex();
2054     fCurrentVertex->SetNContributors(-3);
2055     fCurrentVertex->SetBC(AliVTrack::kTOFBCNA);
2056   }
2057   fCurrentVertex->SetID(-1);
2058   //
2059   // add pileup vertices
2060   int nadd = 0;
2061   for (int iv0=0;iv0<nFND;iv0++) {
2062     int iv = indCont[iv0];
2063     AliESDVertex* fndI = (AliESDVertex*)fMVVertices->At(iv);
2064     if (!fndI) continue;
2065     fndI->SetID(++nadd);
2066     esdEv->AddPileupVertexTracks(fndI);
2067   }
2068   //
2069   fMVVertices->Delete();
2070   //
2071 }
2072