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