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