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()->GetXYZ(vtx);
 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(fCurrentVertex);
 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