- 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 *tT=0;
- // **** rotation **********************
- {
- fAlpha = alp;
- fX = x*ca + p0*sa;
- fP0= -x*sa + p0*ca;
- fP2= sf*ca - cf*sa;
-
- TMatrixD cC(5,5);
- cC(0,0)=c00;
- cC(1,0)=c10; cC(1,1)=c11;
- cC(2,0)=c20; cC(2,1)=c21; cC(2,2)=c22;
- cC(3,0)=c30; cC(3,1)=c31; cC(3,2)=c32; cC(3,3)=c33;
- cC(4,0)=c40; cC(4,1)=c41; cC(4,2)=c42; cC(4,3)=c43; cC(4,4)=c44;
- cC(0,1)=cC(1,0);
- cC(0,2)=cC(2,0); cC(1,2)=cC(2,1);
- cC(0,3)=cC(3,0); cC(1,3)=cC(3,1); cC(2,3)=cC(3,2);
- cC(0,4)=cC(4,0); cC(1,4)=cC(4,1); cC(2,4)=cC(4,2); cC(3,4)=cC(4,3);
-
- TMatrixD mF(6,5);
- mF(0,0)=sa;
- mF(1,0)=ca;
- mF(2,1)=mF(4,3)=mF(5,4)=1;
- mF(3,2)=ca + sf/cf*sa;
-
- TMatrixD tmp(cC,TMatrixD::kMult,TMatrixD(TMatrixD::kTransposed, mF));
- tT=new TMatrixD(mF,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 mF(5,6);
- mF(0,1)=mF(1,2)=mF(2,3)=mF(3,4)=mF(4,5)=1;
- mF(0,3)=dx/(r1+r2)*(2+(f1+f2)*(f2/r2+f1/r1)/(r1+r2));
- mF(0,5)=dx*dx/(r1+r2)*(1+(f1+f2)*f2/(r1+r2));
- mF(1,3)=dx*fP3/(f1*r2 + f2*r1)*(2-(f1+f2)*(r2-f1*f2/r2+r1-f2*f1/r1)/(f1*r2 + f2*r1));
- mF(1,4)=dx*(f1+f2)/(f1*r2 + f2*r1);
- mF(1,5)=dx*dx*fP3/(f1*r2 + f2*r1)*(1-(f1+f2)*(-f1*f2/r2+r1)/(f1*r2 + f2*r1));
- mF(2,5)=dx;
- mF(0,0)=-1/(r1+r2)*((f1+f2)+dx*fP4*(1+(f1+f2)/(r1+r2)*f2/r2));
- mF(1,0)=-fP3/(f1*r2 + f2*r1)*((f1+f2)+dx*fP4*(1+(f1+f2)/(f1*r2 + f2*r1)*(f1*f2/r2-r1)));
- mF(2,0)=-fP4;
-
- TMatrixD tmp(*tT,TMatrixD::kMult,TMatrixD(TMatrixD::kTransposed, mF));
- delete tT;
- TMatrixD cC(mF,TMatrixD::kMult,tmp);
-
- fC00=cC(0,0);
- fC10=cC(1,0); fC11=cC(1,1);
- fC20=cC(2,0); fC21=cC(2,1); fC22=cC(2,2);
- fC30=cC(3,0); fC31=cC(3,1); fC32=cC(3,2); fC33=cC(3,3);
- fC40=cC(4,0); fC41=cC(4,1); fC42=cC(4,2); fC43=cC(4,3); fC44=cC(4,4);