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