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