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