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