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