AnalyticalJacobians/src/JacobianCartesianToLocal.cc

0001 #include "TrackingTools/AnalyticalJacobians/interface/JacobianCartesianToLocal.h"
0002 #include "DataFormats/GeometrySurface/interface/Surface.h"
0003 #include "TrackingTools/TrajectoryParametrization/interface/LocalTrajectoryParameters.h"
0004
0005 JacobianCartesianToLocal::JacobianCartesianToLocal(const Surface& surface,
0006                                                    const LocalTrajectoryParameters& localParameters)
0007     : theJacobian() {
0008   //LocalCoordinates 1 = (q/|p|, dx/dz, dy/dz, x, y)
0009   //LocalCoordinates 2 = (x, y, z, px, py, pz)
0010   //Transformation: q/|p| = q/sqrt(px*px + py*py + pz*pz)
0011   //                dx/dz = px/pz
0012   //                dy/dz = py/pz
0013   //                x     = x
0014   //                y     = y
0015   LocalVector plocal = localParameters.momentum();
0016   double px = plocal.x(), py = plocal.y(), pz = plocal.z();
0017   double p = plocal.mag();
0018   TrackCharge q = localParameters.charge();
0019   // for neutrals: qbp is 1/p instead of q/p -
0020   //   equivalent to charge 1
0021   if (q == 0)
0022     q = 1;
0023   //Jacobian theJacobian( (q/|p|, dx/dz, dy/dz, x, y) = f(x, y, z, px, py, pz) )
0024   theJacobian(0, 3) = -q * px / (p * p * p);
0025   theJacobian(0, 4) = -q * py / (p * p * p);
0026   theJacobian(0, 5) = -q * pz / (p * p * p);
0027   if (fabs(pz) > 0) {
0028     theJacobian(1, 3) = 1. / pz;
0029     theJacobian(1, 5) = -px / (pz * pz);
0030     theJacobian(2, 4) = 1. / pz;
0031     theJacobian(2, 5) = -py / (pz * pz);
0032   }
0033   theJacobian(3, 0) = 1.;
0034   theJacobian(4, 1) = 1.;
0035
0036   /*
0037   LocalVector l1 = surface.toLocal(GlobalVector(1., 0., 0.));
0038   LocalVector l2 = surface.toLocal(GlobalVector(0., 1., 0.));
0039   LocalVector l3 = surface.toLocal(GlobalVector(0., 0., 1.));
0040   AlgebraicMatrix33 Rsub;
0041   Rsub(0,0) = l1.x(); Rsub(0,1) = l2.x(); Rsub(0,2) = l3.x();
0042   Rsub(1,0) = l1.y(); Rsub(1,1) = l2.y(); Rsub(1,2) = l3.y();
0043   Rsub(2,0) = l1.z(); Rsub(2,1) = l2.z(); Rsub(2,2) = l3.z();
0044   */
0045
0046   AlgebraicMatrix33 Rsub;
0047   // need to be copied anhhow to go from float to double...
0048   Surface::RotationType const& rot = surface.rotation();
0049   Rsub(0, 0) = rot.xx();
0050   Rsub(0, 1) = rot.xy();
0051   Rsub(0, 2) = rot.xz();
0052   Rsub(1, 0) = rot.yx();
0053   Rsub(1, 1) = rot.yy();
0054   Rsub(1, 2) = rot.yz();
0055   Rsub(2, 0) = rot.zx();
0056   Rsub(2, 1) = rot.zy();
0057   Rsub(2, 2) = rot.zz();
0058
0059   AlgebraicMatrix66 R;
0060   R.Place_at(Rsub, 0, 0);
0061   R.Place_at(Rsub, 3, 3);
0062   theJacobian = theJacobian * R;
0063   //dbg::dbg_trace(1,"Ca2L", localParameters.vector(),l1,l2,l3,theJacobian);
0064 }