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;
 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       if(normdistx < 3. && normdisty < 3. &&
 612          (sigmaCurr[0]+sigmaCurr[1])<(TMath::Sqrt(fNominalCov[0])+TMath::Sqrt(fNominalCov[2]))) {
 613         propagateOK = track->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
 614       } else {
 615         propagateOK = track->PropagateToDCA(initVertex,GetFieldkG(),100.,d0z0,covd0z0);
 616       }
 617     }
 618
 619     if(!propagateOK) {
 620       AliDebug(1,"     rejected");
 621       delete track; continue;
 622     }
 623
 624     sigmad0 = TMath::Sqrt(covd0z0[0]);
 625     maxd0rphi = fNSigma*sigmad0;
 626     if(optImpParCut==1) maxd0rphi *= 5.;
 627     maxd0rphi = TMath::Min(maxd0rphi,fFiducialR);
 628     //sigmad0z0 = TMath::Sqrt(covd0z0[0]+covd0z0[2]); // for future improvement
 629     //maxd0z0 = 10.*fNSigma*sigmad0z0;
 630
 631     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));
 632
 633
 634     //---- track selection based on impact parameters ----//
 635
 636     // always reject tracks outside fiducial volume
 637     if(TMath::Abs(d0z0[0])>fFiducialR || TMath::Abs(d0z0[1])>fFiducialZ) {
 638       AliDebug(1,"     rejected");
 639       delete track; continue;
 640     }
 641
 642     // during iterations 1 and 2 , reject tracks with d0rphi > maxd0rphi
 643     if(optImpParCut>0 && TMath::Abs(d0z0[0])>maxd0rphi) {
 644       AliDebug(1,"     rejected");
 645       delete track; continue;
 646     }
 647
 648     // if fConstraint=kFALSE, during iterations 1 and 2 ||
 649     // if fConstraint=kTRUE, during iteration 2,
 650     // select tracks with d0oplusz0 < fMaxd0z0
 651     if((!fConstraint && optImpParCut>0 && fVert.GetNContributors()>0) ||
 652        ( fConstraint && optImpParCut==2)) {
 653       if(nTrksOrig>=3 &&
 654          TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1])>fMaxd0z0) {
 655         AliDebug(1,"     rejected");
 656         delete track; continue;
 657       }
 658     }
 659
 660     // track passed all selections
 661     fTrkArraySel.AddLast(track);
 662     fIdSel[nTrksSel] = idOrig[i];
 663     nTrksSel++;
 664   } // end loop on tracks
 665
 666   delete initVertex;
 667
 668   return nTrksSel;
 669 }
 670 //----------------------------------------------------------------------------
 671 Bool_t AliVertexerTracks::PropagateTrackTo(AliExternalTrackParam *track,
 672                                            Double_t xToGo) {
 673   //----------------------------------------------------------------
 674   // COPIED from AliTracker
 675   //
 676   // Propagates the track to the plane X=xk (cm).
 677   //
 678   //  Origin: Marian Ivanov,  Marian.Ivanov@cern.ch
 679   //----------------------------------------------------------------
 680
 681   const Double_t kEpsilon = 0.00001;
 682   Double_t xpos = track->GetX();
 683   Double_t dir = (xpos<xToGo) ? 1. : -1.;
 684   Double_t maxStep = 5;
 685   Double_t maxSnp = 0.8;
 686   //
 687   while ( (xToGo-xpos)*dir > kEpsilon){
 688     Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
 689     Double_t x    = xpos+step;
 690     Double_t xyz0[3],xyz1[3];
 691     track->GetXYZ(xyz0);   //starting global position
 692
 693     if(!track->GetXYZAt(x,GetFieldkG(),xyz1)) return kFALSE;   // no prolongation
 694     xyz1[2]+=kEpsilon; // waiting for bug correction in geo
 695
 696     if(TMath::Abs(track->GetSnpAt(x,GetFieldkG())) >= maxSnp) return kFALSE;
 697     if(!track->PropagateTo(x,GetFieldkG()))  return kFALSE;
 698
 699     if(TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
 700     track->GetXYZ(xyz0);   // global position
 701     Double_t alphan = TMath::ATan2(xyz0[1], xyz0[0]);
 702     //
 703     Double_t ca=TMath::Cos(alphan-track->GetAlpha()),
 704       sa=TMath::Sin(alphan-track->GetAlpha());
 705     Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf));
 706     Double_t sinNew =  sf*ca - cf*sa;
 707     if(TMath::Abs(sinNew) >= maxSnp) return kFALSE;
 708     if(!track->Rotate(alphan)) return kFALSE;
 709
 710     xpos = track->GetX();
 711   }
 712   return kTRUE;
 713 }
 714 //---------------------------------------------------------------------------
 715 AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
 716                                                         TObjArray *trkArray,
 717                                                         UShort_t *id,
 718                                                         Float_t *diamondxy) const
 719 {
 720 //
 721 // Removes tracks in trksTree from fit of inVtx
 722 //
 723
 724   if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraint")) {
 725     printf("ERROR: primary vertex has no constraint: cannot remove tracks");
 726     return 0x0;
 727   }
 728   if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
 729     printf("WARNING: result of tracks' removal will be only approximately correct");
 730
 731   TMatrixD rv(3,1);
 732   rv(0,0) = inVtx->GetXv();
 733   rv(1,0) = inVtx->GetYv();
 734   rv(2,0) = inVtx->GetZv();
 735   TMatrixD vV(3,3);
 736   Double_t cov[6];
 737   inVtx->GetCovMatrix(cov);
 738   vV(0,0) = cov[0];
 739   vV(0,1) = cov[1]; vV(1,0) = cov[1];
 740   vV(1,1) = cov[2];
 741   vV(0,2) = cov[3]; vV(2,0) = cov[3];
 742   vV(1,2) = cov[4]; vV(2,1) = cov[4];
 743   vV(2,2) = cov[5];
 744
 745   TMatrixD sumWi(TMatrixD::kInverted,vV);
 746   TMatrixD sumWiri(sumWi,TMatrixD::kMult,rv);
 747
 748   Int_t nUsedTrks = inVtx->GetNIndices();
 749   Double_t chi2 = inVtx->GetChi2();
 750
 751   AliExternalTrackParam *track = 0;
 752   Int_t ntrks = (Int_t)trkArray->GetEntriesFast();
 753   for(Int_t i=0;i<ntrks;i++) {
 754     track = (AliExternalTrackParam*)trkArray->At(i);
 755     if(!inVtx->UsesTrack(id[i])) {
 756       printf("track %d was not used in vertex fit",id[i]);
 757       continue;
 758     }
 759     Double_t alpha = track->GetAlpha();
 760     Double_t xl = diamondxy[0]*TMath::Cos(alpha)+diamondxy[1]*TMath::Sin(alpha);
 761     track->PropagateTo(xl,GetFieldkG());
 762     // vector of track global coordinates
 763     TMatrixD ri(3,1);
 764     // covariance matrix of ri
 765     TMatrixD wWi(3,3);
 766
 767     // get space point from track
 768     if(!TrackToPoint(track,ri,wWi)) continue;
 769
 770     TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
 771
 772     sumWi -= wWi;
 773     sumWiri -= wWiri;
 774
 775     // track contribution to chi2
 776     TMatrixD deltar = rv; deltar -= ri;
 777     TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
 778     Double_t chi2i = deltar(0,0)*wWideltar(0,0)+
 779                      deltar(1,0)*wWideltar(1,0)+
 780                      deltar(2,0)*wWideltar(2,0);
 781     // remove from total chi2
 782     chi2 -= chi2i;
 783
 784     nUsedTrks--;
 785     if(nUsedTrks<2) {
 786       printf("Trying to remove too many tracks!");
 787       return 0x0;
 788     }
 789   }
 790
 791   TMatrixD rvnew(3,1);
 792   TMatrixD vVnew(3,3);
 793
 794   // new inverted of weights matrix
 795   TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
 796   vVnew = invsumWi;
 797   // new position of primary vertex
 798   rvnew.Mult(vVnew,sumWiri);
 799
 800   Double_t position[3];
 801   position[0] = rvnew(0,0);
 802   position[1] = rvnew(1,0);
 803   position[2] = rvnew(2,0);
 804   cov[0] = vVnew(0,0);
 805   cov[1] = vVnew(0,1);
 806   cov[2] = vVnew(1,1);
 807   cov[3] = vVnew(0,2);
 808   cov[4] = vVnew(1,2);
 809   cov[5] = vVnew(2,2);
 810
 811   // store data in the vertex object
 812   AliESDVertex *outVtx = new AliESDVertex(position,cov,chi2,nUsedTrks+1); // the +1 is for the constraint
 813   outVtx->SetTitle(inVtx->GetTitle());
 814   UShort_t *inindices = inVtx->GetIndices();
 815   Int_t nIndices = nUsedTrks;
 816   UShort_t *outindices = new UShort_t[nIndices];
 817   Int_t j=0;
 818   for(Int_t k=0; k<inVtx->GetNIndices(); k++) {
 819     Bool_t copyindex=kTRUE;
 820     for(Int_t l=0; l<ntrks; l++) {
 821       if(inindices[k]==id[l]) copyindex=kFALSE;
 822     }
 823     if(copyindex) {
 824       outindices[j] = inindices[k]; j++;
 825     }
 826   }
 827   outVtx->SetIndices(nIndices,outindices);
 828   if (outindices) delete outindices;
 829
 830   /*
 831     printf("Vertex before removing tracks:");
 832     inVtx->PrintStatus();
 833     inVtx->PrintIndices();
 834     printf("Vertex after removing tracks:");
 835     outVtx->PrintStatus();
 836     outVtx->PrintIndices();
 837   */
 838
 839   return outVtx;
 840 }
 841 //---------------------------------------------------------------------------
 842 AliESDVertex* AliVertexerTracks::RemoveConstraintFromVertex(AliESDVertex *inVtx,
 843                                                      Float_t *diamondxyz,
 844                                                      Float_t *diamondcov) const
 845 {
 846 //
 847 // Removes diamond constraint from fit of inVtx
 848 //
 849
 850   if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraint")) {
 851     printf("ERROR: primary vertex has no constraint: cannot remove it\n");
 852     return 0x0;
 853   }
 854   if(inVtx->GetNContributors()<3) {
 855     printf("ERROR: primary vertex has less than 2 tracks: cannot remove contraint\n");
 856     return 0x0;
 857   }
 858
 859   // diamond constraint
 860   TMatrixD vVb(3,3);
 861   vVb(0,0) = diamondcov[0];
 862   vVb(0,1) = diamondcov[1];
 863   vVb(0,2) = 0.;
 864   vVb(1,0) = diamondcov[1];
 865   vVb(1,1) = diamondcov[2];
 866   vVb(1,2) = 0.;
 867   vVb(2,0) = 0.;
 868   vVb(2,1) = 0.;
 869   vVb(2,2) = diamondcov[5];
 870   TMatrixD vVbInv(TMatrixD::kInverted,vVb);
 871   TMatrixD rb(3,1);
 872   rb(0,0) = diamondxyz[0];
 873   rb(1,0) = diamondxyz[1];
 874   rb(2,0) = diamondxyz[2];
 875   TMatrixD vVbInvrb(vVbInv,TMatrixD::kMult,rb);
 876
 877   // input vertex
 878   TMatrixD rv(3,1);
 879   rv(0,0) = inVtx->GetXv();
 880   rv(1,0) = inVtx->GetYv();
 881   rv(2,0) = inVtx->GetZv();
 882   TMatrixD vV(3,3);
 883   Double_t cov[6];
 884   inVtx->GetCovMatrix(cov);
 885   vV(0,0) = cov[0];
 886   vV(0,1) = cov[1]; vV(1,0) = cov[1];
 887   vV(1,1) = cov[2];
 888   vV(0,2) = cov[3]; vV(2,0) = cov[3];
 889   vV(1,2) = cov[4]; vV(2,1) = cov[4];
 890   vV(2,2) = cov[5];
 891   TMatrixD vVInv(TMatrixD::kInverted,vV);
 892   TMatrixD vVInvrv(vVInv,TMatrixD::kMult,rv);
 893
 894
 895   TMatrixD sumWi = vVInv - vVbInv;
 896
 897
 898   TMatrixD sumWiri = vVInvrv - vVbInvrb;
 899
 900   TMatrixD rvnew(3,1);
 901   TMatrixD vVnew(3,3);
 902
 903   // new inverted of weights matrix
 904   TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
 905   vVnew = invsumWi;
 906   // new position of primary vertex
 907   rvnew.Mult(vVnew,sumWiri);
 908
 909   Double_t position[3];
 910   position[0] = rvnew(0,0);
 911   position[1] = rvnew(1,0);
 912   position[2] = rvnew(2,0);
 913   cov[0] = vVnew(0,0);
 914   cov[1] = vVnew(0,1);
 915   cov[2] = vVnew(1,1);
 916   cov[3] = vVnew(0,2);
 917   cov[4] = vVnew(1,2);
 918   cov[5] = vVnew(2,2);
 919
 920
 921   Double_t chi2 = inVtx->GetChi2();
 922
 923   // diamond constribution to chi2
 924   TMatrixD deltar = rv; deltar -= rb;
 925   TMatrixD vVbInvdeltar(vVbInv,TMatrixD::kMult,deltar);
 926   Double_t chi2b = deltar(0,0)*vVbInvdeltar(0,0)+
 927                    deltar(1,0)*vVbInvdeltar(1,0)+
 928                    deltar(2,0)*vVbInvdeltar(2,0);
 929   // remove from total chi2
 930   chi2 -= chi2b;
 931
 932   // store data in the vertex object
 933   AliESDVertex *outVtx = new AliESDVertex(position,cov,chi2,inVtx->GetNContributors()-1);
 934   outVtx->SetTitle("VertexerTracksNoConstraint");
 935   UShort_t *inindices = inVtx->GetIndices();
 936   Int_t nIndices = inVtx->GetNIndices();
 937   outVtx->SetIndices(nIndices,inindices);
 938
 939   return outVtx;
 940 }
 941 //---------------------------------------------------------------------------
 942 void AliVertexerTracks::SetCuts(Double_t *cuts)
 943 {
 944 //
 945 //  Cut values
 946 //
 947   SetDCAcut(cuts[0]);
 948   SetDCAcutIter0(cuts[1]);
 949   SetMaxd0z0(cuts[2]);
 950   if(fMode==0 && cuts[3]<0) SetITSrefitNotRequired();
 951   SetMinClusters((Int_t)(TMath::Abs(cuts[3])));
 952   SetMinTracks((Int_t)(cuts[4]));
 953   SetNSigmad0(cuts[5]);
 954   SetMinDetFitter(cuts[6]);
 955   SetMaxTgl(cuts[7]);
 956   SetFiducialRZ(cuts[8],cuts[9]);
 957   fAlgo=(Int_t)(cuts[10]);
 958   fAlgoIter0=(Int_t)(cuts[11]);
 959
 960   return;
 961 }
 962 //---------------------------------------------------------------------------
 963 void AliVertexerTracks::SetITSMode(Double_t dcacut,
 964                                    Double_t dcacutIter0,
 965                                    Double_t maxd0z0,
 966                                    Int_t minCls,
 967                                    Int_t mintrks,
 968                                    Double_t nsigma,
 969                                    Double_t mindetfitter,
 970                                    Double_t maxtgl,
 971                                    Double_t fidR,
 972                                    Double_t fidZ,
 973                                    Int_t finderAlgo,
 974                                    Int_t finderAlgoIter0)
 975 {
 976 //
 977 //  Cut values for ITS mode
 978 //
 979   fMode = 0;
 980   if(minCls>0) {
 981     SetITSrefitRequired();
 982   } else {
 983     SetITSrefitNotRequired();
 984   }
 985   SetDCAcut(dcacut);
 986   SetDCAcutIter0(dcacutIter0);
 987   SetMaxd0z0(maxd0z0);
 988   SetMinClusters(TMath::Abs(minCls));
 989   SetMinTracks(mintrks);
 990   SetNSigmad0(nsigma);
 991   SetMinDetFitter(mindetfitter);
 992   SetMaxTgl(maxtgl);
 993   SetFiducialRZ(fidR,fidZ);
 994   fAlgo=finderAlgo;
 995   fAlgoIter0=finderAlgoIter0;
 996
 997   return;
 998 }
 999 //---------------------------------------------------------------------------
