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