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