- Three classes by MinJung Kweon AliHFEpriVtx, AliHFEsecVtx and AliHFEmcQA for primar...
[u/mrichter/AliRoot.git] / 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) 
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 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);
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 void AliVertexerTracks::SetCuts(Double_t *cuts) 
843 {
844 //
845 //  Cut values
846 //
847   SetDCAcut(cuts[0]);
848   SetDCAcutIter0(cuts[1]);
849   SetMaxd0z0(cuts[2]);
850   if(fMode==0 && cuts[3]<0) SetITSrefitNotRequired();
851   SetMinClusters((Int_t)(TMath::Abs(cuts[3])));
852   SetMinTracks((Int_t)(cuts[4]));
853   SetNSigmad0(cuts[5]);
854   SetMinDetFitter(cuts[6]);
855   SetMaxTgl(cuts[7]);
856   SetFiducialRZ(cuts[8],cuts[9]);
857   fAlgo=(Int_t)(cuts[10]);
858   fAlgoIter0=(Int_t)(cuts[11]);
859
860   return;
861 }
862 //---------------------------------------------------------------------------
863 void AliVertexerTracks::SetITSMode(Double_t dcacut,
864                                    Double_t dcacutIter0,
865                                    Double_t maxd0z0,
866                                    Int_t minCls,
867                                    Int_t mintrks,
868                                    Double_t nsigma,
869                                    Double_t mindetfitter,
870                                    Double_t maxtgl,
871                                    Double_t fidR,
872                                    Double_t fidZ,
873                                    Int_t finderAlgo,
874                                    Int_t finderAlgoIter0)
875 {
876 //
877 //  Cut values for ITS mode
878 //
879   fMode = 0;
880   if(minCls>0) {
881     SetITSrefitRequired();
882   } else {
883     SetITSrefitNotRequired();
884   }
885   SetDCAcut(dcacut);
886   SetDCAcutIter0(dcacutIter0);
887   SetMaxd0z0(maxd0z0);
888   SetMinClusters(TMath::Abs(minCls));
889   SetMinTracks(mintrks);
890   SetNSigmad0(nsigma);
891   SetMinDetFitter(mindetfitter);
892   SetMaxTgl(maxtgl);
893   SetFiducialRZ(fidR,fidZ);
894   fAlgo=finderAlgo;
895   fAlgoIter0=finderAlgoIter0;
896
897   return; 
898 }
899 //---------------------------------------------------------------------------
900 void AliVertexerTracks::SetTPCMode(Double_t dcacut,
901                                    Double_t dcacutIter0,
902                                    Double_t maxd0z0,
903                                    Int_t minCls,
904                                    Int_t mintrks,
905                                    Double_t nsigma,
906                                    Double_t mindetfitter,
907                                    Double_t maxtgl,
908                                    Double_t fidR,
909                                    Double_t fidZ,
910                                    Int_t finderAlgo,
911                                    Int_t finderAlgoIter0) 
912 {
913 //
914 //  Cut values for TPC mode
915 //
916   fMode = 1;
917   SetITSrefitNotRequired();
918   SetDCAcut(dcacut);
919   SetDCAcutIter0(dcacutIter0);
920   SetMaxd0z0(maxd0z0);
921   SetMinClusters(minCls);
922   SetMinTracks(mintrks);
923   SetNSigmad0(nsigma);
924   SetMinDetFitter(mindetfitter);
925   SetMaxTgl(maxtgl);
926   SetFiducialRZ(fidR,fidZ);
927   fAlgo=finderAlgo;
928   fAlgoIter0=finderAlgoIter0;
929
930   return; 
931 }
932 //---------------------------------------------------------------------------
933 void AliVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped) 
934 {
935 //
936 // Mark the tracks not to be used in the vertex reconstruction.
937 // Tracks are identified by AliESDtrack::GetID()
938 //
939   fNTrksToSkip = n;  fTrksToSkip = new Int_t[n]; 
940   for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i]; 
941   return; 
942 }
943 //---------------------------------------------------------------------------
944 void  AliVertexerTracks::SetVtxStart(AliESDVertex *vtx) 
945
946 //
947 // Set initial vertex knowledge
948 //
949   vtx->GetXYZ(fNominalPos);
950   vtx->GetCovMatrix(fNominalCov);
951   SetConstraintOn();
952   return; 
953 }
954 //---------------------------------------------------------------------------
955 void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
956 {
957   AliExternalTrackParam *track1;
958   const Int_t knacc = (Int_t)fTrkArraySel.GetEntriesFast();
959   AliStrLine **linarray = new AliStrLine* [knacc];
960   for(Int_t i=0; i<knacc; i++){
961     track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
962     Double_t alpha=track1->GetAlpha();
963     Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
964     Double_t pos[3],dir[3],sigmasq[3]; 
965     track1->GetXYZAt(mindist,GetFieldkG(),pos);
966     track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
967     sigmasq[0]=TMath::Sin(alpha)*TMath::Sin(alpha)*track1->GetSigmaY2();
968     sigmasq[1]=TMath::Cos(alpha)*TMath::Cos(alpha)*track1->GetSigmaY2();
969     sigmasq[2]=track1->GetSigmaZ2();
970     TMatrixD ri(3,1);
971     TMatrixD wWi(3,3);
972     if(!TrackToPoint(track1,ri,wWi)) {optUseWeights=kFALSE;printf("WARNING\n");}
973     Double_t wmat[9];
974     Int_t iel=0;
975     for(Int_t ia=0;ia<3;ia++){
976       for(Int_t ib=0;ib<3;ib++){
977         wmat[iel]=wWi(ia,ib);
978         iel++;
979       }    
980     }
981     linarray[i] = new AliStrLine(pos,sigmasq,wmat,dir);
982   }
983   fVert=TrackletVertexFinder(linarray,knacc,optUseWeights);
984   for(Int_t i=0; i<knacc; i++) delete linarray[i];
985   delete [] linarray;
986 }
987 //---------------------------------------------------------------------------
988 AliESDVertex AliVertexerTracks::TrackletVertexFinder(TClonesArray *lines, Int_t optUseWeights)
989 {
990   // Calculate the point at minimum distance to prepared tracks (TClonesArray)
991   const Int_t knacc = (Int_t)lines->GetEntriesFast();
992   AliStrLine** lines2 = new AliStrLine* [knacc];
993   for(Int_t i=0; i<knacc; i++){
994     lines2[i]= (AliStrLine*)lines->At(i);
995   }
996   AliESDVertex vert = TrackletVertexFinder(lines2,knacc,optUseWeights); 
997   delete [] lines2;
998   return vert;
999 }
1000
1001 //---------------------------------------------------------------------------
1002 AliESDVertex AliVertexerTracks::TrackletVertexFinder(AliStrLine **lines, const Int_t knacc, Int_t optUseWeights)
1003 {
1004   // Calculate the point at minimum distance to prepared tracks (array of AliStrLine) 
1005
1006   Double_t initPos[3]={0.,0.,0.};
1007
1008   Double_t (*vectP0)[3]=new Double_t [knacc][3];
1009   Double_t (*vectP1)[3]=new Double_t [knacc][3];
1010   
1011   Double_t sum[3][3];
1012   Double_t dsum[3]={0,0,0};
1013   TMatrixD sumWi(3,3);
1014   for(Int_t i=0;i<3;i++){
1015     for(Int_t j=0;j<3;j++){
1016       sum[i][j]=0;
1017       sumWi(i,j)=0.;
1018     }
1019   }
1020
1021   for(Int_t i=0; i<knacc; i++){
1022     AliStrLine *line1 = lines[i]; 
1023     Double_t p0[3],cd[3],sigmasq[3];
1024     Double_t wmat[9];
1025     if(!line1) printf("ERROR %d %d\n",i,knacc);
1026     line1->GetP0(p0);
1027     line1->GetCd(cd);
1028     line1->GetSigma2P0(sigmasq);
1029     line1->GetWMatrix(wmat);
1030     TMatrixD wWi(3,3);
1031     Int_t iel=0;
1032     for(Int_t ia=0;ia<3;ia++){
1033       for(Int_t ib=0;ib<3;ib++){
1034         wWi(ia,ib)=wmat[iel];
1035         iel++;
1036       }    
1037     }
1038
1039     sumWi+=wWi;
1040
1041     Double_t p1[3]={p0[0]+cd[0],p0[1]+cd[1],p0[2]+cd[2]};
1042     vectP0[i][0]=p0[0];
1043     vectP0[i][1]=p0[1];
1044     vectP0[i][2]=p0[2];
1045     vectP1[i][0]=p1[0];
1046     vectP1[i][1]=p1[1];
1047     vectP1[i][2]=p1[2];
1048     
1049     Double_t matr[3][3];
1050     Double_t dknow[3];
1051     if(optUseWeights==0)GetStrLinDerivMatrix(p0,p1,matr,dknow);
1052     else GetStrLinDerivMatrix(p0,p1,sigmasq,matr,dknow);
1053
1054
1055     for(Int_t iii=0;iii<3;iii++){
1056       dsum[iii]+=dknow[iii]; 
1057       for(Int_t lj=0;lj<3;lj++) sum[iii][lj]+=matr[iii][lj];
1058     }
1059   }
1060  
1061   TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
1062   Double_t covmatrix[6];
1063   covmatrix[0] = invsumWi(0,0);
1064   covmatrix[1] = invsumWi(0,1);
1065   covmatrix[2] = invsumWi(1,1);
1066   covmatrix[3] = invsumWi(0,2);
1067   covmatrix[4] = invsumWi(1,2);
1068   covmatrix[5] = invsumWi(2,2);
1069
1070   Double_t vett[3][3];
1071   Double_t det=GetDeterminant3X3(sum);
1072   Double_t sigma=0;
1073   
1074   if(det!=0){
1075     for(Int_t zz=0;zz<3;zz++){
1076       for(Int_t ww=0;ww<3;ww++){
1077         for(Int_t kk=0;kk<3;kk++) vett[ww][kk]=sum[ww][kk];
1078       }
1079       for(Int_t kk=0;kk<3;kk++) vett[kk][zz]=dsum[kk];
1080       initPos[zz]=GetDeterminant3X3(vett)/det;
1081     }
1082
1083
1084     for(Int_t i=0; i<knacc; i++){
1085       Double_t p0[3]={0,0,0},p1[3]={0,0,0};
1086       for(Int_t ii=0;ii<3;ii++){
1087         p0[ii]=vectP0[i][ii];
1088         p1[ii]=vectP1[i][ii];
1089       }
1090       sigma+=GetStrLinMinDist(p0,p1,initPos);
1091     }
1092
1093     if(sigma>0.) {sigma=TMath::Sqrt(sigma);}else{sigma=999;}
1094   }else{
1095     sigma=999;
1096   }
1097   AliESDVertex theVert(initPos,covmatrix,99999.,knacc);
1098   theVert.SetDispersion(sigma);
1099   delete vectP0;
1100   delete vectP1;
1101   return theVert;
1102 }
1103
1104 //---------------------------------------------------------------------------
1105 Bool_t AliVertexerTracks::TrackToPoint(AliExternalTrackParam *t,
1106                                        TMatrixD &ri,TMatrixD &wWi,
1107                                        Bool_t uUi3by3) const 
1108 {
1109 //
1110 // Extract from the AliExternalTrackParam the global coordinates ri and covariance matrix
1111 // wWi of the space point that it represents (to be used in VertexFitter())
1112 //
1113
1114   
1115   Double_t rotAngle = t->GetAlpha();
1116   if(rotAngle<0.) rotAngle += 2.*TMath::Pi();
1117   Double_t cosRot = TMath::Cos(rotAngle);
1118   Double_t sinRot = TMath::Sin(rotAngle);
1119
1120   ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot;
1121   ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot;
1122   ri(2,0) = t->GetZ();
1123
1124
1125
1126   if(!uUi3by3) {
1127     // matrix to go from global (x,y,z) to local (y,z);
1128     TMatrixD qQi(2,3);
1129     qQi(0,0) = -sinRot;
1130     qQi(0,1) = cosRot;
1131     qQi(0,2) = 0.;
1132     qQi(1,0) = 0.;
1133     qQi(1,1) = 0.;
1134     qQi(1,2) = 1.;
1135
1136     // covariance matrix of local (y,z) - inverted
1137     TMatrixD uUi(2,2);
1138     uUi(0,0) = t->GetSigmaY2();
1139     uUi(0,1) = t->GetSigmaZY();
1140     uUi(1,0) = t->GetSigmaZY();
1141     uUi(1,1) = t->GetSigmaZ2();
1142     //printf(" Ui :");
1143     //printf(" %f   %f",uUi(0,0),uUi(0,1));
1144     //printf(" %f   %f",uUi(1,0),uUi(1,1));
1145
1146     if(uUi.Determinant() <= 0.) return kFALSE;
1147     TMatrixD uUiInv(TMatrixD::kInverted,uUi);
1148   
1149     // weights matrix: wWi = qQiT * uUiInv * qQi
1150     TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
1151     TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
1152     wWi = m;
1153   } else {
1154     if(fVert.GetNContributors()<1) AliFatal("Vertex from finder is empty");
1155     // matrix to go from global (x,y,z) to local (x,y,z);
1156     TMatrixD qQi(3,3);
1157     qQi(0,0) = cosRot;
1158     qQi(0,1) = sinRot;
1159     qQi(0,2) = 0.;
1160     qQi(1,0) = -sinRot;
1161     qQi(1,1) = cosRot;
1162     qQi(1,2) = 0.;
1163     qQi(2,0) = 0.;
1164     qQi(2,1) = 0.;
1165     qQi(2,2) = 1.;
1166    
1167     // covariance of fVert along the track  
1168     Double_t p[3],pt,ptot;
1169     t->GetPxPyPz(p);
1170     pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
1171     ptot = TMath::Sqrt(pt*pt+p[2]*p[2]);
1172     Double_t cphi = p[0]/pt;               //cos(phi)=px/pt
1173     Double_t sphi = p[1]/pt;               //sin(phi)=py/pt
1174     Double_t clambda = pt/ptot;            //cos(lambda)=pt/ptot
1175     Double_t slambda = p[2]/ptot;            //sin(lambda)=pz/ptot
1176     Double_t covfVert[6];
1177     fVert.GetCovMatrix(covfVert);
1178     Double_t covfVertalongt = 
1179        covfVert[0]*cphi*cphi*clambda*clambda 
1180       +covfVert[1]*2.*cphi*sphi*clambda*clambda
1181       +covfVert[3]*2.*cphi*clambda*slambda 
1182       +covfVert[2]*sphi*sphi*clambda*clambda 
1183       +covfVert[4]*2.*sphi*clambda*slambda 
1184       +covfVert[5]*slambda*slambda; 
1185     // covariance matrix of local (x,y,z) - inverted
1186     TMatrixD uUi(3,3);
1187     uUi(0,0) = covfVertalongt * fnSigmaForUi00 * fnSigmaForUi00;
1188     AliDebug(1,Form("=====> sqrtUi00 cm  %f",TMath::Sqrt(uUi(0,0))));
1189     uUi(0,1) = 0.;
1190     uUi(0,2) = 0.;
1191     uUi(1,0) = 0.;
1192     uUi(1,1) = t->GetSigmaY2();
1193     uUi(1,2) = t->GetSigmaZY();
1194     uUi(2,0) = 0.;
1195     uUi(2,1) = t->GetSigmaZY();
1196     uUi(2,2) = t->GetSigmaZ2();
1197     //printf(" Ui :\n");
1198     //printf(" %f   %f\n",uUi(0,0),uUi(0,1));
1199     //printf(" %f   %f\n",uUi(1,0),uUi(1,1));
1200   
1201     if(uUi.Determinant() <= 0.) return kFALSE;
1202     TMatrixD uUiInv(TMatrixD::kInverted,uUi);
1203   
1204     // weights matrix: wWi = qQiT * uUiInv * qQi
1205     TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
1206     TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
1207     wWi = m;
1208   }
1209
1210
1211   return kTRUE;
1212
1213 //---------------------------------------------------------------------------
1214 void AliVertexerTracks::TooFewTracks() 
1215 {
1216 //
1217 // When the number of tracks is < fMinTracks,
1218 // deal with vertices not found and prepare to exit
1219 //
1220   AliDebug(1,"TooFewTracks");
1221
1222   Double_t pos[3],err[3];
1223   pos[0] = fNominalPos[0];
1224   err[0] = TMath::Sqrt(fNominalCov[0]);
1225   pos[1] = fNominalPos[1];
1226   err[1] = TMath::Sqrt(fNominalCov[2]);
1227   pos[2] = fNominalPos[2];
1228   err[2] = TMath::Sqrt(fNominalCov[5]);
1229   Int_t    ncontr = (err[0]>1. ? -1 : -3);
1230   if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
1231   fCurrentVertex = new AliESDVertex(pos,err);
1232   fCurrentVertex->SetNContributors(ncontr);
1233
1234   if(fConstraint) {
1235     fCurrentVertex->SetTitle("VertexerTracksWithConstraint");
1236   } else {
1237     fCurrentVertex->SetTitle("VertexerTracksNoConstraint");
1238   }
1239
1240   if(!fTrkArraySel.IsEmpty()) fTrkArraySel.Delete(); 
1241   if(fIdSel) {delete fIdSel; fIdSel=NULL;}
1242   if(fTrksToSkip) {delete fTrksToSkip; fTrksToSkip=NULL;}
1243
1244   return;
1245 }
1246 //---------------------------------------------------------------------------
1247 void AliVertexerTracks::VertexFinder(Int_t optUseWeights) 
1248 {
1249
1250   // Get estimate of vertex position in (x,y) from tracks DCA
1251  
1252   Double_t initPos[3];
1253   initPos[2] = 0.;
1254   for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
1255   Int_t nacc = (Int_t)fTrkArraySel.GetEntriesFast();
1256   Double_t aver[3]={0.,0.,0.};
1257   Double_t aversq[3]={0.,0.,0.};
1258   Double_t sigmasq[3]={0.,0.,0.};
1259   Double_t sigma=0;
1260   Int_t ncombi = 0;
1261   AliExternalTrackParam *track1;
1262   AliExternalTrackParam *track2;
1263   Double_t pos[3],dir[3]; 
1264   Double_t alpha,mindist;
1265
1266   for(Int_t i=0; i<nacc; i++){
1267     track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
1268     alpha=track1->GetAlpha();
1269     mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
1270     track1->GetXYZAt(mindist,GetFieldkG(),pos);
1271     track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
1272     AliStrLine *line1 = new AliStrLine(pos,dir); 
1273
1274    //    AliStrLine *line1 = new AliStrLine();
1275    //    track1->ApproximateHelixWithLine(mindist,GetFieldkG(),line1);
1276    
1277     for(Int_t j=i+1; j<nacc; j++){
1278       track2 = (AliExternalTrackParam*)fTrkArraySel.At(j);
1279       alpha=track2->GetAlpha();
1280       mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
1281       track2->GetXYZAt(mindist,GetFieldkG(),pos);
1282       track2->GetPxPyPzAt(mindist,GetFieldkG(),dir);
1283       AliStrLine *line2 = new AliStrLine(pos,dir); 
1284     //      AliStrLine *line2 = new AliStrLine();
1285     //  track2->ApproximateHelixWithLine(mindist,GetFieldkG(),line2);
1286       Double_t distCA=line2->GetDCA(line1);
1287       //printf("%d   %d   %f\n",i,j,distCA);
1288        if(fDCAcut<=0 || (fDCAcut>0&&distCA<fDCAcut)){
1289         Double_t pnt1[3],pnt2[3],crosspoint[3];
1290
1291         if(optUseWeights<=0){
1292           Int_t retcode = line2->Cross(line1,crosspoint);
1293           if(retcode>=0){
1294             ncombi++;
1295             for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
1296             for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
1297           }
1298         }
1299         if(optUseWeights>0){
1300           Int_t retcode = line1->CrossPoints(line2,pnt1,pnt2);
1301           if(retcode>=0){
1302             Double_t cs, sn;
1303             alpha=track1->GetAlpha();
1304             cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);   
1305             Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
1306             alpha=track2->GetAlpha();
1307             cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
1308             Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
1309             Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
1310             Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
1311             Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
1312             Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
1313             crosspoint[0]=wx1*pnt1[0] + wx2*pnt2[0]; 
1314             crosspoint[1]=wy1*pnt1[1] + wy2*pnt2[1]; 
1315             crosspoint[2]=wz1*pnt1[2] + wz2*pnt2[2];
1316           
1317             ncombi++;
1318             for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
1319             for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
1320           }
1321         }
1322       }
1323       delete line2;
1324     }
1325     delete line1;
1326   }
1327   if(ncombi>0){
1328     for(Int_t jj=0;jj<3;jj++){
1329       initPos[jj] = aver[jj]/ncombi;
1330       //printf("%f\n",initPos[jj]);
1331       aversq[jj]/=ncombi;
1332       sigmasq[jj]=aversq[jj]-initPos[jj]*initPos[jj];
1333       sigma+=sigmasq[jj];
1334     }
1335     sigma=TMath::Sqrt(TMath::Abs(sigma));
1336   }
1337   else {
1338     Warning("VertexFinder","Finder did not succed");
1339     sigma=999;
1340   }
1341   fVert.SetXYZ(initPos);
1342   fVert.SetDispersion(sigma);
1343   fVert.SetNContributors(ncombi);
1344 }
1345 //---------------------------------------------------------------------------
1346 void AliVertexerTracks::VertexFitter() 
1347 {
1348 //
1349 // The optimal estimate of the vertex position is given by a "weighted 
1350 // average of tracks positions".
1351 // Original method: V. Karimaki, CMS Note 97/0051
1352 //
1353   Bool_t useConstraint = fConstraint;
1354   Double_t initPos[3];
1355   if(!fOnlyFitter) {
1356     fVert.GetXYZ(initPos);
1357   } else {
1358     initPos[0]=fNominalPos[0];
1359     initPos[1]=fNominalPos[1];
1360     initPos[2]=fNominalPos[2];
1361   }
1362
1363   Int_t nTrksSel = (Int_t)fTrkArraySel.GetEntries();
1364   if(nTrksSel==1) useConstraint=kTRUE;
1365   AliDebug(1,"--- VertexFitter(): start");
1366   AliDebug(1,Form(" Number of tracks in array: %d\n",nTrksSel));
1367   AliDebug(1,Form(" Minimum # tracks required in fit: %d\n",fMinTracks));
1368   AliDebug(1,Form(" Vertex position after finder: %f,%f,%f\n",initPos[0],initPos[1],initPos[2]));
1369   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]))); 
1370
1371   // special treatment for few-tracks fits (e.g. secondary vertices)
1372   Bool_t uUi3by3 = kFALSE; if(nTrksSel<5 && !useConstraint) uUi3by3 = kTRUE;
1373
1374   Int_t i,j,k,step=0;
1375   TMatrixD rv(3,1);
1376   TMatrixD vV(3,3);
1377   rv(0,0) = initPos[0];
1378   rv(1,0) = initPos[1];
1379   rv(2,0) = 0.;
1380   Double_t xlStart,alpha;
1381   Int_t nTrksUsed;
1382   Double_t chi2,chi2i,chi2b;
1383   AliExternalTrackParam *t = 0;
1384   Int_t failed = 0;
1385
1386   // initial vertex covariance matrix
1387   TMatrixD vVb(3,3);
1388   vVb(0,0) = fNominalCov[0];
1389   vVb(0,1) = fNominalCov[1];
1390   vVb(0,2) = 0.;
1391   vVb(1,0) = fNominalCov[1];
1392   vVb(1,1) = fNominalCov[2];
1393   vVb(1,2) = 0.;
1394   vVb(2,0) = 0.;
1395   vVb(2,1) = 0.;
1396   vVb(2,2) = fNominalCov[5];
1397   TMatrixD vVbInv(TMatrixD::kInverted,vVb);
1398   TMatrixD rb(3,1);
1399   rb(0,0) = fNominalPos[0];
1400   rb(1,0) = fNominalPos[1];
1401   rb(2,0) = fNominalPos[2];
1402   TMatrixD vVbInvrb(vVbInv,TMatrixD::kMult,rb);
1403
1404
1405   // 2 steps:
1406   // 1st - estimate of vtx using all tracks
1407   // 2nd - estimate of global chi2
1408   for(step=0; step<2; step++) {
1409     AliDebug(1,Form(" step = %d\n",step));
1410     chi2 = 0.;
1411     nTrksUsed = 0;
1412
1413     if(step==1) { initPos[0]=rv(0,0); initPos[0]=rv(1,0); }
1414
1415     TMatrixD sumWiri(3,1);
1416     TMatrixD sumWi(3,3);
1417     for(i=0; i<3; i++) {
1418       sumWiri(i,0) = 0.;
1419       for(j=0; j<3; j++) sumWi(j,i) = 0.;
1420     }
1421
1422     // mean vertex constraint
1423     if(useConstraint) {
1424       for(i=0;i<3;i++) {
1425         sumWiri(i,0) += vVbInvrb(i,0);
1426         for(k=0;k<3;k++) sumWi(i,k) += vVbInv(i,k);
1427       }
1428       // chi2
1429       TMatrixD deltar = rv; deltar -= rb;
1430       TMatrixD vVbInvdeltar(vVbInv,TMatrixD::kMult,deltar);
1431       chi2b = deltar(0,0)*vVbInvdeltar(0,0)+
1432               deltar(1,0)*vVbInvdeltar(1,0)+
1433               deltar(2,0)*vVbInvdeltar(2,0);
1434       chi2 += chi2b;
1435     }
1436
1437     // loop on tracks  
1438     for(k=0; k<nTrksSel; k++) {
1439       // get track from track array
1440       t = (AliExternalTrackParam*)fTrkArraySel.At(k);
1441       alpha = t->GetAlpha();
1442       xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha);
1443       // to vtxSeed (from finder)
1444       t->PropagateTo(xlStart,GetFieldkG());   
1445  
1446       // vector of track global coordinates
1447       TMatrixD ri(3,1);
1448       // covariance matrix of ri
1449       TMatrixD wWi(3,3);
1450
1451       // get space point from track
1452       if(!TrackToPoint(t,ri,wWi,uUi3by3)) continue;
1453       TMatrixD wWiri(wWi,TMatrixD::kMult,ri); 
1454
1455       // track chi2
1456       TMatrixD deltar = rv; deltar -= ri;
1457       TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
1458       chi2i = deltar(0,0)*wWideltar(0,0)+
1459               deltar(1,0)*wWideltar(1,0)+
1460               deltar(2,0)*wWideltar(2,0);
1461
1462       // add to total chi2
1463       chi2 += chi2i;
1464
1465       sumWiri += wWiri;
1466       sumWi   += wWi;
1467
1468       nTrksUsed++;
1469     } // end loop on tracks
1470
1471     if(nTrksUsed < fMinTracks) {
1472       failed=1;
1473       continue;
1474     }
1475
1476     Double_t determinant = sumWi.Determinant();
1477     if(determinant < fMinDetFitter)  { 
1478       AliDebug(1,Form("det(V) = %f (<%f)\n",determinant,fMinDetFitter));       
1479       failed=1;
1480       continue;
1481     }
1482
1483     if(step==0) { 
1484       // inverted of weights matrix
1485       TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
1486       vV = invsumWi;
1487       // position of primary vertex
1488       rv.Mult(vV,sumWiri);
1489     }
1490   } // end loop on the 2 steps
1491
1492   if(failed) { 
1493     TooFewTracks();
1494     return; 
1495   }
1496
1497   Double_t position[3];
1498   position[0] = rv(0,0);
1499   position[1] = rv(1,0);
1500   position[2] = rv(2,0);
1501   Double_t covmatrix[6];
1502   covmatrix[0] = vV(0,0);
1503   covmatrix[1] = vV(0,1);
1504   covmatrix[2] = vV(1,1);
1505   covmatrix[3] = vV(0,2);
1506   covmatrix[4] = vV(1,2);
1507   covmatrix[5] = vV(2,2);
1508   
1509   // for correct chi2/ndf, count constraint as additional "track"
1510   if(fConstraint) nTrksUsed++;
1511   // store data in the vertex object
1512   if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
1513   fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);
1514
1515   AliDebug(1," Vertex after fit:");
1516   AliDebug(1,Form("xyz: %f %f %f; nc %d",fCurrentVertex->GetXv(),fCurrentVertex->GetYv(),fCurrentVertex->GetZv(),fCurrentVertex->GetNContributors()));
1517   AliDebug(1,"--- VertexFitter(): finish\n");
1518  
1519
1520   return;
1521 }
1522 //----------------------------------------------------------------------------
1523 AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray,
1524                                                          UShort_t *id,
1525                                                          Bool_t optUseFitter,
1526                                                          Bool_t optPropagate) 
1527 {
1528 //
1529 // Return vertex from tracks (AliExternalTrackParam) in array
1530 //
1531   fCurrentVertex = 0;
1532   SetConstraintOff();
1533
1534   // get tracks and propagate them to initial vertex position
1535   fIdSel = new UShort_t[(Int_t)trkArray->GetEntriesFast()];
1536   Int_t nTrksSel = PrepareTracks(*trkArray,id,0);
1537   if(nTrksSel <  TMath::Max(2,fMinTracks)) {
1538     TooFewTracks();
1539     return fCurrentVertex;
1540   }
1541  
1542   // vertex finder
1543   switch (fAlgo) {
1544     case 1: StrLinVertexFinderMinDist(1); break;
1545     case 2: StrLinVertexFinderMinDist(0); break;
1546     case 3: HelixVertexFinder();          break;
1547     case 4: VertexFinder(1);              break;
1548     case 5: VertexFinder(0);              break;
1549     default: printf("Wrong algorithm\n"); break;  
1550   }
1551   AliDebug(1," Vertex finding completed\n");
1552
1553   // vertex fitter
1554   if(optUseFitter) {
1555     VertexFitter();
1556   } else {
1557     Double_t position[3]={fVert.GetXv(),fVert.GetYv(),fVert.GetZv()};
1558     Double_t covmatrix[6];
1559     fVert.GetCovMatrix(covmatrix);
1560     Double_t chi2=99999.;
1561     Int_t    nTrksUsed=fVert.GetNContributors();
1562     fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);    
1563   }
1564   fCurrentVertex->SetDispersion(fVert.GetDispersion());
1565
1566
1567   // set indices of used tracks and propagate track to found vertex
1568   UShort_t *indices = 0;
1569   Double_t d0z0[2],covd0z0[3];
1570   AliExternalTrackParam *t = 0;
1571   if(fCurrentVertex->GetNContributors()>0) {
1572     indices = new UShort_t[fTrkArraySel.GetEntriesFast()];
1573     for(Int_t jj=0; jj<(Int_t)fTrkArraySel.GetEntriesFast(); jj++) {
1574       indices[jj] = fIdSel[jj];
1575       t = (AliExternalTrackParam*)fTrkArraySel.At(jj);
1576       if(optPropagate && optUseFitter) {
1577         if(TMath::Sqrt(fCurrentVertex->GetXv()*fCurrentVertex->GetXv()+fCurrentVertex->GetYv()*fCurrentVertex->GetYv())<3.) {
1578           t->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
1579           AliDebug(1,Form("Track %d propagated to found vertex",jj));
1580         } else {
1581           AliWarning("Found vertex outside beam pipe!");
1582         }
1583       }
1584     }
1585     fCurrentVertex->SetIndices(fCurrentVertex->GetNContributors(),indices);
1586   }
1587
1588   // clean up
1589   if (indices) {delete indices; indices=NULL;}
1590   delete fIdSel; fIdSel=NULL;
1591   fTrkArraySel.Delete();
1592   
1593   return fCurrentVertex;
1594 }
1595 //----------------------------------------------------------------------------
1596 AliESDVertex* AliVertexerTracks::VertexForSelectedESDTracks(TObjArray *trkArray,
1597                                                          Bool_t optUseFitter,
1598                                                          Bool_t optPropagate) 
1599 {
1600 //
1601 // Return vertex from array of ESD tracks
1602 //
1603
1604   Int_t nTrks = (Int_t)trkArray->GetEntriesFast();
1605   UShort_t *id = new UShort_t[nTrks];
1606
1607   AliESDtrack *esdt = 0;
1608   for(Int_t i=0; i<nTrks; i++){
1609     esdt = (AliESDtrack*)trkArray->At(i);
1610     id[i] = (UShort_t)esdt->GetID();
1611   }
1612     
1613   VertexForSelectedTracks(trkArray,id,optUseFitter,optPropagate);
1614
1615   delete id; id=NULL;
1616
1617   return fCurrentVertex;
1618 }
1619 //--------------------------------------------------------------------------