1 /**************************************************************************
2 * Copyright(c) 1998-2003, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
16 //-----------------------------------------------------------------
17 // Implementation of the vertexer from ESD tracks
19 // Origin: AliITSVertexerTracks
21 // andrea.dainese@pd.infn.it
23 // massimo.masera@to.infn.it
24 // Moved to STEER and adapted to ESD tracks:
25 // F.Prino, Torino, prino@to.infn.it
26 //-----------------------------------------------------------------
28 //---- Root headers --------
31 //---- AliRoot headers -----
32 #include "AliStrLine.h"
33 #include "AliVertexerTracks.h"
35 #include "AliESDtrack.h"
37 ClassImp(AliVertexerTracks)
40 //----------------------------------------------------------------------------
41 AliVertexerTracks::AliVertexerTracks():
58 // Default constructor
66 //-----------------------------------------------------------------------------
67 AliVertexerTracks::AliVertexerTracks(Double_t xStart, Double_t yStart):
84 // Standard constructor
86 SetVtxStart(xStart,yStart);
92 //-----------------------------------------------------------------------------
93 AliVertexerTracks::~AliVertexerTracks() {
95 // The objects pointed by the following pointer are not owned
96 // by this class and are not deleted
100 //----------------------------------------------------------------------------
101 AliESDVertex* AliVertexerTracks::FindPrimaryVertex(const AliESD *esdEvent)
104 // Primary vertex for current ESD event
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)
108 // All ESD tracks with inside the beam pipe are then propagated to found vertex
112 Int_t nTrksTot = (Int_t)esdEvent->GetNumberOfTracks();
113 TTree *trkTree = new TTree("TreeT","tracks");
114 AliESDtrack *esdTrack = 0;
115 trkTree->Branch("tracks","AliESDtrack",&esdTrack);
118 for(Int_t i=0; i<nTrksTot; i++) {
119 // check tracks to skip
121 for(Int_t j=0; j<fNTrksToSkip; j++) {
122 if(i==fTrksToSkip[j]) {
123 if(fDebug) printf("skipping track: %d\n",i);
127 AliESDtrack *et = esdEvent->GetTrack(i);
128 esdTrack = new AliESDtrack(*et);
129 if(skipThis) {delete esdTrack;continue;}
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;}
142 // propagate tracks to initVertex
143 // preselect them (reject for |d0|>5*fNSigma*sigma w.r.t. initVertex)
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");
149 TooFewTracks(esdEvent);
150 if(fDebug) fCurrentVertex->PrintStatus();
153 return fCurrentVertex;
157 if(fDebug) printf("Just one track\n");
158 OneTrackVertFinder();
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;
170 if(fDebug) printf(" vertex finding completed\n");
174 if(fDebug) printf(" vertex fit completed\n");
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)
180 nTrksPrep = PrepareTracks(*trkTree,2);
182 if(fDebug) printf(" tracks prepared - step 2: %d\n",nTrksPrep);
183 if(nTrksPrep < fMinTracks) {
184 if(fDebug) printf("TooFewTracks\n");
185 TooFewTracks(esdEvent);
186 if(fDebug) fCurrentVertex->PrintStatus();
188 return fCurrentVertex;
193 if(fDebug) printf("Just one track\n");
194 OneTrackVertFinder();
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;
205 if(fDebug) printf(" vertex finding completed\n");
209 if(fDebug) printf(" vertex fit completed\n");
212 // take true pos from ESD
214 esdEvent->GetVertex()->GetTruePos(tp);
215 fCurrentVertex->SetTruePos(tp);
216 if(fNominalCov[0]>1.) {
217 fCurrentVertex->SetTitle("VertexerTracksNoConstraint");
219 fCurrentVertex->SetTitle("VertexerTracksWithConstraint");
223 // propagate tracks to found vertex
224 if(TMath::Sqrt(fCurrentVertex->GetXv()*fCurrentVertex->GetXv()+fCurrentVertex->GetYv()*fCurrentVertex->GetYv())<3.) {
225 for(Int_t ii=0; ii<nTrksTot; ii++) {
226 AliESDtrack *et = esdEvent->GetTrack(ii);
227 if(!(et->GetStatus()&AliESDtrack::kITSin)) continue;
228 if(et->GetX()>3.) continue;
229 et->RelateToVertex(fCurrentVertex,GetField(),100.);
232 AliWarning("Found vertex outside beam pipe!");
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();
244 fCurrentVertex->SetIndices(fCurrentVertex->GetNContributors(),indices);
249 if(fTrksToSkip) delete [] fTrksToSkip;
252 if(fDebug) fCurrentVertex->PrintStatus();
253 if(fDebug) fCurrentVertex->PrintIndices();
255 return fCurrentVertex;
257 //------------------------------------------------------------------------
258 Double_t AliVertexerTracks::GetDeterminant3X3(Double_t matr[][3]){
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];
263 //-------------------------------------------------------------------------
264 void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t (*m)[3],Double_t *d){
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;
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;
285 //--------------------------------------------------------------------------
286 void AliVertexerTracks::GetStrLinDerivMatrix(Double_t *p0,Double_t *p1,Double_t *sigmasq,Double_t (*m)[3],Double_t *d){
288 Double_t x12=p1[0]-p0[0];
289 Double_t y12=p1[1]-p0[1];
290 Double_t z12=p1[2]-p0[2];
292 Double_t den= x12*x12*sigmasq[1]*sigmasq[2]+y12*y12*sigmasq[0]*sigmasq[2]+z12*z12*sigmasq[0]*sigmasq[1];
294 Double_t kk= 2*(x12*x12/sigmasq[0]+y12*y12/sigmasq[1]+z12*z12/sigmasq[2]);
297 cc[0]=-x12/sigmasq[0];
298 cc[1]=-y12/sigmasq[1];
299 cc[2]=-z12/sigmasq[2];
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;
303 Double_t ss= -p0[0]*cc[0]-p0[1]*cc[1]-p0[2]*cc[2];
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;
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];
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];
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];
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];
327 //--------------------------------------------------------------------------
328 Double_t AliVertexerTracks::GetStrLinMinDist(Double_t *p0,Double_t *p1,Double_t *x0){
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);
338 //---------------------------------------------------------------------------
339 void 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());
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.};
357 Int_t retcode = zeta->Cross(line1,crosspoint);
359 sigma=line1->GetDistFromPoint(crosspoint);
364 fVert.SetXYZ(crosspoint);
365 fVert.SetDispersion(sigma);
366 fVert.SetNContributors(nContrib);
368 //---------------------------------------------------------------------------
369 void AliVertexerTracks::HelixVertexFinder() {
371 // Get estimate of vertex position in (x,y) from tracks DCA
376 for(Int_t i=0;i<2;i++)initPos[i]=fNominalPos[i];
377 Double_t field=GetField();
379 Int_t nacc = (Int_t)fTrkArray.GetEntriesFast();
381 Double_t aver[3]={0.,0.,0.};
382 Double_t averquad[3]={0.,0.,0.};
383 Double_t sigmaquad[3]={0.,0.,0.};
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);
396 for(Int_t j=i+1; j<nacc; j++){
397 track2 = (AliESDtrack*)fTrkArray.At(j);
399 distCA=track2->PropagateToDCA(track1,field);
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;
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;
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;
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]);
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];
438 sigma=TMath::Sqrt(TMath::Abs(sigma));
441 Warning("HelixVertexFinder","Finder did not succed");
444 fVert.SetXYZ(initPos);
445 fVert.SetDispersion(sigma);
446 fVert.SetNContributors(ncombi);
448 //----------------------------------------------------------------------------
449 Int_t AliVertexerTracks::PrepareTracks(TTree &trkTree,Int_t OptImpParCut) {
451 // Propagate tracks to initial vertex position and store them in a TObjArray
453 Double_t maxd0rphi = 3.;
455 Double_t sigmaCurr[3];
456 Float_t d0z0[2],covd0z0[3];
458 Double_t field=GetField();
460 AliESDVertex *initVertex = new AliESDVertex(fNominalPos,fNominalCov,1,1);
462 Int_t nEntries = (Int_t)trkTree.GetEntries();
463 if(!fTrkArray.IsEmpty()) fTrkArray.Delete();
466 printf(" PrepareTracks()\n");
469 for(Int_t i=0; i<nEntries; i++) {
470 AliESDtrack *track = new AliESDtrack;
471 trkTree.SetBranchAddress("tracks",&track);
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.);
482 track->RelateToVertex(initVertex,field,100.);
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.;
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;
498 fTrkArray.AddLast(track);
506 //---------------------------------------------------------------------------
507 void AliVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped) {
509 // Mark the tracks not ot be used in the vertex finding
512 fTrksToSkip = new Int_t[n];
513 for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i];
516 //---------------------------------------------------------------------------
517 void AliVertexerTracks::SetVtxStart(AliESDVertex *vtx) {
519 // Set initial vertex knowledge
521 vtx->GetXYZ(fNominalPos);
522 vtx->GetCovMatrix(fNominalCov);
525 //---------------------------------------------------------------------------
526 void AliVertexerTracks::StrLinVertexFinderMinDist(Int_t optUseWeights){
528 // Calculate the point at minimum distance to prepared tracks
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();
538 Double_t (*vectP0)[3]=new Double_t [knacc][3];
539 Double_t (*vectP1)[3]=new Double_t [knacc][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];
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);
556 Double_t p0[3],cd[3];
559 Double_t p1[3]={p0[0]+cd[0],p0[1]+cd[1],p0[2]+cd[2]};
569 if(optUseWeights==0)GetStrLinDerivMatrix(p0,p1,matr,dknow);
570 if(optUseWeights==1){
572 sigmasq[0]=track1->GetSigmaY2();
573 sigmasq[1]=track1->GetSigmaY2();
574 sigmasq[2]=track1->GetSigmaZ2();
575 GetStrLinDerivMatrix(p0,p1,sigmasq,matr,dknow);
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];
586 Double_t det=GetDeterminant3X3(sum);
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];
593 for(Int_t kk=0;kk<3;kk++) vett[kk][zz]=dsum[kk];
594 initPos[zz]=GetDeterminant3X3(vett)/det;
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];
604 sigma+=GetStrLinMinDist(p0,p1,initPos);
607 sigma=TMath::Sqrt(sigma);
609 Warning("StrLinVertexFinderMinDist","Finder did not succed");
614 fVert.SetXYZ(initPos);
615 fVert.SetDispersion(sigma);
616 fVert.SetNContributors(knacc);
618 //---------------------------------------------------------------------------
619 void AliVertexerTracks::TooFewTracks(const AliESD* esdEvent) {
621 // When the number of tracks is < fMinTracks
624 // deal with vertices not found
625 Double_t pos[3],err[3];
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)
642 fCurrentVertex = new AliESDVertex(pos,err);
643 fCurrentVertex->SetNContributors(ncontr);
646 esdEvent->GetVertex()->GetTruePos(tp);
647 fCurrentVertex->SetTruePos(tp);
648 fCurrentVertex->SetTitle("VertexerTracksWithConstraint");
649 if(ncontr==-1||ncontr==-2)
650 fCurrentVertex->SetTitle("VertexerTracksNoConstraint");
654 //---------------------------------------------------------------------------
655 void AliVertexerTracks::VertexFinder(Int_t optUseWeights) {
657 // Get estimate of vertex position in (x,y) from tracks DCA
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.};
670 Double_t pos[3],dir[3];
671 Double_t alpha,mindist;
672 Double_t field=GetField();
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);
682 // AliStrLine *line1 = new AliStrLine();
683 // track1->ApproximateHelixWithLine(mindist,field,line1);
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];
698 if(optUseWeights<=0){
699 Int_t retcode = line2->Cross(line1,crosspoint);
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]);
707 Int_t retcode = line1->CrossPoints(line2,pnt1,pnt2);
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];
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]);
735 for(Int_t jj=0;jj<3;jj++){
736 initPos[jj] = aver[jj]/ncombi;
738 sigmasq[jj]=aversq[jj]-initPos[jj]*initPos[jj];
741 sigma=TMath::Sqrt(TMath::Abs(sigma));
744 Warning("VertexFinder","Finder did not succed");
747 fVert.SetXYZ(initPos);
748 fVert.SetDispersion(sigma);
749 fVert.SetNContributors(ncombi);
751 //---------------------------------------------------------------------------
752 void AliVertexerTracks::VertexFitter(Bool_t useNominalVtx) {
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
759 fVert.GetXYZ(initPos);
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]);
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]));
771 rv(0,0) = initPos[0];
772 rv(1,0) = initPos[1];
774 Double_t xlStart,alpha;
776 Double_t cosRot,sinRot;
780 Int_t arrEntries = (Int_t)fTrkArray.GetEntries();
784 // initial vertex covariance matrix
786 vVb(0,0) = fNominalCov[0];
787 vVb(0,1) = fNominalCov[1];
789 vVb(1,0) = fNominalCov[1];
790 vVb(1,1) = fNominalCov[2];
794 vVb(2,2) = fNominalCov[5];
795 TMatrixD vVbInv(TMatrixD::kInverted,vVb);
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);
804 // 1st - estimate of vtx using all tracks
805 // 2nd - estimate of global chi2
806 for(step=0; step<2; step++) {
807 if(fDebug) printf(" step = %d\n",step);
811 TMatrixD sumWiri(3,1);
815 for(j=0; j<3; j++) sumWi(j,i) = 0.;
821 sumWiri(i,0) += vVbInvrb(i,0);
822 for(k=0;k<3;k++) sumWi(i,k) += vVbInv(i,k);
828 for(k=0; k<arrEntries; k++) {
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);
833 t->AliExternalTrackParam::PropagateTo(xlStart,AliTracker::GetBz()); // to vtxSeed
835 if(alpha<0.) rotAngle += 2.*TMath::Pi();
836 cosRot = TMath::Cos(rotAngle);
837 sinRot = TMath::Sin(rotAngle);
839 // vector of track global coordinates
841 ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot;
842 ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot;
845 // matrix to go from global (x,y,z) to local (y,z);
854 // covariance matrix of local (y,z) - inverted
856 t->GetExternalCovariance(cc);
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);
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);
879 TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
885 } // end loop on tracks
887 if(nUsedTrks < fMinTracks) {
892 Double_t determinant = sumWi.Determinant();
893 //cerr<<" determinant: "<<determinant<<endl;
894 if(determinant < 100.) {
895 printf("det(V) = 0\n");
900 // inverted of weights matrix
901 TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
904 // position of primary vertex
907 } // end loop on the 2 steps
912 if(fDebug) printf("TooFewTracks\n");
913 fCurrentVertex = new AliESDVertex(0.,0.,-1);
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);
929 // store data in the vertex object
930 fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nUsedTrks);
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();
940 //----------------------------------------------------------------------------
941 AliVertex* AliVertexerTracks::VertexForSelectedTracks(TTree *trkTree) {
943 // Return vertex from tracks in trkTree
946 // get tracks and propagate them to initial vertex position
947 Int_t nTrksPrep = PrepareTracks(*trkTree,1);
948 if(nTrksPrep < TMath::Max(2,fMinTracks) ) {
949 if(fDebug) printf("TooFewTracks\n");
950 Double_t vtx[3]={0,0,0};
952 fVert.SetDispersion(999);
953 fVert.SetNContributors(-5);
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);
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();
975 fVert.SetIndices(fVert.GetNContributors(),indices);
982 //----------------------------------------------------------------------------
983 AliVertex* AliVertexerTracks::VertexForSelectedTracks(TObjArray *trkArray) {
985 // Return vertex from array of tracks
988 // get tracks and propagate them to initial vertex position
989 Int_t nTrks = trkArray->GetEntriesFast();
990 if(nTrks < TMath::Max(2,fMinTracks) ) {
991 if(fDebug) printf("TooFewTracks\n");
992 Double_t vtx[3]={0,0,0};
994 fVert.SetDispersion(999);
995 fVert.SetNContributors(-5);
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);
1007 AliVertex *vtx = VertexForSelectedTracks(trkTree);
1011 //--------------------------------------------------------------------------