]>
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; } | |
70d3d424 | 350 | try { |
351 | itstrack = new AliITStrackV2(*esdTrack); | |
352 | } | |
353 | catch (const Char_t *msg) { | |
354 | Warning("FindVertexForCurrentEvent",msg); | |
355 | delete esdTrack; | |
356 | continue; | |
357 | } | |
358 | ||
11ba84a4 | 359 | trkTree->Fill(); |
360 | itstrack = 0; | |
361 | delete esdTrack; | |
362 | } | |
363 | delete itstrack; | |
364 | ||
365 | // preselect tracks and propagate them to initial vertex position | |
366 | Int_t nTrks = PrepareTracks(*trkTree); | |
367 | delete trkTree; | |
368 | if(fDebug) printf(" tracks prepared: %d\n",nTrks); | |
369 | if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; } | |
cab6ff9b | 370 | |
11ba84a4 | 371 | // Set initial vertex position from ESD |
2257f27e | 372 | esdEvent->GetVertex()->GetXYZ(vtx); |
11ba84a4 | 373 | SetVtxStart(vtx[0],vtx[1]); |
cab6ff9b | 374 | |
11ba84a4 | 375 | // VERTEX FINDER |
376 | VertexFinder(); | |
cab6ff9b | 377 | |
11ba84a4 | 378 | // VERTEX FITTER |
379 | ComputeMaxChi2PerTrack(nTrks); | |
380 | VertexFitter(); | |
381 | if(fDebug) printf(" vertex fit completed\n"); | |
cab6ff9b | 382 | |
11ba84a4 | 383 | // store vertex information in ESD |
384 | fCurrentVertex->GetXYZ(vtx); | |
385 | fCurrentVertex->GetCovMatrix(cvtx); | |
70d3d424 | 386 | |
387 | Double_t tp[3]; | |
388 | esdEvent->GetVertex()->GetTruePos(tp); | |
389 | fCurrentVertex->SetTruePos(tp); | |
390 | ||
2257f27e | 391 | esdEvent->SetVertex(fCurrentVertex); |
cab6ff9b | 392 | |
11ba84a4 | 393 | cout<<"Vertex: "<<vtx[0]<<", "<<vtx[1]<<", "<<vtx[2]<<endl; |
394 | return fCurrentVertex; | |
cab6ff9b | 395 | } |
396 | //--------------------------------------------------------------------------- | |
11ba84a4 | 397 | void AliITSVertexerTracks::SetSkipTracks(Int_t n,Int_t *skipped) { |
cab6ff9b | 398 | // |
11ba84a4 | 399 | // Mark the tracks not ot be used in the vertex finding |
cab6ff9b | 400 | // |
11ba84a4 | 401 | fNTrksToSkip = n; |
402 | fTrksToSkip = new Int_t[n]; | |
403 | for(Int_t i=0;i<n;i++) fTrksToSkip[i] = skipped[i]; | |
404 | return; | |
cab6ff9b | 405 | } |
406 | //--------------------------------------------------------------------------- | |
407 | void AliITSVertexerTracks::TooFewTracks() { | |
408 | // | |
409 | // When the number of tracks is < fMinTracks the vertex is set to (0,0,0) | |
410 | // and the number of tracks to -1 | |
411 | // | |
d681bb2d | 412 | fCurrentVertex = new AliESDVertex(0.,0.,-1); |
cab6ff9b | 413 | return; |
414 | } | |
415 | //--------------------------------------------------------------------------- | |
416 | void AliITSVertexerTracks::VertexFinder() { | |
417 | ||
418 | // Get estimate of vertex position in (x,y) from tracks DCA | |
419 | // Then this estimate is stored to the data member fInitPos | |
420 | // (previous values are overwritten) | |
421 | ||
422 | ||
423 | ||
424 | /* | |
425 | ******* TEMPORARY!!! FOR TEST ONLY!!! ********************************** | |
426 | ||
427 | fInitPos[0] = fNominalPos[0]+gRandom->Gaus(0.,0.0100); // 100 micron gaussian smearing | |
428 | fInitPos[1] = fNominalPos[1]+gRandom->Gaus(0.,0.0100); // 100 micron gaussian smearing | |
429 | */ | |
430 | ||
431 | fInitPos[2] = 0.; | |
432 | for(Int_t i=0;i<2;i++)fInitPos[i]=fNominalPos[i]; | |
433 | ||
434 | Int_t nacc = (Int_t)fTrkArray.GetEntriesFast(); | |
435 | ||
436 | Double_t aver[3]={0.,0.,0.}; | |
437 | Int_t ncombi = 0; | |
438 | AliITStrackV2 *track1; | |
439 | AliITStrackV2 *track2; | |
96bd751d | 440 | AliITSStrLine line1; |
441 | AliITSStrLine line2; | |
cab6ff9b | 442 | for(Int_t i=0; i<nacc; i++){ |
443 | track1 = (AliITStrackV2*)fTrkArray.At(i); | |
444 | if(fDebug>5){ | |
445 | Double_t xv,par[5]; | |
446 | track1->GetExternalParameters(xv,par); | |
447 | cout<<"Track in position "<<i<<" xr= "<<xv<<endl; | |
448 | for(Int_t ii=0;ii<5;ii++)cout<<par[ii]<<" "; | |
449 | cout<<endl; | |
450 | } | |
451 | Double_t mom1[3]; | |
452 | Double_t alpha = track1->GetAlpha(); | |
453 | Double_t azim = TMath::ASin(track1->GetSnp())+alpha; | |
454 | Double_t theta = TMath::Pi()/2. - TMath::ATan(track1->GetTgl()); | |
455 | mom1[0] = TMath::Sin(theta)*TMath::Cos(azim); | |
456 | mom1[1] = TMath::Sin(theta)*TMath::Sin(azim); | |
457 | mom1[2] = TMath::Cos(theta); | |
458 | ||
459 | Double_t pos1[3]; | |
460 | Double_t mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1]; | |
461 | track1->GetGlobalXYZat(mindist,pos1[0],pos1[1],pos1[2]); | |
96bd751d | 462 | line1.SetP0(pos1); |
463 | line1.SetCd(mom1); | |
cab6ff9b | 464 | for(Int_t j=i+1; j<nacc; j++){ |
465 | track2 = (AliITStrackV2*)fTrkArray.At(j); | |
466 | Double_t mom2[3]; | |
467 | alpha = track2->GetAlpha(); | |
468 | azim = TMath::ASin(track2->GetSnp())+alpha; | |
469 | theta = TMath::Pi()/2. - TMath::ATan(track2->GetTgl()); | |
470 | mom2[0] = TMath::Sin(theta)*TMath::Cos(azim); | |
471 | mom2[1] = TMath::Sin(theta)*TMath::Sin(azim); | |
472 | mom2[2] = TMath::Cos(theta); | |
473 | Double_t pos2[3]; | |
474 | mindist = TMath::Cos(alpha)*fNominalPos[0]+TMath::Sin(alpha)*fNominalPos[1]; | |
475 | track2->GetGlobalXYZat(mindist,pos2[0],pos2[1],pos2[2]); | |
96bd751d | 476 | line2.SetP0(pos2); |
477 | line2.SetCd(mom2); | |
cab6ff9b | 478 | Double_t crosspoint[3]; |
96bd751d | 479 | Int_t retcode = line2.Cross(&line1,crosspoint); |
cab6ff9b | 480 | if(retcode<0){ |
481 | if(fDebug>10)cout<<" i= "<<i<<", j= "<<j<<endl; | |
482 | if(fDebug>10)cout<<"bad intersection\n"; | |
96bd751d | 483 | line1.PrintStatus(); |
484 | line2.PrintStatus(); | |
cab6ff9b | 485 | } |
486 | else { | |
487 | ncombi++; | |
488 | for(Int_t jj=0;jj<3;jj++)aver[jj]+=crosspoint[jj]; | |
489 | if(fDebug>10)cout<<" i= "<<i<<", j= "<<j<<endl; | |
490 | if(fDebug>10)cout<<"\n Cross point: "; | |
491 | if(fDebug>10)cout<<crosspoint[0]<<" "<<crosspoint[1]<<" "<<crosspoint[2]<<endl; | |
492 | } | |
cab6ff9b | 493 | } |
cab6ff9b | 494 | } |
495 | if(ncombi>0){ | |
496 | for(Int_t jj=0;jj<3;jj++)fInitPos[jj] = aver[jj]/ncombi; | |
497 | } | |
498 | else { | |
499 | Warning("VertexFinder","Finder did not succed"); | |
500 | } | |
501 | ||
502 | ||
503 | //************************************************************************ | |
504 | return; | |
505 | ||
506 | } | |
cab6ff9b | 507 | //--------------------------------------------------------------------------- |
508 | void AliITSVertexerTracks::VertexFitter() { | |
509 | // | |
510 | // The optimal estimate of the vertex position is given by a "weighted | |
511 | // average of tracks positions" | |
512 | // Original method: CMS Note 97/0051 | |
513 | // | |
514 | if(fDebug) { | |
515 | printf(" VertexFitter(): start\n"); | |
516 | PrintStatus(); | |
517 | } | |
518 | ||
519 | ||
520 | Int_t i,j,k,step=0; | |
521 | TMatrixD rv(3,1); | |
a48d1ed3 | 522 | TMatrixD vV(3,3); |
cab6ff9b | 523 | rv(0,0) = fInitPos[0]; |
524 | rv(1,0) = fInitPos[1]; | |
525 | rv(2,0) = 0.; | |
526 | Double_t xlStart,alpha; | |
527 | Double_t rotAngle; | |
528 | Double_t cosRot,sinRot; | |
529 | Double_t cc[15]; | |
530 | Int_t nUsedTrks; | |
531 | Double_t chi2,chi2i; | |
532 | Int_t arrEntries = (Int_t)fTrkArray.GetEntries(); | |
533 | AliITStrackV2 *t = 0; | |
534 | Int_t failed = 0; | |
535 | ||
536 | Int_t *skipTrack = new Int_t[arrEntries]; | |
537 | for(i=0; i<arrEntries; i++) skipTrack[i]=0; | |
538 | ||
539 | // 3 steps: | |
540 | // 1st - first estimate of vtx using all tracks | |
541 | // 2nd - apply cut on chi2 max per track | |
542 | // 3rd - estimate of global chi2 | |
543 | for(step=0; step<3; step++) { | |
544 | if(fDebug) printf(" step = %d\n",step); | |
545 | chi2 = 0.; | |
546 | nUsedTrks = 0; | |
547 | ||
a48d1ed3 | 548 | TMatrixD sumWiri(3,1); |
549 | TMatrixD sumWi(3,3); | |
cab6ff9b | 550 | for(i=0; i<3; i++) { |
a48d1ed3 | 551 | sumWiri(i,0) = 0.; |
552 | for(j=0; j<3; j++) sumWi(j,i) = 0.; | |
cab6ff9b | 553 | } |
554 | ||
555 | // loop on tracks | |
556 | for(k=0; k<arrEntries; k++) { | |
557 | if(skipTrack[k]) continue; | |
558 | // get track from track array | |
559 | t = (AliITStrackV2*)fTrkArray.At(k); | |
560 | alpha = t->GetAlpha(); | |
561 | xlStart = fInitPos[0]*TMath::Cos(alpha)+fInitPos[1]*TMath::Sin(alpha); | |
562 | t->PropagateTo(xlStart,0.,0.); // to vtxSeed | |
11ba84a4 | 563 | rotAngle = alpha; |
cab6ff9b | 564 | if(alpha<0.) rotAngle += 2.*TMath::Pi(); |
565 | cosRot = TMath::Cos(rotAngle); | |
566 | sinRot = TMath::Sin(rotAngle); | |
567 | ||
568 | // vector of track global coordinates | |
569 | TMatrixD ri(3,1); | |
570 | ri(0,0) = t->GetX()*cosRot-t->GetY()*sinRot; | |
571 | ri(1,0) = t->GetX()*sinRot+t->GetY()*cosRot; | |
572 | ri(2,0) = t->GetZ(); | |
573 | ||
574 | // matrix to go from global (x,y,z) to local (y,z); | |
a48d1ed3 | 575 | TMatrixD qQi(2,3); |
576 | qQi(0,0) = -sinRot; | |
577 | qQi(0,1) = cosRot; | |
578 | qQi(0,2) = 0.; | |
579 | qQi(1,0) = 0.; | |
580 | qQi(1,1) = 0.; | |
581 | qQi(1,2) = 1.; | |
cab6ff9b | 582 | |
583 | // covariance matrix of local (y,z) - inverted | |
a48d1ed3 | 584 | TMatrixD uUi(2,2); |
cab6ff9b | 585 | t->GetExternalCovariance(cc); |
a48d1ed3 | 586 | uUi(0,0) = cc[0]; |
587 | uUi(0,1) = cc[1]; | |
588 | uUi(1,0) = cc[1]; | |
589 | uUi(1,1) = cc[2]; | |
cab6ff9b | 590 | |
a48d1ed3 | 591 | // weights matrix: wWi = qQiT * uUiInv * qQi |
592 | if(uUi.Determinant() <= 0.) continue; | |
593 | TMatrixD uUiInv(TMatrixD::kInverted,uUi); | |
594 | TMatrixD uUiInvQi(uUiInv,TMatrixD::kMult,qQi); | |
595 | TMatrixD wWi(qQi,TMatrixD::kTransposeMult,uUiInvQi); | |
cab6ff9b | 596 | |
597 | // track chi2 | |
598 | TMatrixD deltar = rv; deltar -= ri; | |
a48d1ed3 | 599 | TMatrixD wWideltar(wWi,TMatrixD::kMult,deltar); |
600 | chi2i = deltar(0,0)*wWideltar(0,0)+ | |
601 | deltar(1,0)*wWideltar(1,0)+ | |
602 | deltar(2,0)*wWideltar(2,0); | |
cab6ff9b | 603 | |
604 | ||
605 | if(step==1 && chi2i > fMaxChi2PerTrack) { | |
606 | skipTrack[k] = 1; | |
607 | continue; | |
608 | } | |
609 | ||
610 | // add to total chi2 | |
611 | chi2 += chi2i; | |
612 | ||
a48d1ed3 | 613 | TMatrixD wWiri(wWi,TMatrixD::kMult,ri); |
cab6ff9b | 614 | |
a48d1ed3 | 615 | sumWiri += wWiri; |
616 | sumWi += wWi; | |
cab6ff9b | 617 | |
618 | nUsedTrks++; | |
619 | } // end loop on tracks | |
620 | ||
621 | if(nUsedTrks < fMinTracks) { | |
622 | failed=1; | |
623 | continue; | |
624 | } | |
625 | ||
a48d1ed3 | 626 | Double_t determinant = sumWi.Determinant(); |
cab6ff9b | 627 | //cerr<<" determinant: "<<determinant<<endl; |
628 | if(determinant < 100.) { | |
629 | printf("det(V) = 0\n"); | |
630 | failed=1; | |
631 | continue; | |
632 | } | |
633 | ||
634 | // inverted of weights matrix | |
a48d1ed3 | 635 | TMatrixD invsumWi(TMatrixD::kInverted,sumWi); |
636 | vV = invsumWi; | |
cab6ff9b | 637 | |
638 | // position of primary vertex | |
a48d1ed3 | 639 | rv.Mult(vV,sumWiri); |
cab6ff9b | 640 | |
641 | } // end loop on the 3 steps | |
642 | ||
643 | delete [] skipTrack; | |
644 | delete t; | |
645 | ||
646 | if(failed) { | |
647 | TooFewTracks(); | |
648 | return; | |
649 | } | |
650 | ||
651 | Double_t position[3]; | |
652 | position[0] = rv(0,0); | |
653 | position[1] = rv(1,0); | |
654 | position[2] = rv(2,0); | |
655 | Double_t covmatrix[6]; | |
a48d1ed3 | 656 | covmatrix[0] = vV(0,0); |
657 | covmatrix[1] = vV(0,1); | |
658 | covmatrix[2] = vV(1,1); | |
659 | covmatrix[3] = vV(0,2); | |
660 | covmatrix[4] = vV(1,2); | |
661 | covmatrix[5] = vV(2,2); | |
cab6ff9b | 662 | |
663 | // store data in the vertex object | |
d681bb2d | 664 | fCurrentVertex = new AliESDVertex(position,covmatrix,chi2,nUsedTrks); |
cab6ff9b | 665 | |
666 | if(fDebug) { | |
667 | printf(" VertexFitter(): finish\n"); | |
668 | printf(" rv = ( %f , %f , %f )\n\n",rv(0,0),rv(1,0),rv(2,0)); | |
669 | fCurrentVertex->PrintStatus(); | |
670 | } | |
671 | ||
672 | return; | |
673 | } | |
674 | //---------------------------------------------------------------------------- | |
d681bb2d | 675 | AliESDVertex *AliITSVertexerTracks::VertexOnTheFly(TTree &trkTree) { |
cab6ff9b | 676 | // |
677 | // Return vertex from tracks in trkTree | |
678 | // | |
679 | if(fCurrentVertex) fCurrentVertex = 0; | |
680 | ||
681 | // get tracks and propagate them to initial vertex position | |
682 | Int_t nTrks = PrepareTracks(*(&trkTree)); | |
683 | if(fDebug) printf(" tracks prepared: %d\n",nTrks); | |
684 | if(nTrks < fMinTracks) { TooFewTracks(); return fCurrentVertex; } | |
685 | ||
686 | // VERTEX FINDER | |
687 | VertexFinder(); | |
688 | ||
689 | // VERTEX FITTER | |
690 | ComputeMaxChi2PerTrack(nTrks); | |
cab6ff9b | 691 | VertexFitter(); |
692 | if(fDebug) printf(" vertex fit completed\n"); | |
693 | ||
694 | return fCurrentVertex; | |
695 | } | |
696 | //---------------------------------------------------------------------------- | |
697 | ||
698 | ||
699 |