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                                            Double_t field,
  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   SetField(field);
  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(Double_t field, TString fn,
  81                                            Double_t xStart,Double_t yStart)
  82                                           :AliITSVertexer(fn) {
  83 //
  84 // Alternative constructor
  85 //
  86   fInFile  = 0;
  87   fOutFile = 0;
  88   SetField(field);
  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     itstrack = new AliITStrackV2(*esdTrack);
 351     trkTree->Fill();
 352     itstrack = 0;
 353     delete esdTrack;
 354   }
 355   delete itstrack;
 356
 357   // preselect tracks and propagate them to initial vertex position
 358   Int_t nTrks = PrepareTracks(*trkTree);
 359   delete trkTree;
 360   if(fDebug) printf(" tracks prepared: %d\n",nTrks);
 361   if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; }
 362
 363   // Set initial vertex position from ESD
 364   esdEvent->GetVertex()->GetXYZ(vtx);
 365   SetVtxStart(vtx[0],vtx[1]);
 366
 367   // VERTEX FINDER
 368   VertexFinder();
 369
 370   // VERTEX FITTER
 371   ComputeMaxChi2PerTrack(nTrks);
 372   VertexFitter();
 373   if(fDebug) printf(" vertex fit completed\n");
 374
 375   // store vertex information in ESD
 376   fCurrentVertex->GetXYZ(vtx);
 377   fCurrentVertex->GetCovMatrix(cvtx);
 378   esdEvent->SetVertex(fCurrentVertex);
 379
 380   cout<<"Vertex: "<<vtx[0]<<", "<<vtx[1]<<", "<<vtx[2]<<endl;
 381   return fCurrentVertex;
 382 }
 383 //---------------------------------------------------------------------------
 384 void AliITSVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped) {
 385 //
 386 // Mark the tracks not ot be used in the vertex finding
 387 //
 388   fNTrksToSkip = n;
 389   fTrksToSkip = new Int_t[n];
 390   for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i];
 391   return;
 392 }
 393 //---------------------------------------------------------------------------
 394 void AliITSVertexerTracks::TooFewTracks() {
 395 //
 396 // When the number of tracks is < fMinTracks the vertex is set to (0,0,0)
 397 // and the number of tracks to -1
 398 //
 399   fCurrentVertex = new AliESDVertex(0.,0.,-1);
 400   return;
 401 }
 402 //---------------------------------------------------------------------------
 403 void AliITSVertexerTracks::VertexFinder() {
 404
 405   // Get estimate of vertex position in (x,y) from tracks DCA
 406   // Then this estimate is stored to the data member fInitPos
 407   // (previous values are overwritten)
 408
 409
 410
 411   /*
 412 ******* TEMPORARY!!! FOR TEST ONLY!!! **********************************
 413
 414 fInitPos[0] = fNominalPos[0]+gRandom->Gaus(0.,0.0100); // 100 micron gaussian smearing
 415 fInitPos[1] = fNominalPos[1]+gRandom->Gaus(0.,0.0100); // 100 micron gaussian smearing
 416   */
 417
 418   fInitPos[2] = 0.;
 419   for(Int_t i=0;i<2;i++)fInitPos[i]=fNominalPos[i];
 420
 421   Int_t nacc = (Int_t)fTrkArray.GetEntriesFast();
 422
 423   Double_t aver[3]={0.,0.,0.};
 424   Int_t ncombi = 0;
 425   AliITStrackV2 *track1;
 426   AliITStrackV2 *track2;
 427   for(Int_t i=0; i<nacc; i++){
 428     track1 = (AliITStrackV2*)fTrkArray.At(i);
 429     if(fDebug>5){
 430       Double_t xv,par[5];
 431       track1->GetExternalParameters(xv,par);
 432       cout<<"Track in position "<<i<<" xr= "<<xv<<endl;
 433       for(Int_t ii=0;ii<5;ii++)cout<<par[ii]<<" ";
 434       cout<<endl;
 435     }
 436     Double_t mom1[3];
 437     Double_t alpha = track1->GetAlpha();
 438     Double_t azim = TMath::ASin(track1->GetSnp())+alpha;
 439     Double_t theta = TMath::Pi()/2. - TMath::ATan(track1->GetTgl());
 440     mom1[0] = TMath::Sin(theta)*TMath::Cos(azim);
 441     mom1[1] = TMath::Sin(theta)*TMath::Sin(azim);
 442     mom1[2] = TMath::Cos(theta);
 443
 444     Double_t pos1[3];
 445     Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
 446     track1->GetGlobalXYZat(mindist,pos1[0],pos1[1],pos1[2]);
 447     AliITSStrLine *line1 = new AliITSStrLine(pos1,mom1);
 448     for(Int_t j=i+1; j<nacc; j++){
 449       track2 = (AliITStrackV2*)fTrkArray.At(j);
 450       Double_t mom2[3];
 451       alpha = track2->GetAlpha();
 452       azim = TMath::ASin(track2->GetSnp())+alpha;
 453       theta = TMath::Pi()/2. - TMath::ATan(track2->GetTgl());
 454       mom2[0] = TMath::Sin(theta)*TMath::Cos(azim);
 455       mom2[1] = TMath::Sin(theta)*TMath::Sin(azim);
 456       mom2[2] = TMath::Cos(theta);
 457       Double_t pos2[3];
 458       mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1];
 459       track2->GetGlobalXYZat(mindist,pos2[0],pos2[1],pos2[2]);
 460       AliITSStrLine *line2 = new AliITSStrLine(pos2,mom2);
 461       Double_t crosspoint[3];
 462       Int_t retcode = line2->Cross(line1,crosspoint);
 463       if(retcode<0){
 464         if(fDebug>10)cout<<" i= "<<i<<",   j= "<<j<<endl;
 465         if(fDebug>10)cout<<"bad intersection\n";
 466         line1->PrintStatus();
 467         line2->PrintStatus();
 468       }
 469       else {
 470         ncombi++;
 471         for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj];
 472         if(fDebug>10)cout<<" i= "<<i<<",   j= "<<j<<endl;
 473         if(fDebug>10)cout<<"\n Cross point: ";
 474         if(fDebug>10)cout<<crosspoint[0]<<" "<<crosspoint[1]<<" "<<crosspoint[2]<<endl;
 475       }
 476       delete line2;
 477     }
 478     delete line1;
 479   }
 480   if(ncombi>0){
 481     for(Int_t jj=0;jj<3;jj++)fInitPos[jj] = aver[jj]/ncombi;
 482   }
 483   else {
 484     Warning("VertexFinder","Finder did not succed");
 485   }
 486
 487
 488   //************************************************************************
 489   return;
 490
 491 }
 492 //---------------------------------------------------------------------------
 493 void AliITSVertexerTracks::VertexFitter() {
 494 //
 495 // The optimal estimate of the vertex position is given by a "weighted
 496 // average of tracks positions"
 497 // Original method: CMS Note 97/0051
 498 //
 499   if(fDebug) {
 500     printf(" VertexFitter(): start\n");
 501     PrintStatus();
 502   }
 503
 504
 505   Int_t i,j,k,step=0;
 506   TMatrixD rv(3,1);
 507   TMatrixD vV(3,3);
 508   rv(0,0) = fInitPos[0];
 509   rv(1,0) = fInitPos[1];
 510   rv(2,0) = 0.;
 511   Double_t xlStart,alpha;
 512   Double_t rotAngle;
 513   Double_t cosRot,sinRot;
 514   Double_t cc[15];
 515   Int_t nUsedTrks;
 516   Double_t chi2,chi2i;
 517   Int_t arrEntries = (Int_t)fTrkArray.GetEntries();
 518   AliITStrackV2 *t = 0;
 519   Int_t failed = 0;
 520
 521   Int_t *skipTrack = new Int_t[arrEntries];
 522   for(i=0; i<arrEntries; i++) skipTrack[i]=0;
 523
 524   // 3 steps:
 525   // 1st - first estimate of vtx using all tracks
 526   // 2nd - apply cut on chi2 max per track
 527   // 3rd - estimate of global chi2
 528   for(step=0; step<3; step++) {
 529     if(fDebug) printf(" step = %d\n",step);
 530     chi2 = 0.;
 531     nUsedTrks = 0;
 532
 533     TMatrixD sumWiri(3,1);
 534     TMatrixD sumWi(3,3);
 535     for(i=0; i<3; i++) {
 536       sumWiri(i,0) = 0.;
 537       for(j=0; j<3; j++) sumWi(j,i) = 0.;
 538     }
 539
 540     // loop on tracks
 541     for(k=0; k<arrEntries; k++) {
 542       if(skipTrack[k]) continue;
 543       // get track from track array
 544       t = (AliITStrackV2*)fTrkArray.At(k);
 545       alpha = t->GetAlpha();
 546       xlStart = fInitPos[0]*TMath::Cos(alpha)+fInitPos[1]*TMath::Sin(alpha);
 547       t->PropagateTo(xlStart,0.,0.);   // to vtxSeed
 548       rotAngle = alpha;
 549       if(alpha<0.) rotAngle += 2.*TMath::Pi();
 550       cosRot = TMath::Cos(rotAngle);
 551       sinRot = TMath::Sin(rotAngle);
 552
 553       // vector of track global coordinates
 554       TMatrixD ri(3,1);
 555       ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot;
 556       ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot;
 557       ri(2,0) = t->GetZ();
 558
 559       // matrix to go from global (x,y,z) to local (y,z);
 560       TMatrixD qQi(2,3);
 561       qQi(0,0) = -sinRot;
 562       qQi(0,1) = cosRot;
 563       qQi(0,2) = 0.;
 564       qQi(1,0) = 0.;
 565       qQi(1,1) = 0.;
 566       qQi(1,2) = 1.;
 567
 568       // covariance matrix of local (y,z) - inverted
 569       TMatrixD uUi(2,2);
 570       t->GetExternalCovariance(cc);
 571       uUi(0,0) = cc[0];
 572       uUi(0,1) = cc[1];
 573       uUi(1,0) = cc[1];
 574       uUi(1,1) = cc[2];
 575
 576       // weights matrix: wWi = qQiT * uUiInv * qQi
 577       if(uUi.Determinant() <= 0.) continue;
 578       TMatrixD uUiInv(TMatrixD::kInverted,uUi);
 579       TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi);
 580       TMatrixD wWi(qQi,TMatrixD::kTransposeMult,uUiInvQi);
 581
 582       // track chi2
 583       TMatrixD deltar = rv; deltar -= ri;
 584       TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar);
 585       chi2i = deltar(0,0)*wWideltar(0,0)+
 586               deltar(1,0)*wWideltar(1,0)+
 587               deltar(2,0)*wWideltar(2,0);
 588
 589
 590       if(step==1 && chi2i > fMaxChi2PerTrack) {
 591         skipTrack[k] = 1;
 592         continue;
 593       }
 594
 595       // add to total chi2
 596       chi2 += chi2i;
 597
 598       TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
 599
 600       sumWiri += wWiri;
 601       sumWi   += wWi;
 602
 603       nUsedTrks++;
 604     } // end loop on tracks
 605
 606     if(nUsedTrks < fMinTracks) {
 607       failed=1;
 608       continue;
 609     }
 610
 611     Double_t determinant = sumWi.Determinant();
 612     //cerr<<" determinant: "<<determinant<<endl;
 613     if(determinant < 100.)  {
 614       printf("det(V) = 0\n");
 615       failed=1;
 616       continue;
 617     }
 618
 619     // inverted of weights matrix
 620     TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
 621     vV = invsumWi;
 622
 623     // position of primary vertex
 624     rv.Mult(vV,sumWiri);
 625
 626   } // end loop on the 3 steps
 627
 628   delete [] skipTrack;
 629   delete t;
 630
 631   if(failed) {
 632     TooFewTracks();
 633     return;
 634   }
 635
 636   Double_t position[3];
 637   position[0] = rv(0,0);
 638   position[1] = rv(1,0);
 639   position[2] = rv(2,0);
 640   Double_t covmatrix[6];
 641   covmatrix[0] = vV(0,0);
 642   covmatrix[1] = vV(0,1);
 643   covmatrix[2] = vV(1,1);
 644   covmatrix[3] = vV(0,2);
 645   covmatrix[4] = vV(1,2);
 646   covmatrix[5] = vV(2,2);
 647
 648   // store data in the vertex object
 649   fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nUsedTrks);
 650
 651   if(fDebug) {
 652     printf(" VertexFitter(): finish\n");
 653     printf(" rv = ( %f , %f , %f )\n\n",rv(0,0),rv(1,0),rv(2,0));
 654     fCurrentVertex->PrintStatus();
 655   }
 656
 657   return;
 658 }
 659 //----------------------------------------------------------------------------
 660 AliESDVertex *AliITSVertexerTracks::VertexOnTheFly(TTree &trkTree) {
 661 //
 662 // Return vertex from tracks in trkTree
 663 //
 664   if(fCurrentVertex) fCurrentVertex = 0;
 665
 666   // get tracks and propagate them to initial vertex position
 667   Int_t nTrks = PrepareTracks(*(&trkTree));
 668   if(fDebug) printf(" tracks prepared: %d\n",nTrks);
 669   if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; }
 670
 671   // VERTEX FINDER
 672   VertexFinder();
 673
 674   // VERTEX FITTER
 675   ComputeMaxChi2PerTrack(nTrks);
 676   VertexFitter();
 677   if(fDebug) printf(" vertex fit completed\n");
 678
 679   return fCurrentVertex;
 680 }
 681 //----------------------------------------------------------------------------
 682
 683
 684