]> git.uio.no Git - u/mrichter/AliRoot.git/blame - ITS/AliITSVertexerTracks.cxx
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Incrementing ClassDefs (Y.Schutz)
[u/mrichter/AliRoot.git] / ITS / AliITSVertexerTracks.cxx
CommitLineData
cab6ff9b 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
a48d1ed3 18// It accepts V2 and ESD tracks
cab6ff9b 19//
a48d1ed3 20// Origin: A.Dainese, Padova,
21// andrea.dainese@pd.infn.it
22// M.Masera, Torino,
23// massimo.masera@to.infn.it
cab6ff9b 24//-----------------------------------------------------------------
25
26//---- standard headers ----
27#include <Riostream.h>
28//---- Root headers --------
11ba84a4 29#include <TKey.h>
cab6ff9b 30#include <TFile.h>
31#include <TTree.h>
cab6ff9b 32#include <TMatrixD.h>
cab6ff9b 33//---- AliRoot headers -----
cab6ff9b 34#include "AliITSStrLine.h"
35#include "AliITStrackV2.h"
d681bb2d 36#include "AliESDVertex.h"
cab6ff9b 37#include "AliITSVertexerTracks.h"
11ba84a4 38#include "AliESD.h"
39#include "AliESDtrack.h"
cab6ff9b 40
41
42ClassImp(AliITSVertexerTracks)
43
44
45//----------------------------------------------------------------------------
46AliITSVertexerTracks::AliITSVertexerTracks():AliITSVertexer() {
47//
48// Default constructor
49//
11ba84a4 50 fInFile = 0;
51 fOutFile = 0;
cab6ff9b 52 SetVtxStart();
53 SetMinTracks();
cab6ff9b 54 fTrksToSkip = 0;
55 fNTrksToSkip = 0;
56 for(Int_t i=0; i<3; i++)fInitPos[i] = 0.;
57}
cab6ff9b 58//----------------------------------------------------------------------------
11ba84a4 59AliITSVertexerTracks::AliITSVertexerTracks(TFile *inFile,TFile *outFile,
c84a5e9e 60 const AliMagF *map,
11ba84a4 61 Int_t fEv,Int_t lEv,
62 Double_t xStart,Double_t yStart) {
cab6ff9b 63//
64// Standard constructor
65//
11ba84a4 66 fCurrentVertex = 0;
67 fInFile = inFile;
68 fOutFile = outFile;
69 SetFirstEvent(fEv);
70 SetLastEvent(lEv);
c84a5e9e 71 SetFieldMap(map);
cab6ff9b 72 SetVtxStart(xStart,yStart);
73 SetMinTracks();
cab6ff9b 74 fTrksToSkip = 0;
75 fNTrksToSkip = 0;
11ba84a4 76 for(Int_t i=0; i<3; i++) fInitPos[i] = 0.;
77 SetDebug();
cab6ff9b 78}
79//----------------------------------------------------------------------------
c84a5e9e 80AliITSVertexerTracks::AliITSVertexerTracks(const AliMagF *map, TString fn,
11ba84a4 81 Double_t xStart,Double_t yStart)
82 :AliITSVertexer(fn) {
83//
84// Alternative constructor
85//
86 fInFile = 0;
87 fOutFile = 0;
c84a5e9e 88 SetFieldMap(map);
11ba84a4 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}
a48d1ed3 95//______________________________________________________________________
96AliITSVertexerTracks::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//______________________________________________________________________
103AliITSVertexerTracks& 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
11ba84a4 110//-----------------------------------------------------------------------------
111AliITSVertexerTracks::~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//-----------------------------------------------------------------------------
cab6ff9b 121Bool_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//---------------------------------------------------------------------------
137void 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//---------------------------------------------------------------------------
147void 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
11ba84a4 155 TDirectory *curdir = 0;
cab6ff9b 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();
11ba84a4 169
170 // write vertex to file
d4221964 171 TString vtxName = "Vertex_";
172 vtxName += ev;
11ba84a4 173 fCurrentVertex->SetName(vtxName.Data());
d4221964 174 fCurrentVertex->SetTitle("VertexerTracks");
11ba84a4 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//---------------------------------------------------------------------------
186void 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;
cab6ff9b 230 } // loop over events
231
232 return;
233}
234//----------------------------------------------------------------------------
235Int_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//----------------------------------------------------------------------------
290void 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);
cab6ff9b 298
299 return;
300}
301//----------------------------------------------------------------------------
d681bb2d 302AliESDVertex* AliITSVertexerTracks::FindVertexForCurrentEvent(Int_t evnumb) {
cab6ff9b 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;
11ba84a4 311 TTree *trkTree=(TTree*)fInFile->Get(treeName.Data());
cab6ff9b 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);
cab6ff9b 326 VertexFitter();
327 if(fDebug) printf(" vertex fit completed\n");
328
cab6ff9b 329 return fCurrentVertex;
330}
cab6ff9b 331//----------------------------------------------------------------------------
d681bb2d 332AliESDVertex* AliITSVertexerTracks::FindVertexForCurrentEvent(AliESD *esdEvent)
11ba84a4 333{
cab6ff9b 334//
11ba84a4 335// Vertex for current ESD event
cab6ff9b 336//
11ba84a4 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);
c84a5e9e 348 if(!esdTrack->GetStatus()&AliESDtrack::kITSin)
11ba84a4 349 { delete esdTrack; continue; }
70d3d424 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
11ba84a4 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; }
cab6ff9b 370
11ba84a4 371 // Set initial vertex position from ESD
2257f27e 372 esdEvent->GetVertex()->GetXYZ(vtx);
11ba84a4 373 SetVtxStart(vtx[0],vtx[1]);
cab6ff9b 374
11ba84a4 375 // VERTEX FINDER
376 VertexFinder();
cab6ff9b 377
11ba84a4 378 // VERTEX FITTER
379 ComputeMaxChi2PerTrack(nTrks);
380 VertexFitter();
381 if(fDebug) printf(" vertex fit completed\n");
cab6ff9b 382
11ba84a4 383 // store vertex information in ESD
384 fCurrentVertex->GetXYZ(vtx);
385 fCurrentVertex->GetCovMatrix(cvtx);
70d3d424 386
387 Double_t tp[3];
388 esdEvent->GetVertex()->GetTruePos(tp);
389 fCurrentVertex->SetTruePos(tp);
390
2257f27e 391 esdEvent->SetVertex(fCurrentVertex);
cab6ff9b 392
11ba84a4 393 cout<<"Vertex: "<<vtx[0]<<", "<<vtx[1]<<", "<<vtx[2]<<endl;
394 return fCurrentVertex;
cab6ff9b 395}
396//---------------------------------------------------------------------------
11ba84a4 397void AliITSVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped) {
cab6ff9b 398//
11ba84a4 399// Mark the tracks not ot be used in the vertex finding
cab6ff9b 400//
11ba84a4 401 fNTrksToSkip = n;
402 fTrksToSkip = new Int_t[n];
403 for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i];
404 return;
cab6ff9b 405}
406//---------------------------------------------------------------------------
407void 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//
d681bb2d 412 fCurrentVertex = new AliESDVertex(0.,0.,-1);
cab6ff9b 413 return;
414}
415//---------------------------------------------------------------------------
416void 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
427fInitPos[0] = fNominalPos[0]+gRandom->Gaus(0.,0.0100); // 100 micron gaussian smearing
428fInitPos[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;
96bd751d 440 AliITSStrLine line1;
441 AliITSStrLine line2;
cab6ff9b 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]);
96bd751d 462 line1.SetP0(pos1);
463 line1.SetCd(mom1);
cab6ff9b 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]);
96bd751d 476 line2.SetP0(pos2);
477 line2.SetCd(mom2);
cab6ff9b 478 Double_t crosspoint[3];
96bd751d 479 Int_t retcode = line2.Cross(&line1,crosspoint);
cab6ff9b 480 if(retcode<0){
481 if(fDebug>10)cout<<" i= "<<i<<", j= "<<j<<endl;
482 if(fDebug>10)cout<<"bad intersection\n";
96bd751d 483 line1.PrintStatus();
484 line2.PrintStatus();
cab6ff9b 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 }
cab6ff9b 493 }
cab6ff9b 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}
cab6ff9b 507//---------------------------------------------------------------------------
508void 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);
a48d1ed3 522 TMatrixD vV(3,3);
cab6ff9b 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
a48d1ed3 548 TMatrixD sumWiri(3,1);
549 TMatrixD sumWi(3,3);
cab6ff9b 550 for(i=0; i<3; i++) {
a48d1ed3 551 sumWiri(i,0) = 0.;
552 for(j=0; j<3; j++) sumWi(j,i) = 0.;
cab6ff9b 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
11ba84a4 563 rotAngle = alpha;
cab6ff9b 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);
a48d1ed3 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.;
cab6ff9b 582
583 // covariance matrix of local (y,z) - inverted
a48d1ed3 584 TMatrixD uUi(2,2);
cab6ff9b 585 t->GetExternalCovariance(cc);
a48d1ed3 586 uUi(0,0) = cc[0];
587 uUi(0,1) = cc[1];
588 uUi(1,0) = cc[1];
589 uUi(1,1) = cc[2];
cab6ff9b 590
a48d1ed3 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);
cab6ff9b 596
597 // track chi2
598 TMatrixD deltar = rv; deltar -= ri;
a48d1ed3 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);
cab6ff9b 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
a48d1ed3 613 TMatrixD wWiri(wWi,TMatrixD::kMult,ri);
cab6ff9b 614
a48d1ed3 615 sumWiri += wWiri;
616 sumWi += wWi;
cab6ff9b 617
618 nUsedTrks++;
619 } // end loop on tracks
620
621 if(nUsedTrks < fMinTracks) {
622 failed=1;
623 continue;
624 }
625
a48d1ed3 626 Double_t determinant = sumWi.Determinant();
cab6ff9b 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
a48d1ed3 635 TMatrixD invsumWi(TMatrixD::kInverted,sumWi);
636 vV = invsumWi;
cab6ff9b 637
638 // position of primary vertex
a48d1ed3 639 rv.Mult(vV,sumWiri);
cab6ff9b 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];
a48d1ed3 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);
cab6ff9b 662
663 // store data in the vertex object
d681bb2d 664 fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nUsedTrks);
cab6ff9b 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//----------------------------------------------------------------------------
d681bb2d 675AliESDVertex *AliITSVertexerTracks::VertexOnTheFly(TTree &trkTree) {
cab6ff9b 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);
cab6ff9b 691 VertexFitter();
692 if(fDebug) printf(" vertex fit completed\n");
693
694 return fCurrentVertex;
695}
696//----------------------------------------------------------------------------
697
698
699