]>
Commit | Line | Data |
---|---|---|
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 -------- | |
29 | #include <TFile.h> | |
30 | #include <TTree.h> | |
cab6ff9b | 31 | #include <TMatrixD.h> |
cab6ff9b | 32 | //---- AliRoot headers ----- |
d681bb2d | 33 | #include "AliESDVertex.h" |
cab6ff9b | 34 | #include "AliITSVertexerTracks.h" |
11ba84a4 | 35 | #include "AliESD.h" |
36 | #include "AliESDtrack.h" | |
2d57349e | 37 | #include "AliVertexerTracks.h" |
cab6ff9b | 38 | |
39 | ||
40 | ClassImp(AliITSVertexerTracks) | |
41 | ||
42 | ||
43 | //---------------------------------------------------------------------------- | |
44 | AliITSVertexerTracks::AliITSVertexerTracks():AliITSVertexer() { | |
45 | // | |
46 | // Default constructor | |
47 | // | |
11ba84a4 | 48 | fInFile = 0; |
49 | fOutFile = 0; | |
cab6ff9b | 50 | SetVtxStart(); |
51 | SetMinTracks(); | |
cab6ff9b | 52 | fTrksToSkip = 0; |
53 | fNTrksToSkip = 0; | |
9c0755ea | 54 | |
cab6ff9b | 55 | } |
cab6ff9b | 56 | //---------------------------------------------------------------------------- |
11ba84a4 | 57 | AliITSVertexerTracks::AliITSVertexerTracks(TFile *inFile,TFile *outFile, |
11ba84a4 | 58 | Int_t fEv,Int_t lEv, |
59 | Double_t xStart,Double_t yStart) { | |
cab6ff9b | 60 | // |
61 | // Standard constructor | |
62 | // | |
11ba84a4 | 63 | fCurrentVertex = 0; |
64 | fInFile = inFile; | |
65 | fOutFile = outFile; | |
66 | SetFirstEvent(fEv); | |
67 | SetLastEvent(lEv); | |
cab6ff9b | 68 | SetVtxStart(xStart,yStart); |
69 | SetMinTracks(); | |
cab6ff9b | 70 | fTrksToSkip = 0; |
71 | fNTrksToSkip = 0; | |
11ba84a4 | 72 | SetDebug(); |
cab6ff9b | 73 | } |
74 | //---------------------------------------------------------------------------- | |
7d0f8548 | 75 | AliITSVertexerTracks::AliITSVertexerTracks(TString fn, |
11ba84a4 | 76 | Double_t xStart,Double_t yStart) |
77 | :AliITSVertexer(fn) { | |
78 | // | |
79 | // Alternative constructor | |
80 | // | |
81 | fInFile = 0; | |
82 | fOutFile = 0; | |
11ba84a4 | 83 | SetVtxStart(xStart,yStart); |
84 | SetMinTracks(); | |
85 | fTrksToSkip = 0; | |
86 | fNTrksToSkip = 0; | |
11ba84a4 | 87 | } |
a48d1ed3 | 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 | ||
11ba84a4 | 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 | //----------------------------------------------------------------------------- | |
cab6ff9b | 114 | Bool_t AliITSVertexerTracks::CheckField() const { |
115 | // | |
116 | // Check if the conv. const. has been set | |
117 | // | |
2d57349e | 118 | Double_t field = AliTracker::GetBz(); // in kG |
cab6ff9b | 119 | |
2d57349e | 120 | if(field<1 || field>6) { |
cab6ff9b | 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 | ||
11ba84a4 | 146 | TDirectory *curdir = 0; |
cab6ff9b | 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 | ||
66e811e2 | 152 | FindPrimaryVertexForCurrentEvent(ev); |
cab6ff9b | 153 | |
154 | if(!fCurrentVertex) { | |
2d57349e | 155 | printf("AliITSVertexerTracks::FindVertices(): no tracks tree for event %d\n",ev); |
cab6ff9b | 156 | continue; |
157 | } | |
158 | ||
159 | if(fDebug) fCurrentVertex->PrintStatus(); | |
11ba84a4 | 160 | |
161 | // write vertex to file | |
d4221964 | 162 | TString vtxName = "Vertex_"; |
163 | vtxName += ev; | |
11ba84a4 | 164 | fCurrentVertex->SetName(vtxName.Data()); |
d4221964 | 165 | fCurrentVertex->SetTitle("VertexerTracks"); |
11ba84a4 | 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 | ||
797b8b68 | 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; | |
cab6ff9b | 234 | } |
235 | //---------------------------------------------------------------------------- | |
236 | Int_t AliITSVertexerTracks::PrepareTracks(TTree &trkTree) { | |
66e811e2 | 237 | // |
238 | // Propagate tracks to initial vertex position and store them in a TObjArray | |
239 | // | |
240 | Double_t maxd0rphi = 3.; | |
2d57349e | 241 | Double_t alpha,xlStart,d0rphi; |
cab6ff9b | 242 | Int_t nTrks = 0; |
243 | Bool_t skipThis; | |
cab6ff9b | 244 | Int_t nEntries = (Int_t)trkTree.GetEntries(); |
245 | ||
2d57349e | 246 | Double_t field=AliTracker::GetBz(); |
247 | ||
cab6ff9b | 248 | if(!fTrkArray.IsEmpty()) fTrkArray.Clear(); |
249 | fTrkArray.Expand(nEntries); | |
250 | ||
251 | if(fDebug) { | |
252 | printf(" PrepareTracks()\n"); | |
9c0755ea | 253 | // trkTree.Print(); |
cab6ff9b | 254 | } |
2d57349e | 255 | |
cab6ff9b | 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 | ||
2d57349e | 267 | AliESDtrack *track = new AliESDtrack; |
268 | trkTree.SetBranchAddress("tracks",&track); | |
cab6ff9b | 269 | trkTree.GetEvent(i); |
270 | ||
2d57349e | 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 | |
cab6ff9b | 275 | |
2d57349e | 276 | // select tracks with d0rphi < maxd0rphi |
9c0755ea | 277 | |
2d57349e | 278 | d0rphi = TMath::Abs(track->GetD(fNominalPos[0],fNominalPos[1],field)); |
279 | if(d0rphi > maxd0rphi) { delete track; continue; } | |
9c0755ea | 280 | |
cab6ff9b | 281 | |
2d57349e | 282 | fTrkArray.AddLast(track); |
66e811e2 | 283 | |
cab6ff9b | 284 | nTrks++; |
66e811e2 | 285 | if(fDebug)cout<<" :-) nTrks, d0rphi "<<nTrks<<" "<<d0rphi<<endl; |
2d57349e | 286 | |
cab6ff9b | 287 | } |
cab6ff9b | 288 | if(fTrksToSkip) delete [] fTrksToSkip; |
2d57349e | 289 | |
cab6ff9b | 290 | return nTrks; |
2d57349e | 291 | } |
cab6ff9b | 292 | //---------------------------------------------------------------------------- |
293 | void AliITSVertexerTracks::PrintStatus() const { | |
294 | // | |
295 | // Print status | |
2d57349e | 296 | // printf(" Initial position (%f,%f)\n",fNominalPos[0],fNominalPos[1]); |
cab6ff9b | 297 | printf(" Number of tracks in array: %d\n",(Int_t)fTrkArray.GetEntriesFast()); |
298 | printf(" Minimum # tracks required in fit: %d\n",fMinTracks); | |
cab6ff9b | 299 | |
300 | return; | |
301 | } | |
302 | //---------------------------------------------------------------------------- | |
2d57349e | 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 | //---------------------------------------------------------------------------- | |
66e811e2 | 346 | AliESDVertex* AliITSVertexerTracks::FindPrimaryVertexForCurrentEvent(Int_t evnumb) { |
cab6ff9b | 347 | // |
348 | // Vertex for current event | |
349 | // | |
350 | fCurrentVertex = 0; | |
351 | ||
352 | // get tree with tracks from input file | |
2d57349e | 353 | fInFile->cd(); |
354 | TTree *esdTree = (TTree*)fInFile->Get("esdTree"); | |
cab6ff9b | 355 | |
2d57349e | 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); | |
cab6ff9b | 364 | } |
cab6ff9b | 365 | //---------------------------------------------------------------------------- |
66e811e2 | 366 | AliESDVertex* AliITSVertexerTracks::FindPrimaryVertexForCurrentEvent(AliESD *esdEvent) |
11ba84a4 | 367 | { |
cab6ff9b | 368 | // |
11ba84a4 | 369 | // Vertex for current ESD event |
cab6ff9b | 370 | // |
11ba84a4 | 371 | fCurrentVertex = 0; |
372 | Double_t vtx[3],cvtx[6]; | |
373 | ||
11ba84a4 | 374 | Int_t entr = (Int_t)esdEvent->GetNumberOfTracks(); |
2d57349e | 375 | TTree *trkTree = new TTree("TreeT","tracks"); |
376 | AliESDtrack *esdTrack = 0; | |
377 | trkTree->Branch("tracks","AliESDtrack",&esdTrack); | |
11ba84a4 | 378 | |
379 | for(Int_t i=0; i<entr; i++) { | |
9c0755ea | 380 | AliESDtrack *et = esdEvent->GetTrack(i); |
381 | esdTrack = new AliESDtrack(*et); | |
2d57349e | 382 | if(!esdTrack->GetStatus()&AliESDtrack::kITSin) continue; |
9c0755ea | 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 | ||
11ba84a4 | 387 | trkTree->Fill(); |
11ba84a4 | 388 | } |
2d57349e | 389 | delete esdTrack; |
11ba84a4 | 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; } | |
cab6ff9b | 396 | |
11ba84a4 | 397 | // Set initial vertex position from ESD |
2257f27e | 398 | esdEvent->GetVertex()->GetXYZ(vtx); |
11ba84a4 | 399 | SetVtxStart(vtx[0],vtx[1]); |
cab6ff9b | 400 | |
11ba84a4 | 401 | // VERTEX FITTER |
402 | ComputeMaxChi2PerTrack(nTrks); | |
403 | VertexFitter(); | |
404 | if(fDebug) printf(" vertex fit completed\n"); | |
cab6ff9b | 405 | |
11ba84a4 | 406 | // store vertex information in ESD |
407 | fCurrentVertex->GetXYZ(vtx); | |
408 | fCurrentVertex->GetCovMatrix(cvtx); | |
70d3d424 | 409 | |
410 | Double_t tp[3]; | |
411 | esdEvent->GetVertex()->GetTruePos(tp); | |
412 | fCurrentVertex->SetTruePos(tp); | |
413 | ||
11ba84a4 | 414 | return fCurrentVertex; |
cab6ff9b | 415 | } |
416 | //--------------------------------------------------------------------------- | |
11ba84a4 | 417 | void AliITSVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped) { |
cab6ff9b | 418 | // |
11ba84a4 | 419 | // Mark the tracks not ot be used in the vertex finding |
cab6ff9b | 420 | // |
11ba84a4 | 421 | fNTrksToSkip = n; |
422 | fTrksToSkip = new Int_t[n]; | |
423 | for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i]; | |
424 | return; | |
cab6ff9b | 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 | // | |
d681bb2d | 432 | fCurrentVertex = new AliESDVertex(0.,0.,-1); |
cab6ff9b | 433 | return; |
434 | } | |
435 | //--------------------------------------------------------------------------- | |
cab6ff9b | 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 | } | |
2d57349e | 446 | AliVertexerTracks *vertexer=new AliVertexerTracks(fNominalPos[0],fNominalPos[1]); |
9c0755ea | 447 | vertexer->SetFinderAlgorithm(1); |
2d57349e | 448 | AliVertex *thevert=(AliVertex*)vertexer->VertexForSelectedTracks(&fTrkArray); |
449 | Double_t initPos[3]; | |
450 | thevert->GetXYZ(initPos); | |
9c0755ea | 451 | // cout<<"Finder: "<<initPos[0]<<"; "<<initPos[1]<<"; "<<initPos[2]<<endl; |
2d57349e | 452 | delete vertexer; |
cab6ff9b | 453 | |
454 | ||
455 | Int_t i,j,k,step=0; | |
456 | TMatrixD rv(3,1); | |
a48d1ed3 | 457 | TMatrixD vV(3,3); |
66e811e2 | 458 | rv(0,0) = initPos[0]; |
459 | rv(1,0) = initPos[1]; | |
cab6ff9b | 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(); | |
2d57349e | 468 | AliESDtrack *t = 0; |
cab6ff9b | 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 | ||
a48d1ed3 | 483 | TMatrixD sumWiri(3,1); |
484 | TMatrixD sumWi(3,3); | |
cab6ff9b | 485 | for(i=0; i<3; i++) { |
a48d1ed3 | 486 | sumWiri(i,0) = 0.; |
487 | for(j=0; j<3; j++) sumWi(j,i) = 0.; | |
cab6ff9b | 488 | } |
489 | ||
490 | // loop on tracks | |
491 | for(k=0; k<arrEntries; k++) { | |
492 | if(skipTrack[k]) continue; | |
493 | // get track from track array | |
2d57349e | 494 | t = (AliESDtrack*)fTrkArray.At(k); |
cab6ff9b | 495 | alpha = t->GetAlpha(); |
66e811e2 | 496 | xlStart = initPos[0]*TMath::Cos(alpha)+initPos[1]*TMath::Sin(alpha); |
2d57349e | 497 | t->PropagateTo(xlStart,AliTracker::GetBz()); // to vtxSeed |
11ba84a4 | 498 | rotAngle = alpha; |
cab6ff9b | 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); | |
a48d1ed3 | 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.; | |
cab6ff9b | 517 | |
518 | // covariance matrix of local (y,z) - inverted | |
a48d1ed3 | 519 | TMatrixD uUi(2,2); |
cab6ff9b | 520 | t->GetExternalCovariance(cc); |
a48d1ed3 | 521 | uUi(0,0) = cc[0]; |
522 | uUi(0,1) = cc[1]; | |
523 | uUi(1,0) = cc[1]; | |
524 | uUi(1,1) = cc[2]; | |
cab6ff9b | 525 | |
a48d1ed3 | 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); | |
cab6ff9b | 531 | |
532 | // track chi2 | |
533 | TMatrixD deltar = rv; deltar -= ri; | |
a48d1ed3 | 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); | |
cab6ff9b | 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 | ||
a48d1ed3 | 548 | TMatrixD wWiri(wWi,TMatrixD::kMult,ri); |
cab6ff9b | 549 | |
a48d1ed3 | 550 | sumWiri += wWiri; |
551 | sumWi += wWi; | |
cab6ff9b | 552 | |
553 | nUsedTrks++; | |
554 | } // end loop on tracks | |
555 | ||
556 | if(nUsedTrks < fMinTracks) { | |
557 | failed=1; | |
558 | continue; | |
559 | } | |
560 | ||
a48d1ed3 | 561 | Double_t determinant = sumWi.Determinant(); |
cab6ff9b | 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 | |
a48d1ed3 | 570 | TMatrixD invsumWi(TMatrixD::kInverted,sumWi); |
571 | vV = invsumWi; | |
cab6ff9b | 572 | |
573 | // position of primary vertex | |
a48d1ed3 | 574 | rv.Mult(vV,sumWiri); |
cab6ff9b | 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]; | |
a48d1ed3 | 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); | |
cab6ff9b | 597 | |
598 | // store data in the vertex object | |
d681bb2d | 599 | fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nUsedTrks); |
cab6ff9b | 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 | //---------------------------------------------------------------------------- | |
d681bb2d | 610 | AliESDVertex *AliITSVertexerTracks::VertexOnTheFly(TTree &trkTree) { |
cab6ff9b | 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 | ||
cab6ff9b | 621 | // VERTEX FITTER |
622 | ComputeMaxChi2PerTrack(nTrks); | |
cab6ff9b | 623 | VertexFitter(); |
624 | if(fDebug) printf(" vertex fit completed\n"); | |
625 | ||
626 | return fCurrentVertex; | |
627 | } | |
628 | //---------------------------------------------------------------------------- | |
629 | ||
630 | ||
631 |