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