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