STEER/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 <TDirectory.h>
  31 #include <TFile.h>
  32 #include <TTree.h>
  33 #include <TMatrixD.h>
  34 //---- AliRoot headers -----
  35 #include "AliLog.h"
  36 #include "AliStrLine.h"
  37 #include "AliExternalTrackParam.h"
  38 #include "AliNeutralTrackParam.h"
  39 #include "AliVEvent.h"
  40 #include "AliVTrack.h"
  41 #include "AliESDEvent.h"
  42 #include "AliESDtrack.h"
  43 #include "AliVertexerTracks.h"
  44
  45 ClassImp(AliVertexerTracks)
  46
  47
  48 //----------------------------------------------------------------------------
  49 AliVertexerTracks::AliVertexerTracks():
  50 TObject(),
  51 fVert(),
  52 fCurrentVertex(0),
  53 fMode(0),
  54 fFieldkG(-999.),
  55 fTrkArraySel(),
  56 fIdSel(0),
  57 fTrksToSkip(0),
  58 fNTrksToSkip(0),
  59 fConstraint(kFALSE),
  60 fOnlyFitter(kFALSE),
  61 fMinTracks(1),
  62 fMinClusters(5),
  63 fDCAcut(0.1),
  64 fDCAcutIter0(0.1),
  65 fNSigma(3.),
  66 fMaxd0z0(0.5),
  67 fMinDetFitter(100.),
  68 fMaxTgl(1000.),
  69 fITSrefit(kTRUE),
  70 fFiducialR(3.),
  71 fFiducialZ(30.),
  72 fnSigmaForUi00(1.5),
  73 fAlgo(1),
  74 fAlgoIter0(4)
  75 {
  76 //
  77 // Default constructor
  78 //
  79   SetVtxStart();
  80   SetVtxStartSigma();
  81 }
  82 //----------------------------------------------------------------------------
  83 AliVertexerTracks::AliVertexerTracks(Double_t fieldkG):
  84 TObject(),
  85 fVert(),
  86 fCurrentVertex(0),
  87 fMode(0),
  88 fFieldkG(fieldkG),
  89 fTrkArraySel(),
  90 fIdSel(0),
  91 fTrksToSkip(0),
  92 fNTrksToSkip(0),
  93 fConstraint(kFALSE),
  94 fOnlyFitter(kFALSE),
  95 fMinTracks(1),
  96 fMinClusters(5),
  97 fDCAcut(0.1),
  98 fDCAcutIter0(0.1),
  99 fNSigma(3.),
 100 fMaxd0z0(0.5),
 101 fMinDetFitter(100.),
 102 fMaxTgl(1000.),
 103 fITSrefit(kTRUE),
 104 fFiducialR(3.),
 105 fFiducialZ(30.),
 106 fnSigmaForUi00(1.5),
 107 fAlgo(1),
 108 fAlgoIter0(4)
 109 {
 110 //
 111 // Standard constructor
 112 //
 113   SetVtxStart();
 114   SetVtxStartSigma();
 115 }
 116 //-----------------------------------------------------------------------------
 117 AliVertexerTracks::~AliVertexerTracks()
 118 {
 119   // Default Destructor
 120   // The objects pointed by the following pointer are not owned
 121   // by this class and are not deleted
 122   fCurrentVertex = 0;
 123   if(fTrksToSkip) { delete fTrksToSkip; fTrksToSkip=NULL; }
 124   if(fIdSel) { delete fIdSel; fIdSel=NULL; }
 125 }
 126 //----------------------------------------------------------------------------
 127 AliESDVertex* AliVertexerTracks::FindPrimaryVertex(AliVEvent *vEvent)
 128 {
 129 //
 130 // Primary vertex for current ESD or AOD event
 131 // (Two iterations:
 132 //  1st with 5*fNSigma*sigma cut w.r.t. to initial vertex
 133 //      + cut on sqrt(d0d0+z0z0) if fConstraint=kFALSE
 134 //  2nd with fNSigma*sigma cut w.r.t. to vertex found in 1st iteration)
 135 //
 136   fCurrentVertex = 0;
 137
 138   TString evtype = vEvent->IsA()->GetName();
 139   Bool_t inputAOD = ((evtype=="AliAODEvent") ? kTRUE : kFALSE);
 140
 141   if(inputAOD && fMode==1) {
 142     printf("Error : AliVertexerTracks: no TPC-only vertex from AOD\n");
 143     TooFewTracks();
 144     return fCurrentVertex;
 145   }
 146
 147   // accept 1-track case only if constraint is available
 148   if(!fConstraint && fMinTracks==1) fMinTracks=2;
 149
 150   // read tracks from AlivEvent
 151   Int_t nTrks = (Int_t)vEvent->GetNumberOfTracks();
 152   if(nTrks<fMinTracks) {
 153     TooFewTracks();
 154     return fCurrentVertex;
 155   }
 156   //
 157
 158   TDirectory * olddir = gDirectory;
 159   TFile *f = 0;
 160   if(nTrks>500) f = new TFile("VertexerTracks.root","recreate");
 161   TObjArray trkArrayOrig(nTrks);
 162   UShort_t *idOrig = new UShort_t[nTrks];
 163
 164   Int_t nTrksOrig=0;
 165   AliExternalTrackParam *t=0;
 166   // loop on tracks
 167   for(Int_t i=0; i<nTrks; i++) {
 168     AliVTrack *track = (AliVTrack*)vEvent->GetTrack(i);
 169     // check tracks to skip
 170     Bool_t skipThis = kFALSE;
 171     for(Int_t j=0; j<fNTrksToSkip; j++) {
 172       if(track->GetID()==fTrksToSkip[j]) {
 173         AliDebug(1,Form("skipping track: %d",i));
 174         skipThis = kTRUE;
 175       }
 176     }
 177     if(skipThis) continue;
 178
 179     // kITSrefit
 180     if(fMode==0 && fITSrefit && !(track->GetStatus()&AliESDtrack::kITSrefit)) continue;
 181
 182     if(!inputAOD) {  // ESD
 183       AliESDtrack* esdt = (AliESDtrack*)track;
 184       if(esdt->GetNcls(fMode) < fMinClusters) continue;
 185       if(fMode==0) {        // ITS mode
 186         Double_t x,p[5],cov[15];
 187         esdt->GetExternalParameters(x,p);
 188         esdt->GetExternalCovariance(cov);
 189         t = new AliExternalTrackParam(x,esdt->GetAlpha(),p,cov);
 190       } else if(fMode==1) { // TPC mode
 191         t = (AliExternalTrackParam*)esdt->GetTPCInnerParam();
 192         if(!t) continue;
 193         Double_t radius = 2.8; //something less than the beam pipe radius
 194         if(!PropagateTrackTo(t,radius)) continue;
 195       }
 196     } else {          // AOD (only ITS mode)
 197       Int_t ncls0=0;
 198       for(Int_t l=0;l<6;l++) if(TESTBIT(track->GetITSClusterMap(),l)) ncls0++;
 199       if(ncls0 < fMinClusters) continue;
 200       t = new AliExternalTrackParam(track);
 201     }
 202     trkArrayOrig.AddLast(t);
 203     idOrig[nTrksOrig]=(UShort_t)track->GetID();
 204     nTrksOrig++;
 205   } // end loop on tracks
 206
 207   // call method that will reconstruct the vertex
 208   FindPrimaryVertex(&trkArrayOrig,idOrig);
 209
 210   if(fMode==0) trkArrayOrig.Delete();
 211   delete idOrig; idOrig=NULL;
 212
 213   if(f) {
 214     f->Close(); delete f; f = NULL;
 215     gSystem->Unlink("VertexerTracks.root");
 216     olddir->cd();
 217   }
 218
 219   // set vertex ID for tracks used in the fit
 220   // (only for ESD)
 221   if(!inputAOD) {
 222     Int_t nIndices = fCurrentVertex->GetNIndices();
 223     UShort_t *indices = fCurrentVertex->GetIndices();
 224     for(Int_t ind=0; ind<nIndices; ind++) {
 225       AliESDtrack *esdt = (AliESDtrack*)vEvent->GetTrack(indices[ind]);
 226       esdt->SetVertexID(-1);
 227     }
 228   }
 229
 230   return fCurrentVertex;
 231 }
 232 //----------------------------------------------------------------------------
 233 AliESDVertex* AliVertexerTracks::FindPrimaryVertex(TObjArray *trkArrayOrig,
 234                                                    UShort_t *idOrig)
 235 {
 236 //
 237 // Primary vertex using the AliExternalTrackParam's in the TObjArray.
 238 // idOrig must contain the track IDs from AliESDtrack::GetID()
 239 // (Two iterations:
 240 //  1st with 5*fNSigma*sigma cut w.r.t. to initial vertex
 241 //      + cut on sqrt(d0d0+z0z0) if fConstraint=kFALSE
 242 //  2nd with fNSigma*sigma cut w.r.t. to vertex found in 1st iteration)
 243 //
 244   fCurrentVertex = 0;
 245
 246   // accept 1-track case only if constraint is available
 247   if(!fConstraint && fMinTracks==1) fMinTracks=2;
 248
 249   // read tracks from array
 250   Int_t nTrksOrig = (Int_t)trkArrayOrig->GetEntriesFast();
 251   AliDebug(1,Form("Initial number of tracks: %d",nTrksOrig));
 252   if(nTrksOrig<fMinTracks) {
 253     AliDebug(1,"TooFewTracks");
 254     TooFewTracks();
 255     return fCurrentVertex;
 256   }
 257
 258   // If fConstraint=kFALSE
 259   // run VertexFinder(1) to get rough estimate of initVertex (x,y)
 260   if(!fConstraint) {
 261     // fill fTrkArraySel, for VertexFinder()
 262     fIdSel = new UShort_t[nTrksOrig];
 263     PrepareTracks(*trkArrayOrig,idOrig,0);
 264     if(fIdSel) { delete fIdSel; fIdSel=NULL; }
 265     Double_t cutsave = fDCAcut;
 266     fDCAcut = fDCAcutIter0;
 267     // vertex finder
 268     switch (fAlgoIter0) {
 269     case 1: StrLinVertexFinderMinDist(1); break;
 270     case 2: StrLinVertexFinderMinDist(0); break;
 271     case 3: HelixVertexFinder();          break;
 272     case 4: VertexFinder(1);              break;
 273     case 5: VertexFinder(0);              break;
 274     default: printf("Wrong algorithm\n"); break;
 275     }
 276     fDCAcut = cutsave;
 277     if(fVert.GetNContributors()>0) {
 278       fVert.GetXYZ(fNominalPos);
 279       fNominalPos[0] = fVert.GetXv();
 280       fNominalPos[1] = fVert.GetYv();
 281       fNominalPos[2] = fVert.GetZv();
 282       AliDebug(1,Form("No mean vertex: VertexFinder gives (%f, %f, %f)",fNominalPos[0],fNominalPos[1],fNominalPos[2]));
 283     } else {
 284       fNominalPos[0] = 0.;
 285       fNominalPos[1] = 0.;
 286       fNominalPos[2] = 0.;
 287       AliDebug(1,"No mean vertex and VertexFinder failed");
 288     }
 289   }
 290
 291   // TWO ITERATIONS:
 292   //
 293   // ITERATION 1
 294   // propagate tracks to fNominalPos vertex
 295   // preselect them:
 296   // if(constraint) reject for |d0|>5*fNSigma*sigma w.r.t. fNominal... vertex
 297   // else  reject for |d0|\oplus|z0| > 5 mm w.r.t. fNominal... vertex
 298   // ITERATION 2
 299   // propagate tracks to best between initVertex and fCurrentVertex
 300   // preselect tracks (reject for |d0|>fNSigma*sigma w.r.t. best
 301   //                   between initVertex and fCurrentVertex)
 302   for(Int_t iter=1; iter<=2; iter++) {
 303     if(fOnlyFitter && iter==1) continue;
 304     if(fIdSel) { delete fIdSel; fIdSel=NULL; }
 305     fIdSel = new UShort_t[nTrksOrig];
 306     Int_t nTrksSel = PrepareTracks(*trkArrayOrig,idOrig,iter);
 307     AliDebug(1,Form("N tracks selected in iteration %d: %d",iter,nTrksSel));
 308     if(nTrksSel < fMinTracks) {
 309       TooFewTracks();
 310       return fCurrentVertex;
 311     }
 312
 313     // vertex finder
 314     if(!fOnlyFitter) {
 315       if(nTrksSel==1) {
 316         AliDebug(1,"Just one track");
 317         OneTrackVertFinder();
 318       } else {
 319         switch (fAlgo) {
 320         case 1: StrLinVertexFinderMinDist(1); break;
 321         case 2: StrLinVertexFinderMinDist(0); break;
 322         case 3: HelixVertexFinder();          break;
 323         case 4: VertexFinder(1);              break;
 324         case 5: VertexFinder(0);              break;
 325         default: printf("Wrong algorithm"); break;
 326         }
 327       }
 328       AliDebug(1," Vertex finding completed");
 329     }
 330
 331     // vertex fitter
 332     VertexFitter();
 333   } // end loop on the two iterations
 334
 335
 336   // set indices of used tracks
 337   UShort_t *indices = 0;
 338   if(fCurrentVertex->GetNContributors()>0) {
 339     Int_t nIndices = (Int_t)fTrkArraySel.GetEntriesFast();
 340     indices = new UShort_t[nIndices];
 341     for(Int_t jj=0; jj<nIndices; jj++)
 342       indices[jj] = fIdSel[jj];
 343     fCurrentVertex->SetIndices(nIndices,indices);
 344   }
 345   if (indices) {delete indices; indices=NULL;}
 346   //
 347
 348   // set vertex title
 349   TString title="VertexerTracksNoConstraint";
 350   if(fConstraint) {
 351     title="VertexerTracksWithConstraint";
 352     if(fOnlyFitter) title.Append("OnlyFitter");
 353   }
 354   fCurrentVertex->SetTitle(title.Data());
 355   //
 356
 357   AliDebug(1,Form("xyz: %f %f %f; nc %d",fCurrentVertex->GetXv(),fCurrentVertex->GetYv(),fCurrentVertex->GetZv(),fCurrentVertex->GetNContributors()));
 358
 359   // clean up
 360   delete fIdSel; fIdSel=NULL;
 361   fTrkArraySel.Delete();
 362   if(fTrksToSkip) { delete fTrksToSkip; fTrksToSkip=NULL; }
 363   //
 364
 365   return fCurrentVertex;
 366 }
 367 //------------------------------------------------------------------------
 368 Double_t AliVertexerTracks::GetDeterminant3X3(Double_t matr[][3])
 369 {
 370   //
 371   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];
 372  return det;
 373 }
 374 //-------------------------------------------------------------------------
 375 void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t (*m)[3],Double_t *d)
 376 {
 377   //
 378   Double_t x12=p0[0]-p1[0];
 379   Double_t y12=p0[1]-p1[1];
 380   Double_t z12=p0[2]-p1[2];
 381   Double_t kk=x12*x12+y12*y12+z12*z12;
 382   m[0][0]=2-2/kk*x12*x12;
 383   m[0][1]=-2/kk*x12*y12;
 384   m[0][2]=-2/kk*x12*z12;
 385   m[1][0]=-2/kk*x12*y12;
 386   m[1][1]=2-2/kk*y12*y12;
 387   m[1][2]=-2/kk*y12*z12;
 388   m[2][0]=-2/kk*x12*z12;
 389   m[2][1]=-2*y12*z12;
 390   m[2][2]=2-2/kk*z12*z12;
 391   d[0]=2*p0[0]-2/kk*p0[0]*x12*x12-2/kk*p0[2]*x12*z12-2/kk*p0[1]*x12*y12;
 392   d[1]=2*p0[1]-2/kk*p0[1]*y12*y12-2/kk*p0[0]*x12*y12-2/kk*p0[2]*z12*y12;
 393   d[2]=2*p0[2]-2/kk*p0[2]*z12*z12-2/kk*p0[0]*x12*z12-2/kk*p0[1]*z12*y12;
 394
 395 }
 396 //--------------------------------------------------------------------------
 397 void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t *sigmasq,Double_t (*m)[3],Double_t *d)
 398 {
 399   //
 400   Double_t x12=p1[0]-p0[0];
 401   Double_t y12=p1[1]-p0[1];
 402   Double_t z12=p1[2]-p0[2];
 403
 404   Double_t den= x12*x12*sigmasq[1]*sigmasq[2]+y12*y12*sigmasq[0]*sigmasq[2]+z12*z12*sigmasq[0]*sigmasq[1];
 405
 406   Double_t kk= 2*(x12*x12/sigmasq[0]+y12*y12/sigmasq[1]+z12*z12/sigmasq[2]);
 407
 408   Double_t cc[3];
 409   cc[0]=-x12/sigmasq[0];
 410   cc[1]=-y12/sigmasq[1];
 411   cc[2]=-z12/sigmasq[2];
 412
 413   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;
 414
 415   Double_t ss= -p0[0]*cc[0]-p0[1]*cc[1]-p0[2]*cc[2];
 416
 417   Double_t aa[3];
 418   aa[0]=x12*sigmasq[1]*sigmasq[2]/den;
 419   aa[1]=y12*sigmasq[0]*sigmasq[2]/den;
 420   aa[2]=z12*sigmasq[0]*sigmasq[1]/den;
 421
 422   m[0][0]=aa[0]*(aa[0]*kk+2*cc[0])+2*cc[0]*aa[0]+2/sigmasq[0];
 423   m[0][1]=aa[1]*(aa[0]*kk+2*cc[0])+2*cc[1]*aa[0];
 424   m[0][2]=aa[2]*(aa[0]*kk+2*cc[0])+2*cc[2]*aa[0];
 425
 426   m[1][0]=aa[0]*(aa[1]*kk+2*cc[1])+2*cc[0]*aa[1];
 427   m[1][1]=aa[1]*(aa[1]*kk+2*cc[1])+2*cc[1]*aa[1]+2/sigmasq[1];
 428   m[1][2]=aa[2]*(aa[1]*kk+2*cc[1])+2*cc[2]*aa[1];
 429
 430   m[2][0]=aa[0]*(aa[2]*kk+2*cc[2])+2*cc[0]*aa[2];
 431   m[2][1]=aa[1]*(aa[2]*kk+2*cc[2])+2*cc[1]*aa[2];
 432   m[2][2]=aa[2]*(aa[2]*kk+2*cc[2])+2*cc[2]*aa[2]+2/sigmasq[2];
 433
 434   d[0]=-ww*(aa[0]*kk+2*cc[0])-2*ss*aa[0]+2*p0[0]/sigmasq[0];
 435   d[1]=-ww*(aa[1]*kk+2*cc[1])-2*ss*aa[1]+2*p0[1]/sigmasq[1];
 436   d[2]=-ww*(aa[2]*kk+2*cc[2])-2*ss*aa[2]+2*p0[2]/sigmasq[2];
 437
 438   }
 439 //--------------------------------------------------------------------------
 440 Double_t AliVertexerTracks::GetStrLinMinDist(Double_t *p0,Double_t *p1,Double_t *x0)
 441 {
 442   //
 443   Double_t x12=p0[0]-p1[0];
 444   Double_t y12=p0[1]-p1[1];
 445   Double_t z12=p0[2]-p1[2];
 446   Double_t x10=p0[0]-x0[0];
 447   Double_t y10=p0[1]-x0[1];
 448   Double_t z10=p0[2]-x0[2];
 449   //  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);
 450   return ((y10*z12-z10*y12)*(y10*z12-z10*y12)+
 451           (z10*x12-x10*z12)*(z10*x12-x10*z12)+
 452           (x10*y12-y10*x12)*(x10*y12-y10*x12))
 453     /(x12*x12+y12*y12+z12*z12);
 454 }
 455 //---------------------------------------------------------------------------
 456 void AliVertexerTracks::OneTrackVertFinder()
 457 {
 458   // find vertex for events with 1 track, using DCA to nominal beam axis
 459   AliDebug(1,Form("Number of prepared tracks =%d - Call OneTrackVertFinder",fTrkArraySel.GetEntries()));
 460   AliExternalTrackParam *track1;
 461   track1 = (AliExternalTrackParam*)fTrkArraySel.At(0);
 462   Double_t alpha=track1->GetAlpha();
 463   Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
 464   Double_t pos[3],dir[3];
 465   track1->GetXYZAt(mindist,GetFieldkG(),pos);
 466   track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
 467   AliStrLine *line1 = new AliStrLine(pos,dir);
 468   Double_t p1[3]={fNominalPos[0],fNominalPos[1],0.};
 469   Double_t p2[3]={fNominalPos[0],fNominalPos[1],10.};
 470   AliStrLine *zeta=new AliStrLine(p1,p2,kTRUE);
 471   Double_t crosspoint[3]={0.,0.,0.};
 472   Double_t sigma=999.;
 473   Int_t nContrib=-1;
 474   Int_t retcode = zeta->Cross(line1,crosspoint);
 475   if(retcode>=0){
 476     sigma=line1->GetDistFromPoint(crosspoint);
 477     nContrib=1;
 478   }
 479   delete zeta;
 480   delete line1;
 481   fVert.SetXYZ(crosspoint);
 482   fVert.SetDispersion(sigma);
 483   fVert.SetNContributors(nContrib);
 484 }
 485 //---------------------------------------------------------------------------
 486 void AliVertexerTracks::HelixVertexFinder()
 487 {
 488   // Get estimate of vertex position in (x,y) from tracks DCA
 489
 490
 491   Double_t initPos[3];
 492   initPos[2] = 0.;
 493   for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
 494
 495   Int_t nacc = (Int_t)fTrkArraySel.GetEntriesFast();
 496
 497   Double_t aver[3]={0.,0.,0.};
 498   Double_t averquad[3]={0.,0.,0.};
 499   Double_t sigmaquad[3]={0.,0.,0.};
 500   Double_t sigma=0;
 501   Int_t ncombi = 0;
 502   AliExternalTrackParam *track1;
 503   AliExternalTrackParam *track2;
 504   Double_t distCA;
 505   Double_t x;
 506   Double_t alpha, cs, sn;
 507   Double_t crosspoint[3];
 508   for(Int_t i=0; i<nacc; i++){
 509     track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
 510
 511
 512     for(Int_t j=i+1; j<nacc; j++){
 513       track2 = (AliExternalTrackParam*)fTrkArraySel.At(j);
 514
 515       distCA=track2->PropagateToDCA(track1,GetFieldkG());
 516       if(fDCAcut<=0 ||(fDCAcut>0&&distCA<fDCAcut)){
 517         x=track1->GetX();
 518         alpha=track1->GetAlpha();
 519         cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
 520         Double_t x1=x*cs - track1->GetY()*sn;
 521         Double_t y1=x*sn + track1->GetY()*cs;
 522         Double_t z1=track1->GetZ();
 523
 524         Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
 525         x=track2->GetX();
 526         alpha=track2->GetAlpha();
 527         cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
 528         Double_t x2=x*cs - track2->GetY()*sn;
 529         Double_t y2=x*sn + track2->GetY()*cs;
 530         Double_t z2=track2->GetZ();
 531         Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
 532         Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
 533         Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
 534         Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
 535         Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
 536         crosspoint[0]=wx1*x1 + wx2*x2;
 537         crosspoint[1]=wy1*y1 + wy2*y2;
 538         crosspoint[2]=wz1*z1 + wz2*z2;
 539
 540         ncombi++;
 541         for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
 542         for(Int_t jj=0;jj<3;jj++)averquad[jj]+=(crosspoint[jj]*crosspoint[jj]);
 543       }
 544     }
 545
 546   }
 547   if(ncombi>0){
 548     for(Int_t jj=0;jj<3;jj++){
 549       initPos[jj] = aver[jj]/ncombi;
 550       averquad[jj]/=ncombi;
 551       sigmaquad[jj]=averquad[jj]-initPos[jj]*initPos[jj];
 552       sigma+=sigmaquad[jj];
 553     }
 554     sigma=TMath::Sqrt(TMath::Abs(sigma));
 555   }
 556   else {
 557     Warning("HelixVertexFinder","Finder did not succed");
 558     sigma=999;
 559   }
 560   fVert.SetXYZ(initPos);
 561   fVert.SetDispersion(sigma);
 562   fVert.SetNContributors(ncombi);
 563 }
 564 //----------------------------------------------------------------------------
 565 Int_t AliVertexerTracks::PrepareTracks(TObjArray &trkArrayOrig,
 566                                        UShort_t *idOrig,
 567                                        Int_t optImpParCut)
 568 {
 569 //
 570 // Propagate tracks to initial vertex position and store them in a TObjArray
 571 //
 572   AliDebug(1," PrepareTracks()");
 573
 574   Int_t nTrksOrig = (Int_t)trkArrayOrig.GetEntriesFast();
 575   Int_t nTrksSel = 0;
 576   Double_t maxd0rphi;
 577   Double_t sigmaCurr[3];
 578   Double_t normdistx,normdisty;
 579   Double_t d0z0[2],covd0z0[3];
 580   Double_t sigmad0;
 581
 582   AliESDVertex *initVertex = new AliESDVertex(fNominalPos,fNominalCov,1,1);
 583
 584   if(!fTrkArraySel.IsEmpty()) fTrkArraySel.Delete();
 585
 586   AliExternalTrackParam *track=0;
 587
 588   // loop on tracks
 589   for(Int_t i=0; i<nTrksOrig; i++) {
 590     AliExternalTrackParam *trackOrig=(AliExternalTrackParam*)trkArrayOrig.At(i);
 591     if(trackOrig->Charge()!=0) { // normal tracks
 592       track = new AliExternalTrackParam(*(AliExternalTrackParam*)trkArrayOrig.At(i));
 593     } else { // neutral tracks (from a V0)
 594       track = new AliNeutralTrackParam(*(AliNeutralTrackParam*)trkArrayOrig.At(i));
 595     }
 596
 597     // tgl cut
 598     if(TMath::Abs(track->GetTgl())>fMaxTgl) {
 599       AliDebug(1,Form(" rejecting track with tgl = %f",track->GetTgl()));
 600       delete track; continue;
 601     }
 602
 603     Bool_t propagateOK = kFALSE, cutond0z0 = kTRUE;
 604     // propagate track to vertex
 605     if(optImpParCut<2 || fOnlyFitter) { // optImpParCut==1 or 0
 606       propagateOK = track->PropagateToDCA(initVertex,GetFieldkG(),100.,d0z0,covd0z0);
 607     } else {              // optImpParCut==2
 608       fCurrentVertex->GetSigmaXYZ(sigmaCurr);
 609       normdistx = TMath::Abs(fCurrentVertex->GetXv()-fNominalPos[0])/TMath::Sqrt(sigmaCurr[0]*sigmaCurr[0]+fNominalCov[0]);
 610       normdisty = TMath::Abs(fCurrentVertex->GetYv()-fNominalPos[1])/TMath::Sqrt(sigmaCurr[1]*sigmaCurr[1]+fNominalCov[2]);
 611       AliDebug(1,Form("normdistx %f  %f    %f",fCurrentVertex->GetXv(),fNominalPos[0],TMath::Sqrt(sigmaCurr[0]*sigmaCurr[0]+fNominalCov[0])));
 612       AliDebug(1,Form("normdisty %f  %f    %f",fCurrentVertex->GetYv(),fNominalPos[1],TMath::Sqrt(sigmaCurr[1]*sigmaCurr[1]+fNominalCov[2])));
 613       AliDebug(1,Form("sigmaCurr %f %f    %f",sigmaCurr[0],sigmaCurr[1],TMath::Sqrt(fNominalCov[0])+TMath::Sqrt(fNominalCov[2])));
 614       if(normdistx < 3. && normdisty < 3. &&
 615          (sigmaCurr[0]+sigmaCurr[1])<(TMath::Sqrt(fNominalCov[0])+TMath::Sqrt(fNominalCov[2]))) {
 616         propagateOK = track->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
 617       } else {
 618         propagateOK = track->PropagateToDCA(initVertex,GetFieldkG(),100.,d0z0,covd0z0);
 619         if(fConstraint) cutond0z0=kFALSE;
 620       }
 621     }
 622
 623     if(!propagateOK) {
 624       AliDebug(1,"     rejected");
 625       delete track; continue;
 626     }
 627
 628     sigmad0 = TMath::Sqrt(covd0z0[0]);
 629     maxd0rphi = fNSigma*sigmad0;
 630     if(optImpParCut==1) maxd0rphi *= 5.;
 631     maxd0rphi = TMath::Min(maxd0rphi,fFiducialR);
 632     //sigmad0z0 = TMath::Sqrt(covd0z0[0]+covd0z0[2]);
 633
 634     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));
 635
 636
 637     //---- track selection based on impact parameters ----//
 638
 639     // always reject tracks outside fiducial volume
 640     if(TMath::Abs(d0z0[0])>fFiducialR || TMath::Abs(d0z0[1])>fFiducialZ) {
 641       AliDebug(1,"     rejected");
 642       delete track; continue;
 643     }
 644
 645     // during iterations 1 and 2 , reject tracks with d0rphi > maxd0rphi
 646     if(optImpParCut>0 && TMath::Abs(d0z0[0])>maxd0rphi) {
 647       AliDebug(1,"     rejected");
 648       delete track; continue;
 649     }
 650
 651     // if fConstraint=kFALSE, during iterations 1 and 2 ||
 652     // if fConstraint=kTRUE, during iteration 2,
 653     // select tracks with d0oplusz0 < fMaxd0z0
 654     if((!fConstraint && optImpParCut>0 && fVert.GetNContributors()>0) ||
 655        ( fConstraint && optImpParCut==2 && cutond0z0)) {
 656       if(nTrksOrig>=3 &&
 657          TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1])>fMaxd0z0) {
 658         AliDebug(1,"     rejected");
 659         delete track; continue;
 660       }
 661     }
 662
 663     // track passed all selections
 664     fTrkArraySel.AddLast(track);
 665     fIdSel[nTrksSel] = idOrig[i];
 666     nTrksSel++;
 667   } // end loop on tracks
 668
 669   delete initVertex;
 670
 671   return nTrksSel;
 672 }
 673 //----------------------------------------------------------------------------
 674 Bool_t AliVertexerTracks::PropagateTrackTo(AliExternalTrackParam *track,
 675                                            Double_t xToGo) {
 676   //----------------------------------------------------------------
 677   // COPIED from AliTracker
 678   //
 679   // Propagates the track to the plane X=xk (cm).
 680   //
 681   //  Origin: Marian Ivanov,  Marian.Ivanov@cern.ch
 682   //----------------------------------------------------------------
 683
 684   const Double_t kEpsilon = 0.00001;
 685   Double_t xpos = track->GetX();
 686   Double_t dir = (xpos<xToGo) ? 1. : -1.;
 687   Double_t maxStep = 5;
 688   Double_t maxSnp = 0.8;
 689   //
 690   while ( (xToGo-xpos)*dir > kEpsilon){
 691     Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
 692     Double_t x    = xpos+step;
 693     Double_t xyz0[3],xyz1[3];
 694     track->GetXYZ(xyz0);   //starting global position
 695
 696     if(!track->GetXYZAt(x,GetFieldkG(),xyz1)) return kFALSE;   // no prolongation
 697     xyz1[2]+=kEpsilon; // waiting for bug correction in geo
 698
 699     if(TMath::Abs(track->GetSnpAt(x,GetFieldkG())) >= maxSnp) return kFALSE;
 700     if(!track->PropagateTo(x,GetFieldkG()))  return kFALSE;
 701
 702     if(TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
 703     track->GetXYZ(xyz0);   // global position
 704     Double_t alphan = TMath::ATan2(xyz0[1], xyz0[0]);
 705     //
 706     Double_t ca=TMath::Cos(alphan-track->GetAlpha()),
 707       sa=TMath::Sin(alphan-track->GetAlpha());
 708     Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf));
 709     Double_t sinNew =  sf*ca - cf*sa;
 710     if(TMath::Abs(sinNew) >= maxSnp) return kFALSE;
 711     if(!track->Rotate(alphan)) return kFALSE;
 712
 713     xpos = track->GetX();
 714   }
 715   return kTRUE;
 716 }
 717 //---------------------------------------------------------------------------
 718 AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
 719                                                         TObjArray *trkArray,
 720                                                         UShort_t *id,
 721                                                         Float_t *diamondxy) const
 722 {
 723 //
 724 // Removes tracks in trksTree from fit of inVtx
 725 //
 726
 727   if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraint")) {
 728     printf("ERROR: primary vertex has no constraint: cannot remove tracks");
 729     return 0x0;
 730   }
 731   if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
 732     printf("WARNING: result of tracks' removal will be only approximately correct");
 733
 734   TMatrixD rv(3,1);
 735   rv(0,0) = inVtx->GetXv();
 736   rv(1,0) = inVtx->GetYv();
 737   rv(2,0) = inVtx->GetZv();
 738   TMatrixD vV(3,3);
 739   Double_t cov[6];
 740   inVtx->GetCovMatrix(cov);
 741   vV(0,0) = cov[0];
 742   vV(0,1) = cov[1]; vV(1,0) = cov[1];
 743   vV(1,1) = cov[2];
 744   vV(0,2) = cov[3]; vV(2,0) = cov[3];
 745   vV(1,2) = cov[4]; vV(2,1) = cov[4];
 746   vV(2,2) = cov[5];
 747
 748   TMatrixD sumWi(TMatrixD::kInverted,vV);
 749   TMatrixD sumWiri(sumWi,TMatrixD::kMult,rv);
 750
 751   Int_t nUsedTrks = inVtx->GetNIndices();
 752   Double_t chi2 = inVtx->GetChi2();
 753
 754   AliExternalTrackParam *track = 0;
 755   Int_t ntrks = (Int_t)trkArray->GetEntriesFast();
 756   for(Int_t i=0;i<ntrks;i++) {
 757     track = (AliExternalTrackParam*)trkArray->At(i);
 758     if(!inVtx->UsesTrack(id[i])) {
 759       printf("track %d was not used in vertex fit",id[i]);
 760       continue;
 761     }
 762     Double_t alpha = track->GetAlpha();
 763     Double_t xl = diamondxy[0]*TMath::Cos(alpha)+diamondxy[1]*TMath::Sin(alpha);
 764     track->PropagateTo(xl,GetFieldkG());
 765     // vector of track global coordinates
 766     TMatrixD ri(3,1);
 767     // covariance matrix of ri
 768     TMatrixD wWi(3,3);
 769
 770     // get space point from track
 771     if(!TrackToPoint(track,ri,wWi)) continue;
 772
 773     TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
 774
 775     sumWi -= wWi;
 776     sumWiri -= wWiri;
 777
 778     // track contribution to chi2
 779     TMatrixD deltar = rv; deltar -= ri;
 780     TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
 781     Double_t chi2i = deltar(0,0)*wWideltar(0,0)+
 782                      deltar(1,0)*wWideltar(1,0)+
 783                      deltar(2,0)*wWideltar(2,0);
 784     // remove from total chi2
 785     chi2 -= chi2i;
 786
 787     nUsedTrks--;
 788     if(nUsedTrks<2) {
 789       printf("Trying to remove too many tracks!");
 790       return 0x0;
 791     }
 792   }
 793
 794   TMatrixD rvnew(3,1);
 795   TMatrixD vVnew(3,3);
 796
 797   // new inverted of weights matrix
 798   TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
 799   vVnew = invsumWi;
 800   // new position of primary vertex
 801   rvnew.Mult(vVnew,sumWiri);
 802
 803   Double_t position[3];
 804   position[0] = rvnew(0,0);
 805   position[1] = rvnew(1,0);
 806   position[2] = rvnew(2,0);
 807   cov[0] = vVnew(0,0);
 808   cov[1] = vVnew(0,1);
 809   cov[2] = vVnew(1,1);
 810   cov[3] = vVnew(0,2);
 811   cov[4] = vVnew(1,2);
 812   cov[5] = vVnew(2,2);
 813
 814   // store data in the vertex object
 815   AliESDVertex *outVtx = new AliESDVertex(position,cov,chi2,nUsedTrks+1); // the +1 is for the constraint
 816   outVtx->SetTitle(inVtx->GetTitle());
 817   UShort_t *inindices = inVtx->GetIndices();
 818   Int_t nIndices = nUsedTrks;
 819   UShort_t *outindices = new UShort_t[nIndices];
 820   Int_t j=0;
 821   for(Int_t k=0; k<inVtx->GetNIndices(); k++) {
 822     Bool_t copyindex=kTRUE;
 823     for(Int_t l=0; l<ntrks; l++) {
 824       if(inindices[k]==id[l]) copyindex=kFALSE;
 825     }
 826     if(copyindex) {
 827       outindices[j] = inindices[k]; j++;
 828     }
 829   }
 830   outVtx->SetIndices(nIndices,outindices);
 831   if (outindices) delete outindices;
 832
 833   /*
 834     printf("Vertex before removing tracks:");
 835     inVtx->PrintStatus();
 836     inVtx->PrintIndices();
 837     printf("Vertex after removing tracks:");
 838     outVtx->PrintStatus();
 839     outVtx->PrintIndices();
 840   */
 841
 842   return outVtx;
 843 }
 844 //---------------------------------------------------------------------------
 845 AliESDVertex* AliVertexerTracks::RemoveConstraintFromVertex(AliESDVertex *inVtx,
 846                                                      Float_t *diamondxyz,
 847                                                      Float_t *diamondcov) const
 848 {
 849 //
 850 // Removes diamond constraint from fit of inVtx
 851 //
 852
 853   if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraint")) {
 854     printf("ERROR: primary vertex has no constraint: cannot remove it\n");
 855     return 0x0;
 856   }
 857   if(inVtx->GetNContributors()<3) {
 858     printf("ERROR: primary vertex has less than 2 tracks: cannot remove contraint\n");
 859     return 0x0;
 860   }
 861
 862   // diamond constraint
 863   TMatrixD vVb(3,3);
 864   vVb(0,0) = diamondcov[0];
 865   vVb(0,1) = diamondcov[1];
 866   vVb(0,2) = 0.;
 867   vVb(1,0) = diamondcov[1];
 868   vVb(1,1) = diamondcov[2];
 869   vVb(1,2) = 0.;
 870   vVb(2,0) = 0.;
 871   vVb(2,1) = 0.;
 872   vVb(2,2) = diamondcov[5];
 873   TMatrixD vVbInv(TMatrixD::kInverted,vVb);
 874   TMatrixD rb(3,1);
 875   rb(0,0) = diamondxyz[0];
 876   rb(1,0) = diamondxyz[1];
 877   rb(2,0) = diamondxyz[2];
 878   TMatrixD vVbInvrb(vVbInv,TMatrixD::kMult,rb);
 879
 880   // input vertex
 881   TMatrixD rv(3,1);
 882   rv(0,0) = inVtx->GetXv();
 883   rv(1,0) = inVtx->GetYv();
 884   rv(2,0) = inVtx->GetZv();
 885   TMatrixD vV(3,3);
 886   Double_t cov[6];
 887   inVtx->GetCovMatrix(cov);
 888   vV(0,0) = cov[0];
 889   vV(0,1) = cov[1]; vV(1,0) = cov[1];
 890   vV(1,1) = cov[2];
 891   vV(0,2) = cov[3]; vV(2,0) = cov[3];
 892   vV(1,2) = cov[4]; vV(2,1) = cov[4];
 893   vV(2,2) = cov[5];
 894   TMatrixD vVInv(TMatrixD::kInverted,vV);
 895   TMatrixD vVInvrv(vVInv,TMatrixD::kMult,rv);
 896
 897
 898   TMatrixD sumWi = vVInv - vVbInv;
 899
 900
 901   TMatrixD sumWiri = vVInvrv - vVbInvrb;
 902
 903   TMatrixD rvnew(3,1);
 904   TMatrixD vVnew(3,3);
 905
 906   // new inverted of weights matrix
 907   TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
 908   vVnew = invsumWi;
 909   // new position of primary vertex
 910   rvnew.Mult(vVnew,sumWiri);
 911
 912   Double_t position[3];
 913   position[0] = rvnew(0,0);
 914   position[1] = rvnew(1,0);
 915   position[2] = rvnew(2,0);
 916   cov[0] = vVnew(0,0);
 917   cov[1] = vVnew(0,1);
 918   cov[2] = vVnew(1,1);
 919   cov[3] = vVnew(0,2);
 920   cov[4] = vVnew(1,2);
 921   cov[5] = vVnew(2,2);
 922
 923
 924   Double_t chi2 = inVtx->GetChi2();
 925
 926   // diamond constribution to chi2
 927   TMatrixD deltar = rv; deltar -= rb;
 928   TMatrixD vVbInvdeltar(vVbInv,TMatrixD::kMult,deltar);
 929   Double_t chi2b = deltar(0,0)*vVbInvdeltar(0,0)+
 930                    deltar(1,0)*vVbInvdeltar(1,0)+
 931                    deltar(2,0)*vVbInvdeltar(2,0);
 932   // remove from total chi2
 933   chi2 -= chi2b;
 934
 935   // store data in the vertex object
 936   AliESDVertex *outVtx = new AliESDVertex(position,cov,chi2,inVtx->GetNContributors()-1);
 937   outVtx->SetTitle("VertexerTracksNoConstraint");
 938   UShort_t *inindices = inVtx->GetIndices();
 939   Int_t nIndices = inVtx->GetNIndices();
 940   outVtx->SetIndices(nIndices,inindices);
 941
 942   return outVtx;
 943 }
 944 //---------------------------------------------------------------------------
 945 void AliVertexerTracks::SetCuts(Double_t *cuts)
 946 {
 947 //
 948 //  Cut values
 949 //
 950   SetDCAcut(cuts[0]);
 951   SetDCAcutIter0(cuts[1]);
 952   SetMaxd0z0(cuts[2]);
 953   if(fMode==0 && cuts[3]<0) SetITSrefitNotRequired();
 954   SetMinClusters((Int_t)(TMath::Abs(cuts[3])));
 955   SetMinTracks((Int_t)(cuts[4]));
 956   SetNSigmad0(cuts[5]);
 957   SetMinDetFitter(cuts[6]);
 958   SetMaxTgl(cuts[7]);
 959   SetFiducialRZ(cuts[8],cuts[9]);
 960   fAlgo=(Int_t)(cuts[10]);
 961   fAlgoIter0=(Int_t)(cuts[11]);
 962
 963   return;
 964 }
 965 //---------------------------------------------------------------------------
 966 void AliVertexerTracks::SetITSMode(Double_t dcacut,
 967                                    Double_t dcacutIter0,
 968                                    Double_t maxd0z0,
 969                                    Int_t minCls,
 970                                    Int_t mintrks,
 971                                    Double_t nsigma,
 972                                    Double_t mindetfitter,
 973                                    Double_t maxtgl,
 974                                    Double_t fidR,
 975                                    Double_t fidZ,
 976                                    Int_t finderAlgo,
 977                                    Int_t finderAlgoIter0)
 978 {
 979 //
 980 //  Cut values for ITS mode
 981 //
 982   fMode = 0;
 983   if(minCls>0) {
 984     SetITSrefitRequired();
 985   } else {
 986     SetITSrefitNotRequired();
 987   }
 988   SetDCAcut(dcacut);
 989   SetDCAcutIter0(dcacutIter0);
 990   SetMaxd0z0(maxd0z0);
 991   SetMinClusters(TMath::Abs(minCls));
 992   SetMinTracks(mintrks);
 993   SetNSigmad0(nsigma);
 994   SetMinDetFitter(mindetfitter);
 995   SetMaxTgl(maxtgl);
 996   SetFiducialRZ(fidR,fidZ);
 997   fAlgo=finderAlgo;
 998   fAlgoIter0=finderAlgoIter0;
 999
1000   return;
1001 }
1002 //---------------------------------------------------------------------------
1003 void AliVertexerTracks::SetTPCMode(Double_t dcacut,
1004                                    Double_t dcacutIter0,
1005                                    Double_t maxd0z0,
1006                                    Int_t minCls,
1007                                    Int_t mintrks,
1008                                    Double_t nsigma,
1009                                    Double_t mindetfitter,
1010                                    Double_t maxtgl,
1011                                    Double_t fidR,
1012                                    Double_t fidZ,
1013                                    Int_t finderAlgo,
1014                                    Int_t finderAlgoIter0)
1015 {
1016 //
1017 //  Cut values for TPC mode
1018 //
1019   fMode = 1;
1020   SetITSrefitNotRequired();
1021   SetDCAcut(dcacut);
1022   SetDCAcutIter0(dcacutIter0);
1023   SetMaxd0z0(maxd0z0);
1024   SetMinClusters(minCls);
1025   SetMinTracks(mintrks);
1026   SetNSigmad0(nsigma);
1027   SetMinDetFitter(mindetfitter);
1028   SetMaxTgl(maxtgl);
1029   SetFiducialRZ(fidR,fidZ);
1030   fAlgo=finderAlgo;
1031   fAlgoIter0=finderAlgoIter0;
1032
1033   return;
1034 }
1035 //---------------------------------------------------------------------------
1036 void AliVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped)
1037 {
1038 //
1039 // Mark the tracks not to be used in the vertex reconstruction.
1040 // Tracks are identified by AliESDtrack::GetID()
1041 //
1042   fNTrksToSkip = n;  fTrksToSkip = new Int_t[n];
1043   for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i];
1044   return;
1045 }
1046 //---------------------------------------------------------------------------
1047 void  AliVertexerTracks::SetVtxStart(AliESDVertex *vtx)
1048 {
1049 //
1050 // Set initial vertex knowledge
1051 //
1052   vtx->GetXYZ(fNominalPos);
1053   vtx->GetCovMatrix(fNominalCov);
1054   SetConstraintOn();
1055   return;
1056 }
1057 //---------------------------------------------------------------------------
1058 void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
1059 {
1060   AliExternalTrackParam *track1;
1061   const Int_t knacc = (Int_t)fTrkArraySel.GetEntriesFast();
1062   AliStrLine **linarray = new AliStrLine* [knacc];
1063   for(Int_t i=0; i<knacc; i++){
1064     track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
1065     Double_t alpha=track1->GetAlpha();
1066     Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
1067     Double_t pos[3],dir[3],sigmasq[3];
1068     track1->GetXYZAt(mindist,GetFieldkG(),pos);
1069     track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
1070     sigmasq[0]=TMath::Sin(alpha)*TMath::Sin(alpha)*track1->GetSigmaY2();
1071     sigmasq[1]=TMath::Cos(alpha)*TMath::Cos(alpha)*track1->GetSigmaY2();
1072     sigmasq[2]=track1->GetSigmaZ2();
1073     TMatrixD ri(3,1);
1074     TMatrixD wWi(3,3);
1075     if(!TrackToPoint(track1,ri,wWi)) {optUseWeights=kFALSE;printf("WARNING\n");}
1076     Double_t wmat[9];
1077     Int_t iel=0;
1078     for(Int_t ia=0;ia<3;ia++){
1079       for(Int_t ib=0;ib<3;ib++){
1080         wmat[iel]=wWi(ia,ib);
1081         iel++;
1082       }
1083     }
1084     linarray[i] = new AliStrLine(pos,sigmasq,wmat,dir);
1085   }
1086   fVert=TrackletVertexFinder(linarray,knacc,optUseWeights);
1087   for(Int_t i=0; i<knacc; i++) delete linarray[i];
1088   delete [] linarray;
1089 }
1090 //---------------------------------------------------------------------------
1091 AliESDVertex AliVertexerTracks::TrackletVertexFinder(TClonesArray *lines, Int_t optUseWeights)
1092 {
1093   // Calculate the point at minimum distance to prepared tracks (TClonesArray)
1094   const Int_t knacc = (Int_t)lines->GetEntriesFast();
1095   AliStrLine** lines2 = new AliStrLine* [knacc];
1096   for(Int_t i=0; i<knacc; i++){
1097     lines2[i]= (AliStrLine*)lines->At(i);
1098   }
1099   AliESDVertex vert = TrackletVertexFinder(lines2,knacc,optUseWeights);
1100   delete [] lines2;
1101   return vert;
1102 }
1103
1104 //---------------------------------------------------------------------------
1105 AliESDVertex AliVertexerTracks::TrackletVertexFinder(AliStrLine **lines, const Int_t knacc, Int_t optUseWeights)
1106 {
1107   // Calculate the point at minimum distance to prepared tracks (array of AliStrLine)
1108
1109   Double_t initPos[3]={0.,0.,0.};
1110
1111   Double_t (*vectP0)[3]=new Double_t [knacc][3];
1112   Double_t (*vectP1)[3]=new Double_t [knacc][3];
1113
1114   Double_t sum[3][3];
1115   Double_t dsum[3]={0,0,0};
1116   TMatrixD sumWi(3,3);
1117   for(Int_t i=0;i<3;i++){
1118     for(Int_t j=0;j<3;j++){
1119       sum[i][j]=0;
1120       sumWi(i,j)=0.;
1121     }
1122   }
1123
1124   for(Int_t i=0; i<knacc; i++){
1125     AliStrLine *line1 = lines[i];
1126     Double_t p0[3],cd[3],sigmasq[3];
1127     Double_t wmat[9];
1128     if(!line1) printf("ERROR %d %d\n",i,knacc);
1129     line1->GetP0(p0);
1130     line1->GetCd(cd);
1131     line1->GetSigma2P0(sigmasq);
1132     line1->GetWMatrix(wmat);
1133     TMatrixD wWi(3,3);
1134     Int_t iel=0;
1135     for(Int_t ia=0;ia<3;ia++){
1136       for(Int_t ib=0;ib<3;ib++){
1137         wWi(ia,ib)=wmat[iel];
1138         iel++;
1139       }
1140     }
1141
1142     sumWi+=wWi;
1143
1144     Double_t p1[3]={p0[0]+cd[0],p0[1]+cd[1],p0[2]+cd[2]};
1145     vectP0[i][0]=p0[0];
1146     vectP0[i][1]=p0[1];
1147     vectP0[i][2]=p0[2];
1148     vectP1[i][0]=p1[0];
1149     vectP1[i][1]=p1[1];
1150     vectP1[i][2]=p1[2];
1151
1152     Double_t matr[3][3];
1153     Double_t dknow[3];
1154     if(optUseWeights==0)GetStrLinDerivMatrix(p0,p1,matr,dknow);
1155     else GetStrLinDerivMatrix(p0,p1,sigmasq,matr,dknow);
1156
1157
1158     for(Int_t iii=0;iii<3;iii++){
1159       dsum[iii]+=dknow[iii];
1160       for(Int_t lj=0;lj<3;lj++) sum[iii][lj]+=matr[iii][lj];
1161     }
1162   }
1163
1164   TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
1165   Double_t covmatrix[6];
1166   covmatrix[0] = invsumWi(0,0);
1167   covmatrix[1] = invsumWi(0,1);
1168   covmatrix[2] = invsumWi(1,1);
1169   covmatrix[3] = invsumWi(0,2);
1170   covmatrix[4] = invsumWi(1,2);
1171   covmatrix[5] = invsumWi(2,2);
1172
1173   Double_t vett[3][3];
1174   Double_t det=GetDeterminant3X3(sum);
1175   Double_t sigma=0;
1176
1177   if(det!=0){
1178     for(Int_t zz=0;zz<3;zz++){
1179       for(Int_t ww=0;ww<3;ww++){
1180         for(Int_t kk=0;kk<3;kk++) vett[ww][kk]=sum[ww][kk];
1181       }
1182       for(Int_t kk=0;kk<3;kk++) vett[kk][zz]=dsum[kk];
1183       initPos[zz]=GetDeterminant3X3(vett)/det;
1184     }
1185
1186
1187     for(Int_t i=0; i<knacc; i++){
1188       Double_t p0[3]={0,0,0},p1[3]={0,0,0};
1189       for(Int_t ii=0;ii<3;ii++){
1190         p0[ii]=vectP0[i][ii];
1191         p1[ii]=vectP1[i][ii];
1192       }
1193       sigma+=GetStrLinMinDist(p0,p1,initPos);
1194     }
1195
1196     if(sigma>0.) {sigma=TMath::Sqrt(sigma);}else{sigma=999;}
1197   }else{
1198     sigma=999;
1199   }
1200   AliESDVertex theVert(initPos,covmatrix,99999.,knacc);
1201   theVert.SetDispersion(sigma);
1202   delete vectP0;
1203   delete vectP1;
1204   return theVert;
1205 }
1206
1207 //---------------------------------------------------------------------------
1208 Bool_t AliVertexerTracks::TrackToPoint(AliExternalTrackParam *t,
1209                                        TMatrixD &ri,TMatrixD &wWi,
1210                                        Bool_t uUi3by3) const
1211 {
1212 //
1213 // Extract from the AliExternalTrackParam the global coordinates ri and covariance matrix
1214 // wWi of the space point that it represents (to be used in VertexFitter())
1215 //
1216
1217
1218   Double_t rotAngle = t->GetAlpha();
1219   if(rotAngle<0.) rotAngle += 2.*TMath::Pi();
1220   Double_t cosRot = TMath::Cos(rotAngle);
1221   Double_t sinRot = TMath::Sin(rotAngle);
1222
1223   ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot;
1224   ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot;
1225   ri(2,0) = t->GetZ();
1226
1227
1228
1229   if(!uUi3by3) {
1230     // matrix to go from global (x,y,z) to local (y,z);
1231     TMatrixD qQi(2,3);
1232     qQi(0,0) = -sinRot;
1233     qQi(0,1) = cosRot;
1234     qQi(0,2) = 0.;
1235     qQi(1,0) = 0.;
1236     qQi(1,1) = 0.;
1237     qQi(1,2) = 1.;
1238
1239     // covariance matrix of local (y,z) - inverted
1240     TMatrixD uUi(2,2);
1241     uUi(0,0) = t->GetSigmaY2();
1242     uUi(0,1) = t->GetSigmaZY();
1243     uUi(1,0) = t->GetSigmaZY();
1244     uUi(1,1) = t->GetSigmaZ2();
1245     //printf(" Ui :");
1246     //printf(" %f   %f",uUi(0,0),uUi(0,1));
1247     //printf(" %f   %f",uUi(1,0),uUi(1,1));
1248
1249     if(uUi.Determinant() <= 0.) return kFALSE;
1250     TMatrixD uUiInv(TMatrixD::kInverted,uUi);
1251
1252     // weights matrix: wWi = qQiT * uUiInv * qQi
1253     TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
1254     TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
1255     wWi = m;
1256   } else {
1257     if(fVert.GetNContributors()<1) AliFatal("Vertex from finder is empty");
1258     // matrix to go from global (x,y,z) to local (x,y,z);
1259     TMatrixD qQi(3,3);
1260     qQi(0,0) = cosRot;
1261     qQi(0,1) = sinRot;
1262     qQi(0,2) = 0.;
1263     qQi(1,0) = -sinRot;
1264     qQi(1,1) = cosRot;
1265     qQi(1,2) = 0.;
1266     qQi(2,0) = 0.;
1267     qQi(2,1) = 0.;
1268     qQi(2,2) = 1.;
1269
1270     // covariance of fVert along the track
1271     Double_t p[3],pt,ptot;
1272     t->GetPxPyPz(p);
1273     pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
1274     ptot = TMath::Sqrt(pt*pt+p[2]*p[2]);
1275     Double_t cphi = p[0]/pt;               //cos(phi)=px/pt
1276     Double_t sphi = p[1]/pt;               //sin(phi)=py/pt
1277     Double_t clambda = pt/ptot;            //cos(lambda)=pt/ptot
1278     Double_t slambda = p[2]/ptot;            //sin(lambda)=pz/ptot
1279     Double_t covfVert[6];
1280     fVert.GetCovMatrix(covfVert);
1281     Double_t covfVertalongt =
1282        covfVert[0]*cphi*cphi*clambda*clambda
1283       +covfVert[1]*2.*cphi*sphi*clambda*clambda
1284       +covfVert[3]*2.*cphi*clambda*slambda
1285       +covfVert[2]*sphi*sphi*clambda*clambda
1286       +covfVert[4]*2.*sphi*clambda*slambda
1287       +covfVert[5]*slambda*slambda;
1288     // covariance matrix of local (x,y,z) - inverted
1289     TMatrixD uUi(3,3);
1290     uUi(0,0) = covfVertalongt * fnSigmaForUi00 * fnSigmaForUi00;
1291     AliDebug(1,Form("=====> sqrtUi00 cm  %f",TMath::Sqrt(uUi(0,0))));
1292     uUi(0,1) = 0.;
1293     uUi(0,2) = 0.;
1294     uUi(1,0) = 0.;
1295     uUi(1,1) = t->GetSigmaY2();
1296     uUi(1,2) = t->GetSigmaZY();
1297     uUi(2,0) = 0.;
1298     uUi(2,1) = t->GetSigmaZY();
1299     uUi(2,2) = t->GetSigmaZ2();
1300     //printf(" Ui :\n");
1301     //printf(" %f   %f\n",uUi(0,0),uUi(0,1));
1302     //printf(" %f   %f\n",uUi(1,0),uUi(1,1));
1303
1304     if(uUi.Determinant() <= 0.) return kFALSE;
1305     TMatrixD uUiInv(TMatrixD::kInverted,uUi);
1306
1307     // weights matrix: wWi = qQiT * uUiInv * qQi
1308     TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
1309     TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
1310     wWi = m;
1311   }
1312
1313
1314   return kTRUE;
1315 }
1316 //---------------------------------------------------------------------------
1317 void AliVertexerTracks::TooFewTracks()
1318 {
1319 //
1320 // When the number of tracks is < fMinTracks,
1321 // deal with vertices not found and prepare to exit
1322 //
1323   AliDebug(1,"TooFewTracks");
1324
1325   Double_t pos[3],err[3];
1326   pos[0] = fNominalPos[0];
1327   err[0] = TMath::Sqrt(fNominalCov[0]);
1328   pos[1] = fNominalPos[1];
1329   err[1] = TMath::Sqrt(fNominalCov[2]);
1330   pos[2] = fNominalPos[2];
1331   err[2] = TMath::Sqrt(fNominalCov[5]);
1332   Int_t    ncontr = (err[0]>1. ? -1 : -3);
1333   if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
1334   fCurrentVertex = new AliESDVertex(pos,err);
1335   fCurrentVertex->SetNContributors(ncontr);
1336
1337   if(fConstraint) {
1338     fCurrentVertex->SetTitle("VertexerTracksWithConstraint");
1339   } else {
1340     fCurrentVertex->SetTitle("VertexerTracksNoConstraint");
1341   }
1342
1343   if(!fTrkArraySel.IsEmpty()) fTrkArraySel.Delete();
1344   if(fIdSel) {delete fIdSel; fIdSel=NULL;}
1345   if(fTrksToSkip) {delete fTrksToSkip; fTrksToSkip=NULL;}
1346
1347   return;
1348 }
1349 //---------------------------------------------------------------------------
1350 void AliVertexerTracks::VertexFinder(Int_t optUseWeights)
1351 {
1352
1353   // Get estimate of vertex position in (x,y) from tracks DCA
1354
1355   Double_t initPos[3];
1356   initPos[2] = 0.;
1357   for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
1358   Int_t nacc = (Int_t)fTrkArraySel.GetEntriesFast();
1359   Double_t aver[3]={0.,0.,0.};
1360   Double_t aversq[3]={0.,0.,0.};
1361   Double_t sigmasq[3]={0.,0.,0.};
1362   Double_t sigma=0;
1363   Int_t ncombi = 0;
1364   AliExternalTrackParam *track1;
1365   AliExternalTrackParam *track2;
1366   Double_t pos[3],dir[3];
1367   Double_t alpha,mindist;
1368
1369   for(Int_t i=0; i<nacc; i++){
1370     track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
1371     alpha=track1->GetAlpha();
1372     mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
1373     track1->GetXYZAt(mindist,GetFieldkG(),pos);
1374     track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
1375     AliStrLine *line1 = new AliStrLine(pos,dir);
1376
1377    //    AliStrLine *line1 = new AliStrLine();
1378    //    track1->ApproximateHelixWithLine(mindist,GetFieldkG(),line1);
1379
1380     for(Int_t j=i+1; j<nacc; j++){
1381       track2 = (AliExternalTrackParam*)fTrkArraySel.At(j);
1382       alpha=track2->GetAlpha();
1383       mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
1384       track2->GetXYZAt(mindist,GetFieldkG(),pos);
1385       track2->GetPxPyPzAt(mindist,GetFieldkG(),dir);
1386       AliStrLine *line2 = new AliStrLine(pos,dir);
1387     //      AliStrLine *line2 = new AliStrLine();
1388     //  track2->ApproximateHelixWithLine(mindist,GetFieldkG(),line2);
1389       Double_t distCA=line2->GetDCA(line1);
1390       //printf("%d   %d   %f\n",i,j,distCA);
1391        if(fDCAcut<=0 || (fDCAcut>0&&distCA<fDCAcut)){
1392         Double_t pnt1[3],pnt2[3],crosspoint[3];
1393
1394         if(optUseWeights<=0){
1395           Int_t retcode = line2->Cross(line1,crosspoint);
1396           if(retcode>=0){
1397             ncombi++;
1398             for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
1399             for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
1400           }
1401         }
1402         if(optUseWeights>0){
1403           Int_t retcode = line1->CrossPoints(line2,pnt1,pnt2);
1404           if(retcode>=0){
1405             Double_t cs, sn;
1406             alpha=track1->GetAlpha();
1407             cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
1408             Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
1409             alpha=track2->GetAlpha();
1410             cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
1411             Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
1412             Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
1413             Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
1414             Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
1415             Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
1416             crosspoint[0]=wx1*pnt1[0] + wx2*pnt2[0];
1417             crosspoint[1]=wy1*pnt1[1] + wy2*pnt2[1];
1418             crosspoint[2]=wz1*pnt1[2] + wz2*pnt2[2];
1419
1420             ncombi++;
1421             for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
1422             for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
1423           }
1424         }
1425       }
1426       delete line2;
1427     }
1428     delete line1;
1429   }
1430   if(ncombi>0){
1431     for(Int_t jj=0;jj<3;jj++){
1432       initPos[jj] = aver[jj]/ncombi;
1433       //printf("%f\n",initPos[jj]);
1434       aversq[jj]/=ncombi;
1435       sigmasq[jj]=aversq[jj]-initPos[jj]*initPos[jj];
1436       sigma+=sigmasq[jj];
1437     }
1438     sigma=TMath::Sqrt(TMath::Abs(sigma));
1439   }
1440   else {
1441     Warning("VertexFinder","Finder did not succed");
1442     sigma=999;
1443   }
1444   fVert.SetXYZ(initPos);
1445   fVert.SetDispersion(sigma);
1446   fVert.SetNContributors(ncombi);
1447 }
1448 //---------------------------------------------------------------------------
1449 void AliVertexerTracks::VertexFitter()
1450 {
1451 //
1452 // The optimal estimate of the vertex position is given by a "weighted
1453 // average of tracks positions".
1454 // Original method: V. Karimaki, CMS Note 97/0051
1455 //
1456   Bool_t useConstraint = fConstraint;
1457   Double_t initPos[3];
1458   if(!fOnlyFitter) {
1459     fVert.GetXYZ(initPos);
1460   } else {
1461     initPos[0]=fNominalPos[0];
1462     initPos[1]=fNominalPos[1];
1463     initPos[2]=fNominalPos[2];
1464   }
1465
1466   Int_t nTrksSel = (Int_t)fTrkArraySel.GetEntries();
1467   if(nTrksSel==1) useConstraint=kTRUE;
1468   AliDebug(1,"--- VertexFitter(): start");
1469   AliDebug(1,Form(" Number of tracks in array: %d\n",nTrksSel));
1470   AliDebug(1,Form(" Minimum # tracks required in fit: %d\n",fMinTracks));
1471   AliDebug(1,Form(" Vertex position after finder: %f,%f,%f\n",initPos[0],initPos[1],initPos[2]));
1472   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])));
1473
1474   // special treatment for few-tracks fits (e.g. secondary vertices)
1475   Bool_t uUi3by3 = kFALSE; if(nTrksSel<5 && !useConstraint) uUi3by3 = kTRUE;
1476
1477   Int_t i,j,k,step=0;
1478   TMatrixD rv(3,1);
1479   TMatrixD vV(3,3);
1480   rv(0,0) = initPos[0];
1481   rv(1,0) = initPos[1];
1482   rv(2,0) = 0.;
1483   Double_t xlStart,alpha;
1484   Int_t nTrksUsed;
1485   Double_t chi2,chi2i,chi2b;
1486   AliExternalTrackParam *t = 0;
1487   Int_t failed = 0;
1488
1489   // initial vertex covariance matrix
1490   TMatrixD vVb(3,3);
1491   vVb(0,0) = fNominalCov[0];
1492   vVb(0,1) = fNominalCov[1];
1493   vVb(0,2) = 0.;
1494   vVb(1,0) = fNominalCov[1];
1495   vVb(1,1) = fNominalCov[2];
1496   vVb(1,2) = 0.;
1497   vVb(2,0) = 0.;
1498   vVb(2,1) = 0.;
1499   vVb(2,2) = fNominalCov[5];
1500   TMatrixD vVbInv(TMatrixD::kInverted,vVb);
1501   TMatrixD rb(3,1);
1502   rb(0,0) = fNominalPos[0];
1503   rb(1,0) = fNominalPos[1];
1504   rb(2,0) = fNominalPos[2];
1505   TMatrixD vVbInvrb(vVbInv,TMatrixD::kMult,rb);
1506
1507
1508   // 2 steps:
1509   // 1st - estimate of vtx using all tracks
1510   // 2nd - estimate of global chi2
1511   for(step=0; step<2; step++) {
1512     AliDebug(1,Form(" step = %d\n",step));
1513     chi2 = 0.;
1514     nTrksUsed = 0;
1515
1516     if(step==1) { initPos[0]=rv(0,0); initPos[0]=rv(1,0); }
1517
1518     TMatrixD sumWiri(3,1);
1519     TMatrixD sumWi(3,3);
1520     for(i=0; i<3; i++) {
1521       sumWiri(i,0) = 0.;
1522       for(j=0; j<3; j++) sumWi(j,i) = 0.;
1523     }
1524
1525     // mean vertex constraint
1526     if(useConstraint) {
1527       for(i=0;i<3;i++) {
1528         sumWiri(i,0) += vVbInvrb(i,0);
1529         for(k=0;k<3;k++) sumWi(i,k) += vVbInv(i,k);
1530       }
1531       // chi2
1532       TMatrixD deltar = rv; deltar -= rb;
1533       TMatrixD vVbInvdeltar(vVbInv,TMatrixD::kMult,deltar);
1534       chi2b = deltar(0,0)*vVbInvdeltar(0,0)+
1535               deltar(1,0)*vVbInvdeltar(1,0)+
1536               deltar(2,0)*vVbInvdeltar(2,0);
1537       chi2 += chi2b;
1538     }
1539
1540     // loop on tracks
1541     for(k=0; k<nTrksSel; k++) {
1542       // get track from track array
1543       t = (AliExternalTrackParam*)fTrkArraySel.At(k);
1544       alpha = t->GetAlpha();
1545       xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha);
1546       // to vtxSeed (from finder)
1547       t->PropagateTo(xlStart,GetFieldkG());
1548
1549       // vector of track global coordinates
1550       TMatrixD ri(3,1);
1551       // covariance matrix of ri
1552       TMatrixD wWi(3,3);
1553
1554       // get space point from track
1555       if(!TrackToPoint(t,ri,wWi,uUi3by3)) continue;
1556       TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
1557
1558       // track chi2
1559       TMatrixD deltar = rv; deltar -= ri;
1560       TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
1561       chi2i = deltar(0,0)*wWideltar(0,0)+
1562               deltar(1,0)*wWideltar(1,0)+
1563               deltar(2,0)*wWideltar(2,0);
1564
1565       // add to total chi2
1566       chi2 += chi2i;
1567
1568       sumWiri += wWiri;
1569       sumWi   += wWi;
1570
1571       nTrksUsed++;
1572     } // end loop on tracks
1573
1574     if(nTrksUsed < fMinTracks) {
1575       failed=1;
1576       continue;
1577     }
1578
1579     Double_t determinant = sumWi.Determinant();
1580     if(determinant < fMinDetFitter)  {
1581       AliDebug(1,Form("det(V) = %f (<%f)\n",determinant,fMinDetFitter));
1582       failed=1;
1583       continue;
1584     }
1585
1586     if(step==0) {
1587       // inverted of weights matrix
1588       TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
1589       vV = invsumWi;
1590       // position of primary vertex
1591       rv.Mult(vV,sumWiri);
1592     }
1593   } // end loop on the 2 steps
1594
1595   if(failed) {
1596     TooFewTracks();
1597     return;
1598   }
1599
1600   Double_t position[3];
1601   position[0] = rv(0,0);
1602   position[1] = rv(1,0);
1603   position[2] = rv(2,0);
1604   Double_t covmatrix[6];
1605   covmatrix[0] = vV(0,0);
1606   covmatrix[1] = vV(0,1);
1607   covmatrix[2] = vV(1,1);
1608   covmatrix[3] = vV(0,2);
1609   covmatrix[4] = vV(1,2);
1610   covmatrix[5] = vV(2,2);
1611
1612   // for correct chi2/ndf, count constraint as additional "track"
1613   if(fConstraint) nTrksUsed++;
1614   // store data in the vertex object
1615   if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
1616   fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);
1617
1618   AliDebug(1," Vertex after fit:");
1619   AliDebug(1,Form("xyz: %f %f %f; nc %d",fCurrentVertex->GetXv(),fCurrentVertex->GetYv(),fCurrentVertex->GetZv(),fCurrentVertex->GetNContributors()));
1620   AliDebug(1,"--- VertexFitter(): finish\n");
1621
1622
1623   return;
1624 }
1625 //----------------------------------------------------------------------------
1626 AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray,
1627                                                          UShort_t *id,
1628                                                          Bool_t optUseFitter,
1629                                                          Bool_t optPropagate,
1630                                                          Bool_t optUseDiamondConstraint)
1631 {
1632 //
1633 // Return vertex from tracks (AliExternalTrackParam) in array
1634 //
1635   fCurrentVertex = 0;
1636   // set optUseDiamondConstraint=TRUE only if you are reconstructing the
1637   // primary vertex!
1638   if(optUseDiamondConstraint) {
1639     SetConstraintOn();
1640   } else {
1641     SetConstraintOff();
1642   }
1643
1644   // get tracks and propagate them to initial vertex position
1645   fIdSel = new UShort_t[(Int_t)trkArray->GetEntriesFast()];
1646   Int_t nTrksSel = PrepareTracks(*trkArray,id,0);
1647   if((!optUseDiamondConstraint && nTrksSel<TMath::Max(2,fMinTracks)) ||
1648      (optUseDiamondConstraint && nTrksSel<1)) {
1649     TooFewTracks();
1650     return fCurrentVertex;
1651   }
1652
1653   // vertex finder
1654   if(nTrksSel==1) {
1655     AliDebug(1,"Just one track");
1656     OneTrackVertFinder();
1657   } else {
1658     switch (fAlgo) {
1659       case 1: StrLinVertexFinderMinDist(1); break;
1660       case 2: StrLinVertexFinderMinDist(0); break;
1661       case 3: HelixVertexFinder();          break;
1662       case 4: VertexFinder(1);              break;
1663       case 5: VertexFinder(0);              break;
1664       default: printf("Wrong algorithm\n"); break;
1665     }
1666   }
1667   AliDebug(1," Vertex finding completed\n");
1668
1669   // vertex fitter
1670   if(optUseFitter) {
1671     VertexFitter();
1672   } else {
1673     Double_t position[3]={fVert.GetXv(),fVert.GetYv(),fVert.GetZv()};
1674     Double_t covmatrix[6];
1675     fVert.GetCovMatrix(covmatrix);
1676     Double_t chi2=99999.;
1677     Int_t    nTrksUsed=fVert.GetNContributors();
1678     fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);
1679   }
1680   fCurrentVertex->SetDispersion(fVert.GetDispersion());
1681
1682
1683   // set indices of used tracks and propagate track to found vertex
1684   UShort_t *indices = 0;
1685   Double_t d0z0[2],covd0z0[3];
1686   AliExternalTrackParam *t = 0;
1687   if(fCurrentVertex->GetNContributors()>0) {
1688     indices = new UShort_t[fTrkArraySel.GetEntriesFast()];
1689     for(Int_t jj=0; jj<(Int_t)fTrkArraySel.GetEntriesFast(); jj++) {
1690       indices[jj] = fIdSel[jj];
1691       t = (AliExternalTrackParam*)fTrkArraySel.At(jj);
1692       if(optPropagate && optUseFitter) {
1693         if(TMath::Sqrt(fCurrentVertex->GetXv()*fCurrentVertex->GetXv()+fCurrentVertex->GetYv()*fCurrentVertex->GetYv())<3.) {
1694           t->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
1695           AliDebug(1,Form("Track %d propagated to found vertex",jj));
1696         } else {
1697           AliWarning("Found vertex outside beam pipe!");
1698         }
1699       }
1700     }
1701     fCurrentVertex->SetIndices(fCurrentVertex->GetNContributors(),indices);
1702   }
1703
1704   // clean up
1705   if (indices) {delete indices; indices=NULL;}
1706   delete fIdSel; fIdSel=NULL;
1707   fTrkArraySel.Delete();
1708
1709   return fCurrentVertex;
1710 }
1711 //----------------------------------------------------------------------------
1712 AliESDVertex* AliVertexerTracks::VertexForSelectedESDTracks(TObjArray *trkArray,
1713                                                             Bool_t optUseFitter,
1714                                                             Bool_t optPropagate,
1715                                                             Bool_t optUseDiamondConstraint)
1716
1717 {
1718 //
1719 // Return vertex from array of ESD tracks
1720 //
1721
1722   Int_t nTrks = (Int_t)trkArray->GetEntriesFast();
1723   UShort_t *id = new UShort_t[nTrks];
1724
1725   AliESDtrack *esdt = 0;
1726   for(Int_t i=0; i<nTrks; i++){
1727     esdt = (AliESDtrack*)trkArray->At(i);
1728     id[i] = (UShort_t)esdt->GetID();
1729   }
1730
1731   VertexForSelectedTracks(trkArray,id,optUseFitter,optPropagate,optUseDiamondConstraint);
1732
1733   delete id; id=NULL;
1734
1735   return fCurrentVertex;
1736 }
1737 //--------------------------------------------------------------------------