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