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