]>
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 | |
18 | // | |
19 | // Origin: A.Dainese, Padova, andrea.dainese@pd.infn.it | |
20 | // M.Masera, Torino, massimo.masera@to.infn.it | |
21 | //----------------------------------------------------------------- | |
22 | ||
23 | //---- standard headers ---- | |
24 | #include <Riostream.h> | |
25 | //---- Root headers -------- | |
11ba84a4 | 26 | #include <TKey.h> |
cab6ff9b | 27 | #include <TFile.h> |
28 | #include <TTree.h> | |
29 | #include <TVector3.h> | |
30 | #include <TMatrixD.h> | |
11ba84a4 | 31 | #include <TRandom.h> |
cab6ff9b | 32 | //---- AliRoot headers ----- |
33 | #include <AliRun.h> | |
34 | #include "AliKalmanTrack.h" | |
35 | #include "AliITSStrLine.h" | |
36 | #include "AliITStrackV2.h" | |
d681bb2d | 37 | #include "AliESDVertex.h" |
cab6ff9b | 38 | #include "AliITSVertexerTracks.h" |
11ba84a4 | 39 | #include "AliESD.h" |
40 | #include "AliESDtrack.h" | |
cab6ff9b | 41 | |
42 | ||
43 | ClassImp(AliITSVertexerTracks) | |
44 | ||
45 | ||
46 | //---------------------------------------------------------------------------- | |
47 | AliITSVertexerTracks::AliITSVertexerTracks():AliITSVertexer() { | |
48 | // | |
49 | // Default constructor | |
50 | // | |
11ba84a4 | 51 | fInFile = 0; |
52 | fOutFile = 0; | |
cab6ff9b | 53 | SetVtxStart(); |
54 | SetMinTracks(); | |
cab6ff9b | 55 | fTrksToSkip = 0; |
56 | fNTrksToSkip = 0; | |
57 | for(Int_t i=0; i<3; i++)fInitPos[i] = 0.; | |
58 | } | |
cab6ff9b | 59 | //---------------------------------------------------------------------------- |
11ba84a4 | 60 | AliITSVertexerTracks::AliITSVertexerTracks(TFile *inFile,TFile *outFile, |
61 | Double_t field, | |
62 | Int_t fEv,Int_t lEv, | |
63 | Double_t xStart,Double_t yStart) { | |
cab6ff9b | 64 | // |
65 | // Standard constructor | |
66 | // | |
11ba84a4 | 67 | fCurrentVertex = 0; |
68 | fInFile = inFile; | |
69 | fOutFile = outFile; | |
70 | SetFirstEvent(fEv); | |
71 | SetLastEvent(lEv); | |
cab6ff9b | 72 | SetField(field); |
73 | SetVtxStart(xStart,yStart); | |
74 | SetMinTracks(); | |
cab6ff9b | 75 | fTrksToSkip = 0; |
76 | fNTrksToSkip = 0; | |
11ba84a4 | 77 | for(Int_t i=0; i<3; i++) fInitPos[i] = 0.; |
78 | SetDebug(); | |
cab6ff9b | 79 | } |
80 | //---------------------------------------------------------------------------- | |
11ba84a4 | 81 | AliITSVertexerTracks::AliITSVertexerTracks(Double_t field, TString fn, |
82 | Double_t xStart,Double_t yStart) | |
83 | :AliITSVertexer(fn) { | |
84 | // | |
85 | // Alternative constructor | |
86 | // | |
87 | fInFile = 0; | |
88 | fOutFile = 0; | |
89 | SetField(field); | |
90 | SetVtxStart(xStart,yStart); | |
91 | SetMinTracks(); | |
92 | fTrksToSkip = 0; | |
93 | fNTrksToSkip = 0; | |
94 | for(Int_t i=0; i<3; i++) fInitPos[i] = 0.; | |
95 | } | |
96 | //----------------------------------------------------------------------------- | |
97 | AliITSVertexerTracks::~AliITSVertexerTracks() { | |
98 | // Default Destructor | |
99 | // The objects poited by the following pointers are not owned | |
100 | // by this class and are not deleted | |
101 | ||
102 | fCurrentVertex = 0; | |
103 | fInFile = 0; | |
104 | fOutFile = 0; | |
105 | } | |
106 | //----------------------------------------------------------------------------- | |
cab6ff9b | 107 | Bool_t AliITSVertexerTracks::CheckField() const { |
108 | // | |
109 | // Check if the conv. const. has been set | |
110 | // | |
111 | AliITStrackV2 t; | |
112 | Double_t cc = t.GetConvConst(); | |
113 | Double_t field = 100./0.299792458/cc; | |
114 | ||
115 | if(field<0.1 || field>0.6) { | |
116 | printf("AliITSVertexerTracks::CheckField():\n ERROR: AliKalmanTrack::fConvConst not set\n Use AliKalmanTrack::SetConvConst() or AliITSVertexerTracks::SetField()\n"); | |
117 | return kFALSE; | |
118 | } | |
119 | printf("AliITSVertexerTracks::CheckField(): Using B = %3.1f T\n",field); | |
120 | return kTRUE; | |
121 | } | |
122 | //--------------------------------------------------------------------------- | |
123 | void AliITSVertexerTracks::ComputeMaxChi2PerTrack(Int_t nTracks) { | |
124 | // | |
125 | // Max. contr. to the chi2 has been tuned as a function of multiplicity | |
126 | // | |
127 | if(nTracks < 7) { fMaxChi2PerTrack = 1.e6; | |
128 | } else { fMaxChi2PerTrack = 100.; } | |
129 | ||
130 | return; | |
131 | } | |
132 | //--------------------------------------------------------------------------- | |
133 | void AliITSVertexerTracks::FindVertices() { | |
134 | // | |
135 | // Vertices for all events from fFirstEvent to fLastEvent | |
136 | // | |
137 | ||
138 | // Check if the conv. const. has been set | |
139 | if(!CheckField()) return; | |
140 | ||
11ba84a4 | 141 | TDirectory *curdir = 0; |
cab6ff9b | 142 | |
143 | // loop over events | |
144 | for(Int_t ev=fFirstEvent; ev<=fLastEvent; ev++) { | |
145 | if(ev % 100 == 0 || fDebug) printf("--- Processing event %d of %d ---\n",ev,fLastEvent); | |
146 | ||
147 | FindVertexForCurrentEvent(ev); | |
148 | ||
149 | if(!fCurrentVertex) { | |
150 | printf("AliITSVertexerTracks::FindVertixes(): no tracks tree for event %d\n",ev); | |
151 | continue; | |
152 | } | |
153 | ||
154 | if(fDebug) fCurrentVertex->PrintStatus(); | |
11ba84a4 | 155 | |
156 | // write vertex to file | |
d4221964 | 157 | TString vtxName = "Vertex_"; |
158 | vtxName += ev; | |
11ba84a4 | 159 | fCurrentVertex->SetName(vtxName.Data()); |
d4221964 | 160 | fCurrentVertex->SetTitle("VertexerTracks"); |
11ba84a4 | 161 | //WriteCurrentVertex(); |
162 | curdir = gDirectory; | |
163 | fOutFile->cd(); | |
164 | fCurrentVertex->Write(); | |
165 | curdir->cd(); | |
166 | fCurrentVertex = 0; | |
167 | } // loop over events | |
168 | ||
169 | return; | |
170 | } | |
171 | //--------------------------------------------------------------------------- | |
172 | void AliITSVertexerTracks::FindVerticesESD() { | |
173 | // | |
174 | // Vertices for all events from fFirstEvent to fLastEvent | |
175 | // | |
176 | ||
177 | // Check if the conv. const. has been set | |
178 | if(!CheckField()) return; | |
179 | ||
180 | TDirectory *curdir = 0; | |
181 | ||
182 | fInFile->cd(); | |
183 | TKey *key=0; | |
184 | TIter next(fInFile->GetListOfKeys()); | |
185 | // loop on events in file | |
186 | while ((key=(TKey*)next())!=0) { | |
187 | AliESD *esdEvent=(AliESD*)key->ReadObj(); | |
188 | if(!esdEvent) { | |
189 | printf("AliITSVertexerTracks::FindVerticesESD(): not an ESD!\n"); | |
190 | return; | |
191 | } | |
192 | Int_t ev = (Int_t)esdEvent->GetEventNumber(); | |
193 | if(ev<fFirstEvent || ev>fLastEvent) { delete esdEvent; continue; } | |
194 | if(ev % 100 == 0 || fDebug) | |
195 | printf("--- Processing event %d of %d ---\n",ev,fLastEvent-fFirstEvent); | |
196 | ||
197 | FindVertexForCurrentEvent(esdEvent); | |
198 | ||
199 | if(!fCurrentVertex) { | |
200 | printf("AliITSVertexerTracks::FindVertixesESD():\n no vertex for event %d\n",ev); | |
201 | continue; | |
202 | } | |
203 | ||
204 | cout<<"VERTICE TROVATO\n"; | |
205 | if(fDebug) fCurrentVertex->PrintStatus(); | |
206 | ||
207 | // write the ESD to file | |
208 | curdir = gDirectory; | |
209 | Char_t ename[100]; | |
210 | sprintf(ename,"%d",ev); | |
211 | fOutFile->cd(); | |
212 | esdEvent->Dump(); | |
213 | esdEvent->Write(ename,TObject::kOverwrite); | |
214 | curdir->cd(); | |
215 | fCurrentVertex = 0; | |
cab6ff9b | 216 | } // loop over events |
217 | ||
218 | return; | |
219 | } | |
220 | //---------------------------------------------------------------------------- | |
221 | Int_t AliITSVertexerTracks::PrepareTracks(TTree &trkTree) { | |
222 | // | |
223 | // Propagate tracks to initial vertex position and store them in a TObjArray | |
224 | // | |
225 | Double_t maxd0rphi = 3.; | |
226 | Double_t alpha,xlStart,d0rphi; | |
227 | Int_t nTrks = 0; | |
228 | Bool_t skipThis; | |
229 | ||
230 | Int_t nEntries = (Int_t)trkTree.GetEntries(); | |
231 | ||
232 | if(!fTrkArray.IsEmpty()) fTrkArray.Clear(); | |
233 | fTrkArray.Expand(nEntries); | |
234 | ||
235 | if(fDebug) { | |
236 | printf(" PrepareTracks()\n"); | |
237 | trkTree.Print(); | |
238 | } | |
239 | ||
240 | for(Int_t i=0; i<nEntries; i++) { | |
241 | // check tracks to skip | |
242 | skipThis = kFALSE; | |
243 | for(Int_t j=0; j<fNTrksToSkip; j++) { | |
244 | if(i==fTrksToSkip[j]) { | |
245 | if(fDebug) printf("skipping track: %d\n",i); | |
246 | skipThis = kTRUE; | |
247 | } | |
248 | } | |
249 | if(skipThis) continue; | |
250 | ||
251 | AliITStrackV2 *itstrack = new AliITStrackV2; | |
252 | trkTree.SetBranchAddress("tracks",&itstrack); | |
253 | trkTree.GetEvent(i); | |
254 | ||
255 | ||
256 | // propagate track to vtxSeed | |
257 | alpha = itstrack->GetAlpha(); | |
258 | xlStart = fNominalPos[0]*TMath::Cos(alpha)+fNominalPos[1]*TMath::Sin(alpha); | |
259 | itstrack->PropagateTo(3.,0.0023,65.19); // to beam pipe (0.8 mm of Be) | |
260 | itstrack->PropagateTo(xlStart,0.,0.); // to vtxSeed | |
261 | ||
262 | // select tracks with d0rphi < maxd0rphi | |
263 | d0rphi = TMath::Abs(itstrack->GetD(fNominalPos[0],fNominalPos[1])); | |
264 | if(d0rphi > maxd0rphi) { delete itstrack; continue; } | |
265 | ||
266 | fTrkArray.AddLast(itstrack); | |
267 | ||
268 | nTrks++; | |
269 | } | |
270 | ||
271 | if(fTrksToSkip) delete [] fTrksToSkip; | |
272 | ||
273 | return nTrks; | |
274 | } | |
275 | //---------------------------------------------------------------------------- | |
276 | void AliITSVertexerTracks::PrintStatus() const { | |
277 | // | |
278 | // Print status | |
279 | // | |
280 | printf(" Initial position (%f,%f)\n",fNominalPos[0],fNominalPos[1]); | |
281 | printf(" Vertex position after vertex finder (%f, %f, %f)\n",fInitPos[0],fInitPos[1],fInitPos[2]); | |
282 | printf(" Number of tracks in array: %d\n",(Int_t)fTrkArray.GetEntriesFast()); | |
283 | printf(" Minimum # tracks required in fit: %d\n",fMinTracks); | |
cab6ff9b | 284 | |
285 | return; | |
286 | } | |
287 | //---------------------------------------------------------------------------- | |
d681bb2d | 288 | AliESDVertex* AliITSVertexerTracks::FindVertexForCurrentEvent(Int_t evnumb) { |
cab6ff9b | 289 | // |
290 | // Vertex for current event | |
291 | // | |
292 | fCurrentVertex = 0; | |
293 | ||
294 | // get tree with tracks from input file | |
295 | TString treeName = "TreeT_ITS_"; | |
296 | treeName += evnumb; | |
11ba84a4 | 297 | TTree *trkTree=(TTree*)fInFile->Get(treeName.Data()); |
cab6ff9b | 298 | if(!trkTree) return fCurrentVertex; |
299 | ||
300 | ||
301 | // get tracks and propagate them to initial vertex position | |
302 | Int_t nTrks = PrepareTracks(*trkTree); | |
303 | delete trkTree; | |
304 | if(fDebug) printf(" tracks prepared: %d\n",nTrks); | |
305 | if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; } | |
306 | ||
307 | // VERTEX FINDER | |
308 | VertexFinder(); | |
309 | ||
310 | // VERTEX FITTER | |
311 | ComputeMaxChi2PerTrack(nTrks); | |
cab6ff9b | 312 | VertexFitter(); |
313 | if(fDebug) printf(" vertex fit completed\n"); | |
314 | ||
cab6ff9b | 315 | return fCurrentVertex; |
316 | } | |
cab6ff9b | 317 | //---------------------------------------------------------------------------- |
d681bb2d | 318 | AliESDVertex* AliITSVertexerTracks::FindVertexForCurrentEvent(AliESD *esdEvent) |
11ba84a4 | 319 | { |
cab6ff9b | 320 | // |
11ba84a4 | 321 | // Vertex for current ESD event |
cab6ff9b | 322 | // |
11ba84a4 | 323 | fCurrentVertex = 0; |
324 | Double_t vtx[3],cvtx[6]; | |
325 | ||
326 | // put tracks reco in ITS in a tree | |
327 | Int_t entr = (Int_t)esdEvent->GetNumberOfTracks(); | |
328 | TTree *trkTree = new TTree("TreeT_ITS","its tracks"); | |
329 | AliITStrackV2 *itstrack = 0; | |
330 | trkTree->Branch("tracks","AliITStrackV2",&itstrack,entr,0); | |
331 | ||
332 | for(Int_t i=0; i<entr; i++) { | |
333 | AliESDtrack *esdTrack = (AliESDtrack*)esdEvent->GetTrack(i); | |
334 | if(!esdTrack->GetStatus()&AliESDtrack::kITSin) | |
335 | { delete esdTrack; continue; } | |
336 | itstrack = new AliITStrackV2(*esdTrack); | |
337 | trkTree->Fill(); | |
338 | itstrack = 0; | |
339 | delete esdTrack; | |
340 | } | |
341 | delete itstrack; | |
342 | ||
343 | // preselect tracks and propagate them to initial vertex position | |
344 | Int_t nTrks = PrepareTracks(*trkTree); | |
345 | delete trkTree; | |
346 | if(fDebug) printf(" tracks prepared: %d\n",nTrks); | |
347 | if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; } | |
cab6ff9b | 348 | |
11ba84a4 | 349 | // Set initial vertex position from ESD |
2257f27e | 350 | esdEvent->GetVertex()->GetXYZ(vtx); |
11ba84a4 | 351 | SetVtxStart(vtx[0],vtx[1]); |
cab6ff9b | 352 | |
11ba84a4 | 353 | // VERTEX FINDER |
354 | VertexFinder(); | |
cab6ff9b | 355 | |
11ba84a4 | 356 | // VERTEX FITTER |
357 | ComputeMaxChi2PerTrack(nTrks); | |
358 | VertexFitter(); | |
359 | if(fDebug) printf(" vertex fit completed\n"); | |
cab6ff9b | 360 | |
11ba84a4 | 361 | // store vertex information in ESD |
362 | fCurrentVertex->GetXYZ(vtx); | |
363 | fCurrentVertex->GetCovMatrix(cvtx); | |
2257f27e | 364 | esdEvent->SetVertex(fCurrentVertex); |
cab6ff9b | 365 | |
11ba84a4 | 366 | cout<<"Vertex: "<<vtx[0]<<", "<<vtx[1]<<", "<<vtx[2]<<endl; |
367 | return fCurrentVertex; | |
cab6ff9b | 368 | } |
369 | //--------------------------------------------------------------------------- | |
11ba84a4 | 370 | void AliITSVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped) { |
cab6ff9b | 371 | // |
11ba84a4 | 372 | // Mark the tracks not ot be used in the vertex finding |
cab6ff9b | 373 | // |
11ba84a4 | 374 | fNTrksToSkip = n; |
375 | fTrksToSkip = new Int_t[n]; | |
376 | for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i]; | |
377 | return; | |
cab6ff9b | 378 | } |
379 | //--------------------------------------------------------------------------- | |
380 | void AliITSVertexerTracks::TooFewTracks() { | |
381 | // | |
382 | // When the number of tracks is < fMinTracks the vertex is set to (0,0,0) | |
383 | // and the number of tracks to -1 | |
384 | // | |
d681bb2d | 385 | fCurrentVertex = new AliESDVertex(0.,0.,-1); |
cab6ff9b | 386 | return; |
387 | } | |
388 | //--------------------------------------------------------------------------- | |
389 | void AliITSVertexerTracks::VertexFinder() { | |
390 | ||
391 | // Get estimate of vertex position in (x,y) from tracks DCA | |
392 | // Then this estimate is stored to the data member fInitPos | |
393 | // (previous values are overwritten) | |
394 | ||
395 | ||
396 | ||
397 | /* | |
398 | ******* TEMPORARY!!! FOR TEST ONLY!!! ********************************** | |
399 | ||
400 | fInitPos[0] = fNominalPos[0]+gRandom->Gaus(0.,0.0100); // 100 micron gaussian smearing | |
401 | fInitPos[1] = fNominalPos[1]+gRandom->Gaus(0.,0.0100); // 100 micron gaussian smearing | |
402 | */ | |
403 | ||
404 | fInitPos[2] = 0.; | |
405 | for(Int_t i=0;i<2;i++)fInitPos[i]=fNominalPos[i]; | |
406 | ||
407 | Int_t nacc = (Int_t)fTrkArray.GetEntriesFast(); | |
408 | ||
409 | Double_t aver[3]={0.,0.,0.}; | |
410 | Int_t ncombi = 0; | |
411 | AliITStrackV2 *track1; | |
412 | AliITStrackV2 *track2; | |
413 | for(Int_t i=0; i<nacc; i++){ | |
414 | track1 = (AliITStrackV2*)fTrkArray.At(i); | |
415 | if(fDebug>5){ | |
416 | Double_t xv,par[5]; | |
417 | track1->GetExternalParameters(xv,par); | |
418 | cout<<"Track in position "<<i<<" xr= "<<xv<<endl; | |
419 | for(Int_t ii=0;ii<5;ii++)cout<<par[ii]<<" "; | |
420 | cout<<endl; | |
421 | } | |
422 | Double_t mom1[3]; | |
423 | Double_t alpha = track1->GetAlpha(); | |
424 | Double_t azim = TMath::ASin(track1->GetSnp())+alpha; | |
425 | Double_t theta = TMath::Pi()/2. - TMath::ATan(track1->GetTgl()); | |
426 | mom1[0] = TMath::Sin(theta)*TMath::Cos(azim); | |
427 | mom1[1] = TMath::Sin(theta)*TMath::Sin(azim); | |
428 | mom1[2] = TMath::Cos(theta); | |
429 | ||
430 | Double_t pos1[3]; | |
431 | Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1]; | |
432 | track1->GetGlobalXYZat(mindist,pos1[0],pos1[1],pos1[2]); | |
433 | AliITSStrLine *line1 = new AliITSStrLine(pos1,mom1); | |
434 | for(Int_t j=i+1; j<nacc; j++){ | |
435 | track2 = (AliITStrackV2*)fTrkArray.At(j); | |
436 | Double_t mom2[3]; | |
437 | alpha = track2->GetAlpha(); | |
438 | azim = TMath::ASin(track2->GetSnp())+alpha; | |
439 | theta = TMath::Pi()/2. - TMath::ATan(track2->GetTgl()); | |
440 | mom2[0] = TMath::Sin(theta)*TMath::Cos(azim); | |
441 | mom2[1] = TMath::Sin(theta)*TMath::Sin(azim); | |
442 | mom2[2] = TMath::Cos(theta); | |
443 | Double_t pos2[3]; | |
444 | mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1]; | |
445 | track2->GetGlobalXYZat(mindist,pos2[0],pos2[1],pos2[2]); | |
446 | AliITSStrLine *line2 = new AliITSStrLine(pos2,mom2); | |
447 | Double_t crosspoint[3]; | |
448 | Int_t retcode = line2->Cross(line1,crosspoint); | |
449 | if(retcode<0){ | |
450 | if(fDebug>10)cout<<" i= "<<i<<", j= "<<j<<endl; | |
451 | if(fDebug>10)cout<<"bad intersection\n"; | |
452 | line1->PrintStatus(); | |
453 | line2->PrintStatus(); | |
454 | } | |
455 | else { | |
456 | ncombi++; | |
457 | for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj]; | |
458 | if(fDebug>10)cout<<" i= "<<i<<", j= "<<j<<endl; | |
459 | if(fDebug>10)cout<<"\n Cross point: "; | |
460 | if(fDebug>10)cout<<crosspoint[0]<<" "<<crosspoint[1]<<" "<<crosspoint[2]<<endl; | |
461 | } | |
462 | delete line2; | |
463 | } | |
464 | delete line1; | |
465 | } | |
466 | if(ncombi>0){ | |
467 | for(Int_t jj=0;jj<3;jj++)fInitPos[jj] = aver[jj]/ncombi; | |
468 | } | |
469 | else { | |
470 | Warning("VertexFinder","Finder did not succed"); | |
471 | } | |
472 | ||
473 | ||
474 | //************************************************************************ | |
475 | return; | |
476 | ||
477 | } | |
cab6ff9b | 478 | //--------------------------------------------------------------------------- |
479 | void AliITSVertexerTracks::VertexFitter() { | |
480 | // | |
481 | // The optimal estimate of the vertex position is given by a "weighted | |
482 | // average of tracks positions" | |
483 | // Original method: CMS Note 97/0051 | |
484 | // | |
485 | if(fDebug) { | |
486 | printf(" VertexFitter(): start\n"); | |
487 | PrintStatus(); | |
488 | } | |
489 | ||
490 | ||
491 | Int_t i,j,k,step=0; | |
492 | TMatrixD rv(3,1); | |
493 | TMatrixD V(3,3); | |
494 | rv(0,0) = fInitPos[0]; | |
495 | rv(1,0) = fInitPos[1]; | |
496 | rv(2,0) = 0.; | |
497 | Double_t xlStart,alpha; | |
498 | Double_t rotAngle; | |
499 | Double_t cosRot,sinRot; | |
500 | Double_t cc[15]; | |
501 | Int_t nUsedTrks; | |
502 | Double_t chi2,chi2i; | |
503 | Int_t arrEntries = (Int_t)fTrkArray.GetEntries(); | |
504 | AliITStrackV2 *t = 0; | |
505 | Int_t failed = 0; | |
506 | ||
507 | Int_t *skipTrack = new Int_t[arrEntries]; | |
508 | for(i=0; i<arrEntries; i++) skipTrack[i]=0; | |
509 | ||
510 | // 3 steps: | |
511 | // 1st - first estimate of vtx using all tracks | |
512 | // 2nd - apply cut on chi2 max per track | |
513 | // 3rd - estimate of global chi2 | |
514 | for(step=0; step<3; step++) { | |
515 | if(fDebug) printf(" step = %d\n",step); | |
516 | chi2 = 0.; | |
517 | nUsedTrks = 0; | |
518 | ||
519 | TMatrixD SumWiri(3,1); | |
520 | TMatrixD SumWi(3,3); | |
521 | for(i=0; i<3; i++) { | |
522 | SumWiri(i,0) = 0.; | |
523 | for(j=0; j<3; j++) SumWi(j,i) = 0.; | |
524 | } | |
525 | ||
526 | // loop on tracks | |
527 | for(k=0; k<arrEntries; k++) { | |
528 | if(skipTrack[k]) continue; | |
529 | // get track from track array | |
530 | t = (AliITStrackV2*)fTrkArray.At(k); | |
531 | alpha = t->GetAlpha(); | |
532 | xlStart = fInitPos[0]*TMath::Cos(alpha)+fInitPos[1]*TMath::Sin(alpha); | |
533 | t->PropagateTo(xlStart,0.,0.); // to vtxSeed | |
11ba84a4 | 534 | rotAngle = alpha; |
cab6ff9b | 535 | if(alpha<0.) rotAngle += 2.*TMath::Pi(); |
536 | cosRot = TMath::Cos(rotAngle); | |
537 | sinRot = TMath::Sin(rotAngle); | |
538 | ||
539 | // vector of track global coordinates | |
540 | TMatrixD ri(3,1); | |
541 | ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot; | |
542 | ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot; | |
543 | ri(2,0) = t->GetZ(); | |
544 | ||
545 | // matrix to go from global (x,y,z) to local (y,z); | |
546 | TMatrixD Qi(2,3); | |
547 | Qi(0,0) = -sinRot; | |
548 | Qi(0,1) = cosRot; | |
549 | Qi(0,2) = 0.; | |
550 | Qi(1,0) = 0.; | |
551 | Qi(1,1) = 0.; | |
552 | Qi(1,2) = 1.; | |
553 | ||
554 | // covariance matrix of local (y,z) - inverted | |
555 | TMatrixD Ui(2,2); | |
556 | t->GetExternalCovariance(cc); | |
557 | Ui(0,0) = cc[0]; | |
558 | Ui(0,1) = cc[1]; | |
559 | Ui(1,0) = cc[1]; | |
560 | Ui(1,1) = cc[2]; | |
561 | ||
562 | // weights matrix: Wi = QiT * UiInv * Qi | |
563 | if(Ui.Determinant() <= 0.) continue; | |
564 | TMatrixD UiInv(TMatrixD::kInverted,Ui); | |
565 | TMatrixD UiInvQi(UiInv,TMatrixD::kMult,Qi); | |
566 | TMatrixD Wi(Qi,TMatrixD::kTransposeMult,UiInvQi); | |
567 | ||
568 | // track chi2 | |
569 | TMatrixD deltar = rv; deltar -= ri; | |
570 | TMatrixD Wideltar(Wi,TMatrixD::kMult,deltar); | |
571 | chi2i = deltar(0,0)*Wideltar(0,0)+ | |
572 | deltar(1,0)*Wideltar(1,0)+ | |
573 | deltar(2,0)*Wideltar(2,0); | |
574 | ||
575 | ||
576 | if(step==1 && chi2i > fMaxChi2PerTrack) { | |
577 | skipTrack[k] = 1; | |
578 | continue; | |
579 | } | |
580 | ||
581 | // add to total chi2 | |
582 | chi2 += chi2i; | |
583 | ||
584 | TMatrixD Wiri(Wi,TMatrixD::kMult,ri); | |
585 | ||
586 | SumWiri += Wiri; | |
587 | SumWi += Wi; | |
588 | ||
589 | nUsedTrks++; | |
590 | } // end loop on tracks | |
591 | ||
592 | if(nUsedTrks < fMinTracks) { | |
593 | failed=1; | |
594 | continue; | |
595 | } | |
596 | ||
597 | Double_t determinant = SumWi.Determinant(); | |
598 | //cerr<<" determinant: "<<determinant<<endl; | |
599 | if(determinant < 100.) { | |
600 | printf("det(V) = 0\n"); | |
601 | failed=1; | |
602 | continue; | |
603 | } | |
604 | ||
605 | // inverted of weights matrix | |
606 | TMatrixD InvSumWi(TMatrixD::kInverted,SumWi); | |
607 | V = InvSumWi; | |
608 | ||
609 | // position of primary vertex | |
610 | rv.Mult(V,SumWiri); | |
611 | ||
612 | } // end loop on the 3 steps | |
613 | ||
614 | delete [] skipTrack; | |
615 | delete t; | |
616 | ||
617 | if(failed) { | |
618 | TooFewTracks(); | |
619 | return; | |
620 | } | |
621 | ||
622 | Double_t position[3]; | |
623 | position[0] = rv(0,0); | |
624 | position[1] = rv(1,0); | |
625 | position[2] = rv(2,0); | |
626 | Double_t covmatrix[6]; | |
627 | covmatrix[0] = V(0,0); | |
628 | covmatrix[1] = V(0,1); | |
629 | covmatrix[2] = V(1,1); | |
630 | covmatrix[3] = V(0,2); | |
631 | covmatrix[4] = V(1,2); | |
632 | covmatrix[5] = V(2,2); | |
633 | ||
634 | // store data in the vertex object | |
d681bb2d | 635 | fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nUsedTrks); |
cab6ff9b | 636 | |
637 | if(fDebug) { | |
638 | printf(" VertexFitter(): finish\n"); | |
639 | printf(" rv = ( %f , %f , %f )\n\n",rv(0,0),rv(1,0),rv(2,0)); | |
640 | fCurrentVertex->PrintStatus(); | |
641 | } | |
642 | ||
643 | return; | |
644 | } | |
645 | //---------------------------------------------------------------------------- | |
d681bb2d | 646 | AliESDVertex *AliITSVertexerTracks::VertexOnTheFly(TTree &trkTree) { |
cab6ff9b | 647 | // |
648 | // Return vertex from tracks in trkTree | |
649 | // | |
650 | if(fCurrentVertex) fCurrentVertex = 0; | |
651 | ||
652 | // get tracks and propagate them to initial vertex position | |
653 | Int_t nTrks = PrepareTracks(*(&trkTree)); | |
654 | if(fDebug) printf(" tracks prepared: %d\n",nTrks); | |
655 | if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; } | |
656 | ||
657 | // VERTEX FINDER | |
658 | VertexFinder(); | |
659 | ||
660 | // VERTEX FITTER | |
661 | ComputeMaxChi2PerTrack(nTrks); | |
cab6ff9b | 662 | VertexFitter(); |
663 | if(fDebug) printf(" vertex fit completed\n"); | |
664 | ||
665 | return fCurrentVertex; | |
666 | } | |
667 | //---------------------------------------------------------------------------- | |
668 | ||
669 | ||
670 |