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 <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   Double_t xv, yv;
 485   Float_t ax, bx, ay, by;
 486   Float_t xr[28], yr[28], zr[28];
 487   Float_t xrc[28], yrc[28], zrc[28];
 488   Int_t chr[28];
 489   char form[1000];
 490
 491   TMatrixD smatrix(2,2);
 492   TMatrixD sums(2,1);
 493   TMatrixD res(2,1);
 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., 0.};
 500
 501   if (fIsMUONTrack) {
 502     if (mtrack->GetMatchTrigger()) {
 503       snprintf(form,1000,"MUONTrack %2d (MT)", fLabel);
 504     } else {
 505       snprintf(form,1000,"MUONTrack %2d     ", fLabel);
 506     }
 507     SetName(form);
 508     SetLineStyle(1);
 509   }
 510
 511   AliMUONTrackParam *trackParam = 0x0;
 512   if (fIsMUONTrack || fIsESDTrack) {
 513     trackParam = mtrack->GetTrackParamAtVertex();
 514     SetPoint(fCount,trackParam->GetNonBendingCoor(),trackParam->GetBendingCoor(),trackParam->GetZ());
 515     fCount++;
 516   }
 517
 518   for (Int_t i = 0; i < 28; i++) {
 519     xr[i]=yr[i]=zr[i]=0.0;
 520     chr[i]=-1;
 521   }
 522
 523   Int_t nTrackHits = mtrack->GetNClusters();
 524   TObjArray* trackParamAtCluster = mtrack->GetTrackParamAtCluster();
 525   for (Int_t iHit = 0; iHit < nTrackHits; iHit++){
 526
 527     trackParam = (AliMUONTrackParam*) trackParamAtCluster->At(iHit);
 528
 529     if (iHit == 0) {
 530       if (IsMUONTrack() || IsESDTrack()) {
 531         pt = TMath::Sqrt(trackParam->Px()*trackParam->Px()+trackParam->Py()*trackParam->Py());
 532         SetLineColor(ColorIndex(pt));
 533       }
 534       pv[0] = trackParam->Px();
 535       pv[1] = trackParam->Py();
 536       pv[2] = trackParam->Pz();
 537       fP.Set(pv);
 538     }
 539
 540     xr[iHit] = trackParam->GetNonBendingCoor();
 541     yr[iHit] = trackParam->GetBendingCoor();
 542     zr[iHit] = trackParam->GetZ();
 543     chr[iHit] = trackParam->GetClusterPtr()->GetChamberId();
 544
 545   }
 546
 547   SetPoint(fCount,xr[0],yr[0],zr[0]);
 548   fCount++;
 549   for (Int_t iHit = 1; iHit < nTrackHits; iHit++) {
 550     if (chr[iHit] > 3 && chr[iHit-1] < 6) Propagate(xr,yr,zr,iHit-1,iHit);
 551     SetPoint(fCount,xr[iHit],yr[iHit],zr[iHit]);
 552     fCount++;
 553   }
 554
 555   if (!fIsMUONTrack && !fIsESDTrack) return;
 556
 557   Int_t nrc = 0;
 558   if (mtrack->GetMatchTrigger()) {
 559
 560     for (Int_t i = 0; i < nTrackHits; i++) {
 561       if (TMath::Abs(zr[i]) > 1000.0) {
 562         //printf("TEveHit %d x %f y %f z %f \n",iHit,xr[i],yr[i],zr[i]);
 563         xrc[nrc] = xr[i];
 564         yrc[nrc] = yr[i];
 565         zrc[nrc] = zr[i];
 566         nrc++;
 567       }
 568     }
 569
 570     if (nrc < 2) return;
 571
 572     // fit x-z
 573     smatrix.Zero();
 574     sums.Zero();
 575     for (Int_t i = 0; i < nrc; i++) {
 576       xv = (Double_t)zrc[i];
 577       yv = (Double_t)xrc[i];
 578       //printf("x-z: xv %f yv %f \n",xv,yv);
 579       smatrix(0,0) += 1.0;
 580       smatrix(1,1) += xv*xv;
 581       smatrix(0,1) += xv;
 582       smatrix(1,0) += xv;
 583       sums(0,0)    += yv;
 584       sums(1,0)    += xv*yv;
 585     }
 586     res = smatrix.Invert() * sums;
 587     ax = res(0,0);
 588     bx = res(1,0);
 589
 590     // fit y-z
 591     smatrix.Zero();
 592     sums.Zero();
 593     for (Int_t i = 0; i < nrc; i++) {
 594       xv = (Double_t)zrc[i];
 595       yv = (Double_t)yrc[i];
 596       //printf("y-z: xv %f yv %f \n",xv,yv);
 597       smatrix(0,0) += 1.0;
 598       smatrix(1,1) += xv*xv;
 599       smatrix(0,1) += xv;
 600       smatrix(1,0) += xv;
 601       sums(0,0)    += yv;
 602       sums(1,0)    += xv*yv;
 603     }
 604     res = smatrix.Invert() * sums;
 605     ay = res(0,0);
 606     by = res(1,0);
 607
 608     Float_t xtc, ytc, ztc;
 609     for (Int_t ii = 0; ii < 4; ii++) {
 610
 611       ztc = zg[ii];
 612       ytc = ay+by*zg[ii];
 613       xtc = ax+bx*zg[ii];
 614
 615       //printf("tc: x %f y %f z %f \n",xtc,ytc,ztc);
 616
 617       SetPoint(fCount,xtc,ytc,ztc);
 618       fCount++;
 619
 620     }
 621
 622   }  // end match trigger
 623
 624 }
 625
 626 //______________________________________________________________________________
 627 void AliEveMUONTrack::MakeMUONTriggerTrack(AliMUONTriggerTrack *mtrack)
 628 {
 629   //
 630   // builds the trigger track from one point and direction
 631   //
 632
 633   Float_t x1   = mtrack->GetX11();
 634   Float_t y1   = mtrack->GetY11();
 635   Float_t thex = mtrack->GetThetax();
 636   Float_t they = mtrack->GetThetay();
 637
 638   Float_t z11 = -1600.0;
 639   Float_t z22 = -1724.0;
 640   Float_t dz  = z22-z11;
 641
 642   Float_t x2 = x1 + dz*TMath::Tan(thex);
 643   Float_t y2 = y1 + dz*TMath::Tan(they);
 644
 645   SetPoint(fCount,x1,y1,z11); fCount++;
 646   SetPoint(fCount,x2,y2,z22); fCount++;
 647
 648   char form[1000];
 649
 650   snprintf(form,1000,"MUONTriggerTrack %2d",mtrack->GetLoTrgNum());
 651   SetName(form);
 652   SetLineStyle(1);
 653
 654 }
 655
 656 //______________________________________________________________________________
 657 void AliEveMUONTrack::MakeESDTrack(AliESDMuonTrack *mtrack)
 658 {
 659   //
 660   // builds the track with dipole field starting from the TParticle
 661   //
 662
 663   fIsESDTrack = kTRUE;
 664
 665   char form[1000];
 666   if (mtrack->GetMatchTrigger()) {
 667     snprintf(form,1000,"ESDTrack %2d (MT)", fLabel);
 668   } else {
 669     snprintf(form,1000,"ESDTrack %2d     ", fLabel);
 670   }
 671   SetName(form);
 672   SetLineStyle(3);
 673   SetLineColor(0);
 674
 675   fTrack = new AliMUONTrack();
 676
 677   // create a simple track from the ESD track
 678   AliMUONESDInterface::ESDToMUON(*mtrack,*fTrack);
 679
 680   // reset track parameters at vertex to the ones at DCA
 681   AliMUONTrackParam paramAtDCA;
 682   AliMUONESDInterface::GetParamAtDCA(*mtrack, paramAtDCA);
 683   fTrack->SetTrackParamAtVertex(&paramAtDCA);
 684
 685   MakeMUONTrack(fTrack);
 686
 687 }
 688
 689 //______________________________________________________________________________
 690 void AliEveMUONTrack::MakeMCTrack(TParticle *part)
 691 {
 692   //
 693   // builds the track with dipole field starting from the TParticle
 694   //
 695
 696   fIsMCTrack = kTRUE;
 697
 698   fPart     = new TParticle(*part);
 699
 700   char form[1000];
 701   snprintf(form,1000,"TEveMCTrack %2d ", fLabel);
 702   SetName(form);
 703   SetLineStyle(2);
 704   SetLineColor(8);
 705
 706   Double_t vect[7], vout[7];
 707   Double_t step = 1.0;
 708
 709   Float_t pv[3];
 710   pv[0] = fPart->Px();
 711   pv[1] = fPart->Py();
 712   pv[2] = fPart->Pz();
 713   fP.Set(pv);
 714
 715   vect[0] = fPart->Vx();
 716   vect[1] = fPart->Vy();
 717   vect[2] = fPart->Vz();
 718   vect[3] = fPart->Px()/fPart->P();
 719   vect[4] = fPart->Py()/fPart->P();
 720   vect[5] = fPart->Pz()/fPart->P();
 721   vect[6] = fPart->P();
 722
 723   TParticlePDG *ppdg = fPart->GetPDG(1);
 724   Int_t charge = (Int_t)(ppdg->Charge()/3.0);
 725
 726   Double_t zMax = -1750.0;
 727   Double_t rMax =   350.0;
 728   Double_t r    =     0.0;
 729
 730   Int_t nSteps = 0;
 731   while ((vect[2] > zMax) && (nSteps < 10000) && (r < rMax)) {
 732     nSteps++;
 733     OneStepRungekutta(charge, step, vect, vout);
 734     SetPoint(fCount,vout[0],vout[1],vout[2]);
 735     fCount++;
 736     for (Int_t i = 0; i < 7; i++) {
 737       vect[i] = vout[i];
 738     }
 739     r = TMath::Sqrt(vect[0]*vect[0]+vect[1]*vect[1]);
 740   }
 741
 742 }
 743
 744 //______________________________________________________________________________
 745 void AliEveMUONTrack::MakeRefTrack(AliMUONTrack *mtrack)
 746 {
 747   //
 748   // builds the track with dipole field starting from the TParticle
 749   //
 750
 751   fIsRefTrack = kTRUE;
 752
 753   char form[1000];
 754   snprintf(form,1000,"RefTrack %2d ", fLabel);
 755   SetName(form);
 756   SetLineStyle(2);
 757   SetLineColor(0);
 758
 759   MakeMUONTrack(mtrack);
 760
 761 }
 762
 763 //______________________________________________________________________________
 764 void AliEveMUONTrack::Propagate(Float_t *xr, Float_t *yr, Float_t *zr, Int_t i1, Int_t i2)
 765 {
 766   //
 767   // propagate in magnetic field between hits of indices i1 and i2
 768   //
 769
 770   Double_t vect[7], vout[7];
 771   Double_t step = 1.0;
 772   Double_t zMax = 0.0;
 773   Int_t  charge =   0;
 774   AliMUONTrackParam *trackParam = 0;
 775   TObjArray *trackParamAtCluster = 0;
 776
 777   if (i2 == 9999) {
 778     zMax = zr[i1]+1.5*step;
 779   } else {
 780     zMax = zr[i2]+1.5*step;
 781   }
 782
 783   trackParamAtCluster = fTrack->GetTrackParamAtCluster();
 784
 785   if (IsMUONTrack() || IsESDTrack() || IsRefTrack()) {
 786     trackParam = (AliMUONTrackParam*)trackParamAtCluster->At(i1);
 787     charge = (Int_t)trackParam->GetCharge();
 788   } else {
 789     return;
 790   }
 791
 792   vect[0] = xr[i1];
 793   vect[1] = yr[i1];
 794   vect[2] = zr[i1];
 795   vect[6] = trackParam->P();
 796   vect[3] = trackParam->Px()/vect[6];
 797   vect[4] = trackParam->Py()/vect[6];
 798   vect[5] = trackParam->Pz()/vect[6];
 799
 800   Int_t nSteps = 0;
 801   while ((vect[2] > zMax) && (nSteps < 10000)) {
 802     nSteps++;
 803     OneStepRungekutta(charge, step, vect, vout);
 804     SetPoint(fCount,vout[0],vout[1],vout[2]);
 805     fCount++;
 806     for (Int_t i = 0; i < 7; i++) {
 807       vect[i] = vout[i];
 808     }
 809   }
 810
 811 }
 812
 813 //______________________________________________________________________________
 814 void AliEveMUONTrack::OneStepRungekutta(Double_t charge, Double_t step,
 815                                   Double_t* vect, Double_t* vout)
 816 {
 817 ///     ******************************************************************
 818 ///     *                                                                *
 819 ///     *  Runge-Kutta method for tracking a particle through a magnetic *
 820 ///     *  field. Uses Nystroem algorithm (See Handbook Nat. Bur. of     *
 821 ///     *  Standards, procedure 25.5.20)                                 *
 822 ///     *                                                                *
 823 ///     *  Input parameters                                              *
 824 ///     *       CHARGE    Particle charge                                *
 825 ///     *       STEP      Step size                                      *
 826 ///     *       VECT      Initial co-ords,direction cosines,momentum     *
 827 ///     *  Output parameters                                             *
 828 ///     *       VOUT      Output co-ords,direction cosines,momentum      *
 829 ///     *  User routine called                                           *
 830 ///     *       CALL GUFLD(X,F)                                          *
 831 ///     *                                                                *
 832 ///     *    ==>Called by : <USER>, GUSWIM                               *
 833 ///     *       Authors    R.Brun, M.Hansroul  *********                 *
 834 ///     *                  V.Perevoztchikov (CUT STEP implementation)    *
 835 ///     *                                                                *
 836 ///     *                                                                *
 837 ///     ******************************************************************
 838
 839     Double_t h2, h4, f[4];
 840     Double_t xyzt[3], a, b, c, ph,ph2;
 841     Double_t secxs[4],secys[4],seczs[4],hxp[3];
 842     Double_t g1, g2, g3, g4, g5, g6, ang2, dxt, dyt, dzt;
 843     Double_t est, at, bt, ct, cba;
 844     Double_t f1, f2, f3, f4, rho, tet, hnorm, hp, rho1, sint, cost;
 845
 846     Double_t x;
 847     Double_t y;
 848     Double_t z;
 849
 850     Double_t xt;
 851     Double_t yt;
 852     Double_t zt;
 853
 854     Double_t maxit = 1992;
 855     Double_t maxcut = 11;
 856
 857     const Double_t kdlt   = 1e-4;
 858     const Double_t kdlt32 = kdlt/32.;
 859     const Double_t kthird = 1./3.;
 860     const Double_t khalf  = 0.5;
 861     const Double_t kec = 2.9979251e-4;
 862
 863     const Double_t kpisqua = 9.86960440109;
 864     const Int_t kix  = 0;
 865     const Int_t kiy  = 1;
 866     const Int_t kiz  = 2;
 867     const Int_t kipx = 3;
 868     const Int_t kipy = 4;
 869     const Int_t kipz = 5;
 870
 871     // *.
 872     // *.    ------------------------------------------------------------------
 873     // *.
 874     // *             this constant is for units cm,gev/c and kgauss
 875     // *
 876     Int_t iter = 0;
 877     Int_t ncut = 0;
 878     for(Int_t j = 0; j < 7; j++)
 879       vout[j] = vect[j];
 880
 881     Double_t  pinv   = kec * charge / vect[6];
 882     Double_t tl = 0.;
 883     Double_t h = step;
 884     Double_t rest;
 885
 886
 887     do {
 888       rest  = step - tl;
 889       if (TMath::Abs(h) > TMath::Abs(rest)) h = rest;
 890       //cmodif: call gufld(vout,f) changed into:
 891       TGeoGlobalMagField::Instance()->Field(vout,f);
 892
 893       // *
 894       // *             start of integration
 895       // *
 896       x      = vout[0];
 897       y      = vout[1];
 898       z      = vout[2];
 899       a      = vout[3];
 900       b      = vout[4];
 901       c      = vout[5];
 902
 903       h2     = khalf * h;
 904       h4     = khalf * h2;
 905       ph     = pinv * h;
 906       ph2    = khalf * ph;
 907       secxs[0] = (b * f[2] - c * f[1]) * ph2;
 908       secys[0] = (c * f[0] - a * f[2]) * ph2;
 909       seczs[0] = (a * f[1] - b * f[0]) * ph2;
 910       ang2 = (secxs[0]*secxs[0] + secys[0]*secys[0] + seczs[0]*seczs[0]);
 911       if (ang2 > kpisqua) break;
 912
 913       dxt    = h2 * a + h4 * secxs[0];
 914       dyt    = h2 * b + h4 * secys[0];
 915       dzt    = h2 * c + h4 * seczs[0];
 916       xt     = x + dxt;
 917       yt     = y + dyt;
 918       zt     = z + dzt;
 919       // *
 920       // *              second intermediate point
 921       // *
 922
 923       est = TMath::Abs(dxt) + TMath::Abs(dyt) + TMath::Abs(dzt);
 924       if (est > h) {
 925         if (ncut++ > maxcut) break;
 926         h *= khalf;
 927         continue;
 928       }
 929
 930       xyzt[0] = xt;
 931       xyzt[1] = yt;
 932       xyzt[2] = zt;
 933
 934       //cmodif: call gufld(xyzt,f) changed into:
 935       TGeoGlobalMagField::Instance()->Field(xyzt,f);
 936
 937       at     = a + secxs[0];
 938       bt     = b + secys[0];
 939       ct     = c + seczs[0];
 940
 941       secxs[1] = (bt * f[2] - ct * f[1]) * ph2;
 942       secys[1] = (ct * f[0] - at * f[2]) * ph2;
 943       seczs[1] = (at * f[1] - bt * f[0]) * ph2;
 944       at     = a + secxs[1];
 945       bt     = b + secys[1];
 946       ct     = c + seczs[1];
 947       secxs[2] = (bt * f[2] - ct * f[1]) * ph2;
 948       secys[2] = (ct * f[0] - at * f[2]) * ph2;
 949       seczs[2] = (at * f[1] - bt * f[0]) * ph2;
 950       dxt    = h * (a + secxs[2]);
 951       dyt    = h * (b + secys[2]);
 952       dzt    = h * (c + seczs[2]);
 953       xt     = x + dxt;
 954       yt     = y + dyt;
 955       zt     = z + dzt;
 956       at     = a + 2.*secxs[2];
 957       bt     = b + 2.*secys[2];
 958       ct     = c + 2.*seczs[2];
 959
 960       est = TMath::Abs(dxt)+TMath::Abs(dyt)+TMath::Abs(dzt);
 961       if (est > 2.*TMath::Abs(h)) {
 962         if (ncut++ > maxcut) break;
 963         h *= khalf;
 964         continue;
 965       }
 966
 967       xyzt[0] = xt;
 968       xyzt[1] = yt;
 969       xyzt[2] = zt;
 970
 971       //cmodif: call gufld(xyzt,f) changed into:
 972       TGeoGlobalMagField::Instance()->Field(xyzt,f);
 973
 974       z      = z + (c + (seczs[0] + seczs[1] + seczs[2]) * kthird) * h;
 975       y      = y + (b + (secys[0] + secys[1] + secys[2]) * kthird) * h;
 976       x      = x + (a + (secxs[0] + secxs[1] + secxs[2]) * kthird) * h;
 977
 978       secxs[3] = (bt*f[2] - ct*f[1])* ph2;
 979       secys[3] = (ct*f[0] - at*f[2])* ph2;
 980       seczs[3] = (at*f[1] - bt*f[0])* ph2;
 981       a      = a+(secxs[0]+secxs[3]+2. * (secxs[1]+secxs[2])) * kthird;
 982       b      = b+(secys[0]+secys[3]+2. * (secys[1]+secys[2])) * kthird;
 983       c      = c+(seczs[0]+seczs[3]+2. * (seczs[1]+seczs[2])) * kthird;
 984
 985       est    = TMath::Abs(secxs[0]+secxs[3] - (secxs[1]+secxs[2]))
 986         + TMath::Abs(secys[0]+secys[3] - (secys[1]+secys[2]))
 987         + TMath::Abs(seczs[0]+seczs[3] - (seczs[1]+seczs[2]));
 988
 989       if (est > kdlt && TMath::Abs(h) > 1.e-4) {
 990         if (ncut++ > maxcut) break;
 991         h *= khalf;
 992         continue;
 993       }
 994
 995       ncut = 0;
 996       // *               if too many iterations, go to helix
 997       if (iter++ > maxit) break;
 998
 999       tl += h;
1000       if (est < kdlt32)
1001         h *= 2.;
1002       cba    = 1./ TMath::Sqrt(a*a + b*b + c*c);
1003       vout[0] = x;
1004       vout[1] = y;
1005       vout[2] = z;
1006       vout[3] = cba*a;
1007       vout[4] = cba*b;
1008       vout[5] = cba*c;
1009       rest = step - tl;
1010       if (step < 0.) rest = -rest;
1011       if (rest < 1.e-5*TMath::Abs(step)) return;
1012
1013     } while(1);
1014
1015     // angle too big, use helix
1016
1017     f1  = f[0];
1018     f2  = f[1];
1019     f3  = f[2];
1020     f4  = TMath::Sqrt(f1*f1+f2*f2+f3*f3);
1021     rho = -f4*pinv;
1022     tet = rho * step;
1023
1024     hnorm = 1./f4;
1025     f1 = f1*hnorm;
1026     f2 = f2*hnorm;
1027     f3 = f3*hnorm;
1028
1029     hxp[0] = f2*vect[kipz] - f3*vect[kipy];
1030     hxp[1] = f3*vect[kipx] - f1*vect[kipz];
1031     hxp[2] = f1*vect[kipy] - f2*vect[kipx];
1032
1033     hp = f1*vect[kipx] + f2*vect[kipy] + f3*vect[kipz];
1034
1035     rho1 = 1./rho;
1036     sint = TMath::Sin(tet);
1037     cost = 2.*TMath::Sin(khalf*tet)*TMath::Sin(khalf*tet);
1038
1039     g1 = sint*rho1;
1040     g2 = cost*rho1;
1041     g3 = (tet-sint) * hp*rho1;
1042     g4 = -cost;
1043     g5 = sint;
1044     g6 = cost * hp;
1045
1046     vout[kix] = vect[kix] + g1*vect[kipx] + g2*hxp[0] + g3*f1;
1047     vout[kiy] = vect[kiy] + g1*vect[kipy] + g2*hxp[1] + g3*f2;
1048     vout[kiz] = vect[kiz] + g1*vect[kipz] + g2*hxp[2] + g3*f3;
1049
1050     vout[kipx] = vect[kipx] + g4*vect[kipx] + g5*hxp[0] + g6*f1;
1051     vout[kipy] = vect[kipy] + g4*vect[kipy] + g5*hxp[1] + g6*f2;
1052     vout[kipz] = vect[kipz] + g4*vect[kipz] + g5*hxp[2] + g6*f3;
1053
1054     return;
1055 }
1056
1057 //______________________________________________________________________________
1058 Int_t AliEveMUONTrack::ColorIndex(Float_t val)
1059 {
1060   //
1061   // returns color index in the palette for a give value
1062   //
1063
1064   Float_t threshold =  0.0;
1065   Float_t maxVal    = 10.0;
1066
1067   Float_t div  = TMath::Max(1, (Int_t)(maxVal - threshold));
1068   Int_t   nCol = gStyle->GetNumberOfColors();
1069   Int_t   cBin = (Int_t) TMath::Nint(nCol*(val - threshold)/div);
1070
1071   return gStyle->GetColorPalette(TMath::Min(nCol - 1, cBin));
1072
1073 }
1074
1075 //==============================================================================
1076 // Temporary AliEveMUONTrackList
1077 //==============================================================================
1078
1079 //______________________________________________________________________________
1080 void AliEveMUONTrackList::HackMomentumLimits(Bool_t recurse)
1081 {
1082   // Find momentum limits from included tracks.
1083
1084   fLimPt = fLimP = 0;
1085
1086   for (List_i i=BeginChildren(); i!=EndChildren(); ++i)
1087   {
1088     TEveTrack* track = dynamic_cast<TEveTrack*>(*i);
1089     if (track)
1090     {
1091       fLimPt = TMath::Max(fLimPt, track->GetMomentum().Perp());
1092       fLimP  = TMath::Max(fLimP,  track->GetMomentum().Mag());
1093     }
1094     if (recurse)
1095       FindMomentumLimits(*i, recurse);
1096   }
1097
1098   fLimPt = RoundMomentumLimit(fLimPt);
1099   fLimP  = RoundMomentumLimit(fLimP);
1100   if (fMaxPt == 0) fMaxPt = fLimPt;
1101   if (fMaxP  == 0) fMaxP  = fLimP;
1102 }