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, Int_t ncuts)
 993 {
 994 //
 995 //  Cut values
 996 //
 997   if (ncuts>0) SetDCAcut(cuts[0]);
 998   if (ncuts>1) SetDCAcutIter0(cuts[1]);
 999   if (ncuts>2) SetMaxd0z0(cuts[2]);
1000   if (ncuts>3) if(fMode==0 && cuts[3]<0) SetITSrefitNotRequired();
1001   if (ncuts>4) SetMinClusters((Int_t)(TMath::Abs(cuts[3])));
1002   if (ncuts>5) SetMinTracks((Int_t)(cuts[4]));
1003   if (ncuts>6) SetNSigmad0(cuts[5]);
1004   if (ncuts>7) SetMinDetFitter(cuts[6]);
1005   if (ncuts>8) SetMaxTgl(cuts[7]);
1006   if (ncuts>9) SetFiducialRZ(cuts[8],cuts[9]);
1007   if (ncuts>10) fAlgo=(Int_t)(cuts[10]);
1008   if (ncuts>11) fAlgoIter0=(Int_t)(cuts[11]);
1009   //
1010   if (ncuts>12) if (cuts[12]>1.)   SetMVTukey2(cuts[12]);
1011   if (ncuts>13) if (cuts[13]>1.)   SetMVSig2Ini(cuts[13]);
1012   if (ncuts>14) if (cuts[14]>0.1)  SetMVMaxSigma2(cuts[14]);
1013   if (ncuts>15) if (cuts[15]>1e-5) SetMVMinSig2Red(cuts[15]);
1014   if (ncuts>16) if (cuts[16]>1e-5) SetMVMinDst(cuts[16]);
1015   if (ncuts>17) if (cuts[17]>0.5)  SetMVScanStep(cuts[17]);
1016   if (ncuts>18) SetMVMaxWghNtr(cuts[18]);
1017   if (ncuts>19) SetMVFinalWBinary(cuts[19]>0);
1018   if (ncuts>20) if (cuts[20]>20.)  SetBCSpacing(int(cuts[20]));
1019   //
1020   if (ncuts>21) {
1021     if (fAlgo==kMultiVertexer) SetSelectOnTOFBunchCrossing(kTRUE,kTRUE);
1022     else                       SetSelectOnTOFBunchCrossing(kFALSE,kTRUE);
1023   }
1024   return;
1025 }
1026 //---------------------------------------------------------------------------
1027 void AliVertexerTracks::SetITSMode(Double_t dcacut,
1028                                    Double_t dcacutIter0,
1029                                    Double_t maxd0z0,
1030                                    Int_t minCls,
1031                                    Int_t mintrks,
1032                                    Double_t nsigma,
1033                                    Double_t mindetfitter,
1034                                    Double_t maxtgl,
1035                                    Double_t fidR,
1036                                    Double_t fidZ,
1037                                    Int_t finderAlgo,
1038                                    Int_t finderAlgoIter0)
1039 {
1040 //
1041 //  Cut values for ITS mode
1042 //
1043   fMode = 0;
1044   if(minCls>0) {
1045     SetITSrefitRequired();
1046   } else {
1047     SetITSrefitNotRequired();
1048   }
1049   SetDCAcut(dcacut);
1050   SetDCAcutIter0(dcacutIter0);
1051   SetMaxd0z0(maxd0z0);
1052   SetMinClusters(TMath::Abs(minCls));
1053   SetMinTracks(mintrks);
1054   SetNSigmad0(nsigma);
1055   SetMinDetFitter(mindetfitter);
1056   SetMaxTgl(maxtgl);
1057   SetFiducialRZ(fidR,fidZ);
1058   fAlgo=finderAlgo;
1059   fAlgoIter0=finderAlgoIter0;
1060
1061   return;
1062 }
1063 //---------------------------------------------------------------------------
1064 void AliVertexerTracks::SetTPCMode(Double_t dcacut,
1065                                    Double_t dcacutIter0,
1066                                    Double_t maxd0z0,
1067                                    Int_t minCls,
1068                                    Int_t mintrks,
1069                                    Double_t nsigma,
1070                                    Double_t mindetfitter,
1071                                    Double_t maxtgl,
1072                                    Double_t fidR,
1073                                    Double_t fidZ,
1074                                    Int_t finderAlgo,
1075                                    Int_t finderAlgoIter0)
1076 {
1077 //
1078 //  Cut values for TPC mode
1079 //
1080   fMode = 1;
1081   SetITSrefitNotRequired();
1082   SetDCAcut(dcacut);
1083   SetDCAcutIter0(dcacutIter0);
1084   SetMaxd0z0(maxd0z0);
1085   SetMinClusters(minCls);
1086   SetMinTracks(mintrks);
1087   SetNSigmad0(nsigma);
1088   SetMinDetFitter(mindetfitter);
1089   SetMaxTgl(maxtgl);
1090   SetFiducialRZ(fidR,fidZ);
1091   fAlgo=finderAlgo;
1092   fAlgoIter0=finderAlgoIter0;
1093
1094   return;
1095 }
1096 //---------------------------------------------------------------------------
1097 void AliVertexerTracks::SetSkipTracks(Int_t n,const Int_t *skipped)
1098 {
1099 //
1100 // Mark the tracks not to be used in the vertex reconstruction.
1101 // Tracks are identified by AliESDtrack::GetID()
1102 //
1103   fNTrksToSkip = n;  fTrksToSkip = new Int_t[n];
1104   for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i];
1105   return;
1106 }
1107 //---------------------------------------------------------------------------
1108 void  AliVertexerTracks::SetVtxStart(AliESDVertex *vtx)
1109 {
1110 //
1111 // Set initial vertex knowledge
1112 //
1113   vtx->GetXYZ(fNominalPos);
1114   vtx->GetCovMatrix(fNominalCov);
1115   SetConstraintOn();
1116   return;
1117 }
1118 //---------------------------------------------------------------------------
1119 void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
1120 {
1121   AliExternalTrackParam *track1;
1122   const Int_t knacc = (Int_t)fTrkArraySel.GetEntriesFast();
1123   AliStrLine **linarray = new AliStrLine* [knacc];
1124   for(Int_t i=0; i<knacc; i++){
1125     track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
1126     Double_t alpha=track1->GetAlpha();
1127     Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
1128     Double_t pos[3],dir[3],sigmasq[3];
1129     track1->GetXYZAt(mindist,GetFieldkG(),pos);
1130     track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
1131     sigmasq[0]=TMath::Sin(alpha)*TMath::Sin(alpha)*track1->GetSigmaY2();
1132     sigmasq[1]=TMath::Cos(alpha)*TMath::Cos(alpha)*track1->GetSigmaY2();
1133     sigmasq[2]=track1->GetSigmaZ2();
1134     TMatrixD ri(3,1);
1135     TMatrixD wWi(3,3);
1136     if(!TrackToPoint(track1,ri,wWi)) {optUseWeights=kFALSE;printf("WARNING\n");}
1137     Double_t wmat[9];
1138     Int_t iel=0;
1139     for(Int_t ia=0;ia<3;ia++){
1140       for(Int_t ib=0;ib<3;ib++){
1141         wmat[iel]=wWi(ia,ib);
1142         iel++;
1143       }
1144     }
1145     linarray[i] = new AliStrLine(pos,sigmasq,wmat,dir);
1146   }
1147   fVert=TrackletVertexFinder(linarray,knacc,optUseWeights);
1148   for(Int_t i=0; i<knacc; i++) delete linarray[i];
1149   delete [] linarray;
1150 }
1151 //---------------------------------------------------------------------------
1152 AliESDVertex AliVertexerTracks::TrackletVertexFinder(const TClonesArray *lines, Int_t optUseWeights)
1153 {
1154   // Calculate the point at minimum distance to prepared tracks (TClonesArray)
1155   const Int_t knacc = (Int_t)lines->GetEntriesFast();
1156   AliStrLine** lines2 = new AliStrLine* [knacc];
1157   for(Int_t i=0; i<knacc; i++){
1158     lines2[i]= (AliStrLine*)lines->At(i);
1159   }
1160   AliESDVertex vert = TrackletVertexFinder(lines2,knacc,optUseWeights);
1161   delete [] lines2;
1162   return vert;
1163 }
1164
1165 //---------------------------------------------------------------------------
1166 AliESDVertex AliVertexerTracks::TrackletVertexFinder(AliStrLine **lines, const Int_t knacc, Int_t optUseWeights)
1167 {
1168   // Calculate the point at minimum distance to prepared tracks (array of AliStrLine)
1169
1170   Double_t initPos[3]={0.,0.,0.};
1171
1172   Double_t (*vectP0)[3]=new Double_t [knacc][3];
1173   Double_t (*vectP1)[3]=new Double_t [knacc][3];
1174
1175   Double_t sum[3][3];
1176   Double_t dsum[3]={0,0,0};
1177   TMatrixD sumWi(3,3);
1178   for(Int_t i=0;i<3;i++){
1179     for(Int_t j=0;j<3;j++){
1180       sum[i][j]=0;
1181       sumWi(i,j)=0.;
1182     }
1183   }
1184
1185   for(Int_t i=0; i<knacc; i++){
1186     AliStrLine *line1 = lines[i];
1187     Double_t p0[3],cd[3],sigmasq[3];
1188     Double_t wmat[9];
1189     if(!line1) {printf("ERROR %d %d\n",i,knacc); continue;}
1190     line1->GetP0(p0);
1191     line1->GetCd(cd);
1192     line1->GetSigma2P0(sigmasq);
1193     line1->GetWMatrix(wmat);
1194     TMatrixD wWi(3,3);
1195     Int_t iel=0;
1196     for(Int_t ia=0;ia<3;ia++){
1197       for(Int_t ib=0;ib<3;ib++){
1198         wWi(ia,ib)=wmat[iel];
1199         iel++;
1200       }
1201     }
1202
1203     sumWi+=wWi;
1204
1205     Double_t p1[3]={p0[0]+cd[0],p0[1]+cd[1],p0[2]+cd[2]};
1206     vectP0[i][0]=p0[0];
1207     vectP0[i][1]=p0[1];
1208     vectP0[i][2]=p0[2];
1209     vectP1[i][0]=p1[0];
1210     vectP1[i][1]=p1[1];
1211     vectP1[i][2]=p1[2];
1212
1213     Double_t matr[3][3];
1214     Double_t dknow[3];
1215     if(optUseWeights==0)GetStrLinDerivMatrix(p0,p1,matr,dknow);
1216     else GetStrLinDerivMatrix(p0,p1,sigmasq,matr,dknow);
1217
1218
1219     for(Int_t iii=0;iii<3;iii++){
1220       dsum[iii]+=dknow[iii];
1221       for(Int_t lj=0;lj<3;lj++) sum[iii][lj]+=matr[iii][lj];
1222     }
1223   }
1224
1225   TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
1226   Double_t covmatrix[6];
1227   covmatrix[0] = invsumWi(0,0);
1228   covmatrix[1] = invsumWi(0,1);
1229   covmatrix[2] = invsumWi(1,1);
1230   covmatrix[3] = invsumWi(0,2);
1231   covmatrix[4] = invsumWi(1,2);
1232   covmatrix[5] = invsumWi(2,2);
1233
1234   Double_t vett[3][3];
1235   Double_t det=GetDeterminant3X3(sum);
1236   Double_t sigma=0;
1237
1238   if(TMath::Abs(det) > kAlmost0){
1239     for(Int_t zz=0;zz<3;zz++){
1240       for(Int_t ww=0;ww<3;ww++){
1241         for(Int_t kk=0;kk<3;kk++) vett[ww][kk]=sum[ww][kk];
1242       }
1243       for(Int_t kk=0;kk<3;kk++) vett[kk][zz]=dsum[kk];
1244       initPos[zz]=GetDeterminant3X3(vett)/det;
1245     }
1246
1247
1248     for(Int_t i=0; i<knacc; i++){
1249       Double_t p0[3]={0,0,0},p1[3]={0,0,0};
1250       for(Int_t ii=0;ii<3;ii++){
1251         p0[ii]=vectP0[i][ii];
1252         p1[ii]=vectP1[i][ii];
1253       }
1254       sigma+=GetStrLinMinDist(p0,p1,initPos);
1255     }
1256
1257     if(sigma>0.) {sigma=TMath::Sqrt(sigma);}else{sigma=999;}
1258   }else{
1259     sigma=999;
1260   }
1261   AliESDVertex theVert(initPos,covmatrix,99999.,knacc);
1262   theVert.SetDispersion(sigma);
1263   delete [] vectP0;
1264   delete [] vectP1;
1265   return theVert;
1266 }
1267
1268 //---------------------------------------------------------------------------
1269 Bool_t AliVertexerTracks::TrackToPoint(AliExternalTrackParam *t,
1270                                        TMatrixD &ri,TMatrixD &wWi,
1271                                        Bool_t uUi3by3) const
1272 {
1273 //
1274 // Extract from the AliExternalTrackParam the global coordinates ri and covariance matrix
1275 // wWi of the space point that it represents (to be used in VertexFitter())
1276 //
1277
1278
1279   Double_t rotAngle = t->GetAlpha();
1280   if(rotAngle<0.) rotAngle += 2.*TMath::Pi();
1281   Double_t cosRot = TMath::Cos(rotAngle);
1282   Double_t sinRot = TMath::Sin(rotAngle);
1283   /*
1284   // RS >>>
1285   Double_t lambda = TMath::ATan(t->GetTgl());
1286   Double_t cosLam   = TMath::Cos(lambda);
1287   Double_t sinLam   = TMath::Sin(lambda);
1288   // RS <<<
1289   */
1290   ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot;
1291   ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot;
1292   ri(2,0) = t->GetZ();
1293
1294   if(!uUi3by3) {
1295     // matrix to go from global (x,y,z) to local (y,z);
1296     TMatrixD qQi(2,3);
1297     qQi(0,0) = -sinRot;
1298     qQi(0,1) = cosRot;
1299     qQi(0,2) = 0.;
1300     //
1301     qQi(1,0) = 0.;
1302     qQi(1,1) = 0.;
1303     qQi(1,2) = 1.;
1304     //
1305     // RS: Added polar inclination
1306     /*
1307     qQi(1,0) = -sinLam*cosRot;
1308     qQi(1,1) = -sinLam*sinRot;
1309     qQi(1,2) = cosLam;
1310     */
1311     // covariance matrix of local (y,z) - inverted
1312     TMatrixD uUi(2,2);
1313     uUi(0,0) = t->GetSigmaY2();
1314     uUi(0,1) = t->GetSigmaZY();
1315     uUi(1,0) = t->GetSigmaZY();
1316     uUi(1,1) = t->GetSigmaZ2();
1317     //printf(" Ui :");
1318     //printf(" %f   %f",uUi(0,0),uUi(0,1));
1319     //printf(" %f   %f",uUi(1,0),uUi(1,1));
1320
1321     if(uUi.Determinant() <= 0.) return kFALSE;
1322     TMatrixD uUiInv(TMatrixD::kInverted,uUi);
1323
1324     // weights matrix: wWi = qQiT * uUiInv * qQi
1325     TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
1326     TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
1327     wWi = m;
1328   } else {
1329     if(fVert.GetNContributors()<1) AliFatal("Vertex from finder is empty");
1330     // matrix to go from global (x,y,z) to local (x,y,z);
1331     TMatrixD qQi(3,3);
1332     qQi(0,0) = cosRot;
1333     qQi(0,1) = sinRot;
1334     qQi(0,2) = 0.;
1335     qQi(1,0) = -sinRot;
1336     qQi(1,1) = cosRot;
1337     qQi(1,2) = 0.;
1338     qQi(2,0) = 0.;
1339     qQi(2,1) = 0.;
1340     qQi(2,2) = 1.;
1341
1342     // covariance of fVert along the track
1343     Double_t p[3],pt,ptot;
1344     t->GetPxPyPz(p);
1345     pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
1346     ptot = TMath::Sqrt(pt*pt+p[2]*p[2]);
1347     Double_t cphi = p[0]/pt;               //cos(phi)=px/pt
1348     Double_t sphi = p[1]/pt;               //sin(phi)=py/pt
1349     Double_t clambda = pt/ptot;            //cos(lambda)=pt/ptot
1350     Double_t slambda = p[2]/ptot;            //sin(lambda)=pz/ptot
1351     Double_t covfVert[6];
1352     fVert.GetCovMatrix(covfVert);
1353     Double_t covfVertalongt =
1354        covfVert[0]*cphi*cphi*clambda*clambda
1355       +covfVert[1]*2.*cphi*sphi*clambda*clambda
1356       +covfVert[3]*2.*cphi*clambda*slambda
1357       +covfVert[2]*sphi*sphi*clambda*clambda
1358       +covfVert[4]*2.*sphi*clambda*slambda
1359       +covfVert[5]*slambda*slambda;
1360     // covariance matrix of local (x,y,z) - inverted
1361     TMatrixD uUi(3,3);
1362     uUi(0,0) = covfVertalongt * fnSigmaForUi00 * fnSigmaForUi00;
1363     AliDebug(1,Form("=====> sqrtUi00 cm  %f",TMath::Sqrt(uUi(0,0))));
1364     uUi(0,1) = 0.;
1365     uUi(0,2) = 0.;
1366     uUi(1,0) = 0.;
1367     uUi(1,1) = t->GetSigmaY2();
1368     uUi(1,2) = t->GetSigmaZY();
1369     uUi(2,0) = 0.;
1370     uUi(2,1) = t->GetSigmaZY();
1371     uUi(2,2) = t->GetSigmaZ2();
1372     //printf(" Ui :\n");
1373     //printf(" %f   %f\n",uUi(0,0),uUi(0,1));
1374     //printf(" %f   %f\n",uUi(1,0),uUi(1,1));
1375
1376     if(uUi.Determinant() <= 0.) return kFALSE;
1377     TMatrixD uUiInv(TMatrixD::kInverted,uUi);
1378
1379     // weights matrix: wWi = qQiT * uUiInv * qQi
1380     TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
1381     TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
1382     wWi = m;
1383   }
1384
1385
1386   return kTRUE;
1387 }
1388 //---------------------------------------------------------------------------
1389 void AliVertexerTracks::TooFewTracks()
1390 {
1391 //
1392 // When the number of tracks is < fMinTracks,
1393 // deal with vertices not found and prepare to exit
1394 //
1395   AliDebug(1,"TooFewTracks");
1396
1397   Double_t pos[3],err[3];
1398   pos[0] = fNominalPos[0];
1399   err[0] = TMath::Sqrt(fNominalCov[0]);
1400   pos[1] = fNominalPos[1];
1401   err[1] = TMath::Sqrt(fNominalCov[2]);
1402   pos[2] = fNominalPos[2];
1403   err[2] = TMath::Sqrt(fNominalCov[5]);
1404   Int_t    ncontr = (err[0]>1. ? -1 : -3);
1405   if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
1406   fCurrentVertex = new AliESDVertex(pos,err);
1407   fCurrentVertex->SetNContributors(ncontr);
1408
1409   if(fConstraint) {
1410     fCurrentVertex->SetTitle("VertexerTracksWithConstraint");
1411   } else {
1412     fCurrentVertex->SetTitle("VertexerTracksNoConstraint");
1413   }
1414
1415   if(!fTrkArraySel.IsEmpty()) fTrkArraySel.Delete();
1416   if(fIdSel) {delete [] fIdSel; fIdSel=NULL;}
1417   if(fTrksToSkip) {delete [] fTrksToSkip; fTrksToSkip=NULL;}
1418
1419   return;
1420 }
1421 //---------------------------------------------------------------------------
1422 void AliVertexerTracks::VertexFinder(Int_t optUseWeights)
1423 {
1424
1425   // Get estimate of vertex position in (x,y) from tracks DCA
1426
1427   Double_t initPos[3];
1428   initPos[2] = 0.;
1429   for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
1430   Int_t nacc = (Int_t)fTrkArraySel.GetEntriesFast();
1431   Double_t aver[3]={0.,0.,0.};
1432   Double_t aversq[3]={0.,0.,0.};
1433   Double_t sigmasq[3]={0.,0.,0.};
1434   Double_t sigma=0;
1435   Int_t ncombi = 0;
1436   AliExternalTrackParam *track1;
1437   AliExternalTrackParam *track2;
1438   Double_t pos[3],dir[3];
1439   Double_t alpha,mindist;
1440
1441   for(Int_t i=0; i<nacc; i++){
1442     track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
1443     alpha=track1->GetAlpha();
1444     mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
1445     track1->GetXYZAt(mindist,GetFieldkG(),pos);
1446     track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
1447     AliStrLine *line1 = new AliStrLine(pos,dir);
1448
1449    //    AliStrLine *line1 = new AliStrLine();
1450    //    track1->ApproximateHelixWithLine(mindist,GetFieldkG(),line1);
1451
1452     for(Int_t j=i+1; j<nacc; j++){
1453       track2 = (AliExternalTrackParam*)fTrkArraySel.At(j);
1454       alpha=track2->GetAlpha();
1455       mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
1456       track2->GetXYZAt(mindist,GetFieldkG(),pos);
1457       track2->GetPxPyPzAt(mindist,GetFieldkG(),dir);
1458       AliStrLine *line2 = new AliStrLine(pos,dir);
1459     //      AliStrLine *line2 = new AliStrLine();
1460     //  track2->ApproximateHelixWithLine(mindist,GetFieldkG(),line2);
1461       Double_t distCA=line2->GetDCA(line1);
1462       //printf("%d   %d   %f\n",i,j,distCA);
1463        if(fDCAcut<=0 || (fDCAcut>0&&distCA<fDCAcut)){
1464         Double_t pnt1[3],pnt2[3],crosspoint[3];
1465
1466         if(optUseWeights<=0){
1467           Int_t retcode = line2->Cross(line1,crosspoint);
1468           if(retcode>=0){
1469             ncombi++;
1470             for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
1471             for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
1472           }
1473         }
1474         if(optUseWeights>0){
1475           Int_t retcode = line1->CrossPoints(line2,pnt1,pnt2);
1476           if(retcode>=0){
1477             Double_t cs, sn;
1478             alpha=track1->GetAlpha();
1479             cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
1480             Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
1481             alpha=track2->GetAlpha();
1482             cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
1483             Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
1484             Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
1485             Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
1486             Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
1487             Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
1488             crosspoint[0]=wx1*pnt1[0] + wx2*pnt2[0];
1489             crosspoint[1]=wy1*pnt1[1] + wy2*pnt2[1];
1490             crosspoint[2]=wz1*pnt1[2] + wz2*pnt2[2];
1491
1492             ncombi++;
1493             for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
1494             for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
1495           }
1496         }
1497       }
1498       delete line2;
1499     }
1500     delete line1;
1501   }
1502   if(ncombi>0){
1503     for(Int_t jj=0;jj<3;jj++){
1504       initPos[jj] = aver[jj]/ncombi;
1505       //printf("%f\n",initPos[jj]);
1506       aversq[jj]/=ncombi;
1507       sigmasq[jj]=aversq[jj]-initPos[jj]*initPos[jj];
1508       sigma+=sigmasq[jj];
1509     }
1510     sigma=TMath::Sqrt(TMath::Abs(sigma));
1511   }
1512   else {
1513     Warning("VertexFinder","Finder did not succed");
1514     sigma=999;
1515   }
1516   fVert.SetXYZ(initPos);
1517   fVert.SetDispersion(sigma);
1518   fVert.SetNContributors(ncombi);
1519 }
1520 //---------------------------------------------------------------------------
1521 void AliVertexerTracks::VertexFitter(Bool_t vfit, Bool_t chiCalc,Int_t useWeights)
1522 {
1523 //
1524 // The optimal estimate of the vertex position is given by a "weighted
1525 // average of tracks positions".
1526 // Original method: V. Karimaki, CMS Note 97/0051
1527 //
1528   const double kTiny = 1e-9;
1529   Bool_t useConstraint = fConstraint;
1530   Double_t initPos[3];
1531   if(!fOnlyFitter) {
1532     fVert.GetXYZ(initPos);
1533   } else {
1534     initPos[0]=fNominalPos[0];
1535     initPos[1]=fNominalPos[1];
1536     initPos[2]=fNominalPos[2];
1537   }
1538
1539   Int_t nTrksSel = (Int_t)fTrkArraySel.GetEntries();
1540   if(nTrksSel==1) useConstraint=kTRUE;
1541   AliDebug(1,Form("--- VertexFitter(): start (%d,%d,%d)",vfit,chiCalc,useWeights));
1542   AliDebug(1,Form(" Number of tracks in array: %d\n",nTrksSel));
1543   AliDebug(1,Form(" Minimum # tracks required in fit: %d\n",fMinTracks));
1544   AliDebug(1,Form(" Vertex position after finder: %f,%f,%f\n",initPos[0],initPos[1],initPos[2]));
1545   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])));
1546
1547   // special treatment for few-tracks fits (e.g. secondary vertices)
1548   Bool_t uUi3by3 = kFALSE; if(nTrksSel<5 && !useConstraint && !useWeights) uUi3by3 = kTRUE;
1549
1550   Int_t i,j,k,step=0;
1551   static TMatrixD rv(3,1);
1552   static TMatrixD vV(3,3);
1553   rv(0,0) = initPos[0];
1554   rv(1,0) = initPos[1];
1555   rv(2,0) = initPos[2];
1556   Double_t xlStart,alpha;
1557   Int_t nTrksUsed = 0;
1558   Double_t chi2=0,chi2i,chi2b;
1559   AliExternalTrackParam *t = 0;
1560   Int_t failed = 0;
1561
1562   // initial vertex covariance matrix
1563   TMatrixD vVb(3,3);
1564   vVb(0,0) = fNominalCov[0];
1565   vVb(0,1) = fNominalCov[1];
1566   vVb(0,2) = 0.;
1567   vVb(1,0) = fNominalCov[1];
1568   vVb(1,1) = fNominalCov[2];
1569   vVb(1,2) = 0.;
1570   vVb(2,0) = 0.;
1571   vVb(2,1) = 0.;
1572   vVb(2,2) = fNominalCov[5];
1573   TMatrixD vVbInv(TMatrixD::kInverted,vVb);
1574   TMatrixD rb(3,1);
1575   rb(0,0) = fNominalPos[0];
1576   rb(1,0) = fNominalPos[1];
1577   rb(2,0) = fNominalPos[2];
1578   TMatrixD vVbInvrb(vVbInv,TMatrixD::kMult,rb);
1579   //
1580   int currBC = fVert.GetBC();
1581   //
1582   // 2 steps:
1583   // 1st - estimate of vtx using all tracks
1584   // 2nd - estimate of global chi2
1585   //
1586   for(step=0; step<2; step++) {
1587     if      (step==0 && !vfit)    continue;
1588     else if (step==1 && !chiCalc) continue;
1589     chi2 = 0.;
1590     nTrksUsed = 0;
1591     fMVWSum = fMVWE2 = 0;
1592     if(step==1) { initPos[0]=rv(0,0); initPos[1]=rv(1,0); initPos[2]=rv(2,0);}
1593     AliDebug(2,Form("Step%d: inipos: %+f %+f %+f MinTr: %d, Sig2:%.2f)",step,initPos[0],initPos[1],initPos[2],fMinTracks,fMVSigma2));
1594     //
1595     TMatrixD sumWiri(3,1);
1596     TMatrixD sumWi(3,3);
1597     for(i=0; i<3; i++) {
1598       sumWiri(i,0) = 0.;
1599       for(j=0; j<3; j++) sumWi(j,i) = 0.;
1600     }
1601
1602     // mean vertex constraint
1603     if(useConstraint) {
1604       for(i=0;i<3;i++) {
1605         sumWiri(i,0) += vVbInvrb(i,0);
1606         for(k=0;k<3;k++) sumWi(i,k) += vVbInv(i,k);
1607       }
1608       // chi2
1609       TMatrixD deltar = rv; deltar -= rb;
1610       TMatrixD vVbInvdeltar(vVbInv,TMatrixD::kMult,deltar);
1611       chi2b = deltar(0,0)*vVbInvdeltar(0,0)+
1612               deltar(1,0)*vVbInvdeltar(1,0)+
1613               deltar(2,0)*vVbInvdeltar(2,0);
1614       chi2 += chi2b;
1615     }
1616
1617     // loop on tracks
1618     for(k=0; k<nTrksSel; k++) {
1619       //
1620       // get track from track array
1621       t = (AliExternalTrackParam*)fTrkArraySel.At(k);
1622       if (useWeights && t->TestBit(kBitUsed)) continue;
1623       //
1624       int tBC = int(t->GetUniqueID()) - kTOFBCShift;    // BC assigned to this track
1625       if (fSelectOnTOFBunchCrossing) {
1626         if (!fKeepAlsoUnflaggedTOFBunchCrossing) continue;   // don't consider tracks with undefined BC
1627         if (currBC!=AliVTrack::kTOFBCNA && tBC!=AliVTrack::kTOFBCNA && tBC!=currBC) continue;  // track does not match to current BCid
1628       }
1629       alpha = t->GetAlpha();
1630       xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha);
1631       // to vtxSeed (from finder)
1632       t->PropagateTo(xlStart,GetFieldkG());
1633
1634       // vector of track global coordinates
1635       TMatrixD ri(3,1);
1636       // covariance matrix of ri
1637       TMatrixD wWi(3,3);
1638
1639       // get space point from track
1640       if(!TrackToPoint(t,ri,wWi,uUi3by3)) continue;
1641
1642       // track chi2
1643       TMatrixD deltar = rv; deltar -= ri;
1644       TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
1645       chi2i = deltar(0,0)*wWideltar(0,0)+
1646               deltar(1,0)*wWideltar(1,0)+
1647               deltar(2,0)*wWideltar(2,0);
1648       //
1649       if (useWeights) {
1650         //double sg = TMath::Sqrt(fMVSigma2);
1651         //double chi2iw = (deltar(0,0)*wWideltar(0,0)+deltar(1,0)*wWideltar(1,0))/sg + deltar(2,0)*wWideltar(2,0)/fMVSigma2;
1652         //double wgh = (1-chi2iw/fMVTukey2);
1653         double chi2iw = chi2i;
1654         double wgh = (1-chi2iw/fMVTukey2/fMVSigma2);
1655
1656         if (wgh<kTiny) wgh = 0;
1657         else if (useWeights==2) wgh = 1.; // use as binary weight
1658         if (step==1) ((AliESDtrack*)t)->SetBit(kBitUsed, wgh>0);
1659         if (wgh<kTiny) continue; // discard the track
1660         wWi *= wgh;  // RS: use weight?
1661         fMVWSum += wgh;
1662         fMVWE2  += wgh*chi2iw;
1663       }
1664       // add to total chi2
1665       if (fSelectOnTOFBunchCrossing && tBC!=AliVTrack::kTOFBCNA && currBC<0) currBC = tBC;
1666       //
1667       chi2 += chi2i;
1668       TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
1669       sumWiri += wWiri;
1670       sumWi   += wWi;
1671
1672       nTrksUsed++;
1673     } // end loop on tracks
1674
1675     if(nTrksUsed < fMinTracks) {
1676       failed=1;
1677       continue;
1678     }
1679
1680     Double_t determinant = sumWi.Determinant();
1681     if(determinant < fMinDetFitter)  {
1682       AliDebug(1,Form("det(V) = %f (<%f)\n",determinant,fMinDetFitter));
1683       failed=1;
1684       continue;
1685     }
1686
1687     if(step==0) {
1688       // inverted of weights matrix
1689       TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
1690       vV = invsumWi;
1691       // position of primary vertex
1692       rv.Mult(vV,sumWiri);
1693     }
1694   } // end loop on the 2 steps
1695
1696   if(failed) {
1697     fVert.SetNContributors(-1);
1698     if (chiCalc) {
1699       TooFewTracks();
1700       if (fCurrentVertex) fVert = *fCurrentVertex;  // RS
1701     }
1702     return;
1703   }
1704   //
1705   Double_t position[3];
1706   position[0] = rv(0,0);
1707   position[1] = rv(1,0);
1708   position[2] = rv(2,0);
1709   Double_t covmatrix[6];
1710   covmatrix[0] = vV(0,0);
1711   covmatrix[1] = vV(0,1);
1712   covmatrix[2] = vV(1,1);
1713   covmatrix[3] = vV(0,2);
1714   covmatrix[4] = vV(1,2);
1715   covmatrix[5] = vV(2,2);
1716
1717   // for correct chi2/ndf, count constraint as additional "track"
1718   if(fConstraint) nTrksUsed++;
1719   //
1720   if (vfit && !chiCalc) { // RS: special mode for multi-vertex finder
1721     fVert.SetXYZ(position);
1722     fVert.SetCovarianceMatrix(covmatrix);
1723     fVert.SetNContributors(nTrksUsed);
1724     return;
1725   }
1726   // store data in the vertex object
1727   if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
1728   fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);
1729   fCurrentVertex->SetBC(currBC);
1730   fVert = *fCurrentVertex;  // RS
1731   AliDebug(1," Vertex after fit:");
1732   AliDebug(1,Form("xyz: %f %f %f; nc %d",fCurrentVertex->GetXv(),fCurrentVertex->GetYv(),fCurrentVertex->GetZv(),fCurrentVertex->GetNContributors()));
1733   AliDebug(1,"--- VertexFitter(): finish\n");
1734
1735
1736   return;
1737 }
1738 //----------------------------------------------------------------------------
1739 AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(const TObjArray *trkArray,
1740                                                          UShort_t *id,
1741                                                          Bool_t optUseFitter,
1742                                                          Bool_t optPropagate,
1743                                                          Bool_t optUseDiamondConstraint)
1744 {
1745 //
1746 // Return vertex from tracks (AliExternalTrackParam) in array
1747 //
1748   fCurrentVertex = 0;
1749   // set optUseDiamondConstraint=TRUE only if you are reconstructing the
1750   // primary vertex!
1751   if(optUseDiamondConstraint) {
1752     SetConstraintOn();
1753   } else {
1754     SetConstraintOff();
1755   }
1756
1757   // get tracks and propagate them to initial vertex position
1758   fIdSel = new UShort_t[(Int_t)trkArray->GetEntriesFast()];
1759   Int_t nTrksSel = PrepareTracks(*trkArray,id,0);
1760   if((!optUseDiamondConstraint && nTrksSel<TMath::Max(2,fMinTracks)) ||
1761      (optUseDiamondConstraint && nTrksSel<1)) {
1762     TooFewTracks();
1763     return fCurrentVertex;
1764   }
1765
1766   // vertex finder
1767   if(nTrksSel==1) {
1768     AliDebug(1,"Just one track");
1769     OneTrackVertFinder();
1770   } else {
1771     switch (fAlgo) {
1772       case 1: StrLinVertexFinderMinDist(1); break;
1773       case 2: StrLinVertexFinderMinDist(0); break;
1774       case 3: HelixVertexFinder();          break;
1775       case 4: VertexFinder(1);              break;
1776       case 5: VertexFinder(0);              break;
1777       default: printf("Wrong algorithm\n"); break;
1778     }
1779   }
1780   AliDebug(1," Vertex finding completed\n");
1781
1782   // vertex fitter
1783   if(optUseFitter) {
1784     VertexFitter();
1785   } else {
1786     Double_t position[3]={fVert.GetXv(),fVert.GetYv(),fVert.GetZv()};
1787     Double_t covmatrix[6];
1788     fVert.GetCovMatrix(covmatrix);
1789     Double_t chi2=99999.;
1790     Int_t    nTrksUsed=fVert.GetNContributors();
1791     fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);
1792   }
1793   fCurrentVertex->SetDispersion(fVert.GetDispersion());
1794
1795
1796   // set indices of used tracks and propagate track to found vertex
1797   UShort_t *indices = 0;
1798   Double_t d0z0[2],covd0z0[3];
1799   AliExternalTrackParam *t = 0;
1800   if(fCurrentVertex->GetNContributors()>0) {
1801     indices = new UShort_t[fTrkArraySel.GetEntriesFast()];
1802     for(Int_t jj=0; jj<(Int_t)fTrkArraySel.GetEntriesFast(); jj++) {
1803       indices[jj] = fIdSel[jj];
1804       t = (AliExternalTrackParam*)fTrkArraySel.At(jj);
1805       if(optPropagate && optUseFitter) {
1806         if(TMath::Sqrt(fCurrentVertex->GetXv()*fCurrentVertex->GetXv()+fCurrentVertex->GetYv()*fCurrentVertex->GetYv())<3.) {
1807           t->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
1808           AliDebug(1,Form("Track %d propagated to found vertex",jj));
1809         } else {
1810           AliWarning("Found vertex outside beam pipe!");
1811         }
1812       }
1813     }
1814     fCurrentVertex->SetIndices(fCurrentVertex->GetNContributors(),indices);
1815   }
1816
1817   // clean up
1818   if (indices) {delete [] indices; indices=NULL;}
1819   delete [] fIdSel; fIdSel=NULL;
1820   fTrkArraySel.Delete();
1821
1822   return fCurrentVertex;
1823 }
1824
1825 //----------------------------------------------------------------------------
1826 AliESDVertex* AliVertexerTracks::VertexForSelectedESDTracks(TObjArray *trkArray,
1827                                                              Bool_t optUseFitter,
1828                                                             Bool_t optPropagate,
1829                                                             Bool_t optUseDiamondConstraint)
1830
1831 {
1832 //
1833 // Return vertex from array of ESD tracks
1834 //
1835
1836   Int_t nTrks = (Int_t)trkArray->GetEntriesFast();
1837   UShort_t *id = new UShort_t[nTrks];
1838
1839   AliESDtrack *esdt = 0;
1840   for(Int_t i=0; i<nTrks; i++){
1841     esdt = (AliESDtrack*)trkArray->At(i);
1842     id[i] = (UShort_t)esdt->GetID();
1843   }
1844
1845   VertexForSelectedTracks(trkArray,id,optUseFitter,optPropagate,optUseDiamondConstraint);
1846
1847   delete [] id; id=NULL;
1848
1849   return fCurrentVertex;
1850 }
1851
1852 //______________________________________________________
1853 Bool_t AliVertexerTracks::FindNextVertexMV()
1854 {
1855   // try to find a new vertex
1856   fMVSigma2 = fMVSig2Ini;
1857   double prevSig2 = fMVSigma2;
1858   double minDst2 = fMVMinDst*fMVMinDst;
1859   const double kSigLimit = 1.0;
1860   const double kSigLimitE = kSigLimit+1e-6;
1861   const double kPushFactor = 0.5;
1862   const int kMaxIter = 20;
1863   double push = kPushFactor;
1864   //
1865   int iter = 0;
1866   double posP[3]={0,0,0},pos[3]={0,0,0};
1867   fVert.GetXYZ(posP);
1868   //
1869
1870   do {
1871     fVert.SetBC(AliVTrack::kTOFBCNA);
1872     VertexFitter(kTRUE,kFALSE,1);
1873     if (fVert.GetNContributors()<fMinTracks) {
1874       AliDebug(3,Form("Failed in iteration %d: No Contributirs",iter));
1875       break;
1876     } // failed
1877     if (fMVWSum>0) fMVSigma2 = TMath::Max(fMVWE2/fMVWSum,kSigLimit);
1878     else {
1879       AliDebug(3,Form("Failed %d, no weithgs",iter));
1880       iter = kMaxIter+1;
1881       break;
1882     } // failed
1883     //
1884     double sigRed = (prevSig2-fMVSigma2)/prevSig2; // reduction of sigma2
1885     //
1886     fVert.GetXYZ(pos);
1887     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]);
1888     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));
1889     if ( (++iter<kMaxIter) && (sigRed<0 || sigRed<fMVMinSig2Red) && fMVSigma2>fMVMaxSigma2) {
1890       fMVSigma2 *= push; // stuck, push little bit
1891       push *= kPushFactor;
1892       if (fMVSigma2<1.) fMVSigma2 = 1.;
1893       AliDebug(3,Form("Pushed sigma2 to %f",fMVSigma2));
1894     }
1895     else if (dst2<minDst2 && ((sigRed<0 || sigRed<fMVMinSig2Red) || fMVSigma2<kSigLimitE)) break;
1896     //
1897     fVert.GetXYZ(posP); // fetch previous vertex position
1898     prevSig2 = fMVSigma2;
1899   } while(iter<kMaxIter);
1900   //
1901   if (fVert.GetNContributors()<0 || iter>kMaxIter || fMVSigma2>fMVMaxSigma2) {
1902     return kFALSE;
1903   }
1904   else {
1905     VertexFitter(kFALSE,kTRUE,fMVFinalWBinary ? 2:1); // final chi2 calculation
1906     int nv = fMVVertices->GetEntries();
1907     // create indices
1908     int ntrk = fTrkArraySel.GetEntries();
1909     int nindices = fCurrentVertex->GetNContributors() - (fConstraint ? 1:0);
1910     UShort_t *indices = 0;
1911     if (nindices>0) indices = new UShort_t[nindices];
1912     int nadded = 0;
1913     for (int itr=0;itr<ntrk;itr++) {
1914       AliExternalTrackParam* t = (AliExternalTrackParam*)fTrkArraySel[itr];
1915       if (t->TestBit(kBitAccounted) || !t->TestBit(kBitUsed)) continue;   // already belongs to some vertex
1916       t->SetBit(kBitAccounted);
1917       indices[nadded++] = fIdSel[itr];
1918     }
1919     if (nadded!=nindices) printf("Mismatch : NInd: %d Nadd: %d\n",nindices,nadded);
1920     fCurrentVertex->SetIndices(nindices,indices);
1921     // set vertex title
1922     TString title="VertexerTracksMVNoConstraint";
1923     if(fConstraint) title="VertexerTracksMVWithConstraint";
1924     fCurrentVertex->SetTitle(title.Data());
1925     fMVVertices->AddLast(fCurrentVertex);
1926     AliDebug(3,Form("Added new vertex #%d NCont:%d XYZ: %f %f %f",nindices,nv,fCurrentVertex->GetX(),fCurrentVertex->GetY(),fCurrentVertex->GetZ()));
1927     if (indices) delete[] indices;
1928     fCurrentVertex = 0; // already attached to fMVVertices
1929     return kTRUE;
1930   }
1931 }
1932
1933 //______________________________________________________
1934 void AliVertexerTracks::FindVerticesMV()
1935 {
1936   // find and fit multiple vertices
1937   //
1938   double step = fMVScanStep>1 ?  fMVScanStep : 1.;
1939   double zmx = 3*TMath::Sqrt(fNominalCov[5]);
1940   double zmn = -zmx;
1941   int nz = TMath::Nint((zmx-zmn)/step); if (nz<1) nz=1;
1942   double dz = (zmx-zmn)/nz;
1943   int izStart=0;
1944   AliDebug(2,Form("%d seeds between %f and %f",nz,zmn+dz/2,zmx+dz/2));
1945   //
1946   if (!fMVVertices) fMVVertices = new TObjArray(10);
1947   fMVVertices->Clear();
1948   //
1949   int ntrLeft = (Int_t)fTrkArraySel.GetEntries();
1950   //
1951   double sig2Scan = fMVSig2Ini;
1952   Bool_t runMore = kTRUE;
1953   int cntWide = 0;
1954   while (runMore) {
1955     fMVSig2Ini = sig2Scan*1e3;  // try wide search
1956     Bool_t found = kFALSE;
1957     cntWide++;
1958     fVert.SetNContributors(-1);
1959     fVert.SetXYZ(fNominalPos);
1960     AliDebug(3,Form("Wide search #%d Z= %f Sigma2=%f",cntWide,fNominalPos[2],fMVSig2Ini));
1961     if (FindNextVertexMV()) { // are there tracks left to consider?
1962       AliESDVertex* vtLast = (AliESDVertex*)fMVVertices->Last();
1963       if (vtLast && vtLast->GetNContributors()>0) ntrLeft -= vtLast->GetNContributors()-(fConstraint ? 1:0);
1964       if (ntrLeft<1) runMore = kFALSE;
1965       found = kTRUE;
1966       continue;
1967     }
1968     // if nothing is found, do narrow sig2ini scan
1969     fMVSig2Ini = sig2Scan;
1970     for (;izStart<nz;izStart++) {
1971       double zSeed = zmn+dz*(izStart+0.5);
1972       AliDebug(3,Form("Seed %d: Z= %f Sigma2=%f",izStart,zSeed,fMVSig2Ini));
1973       fVert.SetNContributors(-1);
1974       fVert.SetXYZ(fNominalPos);
1975       fVert.SetZv(zSeed);
1976       if (FindNextVertexMV()) { // are there tracks left to consider?
1977         AliESDVertex* vtLast = (AliESDVertex*)fMVVertices->Last();
1978         if (vtLast && vtLast->GetNContributors()>0) ntrLeft -= vtLast->GetNContributors()-(fConstraint ? 1:0);
1979         if (ntrLeft<1) runMore = kFALSE;
1980         found = kTRUE;
1981         break;
1982       }
1983     }
1984     runMore = found; // if nothing was found, no need for new iteration
1985   }
1986   fMVSig2Ini = sig2Scan;
1987   int nvFound = fMVVertices->GetEntriesFast();
1988   AliDebug(2,Form("Number of found vertices: %d",nvFound));
1989   if (nvFound<1) TooFewTracks();
1990   if (AliLog::GetGlobalDebugLevel()>0) fMVVertices->Print();
1991   //
1992 }
1993
1994 //______________________________________________________
1995 void AliVertexerTracks::AnalyzePileUp(AliESDEvent* esdEv)
1996 {
1997   // if multiple vertices are found, try to find the primary one and attach it as fCurrentVertex
1998   // then attach pile-up vertices directly to esdEv
1999   //
2000   int nFND = (fMVVertices && fMVVertices->GetEntriesFast()) ? fMVVertices->GetEntriesFast() : 0;
2001   if (nFND<1) { if (!fCurrentVertex) TooFewTracks(); return;} // no multiple vertices
2002   //
2003   int indCont[nFND];
2004   int nIndx[nFND];
2005   for (int iv=0;iv<nFND;iv++) {
2006     AliESDVertex* fnd = (AliESDVertex*)fMVVertices->At(iv);
2007     indCont[iv] = iv;
2008     nIndx[iv]   = fnd->GetNIndices();
2009   }
2010   TMath::Sort(nFND, nIndx, indCont, kTRUE); // sort in decreasing order of Nindices
2011   double dists[nFND];
2012   int    distOrd[nFND],indx[nFND];
2013   for (int iv=0;iv<nFND;iv++) {
2014     AliESDVertex* fndI = (AliESDVertex*)fMVVertices->At(indCont[iv]);
2015     if (fndI->GetStatus()<1) continue; // discarded
2016     int ncomp = 0;
2017     for (int jv=iv+1;jv<nFND;jv++) {
2018       AliESDVertex* fndJ = (AliESDVertex*)fMVVertices->At(indCont[jv]);
2019       if (fndJ->GetStatus()<1) continue;
2020       dists[ncomp] = fndI->GetWDist(fndJ)*fndJ->GetNIndices();
2021       distOrd[ncomp] = indCont[jv];
2022       indx[ncomp]  = ncomp;
2023       ncomp++;
2024     }
2025     if (ncomp<1) continue;
2026     TMath::Sort(ncomp, dists, indx, kFALSE); // sort in increasing distance order
2027     for (int jv0=0;jv0<ncomp;jv0++) {
2028       int jv = distOrd[indx[jv0]];
2029       AliESDVertex* fndJ = (AliESDVertex*)fMVVertices->At(jv);
2030       if (dists[indx[jv0]]<fMVMaxWghNtr) { // candidate for split vertex
2031         //before eliminating the small close vertex, check if we should transfere its BC to largest one
2032         if (fndJ->GetBC()!=AliVTrack::kTOFBCNA && fndI->GetBC()==AliVTrack::kTOFBCNA) fndI->SetBC(fndJ->GetBC());
2033         //
2034         // leave the vertex only if both vertices have definit but different BC's, then this is not splitting.
2035         if ( (fndJ->GetBC()==fndI->GetBC()) || (fndJ->GetBC()==AliVTrack::kTOFBCNA)) fndJ->SetNContributors(-fndJ->GetNContributors());
2036       }
2037     }
2038   }
2039   //
2040   // select as a primary the largest vertex with BC=0, or the largest with BC non-ID
2041   int primBC0=-1,primNoBC=-1;
2042   for (int iv0=0;iv0<nFND;iv0++) {
2043     int iv = indCont[iv0];
2044     AliESDVertex* fndI = (AliESDVertex*)fMVVertices->At(iv);
2045     if (!fndI) continue;
2046     if (fndI->GetStatus()<1) {fMVVertices->RemoveAt(iv); delete fndI; continue;}
2047     if (primBC0<0  && fndI->GetBC()==0) primBC0 = iv;
2048     if (primNoBC<0 && fndI->GetBC()==AliVTrack::kTOFBCNA)  primNoBC = iv;
2049   }
2050   //
2051   if (primBC0>=0) fCurrentVertex = (AliESDVertex*)fMVVertices->At(primBC0);
2052   if (!fCurrentVertex && primNoBC>=0) fCurrentVertex = (AliESDVertex*)fMVVertices->At(primNoBC);
2053   if (fCurrentVertex) fMVVertices->Remove(fCurrentVertex);
2054   else {  // all vertices have BC>0, no primary vertex
2055     fCurrentVertex = new AliESDVertex();
2056     fCurrentVertex->SetNContributors(-3);
2057     fCurrentVertex->SetBC(AliVTrack::kTOFBCNA);
2058   }
2059   fCurrentVertex->SetID(-1);
2060   //
2061   // add pileup vertices
2062   int nadd = 0;
2063   for (int iv0=0;iv0<nFND;iv0++) {
2064     int iv = indCont[iv0];
2065     AliESDVertex* fndI = (AliESDVertex*)fMVVertices->At(iv);
2066     if (!fndI) continue;
2067     fndI->SetID(++nadd);
2068     esdEv->AddPileupVertexTracks(fndI);
2069   }
2070   //
2071   fMVVertices->Delete();
2072   //
2073 }
2074