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