AnalyticalJacobians/src/JacobianLocalToCartesian.cc

0001 #include "TrackingTools/AnalyticalJacobians/interface/JacobianLocalToCartesian.h"
0002 #include "DataFormats/GeometrySurface/interface/Surface.h"
0003 #include "DataFormats/TrajectoryState/interface/LocalTrajectoryParameters.h"
0004
0005 JacobianLocalToCartesian::JacobianLocalToCartesian(const Surface& surface,
0006                                                    const LocalTrajectoryParameters& localParameters)
0007     : theJacobian() {
0008   //LocalCoordinates 1 = (x, y, z, px, py, pz)
0009   //LocalCoordinates 2 = (q/|p|, dx/dz, dy/dz, x, y)
0010   //Transformation: x  = x
0011   //                y  = y
0012   //          px = (q/(q/|p|)) * (dx/dz) *  sqrt(1./(1.+(dx/dz)^2+(dy/dz)^2))
0013   //          py = (q/(q/|p|)) * (dy/dz) * sqrt(1./(1.+(dx/dz)^2+(dy/dz)^2))
0014   //          pz = (q/(q/|p|)) * sqrt(1./(1.+(dx/dz)^2+(dy/dz)^2))
0015   //Jacobian J((x, y, px, py, pz)  = f(q/|p|, dx/dz, dy/dz, x, y))
0016
0017   // AlgebraicVector5 localTrackParams = localParameters.mixedFormatVector();
0018   double qbp = localParameters.qbp();
0019   double dxdz = localParameters.dxdz();
0020   double dydz = localParameters.dydz();
0021   TrackCharge iq = localParameters.charge();
0022   // for neutrals: qbp is 1/p instead of q/p -
0023   //   equivalent to charge 1
0024   if (iq == 0)
0025     iq = 1;
0026   double pzSign = localParameters.pzSign();
0027   double q = iq * pzSign;
0028   double sqr = sqrt(dxdz * dxdz + dydz * dydz + 1);
0029   double den = -q / (sqr * sqr * sqr * qbp);
0030
0031   // no difference between local and data member
0032   AlgebraicMatrix65& lJacobian = theJacobian;
0033   lJacobian(0, 3) = 1.;
0034   lJacobian(1, 4) = 1.;
0035   lJacobian(3, 0) = (dxdz * (-q / (sqr * qbp * qbp)));
0036   lJacobian(3, 1) = (q / (sqr * qbp) + (den * dxdz * dxdz));
0037   lJacobian(3, 2) = ((den * dxdz * dydz));
0038   lJacobian(4, 0) = (dydz * (-q / (sqr * qbp * qbp)));
0039   lJacobian(4, 1) = ((den * dxdz * dydz));
0040   lJacobian(4, 2) = (q / (sqr * qbp) + (den * dydz * dydz));
0041   lJacobian(5, 0) = (-q / (sqr * qbp * qbp));
0042   lJacobian(5, 1) = ((den * dxdz));
0043   lJacobian(5, 2) = ((den * dydz));
0044
0045   /*
0046   GlobalVector g1 = surface.toGlobal(LocalVector(1., 0., 0.));
0047   GlobalVector g2 = surface.toGlobal(LocalVector(0., 1., 0.));
0048   GlobalVector g3 = surface.toGlobal(LocalVector(0., 0., 1.));
0049   */
0050   AlgebraicMatrix33 Rsub;
0051   /*
0052   Rsub(0,0) = g1.x(); Rsub(0,1) = g2.x(); Rsub(0,2) = g3.x();
0053   Rsub(1,0) = g1.y(); Rsub(1,1) = g2.y(); Rsub(1,2) = g3.y();
0054   Rsub(2,0) = g1.z(); Rsub(2,1) = g2.z(); Rsub(2,2) = g3.z();
0055   */
0056   // need to be copied anhhow to go from float to double...
0057   Surface::RotationType const& rot = surface.rotation();
0058   Rsub(0, 0) = rot.xx();
0059   Rsub(0, 1) = rot.yx();
0060   Rsub(0, 2) = rot.zx();
0061   Rsub(1, 0) = rot.xy();
0062   Rsub(1, 1) = rot.yy();
0063   Rsub(1, 2) = rot.zy();
0064   Rsub(2, 0) = rot.xz();
0065   Rsub(2, 1) = rot.yz();
0066   Rsub(2, 2) = rot.zz();
0067
0068   AlgebraicMatrix66 R;
0069   R.Place_at(Rsub, 0, 0);
0070   R.Place_at(Rsub, 3, 3);
0071   theJacobian = R * lJacobian;
0072   //dbg::dbg_trace(1,"Loc2Ca", localParameters.vector(),g1,g2,g3,theJacobian);
0073 }