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