1000 void AliVertexerTracks::SetTPCMode(Double_t dcacut,
1001                                    Double_t dcacutIter0,
1002                                    Double_t maxd0z0,
1003                                    Int_t minCls,
1004                                    Int_t mintrks,
1005                                    Double_t nsigma,
1006                                    Double_t mindetfitter,
1007                                    Double_t maxtgl,
1008                                    Double_t fidR,
1009                                    Double_t fidZ,
1010                                    Int_t finderAlgo,
1011                                    Int_t finderAlgoIter0)
1012 {
1013 //
1014 //  Cut values for TPC mode
1015 //
1016   fMode = 1;
1017   SetITSrefitNotRequired();
1018   SetDCAcut(dcacut);
1019   SetDCAcutIter0(dcacutIter0);
1020   SetMaxd0z0(maxd0z0);
1021   SetMinClusters(minCls);
1022   SetMinTracks(mintrks);
1023   SetNSigmad0(nsigma);
1024   SetMinDetFitter(mindetfitter);
1025   SetMaxTgl(maxtgl);
1026   SetFiducialRZ(fidR,fidZ);
1027   fAlgo=finderAlgo;
1028   fAlgoIter0=finderAlgoIter0;
1029
1030   return;
1031 }
1032 //---------------------------------------------------------------------------
1033 void AliVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped)
1034 {
1035 //
1036 // Mark the tracks not to be used in the vertex reconstruction.
1037 // Tracks are identified by AliESDtrack::GetID()
1038 //
1039   fNTrksToSkip = n;  fTrksToSkip = new Int_t[n];
1040   for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i];
1041   return;
1042 }
1043 //---------------------------------------------------------------------------
1044 void  AliVertexerTracks::SetVtxStart(AliESDVertex *vtx)
1045 {
1046 //
1047 // Set initial vertex knowledge
1048 //
1049   vtx->GetXYZ(fNominalPos);
1050   vtx->GetCovMatrix(fNominalCov);
1051   SetConstraintOn();
1052   return;
1053 }
1054 //---------------------------------------------------------------------------
1055 void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
1056 {
1057   AliExternalTrackParam *track1;
1058   const Int_t knacc = (Int_t)fTrkArraySel.GetEntriesFast();
1059   AliStrLine **linarray = new AliStrLine* [knacc];
1060   for(Int_t i=0; i<knacc; i++){
1061     track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
1062     Double_t alpha=track1->GetAlpha();
1063     Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
1064     Double_t pos[3],dir[3],sigmasq[3];
1065     track1->GetXYZAt(mindist,GetFieldkG(),pos);
1066     track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
1067     sigmasq[0]=TMath::Sin(alpha)*TMath::Sin(alpha)*track1->GetSigmaY2();
1068     sigmasq[1]=TMath::Cos(alpha)*TMath::Cos(alpha)*track1->GetSigmaY2();
1069     sigmasq[2]=track1->GetSigmaZ2();
1070     TMatrixD ri(3,1);
1071     TMatrixD wWi(3,3);
1072     if(!TrackToPoint(track1,ri,wWi)) {optUseWeights=kFALSE;printf("WARNING\n");}
1073     Double_t wmat[9];
1074     Int_t iel=0;
1075     for(Int_t ia=0;ia<3;ia++){
1076       for(Int_t ib=0;ib<3;ib++){
1077         wmat[iel]=wWi(ia,ib);
1078         iel++;
1079       }
1080     }
1081     linarray[i] = new AliStrLine(pos,sigmasq,wmat,dir);
1082   }
1083   fVert=TrackletVertexFinder(linarray,knacc,optUseWeights);
1084   for(Int_t i=0; i<knacc; i++) delete linarray[i];
1085   delete [] linarray;
1086 }
1087 //---------------------------------------------------------------------------
1088 AliESDVertex AliVertexerTracks::TrackletVertexFinder(TClonesArray *lines, Int_t optUseWeights)
1089 {
1090   // Calculate the point at minimum distance to prepared tracks (TClonesArray)
1091   const Int_t knacc = (Int_t)lines->GetEntriesFast();
1092   AliStrLine** lines2 = new AliStrLine* [knacc];
1093   for(Int_t i=0; i<knacc; i++){
1094     lines2[i]= (AliStrLine*)lines->At(i);
1095   }
1096   AliESDVertex vert = TrackletVertexFinder(lines2,knacc,optUseWeights);
1097   delete [] lines2;
1098   return vert;
1099 }
1100
1101 //---------------------------------------------------------------------------
1102 AliESDVertex AliVertexerTracks::TrackletVertexFinder(AliStrLine **lines, const Int_t knacc, Int_t optUseWeights)
1103 {
1104   // Calculate the point at minimum distance to prepared tracks (array of AliStrLine)
1105
1106   Double_t initPos[3]={0.,0.,0.};
1107
1108   Double_t (*vectP0)[3]=new Double_t [knacc][3];
1109   Double_t (*vectP1)[3]=new Double_t [knacc][3];
1110
1111   Double_t sum[3][3];
1112   Double_t dsum[3]={0,0,0};
1113   TMatrixD sumWi(3,3);
1114   for(Int_t i=0;i<3;i++){
1115     for(Int_t j=0;j<3;j++){
1116       sum[i][j]=0;
1117       sumWi(i,j)=0.;
1118     }
1119   }
1120
1121   for(Int_t i=0; i<knacc; i++){
1122     AliStrLine *line1 = lines[i];
1123     Double_t p0[3],cd[3],sigmasq[3];
1124     Double_t wmat[9];
1125     if(!line1) printf("ERROR %d %d\n",i,knacc);
1126     line1->GetP0(p0);
1127     line1->GetCd(cd);
1128     line1->GetSigma2P0(sigmasq);
1129     line1->GetWMatrix(wmat);
1130     TMatrixD wWi(3,3);
1131     Int_t iel=0;
1132     for(Int_t ia=0;ia<3;ia++){
1133       for(Int_t ib=0;ib<3;ib++){
1134         wWi(ia,ib)=wmat[iel];
1135         iel++;
1136       }
1137     }
1138
1139     sumWi+=wWi;
1140
1141     Double_t p1[3]={p0[0]+cd[0],p0[1]+cd[1],p0[2]+cd[2]};
1142     vectP0[i][0]=p0[0];
1143     vectP0[i][1]=p0[1];
1144     vectP0[i][2]=p0[2];
1145     vectP1[i][0]=p1[0];
1146     vectP1[i][1]=p1[1];
1147     vectP1[i][2]=p1[2];
1148
1149     Double_t matr[3][3];
1150     Double_t dknow[3];
1151     if(optUseWeights==0)GetStrLinDerivMatrix(p0,p1,matr,dknow);
1152     else GetStrLinDerivMatrix(p0,p1,sigmasq,matr,dknow);
1153
1154
1155     for(Int_t iii=0;iii<3;iii++){
1156       dsum[iii]+=dknow[iii];
1157       for(Int_t lj=0;lj<3;lj++) sum[iii][lj]+=matr[iii][lj];
1158     }
1159   }
1160
1161   TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
1162   Double_t covmatrix[6];
1163   covmatrix[0] = invsumWi(0,0);
1164   covmatrix[1] = invsumWi(0,1);
1165   covmatrix[2] = invsumWi(1,1);
1166   covmatrix[3] = invsumWi(0,2);
1167   covmatrix[4] = invsumWi(1,2);
1168   covmatrix[5] = invsumWi(2,2);
1169
1170   Double_t vett[3][3];
1171   Double_t det=GetDeterminant3X3(sum);
1172   Double_t sigma=0;
1173
1174   if(det!=0){
1175     for(Int_t zz=0;zz<3;zz++){
1176       for(Int_t ww=0;ww<3;ww++){
1177         for(Int_t kk=0;kk<3;kk++) vett[ww][kk]=sum[ww][kk];
1178       }
1179       for(Int_t kk=0;kk<3;kk++) vett[kk][zz]=dsum[kk];
1180       initPos[zz]=GetDeterminant3X3(vett)/det;
1181     }
1182
1183
1184     for(Int_t i=0; i<knacc; i++){
1185       Double_t p0[3]={0,0,0},p1[3]={0,0,0};
1186       for(Int_t ii=0;ii<3;ii++){
1187         p0[ii]=vectP0[i][ii];
1188         p1[ii]=vectP1[i][ii];
1189       }
1190       sigma+=GetStrLinMinDist(p0,p1,initPos);
1191     }
1192
1193     if(sigma>0.) {sigma=TMath::Sqrt(sigma);}else{sigma=999;}
1194   }else{
1195     sigma=999;
1196   }
1197   AliESDVertex theVert(initPos,covmatrix,99999.,knacc);
1198   theVert.SetDispersion(sigma);
1199   delete vectP0;
1200   delete vectP1;
1201   return theVert;
1202 }
1203
1204 //---------------------------------------------------------------------------
1205 Bool_t AliVertexerTracks::TrackToPoint(AliExternalTrackParam *t,
1206                                        TMatrixD &ri,TMatrixD &wWi,
1207                                        Bool_t uUi3by3) const
1208 {
1209 //
1210 // Extract from the AliExternalTrackParam the global coordinates ri and covariance matrix
1211 // wWi of the space point that it represents (to be used in VertexFitter())
1212 //
1213
1214
1215   Double_t rotAngle = t->GetAlpha();
1216   if(rotAngle<0.) rotAngle += 2.*TMath::Pi();
1217   Double_t cosRot = TMath::Cos(rotAngle);
1218   Double_t sinRot = TMath::Sin(rotAngle);
1219
1220   ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot;
1221   ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot;
1222   ri(2,0) = t->GetZ();
1223
1224
1225
1226   if(!uUi3by3) {
1227     // matrix to go from global (x,y,z) to local (y,z);
1228     TMatrixD qQi(2,3);
1229     qQi(0,0) = -sinRot;
1230     qQi(0,1) = cosRot;
1231     qQi(0,2) = 0.;
1232     qQi(1,0) = 0.;
1233     qQi(1,1) = 0.;
1234     qQi(1,2) = 1.;
1235
1236     // covariance matrix of local (y,z) - inverted
1237     TMatrixD uUi(2,2);
1238     uUi(0,0) = t->GetSigmaY2();
1239     uUi(0,1) = t->GetSigmaZY();
1240     uUi(1,0) = t->GetSigmaZY();
1241     uUi(1,1) = t->GetSigmaZ2();
1242     //printf(" Ui :");
1243     //printf(" %f   %f",uUi(0,0),uUi(0,1));
1244     //printf(" %f   %f",uUi(1,0),uUi(1,1));
1245
1246     if(uUi.Determinant() <= 0.) return kFALSE;
1247     TMatrixD uUiInv(TMatrixD::kInverted,uUi);
1248
1249     // weights matrix: wWi = qQiT * uUiInv * qQi
1250     TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
1251     TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
1252     wWi = m;
1253   } else {
1254     if(fVert.GetNContributors()<1) AliFatal("Vertex from finder is empty");
1255     // matrix to go from global (x,y,z) to local (x,y,z);
1256     TMatrixD qQi(3,3);
1257     qQi(0,0) = cosRot;
1258     qQi(0,1) = sinRot;
1259     qQi(0,2) = 0.;
1260     qQi(1,0) = -sinRot;
1261     qQi(1,1) = cosRot;
1262     qQi(1,2) = 0.;
1263     qQi(2,0) = 0.;
1264     qQi(2,1) = 0.;
1265     qQi(2,2) = 1.;
1266
1267     // covariance of fVert along the track
1268     Double_t p[3],pt,ptot;
1269     t->GetPxPyPz(p);
1270     pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
1271     ptot = TMath::Sqrt(pt*pt+p[2]*p[2]);
1272     Double_t cphi = p[0]/pt;               //cos(phi)=px/pt
1273     Double_t sphi = p[1]/pt;               //sin(phi)=py/pt
1274     Double_t clambda = pt/ptot;            //cos(lambda)=pt/ptot
1275     Double_t slambda = p[2]/ptot;            //sin(lambda)=pz/ptot
1276     Double_t covfVert[6];
1277     fVert.GetCovMatrix(covfVert);
1278     Double_t covfVertalongt =
1279        covfVert[0]*cphi*cphi*clambda*clambda
1280       +covfVert[1]*2.*cphi*sphi*clambda*clambda
1281       +covfVert[3]*2.*cphi*clambda*slambda
1282       +covfVert[2]*sphi*sphi*clambda*clambda
1283       +covfVert[4]*2.*sphi*clambda*slambda
1284       +covfVert[5]*slambda*slambda;
1285     // covariance matrix of local (x,y,z) - inverted
1286     TMatrixD uUi(3,3);
1287     uUi(0,0) = covfVertalongt * fnSigmaForUi00 * fnSigmaForUi00;
1288     AliDebug(1,Form("=====> sqrtUi00 cm  %f",TMath::Sqrt(uUi(0,0))));
1289     uUi(0,1) = 0.;
1290     uUi(0,2) = 0.;
1291     uUi(1,0) = 0.;
1292     uUi(1,1) = t->GetSigmaY2();
1293     uUi(1,2) = t->GetSigmaZY();
1294     uUi(2,0) = 0.;
1295     uUi(2,1) = t->GetSigmaZY();
1296     uUi(2,2) = t->GetSigmaZ2();
1297     //printf(" Ui :\n");
1298     //printf(" %f   %f\n",uUi(0,0),uUi(0,1));
1299     //printf(" %f   %f\n",uUi(1,0),uUi(1,1));
1300
1301     if(uUi.Determinant() <= 0.) return kFALSE;
1302     TMatrixD uUiInv(TMatrixD::kInverted,uUi);
1303
1304     // weights matrix: wWi = qQiT * uUiInv * qQi
1305     TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
1306     TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
1307     wWi = m;
1308   }
1309
1310
1311   return kTRUE;
1312 }
1313 //---------------------------------------------------------------------------
1314 void AliVertexerTracks::TooFewTracks()
1315 {
1316 //
1317 // When the number of tracks is < fMinTracks,
1318 // deal with vertices not found and prepare to exit
1319 //
1320   AliDebug(1,"TooFewTracks");
1321
1322   Double_t pos[3],err[3];
1323   pos[0] = fNominalPos[0];
1324   err[0] = TMath::Sqrt(fNominalCov[0]);
1325   pos[1] = fNominalPos[1];
1326   err[1] = TMath::Sqrt(fNominalCov[2]);
1327   pos[2] = fNominalPos[2];
1328   err[2] = TMath::Sqrt(fNominalCov[5]);
1329   Int_t    ncontr = (err[0]>1. ? -1 : -3);
1330   if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
1331   fCurrentVertex = new AliESDVertex(pos,err);
1332   fCurrentVertex->SetNContributors(ncontr);
1333
1334   if(fConstraint) {
1335     fCurrentVertex->SetTitle("VertexerTracksWithConstraint");
1336   } else {
1337     fCurrentVertex->SetTitle("VertexerTracksNoConstraint");
1338   }
1339
1340   if(!fTrkArraySel.IsEmpty()) fTrkArraySel.Delete();
1341   if(fIdSel) {delete fIdSel; fIdSel=NULL;}
1342   if(fTrksToSkip) {delete fTrksToSkip; fTrksToSkip=NULL;}
1343
1344   return;
1345 }
1346 //---------------------------------------------------------------------------
1347 void AliVertexerTracks::VertexFinder(Int_t optUseWeights)
1348 {
1349
1350   // Get estimate of vertex position in (x,y) from tracks DCA
1351
1352   Double_t initPos[3];
1353   initPos[2] = 0.;
1354   for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
1355   Int_t nacc = (Int_t)fTrkArraySel.GetEntriesFast();
1356   Double_t aver[3]={0.,0.,0.};
1357   Double_t aversq[3]={0.,0.,0.};
1358   Double_t sigmasq[3]={0.,0.,0.};
1359   Double_t sigma=0;
1360   Int_t ncombi = 0;
1361   AliExternalTrackParam *track1;
1362   AliExternalTrackParam *track2;
1363   Double_t pos[3],dir[3];
1364   Double_t alpha,mindist;
1365
1366   for(Int_t i=0; i<nacc; i++){
1367     track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
1368     alpha=track1->GetAlpha();
1369     mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
1370     track1->GetXYZAt(mindist,GetFieldkG(),pos);
1371     track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
1372     AliStrLine *line1 = new AliStrLine(pos,dir);
1373
1374    //    AliStrLine *line1 = new AliStrLine();
1375    //    track1->ApproximateHelixWithLine(mindist,GetFieldkG(),line1);
1376
1377     for(Int_t j=i+1; j<nacc; j++){
1378       track2 = (AliExternalTrackParam*)fTrkArraySel.At(j);
1379       alpha=track2->GetAlpha();
1380       mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
1381       track2->GetXYZAt(mindist,GetFieldkG(),pos);
1382       track2->GetPxPyPzAt(mindist,GetFieldkG(),dir);
1383       AliStrLine *line2 = new AliStrLine(pos,dir);
1384     //      AliStrLine *line2 = new AliStrLine();
1385     //  track2->ApproximateHelixWithLine(mindist,GetFieldkG(),line2);
1386       Double_t distCA=line2->GetDCA(line1);
1387       //printf("%d   %d   %f\n",i,j,distCA);
1388        if(fDCAcut<=0 || (fDCAcut>0&&distCA<fDCAcut)){
1389         Double_t pnt1[3],pnt2[3],crosspoint[3];
1390
1391         if(optUseWeights<=0){
1392           Int_t retcode = line2->Cross(line1,crosspoint);
1393           if(retcode>=0){
1394             ncombi++;
1395             for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
1396             for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
1397           }
1398         }
1399         if(optUseWeights>0){
1400           Int_t retcode = line1->CrossPoints(line2,pnt1,pnt2);
1401           if(retcode>=0){
1402             Double_t cs, sn;
1403             alpha=track1->GetAlpha();
1404             cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
1405             Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
1406             alpha=track2->GetAlpha();
1407             cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
1408             Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
1409             Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
1410             Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
1411             Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
1412             Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
1413             crosspoint[0]=wx1*pnt1[0] + wx2*pnt2[0];
1414             crosspoint[1]=wy1*pnt1[1] + wy2*pnt2[1];
1415             crosspoint[2]=wz1*pnt1[2] + wz2*pnt2[2];
1416
1417             ncombi++;
1418             for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
1419             for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
1420           }
1421         }
1422       }
1423       delete line2;
1424     }
1425     delete line1;
1426   }
1427   if(ncombi>0){
1428     for(Int_t jj=0;jj<3;jj++){
1429       initPos[jj] = aver[jj]/ncombi;
1430       //printf("%f\n",initPos[jj]);
1431       aversq[jj]/=ncombi;
1432       sigmasq[jj]=aversq[jj]-initPos[jj]*initPos[jj];
1433       sigma+=sigmasq[jj];
1434     }
1435     sigma=TMath::Sqrt(TMath::Abs(sigma));
1436   }
1437   else {
1438     Warning("VertexFinder","Finder did not succed");
1439     sigma=999;
1440   }
1441   fVert.SetXYZ(initPos);
1442   fVert.SetDispersion(sigma);
1443   fVert.SetNContributors(ncombi);
1444 }
1445 //---------------------------------------------------------------------------
1446 void AliVertexerTracks::VertexFitter()
1447 {
1448 //
1449 // The optimal estimate of the vertex position is given by a "weighted
1450 // average of tracks positions".
1451 // Original method: V. Karimaki, CMS Note 97/0051
1452 //
1453   Bool_t useConstraint = fConstraint;
1454   Double_t initPos[3];
1455   if(!fOnlyFitter) {
1456     fVert.GetXYZ(initPos);
1457   } else {
1458     initPos[0]=fNominalPos[0];
1459     initPos[1]=fNominalPos[1];
1460     initPos[2]=fNominalPos[2];
1461   }
1462
1463   Int_t nTrksSel = (Int_t)fTrkArraySel.GetEntries();
1464   if(nTrksSel==1) useConstraint=kTRUE;
1465   AliDebug(1,"--- VertexFitter(): start");
1466   AliDebug(1,Form(" Number of tracks in array: %d\n",nTrksSel));
1467   AliDebug(1,Form(" Minimum # tracks required in fit: %d\n",fMinTracks));
1468   AliDebug(1,Form(" Vertex position after finder: %f,%f,%f\n",initPos[0],initPos[1],initPos[2]));
1469   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])));
1470
1471   // special treatment for few-tracks fits (e.g. secondary vertices)
1472   Bool_t uUi3by3 = kFALSE; if(nTrksSel<5 && !useConstraint) uUi3by3 = kTRUE;
1473
1474   Int_t i,j,k,step=0;
1475   TMatrixD rv(3,1);
1476   TMatrixD vV(3,3);
1477   rv(0,0) = initPos[0];
1478   rv(1,0) = initPos[1];
1479   rv(2,0) = 0.;
1480   Double_t xlStart,alpha;
1481   Int_t nTrksUsed;
1482   Double_t chi2,chi2i,chi2b;
1483   AliExternalTrackParam *t = 0;
1484   Int_t failed = 0;
1485
1486   // initial vertex covariance matrix
1487   TMatrixD vVb(3,3);
1488   vVb(0,0) = fNominalCov[0];
1489   vVb(0,1) = fNominalCov[1];
1490   vVb(0,2) = 0.;
1491   vVb(1,0) = fNominalCov[1];
1492   vVb(1,1) = fNominalCov[2];
1493   vVb(1,2) = 0.;
1494   vVb(2,0) = 0.;
1495   vVb(2,1) = 0.;
1496   vVb(2,2) = fNominalCov[5];
1497   TMatrixD vVbInv(TMatrixD::kInverted,vVb);
1498   TMatrixD rb(3,1);
1499   rb(0,0) = fNominalPos[0];
1500   rb(1,0) = fNominalPos[1];
1501   rb(2,0) = fNominalPos[2];
1502   TMatrixD vVbInvrb(vVbInv,TMatrixD::kMult,rb);
1503
1504
1505   // 2 steps:
1506   // 1st - estimate of vtx using all tracks
1507   // 2nd - estimate of global chi2
1508   for(step=0; step<2; step++) {
1509     AliDebug(1,Form(" step = %d\n",step));
1510     chi2 = 0.;
1511     nTrksUsed = 0;
1512
1513     if(step==1) { initPos[0]=rv(0,0); initPos[0]=rv(1,0); }
1514
1515     TMatrixD sumWiri(3,1);
1516     TMatrixD sumWi(3,3);
1517     for(i=0; i<3; i++) {
1518       sumWiri(i,0) = 0.;
1519       for(j=0; j<3; j++) sumWi(j,i) = 0.;
1520     }
1521
1522     // mean vertex constraint
1523     if(useConstraint) {
1524       for(i=0;i<3;i++) {
1525         sumWiri(i,0) += vVbInvrb(i,0);
1526         for(k=0;k<3;k++) sumWi(i,k) += vVbInv(i,k);
1527       }
1528       // chi2
1529       TMatrixD deltar = rv; deltar -= rb;
1530       TMatrixD vVbInvdeltar(vVbInv,TMatrixD::kMult,deltar);
1531       chi2b = deltar(0,0)*vVbInvdeltar(0,0)+
1532               deltar(1,0)*vVbInvdeltar(1,0)+
1533               deltar(2,0)*vVbInvdeltar(2,0);
1534       chi2 += chi2b;
1535     }
1536
1537     // loop on tracks
1538     for(k=0; k<nTrksSel; k++) {
1539       // get track from track array
1540       t = (AliExternalTrackParam*)fTrkArraySel.At(k);
1541       alpha = t->GetAlpha();
1542       xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha);
1543       // to vtxSeed (from finder)
1544       t->PropagateTo(xlStart,GetFieldkG());
1545
1546       // vector of track global coordinates
1547       TMatrixD ri(3,1);
1548       // covariance matrix of ri
1549       TMatrixD wWi(3,3);
1550
1551       // get space point from track
1552       if(!TrackToPoint(t,ri,wWi,uUi3by3)) continue;
1553       TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
1554
1555       // track chi2
1556       TMatrixD deltar = rv; deltar -= ri;
1557       TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
1558       chi2i = deltar(0,0)*wWideltar(0,0)+
1559               deltar(1,0)*wWideltar(1,0)+
1560               deltar(2,0)*wWideltar(2,0);
1561
1562       // add to total chi2
1563       chi2 += chi2i;
1564
1565       sumWiri += wWiri;
1566       sumWi   += wWi;
1567
1568       nTrksUsed++;
1569     } // end loop on tracks
1570
1571     if(nTrksUsed < fMinTracks) {
1572       failed=1;
1573       continue;
1574     }
1575
1576     Double_t determinant = sumWi.Determinant();
1577     if(determinant < fMinDetFitter)  {
1578       AliDebug(1,Form("det(V) = %f (<%f)\n",determinant,fMinDetFitter));
1579       failed=1;
1580       continue;
1581     }
1582
1583     if(step==0) {
1584       // inverted of weights matrix
1585       TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
1586       vV = invsumWi;
1587       // position of primary vertex
1588       rv.Mult(vV,sumWiri);
1589     }
1590   } // end loop on the 2 steps
1591
1592   if(failed) {
1593     TooFewTracks();
1594     return;
1595   }
1596
1597   Double_t position[3];
1598   position[0] = rv(0,0);
1599   position[1] = rv(1,0);
1600   position[2] = rv(2,0);
1601   Double_t covmatrix[6];
1602   covmatrix[0] = vV(0,0);
1603   covmatrix[1] = vV(0,1);
1604   covmatrix[2] = vV(1,1);
1605   covmatrix[3] = vV(0,2);
1606   covmatrix[4] = vV(1,2);
1607   covmatrix[5] = vV(2,2);
1608
1609   // for correct chi2/ndf, count constraint as additional "track"
1610   if(fConstraint) nTrksUsed++;
1611   // store data in the vertex object
1612   if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
1613   fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);
1614
1615   AliDebug(1," Vertex after fit:");
1616   AliDebug(1,Form("xyz: %f %f %f; nc %d",fCurrentVertex->GetXv(),fCurrentVertex->GetYv(),fCurrentVertex->GetZv(),fCurrentVertex->GetNContributors()));
1617   AliDebug(1,"--- VertexFitter(): finish\n");
1618
1619
1620   return;
1621 }
1622 //----------------------------------------------------------------------------
1623 AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray,
1624                                                          UShort_t *id,
1625                                                          Bool_t optUseFitter,
1626                                                          Bool_t optPropagate)
1627 {
1628 //
1629 // Return vertex from tracks (AliExternalTrackParam) in array
1630 //
1631   fCurrentVertex = 0;
1632   SetConstraintOff();
1633
1634   // get tracks and propagate them to initial vertex position
1635   fIdSel = new UShort_t[(Int_t)trkArray->GetEntriesFast()];
1636   Int_t nTrksSel = PrepareTracks(*trkArray,id,0);
1637   if(nTrksSel <  TMath::Max(2,fMinTracks)) {
1638     TooFewTracks();
1639     return fCurrentVertex;
1640   }
1641
1642   // vertex finder
1643   switch (fAlgo) {
1644     case 1: StrLinVertexFinderMinDist(1); break;
1645     case 2: StrLinVertexFinderMinDist(0); break;
1646     case 3: HelixVertexFinder();          break;
1647     case 4: VertexFinder(1);              break;
1648     case 5: VertexFinder(0);              break;
1649     default: printf("Wrong algorithm\n"); break;
1650   }
1651   AliDebug(1," Vertex finding completed\n");
1652
1653   // vertex fitter
1654   if(optUseFitter) {
1655     VertexFitter();
1656   } else {
1657     Double_t position[3]={fVert.GetXv(),fVert.GetYv(),fVert.GetZv()};
1658     Double_t covmatrix[6];
1659     fVert.GetCovMatrix(covmatrix);
1660     Double_t chi2=99999.;
1661     Int_t    nTrksUsed=fVert.GetNContributors();
1662     fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);
1663   }
1664   fCurrentVertex->SetDispersion(fVert.GetDispersion());
1665
1666
1667   // set indices of used tracks and propagate track to found vertex
1668   UShort_t *indices = 0;
1669   Double_t d0z0[2],covd0z0[3];
1670   AliExternalTrackParam *t = 0;
1671   if(fCurrentVertex->GetNContributors()>0) {
1672     indices = new UShort_t[fTrkArraySel.GetEntriesFast()];
1673     for(Int_t jj=0; jj<(Int_t)fTrkArraySel.GetEntriesFast(); jj++) {
1674       indices[jj] = fIdSel[jj];
1675       t = (AliExternalTrackParam*)fTrkArraySel.At(jj);
1676       if(optPropagate && optUseFitter) {
1677         if(TMath::Sqrt(fCurrentVertex->GetXv()*fCurrentVertex->GetXv()+fCurrentVertex->GetYv()*fCurrentVertex->GetYv())<3.) {
1678           t->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
1679           AliDebug(1,Form("Track %d propagated to found vertex",jj));
1680         } else {
1681           AliWarning("Found vertex outside beam pipe!");
1682         }
1683       }
1684     }
1685     fCurrentVertex->SetIndices(fCurrentVertex->GetNContributors(),indices);
1686   }
1687
1688   // clean up
1689   if (indices) {delete indices; indices=NULL;}
1690   delete fIdSel; fIdSel=NULL;
1691   fTrkArraySel.Delete();
1692
1693   return fCurrentVertex;
1694 }
1695 //----------------------------------------------------------------------------
1696 AliESDVertex* AliVertexerTracks::VertexForSelectedESDTracks(TObjArray *trkArray,
1697                                                          Bool_t optUseFitter,
1698                                                          Bool_t optPropagate)
1699 {
1700 //
1701 // Return vertex from array of ESD tracks
1702 //
1703
1704   Int_t nTrks = (Int_t)trkArray->GetEntriesFast();
1705   UShort_t *id = new UShort_t[nTrks];
1706
1707   AliESDtrack *esdt = 0;
1708   for(Int_t i=0; i<nTrks; i++){
1709     esdt = (AliESDtrack*)trkArray->At(i);
1710     id[i] = (UShort_t)esdt->GetID();
1711   }
1712
1713   VertexForSelectedTracks(trkArray,id,optUseFitter,optPropagate);
1714
1715   delete id; id=NULL;
1716
1717   return fCurrentVertex;
1718 }
1719 //--------------------------------------------------------------------------