]>
Commit | Line | Data |
---|---|---|
d810d0de | 1 | // $Id$ |
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> |
18 | #include <AliTrackReference.h> | |
af885e0f | 19 | #include <AliESDEvent.h> |
82e1eece | 20 | #include <AliESDVertex.h> |
eadce74d | 21 | #include <AliRunLoader.h> |
22 | #include <AliRun.h> | |
e19f3684 | 23 | |
24 | #include <AliMUONTrack.h> | |
eadce74d | 25 | #include <AliMUONTriggerTrack.h> |
e19f3684 | 26 | #include <AliMUONTrackParam.h> |
27 | #include <AliMUONConstants.h> | |
103e6575 | 28 | #include <AliMUONESDInterface.h> |
e19f3684 | 29 | |
30 | #include <TClonesArray.h> | |
31 | #include <TMath.h> | |
32 | #include <TMatrixD.h> | |
33 | #include <TStyle.h> | |
f7e76528 | 34 | #include <TROOT.h> |
eadce74d | 35 | #include <TParticle.h> |
36 | #include <TParticlePDG.h> | |
37 | ||
38 | #include <Riostream.h> | |
d810d0de | 39 | |
e19f3684 | 40 | |
57ffa5fb | 41 | //______________________________________________________________________________ |
d810d0de | 42 | // AliEveMUONTrack |
84aff7a4 | 43 | // Produce TEveUtil:TEveTrack from AliMUONTrack with dipole field model |
e19f3684 | 44 | |
d810d0de | 45 | ClassImp(AliEveMUONTrack) |
e19f3684 | 46 | |
d810d0de | 47 | AliMagF* AliEveMUONTrack::fFieldMap = 0; |
eadce74d | 48 | |
57ffa5fb | 49 | //______________________________________________________________________________ |
d810d0de | 50 | AliEveMUONTrack::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 | ||
d810d0de | 65 | fFieldMap = AliEveEventManager::AssertMagField(); |
eadce74d | 66 | |
e19f3684 | 67 | } |
68 | ||
57ffa5fb | 69 | //______________________________________________________________________________ |
d810d0de | 70 | AliEveMUONTrack::~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 | 82 | void 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 | 133 | void AliEveMUONTrack::PrintMUONTrackInfo() |
eadce74d | 134 | { |
135 | // | |
136 | // information about the reconstructed/reference track; at hits and at vertex | |
137 | // | |
138 | ||
139 | Double_t RADDEG = 180.0/TMath::Pi(); | |
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) << | |
51346b82 | 192 | mtp->GetBendingSlope()*RADDEG << " " << |
eadce74d | 193 | |
194 | setw(8) << setprecision(3) << | |
51346b82 | 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) << | |
241 | mtp->GetBendingSlope()*RADDEG << " " << | |
242 | ||
eadce74d | 243 | setw(8) << setprecision(3) << |
51346b82 | 244 | mtp->GetNonBendingSlope()*RADDEG << " " << |
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 | 279 | void AliEveMUONTrack::PrintMUONTriggerTrackInfo() |
eadce74d | 280 | { |
281 | // | |
282 | // information about the trigger track | |
283 | // | |
284 | ||
285 | // Double_t RADDEG = 180.0/TMath::Pi(); | |
286 | ||
287 | } | |
288 | ||
57ffa5fb | 289 | //______________________________________________________________________________ |
d810d0de | 290 | void AliEveMUONTrack::PrintESDTrackInfo() |
eadce74d | 291 | { |
292 | // | |
293 | // information about the reconstructed ESD track at vertex | |
294 | // | |
295 | ||
296 | Double_t RADDEG = 180.0/TMath::Pi(); | |
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) << |
51346b82 | 312 | mtp->GetBendingSlope()*RADDEG << " " << |
313 | ||
eadce74d | 314 | setw(8) << setprecision(3) << |
51346b82 | 315 | mtp->GetNonBendingSlope()*RADDEG << " " << |
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 | 396 | void 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 | 426 | void 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 | 468 | void 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 | 749 | void 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 | 779 | void 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 | 838 | void 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 | 893 | void 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 | 912 | void 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 | 965 | void 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 | ||
eadce74d | 976 | if (fFieldMap) { |
977 | fFieldMap->Field(x,b); | |
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 | 996 | void 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 | 1241 | Int_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 | } |