File indexing completed on 2023-03-17 10:50:00
0001 #include "DataFormats/GeometryVector/interface/TrackingJacobians.h"
0002
0003
0004
0005 AlgebraicMatrix65 jacobianCurvilinearToCartesian(const GlobalVector& momentum, int q) {
0006 AlgebraicMatrix65 theJacobian;
0007 GlobalVector xt = momentum;
0008 GlobalVector yt(-xt.y(), xt.x(), 0.);
0009 GlobalVector zt = xt.cross(yt);
0010
0011 const GlobalVector& pvec = momentum;
0012 double pt = pvec.perp();
0013
0014
0015 if (q == 0)
0016 q = 1;
0017
0018 xt = xt.unit();
0019 if (fabs(pt) > 0) {
0020 yt = yt.unit();
0021 zt = zt.unit();
0022 }
0023
0024 AlgebraicMatrix66 R;
0025 R(0, 0) = xt.x();
0026 R(0, 1) = yt.x();
0027 R(0, 2) = zt.x();
0028 R(1, 0) = xt.y();
0029 R(1, 1) = yt.y();
0030 R(1, 2) = zt.y();
0031 R(2, 0) = xt.z();
0032 R(2, 1) = yt.z();
0033 R(2, 2) = zt.z();
0034 R(3, 3) = 1.;
0035 R(4, 4) = 1.;
0036 R(5, 5) = 1.;
0037
0038 double p = pvec.mag(), p2 = p * p;
0039 double sinlambda = pvec.z() / p, coslambda = pt / p;
0040 double sinphi = pvec.y() / pt, cosphi = pvec.x() / pt;
0041
0042 theJacobian(1, 3) = 1.;
0043 theJacobian(2, 4) = 1.;
0044 theJacobian(3, 0) = -q * p2 * coslambda * cosphi;
0045 theJacobian(3, 1) = -p * sinlambda * cosphi;
0046 theJacobian(3, 2) = -p * coslambda * sinphi;
0047 theJacobian(4, 0) = -q * p2 * coslambda * sinphi;
0048 theJacobian(4, 1) = -p * sinlambda * sinphi;
0049 theJacobian(4, 2) = p * coslambda * cosphi;
0050 theJacobian(5, 0) = -q * p2 * sinlambda;
0051 theJacobian(5, 1) = p * coslambda;
0052 theJacobian(5, 2) = 0.;
0053
0054
0055
0056 theJacobian = R * theJacobian;
0057
0058 return theJacobian;
0059 }
0060
0061 AlgebraicMatrix56 jacobianCartesianToCurvilinear(const GlobalVector& momentum, int q) {
0062 AlgebraicMatrix56 theJacobian;
0063 GlobalVector xt = momentum;
0064 GlobalVector yt(-xt.y(), xt.x(), 0.);
0065 GlobalVector zt = xt.cross(yt);
0066 const GlobalVector& pvec = momentum;
0067 double pt = pvec.perp(), p = pvec.mag();
0068 double px = pvec.x(), py = pvec.y(), pz = pvec.z();
0069 double pt2 = pt * pt, p2 = p * p, p3 = p * p * p;
0070
0071
0072 if (q == 0)
0073 q = 1;
0074 xt = xt.unit();
0075 if (fabs(pt) > 0) {
0076 yt = yt.unit();
0077 zt = zt.unit();
0078 }
0079
0080 AlgebraicMatrix66 R;
0081 R(0, 0) = xt.x();
0082 R(0, 1) = xt.y();
0083 R(0, 2) = xt.z();
0084 R(1, 0) = yt.x();
0085 R(1, 1) = yt.y();
0086 R(1, 2) = yt.z();
0087 R(2, 0) = zt.x();
0088 R(2, 1) = zt.y();
0089 R(2, 2) = zt.z();
0090 R(3, 3) = 1.;
0091 R(4, 4) = 1.;
0092 R(5, 5) = 1.;
0093
0094 theJacobian(0, 3) = -q * px / p3;
0095 theJacobian(0, 4) = -q * py / p3;
0096 theJacobian(0, 5) = -q * pz / p3;
0097 if (fabs(pt) > 0) {
0098
0099 theJacobian(1, 3) = -(px * pz) / (pt * p2);
0100 theJacobian(1, 4) = -(py * pz) / (pt * p2);
0101 theJacobian(1, 5) = pt / p2;
0102 theJacobian(2, 3) = -py / pt2;
0103 theJacobian(2, 4) = px / pt2;
0104 theJacobian(2, 5) = 0.;
0105 }
0106 theJacobian(3, 1) = 1.;
0107 theJacobian(4, 2) = 1.;
0108 theJacobian = theJacobian * R;
0109
0110 return theJacobian;
0111 }