]>
Commit | Line | Data |
---|---|---|
1 | ||
2 | // Main authors: Matevz Tadel & Alja Mrak-Tadel & Bogdan Vulpescu: 2006, 2007 | |
3 | ||
4 | /************************************************************************** | |
5 | * Copyright(c) 1998-2008, ALICE Experiment at CERN, all rights reserved. * | |
6 | * See http://aliceinfo.cern.ch/Offline/AliRoot/License.html for * | |
7 | * full copyright notice. * | |
8 | **************************************************************************/ | |
9 | ||
10 | #include "AliEveMUONTrack.h" | |
11 | ||
12 | #include <EveBase/AliEveEventManager.h> | |
13 | ||
14 | #include <AliMagF.h> | |
15 | #include <AliLog.h> | |
16 | #include <AliESDMuonTrack.h> | |
17 | #include <AliESDEvent.h> | |
18 | #include <AliESDVertex.h> | |
19 | #include <AliRunLoader.h> | |
20 | #include <AliRun.h> | |
21 | ||
22 | #include <AliMUONTrack.h> | |
23 | #include <AliMUONTriggerTrack.h> | |
24 | #include <AliMUONTrackParam.h> | |
25 | #include <AliMUONConstants.h> | |
26 | #include <AliMUONESDInterface.h> | |
27 | #include <AliMUONVCluster.h> | |
28 | ||
29 | #include <TClonesArray.h> | |
30 | #include <TGeoGlobalMagField.h> | |
31 | #include <TParticle.h> | |
32 | #include <TParticlePDG.h> | |
33 | #include <TROOT.h> | |
34 | #include <TStyle.h> | |
35 | ||
36 | #include <Riostream.h> | |
37 | ||
38 | //============================================================================== | |
39 | // AliEveMUONTrack | |
40 | //============================================================================== | |
41 | //============================================================================== | |
42 | ||
43 | //______________________________________________________________________________ | |
44 | // Produce TEveUtil:TEveTrack from AliMUONTrack with dipole field model | |
45 | ||
46 | ClassImp(AliEveMUONTrack) | |
47 | ||
48 | ||
49 | //______________________________________________________________________________ | |
50 | AliEveMUONTrack::AliEveMUONTrack(TEveRecTrack* t, TEveTrackPropagator* rs) : | |
51 | TEveTrack(t,rs), | |
52 | fTrack(0), | |
53 | fPart(0), | |
54 | fCount(0), | |
55 | fIsMUONTrack(kFALSE), | |
56 | fIsMUONTriggerTrack(kFALSE), | |
57 | fIsESDTrack(kFALSE), | |
58 | fIsMCTrack(kFALSE), | |
59 | fIsRefTrack(kFALSE) | |
60 | { | |
61 | // | |
62 | // constructor | |
63 | // | |
64 | } | |
65 | ||
66 | //______________________________________________________________________________ | |
67 | AliEveMUONTrack::~AliEveMUONTrack() | |
68 | { | |
69 | // | |
70 | // destructor | |
71 | // | |
72 | ||
73 | if (fIsRefTrack || fIsESDTrack) delete fTrack; | |
74 | if (fIsMCTrack) delete fPart; | |
75 | ||
76 | } | |
77 | ||
78 | //______________________________________________________________________________ | |
79 | void AliEveMUONTrack::PrintMCTrackInfo() | |
80 | { | |
81 | // | |
82 | // information about the MC particle | |
83 | // | |
84 | ||
85 | Float_t pt, p; | |
86 | ||
87 | if (!fPart) { | |
88 | cout << " ! no particle ..." << endl; | |
89 | return; | |
90 | } | |
91 | ||
92 | cout << endl; | |
93 | cout << " MC track parameters at vertex" << endl; | |
94 | cout << " -------------------------------------------------------------------------------------" << endl; | |
95 | cout << " PDG code Vx Vy Vz Px Py Pz " << endl; | |
96 | ||
97 | cout << " " << | |
98 | setw(8) << setprecision(0) << | |
99 | fPart->GetPdgCode() << " " << | |
100 | setw(8) << setprecision(3) << | |
101 | fPart->Vx() << " " << | |
102 | setw(8) << setprecision(3) << | |
103 | fPart->Vy() << " " << | |
104 | setw(8) << setprecision(3) << | |
105 | fPart->Vz() << " " << | |
106 | setw(8) << setprecision(3) << | |
107 | fPart->Px() << " " << | |
108 | setw(8) << setprecision(3) << | |
109 | fPart->Py() << " " << | |
110 | setw(8) << setprecision(4) << | |
111 | fPart->Pz() << " " << | |
112 | ||
113 | endl; | |
114 | ||
115 | pt = TMath::Sqrt(fPart->Px()*fPart->Px()+fPart->Py()*fPart->Py()); | |
116 | p = TMath::Sqrt(fPart->Px()*fPart->Px()+fPart->Py()*fPart->Py()+fPart->Pz()*fPart->Pz()); | |
117 | ||
118 | cout << endl; | |
119 | cout << " Pt = " << | |
120 | setw(8) << setprecision(3) << | |
121 | pt << " GeV/c" << endl; | |
122 | ||
123 | cout << " P = " << | |
124 | setw(8) << setprecision(4) << | |
125 | p << " GeV/c" << endl; | |
126 | ||
127 | } | |
128 | ||
129 | //______________________________________________________________________________ | |
130 | void AliEveMUONTrack::PrintMUONTrackInfo() | |
131 | { | |
132 | // | |
133 | // information about the reconstructed/reference track; at hits and at vertex | |
134 | // | |
135 | ||
136 | Double_t radDeg = 180.0/TMath::Pi(); | |
137 | ||
138 | Int_t nparam; | |
139 | Float_t pt, bc, nbc, zc; | |
140 | AliMUONTrackParam *mtp; | |
141 | TClonesArray *trackParamAtCluster; | |
142 | ||
143 | if (!fTrack) { | |
144 | cout << " ! no reconstructed track ..." << endl; | |
145 | return; | |
146 | } | |
147 | ||
148 | if (fIsMUONTrack) { | |
149 | cout << endl; | |
150 | cout << " TEveTrack number " << fLabel << endl; | |
151 | cout << " ---------------------------------------------------------------------------------------------------------------------------------" << endl; | |
152 | cout << endl; | |
153 | cout << " Number of clusters " << fTrack->GetNClusters() << endl; | |
154 | cout << " Match to trigger " << fTrack->GetMatchTrigger() << endl; | |
155 | if (fTrack->GetMatchTrigger()) { | |
156 | cout << " Chi2 tracking-trigger " << fTrack->GetChi2MatchTrigger() << endl; | |
157 | cout << " Local trigger number " << fTrack->GetLoTrgNum() << endl; | |
158 | } | |
159 | } | |
160 | ||
161 | if (fIsRefTrack) { | |
162 | cout << endl; | |
163 | cout << " TEveTrack reference number " << fLabel << endl; | |
164 | cout << " ---------------------------------------------------------------------------------------------------------------------------------" << endl; | |
165 | cout << endl; | |
166 | cout << " Number of clusters " << fTrack->GetNClusters() << endl; | |
167 | } | |
168 | ||
169 | trackParamAtCluster = fTrack->GetTrackParamAtCluster(); | |
170 | nparam = trackParamAtCluster->GetEntries(); | |
171 | ||
172 | cout << endl; | |
173 | cout << " trackParamAtCluster entries " << nparam << "" << endl; | |
174 | cout << " ---------------------------------------------------------------------------------------------------------------------------------" << endl; | |
175 | cout << " Number InvBendMom BendSlope NonBendSlope BendCoord NonBendCoord Z Px Py Pz P" << endl; | |
176 | ||
177 | for (Int_t i = 0; i < nparam; i++) { | |
178 | ||
179 | mtp = (AliMUONTrackParam*)trackParamAtCluster->At(i); | |
180 | ||
181 | cout << | |
182 | setw(9)<< setprecision(3) << | |
183 | i << " " << | |
184 | ||
185 | setw(8) << setprecision(3) << | |
186 | mtp->GetInverseBendingMomentum() << " " << | |
187 | ||
188 | setw(8) << setprecision(3) << | |
189 | mtp->GetBendingSlope()*radDeg << " " << | |
190 | ||
191 | setw(8) << setprecision(3) << | |
192 | mtp->GetNonBendingSlope()*radDeg << " " << | |
193 | ||
194 | setw(8) << setprecision(4) << | |
195 | mtp->GetBendingCoor() << " " << | |
196 | ||
197 | setw(8) << setprecision(4) << | |
198 | mtp->GetNonBendingCoor() << " " << | |
199 | ||
200 | setw(10) << setprecision(6) << | |
201 | mtp->GetZ() << " " << | |
202 | ||
203 | setw(8) << setprecision(4) << | |
204 | mtp->Px() << " " << | |
205 | ||
206 | setw(8) << setprecision(4) << | |
207 | mtp->Py() << " " << | |
208 | ||
209 | setw(8) << setprecision(4) << | |
210 | mtp->Pz() << " " << | |
211 | ||
212 | setw(8) << setprecision(4) << | |
213 | mtp->P() << " " << | |
214 | ||
215 | endl; | |
216 | ||
217 | } | |
218 | ||
219 | cout << endl; | |
220 | cout << " TEveTrack parameters at vertex" << endl; | |
221 | cout << " --------------------------------------------------------------------------------------------------------------------" << endl; | |
222 | cout << " InvBendMom BendSlope NonBendSlope BendCoord NonBendCoord Z Px Py Pz P" << endl; | |
223 | ||
224 | mtp = (AliMUONTrackParam*)fTrack->GetTrackParamAtVertex(); | |
225 | ||
226 | bc = mtp->GetBendingCoor(); | |
227 | nbc = mtp->GetNonBendingCoor(); | |
228 | zc = mtp->GetZ(); | |
229 | if (bc < 0.001) bc = 0.0; | |
230 | if (nbc < 0.001) nbc = 0.0; | |
231 | if (zc < 0.001) zc = 0.0; | |
232 | ||
233 | cout << " " << | |
234 | setw(8) << setprecision(3) << | |
235 | mtp->GetInverseBendingMomentum() << " " << | |
236 | ||
237 | setw(8) << setprecision(3) << | |
238 | mtp->GetBendingSlope()*radDeg << " " << | |
239 | ||
240 | setw(8) << setprecision(3) << | |
241 | mtp->GetNonBendingSlope()*radDeg << " " << | |
242 | ||
243 | setw(8) << setprecision(4) << | |
244 | bc << " " << | |
245 | ||
246 | setw(8) << setprecision(4) << | |
247 | nbc << " " << | |
248 | ||
249 | setw(10) << setprecision(6) << | |
250 | zc << " " << | |
251 | ||
252 | setw(8) << setprecision(4) << | |
253 | mtp->Px() << " " << | |
254 | ||
255 | setw(8) << setprecision(4) << | |
256 | mtp->Py() << " " << | |
257 | ||
258 | setw(8) << setprecision(4) << | |
259 | mtp->Pz() << " " << | |
260 | ||
261 | setw(8) << setprecision(4) << | |
262 | mtp->P() << " " << | |
263 | ||
264 | endl; | |
265 | ||
266 | pt = TMath::Sqrt(mtp->Px()*mtp->Px()+mtp->Py()*mtp->Py()); | |
267 | ||
268 | cout << endl; | |
269 | cout << " Pt = " << | |
270 | setw(8) << setprecision(3) << | |
271 | pt << " GeV/c" << endl; | |
272 | ||
273 | } | |
274 | ||
275 | //______________________________________________________________________________ | |
276 | void AliEveMUONTrack::PrintMUONTriggerTrackInfo() | |
277 | { | |
278 | // | |
279 | // information about the trigger track | |
280 | // | |
281 | ||
282 | // Double_t radDeg = 180.0/TMath::Pi(); | |
283 | ||
284 | } | |
285 | ||
286 | //______________________________________________________________________________ | |
287 | void AliEveMUONTrack::PrintESDTrackInfo() | |
288 | { | |
289 | // | |
290 | // information about the reconstructed ESD track at vertex | |
291 | // | |
292 | ||
293 | Double_t radDeg = 180.0/TMath::Pi(); | |
294 | Float_t pt; | |
295 | ||
296 | AliMUONTrackParam *mtp = (AliMUONTrackParam*)fTrack->GetTrackParamAtVertex(); | |
297 | ||
298 | cout << endl; | |
299 | cout << " ESD muon track " << endl; | |
300 | cout << " -----------------------------------------------------------------------------------------------------------" << endl; | |
301 | cout << " InvBendMom BendSlope NonBendSlope BendCoord NonBendCoord Z Px Py Pz" << endl; | |
302 | ||
303 | cout << " " << | |
304 | ||
305 | setw(8) << setprecision(4) << | |
306 | mtp->GetInverseBendingMomentum() << " " << | |
307 | ||
308 | setw(8) << setprecision(3) << | |
309 | mtp->GetBendingSlope()*radDeg << " " << | |
310 | ||
311 | setw(8) << setprecision(3) << | |
312 | mtp->GetNonBendingSlope()*radDeg << " " << | |
313 | ||
314 | setw(8) << setprecision(4) << | |
315 | mtp->GetBendingCoor() << " " << | |
316 | ||
317 | setw(8) << setprecision(4) << | |
318 | mtp->GetNonBendingCoor() << " " << | |
319 | ||
320 | setw(10) << setprecision(6) << | |
321 | mtp->GetZ() << " " << | |
322 | ||
323 | setw(8) << setprecision(3) << | |
324 | mtp->Px() << " " << | |
325 | ||
326 | setw(8) << setprecision(3) << | |
327 | mtp->Py() << " " << | |
328 | ||
329 | setw(8) << setprecision(3) << | |
330 | mtp->Pz() << " " << | |
331 | ||
332 | endl; | |
333 | ||
334 | pt = TMath::Sqrt(mtp->Px()*mtp->Px()+mtp->Py()*mtp->Py()); | |
335 | ||
336 | cout << endl; | |
337 | cout << " Pt = " << | |
338 | setw(8) << setprecision(3) << | |
339 | pt << " GeV/c" << endl; | |
340 | ||
341 | cout << " P = " << | |
342 | setw(8) << setprecision(4) << | |
343 | mtp->P() << " GeV/c" << endl; | |
344 | ||
345 | AliESDEvent* esd = AliEveEventManager::AssertESD(); | |
346 | ||
347 | Double_t spdVertexX = 0; | |
348 | Double_t spdVertexY = 0; | |
349 | Double_t spdVertexZ = 0; | |
350 | Double_t esdVertexX = 0; | |
351 | Double_t esdVertexY = 0; | |
352 | Double_t esdVertexZ = 0; | |
353 | ||
354 | AliESDVertex* spdVertex = (AliESDVertex*) esd->GetVertex(); | |
355 | if (spdVertex->GetNContributors()) { | |
356 | spdVertexZ = spdVertex->GetZv(); | |
357 | spdVertexY = spdVertex->GetYv(); | |
358 | spdVertexX = spdVertex->GetXv(); | |
359 | } | |
360 | ||
361 | AliESDVertex* esdVertex = (AliESDVertex*) esd->GetPrimaryVertex(); | |
362 | if (esdVertex->GetNContributors()) { | |
363 | esdVertexZ = esdVertex->GetZv(); | |
364 | esdVertexY = esdVertex->GetYv(); | |
365 | esdVertexX = esdVertex->GetXv(); | |
366 | } | |
367 | ||
368 | Float_t t0v = esd->GetT0zVertex(); | |
369 | ||
370 | cout << endl; | |
371 | cout << endl; | |
372 | cout << "External vertex SPD: " << | |
373 | setw(3) << | |
374 | spdVertex->GetNContributors() << " " << | |
375 | setw(8) << setprecision(3) << | |
376 | spdVertexX << " " << | |
377 | spdVertexY << " " << | |
378 | spdVertexZ << " " << endl; | |
379 | cout << "External vertex ESD: " << | |
380 | setw(3) << | |
381 | esdVertex->GetNContributors() << " " << | |
382 | setw(8) << setprecision(3) << | |
383 | esdVertexX << " " << | |
384 | esdVertexY << " " << | |
385 | esdVertexZ << " " << endl; | |
386 | cout << "External vertex T0: " << | |
387 | setw(8) << setprecision(3) << | |
388 | t0v << " " << endl; | |
389 | ||
390 | } | |
391 | ||
392 | //______________________________________________________________________________ | |
393 | void AliEveMUONTrack::MUONTrackInfo() | |
394 | { | |
395 | // | |
396 | // MENU function | |
397 | // | |
398 | ||
399 | if (fIsMCTrack) { | |
400 | PrintMCTrackInfo(); | |
401 | } | |
402 | ||
403 | if (fIsMUONTrack || fIsRefTrack) { | |
404 | PrintMUONTrackInfo(); | |
405 | } | |
406 | ||
407 | if (fIsESDTrack) { | |
408 | PrintESDTrackInfo(); | |
409 | } | |
410 | ||
411 | if (fIsMUONTriggerTrack) { | |
412 | PrintMUONTriggerTrackInfo(); | |
413 | } | |
414 | ||
415 | cout << endl; | |
416 | cout << endl; | |
417 | cout << endl; | |
418 | cout << " (slopes [deg], coord [cm], p [GeV/c])" << endl; | |
419 | ||
420 | } | |
421 | ||
422 | //______________________________________________________________________________ | |
423 | void AliEveMUONTrack::MUONTriggerInfo() | |
424 | { | |
425 | // | |
426 | // MENU function | |
427 | // | |
428 | ||
429 | if (fIsMUONTrack) { | |
430 | TEveUtil::TEveUtil::LoadMacro("MUON_trigger_info.C"); | |
431 | gROOT->ProcessLine(Form("MUON_trigger_info(%d);", fLabel)); | |
432 | } | |
433 | if (fIsRefTrack) { | |
434 | cout << "This is a reference track!" << endl; | |
435 | } | |
436 | if (fIsMCTrack) { | |
437 | cout << "This is a Monte-Carlo track!" << endl; | |
438 | } | |
439 | if (fIsESDTrack) { | |
440 | ||
441 | AliESDEvent* esd = AliEveEventManager::AssertESD(); | |
442 | ULong64_t triggerMask = esd->GetTriggerMask(); | |
443 | ||
444 | cout << endl; | |
445 | cout << ">>>>>#########################################################################################################################" << endl; | |
446 | cout << endl; | |
447 | ||
448 | cout << " ESD track trigger info" << endl; | |
449 | cout << " -----------------------------------------------------" << endl; | |
450 | cout << endl; | |
451 | ||
452 | cout << " Match to trigger " << fTrack->GetMatchTrigger() << endl; | |
453 | cout << endl; | |
454 | cout << " ESD trigger mask = " << triggerMask << endl; | |
455 | ||
456 | cout << endl; | |
457 | cout << "#########################################################################################################################<<<<<" << endl; | |
458 | cout << endl; | |
459 | ||
460 | } | |
461 | ||
462 | } | |
463 | ||
464 | //______________________________________________________________________________ | |
465 | void AliEveMUONTrack::MakeMUONTrack(AliMUONTrack *mtrack) | |
466 | { | |
467 | // | |
468 | // builds the track with dipole field | |
469 | // | |
470 | ||
471 | if (!fIsESDTrack) { | |
472 | if (!fIsRefTrack) { | |
473 | fIsMUONTrack = kTRUE; | |
474 | fTrack = mtrack; | |
475 | } else { | |
476 | fTrack = new AliMUONTrack(*mtrack); | |
477 | } | |
478 | } | |
479 | ||
480 | Double_t xv, yv; | |
481 | Float_t ax, bx, ay, by; | |
482 | Float_t xr[28], yr[28], zr[28]; | |
483 | Float_t xrc[28], yrc[28], zrc[28]; | |
484 | Int_t chr[28]; | |
485 | char form[1000]; | |
486 | ||
487 | TMatrixD smatrix(2,2); | |
488 | TMatrixD sums(2,1); | |
489 | TMatrixD res(2,1); | |
490 | ||
491 | // middle z between the two detector planes of the trigger chambers | |
492 | Float_t zg[4] = { -1603.5, -1620.5, -1703.5, -1720.5 }; | |
493 | ||
494 | Float_t pt = 0.0; | |
495 | Float_t pv[3] = {0., 0., 0.}; | |
496 | ||
497 | if (fIsMUONTrack) { | |
498 | if (mtrack->GetMatchTrigger()) { | |
499 | sprintf(form,"MUONTrack %2d (MT)", fLabel); | |
500 | } else { | |
501 | sprintf(form,"MUONTrack %2d ", fLabel); | |
502 | } | |
503 | SetName(form); | |
504 | SetLineStyle(1); | |
505 | } | |
506 | ||
507 | AliMUONTrackParam *trackParam = 0x0; | |
508 | if (fIsMUONTrack || fIsESDTrack) { | |
509 | trackParam = mtrack->GetTrackParamAtVertex(); | |
510 | SetPoint(fCount,trackParam->GetNonBendingCoor(),trackParam->GetBendingCoor(),trackParam->GetZ()); | |
511 | fCount++; | |
512 | } | |
513 | ||
514 | for (Int_t i = 0; i < 28; i++) { | |
515 | xr[i]=yr[i]=zr[i]=0.0; | |
516 | chr[i]=-1; | |
517 | } | |
518 | ||
519 | Int_t nTrackHits = mtrack->GetNClusters(); | |
520 | TClonesArray* trackParamAtCluster = mtrack->GetTrackParamAtCluster(); | |
521 | for (Int_t iHit = 0; iHit < nTrackHits; iHit++){ | |
522 | ||
523 | trackParam = (AliMUONTrackParam*) trackParamAtCluster->At(iHit); | |
524 | ||
525 | if (iHit == 0) { | |
526 | if (IsMUONTrack() || IsESDTrack()) { | |
527 | pt = TMath::Sqrt(trackParam->Px()*trackParam->Px()+trackParam->Py()*trackParam->Py()); | |
528 | SetLineColor(ColorIndex(pt)); | |
529 | } | |
530 | pv[0] = trackParam->Px(); | |
531 | pv[1] = trackParam->Py(); | |
532 | pv[2] = trackParam->Pz(); | |
533 | fP.Set(pv); | |
534 | } | |
535 | ||
536 | xr[iHit] = trackParam->GetNonBendingCoor(); | |
537 | yr[iHit] = trackParam->GetBendingCoor(); | |
538 | zr[iHit] = trackParam->GetZ(); | |
539 | chr[iHit] = trackParam->GetClusterPtr()->GetChamberId(); | |
540 | ||
541 | } | |
542 | ||
543 | SetPoint(fCount,xr[0],yr[0],zr[0]); | |
544 | fCount++; | |
545 | for (Int_t iHit = 1; iHit < nTrackHits; iHit++) { | |
546 | if (chr[iHit] > 3 && chr[iHit-1] < 6) Propagate(xr,yr,zr,iHit-1,iHit); | |
547 | SetPoint(fCount,xr[iHit],yr[iHit],zr[iHit]); | |
548 | fCount++; | |
549 | } | |
550 | ||
551 | if (!fIsMUONTrack && !fIsESDTrack) return; | |
552 | ||
553 | Int_t nrc = 0; | |
554 | if (mtrack->GetMatchTrigger()) { | |
555 | ||
556 | for (Int_t i = 0; i < nTrackHits; i++) { | |
557 | if (TMath::Abs(zr[i]) > 1000.0) { | |
558 | //printf("TEveHit %d x %f y %f z %f \n",iHit,xr[i],yr[i],zr[i]); | |
559 | xrc[nrc] = xr[i]; | |
560 | yrc[nrc] = yr[i]; | |
561 | zrc[nrc] = zr[i]; | |
562 | nrc++; | |
563 | } | |
564 | } | |
565 | ||
566 | if (nrc < 2) return; | |
567 | ||
568 | // fit x-z | |
569 | smatrix.Zero(); | |
570 | sums.Zero(); | |
571 | for (Int_t i = 0; i < nrc; i++) { | |
572 | xv = (Double_t)zrc[i]; | |
573 | yv = (Double_t)xrc[i]; | |
574 | //printf("x-z: xv %f yv %f \n",xv,yv); | |
575 | smatrix(0,0) += 1.0; | |
576 | smatrix(1,1) += xv*xv; | |
577 | smatrix(0,1) += xv; | |
578 | smatrix(1,0) += xv; | |
579 | sums(0,0) += yv; | |
580 | sums(1,0) += xv*yv; | |
581 | } | |
582 | res = smatrix.Invert() * sums; | |
583 | ax = res(0,0); | |
584 | bx = res(1,0); | |
585 | ||
586 | // fit y-z | |
587 | smatrix.Zero(); | |
588 | sums.Zero(); | |
589 | for (Int_t i = 0; i < nrc; i++) { | |
590 | xv = (Double_t)zrc[i]; | |
591 | yv = (Double_t)yrc[i]; | |
592 | //printf("y-z: xv %f yv %f \n",xv,yv); | |
593 | smatrix(0,0) += 1.0; | |
594 | smatrix(1,1) += xv*xv; | |
595 | smatrix(0,1) += xv; | |
596 | smatrix(1,0) += xv; | |
597 | sums(0,0) += yv; | |
598 | sums(1,0) += xv*yv; | |
599 | } | |
600 | res = smatrix.Invert() * sums; | |
601 | ay = res(0,0); | |
602 | by = res(1,0); | |
603 | ||
604 | Float_t xtc, ytc, ztc; | |
605 | for (Int_t ii = 0; ii < 4; ii++) { | |
606 | ||
607 | ztc = zg[ii]; | |
608 | ytc = ay+by*zg[ii]; | |
609 | xtc = ax+bx*zg[ii]; | |
610 | ||
611 | //printf("tc: x %f y %f z %f \n",xtc,ytc,ztc); | |
612 | ||
613 | SetPoint(fCount,xtc,ytc,ztc); | |
614 | fCount++; | |
615 | ||
616 | } | |
617 | ||
618 | } // end match trigger | |
619 | ||
620 | } | |
621 | ||
622 | //______________________________________________________________________________ | |
623 | void AliEveMUONTrack::MakeMUONTriggerTrack(AliMUONTriggerTrack *mtrack) | |
624 | { | |
625 | // | |
626 | // builds the trigger track from one point and direction | |
627 | // | |
628 | ||
629 | Float_t x1 = mtrack->GetX11(); | |
630 | Float_t y1 = mtrack->GetY11(); | |
631 | Float_t thex = mtrack->GetThetax(); | |
632 | Float_t they = mtrack->GetThetay(); | |
633 | ||
634 | Float_t z11 = -1600.0; | |
635 | Float_t z22 = -1724.0; | |
636 | Float_t dz = z22-z11; | |
637 | ||
638 | Float_t x2 = x1 + dz*TMath::Tan(thex); | |
639 | Float_t y2 = y1 + dz*TMath::Tan(they); | |
640 | ||
641 | SetPoint(fCount,x1,y1,z11); fCount++; | |
642 | SetPoint(fCount,x2,y2,z22); fCount++; | |
643 | ||
644 | char form[1000]; | |
645 | ||
646 | sprintf(form,"MUONTriggerTrack %2d",mtrack->GetLoTrgNum()); | |
647 | SetName(form); | |
648 | SetLineStyle(1); | |
649 | ||
650 | } | |
651 | ||
652 | //______________________________________________________________________________ | |
653 | void AliEveMUONTrack::MakeESDTrack(AliESDMuonTrack *mtrack) | |
654 | { | |
655 | // | |
656 | // builds the track with dipole field starting from the TParticle | |
657 | // | |
658 | ||
659 | fIsESDTrack = kTRUE; | |
660 | ||
661 | char form[1000]; | |
662 | if (mtrack->GetMatchTrigger()) { | |
663 | sprintf(form,"ESDTrack %2d (MT)", fLabel); | |
664 | } else { | |
665 | sprintf(form,"ESDTrack %2d ", fLabel); | |
666 | } | |
667 | SetName(form); | |
668 | SetLineStyle(3); | |
669 | SetLineColor(0); | |
670 | ||
671 | fTrack = new AliMUONTrack(); | |
672 | ||
673 | // create a simple track from the ESD track | |
674 | AliMUONESDInterface::ESDToMUON(*mtrack,*fTrack); | |
675 | ||
676 | // reset track parameters at vertex to the ones at DCA | |
677 | AliMUONTrackParam paramAtDCA; | |
678 | AliMUONESDInterface::GetParamAtDCA(*mtrack, paramAtDCA); | |
679 | fTrack->SetTrackParamAtVertex(¶mAtDCA); | |
680 | ||
681 | MakeMUONTrack(fTrack); | |
682 | ||
683 | } | |
684 | ||
685 | //______________________________________________________________________________ | |
686 | void AliEveMUONTrack::MakeMCTrack(TParticle *part) | |
687 | { | |
688 | // | |
689 | // builds the track with dipole field starting from the TParticle | |
690 | // | |
691 | ||
692 | fIsMCTrack = kTRUE; | |
693 | ||
694 | fPart = new TParticle(*part); | |
695 | ||
696 | char form[1000]; | |
697 | sprintf(form,"TEveMCTrack %2d ", fLabel); | |
698 | SetName(form); | |
699 | SetLineStyle(2); | |
700 | SetLineColor(8); | |
701 | ||
702 | Double_t vect[7], vout[7]; | |
703 | Double_t step = 1.0; | |
704 | ||
705 | Float_t pv[3]; | |
706 | pv[0] = fPart->Px(); | |
707 | pv[1] = fPart->Py(); | |
708 | pv[2] = fPart->Pz(); | |
709 | fP.Set(pv); | |
710 | ||
711 | vect[0] = fPart->Vx(); | |
712 | vect[1] = fPart->Vy(); | |
713 | vect[2] = fPart->Vz(); | |
714 | vect[3] = fPart->Px()/fPart->P(); | |
715 | vect[4] = fPart->Py()/fPart->P(); | |
716 | vect[5] = fPart->Pz()/fPart->P(); | |
717 | vect[6] = fPart->P(); | |
718 | ||
719 | TParticlePDG *ppdg = fPart->GetPDG(1); | |
720 | Int_t charge = (Int_t)(ppdg->Charge()/3.0); | |
721 | ||
722 | Double_t zMax = -1750.0; | |
723 | Double_t rMax = 350.0; | |
724 | Double_t r = 0.0; | |
725 | ||
726 | Int_t nSteps = 0; | |
727 | while ((vect[2] > zMax) && (nSteps < 10000) && (r < rMax)) { | |
728 | nSteps++; | |
729 | OneStepRungekutta(charge, step, vect, vout); | |
730 | SetPoint(fCount,vout[0],vout[1],vout[2]); | |
731 | fCount++; | |
732 | for (Int_t i = 0; i < 7; i++) { | |
733 | vect[i] = vout[i]; | |
734 | } | |
735 | r = TMath::Sqrt(vect[0]*vect[0]+vect[1]*vect[1]); | |
736 | } | |
737 | ||
738 | } | |
739 | ||
740 | //______________________________________________________________________________ | |
741 | void AliEveMUONTrack::MakeRefTrack(AliMUONTrack *mtrack) | |
742 | { | |
743 | // | |
744 | // builds the track with dipole field starting from the TParticle | |
745 | // | |
746 | ||
747 | fIsRefTrack = kTRUE; | |
748 | ||
749 | char form[1000]; | |
750 | sprintf(form,"RefTrack %2d ", fLabel); | |
751 | SetName(form); | |
752 | SetLineStyle(2); | |
753 | SetLineColor(0); | |
754 | ||
755 | MakeMUONTrack(mtrack); | |
756 | ||
757 | } | |
758 | ||
759 | //______________________________________________________________________________ | |
760 | void AliEveMUONTrack::Propagate(Float_t *xr, Float_t *yr, Float_t *zr, Int_t i1, Int_t i2) | |
761 | { | |
762 | // | |
763 | // propagate in magnetic field between hits of indices i1 and i2 | |
764 | // | |
765 | ||
766 | Double_t vect[7], vout[7]; | |
767 | Double_t step = 1.0; | |
768 | Double_t zMax = 0.0; | |
769 | Int_t charge = 0; | |
770 | AliMUONTrackParam *trackParam = 0; | |
771 | TClonesArray *trackParamAtCluster = 0; | |
772 | ||
773 | if (i2 == 9999) { | |
774 | zMax = zr[i1]+1.5*step; | |
775 | } else { | |
776 | zMax = zr[i2]+1.5*step; | |
777 | } | |
778 | ||
779 | trackParamAtCluster = fTrack->GetTrackParamAtCluster(); | |
780 | ||
781 | if (IsMUONTrack() || IsESDTrack()) { | |
782 | trackParam = (AliMUONTrackParam*)trackParamAtCluster->At(i1); | |
783 | charge = (Int_t)trackParam->GetCharge(); | |
784 | } | |
785 | if (IsRefTrack()) { | |
786 | trackParam = fTrack->GetTrackParamAtVertex(); | |
787 | charge = (Int_t)trackParam->GetCharge(); | |
788 | trackParam = (AliMUONTrackParam*)trackParamAtCluster->At(i1); | |
789 | } | |
790 | ||
791 | vect[0] = xr[i1]; | |
792 | vect[1] = yr[i1]; | |
793 | vect[2] = zr[i1]; | |
794 | vect[6] = trackParam->P(); | |
795 | vect[3] = trackParam->Px()/vect[6]; | |
796 | vect[4] = trackParam->Py()/vect[6]; | |
797 | vect[5] = trackParam->Pz()/vect[6]; | |
798 | ||
799 | Int_t nSteps = 0; | |
800 | while ((vect[2] > zMax) && (nSteps < 10000)) { | |
801 | nSteps++; | |
802 | OneStepRungekutta(charge, step, vect, vout); | |
803 | SetPoint(fCount,vout[0],vout[1],vout[2]); | |
804 | fCount++; | |
805 | for (Int_t i = 0; i < 7; i++) { | |
806 | vect[i] = vout[i]; | |
807 | } | |
808 | } | |
809 | ||
810 | } | |
811 | ||
812 | //______________________________________________________________________________ | |
813 | void AliEveMUONTrack::OneStepRungekutta(Double_t charge, Double_t step, | |
814 | Double_t* vect, Double_t* vout) | |
815 | { | |
816 | /// ****************************************************************** | |
817 | /// * * | |
818 | /// * Runge-Kutta method for tracking a particle through a magnetic * | |
819 | /// * field. Uses Nystroem algorithm (See Handbook Nat. Bur. of * | |
820 | /// * Standards, procedure 25.5.20) * | |
821 | /// * * | |
822 | /// * Input parameters * | |
823 | /// * CHARGE Particle charge * | |
824 | /// * STEP Step size * | |
825 | /// * VECT Initial co-ords,direction cosines,momentum * | |
826 | /// * Output parameters * | |
827 | /// * VOUT Output co-ords,direction cosines,momentum * | |
828 | /// * User routine called * | |
829 | /// * CALL GUFLD(X,F) * | |
830 | /// * * | |
831 | /// * ==>Called by : <USER>, GUSWIM * | |
832 | /// * Authors R.Brun, M.Hansroul ********* * | |
833 | /// * V.Perevoztchikov (CUT STEP implementation) * | |
834 | /// * * | |
835 | /// * * | |
836 | /// ****************************************************************** | |
837 | ||
838 | Double_t h2, h4, f[4]; | |
839 | Double_t xyzt[3], a, b, c, ph,ph2; | |
840 | Double_t secxs[4],secys[4],seczs[4],hxp[3]; | |
841 | Double_t g1, g2, g3, g4, g5, g6, ang2, dxt, dyt, dzt; | |
842 | Double_t est, at, bt, ct, cba; | |
843 | Double_t f1, f2, f3, f4, rho, tet, hnorm, hp, rho1, sint, cost; | |
844 | ||
845 | Double_t x; | |
846 | Double_t y; | |
847 | Double_t z; | |
848 | ||
849 | Double_t xt; | |
850 | Double_t yt; | |
851 | Double_t zt; | |
852 | ||
853 | Double_t maxit = 1992; | |
854 | Double_t maxcut = 11; | |
855 | ||
856 | const Double_t kdlt = 1e-4; | |
857 | const Double_t kdlt32 = kdlt/32.; | |
858 | const Double_t kthird = 1./3.; | |
859 | const Double_t khalf = 0.5; | |
860 | const Double_t kec = 2.9979251e-4; | |
861 | ||
862 | const Double_t kpisqua = 9.86960440109; | |
863 | const Int_t kix = 0; | |
864 | const Int_t kiy = 1; | |
865 | const Int_t kiz = 2; | |
866 | const Int_t kipx = 3; | |
867 | const Int_t kipy = 4; | |
868 | const Int_t kipz = 5; | |
869 | ||
870 | // *. | |
871 | // *. ------------------------------------------------------------------ | |
872 | // *. | |
873 | // * this constant is for units cm,gev/c and kgauss | |
874 | // * | |
875 | Int_t iter = 0; | |
876 | Int_t ncut = 0; | |
877 | for(Int_t j = 0; j < 7; j++) | |
878 | vout[j] = vect[j]; | |
879 | ||
880 | Double_t pinv = kec * charge / vect[6]; | |
881 | Double_t tl = 0.; | |
882 | Double_t h = step; | |
883 | Double_t rest; | |
884 | ||
885 | ||
886 | do { | |
887 | rest = step - tl; | |
888 | if (TMath::Abs(h) > TMath::Abs(rest)) h = rest; | |
889 | //cmodif: call gufld(vout,f) changed into: | |
890 | TGeoGlobalMagField::Instance()->Field(vout,f); | |
891 | ||
892 | // * | |
893 | // * start of integration | |
894 | // * | |
895 | x = vout[0]; | |
896 | y = vout[1]; | |
897 | z = vout[2]; | |
898 | a = vout[3]; | |
899 | b = vout[4]; | |
900 | c = vout[5]; | |
901 | ||
902 | h2 = khalf * h; | |
903 | h4 = khalf * h2; | |
904 | ph = pinv * h; | |
905 | ph2 = khalf * ph; | |
906 | secxs[0] = (b * f[2] - c * f[1]) * ph2; | |
907 | secys[0] = (c * f[0] - a * f[2]) * ph2; | |
908 | seczs[0] = (a * f[1] - b * f[0]) * ph2; | |
909 | ang2 = (secxs[0]*secxs[0] + secys[0]*secys[0] + seczs[0]*seczs[0]); | |
910 | if (ang2 > kpisqua) break; | |
911 | ||
912 | dxt = h2 * a + h4 * secxs[0]; | |
913 | dyt = h2 * b + h4 * secys[0]; | |
914 | dzt = h2 * c + h4 * seczs[0]; | |
915 | xt = x + dxt; | |
916 | yt = y + dyt; | |
917 | zt = z + dzt; | |
918 | // * | |
919 | // * second intermediate point | |
920 | // * | |
921 | ||
922 | est = TMath::Abs(dxt) + TMath::Abs(dyt) + TMath::Abs(dzt); | |
923 | if (est > h) { | |
924 | if (ncut++ > maxcut) break; | |
925 | h *= khalf; | |
926 | continue; | |
927 | } | |
928 | ||
929 | xyzt[0] = xt; | |
930 | xyzt[1] = yt; | |
931 | xyzt[2] = zt; | |
932 | ||
933 | //cmodif: call gufld(xyzt,f) changed into: | |
934 | TGeoGlobalMagField::Instance()->Field(xyzt,f); | |
935 | ||
936 | at = a + secxs[0]; | |
937 | bt = b + secys[0]; | |
938 | ct = c + seczs[0]; | |
939 | ||
940 | secxs[1] = (bt * f[2] - ct * f[1]) * ph2; | |
941 | secys[1] = (ct * f[0] - at * f[2]) * ph2; | |
942 | seczs[1] = (at * f[1] - bt * f[0]) * ph2; | |
943 | at = a + secxs[1]; | |
944 | bt = b + secys[1]; | |
945 | ct = c + seczs[1]; | |
946 | secxs[2] = (bt * f[2] - ct * f[1]) * ph2; | |
947 | secys[2] = (ct * f[0] - at * f[2]) * ph2; | |
948 | seczs[2] = (at * f[1] - bt * f[0]) * ph2; | |
949 | dxt = h * (a + secxs[2]); | |
950 | dyt = h * (b + secys[2]); | |
951 | dzt = h * (c + seczs[2]); | |
952 | xt = x + dxt; | |
953 | yt = y + dyt; | |
954 | zt = z + dzt; | |
955 | at = a + 2.*secxs[2]; | |
956 | bt = b + 2.*secys[2]; | |
957 | ct = c + 2.*seczs[2]; | |
958 | ||
959 | est = TMath::Abs(dxt)+TMath::Abs(dyt)+TMath::Abs(dzt); | |
960 | if (est > 2.*TMath::Abs(h)) { | |
961 | if (ncut++ > maxcut) break; | |
962 | h *= khalf; | |
963 | continue; | |
964 | } | |
965 | ||
966 | xyzt[0] = xt; | |
967 | xyzt[1] = yt; | |
968 | xyzt[2] = zt; | |
969 | ||
970 | //cmodif: call gufld(xyzt,f) changed into: | |
971 | TGeoGlobalMagField::Instance()->Field(xyzt,f); | |
972 | ||
973 | z = z + (c + (seczs[0] + seczs[1] + seczs[2]) * kthird) * h; | |
974 | y = y + (b + (secys[0] + secys[1] + secys[2]) * kthird) * h; | |
975 | x = x + (a + (secxs[0] + secxs[1] + secxs[2]) * kthird) * h; | |
976 | ||
977 | secxs[3] = (bt*f[2] - ct*f[1])* ph2; | |
978 | secys[3] = (ct*f[0] - at*f[2])* ph2; | |
979 | seczs[3] = (at*f[1] - bt*f[0])* ph2; | |
980 | a = a+(secxs[0]+secxs[3]+2. * (secxs[1]+secxs[2])) * kthird; | |
981 | b = b+(secys[0]+secys[3]+2. * (secys[1]+secys[2])) * kthird; | |
982 | c = c+(seczs[0]+seczs[3]+2. * (seczs[1]+seczs[2])) * kthird; | |
983 | ||
984 | est = TMath::Abs(secxs[0]+secxs[3] - (secxs[1]+secxs[2])) | |
985 | + TMath::Abs(secys[0]+secys[3] - (secys[1]+secys[2])) | |
986 | + TMath::Abs(seczs[0]+seczs[3] - (seczs[1]+seczs[2])); | |
987 | ||
988 | if (est > kdlt && TMath::Abs(h) > 1.e-4) { | |
989 | if (ncut++ > maxcut) break; | |
990 | h *= khalf; | |
991 | continue; | |
992 | } | |
993 | ||
994 | ncut = 0; | |
995 | // * if too many iterations, go to helix | |
996 | if (iter++ > maxit) break; | |
997 | ||
998 | tl += h; | |
999 | if (est < kdlt32) | |
1000 | h *= 2.; | |
1001 | cba = 1./ TMath::Sqrt(a*a + b*b + c*c); | |
1002 | vout[0] = x; | |
1003 | vout[1] = y; | |
1004 | vout[2] = z; | |
1005 | vout[3] = cba*a; | |
1006 | vout[4] = cba*b; | |
1007 | vout[5] = cba*c; | |
1008 | rest = step - tl; | |
1009 | if (step < 0.) rest = -rest; | |
1010 | if (rest < 1.e-5*TMath::Abs(step)) return; | |
1011 | ||
1012 | } while(1); | |
1013 | ||
1014 | // angle too big, use helix | |
1015 | ||
1016 | f1 = f[0]; | |
1017 | f2 = f[1]; | |
1018 | f3 = f[2]; | |
1019 | f4 = TMath::Sqrt(f1*f1+f2*f2+f3*f3); | |
1020 | rho = -f4*pinv; | |
1021 | tet = rho * step; | |
1022 | ||
1023 | hnorm = 1./f4; | |
1024 | f1 = f1*hnorm; | |
1025 | f2 = f2*hnorm; | |
1026 | f3 = f3*hnorm; | |
1027 | ||
1028 | hxp[0] = f2*vect[kipz] - f3*vect[kipy]; | |
1029 | hxp[1] = f3*vect[kipx] - f1*vect[kipz]; | |
1030 | hxp[2] = f1*vect[kipy] - f2*vect[kipx]; | |
1031 | ||
1032 | hp = f1*vect[kipx] + f2*vect[kipy] + f3*vect[kipz]; | |
1033 | ||
1034 | rho1 = 1./rho; | |
1035 | sint = TMath::Sin(tet); | |
1036 | cost = 2.*TMath::Sin(khalf*tet)*TMath::Sin(khalf*tet); | |
1037 | ||
1038 | g1 = sint*rho1; | |
1039 | g2 = cost*rho1; | |
1040 | g3 = (tet-sint) * hp*rho1; | |
1041 | g4 = -cost; | |
1042 | g5 = sint; | |
1043 | g6 = cost * hp; | |
1044 | ||
1045 | vout[kix] = vect[kix] + g1*vect[kipx] + g2*hxp[0] + g3*f1; | |
1046 | vout[kiy] = vect[kiy] + g1*vect[kipy] + g2*hxp[1] + g3*f2; | |
1047 | vout[kiz] = vect[kiz] + g1*vect[kipz] + g2*hxp[2] + g3*f3; | |
1048 | ||
1049 | vout[kipx] = vect[kipx] + g4*vect[kipx] + g5*hxp[0] + g6*f1; | |
1050 | vout[kipy] = vect[kipy] + g4*vect[kipy] + g5*hxp[1] + g6*f2; | |
1051 | vout[kipz] = vect[kipz] + g4*vect[kipz] + g5*hxp[2] + g6*f3; | |
1052 | ||
1053 | return; | |
1054 | } | |
1055 | ||
1056 | //______________________________________________________________________________ | |
1057 | Int_t AliEveMUONTrack::ColorIndex(Float_t val) | |
1058 | { | |
1059 | // | |
1060 | // returns color index in the palette for a give value | |
1061 | // | |
1062 | ||
1063 | Float_t threshold = 0.0; | |
1064 | Float_t maxVal = 10.0; | |
1065 | ||
1066 | Float_t div = TMath::Max(1, (Int_t)(maxVal - threshold)); | |
1067 | Int_t nCol = gStyle->GetNumberOfColors(); | |
1068 | Int_t cBin = (Int_t) TMath::Nint(nCol*(val - threshold)/div); | |
1069 | ||
1070 | return gStyle->GetColorPalette(TMath::Min(nCol - 1, cBin)); | |
1071 | ||
1072 | } | |
1073 | ||
1074 | //============================================================================== | |
1075 | // Temporary AliEveMUONTrackList | |
1076 | //============================================================================== | |
1077 | ||
1078 | //______________________________________________________________________________ | |
1079 | void AliEveMUONTrackList::HackMomentumLimits(Bool_t recurse) | |
1080 | { | |
1081 | // Find momentum limits from included tracks. | |
1082 | ||
1083 | fLimPt = fLimP = 0; | |
1084 | ||
1085 | for (List_i i=BeginChildren(); i!=EndChildren(); ++i) | |
1086 | { | |
1087 | TEveTrack* track = dynamic_cast<TEveTrack*>(*i); | |
1088 | if (track) | |
1089 | { | |
1090 | fLimPt = TMath::Max(fLimPt, track->GetMomentum().Perp()); | |
1091 | fLimP = TMath::Max(fLimP, track->GetMomentum().Mag()); | |
1092 | } | |
1093 | if (recurse) | |
1094 | FindMomentumLimits(*i, recurse); | |
1095 | } | |
1096 | ||
1097 | fLimPt = RoundMomentumLimit(fLimPt); | |
1098 | fLimP = RoundMomentumLimit(fLimP); | |
1099 | if (fMaxPt == 0) fMaxPt = fLimPt; | |
1100 | if (fMaxP == 0) fMaxP = fLimP; | |
1101 | } |