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