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: {AliFatal(Form("Wrong seeder algorithm %d",fAlgoIter0));} 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: {AliFatal(Form("Wrong vertexer algorithm %d",fAlgo));} 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; break;}
 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)) {
1153       optUseWeights=kFALSE;
1154       AliDebug(1,"WARNING: TrackToPoint failed");
1155     }
1156     Double_t wmat[9];
1157     Int_t iel=0;
1158     for(Int_t ia=0;ia<3;ia++){
1159       for(Int_t ib=0;ib<3;ib++){
1160         wmat[iel]=wWi(ia,ib);
1161         iel++;
1162       }
1163     }
1164     linarray[i] = new AliStrLine(pos,sigmasq,wmat,dir);
1165   }
1166   fVert=TrackletVertexFinder(linarray,knacc,optUseWeights);
1167   for(Int_t i=0; i<knacc; i++) delete linarray[i];
1168   delete [] linarray;
1169 }
1170 //---------------------------------------------------------------------------
1171 AliESDVertex AliVertexerTracks::TrackletVertexFinder(const TClonesArray *lines, Int_t optUseWeights)
1172 {
1173   // Calculate the point at minimum distance to prepared tracks (TClonesArray)
1174   const Int_t knacc = (Int_t)lines->GetEntriesFast();
1175   AliStrLine** lines2 = new AliStrLine* [knacc];
1176   for(Int_t i=0; i<knacc; i++){
1177     lines2[i]= (AliStrLine*)lines->At(i);
1178   }
1179   AliESDVertex vert = TrackletVertexFinder(lines2,knacc,optUseWeights);
1180   delete [] lines2;
1181   return vert;
1182 }
1183
1184 //---------------------------------------------------------------------------
1185 AliESDVertex AliVertexerTracks::TrackletVertexFinder(AliStrLine **lines, const Int_t knacc, Int_t optUseWeights)
1186 {
1187   // Calculate the point at minimum distance to prepared tracks (array of AliStrLine)
1188
1189   Double_t initPos[3]={0.,0.,0.};
1190
1191   Double_t (*vectP0)[3]=new Double_t [knacc][3];
1192   Double_t (*vectP1)[3]=new Double_t [knacc][3];
1193
1194   Double_t sum[3][3];
1195   Double_t dsum[3]={0,0,0};
1196   TMatrixD sumWi(3,3);
1197   for(Int_t i=0;i<3;i++){
1198     for(Int_t j=0;j<3;j++){
1199       sum[i][j]=0;
1200       sumWi(i,j)=0.;
1201     }
1202   }
1203
1204   for(Int_t i=0; i<knacc; i++){
1205     AliStrLine *line1 = lines[i];
1206     Double_t p0[3],cd[3],sigmasq[3];
1207     Double_t wmat[9];
1208     if(!line1) {printf("ERROR %d %d\n",i,knacc); continue;}
1209     line1->GetP0(p0);
1210     line1->GetCd(cd);
1211     line1->GetSigma2P0(sigmasq);
1212     line1->GetWMatrix(wmat);
1213     TMatrixD wWi(3,3);
1214     Int_t iel=0;
1215     for(Int_t ia=0;ia<3;ia++){
1216       for(Int_t ib=0;ib<3;ib++){
1217         wWi(ia,ib)=wmat[iel];
1218         iel++;
1219       }
1220     }
1221
1222     sumWi+=wWi;
1223
1224     Double_t p1[3]={p0[0]+cd[0],p0[1]+cd[1],p0[2]+cd[2]};
1225     vectP0[i][0]=p0[0];
1226     vectP0[i][1]=p0[1];
1227     vectP0[i][2]=p0[2];
1228     vectP1[i][0]=p1[0];
1229     vectP1[i][1]=p1[1];
1230     vectP1[i][2]=p1[2];
1231
1232     Double_t matr[3][3];
1233     Double_t dknow[3];
1234     if(optUseWeights==0)GetStrLinDerivMatrix(p0,p1,matr,dknow);
1235     else GetStrLinDerivMatrix(p0,p1,sigmasq,matr,dknow);
1236
1237
1238     for(Int_t iii=0;iii<3;iii++){
1239       dsum[iii]+=dknow[iii];
1240       for(Int_t lj=0;lj<3;lj++) sum[iii][lj]+=matr[iii][lj];
1241     }
1242   }
1243
1244   TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
1245   Double_t covmatrix[6];
1246   covmatrix[0] = invsumWi(0,0);
1247   covmatrix[1] = invsumWi(0,1);
1248   covmatrix[2] = invsumWi(1,1);
1249   covmatrix[3] = invsumWi(0,2);
1250   covmatrix[4] = invsumWi(1,2);
1251   covmatrix[5] = invsumWi(2,2);
1252
1253   Double_t vett[3][3];
1254   Double_t det=GetDeterminant3X3(sum);
1255   Double_t sigma=0;
1256
1257   if(TMath::Abs(det) > kAlmost0){
1258     for(Int_t zz=0;zz<3;zz++){
1259       for(Int_t ww=0;ww<3;ww++){
1260         for(Int_t kk=0;kk<3;kk++) vett[ww][kk]=sum[ww][kk];
1261       }
1262       for(Int_t kk=0;kk<3;kk++) vett[kk][zz]=dsum[kk];
1263       initPos[zz]=GetDeterminant3X3(vett)/det;
1264     }
1265
1266
1267     for(Int_t i=0; i<knacc; i++){
1268       Double_t p0[3]={0,0,0},p1[3]={0,0,0};
1269       for(Int_t ii=0;ii<3;ii++){
1270         p0[ii]=vectP0[i][ii];
1271         p1[ii]=vectP1[i][ii];
1272       }
1273       sigma+=GetStrLinMinDist(p0,p1,initPos);
1274     }
1275
1276     if(sigma>0.) {sigma=TMath::Sqrt(sigma);}else{sigma=999;}
1277   }else{
1278     sigma=999;
1279   }
1280   AliESDVertex theVert(initPos,covmatrix,99999.,knacc);
1281   theVert.SetDispersion(sigma);
1282   delete [] vectP0;
1283   delete [] vectP1;
1284   return theVert;
1285 }
1286
1287 //---------------------------------------------------------------------------
1288 Bool_t AliVertexerTracks::TrackToPoint(AliExternalTrackParam *t,
1289                                        TMatrixD &ri,TMatrixD &wWi,
1290                                        Bool_t uUi3by3) const
1291 {
1292 //
1293 // Extract from the AliExternalTrackParam the global coordinates ri and covariance matrix
1294 // wWi of the space point that it represents (to be used in VertexFitter())
1295 //
1296
1297
1298   Double_t rotAngle = t->GetAlpha();
1299   if(rotAngle<0.) rotAngle += 2.*TMath::Pi();
1300   Double_t cosRot = TMath::Cos(rotAngle);
1301   Double_t sinRot = TMath::Sin(rotAngle);
1302   /*
1303   // RS >>>
1304   Double_t lambda = TMath::ATan(t->GetTgl());
1305   Double_t cosLam   = TMath::Cos(lambda);
1306   Double_t sinLam   = TMath::Sin(lambda);
1307   // RS <<<
1308   */
1309   ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot;
1310   ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot;
1311   ri(2,0) = t->GetZ();
1312
1313   if(!uUi3by3) {
1314     // matrix to go from global (x,y,z) to local (y,z);
1315     TMatrixD qQi(2,3);
1316     qQi(0,0) = -sinRot;
1317     qQi(0,1) = cosRot;
1318     qQi(0,2) = 0.;
1319     //
1320     qQi(1,0) = 0.;
1321     qQi(1,1) = 0.;
1322     qQi(1,2) = 1.;
1323     //
1324     // RS: Added polar inclination
1325     /*
1326     qQi(1,0) = -sinLam*cosRot;
1327     qQi(1,1) = -sinLam*sinRot;
1328     qQi(1,2) = cosLam;
1329     */
1330     // covariance matrix of local (y,z) - inverted
1331     TMatrixD uUi(2,2);
1332     uUi(0,0) = t->GetSigmaY2();
1333     uUi(0,1) = t->GetSigmaZY();
1334     uUi(1,0) = t->GetSigmaZY();
1335     uUi(1,1) = t->GetSigmaZ2();
1336     //printf(" Ui :");
1337     //printf(" %f   %f",uUi(0,0),uUi(0,1));
1338     //printf(" %f   %f",uUi(1,0),uUi(1,1));
1339
1340     if(uUi.Determinant() <= 0.) return kFALSE;
1341     TMatrixD uUiInv(TMatrixD::kInverted,uUi);
1342
1343     // weights matrix: wWi = qQiT * uUiInv * qQi
1344     TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
1345     TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
1346     wWi = m;
1347   } else {
1348     if(fVert.GetNContributors()<1) AliFatal("Vertex from finder is empty");
1349     // matrix to go from global (x,y,z) to local (x,y,z);
1350     TMatrixD qQi(3,3);
1351     qQi(0,0) = cosRot;
1352     qQi(0,1) = sinRot;
1353     qQi(0,2) = 0.;
1354     qQi(1,0) = -sinRot;
1355     qQi(1,1) = cosRot;
1356     qQi(1,2) = 0.;
1357     qQi(2,0) = 0.;
1358     qQi(2,1) = 0.;
1359     qQi(2,2) = 1.;
1360
1361     // covariance of fVert along the track
1362     Double_t p[3],pt,ptot;
1363     t->GetPxPyPz(p);
1364     pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
1365     ptot = TMath::Sqrt(pt*pt+p[2]*p[2]);
1366     Double_t cphi = p[0]/pt;               //cos(phi)=px/pt
1367     Double_t sphi = p[1]/pt;               //sin(phi)=py/pt
1368     Double_t clambda = pt/ptot;            //cos(lambda)=pt/ptot
1369     Double_t slambda = p[2]/ptot;            //sin(lambda)=pz/ptot
1370     Double_t covfVert[6];
1371     fVert.GetCovMatrix(covfVert);
1372     Double_t covfVertalongt =
1373        covfVert[0]*cphi*cphi*clambda*clambda
1374       +covfVert[1]*2.*cphi*sphi*clambda*clambda
1375       +covfVert[3]*2.*cphi*clambda*slambda
1376       +covfVert[2]*sphi*sphi*clambda*clambda
1377       +covfVert[4]*2.*sphi*clambda*slambda
1378       +covfVert[5]*slambda*slambda;
1379     // covariance matrix of local (x,y,z) - inverted
1380     TMatrixD uUi(3,3);
1381     uUi(0,0) = covfVertalongt * fnSigmaForUi00 * fnSigmaForUi00;
1382     AliDebug(1,Form("=====> sqrtUi00 cm  %f",TMath::Sqrt(uUi(0,0))));
1383     uUi(0,1) = 0.;
1384     uUi(0,2) = 0.;
1385     uUi(1,0) = 0.;
1386     uUi(1,1) = t->GetSigmaY2();
1387     uUi(1,2) = t->GetSigmaZY();
1388     uUi(2,0) = 0.;
1389     uUi(2,1) = t->GetSigmaZY();
1390     uUi(2,2) = t->GetSigmaZ2();
1391     //printf(" Ui :\n");
1392     //printf(" %f   %f\n",uUi(0,0),uUi(0,1));
1393     //printf(" %f   %f\n",uUi(1,0),uUi(1,1));
1394
1395     if(uUi.Determinant() <= 0.) return kFALSE;
1396     TMatrixD uUiInv(TMatrixD::kInverted,uUi);
1397
1398     // weights matrix: wWi = qQiT * uUiInv * qQi
1399     TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
1400     TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
1401     wWi = m;
1402   }
1403
1404
1405   return kTRUE;
1406 }
1407 //---------------------------------------------------------------------------
1408 void AliVertexerTracks::TooFewTracks()
1409 {
1410 //
1411 // When the number of tracks is < fMinTracks,
1412 // deal with vertices not found and prepare to exit
1413 //
1414   AliDebug(1,"TooFewTracks");
1415
1416   Double_t pos[3],err[3];
1417   pos[0] = fNominalPos[0];
1418   err[0] = TMath::Sqrt(fNominalCov[0]);
1419   pos[1] = fNominalPos[1];
1420   err[1] = TMath::Sqrt(fNominalCov[2]);
1421   pos[2] = fNominalPos[2];
1422   err[2] = TMath::Sqrt(fNominalCov[5]);
1423   Int_t    ncontr = (err[0]>1. ? -1 : -3);
1424   if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
1425   fCurrentVertex = new AliESDVertex(pos,err);
1426   fCurrentVertex->SetNContributors(ncontr);
1427
1428   if(fConstraint) {
1429     fCurrentVertex->SetTitle("VertexerTracksWithConstraint");
1430   } else {
1431     fCurrentVertex->SetTitle("VertexerTracksNoConstraint");
1432   }
1433
1434   if(!fTrkArraySel.IsEmpty()) fTrkArraySel.Delete();
1435   if(fIdSel) {delete [] fIdSel; fIdSel=NULL;}
1436   if(fTrksToSkip) {delete [] fTrksToSkip; fTrksToSkip=NULL;}
1437
1438   return;
1439 }
1440 //---------------------------------------------------------------------------
1441 void AliVertexerTracks::VertexFinder(Int_t optUseWeights)
1442 {
1443
1444   // Get estimate of vertex position in (x,y) from tracks DCA
1445
1446   Double_t initPos[3];
1447   initPos[2] = 0.;
1448   for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
1449   Int_t nacc = (Int_t)fTrkArraySel.GetEntriesFast();
1450   Double_t aver[3]={0.,0.,0.};
1451   Double_t aversq[3]={0.,0.,0.};
1452   Double_t sigmasq[3]={0.,0.,0.};
1453   Double_t sigma=0;
1454   Int_t ncombi = 0;
1455   AliExternalTrackParam *track1;
1456   AliExternalTrackParam *track2;
1457   Double_t pos[3],dir[3];
1458   Double_t alpha,mindist;
1459
1460   for(Int_t i=0; i<nacc; i++){
1461     track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
1462     alpha=track1->GetAlpha();
1463     mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
1464     track1->GetXYZAt(mindist,GetFieldkG(),pos);
1465     track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
1466     AliStrLine *line1 = new AliStrLine(pos,dir);
1467
1468    //    AliStrLine *line1 = new AliStrLine();
1469    //    track1->ApproximateHelixWithLine(mindist,GetFieldkG(),line1);
1470
1471     for(Int_t j=i+1; j<nacc; j++){
1472       track2 = (AliExternalTrackParam*)fTrkArraySel.At(j);
1473       alpha=track2->GetAlpha();
1474       mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
1475       track2->GetXYZAt(mindist,GetFieldkG(),pos);
1476       track2->GetPxPyPzAt(mindist,GetFieldkG(),dir);
1477       AliStrLine *line2 = new AliStrLine(pos,dir);
1478     //      AliStrLine *line2 = new AliStrLine();
1479     //  track2->ApproximateHelixWithLine(mindist,GetFieldkG(),line2);
1480       Double_t distCA=line2->GetDCA(line1);
1481       //printf("%d   %d   %f\n",i,j,distCA);
1482        if(fDCAcut<=0 || (fDCAcut>0&&distCA<fDCAcut)){
1483         Double_t pnt1[3],pnt2[3],crosspoint[3];
1484
1485         if(optUseWeights<=0){
1486           Int_t retcode = line2->Cross(line1,crosspoint);
1487           if(retcode>=0){
1488             ncombi++;
1489             for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
1490             for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
1491           }
1492         }
1493         if(optUseWeights>0){
1494           Int_t retcode = line1->CrossPoints(line2,pnt1,pnt2);
1495           if(retcode>=0){
1496             Double_t cs, sn;
1497             alpha=track1->GetAlpha();
1498             cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
1499             Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
1500             alpha=track2->GetAlpha();
1501             cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
1502             Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
1503             Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
1504             Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
1505             Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
1506             Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
1507             crosspoint[0]=wx1*pnt1[0] + wx2*pnt2[0];
1508             crosspoint[1]=wy1*pnt1[1] + wy2*pnt2[1];
1509             crosspoint[2]=wz1*pnt1[2] + wz2*pnt2[2];
1510
1511             ncombi++;
1512             for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
1513             for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
1514           }
1515         }
1516       }
1517       delete line2;
1518     }
1519     delete line1;
1520   }
1521   if(ncombi>0){
1522     for(Int_t jj=0;jj<3;jj++){
1523       initPos[jj] = aver[jj]/ncombi;
1524       //printf("%f\n",initPos[jj]);
1525       aversq[jj]/=ncombi;
1526       sigmasq[jj]=aversq[jj]-initPos[jj]*initPos[jj];
1527       sigma+=sigmasq[jj];
1528     }
1529     sigma=TMath::Sqrt(TMath::Abs(sigma));
1530   }
1531   else {
1532     Warning("VertexFinder","Finder did not succed");
1533     sigma=999;
1534   }
1535   fVert.SetXYZ(initPos);
1536   fVert.SetDispersion(sigma);
1537   fVert.SetNContributors(ncombi);
1538 }
1539 //---------------------------------------------------------------------------
1540 void AliVertexerTracks::VertexFitter(Bool_t vfit, Bool_t chiCalc,Int_t useWeights)
1541 {
1542 //
1543 // The optimal estimate of the vertex position is given by a "weighted
1544 // average of tracks positions".
1545 // Original method: V. Karimaki, CMS Note 97/0051
1546 //
1547   const double kTiny = 1e-9;
1548   Bool_t useConstraint = fConstraint;
1549   Double_t initPos[3];
1550   if(!fOnlyFitter) {
1551     fVert.GetXYZ(initPos);
1552   } else {
1553     initPos[0]=fNominalPos[0];
1554     initPos[1]=fNominalPos[1];
1555     initPos[2]=fNominalPos[2];
1556   }
1557
1558   Int_t nTrksSel = (Int_t)fTrkArraySel.GetEntries();
1559   if(nTrksSel==1) useConstraint=kTRUE;
1560   AliDebug(1,Form("--- VertexFitter(): start (%d,%d,%d)",vfit,chiCalc,useWeights));
1561   AliDebug(1,Form(" Number of tracks in array: %d\n",nTrksSel));
1562   AliDebug(1,Form(" Minimum # tracks required in fit: %d\n",fMinTracks));
1563   AliDebug(1,Form(" Vertex position after finder: %f,%f,%f\n",initPos[0],initPos[1],initPos[2]));
1564   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])));
1565
1566   // special treatment for few-tracks fits (e.g. secondary vertices)
1567   Bool_t uUi3by3 = kFALSE; if(nTrksSel<5 && !useConstraint && !useWeights) uUi3by3 = kTRUE;
1568
1569   Int_t i,j,k,step=0;
1570   static TMatrixD rv(3,1);
1571   static TMatrixD vV(3,3);
1572   rv(0,0) = initPos[0];
1573   rv(1,0) = initPos[1];
1574   rv(2,0) = initPos[2];
1575   Double_t xlStart,alpha;
1576   Int_t nTrksUsed = 0;
1577   Double_t chi2=0,chi2i,chi2b;
1578   AliExternalTrackParam *t = 0;
1579   Int_t failed = 0;
1580
1581   // initial vertex covariance matrix
1582   TMatrixD vVb(3,3);
1583   vVb(0,0) = fNominalCov[0];
1584   vVb(0,1) = fNominalCov[1];
1585   vVb(0,2) = 0.;
1586   vVb(1,0) = fNominalCov[1];
1587   vVb(1,1) = fNominalCov[2];
1588   vVb(1,2) = 0.;
1589   vVb(2,0) = 0.;
1590   vVb(2,1) = 0.;
1591   vVb(2,2) = fNominalCov[5];
1592   TMatrixD vVbInv(TMatrixD::kInverted,vVb);
1593   TMatrixD rb(3,1);
1594   rb(0,0) = fNominalPos[0];
1595   rb(1,0) = fNominalPos[1];
1596   rb(2,0) = fNominalPos[2];
1597   TMatrixD vVbInvrb(vVbInv,TMatrixD::kMult,rb);
1598   //
1599   int currBC = fVert.GetBC();
1600   //
1601   // 2 steps:
1602   // 1st - estimate of vtx using all tracks
1603   // 2nd - estimate of global chi2
1604   //
1605   for(step=0; step<2; step++) {
1606     if      (step==0 && !vfit)    continue;
1607     else if (step==1 && !chiCalc) continue;
1608     chi2 = 0.;
1609     nTrksUsed = 0;
1610     fMVWSum = fMVWE2 = 0;
1611     if(step==1) { initPos[0]=rv(0,0); initPos[1]=rv(1,0); initPos[2]=rv(2,0);}
1612     AliDebug(2,Form("Step%d: inipos: %+f %+f %+f MinTr: %d, Sig2:%.2f)",step,initPos[0],initPos[1],initPos[2],fMinTracks,fMVSigma2));
1613     //
1614     TMatrixD sumWiri(3,1);
1615     TMatrixD sumWi(3,3);
1616     for(i=0; i<3; i++) {
1617       sumWiri(i,0) = 0.;
1618       for(j=0; j<3; j++) sumWi(j,i) = 0.;
1619     }
1620
1621     // mean vertex constraint
1622     if(useConstraint) {
1623       for(i=0;i<3;i++) {
1624         sumWiri(i,0) += vVbInvrb(i,0);
1625         for(k=0;k<3;k++) sumWi(i,k) += vVbInv(i,k);
1626       }
1627       // chi2
1628       TMatrixD deltar = rv; deltar -= rb;
1629       TMatrixD vVbInvdeltar(vVbInv,TMatrixD::kMult,deltar);
1630       chi2b = deltar(0,0)*vVbInvdeltar(0,0)+
1631               deltar(1,0)*vVbInvdeltar(1,0)+
1632               deltar(2,0)*vVbInvdeltar(2,0);
1633       chi2 += chi2b;
1634     }
1635
1636     // loop on tracks
1637     for(k=0; k<nTrksSel; k++) {
1638       //
1639       // get track from track array
1640       t = (AliExternalTrackParam*)fTrkArraySel.At(k);
1641       if (useWeights && t->TestBit(kBitUsed)) continue;
1642       //
1643       int tBC = int(t->GetUniqueID()) - kTOFBCShift;    // BC assigned to this track
1644       if (fSelectOnTOFBunchCrossing) {
1645         if (!fKeepAlsoUnflaggedTOFBunchCrossing) continue;   // don't consider tracks with undefined BC
1646         if (currBC!=AliVTrack::kTOFBCNA && tBC!=AliVTrack::kTOFBCNA && tBC!=currBC) continue;  // track does not match to current BCid
1647       }
1648       alpha = t->GetAlpha();
1649       xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha);
1650       // to vtxSeed (from finder)
1651       t->PropagateTo(xlStart,GetFieldkG());
1652
1653       // vector of track global coordinates
1654       TMatrixD ri(3,1);
1655       // covariance matrix of ri
1656       TMatrixD wWi(3,3);
1657
1658       // get space point from track
1659       if(!TrackToPoint(t,ri,wWi,uUi3by3)) continue;
1660
1661       // track chi2
1662       TMatrixD deltar = rv; deltar -= ri;
1663       TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
1664       chi2i = deltar(0,0)*wWideltar(0,0)+
1665               deltar(1,0)*wWideltar(1,0)+
1666               deltar(2,0)*wWideltar(2,0);
1667       //
1668       if (useWeights) {
1669         //double sg = TMath::Sqrt(fMVSigma2);
1670         //double chi2iw = (deltar(0,0)*wWideltar(0,0)+deltar(1,0)*wWideltar(1,0))/sg + deltar(2,0)*wWideltar(2,0)/fMVSigma2;
1671         //double wgh = (1-chi2iw/fMVTukey2);
1672         double chi2iw = chi2i;
1673         double wgh = (1-chi2iw/fMVTukey2/fMVSigma2);
1674
1675         if (wgh<kTiny) wgh = 0;
1676         else if (useWeights==2) wgh = 1.; // use as binary weight
1677         if (step==1) ((AliESDtrack*)t)->SetBit(kBitUsed, wgh>0);
1678         if (wgh<kTiny) continue; // discard the track
1679         wWi *= wgh;  // RS: use weight?
1680         fMVWSum += wgh;
1681         fMVWE2  += wgh*chi2iw;
1682       }
1683       // add to total chi2
1684       if (fSelectOnTOFBunchCrossing && tBC!=AliVTrack::kTOFBCNA && currBC<0) currBC = tBC;
1685       //
1686       chi2 += chi2i;
1687       TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
1688       sumWiri += wWiri;
1689       sumWi   += wWi;
1690
1691       nTrksUsed++;
1692     } // end loop on tracks
1693
1694     if(nTrksUsed < fMinTracks) {
1695       failed=1;
1696       continue;
1697     }
1698
1699     Double_t determinant = sumWi.Determinant();
1700     if(determinant < fMinDetFitter)  {
1701       AliDebug(1,Form("det(V) = %f (<%f)\n",determinant,fMinDetFitter));
1702       failed=1;
1703       continue;
1704     }
1705
1706     if(step==0) {
1707       // inverted of weights matrix
1708       TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
1709       vV = invsumWi;
1710       // position of primary vertex
1711       rv.Mult(vV,sumWiri);
1712     }
1713   } // end loop on the 2 steps
1714
1715   if(failed) {
1716     fVert.SetNContributors(-1);
1717     if (chiCalc) {
1718       TooFewTracks();
1719       if (fCurrentVertex) fVert = *fCurrentVertex;  // RS
1720     }
1721     return;
1722   }
1723   //
1724   Double_t position[3];
1725   position[0] = rv(0,0);
1726   position[1] = rv(1,0);
1727   position[2] = rv(2,0);
1728   Double_t covmatrix[6];
1729   covmatrix[0] = vV(0,0);
1730   covmatrix[1] = vV(0,1);
1731   covmatrix[2] = vV(1,1);
1732   covmatrix[3] = vV(0,2);
1733   covmatrix[4] = vV(1,2);
1734   covmatrix[5] = vV(2,2);
1735
1736   // for correct chi2/ndf, count constraint as additional "track"
1737   if(fConstraint) nTrksUsed++;
1738   //
1739   if (vfit && !chiCalc) { // RS: special mode for multi-vertex finder
1740     fVert.SetXYZ(position);
1741     fVert.SetCovarianceMatrix(covmatrix);
1742     fVert.SetNContributors(nTrksUsed);
1743     return;
1744   }
1745   // store data in the vertex object
1746   if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
1747   fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);
1748   fCurrentVertex->SetBC(currBC);
1749   fVert = *fCurrentVertex;  // RS
1750   AliDebug(1," Vertex after fit:");
1751   AliDebug(1,Form("xyz: %f %f %f; nc %d",fCurrentVertex->GetXv(),fCurrentVertex->GetYv(),fCurrentVertex->GetZv(),fCurrentVertex->GetNContributors()));
1752   AliDebug(1,"--- VertexFitter(): finish\n");
1753
1754
1755   return;
1756 }
1757 //----------------------------------------------------------------------------
1758 AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(const TObjArray *trkArray,
1759                                                          UShort_t *id,
1760                                                          Bool_t optUseFitter,
1761                                                          Bool_t optPropagate,
1762                                                          Bool_t optUseDiamondConstraint)
1763 {
1764 //
1765 // Return vertex from tracks (AliExternalTrackParam) in array
1766 //
1767   fCurrentVertex = 0;
1768   // set optUseDiamondConstraint=TRUE only if you are reconstructing the
1769   // primary vertex!
1770   if(optUseDiamondConstraint) {
1771     SetConstraintOn();
1772   } else {
1773     SetConstraintOff();
1774   }
1775
1776   // get tracks and propagate them to initial vertex position
1777   fIdSel = new UShort_t[(Int_t)trkArray->GetEntriesFast()];
1778   Int_t nTrksSel = PrepareTracks(*trkArray,id,0);
1779   if((!optUseDiamondConstraint && nTrksSel<TMath::Max(2,fMinTracks)) ||
1780      (optUseDiamondConstraint && nTrksSel<1)) {
1781     TooFewTracks();
1782     return fCurrentVertex;
1783   }
1784
1785   // vertex finder
1786   if(nTrksSel==1) {
1787     AliDebug(1,"Just one track");
1788     OneTrackVertFinder();
1789   } else {
1790     switch (fAlgo) {
1791       case 1: StrLinVertexFinderMinDist(1); break;
1792       case 2: StrLinVertexFinderMinDist(0); break;
1793       case 3: HelixVertexFinder();          break;
1794       case 4: VertexFinder(1);              break;
1795       case 5: VertexFinder(0);              break;
1796       default: printf("Wrong algorithm\n"); break;
1797     }
1798   }
1799   AliDebug(1," Vertex finding completed\n");
1800   Double_t vdispersion=fVert.GetDispersion();
1801
1802   // vertex fitter
1803   if(optUseFitter) {
1804     VertexFitter();
1805   } else {
1806     Double_t position[3]={fVert.GetXv(),fVert.GetYv(),fVert.GetZv()};
1807     Double_t covmatrix[6];
1808     fVert.GetCovMatrix(covmatrix);
1809     Double_t chi2=99999.;
1810     Int_t    nTrksUsed=fVert.GetNContributors();
1811     fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);
1812   }
1813   fCurrentVertex->SetDispersion(vdispersion);
1814
1815
1816   // set indices of used tracks and propagate track to found vertex
1817   UShort_t *indices = 0;
1818   Double_t d0z0[2],covd0z0[3];
1819   AliExternalTrackParam *t = 0;
1820   if(fCurrentVertex->GetNContributors()>0) {
1821     indices = new UShort_t[fTrkArraySel.GetEntriesFast()];
1822     for(Int_t jj=0; jj<(Int_t)fTrkArraySel.GetEntriesFast(); jj++) {
1823       indices[jj] = fIdSel[jj];
1824       t = (AliExternalTrackParam*)fTrkArraySel.At(jj);
1825       if(optPropagate && optUseFitter) {
1826         if(TMath::Sqrt(fCurrentVertex->GetXv()*fCurrentVertex->GetXv()+fCurrentVertex->GetYv()*fCurrentVertex->GetYv())<3.) {
1827           t->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
1828           AliDebug(1,Form("Track %d propagated to found vertex",jj));
1829         } else {
1830           AliWarning("Found vertex outside beam pipe!");
1831         }
1832       }
1833     }
1834     fCurrentVertex->SetIndices(fCurrentVertex->GetNContributors(),indices);
1835   }
1836
1837   // clean up
1838   if (indices) {delete [] indices; indices=NULL;}
1839   delete [] fIdSel; fIdSel=NULL;
1840   fTrkArraySel.Delete();
1841
1842   return fCurrentVertex;
1843 }
1844
1845 //----------------------------------------------------------------------------
1846 AliESDVertex* AliVertexerTracks::VertexForSelectedESDTracks(TObjArray *trkArray,
1847                                                              Bool_t optUseFitter,
1848                                                             Bool_t optPropagate,
1849                                                             Bool_t optUseDiamondConstraint)
1850
1851 {
1852 //
1853 // Return vertex from array of ESD tracks
1854 //
1855
1856   Int_t nTrks = (Int_t)trkArray->GetEntriesFast();
1857   UShort_t *id = new UShort_t[nTrks];
1858
1859   AliESDtrack *esdt = 0;
1860   for(Int_t i=0; i<nTrks; i++){
1861     esdt = (AliESDtrack*)trkArray->At(i);
1862     id[i] = (UShort_t)esdt->GetID();
1863   }
1864
1865   VertexForSelectedTracks(trkArray,id,optUseFitter,optPropagate,optUseDiamondConstraint);
1866
1867   delete [] id; id=NULL;
1868
1869   return fCurrentVertex;
1870 }
1871
1872 //______________________________________________________
1873 Bool_t AliVertexerTracks::FindNextVertexMV()
1874 {
1875   // try to find a new vertex
1876   fMVSigma2 = fMVSig2Ini;
1877   double prevSig2 = fMVSigma2;
1878   double minDst2 = fMVMinDst*fMVMinDst;
1879   const double kSigLimit = 1.0;
1880   const double kSigLimitE = kSigLimit+1e-6;
1881   const double kPushFactor = 0.5;
1882   const int kMaxIter = 20;
1883   double push = kPushFactor;
1884   //
1885   int iter = 0;
1886   double posP[3]={0,0,0},pos[3]={0,0,0};
1887   fVert.GetXYZ(posP);
1888   //
1889
1890   do {
1891     fVert.SetBC(AliVTrack::kTOFBCNA);
1892     VertexFitter(kTRUE,kFALSE,1);
1893     if (fVert.GetNContributors()<fMinTracks) {
1894       AliDebug(3,Form("Failed in iteration %d: No Contributirs",iter));
1895       break;
1896     } // failed
1897     if (fMVWSum>0) fMVSigma2 = TMath::Max(fMVWE2/fMVWSum,kSigLimit);
1898     else {
1899       AliDebug(3,Form("Failed %d, no weithgs",iter));
1900       iter = kMaxIter+1;
1901       break;
1902     } // failed
1903     //
1904     double sigRed = (prevSig2-fMVSigma2)/prevSig2; // reduction of sigma2
1905     //
1906     fVert.GetXYZ(pos);
1907     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]);
1908     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));
1909     if ( (++iter<kMaxIter) && (sigRed<0 || sigRed<fMVMinSig2Red) && fMVSigma2>fMVMaxSigma2) {
1910       fMVSigma2 *= push; // stuck, push little bit
1911       push *= kPushFactor;
1912       if (fMVSigma2<1.) fMVSigma2 = 1.;
1913       AliDebug(3,Form("Pushed sigma2 to %f",fMVSigma2));
1914     }
1915     else if (dst2<minDst2 && ((sigRed<0 || sigRed<fMVMinSig2Red) || fMVSigma2<kSigLimitE)) break;
1916     //
1917     fVert.GetXYZ(posP); // fetch previous vertex position
1918     prevSig2 = fMVSigma2;
1919   } while(iter<kMaxIter);
1920   //
1921   if (fVert.GetNContributors()<0 || iter>kMaxIter || fMVSigma2>fMVMaxSigma2) {
1922     return kFALSE;
1923   }
1924   else {
1925     VertexFitter(kFALSE,kTRUE,fMVFinalWBinary ? 2:1); // final chi2 calculation
1926     int nv = fMVVertices->GetEntries();
1927     // create indices
1928     int ntrk = fTrkArraySel.GetEntries();
1929     int nindices = fCurrentVertex->GetNContributors() - (fConstraint ? 1:0);
1930     UShort_t *indices = 0;
1931     if (nindices>0) indices = new UShort_t[nindices];
1932     int nadded = 0;
1933     for (int itr=0;itr<ntrk;itr++) {
1934       AliExternalTrackParam* t = (AliExternalTrackParam*)fTrkArraySel[itr];
1935       if (t->TestBit(kBitAccounted) || !t->TestBit(kBitUsed)) continue;   // already belongs to some vertex
1936       t->SetBit(kBitAccounted);
1937       indices[nadded++] = fIdSel[itr];
1938     }
1939     if (nadded!=nindices) printf("Mismatch : NInd: %d Nadd: %d\n",nindices,nadded);
1940     fCurrentVertex->SetIndices(nindices,indices);
1941     // set vertex title
1942     TString title="VertexerTracksMVNoConstraint";
1943     if(fConstraint) title="VertexerTracksMVWithConstraint";
1944     fCurrentVertex->SetTitle(title.Data());
1945     fMVVertices->AddLast(fCurrentVertex);
1946     AliDebug(3,Form("Added new vertex #%d NCont:%d XYZ: %f %f %f",nindices,nv,fCurrentVertex->GetX(),fCurrentVertex->GetY(),fCurrentVertex->GetZ()));
1947     if (indices) delete[] indices;
1948     fCurrentVertex = 0; // already attached to fMVVertices
1949     return kTRUE;
1950   }
1951 }
1952
1953 //______________________________________________________
1954 void AliVertexerTracks::FindVerticesMV()
1955 {
1956   // find and fit multiple vertices
1957   //
1958   double step = fMVScanStep>1 ?  fMVScanStep : 1.;
1959   double zmx = 3*TMath::Sqrt(fNominalCov[5]);
1960   double zmn = -zmx;
1961   int nz = TMath::Nint((zmx-zmn)/step); if (nz<1) nz=1;
1962   double dz = (zmx-zmn)/nz;
1963   int izStart=0;
1964   AliDebug(2,Form("%d seeds between %f and %f",nz,zmn+dz/2,zmx+dz/2));
1965   //
1966   if (!fMVVertices) fMVVertices = new TObjArray(10);
1967   fMVVertices->Clear();
1968   //
1969   int ntrLeft = (Int_t)fTrkArraySel.GetEntries();
1970   //
1971   double sig2Scan = fMVSig2Ini;
1972   Bool_t runMore = kTRUE;
1973   int cntWide = 0;
1974   while (runMore) {
1975     fMVSig2Ini = sig2Scan*1e3;  // try wide search
1976     Bool_t found = kFALSE;
1977     cntWide++;
1978     fVert.SetNContributors(-1);
1979     fVert.SetXYZ(fNominalPos);
1980     AliDebug(3,Form("Wide search #%d Z= %f Sigma2=%f",cntWide,fNominalPos[2],fMVSig2Ini));
1981     if (FindNextVertexMV()) { // are there tracks left to consider?
1982       AliESDVertex* vtLast = (AliESDVertex*)fMVVertices->Last();
1983       if (vtLast && vtLast->GetNContributors()>0) ntrLeft -= vtLast->GetNContributors()-(fConstraint ? 1:0);
1984       if (ntrLeft<1) runMore = kFALSE;
1985       found = kTRUE;
1986       continue;
1987     }
1988     // if nothing is found, do narrow sig2ini scan
1989     fMVSig2Ini = sig2Scan;
1990     for (;izStart<nz;izStart++) {
1991       double zSeed = zmn+dz*(izStart+0.5);
1992       AliDebug(3,Form("Seed %d: Z= %f Sigma2=%f",izStart,zSeed,fMVSig2Ini));
1993       fVert.SetNContributors(-1);
1994       fVert.SetXYZ(fNominalPos);
1995       fVert.SetZv(zSeed);
1996       if (FindNextVertexMV()) { // are there tracks left to consider?
1997         AliESDVertex* vtLast = (AliESDVertex*)fMVVertices->Last();
1998         if (vtLast && vtLast->GetNContributors()>0) ntrLeft -= vtLast->GetNContributors()-(fConstraint ? 1:0);
1999         if (ntrLeft<1) runMore = kFALSE;
2000         found = kTRUE;
2001         break;
2002       }
2003     }
2004     runMore = found; // if nothing was found, no need for new iteration
2005   }
2006   fMVSig2Ini = sig2Scan;
2007   int nvFound = fMVVertices->GetEntriesFast();
2008   AliDebug(2,Form("Number of found vertices: %d",nvFound));
2009   if (nvFound<1) TooFewTracks();
2010   if (AliLog::GetGlobalDebugLevel()>0) fMVVertices->Print();
2011   //
2012 }
2013
2014 //______________________________________________________
2015 void AliVertexerTracks::AnalyzePileUp(AliESDEvent* esdEv)
2016 {
2017   // if multiple vertices are found, try to find the primary one and attach it as fCurrentVertex
2018   // then attach pile-up vertices directly to esdEv
2019   //
2020   int nFND = (fMVVertices && fMVVertices->GetEntriesFast()) ? fMVVertices->GetEntriesFast() : 0;
2021   if (nFND<1) { if (!fCurrentVertex) TooFewTracks(); return;} // no multiple vertices
2022   //
2023   int indCont[nFND];
2024   int nIndx[nFND];
2025   for (int iv=0;iv<nFND;iv++) {
2026     AliESDVertex* fnd = (AliESDVertex*)fMVVertices->At(iv);
2027     indCont[iv] = iv;
2028     nIndx[iv]   = fnd->GetNIndices();
2029   }
2030   TMath::Sort(nFND, nIndx, indCont, kTRUE); // sort in decreasing order of Nindices
2031   double dists[nFND];
2032   int    distOrd[nFND],indx[nFND];
2033   for (int iv=0;iv<nFND;iv++) {
2034     AliESDVertex* fndI = (AliESDVertex*)fMVVertices->At(indCont[iv]);
2035     if (fndI->GetStatus()<1) continue; // discarded
2036     int ncomp = 0;
2037     for (int jv=iv+1;jv<nFND;jv++) {
2038       AliESDVertex* fndJ = (AliESDVertex*)fMVVertices->At(indCont[jv]);
2039       if (fndJ->GetStatus()<1) continue;
2040       dists[ncomp] = fndI->GetWDist(fndJ)*fndJ->GetNIndices();
2041       distOrd[ncomp] = indCont[jv];
2042       indx[ncomp]  = ncomp;
2043       ncomp++;
2044     }
2045     if (ncomp<1) continue;
2046     TMath::Sort(ncomp, dists, indx, kFALSE); // sort in increasing distance order
2047     for (int jv0=0;jv0<ncomp;jv0++) {
2048       int jv = distOrd[indx[jv0]];
2049       AliESDVertex* fndJ = (AliESDVertex*)fMVVertices->At(jv);
2050       if (dists[indx[jv0]]<fMVMaxWghNtr) { // candidate for split vertex
2051         //before eliminating the small close vertex, check if we should transfere its BC to largest one
2052         if (fndJ->GetBC()!=AliVTrack::kTOFBCNA && fndI->GetBC()==AliVTrack::kTOFBCNA) fndI->SetBC(fndJ->GetBC());
2053         //
2054         // leave the vertex only if both vertices have definit but different BC's, then this is not splitting.
2055         if ( (fndJ->GetBC()==fndI->GetBC()) || (fndJ->GetBC()==AliVTrack::kTOFBCNA)) fndJ->SetNContributors(-fndJ->GetNContributors());
2056       }
2057     }
2058   }
2059   //
2060   // select as a primary the largest vertex with BC=0, or the largest with BC non-ID
2061   int primBC0=-1,primNoBC=-1;
2062   for (int iv0=0;iv0<nFND;iv0++) {
2063     int iv = indCont[iv0];
2064     AliESDVertex* fndI = (AliESDVertex*)fMVVertices->At(iv);
2065     if (!fndI) continue;
2066     if (fndI->GetStatus()<1) {fMVVertices->RemoveAt(iv); delete fndI; continue;}
2067     if (primBC0<0  && fndI->GetBC()==0) primBC0 = iv;
2068     if (primNoBC<0 && fndI->GetBC()==AliVTrack::kTOFBCNA)  primNoBC = iv;
2069   }
2070   //
2071   if (primBC0>=0) fCurrentVertex = (AliESDVertex*)fMVVertices->At(primBC0);
2072   if (!fCurrentVertex && primNoBC>=0) fCurrentVertex = (AliESDVertex*)fMVVertices->At(primNoBC);
2073   if (fCurrentVertex) fMVVertices->Remove(fCurrentVertex);
2074   else {  // all vertices have BC>0, no primary vertex
2075     fCurrentVertex = new AliESDVertex();
2076     fCurrentVertex->SetNContributors(-3);
2077     fCurrentVertex->SetBC(AliVTrack::kTOFBCNA);
2078   }
2079   fCurrentVertex->SetID(-1);
2080   //
2081   // add pileup vertices
2082   int nadd = 0;
2083   for (int iv0=0;iv0<nFND;iv0++) {
2084     int iv = indCont[iv0];
2085     AliESDVertex* fndI = (AliESDVertex*)fMVVertices->At(iv);
2086     if (!fndI) continue;
2087     fndI->SetID(++nadd);
2088     esdEv->AddPileupVertexTracks(fndI);
2089   }
2090   //
2091   fMVVertices->Delete();
2092   //
2093 }
2094
2095 //______________________________________________________
2096 void AliVertexerTracks::FindAllVertices(Int_t nTrksOrig,
2097                                         const TObjArray *trkArrayOrig,
2098                                         Double_t* zTr,
2099                                         Double_t* err2zTr,
2100                                         UShort_t* idOrig){
2101
2102   // clusterize tracks using z coordinates of intersection with beam axis
2103   // and compute all vertices
2104
2105   UShort_t* posOrig=new UShort_t[nTrksOrig];
2106   for(Int_t iTr=0; iTr<nTrksOrig; iTr++) posOrig[iTr]=iTr;
2107
2108
2109   // sort points along Z
2110   AliDebug(1,Form("Sort points along Z, used tracks %d",nTrksOrig));
2111   for(Int_t iTr1=0; iTr1<nTrksOrig; iTr1++){
2112     for(Int_t iTr2=iTr1+1; iTr2<nTrksOrig; iTr2++){
2113       if(zTr[iTr1]>zTr[iTr2]){
2114         Double_t tmpz=zTr[iTr2];
2115         Double_t tmperr=err2zTr[iTr2];
2116         UShort_t tmppos=posOrig[iTr2];
2117         UShort_t tmpid=idOrig[iTr2];
2118         zTr[iTr2]=zTr[iTr1];
2119         err2zTr[iTr2]=err2zTr[iTr1];
2120         posOrig[iTr2]=posOrig[iTr1];
2121         idOrig[iTr2]=idOrig[iTr1];
2122         zTr[iTr1]=tmpz;
2123         err2zTr[iTr1]=tmperr;
2124         idOrig[iTr1]=tmpid;
2125         posOrig[iTr1]=tmppos;
2126       }
2127     }
2128   }
2129
2130   // clusterize
2131   Int_t nClusters=0;
2132   Int_t* firstTr=new Int_t[nTrksOrig];
2133   Int_t* lastTr=new Int_t[nTrksOrig];
2134
2135   firstTr[0]=0;
2136   for(Int_t iTr=0; iTr<nTrksOrig-1; iTr++){
2137     Double_t distz=zTr[iTr+1]-zTr[iTr];
2138     Double_t errdistz=TMath::Sqrt(err2zTr[iTr+1]+err2zTr[iTr]);
2139     if(errdistz<=0.000001) errdistz=0.000001;
2140     if(distz>fDeltaZCutForCluster || (distz/errdistz)>fnSigmaZCutForCluster){
2141       lastTr[nClusters]=iTr;
2142       firstTr[nClusters+1]=iTr+1;
2143       nClusters++;
2144     }
2145   }
2146   lastTr[nClusters]=nTrksOrig-1;
2147
2148   // Compute vertices
2149   AliDebug(1,Form("Number of found clusters %d",nClusters+1));
2150   Int_t nFoundVertices=0;
2151
2152   if (!fMVVertices) fMVVertices = new TObjArray(nClusters+1);
2153
2154   fMVVertices->Clear();
2155   TObjArray cluTrackArr(nTrksOrig);
2156   UShort_t *idTrkClu=new UShort_t[nTrksOrig];
2157   //  Int_t maxContr=0;
2158   //  Int_t maxPos=-1;
2159
2160   for(Int_t iClu=0; iClu<=nClusters; iClu++){
2161     Int_t nCluTracks=lastTr[iClu]-firstTr[iClu]+1;
2162     cluTrackArr.Clear();
2163     AliDebug(1,Form(" Vertex #%d tracks %d first tr %d  last track %d",iClu,nCluTracks,firstTr[iClu],lastTr[iClu]));
2164     Int_t nSelTr=0;
2165     for(Int_t iTr=firstTr[iClu]; iTr<=lastTr[iClu]; iTr++){
2166       AliExternalTrackParam* t=(AliExternalTrackParam*)trkArrayOrig->At(posOrig[iTr]);
2167       if(TMath::Abs(t->GetZ()-zTr[iTr])>0.00001){
2168         AliError(Form("Clu %d Track %d zTrack=%f  zVec=%f\n",iClu,iTr,t->GetZ(),zTr[iTr]));
2169       }
2170       cluTrackArr.AddAt(t,nSelTr);
2171       idTrkClu[nSelTr]=idOrig[iTr];
2172       AliDebug(1,Form("   Add track %d: id %d, z=%f",iTr,idOrig[iTr],zTr[iTr]));
2173       nSelTr++;
2174     }
2175     AliESDVertex* vert=FindPrimaryVertex(&cluTrackArr,idTrkClu);
2176     AliDebug(1,Form("Found vertex in z=%f with %d contributors",vert->GetZv(),
2177                  vert->GetNContributors()));
2178
2179     fCurrentVertex=0;
2180     if(vert->GetNContributors()>0){
2181       nFoundVertices++;
2182       fMVVertices->AddLast(vert);
2183     }
2184     //    if(vert->GetNContributors()>maxContr){
2185     //      maxContr=vert->GetNContributors();
2186     //      maxPos=nFoundVertices-1;
2187     //    }
2188   }
2189
2190   AliDebug(1,Form("Number of found vertices %d (%d)",nFoundVertices,fMVVertices->GetEntriesFast()));
2191   // if(maxPos>=0 && maxContr>0){
2192   //   AliESDVertex* vtxMax=(AliESDVertex*)fMVVertices->At(maxPos);
2193   //   if(fCurrentVertex) delete fCurrentVertex;
2194   //   fCurrentVertex=new AliESDVertex(*vtxMax);
2195   // }
2196
2197   delete [] firstTr;
2198   delete [] lastTr;
2199   delete [] idTrkClu;
2200   delete [] posOrig;
2201
2202   return;
2203
2204 }