EVE/EveDet/AliEveMUONTrack.cxx

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