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