A few improvements, a small bug fix, and removal of the dependency on AliTracker...
[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 ESD 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 <TTree.h>
30 #include <TMatrixD.h>
31 //---- AliRoot headers -----
32 #include "AliStrLine.h"
33 #include "AliVertexerTracks.h"
34 #include "AliESD.h"
35 #include "AliESDtrack.h"
36
37 ClassImp(AliVertexerTracks)
38
39
40 //----------------------------------------------------------------------------
41 AliVertexerTracks::AliVertexerTracks():
42 TObject(),
43 fVert(),
44 fCurrentVertex(0),
45 fFieldkG(-999.),
46 fConstraint(kFALSE),
47 fOnlyFitter(kFALSE),
48 fMinTracks(1),
49 fMinITSClusters(5),
50 fTrkArray(),
51 fTrksToSkip(0),
52 fNTrksToSkip(0),
53 fDCAcut(0),
54 fAlgo(1),
55 fNSigma(3),
56 fMaxd0z0(0.5),
57 fITSin(kTRUE),
58 fITSrefit(kTRUE),
59 fDebug(0)
60 {
61 //
62 // Default constructor
63 //
64   SetVtxStart();
65   SetVtxStartSigma();
66   SetMinTracks();
67   SetMinITSClusters();
68   SetNSigmad0(); 
69   SetMaxd0z0(); 
70 }
71 //----------------------------------------------------------------------------
72 AliVertexerTracks::AliVertexerTracks(Double_t fieldkG):
73 TObject(),
74 fVert(),
75 fCurrentVertex(0),
76 fFieldkG(-999.),
77 fConstraint(kFALSE),
78 fOnlyFitter(kFALSE),
79 fMinTracks(1),
80 fMinITSClusters(5),
81 fTrkArray(),
82 fTrksToSkip(0),
83 fNTrksToSkip(0),
84 fDCAcut(0),
85 fAlgo(1),
86 fNSigma(3),
87 fMaxd0z0(0.5),
88 fITSin(kTRUE),
89 fITSrefit(kTRUE),
90 fDebug(0)
91 {
92 //
93 // Standard constructor
94 //
95   SetVtxStart();
96   SetVtxStartSigma();
97   SetMinTracks();
98   SetMinITSClusters();
99   SetNSigmad0(); 
100   SetMaxd0z0();
101   SetFieldkG(fieldkG);
102 }
103 //-----------------------------------------------------------------------------
104 AliVertexerTracks::~AliVertexerTracks() 
105 {
106   // Default Destructor
107   // The objects pointed by the following pointer are not owned
108   // by this class and are not deleted
109   fCurrentVertex = 0;
110   delete[] fTrksToSkip;
111 }
112 //----------------------------------------------------------------------------
113 AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliESD *esdEvent)
114 {
115 //
116 // Primary vertex for current ESD event
117 // (Two iterations: 
118 //  1st with 5*fNSigma*sigma cut w.r.t. to initial vertex
119 //      + cut on sqrt(d0d0+z0z0) if fConstraint=kFALSE  
120 //  2nd with fNSigma*sigma cut w.r.t. to vertex found in 1st iteration) 
121 // All ESD tracks with inside the beam pipe are then propagated to found vertex
122 //
123   fCurrentVertex = 0;
124
125   // accept 1-track case only if constraint is available
126   if(!fConstraint && fMinTracks==1) fMinTracks=2;
127
128   // get Bz from ESD
129   //  SetFieldkG(esdEvent->GetMagneticField());
130
131   // read tracks from ESD
132   Int_t nTrksTot = (Int_t)esdEvent->GetNumberOfTracks();
133   if (nTrksTot<=0) {
134     if(fDebug) printf("TooFewTracks\n");
135     TooFewTracks(esdEvent);
136     return fCurrentVertex;
137   } 
138
139   TTree *trkTree = new TTree("TreeT","tracks");
140   AliESDtrack *esdTrack = 0;
141   trkTree->Branch("tracks","AliESDtrack",&esdTrack);
142
143   Bool_t   skipThis;
144   for(Int_t i=0; i<nTrksTot; i++) {
145     AliESDtrack *et = esdEvent->GetTrack(i);
146     esdTrack = new AliESDtrack(*et);
147     // check tracks to skip
148     skipThis = kFALSE;
149     for(Int_t j=0; j<fNTrksToSkip; j++) { 
150       if(et->GetID()==fTrksToSkip[j]) {
151         if(fDebug) printf("skipping track: %d\n",i);
152         skipThis = kTRUE;
153       }
154     }
155     if(skipThis) {delete esdTrack;continue;}
156     if(fITSin) {
157       if(!(esdTrack->GetStatus()&AliESDtrack::kITSin)) {delete esdTrack;continue;}
158       if(fITSrefit && !(esdTrack->GetStatus()&AliESDtrack::kITSrefit)) {delete esdTrack;continue;}
159       Int_t nclus=esdTrack->GetNcls(0); // check number of clusters in ITS
160       if(nclus<fMinITSClusters) {delete esdTrack;continue;}
161     }
162     trkTree->Fill();
163     delete esdTrack;
164   }
165   
166   // If fConstraint=kFALSE
167   // run VertexFinder(1) to get rough estimate of initVertex (x,y)
168   if(!fConstraint) {
169     // fill fTrkArray, for VertexFinder()
170     if(!fTrkArray.IsEmpty()) fTrkArray.Delete();
171     PrepareTracks(*trkTree,0);
172     Double_t cutsave = fDCAcut;  fDCAcut = 0.1; // 1 mm
173     VertexFinder(1); // using weights, cutting dca < fDCAcut
174     fDCAcut = cutsave;
175     fTrkArray.Delete();
176     if(fVert.GetNContributors()>0) {
177       fVert.GetXYZ(fNominalPos);
178       fNominalPos[0] = fVert.GetXv();
179       fNominalPos[1] = fVert.GetYv();
180       fNominalPos[2] = fVert.GetZv();
181       if(fDebug) printf("No mean vertex: VertexFinder gives (%f, %f, %f)\n",fNominalPos[0],fNominalPos[1],fNominalPos[2]);
182     } else {
183       fNominalPos[0] = 0.;
184       fNominalPos[1] = 0.;
185       fNominalPos[2] = 0.;
186       if(fDebug) printf("No mean vertex and VertexFinder failed\n");
187     }
188   }
189   
190   // TWO ITERATIONS:
191   //
192   // ITERATION 1
193   // propagate tracks to fNominalPos vertex
194   // preselect them:
195   // if(constraint) reject for |d0|>5*fNSigma*sigma w.r.t. fNominal... vertex
196   // else  reject for |d0|\oplus|z0| > 5 mm w.r.t. fNominal... vertex
197   // ITERATION 2
198   // propagate tracks to best between initVertex and fCurrentVertex
199   // preselect tracks (reject for |d0|>fNSigma*sigma w.r.t. best 
200   //                   between initVertex and fCurrentVertex) 
201   for(Int_t iter=0; iter<2; iter++) {
202     if(fOnlyFitter && iter==0) continue; 
203     Int_t nTrksPrep = PrepareTracks(*trkTree,iter+1);
204     if(fDebug) printf(" tracks prepared - iteration %d: %d\n",iter+1,nTrksPrep);
205     if(nTrksPrep < fMinTracks) {
206       if(fDebug) printf("TooFewTracks\n");
207       TooFewTracks(esdEvent);
208       if(fDebug) fCurrentVertex->PrintStatus();
209       fTrkArray.Delete();
210       delete trkTree;
211       return fCurrentVertex; 
212     }
213
214     // vertex finder
215     if(!fOnlyFitter) {
216       if(nTrksPrep==1){
217         if(fDebug) printf("Just one track\n");
218         OneTrackVertFinder();
219       }else{
220         switch (fAlgo) {
221         case 1: StrLinVertexFinderMinDist(1); break;
222         case 2: StrLinVertexFinderMinDist(0); break;
223         case 3: HelixVertexFinder();          break;
224         case 4: VertexFinder(1);              break;
225         case 5: VertexFinder(0);              break;
226         default: printf("Wrong algorithm\n"); break;  
227         }
228       }
229       if(fDebug) printf(" Vertex finding completed\n");
230     }
231
232     // vertex fitter
233     VertexFitter(fConstraint);
234     if(fDebug) printf(" Vertex fit completed\n");
235     if(iter==0) fTrkArray.Delete();
236   } // end loop on the two iterations
237
238
239   // take true pos from SPD vertex in ESD and write it in tracks' vertex
240   Double_t tp[3];
241   esdEvent->GetVertex()->GetTruePos(tp);
242   fCurrentVertex->SetTruePos(tp);
243
244   if(fConstraint) {
245     if(fOnlyFitter) {
246       fCurrentVertex->SetTitle("VertexerTracksWithConstraintOnlyFitter");
247     } else {
248       fCurrentVertex->SetTitle("VertexerTracksWithConstraint");
249     }
250   } else {
251     fCurrentVertex->SetTitle("VertexerTracksNoConstraint");
252   }
253
254     
255   // set indices of used tracks
256   UShort_t *indices = 0;
257   AliESDtrack *ett = 0;
258   if(fCurrentVertex->GetNContributors()>0) {
259     indices = new UShort_t[fCurrentVertex->GetNContributors()];
260     for(Int_t jj=0;jj<(Int_t)fTrkArray.GetEntriesFast();jj++) {
261       ett = (AliESDtrack*)fTrkArray.At(jj);
262       indices[jj] = (UShort_t)ett->GetID();
263     }
264     fCurrentVertex->SetIndices(fCurrentVertex->GetNContributors(),indices);
265   }
266   delete [] indices;
267
268   delete trkTree;
269   
270   fTrkArray.Delete();
271
272   if(fTrksToSkip) delete [] fTrksToSkip;
273
274
275   if(fDebug) fCurrentVertex->PrintStatus();
276   if(fDebug) fCurrentVertex->PrintIndices();
277
278   
279   return fCurrentVertex;
280 }
281 //------------------------------------------------------------------------
282 Double_t AliVertexerTracks::GetDeterminant3X3(Double_t matr[][3])
283 {
284   //
285   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];
286  return det;
287 }
288 //-------------------------------------------------------------------------
289 void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t (*m)[3],Double_t *d)
290 {
291   //
292   Double_t x12=p0[0]-p1[0];
293   Double_t y12=p0[1]-p1[1];
294   Double_t z12=p0[2]-p1[2];
295   Double_t kk=x12*x12+y12*y12+z12*z12;
296   m[0][0]=2-2/kk*x12*x12;
297   m[0][1]=-2/kk*x12*y12;
298   m[0][2]=-2/kk*x12*z12;
299   m[1][0]=-2/kk*x12*y12;
300   m[1][1]=2-2/kk*y12*y12;
301   m[1][2]=-2/kk*y12*z12;
302   m[2][0]=-2/kk*x12*z12;
303   m[2][1]=-2*y12*z12;
304   m[2][2]=2-2/kk*z12*z12;
305   d[0]=2*p0[0]-2/kk*p0[0]*x12*x12-2/kk*p0[2]*x12*z12-2/kk*p0[1]*x12*y12;
306   d[1]=2*p0[1]-2/kk*p0[1]*y12*y12-2/kk*p0[0]*x12*y12-2/kk*p0[2]*z12*y12;
307   d[2]=2*p0[2]-2/kk*p0[2]*z12*z12-2/kk*p0[0]*x12*z12-2/kk*p0[1]*z12*y12;
308
309 }
310 //--------------------------------------------------------------------------  
311 void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t *sigmasq,Double_t (*m)[3],Double_t *d)
312 {
313   //
314   Double_t x12=p1[0]-p0[0];
315   Double_t y12=p1[1]-p0[1];
316   Double_t z12=p1[2]-p0[2];
317
318   Double_t den= x12*x12*sigmasq[1]*sigmasq[2]+y12*y12*sigmasq[0]*sigmasq[2]+z12*z12*sigmasq[0]*sigmasq[1];
319
320   Double_t kk= 2*(x12*x12/sigmasq[0]+y12*y12/sigmasq[1]+z12*z12/sigmasq[2]);
321
322   Double_t cc[3];
323   cc[0]=-x12/sigmasq[0];
324   cc[1]=-y12/sigmasq[1];
325   cc[2]=-z12/sigmasq[2];
326
327   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;
328
329   Double_t ss= -p0[0]*cc[0]-p0[1]*cc[1]-p0[2]*cc[2];
330
331   Double_t aa[3];
332   aa[0]=x12*sigmasq[1]*sigmasq[2]/den;
333   aa[1]=y12*sigmasq[0]*sigmasq[2]/den;
334   aa[2]=z12*sigmasq[0]*sigmasq[1]/den;
335
336   m[0][0]=aa[0]*(aa[0]*kk+2*cc[0])+2*cc[0]*aa[0]+2/sigmasq[0];
337   m[0][1]=aa[1]*(aa[0]*kk+2*cc[0])+2*cc[1]*aa[0];
338   m[0][2]=aa[2]*(aa[0]*kk+2*cc[0])+2*cc[2]*aa[0];
339
340   m[1][0]=aa[0]*(aa[1]*kk+2*cc[1])+2*cc[0]*aa[1];
341   m[1][1]=aa[1]*(aa[1]*kk+2*cc[1])+2*cc[1]*aa[1]+2/sigmasq[1];
342   m[1][2]=aa[2]*(aa[1]*kk+2*cc[1])+2*cc[2]*aa[1];
343
344   m[2][0]=aa[0]*(aa[2]*kk+2*cc[2])+2*cc[0]*aa[2];
345   m[2][1]=aa[1]*(aa[2]*kk+2*cc[2])+2*cc[1]*aa[2];
346   m[2][2]=aa[2]*(aa[2]*kk+2*cc[2])+2*cc[2]*aa[2]+2/sigmasq[2];
347
348   d[0]=-ww*(aa[0]*kk+2*cc[0])-2*ss*aa[0]+2*p0[0]/sigmasq[0];
349   d[1]=-ww*(aa[1]*kk+2*cc[1])-2*ss*aa[1]+2*p0[1]/sigmasq[1];
350   d[2]=-ww*(aa[2]*kk+2*cc[2])-2*ss*aa[2]+2*p0[2]/sigmasq[2];
351
352   }
353 //--------------------------------------------------------------------------   
354 Double_t AliVertexerTracks::GetStrLinMinDist(Double_t *p0,Double_t *p1,Double_t *x0)
355 {
356   //
357   Double_t x12=p0[0]-p1[0];
358   Double_t y12=p0[1]-p1[1];
359   Double_t z12=p0[2]-p1[2];
360   Double_t x10=p0[0]-x0[0];
361   Double_t y10=p0[1]-x0[1];
362   Double_t z10=p0[2]-x0[2];
363   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);
364 }
365 //---------------------------------------------------------------------------
366 void AliVertexerTracks::OneTrackVertFinder() 
367 {
368   // find vertex for events with 1 track, using DCA to nominal beam axis
369   if(fDebug) printf("Number of prepared tracks =%d - Call OneTrackVertFinder",fTrkArray.GetEntries());
370   AliESDtrack *track1;
371   track1 = (AliESDtrack*)fTrkArray.At(0);
372   Double_t field=GetFieldkG();
373   Double_t alpha=track1->GetAlpha();
374   Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
375   Double_t pos[3],dir[3]; 
376   track1->GetXYZAt(mindist,field,pos);
377   track1->GetPxPyPzAt(mindist,field,dir);
378   AliStrLine *line1 = new AliStrLine(pos,dir);
379   Double_t p1[3]={fNominalPos[0],fNominalPos[1],0.}; 
380   Double_t p2[3]={fNominalPos[0],fNominalPos[1],10.}; 
381   AliStrLine *zeta=new AliStrLine(p1,p2,kTRUE);
382   Double_t crosspoint[3]={0.,0.,0.};
383   Double_t sigma=999.;
384   Int_t nContrib=-1;
385   Int_t retcode = zeta->Cross(line1,crosspoint);
386   if(retcode>=0){
387     sigma=line1->GetDistFromPoint(crosspoint);
388     nContrib=1;
389   }
390   delete zeta;
391   delete line1;
392   fVert.SetXYZ(crosspoint);
393   fVert.SetDispersion(sigma);
394   fVert.SetNContributors(nContrib);  
395 }
396 //---------------------------------------------------------------------------
397 void AliVertexerTracks::HelixVertexFinder() 
398 {
399   // Get estimate of vertex position in (x,y) from tracks DCA
400
401
402   Double_t initPos[3];
403   initPos[2] = 0.;
404   for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
405   Double_t field=GetFieldkG();
406
407   Int_t nacc = (Int_t)fTrkArray.GetEntriesFast();
408
409   Double_t aver[3]={0.,0.,0.};
410   Double_t averquad[3]={0.,0.,0.};
411   Double_t sigmaquad[3]={0.,0.,0.};
412   Double_t sigma=0;
413   Int_t ncombi = 0;
414   AliESDtrack *track1;
415   AliESDtrack *track2;
416   Double_t distCA;
417   Double_t x, par[5];
418   Double_t alpha, cs, sn;
419   Double_t crosspoint[3];
420   for(Int_t i=0; i<nacc; i++){
421     track1 = (AliESDtrack*)fTrkArray.At(i);
422     
423
424     for(Int_t j=i+1; j<nacc; j++){
425       track2 = (AliESDtrack*)fTrkArray.At(j);
426
427       distCA=track2->PropagateToDCA(track1,field);
428       if(fDCAcut<=0 ||(fDCAcut>0&&distCA<fDCAcut)){
429         track1->GetExternalParameters(x,par);
430         alpha=track1->GetAlpha();
431         cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
432         Double_t x1=x*cs - par[0]*sn;
433         Double_t y1=x*sn + par[0]*cs;
434         Double_t z1=par[1];
435         Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2(); 
436         track2->GetExternalParameters(x,par);
437         alpha=track2->GetAlpha();
438         cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
439         Double_t x2=x*cs - par[0]*sn;
440         Double_t y2=x*sn + par[0]*cs;
441         Double_t z2=par[1];
442         Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
443         Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
444         Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
445         Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
446         Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
447         crosspoint[0]=wx1*x1 + wx2*x2; 
448         crosspoint[1]=wy1*y1 + wy2*y2; 
449         crosspoint[2]=wz1*z1 + wz2*z2;
450
451         ncombi++;
452         for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
453         for(Int_t jj=0;jj<3;jj++)averquad[jj]+=(crosspoint[jj]*crosspoint[jj]);
454       }
455     }
456       
457   }
458   if(ncombi>0){
459     for(Int_t jj=0;jj<3;jj++){
460       initPos[jj] = aver[jj]/ncombi;
461       averquad[jj]/=ncombi;
462       sigmaquad[jj]=averquad[jj]-initPos[jj]*initPos[jj];
463       sigma+=sigmaquad[jj];
464     }
465     sigma=TMath::Sqrt(TMath::Abs(sigma));
466   }
467   else {
468     Warning("HelixVertexFinder","Finder did not succed");
469     sigma=999;
470   }
471   fVert.SetXYZ(initPos);
472   fVert.SetDispersion(sigma);
473   fVert.SetNContributors(ncombi);
474 }
475 //----------------------------------------------------------------------------
476 Int_t AliVertexerTracks::PrepareTracks(TTree &trkTree,Int_t optImpParCut) 
477 {
478 //
479 // Propagate tracks to initial vertex position and store them in a TObjArray
480 //
481   Double_t maxd0rphi; 
482   Double_t maxd0z0 = fMaxd0z0; // default is 5 mm  
483   Int_t    nTrks    = 0;
484   Double_t sigmaCurr[3];
485   Double_t normdistx,normdisty;
486   Float_t  d0z0[2],covd0z0[3]; 
487   Double_t sigma;
488   Double_t field=GetFieldkG();
489
490   AliESDVertex *initVertex = new AliESDVertex(fNominalPos,fNominalCov,1,1);
491
492   Int_t    nEntries = (Int_t)trkTree.GetEntries();
493   if(!fTrkArray.IsEmpty()) fTrkArray.Delete();
494
495   if(fDebug) {
496     printf(" PrepareTracks()\n");
497   }
498
499   for(Int_t i=0; i<nEntries; i++) {
500     AliESDtrack *track = new AliESDtrack; 
501     trkTree.SetBranchAddress("tracks",&track);
502     trkTree.GetEvent(i);
503
504     // propagate track to vertex
505     if(optImpParCut<=1 || fOnlyFitter) { // optImpParCut==1 or 0
506       track->RelateToVertex(initVertex,field,100.);
507     } else {              // optImpParCut==2
508       fCurrentVertex->GetSigmaXYZ(sigmaCurr);
509       normdistx = TMath::Abs(fCurrentVertex->GetXv()-fNominalPos[0])/TMath::Sqrt(sigmaCurr[0]*sigmaCurr[0]+fNominalCov[0]);
510       normdisty = TMath::Abs(fCurrentVertex->GetYv()-fNominalPos[1])/TMath::Sqrt(sigmaCurr[1]*sigmaCurr[1]+fNominalCov[2]);
511       if(normdistx < 3. && normdisty < 3. &&
512          (sigmaCurr[0]+sigmaCurr[1])<(TMath::Sqrt(fNominalCov[0])+TMath::Sqrt(fNominalCov[2]))) {
513         track->RelateToVertex(fCurrentVertex,field,100.);
514       } else {
515         track->RelateToVertex(initVertex,field,100.);
516       }
517     }
518
519     track->GetImpactParameters(d0z0,covd0z0);
520     sigma = TMath::Sqrt(covd0z0[0]);
521     maxd0rphi = fNSigma*sigma;
522     if(optImpParCut==1) maxd0rphi *= 5.;
523
524
525
526     if(fDebug) printf("trk %d; lab %d; |d0| = %f;  d0 cut = %f; |z0| = %f; |d0|oplus|z0| = %f; d0z0 cut = %f\n",i,track->GetLabel(),TMath::Abs(d0z0[0]),maxd0rphi,TMath::Abs(d0z0[1]),TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1]),maxd0z0);
527
528     // during iterations 1 and 2, if fConstraint=kFALSE,
529     // select tracks with d0oplusz0 < maxd0z0
530     if(optImpParCut>=1 && !fConstraint && nEntries>=3 && 
531        fVert.GetNContributors()>0) {
532       if(TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1]) > maxd0z0) { 
533         if(fDebug) printf("     rejected\n");
534         delete track; continue; 
535       }
536     }
537
538     // select tracks with d0rphi < maxd0rphi
539     if(optImpParCut>0 && TMath::Abs(d0z0[0]) > maxd0rphi) { 
540       if(fDebug) printf("     rejected\n");
541       delete track; continue; 
542     }
543
544     fTrkArray.AddLast(track);
545     nTrks++; 
546   }
547
548   delete initVertex;
549
550   return nTrks;
551
552 //---------------------------------------------------------------------------
553 AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
554                                                         TTree *trksTree,
555                                                         Float_t *diamondxy) 
556 {
557 //
558 // Removes tracks in trksTree from fit of inVtx
559 //
560
561   if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
562     printf("WARNING: result of tracks' removal will be only approximately correct\n");
563
564   TMatrixD rv(3,1);
565   rv(0,0) = inVtx->GetXv();
566   rv(1,0) = inVtx->GetYv();
567   rv(2,0) = inVtx->GetZv();
568   TMatrixD vV(3,3);
569   Double_t cov[6];
570   inVtx->GetCovMatrix(cov);
571   vV(0,0) = cov[0];
572   vV(0,1) = cov[1]; vV(1,0) = cov[1];
573   vV(1,1) = cov[2];
574   vV(0,2) = cov[3]; vV(2,0) = cov[3];
575   vV(1,2) = cov[4]; vV(2,1) = cov[4]; 
576   vV(2,2) = cov[5];
577
578   TMatrixD sumWi(TMatrixD::kInverted,vV);
579   TMatrixD sumWiri(sumWi,TMatrixD::kMult,rv);
580
581   Int_t nUsedTrks = inVtx->GetNContributors();
582   Double_t chi2 = inVtx->GetChi2();
583
584   AliESDtrack *track = 0;
585   trksTree->SetBranchAddress("tracks",&track);
586   Int_t ntrks = trksTree->GetEntries();
587   for(Int_t i=0;i<ntrks;i++) {
588     trksTree->GetEvent(i);
589     if(!inVtx->UsesTrack(track->GetID())) {
590       printf("track %d was not used in vertex fit\n",track->GetID());
591       continue;
592     }
593     Double_t alpha = track->GetAlpha();
594     Double_t xl = diamondxy[0]*TMath::Cos(alpha)+diamondxy[1]*TMath::Sin(alpha);
595     track->AliExternalTrackParam::PropagateTo(xl,GetFieldkG()); 
596     // vector of track global coordinates
597     TMatrixD ri(3,1);
598     // covariance matrix of ri
599     TMatrixD wWi(3,3);
600     
601     // get space point from track
602     if(!TrackToPoint(track,ri,wWi)) continue;
603
604     TMatrixD wWiri(wWi,TMatrixD::kMult,ri); 
605
606     sumWi -= wWi;
607     sumWiri -= wWiri;
608
609     // track chi2
610     TMatrixD deltar = rv; deltar -= ri;
611     TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
612     Double_t chi2i = deltar(0,0)*wWideltar(0,0)+
613                      deltar(1,0)*wWideltar(1,0)+
614                      deltar(2,0)*wWideltar(2,0);
615     // remove from total chi2
616     chi2 -= chi2i;
617
618     nUsedTrks--;
619     if(nUsedTrks<2) {
620       printf("Trying to remove too many tracks!\n");
621       return 0x0;
622     }
623   }
624
625   TMatrixD rvnew(3,1);
626   TMatrixD vVnew(3,3);
627
628   // new inverted of weights matrix
629   TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
630   vVnew = invsumWi;
631   // new position of primary vertex
632   rvnew.Mult(vVnew,sumWiri);
633
634   Double_t position[3];
635   position[0] = rvnew(0,0);
636   position[1] = rvnew(1,0);
637   position[2] = rvnew(2,0);
638   cov[0] = vVnew(0,0);
639   cov[1] = vVnew(0,1);
640   cov[2] = vVnew(1,1);
641   cov[3] = vVnew(0,2);
642   cov[4] = vVnew(1,2);
643   cov[5] = vVnew(2,2);
644   
645   // store data in the vertex object
646   AliESDVertex *outVtx = new AliESDVertex(position,cov,chi2,nUsedTrks);
647   outVtx->SetTitle(inVtx->GetTitle());
648   Double_t tp[3];
649   inVtx->GetTruePos(tp);
650   outVtx->SetTruePos(tp);
651   UShort_t *inindices = inVtx->GetIndices();
652   UShort_t *outindices = new UShort_t[outVtx->GetNContributors()];
653   Int_t j=0;
654   Bool_t copyindex;
655   for(Int_t k=0; k<inVtx->GetNIndices(); k++) {
656     copyindex=kTRUE;
657     for(Int_t l=0; l<ntrks; l++) {
658       trksTree->GetEvent(l);
659       if(inindices[k]==track->GetID()) copyindex=kFALSE;
660     }
661     if(copyindex) {
662       outindices[j] = inindices[k]; j++;
663     }
664   }
665   outVtx->SetIndices(outVtx->GetNContributors(),outindices);
666   delete [] outindices;
667
668   if(fDebug) {
669     printf("Vertex before removing tracks:\n");
670     inVtx->PrintStatus();
671     inVtx->PrintIndices();
672     printf("Vertex after removing tracks:\n");
673     outVtx->PrintStatus();
674     outVtx->PrintIndices();
675   }
676
677   return outVtx;
678 }
679 //---------------------------------------------------------------------------
680 void AliVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped) 
681 {
682 //
683 // Mark the tracks not to be used in the vertex reconstruction.
684 // Tracks are identified by AliESDtrack::GetID()
685 //
686   fNTrksToSkip = n;  fTrksToSkip = new Int_t[n]; 
687   for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i]; 
688   return; 
689 }
690 //---------------------------------------------------------------------------
691 void  AliVertexerTracks::SetVtxStart(AliESDVertex *vtx) 
692
693 //
694 // Set initial vertex knowledge
695 //
696   vtx->GetXYZ(fNominalPos);
697   vtx->GetCovMatrix(fNominalCov);
698   SetConstraintOn();
699   return; 
700 }
701 //---------------------------------------------------------------------------
702 void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
703 {
704   AliESDtrack *track1;
705   Double_t field=GetFieldkG();
706   const Int_t knacc = (Int_t)fTrkArray.GetEntriesFast();
707   TClonesArray *linarray = new TClonesArray("AliStrLine",1000);
708   TClonesArray &lines = *linarray;
709   for(Int_t i=0; i<knacc; i++){
710     track1 = (AliESDtrack*)fTrkArray.At(i);
711     Double_t alpha=track1->GetAlpha();
712     Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
713     Double_t pos[3],dir[3],sigmasq[3]; 
714     track1->GetXYZAt(mindist,field,pos);
715     track1->GetPxPyPzAt(mindist,field,dir);
716     sigmasq[0]=TMath::Sin(alpha)*TMath::Sin(alpha)*track1->GetSigmaY2();
717     sigmasq[1]=TMath::Cos(alpha)*TMath::Cos(alpha)*track1->GetSigmaY2();
718     sigmasq[2]=track1->GetSigmaZ2();
719     new(lines[i]) AliStrLine(pos,sigmasq,dir); 
720   }
721   fVert=TrackletVertexFinder(linarray,optUseWeights);
722   linarray->Delete();
723   delete linarray;
724 }
725 //---------------------------------------------------------------------------
726 AliVertex AliVertexerTracks::TrackletVertexFinder(TClonesArray *lines, Int_t optUseWeights)
727 {
728   // Calculate the point at minimum distance to prepared tracks 
729   
730   const Int_t knacc = (Int_t)lines->GetEntriesFast();
731   Double_t initPos[3]={0.,0.,0.};
732
733   Double_t (*vectP0)[3]=new Double_t [knacc][3];
734   Double_t (*vectP1)[3]=new Double_t [knacc][3];
735   
736   Double_t sum[3][3];
737   Double_t dsum[3]={0,0,0};
738   for(Int_t i=0;i<3;i++)
739     for(Int_t j=0;j<3;j++)sum[i][j]=0;
740   for(Int_t i=0; i<knacc; i++){
741     AliStrLine* line1 = (AliStrLine*)lines->At(i); 
742     Double_t p0[3],cd[3],sigmasq[3];
743     line1->GetP0(p0);
744     line1->GetCd(cd);
745     line1->GetSigma2P0(sigmasq);
746     Double_t p1[3]={p0[0]+cd[0],p0[1]+cd[1],p0[2]+cd[2]};
747     vectP0[i][0]=p0[0];
748     vectP0[i][1]=p0[1];
749     vectP0[i][2]=p0[2];
750     vectP1[i][0]=p1[0];
751     vectP1[i][1]=p1[1];
752     vectP1[i][2]=p1[2];
753     
754     Double_t matr[3][3];
755     Double_t dknow[3];
756     if(optUseWeights==0)GetStrLinDerivMatrix(p0,p1,matr,dknow);
757     else GetStrLinDerivMatrix(p0,p1,sigmasq,matr,dknow);
758
759
760     for(Int_t iii=0;iii<3;iii++){
761       dsum[iii]+=dknow[iii]; 
762       for(Int_t lj=0;lj<3;lj++) sum[iii][lj]+=matr[iii][lj];
763     }
764   }
765   
766   Double_t vett[3][3];
767   Double_t det=GetDeterminant3X3(sum);
768   Double_t sigma=0;
769   
770   if(det!=0){
771     for(Int_t zz=0;zz<3;zz++){
772       for(Int_t ww=0;ww<3;ww++){
773         for(Int_t kk=0;kk<3;kk++) vett[ww][kk]=sum[ww][kk];
774       }
775       for(Int_t kk=0;kk<3;kk++) vett[kk][zz]=dsum[kk];
776       initPos[zz]=GetDeterminant3X3(vett)/det;
777     }
778
779
780     for(Int_t i=0; i<knacc; i++){
781       Double_t p0[3]={0,0,0},p1[3]={0,0,0};
782       for(Int_t ii=0;ii<3;ii++){
783         p0[ii]=vectP0[i][ii];
784         p1[ii]=vectP1[i][ii];
785       }
786       sigma+=GetStrLinMinDist(p0,p1,initPos);
787     }
788
789     sigma=TMath::Sqrt(sigma);
790   }else{
791     sigma=999;
792   }
793   AliVertex theVert(initPos,sigma,knacc);
794   delete vectP0;
795   delete vectP1;
796   return theVert;
797 }
798 //---------------------------------------------------------------------------
799 Bool_t AliVertexerTracks::TrackToPoint(AliESDtrack *t,
800                                        TMatrixD &ri,TMatrixD &wWi) const 
801 {
802 //
803 // Extract from the AliESDtrack the global coordinates ri and covariance matrix
804 // wWi of the space point that it represents (to be used in VertexFitter())
805 //
806
807   
808   Double_t rotAngle = t->GetAlpha();
809   if(rotAngle<0.) rotAngle += 2.*TMath::Pi();
810   Double_t cosRot = TMath::Cos(rotAngle);
811   Double_t sinRot = TMath::Sin(rotAngle);
812
813   ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot;
814   ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot;
815   ri(2,0) = t->GetZ();
816
817   // matrix to go from global (x,y,z) to local (y,z);
818   TMatrixD qQi(2,3);
819   qQi(0,0) = -sinRot;
820   qQi(0,1) = cosRot;
821   qQi(0,2) = 0.;
822   qQi(1,0) = 0.;
823   qQi(1,1) = 0.;
824   qQi(1,2) = 1.;
825
826   // covariance matrix of local (y,z) - inverted
827   TMatrixD uUi(2,2);
828   Double_t cc[15];
829   t->GetExternalCovariance(cc);
830   uUi(0,0) = cc[0];
831   uUi(0,1) = cc[1];
832   uUi(1,0) = cc[1];
833   uUi(1,1) = cc[2];
834   if(uUi.Determinant() <= 0.) return kFALSE;
835   TMatrixD uUiInv(TMatrixD::kInverted,uUi);
836   
837   // weights matrix: wWi = qQiT * uUiInv * qQi
838   TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
839   TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
840   wWi = m;
841
842   return kTRUE;
843
844 //---------------------------------------------------------------------------
845 void AliVertexerTracks::TooFewTracks(const AliESD* esdEvent) 
846 {
847 //
848 // When the number of tracks is < fMinTracks
849 //
850
851   // deal with vertices not found
852   Double_t pos[3],err[3];
853   Int_t    ncontr=0;
854   pos[0] = fNominalPos[0];
855   err[0] = TMath::Sqrt(fNominalCov[0]);
856   pos[1] = fNominalPos[1];
857   err[1] = TMath::Sqrt(fNominalCov[2]);
858   pos[2] = esdEvent->GetVertex()->GetZv();
859   err[2] = esdEvent->GetVertex()->GetZRes();
860   if(err[0]>1. && esdEvent->GetVertex()->GetNContributors()<=0) 
861     ncontr = -1; // (x,y,z) = (0,0,0) 
862   if(err[0]>1. && esdEvent->GetVertex()->GetNContributors()>0) 
863     ncontr = -2; // (x,y,z) = (0,0,z_fromSPD) 
864   if(err[0]<1. && esdEvent->GetVertex()->GetNContributors()<=0) 
865     ncontr = -3; // (x,y,z) = (x_mean,y_mean,0)
866   if(err[0]<1. && esdEvent->GetVertex()->GetNContributors()>0) 
867     ncontr = -4; // (x,y,z) = (x_mean,y_mean,z_fromSPD)
868   fCurrentVertex = 0;
869   fCurrentVertex = new AliESDVertex(pos,err);
870   fCurrentVertex->SetNContributors(ncontr);
871
872   Double_t tp[3];
873   esdEvent->GetVertex()->GetTruePos(tp);
874   fCurrentVertex->SetTruePos(tp);
875   if(fConstraint) {
876     fCurrentVertex->SetTitle("VertexerTracksWithConstraint");
877   } else {
878     fCurrentVertex->SetTitle("VertexerTracksNoConstraint");
879   }
880
881   return;
882 }
883 //---------------------------------------------------------------------------
884 void AliVertexerTracks::VertexFinder(Int_t optUseWeights) 
885 {
886
887   // Get estimate of vertex position in (x,y) from tracks DCA
888  
889   Double_t initPos[3];
890   initPos[2] = 0.;
891   for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
892   Int_t nacc = (Int_t)fTrkArray.GetEntriesFast();
893   Double_t aver[3]={0.,0.,0.};
894   Double_t aversq[3]={0.,0.,0.};
895   Double_t sigmasq[3]={0.,0.,0.};
896   Double_t sigma=0;
897   Int_t ncombi = 0;
898   AliESDtrack *track1;
899   AliESDtrack *track2;
900   Double_t pos[3],dir[3]; 
901   Double_t alpha,mindist;
902   Double_t field=GetFieldkG();
903
904   for(Int_t i=0; i<nacc; i++){
905     track1 = (AliESDtrack*)fTrkArray.At(i);
906     alpha=track1->GetAlpha();
907     mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
908     track1->GetXYZAt(mindist,field,pos);
909     track1->GetPxPyPzAt(mindist,field,dir);
910     AliStrLine *line1 = new AliStrLine(pos,dir); 
911
912    //    AliStrLine *line1 = new AliStrLine();
913    //    track1->ApproximateHelixWithLine(mindist,field,line1);
914    
915     for(Int_t j=i+1; j<nacc; j++){
916       track2 = (AliESDtrack*)fTrkArray.At(j);
917       alpha=track2->GetAlpha();
918       mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
919       track2->GetXYZAt(mindist,field,pos);
920       track2->GetPxPyPzAt(mindist,field,dir);
921       AliStrLine *line2 = new AliStrLine(pos,dir); 
922     //      AliStrLine *line2 = new AliStrLine();
923     //  track2->ApproximateHelixWithLine(mindist,field,line2);
924       Double_t distCA=line2->GetDCA(line1);
925       //printf("%d   %d   %f\n",i,j,distCA);
926        if(fDCAcut<=0 || (fDCAcut>0&&distCA<fDCAcut)){
927         Double_t pnt1[3],pnt2[3],crosspoint[3];
928
929         if(optUseWeights<=0){
930           Int_t retcode = line2->Cross(line1,crosspoint);
931           if(retcode>=0){
932             ncombi++;
933             for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
934             for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
935           }
936         }
937         if(optUseWeights>0){
938           Int_t retcode = line1->CrossPoints(line2,pnt1,pnt2);
939           if(retcode>=0){
940             Double_t cs, sn;
941             alpha=track1->GetAlpha();
942             cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);   
943             Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
944             alpha=track2->GetAlpha();
945             cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
946             Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
947             Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
948             Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
949             Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
950             Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
951             crosspoint[0]=wx1*pnt1[0] + wx2*pnt2[0]; 
952             crosspoint[1]=wy1*pnt1[1] + wy2*pnt2[1]; 
953             crosspoint[2]=wz1*pnt1[2] + wz2*pnt2[2];
954           
955             ncombi++;
956             for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
957             for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
958           }
959         }
960       }
961       delete line2;
962     }
963     delete line1;
964   }
965   if(ncombi>0){
966     for(Int_t jj=0;jj<3;jj++){
967       initPos[jj] = aver[jj]/ncombi;
968       //printf("%f\n",initPos[jj]);
969       aversq[jj]/=ncombi;
970       sigmasq[jj]=aversq[jj]-initPos[jj]*initPos[jj];
971       sigma+=sigmasq[jj];
972     }
973     sigma=TMath::Sqrt(TMath::Abs(sigma));
974   }
975   else {
976     Warning("VertexFinder","Finder did not succed");
977     sigma=999;
978   }
979   fVert.SetXYZ(initPos);
980   fVert.SetDispersion(sigma);
981   fVert.SetNContributors(ncombi);
982 }
983 //---------------------------------------------------------------------------
984 void AliVertexerTracks::VertexFitter(Bool_t useConstraint) 
985 {
986 //
987 // The optimal estimate of the vertex position is given by a "weighted 
988 // average of tracks positions"
989 // Original method: V. Karimaki, CMS Note 97/0051
990 //
991   Double_t initPos[3];
992   fVert.GetXYZ(initPos);
993   Int_t arrEntries = (Int_t)fTrkArray.GetEntries();
994   if(arrEntries==1) useConstraint=kTRUE;
995   if(fDebug) { 
996     printf(" VertexFitter(): start\n");
997     printf(" Number of tracks in array: %d\n",(Int_t)fTrkArray.GetEntriesFast());
998     printf(" Minimum # tracks required in fit: %d\n",fMinTracks);
999     printf("Vertex position after finder: %f,%f,%f\n",initPos[0],initPos[1],initPos[2]);
1000     if(useConstraint) printf(" This vertex will be used in fit: (%f+-%f,%f+-%f)\n",fNominalPos[0],TMath::Sqrt(fNominalCov[0]),fNominalPos[1],TMath::Sqrt(fNominalCov[2])); 
1001   }
1002
1003   Int_t i,j,k,step=0;
1004   TMatrixD rv(3,1);
1005   TMatrixD vV(3,3);
1006   rv(0,0) = initPos[0];
1007   rv(1,0) = initPos[1];
1008   rv(2,0) = 0.;
1009   Double_t xlStart,alpha;
1010   Int_t nUsedTrks;
1011   Double_t chi2,chi2i,chi2b;
1012   AliESDtrack *t = 0;
1013   Int_t failed = 0;
1014
1015   // initial vertex covariance matrix
1016   TMatrixD vVb(3,3);
1017   vVb(0,0) = fNominalCov[0];
1018   vVb(0,1) = fNominalCov[1];
1019   vVb(0,2) = 0.;
1020   vVb(1,0) = fNominalCov[1];
1021   vVb(1,1) = fNominalCov[2];
1022   vVb(1,2) = 0.;
1023   vVb(2,0) = 0.;
1024   vVb(2,1) = 0.;
1025   vVb(2,2) = fNominalCov[5];
1026   TMatrixD vVbInv(TMatrixD::kInverted,vVb);
1027   TMatrixD rb(3,1);
1028   rb(0,0) = fNominalPos[0];
1029   rb(1,0) = fNominalPos[1];
1030   rb(2,0) = fNominalPos[2];
1031   TMatrixD vVbInvrb(vVbInv,TMatrixD::kMult,rb);
1032
1033
1034   // 2 steps:
1035   // 1st - estimate of vtx using all tracks
1036   // 2nd - estimate of global chi2
1037   for(step=0; step<2; step++) {
1038     if(fDebug) printf(" step = %d\n",step);
1039     chi2 = 0.;
1040     nUsedTrks = 0;
1041
1042     TMatrixD sumWiri(3,1);
1043     TMatrixD sumWi(3,3);
1044     for(i=0; i<3; i++) {
1045       sumWiri(i,0) = 0.;
1046       for(j=0; j<3; j++) sumWi(j,i) = 0.;
1047     }
1048
1049     // mean vertex constraint
1050     if(useConstraint) {
1051       for(i=0;i<3;i++) {
1052         sumWiri(i,0) += vVbInvrb(i,0);
1053         for(k=0;k<3;k++) sumWi(i,k) += vVbInv(i,k);
1054       }
1055       // chi2
1056       TMatrixD deltar = rv; deltar -= rb;
1057       TMatrixD vVbInvdeltar(vVbInv,TMatrixD::kMult,deltar);
1058       chi2b = deltar(0,0)*vVbInvdeltar(0,0)+
1059               deltar(1,0)*vVbInvdeltar(1,0)+
1060               deltar(2,0)*vVbInvdeltar(2,0);
1061       chi2 += chi2b;
1062     }
1063
1064
1065     // loop on tracks  
1066     for(k=0; k<arrEntries; k++) {
1067       // get track from track array
1068       t = (AliESDtrack*)fTrkArray.At(k);
1069       alpha = t->GetAlpha();
1070       xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha);
1071       // to vtxSeed (from finder)
1072       t->AliExternalTrackParam::PropagateTo(xlStart,GetFieldkG());   
1073  
1074
1075       // vector of track global coordinates
1076       TMatrixD ri(3,1);
1077       // covariance matrix of ri
1078       TMatrixD wWi(3,3);
1079
1080       // get space point from track
1081       if(!TrackToPoint(t,ri,wWi)) continue;
1082
1083       TMatrixD wWiri(wWi,TMatrixD::kMult,ri); 
1084
1085       // track chi2
1086       TMatrixD deltar = rv; deltar -= ri;
1087       TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
1088       chi2i = deltar(0,0)*wWideltar(0,0)+
1089               deltar(1,0)*wWideltar(1,0)+
1090               deltar(2,0)*wWideltar(2,0);
1091
1092       // add to total chi2
1093       chi2 += chi2i;
1094
1095       sumWiri += wWiri;
1096       sumWi   += wWi;
1097
1098       nUsedTrks++;
1099     } // end loop on tracks
1100
1101     if(nUsedTrks < fMinTracks) {
1102       failed=1;
1103       continue;
1104     }
1105
1106     Double_t determinant = sumWi.Determinant();
1107     //cerr<<" determinant: "<<determinant<<endl;
1108     if(determinant < 100.)  { 
1109       printf("det(V) = 0\n");       
1110       failed=1;
1111       continue;
1112     }
1113
1114     // inverted of weights matrix
1115     TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
1116     vV = invsumWi;
1117      
1118     // position of primary vertex
1119     rv.Mult(vV,sumWiri);
1120
1121   } // end loop on the 2 steps
1122
1123   // delete t;
1124
1125   if(failed) { 
1126     if(fDebug) printf("TooFewTracks\n");
1127     fCurrentVertex = new AliESDVertex(0.,0.,-1);
1128     return; 
1129   }
1130
1131   Double_t position[3];
1132   position[0] = rv(0,0);
1133   position[1] = rv(1,0);
1134   position[2] = rv(2,0);
1135   Double_t covmatrix[6];
1136   covmatrix[0] = vV(0,0);
1137   covmatrix[1] = vV(0,1);
1138   covmatrix[2] = vV(1,1);
1139   covmatrix[3] = vV(0,2);
1140   covmatrix[4] = vV(1,2);
1141   covmatrix[5] = vV(2,2);
1142   
1143   // for correct chi2/ndf, count constraint as additional "track"
1144   if(fConstraint) nUsedTrks++;
1145   // store data in the vertex object
1146   fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nUsedTrks);
1147
1148   if(fDebug) {
1149     printf(" VertexFitter(): finish\n");
1150     printf(" rv = ( %f , %f , %f )\n\n",rv(0,0),rv(1,0),rv(2,0));
1151     fCurrentVertex->PrintStatus();
1152   }
1153
1154   return;
1155 }
1156 //----------------------------------------------------------------------------
1157 AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(TTree *trkTree,Bool_t optUseFitter,Bool_t optPropagate) 
1158 {
1159 //
1160 // Return vertex from tracks in trkTree
1161 //
1162
1163   SetConstraintOff();
1164
1165   // get tracks and propagate them to initial vertex position
1166   Int_t nTrksPrep = PrepareTracks(*trkTree,0);
1167   if(nTrksPrep <  TMath::Max(2,fMinTracks) ) {
1168     if(fDebug) printf("TooFewTracks\n");
1169     fCurrentVertex = new AliESDVertex(0.,0.,-1);
1170     return fCurrentVertex;
1171   }
1172  
1173   switch (fAlgo) {
1174     case 1: StrLinVertexFinderMinDist(1); break;
1175     case 2: StrLinVertexFinderMinDist(0); break;
1176     case 3: HelixVertexFinder();          break;
1177     case 4: VertexFinder(1);              break;
1178     case 5: VertexFinder(0);              break;
1179     default: printf("Wrong algorithm\n"); break;  
1180   }
1181   
1182   if(fDebug) printf(" Vertex finding completed\n");
1183
1184   // vertex fitter
1185   if(optUseFitter){
1186     VertexFitter(fConstraint);
1187     if(fDebug) printf(" Vertex fit completed\n");
1188   }else{
1189     Double_t position[3]={fVert.GetXv(),fVert.GetYv(),fVert.GetZv()};
1190     Double_t covmatrix[6]={0.,0.,0.,0.,0.,0.};
1191     Double_t chi2=99999.;
1192     Int_t    nUsedTrks=fVert.GetNContributors();
1193     fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nUsedTrks);    
1194   }
1195   fCurrentVertex->SetDispersion(fVert.GetDispersion());
1196
1197
1198   // set indices of used tracks and propagate track to found vertex
1199   UShort_t *indices = 0;
1200   AliESDtrack *eta = 0;
1201   if(fCurrentVertex->GetNContributors()>0) {
1202     indices = new UShort_t[fCurrentVertex->GetNContributors()];
1203     for(Int_t jj=0;jj<(Int_t)fTrkArray.GetEntriesFast();jj++) {
1204       eta = (AliESDtrack*)fTrkArray.At(jj);
1205       indices[jj] = (UShort_t)eta->GetID();
1206       if(optPropagate&&optUseFitter){
1207         if(TMath::Sqrt(fCurrentVertex->GetXv()*fCurrentVertex->GetXv()+fCurrentVertex->GetYv()*fCurrentVertex->GetYv())<3.) {
1208           eta->RelateToVertex(fCurrentVertex,GetFieldkG(),100.);
1209           if(fDebug) printf("Track %d propagated to found vertex\n",jj);
1210         }else{
1211           AliWarning("Found vertex outside beam pipe!");
1212         }
1213       }
1214     }
1215     fCurrentVertex->SetIndices(fCurrentVertex->GetNContributors(),indices);
1216   }
1217   delete [] indices;
1218   fTrkArray.Delete();
1219   
1220   return fCurrentVertex;
1221 }
1222 //----------------------------------------------------------------------------
1223 AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray,Bool_t optUseFitter, Bool_t optPropagate) 
1224 {
1225 //
1226 // Return vertex from array of tracks
1227 //
1228
1229   // get tracks and propagate them to initial vertex position
1230   Int_t nTrks = trkArray->GetEntriesFast();
1231   if(nTrks < TMath::Max(2,fMinTracks) ) {
1232     if(fDebug) printf("TooFewTracks\n");
1233     fCurrentVertex = new AliESDVertex(0.,0.,-1);
1234     return fCurrentVertex;
1235   }
1236   TTree *trkTree = new TTree("TreeT","tracks");
1237   AliESDtrack *esdTrack = 0;
1238   trkTree->Branch("tracks","AliESDtrack",&esdTrack);
1239   for(Int_t i=0; i<nTrks; i++){
1240     esdTrack = (AliESDtrack*)trkArray->At(i);
1241     trkTree->Fill();
1242   }
1243     
1244   AliESDVertex *vtx =  VertexForSelectedTracks(trkTree,optUseFitter,optPropagate);
1245   delete trkTree;
1246   return vtx;
1247 }
1248 //--------------------------------------------------------------------------