]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STEER/AliVertexerTracks.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
further update on the HLT fast track finding
[u/mrichter/AliRoot.git] / STEER / AliVertexerTracks.cxx
CommitLineData
2d57349e 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//-----------------------------------------------------------------
6d8df534 17// Implementation of the vertexer from tracks
2d57349e 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:
6d8df534 25// F.Prino, Torino, prino@to.infn.it
2d57349e 26//-----------------------------------------------------------------
27
28//---- Root headers --------
a7ede274 29#include <TSystem.h>
30#include <TDirectory.h>
31#include <TFile.h>
2d57349e 32#include <TTree.h>
ec1be5d5 33#include <TMatrixD.h>
2d57349e 34//---- AliRoot headers -----
35#include "AliStrLine.h"
6d8df534 36#include "AliExternalTrackParam.h"
af885e0f 37#include "AliESDEvent.h"
2d57349e 38#include "AliESDtrack.h"
6d8df534 39#include "AliVertexerTracks.h"
2d57349e 40
41ClassImp(AliVertexerTracks)
42
43
44//----------------------------------------------------------------------------
45AliVertexerTracks::AliVertexerTracks():
07680cae 46TObject(),
47fVert(),
48fCurrentVertex(0),
f09c879d 49fMode(0),
dc3719f3 50fFieldkG(-999.),
6d8df534 51fTrkArraySel(),
52fIdSel(0),
07680cae 53fTrksToSkip(0),
54fNTrksToSkip(0),
f09c879d 55fConstraint(kFALSE),
56fOnlyFitter(kFALSE),
57fMinTracks(1),
58fMinITSClusters(5),
59fDCAcut(0.1),
60fDCAcutIter0(0.1),
61fNSigma(3.),
fee9ef0d 62fMaxd0z0(0.5),
f09c879d 63fMinDetFitter(100.),
64fMaxTgl(1000.),
50ff0bcd 65fITSrefit(kTRUE),
3b768a1d 66fnSigmaForUi00(1.5),
f09c879d 67fDebug(0),
68fAlgo(1)
2d57349e 69{
70//
71// Default constructor
72//
3b768a1d 73 SetVtxStart();
07680cae 74 SetVtxStartSigma();
f09c879d 75 SetITSMode();
2d57349e 76}
dc3719f3 77//----------------------------------------------------------------------------
78AliVertexerTracks::AliVertexerTracks(Double_t fieldkG):
07680cae 79TObject(),
80fVert(),
81fCurrentVertex(0),
f09c879d 82fMode(0),
83fFieldkG(fieldkG),
6d8df534 84fTrkArraySel(),
85fIdSel(0),
07680cae 86fTrksToSkip(0),
87fNTrksToSkip(0),
f09c879d 88fConstraint(kFALSE),
89fOnlyFitter(kFALSE),
90fMinTracks(1),
91fMinITSClusters(5),
92fDCAcut(0.1),
93fDCAcutIter0(0.1),
94fNSigma(3.),
fee9ef0d 95fMaxd0z0(0.5),
f09c879d 96fMinDetFitter(100.),
97fMaxTgl(1000.),
50ff0bcd 98fITSrefit(kTRUE),
3b768a1d 99fnSigmaForUi00(1.5),
f09c879d 100fDebug(0),
101fAlgo(1)
2d57349e 102{
103//
104// Standard constructor
105//
3b768a1d 106 SetVtxStart();
07680cae 107 SetVtxStartSigma();
f09c879d 108 SetITSMode();
2d57349e 109}
110//-----------------------------------------------------------------------------
fee9ef0d 111AliVertexerTracks::~AliVertexerTracks()
112{
2d57349e 113 // Default Destructor
eae1677e 114 // The objects pointed by the following pointer are not owned
2d57349e 115 // by this class and are not deleted
ec1be5d5 116 fCurrentVertex = 0;
146c29df 117 if(fTrksToSkip) { delete [] fTrksToSkip; fTrksToSkip=NULL; }
6d8df534 118 if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
ec1be5d5 119}
c5e3e5d1 120//----------------------------------------------------------------------------
af885e0f 121AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliESDEvent *esdEvent)
ec1be5d5 122{
c5e3e5d1 123//
ec1be5d5 124// Primary vertex for current ESD event
eae1677e 125// (Two iterations:
fee9ef0d 126// 1st with 5*fNSigma*sigma cut w.r.t. to initial vertex
127// + cut on sqrt(d0d0+z0z0) if fConstraint=kFALSE
128// 2nd with fNSigma*sigma cut w.r.t. to vertex found in 1st iteration)
07680cae 129//
130 fCurrentVertex = 0;
131
dc3719f3 132 // accept 1-track case only if constraint is available
133 if(!fConstraint && fMinTracks==1) fMinTracks=2;
134
fee9ef0d 135 // read tracks from ESD
6d8df534 136 Int_t nTrks = (Int_t)esdEvent->GetNumberOfTracks();
f09c879d 137 if(nTrks<fMinTracks) {
6d8df534 138 TooFewTracks();
80d36431 139 return fCurrentVertex;
140 }
6d8df534 141 //
713e5fbf 142
a7ede274 143 TDirectory * olddir = gDirectory;
f09c879d 144 TFile *f = 0;
145 if(nTrks>500) f = new TFile("VertexerTracks.root","recreate");
6d8df534 146 TObjArray trkArrayOrig(nTrks);
147 UShort_t *idOrig = new UShort_t[nTrks];
07680cae 148
6d8df534 149 Int_t nTrksOrig=0;
f09c879d 150 AliExternalTrackParam *t=0;
151 // loop on ESD tracks
6d8df534 152 for(Int_t i=0; i<nTrks; i++) {
153 AliESDtrack *esdt = esdEvent->GetTrack(i);
eae1677e 154 // check tracks to skip
6d8df534 155 Bool_t skipThis = kFALSE;
eae1677e 156 for(Int_t j=0; j<fNTrksToSkip; j++) {
6d8df534 157 if(esdt->GetID()==fTrksToSkip[j]) {
eae1677e 158 if(fDebug) printf("skipping track: %d\n",i);
159 skipThis = kTRUE;
160 }
161 }
a7ede274 162 if(skipThis) continue;
f09c879d 163
164 if(fMode==0) { // ITS mode
165 if(fITSrefit && !(esdt->GetStatus()&AliESDtrack::kITSrefit)) continue;
6d8df534 166 // check number of clusters in ITS
167 if(esdt->GetNcls(0)<fMinITSClusters) continue;
f09c879d 168 Double_t x,p[5],cov[15];
169 esdt->GetExternalParameters(x,p);
170 esdt->GetExternalCovariance(cov);
171 t = new AliExternalTrackParam(x,esdt->GetAlpha(),p,cov);
172 } else if(fMode==1) { // TPC mode
173 t = (AliExternalTrackParam*)esdt->GetTPCInnerParam();
174 if(!t) continue;
175 Double_t radius = 2.8; //something less than the beam pipe radius
176 if(!PropagateTrackTo(t,radius)) continue;
50ff0bcd 177 }
6d8df534 178 trkArrayOrig.AddLast(t);
179 idOrig[nTrksOrig]=(UShort_t)esdt->GetID();
180 nTrksOrig++;
f09c879d 181 } // end loop on ESD tracks
fee9ef0d 182
f09c879d 183 // call method that will reconstruct the vertex
6d8df534 184 FindPrimaryVertex(&trkArrayOrig,idOrig);
185
f09c879d 186 if(fMode==0) trkArrayOrig.Delete();
6d8df534 187 delete [] idOrig; idOrig=NULL;
f09c879d 188
189 if(f) {
190 f->Close(); delete f; f = NULL;
191 gSystem->Unlink("VertexerTracks.root");
192 olddir->cd();
193 }
6d8df534 194
195 return fCurrentVertex;
196}
197//----------------------------------------------------------------------------
198AliESDVertex* AliVertexerTracks::FindPrimaryVertex(TObjArray *trkArrayOrig,
199 UShort_t *idOrig)
200{
201//
202// Primary vertex using the AliExternalTrackParam's in the TObjArray.
203// idOrig must contain the track IDs from AliESDtrack::GetID()
204// (Two iterations:
205// 1st with 5*fNSigma*sigma cut w.r.t. to initial vertex
206// + cut on sqrt(d0d0+z0z0) if fConstraint=kFALSE
207// 2nd with fNSigma*sigma cut w.r.t. to vertex found in 1st iteration)
208//
209 fCurrentVertex = 0;
210
211 // accept 1-track case only if constraint is available
212 if(!fConstraint && fMinTracks==1) fMinTracks=2;
213
f09c879d 214 // read tracks from array
6d8df534 215 Int_t nTrksOrig = (Int_t)trkArrayOrig->GetEntriesFast();
43c9dae1 216 if(fDebug) printf("Initial number of tracks: %d\n",nTrksOrig);
f09c879d 217 if(nTrksOrig<fMinTracks) {
6d8df534 218 if(fDebug) printf("TooFewTracks\n");
219 TooFewTracks();
220 return fCurrentVertex;
221 }
222
fee9ef0d 223 // If fConstraint=kFALSE
224 // run VertexFinder(1) to get rough estimate of initVertex (x,y)
225 if(!fConstraint) {
6d8df534 226 // fill fTrkArraySel, for VertexFinder()
227 fIdSel = new UShort_t[nTrksOrig];
228 PrepareTracks(*trkArrayOrig,idOrig,0);
229 if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
230 Double_t cutsave = fDCAcut;
f09c879d 231 fDCAcut = fDCAcutIter0;
232 VertexFinder(1); // using weights, cutting dca < fDCAcutIter0
fee9ef0d 233 fDCAcut = cutsave;
fee9ef0d 234 if(fVert.GetNContributors()>0) {
235 fVert.GetXYZ(fNominalPos);
236 fNominalPos[0] = fVert.GetXv();
237 fNominalPos[1] = fVert.GetYv();
238 fNominalPos[2] = fVert.GetZv();
239 if(fDebug) printf("No mean vertex: VertexFinder gives (%f, %f, %f)\n",fNominalPos[0],fNominalPos[1],fNominalPos[2]);
240 } else {
241 fNominalPos[0] = 0.;
242 fNominalPos[1] = 0.;
243 fNominalPos[2] = 0.;
244 if(fDebug) printf("No mean vertex and VertexFinder failed\n");
0c69be49 245 }
07680cae 246 }
fee9ef0d 247
248 // TWO ITERATIONS:
249 //
250 // ITERATION 1
251 // propagate tracks to fNominalPos vertex
252 // preselect them:
253 // if(constraint) reject for |d0|>5*fNSigma*sigma w.r.t. fNominal... vertex
254 // else reject for |d0|\oplus|z0| > 5 mm w.r.t. fNominal... vertex
07680cae 255 // ITERATION 2
256 // propagate tracks to best between initVertex and fCurrentVertex
257 // preselect tracks (reject for |d0|>fNSigma*sigma w.r.t. best
258 // between initVertex and fCurrentVertex)
6d8df534 259 for(Int_t iter=1; iter<=2; iter++) {
260 if(fOnlyFitter && iter==1) continue;
261 if(fIdSel) { delete [] fIdSel; fIdSel=NULL; }
262 fIdSel = new UShort_t[nTrksOrig];
263 Int_t nTrksSel = PrepareTracks(*trkArrayOrig,idOrig,iter);
264 if(fDebug) printf("N tracks selected in iteration %d: %d\n",iter,nTrksSel);
265 if(nTrksSel < fMinTracks) {
266 TooFewTracks();
fee9ef0d 267 return fCurrentVertex;
268 }
07680cae 269
fee9ef0d 270 // vertex finder
271 if(!fOnlyFitter) {
6d8df534 272 if(nTrksSel==1) {
fee9ef0d 273 if(fDebug) printf("Just one track\n");
274 OneTrackVertFinder();
6d8df534 275 } else {
fee9ef0d 276 switch (fAlgo) {
277 case 1: StrLinVertexFinderMinDist(1); break;
278 case 2: StrLinVertexFinderMinDist(0); break;
279 case 3: HelixVertexFinder(); break;
280 case 4: VertexFinder(1); break;
281 case 5: VertexFinder(0); break;
282 default: printf("Wrong algorithm\n"); break;
283 }
284 }
285 if(fDebug) printf(" Vertex finding completed\n");
0c69be49 286 }
07680cae 287
fee9ef0d 288 // vertex fitter
6d8df534 289 VertexFitter();
fee9ef0d 290 } // end loop on the two iterations
07680cae 291
fee9ef0d 292
50ff0bcd 293 // set indices of used tracks
294 UShort_t *indices = 0;
50ff0bcd 295 if(fCurrentVertex->GetNContributors()>0) {
4d023a8f 296 Int_t nIndices = (Int_t)fTrkArraySel.GetEntriesFast();
297 indices = new UShort_t[nIndices];
298 for(Int_t jj=0; jj<nIndices; jj++)
6d8df534 299 indices[jj] = fIdSel[jj];
4d023a8f 300 fCurrentVertex->SetIndices(nIndices,indices);
50ff0bcd 301 }
6d8df534 302 delete [] indices; indices=NULL;
303 //
fee9ef0d 304
6d8df534 305 // set vertex title
306 TString title="VertexerTracksNoConstraint";
307 if(fConstraint) {
308 title="VertexerTracksWithConstraint";
309 if(fOnlyFitter) title.Append("OnlyFitter");
310 }
311 fCurrentVertex->SetTitle(title.Data());
312 //
a7ede274 313
6d8df534 314 if(fDebug) {
315 fCurrentVertex->PrintStatus();
316 fCurrentVertex->PrintIndices();
317 }
50ff0bcd 318
6d8df534 319 // clean up
320 delete [] fIdSel; fIdSel=NULL;
321 fTrkArraySel.Delete();
146c29df 322 if(fTrksToSkip) { delete [] fTrksToSkip; fTrksToSkip=NULL; }
6d8df534 323 //
fee9ef0d 324
07680cae 325 return fCurrentVertex;
326}
50ff0bcd 327//------------------------------------------------------------------------
fee9ef0d 328Double_t AliVertexerTracks::GetDeterminant3X3(Double_t matr[][3])
329{
50ff0bcd 330 //
331 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];
332 return det;
333}
334//-------------------------------------------------------------------------
fee9ef0d 335void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t (*m)[3],Double_t *d)
336{
50ff0bcd 337 //
338 Double_t x12=p0[0]-p1[0];
339 Double_t y12=p0[1]-p1[1];
340 Double_t z12=p0[2]-p1[2];
341 Double_t kk=x12*x12+y12*y12+z12*z12;
342 m[0][0]=2-2/kk*x12*x12;
343 m[0][1]=-2/kk*x12*y12;
344 m[0][2]=-2/kk*x12*z12;
345 m[1][0]=-2/kk*x12*y12;
346 m[1][1]=2-2/kk*y12*y12;
347 m[1][2]=-2/kk*y12*z12;
348 m[2][0]=-2/kk*x12*z12;
349 m[2][1]=-2*y12*z12;
350 m[2][2]=2-2/kk*z12*z12;
351 d[0]=2*p0[0]-2/kk*p0[0]*x12*x12-2/kk*p0[2]*x12*z12-2/kk*p0[1]*x12*y12;
352 d[1]=2*p0[1]-2/kk*p0[1]*y12*y12-2/kk*p0[0]*x12*y12-2/kk*p0[2]*z12*y12;
353 d[2]=2*p0[2]-2/kk*p0[2]*z12*z12-2/kk*p0[0]*x12*z12-2/kk*p0[1]*z12*y12;
c5e3e5d1 354
50ff0bcd 355}
356//--------------------------------------------------------------------------
fee9ef0d 357void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t *sigmasq,Double_t (*m)[3],Double_t *d)
358{
50ff0bcd 359 //
360 Double_t x12=p1[0]-p0[0];
361 Double_t y12=p1[1]-p0[1];
362 Double_t z12=p1[2]-p0[2];
c5e3e5d1 363
50ff0bcd 364 Double_t den= x12*x12*sigmasq[1]*sigmasq[2]+y12*y12*sigmasq[0]*sigmasq[2]+z12*z12*sigmasq[0]*sigmasq[1];
c5e3e5d1 365
50ff0bcd 366 Double_t kk= 2*(x12*x12/sigmasq[0]+y12*y12/sigmasq[1]+z12*z12/sigmasq[2]);
c5e3e5d1 367
50ff0bcd 368 Double_t cc[3];
369 cc[0]=-x12/sigmasq[0];
370 cc[1]=-y12/sigmasq[1];
371 cc[2]=-z12/sigmasq[2];
ec1be5d5 372
50ff0bcd 373 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;
ec1be5d5 374
50ff0bcd 375 Double_t ss= -p0[0]*cc[0]-p0[1]*cc[1]-p0[2]*cc[2];
2d57349e 376
50ff0bcd 377 Double_t aa[3];
378 aa[0]=x12*sigmasq[1]*sigmasq[2]/den;
379 aa[1]=y12*sigmasq[0]*sigmasq[2]/den;
380 aa[2]=z12*sigmasq[0]*sigmasq[1]/den;
2d57349e 381
50ff0bcd 382 m[0][0]=aa[0]*(aa[0]*kk+2*cc[0])+2*cc[0]*aa[0]+2/sigmasq[0];
383 m[0][1]=aa[1]*(aa[0]*kk+2*cc[0])+2*cc[1]*aa[0];
384 m[0][2]=aa[2]*(aa[0]*kk+2*cc[0])+2*cc[2]*aa[0];
2d57349e 385
50ff0bcd 386 m[1][0]=aa[0]*(aa[1]*kk+2*cc[1])+2*cc[0]*aa[1];
387 m[1][1]=aa[1]*(aa[1]*kk+2*cc[1])+2*cc[1]*aa[1]+2/sigmasq[1];
388 m[1][2]=aa[2]*(aa[1]*kk+2*cc[1])+2*cc[2]*aa[1];
2d57349e 389
50ff0bcd 390 m[2][0]=aa[0]*(aa[2]*kk+2*cc[2])+2*cc[0]*aa[2];
391 m[2][1]=aa[1]*(aa[2]*kk+2*cc[2])+2*cc[1]*aa[2];
392 m[2][2]=aa[2]*(aa[2]*kk+2*cc[2])+2*cc[2]*aa[2]+2/sigmasq[2];
08df6187 393
50ff0bcd 394 d[0]=-ww*(aa[0]*kk+2*cc[0])-2*ss*aa[0]+2*p0[0]/sigmasq[0];
395 d[1]=-ww*(aa[1]*kk+2*cc[1])-2*ss*aa[1]+2*p0[1]/sigmasq[1];
396 d[2]=-ww*(aa[2]*kk+2*cc[2])-2*ss*aa[2]+2*p0[2]/sigmasq[2];
2d57349e 397
2d57349e 398 }
50ff0bcd 399//--------------------------------------------------------------------------
fee9ef0d 400Double_t AliVertexerTracks::GetStrLinMinDist(Double_t *p0,Double_t *p1,Double_t *x0)
401{
50ff0bcd 402 //
403 Double_t x12=p0[0]-p1[0];
404 Double_t y12=p0[1]-p1[1];
405 Double_t z12=p0[2]-p1[2];
406 Double_t x10=p0[0]-x0[0];
407 Double_t y10=p0[1]-x0[1];
408 Double_t z10=p0[2]-x0[2];
3446476f 409 // 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);
3446476f 410 return ((y10*z12-z10*y12)*(y10*z12-z10*y12)+
411 (z10*x12-x10*z12)*(z10*x12-x10*z12)+
412 (x10*y12-y10*x12)*(x10*y12-y10*x12))
413 /(x12*x12+y12*y12+z12*z12);
2d57349e 414}
415//---------------------------------------------------------------------------
fee9ef0d 416void AliVertexerTracks::OneTrackVertFinder()
417{
0c69be49 418 // find vertex for events with 1 track, using DCA to nominal beam axis
6d8df534 419 if(fDebug) printf("Number of prepared tracks =%d - Call OneTrackVertFinder",fTrkArraySel.GetEntries());
420 AliExternalTrackParam *track1;
421 track1 = (AliExternalTrackParam*)fTrkArraySel.At(0);
0c69be49 422 Double_t alpha=track1->GetAlpha();
423 Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
424 Double_t pos[3],dir[3];
f09c879d 425 track1->GetXYZAt(mindist,GetFieldkG(),pos);
426 track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
0c69be49 427 AliStrLine *line1 = new AliStrLine(pos,dir);
428 Double_t p1[3]={fNominalPos[0],fNominalPos[1],0.};
429 Double_t p2[3]={fNominalPos[0],fNominalPos[1],10.};
430 AliStrLine *zeta=new AliStrLine(p1,p2,kTRUE);
431 Double_t crosspoint[3]={0.,0.,0.};
432 Double_t sigma=999.;
433 Int_t nContrib=-1;
434 Int_t retcode = zeta->Cross(line1,crosspoint);
435 if(retcode>=0){
436 sigma=line1->GetDistFromPoint(crosspoint);
437 nContrib=1;
438 }
439 delete zeta;
440 delete line1;
441 fVert.SetXYZ(crosspoint);
442 fVert.SetDispersion(sigma);
443 fVert.SetNContributors(nContrib);
444}
445//---------------------------------------------------------------------------
fee9ef0d 446void AliVertexerTracks::HelixVertexFinder()
447{
2d57349e 448 // Get estimate of vertex position in (x,y) from tracks DCA
449
450
451 Double_t initPos[3];
452 initPos[2] = 0.;
453 for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
2d57349e 454
6d8df534 455 Int_t nacc = (Int_t)fTrkArraySel.GetEntriesFast();
2d57349e 456
457 Double_t aver[3]={0.,0.,0.};
458 Double_t averquad[3]={0.,0.,0.};
459 Double_t sigmaquad[3]={0.,0.,0.};
460 Double_t sigma=0;
461 Int_t ncombi = 0;
6d8df534 462 AliExternalTrackParam *track1;
463 AliExternalTrackParam *track2;
2d57349e 464 Double_t distCA;
6d8df534 465 Double_t x;
2d57349e 466 Double_t alpha, cs, sn;
467 Double_t crosspoint[3];
468 for(Int_t i=0; i<nacc; i++){
6d8df534 469 track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
2d57349e 470
471
472 for(Int_t j=i+1; j<nacc; j++){
6d8df534 473 track2 = (AliExternalTrackParam*)fTrkArraySel.At(j);
2d57349e 474
f09c879d 475 distCA=track2->PropagateToDCA(track1,GetFieldkG());
2d57349e 476 if(fDCAcut<=0 ||(fDCAcut>0&&distCA<fDCAcut)){
6d8df534 477 x=track1->GetX();
2d57349e 478 alpha=track1->GetAlpha();
479 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
6d8df534 480 Double_t x1=x*cs - track1->GetY()*sn;
481 Double_t y1=x*sn + track1->GetY()*cs;
482 Double_t z1=track1->GetZ();
483
2d57349e 484 Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
6d8df534 485 x=track2->GetX();
2d57349e 486 alpha=track2->GetAlpha();
487 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
6d8df534 488 Double_t x2=x*cs - track2->GetY()*sn;
489 Double_t y2=x*sn + track2->GetY()*cs;
490 Double_t z2=track2->GetZ();
2d57349e 491 Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
492 Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
493 Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
494 Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
495 Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
496 crosspoint[0]=wx1*x1 + wx2*x2;
497 crosspoint[1]=wy1*y1 + wy2*y2;
498 crosspoint[2]=wz1*z1 + wz2*z2;
499
500 ncombi++;
501 for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
502 for(Int_t jj=0;jj<3;jj++)averquad[jj]+=(crosspoint[jj]*crosspoint[jj]);
503 }
504 }
505
506 }
507 if(ncombi>0){
508 for(Int_t jj=0;jj<3;jj++){
509 initPos[jj] = aver[jj]/ncombi;
510 averquad[jj]/=ncombi;
511 sigmaquad[jj]=averquad[jj]-initPos[jj]*initPos[jj];
512 sigma+=sigmaquad[jj];
513 }
514 sigma=TMath::Sqrt(TMath::Abs(sigma));
515 }
516 else {
517 Warning("HelixVertexFinder","Finder did not succed");
518 sigma=999;
519 }
520 fVert.SetXYZ(initPos);
521 fVert.SetDispersion(sigma);
522 fVert.SetNContributors(ncombi);
523}
50ff0bcd 524//----------------------------------------------------------------------------
6d8df534 525Int_t AliVertexerTracks::PrepareTracks(TObjArray &trkArrayOrig,
526 UShort_t *idOrig,
527 Int_t optImpParCut)
fee9ef0d 528{
50ff0bcd 529//
530// Propagate tracks to initial vertex position and store them in a TObjArray
531//
6d8df534 532 if(fDebug) printf(" PrepareTracks()\n");
533
534 Int_t nTrksOrig = (Int_t)trkArrayOrig.GetEntriesFast();
535 Int_t nTrksSel = 0;
fee9ef0d 536 Double_t maxd0rphi;
f09c879d 537 Double_t fiducialR = 3.0, fiducialZ = 20.0; // pipe and 1.5xSPD
50ff0bcd 538 Double_t sigmaCurr[3];
fee9ef0d 539 Double_t normdistx,normdisty;
6d8df534 540 Double_t d0z0[2],covd0z0[3];
43c9dae1 541 Double_t sigmad0;
50ff0bcd 542
543 AliESDVertex *initVertex = new AliESDVertex(fNominalPos,fNominalCov,1,1);
544
6d8df534 545 if(!fTrkArraySel.IsEmpty()) fTrkArraySel.Delete();
50ff0bcd 546
6d8df534 547 AliExternalTrackParam *track=0;
50ff0bcd 548
f09c879d 549 // loop on tracks
6d8df534 550 for(Int_t i=0; i<nTrksOrig; i++) {
551 track = new AliExternalTrackParam(*(AliExternalTrackParam*)trkArrayOrig.At(i));
f09c879d 552 // tgl cut
553 if(TMath::Abs(track->GetTgl())>fMaxTgl) {
43c9dae1 554 if(fDebug) printf(" rejecting track with tgl = %f\n",track->GetTgl());
f09c879d 555 delete track; continue;
6d8df534 556 }
50ff0bcd 557
558 // propagate track to vertex
43c9dae1 559 if(optImpParCut<2 || fOnlyFitter) { // optImpParCut==1 or 0
f09c879d 560 track->PropagateToDCA(initVertex,GetFieldkG(),100.,d0z0,covd0z0);
fee9ef0d 561 } else { // optImpParCut==2
50ff0bcd 562 fCurrentVertex->GetSigmaXYZ(sigmaCurr);
fee9ef0d 563 normdistx = TMath::Abs(fCurrentVertex->GetXv()-fNominalPos[0])/TMath::Sqrt(sigmaCurr[0]*sigmaCurr[0]+fNominalCov[0]);
80d36431 564 normdisty = TMath::Abs(fCurrentVertex->GetYv()-fNominalPos[1])/TMath::Sqrt(sigmaCurr[1]*sigmaCurr[1]+fNominalCov[2]);
fee9ef0d 565 if(normdistx < 3. && normdisty < 3. &&
80d36431 566 (sigmaCurr[0]+sigmaCurr[1])<(TMath::Sqrt(fNominalCov[0])+TMath::Sqrt(fNominalCov[2]))) {
f09c879d 567 track->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
50ff0bcd 568 } else {
f09c879d 569 track->PropagateToDCA(initVertex,GetFieldkG(),100.,d0z0,covd0z0);
50ff0bcd 570 }
571 }
572
43c9dae1 573 sigmad0 = TMath::Sqrt(covd0z0[0]);
574 maxd0rphi = fNSigma*sigmad0;
fee9ef0d 575 if(optImpParCut==1) maxd0rphi *= 5.;
f09c879d 576 maxd0rphi = TMath::Min(maxd0rphi,fiducialR);
43c9dae1 577 //sigmad0z0 = TMath::Sqrt(covd0z0[0]+covd0z0[2]); // for future improvement
578 //maxd0z0 = 10.*fNSigma*sigmad0z0;
fee9ef0d 579
f09c879d 580 if(fDebug) printf("trk %d; id %d; |d0| = %f; d0 cut = %f; |z0| = %f; |d0|oplus|z0| = %f; d0z0 cut = %f\n",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);
581
582
583 //---- track selection based on impact parameters ----//
584
585 // always reject tracks outside fiducial volume
586 if(TMath::Abs(d0z0[0])>fiducialR || TMath::Abs(d0z0[1])>fiducialZ) {
587 if(fDebug) printf(" rejected\n");
588 delete track; continue;
589 }
590
591 // during iterations 1 and 2 , reject tracks with d0rphi > maxd0rphi
592 if(optImpParCut>0 && TMath::Abs(d0z0[0])>maxd0rphi) {
593 if(fDebug) printf(" rejected\n");
594 delete track; continue;
595 }
fee9ef0d 596
43c9dae1 597 // if fConstraint=kFALSE, during iterations 1 and 2 ||
598 // if fConstraint=kTRUE, during iteration 2,
f09c879d 599 // select tracks with d0oplusz0 < fMaxd0z0
43c9dae1 600 if((!fConstraint && optImpParCut>0 && fVert.GetNContributors()>0) ||
601 ( fConstraint && optImpParCut==2)) {
602 if(nTrksOrig>=3 &&
f09c879d 603 TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1])>fMaxd0z0) {
fee9ef0d 604 if(fDebug) printf(" rejected\n");
605 delete track; continue;
606 }
607 }
43c9dae1 608
f09c879d 609 // track passed all selections
6d8df534 610 fTrkArraySel.AddLast(track);
611 fIdSel[nTrksSel] = idOrig[i];
612 nTrksSel++;
f09c879d 613 } // end loop on tracks
50ff0bcd 614
615 delete initVertex;
616
6d8df534 617 return nTrksSel;
50ff0bcd 618}
f09c879d 619//----------------------------------------------------------------------------
620Bool_t AliVertexerTracks::PropagateTrackTo(AliExternalTrackParam *track,
621 Double_t xToGo) {
622 //----------------------------------------------------------------
623 // COPIED from AliTracker
624 //
625 // Propagates the track to the plane X=xk (cm).
626 //
627 // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
628 //----------------------------------------------------------------
629
630 const Double_t kEpsilon = 0.00001;
631 Double_t xpos = track->GetX();
632 Double_t dir = (xpos<xToGo) ? 1. : -1.;
633 Double_t maxStep = 5;
634 Double_t maxSnp = 0.8;
635 //
636 while ( (xToGo-xpos)*dir > kEpsilon){
637 Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
638 Double_t x = xpos+step;
639 Double_t xyz0[3],xyz1[3];
640 track->GetXYZ(xyz0); //starting global position
641
642 if(!track->GetXYZAt(x,GetFieldkG(),xyz1)) return kFALSE; // no prolongation
643 xyz1[2]+=kEpsilon; // waiting for bug correction in geo
644
645 if(TMath::Abs(track->GetSnpAt(x,GetFieldkG())) >= maxSnp) return kFALSE;
646 if(!track->PropagateTo(x,GetFieldkG())) return kFALSE;
647
648 if(TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
649 track->GetXYZ(xyz0); // global position
650 Double_t alphan = TMath::ATan2(xyz0[1], xyz0[0]);
651 //
652 Double_t ca=TMath::Cos(alphan-track->GetAlpha()),
653 sa=TMath::Sin(alphan-track->GetAlpha());
654 Double_t sf=track->GetSnp(), cf=TMath::Sqrt(1.- sf*sf);
655 Double_t sinNew = sf*ca - cf*sa;
656 if(TMath::Abs(sinNew) >= maxSnp) return kFALSE;
657 if(!track->Rotate(alphan)) return kFALSE;
658
659 xpos = track->GetX();
660 }
661 return kTRUE;
662}
50ff0bcd 663//---------------------------------------------------------------------------
fee9ef0d 664AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
6d8df534 665 TObjArray *trkArray,
666 UShort_t *id,
fee9ef0d 667 Float_t *diamondxy)
668{
50ff0bcd 669//
fee9ef0d 670// Removes tracks in trksTree from fit of inVtx
50ff0bcd 671//
fee9ef0d 672
146c29df 673 if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraint")) {
674 printf("ERROR: primary vertex has no constraint: cannot remove tracks\n");
675 return 0x0;
676 }
fee9ef0d 677 if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
678 printf("WARNING: result of tracks' removal will be only approximately correct\n");
679
680 TMatrixD rv(3,1);
681 rv(0,0) = inVtx->GetXv();
682 rv(1,0) = inVtx->GetYv();
683 rv(2,0) = inVtx->GetZv();
684 TMatrixD vV(3,3);
685 Double_t cov[6];
686 inVtx->GetCovMatrix(cov);
687 vV(0,0) = cov[0];
688 vV(0,1) = cov[1]; vV(1,0) = cov[1];
689 vV(1,1) = cov[2];
690 vV(0,2) = cov[3]; vV(2,0) = cov[3];
691 vV(1,2) = cov[4]; vV(2,1) = cov[4];
692 vV(2,2) = cov[5];
693
694 TMatrixD sumWi(TMatrixD::kInverted,vV);
695 TMatrixD sumWiri(sumWi,TMatrixD::kMult,rv);
696
dad616ce 697 Int_t nUsedTrks = inVtx->GetNIndices();
fee9ef0d 698 Double_t chi2 = inVtx->GetChi2();
699
6d8df534 700 AliExternalTrackParam *track = 0;
701 Int_t ntrks = (Int_t)trkArray->GetEntriesFast();
fee9ef0d 702 for(Int_t i=0;i<ntrks;i++) {
6d8df534 703 track = (AliExternalTrackParam*)trkArray->At(i);
704 if(!inVtx->UsesTrack(id[i])) {
705 printf("track %d was not used in vertex fit\n",id[i]);
fee9ef0d 706 continue;
707 }
708 Double_t alpha = track->GetAlpha();
709 Double_t xl = diamondxy[0]*TMath::Cos(alpha)+diamondxy[1]*TMath::Sin(alpha);
6d8df534 710 track->PropagateTo(xl,GetFieldkG());
fee9ef0d 711 // vector of track global coordinates
712 TMatrixD ri(3,1);
713 // covariance matrix of ri
714 TMatrixD wWi(3,3);
715
716 // get space point from track
717 if(!TrackToPoint(track,ri,wWi)) continue;
718
719 TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
720
721 sumWi -= wWi;
722 sumWiri -= wWiri;
723
724 // track chi2
725 TMatrixD deltar = rv; deltar -= ri;
726 TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
727 Double_t chi2i = deltar(0,0)*wWideltar(0,0)+
728 deltar(1,0)*wWideltar(1,0)+
729 deltar(2,0)*wWideltar(2,0);
730 // remove from total chi2
731 chi2 -= chi2i;
732
733 nUsedTrks--;
734 if(nUsedTrks<2) {
735 printf("Trying to remove too many tracks!\n");
736 return 0x0;
737 }
738 }
739
740 TMatrixD rvnew(3,1);
741 TMatrixD vVnew(3,3);
742
743 // new inverted of weights matrix
744 TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
745 vVnew = invsumWi;
746 // new position of primary vertex
747 rvnew.Mult(vVnew,sumWiri);
748
749 Double_t position[3];
750 position[0] = rvnew(0,0);
751 position[1] = rvnew(1,0);
752 position[2] = rvnew(2,0);
753 cov[0] = vVnew(0,0);
754 cov[1] = vVnew(0,1);
755 cov[2] = vVnew(1,1);
756 cov[3] = vVnew(0,2);
757 cov[4] = vVnew(1,2);
758 cov[5] = vVnew(2,2);
759
760 // store data in the vertex object
761 AliESDVertex *outVtx = new AliESDVertex(position,cov,chi2,nUsedTrks);
762 outVtx->SetTitle(inVtx->GetTitle());
fee9ef0d 763 UShort_t *inindices = inVtx->GetIndices();
dad616ce 764 Int_t nIndices = nUsedTrks;
4d023a8f 765 UShort_t *outindices = new UShort_t[nIndices];
fee9ef0d 766 Int_t j=0;
fee9ef0d 767 for(Int_t k=0; k<inVtx->GetNIndices(); k++) {
6d8df534 768 Bool_t copyindex=kTRUE;
fee9ef0d 769 for(Int_t l=0; l<ntrks; l++) {
6d8df534 770 if(inindices[k]==id[l]) copyindex=kFALSE;
fee9ef0d 771 }
772 if(copyindex) {
773 outindices[j] = inindices[k]; j++;
774 }
775 }
4d023a8f 776 outVtx->SetIndices(nIndices,outindices);
fee9ef0d 777 delete [] outindices;
778
779 if(fDebug) {
780 printf("Vertex before removing tracks:\n");
781 inVtx->PrintStatus();
782 inVtx->PrintIndices();
783 printf("Vertex after removing tracks:\n");
784 outVtx->PrintStatus();
785 outVtx->PrintIndices();
786 }
787
788 return outVtx;
789}
790//---------------------------------------------------------------------------
f09c879d 791void AliVertexerTracks::SetITSMode(Double_t dcacut,
792 Double_t dcacutIter0,
793 Double_t maxd0z0,
794 Int_t minITScls,
795 Int_t mintrks,
796 Double_t nsigma,
797 Double_t mindetfitter,
798 Double_t maxtgl)
799{
800//
801// Cut values for ITS mode
802//
803 fMode = 0;
804 SetITSrefitRequired();
805 SetDCAcut(dcacut);
806 SetDCAcutIter0(dcacutIter0);
807 SetMaxd0z0(maxd0z0);
808 SetMinITSClusters(minITScls);
809 SetMinTracks(mintrks);
810 SetNSigmad0(nsigma);
811 SetMinDetFitter(mindetfitter);
812 SetMaxTgl(maxtgl);
813
814 return;
815}
816//---------------------------------------------------------------------------
817void AliVertexerTracks::SetTPCMode(Double_t dcacut,
818 Double_t dcacutIter0,
819 Double_t maxd0z0,
820 Int_t minITScls,
821 Int_t mintrks,
822 Double_t nsigma,
823 Double_t mindetfitter,
824 Double_t maxtgl)
825{
826//
827// Cut values for TPC mode
828//
829 fMode = 1;
830 SetITSrefitNotRequired();
831 SetDCAcut(dcacut);
832 SetDCAcutIter0(dcacutIter0);
833 SetMaxd0z0(maxd0z0);
834 SetMinITSClusters(minITScls);
835 SetMinTracks(mintrks);
836 SetNSigmad0(nsigma);
837 SetMinDetFitter(mindetfitter);
838 SetMaxTgl(maxtgl);
839
840 return;
841}
842//---------------------------------------------------------------------------
fee9ef0d 843void AliVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped)
844{
845//
846// Mark the tracks not to be used in the vertex reconstruction.
847// Tracks are identified by AliESDtrack::GetID()
848//
849 fNTrksToSkip = n; fTrksToSkip = new Int_t[n];
50ff0bcd 850 for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i];
851 return;
852}
853//---------------------------------------------------------------------------
fee9ef0d 854void AliVertexerTracks::SetVtxStart(AliESDVertex *vtx)
855{
50ff0bcd 856//
857// Set initial vertex knowledge
858//
859 vtx->GetXYZ(fNominalPos);
860 vtx->GetCovMatrix(fNominalCov);
fee9ef0d 861 SetConstraintOn();
50ff0bcd 862 return;
863}
2d57349e 864//---------------------------------------------------------------------------
fee9ef0d 865void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
dc3719f3 866{
6d8df534 867 AliExternalTrackParam *track1;
6d8df534 868 const Int_t knacc = (Int_t)fTrkArraySel.GetEntriesFast();
f54bad66 869 static TClonesArray linarray("AliStrLine",knacc);
dc3719f3 870 for(Int_t i=0; i<knacc; i++){
6d8df534 871 track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
dc3719f3 872 Double_t alpha=track1->GetAlpha();
873 Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
874 Double_t pos[3],dir[3],sigmasq[3];
f09c879d 875 track1->GetXYZAt(mindist,GetFieldkG(),pos);
876 track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
dc3719f3 877 sigmasq[0]=TMath::Sin(alpha)*TMath::Sin(alpha)*track1->GetSigmaY2();
878 sigmasq[1]=TMath::Cos(alpha)*TMath::Cos(alpha)*track1->GetSigmaY2();
879 sigmasq[2]=track1->GetSigmaZ2();
146c29df 880 TMatrixD ri(3,1);
881 TMatrixD wWi(3,3);
882 if(!TrackToPoint(track1,ri,wWi)) continue;
883 Double_t wmat[9];
884 Int_t iel=0;
885 for(Int_t ia=0;ia<3;ia++){
886 for(Int_t ib=0;ib<3;ib++){
887 wmat[iel]=wWi(ia,ib);
888 iel++;
889 }
890 }
f54bad66 891 new(linarray[i]) AliStrLine(pos,sigmasq,wmat,dir);
dc3719f3 892 }
f54bad66 893 fVert=TrackletVertexFinder(&linarray,optUseWeights);
894 linarray.Clear("C");
dc3719f3 895}
896//---------------------------------------------------------------------------
146c29df 897AliESDVertex AliVertexerTracks::TrackletVertexFinder(TClonesArray *lines, Int_t optUseWeights)
fee9ef0d 898{
2d57349e 899 // Calculate the point at minimum distance to prepared tracks
900
dc3719f3 901 const Int_t knacc = (Int_t)lines->GetEntriesFast();
902 Double_t initPos[3]={0.,0.,0.};
2d57349e 903
2d57349e 904 Double_t (*vectP0)[3]=new Double_t [knacc][3];
905 Double_t (*vectP1)[3]=new Double_t [knacc][3];
906
907 Double_t sum[3][3];
908 Double_t dsum[3]={0,0,0};
146c29df 909 TMatrixD sumWi(3,3);
910 for(Int_t i=0;i<3;i++){
911 for(Int_t j=0;j<3;j++){
912 sum[i][j]=0;
913 sumWi(i,j)=0.;
914 }
915 }
916
2d57349e 917 for(Int_t i=0; i<knacc; i++){
f54bad66 918 AliStrLine* line1 = (AliStrLine*)lines->At(i);
dc3719f3 919 Double_t p0[3],cd[3],sigmasq[3];
146c29df 920 Double_t wmat[9];
2d57349e 921 line1->GetP0(p0);
922 line1->GetCd(cd);
dc3719f3 923 line1->GetSigma2P0(sigmasq);
146c29df 924 line1->GetWMatrix(wmat);
925 TMatrixD wWi(3,3);
926 Int_t iel=0;
927 for(Int_t ia=0;ia<3;ia++){
928 for(Int_t ib=0;ib<3;ib++){
929 wWi(ia,ib)=wmat[iel];
930 iel++;
931 }
932 }
933
934 sumWi+=wWi;
935
2d57349e 936 Double_t p1[3]={p0[0]+cd[0],p0[1]+cd[1],p0[2]+cd[2]};
937 vectP0[i][0]=p0[0];
938 vectP0[i][1]=p0[1];
939 vectP0[i][2]=p0[2];
940 vectP1[i][0]=p1[0];
941 vectP1[i][1]=p1[1];
942 vectP1[i][2]=p1[2];
943
944 Double_t matr[3][3];
945 Double_t dknow[3];
946 if(optUseWeights==0)GetStrLinDerivMatrix(p0,p1,matr,dknow);
dc3719f3 947 else GetStrLinDerivMatrix(p0,p1,sigmasq,matr,dknow);
948
2d57349e 949
950 for(Int_t iii=0;iii<3;iii++){
951 dsum[iii]+=dknow[iii];
952 for(Int_t lj=0;lj<3;lj++) sum[iii][lj]+=matr[iii][lj];
953 }
2d57349e 954 }
146c29df 955
956 TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
957 Double_t covmatrix[6];
958 covmatrix[0] = invsumWi(0,0);
959 covmatrix[1] = invsumWi(0,1);
960 covmatrix[2] = invsumWi(1,1);
961 covmatrix[3] = invsumWi(0,2);
962 covmatrix[4] = invsumWi(1,2);
963 covmatrix[5] = invsumWi(2,2);
964
2d57349e 965 Double_t vett[3][3];
966 Double_t det=GetDeterminant3X3(sum);
dc3719f3 967 Double_t sigma=0;
2d57349e 968
fee9ef0d 969 if(det!=0){
970 for(Int_t zz=0;zz<3;zz++){
971 for(Int_t ww=0;ww<3;ww++){
972 for(Int_t kk=0;kk<3;kk++) vett[ww][kk]=sum[ww][kk];
973 }
974 for(Int_t kk=0;kk<3;kk++) vett[kk][zz]=dsum[kk];
975 initPos[zz]=GetDeterminant3X3(vett)/det;
976 }
977
978
979 for(Int_t i=0; i<knacc; i++){
980 Double_t p0[3]={0,0,0},p1[3]={0,0,0};
981 for(Int_t ii=0;ii<3;ii++){
982 p0[ii]=vectP0[i][ii];
983 p1[ii]=vectP1[i][ii];
984 }
985 sigma+=GetStrLinMinDist(p0,p1,initPos);
986 }
987
f09c879d 988 if(sigma>0.) {sigma=TMath::Sqrt(sigma);}else{sigma=999;}
fee9ef0d 989 }else{
50ff0bcd 990 sigma=999;
991 }
146c29df 992 AliESDVertex theVert(initPos,covmatrix,99999.,knacc);
993 theVert.SetDispersion(sigma);
50ff0bcd 994 delete vectP0;
995 delete vectP1;
dc3719f3 996 return theVert;
50ff0bcd 997}
998//---------------------------------------------------------------------------
6d8df534 999Bool_t AliVertexerTracks::TrackToPoint(AliExternalTrackParam *t,
0bf2e2b4 1000 TMatrixD &ri,TMatrixD &wWi,
1001 Bool_t uUi3by3) const
fee9ef0d 1002{
1003//
6d8df534 1004// Extract from the AliExternalTrackParam the global coordinates ri and covariance matrix
fee9ef0d 1005// wWi of the space point that it represents (to be used in VertexFitter())
1006//
1007
1008
1009 Double_t rotAngle = t->GetAlpha();
1010 if(rotAngle<0.) rotAngle += 2.*TMath::Pi();
1011 Double_t cosRot = TMath::Cos(rotAngle);
1012 Double_t sinRot = TMath::Sin(rotAngle);
1013
1014 ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot;
1015 ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot;
1016 ri(2,0) = t->GetZ();
1017
0bf2e2b4 1018
1019
1020 if(!uUi3by3) {
1021 // matrix to go from global (x,y,z) to local (y,z);
1022 TMatrixD qQi(2,3);
1023 qQi(0,0) = -sinRot;
1024 qQi(0,1) = cosRot;
1025 qQi(0,2) = 0.;
1026 qQi(1,0) = 0.;
1027 qQi(1,1) = 0.;
1028 qQi(1,2) = 1.;
1029
1030 // covariance matrix of local (y,z) - inverted
0bf2e2b4 1031 TMatrixD uUi(2,2);
6d8df534 1032 uUi(0,0) = t->GetSigmaY2();
1033 uUi(0,1) = t->GetSigmaZY();
1034 uUi(1,0) = t->GetSigmaZY();
1035 uUi(1,1) = t->GetSigmaZ2();
0bf2e2b4 1036 //printf(" Ui :\n");
1037 //printf(" %f %f\n",uUi(0,0),uUi(0,1));
1038 //printf(" %f %f\n",uUi(1,0),uUi(1,1));
1039
1040 if(uUi.Determinant() <= 0.) return kFALSE;
1041 TMatrixD uUiInv(TMatrixD::kInverted,uUi);
1042
1043 // weights matrix: wWi = qQiT * uUiInv * qQi
1044 TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
1045 TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
1046 wWi = m;
1047 } else {
1048 if(fVert.GetNContributors()<1) AliFatal("Vertex from finder is empty");
1049 // matrix to go from global (x,y,z) to local (x,y,z);
1050 TMatrixD qQi(3,3);
1051 qQi(0,0) = cosRot;
1052 qQi(0,1) = sinRot;
1053 qQi(0,2) = 0.;
1054 qQi(1,0) = -sinRot;
1055 qQi(1,1) = cosRot;
1056 qQi(1,2) = 0.;
1057 qQi(2,0) = 0.;
1058 qQi(2,1) = 0.;
1059 qQi(2,2) = 1.;
1060
1061 // covariance of fVert along the track
1062 Double_t p[3],pt,ptot;
1063 t->GetPxPyPz(p);
1064 pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
1065 ptot = TMath::Sqrt(pt*pt+p[2]*p[2]);
3b768a1d 1066 Double_t cphi = p[0]/pt; //cos(phi)=px/pt
1067 Double_t sphi = p[1]/pt; //sin(phi)=py/pt
1068 Double_t clambda = pt/ptot; //cos(lambda)=pt/ptot
1069 Double_t slambda = p[2]/ptot; //sin(lambda)=pz/ptot
0bf2e2b4 1070 Double_t covfVert[6];
1071 fVert.GetCovMatrix(covfVert);
1072 Double_t covfVertalongt =
3b768a1d 1073 covfVert[0]*cphi*cphi*clambda*clambda
1074 +covfVert[1]*2.*cphi*sphi*clambda*clambda
1075 +covfVert[3]*2.*cphi*clambda*slambda
1076 +covfVert[2]*sphi*sphi*clambda*clambda
1077 +covfVert[4]*2.*sphi*clambda*slambda
1078 +covfVert[5]*slambda*slambda;
0bf2e2b4 1079 // covariance matrix of local (x,y,z) - inverted
1080 TMatrixD uUi(3,3);
3b768a1d 1081 uUi(0,0) = covfVertalongt * fnSigmaForUi00 * fnSigmaForUi00;
0bf2e2b4 1082 if(fDebug) printf("=====> sqrtUi00 cm %f\n",TMath::Sqrt(uUi(0,0)));
1083 uUi(0,1) = 0.;
1084 uUi(0,2) = 0.;
1085 uUi(1,0) = 0.;
6d8df534 1086 uUi(1,1) = t->GetSigmaY2();
1087 uUi(1,2) = t->GetSigmaZY();
0bf2e2b4 1088 uUi(2,0) = 0.;
6d8df534 1089 uUi(2,1) = t->GetSigmaZY();
1090 uUi(2,2) = t->GetSigmaZ2();
0bf2e2b4 1091 //printf(" Ui :\n");
1092 //printf(" %f %f\n",uUi(0,0),uUi(0,1));
1093 //printf(" %f %f\n",uUi(1,0),uUi(1,1));
fee9ef0d 1094
0bf2e2b4 1095 if(uUi.Determinant() <= 0.) return kFALSE;
1096 TMatrixD uUiInv(TMatrixD::kInverted,uUi);
1097
1098 // weights matrix: wWi = qQiT * uUiInv * qQi
1099 TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
1100 TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
1101 wWi = m;
1102 }
1103
fee9ef0d 1104
1105 return kTRUE;
1106}
1107//---------------------------------------------------------------------------
6d8df534 1108void AliVertexerTracks::TooFewTracks()
fee9ef0d 1109{
50ff0bcd 1110//
6d8df534 1111// When the number of tracks is < fMinTracks,
1112// deal with vertices not found and prepare to exit
50ff0bcd 1113//
6d8df534 1114 if(fDebug) printf("TooFewTracks\n");
2d57349e 1115
50ff0bcd 1116 Double_t pos[3],err[3];
50ff0bcd 1117 pos[0] = fNominalPos[0];
1118 err[0] = TMath::Sqrt(fNominalCov[0]);
1119 pos[1] = fNominalPos[1];
1120 err[1] = TMath::Sqrt(fNominalCov[2]);
6d8df534 1121 pos[2] = fNominalPos[2];
1122 err[2] = TMath::Sqrt(fNominalCov[5]);
1123 Int_t ncontr = (err[0]>1. ? -1 : -3);
50ff0bcd 1124 fCurrentVertex = 0;
1125 fCurrentVertex = new AliESDVertex(pos,err);
1126 fCurrentVertex->SetNContributors(ncontr);
2d57349e 1127
fee9ef0d 1128 if(fConstraint) {
1129 fCurrentVertex->SetTitle("VertexerTracksWithConstraint");
1130 } else {
50ff0bcd 1131 fCurrentVertex->SetTitle("VertexerTracksNoConstraint");
fee9ef0d 1132 }
2d57349e 1133
6d8df534 1134 if(!fTrkArraySel.IsEmpty()) fTrkArraySel.Delete();
1135 if(fIdSel) {delete [] fIdSel; fIdSel=NULL;}
1136 if(fTrksToSkip) {delete [] fTrksToSkip; fTrksToSkip=NULL;}
1137
50ff0bcd 1138 return;
1139}
1140//---------------------------------------------------------------------------
fee9ef0d 1141void AliVertexerTracks::VertexFinder(Int_t optUseWeights)
1142{
2d57349e 1143
50ff0bcd 1144 // Get estimate of vertex position in (x,y) from tracks DCA
1145
1146 Double_t initPos[3];
1147 initPos[2] = 0.;
1148 for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
6d8df534 1149 Int_t nacc = (Int_t)fTrkArraySel.GetEntriesFast();
50ff0bcd 1150 Double_t aver[3]={0.,0.,0.};
1151 Double_t aversq[3]={0.,0.,0.};
1152 Double_t sigmasq[3]={0.,0.,0.};
1153 Double_t sigma=0;
1154 Int_t ncombi = 0;
6d8df534 1155 AliExternalTrackParam *track1;
1156 AliExternalTrackParam *track2;
50ff0bcd 1157 Double_t pos[3],dir[3];
1158 Double_t alpha,mindist;
2d57349e 1159
50ff0bcd 1160 for(Int_t i=0; i<nacc; i++){
6d8df534 1161 track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
50ff0bcd 1162 alpha=track1->GetAlpha();
1163 mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
f09c879d 1164 track1->GetXYZAt(mindist,GetFieldkG(),pos);
1165 track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
50ff0bcd 1166 AliStrLine *line1 = new AliStrLine(pos,dir);
2d57349e 1167
50ff0bcd 1168 // AliStrLine *line1 = new AliStrLine();
f09c879d 1169 // track1->ApproximateHelixWithLine(mindist,GetFieldkG(),line1);
50ff0bcd 1170
1171 for(Int_t j=i+1; j<nacc; j++){
6d8df534 1172 track2 = (AliExternalTrackParam*)fTrkArraySel.At(j);
50ff0bcd 1173 alpha=track2->GetAlpha();
1174 mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
f09c879d 1175 track2->GetXYZAt(mindist,GetFieldkG(),pos);
1176 track2->GetPxPyPzAt(mindist,GetFieldkG(),dir);
50ff0bcd 1177 AliStrLine *line2 = new AliStrLine(pos,dir);
1178 // AliStrLine *line2 = new AliStrLine();
f09c879d 1179 // track2->ApproximateHelixWithLine(mindist,GetFieldkG(),line2);
50ff0bcd 1180 Double_t distCA=line2->GetDCA(line1);
fee9ef0d 1181 //printf("%d %d %f\n",i,j,distCA);
1182 if(fDCAcut<=0 || (fDCAcut>0&&distCA<fDCAcut)){
50ff0bcd 1183 Double_t pnt1[3],pnt2[3],crosspoint[3];
2d57349e 1184
50ff0bcd 1185 if(optUseWeights<=0){
1186 Int_t retcode = line2->Cross(line1,crosspoint);
1187 if(retcode>=0){
1188 ncombi++;
1189 for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
1190 for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
1191 }
1192 }
1193 if(optUseWeights>0){
1194 Int_t retcode = line1->CrossPoints(line2,pnt1,pnt2);
1195 if(retcode>=0){
5c03260f 1196 Double_t cs, sn;
50ff0bcd 1197 alpha=track1->GetAlpha();
1198 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
1199 Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
1200 alpha=track2->GetAlpha();
1201 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
1202 Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
1203 Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
1204 Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
1205 Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
1206 Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
1207 crosspoint[0]=wx1*pnt1[0] + wx2*pnt2[0];
1208 crosspoint[1]=wy1*pnt1[1] + wy2*pnt2[1];
1209 crosspoint[2]=wz1*pnt1[2] + wz2*pnt2[2];
1210
1211 ncombi++;
1212 for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
1213 for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
1214 }
1215 }
1216 }
1217 delete line2;
1218 }
1219 delete line1;
71d84967 1220 }
50ff0bcd 1221 if(ncombi>0){
1222 for(Int_t jj=0;jj<3;jj++){
1223 initPos[jj] = aver[jj]/ncombi;
fee9ef0d 1224 //printf("%f\n",initPos[jj]);
50ff0bcd 1225 aversq[jj]/=ncombi;
1226 sigmasq[jj]=aversq[jj]-initPos[jj]*initPos[jj];
1227 sigma+=sigmasq[jj];
1228 }
1229 sigma=TMath::Sqrt(TMath::Abs(sigma));
1230 }
1231 else {
1232 Warning("VertexFinder","Finder did not succed");
1233 sigma=999;
1234 }
1235 fVert.SetXYZ(initPos);
1236 fVert.SetDispersion(sigma);
1237 fVert.SetNContributors(ncombi);
2d57349e 1238}
ec1be5d5 1239//---------------------------------------------------------------------------
6d8df534 1240void AliVertexerTracks::VertexFitter()
fee9ef0d 1241{
ec1be5d5 1242//
1243// The optimal estimate of the vertex position is given by a "weighted
0bf2e2b4 1244// average of tracks positions".
fee9ef0d 1245// Original method: V. Karimaki, CMS Note 97/0051
ec1be5d5 1246//
6d8df534 1247 Bool_t useConstraint = fConstraint;
ec1be5d5 1248 Double_t initPos[3];
0bf2e2b4 1249 if(!fOnlyFitter) {
1250 fVert.GetXYZ(initPos);
1251 } else {
1252 initPos[0]=fNominalPos[0];
1253 initPos[1]=fNominalPos[1];
1254 initPos[2]=fNominalPos[2];
1255 }
1256
6d8df534 1257 Int_t nTrksSel = (Int_t)fTrkArraySel.GetEntries();
1258 if(nTrksSel==1) useConstraint=kTRUE;
ec1be5d5 1259 if(fDebug) {
0bf2e2b4 1260 printf("--- VertexFitter(): start\n");
6d8df534 1261 printf(" Number of tracks in array: %d\n",nTrksSel);
ec1be5d5 1262 printf(" Minimum # tracks required in fit: %d\n",fMinTracks);
3b768a1d 1263 printf(" Vertex position after finder: %f,%f,%f\n",initPos[0],initPos[1],initPos[2]);
fee9ef0d 1264 if(useConstraint) printf(" This vertex will be used in fit: (%f+-%f,%f+-%f)\n",fNominalPos[0],TMath::Sqrt(fNominalCov[0]),fNominalPos[1],TMath::Sqrt(fNominalCov[2]));
ec1be5d5 1265 }
1266
0bf2e2b4 1267 // special treatment for few-tracks fits (e.g. secondary vertices)
6d8df534 1268 Bool_t uUi3by3 = kFALSE; if(nTrksSel<5 && !useConstraint) uUi3by3 = kTRUE;
0bf2e2b4 1269
ec1be5d5 1270 Int_t i,j,k,step=0;
1271 TMatrixD rv(3,1);
1272 TMatrixD vV(3,3);
1273 rv(0,0) = initPos[0];
1274 rv(1,0) = initPos[1];
1275 rv(2,0) = 0.;
1276 Double_t xlStart,alpha;
6d8df534 1277 Int_t nTrksUsed;
fee9ef0d 1278 Double_t chi2,chi2i,chi2b;
6d8df534 1279 AliExternalTrackParam *t = 0;
ec1be5d5 1280 Int_t failed = 0;
1281
50ff0bcd 1282 // initial vertex covariance matrix
1283 TMatrixD vVb(3,3);
1284 vVb(0,0) = fNominalCov[0];
1285 vVb(0,1) = fNominalCov[1];
1286 vVb(0,2) = 0.;
1287 vVb(1,0) = fNominalCov[1];
1288 vVb(1,1) = fNominalCov[2];
1289 vVb(1,2) = 0.;
1290 vVb(2,0) = 0.;
1291 vVb(2,1) = 0.;
1292 vVb(2,2) = fNominalCov[5];
1293 TMatrixD vVbInv(TMatrixD::kInverted,vVb);
1294 TMatrixD rb(3,1);
1295 rb(0,0) = fNominalPos[0];
1296 rb(1,0) = fNominalPos[1];
1297 rb(2,0) = fNominalPos[2];
1298 TMatrixD vVbInvrb(vVbInv,TMatrixD::kMult,rb);
1299
1300
eae1677e 1301 // 2 steps:
1302 // 1st - estimate of vtx using all tracks
1303 // 2nd - estimate of global chi2
1304 for(step=0; step<2; step++) {
ec1be5d5 1305 if(fDebug) printf(" step = %d\n",step);
1306 chi2 = 0.;
6d8df534 1307 nTrksUsed = 0;
ec1be5d5 1308
0bf2e2b4 1309 if(step==1) { initPos[0]=rv(0,0); initPos[0]=rv(1,0); }
1310
ec1be5d5 1311 TMatrixD sumWiri(3,1);
1312 TMatrixD sumWi(3,3);
1313 for(i=0; i<3; i++) {
1314 sumWiri(i,0) = 0.;
1315 for(j=0; j<3; j++) sumWi(j,i) = 0.;
1316 }
1317
fee9ef0d 1318 // mean vertex constraint
1319 if(useConstraint) {
50ff0bcd 1320 for(i=0;i<3;i++) {
1321 sumWiri(i,0) += vVbInvrb(i,0);
1322 for(k=0;k<3;k++) sumWi(i,k) += vVbInv(i,k);
07680cae 1323 }
fee9ef0d 1324 // chi2
1325 TMatrixD deltar = rv; deltar -= rb;
1326 TMatrixD vVbInvdeltar(vVbInv,TMatrixD::kMult,deltar);
1327 chi2b = deltar(0,0)*vVbInvdeltar(0,0)+
1328 deltar(1,0)*vVbInvdeltar(1,0)+
1329 deltar(2,0)*vVbInvdeltar(2,0);
1330 chi2 += chi2b;
07680cae 1331 }
1332
ec1be5d5 1333 // loop on tracks
6d8df534 1334 for(k=0; k<nTrksSel; k++) {
ec1be5d5 1335 // get track from track array
6d8df534 1336 t = (AliExternalTrackParam*)fTrkArraySel.At(k);
ec1be5d5 1337 alpha = t->GetAlpha();
1338 xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha);
fee9ef0d 1339 // to vtxSeed (from finder)
6d8df534 1340 t->PropagateTo(xlStart,GetFieldkG());
fee9ef0d 1341
ec1be5d5 1342 // vector of track global coordinates
1343 TMatrixD ri(3,1);
fee9ef0d 1344 // covariance matrix of ri
1345 TMatrixD wWi(3,3);
1346
1347 // get space point from track
0bf2e2b4 1348 if(!TrackToPoint(t,ri,wWi,uUi3by3)) continue;
fee9ef0d 1349 TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
ec1be5d5 1350
1351 // track chi2
1352 TMatrixD deltar = rv; deltar -= ri;
1353 TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
1354 chi2i = deltar(0,0)*wWideltar(0,0)+
1355 deltar(1,0)*wWideltar(1,0)+
1356 deltar(2,0)*wWideltar(2,0);
1357
ec1be5d5 1358 // add to total chi2
1359 chi2 += chi2i;
1360
ec1be5d5 1361 sumWiri += wWiri;
1362 sumWi += wWi;
1363
6d8df534 1364 nTrksUsed++;
ec1be5d5 1365 } // end loop on tracks
1366
6d8df534 1367 if(nTrksUsed < fMinTracks) {
ec1be5d5 1368 failed=1;
1369 continue;
1370 }
1371
1372 Double_t determinant = sumWi.Determinant();
f09c879d 1373 if(determinant < fMinDetFitter) {
1374 printf("det(V) = %f (<%f)\n",determinant,fMinDetFitter);
ec1be5d5 1375 failed=1;
1376 continue;
1377 }
1378
0bf2e2b4 1379 if(step==0) {
1380 // inverted of weights matrix
1381 TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
1382 vV = invsumWi;
1383 // position of primary vertex
1384 rv.Mult(vV,sumWiri);
1385 }
eae1677e 1386 } // end loop on the 2 steps
ec1be5d5 1387
ec1be5d5 1388 if(failed) {
6d8df534 1389 TooFewTracks();
ec1be5d5 1390 return;
1391 }
1392
1393 Double_t position[3];
1394 position[0] = rv(0,0);
1395 position[1] = rv(1,0);
1396 position[2] = rv(2,0);
1397 Double_t covmatrix[6];
1398 covmatrix[0] = vV(0,0);
1399 covmatrix[1] = vV(0,1);
1400 covmatrix[2] = vV(1,1);
1401 covmatrix[3] = vV(0,2);
1402 covmatrix[4] = vV(1,2);
1403 covmatrix[5] = vV(2,2);
1404
dc3719f3 1405 // for correct chi2/ndf, count constraint as additional "track"
6d8df534 1406 if(fConstraint) nTrksUsed++;
ec1be5d5 1407 // store data in the vertex object
6d8df534 1408 fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);
ec1be5d5 1409
1410 if(fDebug) {
0bf2e2b4 1411 printf(" Vertex after fit:\n");
ec1be5d5 1412 fCurrentVertex->PrintStatus();
0bf2e2b4 1413 printf("--- VertexFitter(): finish\n");
ec1be5d5 1414 }
1415
1416 return;
1417}
50ff0bcd 1418//----------------------------------------------------------------------------
6d8df534 1419AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray,
1420 UShort_t *id,
1421 Bool_t optUseFitter,
1422 Bool_t optPropagate)
fee9ef0d 1423{
eae1677e 1424//
6d8df534 1425// Return vertex from tracks (AliExternalTrackParam) in array
eae1677e 1426//
1427
fee9ef0d 1428 SetConstraintOff();
1429
50ff0bcd 1430 // get tracks and propagate them to initial vertex position
6d8df534 1431 fIdSel = new UShort_t[(Int_t)trkArray->GetEntriesFast()];
1432 Int_t nTrksSel = PrepareTracks(*trkArray,id,0);
1433 if(nTrksSel < TMath::Max(2,fMinTracks)) {
1434 TooFewTracks();
fee9ef0d 1435 return fCurrentVertex;
50ff0bcd 1436 }
1437
6d8df534 1438 // vertex finder
fee9ef0d 1439 switch (fAlgo) {
1440 case 1: StrLinVertexFinderMinDist(1); break;
1441 case 2: StrLinVertexFinderMinDist(0); break;
1442 case 3: HelixVertexFinder(); break;
1443 case 4: VertexFinder(1); break;
1444 case 5: VertexFinder(0); break;
1445 default: printf("Wrong algorithm\n"); break;
1446 }
fee9ef0d 1447 if(fDebug) printf(" Vertex finding completed\n");
1448
1449 // vertex fitter
6d8df534 1450 if(optUseFitter) {
1451 VertexFitter();
1452 } else {
fee9ef0d 1453 Double_t position[3]={fVert.GetXv(),fVert.GetYv(),fVert.GetZv()};
146c29df 1454 Double_t covmatrix[6];
1455 fVert.GetCovMatrix(covmatrix);
fee9ef0d 1456 Double_t chi2=99999.;
6d8df534 1457 Int_t nTrksUsed=fVert.GetNContributors();
1458 fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);
fee9ef0d 1459 }
1460 fCurrentVertex->SetDispersion(fVert.GetDispersion());
1461
1462
1463 // set indices of used tracks and propagate track to found vertex
50ff0bcd 1464 UShort_t *indices = 0;
6d8df534 1465 Double_t d0z0[2],covd0z0[3];
1466 AliExternalTrackParam *t = 0;
fee9ef0d 1467 if(fCurrentVertex->GetNContributors()>0) {
1468 indices = new UShort_t[fCurrentVertex->GetNContributors()];
6d8df534 1469 for(Int_t jj=0; jj<(Int_t)fTrkArraySel.GetEntriesFast(); jj++) {
1470 indices[jj] = fIdSel[jj];
1471 t = (AliExternalTrackParam*)fTrkArraySel.At(jj);
1472 if(optPropagate && optUseFitter) {
fee9ef0d 1473 if(TMath::Sqrt(fCurrentVertex->GetXv()*fCurrentVertex->GetXv()+fCurrentVertex->GetYv()*fCurrentVertex->GetYv())<3.) {
6d8df534 1474 t->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
fee9ef0d 1475 if(fDebug) printf("Track %d propagated to found vertex\n",jj);
6d8df534 1476 } else {
fee9ef0d 1477 AliWarning("Found vertex outside beam pipe!");
1478 }
1479 }
50ff0bcd 1480 }
fee9ef0d 1481 fCurrentVertex->SetIndices(fCurrentVertex->GetNContributors(),indices);
50ff0bcd 1482 }
6d8df534 1483
1484 // clean up
1485 delete [] indices; indices=NULL;
1486 delete [] fIdSel; fIdSel=NULL;
1487 fTrkArraySel.Delete();
50ff0bcd 1488
fee9ef0d 1489 return fCurrentVertex;
eae1677e 1490}
50ff0bcd 1491//----------------------------------------------------------------------------
6d8df534 1492AliESDVertex* AliVertexerTracks::VertexForSelectedESDTracks(TObjArray *trkArray,
1493 Bool_t optUseFitter,
1494 Bool_t optPropagate)
fee9ef0d 1495{
50ff0bcd 1496//
6d8df534 1497// Return vertex from array of ESD tracks
50ff0bcd 1498//
5834e9a1 1499
6d8df534 1500 Int_t nTrks = (Int_t)trkArray->GetEntriesFast();
1501 UShort_t *id = new UShort_t[nTrks];
1502
1503 AliESDtrack *esdt = 0;
50ff0bcd 1504 for(Int_t i=0; i<nTrks; i++){
6d8df534 1505 esdt = (AliESDtrack*)trkArray->At(i);
1506 id[i] = (UShort_t)esdt->GetID();
50ff0bcd 1507 }
1508
6d8df534 1509 VertexForSelectedTracks(trkArray,id,optUseFitter,optPropagate);
1510
1511 delete [] id; id=NULL;
1512
1513 return fCurrentVertex;
50ff0bcd 1514}
1515//--------------------------------------------------------------------------