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