a9e2aefa |
1 | /************************************************************************** |
2 | * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. * |
3 | * * |
4 | * Author: The ALICE Off-line Project. * |
5 | * Contributors are mentioned in the code where appropriate. * |
6 | * * |
7 | * Permission to use, copy, modify and distribute this software and its * |
8 | * documentation strictly for non-commercial purposes is hereby granted * |
9 | * without fee, provided that the above copyright notice appears in all * |
10 | * copies and that both the copyright notice and this permission notice * |
11 | * appear in the supporting documentation. The authors make no claims * |
12 | * about the suitability of this software for any purpose. It is * |
13 | * provided "as is" without express or implied warranty. * |
14 | **************************************************************************/ |
15 | |
88cb7938 |
16 | /* $Id$ */ |
a9e2aefa |
17 | |
3831f268 |
18 | /////////////////////////////////////////////////// |
19 | // |
20 | // Track parameters |
21 | // in |
22 | // ALICE |
23 | // dimuon |
24 | // spectrometer |
a9e2aefa |
25 | // |
3831f268 |
26 | /////////////////////////////////////////////////// |
a9e2aefa |
27 | |
70479d0e |
28 | #include <Riostream.h> |
3831f268 |
29 | #include "AliMUON.h" |
a9e2aefa |
30 | #include "AliMUONTrackParam.h" |
3831f268 |
31 | #include "AliMUONChamber.h" |
a9e2aefa |
32 | #include "AliRun.h" |
94de3818 |
33 | #include "AliMagF.h" |
8c343c7c |
34 | #include "AliLog.h" |
a9e2aefa |
35 | |
36 | ClassImp(AliMUONTrackParam) // Class implementation in ROOT context |
37 | |
61adb9bd |
38 | //_________________________________________________________________________ |
30178c30 |
39 | AliMUONTrackParam::AliMUONTrackParam() |
40 | : TObject() |
41 | { |
42 | // Constructor |
43 | |
44 | fInverseBendingMomentum = 0; |
45 | fBendingSlope = 0; |
46 | fNonBendingSlope = 0; |
47 | fZ = 0; |
48 | fBendingCoor = 0; |
49 | fNonBendingCoor = 0; |
50 | } |
61adb9bd |
51 | |
30178c30 |
52 | //_________________________________________________________________________ |
53 | AliMUONTrackParam& |
54 | AliMUONTrackParam::operator=(const AliMUONTrackParam& theMUONTrackParam) |
61adb9bd |
55 | { |
30178c30 |
56 | if (this == &theMUONTrackParam) |
61adb9bd |
57 | return *this; |
58 | |
30178c30 |
59 | // base class assignement |
60 | TObject::operator=(theMUONTrackParam); |
61 | |
62 | fInverseBendingMomentum = theMUONTrackParam.fInverseBendingMomentum; |
63 | fBendingSlope = theMUONTrackParam.fBendingSlope; |
64 | fNonBendingSlope = theMUONTrackParam.fNonBendingSlope; |
65 | fZ = theMUONTrackParam.fZ; |
66 | fBendingCoor = theMUONTrackParam.fBendingCoor; |
67 | fNonBendingCoor = theMUONTrackParam.fNonBendingCoor; |
61adb9bd |
68 | |
69 | return *this; |
70 | } |
71 | //_________________________________________________________________________ |
30178c30 |
72 | AliMUONTrackParam::AliMUONTrackParam(const AliMUONTrackParam& theMUONTrackParam) |
73 | : TObject(theMUONTrackParam) |
61adb9bd |
74 | { |
30178c30 |
75 | fInverseBendingMomentum = theMUONTrackParam.fInverseBendingMomentum; |
76 | fBendingSlope = theMUONTrackParam.fBendingSlope; |
77 | fNonBendingSlope = theMUONTrackParam.fNonBendingSlope; |
78 | fZ = theMUONTrackParam.fZ; |
79 | fBendingCoor = theMUONTrackParam.fBendingCoor; |
80 | fNonBendingCoor = theMUONTrackParam.fNonBendingCoor; |
61adb9bd |
81 | } |
a9e2aefa |
82 | |
a9e2aefa |
83 | //__________________________________________________________________________ |
84 | void AliMUONTrackParam::ExtrapToZ(Double_t Z) |
85 | { |
86 | // Track parameter extrapolation to the plane at "Z". |
87 | // On return, the track parameters resulting from the extrapolation |
88 | // replace the current track parameters. |
a9e2aefa |
89 | if (this->fZ == Z) return; // nothing to be done if same Z |
90 | Double_t forwardBackward; // +1 if forward, -1 if backward |
5b64e914 |
91 | if (Z < this->fZ) forwardBackward = 1.0; // spectro. z<0 |
a9e2aefa |
92 | else forwardBackward = -1.0; |
a6f03ddb |
93 | Double_t vGeant3[7], vGeant3New[7]; // 7 in parameter ???? |
a9e2aefa |
94 | Int_t iGeant3, stepNumber; |
95 | Int_t maxStepNumber = 5000; // in parameter ???? |
96 | // For safety: return kTRUE or kFALSE ???? |
a6f03ddb |
97 | // Parameter vector for calling EXTRAP_ONESTEP |
a9e2aefa |
98 | SetGeant3Parameters(vGeant3, forwardBackward); |
956019b6 |
99 | // sign of charge (sign of fInverseBendingMomentum if forward motion) |
a6f03ddb |
100 | // must be changed if backward extrapolation |
956019b6 |
101 | Double_t chargeExtrap = forwardBackward * |
102 | TMath::Sign(Double_t(1.0), this->fInverseBendingMomentum); |
a9e2aefa |
103 | Double_t stepLength = 6.0; // in parameter ???? |
104 | // Extrapolation loop |
105 | stepNumber = 0; |
5b64e914 |
106 | while (((-forwardBackward * (vGeant3[2] - Z)) <= 0.0) && // spectro. z<0 |
a9e2aefa |
107 | (stepNumber < maxStepNumber)) { |
108 | stepNumber++; |
a6f03ddb |
109 | // Option for switching between helix and Runge-Kutta ???? |
4d03a78e |
110 | //ExtrapOneStepRungekutta(chargeExtrap, stepLength, vGeant3, vGeant3New); |
111 | ExtrapOneStepHelix(chargeExtrap, stepLength, vGeant3, vGeant3New); |
5b64e914 |
112 | if ((-forwardBackward * (vGeant3New[2] - Z)) > 0.0) break; // one is beyond Z spectro. z<0 |
a9e2aefa |
113 | // better use TArray ???? |
114 | for (iGeant3 = 0; iGeant3 < 7; iGeant3++) |
115 | {vGeant3[iGeant3] = vGeant3New[iGeant3];} |
116 | } |
117 | // check maxStepNumber ???? |
a9e2aefa |
118 | // Interpolation back to exact Z (2nd order) |
119 | // should be in function ???? using TArray ???? |
120 | Double_t dZ12 = vGeant3New[2] - vGeant3[2]; // 1->2 |
121 | Double_t dZ1i = Z - vGeant3[2]; // 1-i |
122 | Double_t dZi2 = vGeant3New[2] - Z; // i->2 |
123 | Double_t xPrime = (vGeant3New[0] - vGeant3[0]) / dZ12; |
124 | Double_t xSecond = |
125 | ((vGeant3New[3] / vGeant3New[5]) - (vGeant3[3] / vGeant3[5])) / dZ12; |
126 | Double_t yPrime = (vGeant3New[1] - vGeant3[1]) / dZ12; |
127 | Double_t ySecond = |
128 | ((vGeant3New[4] / vGeant3New[5]) - (vGeant3[4] / vGeant3[5])) / dZ12; |
129 | vGeant3[0] = vGeant3[0] + xPrime * dZ1i - 0.5 * xSecond * dZ1i * dZi2; // X |
130 | vGeant3[1] = vGeant3[1] + yPrime * dZ1i - 0.5 * ySecond * dZ1i * dZi2; // Y |
131 | vGeant3[2] = Z; // Z |
132 | Double_t xPrimeI = xPrime - 0.5 * xSecond * (dZi2 - dZ1i); |
133 | Double_t yPrimeI = yPrime - 0.5 * ySecond * (dZi2 - dZ1i); |
956019b6 |
134 | // (PX, PY, PZ)/PTOT assuming forward motion |
a9e2aefa |
135 | vGeant3[5] = |
136 | 1.0 / TMath::Sqrt(1.0 + xPrimeI * xPrimeI + yPrimeI * yPrimeI); // PZ/PTOT |
137 | vGeant3[3] = xPrimeI * vGeant3[5]; // PX/PTOT |
138 | vGeant3[4] = yPrimeI * vGeant3[5]; // PY/PTOT |
956019b6 |
139 | // Track parameters from Geant3 parameters, |
140 | // with charge back for forward motion |
141 | GetFromGeant3Parameters(vGeant3, chargeExtrap * forwardBackward); |
a9e2aefa |
142 | } |
143 | |
144 | //__________________________________________________________________________ |
145 | void AliMUONTrackParam::SetGeant3Parameters(Double_t *VGeant3, Double_t ForwardBackward) |
146 | { |
147 | // Set vector of Geant3 parameters pointed to by "VGeant3" |
148 | // from track parameters in current AliMUONTrackParam. |
149 | // Since AliMUONTrackParam is only geometry, one uses "ForwardBackward" |
150 | // to know whether the particle is going forward (+1) or backward (-1). |
151 | VGeant3[0] = this->fNonBendingCoor; // X |
152 | VGeant3[1] = this->fBendingCoor; // Y |
153 | VGeant3[2] = this->fZ; // Z |
154 | Double_t pYZ = TMath::Abs(1.0 / this->fInverseBendingMomentum); |
155 | Double_t pZ = |
156 | pYZ / TMath::Sqrt(1.0 + this->fBendingSlope * this->fBendingSlope); |
157 | VGeant3[6] = |
158 | TMath::Sqrt(pYZ * pYZ + |
159 | pZ * pZ * this->fNonBendingSlope * this->fNonBendingSlope); // PTOT |
5b64e914 |
160 | VGeant3[5] = -ForwardBackward * pZ / VGeant3[6]; // PZ/PTOT spectro. z<0 |
a9e2aefa |
161 | VGeant3[3] = this->fNonBendingSlope * VGeant3[5]; // PX/PTOT |
162 | VGeant3[4] = this->fBendingSlope * VGeant3[5]; // PY/PTOT |
163 | } |
164 | |
165 | //__________________________________________________________________________ |
166 | void AliMUONTrackParam::GetFromGeant3Parameters(Double_t *VGeant3, Double_t Charge) |
167 | { |
168 | // Get track parameters in current AliMUONTrackParam |
956019b6 |
169 | // from Geant3 parameters pointed to by "VGeant3", |
170 | // assumed to be calculated for forward motion in Z. |
a9e2aefa |
171 | // "InverseBendingMomentum" is signed with "Charge". |
172 | this->fNonBendingCoor = VGeant3[0]; // X |
173 | this->fBendingCoor = VGeant3[1]; // Y |
174 | this->fZ = VGeant3[2]; // Z |
175 | Double_t pYZ = VGeant3[6] * TMath::Sqrt(1.0 - VGeant3[3] * VGeant3[3]); |
176 | this->fInverseBendingMomentum = Charge / pYZ; |
177 | this->fBendingSlope = VGeant3[4] / VGeant3[5]; |
178 | this->fNonBendingSlope = VGeant3[3] / VGeant3[5]; |
179 | } |
180 | |
181 | //__________________________________________________________________________ |
182 | void AliMUONTrackParam::ExtrapToStation(Int_t Station, AliMUONTrackParam *TrackParam) |
183 | { |
184 | // Track parameters extrapolated from current track parameters ("this") |
185 | // to both chambers of the station(0..) "Station" |
186 | // are returned in the array (dimension 2) of track parameters |
187 | // pointed to by "TrackParam" (index 0 and 1 for first and second chambers). |
188 | Double_t extZ[2], z1, z2; |
ecfa008b |
189 | Int_t i1 = -1, i2 = -1; // = -1 to avoid compilation warnings |
a9e2aefa |
190 | AliMUON *pMUON = (AliMUON*) gAlice->GetModule("MUON"); // necessary ???? |
191 | // range of Station to be checked ???? |
192 | z1 = (&(pMUON->Chamber(2 * Station)))->Z(); // Z of first chamber |
193 | z2 = (&(pMUON->Chamber(2 * Station + 1)))->Z(); // Z of second chamber |
194 | // First and second Z to extrapolate at |
195 | if ((z1 > this->fZ) && (z2 > this->fZ)) {i1 = 0; i2 = 1;} |
196 | else if ((z1 < this->fZ) && (z2 < this->fZ)) {i1 = 1; i2 = 0;} |
197 | else { |
8c343c7c |
198 | AliError(Form("Starting Z (%f) in between z1 (%f) and z2 (%f) of station(0..)%d",this->fZ,z1,z2,Station)); |
199 | // cout << "ERROR in AliMUONTrackParam::CreateExtrapSegmentInStation" << endl; |
200 | // cout << "Starting Z (" << this->fZ << ") in between z1 (" << z1 << |
201 | // ") and z2 (" << z2 << ") of station(0..) " << Station << endl; |
a9e2aefa |
202 | } |
203 | extZ[i1] = z1; |
204 | extZ[i2] = z2; |
205 | // copy of track parameters |
206 | TrackParam[i1] = *this; |
207 | // first extrapolation |
208 | (&(TrackParam[i1]))->ExtrapToZ(extZ[0]); |
209 | TrackParam[i2] = TrackParam[i1]; |
210 | // second extrapolation |
211 | (&(TrackParam[i2]))->ExtrapToZ(extZ[1]); |
212 | return; |
213 | } |
214 | |
04b5ea16 |
215 | //__________________________________________________________________________ |
889a0215 |
216 | void AliMUONTrackParam::ExtrapToVertex(Double_t xVtx, Double_t yVtx, Double_t zVtx) |
04b5ea16 |
217 | { |
218 | // Extrapolation to the vertex. |
219 | // Returns the track parameters resulting from the extrapolation, |
220 | // in the current TrackParam. |
956019b6 |
221 | // Changes parameters according to Branson correction through the absorber |
04b5ea16 |
222 | |
5b64e914 |
223 | Double_t zAbsorber = -503.0; // to be coherent with the Geant absorber geometry !!!! |
224 | // spectro. (z<0) |
04b5ea16 |
225 | // Extrapolates track parameters upstream to the "Z" end of the front absorber |
b45fd22b |
226 | ExtrapToZ(zAbsorber); // !!! |
5b64e914 |
227 | // Makes Branson correction (multiple scattering + energy loss) |
889a0215 |
228 | BransonCorrection(xVtx,yVtx,zVtx); |
5b64e914 |
229 | // Makes a simple magnetic field correction through the absorber |
b45fd22b |
230 | FieldCorrection(zAbsorber); |
04b5ea16 |
231 | } |
232 | |
43af2cb6 |
233 | |
234 | // Keep this version for future developments |
04b5ea16 |
235 | //__________________________________________________________________________ |
43af2cb6 |
236 | // void AliMUONTrackParam::BransonCorrection() |
237 | // { |
238 | // // Branson correction of track parameters |
239 | // // the entry parameters have to be calculated at the end of the absorber |
240 | // Double_t zEndAbsorber, zBP, xBP, yBP; |
241 | // Double_t pYZ, pX, pY, pZ, pTotal, xEndAbsorber, yEndAbsorber, radiusEndAbsorber2, pT, theta; |
242 | // Int_t sign; |
243 | // // Would it be possible to calculate all that from Geant configuration ???? |
244 | // // and to get the Branson parameters from a function in ABSO module ???? |
245 | // // with an eventual contribution from other detectors like START ???? |
246 | // // Radiation lengths outer part theta > 3 degres |
247 | // static Double_t x01[9] = { 18.8, // C (cm) |
248 | // 10.397, // Concrete (cm) |
249 | // 0.56, // Plomb (cm) |
250 | // 47.26, // Polyethylene (cm) |
251 | // 0.56, // Plomb (cm) |
252 | // 47.26, // Polyethylene (cm) |
253 | // 0.56, // Plomb (cm) |
254 | // 47.26, // Polyethylene (cm) |
255 | // 0.56 }; // Plomb (cm) |
256 | // // inner part theta < 3 degres |
257 | // static Double_t x02[3] = { 18.8, // C (cm) |
258 | // 10.397, // Concrete (cm) |
259 | // 0.35 }; // W (cm) |
260 | // // z positions of the materials inside the absober outer part theta > 3 degres |
261 | // static Double_t z1[10] = { 90, 315, 467, 472, 477, 482, 487, 492, 497, 502 }; |
262 | // // inner part theta < 3 degres |
263 | // static Double_t z2[4] = { 90, 315, 467, 503 }; |
264 | // static Bool_t first = kTRUE; |
265 | // static Double_t zBP1, zBP2, rLimit; |
266 | // // Calculates z positions of the Branson's planes: zBP1 for outer part and zBP2 for inner part (only at the first call) |
267 | // if (first) { |
268 | // first = kFALSE; |
269 | // Double_t aNBP = 0.0; |
270 | // Double_t aDBP = 0.0; |
271 | // Int_t iBound; |
272 | |
273 | // for (iBound = 0; iBound < 9; iBound++) { |
274 | // aNBP = aNBP + |
275 | // (z1[iBound+1] * z1[iBound+1] * z1[iBound+1] - |
276 | // z1[iBound] * z1[iBound] * z1[iBound] ) / x01[iBound]; |
277 | // aDBP = aDBP + |
278 | // (z1[iBound+1] * z1[iBound+1] - z1[iBound] * z1[iBound] ) / x01[iBound]; |
279 | // } |
280 | // zBP1 = (2.0 * aNBP) / (3.0 * aDBP); |
281 | // aNBP = 0.0; |
282 | // aDBP = 0.0; |
283 | // for (iBound = 0; iBound < 3; iBound++) { |
284 | // aNBP = aNBP + |
285 | // (z2[iBound+1] * z2[iBound+1] * z2[iBound+1] - |
286 | // z2[iBound] * z2[iBound ] * z2[iBound] ) / x02[iBound]; |
287 | // aDBP = aDBP + |
288 | // (z2[iBound+1] * z2[iBound+1] - z2[iBound] * z2[iBound]) / x02[iBound]; |
289 | // } |
290 | // zBP2 = (2.0 * aNBP) / (3.0 * aDBP); |
291 | // rLimit = z2[3] * TMath::Tan(3.0 * (TMath::Pi()) / 180.); |
292 | // } |
293 | |
294 | // pYZ = TMath::Abs(1.0 / fInverseBendingMomentum); |
295 | // sign = 1; |
296 | // if (fInverseBendingMomentum < 0) sign = -1; |
297 | // pZ = pYZ / (TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope)); |
298 | // pX = pZ * fNonBendingSlope; |
299 | // pY = pZ * fBendingSlope; |
300 | // pTotal = TMath::Sqrt(pYZ *pYZ + pX * pX); |
301 | // xEndAbsorber = fNonBendingCoor; |
302 | // yEndAbsorber = fBendingCoor; |
303 | // radiusEndAbsorber2 = xEndAbsorber * xEndAbsorber + yEndAbsorber * yEndAbsorber; |
304 | |
305 | // if (radiusEndAbsorber2 > rLimit*rLimit) { |
306 | // zEndAbsorber = z1[9]; |
307 | // zBP = zBP1; |
308 | // } else { |
309 | // zEndAbsorber = z2[3]; |
310 | // zBP = zBP2; |
311 | // } |
312 | |
313 | // xBP = xEndAbsorber - (pX / pZ) * (zEndAbsorber - zBP); |
314 | // yBP = yEndAbsorber - (pY / pZ) * (zEndAbsorber - zBP); |
315 | |
316 | // // new parameters after Branson and energy loss corrections |
317 | // pZ = pTotal * zBP / TMath::Sqrt(xBP * xBP + yBP * yBP + zBP * zBP); |
318 | // pX = pZ * xBP / zBP; |
319 | // pY = pZ * yBP / zBP; |
320 | // fBendingSlope = pY / pZ; |
321 | // fNonBendingSlope = pX / pZ; |
322 | |
323 | // pT = TMath::Sqrt(pX * pX + pY * pY); |
324 | // theta = TMath::ATan2(pT, pZ); |
325 | // pTotal = |
326 | // TotalMomentumEnergyLoss(rLimit, pTotal, theta, xEndAbsorber, yEndAbsorber); |
327 | |
328 | // fInverseBendingMomentum = (sign / pTotal) * |
329 | // TMath::Sqrt(1.0 + |
330 | // fBendingSlope * fBendingSlope + |
331 | // fNonBendingSlope * fNonBendingSlope) / |
332 | // TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope); |
333 | |
334 | // // vertex position at (0,0,0) |
335 | // // should be taken from vertex measurement ??? |
336 | // fBendingCoor = 0.0; |
337 | // fNonBendingCoor = 0; |
338 | // fZ= 0; |
339 | // } |
340 | |
889a0215 |
341 | void AliMUONTrackParam::BransonCorrection(Double_t xVtx,Double_t yVtx,Double_t zVtx) |
04b5ea16 |
342 | { |
343 | // Branson correction of track parameters |
344 | // the entry parameters have to be calculated at the end of the absorber |
43af2cb6 |
345 | // simplified version: the z positions of Branson's planes are no longer calculated |
346 | // but are given as inputs. One can use the macros MUONTestAbso.C and DrawTestAbso.C |
347 | // to test this correction. |
04b5ea16 |
348 | // Would it be possible to calculate all that from Geant configuration ???? |
956019b6 |
349 | // and to get the Branson parameters from a function in ABSO module ???? |
350 | // with an eventual contribution from other detectors like START ???? |
889a0215 |
351 | //change to take into account the vertex postition (real, reconstruct,....) |
352 | |
43af2cb6 |
353 | Double_t zBP, xBP, yBP; |
354 | Double_t pYZ, pX, pY, pZ, pTotal, xEndAbsorber, yEndAbsorber, radiusEndAbsorber2, pT, theta; |
355 | Int_t sign; |
04b5ea16 |
356 | static Bool_t first = kTRUE; |
b45fd22b |
357 | static Double_t zBP1, zBP2, rLimit, thetaLimit, zEndAbsorber; |
43af2cb6 |
358 | // zBP1 for outer part and zBP2 for inner part (only at the first call) |
04b5ea16 |
359 | if (first) { |
360 | first = kFALSE; |
43af2cb6 |
361 | |
5b64e914 |
362 | zEndAbsorber = -503; // spectro (z<0) |
b45fd22b |
363 | thetaLimit = 3.0 * (TMath::Pi()) / 180.; |
5b64e914 |
364 | rLimit = TMath::Abs(zEndAbsorber) * TMath::Tan(thetaLimit); |
365 | zBP1 = -450; // values close to those calculated with EvalAbso.C |
366 | zBP2 = -480; |
04b5ea16 |
367 | } |
368 | |
369 | pYZ = TMath::Abs(1.0 / fInverseBendingMomentum); |
370 | sign = 1; |
b8dc484b |
371 | if (fInverseBendingMomentum < 0) sign = -1; |
372 | pZ = Pz(); |
373 | pX = Px(); |
374 | pY = Py(); |
04b5ea16 |
375 | pTotal = TMath::Sqrt(pYZ *pYZ + pX * pX); |
376 | xEndAbsorber = fNonBendingCoor; |
377 | yEndAbsorber = fBendingCoor; |
378 | radiusEndAbsorber2 = xEndAbsorber * xEndAbsorber + yEndAbsorber * yEndAbsorber; |
379 | |
380 | if (radiusEndAbsorber2 > rLimit*rLimit) { |
04b5ea16 |
381 | zBP = zBP1; |
382 | } else { |
04b5ea16 |
383 | zBP = zBP2; |
384 | } |
385 | |
386 | xBP = xEndAbsorber - (pX / pZ) * (zEndAbsorber - zBP); |
387 | yBP = yEndAbsorber - (pY / pZ) * (zEndAbsorber - zBP); |
388 | |
389 | // new parameters after Branson and energy loss corrections |
b45fd22b |
390 | // Float_t zSmear = zBP - gRandom->Gaus(0.,2.); // !!! possible smearing of Z vertex position |
889a0215 |
391 | |
392 | Float_t zSmear = zBP ; |
b45fd22b |
393 | |
889a0215 |
394 | pZ = pTotal * (zSmear-zVtx) / TMath::Sqrt((xBP-xVtx) * (xBP-xVtx) + (yBP-yVtx) * (yBP-yVtx) +( zSmear-zVtx) * (zSmear-zVtx) ); |
395 | pX = pZ * (xBP - xVtx)/ (zSmear-zVtx); |
396 | pY = pZ * (yBP - yVtx) / (zSmear-zVtx); |
04b5ea16 |
397 | fBendingSlope = pY / pZ; |
398 | fNonBendingSlope = pX / pZ; |
5b64e914 |
399 | |
04b5ea16 |
400 | |
401 | pT = TMath::Sqrt(pX * pX + pY * pY); |
5b64e914 |
402 | theta = TMath::ATan2(pT, TMath::Abs(pZ)); |
b45fd22b |
403 | pTotal = TotalMomentumEnergyLoss(thetaLimit, pTotal, theta); |
04b5ea16 |
404 | |
405 | fInverseBendingMomentum = (sign / pTotal) * |
406 | TMath::Sqrt(1.0 + |
407 | fBendingSlope * fBendingSlope + |
408 | fNonBendingSlope * fNonBendingSlope) / |
409 | TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope); |
410 | |
411 | // vertex position at (0,0,0) |
412 | // should be taken from vertex measurement ??? |
889a0215 |
413 | |
414 | fBendingCoor = xVtx; |
415 | fNonBendingCoor = yVtx; |
416 | fZ= zVtx; |
417 | |
04b5ea16 |
418 | } |
b45fd22b |
419 | |
04b5ea16 |
420 | //__________________________________________________________________________ |
b45fd22b |
421 | Double_t AliMUONTrackParam::TotalMomentumEnergyLoss(Double_t thetaLimit, Double_t pTotal, Double_t theta) |
04b5ea16 |
422 | { |
423 | // Returns the total momentum corrected from energy loss in the front absorber |
43af2cb6 |
424 | // One can use the macros MUONTestAbso.C and DrawTestAbso.C |
425 | // to test this correction. |
b45fd22b |
426 | // Momentum energy loss behaviour evaluated with the simulation of single muons (april 2002) |
04b5ea16 |
427 | Double_t deltaP, pTotalCorrected; |
428 | |
b45fd22b |
429 | // Parametrization to be redone according to change of absorber material ???? |
956019b6 |
430 | // See remark in function BransonCorrection !!!! |
04b5ea16 |
431 | // The name is not so good, and there are many arguments !!!! |
b45fd22b |
432 | if (theta < thetaLimit ) { |
433 | if (pTotal < 20) { |
434 | deltaP = 2.5938 + 0.0570 * pTotal - 0.001151 * pTotal * pTotal; |
04b5ea16 |
435 | } else { |
b45fd22b |
436 | deltaP = 3.0714 + 0.011767 *pTotal; |
04b5ea16 |
437 | } |
438 | } else { |
b45fd22b |
439 | if (pTotal < 20) { |
440 | deltaP = 2.1207 + 0.05478 * pTotal - 0.00145079 * pTotal * pTotal; |
04b5ea16 |
441 | } else { |
b45fd22b |
442 | deltaP = 2.6069 + 0.0051705 * pTotal; |
04b5ea16 |
443 | } |
444 | } |
445 | pTotalCorrected = pTotal + deltaP / TMath::Cos(theta); |
446 | return pTotalCorrected; |
447 | } |
448 | |
b45fd22b |
449 | //__________________________________________________________________________ |
450 | void AliMUONTrackParam::FieldCorrection(Double_t Z) |
451 | { |
452 | // |
453 | // Correction of the effect of the magnetic field in the absorber |
454 | // Assume a constant field along Z axis. |
455 | |
456 | Float_t b[3],x[3]; |
457 | Double_t bZ; |
458 | Double_t pYZ,pX,pY,pZ,pT; |
459 | Double_t pXNew,pYNew; |
460 | Double_t c; |
461 | |
462 | pYZ = TMath::Abs(1.0 / fInverseBendingMomentum); |
463 | c = TMath::Sign(1.0,fInverseBendingMomentum); // particle charge |
464 | |
b8dc484b |
465 | pZ = Pz(); |
466 | pX = Px(); |
467 | pY = Py(); |
b45fd22b |
468 | pT = TMath::Sqrt(pX*pX+pY*pY); |
469 | |
5b64e914 |
470 | if (TMath::Abs(pZ) <= 0) return; |
b45fd22b |
471 | x[2] = Z/2; |
472 | x[0] = x[2]*fNonBendingSlope; |
473 | x[1] = x[2]*fBendingSlope; |
474 | |
475 | // Take magn. field value at position x. |
476 | gAlice->Field()->Field(x, b); |
477 | bZ = b[2]; |
478 | |
479 | // Transverse momentum rotation |
480 | // Parameterized with the study of DeltaPhi = phiReco - phiGen as a function of pZ. |
5b64e914 |
481 | Double_t phiShift = c*0.436*0.0003*bZ*Z/pZ; |
b45fd22b |
482 | // Rotate momentum around Z axis. |
483 | pXNew = pX*TMath::Cos(phiShift) - pY*TMath::Sin(phiShift); |
484 | pYNew = pX*TMath::Sin(phiShift) + pY*TMath::Cos(phiShift); |
485 | |
486 | fBendingSlope = pYNew / pZ; |
487 | fNonBendingSlope = pXNew / pZ; |
488 | |
489 | fInverseBendingMomentum = c / TMath::Sqrt(pYNew*pYNew+pZ*pZ); |
490 | |
b8dc484b |
491 | } |
492 | //__________________________________________________________________________ |
493 | Double_t AliMUONTrackParam::Px() |
494 | { |
495 | // return px from track paramaters |
496 | Double_t pYZ, pZ, pX; |
497 | pYZ = 0; |
498 | if ( TMath::Abs(fInverseBendingMomentum) > 0 ) |
499 | pYZ = TMath::Abs(1.0 / fInverseBendingMomentum); |
500 | pZ = -pYZ / (TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope)); // spectro. (z<0) |
501 | pX = pZ * fNonBendingSlope; |
502 | return pX; |
503 | } |
504 | //__________________________________________________________________________ |
505 | Double_t AliMUONTrackParam::Py() |
506 | { |
507 | // return px from track paramaters |
508 | Double_t pYZ, pZ, pY; |
509 | pYZ = 0; |
510 | if ( TMath::Abs(fInverseBendingMomentum) > 0 ) |
511 | pYZ = TMath::Abs(1.0 / fInverseBendingMomentum); |
512 | pZ = -pYZ / (TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope)); // spectro. (z<0) |
513 | pY = pZ * fBendingSlope; |
514 | return pY; |
515 | } |
516 | //__________________________________________________________________________ |
517 | Double_t AliMUONTrackParam::Pz() |
518 | { |
519 | // return px from track paramaters |
520 | Double_t pYZ, pZ; |
521 | pYZ = 0; |
522 | if ( TMath::Abs(fInverseBendingMomentum) > 0 ) |
523 | pYZ = TMath::Abs(1.0 / fInverseBendingMomentum); |
524 | pZ = -pYZ / (TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope)); // spectro. (z<0) |
525 | return pZ; |
526 | } |
527 | //__________________________________________________________________________ |
528 | Double_t AliMUONTrackParam::P() |
529 | { |
530 | // return p from track paramaters |
531 | Double_t pYZ, pZ, p; |
532 | pYZ = 0; |
533 | if ( TMath::Abs(fInverseBendingMomentum) > 0 ) |
534 | pYZ = TMath::Abs(1.0 / fInverseBendingMomentum); |
535 | pZ = -pYZ / (TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope)); // spectro. (z<0) |
536 | p = TMath::Abs(pZ) * |
537 | TMath::Sqrt(1.0 + fBendingSlope * fBendingSlope + fNonBendingSlope * fNonBendingSlope); |
538 | return p; |
539 | |
b45fd22b |
540 | } |
4d03a78e |
541 | //__________________________________________________________________________ |
542 | void AliMUONTrackParam::ExtrapOneStepHelix(Double_t charge, Double_t step, |
543 | Double_t *vect, Double_t *vout) |
544 | { |
545 | // ****************************************************************** |
546 | // * * |
547 | // * Performs the tracking of one step in a magnetic field * |
548 | // * The trajectory is assumed to be a helix in a constant field * |
549 | // * taken at the mid point of the step. * |
550 | // * Parameters: * |
551 | // * input * |
552 | // * STEP =arc length of the step asked * |
553 | // * VECT =input vector (position,direction cos and momentum) * |
554 | // * CHARGE= electric charge of the particle * |
555 | // * output * |
556 | // * VOUT = same as VECT after completion of the step * |
557 | // * * |
558 | // * ==>Called by : <USER>, GUSWIM * |
559 | // * Author m.hansroul ********* * |
560 | // * modified s.egli, s.v.levonian * |
561 | // * modified v.perevoztchikov |
562 | // * * |
563 | // ****************************************************************** |
564 | // |
565 | |
566 | // modif: everything in double precision |
567 | |
568 | Double_t xyz[3], h[4], hxp[3]; |
569 | Double_t h2xy, hp, rho, tet; |
570 | Double_t sint, sintt, tsint, cos1t; |
571 | Double_t f1, f2, f3, f4, f5, f6; |
572 | |
573 | const Int_t ix = 0; |
574 | const Int_t iy = 1; |
575 | const Int_t iz = 2; |
576 | const Int_t ipx = 3; |
577 | const Int_t ipy = 4; |
578 | const Int_t ipz = 5; |
579 | const Int_t ipp = 6; |
580 | |
581 | const Double_t ec = 2.9979251e-4; |
582 | // |
583 | // ------------------------------------------------------------------ |
584 | // |
585 | // units are kgauss,centimeters,gev/c |
586 | // |
587 | vout[ipp] = vect[ipp]; |
588 | if (TMath::Abs(charge) < 0.00001) { |
589 | for (Int_t i = 0; i < 3; i++) { |
590 | vout[i] = vect[i] + step * vect[i+3]; |
591 | vout[i+3] = vect[i+3]; |
592 | } |
593 | return; |
594 | } |
595 | xyz[0] = vect[ix] + 0.5 * step * vect[ipx]; |
596 | xyz[1] = vect[iy] + 0.5 * step * vect[ipy]; |
597 | xyz[2] = vect[iz] + 0.5 * step * vect[ipz]; |
598 | |
599 | //cmodif: call gufld (xyz, h) changed into: |
600 | GetField (xyz, h); |
601 | |
602 | h2xy = h[0]*h[0] + h[1]*h[1]; |
603 | h[3] = h[2]*h[2]+ h2xy; |
604 | if (h[3] < 1.e-12) { |
605 | for (Int_t i = 0; i < 3; i++) { |
606 | vout[i] = vect[i] + step * vect[i+3]; |
607 | vout[i+3] = vect[i+3]; |
608 | } |
609 | return; |
610 | } |
611 | if (h2xy < 1.e-12*h[3]) { |
612 | ExtrapOneStepHelix3(charge*h[2], step, vect, vout); |
613 | return; |
614 | } |
615 | h[3] = TMath::Sqrt(h[3]); |
616 | h[0] /= h[3]; |
617 | h[1] /= h[3]; |
618 | h[2] /= h[3]; |
619 | h[3] *= ec; |
620 | |
621 | hxp[0] = h[1]*vect[ipz] - h[2]*vect[ipy]; |
622 | hxp[1] = h[2]*vect[ipx] - h[0]*vect[ipz]; |
623 | hxp[2] = h[0]*vect[ipy] - h[1]*vect[ipx]; |
624 | |
625 | hp = h[0]*vect[ipx] + h[1]*vect[ipy] + h[2]*vect[ipz]; |
626 | |
627 | rho = -charge*h[3]/vect[ipp]; |
628 | tet = rho * step; |
629 | |
630 | if (TMath::Abs(tet) > 0.15) { |
631 | sint = TMath::Sin(tet); |
632 | sintt = (sint/tet); |
633 | tsint = (tet-sint)/tet; |
634 | cos1t = 2.*(TMath::Sin(0.5*tet))*(TMath::Sin(0.5*tet))/tet; |
635 | } else { |
636 | tsint = tet*tet/36.; |
637 | sintt = (1. - tsint); |
638 | sint = tet*sintt; |
639 | cos1t = 0.5*tet; |
640 | } |
641 | |
642 | f1 = step * sintt; |
643 | f2 = step * cos1t; |
644 | f3 = step * tsint * hp; |
645 | f4 = -tet*cos1t; |
646 | f5 = sint; |
647 | f6 = tet * cos1t * hp; |
648 | |
649 | vout[ix] = vect[ix] + f1*vect[ipx] + f2*hxp[0] + f3*h[0]; |
650 | vout[iy] = vect[iy] + f1*vect[ipy] + f2*hxp[1] + f3*h[1]; |
651 | vout[iz] = vect[iz] + f1*vect[ipz] + f2*hxp[2] + f3*h[2]; |
652 | |
653 | vout[ipx] = vect[ipx] + f4*vect[ipx] + f5*hxp[0] + f6*h[0]; |
654 | vout[ipy] = vect[ipy] + f4*vect[ipy] + f5*hxp[1] + f6*h[1]; |
655 | vout[ipz] = vect[ipz] + f4*vect[ipz] + f5*hxp[2] + f6*h[2]; |
656 | |
657 | return; |
658 | } |
659 | |
660 | //__________________________________________________________________________ |
661 | void AliMUONTrackParam::ExtrapOneStepHelix3(Double_t field, Double_t step, |
662 | Double_t *vect, Double_t *vout) |
663 | { |
664 | // |
665 | // ****************************************************************** |
666 | // * * |
667 | // * Tracking routine in a constant field oriented * |
668 | // * along axis 3 * |
669 | // * Tracking is performed with a conventional * |
670 | // * helix step method * |
671 | // * * |
672 | // * ==>Called by : <USER>, GUSWIM * |
673 | // * Authors R.Brun, M.Hansroul ********* * |
674 | // * Rewritten V.Perevoztchikov |
675 | // * * |
676 | // ****************************************************************** |
677 | // |
678 | |
679 | Double_t hxp[3]; |
680 | Double_t h4, hp, rho, tet; |
681 | Double_t sint, sintt, tsint, cos1t; |
682 | Double_t f1, f2, f3, f4, f5, f6; |
683 | |
684 | const Int_t ix = 0; |
685 | const Int_t iy = 1; |
686 | const Int_t iz = 2; |
687 | const Int_t ipx = 3; |
688 | const Int_t ipy = 4; |
689 | const Int_t ipz = 5; |
690 | const Int_t ipp = 6; |
691 | |
692 | const Double_t ec = 2.9979251e-4; |
693 | |
694 | // |
695 | // ------------------------------------------------------------------ |
696 | // |
697 | // units are kgauss,centimeters,gev/c |
698 | // |
699 | vout[ipp] = vect[ipp]; |
700 | h4 = field * ec; |
701 | |
702 | hxp[0] = - vect[ipy]; |
703 | hxp[1] = + vect[ipx]; |
704 | |
705 | hp = vect[ipz]; |
706 | |
707 | rho = -h4/vect[ipp]; |
708 | tet = rho * step; |
709 | if (TMath::Abs(tet) > 0.15) { |
710 | sint = TMath::Sin(tet); |
711 | sintt = (sint/tet); |
712 | tsint = (tet-sint)/tet; |
713 | cos1t = 2.* TMath::Sin(0.5*tet) * TMath::Sin(0.5*tet)/tet; |
714 | } else { |
715 | tsint = tet*tet/36.; |
716 | sintt = (1. - tsint); |
717 | sint = tet*sintt; |
718 | cos1t = 0.5*tet; |
719 | } |
720 | |
721 | f1 = step * sintt; |
722 | f2 = step * cos1t; |
723 | f3 = step * tsint * hp; |
724 | f4 = -tet*cos1t; |
725 | f5 = sint; |
726 | f6 = tet * cos1t * hp; |
727 | |
728 | vout[ix] = vect[ix] + f1*vect[ipx] + f2*hxp[0]; |
729 | vout[iy] = vect[iy] + f1*vect[ipy] + f2*hxp[1]; |
730 | vout[iz] = vect[iz] + f1*vect[ipz] + f3; |
731 | |
732 | vout[ipx] = vect[ipx] + f4*vect[ipx] + f5*hxp[0]; |
733 | vout[ipy] = vect[ipy] + f4*vect[ipy] + f5*hxp[1]; |
734 | vout[ipz] = vect[ipz] + f4*vect[ipz] + f6; |
735 | |
736 | return; |
737 | } |
738 | //__________________________________________________________________________ |
739 | void AliMUONTrackParam::ExtrapOneStepRungekutta(Double_t charge, Double_t step, |
740 | Double_t* vect, Double_t* vout) |
741 | { |
742 | // |
743 | // ****************************************************************** |
744 | // * * |
745 | // * Runge-Kutta method for tracking a particle through a magnetic * |
746 | // * field. Uses Nystroem algorithm (See Handbook Nat. Bur. of * |
747 | // * Standards, procedure 25.5.20) * |
748 | // * * |
749 | // * Input parameters * |
750 | // * CHARGE Particle charge * |
751 | // * STEP Step size * |
752 | // * VECT Initial co-ords,direction cosines,momentum * |
753 | // * Output parameters * |
754 | // * VOUT Output co-ords,direction cosines,momentum * |
755 | // * User routine called * |
756 | // * CALL GUFLD(X,F) * |
757 | // * * |
758 | // * ==>Called by : <USER>, GUSWIM * |
759 | // * Authors R.Brun, M.Hansroul ********* * |
760 | // * V.Perevoztchikov (CUT STEP implementation) * |
761 | // * * |
762 | // * * |
763 | // ****************************************************************** |
764 | // |
765 | |
766 | Double_t h2, h4, f[4]; |
767 | Double_t xyzt[3], a, b, c, ph,ph2; |
768 | Double_t secxs[4],secys[4],seczs[4],hxp[3]; |
769 | Double_t g1, g2, g3, g4, g5, g6, ang2, dxt, dyt, dzt; |
770 | Double_t est, at, bt, ct, cba; |
771 | Double_t f1, f2, f3, f4, rho, tet, hnorm, hp, rho1, sint, cost; |
772 | |
773 | Double_t x; |
774 | Double_t y; |
775 | Double_t z; |
776 | |
777 | Double_t xt; |
778 | Double_t yt; |
779 | Double_t zt; |
780 | |
781 | Double_t maxit = 1992; |
782 | Double_t maxcut = 11; |
783 | |
784 | const Double_t dlt = 1e-4; |
785 | const Double_t dlt32 = dlt/32.; |
786 | const Double_t third = 1./3.; |
787 | const Double_t half = 0.5; |
788 | const Double_t ec = 2.9979251e-4; |
789 | |
790 | const Double_t pisqua = 9.86960440109; |
791 | const Int_t ix = 0; |
792 | const Int_t iy = 1; |
793 | const Int_t iz = 2; |
794 | const Int_t ipx = 3; |
795 | const Int_t ipy = 4; |
796 | const Int_t ipz = 5; |
797 | |
798 | // *. |
799 | // *. ------------------------------------------------------------------ |
800 | // *. |
801 | // * this constant is for units cm,gev/c and kgauss |
802 | // * |
803 | Int_t iter = 0; |
804 | Int_t ncut = 0; |
805 | for(Int_t j = 0; j < 7; j++) |
806 | vout[j] = vect[j]; |
807 | |
808 | Double_t pinv = ec * charge / vect[6]; |
809 | Double_t tl = 0.; |
810 | Double_t h = step; |
811 | Double_t rest; |
812 | |
813 | |
814 | do { |
815 | rest = step - tl; |
816 | if (TMath::Abs(h) > TMath::Abs(rest)) h = rest; |
817 | //cmodif: call gufld(vout,f) changed into: |
818 | |
819 | GetField(vout,f); |
820 | |
821 | // * |
822 | // * start of integration |
823 | // * |
824 | x = vout[0]; |
825 | y = vout[1]; |
826 | z = vout[2]; |
827 | a = vout[3]; |
828 | b = vout[4]; |
829 | c = vout[5]; |
830 | |
831 | h2 = half * h; |
832 | h4 = half * h2; |
833 | ph = pinv * h; |
834 | ph2 = half * ph; |
835 | secxs[0] = (b * f[2] - c * f[1]) * ph2; |
836 | secys[0] = (c * f[0] - a * f[2]) * ph2; |
837 | seczs[0] = (a * f[1] - b * f[0]) * ph2; |
838 | ang2 = (secxs[0]*secxs[0] + secys[0]*secys[0] + seczs[0]*seczs[0]); |
839 | if (ang2 > pisqua) break; |
840 | |
841 | dxt = h2 * a + h4 * secxs[0]; |
842 | dyt = h2 * b + h4 * secys[0]; |
843 | dzt = h2 * c + h4 * seczs[0]; |
844 | xt = x + dxt; |
845 | yt = y + dyt; |
846 | zt = z + dzt; |
847 | // * |
848 | // * second intermediate point |
849 | // * |
850 | |
851 | est = TMath::Abs(dxt) + TMath::Abs(dyt) + TMath::Abs(dzt); |
852 | if (est > h) { |
853 | if (ncut++ > maxcut) break; |
854 | h *= half; |
855 | continue; |
856 | } |
857 | |
858 | xyzt[0] = xt; |
859 | xyzt[1] = yt; |
860 | xyzt[2] = zt; |
861 | |
862 | //cmodif: call gufld(xyzt,f) changed into: |
863 | GetField(xyzt,f); |
864 | |
865 | at = a + secxs[0]; |
866 | bt = b + secys[0]; |
867 | ct = c + seczs[0]; |
868 | |
869 | secxs[1] = (bt * f[2] - ct * f[1]) * ph2; |
870 | secys[1] = (ct * f[0] - at * f[2]) * ph2; |
871 | seczs[1] = (at * f[1] - bt * f[0]) * ph2; |
872 | at = a + secxs[1]; |
873 | bt = b + secys[1]; |
874 | ct = c + seczs[1]; |
875 | secxs[2] = (bt * f[2] - ct * f[1]) * ph2; |
876 | secys[2] = (ct * f[0] - at * f[2]) * ph2; |
877 | seczs[2] = (at * f[1] - bt * f[0]) * ph2; |
878 | dxt = h * (a + secxs[2]); |
879 | dyt = h * (b + secys[2]); |
880 | dzt = h * (c + seczs[2]); |
881 | xt = x + dxt; |
882 | yt = y + dyt; |
883 | zt = z + dzt; |
884 | at = a + 2.*secxs[2]; |
885 | bt = b + 2.*secys[2]; |
886 | ct = c + 2.*seczs[2]; |
887 | |
888 | est = TMath::Abs(dxt)+TMath::Abs(dyt)+TMath::Abs(dzt); |
889 | if (est > 2.*TMath::Abs(h)) { |
890 | if (ncut++ > maxcut) break; |
891 | h *= half; |
892 | continue; |
893 | } |
894 | |
895 | xyzt[0] = xt; |
896 | xyzt[1] = yt; |
897 | xyzt[2] = zt; |
898 | |
899 | //cmodif: call gufld(xyzt,f) changed into: |
900 | GetField(xyzt,f); |
901 | |
902 | z = z + (c + (seczs[0] + seczs[1] + seczs[2]) * third) * h; |
903 | y = y + (b + (secys[0] + secys[1] + secys[2]) * third) * h; |
904 | x = x + (a + (secxs[0] + secxs[1] + secxs[2]) * third) * h; |
905 | |
906 | secxs[3] = (bt*f[2] - ct*f[1])* ph2; |
907 | secys[3] = (ct*f[0] - at*f[2])* ph2; |
908 | seczs[3] = (at*f[1] - bt*f[0])* ph2; |
909 | a = a+(secxs[0]+secxs[3]+2. * (secxs[1]+secxs[2])) * third; |
910 | b = b+(secys[0]+secys[3]+2. * (secys[1]+secys[2])) * third; |
911 | c = c+(seczs[0]+seczs[3]+2. * (seczs[1]+seczs[2])) * third; |
912 | |
913 | est = TMath::Abs(secxs[0]+secxs[3] - (secxs[1]+secxs[2])) |
914 | + TMath::Abs(secys[0]+secys[3] - (secys[1]+secys[2])) |
915 | + TMath::Abs(seczs[0]+seczs[3] - (seczs[1]+seczs[2])); |
916 | |
917 | if (est > dlt && TMath::Abs(h) > 1.e-4) { |
918 | if (ncut++ > maxcut) break; |
919 | h *= half; |
920 | continue; |
921 | } |
922 | |
923 | ncut = 0; |
924 | // * if too many iterations, go to helix |
925 | if (iter++ > maxit) break; |
926 | |
927 | tl += h; |
928 | if (est < dlt32) |
929 | h *= 2.; |
930 | cba = 1./ TMath::Sqrt(a*a + b*b + c*c); |
931 | vout[0] = x; |
932 | vout[1] = y; |
933 | vout[2] = z; |
934 | vout[3] = cba*a; |
935 | vout[4] = cba*b; |
936 | vout[5] = cba*c; |
937 | rest = step - tl; |
938 | if (step < 0.) rest = -rest; |
939 | if (rest < 1.e-5*TMath::Abs(step)) return; |
940 | |
941 | } while(1); |
942 | |
943 | // angle too big, use helix |
944 | |
945 | f1 = f[0]; |
946 | f2 = f[1]; |
947 | f3 = f[2]; |
948 | f4 = TMath::Sqrt(f1*f1+f2*f2+f3*f3); |
949 | rho = -f4*pinv; |
950 | tet = rho * step; |
951 | |
952 | hnorm = 1./f4; |
953 | f1 = f1*hnorm; |
954 | f2 = f2*hnorm; |
955 | f3 = f3*hnorm; |
956 | |
957 | hxp[0] = f2*vect[ipz] - f3*vect[ipy]; |
958 | hxp[1] = f3*vect[ipx] - f1*vect[ipz]; |
959 | hxp[2] = f1*vect[ipy] - f2*vect[ipx]; |
960 | |
961 | hp = f1*vect[ipx] + f2*vect[ipy] + f3*vect[ipz]; |
962 | |
963 | rho1 = 1./rho; |
964 | sint = TMath::Sin(tet); |
965 | cost = 2.*TMath::Sin(half*tet)*TMath::Sin(half*tet); |
966 | |
967 | g1 = sint*rho1; |
968 | g2 = cost*rho1; |
969 | g3 = (tet-sint) * hp*rho1; |
970 | g4 = -cost; |
971 | g5 = sint; |
972 | g6 = cost * hp; |
973 | |
974 | vout[ix] = vect[ix] + g1*vect[ipx] + g2*hxp[0] + g3*f1; |
975 | vout[iy] = vect[iy] + g1*vect[ipy] + g2*hxp[1] + g3*f2; |
976 | vout[iz] = vect[iz] + g1*vect[ipz] + g2*hxp[2] + g3*f3; |
977 | |
978 | vout[ipx] = vect[ipx] + g4*vect[ipx] + g5*hxp[0] + g6*f1; |
979 | vout[ipy] = vect[ipy] + g4*vect[ipy] + g5*hxp[1] + g6*f2; |
980 | vout[ipz] = vect[ipz] + g4*vect[ipz] + g5*hxp[2] + g6*f3; |
981 | |
982 | return; |
983 | } |
984 | //___________________________________________________________ |
985 | void AliMUONTrackParam::GetField(Double_t *Position, Double_t *Field) |
986 | { |
987 | // interface to "gAlice->Field()->Field" for arguments in double precision |
988 | |
989 | Float_t x[3], b[3]; |
990 | |
991 | x[0] = Position[0]; x[1] = Position[1]; x[2] = Position[2]; |
992 | |
993 | gAlice->Field()->Field(x, b); |
994 | Field[0] = b[0]; Field[1] = b[1]; Field[2] = b[2]; |
995 | |
996 | return; |
997 | } |