EVE/Alieve/AliEveMUONTrack.cxx

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