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