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