]> 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
Updating filename
[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
18eba197 561 Bool_t propagateOK = kFALSE;
50ff0bcd 562 // propagate track to vertex
43c9dae1 563 if(optImpParCut<2 || fOnlyFitter) { // optImpParCut==1 or 0
18eba197 564 propagateOK = track->PropagateToDCA(initVertex,GetFieldkG(),100.,d0z0,covd0z0);
fee9ef0d 565 } else { // optImpParCut==2
50ff0bcd 566 fCurrentVertex->GetSigmaXYZ(sigmaCurr);
fee9ef0d 567 normdistx = TMath::Abs(fCurrentVertex->GetXv()-fNominalPos[0])/TMath::Sqrt(sigmaCurr[0]*sigmaCurr[0]+fNominalCov[0]);
80d36431 568 normdisty = TMath::Abs(fCurrentVertex->GetYv()-fNominalPos[1])/TMath::Sqrt(sigmaCurr[1]*sigmaCurr[1]+fNominalCov[2]);
fee9ef0d 569 if(normdistx < 3. && normdisty < 3. &&
80d36431 570 (sigmaCurr[0]+sigmaCurr[1])<(TMath::Sqrt(fNominalCov[0])+TMath::Sqrt(fNominalCov[2]))) {
18eba197 571 propagateOK = track->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
50ff0bcd 572 } else {
18eba197 573 propagateOK = track->PropagateToDCA(initVertex,GetFieldkG(),100.,d0z0,covd0z0);
50ff0bcd 574 }
575 }
576
18eba197 577 if(!propagateOK) {
578 if(fDebug) printf(" rejected\n");
579 delete track; continue;
580 }
581
43c9dae1 582 sigmad0 = TMath::Sqrt(covd0z0[0]);
583 maxd0rphi = fNSigma*sigmad0;
fee9ef0d 584 if(optImpParCut==1) maxd0rphi *= 5.;
6b29d399 585 maxd0rphi = TMath::Min(maxd0rphi,fFiducialR);
43c9dae1 586 //sigmad0z0 = TMath::Sqrt(covd0z0[0]+covd0z0[2]); // for future improvement
587 //maxd0z0 = 10.*fNSigma*sigmad0z0;
fee9ef0d 588
f09c879d 589 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);
590
591
592 //---- track selection based on impact parameters ----//
593
594 // always reject tracks outside fiducial volume
6b29d399 595 if(TMath::Abs(d0z0[0])>fFiducialR || TMath::Abs(d0z0[1])>fFiducialZ) {
f09c879d 596 if(fDebug) printf(" rejected\n");
597 delete track; continue;
598 }
599
600 // during iterations 1 and 2 , reject tracks with d0rphi > maxd0rphi
601 if(optImpParCut>0 && TMath::Abs(d0z0[0])>maxd0rphi) {
602 if(fDebug) printf(" rejected\n");
603 delete track; continue;
604 }
fee9ef0d 605
43c9dae1 606 // if fConstraint=kFALSE, during iterations 1 and 2 ||
607 // if fConstraint=kTRUE, during iteration 2,
f09c879d 608 // select tracks with d0oplusz0 < fMaxd0z0
43c9dae1 609 if((!fConstraint && optImpParCut>0 && fVert.GetNContributors()>0) ||
610 ( fConstraint && optImpParCut==2)) {
611 if(nTrksOrig>=3 &&
f09c879d 612 TMath::Sqrt(d0z0[0]*d0z0[0]+d0z0[1]*d0z0[1])>fMaxd0z0) {
fee9ef0d 613 if(fDebug) printf(" rejected\n");
614 delete track; continue;
615 }
616 }
43c9dae1 617
f09c879d 618 // track passed all selections
6d8df534 619 fTrkArraySel.AddLast(track);
620 fIdSel[nTrksSel] = idOrig[i];
621 nTrksSel++;
f09c879d 622 } // end loop on tracks
50ff0bcd 623
624 delete initVertex;
625
6d8df534 626 return nTrksSel;
50ff0bcd 627}
f09c879d 628//----------------------------------------------------------------------------
629Bool_t AliVertexerTracks::PropagateTrackTo(AliExternalTrackParam *track,
630 Double_t xToGo) {
631 //----------------------------------------------------------------
632 // COPIED from AliTracker
633 //
634 // Propagates the track to the plane X=xk (cm).
635 //
636 // Origin: Marian Ivanov, Marian.Ivanov@cern.ch
637 //----------------------------------------------------------------
638
639 const Double_t kEpsilon = 0.00001;
640 Double_t xpos = track->GetX();
641 Double_t dir = (xpos<xToGo) ? 1. : -1.;
642 Double_t maxStep = 5;
643 Double_t maxSnp = 0.8;
644 //
645 while ( (xToGo-xpos)*dir > kEpsilon){
646 Double_t step = dir*TMath::Min(TMath::Abs(xToGo-xpos), maxStep);
647 Double_t x = xpos+step;
648 Double_t xyz0[3],xyz1[3];
649 track->GetXYZ(xyz0); //starting global position
650
651 if(!track->GetXYZAt(x,GetFieldkG(),xyz1)) return kFALSE; // no prolongation
652 xyz1[2]+=kEpsilon; // waiting for bug correction in geo
653
654 if(TMath::Abs(track->GetSnpAt(x,GetFieldkG())) >= maxSnp) return kFALSE;
655 if(!track->PropagateTo(x,GetFieldkG())) return kFALSE;
656
657 if(TMath::Abs(track->GetSnp()) >= maxSnp) return kFALSE;
658 track->GetXYZ(xyz0); // global position
659 Double_t alphan = TMath::ATan2(xyz0[1], xyz0[0]);
660 //
661 Double_t ca=TMath::Cos(alphan-track->GetAlpha()),
662 sa=TMath::Sin(alphan-track->GetAlpha());
663 Double_t sf=track->GetSnp(), cf=TMath::Sqrt(1.- sf*sf);
664 Double_t sinNew = sf*ca - cf*sa;
665 if(TMath::Abs(sinNew) >= maxSnp) return kFALSE;
666 if(!track->Rotate(alphan)) return kFALSE;
667
668 xpos = track->GetX();
669 }
670 return kTRUE;
671}
50ff0bcd 672//---------------------------------------------------------------------------
fee9ef0d 673AliESDVertex* AliVertexerTracks::RemoveTracksFromVertex(AliESDVertex *inVtx,
6d8df534 674 TObjArray *trkArray,
675 UShort_t *id,
fee9ef0d 676 Float_t *diamondxy)
677{
50ff0bcd 678//
fee9ef0d 679// Removes tracks in trksTree from fit of inVtx
50ff0bcd 680//
fee9ef0d 681
146c29df 682 if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraint")) {
683 printf("ERROR: primary vertex has no constraint: cannot remove tracks\n");
684 return 0x0;
685 }
fee9ef0d 686 if(!strstr(inVtx->GetTitle(),"VertexerTracksWithConstraintOnlyFitter"))
687 printf("WARNING: result of tracks' removal will be only approximately correct\n");
688
689 TMatrixD rv(3,1);
690 rv(0,0) = inVtx->GetXv();
691 rv(1,0) = inVtx->GetYv();
692 rv(2,0) = inVtx->GetZv();
693 TMatrixD vV(3,3);
694 Double_t cov[6];
695 inVtx->GetCovMatrix(cov);
696 vV(0,0) = cov[0];
697 vV(0,1) = cov[1]; vV(1,0) = cov[1];
698 vV(1,1) = cov[2];
699 vV(0,2) = cov[3]; vV(2,0) = cov[3];
700 vV(1,2) = cov[4]; vV(2,1) = cov[4];
701 vV(2,2) = cov[5];
702
703 TMatrixD sumWi(TMatrixD::kInverted,vV);
704 TMatrixD sumWiri(sumWi,TMatrixD::kMult,rv);
705
dad616ce 706 Int_t nUsedTrks = inVtx->GetNIndices();
fee9ef0d 707 Double_t chi2 = inVtx->GetChi2();
708
6d8df534 709 AliExternalTrackParam *track = 0;
710 Int_t ntrks = (Int_t)trkArray->GetEntriesFast();
fee9ef0d 711 for(Int_t i=0;i<ntrks;i++) {
6d8df534 712 track = (AliExternalTrackParam*)trkArray->At(i);
713 if(!inVtx->UsesTrack(id[i])) {
714 printf("track %d was not used in vertex fit\n",id[i]);
fee9ef0d 715 continue;
716 }
717 Double_t alpha = track->GetAlpha();
718 Double_t xl = diamondxy[0]*TMath::Cos(alpha)+diamondxy[1]*TMath::Sin(alpha);
6d8df534 719 track->PropagateTo(xl,GetFieldkG());
fee9ef0d 720 // vector of track global coordinates
721 TMatrixD ri(3,1);
722 // covariance matrix of ri
723 TMatrixD wWi(3,3);
724
725 // get space point from track
726 if(!TrackToPoint(track,ri,wWi)) continue;
727
728 TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
729
730 sumWi -= wWi;
731 sumWiri -= wWiri;
732
733 // track chi2
734 TMatrixD deltar = rv; deltar -= ri;
735 TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
736 Double_t chi2i = deltar(0,0)*wWideltar(0,0)+
737 deltar(1,0)*wWideltar(1,0)+
738 deltar(2,0)*wWideltar(2,0);
739 // remove from total chi2
740 chi2 -= chi2i;
741
742 nUsedTrks--;
743 if(nUsedTrks<2) {
744 printf("Trying to remove too many tracks!\n");
745 return 0x0;
746 }
747 }
748
749 TMatrixD rvnew(3,1);
750 TMatrixD vVnew(3,3);
751
752 // new inverted of weights matrix
753 TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
754 vVnew = invsumWi;
755 // new position of primary vertex
756 rvnew.Mult(vVnew,sumWiri);
757
758 Double_t position[3];
759 position[0] = rvnew(0,0);
760 position[1] = rvnew(1,0);
761 position[2] = rvnew(2,0);
762 cov[0] = vVnew(0,0);
763 cov[1] = vVnew(0,1);
764 cov[2] = vVnew(1,1);
765 cov[3] = vVnew(0,2);
766 cov[4] = vVnew(1,2);
767 cov[5] = vVnew(2,2);
768
769 // store data in the vertex object
770 AliESDVertex *outVtx = new AliESDVertex(position,cov,chi2,nUsedTrks);
771 outVtx->SetTitle(inVtx->GetTitle());
fee9ef0d 772 UShort_t *inindices = inVtx->GetIndices();
dad616ce 773 Int_t nIndices = nUsedTrks;
4d023a8f 774 UShort_t *outindices = new UShort_t[nIndices];
fee9ef0d 775 Int_t j=0;
fee9ef0d 776 for(Int_t k=0; k<inVtx->GetNIndices(); k++) {
6d8df534 777 Bool_t copyindex=kTRUE;
fee9ef0d 778 for(Int_t l=0; l<ntrks; l++) {
6d8df534 779 if(inindices[k]==id[l]) copyindex=kFALSE;
fee9ef0d 780 }
781 if(copyindex) {
782 outindices[j] = inindices[k]; j++;
783 }
784 }
4d023a8f 785 outVtx->SetIndices(nIndices,outindices);
fee9ef0d 786 delete [] outindices;
787
788 if(fDebug) {
789 printf("Vertex before removing tracks:\n");
790 inVtx->PrintStatus();
791 inVtx->PrintIndices();
792 printf("Vertex after removing tracks:\n");
793 outVtx->PrintStatus();
794 outVtx->PrintIndices();
795 }
796
797 return outVtx;
798}
799//---------------------------------------------------------------------------
a00021a7 800void AliVertexerTracks::SetCuts(Double_t *cuts)
801{
802//
803// Cut values
804//
805 SetDCAcut(cuts[0]);
806 SetDCAcutIter0(cuts[1]);
807 SetMaxd0z0(cuts[2]);
5b78b791 808 if(fMode==0 && cuts[3]<0) SetITSrefitNotRequired();
809 SetMinClusters((Int_t)(TMath::Abs(cuts[3])));
a00021a7 810 SetMinTracks((Int_t)(cuts[4]));
811 SetNSigmad0(cuts[5]);
812 SetMinDetFitter(cuts[6]);
813 SetMaxTgl(cuts[7]);
814 SetFiducialRZ(cuts[8],cuts[9]);
815
816 return;
817}
818//---------------------------------------------------------------------------
f09c879d 819void AliVertexerTracks::SetITSMode(Double_t dcacut,
820 Double_t dcacutIter0,
821 Double_t maxd0z0,
6b29d399 822 Int_t minCls,
f09c879d 823 Int_t mintrks,
824 Double_t nsigma,
825 Double_t mindetfitter,
6b29d399 826 Double_t maxtgl,
827 Double_t fidR,
828 Double_t fidZ)
f09c879d 829{
830//
831// Cut values for ITS mode
832//
833 fMode = 0;
5b78b791 834 if(minCls>0) {
835 SetITSrefitRequired();
836 } else {
837 SetITSrefitNotRequired();
838 }
f09c879d 839 SetDCAcut(dcacut);
840 SetDCAcutIter0(dcacutIter0);
841 SetMaxd0z0(maxd0z0);
5b78b791 842 SetMinClusters(TMath::Abs(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//---------------------------------------------------------------------------
852void AliVertexerTracks::SetTPCMode(Double_t dcacut,
853 Double_t dcacutIter0,
854 Double_t maxd0z0,
6b29d399 855 Int_t minCls,
f09c879d 856 Int_t mintrks,
857 Double_t nsigma,
858 Double_t mindetfitter,
6b29d399 859 Double_t maxtgl,
860 Double_t fidR,
861 Double_t fidZ)
f09c879d 862{
863//
864// Cut values for TPC mode
865//
866 fMode = 1;
867 SetITSrefitNotRequired();
868 SetDCAcut(dcacut);
869 SetDCAcutIter0(dcacutIter0);
870 SetMaxd0z0(maxd0z0);
6b29d399 871 SetMinClusters(minCls);
f09c879d 872 SetMinTracks(mintrks);
873 SetNSigmad0(nsigma);
874 SetMinDetFitter(mindetfitter);
875 SetMaxTgl(maxtgl);
6b29d399 876 SetFiducialRZ(fidR,fidZ);
f09c879d 877
878 return;
879}
880//---------------------------------------------------------------------------
fee9ef0d 881void AliVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped)
882{
883//
884// Mark the tracks not to be used in the vertex reconstruction.
885// Tracks are identified by AliESDtrack::GetID()
886//
887 fNTrksToSkip = n; fTrksToSkip = new Int_t[n];
50ff0bcd 888 for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i];
889 return;
890}
891//---------------------------------------------------------------------------
fee9ef0d 892void AliVertexerTracks::SetVtxStart(AliESDVertex *vtx)
893{
50ff0bcd 894//
895// Set initial vertex knowledge
896//
897 vtx->GetXYZ(fNominalPos);
898 vtx->GetCovMatrix(fNominalCov);
fee9ef0d 899 SetConstraintOn();
50ff0bcd 900 return;
901}
2d57349e 902//---------------------------------------------------------------------------
fee9ef0d 903void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights)
dc3719f3 904{
6d8df534 905 AliExternalTrackParam *track1;
6d8df534 906 const Int_t knacc = (Int_t)fTrkArraySel.GetEntriesFast();
f54bad66 907 static TClonesArray linarray("AliStrLine",knacc);
dc3719f3 908 for(Int_t i=0; i<knacc; i++){
6d8df534 909 track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
dc3719f3 910 Double_t alpha=track1->GetAlpha();
911 Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
912 Double_t pos[3],dir[3],sigmasq[3];
f09c879d 913 track1->GetXYZAt(mindist,GetFieldkG(),pos);
914 track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
dc3719f3 915 sigmasq[0]=TMath::Sin(alpha)*TMath::Sin(alpha)*track1->GetSigmaY2();
916 sigmasq[1]=TMath::Cos(alpha)*TMath::Cos(alpha)*track1->GetSigmaY2();
917 sigmasq[2]=track1->GetSigmaZ2();
146c29df 918 TMatrixD ri(3,1);
919 TMatrixD wWi(3,3);
920 if(!TrackToPoint(track1,ri,wWi)) continue;
921 Double_t wmat[9];
922 Int_t iel=0;
923 for(Int_t ia=0;ia<3;ia++){
924 for(Int_t ib=0;ib<3;ib++){
925 wmat[iel]=wWi(ia,ib);
926 iel++;
927 }
928 }
f54bad66 929 new(linarray[i]) AliStrLine(pos,sigmasq,wmat,dir);
dc3719f3 930 }
f54bad66 931 fVert=TrackletVertexFinder(&linarray,optUseWeights);
932 linarray.Clear("C");
dc3719f3 933}
934//---------------------------------------------------------------------------
146c29df 935AliESDVertex AliVertexerTracks::TrackletVertexFinder(TClonesArray *lines, Int_t optUseWeights)
fee9ef0d 936{
2d57349e 937 // Calculate the point at minimum distance to prepared tracks
938
dc3719f3 939 const Int_t knacc = (Int_t)lines->GetEntriesFast();
940 Double_t initPos[3]={0.,0.,0.};
2d57349e 941
2d57349e 942 Double_t (*vectP0)[3]=new Double_t [knacc][3];
943 Double_t (*vectP1)[3]=new Double_t [knacc][3];
944
945 Double_t sum[3][3];
946 Double_t dsum[3]={0,0,0};
146c29df 947 TMatrixD sumWi(3,3);
948 for(Int_t i=0;i<3;i++){
949 for(Int_t j=0;j<3;j++){
950 sum[i][j]=0;
951 sumWi(i,j)=0.;
952 }
953 }
954
2d57349e 955 for(Int_t i=0; i<knacc; i++){
f54bad66 956 AliStrLine* line1 = (AliStrLine*)lines->At(i);
dc3719f3 957 Double_t p0[3],cd[3],sigmasq[3];
146c29df 958 Double_t wmat[9];
2d57349e 959 line1->GetP0(p0);
960 line1->GetCd(cd);
dc3719f3 961 line1->GetSigma2P0(sigmasq);
146c29df 962 line1->GetWMatrix(wmat);
963 TMatrixD wWi(3,3);
964 Int_t iel=0;
965 for(Int_t ia=0;ia<3;ia++){
966 for(Int_t ib=0;ib<3;ib++){
967 wWi(ia,ib)=wmat[iel];
968 iel++;
969 }
970 }
971
972 sumWi+=wWi;
973
2d57349e 974 Double_t p1[3]={p0[0]+cd[0],p0[1]+cd[1],p0[2]+cd[2]};
975 vectP0[i][0]=p0[0];
976 vectP0[i][1]=p0[1];
977 vectP0[i][2]=p0[2];
978 vectP1[i][0]=p1[0];
979 vectP1[i][1]=p1[1];
980 vectP1[i][2]=p1[2];
981
982 Double_t matr[3][3];
983 Double_t dknow[3];
984 if(optUseWeights==0)GetStrLinDerivMatrix(p0,p1,matr,dknow);
dc3719f3 985 else GetStrLinDerivMatrix(p0,p1,sigmasq,matr,dknow);
986
2d57349e 987
988 for(Int_t iii=0;iii<3;iii++){
989 dsum[iii]+=dknow[iii];
990 for(Int_t lj=0;lj<3;lj++) sum[iii][lj]+=matr[iii][lj];
991 }
2d57349e 992 }
146c29df 993
994 TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
995 Double_t covmatrix[6];
996 covmatrix[0] = invsumWi(0,0);
997 covmatrix[1] = invsumWi(0,1);
998 covmatrix[2] = invsumWi(1,1);
999 covmatrix[3] = invsumWi(0,2);
1000 covmatrix[4] = invsumWi(1,2);
1001 covmatrix[5] = invsumWi(2,2);
1002
2d57349e 1003 Double_t vett[3][3];
1004 Double_t det=GetDeterminant3X3(sum);
dc3719f3 1005 Double_t sigma=0;
2d57349e 1006
fee9ef0d 1007 if(det!=0){
1008 for(Int_t zz=0;zz<3;zz++){
1009 for(Int_t ww=0;ww<3;ww++){
1010 for(Int_t kk=0;kk<3;kk++) vett[ww][kk]=sum[ww][kk];
1011 }
1012 for(Int_t kk=0;kk<3;kk++) vett[kk][zz]=dsum[kk];
1013 initPos[zz]=GetDeterminant3X3(vett)/det;
1014 }
1015
1016
1017 for(Int_t i=0; i<knacc; i++){
1018 Double_t p0[3]={0,0,0},p1[3]={0,0,0};
1019 for(Int_t ii=0;ii<3;ii++){
1020 p0[ii]=vectP0[i][ii];
1021 p1[ii]=vectP1[i][ii];
1022 }
1023 sigma+=GetStrLinMinDist(p0,p1,initPos);
1024 }
1025
f09c879d 1026 if(sigma>0.) {sigma=TMath::Sqrt(sigma);}else{sigma=999;}
fee9ef0d 1027 }else{
50ff0bcd 1028 sigma=999;
1029 }
146c29df 1030 AliESDVertex theVert(initPos,covmatrix,99999.,knacc);
1031 theVert.SetDispersion(sigma);
50ff0bcd 1032 delete vectP0;
1033 delete vectP1;
dc3719f3 1034 return theVert;
50ff0bcd 1035}
1036//---------------------------------------------------------------------------
6d8df534 1037Bool_t AliVertexerTracks::TrackToPoint(AliExternalTrackParam *t,
0bf2e2b4 1038 TMatrixD &ri,TMatrixD &wWi,
1039 Bool_t uUi3by3) const
fee9ef0d 1040{
1041//
6d8df534 1042// Extract from the AliExternalTrackParam the global coordinates ri and covariance matrix
fee9ef0d 1043// wWi of the space point that it represents (to be used in VertexFitter())
1044//
1045
1046
1047 Double_t rotAngle = t->GetAlpha();
1048 if(rotAngle<0.) rotAngle += 2.*TMath::Pi();
1049 Double_t cosRot = TMath::Cos(rotAngle);
1050 Double_t sinRot = TMath::Sin(rotAngle);
1051
1052 ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot;
1053 ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot;
1054 ri(2,0) = t->GetZ();
1055
0bf2e2b4 1056
1057
1058 if(!uUi3by3) {
1059 // matrix to go from global (x,y,z) to local (y,z);
1060 TMatrixD qQi(2,3);
1061 qQi(0,0) = -sinRot;
1062 qQi(0,1) = cosRot;
1063 qQi(0,2) = 0.;
1064 qQi(1,0) = 0.;
1065 qQi(1,1) = 0.;
1066 qQi(1,2) = 1.;
1067
1068 // covariance matrix of local (y,z) - inverted
0bf2e2b4 1069 TMatrixD uUi(2,2);
6d8df534 1070 uUi(0,0) = t->GetSigmaY2();
1071 uUi(0,1) = t->GetSigmaZY();
1072 uUi(1,0) = t->GetSigmaZY();
1073 uUi(1,1) = t->GetSigmaZ2();
0bf2e2b4 1074 //printf(" Ui :\n");
1075 //printf(" %f %f\n",uUi(0,0),uUi(0,1));
1076 //printf(" %f %f\n",uUi(1,0),uUi(1,1));
1077
1078 if(uUi.Determinant() <= 0.) return kFALSE;
1079 TMatrixD uUiInv(TMatrixD::kInverted,uUi);
1080
1081 // weights matrix: wWi = qQiT * uUiInv * qQi
1082 TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
1083 TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
1084 wWi = m;
1085 } else {
1086 if(fVert.GetNContributors()<1) AliFatal("Vertex from finder is empty");
1087 // matrix to go from global (x,y,z) to local (x,y,z);
1088 TMatrixD qQi(3,3);
1089 qQi(0,0) = cosRot;
1090 qQi(0,1) = sinRot;
1091 qQi(0,2) = 0.;
1092 qQi(1,0) = -sinRot;
1093 qQi(1,1) = cosRot;
1094 qQi(1,2) = 0.;
1095 qQi(2,0) = 0.;
1096 qQi(2,1) = 0.;
1097 qQi(2,2) = 1.;
1098
1099 // covariance of fVert along the track
1100 Double_t p[3],pt,ptot;
1101 t->GetPxPyPz(p);
1102 pt = TMath::Sqrt(p[0]*p[0]+p[1]*p[1]);
1103 ptot = TMath::Sqrt(pt*pt+p[2]*p[2]);
3b768a1d 1104 Double_t cphi = p[0]/pt; //cos(phi)=px/pt
1105 Double_t sphi = p[1]/pt; //sin(phi)=py/pt
1106 Double_t clambda = pt/ptot; //cos(lambda)=pt/ptot
1107 Double_t slambda = p[2]/ptot; //sin(lambda)=pz/ptot
0bf2e2b4 1108 Double_t covfVert[6];
1109 fVert.GetCovMatrix(covfVert);
1110 Double_t covfVertalongt =
3b768a1d 1111 covfVert[0]*cphi*cphi*clambda*clambda
1112 +covfVert[1]*2.*cphi*sphi*clambda*clambda
1113 +covfVert[3]*2.*cphi*clambda*slambda
1114 +covfVert[2]*sphi*sphi*clambda*clambda
1115 +covfVert[4]*2.*sphi*clambda*slambda
1116 +covfVert[5]*slambda*slambda;
0bf2e2b4 1117 // covariance matrix of local (x,y,z) - inverted
1118 TMatrixD uUi(3,3);
3b768a1d 1119 uUi(0,0) = covfVertalongt * fnSigmaForUi00 * fnSigmaForUi00;
0bf2e2b4 1120 if(fDebug) printf("=====> sqrtUi00 cm %f\n",TMath::Sqrt(uUi(0,0)));
1121 uUi(0,1) = 0.;
1122 uUi(0,2) = 0.;
1123 uUi(1,0) = 0.;
6d8df534 1124 uUi(1,1) = t->GetSigmaY2();
1125 uUi(1,2) = t->GetSigmaZY();
0bf2e2b4 1126 uUi(2,0) = 0.;
6d8df534 1127 uUi(2,1) = t->GetSigmaZY();
1128 uUi(2,2) = t->GetSigmaZ2();
0bf2e2b4 1129 //printf(" Ui :\n");
1130 //printf(" %f %f\n",uUi(0,0),uUi(0,1));
1131 //printf(" %f %f\n",uUi(1,0),uUi(1,1));
fee9ef0d 1132
0bf2e2b4 1133 if(uUi.Determinant() <= 0.) return kFALSE;
1134 TMatrixD uUiInv(TMatrixD::kInverted,uUi);
1135
1136 // weights matrix: wWi = qQiT * uUiInv * qQi
1137 TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
1138 TMatrixD m(qQi,TMatrixD::kTransposeMult,uUiInvQi);
1139 wWi = m;
1140 }
1141
fee9ef0d 1142
1143 return kTRUE;
1144}
1145//---------------------------------------------------------------------------
6d8df534 1146void AliVertexerTracks::TooFewTracks()
fee9ef0d 1147{
50ff0bcd 1148//
6d8df534 1149// When the number of tracks is < fMinTracks,
1150// deal with vertices not found and prepare to exit
50ff0bcd 1151//
6d8df534 1152 if(fDebug) printf("TooFewTracks\n");
2d57349e 1153
50ff0bcd 1154 Double_t pos[3],err[3];
50ff0bcd 1155 pos[0] = fNominalPos[0];
1156 err[0] = TMath::Sqrt(fNominalCov[0]);
1157 pos[1] = fNominalPos[1];
1158 err[1] = TMath::Sqrt(fNominalCov[2]);
6d8df534 1159 pos[2] = fNominalPos[2];
1160 err[2] = TMath::Sqrt(fNominalCov[5]);
1161 Int_t ncontr = (err[0]>1. ? -1 : -3);
5b78b791 1162 if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
50ff0bcd 1163 fCurrentVertex = new AliESDVertex(pos,err);
1164 fCurrentVertex->SetNContributors(ncontr);
2d57349e 1165
fee9ef0d 1166 if(fConstraint) {
1167 fCurrentVertex->SetTitle("VertexerTracksWithConstraint");
1168 } else {
50ff0bcd 1169 fCurrentVertex->SetTitle("VertexerTracksNoConstraint");
fee9ef0d 1170 }
2d57349e 1171
6d8df534 1172 if(!fTrkArraySel.IsEmpty()) fTrkArraySel.Delete();
1173 if(fIdSel) {delete [] fIdSel; fIdSel=NULL;}
1174 if(fTrksToSkip) {delete [] fTrksToSkip; fTrksToSkip=NULL;}
1175
50ff0bcd 1176 return;
1177}
1178//---------------------------------------------------------------------------
fee9ef0d 1179void AliVertexerTracks::VertexFinder(Int_t optUseWeights)
1180{
2d57349e 1181
50ff0bcd 1182 // Get estimate of vertex position in (x,y) from tracks DCA
1183
1184 Double_t initPos[3];
1185 initPos[2] = 0.;
1186 for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
6d8df534 1187 Int_t nacc = (Int_t)fTrkArraySel.GetEntriesFast();
50ff0bcd 1188 Double_t aver[3]={0.,0.,0.};
1189 Double_t aversq[3]={0.,0.,0.};
1190 Double_t sigmasq[3]={0.,0.,0.};
1191 Double_t sigma=0;
1192 Int_t ncombi = 0;
6d8df534 1193 AliExternalTrackParam *track1;
1194 AliExternalTrackParam *track2;
50ff0bcd 1195 Double_t pos[3],dir[3];
1196 Double_t alpha,mindist;
2d57349e 1197
50ff0bcd 1198 for(Int_t i=0; i<nacc; i++){
6d8df534 1199 track1 = (AliExternalTrackParam*)fTrkArraySel.At(i);
50ff0bcd 1200 alpha=track1->GetAlpha();
1201 mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
f09c879d 1202 track1->GetXYZAt(mindist,GetFieldkG(),pos);
1203 track1->GetPxPyPzAt(mindist,GetFieldkG(),dir);
50ff0bcd 1204 AliStrLine *line1 = new AliStrLine(pos,dir);
2d57349e 1205
50ff0bcd 1206 // AliStrLine *line1 = new AliStrLine();
f09c879d 1207 // track1->ApproximateHelixWithLine(mindist,GetFieldkG(),line1);
50ff0bcd 1208
1209 for(Int_t j=i+1; j<nacc; j++){
6d8df534 1210 track2 = (AliExternalTrackParam*)fTrkArraySel.At(j);
50ff0bcd 1211 alpha=track2->GetAlpha();
1212 mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
f09c879d 1213 track2->GetXYZAt(mindist,GetFieldkG(),pos);
1214 track2->GetPxPyPzAt(mindist,GetFieldkG(),dir);
50ff0bcd 1215 AliStrLine *line2 = new AliStrLine(pos,dir);
1216 // AliStrLine *line2 = new AliStrLine();
f09c879d 1217 // track2->ApproximateHelixWithLine(mindist,GetFieldkG(),line2);
50ff0bcd 1218 Double_t distCA=line2->GetDCA(line1);
fee9ef0d 1219 //printf("%d %d %f\n",i,j,distCA);
1220 if(fDCAcut<=0 || (fDCAcut>0&&distCA<fDCAcut)){
50ff0bcd 1221 Double_t pnt1[3],pnt2[3],crosspoint[3];
2d57349e 1222
50ff0bcd 1223 if(optUseWeights<=0){
1224 Int_t retcode = line2->Cross(line1,crosspoint);
1225 if(retcode>=0){
1226 ncombi++;
1227 for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
1228 for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
1229 }
1230 }
1231 if(optUseWeights>0){
1232 Int_t retcode = line1->CrossPoints(line2,pnt1,pnt2);
1233 if(retcode>=0){
5c03260f 1234 Double_t cs, sn;
50ff0bcd 1235 alpha=track1->GetAlpha();
1236 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
1237 Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
1238 alpha=track2->GetAlpha();
1239 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
1240 Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
1241 Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
1242 Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
1243 Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
1244 Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
1245 crosspoint[0]=wx1*pnt1[0] + wx2*pnt2[0];
1246 crosspoint[1]=wy1*pnt1[1] + wy2*pnt2[1];
1247 crosspoint[2]=wz1*pnt1[2] + wz2*pnt2[2];
1248
1249 ncombi++;
1250 for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
1251 for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
1252 }
1253 }
1254 }
1255 delete line2;
1256 }
1257 delete line1;
71d84967 1258 }
50ff0bcd 1259 if(ncombi>0){
1260 for(Int_t jj=0;jj<3;jj++){
1261 initPos[jj] = aver[jj]/ncombi;
fee9ef0d 1262 //printf("%f\n",initPos[jj]);
50ff0bcd 1263 aversq[jj]/=ncombi;
1264 sigmasq[jj]=aversq[jj]-initPos[jj]*initPos[jj];
1265 sigma+=sigmasq[jj];
1266 }
1267 sigma=TMath::Sqrt(TMath::Abs(sigma));
1268 }
1269 else {
1270 Warning("VertexFinder","Finder did not succed");
1271 sigma=999;
1272 }
1273 fVert.SetXYZ(initPos);
1274 fVert.SetDispersion(sigma);
1275 fVert.SetNContributors(ncombi);
2d57349e 1276}
ec1be5d5 1277//---------------------------------------------------------------------------
6d8df534 1278void AliVertexerTracks::VertexFitter()
fee9ef0d 1279{
ec1be5d5 1280//
1281// The optimal estimate of the vertex position is given by a "weighted
0bf2e2b4 1282// average of tracks positions".
fee9ef0d 1283// Original method: V. Karimaki, CMS Note 97/0051
ec1be5d5 1284//
6d8df534 1285 Bool_t useConstraint = fConstraint;
ec1be5d5 1286 Double_t initPos[3];
0bf2e2b4 1287 if(!fOnlyFitter) {
1288 fVert.GetXYZ(initPos);
1289 } else {
1290 initPos[0]=fNominalPos[0];
1291 initPos[1]=fNominalPos[1];
1292 initPos[2]=fNominalPos[2];
1293 }
1294
6d8df534 1295 Int_t nTrksSel = (Int_t)fTrkArraySel.GetEntries();
1296 if(nTrksSel==1) useConstraint=kTRUE;
ec1be5d5 1297 if(fDebug) {
0bf2e2b4 1298 printf("--- VertexFitter(): start\n");
6d8df534 1299 printf(" Number of tracks in array: %d\n",nTrksSel);
ec1be5d5 1300 printf(" Minimum # tracks required in fit: %d\n",fMinTracks);
3b768a1d 1301 printf(" Vertex position after finder: %f,%f,%f\n",initPos[0],initPos[1],initPos[2]);
fee9ef0d 1302 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 1303 }
1304
0bf2e2b4 1305 // special treatment for few-tracks fits (e.g. secondary vertices)
6d8df534 1306 Bool_t uUi3by3 = kFALSE; if(nTrksSel<5 && !useConstraint) uUi3by3 = kTRUE;
0bf2e2b4 1307
ec1be5d5 1308 Int_t i,j,k,step=0;
1309 TMatrixD rv(3,1);
1310 TMatrixD vV(3,3);
1311 rv(0,0) = initPos[0];
1312 rv(1,0) = initPos[1];
1313 rv(2,0) = 0.;
1314 Double_t xlStart,alpha;
6d8df534 1315 Int_t nTrksUsed;
fee9ef0d 1316 Double_t chi2,chi2i,chi2b;
6d8df534 1317 AliExternalTrackParam *t = 0;
ec1be5d5 1318 Int_t failed = 0;
1319
50ff0bcd 1320 // initial vertex covariance matrix
1321 TMatrixD vVb(3,3);
1322 vVb(0,0) = fNominalCov[0];
1323 vVb(0,1) = fNominalCov[1];
1324 vVb(0,2) = 0.;
1325 vVb(1,0) = fNominalCov[1];
1326 vVb(1,1) = fNominalCov[2];
1327 vVb(1,2) = 0.;
1328 vVb(2,0) = 0.;
1329 vVb(2,1) = 0.;
1330 vVb(2,2) = fNominalCov[5];
1331 TMatrixD vVbInv(TMatrixD::kInverted,vVb);
1332 TMatrixD rb(3,1);
1333 rb(0,0) = fNominalPos[0];
1334 rb(1,0) = fNominalPos[1];
1335 rb(2,0) = fNominalPos[2];
1336 TMatrixD vVbInvrb(vVbInv,TMatrixD::kMult,rb);
1337
1338
eae1677e 1339 // 2 steps:
1340 // 1st - estimate of vtx using all tracks
1341 // 2nd - estimate of global chi2
1342 for(step=0; step<2; step++) {
ec1be5d5 1343 if(fDebug) printf(" step = %d\n",step);
1344 chi2 = 0.;
6d8df534 1345 nTrksUsed = 0;
ec1be5d5 1346
0bf2e2b4 1347 if(step==1) { initPos[0]=rv(0,0); initPos[0]=rv(1,0); }
1348
ec1be5d5 1349 TMatrixD sumWiri(3,1);
1350 TMatrixD sumWi(3,3);
1351 for(i=0; i<3; i++) {
1352 sumWiri(i,0) = 0.;
1353 for(j=0; j<3; j++) sumWi(j,i) = 0.;
1354 }
1355
fee9ef0d 1356 // mean vertex constraint
1357 if(useConstraint) {
50ff0bcd 1358 for(i=0;i<3;i++) {
1359 sumWiri(i,0) += vVbInvrb(i,0);
1360 for(k=0;k<3;k++) sumWi(i,k) += vVbInv(i,k);
07680cae 1361 }
fee9ef0d 1362 // chi2
1363 TMatrixD deltar = rv; deltar -= rb;
1364 TMatrixD vVbInvdeltar(vVbInv,TMatrixD::kMult,deltar);
1365 chi2b = deltar(0,0)*vVbInvdeltar(0,0)+
1366 deltar(1,0)*vVbInvdeltar(1,0)+
1367 deltar(2,0)*vVbInvdeltar(2,0);
1368 chi2 += chi2b;
07680cae 1369 }
1370
ec1be5d5 1371 // loop on tracks
6d8df534 1372 for(k=0; k<nTrksSel; k++) {
ec1be5d5 1373 // get track from track array
6d8df534 1374 t = (AliExternalTrackParam*)fTrkArraySel.At(k);
ec1be5d5 1375 alpha = t->GetAlpha();
1376 xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha);
fee9ef0d 1377 // to vtxSeed (from finder)
6d8df534 1378 t->PropagateTo(xlStart,GetFieldkG());
fee9ef0d 1379
ec1be5d5 1380 // vector of track global coordinates
1381 TMatrixD ri(3,1);
fee9ef0d 1382 // covariance matrix of ri
1383 TMatrixD wWi(3,3);
1384
1385 // get space point from track
0bf2e2b4 1386 if(!TrackToPoint(t,ri,wWi,uUi3by3)) continue;
fee9ef0d 1387 TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
ec1be5d5 1388
1389 // track chi2
1390 TMatrixD deltar = rv; deltar -= ri;
1391 TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
1392 chi2i = deltar(0,0)*wWideltar(0,0)+
1393 deltar(1,0)*wWideltar(1,0)+
1394 deltar(2,0)*wWideltar(2,0);
1395
ec1be5d5 1396 // add to total chi2
1397 chi2 += chi2i;
1398
ec1be5d5 1399 sumWiri += wWiri;
1400 sumWi += wWi;
1401
6d8df534 1402 nTrksUsed++;
ec1be5d5 1403 } // end loop on tracks
1404
6d8df534 1405 if(nTrksUsed < fMinTracks) {
ec1be5d5 1406 failed=1;
1407 continue;
1408 }
1409
1410 Double_t determinant = sumWi.Determinant();
f09c879d 1411 if(determinant < fMinDetFitter) {
1412 printf("det(V) = %f (<%f)\n",determinant,fMinDetFitter);
ec1be5d5 1413 failed=1;
1414 continue;
1415 }
1416
0bf2e2b4 1417 if(step==0) {
1418 // inverted of weights matrix
1419 TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
1420 vV = invsumWi;
1421 // position of primary vertex
1422 rv.Mult(vV,sumWiri);
1423 }
eae1677e 1424 } // end loop on the 2 steps
ec1be5d5 1425
ec1be5d5 1426 if(failed) {
6d8df534 1427 TooFewTracks();
ec1be5d5 1428 return;
1429 }
1430
1431 Double_t position[3];
1432 position[0] = rv(0,0);
1433 position[1] = rv(1,0);
1434 position[2] = rv(2,0);
1435 Double_t covmatrix[6];
1436 covmatrix[0] = vV(0,0);
1437 covmatrix[1] = vV(0,1);
1438 covmatrix[2] = vV(1,1);
1439 covmatrix[3] = vV(0,2);
1440 covmatrix[4] = vV(1,2);
1441 covmatrix[5] = vV(2,2);
1442
dc3719f3 1443 // for correct chi2/ndf, count constraint as additional "track"
6d8df534 1444 if(fConstraint) nTrksUsed++;
ec1be5d5 1445 // store data in the vertex object
5b78b791 1446 if(fCurrentVertex) { delete fCurrentVertex; fCurrentVertex=0; }
6d8df534 1447 fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);
ec1be5d5 1448
1449 if(fDebug) {
0bf2e2b4 1450 printf(" Vertex after fit:\n");
ec1be5d5 1451 fCurrentVertex->PrintStatus();
0bf2e2b4 1452 printf("--- VertexFitter(): finish\n");
ec1be5d5 1453 }
1454
1455 return;
1456}
50ff0bcd 1457//----------------------------------------------------------------------------
6d8df534 1458AliESDVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray,
1459 UShort_t *id,
1460 Bool_t optUseFitter,
1461 Bool_t optPropagate)
fee9ef0d 1462{
eae1677e 1463//
6d8df534 1464// Return vertex from tracks (AliExternalTrackParam) in array
eae1677e 1465//
5b78b791 1466 fCurrentVertex = 0;
fee9ef0d 1467 SetConstraintOff();
1468
50ff0bcd 1469 // get tracks and propagate them to initial vertex position
6d8df534 1470 fIdSel = new UShort_t[(Int_t)trkArray->GetEntriesFast()];
1471 Int_t nTrksSel = PrepareTracks(*trkArray,id,0);
1472 if(nTrksSel < TMath::Max(2,fMinTracks)) {
1473 TooFewTracks();
fee9ef0d 1474 return fCurrentVertex;
50ff0bcd 1475 }
1476
6d8df534 1477 // vertex finder
fee9ef0d 1478 switch (fAlgo) {
1479 case 1: StrLinVertexFinderMinDist(1); break;
1480 case 2: StrLinVertexFinderMinDist(0); break;
1481 case 3: HelixVertexFinder(); break;
1482 case 4: VertexFinder(1); break;
1483 case 5: VertexFinder(0); break;
1484 default: printf("Wrong algorithm\n"); break;
1485 }
fee9ef0d 1486 if(fDebug) printf(" Vertex finding completed\n");
1487
1488 // vertex fitter
6d8df534 1489 if(optUseFitter) {
1490 VertexFitter();
1491 } else {
fee9ef0d 1492 Double_t position[3]={fVert.GetXv(),fVert.GetYv(),fVert.GetZv()};
146c29df 1493 Double_t covmatrix[6];
1494 fVert.GetCovMatrix(covmatrix);
fee9ef0d 1495 Double_t chi2=99999.;
6d8df534 1496 Int_t nTrksUsed=fVert.GetNContributors();
1497 fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nTrksUsed);
fee9ef0d 1498 }
1499 fCurrentVertex->SetDispersion(fVert.GetDispersion());
1500
1501
1502 // set indices of used tracks and propagate track to found vertex
50ff0bcd 1503 UShort_t *indices = 0;
6d8df534 1504 Double_t d0z0[2],covd0z0[3];
1505 AliExternalTrackParam *t = 0;
fee9ef0d 1506 if(fCurrentVertex->GetNContributors()>0) {
1507 indices = new UShort_t[fCurrentVertex->GetNContributors()];
6d8df534 1508 for(Int_t jj=0; jj<(Int_t)fTrkArraySel.GetEntriesFast(); jj++) {
1509 indices[jj] = fIdSel[jj];
1510 t = (AliExternalTrackParam*)fTrkArraySel.At(jj);
1511 if(optPropagate && optUseFitter) {
fee9ef0d 1512 if(TMath::Sqrt(fCurrentVertex->GetXv()*fCurrentVertex->GetXv()+fCurrentVertex->GetYv()*fCurrentVertex->GetYv())<3.) {
6d8df534 1513 t->PropagateToDCA(fCurrentVertex,GetFieldkG(),100.,d0z0,covd0z0);
fee9ef0d 1514 if(fDebug) printf("Track %d propagated to found vertex\n",jj);
6d8df534 1515 } else {
fee9ef0d 1516 AliWarning("Found vertex outside beam pipe!");
1517 }
1518 }
50ff0bcd 1519 }
fee9ef0d 1520 fCurrentVertex->SetIndices(fCurrentVertex->GetNContributors(),indices);
50ff0bcd 1521 }
6d8df534 1522
1523 // clean up
1524 delete [] indices; indices=NULL;
1525 delete [] fIdSel; fIdSel=NULL;
1526 fTrkArraySel.Delete();
50ff0bcd 1527
fee9ef0d 1528 return fCurrentVertex;
eae1677e 1529}
50ff0bcd 1530//----------------------------------------------------------------------------
6d8df534 1531AliESDVertex* AliVertexerTracks::VertexForSelectedESDTracks(TObjArray *trkArray,
1532 Bool_t optUseFitter,
1533 Bool_t optPropagate)
fee9ef0d 1534{
50ff0bcd 1535//
6d8df534 1536// Return vertex from array of ESD tracks
50ff0bcd 1537//
5834e9a1 1538
6d8df534 1539 Int_t nTrks = (Int_t)trkArray->GetEntriesFast();
1540 UShort_t *id = new UShort_t[nTrks];
1541
1542 AliESDtrack *esdt = 0;
50ff0bcd 1543 for(Int_t i=0; i<nTrks; i++){
6d8df534 1544 esdt = (AliESDtrack*)trkArray->At(i);
1545 id[i] = (UShort_t)esdt->GetID();
50ff0bcd 1546 }
1547
6d8df534 1548 VertexForSelectedTracks(trkArray,id,optUseFitter,optPropagate);
1549
1550 delete [] id; id=NULL;
1551
1552 return fCurrentVertex;
50ff0bcd 1553}
1554//--------------------------------------------------------------------------