GeometryVector/src/TrackingJacobians.cc

0001 #include "DataFormats/GeometryVector/interface/TrackingJacobians.h"
0002
0003 // Code moved from TrackingTools/AnalyticalJacobians
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   // for neutrals: qbp is 1/p instead of q/p -
0014   //   equivalent to charge 1
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   //ErrorPropagation:
0055   //    C(Cart) = R(6*6) * J(6*5) * C(Curvi) * J(5*6)_T * R(6*6)_T
0056   theJacobian = R * theJacobian;
0057   //dbg::dbg_trace(1,"Cu2Ca", globalParameters.vector(),theJacobian);
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   // for neutrals: qbp is 1/p instead of q/p -
0071   //   equivalent to charge 1
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     //theJacobian(1,3) = (px*pz)/(pt*p2); theJacobian(1,4) = (py*pz)/(pt*p2); theJacobian(1,5) = -pt/p2; //wrong sign
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   //dbg::dbg_trace(1,"Ca2Cu", globalParameters.vector(),theJacobian);
0110   return theJacobian;
0111 }