]> 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
EMCAL geometry can be created independently form anything now
[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//-----------------------------------------------------------------
17// Implementation of the vertexer from ESD tracks
18//
19// Origin: AliITSVertexerTracks
20// A.Dainese, Padova,
21// andrea.dainese@pd.infn.it
22// M.Masera, Torino,
23// massimo.masera@to.infn.it
24// Moved to STEER and adapted to ESD tracks:
25// F.Prino, Torino, prino@to.infn.it
26//-----------------------------------------------------------------
27
28//---- Root headers --------
29#include <TTree.h>
ec1be5d5 30#include <TMatrixD.h>
2d57349e 31//---- AliRoot headers -----
32#include "AliStrLine.h"
33#include "AliVertexerTracks.h"
c5e3e5d1 34#include "AliESD.h"
2d57349e 35#include "AliESDtrack.h"
36
37ClassImp(AliVertexerTracks)
38
39
40//----------------------------------------------------------------------------
41AliVertexerTracks::AliVertexerTracks():
07680cae 42TObject(),
43fVert(),
44fCurrentVertex(0),
0c69be49 45fMinTracks(1),
fe12e09c 46fMinITSClusters(5),
07680cae 47fTrkArray(),
48fTrksToSkip(0),
49fNTrksToSkip(0),
50fDCAcut(0),
51fAlgo(1),
fe12e09c 52fNSigma(3),
50ff0bcd 53fITSin(kTRUE),
54fITSrefit(kTRUE),
07680cae 55fDebug(0)
2d57349e 56{
57//
58// Default constructor
59//
60 SetVtxStart();
07680cae 61 SetVtxStartSigma();
2d57349e 62 SetMinTracks();
07680cae 63 SetMinITSClusters();
64 SetNSigmad0();
2d57349e 65}
66//-----------------------------------------------------------------------------
67AliVertexerTracks::AliVertexerTracks(Double_t xStart, Double_t yStart):
07680cae 68TObject(),
69fVert(),
70fCurrentVertex(0),
0c69be49 71fMinTracks(1),
fe12e09c 72fMinITSClusters(5),
07680cae 73fTrkArray(),
74fTrksToSkip(0),
75fNTrksToSkip(0),
76fDCAcut(0),
77fAlgo(1),
fe12e09c 78fNSigma(3),
50ff0bcd 79fITSin(kTRUE),
80fITSrefit(kTRUE),
07680cae 81fDebug(0)
2d57349e 82{
83//
84// Standard constructor
85//
86 SetVtxStart(xStart,yStart);
07680cae 87 SetVtxStartSigma();
2d57349e 88 SetMinTracks();
07680cae 89 SetMinITSClusters();
90 SetNSigmad0();
2d57349e 91}
92//-----------------------------------------------------------------------------
93AliVertexerTracks::~AliVertexerTracks() {
94 // Default Destructor
eae1677e 95 // The objects pointed by the following pointer are not owned
2d57349e 96 // by this class and are not deleted
ec1be5d5 97 fCurrentVertex = 0;
07680cae 98 delete[] fTrksToSkip;
ec1be5d5 99}
c5e3e5d1 100//----------------------------------------------------------------------------
ec1be5d5 101AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliESD *esdEvent)
102{
c5e3e5d1 103//
ec1be5d5 104// Primary vertex for current ESD event
eae1677e 105// (Two iterations:
106// 1st with 5*fNSigma*sigma(pt) cut w.r.t. to initial vertex;
107// 2nd with fNSigma*sigma(pt) cut w.r.t. to vertex found in 1st iteration)
318aecb9 108// All ESD tracks with inside the beam pipe are then propagated to found vertex
07680cae 109//
110 fCurrentVertex = 0;
111
50ff0bcd 112 Int_t nTrksTot = (Int_t)esdEvent->GetNumberOfTracks();
07680cae 113 TTree *trkTree = new TTree("TreeT","tracks");
114 AliESDtrack *esdTrack = 0;
115 trkTree->Branch("tracks","AliESDtrack",&esdTrack);
116
eae1677e 117 Bool_t skipThis;
50ff0bcd 118 for(Int_t i=0; i<nTrksTot; i++) {
eae1677e 119 // check tracks to skip
120 skipThis = kFALSE;
121 for(Int_t j=0; j<fNTrksToSkip; j++) {
122 if(i==fTrksToSkip[j]) {
123 if(fDebug) printf("skipping track: %d\n",i);
124 skipThis = kTRUE;
125 }
126 }
07680cae 127 AliESDtrack *et = esdEvent->GetTrack(i);
128 esdTrack = new AliESDtrack(*et);
50ff0bcd 129 if(skipThis) {delete esdTrack;continue;}
130 if(fITSin) {
131 if(!(esdTrack->GetStatus()&AliESDtrack::kITSin)) {delete esdTrack;continue;}
132 if(fITSrefit && !(esdTrack->GetStatus()&AliESDtrack::kITSrefit)) {delete esdTrack;continue;}
133 Int_t nclus=esdTrack->GetNcls(0); // check number of clusters in ITS
134 if(nclus<fMinITSClusters) {delete esdTrack;continue;}
135 }
07680cae 136 trkTree->Fill();
5834e9a1 137 delete esdTrack;
07680cae 138 }
318aecb9 139
07680cae 140
141 // ITERATION 1
142 // propagate tracks to initVertex
eae1677e 143 // preselect them (reject for |d0|>5*fNSigma*sigma w.r.t. initVertex)
50ff0bcd 144 Int_t nTrksPrep;
145 nTrksPrep = PrepareTracks(*trkTree,1);
146 if(fDebug) printf(" tracks prepared - step 1: %d\n",nTrksPrep);
147 if(nTrksPrep < fMinTracks) {
148 if(fDebug) printf("TooFewTracks\n");
eae1677e 149 TooFewTracks(esdEvent);
150 if(fDebug) fCurrentVertex->PrintStatus();
5834e9a1 151 fTrkArray.Delete();
152 delete trkTree;
07680cae 153 return fCurrentVertex;
154 }
155
0c69be49 156 if(nTrksPrep==1){
157 if(fDebug) printf("Just one track\n");
158 OneTrackVertFinder();
159 }else{
07680cae 160 // vertex finder
0c69be49 161 switch (fAlgo) {
07680cae 162 case 1: StrLinVertexFinderMinDist(1); break;
163 case 2: StrLinVertexFinderMinDist(0); break;
164 case 3: HelixVertexFinder(); break;
165 case 4: VertexFinder(1); break;
166 case 5: VertexFinder(0); break;
167 default: printf("Wrong algorithm\n"); break;
0c69be49 168 }
07680cae 169 }
170 if(fDebug) printf(" vertex finding completed\n");
171
172 // vertex fitter
eae1677e 173 VertexFitter(kTRUE);
07680cae 174 if(fDebug) printf(" vertex fit completed\n");
5834e9a1 175 fTrkArray.Delete();
07680cae 176 // ITERATION 2
177 // propagate tracks to best between initVertex and fCurrentVertex
178 // preselect tracks (reject for |d0|>fNSigma*sigma w.r.t. best
179 // between initVertex and fCurrentVertex)
50ff0bcd 180 nTrksPrep = PrepareTracks(*trkTree,2);
07680cae 181 delete trkTree;
50ff0bcd 182 if(fDebug) printf(" tracks prepared - step 2: %d\n",nTrksPrep);
183 if(nTrksPrep < fMinTracks) {
184 if(fDebug) printf("TooFewTracks\n");
eae1677e 185 TooFewTracks(esdEvent);
186 if(fDebug) fCurrentVertex->PrintStatus();
5834e9a1 187 fTrkArray.Delete();
07680cae 188 return fCurrentVertex;
189 }
190
191 // vertex finder
0c69be49 192 if(nTrksPrep==1){
193 if(fDebug) printf("Just one track\n");
194 OneTrackVertFinder();
195 }else{
196 switch (fAlgo) {
07680cae 197 case 1: StrLinVertexFinderMinDist(1); break;
198 case 2: StrLinVertexFinderMinDist(0); break;
199 case 3: HelixVertexFinder(); break;
200 case 4: VertexFinder(1); break;
201 case 5: VertexFinder(0); break;
202 default: printf("Wrong algorithm\n"); break;
0c69be49 203 }
07680cae 204 }
205 if(fDebug) printf(" vertex finding completed\n");
206
207 // fitter
eae1677e 208 VertexFitter(kTRUE);
07680cae 209 if(fDebug) printf(" vertex fit completed\n");
210
eae1677e 211
eae1677e 212 // take true pos from ESD
07680cae 213 Double_t tp[3];
214 esdEvent->GetVertex()->GetTruePos(tp);
215 fCurrentVertex->SetTruePos(tp);
50ff0bcd 216 if(fNominalCov[0]>1.) {
eae1677e 217 fCurrentVertex->SetTitle("VertexerTracksNoConstraint");
218 } else {
219 fCurrentVertex->SetTitle("VertexerTracksWithConstraint");
220 }
221
318aecb9 222
223 // propagate tracks to found vertex
224 if(TMath::Sqrt(fCurrentVertex->GetXv()*fCurrentVertex->GetXv()+fCurrentVertex->GetYv()*fCurrentVertex->GetYv())<3.) {
50ff0bcd 225 for(Int_t ii=0; ii<nTrksTot; ii++) {
318aecb9 226 AliESDtrack *et = esdEvent->GetTrack(ii);
50ff0bcd 227 if(!(et->GetStatus()&AliESDtrack::kITSin)) continue;
318aecb9 228 if(et->GetX()>3.) continue;
229 et->RelateToVertex(fCurrentVertex,GetField(),100.);
230 }
231 } else {
232 AliWarning("Found vertex outside beam pipe!");
233 }
50ff0bcd 234
235 // set indices of used tracks
236 UShort_t *indices = 0;
237 AliESDtrack *ett = 0;
238 if(fCurrentVertex->GetNContributors()>0) {
239 indices = new UShort_t[fCurrentVertex->GetNContributors()];
240 for(Int_t jj=0;jj<(Int_t)fTrkArray.GetEntriesFast();jj++) {
241 ett = (AliESDtrack*)fTrkArray.At(jj);
242 indices[jj] = (UShort_t)ett->GetID();
243 }
244 fCurrentVertex->SetIndices(fCurrentVertex->GetNContributors(),indices);
245 }
246 delete [] indices;
247 fTrkArray.Delete();
248
249 if(fTrksToSkip) delete [] fTrksToSkip;
250
251
252 if(fDebug) fCurrentVertex->PrintStatus();
253 if(fDebug) fCurrentVertex->PrintIndices();
318aecb9 254
07680cae 255 return fCurrentVertex;
256}
50ff0bcd 257//------------------------------------------------------------------------
258Double_t AliVertexerTracks::GetDeterminant3X3(Double_t matr[][3]){
259 //
260 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];
261 return det;
262}
263//-------------------------------------------------------------------------
264void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t (*m)[3],Double_t *d){
ec1be5d5 265
50ff0bcd 266 //
267 Double_t x12=p0[0]-p1[0];
268 Double_t y12=p0[1]-p1[1];
269 Double_t z12=p0[2]-p1[2];
270 Double_t kk=x12*x12+y12*y12+z12*z12;
271 m[0][0]=2-2/kk*x12*x12;
272 m[0][1]=-2/kk*x12*y12;
273 m[0][2]=-2/kk*x12*z12;
274 m[1][0]=-2/kk*x12*y12;
275 m[1][1]=2-2/kk*y12*y12;
276 m[1][2]=-2/kk*y12*z12;
277 m[2][0]=-2/kk*x12*z12;
278 m[2][1]=-2*y12*z12;
279 m[2][2]=2-2/kk*z12*z12;
280 d[0]=2*p0[0]-2/kk*p0[0]*x12*x12-2/kk*p0[2]*x12*z12-2/kk*p0[1]*x12*y12;
281 d[1]=2*p0[1]-2/kk*p0[1]*y12*y12-2/kk*p0[0]*x12*y12-2/kk*p0[2]*z12*y12;
282 d[2]=2*p0[2]-2/kk*p0[2]*z12*z12-2/kk*p0[0]*x12*z12-2/kk*p0[1]*z12*y12;
c5e3e5d1 283
50ff0bcd 284}
285//--------------------------------------------------------------------------
286void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t *sigmasq,Double_t (*m)[3],Double_t *d){
287 //
288 Double_t x12=p1[0]-p0[0];
289 Double_t y12=p1[1]-p0[1];
290 Double_t z12=p1[2]-p0[2];
c5e3e5d1 291
50ff0bcd 292 Double_t den= x12*x12*sigmasq[1]*sigmasq[2]+y12*y12*sigmasq[0]*sigmasq[2]+z12*z12*sigmasq[0]*sigmasq[1];
c5e3e5d1 293
50ff0bcd 294 Double_t kk= 2*(x12*x12/sigmasq[0]+y12*y12/sigmasq[1]+z12*z12/sigmasq[2]);
c5e3e5d1 295
50ff0bcd 296 Double_t cc[3];
297 cc[0]=-x12/sigmasq[0];
298 cc[1]=-y12/sigmasq[1];
299 cc[2]=-z12/sigmasq[2];
ec1be5d5 300
50ff0bcd 301 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 302
50ff0bcd 303 Double_t ss= -p0[0]*cc[0]-p0[1]*cc[1]-p0[2]*cc[2];
2d57349e 304
50ff0bcd 305 Double_t aa[3];
306 aa[0]=x12*sigmasq[1]*sigmasq[2]/den;
307 aa[1]=y12*sigmasq[0]*sigmasq[2]/den;
308 aa[2]=z12*sigmasq[0]*sigmasq[1]/den;
2d57349e 309
50ff0bcd 310 m[0][0]=aa[0]*(aa[0]*kk+2*cc[0])+2*cc[0]*aa[0]+2/sigmasq[0];
311 m[0][1]=aa[1]*(aa[0]*kk+2*cc[0])+2*cc[1]*aa[0];
312 m[0][2]=aa[2]*(aa[0]*kk+2*cc[0])+2*cc[2]*aa[0];
2d57349e 313
50ff0bcd 314 m[1][0]=aa[0]*(aa[1]*kk+2*cc[1])+2*cc[0]*aa[1];
315 m[1][1]=aa[1]*(aa[1]*kk+2*cc[1])+2*cc[1]*aa[1]+2/sigmasq[1];
316 m[1][2]=aa[2]*(aa[1]*kk+2*cc[1])+2*cc[2]*aa[1];
2d57349e 317
50ff0bcd 318 m[2][0]=aa[0]*(aa[2]*kk+2*cc[2])+2*cc[0]*aa[2];
319 m[2][1]=aa[1]*(aa[2]*kk+2*cc[2])+2*cc[1]*aa[2];
320 m[2][2]=aa[2]*(aa[2]*kk+2*cc[2])+2*cc[2]*aa[2]+2/sigmasq[2];
08df6187 321
50ff0bcd 322 d[0]=-ww*(aa[0]*kk+2*cc[0])-2*ss*aa[0]+2*p0[0]/sigmasq[0];
323 d[1]=-ww*(aa[1]*kk+2*cc[1])-2*ss*aa[1]+2*p0[1]/sigmasq[1];
324 d[2]=-ww*(aa[2]*kk+2*cc[2])-2*ss*aa[2]+2*p0[2]/sigmasq[2];
2d57349e 325
2d57349e 326 }
50ff0bcd 327//--------------------------------------------------------------------------
328Double_t AliVertexerTracks::GetStrLinMinDist(Double_t *p0,Double_t *p1,Double_t *x0){
329 //
330 Double_t x12=p0[0]-p1[0];
331 Double_t y12=p0[1]-p1[1];
332 Double_t z12=p0[2]-p1[2];
333 Double_t x10=p0[0]-x0[0];
334 Double_t y10=p0[1]-x0[1];
335 Double_t z10=p0[2]-x0[2];
336 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);
2d57349e 337}
338//---------------------------------------------------------------------------
0c69be49 339void AliVertexerTracks::OneTrackVertFinder() {
340 // find vertex for events with 1 track, using DCA to nominal beam axis
341 if(fDebug) printf("Number of prepared tracks =%d - Call OneTrackVertFinder",fTrkArray.GetEntries());
342 AliESDtrack *track1;
343 track1 = (AliESDtrack*)fTrkArray.At(0);
344 Double_t field=GetField();
345 Double_t alpha=track1->GetAlpha();
346 Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
347 Double_t pos[3],dir[3];
348 track1->GetXYZAt(mindist,field,pos);
349 track1->GetPxPyPzAt(mindist,field,dir);
350 AliStrLine *line1 = new AliStrLine(pos,dir);
351 Double_t p1[3]={fNominalPos[0],fNominalPos[1],0.};
352 Double_t p2[3]={fNominalPos[0],fNominalPos[1],10.};
353 AliStrLine *zeta=new AliStrLine(p1,p2,kTRUE);
354 Double_t crosspoint[3]={0.,0.,0.};
355 Double_t sigma=999.;
356 Int_t nContrib=-1;
357 Int_t retcode = zeta->Cross(line1,crosspoint);
358 if(retcode>=0){
359 sigma=line1->GetDistFromPoint(crosspoint);
360 nContrib=1;
361 }
362 delete zeta;
363 delete line1;
364 fVert.SetXYZ(crosspoint);
365 fVert.SetDispersion(sigma);
366 fVert.SetNContributors(nContrib);
367}
368//---------------------------------------------------------------------------
2d57349e 369void AliVertexerTracks::HelixVertexFinder() {
370
371 // Get estimate of vertex position in (x,y) from tracks DCA
372
373
374 Double_t initPos[3];
375 initPos[2] = 0.;
376 for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
377 Double_t field=GetField();
378
379 Int_t nacc = (Int_t)fTrkArray.GetEntriesFast();
380
381 Double_t aver[3]={0.,0.,0.};
382 Double_t averquad[3]={0.,0.,0.};
383 Double_t sigmaquad[3]={0.,0.,0.};
384 Double_t sigma=0;
385 Int_t ncombi = 0;
386 AliESDtrack *track1;
387 AliESDtrack *track2;
388 Double_t distCA;
389 Double_t x, par[5];
390 Double_t alpha, cs, sn;
391 Double_t crosspoint[3];
392 for(Int_t i=0; i<nacc; i++){
393 track1 = (AliESDtrack*)fTrkArray.At(i);
394
395
396 for(Int_t j=i+1; j<nacc; j++){
397 track2 = (AliESDtrack*)fTrkArray.At(j);
398
399 distCA=track2->PropagateToDCA(track1,field);
400
401 if(fDCAcut<=0 ||(fDCAcut>0&&distCA<fDCAcut)){
402 track1->GetExternalParameters(x,par);
403 alpha=track1->GetAlpha();
404 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
405 Double_t x1=x*cs - par[0]*sn;
406 Double_t y1=x*sn + par[0]*cs;
407 Double_t z1=par[1];
408 Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
409 track2->GetExternalParameters(x,par);
410 alpha=track2->GetAlpha();
411 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
412 Double_t x2=x*cs - par[0]*sn;
413 Double_t y2=x*sn + par[0]*cs;
414 Double_t z2=par[1];
415 Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
416 Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
417 Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
418 Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
419 Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
420 crosspoint[0]=wx1*x1 + wx2*x2;
421 crosspoint[1]=wy1*y1 + wy2*y2;
422 crosspoint[2]=wz1*z1 + wz2*z2;
423
424 ncombi++;
425 for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
426 for(Int_t jj=0;jj<3;jj++)averquad[jj]+=(crosspoint[jj]*crosspoint[jj]);
427 }
428 }
429
430 }
431 if(ncombi>0){
432 for(Int_t jj=0;jj<3;jj++){
433 initPos[jj] = aver[jj]/ncombi;
434 averquad[jj]/=ncombi;
435 sigmaquad[jj]=averquad[jj]-initPos[jj]*initPos[jj];
436 sigma+=sigmaquad[jj];
437 }
438 sigma=TMath::Sqrt(TMath::Abs(sigma));
439 }
440 else {
441 Warning("HelixVertexFinder","Finder did not succed");
442 sigma=999;
443 }
444 fVert.SetXYZ(initPos);
445 fVert.SetDispersion(sigma);
446 fVert.SetNContributors(ncombi);
447}
50ff0bcd 448//----------------------------------------------------------------------------
449Int_t AliVertexerTracks::PrepareTracks(TTree &trkTree,Int_t OptImpParCut) {
450//
451// Propagate tracks to initial vertex position and store them in a TObjArray
452//
453 Double_t maxd0rphi = 3.;
454 Int_t nTrks = 0;
455 Double_t sigmaCurr[3];
456 Float_t d0z0[2],covd0z0[3];
457 Double_t sigma;
458 Double_t field=GetField();
459
460 AliESDVertex *initVertex = new AliESDVertex(fNominalPos,fNominalCov,1,1);
461
462 Int_t nEntries = (Int_t)trkTree.GetEntries();
463 if(!fTrkArray.IsEmpty()) fTrkArray.Delete();
464
465 if(fDebug) {
466 printf(" PrepareTracks()\n");
467 }
468
469 for(Int_t i=0; i<nEntries; i++) {
470 AliESDtrack *track = new AliESDtrack;
471 trkTree.SetBranchAddress("tracks",&track);
472 trkTree.GetEvent(i);
473
474 // propagate track to vertex
475 if(OptImpParCut==1) { // OptImpParCut==1
476 track->RelateToVertex(initVertex,field,100.);
477 } else { // OptImpParCut==2
478 fCurrentVertex->GetSigmaXYZ(sigmaCurr);
479 if((sigmaCurr[0]+sigmaCurr[1])<(TMath::Sqrt(fNominalCov[0])+TMath::Sqrt(fNominalCov[2]))) {
480 track->RelateToVertex(fCurrentVertex,field,100.);
481 } else {
482 track->RelateToVertex(initVertex,field,100.);
483 }
484 }
485
486 // select tracks with d0rphi < maxd0rphi
487 track->GetImpactParameters(d0z0,covd0z0);
488 sigma = TMath::Sqrt(covd0z0[0]);
489 maxd0rphi = fNSigma*sigma;
490 if(OptImpParCut==1) maxd0rphi *= 5.;
491
492 if(fDebug) printf("trk %d; lab %d; |d0| = %f; cut = %f\n",i,track->GetLabel(),TMath::Abs(d0z0[0]),maxd0rphi);
493 if(TMath::Abs(d0z0[0]) > maxd0rphi) {
494 if(fDebug) printf(" rejected\n");
495 delete track; continue;
496 }
497
498 fTrkArray.AddLast(track);
499 nTrks++;
500 }
501
502 delete initVertex;
503
504 return nTrks;
505}
506//---------------------------------------------------------------------------
507void AliVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped) {
508//
509// Mark the tracks not ot be used in the vertex finding
510//
511 fNTrksToSkip = n;
512 fTrksToSkip = new Int_t[n];
513 for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i];
514 return;
515}
516//---------------------------------------------------------------------------
517void AliVertexerTracks::SetVtxStart(AliESDVertex *vtx) {
518//
519// Set initial vertex knowledge
520//
521 vtx->GetXYZ(fNominalPos);
522 vtx->GetCovMatrix(fNominalCov);
523 return;
524}
2d57349e 525//---------------------------------------------------------------------------
526void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights){
527
528 // Calculate the point at minimum distance to prepared tracks
529
530 Double_t initPos[3];
531 initPos[2] = 0.;
532 Double_t sigma=0;
533 for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
534 const Int_t knacc = (Int_t)fTrkArray.GetEntriesFast();
535 Double_t field=GetField();
536
537 AliESDtrack *track1;
538 Double_t (*vectP0)[3]=new Double_t [knacc][3];
539 Double_t (*vectP1)[3]=new Double_t [knacc][3];
540
541 Double_t sum[3][3];
542 Double_t dsum[3]={0,0,0};
543 for(Int_t i=0;i<3;i++)
544 for(Int_t j=0;j<3;j++)sum[i][j]=0;
545 for(Int_t i=0; i<knacc; i++){
546 track1 = (AliESDtrack*)fTrkArray.At(i);
547 Double_t alpha=track1->GetAlpha();
548 Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
ec1be5d5 549 Double_t pos[3],dir[3];
550 track1->GetXYZAt(mindist,field,pos);
551 track1->GetPxPyPzAt(mindist,field,dir);
552 AliStrLine *line1 = new AliStrLine(pos,dir);
553 // AliStrLine *line1 = new AliStrLine();
554 // track1->ApproximateHelixWithLine(mindist,field,line1);
2d57349e 555
556 Double_t p0[3],cd[3];
557 line1->GetP0(p0);
558 line1->GetCd(cd);
559 Double_t p1[3]={p0[0]+cd[0],p0[1]+cd[1],p0[2]+cd[2]};
560 vectP0[i][0]=p0[0];
561 vectP0[i][1]=p0[1];
562 vectP0[i][2]=p0[2];
563 vectP1[i][0]=p1[0];
564 vectP1[i][1]=p1[1];
565 vectP1[i][2]=p1[2];
566
567 Double_t matr[3][3];
568 Double_t dknow[3];
569 if(optUseWeights==0)GetStrLinDerivMatrix(p0,p1,matr,dknow);
570 if(optUseWeights==1){
571 Double_t sigmasq[3];
572 sigmasq[0]=track1->GetSigmaY2();
573 sigmasq[1]=track1->GetSigmaY2();
574 sigmasq[2]=track1->GetSigmaZ2();
575 GetStrLinDerivMatrix(p0,p1,sigmasq,matr,dknow);
576 }
577
578 for(Int_t iii=0;iii<3;iii++){
579 dsum[iii]+=dknow[iii];
580 for(Int_t lj=0;lj<3;lj++) sum[iii][lj]+=matr[iii][lj];
581 }
582 delete line1;
583 }
584
585 Double_t vett[3][3];
586 Double_t det=GetDeterminant3X3(sum);
587
588 if(det!=0){
589 for(Int_t zz=0;zz<3;zz++){
590 for(Int_t ww=0;ww<3;ww++){
591 for(Int_t kk=0;kk<3;kk++) vett[ww][kk]=sum[ww][kk];
592 }
593 for(Int_t kk=0;kk<3;kk++) vett[kk][zz]=dsum[kk];
594 initPos[zz]=GetDeterminant3X3(vett)/det;
595 }
596
597
50ff0bcd 598 for(Int_t i=0; i<knacc; i++){
599 Double_t p0[3]={0,0,0},p1[3]={0,0,0};
600 for(Int_t ii=0;ii<3;ii++){
601 p0[ii]=vectP0[i][ii];
602 p1[ii]=vectP1[i][ii];
603 }
604 sigma+=GetStrLinMinDist(p0,p1,initPos);
605 }
2d57349e 606
50ff0bcd 607 sigma=TMath::Sqrt(sigma);
608 }else{
609 Warning("StrLinVertexFinderMinDist","Finder did not succed");
610 sigma=999;
611 }
612 delete vectP0;
613 delete vectP1;
614 fVert.SetXYZ(initPos);
615 fVert.SetDispersion(sigma);
616 fVert.SetNContributors(knacc);
617}
618//---------------------------------------------------------------------------
619void AliVertexerTracks::TooFewTracks(const AliESD* esdEvent) {
620//
621// When the number of tracks is < fMinTracks
622//
2d57349e 623
50ff0bcd 624 // deal with vertices not found
625 Double_t pos[3],err[3];
626 Int_t ncontr=0;
627 pos[0] = fNominalPos[0];
628 err[0] = TMath::Sqrt(fNominalCov[0]);
629 pos[1] = fNominalPos[1];
630 err[1] = TMath::Sqrt(fNominalCov[2]);
631 pos[2] = esdEvent->GetVertex()->GetZv();
632 err[2] = esdEvent->GetVertex()->GetZRes();
633 if(err[0]>1. && esdEvent->GetVertex()->GetNContributors()<=0)
634 ncontr = -1; // (x,y,z) = (0,0,0)
635 if(err[0]>1. && esdEvent->GetVertex()->GetNContributors()>0)
636 ncontr = -2; // (x,y,z) = (0,0,z_fromSPD)
637 if(err[0]<1. && esdEvent->GetVertex()->GetNContributors()<=0)
638 ncontr = -3; // (x,y,z) = (x_mean,y_mean,0)
639 if(err[0]<1. && esdEvent->GetVertex()->GetNContributors()>0)
640 ncontr = -4; // (x,y,z) = (x_mean,y_mean,z_fromSPD)
641 fCurrentVertex = 0;
642 fCurrentVertex = new AliESDVertex(pos,err);
643 fCurrentVertex->SetNContributors(ncontr);
2d57349e 644
50ff0bcd 645 Double_t tp[3];
646 esdEvent->GetVertex()->GetTruePos(tp);
647 fCurrentVertex->SetTruePos(tp);
648 fCurrentVertex->SetTitle("VertexerTracksWithConstraint");
649 if(ncontr==-1||ncontr==-2)
650 fCurrentVertex->SetTitle("VertexerTracksNoConstraint");
2d57349e 651
50ff0bcd 652 return;
653}
654//---------------------------------------------------------------------------
655void AliVertexerTracks::VertexFinder(Int_t optUseWeights) {
2d57349e 656
50ff0bcd 657 // Get estimate of vertex position in (x,y) from tracks DCA
658
659 Double_t initPos[3];
660 initPos[2] = 0.;
661 for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
662 Int_t nacc = (Int_t)fTrkArray.GetEntriesFast();
663 Double_t aver[3]={0.,0.,0.};
664 Double_t aversq[3]={0.,0.,0.};
665 Double_t sigmasq[3]={0.,0.,0.};
666 Double_t sigma=0;
667 Int_t ncombi = 0;
668 AliESDtrack *track1;
669 AliESDtrack *track2;
670 Double_t pos[3],dir[3];
671 Double_t alpha,mindist;
672 Double_t field=GetField();
2d57349e 673
50ff0bcd 674 for(Int_t i=0; i<nacc; i++){
675 track1 = (AliESDtrack*)fTrkArray.At(i);
676 alpha=track1->GetAlpha();
677 mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
678 track1->GetXYZAt(mindist,field,pos);
679 track1->GetPxPyPzAt(mindist,field,dir);
680 AliStrLine *line1 = new AliStrLine(pos,dir);
2d57349e 681
50ff0bcd 682 // AliStrLine *line1 = new AliStrLine();
683 // track1->ApproximateHelixWithLine(mindist,field,line1);
684
685 for(Int_t j=i+1; j<nacc; j++){
686 track2 = (AliESDtrack*)fTrkArray.At(j);
687 alpha=track2->GetAlpha();
688 mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
689 track2->GetXYZAt(mindist,field,pos);
690 track2->GetPxPyPzAt(mindist,field,dir);
691 AliStrLine *line2 = new AliStrLine(pos,dir);
692 // AliStrLine *line2 = new AliStrLine();
693 // track2->ApproximateHelixWithLine(mindist,field,line2);
694 Double_t distCA=line2->GetDCA(line1);
695 if(fDCAcut<=0 || (fDCAcut>0&&distCA<fDCAcut)){
696 Double_t pnt1[3],pnt2[3],crosspoint[3];
2d57349e 697
50ff0bcd 698 if(optUseWeights<=0){
699 Int_t retcode = line2->Cross(line1,crosspoint);
700 if(retcode>=0){
701 ncombi++;
702 for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
703 for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
704 }
705 }
706 if(optUseWeights>0){
707 Int_t retcode = line1->CrossPoints(line2,pnt1,pnt2);
708 if(retcode>=0){
709 Double_t alpha, cs, sn;
710 alpha=track1->GetAlpha();
711 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
712 Double_t sx1=sn*sn*track1->GetSigmaY2(), sy1=cs*cs*track1->GetSigmaY2();
713 alpha=track2->GetAlpha();
714 cs=TMath::Cos(alpha); sn=TMath::Sin(alpha);
715 Double_t sx2=sn*sn*track2->GetSigmaY2(), sy2=cs*cs*track2->GetSigmaY2();
716 Double_t sz1=track1->GetSigmaZ2(), sz2=track2->GetSigmaZ2();
717 Double_t wx1=sx2/(sx1+sx2), wx2=1.- wx1;
718 Double_t wy1=sy2/(sy1+sy2), wy2=1.- wy1;
719 Double_t wz1=sz2/(sz1+sz2), wz2=1.- wz1;
720 crosspoint[0]=wx1*pnt1[0] + wx2*pnt2[0];
721 crosspoint[1]=wy1*pnt1[1] + wy2*pnt2[1];
722 crosspoint[2]=wz1*pnt1[2] + wz2*pnt2[2];
723
724 ncombi++;
725 for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
726 for(Int_t jj=0;jj<3;jj++)aversq[jj]+=(crosspoint[jj]*crosspoint[jj]);
727 }
728 }
729 }
730 delete line2;
731 }
732 delete line1;
71d84967 733 }
50ff0bcd 734 if(ncombi>0){
735 for(Int_t jj=0;jj<3;jj++){
736 initPos[jj] = aver[jj]/ncombi;
737 aversq[jj]/=ncombi;
738 sigmasq[jj]=aversq[jj]-initPos[jj]*initPos[jj];
739 sigma+=sigmasq[jj];
740 }
741 sigma=TMath::Sqrt(TMath::Abs(sigma));
742 }
743 else {
744 Warning("VertexFinder","Finder did not succed");
745 sigma=999;
746 }
747 fVert.SetXYZ(initPos);
748 fVert.SetDispersion(sigma);
749 fVert.SetNContributors(ncombi);
2d57349e 750}
ec1be5d5 751//---------------------------------------------------------------------------
07680cae 752void AliVertexerTracks::VertexFitter(Bool_t useNominalVtx) {
ec1be5d5 753//
754// The optimal estimate of the vertex position is given by a "weighted
755// average of tracks positions"
756// Original method: CMS Note 97/0051
757//
758 Double_t initPos[3];
759 fVert.GetXYZ(initPos);
760 if(fDebug) {
761 printf(" VertexFitter(): start\n");
762 printf(" Number of tracks in array: %d\n",(Int_t)fTrkArray.GetEntriesFast());
763 printf(" Minimum # tracks required in fit: %d\n",fMinTracks);
764 printf("Vertex position after finder: %f,%f,%f\n",initPos[0],initPos[1],initPos[2]);
50ff0bcd 765 if(useNominalVtx) 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 766 }
767
768 Int_t i,j,k,step=0;
769 TMatrixD rv(3,1);
770 TMatrixD vV(3,3);
771 rv(0,0) = initPos[0];
772 rv(1,0) = initPos[1];
773 rv(2,0) = 0.;
774 Double_t xlStart,alpha;
775 Double_t rotAngle;
776 Double_t cosRot,sinRot;
777 Double_t cc[15];
778 Int_t nUsedTrks;
779 Double_t chi2,chi2i;
780 Int_t arrEntries = (Int_t)fTrkArray.GetEntries();
781 AliESDtrack *t = 0;
782 Int_t failed = 0;
783
50ff0bcd 784 // initial vertex covariance matrix
785 TMatrixD vVb(3,3);
786 vVb(0,0) = fNominalCov[0];
787 vVb(0,1) = fNominalCov[1];
788 vVb(0,2) = 0.;
789 vVb(1,0) = fNominalCov[1];
790 vVb(1,1) = fNominalCov[2];
791 vVb(1,2) = 0.;
792 vVb(2,0) = 0.;
793 vVb(2,1) = 0.;
794 vVb(2,2) = fNominalCov[5];
795 TMatrixD vVbInv(TMatrixD::kInverted,vVb);
796 TMatrixD rb(3,1);
797 rb(0,0) = fNominalPos[0];
798 rb(1,0) = fNominalPos[1];
799 rb(2,0) = fNominalPos[2];
800 TMatrixD vVbInvrb(vVbInv,TMatrixD::kMult,rb);
801
802
eae1677e 803 // 2 steps:
804 // 1st - estimate of vtx using all tracks
805 // 2nd - estimate of global chi2
806 for(step=0; step<2; step++) {
ec1be5d5 807 if(fDebug) printf(" step = %d\n",step);
808 chi2 = 0.;
809 nUsedTrks = 0;
810
811 TMatrixD sumWiri(3,1);
812 TMatrixD sumWi(3,3);
813 for(i=0; i<3; i++) {
814 sumWiri(i,0) = 0.;
815 for(j=0; j<3; j++) sumWi(j,i) = 0.;
816 }
817
50ff0bcd 818
07680cae 819 if(useNominalVtx) {
50ff0bcd 820 for(i=0;i<3;i++) {
821 sumWiri(i,0) += vVbInvrb(i,0);
822 for(k=0;k<3;k++) sumWi(i,k) += vVbInv(i,k);
07680cae 823 }
824 }
825
826
ec1be5d5 827 // loop on tracks
828 for(k=0; k<arrEntries; k++) {
ec1be5d5 829 // get track from track array
830 t = (AliESDtrack*)fTrkArray.At(k);
831 alpha = t->GetAlpha();
832 xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha);
6c94f330 833 t->AliExternalTrackParam::PropagateTo(xlStart,AliTracker::GetBz()); // to vtxSeed
ec1be5d5 834 rotAngle = alpha;
835 if(alpha<0.) rotAngle += 2.*TMath::Pi();
836 cosRot = TMath::Cos(rotAngle);
837 sinRot = TMath::Sin(rotAngle);
838
839 // vector of track global coordinates
840 TMatrixD ri(3,1);
841 ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot;
842 ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot;
843 ri(2,0) = t->GetZ();
844
845 // matrix to go from global (x,y,z) to local (y,z);
846 TMatrixD qQi(2,3);
847 qQi(0,0) = -sinRot;
848 qQi(0,1) = cosRot;
849 qQi(0,2) = 0.;
850 qQi(1,0) = 0.;
851 qQi(1,1) = 0.;
852 qQi(1,2) = 1.;
853
854 // covariance matrix of local (y,z) - inverted
855 TMatrixD uUi(2,2);
856 t->GetExternalCovariance(cc);
857 uUi(0,0) = cc[0];
858 uUi(0,1) = cc[1];
859 uUi(1,0) = cc[1];
860 uUi(1,1) = cc[2];
861
862 // weights matrix: wWi = qQiT * uUiInv * qQi
863 if(uUi.Determinant() <= 0.) continue;
864 TMatrixD uUiInv(TMatrixD::kInverted,uUi);
865 TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
866 TMatrixD wWi(qQi,TMatrixD::kTransposeMult,uUiInvQi);
867
868 // track chi2
869 TMatrixD deltar = rv; deltar -= ri;
870 TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
871 chi2i = deltar(0,0)*wWideltar(0,0)+
872 deltar(1,0)*wWideltar(1,0)+
873 deltar(2,0)*wWideltar(2,0);
874
875
ec1be5d5 876 // add to total chi2
877 chi2 += chi2i;
878
879 TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
880
881 sumWiri += wWiri;
882 sumWi += wWi;
883
884 nUsedTrks++;
885 } // end loop on tracks
886
887 if(nUsedTrks < fMinTracks) {
888 failed=1;
889 continue;
890 }
891
892 Double_t determinant = sumWi.Determinant();
893 //cerr<<" determinant: "<<determinant<<endl;
894 if(determinant < 100.) {
895 printf("det(V) = 0\n");
896 failed=1;
897 continue;
898 }
899
900 // inverted of weights matrix
901 TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
902 vV = invsumWi;
903
904 // position of primary vertex
905 rv.Mult(vV,sumWiri);
906
eae1677e 907 } // end loop on the 2 steps
ec1be5d5 908
5834e9a1 909 // delete t;
ec1be5d5 910
911 if(failed) {
50ff0bcd 912 if(fDebug) printf("TooFewTracks\n");
ec1be5d5 913 fCurrentVertex = new AliESDVertex(0.,0.,-1);
914 return;
915 }
916
917 Double_t position[3];
918 position[0] = rv(0,0);
919 position[1] = rv(1,0);
920 position[2] = rv(2,0);
921 Double_t covmatrix[6];
922 covmatrix[0] = vV(0,0);
923 covmatrix[1] = vV(0,1);
924 covmatrix[2] = vV(1,1);
925 covmatrix[3] = vV(0,2);
926 covmatrix[4] = vV(1,2);
927 covmatrix[5] = vV(2,2);
928
929 // store data in the vertex object
930 fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nUsedTrks);
931
932 if(fDebug) {
933 printf(" VertexFitter(): finish\n");
934 printf(" rv = ( %f , %f , %f )\n\n",rv(0,0),rv(1,0),rv(2,0));
935 fCurrentVertex->PrintStatus();
936 }
937
938 return;
939}
50ff0bcd 940//----------------------------------------------------------------------------
941AliVertex* AliVertexerTracks::VertexForSelectedTracks(TTree *trkTree) {
eae1677e 942//
50ff0bcd 943// Return vertex from tracks in trkTree
eae1677e 944//
945
50ff0bcd 946 // get tracks and propagate them to initial vertex position
947 Int_t nTrksPrep = PrepareTracks(*trkTree,1);
0c69be49 948 if(nTrksPrep < TMath::Max(2,fMinTracks) ) {
50ff0bcd 949 if(fDebug) printf("TooFewTracks\n");
950 Double_t vtx[3]={0,0,0};
951 fVert.SetXYZ(vtx);
952 fVert.SetDispersion(999);
953 fVert.SetNContributors(-5);
954 fTrkArray.Delete();
955 return &fVert;
956 }
957
958 // Set initial vertex position from ESD
959 if(fAlgo==1) StrLinVertexFinderMinDist(1);
960 if(fAlgo==2) StrLinVertexFinderMinDist(0);
961 if(fAlgo==3) HelixVertexFinder();
962 if(fAlgo==4) VertexFinder(1);
963 if(fAlgo==5) VertexFinder(0);
eae1677e 964
50ff0bcd 965
966 // set indices of used tracks
967 UShort_t *indices = 0;
968 AliESDtrack *eta = 0;
969 if(fVert.GetNContributors()>0) {
970 indices = new UShort_t[fVert.GetNContributors()];
971 for(Int_t jj=0;jj<(Int_t)fTrkArray.GetEntriesFast();jj++) {
972 eta = (AliESDtrack*)fTrkArray.At(jj);
973 indices[jj] = (UShort_t)eta->GetID();
974 }
975 fVert.SetIndices(fVert.GetNContributors(),indices);
976 }
977 delete [] indices;
978 fTrkArray.Delete();
979
980 return &fVert;
eae1677e 981}
50ff0bcd 982//----------------------------------------------------------------------------
983AliVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray) {
984//
985// Return vertex from array of tracks
986//
5834e9a1 987
50ff0bcd 988 // get tracks and propagate them to initial vertex position
989 Int_t nTrks = trkArray->GetEntriesFast();
0c69be49 990 if(nTrks < TMath::Max(2,fMinTracks) ) {
50ff0bcd 991 if(fDebug) printf("TooFewTracks\n");
992 Double_t vtx[3]={0,0,0};
993 fVert.SetXYZ(vtx);
994 fVert.SetDispersion(999);
995 fVert.SetNContributors(-5);
996 fTrkArray.Delete();
997 return &fVert;
998 }
999 TTree *trkTree = new TTree("TreeT","tracks");
1000 AliESDtrack *esdTrack = 0;
1001 trkTree->Branch("tracks","AliESDtrack",&esdTrack);
1002 for(Int_t i=0; i<nTrks; i++){
1003 esdTrack = (AliESDtrack*)trkArray->At(i);
1004 trkTree->Fill();
1005 }
1006
1007 AliVertex *vtx = VertexForSelectedTracks(trkTree);
1008 delete trkTree;
1009 return vtx;
1010}
1011//--------------------------------------------------------------------------