- Double_t r2 = TMath::Sqrt(1.0 - 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;
- Double_t cc = c1+c2;
- Double_t 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;
- Double_t b01 = f12*fCey + f14*fCcy + f13*fCty;
- Double_t b10 = f02*fCez + f04*fCcz;
- Double_t b11 = f12*fCez + f14*fCcz + f13*fCtz;
- Double_t b20 = f02*fCee + f04*fCce;
- Double_t b21 = f12*fCee + f14*fCce + f13*fCte;
- Double_t b30 = f02*fCte + f04*fCct;
- Double_t b31 = f12*fCte + f14*fCct + f13*fCtt;
- Double_t b40 = f02*fCce + f04*fCcc;
- Double_t b41 = f12*fCce + f14*fCcc + f13*fCct;
-
- // a = f*b = f*C*ft
- Double_t a00 = f02*b20 + f04*b40;
- Double_t a01 = f02*b21 + f04*b41;
- Double_t a11 = f12*b21 + f14*b41 + f13*b31;
-
- // F*C*Ft = C + (a + b + bt)
- fCyy += a00 + 2.0*b00;
- fCzy += a01 + b01 + b10;
- fCey += b20;
- fCty += b30;
- fCcy += b40;
- fCzz += a11 + 2.0*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
- // What is 14.1 ????
- Double_t d = TMath::Sqrt((x1-fX)*(x1-fX) + (y1-fY)*(y1-fY) + (z1-fZ)*(z1-fZ));
- Double_t p2 = (1.0 + 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;
- Double_t ez = fT;
- Double_t xz = fC*ez;
- Double_t zz1 = ez*ez + 1.0;
- Double_t xy = fE + ey;
-
- fCee += (2.0*ey*ez*ez*fE + 1.0 - 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;
-
- // Energy losses
- // What is 5940.0 ???? and 0.153e-3 ????
- if ((5940.0*beta2 / (1.0 - beta2 + 1e-10) - beta2) < 0.0) {
- return 0;
- }
- Double_t dE = 0.153e-3/beta2 * (TMath::Log(5940.0*beta2 / (1.0 - beta2 + 1e-10)) - beta2)
- * d * rho;
- Float_t budget = d * rho;
- fBudget[0] +=budget;
-
- // Suspicious part - think about it ????
- Double_t kinE = TMath::Sqrt(p2);
- if (dE > 0.8 * kinE) {
- dE = 0.8 * kinE;
- }
- if (dE < 0.0) {
- dE = 0.0; // Not valid region for Bethe bloch
- }
- fDE += dE;
- if (x1 < x2) {
- dE = -dE;
- }
- cc = fC;
- fC *= (1.0 - TMath::Sqrt(p2 + GetMass()*GetMass()) / p2 * dE);
- fE += fX * (fC - cc);
-
- // Energy loss fluctuation
- // Why 0.07 ????
- Double_t sigmade = 0.07 * TMath::Sqrt(TMath::Abs(dE));
- Double_t sigmac2 = sigmade*sigmade * fC*fC * (p2 + GetMass()*GetMass()) / (p2*p2);
- fCcc += sigmac2;
- fCee += fX*fX * sigmac2;
-
- // Track time measurement
- if (x1 < x2) {
- if (IsStartedTimeIntegral()) {
- Double_t l2 = TMath::Sqrt((fX-oldX)*(fX-oldX)
- + (fY-oldY)*(fY-oldY)
- + (fZ-oldZ)*(fZ-oldZ));
- if (TMath::Abs(l2*fC) > 0.0001){
- // <ake correction for curvature if neccesary
- l2 = 0.5 * TMath::Sqrt((fX-oldX)*(fX-oldX)
- + (fY-oldY)*(fY-oldY));
- l2 = 2.0 * TMath::ASin(l2 * fC) / fC;
- l2 = TMath::Sqrt(l2*l2 + (fZ-oldZ)*(fZ-oldZ));