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