]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AliITSVertexerTracks.cxx
Transient pointer to AliRunDigitizer
[u/mrichter/AliRoot.git] / ITS / AliITSVertexerTracks.cxx
CommitLineData
cab6ff9b 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 tracks
18//
19// Origin: A.Dainese, Padova, andrea.dainese@pd.infn.it
20// M.Masera, Torino, massimo.masera@to.infn.it
21//-----------------------------------------------------------------
22
23//---- standard headers ----
24#include <Riostream.h>
25//---- Root headers --------
26#include <TFile.h>
27#include <TTree.h>
28#include <TVector3.h>
29#include <TMatrixD.h>
cab6ff9b 30#include <TRandom.h> // TEMPORARY !!!!!!!
31//---- AliRoot headers -----
32#include <AliRun.h>
33#include "AliKalmanTrack.h"
34#include "AliITSStrLine.h"
35#include "AliITStrackV2.h"
36#include "AliITSVertex.h"
37#include "AliITSVertexerTracks.h"
38
39
40ClassImp(AliITSVertexerTracks)
41
42
43//----------------------------------------------------------------------------
44AliITSVertexerTracks::AliITSVertexerTracks():AliITSVertexer() {
45//
46// Default constructor
47//
48
49 SetVtxStart();
50 SetMinTracks();
51 SetUseThrustFrame();
52 SetPhiThrust();
53 fTrksToSkip = 0;
54 fNTrksToSkip = 0;
55 for(Int_t i=0; i<3; i++)fInitPos[i] = 0.;
56}
57
58//----------------------------------------------------------------------------
88cb7938 59AliITSVertexerTracks::AliITSVertexerTracks(Double_t field, TString fn,
cab6ff9b 60 Double_t xStart,Double_t yStart,
61 Int_t useThFr)
88cb7938 62 :AliITSVertexer(fn) {
cab6ff9b 63//
64// Standard constructor
65//
66
67 SetField(field);
68 SetVtxStart(xStart,yStart);
69 SetMinTracks();
70 SetUseThrustFrame(useThFr);
71 SetPhiThrust();
72 fTrksToSkip = 0;
73 fNTrksToSkip = 0;
74 for(Int_t i=0; i<3; i++)fInitPos[i] = 0.;
75}
76//----------------------------------------------------------------------------
77Bool_t AliITSVertexerTracks::CheckField() const {
78//
79// Check if the conv. const. has been set
80//
81 AliITStrackV2 t;
82 Double_t cc = t.GetConvConst();
83 Double_t field = 100./0.299792458/cc;
84
85 if(field<0.1 || field>0.6) {
86 printf("AliITSVertexerTracks::CheckField():\n ERROR: AliKalmanTrack::fConvConst not set\n Use AliKalmanTrack::SetConvConst() or AliITSVertexerTracks::SetField()\n");
87 return kFALSE;
88 }
89 printf("AliITSVertexerTracks::CheckField(): Using B = %3.1f T\n",field);
90 return kTRUE;
91}
92//---------------------------------------------------------------------------
93void AliITSVertexerTracks::ComputeMaxChi2PerTrack(Int_t nTracks) {
94//
95// Max. contr. to the chi2 has been tuned as a function of multiplicity
96//
97 if(nTracks < 7) { fMaxChi2PerTrack = 1.e6;
98 } else { fMaxChi2PerTrack = 100.; }
99
100 return;
101}
102//---------------------------------------------------------------------------
103void AliITSVertexerTracks::FindVertices() {
104//
105// Vertices for all events from fFirstEvent to fLastEvent
106//
107
108 // Check if the conv. const. has been set
109 if(!CheckField()) return;
110
111
112 // loop over events
113 for(Int_t ev=fFirstEvent; ev<=fLastEvent; ev++) {
114 if(ev % 100 == 0 || fDebug) printf("--- Processing event %d of %d ---\n",ev,fLastEvent);
115
116 FindVertexForCurrentEvent(ev);
117
118 if(!fCurrentVertex) {
119 printf("AliITSVertexerTracks::FindVertixes(): no tracks tree for event %d\n",ev);
120 continue;
121 }
122
123 if(fDebug) fCurrentVertex->PrintStatus();
d4221964 124 TString vtxName = "Vertex_";
125 vtxName += ev;
88cb7938 126 // fCurrentVertex->SetName(vtxName.Data());
d4221964 127 fCurrentVertex->SetTitle("VertexerTracks");
cab6ff9b 128 WriteCurrentVertex();
129 } // loop over events
130
131 return;
132}
133//----------------------------------------------------------------------------
134Int_t AliITSVertexerTracks::PrepareTracks(TTree &trkTree) {
135//
136// Propagate tracks to initial vertex position and store them in a TObjArray
137//
138 Double_t maxd0rphi = 3.;
139 Double_t alpha,xlStart,d0rphi;
140 Int_t nTrks = 0;
141 Bool_t skipThis;
142
143 Int_t nEntries = (Int_t)trkTree.GetEntries();
144
145 if(!fTrkArray.IsEmpty()) fTrkArray.Clear();
146 fTrkArray.Expand(nEntries);
147
148 if(fDebug) {
149 printf(" PrepareTracks()\n");
150 trkTree.Print();
151 }
152
153 for(Int_t i=0; i<nEntries; i++) {
154 // check tracks to skip
155 skipThis = kFALSE;
156 for(Int_t j=0; j<fNTrksToSkip; j++) {
157 if(i==fTrksToSkip[j]) {
158 if(fDebug) printf("skipping track: %d\n",i);
159 skipThis = kTRUE;
160 }
161 }
162 if(skipThis) continue;
163
164 AliITStrackV2 *itstrack = new AliITStrackV2;
165 trkTree.SetBranchAddress("tracks",&itstrack);
166 trkTree.GetEvent(i);
167
168
169 // propagate track to vtxSeed
170 alpha = itstrack->GetAlpha();
171 xlStart = fNominalPos[0]*TMath::Cos(alpha)+fNominalPos[1]*TMath::Sin(alpha);
172 itstrack->PropagateTo(3.,0.0023,65.19); // to beam pipe (0.8 mm of Be)
173 itstrack->PropagateTo(xlStart,0.,0.); // to vtxSeed
174
175 // select tracks with d0rphi < maxd0rphi
176 d0rphi = TMath::Abs(itstrack->GetD(fNominalPos[0],fNominalPos[1]));
177 if(d0rphi > maxd0rphi) { delete itstrack; continue; }
178
179 fTrkArray.AddLast(itstrack);
180
181 nTrks++;
182 }
183
184 if(fTrksToSkip) delete [] fTrksToSkip;
185
186 return nTrks;
187}
188//----------------------------------------------------------------------------
189void AliITSVertexerTracks::PrintStatus() const {
190//
191// Print status
192//
193 printf(" Initial position (%f,%f)\n",fNominalPos[0],fNominalPos[1]);
194 printf(" Vertex position after vertex finder (%f, %f, %f)\n",fInitPos[0],fInitPos[1],fInitPos[2]);
195 printf(" Number of tracks in array: %d\n",(Int_t)fTrkArray.GetEntriesFast());
196 printf(" Minimum # tracks required in fit: %d\n",fMinTracks);
197 printf(" Using Thrust Frame: %d fPhiThrust = %f\n",fUseThrustFrame,fPhiThrust);
198
199 return;
200}
201//----------------------------------------------------------------------------
202AliITSVertex* AliITSVertexerTracks::FindVertexForCurrentEvent(Int_t evnumb) {
203//
204// Vertex for current event
205//
206 fCurrentVertex = 0;
207
208 // get tree with tracks from input file
209 TString treeName = "TreeT_ITS_";
210 treeName += evnumb;
88cb7938 211 // TTree *trkTree=(TTree*)fInFile->Get(treeName.Data()); masera
088e0b8d 212 TTree *trkTree=0;
cab6ff9b 213 if(!trkTree) return fCurrentVertex;
214
215
216 // get tracks and propagate them to initial vertex position
217 Int_t nTrks = PrepareTracks(*trkTree);
218 delete trkTree;
219 if(fDebug) printf(" tracks prepared: %d\n",nTrks);
220 if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; }
221
222 // VERTEX FINDER
223 VertexFinder();
224
225 // VERTEX FITTER
226 ComputeMaxChi2PerTrack(nTrks);
227 if(fUseThrustFrame) ThrustFinderXY();
228 if(fDebug) printf(" thrust found: phi = %f\n",fPhiThrust);
229 VertexFitter();
230 if(fDebug) printf(" vertex fit completed\n");
231
232 TString vtxName;
d4221964 233 vtxName = "Vertex_";
cab6ff9b 234 vtxName += evnumb;
88cb7938 235 // fCurrentVertex->SetName(vtxName.Data());
cab6ff9b 236 return fCurrentVertex;
237}
238//---------------------------------------------------------------------------
239void AliITSVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped) {
240//
241// Mark the tracks not ot be used in the vertex finding
242//
243 fNTrksToSkip = n;
244 fTrksToSkip = new Int_t[n];
245 for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i];
246 return;
247}
248//----------------------------------------------------------------------------
249Double_t AliITSVertexerTracks::SumPl(TTree &momTree,Double_t phi) const {
250//
251// Function to be maximized for thrust determination
252//
253 TVector3 u(1.,1.,0);
254 u.SetMag(1.);
255 u.SetPhi(phi);
256 Double_t pl;
257
258 Double_t sum = 0.;
259
260 TVector3* mom=0;
261 momTree.SetBranchAddress("momenta",&mom);
262 Int_t entries = (Int_t)momTree.GetEntries();
263
264 for(Int_t i=0; i<entries; i++) {
265 momTree.GetEvent(i);
266 pl = mom->Dot(u);
267 if(pl>0.) sum += pl;
268 }
269
270 delete mom;
271
272 return sum;
273}
274//---------------------------------------------------------------------------
275void AliITSVertexerTracks::ThrustFinderXY() {
276//
277// This function looks for the thrust direction, \vec{u}, in the (x,y) plane.
278// The following function is maximized:
279// \Sum_{\vec{p}\cdot\vec{u}} \vec{p}\cdot\vec{u} / \Sum |\vec{p}|
280// where \vec{p} = (p_x,p_y)
281//
282 Double_t pt,alpha,phi;
283 Double_t totPt;
284
285 // tree for thrust determination
286 TVector3 *ioMom = new TVector3;
287 TTree *t = new TTree("Tree_Momenta","Tree with momenta");
288 t->Branch("momenta","TVector3",&ioMom);
289 totPt = 0.;
290
291 AliITStrackV2 *itstrack = 0;
292 Int_t arrEntries = (Int_t)fTrkArray.GetEntries();
293
294 // loop on tracks
295 for(Int_t i=0; i<arrEntries; i++) {
296 itstrack = (AliITStrackV2*)fTrkArray.At(i);
297 // momentum of the track at the vertex
298 pt = 1./TMath::Abs(itstrack->Get1Pt());
299 alpha = itstrack->GetAlpha();
300 phi = alpha+TMath::ASin(itstrack->GetSnp());
301 ioMom->SetX(pt*TMath::Cos(phi));
302 ioMom->SetY(pt*TMath::Sin(phi));
303 ioMom->SetZ(0.);
304
305 totPt += ioMom->Pt();
306 t->Fill();
307 } // end loop on tracks
308
309 Double_t tValue=0.,tPhi=0.;
310 Double_t maxSumPl = 0.;
311 Double_t thisSumPl;
312 Double_t dPhi;
313 Int_t nSteps,iStep;
314
315 phi = 0.;
316 nSteps = 100;
317 dPhi = 2.*TMath::Pi()/(Double_t)nSteps;
318
319 for(iStep=0; iStep<nSteps; iStep++) {
320 phi += dPhi;
321 thisSumPl = SumPl(*t,phi);
322 if(thisSumPl > maxSumPl) {
323 maxSumPl = thisSumPl;
324 tPhi = phi;
325 }
326 }
327
328 phi = tPhi-dPhi;
329 nSteps = 10;
330 dPhi /= 5.;
331
332 for(iStep=0; iStep<nSteps; iStep++) {
333 phi += dPhi;
334 thisSumPl = SumPl(*t,phi);
335 if(thisSumPl > maxSumPl) {
336 maxSumPl = thisSumPl;
337 tPhi = phi;
338 }
339 }
340
341 tValue = 2.*maxSumPl/totPt;
342 if(tPhi<0.) tPhi += 2.*TMath::Pi();
343 if(tPhi>2.*TMath::Pi()) tPhi -= 2.*TMath::Pi();
344
345 SetPhiThrust(tPhi);
346
347 delete t;
348 delete ioMom;
349
350 return;
351}
352//---------------------------------------------------------------------------
353void AliITSVertexerTracks::TooFewTracks() {
354//
355// When the number of tracks is < fMinTracks the vertex is set to (0,0,0)
356// and the number of tracks to -1
357//
358 fCurrentVertex = new AliITSVertex(0.,0.,-1);
359 return;
360}
361//---------------------------------------------------------------------------
362void AliITSVertexerTracks::VertexFinder() {
363
364 // Get estimate of vertex position in (x,y) from tracks DCA
365 // Then this estimate is stored to the data member fInitPos
366 // (previous values are overwritten)
367
368
369
370 /*
371******* TEMPORARY!!! FOR TEST ONLY!!! **********************************
372
373fInitPos[0] = fNominalPos[0]+gRandom->Gaus(0.,0.0100); // 100 micron gaussian smearing
374fInitPos[1] = fNominalPos[1]+gRandom->Gaus(0.,0.0100); // 100 micron gaussian smearing
375 */
376
377 fInitPos[2] = 0.;
378 for(Int_t i=0;i<2;i++)fInitPos[i]=fNominalPos[i];
379
380 Int_t nacc = (Int_t)fTrkArray.GetEntriesFast();
381
382 Double_t aver[3]={0.,0.,0.};
383 Int_t ncombi = 0;
384 AliITStrackV2 *track1;
385 AliITStrackV2 *track2;
386 for(Int_t i=0; i<nacc; i++){
387 track1 = (AliITStrackV2*)fTrkArray.At(i);
388 if(fDebug>5){
389 Double_t xv,par[5];
390 track1->GetExternalParameters(xv,par);
391 cout<<"Track in position "<<i<<" xr= "<<xv<<endl;
392 for(Int_t ii=0;ii<5;ii++)cout<<par[ii]<<" ";
393 cout<<endl;
394 }
395 Double_t mom1[3];
396 Double_t alpha = track1->GetAlpha();
397 Double_t azim = TMath::ASin(track1->GetSnp())+alpha;
398 Double_t theta = TMath::Pi()/2. - TMath::ATan(track1->GetTgl());
399 mom1[0] = TMath::Sin(theta)*TMath::Cos(azim);
400 mom1[1] = TMath::Sin(theta)*TMath::Sin(azim);
401 mom1[2] = TMath::Cos(theta);
402
403 Double_t pos1[3];
404 Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
405 track1->GetGlobalXYZat(mindist,pos1[0],pos1[1],pos1[2]);
406 AliITSStrLine *line1 = new AliITSStrLine(pos1,mom1);
407 for(Int_t j=i+1; j<nacc; j++){
408 track2 = (AliITStrackV2*)fTrkArray.At(j);
409 Double_t mom2[3];
410 alpha = track2->GetAlpha();
411 azim = TMath::ASin(track2->GetSnp())+alpha;
412 theta = TMath::Pi()/2. - TMath::ATan(track2->GetTgl());
413 mom2[0] = TMath::Sin(theta)*TMath::Cos(azim);
414 mom2[1] = TMath::Sin(theta)*TMath::Sin(azim);
415 mom2[2] = TMath::Cos(theta);
416 Double_t pos2[3];
417 mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
418 track2->GetGlobalXYZat(mindist,pos2[0],pos2[1],pos2[2]);
419 AliITSStrLine *line2 = new AliITSStrLine(pos2,mom2);
420 Double_t crosspoint[3];
421 Int_t retcode = line2->Cross(line1,crosspoint);
422 if(retcode<0){
423 if(fDebug>10)cout<<" i= "<<i<<", j= "<<j<<endl;
424 if(fDebug>10)cout<<"bad intersection\n";
425 line1->PrintStatus();
426 line2->PrintStatus();
427 }
428 else {
429 ncombi++;
430 for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
431 if(fDebug>10)cout<<" i= "<<i<<", j= "<<j<<endl;
432 if(fDebug>10)cout<<"\n Cross point: ";
433 if(fDebug>10)cout<<crosspoint[0]<<" "<<crosspoint[1]<<" "<<crosspoint[2]<<endl;
434 }
435 delete line2;
436 }
437 delete line1;
438 }
439 if(ncombi>0){
440 for(Int_t jj=0;jj<3;jj++)fInitPos[jj] = aver[jj]/ncombi;
441 }
442 else {
443 Warning("VertexFinder","Finder did not succed");
444 }
445
446
447 //************************************************************************
448 return;
449
450}
451
452//---------------------------------------------------------------------------
453void AliITSVertexerTracks::VertexFitter() {
454//
455// The optimal estimate of the vertex position is given by a "weighted
456// average of tracks positions"
457// Original method: CMS Note 97/0051
458//
459 if(fDebug) {
460 printf(" VertexFitter(): start\n");
461 PrintStatus();
462 }
463
464
465 Int_t i,j,k,step=0;
466 TMatrixD rv(3,1);
467 TMatrixD V(3,3);
468 rv(0,0) = fInitPos[0];
469 rv(1,0) = fInitPos[1];
470 rv(2,0) = 0.;
471 Double_t xlStart,alpha;
472 Double_t rotAngle;
473 Double_t cosRot,sinRot;
474 Double_t cc[15];
475 Int_t nUsedTrks;
476 Double_t chi2,chi2i;
477 Int_t arrEntries = (Int_t)fTrkArray.GetEntries();
478 AliITStrackV2 *t = 0;
479 Int_t failed = 0;
480
481 Int_t *skipTrack = new Int_t[arrEntries];
482 for(i=0; i<arrEntries; i++) skipTrack[i]=0;
483
484 // 3 steps:
485 // 1st - first estimate of vtx using all tracks
486 // 2nd - apply cut on chi2 max per track
487 // 3rd - estimate of global chi2
488 for(step=0; step<3; step++) {
489 if(fDebug) printf(" step = %d\n",step);
490 chi2 = 0.;
491 nUsedTrks = 0;
492
493 TMatrixD SumWiri(3,1);
494 TMatrixD SumWi(3,3);
495 for(i=0; i<3; i++) {
496 SumWiri(i,0) = 0.;
497 for(j=0; j<3; j++) SumWi(j,i) = 0.;
498 }
499
500 // loop on tracks
501 for(k=0; k<arrEntries; k++) {
502 if(skipTrack[k]) continue;
503 // get track from track array
504 t = (AliITStrackV2*)fTrkArray.At(k);
505 alpha = t->GetAlpha();
506 xlStart = fInitPos[0]*TMath::Cos(alpha)+fInitPos[1]*TMath::Sin(alpha);
507 t->PropagateTo(xlStart,0.,0.); // to vtxSeed
508 rotAngle = alpha-fPhiThrust;
509 if(alpha<0.) rotAngle += 2.*TMath::Pi();
510 cosRot = TMath::Cos(rotAngle);
511 sinRot = TMath::Sin(rotAngle);
512
513 // vector of track global coordinates
514 TMatrixD ri(3,1);
515 ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot;
516 ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot;
517 ri(2,0) = t->GetZ();
518
519 // matrix to go from global (x,y,z) to local (y,z);
520 TMatrixD Qi(2,3);
521 Qi(0,0) = -sinRot;
522 Qi(0,1) = cosRot;
523 Qi(0,2) = 0.;
524 Qi(1,0) = 0.;
525 Qi(1,1) = 0.;
526 Qi(1,2) = 1.;
527
528 // covariance matrix of local (y,z) - inverted
529 TMatrixD Ui(2,2);
530 t->GetExternalCovariance(cc);
531 Ui(0,0) = cc[0];
532 Ui(0,1) = cc[1];
533 Ui(1,0) = cc[1];
534 Ui(1,1) = cc[2];
535
536 // weights matrix: Wi = QiT * UiInv * Qi
537 if(Ui.Determinant() <= 0.) continue;
538 TMatrixD UiInv(TMatrixD::kInverted,Ui);
539 TMatrixD UiInvQi(UiInv,TMatrixD::kMult,Qi);
540 TMatrixD Wi(Qi,TMatrixD::kTransposeMult,UiInvQi);
541
542 // track chi2
543 TMatrixD deltar = rv; deltar -= ri;
544 TMatrixD Wideltar(Wi,TMatrixD::kMult,deltar);
545 chi2i = deltar(0,0)*Wideltar(0,0)+
546 deltar(1,0)*Wideltar(1,0)+
547 deltar(2,0)*Wideltar(2,0);
548
549
550 if(step==1 && chi2i > fMaxChi2PerTrack) {
551 skipTrack[k] = 1;
552 continue;
553 }
554
555 // add to total chi2
556 chi2 += chi2i;
557
558 TMatrixD Wiri(Wi,TMatrixD::kMult,ri);
559
560 SumWiri += Wiri;
561 SumWi += Wi;
562
563 nUsedTrks++;
564 } // end loop on tracks
565
566 if(nUsedTrks < fMinTracks) {
567 failed=1;
568 continue;
569 }
570
571 Double_t determinant = SumWi.Determinant();
572 //cerr<<" determinant: "<<determinant<<endl;
573 if(determinant < 100.) {
574 printf("det(V) = 0\n");
575 failed=1;
576 continue;
577 }
578
579 // inverted of weights matrix
580 TMatrixD InvSumWi(TMatrixD::kInverted,SumWi);
581 V = InvSumWi;
582
583 // position of primary vertex
584 rv.Mult(V,SumWiri);
585
586 } // end loop on the 3 steps
587
588 delete [] skipTrack;
589 delete t;
590
591 if(failed) {
592 TooFewTracks();
593 return;
594 }
595
596 Double_t position[3];
597 position[0] = rv(0,0);
598 position[1] = rv(1,0);
599 position[2] = rv(2,0);
600 Double_t covmatrix[6];
601 covmatrix[0] = V(0,0);
602 covmatrix[1] = V(0,1);
603 covmatrix[2] = V(1,1);
604 covmatrix[3] = V(0,2);
605 covmatrix[4] = V(1,2);
606 covmatrix[5] = V(2,2);
607
608 // store data in the vertex object
609 fCurrentVertex = new AliITSVertex(fPhiThrust,position,covmatrix,chi2,nUsedTrks);
610
611 if(fDebug) {
612 printf(" VertexFitter(): finish\n");
613 printf(" rv = ( %f , %f , %f )\n\n",rv(0,0),rv(1,0),rv(2,0));
614 fCurrentVertex->PrintStatus();
615 }
616
617 return;
618}
619//----------------------------------------------------------------------------
620AliITSVertex *AliITSVertexerTracks::VertexOnTheFly(TTree &trkTree) {
621//
622// Return vertex from tracks in trkTree
623//
624 if(fCurrentVertex) fCurrentVertex = 0;
625
626 // get tracks and propagate them to initial vertex position
627 Int_t nTrks = PrepareTracks(*(&trkTree));
628 if(fDebug) printf(" tracks prepared: %d\n",nTrks);
629 if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; }
630
631 // VERTEX FINDER
632 VertexFinder();
633
634 // VERTEX FITTER
635 ComputeMaxChi2PerTrack(nTrks);
636 if(fUseThrustFrame) ThrustFinderXY();
637 if(fDebug) printf(" thrust found: phi = %f\n",fPhiThrust);
638 VertexFitter();
639 if(fDebug) printf(" vertex fit completed\n");
640
641 return fCurrentVertex;
642}
643//----------------------------------------------------------------------------
644
645
646