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