73d450a46fb827a15f83064fa0e349731bef5dde
[u/mrichter/AliRoot.git] / STEER / ESD / AliVertexerTracks.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 //-----------------------------------------------------------------
17 //    Implementation of the vertexer from tracks
18 //
19 // Origin: AliITSVertexerTracks
20 //         A.Dainese, Padova, 
21 //         andrea.dainese@pd.infn.it
22 //         M.Masera,  Torino, 
23 //         massimo.masera@to.infn.it
24 // Moved to STEER and adapted to ESD tracks: 
25 //         F.Prino,  Torino, prino@to.infn.it
26 //-----------------------------------------------------------------
27
28 //---- Root headers --------
29 #include <TSystem.h>
30 #include <TClonesArray.h>
31 #include <TDirectory.h>
32 #include <TFile.h>
33 //---- AliRoot headers -----
34 #include "AliStrLine.h"
35 #include "AliExternalTrackParam.h"
36 #include "AliNeutralTrackParam.h"
37 #include "AliVEvent.h"
38 #include "AliVTrack.h"
39 #include "AliESDtrack.h"
40 #include "AliESDEvent.h"
41 #include "AliVertexerTracks.h"
42
43
44 ClassImp(AliVertexerTracks)
45
46
47 //----------------------------------------------------------------------------
48 AliVertexerTracks::AliVertexerTracks():
49 TObject(),
50 fVert(),
51 fCurrentVertex(0),
52 fMode(0),
53 fFieldkG(-999.),
54 fTrkArraySel(),
55 fIdSel(0),
56 fTrksToSkip(0),
57 fNTrksToSkip(0),
58 fConstraint(kFALSE),
59 fOnlyFitter(kFALSE),
60 fMinTracks(1),
61 fMinClusters(5),
62 fDCAcut(0.1),
63 fDCAcutIter0(0.1),
64 fNSigma(3.),
65 fMaxd0z0(0.5),
66 fMinDetFitter(100.),
67 fMaxTgl(1000.),
68 fITSrefit(kTRUE),
69 fITSpureSA(kFALSE),
70 fFiducialR(3.),
71 fFiducialZ(30.),
72 fnSigmaForUi00(1.5),
73 fAlgo(1),
74 fAlgoIter0(4),
75 fSelectOnTOFBunchCrossing(kFALSE), 
76 fKeepAlsoUnflaggedTOFBunchCrossing(kTRUE),
77 fMVWSum(0),
78 fMVWE2(0),
79 fMVTukey2(6.),
80 fMVSigma2(1.),
81 fMVSig2Ini(1.e3),
82 fMVMaxSigma2(3.),
83 fMVMinSig2Red(0.005),
84 fMVMinDst(10.e-4),
85 fMVScanStep(3.),
86 fMVMaxWghNtr(10.),
87 fMVFinalWBinary(kTRUE),
88 fBCSpacing(50),
89 fMVVertices(0)
90 {
91 //
92 // Default constructor
93 //
94   SetVtxStart();
95   SetVtxStartSigma();
96 }
97 //----------------------------------------------------------------------------
98 AliVertexerTracks::AliVertexerTracks(Double_t fieldkG):
99 TObject(),
100 fVert(),
101 fCurrentVertex(0),
102 fMode(0),
103 fFieldkG(fieldkG),
104 fTrkArraySel(),
105 fIdSel(0),
106 fTrksToSkip(0),
107 fNTrksToSkip(0),
108 fConstraint(kFALSE),
109 fOnlyFitter(kFALSE),
110 fMinTracks(1),
111 fMinClusters(5),
112 fDCAcut(0.1),
113 fDCAcutIter0(0.1),
114 fNSigma(3.),
115 fMaxd0z0(0.5),
116 fMinDetFitter(100.),
117 fMaxTgl(1000.),
118 fITSrefit(kTRUE),
119 fITSpureSA(kFALSE),
120 fFiducialR(3.),
121 fFiducialZ(30.),
122 fnSigmaForUi00(1.5),
123 fAlgo(1),
124 fAlgoIter0(4),
125 fSelectOnTOFBunchCrossing(kFALSE), 
126 fKeepAlsoUnflaggedTOFBunchCrossing(kTRUE),
127 fMVWSum(0),
128 fMVWE2(0),
129 fMVTukey2(6.),
130 fMVSigma2(1.),
131 fMVSig2Ini(1.e3),
132 fMVMaxSigma2(3.),
133 fMVMinSig2Red(0.005),
134 fMVMinDst(10.e-4),
135 fMVScanStep(3.),
136 fMVMaxWghNtr(10.),
137 fMVFinalWBinary(kTRUE),
138 fBCSpacing(50),
139 fMVVertices(0)
140 {
141 //
142 // Standard constructor
143 //
144   SetVtxStart();
145   SetVtxStartSigma();
146 }
147 //-----------------------------------------------------------------------------
148 AliVertexerTracks::~AliVertexerTracks() 
149 {
150   // Default Destructor
151   // The objects pointed by the following pointer are not owned
152   // by this class and are not deleted
153   fCurrentVertex = 0;
154   if(fTrksToSkip) { delete [] fTrksToSkip; fTrksToSkip=NULL; fNTrksToSkip=0;}
155   if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
156   if(fMVVertices) delete fMVVertices;
157 }
158
159 //----------------------------------------------------------------------------
160 AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliVEvent *vEvent)
161 {
162 //
163 // Primary vertex for current ESD or AOD event
164 // (Two iterations: 
165 //  1st with 5*fNSigma*sigma cut w.r.t. to initial vertex
166 //      + cut on sqrt(d0d0+z0z0) if fConstraint=kFALSE  
167 //  2nd with fNSigma*sigma cut w.r.t. to vertex found in 1st iteration) 
168 //
169   fCurrentVertex = 0;
170   TString evtype = vEvent->IsA()->GetName();
171   Bool_t inputAOD = ((evtype=="AliAODEvent") ? kTRUE : kFALSE);
172
173   if(inputAOD && fMode==1) {
174     printf("Error : AliVertexerTracks: no TPC-only vertex from AOD\n"); 
175     TooFewTracks(); 
176     return fCurrentVertex;
177   }
178
179   // accept 1-track case only if constraint is available
180   if(!fConstraint && fMinTracks==1) fMinTracks=2;
181
182   // read tracks from AlivEvent
183   Int_t nTrks = (Int_t)vEvent->GetNumberOfTracks();
184   if(nTrks<fMinTracks) {
185     TooFewTracks();
186     return fCurrentVertex;
187   } 
188   //
189   int bcRound = fBCSpacing/25;   // profit from larger than 25ns spacing and set correct BC
190   TDirectory * olddir = gDirectory;
191   TFile *f = 0;
192   if(nTrks>500) f = new TFile("VertexerTracks.root","recreate");
193   TObjArray trkArrayOrig(nTrks);
194   UShort_t *idOrig = new UShort_t[nTrks];
195
196   Int_t nTrksOrig=0;
197   AliExternalTrackParam *t=0;
198   // loop on tracks
199   for(Int_t i=0; i<nTrks; i++) {
200     AliVTrack *track = (AliVTrack*)vEvent->GetTrack(i);
201     // check tracks to skip
202     Bool_t skipThis = kFALSE;
203     for(Int_t j=0; j<fNTrksToSkip; j++) { 
204       if(track->GetID()==fTrksToSkip[j]) {
205         AliDebug(1,Form("skipping track: %d",i));
206         skipThis = kTRUE;
207       }
208     }
209     if(skipThis) continue;
210
211     // skip pure ITS SA tracks (default)
212     if(!fITSpureSA && (track->GetStatus()&AliESDtrack::kITSpureSA)) continue;
213     // or use only pure ITS SA tracks
214     if(fITSpureSA && !(track->GetStatus()&AliESDtrack::kITSpureSA)) continue;
215
216     // kITSrefit
217     if(fMode==0 && fITSrefit && !(track->GetStatus()&AliESDtrack::kITSrefit)) continue;
218
219     if(!inputAOD) {  // ESD
220       AliESDtrack* esdt = (AliESDtrack*)track;
221       if(esdt->GetNcls(fMode) < fMinClusters) continue;
222       if(fMode==0) {        // ITS mode
223         Double_t x,p[5],cov[15];
224         esdt->GetExternalParameters(x,p);
225         esdt->GetExternalCovariance(cov);
226         t = new AliExternalTrackParam(x,esdt->GetAlpha(),p,cov);
227       } else if(fMode==1) { // TPC mode
228         t = (AliExternalTrackParam*)esdt->GetTPCInnerParam();
229         if(!t) continue;
230         Double_t radius = 2.8; //something less than the beam pipe radius
231         if(!PropagateTrackTo(t,radius)) continue;
232       }
233     } else {          // AOD (only ITS mode)
234       Int_t ncls0=0;
235       for(Int_t l=0;l<6;l++) if(TESTBIT(track->GetITSClusterMap(),l)) ncls0++;
236       if(ncls0 < fMinClusters) continue;
237       t = new AliExternalTrackParam(); t->CopyFromVTrack(track);
238     }
239
240     // use TOF info about bunch crossing
241     if(fSelectOnTOFBunchCrossing) {
242       int bc = track->GetTOFBunchCrossing(fFieldkG);
243       if (bc>AliVTrack::kTOFBCNA) bc /= bcRound;
244       t->SetUniqueID(UInt_t(bc + kTOFBCShift));
245     }
246     //
247     trkArrayOrig.AddLast(t);
248     idOrig[nTrksOrig]=(UShort_t)track->GetID();
249     nTrksOrig++;
250   } // end loop on tracks
251   
252   // call method that will reconstruct the vertex
253   FindPrimaryVertex(&trkArrayOrig,idOrig);
254
255   if(!inputAOD) AnalyzePileUp((AliESDEvent*)vEvent);
256
257   if(fMode==0) trkArrayOrig.Delete();
258   delete [] idOrig; idOrig=NULL;
259
260   if(f) {
261     f->Close(); delete f; f = NULL;
262     gSystem->Unlink("VertexerTracks.root");
263     olddir->cd();
264   }
265
266   // set vertex ID for tracks used in the fit
267   // (only for ESD)
268   if(!inputAOD && fCurrentVertex) {
269     Int_t nIndices = fCurrentVertex->GetNIndices();
270     UShort_t *indices = fCurrentVertex->GetIndices();
271     for(Int_t ind=0; ind<nIndices; ind++) {
272       AliESDtrack *esdt = (AliESDtrack*)vEvent->GetTrack(indices[ind]);
273       esdt->SetVertexID(-1);
274     }
275   }
276
277   return fCurrentVertex;
278 }
279 //----------------------------------------------------------------------------
280 AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const TObjArray *trkArrayOrig,
281                                                    UShort_t *idOrig)
282 {
283 //
284 // Primary vertex using the AliExternalTrackParam's in the TObjArray.
285 // idOrig must contain the track IDs from AliESDtrack::GetID()
286 // (Two iterations: 
287 //  1st with 5*fNSigma*sigma cut w.r.t. to initial vertex
288 //      + cut on sqrt(d0d0+z0z0) if fConstraint=kFALSE  
289 //  2nd with fNSigma*sigma cut w.r.t. to vertex found in 1st iteration) 
290 //
291   fCurrentVertex = 0;
292   // accept 1-track case only if constraint is available
293   if(!fConstraint && fMinTracks==1) fMinTracks=2;
294
295   // read tracks from array
296   Int_t nTrksOrig = (Int_t)trkArrayOrig->GetEntriesFast();
297   AliDebug(1,Form("Initial number of tracks: %d",nTrksOrig));
298   if(nTrksOrig<fMinTracks) {
299     AliDebug(1,"TooFewTracks");
300     TooFewTracks();
301     return fCurrentVertex;
302   } 
303
304   // If fConstraint=kFALSE
305   // run VertexFinder(1) to get rough estimate of initVertex (x,y)
306   if(!fConstraint) {
307     // fill fTrkArraySel, for VertexFinder()
308     fIdSel = new UShort_t[nTrksOrig];
309     PrepareTracks(*trkArrayOrig,idOrig,0);
310     if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
311     Double_t cutsave = fDCAcut;  
312     fDCAcut = fDCAcutIter0;
313     // vertex finder
314     switch (fAlgoIter0) {
315     case 1: StrLinVertexFinderMinDist(1); break;
316     case 2: StrLinVertexFinderMinDist(0); break;
317     case 3: HelixVertexFinder();          break;
318     case 4: VertexFinder(1);              break;
319     case 5: VertexFinder(0);              break;
320     default: printf("Wrong algorithm\n"); break;  
321     }
322     fDCAcut = cutsave;
323     if(fVert.GetNContributors()>0) {
324       fVert.GetXYZ(fNominalPos);
325       fNominalPos[0] = fVert.GetXv();
326       fNominalPos[1] = fVert.GetYv();
327       fNominalPos[2] = fVert.GetZv();
328       AliDebug(1,Form("No mean vertex: VertexFinder gives (%f, %f, %f)",fNominalPos[0],fNominalPos[1],fNominalPos[2]));
329     } else {
330       fNominalPos[0] = 0.;
331       fNominalPos[1] = 0.;
332       fNominalPos[2] = 0.;
333       AliDebug(1,"No mean vertex and VertexFinder failed");
334     }
335   }
336   
337   // TWO ITERATIONS:
338   //
339   // ITERATION 1
340   // propagate tracks to fNominalPos vertex
341   // preselect them:
342   // if(constraint) reject for |d0|>5*fNSigma*sigma w.r.t. fNominal... vertex
343   // else  reject for |d0|\oplus|z0| > 5 mm w.r.t. fNominal... vertex
344   // ITERATION 2
345   // propagate tracks to best between initVertex and fCurrentVertex
346   // preselect tracks (reject for |d0|>fNSigma*sigma w.r.t. best 
347   //                   between initVertex and fCurrentVertex) 
348   Bool_t multiMode = kFALSE;
349   for(Int_t iter=1; iter<=2; iter++) {
350     if (fAlgo==kMultiVertexer && iter==2) break; // multivertexer does not need 2 iterations
351     if(fOnlyFitter && iter==1) continue; 
352     if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
353     fIdSel = new UShort_t[nTrksOrig];
354     Int_t nTrksSel = PrepareTracks(*trkArrayOrig,idOrig,iter);
355     AliDebug(1,Form("N tracks selected in iteration %d: %d",iter,nTrksSel));
356     if(nTrksSel < fMinTracks) {
357       TooFewTracks();
358       return fCurrentVertex; 
359     }
360
361     // vertex finder
362     if(!fOnlyFitter) {
363       if(nTrksSel==1 && fAlgo!=kMultiVertexer) {
364         AliDebug(1,"Just one track");
365         OneTrackVertFinder();
366       } else {
367         switch (fAlgo) {
368         case kStrLinVertexFinderMinDist1: StrLinVertexFinderMinDist(1); break;
369         case kStrLinVertexFinderMinDist0: StrLinVertexFinderMinDist(0); break;
370         case kHelixVertexFinder         : HelixVertexFinder();          break;
371         case kVertexFinder1             : VertexFinder(1);              break;
372         case kVertexFinder0             : VertexFinder(0);              break;
373         case kMultiVertexer             : FindVerticesMV(); multiMode = kTRUE; break;
374         default: printf("Wrong algorithm"); break;  
375         }
376       }
377       AliDebug(1," Vertex finding completed");
378     }
379     if (multiMode) break; // // multivertexer does not need 2nd iteration
380     // vertex fitter
381     VertexFitter();
382   } // end loop on the two iterations
383
384   if (!multiMode || fMVVertices->GetEntries()==0) { // in multi-vertex mode this is already done for found vertices
385     // set indices of used tracks
386     UShort_t *indices = 0;
387     if(fCurrentVertex->GetNContributors()>0) {
388       Int_t nIndices = (Int_t)fTrkArraySel.GetEntriesFast();
389       indices = new UShort_t[nIndices];
390       for(Int_t jj=0; jj<nIndices; jj++)
391         indices[jj] = fIdSel[jj];
392       fCurrentVertex->SetIndices(nIndices,indices);
393     }
394     if (indices) {delete [] indices; indices=NULL;}
395     //
396     // set vertex title
397     TString title="VertexerTracksNoConstraint";
398     if(fConstraint) {
399       title="VertexerTracksWithConstraint";
400       if(fOnlyFitter) title.Append("OnlyFitter");
401     }
402     fCurrentVertex->SetTitle(title.Data());
403     //
404     AliDebug(1,Form("xyz: %f %f %f; nc %d",fCurrentVertex->GetXv(),fCurrentVertex->GetYv(),fCurrentVertex->GetZv(),fCurrentVertex->GetNContributors()));
405   }
406   // clean up
407   delete [] fIdSel; fIdSel=NULL;
408   fTrkArraySel.Delete();
409   if(fTrksToSkip) { delete [] fTrksToSkip; fTrksToSkip=NULL; }
410   //
411   
412   return fCurrentVertex;
413 }
414 //------------------------------------------------------------------------
415 Double_t AliVertexerTracks::GetDeterminant3X3(Double_t matr[][3])
416 {
417   //
418   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];
419  return det;
420 }
421 //-------------------------------------------------------------------------
422 void AliVertexerTracks::GetStrLinDerivMatrix(const Double_t *p0,const Double_t *p1,Double_t (*m)[3],Double_t *d)
423 {
424   //
425   Double_t x12=p0[0]-p1[0];
426   Double_t y12=p0[1]-p1[1];
427   Double_t z12=p0[2]-p1[2];
428   Double_t kk=x12*x12+y12*y12+z12*z12;
429   m[0][0]=2-2/kk*x12*x12;
430   m[0][1]=-2/kk*x12*y12;
431   m[0][2]=-2/kk*x12*z12;
432   m[1][0]=-2/kk*x12*y12;
433   m[1][1]=2-2/kk*y12*y12;
434   m[1][2]=-2/kk*y12*z12;
435   m[2][0]=-2/kk*x12*z12;
436   m[2][1]=-2*y12*z12;
437   m[2][2]=2-2/kk*z12*z12;
438   d[0]=2*p0[0]-2/kk*p0[0]*x12*x12-2/kk*p0[2]*x12*z12-2/kk*p0[1]*x12*y12;
439   d[1]=2*p0[1]-2/kk*p0[1]*y12*y12-2/kk*p0[0]*x12*y12-2/kk*p0[2]*z12*y12;
440   d[2]=2*p0[2]-2/kk*p0[2]*z12*z12-2/kk*p0[0]*x12*z12-2/kk*p0[1]*z12*y12;
441
442 }
443 //--------------------------------------------------------------------------  
444 void AliVertexerTracks::GetStrLinDerivMatrix(const Double_t *p0,const Double_t *p1,const Double_t *sigmasq,Double_t (*m)[3],Double_t *d)
445 {
446   //
447   Double_t x12=p1[0]-p0[0];
448   Double_t y12=p1[1]-p0[1];
449   Double_t z12=p1[2]-p0[2];
450
451   Double_t den= x12*x12*sigmasq[1]*sigmasq[2]+y12*y12*sigmasq[0]*sigmasq[2]+z12*z12*sigmasq[0]*sigmasq[1];
452
453   Double_t kk= 2*(x12*x12/sigmasq[0]+y12*y12/sigmasq[1]+z12*z12/sigmasq[2]);
454
455   Double_t cc[3];
456   cc[0]=-x12/sigmasq[0];
457   cc[1]=-y12/sigmasq[1];
458   cc[2]=-z12/sigmasq[2];
459
460   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;
461
462   Double_t ss= -p0[0]*cc[0]-p0[1]*cc[1]-p0[2]*cc[2];
463
464   Double_t aa[3];
465   aa[0]=x12*sigmasq[1]*sigmasq[2]/den;
466   aa[1]=y12*sigmasq[0]*sigmasq[2]/den;
467   aa[2]=z12*sigmasq[0]*sigmasq[1]/den;
468
469   m[0][0]=aa[0]*(aa[0]*kk+2*cc[0])+2*cc[0]*aa[0]+2/sigmasq[0];
470   m[0][1]=aa[1]*(aa[0]*kk+2*cc[0])+2*cc[1]*aa[0];
471   m[0][2]=aa[2]*(aa[0]*kk+2*cc[0])+2*cc[2]*aa[0];
472
473   m[1][0]=aa[0]*(aa[1]*kk+2*cc[1])+2*cc[0]*aa[1];
474   m[1][1]=aa[1]*(aa[1]*kk+2*cc[1])+2*cc[1]*aa[1]+2/sigmasq[1];
475   m[1][2]=aa[2]*(aa[1]*kk+2*cc[1])+2*cc[2]*aa[1];
476
477   m[2][0]=aa[0]*(aa[2]*kk+2*cc[2])+2*cc[0]*aa[2];
478   m[2][1]=aa[1]*(aa[2]*kk+2*cc[2])+2*cc[1]*aa[2];
479   m[2][2]=aa[2]*(aa[2]*kk+2*cc[2])+2*cc[2]*aa[2]+2/sigmasq[2];
480
481   d[0]=-ww*(aa[0]*kk+2*cc[0])-2*ss*aa[0]+2*p0[0]/sigmasq[0];
482   d[1]=-ww*(aa[1]*kk+2*cc[1])-2*ss*aa[1]+2*p0[1]/sigmasq[1];
483   d[2]=-ww*(aa[2]*kk+2*cc[2])-2*ss*aa[2]+2*p0[2]/sigmasq[2];
484
485   }
486 //--------------------------------------------------------------------------   
487 Double_t AliVertexerTracks::GetStrLinMinDist(const Double_t *p0,const Double_t *p1,const Double_t *x0)
488 {
489   //
490   Double_t x12=p0[0]-p1[0];
491   Double_t y12=p0[1]-p1[1];
492   Double_t z12=p0[2]-p1[2];
493   Double_t x10=p0[0]-x0[0];
494   Double_t y10=p0[1]-x0[1];
495   Double_t z10=p0[2]-x0[2];
496   //  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);
497   return ((y10*z12-z10*y12)*(y10*z12-z10*y12)+
498           (z10*x12-x10*z12)*(z10*x12-x10*z12)+
499           (x10*y12-y10*x12)*(x10*y12-y10*x12))
500     /(x12*x12+y12*y12+z12*z12);
501 }
502 //---------------------------------------------------------------------------
503 void AliVertexerTracks::OneTrackVertFinder() 
504 {
505   // find vertex for events with 1 track, using DCA to nominal beam axis
506   AliDebug(1,Form("Number of prepared tracks =%d - Call OneTrackVertFinder",fTrkArraySel.GetEntries()));
507   AliExternalTrackParam *track1;
508   track1 = (AliExternalTrackParam*)fTrkArraySel.At(0);
509   Double_t alpha=track1->GetAlpha();
510   Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
511   Double_t pos[3],dir[3]; 
512   track1->GetXYZAt(mindist,GetFieldkG(),pos);
513   track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
514   AliStrLine *line1 = new AliStrLine(pos,dir);
515   Double_t p1[3]={fNominalPos[0],fNominalPos[1],0.}; 
516   Double_t p2[3]={fNominalPos[0],fNominalPos[1],10.}; 
517   AliStrLine *zeta=new AliStrLine(p1,p2,kTRUE);
518   Double_t crosspoint[3]={0.,0.,0.};
519   Double_t sigma=999.;
520   Int_t nContrib=-1;
521   Int_t retcode = zeta->Cross(line1,crosspoint);
522   if(retcode>=0){
523     sigma=line1->GetDistFromPoint(crosspoint);
524     nContrib=1;
525   }
526   delete zeta;
527   delete line1;
528   fVert.SetXYZ(crosspoint);
529   fVert.SetDispersion(sigma);
530   fVert.SetNContributors(nContrib);  
531 }
532 //---------------------------------------------------------------------------
533 void AliVertexerTracks::HelixVertexFinder() 
534 {
535   // Get estimate of vertex position in (x,y) from tracks DCA
536
537
538   Double_t initPos[3];
539   initPos[2] = 0.;
540   for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
541
542   Int_t nacc = (Int_t)fTrkArraySel.GetEntriesFast();
543
544   Double_t aver[3]={0.,0.,0.};
545   Double_t averquad[3]={0.,0.,0.};
546   Double_t sigmaquad[3]={0.,0.,0.};
547   Double_t sigma=0;
548   Int_t ncombi = 0;
549   AliExternalTrackParam *track1;
550   AliExternalTrackParam *track2;
551   Double_t distCA;
552   Double_t x;
553   Double_t alpha, cs, sn;
554   Double_t crosspoint[3];
555   for(Int_t i=0; i<nacc; i++){
556     track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
557     
558
559     for(Int_t j=i+1; j<nacc; j++){
560       track2 = (AliExternalTrackParam*)fTrkArraySel.At(j);
561
562       distCA=track2->PropagateToDCA(track1,GetFieldkG());
563       if(fDCAcut<=0 ||(fDCAcut>0&&distCA<fDCAcut)){
564         x=track1->GetX();
565         alpha=track1->GetAlpha();
566         cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
567         Double_t x1=x*cs - track1->GetY()*sn;
568         Double_t y1=x*sn + track1->GetY()*cs;
569         Double_t z1=track1->GetZ();
570         
571         Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2(); 
572         x=track2->GetX();
573         alpha=track2->GetAlpha();
574         cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
575         Double_t x2=x*cs - track2->GetY()*sn;
576         Double_t y2=x*sn + track2->GetY()*cs;
577         Double_t z2=track2->GetZ();
578         Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
579         Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
580         Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
581         Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
582         Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
583         crosspoint[0]=wx1*x1 + wx2*x2; 
584         crosspoint[1]=wy1*y1 + wy2*y2; 
585         crosspoint[2]=wz1*z1 + wz2*z2;
586
587         ncombi++;
588         for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
589         for(Int_t jj=0;jj<3;jj++)averquad[jj]+=(crosspoint[jj]*crosspoint[jj]);
590       }
591     }
592       
593   }
594   if(ncombi>0){
595     for(Int_t jj=0;jj<3;jj++){
596       initPos[jj] = aver[jj]/ncombi;
597       averquad[jj]/=ncombi;
598       sigmaquad[jj]=averquad[jj]-initPos[jj]*initPos[jj];
599       sigma+=sigmaquad[jj];
600     }
601     sigma=TMath::Sqrt(TMath::Abs(sigma));
602   }
603   else {
604     Warning("HelixVertexFinder","Finder did not succed");
605     sigma=999;
606   }
607   fVert.SetXYZ(initPos);
608   fVert.SetDispersion(sigma);
609   fVert.SetNContributors(ncombi);
610 }
611 //----------------------------------------------------------------------------
612 Int_t AliVertexerTracks::PrepareTracks(const TObjArray &trkArrayOrig,
613                                        const UShort_t *idOrig,
614                                        Int_t optImpParCut) 
615 {
616 //
617 // Propagate tracks to initial vertex position and store them in a TObjArray
618 //
619   AliDebug(1," PrepareTracks()");
620
621   Int_t nTrksOrig = (Int_t)trkArrayOrig.GetEntriesFast();
622   Int_t nTrksSel = 0;
623   Double_t maxd0rphi; 
624   Double_t sigmaCurr[3];
625   Double_t normdistx,normdisty;
626   Double_t d0z0[2],covd0z0[3]; 
627   Double_t sigmad0;
628
629   AliESDVertex *initVertex = new AliESDVertex(fNominalPos,fNominalCov,1,1);
630
631   if(!fTrkArraySel.IsEmpty()) fTrkArraySel.Delete();
632
633   AliExternalTrackParam *track=0;
634
635   // loop on tracks
636   for(Int_t i=0; i<nTrksOrig; i++) {
637     AliExternalTrackParam *trackOrig=(AliExternalTrackParam*)trkArrayOrig.At(i);
638     if(trackOrig->Charge()!=0) { // normal tracks
639       track = new AliExternalTrackParam(*(AliExternalTrackParam*)trkArrayOrig.At(i));
640     } else { // neutral tracks (from a V0)
641       track = new AliNeutralTrackParam(*(AliNeutralTrackParam*)trkArrayOrig.At(i));
642     }
643     track->SetUniqueID(trackOrig->GetUniqueID());
644     // tgl cut
645     if(TMath::Abs(track->GetTgl())>fMaxTgl) {
646       AliDebug(1,Form(" rejecting track with tgl = %f",track->GetTgl()));
647       delete track; continue;
648     }
649
650     Bool_t propagateOK = kFALSE, cutond0z0 = kTRUE;
651     // propagate track to vertex
652     if(optImpParCut<2 || fOnlyFitter) { // optImpParCut==1 or 0
653       propagateOK = track->PropagateToDCA(initVertex,GetFieldkG(),100.,d0z0,covd0z0);
654     } else {              // optImpParCut==2
655       fCurrentVertex->GetSigmaXYZ(sigmaCurr);
656       normdistx = TMath::Abs(fCurrentVertex->GetXv()-fNominalPos[0])/TMath::Sqrt(sigmaCurr[0]*sigmaCurr[0]+fNominalCov[0]);
657       normdisty = TMath::Abs(fCurrentVertex->GetYv()-fNominalPos[1])/TMath::Sqrt(sigmaCurr[1]*sigmaCurr[1]+fNominalCov[2]);
658       AliDebug(1,Form("normdistx %f  %f    %f",fCurrentVertex->GetXv(),fNominalPos[0],TMath::Sqrt(sigmaCurr[0]*sigmaCurr[0]+fNominalCov[0])));
659       AliDebug(1,Form("normdisty %f  %f    %f",fCurrentVertex->GetYv(),fNominalPos[1],TMath::Sqrt(sigmaCurr[1]*sigmaCurr[1]+fNominalCov[2])));
660       AliDebug(1,Form("sigmaCurr %f %f    %f",sigmaCurr[0],sigmaCurr[1],TMath::Sqrt(fNominalCov[0])+TMath::Sqrt(fNominalCov[2])));
661       if(normdistx < 3. && normdisty < 3. &&
662          (sigmaCurr[0]+sigmaCurr[1])<(TMath::Sqrt(fNominalCov[0])+TMath::Sqrt(fNominalCov[2]))) {
663         propagateOK = track->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
664       } else {
665         propagateOK = track->PropagateToDCA(initVertex,GetFieldkG(),100.,d0z0,covd0z0);
666         if(fConstraint) cutond0z0=kFALSE;
667       }
668     }
669
670     if(!propagateOK) { 
671       AliDebug(1,"     rejected");
672       delete track; continue; 
673     }
674
675     sigmad0 = TMath::Sqrt(covd0z0[0]);
676     maxd0rphi = fNSigma*sigmad0;
677     if(optImpParCut==1) maxd0rphi *= 5.;
678     maxd0rphi = TMath::Min(maxd0rphi,fFiducialR); 
679     //sigmad0z0 = TMath::Sqrt(covd0z0[0]+covd0z0[2]);
680
681     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));
682
683
684     //---- track selection based on impact parameters ----//
685
686     // always reject tracks outside fiducial volume
687     if(TMath::Abs(d0z0[0])>fFiducialR || TMath::Abs(d0z0[1])>fFiducialZ) { 
688       AliDebug(1,"     rejected");
689       delete track; continue; 
690     }
691
692     // during iterations 1 and 2 , reject tracks with d0rphi > maxd0rphi
693     if(optImpParCut>0 && TMath::Abs(d0z0[0])>maxd0rphi) { 
694       AliDebug(1,"     rejected");
695       delete track; continue; 
696     }
697
698     // if fConstraint=kFALSE, during iterations 1 and 2 ||
699     // if fConstraint=kTRUE, during iteration 2,
700     // select tracks with d0oplusz0 < fMaxd0z0
701     if((!fConstraint && optImpParCut>0 && fVert.GetNContributors()>0) ||
702        ( fConstraint && optImpParCut==2 && cutond0z0)) {
703       if(nTrksOrig>=3 && 
704          TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1])>fMaxd0z0) { 
705         AliDebug(1,"     rejected");
706         delete track; continue; 
707       }
708     }
709     
710     // track passed all selections
711     fTrkArraySel.AddLast(track);
712     fIdSel[nTrksSel] = idOrig[i];
713     nTrksSel++; 
714   } // end loop on tracks
715
716   delete initVertex;
717
718   return nTrksSel;
719
720 //----------------------------------------------------------------------------
721 Bool_t AliVertexerTracks::PropagateTrackTo(AliExternalTrackParam *track,
722                                            Double_t xToGo) {
723   //----------------------------------------------------------------
724   // COPIED from AliTracker
725   //
726   // Propagates the track to the plane X=xk (cm).
727   //
728   //  Origin: Marian Ivanov,  Marian.Ivanov@cern.ch
729   //----------------------------------------------------------------
730
731   const Double_t kEpsilon = 0.00001;
732   Double_t xpos = track->GetX();
733   Double_t dir = (xpos<xToGo) ? 1. : -1.;
734   Double_t maxStep = 5;
735   Double_t maxSnp = 0.8;
736   //
737   while ( (xToGo-xpos)*dir > kEpsilon){
738     Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
739     Double_t x    = xpos+step;
740     Double_t xyz0[3],xyz1[3];
741     track->GetXYZ(xyz0);   //starting global position
742
743     if(!track->GetXYZAt(x,GetFieldkG(),xyz1)) return kFALSE;   // no prolongation
744     xyz1[2]+=kEpsilon; // waiting for bug correction in geo
745
746     if(TMath::Abs(track->GetSnpAt(x,GetFieldkG())) >= maxSnp) return kFALSE;
747     if(!track->PropagateTo(x,GetFieldkG()))  return kFALSE;
748
749     if(TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
750     track->GetXYZ(xyz0);   // global position
751     Double_t alphan = TMath::ATan2(xyz0[1], xyz0[0]); 
752     //
753     Double_t ca=TMath::Cos(alphan-track->GetAlpha()), 
754       sa=TMath::Sin(alphan-track->GetAlpha());
755     Double_t sf=track->GetSnp(), cf=TMath::Sqrt((1.-sf)*(1.+sf));
756     Double_t sinNew =  sf*ca - cf*sa;
757     if(TMath::Abs(sinNew) >= maxSnp) return kFALSE;
758     if(!track->Rotate(alphan)) return kFALSE;
759  
760     xpos = track->GetX();
761   }
762   return kTRUE;
763 }
764 //---------------------------------------------------------------------------
765 AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
766                                                         const TObjArray *trkArray,
767                                                         UShort_t *id,
768                                                         const Float_t *diamondxy) const
769 {
770 //
771 // Removes tracks in trksTree from fit of inVtx
772 //
773
774   if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraint")) {
775     printf("ERROR: primary vertex has no constraint: cannot remove tracks");
776     return 0x0;
777   }
778   if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
779     printf("WARNING: result of tracks' removal will be only approximately correct");
780
781   TMatrixD rv(3,1);
782   rv(0,0) = inVtx->GetXv();
783   rv(1,0) = inVtx->GetYv();
784   rv(2,0) = inVtx->GetZv();
785   TMatrixD vV(3,3);
786   Double_t cov[6];
787   inVtx->GetCovMatrix(cov);
788   vV(0,0) = cov[0];
789   vV(0,1) = cov[1]; vV(1,0) = cov[1];
790   vV(1,1) = cov[2];
791   vV(0,2) = cov[3]; vV(2,0) = cov[3];
792   vV(1,2) = cov[4]; vV(2,1) = cov[4]; 
793   vV(2,2) = cov[5];
794
795   TMatrixD sumWi(TMatrixD::kInverted,vV);
796   TMatrixD sumWiri(sumWi,TMatrixD::kMult,rv);
797
798   Int_t nUsedTrks = inVtx->GetNIndices();
799   Double_t chi2 = inVtx->GetChi2();
800
801   AliExternalTrackParam *track = 0;
802   Int_t ntrks = (Int_t)trkArray->GetEntriesFast();
803   for(Int_t i=0;i<ntrks;i++) {
804     track = (AliExternalTrackParam*)trkArray->At(i);
805     if(!inVtx->UsesTrack(id[i])) {
806       printf("track %d was not used in vertex fit",id[i]);
807       continue;
808     }
809     Double_t alpha = track->GetAlpha();
810     Double_t xl = diamondxy[0]*TMath::Cos(alpha)+diamondxy[1]*TMath::Sin(alpha);
811     track->PropagateTo(xl,GetFieldkG()); 
812     // vector of track global coordinates
813     TMatrixD ri(3,1);
814     // covariance matrix of ri
815     TMatrixD wWi(3,3);
816     
817     // get space point from track
818     if(!TrackToPoint(track,ri,wWi)) continue;
819
820     TMatrixD wWiri(wWi,TMatrixD::kMult,ri); 
821
822     sumWi -= wWi;
823     sumWiri -= wWiri;
824
825     // track contribution to chi2
826     TMatrixD deltar = rv; deltar -= ri;
827     TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
828     Double_t chi2i = deltar(0,0)*wWideltar(0,0)+
829                      deltar(1,0)*wWideltar(1,0)+
830                      deltar(2,0)*wWideltar(2,0);
831     // remove from total chi2
832     chi2 -= chi2i;
833
834     nUsedTrks--;
835     if(nUsedTrks<2) {
836       printf("Trying to remove too many tracks!");
837       return 0x0;
838     }
839   }
840
841   TMatrixD rvnew(3,1);
842   TMatrixD vVnew(3,3);
843
844   // new inverted of weights matrix
845   TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
846   vVnew = invsumWi;
847   // new position of primary vertex
848   rvnew.Mult(vVnew,sumWiri);
849
850   Double_t position[3];
851   position[0] = rvnew(0,0);
852   position[1] = rvnew(1,0);
853   position[2] = rvnew(2,0);
854   cov[0] = vVnew(0,0);
855   cov[1] = vVnew(0,1);
856   cov[2] = vVnew(1,1);
857   cov[3] = vVnew(0,2);
858   cov[4] = vVnew(1,2);
859   cov[5] = vVnew(2,2);
860   
861   // store data in the vertex object
862   AliESDVertex *outVtx = new AliESDVertex(position,cov,chi2,nUsedTrks+1); // the +1 is for the constraint
863   outVtx->SetTitle(inVtx->GetTitle());
864   UShort_t *inindices = inVtx->GetIndices();
865   Int_t nIndices = nUsedTrks;
866   UShort_t *outindices = new UShort_t[nIndices];
867   Int_t j=0;
868   for(Int_t k=0; k<inVtx->GetNIndices(); k++) {
869     Bool_t copyindex=kTRUE;
870     for(Int_t l=0; l<ntrks; l++) {
871       if(inindices[k]==id[l]) copyindex=kFALSE;
872     }
873     if(copyindex) {
874       outindices[j] = inindices[k]; j++;
875     }
876   }
877   outVtx->SetIndices(nIndices,outindices);
878   if (outindices) delete [] outindices;
879
880   /*
881     printf("Vertex before removing tracks:");
882     inVtx->PrintStatus();
883     inVtx->PrintIndices();
884     printf("Vertex after removing tracks:");
885     outVtx->PrintStatus();
886     outVtx->PrintIndices();
887   */
888
889   return outVtx;
890 }
891 //---------------------------------------------------------------------------
892 AliESDVertex* AliVertexerTracks::RemoveConstraintFromVertex(AliESDVertex *inVtx,
893                                                      Float_t *diamondxyz,
894                                                      Float_t *diamondcov) const
895 {
896 //
897 // Removes diamond constraint from fit of inVtx
898 //
899
900   if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraint")) {
901     printf("ERROR: primary vertex has no constraint: cannot remove it\n");
902     return 0x0;
903   }
904   if(inVtx->GetNContributors()<3) {
905     printf("ERROR: primary vertex has less than 2 tracks: cannot remove contraint\n");
906     return 0x0;
907   }
908
909   // diamond constraint 
910   TMatrixD vVb(3,3);
911   vVb(0,0) = diamondcov[0];
912   vVb(0,1) = diamondcov[1];
913   vVb(0,2) = 0.;
914   vVb(1,0) = diamondcov[1];
915   vVb(1,1) = diamondcov[2];
916   vVb(1,2) = 0.;
917   vVb(2,0) = 0.;
918   vVb(2,1) = 0.;
919   vVb(2,2) = diamondcov[5];
920   TMatrixD vVbInv(TMatrixD::kInverted,vVb);
921   TMatrixD rb(3,1);
922   rb(0,0) = diamondxyz[0];
923   rb(1,0) = diamondxyz[1];
924   rb(2,0) = diamondxyz[2];
925   TMatrixD vVbInvrb(vVbInv,TMatrixD::kMult,rb);
926
927   // input vertex
928   TMatrixD rv(3,1);
929   rv(0,0) = inVtx->GetXv();
930   rv(1,0) = inVtx->GetYv();
931   rv(2,0) = inVtx->GetZv();
932   TMatrixD vV(3,3);
933   Double_t cov[6];
934   inVtx->GetCovMatrix(cov);
935   vV(0,0) = cov[0];
936   vV(0,1) = cov[1]; vV(1,0) = cov[1];
937   vV(1,1) = cov[2];
938   vV(0,2) = cov[3]; vV(2,0) = cov[3];
939   vV(1,2) = cov[4]; vV(2,1) = cov[4]; 
940   vV(2,2) = cov[5];
941   TMatrixD vVInv(TMatrixD::kInverted,vV);
942   TMatrixD vVInvrv(vVInv,TMatrixD::kMult,rv);
943
944
945   TMatrixD sumWi = vVInv - vVbInv;
946
947
948   TMatrixD sumWiri = vVInvrv - vVbInvrb;
949
950   TMatrixD rvnew(3,1);
951   TMatrixD vVnew(3,3);
952
953   // new inverted of weights matrix
954   TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
955   vVnew = invsumWi;
956   // new position of primary vertex
957   rvnew.Mult(vVnew,sumWiri);
958
959   Double_t position[3];
960   position[0] = rvnew(0,0);
961   position[1] = rvnew(1,0);
962   position[2] = rvnew(2,0);
963   cov[0] = vVnew(0,0);
964   cov[1] = vVnew(0,1);
965   cov[2] = vVnew(1,1);
966   cov[3] = vVnew(0,2);
967   cov[4] = vVnew(1,2);
968   cov[5] = vVnew(2,2);
969
970
971   Double_t chi2 = inVtx->GetChi2();
972
973   // diamond constribution to chi2
974   TMatrixD deltar = rv; deltar -= rb;
975   TMatrixD vVbInvdeltar(vVbInv,TMatrixD::kMult,deltar);
976   Double_t chi2b = deltar(0,0)*vVbInvdeltar(0,0)+
977                    deltar(1,0)*vVbInvdeltar(1,0)+
978                    deltar(2,0)*vVbInvdeltar(2,0);
979   // remove from total chi2
980   chi2 -= chi2b;
981
982   // store data in the vertex object
983   AliESDVertex *outVtx = new AliESDVertex(position,cov,chi2,inVtx->GetNContributors()-1);
984   outVtx->SetTitle("VertexerTracksNoConstraint");
985   UShort_t *inindices = inVtx->GetIndices();
986   Int_t nIndices = inVtx->GetNIndices();
987   outVtx->SetIndices(nIndices,inindices);
988
989   return outVtx;
990 }
991 //---------------------------------------------------------------------------
992 void AliVertexerTracks::SetCuts(Double_t *cuts, Int_t ncuts) 
993 {
994 //
995 //  Cut values
996 //
997   if (ncuts>0) SetDCAcut(cuts[0]);
998   if (ncuts>1) SetDCAcutIter0(cuts[1]);
999   if (ncuts>2) SetMaxd0z0(cuts[2]);
1000   if (ncuts>3) if(fMode==0 && cuts[3]<0) SetITSrefitNotRequired();
1001   if (ncuts>4) SetMinClusters((Int_t)(TMath::Abs(cuts[3])));
1002   if (ncuts>5) SetMinTracks((Int_t)(cuts[4]));
1003   if (ncuts>6) SetNSigmad0(cuts[5]);
1004   if (ncuts>7) SetMinDetFitter(cuts[6]);
1005   if (ncuts>8) SetMaxTgl(cuts[7]);
1006   if (ncuts>9) SetFiducialRZ(cuts[8],cuts[9]);
1007   if (ncuts>10) fAlgo=(Int_t)(cuts[10]);
1008   if (ncuts>11) fAlgoIter0=(Int_t)(cuts[11]);
1009   //
1010   if (ncuts>12) if (cuts[12]>1.)   SetMVTukey2(cuts[12]);
1011   if (ncuts>13) if (cuts[13]>1.)   SetMVSig2Ini(cuts[13]);
1012   if (ncuts>14) if (cuts[14]>0.1)  SetMVMaxSigma2(cuts[14]);
1013   if (ncuts>15) if (cuts[15]>1e-5) SetMVMinSig2Red(cuts[15]);
1014   if (ncuts>16) if (cuts[16]>1e-5) SetMVMinDst(cuts[16]);
1015   if (ncuts>17) if (cuts[17]>0.5)  SetMVScanStep(cuts[17]);
1016   if (ncuts>18) SetMVMaxWghNtr(cuts[18]);
1017   if (ncuts>19) SetMVFinalWBinary(cuts[19]>0);
1018   if (ncuts>20) if (cuts[20]>20.)  SetBCSpacing(int(cuts[20]));
1019   //
1020   if (fAlgo==kMultiVertexer) SetSelectOnTOFBunchCrossing(kTRUE,kTRUE);
1021   else                       SetSelectOnTOFBunchCrossing(kFALSE,kTRUE);
1022   //
1023   return;
1024 }
1025 //---------------------------------------------------------------------------
1026 void AliVertexerTracks::SetITSMode(Double_t dcacut,
1027                                    Double_t dcacutIter0,
1028                                    Double_t maxd0z0,
1029                                    Int_t minCls,
1030                                    Int_t mintrks,
1031                                    Double_t nsigma,
1032                                    Double_t mindetfitter,
1033                                    Double_t maxtgl,
1034                                    Double_t fidR,
1035                                    Double_t fidZ,
1036                                    Int_t finderAlgo,
1037                                    Int_t finderAlgoIter0)
1038 {
1039 //
1040 //  Cut values for ITS mode
1041 //
1042   fMode = 0;
1043   if(minCls>0) {
1044     SetITSrefitRequired();
1045   } else {
1046     SetITSrefitNotRequired();
1047   }
1048   SetDCAcut(dcacut);
1049   SetDCAcutIter0(dcacutIter0);
1050   SetMaxd0z0(maxd0z0);
1051   SetMinClusters(TMath::Abs(minCls));
1052   SetMinTracks(mintrks);
1053   SetNSigmad0(nsigma);
1054   SetMinDetFitter(mindetfitter);
1055   SetMaxTgl(maxtgl);
1056   SetFiducialRZ(fidR,fidZ);
1057   fAlgo=finderAlgo;
1058   fAlgoIter0=finderAlgoIter0;
1059
1060   return; 
1061 }
1062 //---------------------------------------------------------------------------
1063 void AliVertexerTracks::SetTPCMode(Double_t dcacut,
1064                                    Double_t dcacutIter0,
1065                                    Double_t maxd0z0,
1066                                    Int_t minCls,
1067                                    Int_t mintrks,
1068                                    Double_t nsigma,
1069                                    Double_t mindetfitter,
1070                                    Double_t maxtgl,
1071                                    Double_t fidR,
1072                                    Double_t fidZ,
1073                                    Int_t finderAlgo,
1074                                    Int_t finderAlgoIter0) 
1075 {
1076 //
1077 //  Cut values for TPC mode
1078 //
1079   fMode = 1;
1080   SetITSrefitNotRequired();
1081   SetDCAcut(dcacut);
1082   SetDCAcutIter0(dcacutIter0);
1083   SetMaxd0z0(maxd0z0);
1084   SetMinClusters(minCls);
1085   SetMinTracks(mintrks);
1086   SetNSigmad0(nsigma);
1087   SetMinDetFitter(mindetfitter);
1088   SetMaxTgl(maxtgl);
1089   SetFiducialRZ(fidR,fidZ);
1090   fAlgo=finderAlgo;
1091   fAlgoIter0=finderAlgoIter0;
1092
1093   return; 
1094 }
1095 //---------------------------------------------------------------------------
1096 void AliVertexerTracks::SetSkipTracks(Int_t n,const Int_t *skipped) 
1097 {
1098 //
1099 // Mark the tracks not to be used in the vertex reconstruction.
1100 // Tracks are identified by AliESDtrack::GetID()
1101 //
1102   fNTrksToSkip = n;  fTrksToSkip = new Int_t[n]; 
1103   for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i]; 
1104   return; 
1105 }
1106 //---------------------------------------------------------------------------
1107 void  AliVertexerTracks::SetVtxStart(AliESDVertex *vtx) 
1108
1109 //
1110 // Set initial vertex knowledge
1111 //
1112   vtx->GetXYZ(fNominalPos);
1113   vtx->GetCovMatrix(fNominalCov);
1114   SetConstraintOn();
1115   return; 
1116 }
1117 //---------------------------------------------------------------------------
1118 void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
1119 {
1120   AliExternalTrackParam *track1;
1121   const Int_t knacc = (Int_t)fTrkArraySel.GetEntriesFast();
1122   AliStrLine **linarray = new AliStrLine* [knacc];
1123   for(Int_t i=0; i<knacc; i++){
1124     track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
1125     Double_t alpha=track1->GetAlpha();
1126     Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
1127     Double_t pos[3],dir[3],sigmasq[3]; 
1128     track1->GetXYZAt(mindist,GetFieldkG(),pos);
1129     track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
1130     sigmasq[0]=TMath::Sin(alpha)*TMath::Sin(alpha)*track1->GetSigmaY2();
1131     sigmasq[1]=TMath::Cos(alpha)*TMath::Cos(alpha)*track1->GetSigmaY2();
1132     sigmasq[2]=track1->GetSigmaZ2();
1133     TMatrixD ri(3,1);
1134     TMatrixD wWi(3,3);
1135     if(!TrackToPoint(track1,ri,wWi)) {optUseWeights=kFALSE;printf("WARNING\n");}
1136     Double_t wmat[9];
1137     Int_t iel=0;
1138     for(Int_t ia=0;ia<3;ia++){
1139       for(Int_t ib=0;ib<3;ib++){
1140         wmat[iel]=wWi(ia,ib);
1141         iel++;
1142       }    
1143     }
1144     linarray[i] = new AliStrLine(pos,sigmasq,wmat,dir);
1145   }
1146   fVert=TrackletVertexFinder(linarray,knacc,optUseWeights);
1147   for(Int_t i=0; i<knacc; i++) delete linarray[i];
1148   delete [] linarray;
1149 }
1150 //---------------------------------------------------------------------------
1151 AliESDVertex AliVertexerTracks::TrackletVertexFinder(const TClonesArray *lines, Int_t optUseWeights)
1152 {
1153   // Calculate the point at minimum distance to prepared tracks (TClonesArray)
1154   const Int_t knacc = (Int_t)lines->GetEntriesFast();
1155   AliStrLine** lines2 = new AliStrLine* [knacc];
1156   for(Int_t i=0; i<knacc; i++){
1157     lines2[i]= (AliStrLine*)lines->At(i);
1158   }
1159   AliESDVertex vert = TrackletVertexFinder(lines2,knacc,optUseWeights); 
1160   delete [] lines2;
1161   return vert;
1162 }
1163
1164 //---------------------------------------------------------------------------
1165 AliESDVertex AliVertexerTracks::TrackletVertexFinder(AliStrLine **lines, const Int_t knacc, Int_t optUseWeights)
1166 {
1167   // Calculate the point at minimum distance to prepared tracks (array of AliStrLine) 
1168
1169   Double_t initPos[3]={0.,0.,0.};
1170
1171   Double_t (*vectP0)[3]=new Double_t [knacc][3];
1172   Double_t (*vectP1)[3]=new Double_t [knacc][3];
1173   
1174   Double_t sum[3][3];
1175   Double_t dsum[3]={0,0,0};
1176   TMatrixD sumWi(3,3);
1177   for(Int_t i=0;i<3;i++){
1178     for(Int_t j=0;j<3;j++){
1179       sum[i][j]=0;
1180       sumWi(i,j)=0.;
1181     }
1182   }
1183
1184   for(Int_t i=0; i<knacc; i++){
1185     AliStrLine *line1 = lines[i]; 
1186     Double_t p0[3],cd[3],sigmasq[3];
1187     Double_t wmat[9];
1188     if(!line1) {printf("ERROR %d %d\n",i,knacc); continue;}
1189     line1->GetP0(p0);
1190     line1->GetCd(cd);
1191     line1->GetSigma2P0(sigmasq);
1192     line1->GetWMatrix(wmat);
1193     TMatrixD wWi(3,3);
1194     Int_t iel=0;
1195     for(Int_t ia=0;ia<3;ia++){
1196       for(Int_t ib=0;ib<3;ib++){
1197         wWi(ia,ib)=wmat[iel];
1198         iel++;
1199       }    
1200     }
1201
1202     sumWi+=wWi;
1203
1204     Double_t p1[3]={p0[0]+cd[0],p0[1]+cd[1],p0[2]+cd[2]};
1205     vectP0[i][0]=p0[0];
1206     vectP0[i][1]=p0[1];
1207     vectP0[i][2]=p0[2];
1208     vectP1[i][0]=p1[0];
1209     vectP1[i][1]=p1[1];
1210     vectP1[i][2]=p1[2];
1211     
1212     Double_t matr[3][3];
1213     Double_t dknow[3];
1214     if(optUseWeights==0)GetStrLinDerivMatrix(p0,p1,matr,dknow);
1215     else GetStrLinDerivMatrix(p0,p1,sigmasq,matr,dknow);
1216
1217
1218     for(Int_t iii=0;iii<3;iii++){
1219       dsum[iii]+=dknow[iii]; 
1220       for(Int_t lj=0;lj<3;lj++) sum[iii][lj]+=matr[iii][lj];
1221     }
1222   }
1223  
1224   TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
1225   Double_t covmatrix[6];
1226   covmatrix[0] = invsumWi(0,0);
1227   covmatrix[1] = invsumWi(0,1);
1228   covmatrix[2] = invsumWi(1,1);
1229   covmatrix[3] = invsumWi(0,2);
1230   covmatrix[4] = invsumWi(1,2);
1231   covmatrix[5] = invsumWi(2,2);
1232
1233   Double_t vett[3][3];
1234   Double_t det=GetDeterminant3X3(sum);
1235   Double_t sigma=0;
1236   
1237   if(TMath::Abs(det) > kAlmost0){
1238     for(Int_t zz=0;zz<3;zz++){
1239       for(Int_t ww=0;ww<3;ww++){
1240         for(Int_t kk=0;kk<3;kk++) vett[ww][kk]=sum[ww][kk];
1241       }
1242       for(Int_t kk=0;kk<3;kk++) vett[kk][zz]=dsum[kk];
1243       initPos[zz]=GetDeterminant3X3(vett)/det;
1244     }
1245
1246
1247     for(Int_t i=0; i<knacc; i++){
1248       Double_t p0[3]={0,0,0},p1[3]={0,0,0};
1249       for(Int_t ii=0;ii<3;ii++){
1250         p0[ii]=vectP0[i][ii];
1251         p1[ii]=vectP1[i][ii];
1252       }
1253       sigma+=GetStrLinMinDist(p0,p1,initPos);
1254     }
1255
1256     if(sigma>0.) {sigma=TMath::Sqrt(sigma);}else{sigma=999;}
1257   }else{
1258     sigma=999;
1259   }
1260   AliESDVertex theVert(initPos,covmatrix,99999.,knacc);
1261   theVert.SetDispersion(sigma);
1262   delete [] vectP0;
1263   delete [] vectP1;
1264   return theVert;
1265 }
1266
1267 //---------------------------------------------------------------------------
1268 Bool_t AliVertexerTracks::TrackToPoint(AliExternalTrackParam *t,
1269                                        TMatrixD &ri,TMatrixD &wWi,
1270                                        Bool_t uUi3by3) const 
1271 {
1272 //
1273 // Extract from the AliExternalTrackParam the global coordinates ri and covariance matrix
1274 // wWi of the space point that it represents (to be used in VertexFitter())
1275 //
1276
1277   
1278   Double_t rotAngle = t->GetAlpha();
1279   if(rotAngle<0.) rotAngle += 2.*TMath::Pi();
1280   Double_t cosRot = TMath::Cos(rotAngle);
1281   Double_t sinRot = TMath::Sin(rotAngle);
1282   /*
1283   // RS >>>
1284   Double_t lambda = TMath::ATan(t->GetTgl());
1285   Double_t cosLam   = TMath::Cos(lambda);
1286   Double_t sinLam   = TMath::Sin(lambda);
1287   // RS <<<
1288   */
1289   ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot;
1290   ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot;
1291   ri(2,0) = t->GetZ();
1292
1293   if(!uUi3by3) {
1294     // matrix to go from global (x,y,z) to local (y,z);
1295     TMatrixD qQi(2,3);
1296     qQi(0,0) = -sinRot;
1297     qQi(0,1) = cosRot;
1298     qQi(0,2) = 0.;
1299     //
1300     qQi(1,0) = 0.;
1301     qQi(1,1) = 0.;
1302     qQi(1,2) = 1.;
1303     //
1304     // RS: Added polar inclination
1305     /*
1306     qQi(1,0) = -sinLam*cosRot;
1307     qQi(1,1) = -sinLam*sinRot;
1308     qQi(1,2) = cosLam;
1309     */
1310     // covariance matrix of local (y,z) - inverted
1311     TMatrixD uUi(2,2);
1312     uUi(0,0) = t->GetSigmaY2();
1313     uUi(0,1) = t->GetSigmaZY();
1314     uUi(1,0) = t->GetSigmaZY();
1315     uUi(1,1) = t->GetSigmaZ2();
1316     //printf(" Ui :");
1317     //printf(" %f   %f",uUi(0,0),uUi(0,1));
1318     //printf(" %f   %f",uUi(1,0),uUi(1,1));
1319
1320     if(uUi.Determinant() <= 0.) return kFALSE;
1321     TMatrixD uUiInv(TMatrixD::kInverted,uUi);
1322   
1323     // weights matrix: wWi = qQiT * uUiInv * qQi
1324     TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
1325     TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
1326     wWi = m;
1327   } else {
1328     if(fVert.GetNContributors()<1) AliFatal("Vertex from finder is empty");
1329     // matrix to go from global (x,y,z) to local (x,y,z);
1330     TMatrixD qQi(3,3);
1331     qQi(0,0) = cosRot;
1332     qQi(0,1) = sinRot;
1333     qQi(0,2) = 0.;
1334     qQi(1,0) = -sinRot;
1335     qQi(1,1) = cosRot;
1336     qQi(1,2) = 0.;
1337     qQi(2,0) = 0.;
1338     qQi(2,1) = 0.;
1339     qQi(2,2) = 1.;
1340    
1341     // covariance of fVert along the track  
1342     Double_t p[3],pt,ptot;
1343     t->GetPxPyPz(p);
1344     pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
1345     ptot = TMath::Sqrt(pt*pt+p[2]*p[2]);
1346     Double_t cphi = p[0]/pt;               //cos(phi)=px/pt
1347     Double_t sphi = p[1]/pt;               //sin(phi)=py/pt
1348     Double_t clambda = pt/ptot;            //cos(lambda)=pt/ptot
1349     Double_t slambda = p[2]/ptot;            //sin(lambda)=pz/ptot
1350     Double_t covfVert[6];
1351     fVert.GetCovMatrix(covfVert);
1352     Double_t covfVertalongt = 
1353        covfVert[0]*cphi*cphi*clambda*clambda 
1354       +covfVert[1]*2.*cphi*sphi*clambda*clambda
1355       +covfVert[3]*2.*cphi*clambda*slambda 
1356       +covfVert[2]*sphi*sphi*clambda*clambda 
1357       +covfVert[4]*2.*sphi*clambda*slambda 
1358       +covfVert[5]*slambda*slambda; 
1359     // covariance matrix of local (x,y,z) - inverted
1360     TMatrixD uUi(3,3);
1361     uUi(0,0) = covfVertalongt * fnSigmaForUi00 * fnSigmaForUi00;
1362     AliDebug(1,Form("=====> sqrtUi00 cm  %f",TMath::Sqrt(uUi(0,0))));
1363     uUi(0,1) = 0.;
1364     uUi(0,2) = 0.;
1365     uUi(1,0) = 0.;
1366     uUi(1,1) = t->GetSigmaY2();
1367     uUi(1,2) = t->GetSigmaZY();
1368     uUi(2,0) = 0.;
1369     uUi(2,1) = t->GetSigmaZY();
1370     uUi(2,2) = t->GetSigmaZ2();
1371     //printf(" Ui :\n");
1372     //printf(" %f   %f\n",uUi(0,0),uUi(0,1));
1373     //printf(" %f   %f\n",uUi(1,0),uUi(1,1));
1374   
1375     if(uUi.Determinant() <= 0.) return kFALSE;
1376     TMatrixD uUiInv(TMatrixD::kInverted,uUi);
1377   
1378     // weights matrix: wWi = qQiT * uUiInv * qQi
1379     TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
1380     TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
1381     wWi = m;
1382   }
1383
1384
1385   return kTRUE;
1386
1387 //---------------------------------------------------------------------------
1388 void AliVertexerTracks::TooFewTracks() 
1389 {
1390 //
1391 // When the number of tracks is < fMinTracks,
1392 // deal with vertices not found and prepare to exit
1393 //
1394   AliDebug(1,"TooFewTracks");
1395
1396   Double_t pos[3],err[3];
1397   pos[0] = fNominalPos[0];
1398   err[0] = TMath::Sqrt(fNominalCov[0]);
1399   pos[1] = fNominalPos[1];
1400   err[1] = TMath::Sqrt(fNominalCov[2]);
1401   pos[2] = fNominalPos[2];
1402   err[2] = TMath::Sqrt(fNominalCov[5]);
1403   Int_t    ncontr = (err[0]>1. ? -1 : -3);
1404   if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
1405   fCurrentVertex = new AliESDVertex(pos,err);
1406   fCurrentVertex->SetNContributors(ncontr);
1407
1408   if(fConstraint) {
1409     fCurrentVertex->SetTitle("VertexerTracksWithConstraint");
1410   } else {
1411     fCurrentVertex->SetTitle("VertexerTracksNoConstraint");
1412   }
1413
1414   if(!fTrkArraySel.IsEmpty()) fTrkArraySel.Delete(); 
1415   if(fIdSel) {delete [] fIdSel; fIdSel=NULL;}
1416   if(fTrksToSkip) {delete [] fTrksToSkip; fTrksToSkip=NULL;}
1417
1418   return;
1419 }
1420 //---------------------------------------------------------------------------
1421 void AliVertexerTracks::VertexFinder(Int_t optUseWeights) 
1422 {
1423
1424   // Get estimate of vertex position in (x,y) from tracks DCA
1425  
1426   Double_t initPos[3];
1427   initPos[2] = 0.;
1428   for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
1429   Int_t nacc = (Int_t)fTrkArraySel.GetEntriesFast();
1430   Double_t aver[3]={0.,0.,0.};
1431   Double_t aversq[3]={0.,0.,0.};
1432   Double_t sigmasq[3]={0.,0.,0.};
1433   Double_t sigma=0;
1434   Int_t ncombi = 0;
1435   AliExternalTrackParam *track1;
1436   AliExternalTrackParam *track2;
1437   Double_t pos[3],dir[3]; 
1438   Double_t alpha,mindist;
1439
1440   for(Int_t i=0; i<nacc; i++){
1441     track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
1442     alpha=track1->GetAlpha();
1443     mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
1444     track1->GetXYZAt(mindist,GetFieldkG(),pos);
1445     track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
1446     AliStrLine *line1 = new AliStrLine(pos,dir); 
1447
1448    //    AliStrLine *line1 = new AliStrLine();
1449    //    track1->ApproximateHelixWithLine(mindist,GetFieldkG(),line1);
1450    
1451     for(Int_t j=i+1; j<nacc; j++){
1452       track2 = (AliExternalTrackParam*)fTrkArraySel.At(j);
1453       alpha=track2->GetAlpha();
1454       mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
1455       track2->GetXYZAt(mindist,GetFieldkG(),pos);
1456       track2->GetPxPyPzAt(mindist,GetFieldkG(),dir);
1457       AliStrLine *line2 = new AliStrLine(pos,dir); 
1458     //      AliStrLine *line2 = new AliStrLine();
1459     //  track2->ApproximateHelixWithLine(mindist,GetFieldkG(),line2);
1460       Double_t distCA=line2->GetDCA(line1);
1461       //printf("%d   %d   %f\n",i,j,distCA);
1462        if(fDCAcut<=0 || (fDCAcut>0&&distCA<fDCAcut)){
1463         Double_t pnt1[3],pnt2[3],crosspoint[3];
1464
1465         if(optUseWeights<=0){
1466           Int_t retcode = line2->Cross(line1,crosspoint);
1467           if(retcode>=0){
1468             ncombi++;
1469             for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
1470             for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
1471           }
1472         }
1473         if(optUseWeights>0){
1474           Int_t retcode = line1->CrossPoints(line2,pnt1,pnt2);
1475           if(retcode>=0){
1476             Double_t cs, sn;
1477             alpha=track1->GetAlpha();
1478             cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);   
1479             Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
1480             alpha=track2->GetAlpha();
1481             cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
1482             Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
1483             Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
1484             Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
1485             Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
1486             Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
1487             crosspoint[0]=wx1*pnt1[0] + wx2*pnt2[0]; 
1488             crosspoint[1]=wy1*pnt1[1] + wy2*pnt2[1]; 
1489             crosspoint[2]=wz1*pnt1[2] + wz2*pnt2[2];
1490           
1491             ncombi++;
1492             for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
1493             for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
1494           }
1495         }
1496       }
1497       delete line2;
1498     }
1499     delete line1;
1500   }
1501   if(ncombi>0){
1502     for(Int_t jj=0;jj<3;jj++){
1503       initPos[jj] = aver[jj]/ncombi;
1504       //printf("%f\n",initPos[jj]);
1505       aversq[jj]/=ncombi;
1506       sigmasq[jj]=aversq[jj]-initPos[jj]*initPos[jj];
1507       sigma+=sigmasq[jj];
1508     }
1509     sigma=TMath::Sqrt(TMath::Abs(sigma));
1510   }
1511   else {
1512     Warning("VertexFinder","Finder did not succed");
1513     sigma=999;
1514   }
1515   fVert.SetXYZ(initPos);
1516   fVert.SetDispersion(sigma);
1517   fVert.SetNContributors(ncombi);
1518 }
1519 //---------------------------------------------------------------------------
1520 void AliVertexerTracks::VertexFitter(Bool_t vfit, Bool_t chiCalc,Int_t useWeights) 
1521 {
1522 //
1523 // The optimal estimate of the vertex position is given by a "weighted 
1524 // average of tracks positions".
1525 // Original method: V. Karimaki, CMS Note 97/0051
1526 //
1527   const double kTiny = 1e-9;
1528   Bool_t useConstraint = fConstraint;
1529   Double_t initPos[3];
1530   if(!fOnlyFitter) {
1531     fVert.GetXYZ(initPos);
1532   } else {
1533     initPos[0]=fNominalPos[0];
1534     initPos[1]=fNominalPos[1];
1535     initPos[2]=fNominalPos[2];
1536   }
1537
1538   Int_t nTrksSel = (Int_t)fTrkArraySel.GetEntries();
1539   if(nTrksSel==1) useConstraint=kTRUE;
1540   AliDebug(1,Form("--- VertexFitter(): start (%d,%d,%d)",vfit,chiCalc,useWeights));
1541   AliDebug(1,Form(" Number of tracks in array: %d\n",nTrksSel));
1542   AliDebug(1,Form(" Minimum # tracks required in fit: %d\n",fMinTracks));
1543   AliDebug(1,Form(" Vertex position after finder: %f,%f,%f\n",initPos[0],initPos[1],initPos[2]));
1544   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]))); 
1545
1546   // special treatment for few-tracks fits (e.g. secondary vertices)
1547   Bool_t uUi3by3 = kFALSE; if(nTrksSel<5 && !useConstraint && !useWeights) uUi3by3 = kTRUE;
1548
1549   Int_t i,j,k,step=0;
1550   static TMatrixD rv(3,1);
1551   static TMatrixD vV(3,3);
1552   rv(0,0) = initPos[0];
1553   rv(1,0) = initPos[1];
1554   rv(2,0) = initPos[2];
1555   Double_t xlStart,alpha;
1556   Int_t nTrksUsed = 0;
1557   Double_t chi2=0,chi2i,chi2b;
1558   AliExternalTrackParam *t = 0;
1559   Int_t failed = 0;
1560
1561   // initial vertex covariance matrix
1562   TMatrixD vVb(3,3);
1563   vVb(0,0) = fNominalCov[0];
1564   vVb(0,1) = fNominalCov[1];
1565   vVb(0,2) = 0.;
1566   vVb(1,0) = fNominalCov[1];
1567   vVb(1,1) = fNominalCov[2];
1568   vVb(1,2) = 0.;
1569   vVb(2,0) = 0.;
1570   vVb(2,1) = 0.;
1571   vVb(2,2) = fNominalCov[5];
1572   TMatrixD vVbInv(TMatrixD::kInverted,vVb);
1573   TMatrixD rb(3,1);
1574   rb(0,0) = fNominalPos[0];
1575   rb(1,0) = fNominalPos[1];
1576   rb(2,0) = fNominalPos[2];
1577   TMatrixD vVbInvrb(vVbInv,TMatrixD::kMult,rb);
1578   //
1579   int currBC = fVert.GetBC();
1580   //
1581   // 2 steps:
1582   // 1st - estimate of vtx using all tracks
1583   // 2nd - estimate of global chi2
1584   //
1585   for(step=0; step<2; step++) {
1586     if      (step==0 && !vfit)    continue;
1587     else if (step==1 && !chiCalc) continue;
1588     chi2 = 0.;
1589     nTrksUsed = 0;
1590     fMVWSum = fMVWE2 = 0;
1591     if(step==1) { initPos[0]=rv(0,0); initPos[1]=rv(1,0); initPos[2]=rv(2,0);}
1592     AliDebug(2,Form("Step%d: inipos: %+f %+f %+f MinTr: %d, Sig2:%.2f)",step,initPos[0],initPos[1],initPos[2],fMinTracks,fMVSigma2));
1593     //
1594     TMatrixD sumWiri(3,1);
1595     TMatrixD sumWi(3,3);
1596     for(i=0; i<3; i++) {
1597       sumWiri(i,0) = 0.;
1598       for(j=0; j<3; j++) sumWi(j,i) = 0.;
1599     }
1600
1601     // mean vertex constraint
1602     if(useConstraint) {
1603       for(i=0;i<3;i++) {
1604         sumWiri(i,0) += vVbInvrb(i,0);
1605         for(k=0;k<3;k++) sumWi(i,k) += vVbInv(i,k);
1606       }
1607       // chi2
1608       TMatrixD deltar = rv; deltar -= rb;
1609       TMatrixD vVbInvdeltar(vVbInv,TMatrixD::kMult,deltar);
1610       chi2b = deltar(0,0)*vVbInvdeltar(0,0)+
1611               deltar(1,0)*vVbInvdeltar(1,0)+
1612               deltar(2,0)*vVbInvdeltar(2,0);
1613       chi2 += chi2b;
1614     }
1615
1616     // loop on tracks  
1617     for(k=0; k<nTrksSel; k++) {
1618       //
1619       // get track from track array
1620       t = (AliExternalTrackParam*)fTrkArraySel.At(k);
1621       if (useWeights && t->TestBit(kBitUsed)) continue;
1622       //      
1623       int tBC = int(t->GetUniqueID()) - kTOFBCShift;    // BC assigned to this track
1624       if (fSelectOnTOFBunchCrossing) {
1625         if (!fKeepAlsoUnflaggedTOFBunchCrossing) continue;   // don't consider tracks with undefined BC 
1626         if (currBC!=AliVTrack::kTOFBCNA && tBC!=AliVTrack::kTOFBCNA && tBC!=currBC) continue;  // track does not match to current BCid
1627       }
1628       alpha = t->GetAlpha();
1629       xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha);
1630       // to vtxSeed (from finder)
1631       t->PropagateTo(xlStart,GetFieldkG());   
1632  
1633       // vector of track global coordinates
1634       TMatrixD ri(3,1);
1635       // covariance matrix of ri
1636       TMatrixD wWi(3,3);
1637
1638       // get space point from track
1639       if(!TrackToPoint(t,ri,wWi,uUi3by3)) continue;
1640
1641       // track chi2
1642       TMatrixD deltar = rv; deltar -= ri;
1643       TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
1644       chi2i = deltar(0,0)*wWideltar(0,0)+
1645               deltar(1,0)*wWideltar(1,0)+
1646               deltar(2,0)*wWideltar(2,0);
1647       //
1648       if (useWeights) {
1649         //double sg = TMath::Sqrt(fMVSigma2);
1650         //double chi2iw = (deltar(0,0)*wWideltar(0,0)+deltar(1,0)*wWideltar(1,0))/sg + deltar(2,0)*wWideltar(2,0)/fMVSigma2;
1651         //double wgh = (1-chi2iw/fMVTukey2); 
1652         double chi2iw = chi2i;
1653         double wgh = (1-chi2iw/fMVTukey2/fMVSigma2); 
1654
1655         if (wgh<kTiny) wgh = 0;
1656         else if (useWeights==2) wgh = 1.; // use as binary weight
1657         if (step==1) ((AliESDtrack*)t)->SetBit(kBitUsed, wgh>0);
1658         if (wgh<kTiny) continue; // discard the track
1659         wWi *= wgh;  // RS: use weight?
1660         fMVWSum += wgh;
1661         fMVWE2  += wgh*chi2iw;
1662       }
1663       // add to total chi2
1664       if (fSelectOnTOFBunchCrossing && tBC!=AliVTrack::kTOFBCNA && currBC<0) currBC = tBC;
1665       //
1666       chi2 += chi2i;
1667       TMatrixD wWiri(wWi,TMatrixD::kMult,ri); 
1668       sumWiri += wWiri;
1669       sumWi   += wWi;
1670
1671       nTrksUsed++;
1672     } // end loop on tracks
1673
1674     if(nTrksUsed < fMinTracks) {
1675       failed=1;
1676       continue;
1677     }
1678
1679     Double_t determinant = sumWi.Determinant();
1680     if(determinant < fMinDetFitter)  { 
1681       AliDebug(1,Form("det(V) = %f (<%f)\n",determinant,fMinDetFitter));       
1682       failed=1;
1683       continue;
1684     }
1685
1686     if(step==0) { 
1687       // inverted of weights matrix
1688       TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
1689       vV = invsumWi;
1690       // position of primary vertex
1691       rv.Mult(vV,sumWiri);
1692     }
1693   } // end loop on the 2 steps
1694
1695   if(failed) { 
1696     fVert.SetNContributors(-1);
1697     if (chiCalc) {
1698       TooFewTracks();
1699       if (fCurrentVertex) fVert = *fCurrentVertex;  // RS
1700     }
1701     return; 
1702   }
1703   //
1704   Double_t position[3];
1705   position[0] = rv(0,0);
1706   position[1] = rv(1,0);
1707   position[2] = rv(2,0);
1708   Double_t covmatrix[6];
1709   covmatrix[0] = vV(0,0);
1710   covmatrix[1] = vV(0,1);
1711   covmatrix[2] = vV(1,1);
1712   covmatrix[3] = vV(0,2);
1713   covmatrix[4] = vV(1,2);
1714   covmatrix[5] = vV(2,2);
1715   
1716   // for correct chi2/ndf, count constraint as additional "track"
1717   if(fConstraint) nTrksUsed++;
1718   //
1719   if (vfit && !chiCalc) { // RS: special mode for multi-vertex finder
1720     fVert.SetXYZ(position);
1721     fVert.SetCovarianceMatrix(covmatrix);
1722     fVert.SetNContributors(nTrksUsed);
1723     return;
1724   } 
1725   // store data in the vertex object
1726   if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
1727   fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);
1728   fCurrentVertex->SetBC(currBC);
1729   fVert = *fCurrentVertex;  // RS
1730   AliDebug(1," Vertex after fit:");
1731   AliDebug(1,Form("xyz: %f %f %f; nc %d",fCurrentVertex->GetXv(),fCurrentVertex->GetYv(),fCurrentVertex->GetZv(),fCurrentVertex->GetNContributors()));
1732   AliDebug(1,"--- VertexFitter(): finish\n");
1733  
1734
1735   return;
1736 }
1737 //----------------------------------------------------------------------------
1738 AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(const TObjArray *trkArray,
1739                                                          UShort_t *id,
1740                                                          Bool_t optUseFitter,
1741                                                          Bool_t optPropagate,
1742                                                          Bool_t optUseDiamondConstraint) 
1743 {
1744 //
1745 // Return vertex from tracks (AliExternalTrackParam) in array
1746 //
1747   fCurrentVertex = 0;
1748   // set optUseDiamondConstraint=TRUE only if you are reconstructing the 
1749   // primary vertex!
1750   if(optUseDiamondConstraint) {
1751     SetConstraintOn();
1752   } else {    
1753     SetConstraintOff();
1754   }
1755
1756   // get tracks and propagate them to initial vertex position
1757   fIdSel = new UShort_t[(Int_t)trkArray->GetEntriesFast()];
1758   Int_t nTrksSel = PrepareTracks(*trkArray,id,0);
1759   if((!optUseDiamondConstraint && nTrksSel<TMath::Max(2,fMinTracks)) ||
1760      (optUseDiamondConstraint && nTrksSel<1)) {
1761     TooFewTracks();
1762     return fCurrentVertex;
1763   }
1764  
1765   // vertex finder
1766   if(nTrksSel==1) {
1767     AliDebug(1,"Just one track");
1768     OneTrackVertFinder();
1769   } else {
1770     switch (fAlgo) {
1771       case 1: StrLinVertexFinderMinDist(1); break;
1772       case 2: StrLinVertexFinderMinDist(0); break;
1773       case 3: HelixVertexFinder();          break;
1774       case 4: VertexFinder(1);              break;
1775       case 5: VertexFinder(0);              break;
1776       default: printf("Wrong algorithm\n"); break;  
1777     }
1778   }
1779   AliDebug(1," Vertex finding completed\n");
1780
1781   // vertex fitter
1782   if(optUseFitter) {
1783     VertexFitter();
1784   } else {
1785     Double_t position[3]={fVert.GetXv(),fVert.GetYv(),fVert.GetZv()};
1786     Double_t covmatrix[6];
1787     fVert.GetCovMatrix(covmatrix);
1788     Double_t chi2=99999.;
1789     Int_t    nTrksUsed=fVert.GetNContributors();
1790     fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);    
1791   }
1792   fCurrentVertex->SetDispersion(fVert.GetDispersion());
1793
1794
1795   // set indices of used tracks and propagate track to found vertex
1796   UShort_t *indices = 0;
1797   Double_t d0z0[2],covd0z0[3];
1798   AliExternalTrackParam *t = 0;
1799   if(fCurrentVertex->GetNContributors()>0) {
1800     indices = new UShort_t[fTrkArraySel.GetEntriesFast()];
1801     for(Int_t jj=0; jj<(Int_t)fTrkArraySel.GetEntriesFast(); jj++) {
1802       indices[jj] = fIdSel[jj];
1803       t = (AliExternalTrackParam*)fTrkArraySel.At(jj);
1804       if(optPropagate && optUseFitter) {
1805         if(TMath::Sqrt(fCurrentVertex->GetXv()*fCurrentVertex->GetXv()+fCurrentVertex->GetYv()*fCurrentVertex->GetYv())<3.) {
1806           t->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
1807           AliDebug(1,Form("Track %d propagated to found vertex",jj));
1808         } else {
1809           AliWarning("Found vertex outside beam pipe!");
1810         }
1811       }
1812     }
1813     fCurrentVertex->SetIndices(fCurrentVertex->GetNContributors(),indices);
1814   }
1815
1816   // clean up
1817   if (indices) {delete [] indices; indices=NULL;}
1818   delete [] fIdSel; fIdSel=NULL;
1819   fTrkArraySel.Delete();
1820   
1821   return fCurrentVertex;
1822 }
1823  
1824 //----------------------------------------------------------------------------
1825 AliESDVertex* AliVertexerTracks::VertexForSelectedESDTracks(TObjArray *trkArray,
1826                                                              Bool_t optUseFitter,
1827                                                             Bool_t optPropagate,
1828                                                             Bool_t optUseDiamondConstraint)
1829
1830 {
1831 //
1832 // Return vertex from array of ESD tracks
1833 //
1834
1835   Int_t nTrks = (Int_t)trkArray->GetEntriesFast();
1836   UShort_t *id = new UShort_t[nTrks];
1837
1838   AliESDtrack *esdt = 0;
1839   for(Int_t i=0; i<nTrks; i++){
1840     esdt = (AliESDtrack*)trkArray->At(i);
1841     id[i] = (UShort_t)esdt->GetID();
1842   }
1843     
1844   VertexForSelectedTracks(trkArray,id,optUseFitter,optPropagate,optUseDiamondConstraint);
1845
1846   delete [] id; id=NULL;
1847
1848   return fCurrentVertex;
1849 }
1850  
1851 //______________________________________________________
1852 Bool_t AliVertexerTracks::FindNextVertexMV()
1853 {
1854   // try to find a new vertex
1855   fMVSigma2 = fMVSig2Ini;
1856   double prevSig2 = fMVSigma2;
1857   double minDst2 = fMVMinDst*fMVMinDst;
1858   const double kSigLimit = 1.0;
1859   const double kSigLimitE = kSigLimit+1e-6;
1860   const double kPushFactor = 0.5;
1861   const int kMaxIter = 20;
1862   double push = kPushFactor;
1863   //
1864   int iter = 0;
1865   double posP[3]={0,0,0},pos[3]={0,0,0};
1866   fVert.GetXYZ(posP);
1867   //
1868   
1869   do {    
1870     fVert.SetBC(AliVTrack::kTOFBCNA);
1871     VertexFitter(kTRUE,kFALSE,1);
1872     if (fVert.GetNContributors()<fMinTracks) {
1873       AliDebug(3,Form("Failed in iteration %d: No Contributirs",iter)); 
1874       break;
1875     } // failed
1876     if (fMVWSum>0) fMVSigma2 = TMath::Max(fMVWE2/fMVWSum,kSigLimit);
1877     else {
1878       AliDebug(3,Form("Failed %d, no weithgs",iter)); 
1879       iter = kMaxIter+1; 
1880       break;
1881     } // failed
1882     //
1883     double sigRed = (prevSig2-fMVSigma2)/prevSig2; // reduction of sigma2
1884     //
1885     fVert.GetXYZ(pos);
1886     double dst2 = (pos[0]-posP[0])*(pos[0]-posP[0])+(pos[1]-posP[1])*(pos[1]-posP[1])+(pos[2]-posP[2])*(pos[2]-posP[2]);
1887     AliDebug(3,Form("It#%2d Vtx: %+f %+f %+f Dst:%f Sig: %f [%.2f/%.2f] SigRed:%f",iter,pos[0],pos[1],pos[2],TMath::Sqrt(dst2),fMVSigma2,fMVWE2,fMVWSum,sigRed));
1888     if ( (++iter<kMaxIter) && (sigRed<0 || sigRed<fMVMinSig2Red) && fMVSigma2>fMVMaxSigma2) {
1889       fMVSigma2 *= push; // stuck, push little bit
1890       push *= kPushFactor;
1891       if (fMVSigma2<1.) fMVSigma2 = 1.; 
1892       AliDebug(3,Form("Pushed sigma2 to %f",fMVSigma2));
1893     }
1894     else if (dst2<minDst2 && ((sigRed<0 || sigRed<fMVMinSig2Red) || fMVSigma2<kSigLimitE)) break;
1895     //
1896     fVert.GetXYZ(posP); // fetch previous vertex position
1897     prevSig2 = fMVSigma2;
1898   } while(iter<kMaxIter);
1899   //
1900   if (fVert.GetNContributors()<0 || iter>kMaxIter || fMVSigma2>fMVMaxSigma2) {
1901     return kFALSE;
1902   }
1903   else {
1904     VertexFitter(kFALSE,kTRUE,fMVFinalWBinary ? 2:1); // final chi2 calculation
1905     int nv = fMVVertices->GetEntries();
1906     // create indices
1907     int ntrk = fTrkArraySel.GetEntries();
1908     int nindices = fCurrentVertex->GetNContributors() - (fConstraint ? 1:0);
1909     UShort_t *indices = 0;
1910     if (nindices>0) indices = new UShort_t[nindices];
1911     int nadded = 0;
1912     for (int itr=0;itr<ntrk;itr++) {
1913       AliExternalTrackParam* t = (AliExternalTrackParam*)fTrkArraySel[itr];
1914       if (t->TestBit(kBitAccounted) || !t->TestBit(kBitUsed)) continue;   // already belongs to some vertex
1915       t->SetBit(kBitAccounted);
1916       indices[nadded++] = fIdSel[itr];
1917     }
1918     if (nadded!=nindices) printf("Mismatch : NInd: %d Nadd: %d\n",nindices,nadded);
1919     fCurrentVertex->SetIndices(nindices,indices);
1920     // set vertex title
1921     TString title="VertexerTracksMVNoConstraint";
1922     if(fConstraint) title="VertexerTracksMVWithConstraint";
1923     fCurrentVertex->SetTitle(title.Data());    
1924     fMVVertices->AddLast(fCurrentVertex);
1925     AliDebug(3,Form("Added new vertex #%d NCont:%d XYZ: %f %f %f",nindices,nv,fCurrentVertex->GetX(),fCurrentVertex->GetY(),fCurrentVertex->GetZ()));
1926     if (indices) delete[] indices;
1927     fCurrentVertex = 0; // already attached to fMVVertices
1928     return kTRUE;
1929   }
1930 }
1931
1932 //______________________________________________________
1933 void AliVertexerTracks::FindVerticesMV()
1934 {
1935   // find and fit multiple vertices
1936   // 
1937   double step = fMVScanStep>1 ?  fMVScanStep : 1.;
1938   double zmx = 3*TMath::Sqrt(fNominalCov[5]);
1939   double zmn = -zmx;
1940   int nz = TMath::Nint((zmx-zmn)/step); if (nz<1) nz=1;
1941   double dz = (zmx-zmn)/nz;
1942   int izStart=0;
1943   AliDebug(2,Form("%d seeds between %f and %f",nz,zmn+dz/2,zmx+dz/2));
1944   //
1945   if (!fMVVertices) fMVVertices = new TObjArray(10);
1946   fMVVertices->Clear();
1947   //
1948   int ntrLeft = (Int_t)fTrkArraySel.GetEntries();
1949   //
1950   double sig2Scan = fMVSig2Ini;
1951   Bool_t runMore = kTRUE;
1952   int cntWide = 0;
1953   while (runMore) {
1954     fMVSig2Ini = sig2Scan*1e3;  // try wide search
1955     Bool_t found = kFALSE;
1956     cntWide++;
1957     fVert.SetNContributors(-1);
1958     fVert.SetXYZ(fNominalPos);
1959     AliDebug(3,Form("Wide search #%d Z= %f Sigma2=%f",cntWide,fNominalPos[2],fMVSig2Ini));
1960     if (FindNextVertexMV()) { // are there tracks left to consider?
1961       AliESDVertex* vtLast = (AliESDVertex*)fMVVertices->Last();
1962       if (vtLast && vtLast->GetNContributors()>0) ntrLeft -= vtLast->GetNContributors()-(fConstraint ? 1:0);
1963       if (ntrLeft<1) runMore = kFALSE; 
1964       found = kTRUE;
1965       continue;
1966     }  
1967     // if nothing is found, do narrow sig2ini scan
1968     fMVSig2Ini = sig2Scan;
1969     for (;izStart<nz;izStart++) {
1970       double zSeed = zmn+dz*(izStart+0.5);
1971       AliDebug(3,Form("Seed %d: Z= %f Sigma2=%f",izStart,zSeed,fMVSig2Ini));
1972       fVert.SetNContributors(-1);
1973       fVert.SetXYZ(fNominalPos);
1974       fVert.SetZv(zSeed);
1975       if (FindNextVertexMV()) { // are there tracks left to consider?
1976         AliESDVertex* vtLast = (AliESDVertex*)fMVVertices->Last();
1977         if (vtLast && vtLast->GetNContributors()>0) ntrLeft -= vtLast->GetNContributors()-(fConstraint ? 1:0);
1978         if (ntrLeft<1) runMore = kFALSE;
1979         found = kTRUE;
1980         break;
1981       }    
1982     }
1983     runMore = found; // if nothing was found, no need for new iteration
1984   }
1985   fMVSig2Ini = sig2Scan;
1986   int nvFound = fMVVertices->GetEntriesFast();
1987   AliDebug(2,Form("Number of found vertices: %d",nvFound));
1988   if (nvFound<1) TooFewTracks();
1989   if (AliLog::GetGlobalDebugLevel()>0) fMVVertices->Print();
1990   //
1991 }
1992
1993 //______________________________________________________
1994 void AliVertexerTracks::AnalyzePileUp(AliESDEvent* esdEv)
1995 {
1996   // if multiple vertices are found, try to find the primary one and attach it as fCurrentVertex
1997   // then attach pile-up vertices directly to esdEv
1998   //
1999   int nFND = (fMVVertices && fMVVertices->GetEntriesFast()) ? fMVVertices->GetEntriesFast() : 0;
2000   if (nFND<1) { if (!fCurrentVertex) TooFewTracks(); return;} // no multiple vertices
2001   //
2002   int indCont[nFND];
2003   int nIndx[nFND];
2004   for (int iv=0;iv<nFND;iv++) {
2005     AliESDVertex* fnd = (AliESDVertex*)fMVVertices->At(iv);
2006     indCont[iv] = iv;
2007     nIndx[iv]   = fnd->GetNIndices();
2008   }
2009   TMath::Sort(nFND, nIndx, indCont, kTRUE); // sort in decreasing order of Nindices
2010   double dists[nFND];
2011   int    distOrd[nFND],indx[nFND];
2012   for (int iv=0;iv<nFND;iv++) {
2013     AliESDVertex* fndI = (AliESDVertex*)fMVVertices->At(indCont[iv]);
2014     if (fndI->GetStatus()<1) continue; // discarded
2015     int ncomp = 0;
2016     for (int jv=iv+1;jv<nFND;jv++) {
2017       AliESDVertex* fndJ = (AliESDVertex*)fMVVertices->At(indCont[jv]);
2018       if (fndJ->GetStatus()<1) continue;
2019       dists[ncomp] = fndI->GetWDist(fndJ)*fndJ->GetNIndices();
2020       distOrd[ncomp] = indCont[jv];
2021       indx[ncomp]  = ncomp;
2022       ncomp++;
2023     }
2024     if (ncomp<1) continue;
2025     TMath::Sort(ncomp, dists, indx, kFALSE); // sort in increasing distance order
2026     for (int jv0=0;jv0<ncomp;jv0++) {
2027       int jv = distOrd[indx[jv0]];
2028       AliESDVertex* fndJ = (AliESDVertex*)fMVVertices->At(jv);
2029       if (dists[indx[jv0]]<fMVMaxWghNtr) { // candidate for split vertex
2030         //before eliminating the small close vertex, check if we should transfere its BC to largest one
2031         if (fndJ->GetBC()!=AliVTrack::kTOFBCNA && fndI->GetBC()==AliVTrack::kTOFBCNA) fndI->SetBC(fndJ->GetBC());
2032         //
2033         // leave the vertex only if both vertices have definit but different BC's, then this is not splitting.
2034         if ( (fndJ->GetBC()==fndI->GetBC()) || (fndJ->GetBC()==AliVTrack::kTOFBCNA)) fndJ->SetNContributors(-fndJ->GetNContributors());
2035       }
2036     }
2037   }
2038   //
2039   // select as a primary the largest vertex with BC=0, or the largest with BC non-ID
2040   int primBC0=-1,primNoBC=-1;
2041   for (int iv0=0;iv0<nFND;iv0++) {
2042     int iv = indCont[iv0];
2043     AliESDVertex* fndI = (AliESDVertex*)fMVVertices->At(iv);
2044     if (!fndI) continue;
2045     if (fndI->GetStatus()<1) {fMVVertices->RemoveAt(iv); delete fndI; continue;}
2046     if (primBC0<0  && fndI->GetBC()==0) primBC0 = iv;
2047     if (primNoBC<0 && fndI->GetBC()==AliVTrack::kTOFBCNA)  primNoBC = iv;
2048   }
2049   //
2050   if (primBC0>=0) fCurrentVertex = (AliESDVertex*)fMVVertices->At(primBC0);
2051   if (!fCurrentVertex && primNoBC>=0) fCurrentVertex = (AliESDVertex*)fMVVertices->At(primNoBC);
2052   if (fCurrentVertex) fMVVertices->Remove(fCurrentVertex);
2053   else {  // all vertices have BC>0, no primary vertex
2054     fCurrentVertex = new AliESDVertex();
2055     fCurrentVertex->SetNContributors(-3);
2056     fCurrentVertex->SetBC(AliVTrack::kTOFBCNA);
2057   }
2058   fCurrentVertex->SetID(-1);
2059   //
2060   // add pileup vertices
2061   int nadd = 0;
2062   for (int iv0=0;iv0<nFND;iv0++) {
2063     int iv = indCont[iv0];
2064     AliESDVertex* fndI = (AliESDVertex*)fMVVertices->At(iv);
2065     if (!fndI) continue;
2066     fndI->SetID(++nadd);
2067     esdEv->AddPileupVertexTracks(fndI);
2068   }
2069   //
2070   fMVVertices->Delete();
2071   //
2072 }
2073