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