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