]>
Commit | Line | Data |
---|---|---|
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 <TFile.h> | |
30 | #include <TTree.h> | |
31 | #include <TMatrixD.h> | |
32 | //---- AliRoot headers ----- | |
33 | #include "AliESDVertex.h" | |
34 | #include "AliITSVertexerTracks.h" | |
35 | #include "AliESD.h" | |
36 | #include "AliESDtrack.h" | |
37 | #include "AliVertexerTracks.h" | |
38 | ||
39 | ||
40 | ClassImp(AliITSVertexerTracks) | |
41 | ||
42 | ||
43 | //---------------------------------------------------------------------------- | |
44 | AliITSVertexerTracks::AliITSVertexerTracks():AliITSVertexer() { | |
45 | // | |
46 | // Default constructor | |
47 | // | |
48 | fInFile = 0; | |
49 | fOutFile = 0; | |
50 | SetVtxStart(); | |
51 | SetMinTracks(); | |
52 | fTrksToSkip = 0; | |
53 | fNTrksToSkip = 0; | |
54 | ||
55 | } | |
56 | //---------------------------------------------------------------------------- | |
57 | AliITSVertexerTracks::AliITSVertexerTracks(TFile *inFile,TFile *outFile, | |
58 | Int_t fEv,Int_t lEv, | |
59 | Double_t xStart,Double_t yStart) { | |
60 | // | |
61 | // Standard constructor | |
62 | // | |
63 | fCurrentVertex = 0; | |
64 | fInFile = inFile; | |
65 | fOutFile = outFile; | |
66 | SetFirstEvent(fEv); | |
67 | SetLastEvent(lEv); | |
68 | SetVtxStart(xStart,yStart); | |
69 | SetMinTracks(); | |
70 | fTrksToSkip = 0; | |
71 | fNTrksToSkip = 0; | |
72 | SetDebug(); | |
73 | } | |
74 | //---------------------------------------------------------------------------- | |
75 | AliITSVertexerTracks::AliITSVertexerTracks(TString fn, | |
76 | Double_t xStart,Double_t yStart) | |
77 | :AliITSVertexer(fn) { | |
78 | // | |
79 | // Alternative constructor | |
80 | // | |
81 | fInFile = 0; | |
82 | fOutFile = 0; | |
83 | SetVtxStart(xStart,yStart); | |
84 | SetMinTracks(); | |
85 | fTrksToSkip = 0; | |
86 | fNTrksToSkip = 0; | |
87 | } | |
88 | //______________________________________________________________________ | |
89 | AliITSVertexerTracks::AliITSVertexerTracks(const AliITSVertexerTracks &vtxr) : AliITSVertexer(vtxr) { | |
90 | // Copy constructor | |
91 | // Copies are not allowed. The method is protected to avoid misuse. | |
92 | Error("AliITSVertexerTracks","Copy constructor not allowed\n"); | |
93 | } | |
94 | ||
95 | //______________________________________________________________________ | |
96 | AliITSVertexerTracks& AliITSVertexerTracks::operator=(const AliITSVertexerTracks& /* vtxr */){ | |
97 | // Assignment operator | |
98 | // Assignment is not allowed. The method is protected to avoid misuse. | |
99 | Error("= operator","Assignment operator not allowed\n"); | |
100 | return *this; | |
101 | } | |
102 | ||
103 | //----------------------------------------------------------------------------- | |
104 | AliITSVertexerTracks::~AliITSVertexerTracks() { | |
105 | // Default Destructor | |
106 | // The objects poited by the following pointers are not owned | |
107 | // by this class and are not deleted | |
108 | ||
109 | fCurrentVertex = 0; | |
110 | fInFile = 0; | |
111 | fOutFile = 0; | |
112 | } | |
113 | //----------------------------------------------------------------------------- | |
114 | Bool_t AliITSVertexerTracks::CheckField() const { | |
115 | // | |
116 | // Check if the conv. const. has been set | |
117 | // | |
118 | Double_t field = AliTracker::GetBz(); // in kG | |
119 | ||
120 | if(field<1 || field>6) { | |
121 | printf("AliITSVertexerTracks::CheckField():\n ERROR: AliKalmanTrack::fConvConst not set\n Use AliKalmanTrack::SetConvConst() or AliITSVertexerTracks::SetField()\n"); | |
122 | return kFALSE; | |
123 | } | |
124 | printf("AliITSVertexerTracks::CheckField(): Using B = %3.1f T\n",field); | |
125 | return kTRUE; | |
126 | } | |
127 | //--------------------------------------------------------------------------- | |
128 | void AliITSVertexerTracks::ComputeMaxChi2PerTrack(Int_t nTracks) { | |
129 | // | |
130 | // Max. contr. to the chi2 has been tuned as a function of multiplicity | |
131 | // | |
132 | if(nTracks < 7) { fMaxChi2PerTrack = 1.e6; | |
133 | } else { fMaxChi2PerTrack = 100.; } | |
134 | ||
135 | return; | |
136 | } | |
137 | //--------------------------------------------------------------------------- | |
138 | void AliITSVertexerTracks::FindVertices() { | |
139 | // | |
140 | // Vertices for all events from fFirstEvent to fLastEvent | |
141 | // | |
142 | ||
143 | // Check if the conv. const. has been set | |
144 | if(!CheckField()) return; | |
145 | ||
146 | TDirectory *curdir = 0; | |
147 | ||
148 | // loop over events | |
149 | for(Int_t ev=fFirstEvent; ev<=fLastEvent; ev++) { | |
150 | if(ev % 100 == 0 || fDebug) printf("--- Processing event %d of %d ---\n",ev,fLastEvent); | |
151 | ||
152 | FindPrimaryVertexForCurrentEvent(ev); | |
153 | ||
154 | if(!fCurrentVertex) { | |
155 | printf("AliITSVertexerTracks::FindVertices(): no tracks tree for event %d\n",ev); | |
156 | continue; | |
157 | } | |
158 | ||
159 | if(fDebug) fCurrentVertex->PrintStatus(); | |
160 | ||
161 | // write vertex to file | |
162 | TString vtxName = "Vertex_"; | |
163 | vtxName += ev; | |
164 | fCurrentVertex->SetName(vtxName.Data()); | |
165 | fCurrentVertex->SetTitle("VertexerTracks"); | |
166 | //WriteCurrentVertex(); | |
167 | curdir = gDirectory; | |
168 | fOutFile->cd(); | |
169 | fCurrentVertex->Write(); | |
170 | curdir->cd(); | |
171 | fCurrentVertex = 0; | |
172 | } // loop over events | |
173 | ||
174 | return; | |
175 | } | |
176 | //--------------------------------------------------------------------------- | |
177 | void AliITSVertexerTracks::FindVerticesESD() { | |
178 | // | |
179 | // Vertices for all events from fFirstEvent to fLastEvent | |
180 | // | |
181 | ||
182 | // Check if the conv. const. has been set | |
183 | if(!CheckField()) return; | |
184 | ||
185 | TDirectory *curdir = 0; | |
186 | ||
187 | fInFile->cd(); | |
188 | TTree *esdTree = (TTree*)fInFile->Get("esdTree"); | |
189 | if(!esdTree) { | |
190 | printf("AliITSVertexerTracks::FindVerticesESD(): no tree in file!\n"); | |
191 | return; | |
192 | } | |
193 | Int_t nev = (Int_t)esdTree->GetEntries(); | |
194 | Int_t ev; | |
195 | // loop on events in tree | |
196 | for(Int_t i=0; i<nev; i++) { | |
197 | AliESD *esdEvent = new AliESD; | |
198 | esdTree->SetBranchAddress("ESD",&esdEvent); | |
199 | if(!esdTree->GetEvent(i)) { | |
200 | printf("AliITSVertexerTracks::FindVerticesESD(): not an ESD!\n"); | |
201 | delete esdEvent; | |
202 | return; | |
203 | } | |
204 | ev = (Int_t)esdEvent->GetEventNumber(); | |
205 | if(ev<fFirstEvent || ev>fLastEvent) { delete esdEvent; continue; } | |
206 | if(ev % 100 == 0 || fDebug) | |
207 | printf("--- Processing event %d of %d ---\n",ev,fLastEvent-fFirstEvent); | |
208 | ||
209 | FindPrimaryVertexForCurrentEvent(esdEvent); | |
210 | ||
211 | if(!fCurrentVertex) { | |
212 | printf("AliITSVertexerTracks::FindVertixesESD():\n no vertex for event %d\n",ev); | |
213 | continue; | |
214 | } | |
215 | ||
216 | if(fDebug) fCurrentVertex->PrintStatus(); | |
217 | ||
218 | // write vertex to file | |
219 | TString vtxName = "Vertex_"; | |
220 | vtxName += ev; | |
221 | fCurrentVertex->SetName(vtxName.Data()); | |
222 | fCurrentVertex->SetTitle("VertexerTracks"); | |
223 | //WriteCurrentVertex(); | |
224 | curdir = gDirectory; | |
225 | fOutFile->cd(); | |
226 | fCurrentVertex->Write(); | |
227 | curdir->cd(); | |
228 | fCurrentVertex = 0; | |
229 | esdEvent = 0; | |
230 | delete esdEvent; | |
231 | } // end loop over events | |
232 | ||
233 | return; | |
234 | } | |
235 | //---------------------------------------------------------------------------- | |
236 | Int_t AliITSVertexerTracks::PrepareTracks(TTree &trkTree) { | |
237 | // | |
238 | // Propagate tracks to initial vertex position and store them in a TObjArray | |
239 | // | |
240 | Double_t maxd0rphi = 3.; | |
241 | Double_t alpha,xlStart,d0rphi; | |
242 | Int_t nTrks = 0; | |
243 | Bool_t skipThis; | |
244 | Int_t nEntries = (Int_t)trkTree.GetEntries(); | |
245 | ||
246 | Double_t field=AliTracker::GetBz(); | |
247 | ||
248 | if(!fTrkArray.IsEmpty()) fTrkArray.Clear(); | |
249 | fTrkArray.Expand(nEntries); | |
250 | ||
251 | if(fDebug) { | |
252 | printf(" PrepareTracks()\n"); | |
253 | // trkTree.Print(); | |
254 | } | |
255 | ||
256 | for(Int_t i=0; i<nEntries; i++) { | |
257 | // check tracks to skip | |
258 | skipThis = kFALSE; | |
259 | for(Int_t j=0; j<fNTrksToSkip; j++) { | |
260 | if(i==fTrksToSkip[j]) { | |
261 | if(fDebug) printf("skipping track: %d\n",i); | |
262 | skipThis = kTRUE; | |
263 | } | |
264 | } | |
265 | if(skipThis) continue; | |
266 | ||
267 | AliESDtrack *track = new AliESDtrack; | |
268 | trkTree.SetBranchAddress("tracks",&track); | |
269 | trkTree.GetEvent(i); | |
270 | ||
271 | // propagate track to vtxSeed | |
272 | alpha = track->GetAlpha(); | |
273 | xlStart = fNominalPos[0]*TMath::Cos(alpha)+fNominalPos[1]*TMath::Sin(alpha); | |
274 | track->PropagateTo(xlStart,field); // to vtxSeed | |
275 | ||
276 | // select tracks with d0rphi < maxd0rphi | |
277 | ||
278 | d0rphi = TMath::Abs(track->GetD(fNominalPos[0],fNominalPos[1],field)); | |
279 | if(d0rphi > maxd0rphi) { delete track; continue; } | |
280 | ||
281 | ||
282 | fTrkArray.AddLast(track); | |
283 | ||
284 | nTrks++; | |
285 | if(fDebug)cout<<" :-) nTrks, d0rphi "<<nTrks<<" "<<d0rphi<<endl; | |
286 | ||
287 | } | |
288 | if(fTrksToSkip) delete [] fTrksToSkip; | |
289 | ||
290 | return nTrks; | |
291 | } | |
292 | //---------------------------------------------------------------------------- | |
293 | void AliITSVertexerTracks::PrintStatus() const { | |
294 | // | |
295 | // Print status | |
296 | // printf(" Initial position (%f,%f)\n",fNominalPos[0],fNominalPos[1]); | |
297 | printf(" Number of tracks in array: %d\n",(Int_t)fTrkArray.GetEntriesFast()); | |
298 | printf(" Minimum # tracks required in fit: %d\n",fMinTracks); | |
299 | ||
300 | return; | |
301 | } | |
302 | //---------------------------------------------------------------------------- | |
303 | AliVertex* AliITSVertexerTracks::VertexForSelectedTracks(AliESD *esdEvent,Int_t nofCand, Int_t *trkPos, Int_t opt){ | |
304 | ||
305 | // | |
306 | // Computes the vertex for selected tracks | |
307 | // trkPos=vector with track positions in ESD | |
308 | // | |
309 | Double_t vtx[3]; | |
310 | esdEvent->GetVertex()->GetXYZ(vtx); | |
311 | TTree *trkTree = new TTree("TreeT","tracks"); | |
312 | AliESDtrack *esdTrack = 0; | |
313 | trkTree->Branch("tracks","AliESDtrack",&esdTrack); | |
314 | for(Int_t i=0; i<nofCand;i++){ | |
315 | esdTrack = (AliESDtrack*)esdEvent->GetTrack(trkPos[i]); | |
316 | ||
317 | if(!esdTrack->GetStatus()&AliESDtrack::kTPCin) continue; | |
318 | if(!esdTrack->GetStatus()&AliESDtrack::kITSin) continue; | |
319 | if(!esdTrack->GetStatus()&AliESDtrack::kITSrefit) continue; | |
320 | ||
321 | Int_t nclus=esdTrack->GetNcls(0); // check number of clusters in ITS | |
322 | if(nclus<6) continue; | |
323 | trkTree->Fill(); | |
324 | } | |
325 | delete esdTrack; | |
326 | Int_t nTrks = PrepareTracks(*trkTree); | |
327 | //delete trkTree;// :-)) | |
328 | if(fDebug) printf(" tracks prepared: %d\n",nTrks); | |
329 | if(nTrks < fMinTracks) { | |
330 | printf("TooFewTracks\n"); | |
331 | AliVertex *theVert=new AliVertex(); | |
332 | theVert->SetDispersion(999); | |
333 | theVert->SetNContributors(-5); | |
334 | return theVert; | |
335 | } | |
336 | ||
337 | AliVertexerTracks *vertexer=new AliVertexerTracks(vtx[0],vtx[1]); | |
338 | vertexer->SetFinderAlgorithm(opt); | |
339 | AliVertex *theVert=(AliVertex*)vertexer->VertexForSelectedTracks(&fTrkArray); | |
340 | // beware: newvt object should be deleted by the caller | |
341 | AliVertex *newvt = new AliVertex(*theVert); | |
342 | delete vertexer; | |
343 | return newvt; | |
344 | } | |
345 | //---------------------------------------------------------------------------- | |
346 | AliESDVertex* AliITSVertexerTracks::FindPrimaryVertexForCurrentEvent(Int_t evnumb) { | |
347 | // | |
348 | // Vertex for current event | |
349 | // | |
350 | fCurrentVertex = 0; | |
351 | ||
352 | // get tree with tracks from input file | |
353 | fInFile->cd(); | |
354 | TTree *esdTree = (TTree*)fInFile->Get("esdTree"); | |
355 | ||
356 | if(!esdTree) { | |
357 | printf("AliITSVertexerTracks::FindPrimaryVertexForCurrentEvent(): no tree in file!\n"); | |
358 | return fCurrentVertex; | |
359 | } | |
360 | AliESD *esdEvent = new AliESD; | |
361 | esdTree->SetBranchAddress("ESD",&esdEvent); | |
362 | esdTree->GetEvent(evnumb); | |
363 | return FindPrimaryVertexForCurrentEvent(esdEvent); | |
364 | } | |
365 | //---------------------------------------------------------------------------- | |
366 | AliESDVertex* AliITSVertexerTracks::FindPrimaryVertexForCurrentEvent(AliESD *esdEvent) | |
367 | { | |
368 | // | |
369 | // Vertex for current ESD event | |
370 | // | |
371 | fCurrentVertex = 0; | |
372 | Double_t vtx[3],cvtx[6]; | |
373 | ||
374 | Int_t entr = (Int_t)esdEvent->GetNumberOfTracks(); | |
375 | TTree *trkTree = new TTree("TreeT","tracks"); | |
376 | AliESDtrack *esdTrack = 0; | |
377 | trkTree->Branch("tracks","AliESDtrack",&esdTrack); | |
378 | ||
379 | for(Int_t i=0; i<entr; i++) { | |
380 | AliESDtrack *et = esdEvent->GetTrack(i); | |
381 | esdTrack = new AliESDtrack(*et); | |
382 | if(!esdTrack->GetStatus()&AliESDtrack::kITSin) continue; | |
383 | if(!esdTrack->GetStatus()&AliESDtrack::kITSrefit) continue; | |
384 | Int_t nclus=esdTrack->GetNcls(0); // check number of clusters in ITS | |
385 | if(nclus<5) continue; | |
386 | ||
387 | trkTree->Fill(); | |
388 | } | |
389 | delete esdTrack; | |
390 | ||
391 | // preselect tracks and propagate them to initial vertex position | |
392 | Int_t nTrks = PrepareTracks(*trkTree); | |
393 | delete trkTree; | |
394 | if(fDebug) printf(" tracks prepared: %d\n",nTrks); | |
395 | if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; } | |
396 | ||
397 | // Set initial vertex position from ESD | |
398 | esdEvent->GetVertex()->GetXYZ(vtx); | |
399 | SetVtxStart(vtx[0],vtx[1]); | |
400 | ||
401 | // VERTEX FITTER | |
402 | ComputeMaxChi2PerTrack(nTrks); | |
403 | VertexFitter(); | |
404 | if(fDebug) printf(" vertex fit completed\n"); | |
405 | ||
406 | // store vertex information in ESD | |
407 | fCurrentVertex->GetXYZ(vtx); | |
408 | fCurrentVertex->GetCovMatrix(cvtx); | |
409 | ||
410 | Double_t tp[3]; | |
411 | esdEvent->GetVertex()->GetTruePos(tp); | |
412 | fCurrentVertex->SetTruePos(tp); | |
413 | ||
414 | return fCurrentVertex; | |
415 | } | |
416 | //--------------------------------------------------------------------------- | |
417 | void AliITSVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped) { | |
418 | // | |
419 | // Mark the tracks not ot be used in the vertex finding | |
420 | // | |
421 | fNTrksToSkip = n; | |
422 | fTrksToSkip = new Int_t[n]; | |
423 | for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i]; | |
424 | return; | |
425 | } | |
426 | //--------------------------------------------------------------------------- | |
427 | void AliITSVertexerTracks::TooFewTracks() { | |
428 | // | |
429 | // When the number of tracks is < fMinTracks the vertex is set to (0,0,0) | |
430 | // and the number of tracks to -1 | |
431 | // | |
432 | fCurrentVertex = new AliESDVertex(0.,0.,-1); | |
433 | return; | |
434 | } | |
435 | //--------------------------------------------------------------------------- | |
436 | void AliITSVertexerTracks::VertexFitter() { | |
437 | // | |
438 | // The optimal estimate of the vertex position is given by a "weighted | |
439 | // average of tracks positions" | |
440 | // Original method: CMS Note 97/0051 | |
441 | // | |
442 | if(fDebug) { | |
443 | printf(" VertexFitter(): start\n"); | |
444 | PrintStatus(); | |
445 | } | |
446 | AliVertexerTracks *vertexer=new AliVertexerTracks(fNominalPos[0],fNominalPos[1]); | |
447 | vertexer->SetFinderAlgorithm(1); | |
448 | AliVertex *thevert=(AliVertex*)vertexer->VertexForSelectedTracks(&fTrkArray); | |
449 | Double_t initPos[3]; | |
450 | thevert->GetXYZ(initPos); | |
451 | // cout<<"Finder: "<<initPos[0]<<"; "<<initPos[1]<<"; "<<initPos[2]<<endl; | |
452 | delete vertexer; | |
453 | ||
454 | ||
455 | Int_t i,j,k,step=0; | |
456 | TMatrixD rv(3,1); | |
457 | TMatrixD vV(3,3); | |
458 | rv(0,0) = initPos[0]; | |
459 | rv(1,0) = initPos[1]; | |
460 | rv(2,0) = 0.; | |
461 | Double_t xlStart,alpha; | |
462 | Double_t rotAngle; | |
463 | Double_t cosRot,sinRot; | |
464 | Double_t cc[15]; | |
465 | Int_t nUsedTrks; | |
466 | Double_t chi2,chi2i; | |
467 | Int_t arrEntries = (Int_t)fTrkArray.GetEntries(); | |
468 | AliESDtrack *t = 0; | |
469 | Int_t failed = 0; | |
470 | ||
471 | Int_t *skipTrack = new Int_t[arrEntries]; | |
472 | for(i=0; i<arrEntries; i++) skipTrack[i]=0; | |
473 | ||
474 | // 3 steps: | |
475 | // 1st - first estimate of vtx using all tracks | |
476 | // 2nd - apply cut on chi2 max per track | |
477 | // 3rd - estimate of global chi2 | |
478 | for(step=0; step<3; step++) { | |
479 | if(fDebug) printf(" step = %d\n",step); | |
480 | chi2 = 0.; | |
481 | nUsedTrks = 0; | |
482 | ||
483 | TMatrixD sumWiri(3,1); | |
484 | TMatrixD sumWi(3,3); | |
485 | for(i=0; i<3; i++) { | |
486 | sumWiri(i,0) = 0.; | |
487 | for(j=0; j<3; j++) sumWi(j,i) = 0.; | |
488 | } | |
489 | ||
490 | // loop on tracks | |
491 | for(k=0; k<arrEntries; k++) { | |
492 | if(skipTrack[k]) continue; | |
493 | // get track from track array | |
494 | t = (AliESDtrack*)fTrkArray.At(k); | |
495 | alpha = t->GetAlpha(); | |
496 | xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha); | |
497 | t->PropagateTo(xlStart,AliTracker::GetBz()); // to vtxSeed | |
498 | rotAngle = alpha; | |
499 | if(alpha<0.) rotAngle += 2.*TMath::Pi(); | |
500 | cosRot = TMath::Cos(rotAngle); | |
501 | sinRot = TMath::Sin(rotAngle); | |
502 | ||
503 | // vector of track global coordinates | |
504 | TMatrixD ri(3,1); | |
505 | ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot; | |
506 | ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot; | |
507 | ri(2,0) = t->GetZ(); | |
508 | ||
509 | // matrix to go from global (x,y,z) to local (y,z); | |
510 | TMatrixD qQi(2,3); | |
511 | qQi(0,0) = -sinRot; | |
512 | qQi(0,1) = cosRot; | |
513 | qQi(0,2) = 0.; | |
514 | qQi(1,0) = 0.; | |
515 | qQi(1,1) = 0.; | |
516 | qQi(1,2) = 1.; | |
517 | ||
518 | // covariance matrix of local (y,z) - inverted | |
519 | TMatrixD uUi(2,2); | |
520 | t->GetExternalCovariance(cc); | |
521 | uUi(0,0) = cc[0]; | |
522 | uUi(0,1) = cc[1]; | |
523 | uUi(1,0) = cc[1]; | |
524 | uUi(1,1) = cc[2]; | |
525 | ||
526 | // weights matrix: wWi = qQiT * uUiInv * qQi | |
527 | if(uUi.Determinant() <= 0.) continue; | |
528 | TMatrixD uUiInv(TMatrixD::kInverted,uUi); | |
529 | TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi); | |
530 | TMatrixD wWi(qQi,TMatrixD::kTransposeMult,uUiInvQi); | |
531 | ||
532 | // track chi2 | |
533 | TMatrixD deltar = rv; deltar -= ri; | |
534 | TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar); | |
535 | chi2i = deltar(0,0)*wWideltar(0,0)+ | |
536 | deltar(1,0)*wWideltar(1,0)+ | |
537 | deltar(2,0)*wWideltar(2,0); | |
538 | ||
539 | ||
540 | if(step==1 && chi2i > fMaxChi2PerTrack) { | |
541 | skipTrack[k] = 1; | |
542 | continue; | |
543 | } | |
544 | ||
545 | // add to total chi2 | |
546 | chi2 += chi2i; | |
547 | ||
548 | TMatrixD wWiri(wWi,TMatrixD::kMult,ri); | |
549 | ||
550 | sumWiri += wWiri; | |
551 | sumWi += wWi; | |
552 | ||
553 | nUsedTrks++; | |
554 | } // end loop on tracks | |
555 | ||
556 | if(nUsedTrks < fMinTracks) { | |
557 | failed=1; | |
558 | continue; | |
559 | } | |
560 | ||
561 | Double_t determinant = sumWi.Determinant(); | |
562 | //cerr<<" determinant: "<<determinant<<endl; | |
563 | if(determinant < 100.) { | |
564 | printf("det(V) = 0\n"); | |
565 | failed=1; | |
566 | continue; | |
567 | } | |
568 | ||
569 | // inverted of weights matrix | |
570 | TMatrixD invsumWi(TMatrixD::kInverted,sumWi); | |
571 | vV = invsumWi; | |
572 | ||
573 | // position of primary vertex | |
574 | rv.Mult(vV,sumWiri); | |
575 | ||
576 | } // end loop on the 3 steps | |
577 | ||
578 | delete [] skipTrack; | |
579 | delete t; | |
580 | ||
581 | if(failed) { | |
582 | TooFewTracks(); | |
583 | return; | |
584 | } | |
585 | ||
586 | Double_t position[3]; | |
587 | position[0] = rv(0,0); | |
588 | position[1] = rv(1,0); | |
589 | position[2] = rv(2,0); | |
590 | Double_t covmatrix[6]; | |
591 | covmatrix[0] = vV(0,0); | |
592 | covmatrix[1] = vV(0,1); | |
593 | covmatrix[2] = vV(1,1); | |
594 | covmatrix[3] = vV(0,2); | |
595 | covmatrix[4] = vV(1,2); | |
596 | covmatrix[5] = vV(2,2); | |
597 | ||
598 | // store data in the vertex object | |
599 | fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nUsedTrks); | |
600 | ||
601 | if(fDebug) { | |
602 | printf(" VertexFitter(): finish\n"); | |
603 | printf(" rv = ( %f , %f , %f )\n\n",rv(0,0),rv(1,0),rv(2,0)); | |
604 | fCurrentVertex->PrintStatus(); | |
605 | } | |
606 | ||
607 | return; | |
608 | } | |
609 | //---------------------------------------------------------------------------- | |
610 | AliESDVertex *AliITSVertexerTracks::VertexOnTheFly(TTree &trkTree) { | |
611 | // | |
612 | // Return vertex from tracks in trkTree | |
613 | // | |
614 | if(fCurrentVertex) fCurrentVertex = 0; | |
615 | ||
616 | // get tracks and propagate them to initial vertex position | |
617 | Int_t nTrks = PrepareTracks(*(&trkTree)); | |
618 | if(fDebug) printf(" tracks prepared: %d\n",nTrks); | |
619 | if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; } | |
620 | ||
621 | // VERTEX FITTER | |
622 | ComputeMaxChi2PerTrack(nTrks); | |
623 | VertexFitter(); | |
624 | if(fDebug) printf(" vertex fit completed\n"); | |
625 | ||
626 | return fCurrentVertex; | |
627 | } | |
628 | //---------------------------------------------------------------------------- | |
629 | ||
630 | ||
631 |