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     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