- Double_t alpha=fAlpha, x=fX;
- Double_t p0=fP0,p1=fP1,p2=fP2,p3=fP3,p4=fP4;
- Double_t c00=fC00;
- Double_t c10=fC10, c11=fC11;
- Double_t c20=fC20, c21=fC21, c22=fC22;
- Double_t c30=fC30, c31=fC31, c32=fC32, c33=fC33;
- Double_t c40=fC40, c41=fC41, c42=fC42, c43=fC43, c44=fC44;
-
- if (alp < -TMath::Pi()) alp += 2*TMath::Pi();
- else if (alp >= TMath::Pi()) alp -= 2*TMath::Pi();
- Double_t ca=TMath::Cos(alp-fAlpha), sa=TMath::Sin(alp-fAlpha);
- Double_t sf=fP2, cf=TMath::Sqrt(1.- fP2*fP2);
-
- TMatrixD *T=0;
- // **** rotation **********************
- {
- fAlpha = alp;
- fX = x*ca + p0*sa;
- fP0= -x*sa + p0*ca;
- fP2= sf*ca - cf*sa;
-
- TMatrixD C(5,5);
- C(0,0)=c00;
- C(1,0)=c10; C(1,1)=c11;
- C(2,0)=c20; C(2,1)=c21; C(2,2)=c22;
- C(3,0)=c30; C(3,1)=c31; C(3,2)=c32; C(3,3)=c33;
- C(4,0)=c40; C(4,1)=c41; C(4,2)=c42; C(4,3)=c43; C(4,4)=c44;
- C(0,1)=C(1,0);
- C(0,2)=C(2,0); C(1,2)=C(2,1);
- C(0,3)=C(3,0); C(1,3)=C(3,1); C(2,3)=C(3,2);
- C(0,4)=C(4,0); C(1,4)=C(4,1); C(2,4)=C(4,2); C(3,4)=C(4,3);
-
- TMatrixD F(6,5);
- F(0,0)=sa;
- F(1,0)=ca;
- F(2,1)=F(4,3)=F(5,4)=1;
- F(3,2)=ca + sf/cf*sa;
-
- TMatrixD tmp(C,TMatrixD::kMult,TMatrixD(TMatrixD::kTransposed, F));
- T=new TMatrixD(F,TMatrixD::kMult,tmp);
- }
-
- // **** translation ******************
- {
- Double_t dx=xk-fX;
- Double_t f1=fP2, f2=f1 + fP4*dx;
- if (TMath::Abs(f2) >= 0.9999) {
- Int_t n=GetNumberOfClusters();
- if (n>kWARN)
- Warning("Propagate","Propagation failed (%d) !\n",n);
- return 0;
- }
- Double_t r1=TMath::Sqrt(1.- f1*f1), r2=TMath::Sqrt(1.- f2*f2);
-
- fX=xk;
- fP0 += dx*(f1+f2)/(r1+r2);
- fP1 += dx*(f1+f2)/(f1*r2 + f2*r1)*fP3;
- fP2 += dx*fP4;
-
- TMatrixD F(5,6);
- F(0,1)=F(1,2)=F(2,3)=F(3,4)=F(4,5)=1;
- F(0,3)=dx/(r1+r2)*(2+(f1+f2)*(f2/r2+f1/r1)/(r1+r2));
- F(0,5)=dx*dx/(r1+r2)*(1+(f1+f2)*f2/(r1+r2));
- F(1,3)=dx*fP3/(f1*r2 + f2*r1)*(2-(f1+f2)*(r2-f1*f2/r2+r1-f2*f1/r1)/(f1*r2 + f2*r1));
- F(1,4)=dx*(f1+f2)/(f1*r2 + f2*r1);
- F(1,5)=dx*dx*fP3/(f1*r2 + f2*r1)*(1-(f1+f2)*(-f1*f2/r2+r1)/(f1*r2 + f2*r1));
- F(2,5)=dx;
- F(0,0)=-1/(r1+r2)*((f1+f2)+dx*fP4*(1+(f1+f2)/(r1+r2)*f2/r2));
- F(1,0)=-fP3/(f1*r2 + f2*r1)*((f1+f2)+dx*fP4*(1+(f1+f2)/(f1*r2 + f2*r1)*(f1*f2/r2-r1)));
- F(2,0)=-fP4;
-
- TMatrixD tmp(*T,TMatrixD::kMult,TMatrixD(TMatrixD::kTransposed, F));
- delete T;
- TMatrixD C(F,TMatrixD::kMult,tmp);
-
- fC00=C(0,0);
- fC10=C(1,0); fC11=C(1,1);
- fC20=C(2,0); fC21=C(2,1); fC22=C(2,2);
- fC30=C(3,0); fC31=C(3,1); fC32=C(3,2); fC33=C(3,3);
- fC40=C(4,0); fC41=C(4,1); fC42=C(4,2); fC43=C(4,3); fC44=C(4,4);