- Double_t x1=fX, x2=x1+(xk-x1), dx=x2-x1, y1=fY, z1=fZ;
- Double_t c1=fC*x1 - fE;
- if((c1*c1) > 1){
- return 0;}
- Double_t r1=sqrt(1.- c1*c1);
- Double_t c2=fC*x2 - fE;
- if((c2*c2) > 1) {
- return 0;
- }
- Double_t r2=sqrt(1.- c2*c2);
-
- fY += dx*(c1+c2)/(r1+r2);
- fZ += dx*(c1+c2)/(c1*r2 + c2*r1)*fT;
-
- //f = F - 1
- Double_t rr=r1+r2, cc=c1+c2, xx=x1+x2;
- Double_t f02=-dx*(2*rr + cc*(c1/r1 + c2/r2))/(rr*rr);
- Double_t f04= dx*(rr*xx + cc*(c1*x1/r1+c2*x2/r2))/(rr*rr);
- Double_t cr=c1*r2+c2*r1;
- Double_t f12=-dx*fT*(2*cr + cc*(c2*c1/r1-r1 + c1*c2/r2-r2))/(cr*cr);
- Double_t f13= dx*cc/cr;
- Double_t f14=dx*fT*(cr*xx-cc*(r1*x2-c2*c1*x1/r1+r2*x1-c1*c2*x2/r2))/(cr*cr);
-
- //b = C*ft
- Double_t b00=f02*fCey + f04*fCcy, b01=f12*fCey + f14*fCcy + f13*fCty;
- Double_t b10=f02*fCez + f04*fCcz, b11=f12*fCez + f14*fCcz + f13*fCtz;
- Double_t b20=f02*fCee + f04*fCce, b21=f12*fCee + f14*fCce + f13*fCte;
- Double_t b30=f02*fCte + f04*fCct, b31=f12*fCte + f14*fCct + f13*fCtt;
- Double_t b40=f02*fCce + f04*fCcc, b41=f12*fCce + f14*fCcc + f13*fCct;
-
- //a = f*b = f*C*ft
- Double_t a00=f02*b20+f04*b40,a01=f02*b21+f04*b41,a11=f12*b21+f14*b41+f13*b31;
-
- //F*C*Ft = C + (a + b + bt)
- fCyy += a00 + 2*b00;
- fCzy += a01 + b01 + b10;
- fCey += b20;
- fCty += b30;
- fCcy += b40;
- fCzz += a11 + 2*b11;
- fCez += b21;
- fCtz += b31;
- fCcz += b41;
-
- fX=x2;
-
- //Change of the magnetic field *************
- SaveLocalConvConst();
- cc=fC;
- fC*=lcc/GetLocalConvConst();
- fE+=fX*(fC-cc);
-
- //Multiple scattering ******************
- Double_t d=sqrt((x1-fX)*(x1-fX)+(y1-fY)*(y1-fY)+(z1-fZ)*(z1-fZ));
- Double_t p2=(1.+ GetTgl()*GetTgl())/(Get1Pt()*Get1Pt());
- Double_t beta2=p2/(p2 + GetMass()*GetMass());
- Double_t theta2=14.1*14.1/(beta2*p2*1e6)*d/x0*rho;
-
- Double_t ey=fC*fX - fE, ez=fT;
- Double_t xz=fC*ez, zz1=ez*ez+1, xy=fE+ey;
-
- fCee += (2*ey*ez*ez*fE+1-ey*ey+ez*ez+fE*fE*ez*ez)*theta2;
- fCte += ez*zz1*xy*theta2;
- fCtt += zz1*zz1*theta2;
- fCce += xz*ez*xy*theta2;
- fCct += xz*zz1*theta2;
- fCcc += xz*xz*theta2;