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