```///    ******************************************************************
///    *                                                                *
///    *  Performs the tracking of one step in a magnetic field         *
///    *  The trajectory is assumed to be a helix in a constant field   *
///    *  taken at the mid point of the step.                           *
///    *  Parameters:                                                   *
///    *   input                                                        *
///    *     STEP =arc length of the step asked                         *
///    *     VECT =input vector (position,direction cos and momentum)   *
///    *     CHARGE=  electric charge of the particle                   *
///    *   output                                                       *
///    *     VOUT = same as VECT after completion of the step           *
///    *                                                                *
///    *    ==>Called by : USER, GUSWIM                               *
///    *       Author    m.hansroul  *********                          *
///    *       modified  s.egli, s.v.levonian                           *
///    *       modified  v.perevoztchikov
///    *                                                                *
///    ******************************************************************
/// ```
// modif: everything in double precision Double_t xyz[3], h[4], hxp[3]; Double_t h2xy, hp, rho, tet; Double_t sint, sintt, tsint, cos1t; Double_t f1, f2, f3, f4, f5, f6; const Int_t kix = 0; const Int_t kiy = 1; const Int_t kiz = 2; const Int_t kipx = 3; const Int_t kipy = 4; const Int_t kipz = 5; const Int_t kipp = 6; const Double_t kec = 2.9979251e-4; // // ------------------------------------------------------------------ // // units are kgauss,centimeters,gev/c // vout[kipp] = vect[kipp]; if (TMath::Abs(charge) < 0.00001) { for (Int_t i = 0; i < 3; i++) { vout[i] = vect[i] + step * vect[i+3]; vout[i+3] = vect[i+3]; } return; } xyz[0] = vect[kix] + 0.5 * step * vect[kipx]; xyz[1] = vect[kiy] + 0.5 * step * vect[kipy]; xyz[2] = vect[kiz] + 0.5 * step * vect[kipz]; //cmodif: call gufld (xyz, h) changed into: TGeoGlobalMagField::Instance()->Field(xyz,h); h2xy = h[0]*h[0] + h[1]*h[1]; h[3] = h[2]*h[2]+ h2xy; if (h[3] < 1.e-12) { for (Int_t i = 0; i < 3; i++) { vout[i] = vect[i] + step * vect[i+3]; vout[i+3] = vect[i+3]; } return; } if (h2xy < 1.e-12*h[3]) { ExtrapOneStepHelix3(charge*h[2], step, vect, vout); return; } h[3] = TMath::Sqrt(h[3]); h[0] /= h[3]; h[1] /= h[3]; h[2] /= h[3]; h[3] *= kec; hxp[0] = h[1]*vect[kipz] - h[2]*vect[kipy]; hxp[1] = h[2]*vect[kipx] - h[0]*vect[kipz]; hxp[2] = h[0]*vect[kipy] - h[1]*vect[kipx]; hp = h[0]*vect[kipx] + h[1]*vect[kipy] + h[2]*vect[kipz]; rho = -charge*h[3]/vect[kipp]; tet = rho * step; if (TMath::Abs(tet) > 0.15) { sint = TMath::Sin(tet); sintt = (sint/tet); tsint = (tet-sint)/tet; cos1t = 2.*(TMath::Sin(0.5*tet))*(TMath::Sin(0.5*tet))/tet; } else { tsint = tet*tet/36.; sintt = (1. - tsint); sint = tet*sintt; cos1t = 0.5*tet; } f1 = step * sintt; f2 = step * cos1t; f3 = step * tsint * hp; f4 = -tet*cos1t; f5 = sint; f6 = tet * cos1t * hp; vout[kix] = vect[kix] + f1*vect[kipx] + f2*hxp[0] + f3*h[0]; vout[kiy] = vect[kiy] + f1*vect[kipy] + f2*hxp[1] + f3*h[1]; vout[kiz] = vect[kiz] + f1*vect[kipz] + f2*hxp[2] + f3*h[2]; vout[kipx] = vect[kipx] + f4*vect[kipx] + f5*hxp[0] + f6*h[0]; vout[kipy] = vect[kipy] + f4*vect[kipy] + f5*hxp[1] + f6*h[1]; vout[kipz] = vect[kipz] + f4*vect[kipz] + f5*hxp[2] + f6*h[2]; return; } //__________________________________________________________________________ void AliMUONTrackExtrap::ExtrapOneStepHelix3(Double_t field, Double_t step, const Double_t *vect, Double_t *vout) { ///
```///	******************************************************************
///	*								 *
///	*	Tracking routine in a constant field oriented		 *
///	*	along axis 3						 *
///	*	Tracking is performed with a conventional		 *
///	*	helix step method					 *
///	*								 *
///	*    ==>Called by : USER, GUSWIM				 *
///	*	Authors    R.Brun, M.Hansroul  *********		 *
///	*	Rewritten  V.Perevoztchikov
///	*								 *
///	******************************************************************
/// ```
Double_t hxp[3]; Double_t h4, hp, rho, tet; Double_t sint, sintt, tsint, cos1t; Double_t f1, f2, f3, f4, f5, f6; const Int_t kix = 0; const Int_t kiy = 1; const Int_t kiz = 2; const Int_t kipx = 3; const Int_t kipy = 4; const Int_t kipz = 5; const Int_t kipp = 6; const Double_t kec = 2.9979251e-4; // // ------------------------------------------------------------------ // // units are kgauss,centimeters,gev/c // vout[kipp] = vect[kipp]; h4 = field * kec; hxp[0] = - vect[kipy]; hxp[1] = + vect[kipx]; hp = vect[kipz]; rho = -h4/vect[kipp]; tet = rho * step; if (TMath::Abs(tet) > 0.15) { sint = TMath::Sin(tet); sintt = (sint/tet); tsint = (tet-sint)/tet; cos1t = 2.* TMath::Sin(0.5*tet) * TMath::Sin(0.5*tet)/tet; } else { tsint = tet*tet/36.; sintt = (1. - tsint); sint = tet*sintt; cos1t = 0.5*tet; } f1 = step * sintt; f2 = step * cos1t; f3 = step * tsint * hp; f4 = -tet*cos1t; f5 = sint; f6 = tet * cos1t * hp; vout[kix] = vect[kix] + f1*vect[kipx] + f2*hxp[0]; vout[kiy] = vect[kiy] + f1*vect[kipy] + f2*hxp[1]; vout[kiz] = vect[kiz] + f1*vect[kipz] + f3; vout[kipx] = vect[kipx] + f4*vect[kipx] + f5*hxp[0]; vout[kipy] = vect[kipy] + f4*vect[kipy] + f5*hxp[1]; vout[kipz] = vect[kipz] + f4*vect[kipz] + f6; return; } //__________________________________________________________________________ Bool_t AliMUONTrackExtrap::ExtrapOneStepRungekutta(Double_t charge, Double_t step, const Double_t* vect, Double_t* vout) { ///
```///	******************************************************************
///	*								 *
///	*  Runge-Kutta method for tracking a particle through a magnetic *
///	*  field. Uses Nystroem algorithm (See Handbook Nat. Bur. of	 *
///	*  Standards, procedure 25.5.20)				 *
///	*								 *
///	*  Input parameters						 *
///	*	CHARGE    Particle charge				 *
///	*	STEP	  Step size					 *
///	*	VECT	  Initial co-ords,direction cosines,momentum	 *
///	*  Output parameters						 *
///	*	VOUT	  Output co-ords,direction cosines,momentum	 *
///	*  User routine called  					 *
///	*	CALL GUFLD(X,F) 					 *
///	*								 *
///	*    ==>Called by : USER, GUSWIM				 *
///	*	Authors    R.Brun, M.Hansroul  *********		 *
///	*		   V.Perevoztchikov (CUT STEP implementation)	 *
///	*								 *
///	*								 *
///	******************************************************************
/// ```
Double_t h2, h4, f[4]; Double_t xyzt[3] = {FLT_MAX, FLT_MAX, FLT_MAX}; Double_t a, b, c, ph,ph2; Double_t secxs[4],secys[4],seczs[4],hxp[3]; Double_t g1, g2, g3, g4, g5, g6, ang2, dxt, dyt, dzt; Double_t est, at, bt, ct, cba; Double_t f1, f2, f3, f4, rho, tet, hnorm, hp, rho1, sint, cost; Double_t x; Double_t y; Double_t z; Double_t xt; Double_t yt; Double_t zt; Double_t maxit = 1992; Double_t maxcut = 11; const Double_t kdlt = 1e-4; const Double_t kdlt32 = kdlt/32.; const Double_t kthird = 1./3.; const Double_t khalf = 0.5; const Double_t kec = 2.9979251e-4; const Double_t kpisqua = 9.86960440109; const Int_t kix = 0; const Int_t kiy = 1; const Int_t kiz = 2; const Int_t kipx = 3; const Int_t kipy = 4; const Int_t kipz = 5; // *. // *. ------------------------------------------------------------------ // *. // * this constant is for units cm,gev/c and kgauss // * Int_t iter = 0; Int_t ncut = 0; for(Int_t j = 0; j < 7; j++) vout[j] = vect[j]; Double_t pinv = kec * charge / vect[6]; Double_t tl = 0.; Double_t h = step; Double_t rest; do { rest = step - tl; if (TMath::Abs(h) > TMath::Abs(rest)) h = rest; //cmodif: call gufld(vout,f) changed into: TGeoGlobalMagField::Instance()->Field(vout,f); // * // * start of integration // * x = vout[0]; y = vout[1]; z = vout[2]; a = vout[3]; b = vout[4]; c = vout[5]; h2 = khalf * h; h4 = khalf * h2; ph = pinv * h; ph2 = khalf * ph; secxs[0] = (b * f[2] - c * f[1]) * ph2; secys[0] = (c * f[0] - a * f[2]) * ph2; seczs[0] = (a * f[1] - b * f[0]) * ph2; ang2 = (secxs[0]*secxs[0] + secys[0]*secys[0] + seczs[0]*seczs[0]); if (ang2 > kpisqua) break; dxt = h2 * a + h4 * secxs[0]; dyt = h2 * b + h4 * secys[0]; dzt = h2 * c + h4 * seczs[0]; xt = x + dxt; yt = y + dyt; zt = z + dzt; // * // * second intermediate point // * est = TMath::Abs(dxt) + TMath::Abs(dyt) + TMath::Abs(dzt); if (est > h) { if (ncut++ > maxcut) break; h *= khalf; continue; } xyzt[0] = xt; xyzt[1] = yt; xyzt[2] = zt; //cmodif: call gufld(xyzt,f) changed into: TGeoGlobalMagField::Instance()->Field(xyzt,f); at = a + secxs[0]; bt = b + secys[0]; ct = c + seczs[0]; secxs[1] = (bt * f[2] - ct * f[1]) * ph2; secys[1] = (ct * f[0] - at * f[2]) * ph2; seczs[1] = (at * f[1] - bt * f[0]) * ph2; at = a + secxs[1]; bt = b + secys[1]; ct = c + seczs[1]; secxs[2] = (bt * f[2] - ct * f[1]) * ph2; secys[2] = (ct * f[0] - at * f[2]) * ph2; seczs[2] = (at * f[1] - bt * f[0]) * ph2; dxt = h * (a + secxs[2]); dyt = h * (b + secys[2]); dzt = h * (c + seczs[2]); xt = x + dxt; yt = y + dyt; zt = z + dzt; at = a + 2.*secxs[2]; bt = b + 2.*secys[2]; ct = c + 2.*seczs[2]; est = TMath::Abs(dxt)+TMath::Abs(dyt)+TMath::Abs(dzt); if (est > 2.*TMath::Abs(h)) { if (ncut++ > maxcut) break; h *= khalf; continue; } xyzt[0] = xt; xyzt[1] = yt; xyzt[2] = zt; //cmodif: call gufld(xyzt,f) changed into: TGeoGlobalMagField::Instance()->Field(xyzt,f); z = z + (c + (seczs[0] + seczs[1] + seczs[2]) * kthird) * h; y = y + (b + (secys[0] + secys[1] + secys[2]) * kthird) * h; x = x + (a + (secxs[0] + secxs[1] + secxs[2]) * kthird) * h; secxs[3] = (bt*f[2] - ct*f[1])* ph2; secys[3] = (ct*f[0] - at*f[2])* ph2; seczs[3] = (at*f[1] - bt*f[0])* ph2; a = a+(secxs[0]+secxs[3]+2. * (secxs[1]+secxs[2])) * kthird; b = b+(secys[0]+secys[3]+2. * (secys[1]+secys[2])) * kthird; c = c+(seczs[0]+seczs[3]+2. * (seczs[1]+seczs[2])) * kthird; est = TMath::Abs(secxs[0]+secxs[3] - (secxs[1]+secxs[2])) + TMath::Abs(secys[0]+secys[3] - (secys[1]+secys[2])) + TMath::Abs(seczs[0]+seczs[3] - (seczs[1]+seczs[2])); if (est > kdlt && TMath::Abs(h) > 1.e-4) { if (ncut++ > maxcut) break; h *= khalf; continue; } ncut = 0; // * if too many iterations, go to helix if (iter++ > maxit) break; tl += h; if (est < kdlt32) h *= 2.; cba = 1./ TMath::Sqrt(a*a + b*b + c*c); vout[0] = x; vout[1] = y; vout[2] = z; vout[3] = cba*a; vout[4] = cba*b; vout[5] = cba*c; rest = step - tl; if (step < 0.) rest = -rest; if (rest < 1.e-5*TMath::Abs(step)) return kTRUE; } while(1); // angle too big, use helix cout<<"W-AliMUONTrackExtrap::ExtrapOneStepRungekutta: Ruge-Kutta failed: switch to helix"<