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