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