EVE/Alieve/MUONTrack.cxx

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