bug fixed
[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 // Class AliMUONTrackExtrap
20 // ------------------------
21 // Tools for track extrapolation in ALICE dimuon spectrometer
22 // Author: Philippe Pillot
23 //-----------------------------------------------------------------------------
24
25 #include "AliMUONTrackExtrap.h" 
26 #include "AliMUONTrackParam.h"
27 #include "AliMUONConstants.h"
28 #include "AliMUONReconstructor.h"
29
30 #include "AliMagF.h"
31 #include "AliExternalTrackParam.h"
32
33 #include <TGeoGlobalMagField.h>
34 #include <TGeoManager.h>
35 #include <TMath.h>
36 #include <TDatabasePDG.h>
37
38 #include <Riostream.h>
39
40 /// \cond CLASSIMP
41 ClassImp(AliMUONTrackExtrap) // Class implementation in ROOT context
42 /// \endcond
43
44 const Double_t AliMUONTrackExtrap::fgkSimpleBPosition = 0.5 * (AliMUONConstants::CoilZ() + AliMUONConstants::YokeZ());
45 const Double_t AliMUONTrackExtrap::fgkSimpleBLength = 0.5 * (AliMUONConstants::CoilL() + AliMUONConstants::YokeL());
46       Double_t AliMUONTrackExtrap::fgSimpleBValue = 0.;
47       Bool_t   AliMUONTrackExtrap::fgFieldON = kFALSE;
48 const Bool_t   AliMUONTrackExtrap::fgkUseHelix = kFALSE;
49 const Int_t    AliMUONTrackExtrap::fgkMaxStepNumber = 5000;
50 const Double_t AliMUONTrackExtrap::fgkHelixStepLength = 6.;
51 const Double_t AliMUONTrackExtrap::fgkRungeKuttaMaxResidue = 0.002;
52
53 //__________________________________________________________________________
54 void AliMUONTrackExtrap::SetField()
55 {
56   /// set field on/off flag;  
57   /// set field at the centre of the dipole
58   const Double_t x[3] = {50.,50.,fgkSimpleBPosition};
59   Double_t b[3] = {0.,0.,0.};
60   TGeoGlobalMagField::Instance()->Field(x,b);
61   fgSimpleBValue = b[0];
62   fgFieldON = (TMath::Abs(fgSimpleBValue) > 1.e-10) ? kTRUE : kFALSE;
63   
64 }
65
66 //__________________________________________________________________________
67 Double_t AliMUONTrackExtrap::GetImpactParamFromBendingMomentum(Double_t bendingMomentum)
68 {
69   /// Returns impact parameter at vertex in bending plane (cm),
70   /// from the signed bending momentum "BendingMomentum" in bending plane (GeV/c),
71   /// using simple values for dipole magnetic field.
72   /// The sign of "BendingMomentum" is the sign of the charge.
73   
74   if (bendingMomentum == 0.) return 1.e10;
75   
76   const Double_t kCorrectionFactor = 1.1; // impact parameter is 10% underestimated
77   
78   return kCorrectionFactor * (-0.0003 * fgSimpleBValue * fgkSimpleBLength * fgkSimpleBPosition / bendingMomentum);
79 }
80
81 //__________________________________________________________________________
82 Double_t 
83 AliMUONTrackExtrap::GetBendingMomentumFromImpactParam(Double_t impactParam)
84 {
85   /// Returns signed bending momentum in bending plane (GeV/c),
86   /// the sign being the sign of the charge for particles moving forward in Z,
87   /// from the impact parameter "ImpactParam" at vertex in bending plane (cm),
88   /// using simple values for dipole magnetic field.
89   
90   if (impactParam == 0.) return 1.e10;
91   
92   const Double_t kCorrectionFactor = 1.1; // bending momentum is 10% underestimated
93   
94   if (fgFieldON) 
95   {
96     return kCorrectionFactor * (-0.0003 * fgSimpleBValue * fgkSimpleBLength * fgkSimpleBPosition / impactParam);
97   }
98   else 
99   {
100     return AliMUONConstants::GetMostProbBendingMomentum();
101   }
102 }
103
104 //__________________________________________________________________________
105 void AliMUONTrackExtrap::LinearExtrapToZ(AliMUONTrackParam* trackParam, Double_t zEnd)
106 {
107   /// Track parameters linearly extrapolated to the plane at "zEnd".
108   /// On return, results from the extrapolation are updated in trackParam.
109   
110   if (trackParam->GetZ() == zEnd) return; // nothing to be done if same z
111   
112   // Compute track parameters
113   Double_t dZ = zEnd - trackParam->GetZ();
114   trackParam->SetNonBendingCoor(trackParam->GetNonBendingCoor() + trackParam->GetNonBendingSlope() * dZ);
115   trackParam->SetBendingCoor(trackParam->GetBendingCoor() + trackParam->GetBendingSlope() * dZ);
116   trackParam->SetZ(zEnd);
117 }
118
119 //__________________________________________________________________________
120 void AliMUONTrackExtrap::LinearExtrapToZCov(AliMUONTrackParam* trackParam, Double_t zEnd, Bool_t updatePropagator)
121 {
122   /// Track parameters and their covariances linearly extrapolated to the plane at "zEnd".
123   /// On return, results from the extrapolation are updated in trackParam.
124   
125   if (trackParam->GetZ() == zEnd) return; // nothing to be done if same z
126   
127   // No need to propagate the covariance matrix if it does not exist
128   if (!trackParam->CovariancesExist()) {
129     cout<<"W-AliMUONTrackExtrap::LinearExtrapToZCov: Covariance matrix does not exist"<<endl;
130     // Extrapolate linearly track parameters to "zEnd"
131     LinearExtrapToZ(trackParam,zEnd);
132     return;
133   }
134   
135   // Compute track parameters
136   Double_t dZ = zEnd - trackParam->GetZ();
137   trackParam->SetNonBendingCoor(trackParam->GetNonBendingCoor() + trackParam->GetNonBendingSlope() * dZ);
138   trackParam->SetBendingCoor(trackParam->GetBendingCoor() + trackParam->GetBendingSlope() * dZ);
139   trackParam->SetZ(zEnd);
140   
141   // Calculate the jacobian related to the track parameters linear extrapolation to "zEnd"
142   TMatrixD jacob(5,5);
143   jacob.UnitMatrix();
144   jacob(0,1) = dZ;
145   jacob(2,3) = dZ;
146   
147   // Extrapolate track parameter covariances to "zEnd"
148   TMatrixD tmp(trackParam->GetCovariances(),TMatrixD::kMultTranspose,jacob);
149   TMatrixD tmp2(jacob,TMatrixD::kMult,tmp);
150   trackParam->SetCovariances(tmp2);
151   
152   // Update the propagator if required
153   if (updatePropagator) trackParam->UpdatePropagator(jacob);
154 }
155
156 //__________________________________________________________________________
157 Bool_t AliMUONTrackExtrap::ExtrapToZ(AliMUONTrackParam* trackParam, Double_t zEnd)
158 {
159   /// Interface to track parameter extrapolation to the plane at "Z" using Helix or Rungekutta algorithm.
160   /// On return, the track parameters resulting from the extrapolation are updated in trackParam.
161   if (!fgFieldON) {
162     AliMUONTrackExtrap::LinearExtrapToZ(trackParam,zEnd);
163     return kTRUE;
164   }
165   else if (fgkUseHelix) return AliMUONTrackExtrap::ExtrapToZHelix(trackParam,zEnd);
166   else return AliMUONTrackExtrap::ExtrapToZRungekutta(trackParam,zEnd);
167 }
168
169 //__________________________________________________________________________
170 Bool_t AliMUONTrackExtrap::ExtrapToZHelix(AliMUONTrackParam* trackParam, Double_t zEnd)
171 {
172   /// Track parameter extrapolation to the plane at "Z" using Helix algorithm.
173   /// On return, the track parameters resulting from the extrapolation are updated in trackParam.
174   if (trackParam->GetZ() == zEnd) return kTRUE; // nothing to be done if same Z
175   Double_t forwardBackward; // +1 if forward, -1 if backward
176   if (zEnd < trackParam->GetZ()) forwardBackward = 1.0; // spectro. z<0 
177   else forwardBackward = -1.0;
178   Double_t v3[7], v3New[7]; // 7 in parameter ????
179   Int_t i3, stepNumber;
180   // For safety: return kTRUE or kFALSE ????
181   // Parameter vector for calling EXTRAP_ONESTEP
182   ConvertTrackParamForExtrap(trackParam, forwardBackward, v3);
183   // sign of charge (sign of fInverseBendingMomentum if forward motion)
184   // must be changed if backward extrapolation
185   Double_t chargeExtrap = forwardBackward * TMath::Sign(Double_t(1.0), trackParam->GetInverseBendingMomentum());
186   // Extrapolation loop
187   stepNumber = 0;
188   while (((-forwardBackward * (v3[2] - zEnd)) <= 0.0) && (stepNumber < fgkMaxStepNumber)) { // spectro. z<0
189     stepNumber++;
190     ExtrapOneStepHelix(chargeExtrap, fgkHelixStepLength, v3, v3New);
191     if ((-forwardBackward * (v3New[2] - zEnd)) > 0.0) break; // one is beyond Z spectro. z<0
192                                                              // better use TArray ????
193     for (i3 = 0; i3 < 7; i3++) {v3[i3] = v3New[i3];}
194   }
195   // check fgkMaxStepNumber ????
196   // Interpolation back to exact Z (2nd order)
197   // should be in function ???? using TArray ????
198   Double_t dZ12 = v3New[2] - v3[2]; // 1->2
199   if (TMath::Abs(dZ12) > 0) {
200     Double_t dZ1i = zEnd - v3[2]; // 1-i
201     Double_t dZi2 = v3New[2] - zEnd; // i->2
202     Double_t xPrime = (v3New[0] - v3[0]) / dZ12;
203     Double_t xSecond = ((v3New[3] / v3New[5]) - (v3[3] / v3[5])) / dZ12;
204     Double_t yPrime = (v3New[1] - v3[1]) / dZ12;
205     Double_t ySecond = ((v3New[4] / v3New[5]) - (v3[4] / v3[5])) / dZ12;
206     v3[0] = v3[0] + xPrime * dZ1i - 0.5 * xSecond * dZ1i * dZi2; // X
207     v3[1] = v3[1] + yPrime * dZ1i - 0.5 * ySecond * dZ1i * dZi2; // Y
208     v3[2] = zEnd; // Z
209     Double_t xPrimeI = xPrime - 0.5 * xSecond * (dZi2 - dZ1i);
210     Double_t yPrimeI = yPrime - 0.5 * ySecond * (dZi2 - dZ1i);
211     // (PX, PY, PZ)/PTOT assuming forward motion
212     v3[5] = 1.0 / TMath::Sqrt(1.0 + xPrimeI * xPrimeI + yPrimeI * yPrimeI); // PZ/PTOT
213     v3[3] = xPrimeI * v3[5]; // PX/PTOT
214     v3[4] = yPrimeI * v3[5]; // PY/PTOT
215   } else {
216     cout<<"W-AliMUONTrackExtrap::ExtrapToZHelix: Extrap. to Z not reached, Z = "<<zEnd<<endl;
217   }
218   // Recover track parameters (charge back for forward motion)
219   RecoverTrackParam(v3, chargeExtrap * forwardBackward, trackParam);
220   return kTRUE;
221 }
222
223 //__________________________________________________________________________
224 Bool_t AliMUONTrackExtrap::ExtrapToZRungekutta(AliMUONTrackParam* trackParam, Double_t zEnd)
225 {
226   /// Track parameter extrapolation to the plane at "Z" using Rungekutta algorithm.
227   /// On return, the track parameters resulting from the extrapolation are updated in trackParam.
228   if (trackParam->GetZ() == zEnd) return kTRUE; // nothing to be done if same Z
229   Double_t forwardBackward; // +1 if forward, -1 if backward
230   if (zEnd < trackParam->GetZ()) forwardBackward = 1.0; // spectro. z<0 
231   else forwardBackward = -1.0;
232   // sign of charge (sign of fInverseBendingMomentum if forward motion)
233   // must be changed if backward extrapolation
234   Double_t chargeExtrap = forwardBackward * TMath::Sign(Double_t(1.0), trackParam->GetInverseBendingMomentum());
235   Double_t v3[7], v3New[7];
236   Double_t dZ, step;
237   Int_t stepNumber = 0;
238   
239   // Extrapolation loop (until within tolerance or the track turn around)
240   Double_t residue = zEnd - trackParam->GetZ();
241   Bool_t uturn = kFALSE;
242   Bool_t trackingFailed = kFALSE;
243   Bool_t tooManyStep = kFALSE;
244   while (TMath::Abs(residue) > fgkRungeKuttaMaxResidue && stepNumber <= fgkMaxStepNumber) {
245     
246     dZ = zEnd - trackParam->GetZ();
247     // step lenght assuming linear trajectory
248     step = dZ * TMath::Sqrt(1.0 + trackParam->GetBendingSlope()*trackParam->GetBendingSlope() +
249                             trackParam->GetNonBendingSlope()*trackParam->GetNonBendingSlope());
250     ConvertTrackParamForExtrap(trackParam, forwardBackward, v3);
251     
252     do { // reduce step lenght while zEnd oversteped
253       if (stepNumber > fgkMaxStepNumber) {
254         cout<<"W-AliMUONTrackExtrap::ExtrapToZRungekutta: Too many trials: "<<stepNumber<<endl;
255         tooManyStep = kTRUE;
256         break;
257       }
258       stepNumber ++;
259       step = TMath::Abs(step);
260       if (!AliMUONTrackExtrap::ExtrapOneStepRungekutta(chargeExtrap,step,v3,v3New)) {
261         trackingFailed = kTRUE;
262         break;
263       }
264       residue = zEnd - v3New[2];
265       step *= dZ/(v3New[2]-trackParam->GetZ());
266     } while (residue*dZ < 0 && TMath::Abs(residue) > fgkRungeKuttaMaxResidue);
267     
268     if (trackingFailed) break;
269     else if (v3New[5]*v3[5] < 0) { // the track turned around
270       cout<<"W-AliMUONTrackExtrap::ExtrapToZRungekutta: The track turned around"<<endl;
271       uturn = kTRUE;
272       break;
273     } else RecoverTrackParam(v3New, chargeExtrap * forwardBackward, trackParam);
274     
275   }
276   
277   // terminate the extropolation with a straight line up to the exact "zEnd" value
278   if (trackingFailed || uturn) {
279     
280     // track ends +-100 meters away in the bending direction
281     dZ = zEnd - v3[2];
282     Double_t bendingSlope = TMath::Sign(1.e4,-fgSimpleBValue*trackParam->GetInverseBendingMomentum()) / dZ;
283     Double_t pZ = TMath::Abs(1. / trackParam->GetInverseBendingMomentum()) / TMath::Sqrt(1.0 + bendingSlope * bendingSlope);
284     Double_t nonBendingSlope = TMath::Sign(TMath::Abs(v3[3]) * v3[6] / pZ, trackParam->GetNonBendingSlope());
285     trackParam->SetNonBendingCoor(trackParam->GetNonBendingCoor() + dZ * nonBendingSlope);
286     trackParam->SetNonBendingSlope(nonBendingSlope);
287     trackParam->SetBendingCoor(trackParam->GetBendingCoor() + dZ * bendingSlope);
288     trackParam->SetBendingSlope(bendingSlope);
289     trackParam->SetZ(zEnd);
290     
291     return kFALSE;
292     
293   } else {
294     
295     // track extrapolated normally
296     trackParam->SetNonBendingCoor(trackParam->GetNonBendingCoor() + residue * trackParam->GetNonBendingSlope());
297     trackParam->SetBendingCoor(trackParam->GetBendingCoor() + residue * trackParam->GetBendingSlope());
298     trackParam->SetZ(zEnd);
299     
300     return !tooManyStep;
301     
302   }
303   
304 }
305
306 //__________________________________________________________________________
307 void AliMUONTrackExtrap::ConvertTrackParamForExtrap(AliMUONTrackParam* trackParam, Double_t forwardBackward, Double_t *v3)
308 {
309   /// Set vector of Geant3 parameters pointed to by "v3" from track parameters in trackParam.
310   /// Since AliMUONTrackParam is only geometry, one uses "forwardBackward"
311   /// to know whether the particle is going forward (+1) or backward (-1).
312   v3[0] = trackParam->GetNonBendingCoor(); // X
313   v3[1] = trackParam->GetBendingCoor(); // Y
314   v3[2] = trackParam->GetZ(); // Z
315   Double_t pYZ = TMath::Abs(1.0 / trackParam->GetInverseBendingMomentum());
316   Double_t pZ = pYZ / TMath::Sqrt(1.0 + trackParam->GetBendingSlope() * trackParam->GetBendingSlope());
317   v3[6] = TMath::Sqrt(pYZ * pYZ + pZ * pZ * trackParam->GetNonBendingSlope() * trackParam->GetNonBendingSlope()); // PTOT
318   v3[5] = -forwardBackward * pZ / v3[6]; // PZ/PTOT spectro. z<0
319   v3[3] = trackParam->GetNonBendingSlope() * v3[5]; // PX/PTOT
320   v3[4] = trackParam->GetBendingSlope() * v3[5]; // PY/PTOT
321 }
322
323 //__________________________________________________________________________
324 void AliMUONTrackExtrap::RecoverTrackParam(Double_t *v3, Double_t charge, AliMUONTrackParam* trackParam)
325 {
326   /// Set track parameters in trackParam from Geant3 parameters pointed to by "v3",
327   /// assumed to be calculated for forward motion in Z.
328   /// "InverseBendingMomentum" is signed with "charge".
329   trackParam->SetNonBendingCoor(v3[0]); // X
330   trackParam->SetBendingCoor(v3[1]); // Y
331   trackParam->SetZ(v3[2]); // Z
332   Double_t pYZ = v3[6] * TMath::Sqrt((1.-v3[3])*(1.+v3[3]));
333   trackParam->SetInverseBendingMomentum(charge/pYZ);
334   trackParam->SetBendingSlope(v3[4]/v3[5]);
335   trackParam->SetNonBendingSlope(v3[3]/v3[5]);
336 }
337
338 //__________________________________________________________________________
339 Bool_t AliMUONTrackExtrap::ExtrapToZCov(AliMUONTrackParam* trackParam, Double_t zEnd, Bool_t updatePropagator)
340 {
341   /// Track parameters and their covariances extrapolated to the plane at "zEnd".
342   /// On return, results from the extrapolation are updated in trackParam.
343   
344   if (trackParam->GetZ() == zEnd) return kTRUE; // nothing to be done if same z
345   
346   if (!fgFieldON) { // linear extrapolation if no magnetic field
347     AliMUONTrackExtrap::LinearExtrapToZCov(trackParam,zEnd,updatePropagator);
348     return kTRUE;
349   }
350   
351   // No need to propagate the covariance matrix if it does not exist
352   if (!trackParam->CovariancesExist()) {
353     cout<<"W-AliMUONTrackExtrap::ExtrapToZCov: Covariance matrix does not exist"<<endl;
354     // Extrapolate track parameters to "zEnd"
355     return ExtrapToZ(trackParam,zEnd);
356   }
357   
358   // Save the actual track parameters
359   AliMUONTrackParam trackParamSave(*trackParam);
360   TMatrixD paramSave(trackParamSave.GetParameters());
361   Double_t zBegin = trackParamSave.GetZ();
362   
363   // Get reference to the parameter covariance matrix
364   const TMatrixD& kParamCov = trackParam->GetCovariances();
365         
366   // Extrapolate track parameters to "zEnd"
367   // Do not update the covariance matrix if the extrapolation failed
368   if (!ExtrapToZ(trackParam,zEnd)) return kFALSE;
369   
370   // Get reference to the extrapolated parameters
371   const TMatrixD& extrapParam = trackParam->GetParameters();
372   
373   // Calculate the jacobian related to the track parameters extrapolation to "zEnd"
374   Bool_t extrapStatus = kTRUE;
375   TMatrixD jacob(5,5);
376   jacob.Zero();
377   TMatrixD dParam(5,1);
378   Double_t direction[5] = {-1.,-1.,1.,1.,-1.};
379   for (Int_t i=0; i<5; i++) {
380     // Skip jacobian calculation for parameters with no associated error
381     if (kParamCov(i,i) <= 0.) continue;
382     
383     // Small variation of parameter i only
384     for (Int_t j=0; j<5; j++) {
385       if (j==i) {
386         dParam(j,0) = TMath::Sqrt(kParamCov(i,i));
387         dParam(j,0) *= TMath::Sign(1.,direction[j]*paramSave(j,0)); // variation always in the same direction
388       } else dParam(j,0) = 0.;
389     }
390     
391     // Set new parameters
392     trackParamSave.SetParameters(paramSave);
393     trackParamSave.AddParameters(dParam);
394     trackParamSave.SetZ(zBegin);
395     
396     // Extrapolate new track parameters to "zEnd"
397     if (!ExtrapToZ(&trackParamSave,zEnd)) {
398       cout<<"W-AliMUONTrackExtrap::ExtrapToZCov: Bad covariance matrix"<<endl;
399       extrapStatus = kFALSE;
400     }
401     
402     // Calculate the jacobian
403     TMatrixD jacobji(trackParamSave.GetParameters(),TMatrixD::kMinus,extrapParam);
404     jacobji *= 1. / dParam(i,0);
405     jacob.SetSub(0,i,jacobji);
406   }
407   
408   // Extrapolate track parameter covariances to "zEnd"
409   TMatrixD tmp(kParamCov,TMatrixD::kMultTranspose,jacob);
410   TMatrixD tmp2(jacob,TMatrixD::kMult,tmp);
411   trackParam->SetCovariances(tmp2);
412   
413   // Update the propagator if required
414   if (updatePropagator) trackParam->UpdatePropagator(jacob);
415   
416   return extrapStatus;
417 }
418
419 //__________________________________________________________________________
420 void AliMUONTrackExtrap::AddMCSEffectInAbsorber(AliMUONTrackParam* param, Double_t signedPathLength, Double_t f0, Double_t f1, Double_t f2)
421 {
422   /// Add to the track parameter covariances the effects of multiple Coulomb scattering
423   /// signedPathLength must have the sign of (zOut - zIn) where all other parameters are assumed to be given at zOut.
424   
425   // absorber related covariance parameters
426   Double_t bendingSlope = param->GetBendingSlope();
427   Double_t nonBendingSlope = param->GetNonBendingSlope();
428   Double_t inverseBendingMomentum = param->GetInverseBendingMomentum();
429   Double_t alpha2 = 0.0136 * 0.0136 * inverseBendingMomentum * inverseBendingMomentum * (1.0 + bendingSlope * bendingSlope) /
430                     (1.0 + bendingSlope *bendingSlope + nonBendingSlope * nonBendingSlope); // velocity = 1
431   Double_t pathLength = TMath::Abs(signedPathLength);
432   Double_t varCoor = alpha2 * (pathLength * pathLength * f0 - 2. * pathLength * f1 + f2);
433   Double_t covCorrSlope = TMath::Sign(1.,signedPathLength) * alpha2 * (pathLength * f0 - f1);
434   Double_t varSlop = alpha2 * f0;
435   
436   // Set MCS covariance matrix
437   TMatrixD newParamCov(param->GetCovariances());
438   // Non bending plane
439   newParamCov(0,0) += varCoor;       newParamCov(0,1) += covCorrSlope;
440   newParamCov(1,0) += covCorrSlope;  newParamCov(1,1) += varSlop;
441   // Bending plane
442   newParamCov(2,2) += varCoor;       newParamCov(2,3) += covCorrSlope;
443   newParamCov(3,2) += covCorrSlope;  newParamCov(3,3) += varSlop;
444   
445   // Set momentum related covariances if B!=0
446   if (fgFieldON) {
447     // compute derivative d(q/Pxy) / dSlopeX and d(q/Pxy) / dSlopeY
448     Double_t dqPxydSlopeX = inverseBendingMomentum * nonBendingSlope / (1. + nonBendingSlope*nonBendingSlope + bendingSlope*bendingSlope);
449     Double_t dqPxydSlopeY = - inverseBendingMomentum * nonBendingSlope*nonBendingSlope * bendingSlope /
450                               (1. + bendingSlope*bendingSlope) / (1. + nonBendingSlope*nonBendingSlope + bendingSlope*bendingSlope);
451     // Inverse bending momentum (due to dependences with bending and non bending slopes)
452     newParamCov(4,0) += dqPxydSlopeX * covCorrSlope; newParamCov(0,4) += dqPxydSlopeX * covCorrSlope;
453     newParamCov(4,1) += dqPxydSlopeX * varSlop;      newParamCov(1,4) += dqPxydSlopeX * varSlop;
454     newParamCov(4,2) += dqPxydSlopeY * covCorrSlope; newParamCov(2,4) += dqPxydSlopeY * covCorrSlope;
455     newParamCov(4,3) += dqPxydSlopeY * varSlop;      newParamCov(3,4) += dqPxydSlopeY * varSlop;
456     newParamCov(4,4) += (dqPxydSlopeX*dqPxydSlopeX + dqPxydSlopeY*dqPxydSlopeY) * varSlop;
457   }
458   
459   // Set new covariances
460   param->SetCovariances(newParamCov);
461 }
462
463 //__________________________________________________________________________
464 void AliMUONTrackExtrap::CorrectMCSEffectInAbsorber(AliMUONTrackParam* param,
465                                                     Double_t xVtx, Double_t yVtx, Double_t zVtx,
466                                                     Double_t errXVtx, Double_t errYVtx,
467                                                     Double_t absZBeg, Double_t pathLength, Double_t f0, Double_t f1, Double_t f2)
468 {
469   /// Correct parameters and corresponding covariances using Branson correction
470   /// - input param are parameters and covariances at the end of absorber
471   /// - output param are parameters and covariances at vertex
472   /// Absorber correction parameters are supposed to be calculated at the current track z-position
473   
474   // Position of the Branson plane (spectro. (z<0))
475   Double_t zB = (f1>0.) ? absZBeg - f2/f1 : 0.;
476   
477   // Add MCS effects to current parameter covariances (spectro. (z<0))
478   AddMCSEffectInAbsorber(param, -pathLength, f0, f1, f2);
479   
480   // Get track parameters and covariances in the Branson plane corrected for magnetic field effect
481   ExtrapToZCov(param,zVtx);
482   LinearExtrapToZCov(param,zB);
483   
484   // compute track parameters at vertex
485   TMatrixD newParam(5,1);
486   newParam(0,0) = xVtx;
487   newParam(1,0) = (param->GetNonBendingCoor() - xVtx) / (zB - zVtx);
488   newParam(2,0) = yVtx;
489   newParam(3,0) = (param->GetBendingCoor() - yVtx) / (zB - zVtx);
490   newParam(4,0) = param->GetCharge() / param->P() *
491                   TMath::Sqrt(1.0 + newParam(1,0)*newParam(1,0) + newParam(3,0)*newParam(3,0)) /
492                   TMath::Sqrt(1.0 + newParam(3,0)*newParam(3,0));
493   
494   // Get covariances in (X, SlopeX, Y, SlopeY, q*PTot) coordinate system
495   TMatrixD paramCovP(param->GetCovariances());
496   Cov2CovP(param->GetParameters(),paramCovP);
497   
498   // Get the covariance matrix in the (XVtx, X, YVtx, Y, q*PTot) coordinate system
499   TMatrixD paramCovVtx(5,5);
500   paramCovVtx.Zero();
501   paramCovVtx(0,0) = errXVtx * errXVtx;
502   paramCovVtx(1,1) = paramCovP(0,0);
503   paramCovVtx(2,2) = errYVtx * errYVtx;
504   paramCovVtx(3,3) = paramCovP(2,2);
505   paramCovVtx(4,4) = paramCovP(4,4);
506   paramCovVtx(1,3) = paramCovP(0,2);
507   paramCovVtx(3,1) = paramCovP(2,0);
508   paramCovVtx(1,4) = paramCovP(0,4);
509   paramCovVtx(4,1) = paramCovP(4,0);
510   paramCovVtx(3,4) = paramCovP(2,4);
511   paramCovVtx(4,3) = paramCovP(4,2);
512   
513   // Jacobian of the transformation (XVtx, X, YVtx, Y, q*PTot) -> (XVtx, SlopeXVtx, YVtx, SlopeYVtx, q*PTotVtx)
514   TMatrixD jacob(5,5);
515   jacob.UnitMatrix();
516   jacob(1,0) = - 1. / (zB - zVtx);
517   jacob(1,1) = 1. / (zB - zVtx);
518   jacob(3,2) = - 1. / (zB - zVtx);
519   jacob(3,3) = 1. / (zB - zVtx);
520   
521   // Compute covariances at vertex in the (XVtx, SlopeXVtx, YVtx, SlopeYVtx, q*PTotVtx) coordinate system
522   TMatrixD tmp(paramCovVtx,TMatrixD::kMultTranspose,jacob);
523   TMatrixD newParamCov(jacob,TMatrixD::kMult,tmp);
524   
525   // Compute covariances at vertex in the (XVtx, SlopeXVtx, YVtx, SlopeYVtx, q/PyzVtx) coordinate system
526   CovP2Cov(newParam,newParamCov);
527   
528   // Set parameters and covariances at vertex
529   param->SetParameters(newParam);
530   param->SetZ(zVtx);
531   param->SetCovariances(newParamCov);
532 }
533
534 //__________________________________________________________________________
535 void AliMUONTrackExtrap::CorrectELossEffectInAbsorber(AliMUONTrackParam* param, Double_t eLoss, Double_t sigmaELoss2)
536 {
537   /// Correct parameters for energy loss and add energy loss fluctuation effect to covariances
538   
539   // Get parameter covariances in (X, SlopeX, Y, SlopeY, q*PTot) coordinate system
540   TMatrixD newParamCov(param->GetCovariances());
541   Cov2CovP(param->GetParameters(),newParamCov);
542   
543   // Compute new parameters corrected for energy loss
544   Double_t muMass = TDatabasePDG::Instance()->GetParticle("mu-")->Mass(); // GeV
545   Double_t p = param->P();
546   Double_t e = TMath::Sqrt(p*p + muMass*muMass);
547   Double_t eCorr = e + eLoss;
548   Double_t pCorr = TMath::Sqrt(eCorr*eCorr - muMass*muMass);
549   Double_t nonBendingSlope = param->GetNonBendingSlope();
550   Double_t bendingSlope = param->GetBendingSlope();
551   param->SetInverseBendingMomentum(param->GetCharge() / pCorr *
552                                    TMath::Sqrt(1.0 + nonBendingSlope*nonBendingSlope + bendingSlope*bendingSlope) /
553                                    TMath::Sqrt(1.0 + bendingSlope*bendingSlope));
554   
555   // Add effects of energy loss fluctuation to covariances
556   newParamCov(4,4) += eCorr * eCorr / pCorr / pCorr * sigmaELoss2;
557   
558   // Get new parameter covariances in (X, SlopeX, Y, SlopeY, q/Pyz) coordinate system
559   CovP2Cov(param->GetParameters(),newParamCov);
560   
561   // Set new parameter covariances
562   param->SetCovariances(newParamCov);
563 }
564
565 //__________________________________________________________________________
566 Bool_t AliMUONTrackExtrap::GetAbsorberCorrectionParam(Double_t trackXYZIn[3], Double_t trackXYZOut[3], Double_t pTotal,
567                                                       Double_t &pathLength, Double_t &f0, Double_t &f1, Double_t &f2,
568                                                       Double_t &meanRho, Double_t &totalELoss, Double_t &sigmaELoss2)
569 {
570   /// Parameters used to correct for Multiple Coulomb Scattering and energy loss in absorber
571   /// Calculated assuming a linear propagation from trackXYZIn to trackXYZOut (order is important)
572   // pathLength: path length between trackXYZIn and trackXYZOut (cm)
573   // f0:         0th moment of z calculated with the inverse radiation-length distribution
574   // f1:         1st moment of z calculated with the inverse radiation-length distribution
575   // f2:         2nd moment of z calculated with the inverse radiation-length distribution
576   // meanRho:    average density of crossed material (g/cm3)
577   // totalELoss: total energy loss in absorber
578   
579   // Reset absorber's parameters
580   pathLength = 0.;
581   f0 = 0.;
582   f1 = 0.;
583   f2 = 0.;
584   meanRho = 0.;
585   totalELoss = 0.;
586   sigmaELoss2 = 0.;
587   
588   // Check whether the geometry is available
589   if (!gGeoManager) {
590     cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: no TGeo"<<endl;
591     return kFALSE;
592   }
593   
594   // Initialize starting point and direction
595   pathLength = TMath::Sqrt((trackXYZOut[0] - trackXYZIn[0])*(trackXYZOut[0] - trackXYZIn[0])+
596                            (trackXYZOut[1] - trackXYZIn[1])*(trackXYZOut[1] - trackXYZIn[1])+
597                            (trackXYZOut[2] - trackXYZIn[2])*(trackXYZOut[2] - trackXYZIn[2]));
598   if (pathLength < TGeoShape::Tolerance()) return kFALSE;
599   Double_t b[3];
600   b[0] = (trackXYZOut[0] - trackXYZIn[0]) / pathLength;
601   b[1] = (trackXYZOut[1] - trackXYZIn[1]) / pathLength;
602   b[2] = (trackXYZOut[2] - trackXYZIn[2]) / pathLength;
603   TGeoNode *currentnode = gGeoManager->InitTrack(trackXYZIn, b);
604   if (!currentnode) {
605     cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: start point out of geometry"<<endl;
606     return kFALSE;
607   }
608   
609   // loop over absorber slices and calculate absorber's parameters
610   Double_t rho = 0.; // material density (g/cm3)
611   Double_t x0 = 0.;  // radiation-length (cm-1)
612   Double_t atomicA = 0.; // A of material
613   Double_t atomicZ = 0.; // Z of material
614   Double_t atomicZoverA = 0.; // Z/A of material
615   Double_t localPathLength = 0;
616   Double_t remainingPathLength = pathLength;
617   Double_t zB = trackXYZIn[2];
618   Double_t zE, dzB, dzE;
619   do {
620     // Get material properties
621     TGeoMaterial *material = currentnode->GetVolume()->GetMedium()->GetMaterial();
622     rho = material->GetDensity();
623     x0 = material->GetRadLen();
624     atomicA = material->GetA();
625     atomicZ = material->GetZ();
626     if(material->IsMixture()){
627       TGeoMixture * mixture = (TGeoMixture*)material;
628       atomicZoverA = 0.;
629       Double_t sum = 0.;
630       for (Int_t iel=0;iel<mixture->GetNelements();iel++){
631         sum  += mixture->GetWmixt()[iel];
632         atomicZoverA += mixture->GetZmixt()[iel]*mixture->GetWmixt()[iel]/mixture->GetAmixt()[iel];
633       }
634       atomicZoverA/=sum;
635     }
636     else atomicZoverA = atomicZ/atomicA;
637     
638     // Get path length within this material
639     gGeoManager->FindNextBoundary(remainingPathLength);
640     localPathLength = gGeoManager->GetStep() + 1.e-6;
641     // Check if boundary within remaining path length. If so, make sure to cross the boundary to prepare the next step
642     if (localPathLength >= remainingPathLength) localPathLength = remainingPathLength;
643     else {
644       currentnode = gGeoManager->Step();
645       if (!currentnode) {
646         cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: navigation failed"<<endl;
647         f0 = f1 = f2 = meanRho = totalELoss = sigmaELoss2 = 0.;
648         return kFALSE;
649       }
650       if (!gGeoManager->IsEntering()) {
651         // make another small step to try to enter in new absorber slice
652         gGeoManager->SetStep(0.001);
653         currentnode = gGeoManager->Step();
654         if (!gGeoManager->IsEntering() || !currentnode) {
655           cout<<"E-AliMUONTrackExtrap::GetAbsorberCorrectionParam: navigation failed"<<endl;
656           f0 = f1 = f2 = meanRho = totalELoss = sigmaELoss2 = 0.;
657           return kFALSE;
658         }
659         localPathLength += 0.001;
660       }
661     }
662     
663     // calculate absorber's parameters
664     zE = b[2] * localPathLength + zB;
665     dzB = zB - trackXYZIn[2];
666     dzE = zE - trackXYZIn[2];
667     f0 += localPathLength / x0;
668     f1 += (dzE*dzE - dzB*dzB) / b[2] / b[2] / x0 / 2.;
669     f2 += (dzE*dzE*dzE - dzB*dzB*dzB) / b[2] / b[2] / b[2] / x0 / 3.;
670     meanRho += localPathLength * rho;
671     totalELoss += BetheBloch(pTotal, localPathLength, rho, atomicZ, atomicZoverA);
672     sigmaELoss2 += EnergyLossFluctuation2(pTotal, localPathLength, rho, atomicZoverA);
673     
674     // prepare next step
675     zB = zE;
676     remainingPathLength -= localPathLength;
677   } while (remainingPathLength > TGeoShape::Tolerance());
678   
679   meanRho /= pathLength;
680   
681   return kTRUE;
682 }
683
684 //__________________________________________________________________________
685 Double_t AliMUONTrackExtrap::GetMCSAngle2(const AliMUONTrackParam& param, Double_t dZ, Double_t x0)
686 {
687   /// Return the angular dispersion square due to multiple Coulomb scattering
688   /// through a material of thickness "dZ" and of radiation length "x0"
689   /// assuming linear propagation and using the small angle approximation.
690   
691   Double_t bendingSlope = param.GetBendingSlope();
692   Double_t nonBendingSlope = param.GetNonBendingSlope();
693   Double_t inverseTotalMomentum2 = param.GetInverseBendingMomentum() * param.GetInverseBendingMomentum() *
694                                    (1.0 + bendingSlope * bendingSlope) /
695                                    (1.0 + bendingSlope *bendingSlope + nonBendingSlope * nonBendingSlope); 
696   // Path length in the material
697   Double_t pathLength = TMath::Abs(dZ) * TMath::Sqrt(1.0 + bendingSlope*bendingSlope + nonBendingSlope*nonBendingSlope);
698   // relativistic velocity
699   Double_t velo = 1.;
700   // Angular dispersion square of the track (variance) in a plane perpendicular to the trajectory
701   Double_t theta02 = 0.0136 / velo * (1 + 0.038 * TMath::Log(pathLength/x0));
702   
703   return theta02 * theta02 * inverseTotalMomentum2 * pathLength / x0;
704 }
705
706 //__________________________________________________________________________
707 void AliMUONTrackExtrap::AddMCSEffect(AliMUONTrackParam *param, Double_t dZ, Double_t x0)
708 {
709   /// Add to the track parameter covariances the effects of multiple Coulomb scattering
710   /// through a material of thickness "Abs(dZ)" and of radiation length "x0"
711   /// assuming linear propagation and using the small angle approximation.
712   /// dZ = zOut - zIn (sign is important) and "param" is assumed to be given zOut.
713   /// If x0 <= 0., assume dZ = pathLength/x0 and consider the material thickness as negligible.
714   
715   Double_t bendingSlope = param->GetBendingSlope();
716   Double_t nonBendingSlope = param->GetNonBendingSlope();
717   Double_t inverseBendingMomentum = param->GetInverseBendingMomentum();
718   Double_t inverseTotalMomentum2 = inverseBendingMomentum * inverseBendingMomentum *
719                                    (1.0 + bendingSlope * bendingSlope) /
720                                    (1.0 + bendingSlope *bendingSlope + nonBendingSlope * nonBendingSlope); 
721   // Path length in the material
722   Double_t signedPathLength = dZ * TMath::Sqrt(1.0 + bendingSlope*bendingSlope + nonBendingSlope*nonBendingSlope);
723   Double_t pathLengthOverX0 = (x0 > 0.) ? TMath::Abs(signedPathLength) / x0 : TMath::Abs(signedPathLength);
724   // relativistic velocity
725   Double_t velo = 1.;
726   // Angular dispersion square of the track (variance) in a plane perpendicular to the trajectory
727   Double_t theta02 = 0.0136 / velo * (1 + 0.038 * TMath::Log(pathLengthOverX0));
728   theta02 *= theta02 * inverseTotalMomentum2 * pathLengthOverX0;
729   
730   Double_t varCoor      = (x0 > 0.) ? signedPathLength * signedPathLength * theta02 / 3. : 0.;
731   Double_t varSlop      = theta02;
732   Double_t covCorrSlope = (x0 > 0.) ? signedPathLength * theta02 / 2. : 0.;
733   
734   // Set MCS covariance matrix
735   TMatrixD newParamCov(param->GetCovariances());
736   // Non bending plane
737   newParamCov(0,0) += varCoor;       newParamCov(0,1) += covCorrSlope;
738   newParamCov(1,0) += covCorrSlope;  newParamCov(1,1) += varSlop;
739   // Bending plane
740   newParamCov(2,2) += varCoor;       newParamCov(2,3) += covCorrSlope;
741   newParamCov(3,2) += covCorrSlope;  newParamCov(3,3) += varSlop;
742   
743   // Set momentum related covariances if B!=0
744   if (fgFieldON) {
745     // compute derivative d(q/Pxy) / dSlopeX and d(q/Pxy) / dSlopeY
746     Double_t dqPxydSlopeX = inverseBendingMomentum * nonBendingSlope / (1. + nonBendingSlope*nonBendingSlope + bendingSlope*bendingSlope);
747     Double_t dqPxydSlopeY = - inverseBendingMomentum * nonBendingSlope*nonBendingSlope * bendingSlope /
748                               (1. + bendingSlope*bendingSlope) / (1. + nonBendingSlope*nonBendingSlope + bendingSlope*bendingSlope);
749     // Inverse bending momentum (due to dependences with bending and non bending slopes)
750     newParamCov(4,0) += dqPxydSlopeX * covCorrSlope; newParamCov(0,4) += dqPxydSlopeX * covCorrSlope;
751     newParamCov(4,1) += dqPxydSlopeX * varSlop;      newParamCov(1,4) += dqPxydSlopeX * varSlop;
752     newParamCov(4,2) += dqPxydSlopeY * covCorrSlope; newParamCov(2,4) += dqPxydSlopeY * covCorrSlope;
753     newParamCov(4,3) += dqPxydSlopeY * varSlop;      newParamCov(3,4) += dqPxydSlopeY * varSlop;
754     newParamCov(4,4) += (dqPxydSlopeX*dqPxydSlopeX + dqPxydSlopeY*dqPxydSlopeY) * varSlop;
755   }
756   
757   // Set new covariances
758   param->SetCovariances(newParamCov);
759 }
760
761 //__________________________________________________________________________
762 void AliMUONTrackExtrap::ExtrapToVertex(AliMUONTrackParam* trackParam,
763                                         Double_t xVtx, Double_t yVtx, Double_t zVtx,
764                                         Double_t errXVtx, Double_t errYVtx,
765                                         Bool_t correctForMCS, Bool_t correctForEnergyLoss)
766 {
767   /// Main method for extrapolation to the vertex:
768   /// Returns the track parameters and covariances resulting from the extrapolation of the current trackParam
769   /// Changes parameters and covariances according to multiple scattering and energy loss corrections:
770   /// if correctForMCS=kTRUE:  compute parameters using Branson correction and add correction resolution to covariances
771   /// if correctForMCS=kFALSE: add parameter dispersion due to MCS in parameter covariances
772   /// if correctForEnergyLoss=kTRUE:  correct parameters for energy loss and add energy loss fluctuation to covariances
773   /// if correctForEnergyLoss=kFALSE: do nothing about energy loss
774   
775   if (trackParam->GetZ() == zVtx) return; // nothing to be done if already at vertex
776   
777   if (trackParam->GetZ() > zVtx) { // spectro. (z<0)
778     cout<<"E-AliMUONTrackExtrap::ExtrapToVertex: Starting Z ("<<trackParam->GetZ()
779         <<") upstream the vertex (zVtx = "<<zVtx<<")"<<endl;
780     return;
781   }
782   
783   // Check the vertex position relatively to the absorber
784   if (zVtx < AliMUONConstants::AbsZBeg() && zVtx > AliMUONConstants::AbsZEnd()) { // spectro. (z<0)
785     cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Ending Z ("<<zVtx
786         <<") inside the front absorber ("<<AliMUONConstants::AbsZBeg()<<","<<AliMUONConstants::AbsZEnd()<<")"<<endl;
787   } else if (zVtx < AliMUONConstants::AbsZEnd() ) { // spectro. (z<0)
788     cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Ending Z ("<<zVtx
789         <<") downstream the front absorber (zAbsorberEnd = "<<AliMUONConstants::AbsZEnd()<<")"<<endl;
790     if (trackParam->CovariancesExist()) ExtrapToZCov(trackParam,zVtx);
791     else ExtrapToZ(trackParam,zVtx);
792     return;
793   }
794   
795   // Check the track position relatively to the absorber and extrapolate track parameters to the end of the absorber if needed
796   if (trackParam->GetZ() > AliMUONConstants::AbsZBeg()) { // spectro. (z<0)
797     cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Starting Z ("<<trackParam->GetZ()
798         <<") upstream the front absorber (zAbsorberBegin = "<<AliMUONConstants::AbsZBeg()<<")"<<endl;
799     if (trackParam->CovariancesExist()) ExtrapToZCov(trackParam,zVtx);
800     else ExtrapToZ(trackParam,zVtx);
801     return;
802   } else if (trackParam->GetZ() > AliMUONConstants::AbsZEnd()) { // spectro. (z<0)
803     cout<<"W-AliMUONTrackExtrap::ExtrapToVertex: Starting Z ("<<trackParam->GetZ()
804         <<") inside the front absorber ("<<AliMUONConstants::AbsZBeg()<<","<<AliMUONConstants::AbsZEnd()<<")"<<endl;
805   } else {
806     if (trackParam->CovariancesExist()) ExtrapToZCov(trackParam,AliMUONConstants::AbsZEnd());
807     else ExtrapToZ(trackParam,AliMUONConstants::AbsZEnd());
808   }
809   
810   // Get absorber correction parameters assuming linear propagation in absorber
811   Double_t trackXYZOut[3];
812   trackXYZOut[0] = trackParam->GetNonBendingCoor();
813   trackXYZOut[1] = trackParam->GetBendingCoor();
814   trackXYZOut[2] = trackParam->GetZ();
815   Double_t trackXYZIn[3];
816   if (correctForMCS) { // assume linear propagation until the vertex
817     trackXYZIn[2] = TMath::Min(zVtx, AliMUONConstants::AbsZBeg()); // spectro. (z<0)
818     trackXYZIn[0] = trackXYZOut[0] + (xVtx - trackXYZOut[0]) / (zVtx - trackXYZOut[2]) * (trackXYZIn[2] - trackXYZOut[2]);
819     trackXYZIn[1] = trackXYZOut[1] + (yVtx - trackXYZOut[1]) / (zVtx - trackXYZOut[2]) * (trackXYZIn[2] - trackXYZOut[2]);
820   } else {
821     AliMUONTrackParam trackParamIn(*trackParam);
822     ExtrapToZ(&trackParamIn, TMath::Min(zVtx, AliMUONConstants::AbsZBeg()));
823     trackXYZIn[0] = trackParamIn.GetNonBendingCoor();
824     trackXYZIn[1] = trackParamIn.GetBendingCoor();
825     trackXYZIn[2] = trackParamIn.GetZ();
826   }
827   Double_t pTot = trackParam->P();
828   Double_t pathLength, f0, f1, f2, meanRho, totalELoss, sigmaELoss2;
829   if (!GetAbsorberCorrectionParam(trackXYZIn,trackXYZOut,pTot,pathLength,f0,f1,f2,meanRho,totalELoss,sigmaELoss2)) {
830     cout<<"E-AliMUONTrackExtrap::ExtrapToVertex: Unable to take into account the absorber effects"<<endl;
831     if (trackParam->CovariancesExist()) ExtrapToZCov(trackParam,zVtx);
832     else ExtrapToZ(trackParam,zVtx);
833     return;
834   }
835   
836   // Compute track parameters and covariances at vertex according to correctForMCS and correctForEnergyLoss flags
837   if (correctForMCS) {
838     
839     if (correctForEnergyLoss) {
840       
841       // Correct for multiple scattering and energy loss
842       CorrectELossEffectInAbsorber(trackParam, 0.5*totalELoss, 0.5*sigmaELoss2);
843       CorrectMCSEffectInAbsorber(trackParam, xVtx, yVtx, zVtx, errXVtx, errYVtx,
844                                  trackXYZIn[2], pathLength, f0, f1, f2);
845       CorrectELossEffectInAbsorber(trackParam, 0.5*totalELoss, 0.5*sigmaELoss2);
846       
847     } else {
848       
849       // Correct for multiple scattering
850       CorrectMCSEffectInAbsorber(trackParam, xVtx, yVtx, zVtx, errXVtx, errYVtx,
851                                  trackXYZIn[2], pathLength, f0, f1, f2);
852     }
853     
854   } else {
855     
856     if (correctForEnergyLoss) {
857       
858       // Correct for energy loss add multiple scattering dispersion in covariance matrix
859       CorrectELossEffectInAbsorber(trackParam, 0.5*totalELoss, 0.5*sigmaELoss2);
860       AddMCSEffectInAbsorber(trackParam, -pathLength, f0, f1, f2); // (spectro. (z<0))
861       ExtrapToZCov(trackParam, trackXYZIn[2]);
862       CorrectELossEffectInAbsorber(trackParam, 0.5*totalELoss, 0.5*sigmaELoss2);
863       ExtrapToZCov(trackParam, zVtx);
864       
865     } else {
866       
867       // add multiple scattering dispersion in covariance matrix
868       AddMCSEffectInAbsorber(trackParam, -pathLength, f0, f1, f2); // (spectro. (z<0))
869       ExtrapToZCov(trackParam, zVtx);
870       
871     }
872     
873   }
874   
875 }
876
877 //__________________________________________________________________________
878 void AliMUONTrackExtrap::ExtrapToVertex(AliMUONTrackParam* trackParam,
879                                         Double_t xVtx, Double_t yVtx, Double_t zVtx,
880                                         Double_t errXVtx, Double_t errYVtx)
881 {
882   /// Extrapolate track parameters to vertex, corrected for multiple scattering and energy loss effects
883   /// Add branson correction resolution and energy loss fluctuation to parameter covariances
884   ExtrapToVertex(trackParam, xVtx, yVtx, zVtx, errXVtx, errYVtx, kTRUE, kTRUE);
885 }
886
887 //__________________________________________________________________________
888 void AliMUONTrackExtrap::ExtrapToVertexWithoutELoss(AliMUONTrackParam* trackParam,
889                                                     Double_t xVtx, Double_t yVtx, Double_t zVtx,
890                                                     Double_t errXVtx, Double_t errYVtx)
891 {
892   /// Extrapolate track parameters to vertex, corrected for multiple scattering effects only
893   /// Add branson correction resolution to parameter covariances
894   ExtrapToVertex(trackParam, xVtx, yVtx, zVtx, errXVtx, errYVtx, kTRUE, kFALSE);
895 }
896
897 //__________________________________________________________________________
898 void AliMUONTrackExtrap::ExtrapToVertexWithoutBranson(AliMUONTrackParam* trackParam, Double_t zVtx)
899 {
900   /// Extrapolate track parameters to vertex, corrected for energy loss effects only
901   /// Add dispersion due to multiple scattering and energy loss fluctuation to parameter covariances
902   ExtrapToVertex(trackParam, 0., 0., zVtx, 0., 0., kFALSE, kTRUE);
903 }
904
905 //__________________________________________________________________________
906 void AliMUONTrackExtrap::ExtrapToVertexUncorrected(AliMUONTrackParam* trackParam, Double_t zVtx)
907 {
908   /// Extrapolate track parameters to vertex without multiple scattering and energy loss corrections
909   /// Add dispersion due to multiple scattering to parameter covariances
910   ExtrapToVertex(trackParam, 0., 0., zVtx, 0., 0., kFALSE, kFALSE);
911 }
912
913 //__________________________________________________________________________
914 Double_t AliMUONTrackExtrap::TotalMomentumEnergyLoss(AliMUONTrackParam* trackParam, Double_t xVtx, Double_t yVtx, Double_t zVtx)
915 {
916   /// Calculate the total momentum energy loss in-between the track position and the vertex assuming a linear propagation
917   
918   if (trackParam->GetZ() == zVtx) return 0.; // nothing to be done if already at vertex
919   
920   // Check whether the geometry is available
921   if (!gGeoManager) {
922     cout<<"E-AliMUONTrackExtrap::TotalMomentumEnergyLoss: no TGeo"<<endl;
923     return 0.;
924   }
925   
926   // Get encountered material correction parameters assuming linear propagation from vertex to the track position
927   Double_t trackXYZOut[3];
928   trackXYZOut[0] = trackParam->GetNonBendingCoor();
929   trackXYZOut[1] = trackParam->GetBendingCoor();
930   trackXYZOut[2] = trackParam->GetZ();
931   Double_t trackXYZIn[3];
932   trackXYZIn[0] = xVtx;
933   trackXYZIn[1] = yVtx;
934   trackXYZIn[2] = zVtx;
935   Double_t pTot = trackParam->P();
936   Double_t pathLength, f0, f1, f2, meanRho, totalELoss, sigmaELoss2;
937   GetAbsorberCorrectionParam(trackXYZIn,trackXYZOut,pTot,pathLength,f0,f1,f2,meanRho,totalELoss,sigmaELoss2);
938   
939   // total momentum corrected for energy loss
940   Double_t muMass = TDatabasePDG::Instance()->GetParticle("mu-")->Mass(); // GeV
941   Double_t e = TMath::Sqrt(pTot*pTot + muMass*muMass);
942   Double_t eCorr = e + totalELoss;
943   Double_t pTotCorr = TMath::Sqrt(eCorr*eCorr - muMass*muMass);
944   
945   return pTotCorr - pTot;
946 }
947
948 //__________________________________________________________________________
949 Double_t AliMUONTrackExtrap::BetheBloch(Double_t pTotal, Double_t pathLength, Double_t rho, Double_t atomicZ, Double_t atomicZoverA)
950 {
951   /// Returns the mean total momentum energy loss of muon with total momentum='pTotal'
952   /// in the absorber layer of lenght='pathLength', density='rho', A='atomicA' and Z='atomicZ'
953   Double_t muMass = TDatabasePDG::Instance()->GetParticle("mu-")->Mass(); // GeV
954   
955   // mean exitation energy (GeV)
956   Double_t i;
957   if (atomicZ < 13) i = (12. * atomicZ + 7.) * 1.e-9;
958   else i = (9.76 * atomicZ + 58.8 * TMath::Power(atomicZ,-0.19)) * 1.e-9;
959   
960   return pathLength * rho * AliExternalTrackParam::BetheBlochGeant(pTotal/muMass, rho, 0.20, 3.00, i, atomicZoverA);
961 }
962
963 //__________________________________________________________________________
964 Double_t AliMUONTrackExtrap::EnergyLossFluctuation2(Double_t pTotal, Double_t pathLength, Double_t rho, Double_t atomicZoverA)
965 {
966   /// Returns the total momentum energy loss fluctuation of muon with total momentum='pTotal'
967   /// in the absorber layer of lenght='pathLength', density='rho', A='atomicA' and Z='atomicZ'
968   Double_t muMass = TDatabasePDG::Instance()->GetParticle("mu-")->Mass(); // GeV
969   //Double_t eMass = 0.510998918e-3; // GeV
970   Double_t k = 0.307075e-3; // GeV.g^-1.cm^2
971   Double_t p2=pTotal*pTotal;
972   Double_t beta2=p2/(p2 + muMass*muMass);
973   
974   Double_t fwhm = 2. * k * rho * pathLength * atomicZoverA / beta2; // FWHM of the energy loss Landau distribution
975   Double_t sigma2 = fwhm * fwhm / (8.*log(2.)); // gaussian: fwmh = 2 * srqt(2*ln(2)) * sigma (i.e. fwmh = 2.35 * sigma)
976   
977   //sigma2 = k * rho * pathLength * atomicZ / atomicA * eMass; // sigma2 of the energy loss gaussian distribution
978   
979   return sigma2;
980 }
981
982 //__________________________________________________________________________
983 void AliMUONTrackExtrap::Cov2CovP(const TMatrixD &param, TMatrixD &cov)
984 {
985   /// change coordinate system: (X, SlopeX, Y, SlopeY, q/Pyz) -> (X, SlopeX, Y, SlopeY, q*PTot)
986   /// parameters (param) are given in the (X, SlopeX, Y, SlopeY, q/Pyz) coordinate system
987   
988   // charge * total momentum
989   Double_t qPTot = TMath::Sqrt(1. + param(1,0)*param(1,0) + param(3,0)*param(3,0)) /
990                    TMath::Sqrt(1. + param(3,0)*param(3,0)) / param(4,0);
991   
992   // Jacobian of the opposite transformation
993   TMatrixD jacob(5,5);
994   jacob.UnitMatrix();
995   jacob(4,1) = qPTot * param(1,0) / (1. + param(1,0)*param(1,0) + param(3,0)*param(3,0));
996   jacob(4,3) = - qPTot * param(1,0) * param(1,0) * param(3,0) /
997                  (1. + param(3,0)*param(3,0)) / (1. + param(1,0)*param(1,0) + param(3,0)*param(3,0));
998   jacob(4,4) = - qPTot / param(4,0);
999   
1000   // compute covariances in new coordinate system
1001   TMatrixD tmp(cov,TMatrixD::kMultTranspose,jacob);
1002   cov.Mult(jacob,tmp);
1003 }
1004
1005 //__________________________________________________________________________
1006 void AliMUONTrackExtrap::CovP2Cov(const TMatrixD &param, TMatrixD &covP)
1007 {
1008   /// change coordinate system: (X, SlopeX, Y, SlopeY, q*PTot) -> (X, SlopeX, Y, SlopeY, q/Pyz)
1009   /// parameters (param) are given in the (X, SlopeX, Y, SlopeY, q/Pyz) coordinate system
1010   
1011   // charge * total momentum
1012   Double_t qPTot = TMath::Sqrt(1. + param(1,0)*param(1,0) + param(3,0)*param(3,0)) /
1013                    TMath::Sqrt(1. + param(3,0)*param(3,0)) / param(4,0);
1014   
1015   // Jacobian of the transformation
1016   TMatrixD jacob(5,5);
1017   jacob.UnitMatrix();
1018   jacob(4,1) = param(4,0) * param(1,0) / (1. + param(1,0)*param(1,0) + param(3,0)*param(3,0));
1019   jacob(4,3) = - param(4,0) * param(1,0) * param(1,0) * param(3,0) /
1020                  (1. + param(3,0)*param(3,0)) / (1. + param(1,0)*param(1,0) + param(3,0)*param(3,0));
1021   jacob(4,4) = - param(4,0) / qPTot;
1022   
1023   // compute covariances in new coordinate system
1024   TMatrixD tmp(covP,TMatrixD::kMultTranspose,jacob);
1025   covP.Mult(jacob,tmp);
1026 }
1027
1028  //__________________________________________________________________________
1029 void AliMUONTrackExtrap::ExtrapOneStepHelix(Double_t charge, Double_t step, const Double_t *vect, Double_t *vout)
1030 {
1031 /// <pre>
1032 ///    ******************************************************************
1033 ///    *                                                                *
1034 ///    *  Performs the tracking of one step in a magnetic field         *
1035 ///    *  The trajectory is assumed to be a helix in a constant field   *
1036 ///    *  taken at the mid point of the step.                           *
1037 ///    *  Parameters:                                                   *
1038 ///    *   input                                                        *
1039 ///    *     STEP =arc length of the step asked                         *
1040 ///    *     VECT =input vector (position,direction cos and momentum)   *
1041 ///    *     CHARGE=  electric charge of the particle                   *
1042 ///    *   output                                                       *
1043 ///    *     VOUT = same as VECT after completion of the step           *
1044 ///    *                                                                *
1045 ///    *    ==>Called by : USER, GUSWIM                               *
1046 ///    *       Author    m.hansroul  *********                          *
1047 ///    *       modified  s.egli, s.v.levonian                           *
1048 ///    *       modified  v.perevoztchikov
1049 ///    *                                                                *
1050 ///    ******************************************************************
1051 /// </pre>
1052
1053 // modif: everything in double precision
1054
1055     Double_t xyz[3], h[4], hxp[3];
1056     Double_t h2xy, hp, rho, tet;
1057     Double_t sint, sintt, tsint, cos1t;
1058     Double_t f1, f2, f3, f4, f5, f6;
1059
1060     const Int_t kix  = 0;
1061     const Int_t kiy  = 1;
1062     const Int_t kiz  = 2;
1063     const Int_t kipx = 3;
1064     const Int_t kipy = 4;
1065     const Int_t kipz = 5;
1066     const Int_t kipp = 6;
1067
1068     const Double_t kec = 2.9979251e-4;
1069     //
1070     //    ------------------------------------------------------------------
1071     //
1072     //       units are kgauss,centimeters,gev/c
1073     //
1074     vout[kipp] = vect[kipp];
1075     if (TMath::Abs(charge) < 0.00001) {
1076       for (Int_t i = 0; i < 3; i++) {
1077         vout[i] = vect[i] + step * vect[i+3];
1078         vout[i+3] = vect[i+3];
1079       }
1080       return;
1081     }
1082     xyz[0]    = vect[kix] + 0.5 * step * vect[kipx];
1083     xyz[1]    = vect[kiy] + 0.5 * step * vect[kipy];
1084     xyz[2]    = vect[kiz] + 0.5 * step * vect[kipz];
1085
1086     //cmodif: call gufld (xyz, h) changed into:
1087     TGeoGlobalMagField::Instance()->Field(xyz,h);
1088  
1089     h2xy = h[0]*h[0] + h[1]*h[1];
1090     h[3] = h[2]*h[2]+ h2xy;
1091     if (h[3] < 1.e-12) {
1092       for (Int_t i = 0; i < 3; i++) {
1093         vout[i] = vect[i] + step * vect[i+3];
1094         vout[i+3] = vect[i+3];
1095       }
1096       return;
1097     }
1098     if (h2xy < 1.e-12*h[3]) {
1099       ExtrapOneStepHelix3(charge*h[2], step, vect, vout);
1100       return;
1101     }
1102     h[3] = TMath::Sqrt(h[3]);
1103     h[0] /= h[3];
1104     h[1] /= h[3];
1105     h[2] /= h[3];
1106     h[3] *= kec;
1107
1108     hxp[0] = h[1]*vect[kipz] - h[2]*vect[kipy];
1109     hxp[1] = h[2]*vect[kipx] - h[0]*vect[kipz];
1110     hxp[2] = h[0]*vect[kipy] - h[1]*vect[kipx];
1111  
1112     hp = h[0]*vect[kipx] + h[1]*vect[kipy] + h[2]*vect[kipz];
1113
1114     rho = -charge*h[3]/vect[kipp];
1115     tet = rho * step;
1116
1117     if (TMath::Abs(tet) > 0.15) {
1118       sint = TMath::Sin(tet);
1119       sintt = (sint/tet);
1120       tsint = (tet-sint)/tet;
1121       cos1t = 2.*(TMath::Sin(0.5*tet))*(TMath::Sin(0.5*tet))/tet;
1122     } else {
1123       tsint = tet*tet/36.;
1124       sintt = (1. - tsint);
1125       sint = tet*sintt;
1126       cos1t = 0.5*tet;
1127     }
1128
1129     f1 = step * sintt;
1130     f2 = step * cos1t;
1131     f3 = step * tsint * hp;
1132     f4 = -tet*cos1t;
1133     f5 = sint;
1134     f6 = tet * cos1t * hp;
1135  
1136     vout[kix] = vect[kix] + f1*vect[kipx] + f2*hxp[0] + f3*h[0];
1137     vout[kiy] = vect[kiy] + f1*vect[kipy] + f2*hxp[1] + f3*h[1];
1138     vout[kiz] = vect[kiz] + f1*vect[kipz] + f2*hxp[2] + f3*h[2];
1139  
1140     vout[kipx] = vect[kipx] + f4*vect[kipx] + f5*hxp[0] + f6*h[0];
1141     vout[kipy] = vect[kipy] + f4*vect[kipy] + f5*hxp[1] + f6*h[1];
1142     vout[kipz] = vect[kipz] + f4*vect[kipz] + f5*hxp[2] + f6*h[2];
1143  
1144     return;
1145 }
1146
1147  //__________________________________________________________________________
1148 void AliMUONTrackExtrap::ExtrapOneStepHelix3(Double_t field, Double_t step, const Double_t *vect, Double_t *vout)
1149 {
1150 /// <pre>
1151 ///     ******************************************************************
1152 ///     *                                                                *
1153 ///     *       Tracking routine in a constant field oriented            *
1154 ///     *       along axis 3                                             *
1155 ///     *       Tracking is performed with a conventional                *
1156 ///     *       helix step method                                        *
1157 ///     *                                                                *
1158 ///     *    ==>Called by : USER, GUSWIM                                 *
1159 ///     *       Authors    R.Brun, M.Hansroul  *********                 *
1160 ///     *       Rewritten  V.Perevoztchikov
1161 ///     *                                                                *
1162 ///     ******************************************************************
1163 /// </pre>
1164
1165     Double_t hxp[3];
1166     Double_t h4, hp, rho, tet;
1167     Double_t sint, sintt, tsint, cos1t;
1168     Double_t f1, f2, f3, f4, f5, f6;
1169
1170     const Int_t kix  = 0;
1171     const Int_t kiy  = 1;
1172     const Int_t kiz  = 2;
1173     const Int_t kipx = 3;
1174     const Int_t kipy = 4;
1175     const Int_t kipz = 5;
1176     const Int_t kipp = 6;
1177
1178     const Double_t kec = 2.9979251e-4;
1179
1180 // 
1181 //     ------------------------------------------------------------------
1182 // 
1183 //       units are kgauss,centimeters,gev/c
1184 // 
1185     vout[kipp] = vect[kipp];
1186     h4 = field * kec;
1187
1188     hxp[0] = - vect[kipy];
1189     hxp[1] = + vect[kipx];
1190  
1191     hp = vect[kipz];
1192
1193     rho = -h4/vect[kipp];
1194     tet = rho * step;
1195     if (TMath::Abs(tet) > 0.15) {
1196       sint = TMath::Sin(tet);
1197       sintt = (sint/tet);
1198       tsint = (tet-sint)/tet;
1199       cos1t = 2.* TMath::Sin(0.5*tet) * TMath::Sin(0.5*tet)/tet;
1200     } else {
1201       tsint = tet*tet/36.;
1202       sintt = (1. - tsint);
1203       sint = tet*sintt;
1204       cos1t = 0.5*tet;
1205     }
1206
1207     f1 = step * sintt;
1208     f2 = step * cos1t;
1209     f3 = step * tsint * hp;
1210     f4 = -tet*cos1t;
1211     f5 = sint;
1212     f6 = tet * cos1t * hp;
1213  
1214     vout[kix] = vect[kix] + f1*vect[kipx] + f2*hxp[0];
1215     vout[kiy] = vect[kiy] + f1*vect[kipy] + f2*hxp[1];
1216     vout[kiz] = vect[kiz] + f1*vect[kipz] + f3;
1217  
1218     vout[kipx] = vect[kipx] + f4*vect[kipx] + f5*hxp[0];
1219     vout[kipy] = vect[kipy] + f4*vect[kipy] + f5*hxp[1];
1220     vout[kipz] = vect[kipz] + f4*vect[kipz] + f6;
1221
1222     return;
1223 }
1224
1225  //__________________________________________________________________________
1226 Bool_t AliMUONTrackExtrap::ExtrapOneStepRungekutta(Double_t charge, Double_t step, const Double_t* vect, Double_t* vout)
1227 {
1228 /// <pre>
1229 ///     ******************************************************************
1230 ///     *                                                                *
1231 ///     *  Runge-Kutta method for tracking a particle through a magnetic *
1232 ///     *  field. Uses Nystroem algorithm (See Handbook Nat. Bur. of     *
1233 ///     *  Standards, procedure 25.5.20)                                 *
1234 ///     *                                                                *
1235 ///     *  Input parameters                                              *
1236 ///     *       CHARGE    Particle charge                                *
1237 ///     *       STEP      Step size                                      *
1238 ///     *       VECT      Initial co-ords,direction cosines,momentum     *
1239 ///     *  Output parameters                                             *
1240 ///     *       VOUT      Output co-ords,direction cosines,momentum      *
1241 ///     *  User routine called                                           *
1242 ///     *       CALL GUFLD(X,F)                                          *
1243 ///     *                                                                *
1244 ///     *    ==>Called by : USER, GUSWIM                                 *
1245 ///     *       Authors    R.Brun, M.Hansroul  *********                 *
1246 ///     *                  V.Perevoztchikov (CUT STEP implementation)    *
1247 ///     *                                                                *
1248 ///     *                                                                *
1249 ///     ******************************************************************
1250 /// </pre>
1251
1252     Double_t h2, h4, f[4];
1253     Double_t xyzt[3], a, b, c, ph,ph2;
1254     Double_t secxs[4],secys[4],seczs[4],hxp[3];
1255     Double_t g1, g2, g3, g4, g5, g6, ang2, dxt, dyt, dzt;
1256     Double_t est, at, bt, ct, cba;
1257     Double_t f1, f2, f3, f4, rho, tet, hnorm, hp, rho1, sint, cost;
1258     
1259     Double_t x;
1260     Double_t y;
1261     Double_t z;
1262     
1263     Double_t xt;
1264     Double_t yt;
1265     Double_t zt;
1266
1267     Double_t maxit = 1992;
1268     Double_t maxcut = 11;
1269
1270     const Double_t kdlt   = 1e-4;
1271     const Double_t kdlt32 = kdlt/32.;
1272     const Double_t kthird = 1./3.;
1273     const Double_t khalf  = 0.5;
1274     const Double_t kec = 2.9979251e-4;
1275
1276     const Double_t kpisqua = 9.86960440109;
1277     const Int_t kix  = 0;
1278     const Int_t kiy  = 1;
1279     const Int_t kiz  = 2;
1280     const Int_t kipx = 3;
1281     const Int_t kipy = 4;
1282     const Int_t kipz = 5;
1283   
1284     // *.
1285     // *.    ------------------------------------------------------------------
1286     // *.
1287     // *             this constant is for units cm,gev/c and kgauss
1288     // *
1289     Int_t iter = 0;
1290     Int_t ncut = 0;
1291     for(Int_t j = 0; j < 7; j++)
1292       vout[j] = vect[j];
1293
1294     Double_t  pinv   = kec * charge / vect[6];
1295     Double_t tl = 0.;
1296     Double_t h = step;
1297     Double_t rest;
1298
1299  
1300     do {
1301       rest  = step - tl;
1302       if (TMath::Abs(h) > TMath::Abs(rest)) h = rest;
1303       //cmodif: call gufld(vout,f) changed into:
1304       TGeoGlobalMagField::Instance()->Field(vout,f);
1305
1306       // *
1307       // *             start of integration
1308       // *
1309       x      = vout[0];
1310       y      = vout[1];
1311       z      = vout[2];
1312       a      = vout[3];
1313       b      = vout[4];
1314       c      = vout[5];
1315
1316       h2     = khalf * h;
1317       h4     = khalf * h2;
1318       ph     = pinv * h;
1319       ph2    = khalf * ph;
1320       secxs[0] = (b * f[2] - c * f[1]) * ph2;
1321       secys[0] = (c * f[0] - a * f[2]) * ph2;
1322       seczs[0] = (a * f[1] - b * f[0]) * ph2;
1323       ang2 = (secxs[0]*secxs[0] + secys[0]*secys[0] + seczs[0]*seczs[0]);
1324       if (ang2 > kpisqua) break;
1325
1326       dxt    = h2 * a + h4 * secxs[0];
1327       dyt    = h2 * b + h4 * secys[0];
1328       dzt    = h2 * c + h4 * seczs[0];
1329       xt     = x + dxt;
1330       yt     = y + dyt;
1331       zt     = z + dzt;
1332       // *
1333       // *              second intermediate point
1334       // *
1335
1336       est = TMath::Abs(dxt) + TMath::Abs(dyt) + TMath::Abs(dzt);
1337       if (est > h) {
1338         if (ncut++ > maxcut) break;
1339         h *= khalf;
1340         continue;
1341       }
1342  
1343       xyzt[0] = xt;
1344       xyzt[1] = yt;
1345       xyzt[2] = zt;
1346
1347       //cmodif: call gufld(xyzt,f) changed into:
1348       TGeoGlobalMagField::Instance()->Field(xyzt,f);
1349
1350       at     = a + secxs[0];
1351       bt     = b + secys[0];
1352       ct     = c + seczs[0];
1353
1354       secxs[1] = (bt * f[2] - ct * f[1]) * ph2;
1355       secys[1] = (ct * f[0] - at * f[2]) * ph2;
1356       seczs[1] = (at * f[1] - bt * f[0]) * ph2;
1357       at     = a + secxs[1];
1358       bt     = b + secys[1];
1359       ct     = c + seczs[1];
1360       secxs[2] = (bt * f[2] - ct * f[1]) * ph2;
1361       secys[2] = (ct * f[0] - at * f[2]) * ph2;
1362       seczs[2] = (at * f[1] - bt * f[0]) * ph2;
1363       dxt    = h * (a + secxs[2]);
1364       dyt    = h * (b + secys[2]);
1365       dzt    = h * (c + seczs[2]);
1366       xt     = x + dxt;
1367       yt     = y + dyt;
1368       zt     = z + dzt;
1369       at     = a + 2.*secxs[2];
1370       bt     = b + 2.*secys[2];
1371       ct     = c + 2.*seczs[2];
1372
1373       est = TMath::Abs(dxt)+TMath::Abs(dyt)+TMath::Abs(dzt);
1374       if (est > 2.*TMath::Abs(h)) {
1375         if (ncut++ > maxcut) break;
1376         h *= khalf;
1377         continue;
1378       }
1379  
1380       xyzt[0] = xt;
1381       xyzt[1] = yt;
1382       xyzt[2] = zt;
1383
1384       //cmodif: call gufld(xyzt,f) changed into:
1385       TGeoGlobalMagField::Instance()->Field(xyzt,f);
1386
1387       z      = z + (c + (seczs[0] + seczs[1] + seczs[2]) * kthird) * h;
1388       y      = y + (b + (secys[0] + secys[1] + secys[2]) * kthird) * h;
1389       x      = x + (a + (secxs[0] + secxs[1] + secxs[2]) * kthird) * h;
1390
1391       secxs[3] = (bt*f[2] - ct*f[1])* ph2;
1392       secys[3] = (ct*f[0] - at*f[2])* ph2;
1393       seczs[3] = (at*f[1] - bt*f[0])* ph2;
1394       a      = a+(secxs[0]+secxs[3]+2. * (secxs[1]+secxs[2])) * kthird;
1395       b      = b+(secys[0]+secys[3]+2. * (secys[1]+secys[2])) * kthird;
1396       c      = c+(seczs[0]+seczs[3]+2. * (seczs[1]+seczs[2])) * kthird;
1397
1398       est    = TMath::Abs(secxs[0]+secxs[3] - (secxs[1]+secxs[2]))
1399         + TMath::Abs(secys[0]+secys[3] - (secys[1]+secys[2]))
1400         + TMath::Abs(seczs[0]+seczs[3] - (seczs[1]+seczs[2]));
1401
1402       if (est > kdlt && TMath::Abs(h) > 1.e-4) {
1403         if (ncut++ > maxcut) break;
1404         h *= khalf;
1405         continue;
1406       }
1407
1408       ncut = 0;
1409       // *               if too many iterations, go to helix
1410       if (iter++ > maxit) break;
1411
1412       tl += h;
1413       if (est < kdlt32) 
1414         h *= 2.;
1415       cba    = 1./ TMath::Sqrt(a*a + b*b + c*c);
1416       vout[0] = x;
1417       vout[1] = y;
1418       vout[2] = z;
1419       vout[3] = cba*a;
1420       vout[4] = cba*b;
1421       vout[5] = cba*c;
1422       rest = step - tl;
1423       if (step < 0.) rest = -rest;
1424       if (rest < 1.e-5*TMath::Abs(step)) return kTRUE;
1425
1426     } while(1);
1427
1428     // angle too big, use helix
1429     cout<<"W-AliMUONTrackExtrap::ExtrapOneStepRungekutta: Ruge-Kutta failed: switch to helix"<<endl;
1430
1431     f1  = f[0];
1432     f2  = f[1];
1433     f3  = f[2];
1434     f4  = TMath::Sqrt(f1*f1+f2*f2+f3*f3);
1435     if (f4 < 1.e-10) {
1436       cout<<"E-AliMUONTrackExtrap::ExtrapOneStepRungekutta: magnetic field at (";
1437       cout<<xyzt[0]<<", "<<xyzt[1]<<", "<<xyzt[2]<<") = "<<f4<<": giving up"<<endl;
1438       return kFALSE;
1439     }
1440     rho = -f4*pinv;
1441     tet = rho * step;
1442  
1443     hnorm = 1./f4;
1444     f1 = f1*hnorm;
1445     f2 = f2*hnorm;
1446     f3 = f3*hnorm;
1447
1448     hxp[0] = f2*vect[kipz] - f3*vect[kipy];
1449     hxp[1] = f3*vect[kipx] - f1*vect[kipz];
1450     hxp[2] = f1*vect[kipy] - f2*vect[kipx];
1451  
1452     hp = f1*vect[kipx] + f2*vect[kipy] + f3*vect[kipz];
1453
1454     rho1 = 1./rho;
1455     sint = TMath::Sin(tet);
1456     cost = 2.*TMath::Sin(khalf*tet)*TMath::Sin(khalf*tet);
1457
1458     g1 = sint*rho1;
1459     g2 = cost*rho1;
1460     g3 = (tet-sint) * hp*rho1;
1461     g4 = -cost;
1462     g5 = sint;
1463     g6 = cost * hp;
1464  
1465     vout[kix] = vect[kix] + g1*vect[kipx] + g2*hxp[0] + g3*f1;
1466     vout[kiy] = vect[kiy] + g1*vect[kipy] + g2*hxp[1] + g3*f2;
1467     vout[kiz] = vect[kiz] + g1*vect[kipz] + g2*hxp[2] + g3*f3;
1468  
1469     vout[kipx] = vect[kipx] + g4*vect[kipx] + g5*hxp[0] + g6*f1;
1470     vout[kipy] = vect[kipy] + g4*vect[kipy] + g5*hxp[1] + g6*f2;
1471     vout[kipz] = vect[kipz] + g4*vect[kipz] + g5*hxp[2] + g6*f3;
1472
1473     return kTRUE;
1474 }
1475