AnalyticalJacobians/src/JacobianCurvilinearToLocal.cc

0001 #include "TrackingTools/AnalyticalJacobians/interface/JacobianCurvilinearToLocal.h"
0002 #include "DataFormats/GeometrySurface/interface/Surface.h"
0003 #include "TrackingTools/TrajectoryParametrization/interface/LocalTrajectoryParameters.h"
0004 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0005
0006 #include "MagneticField/Engine/interface/MagneticField.h"
0007
0008 JacobianCurvilinearToLocal::JacobianCurvilinearToLocal(const Surface& surface,
0009                                                        const LocalTrajectoryParameters& localParameters,
0010                                                        const MagneticField& magField)
0011     : theJacobian(ROOT::Math::SMatrixNoInit()) {
0012   GlobalPoint x = surface.toGlobal(localParameters.position());
0013   GlobalVector h = magField.inInverseGeV(x);
0014   GlobalVector qh = h * localParameters.signedInverseMomentum();  // changed sign
0015
0016   LocalVector tnl = localParameters.direction();
0017   GlobalVector tn = surface.toGlobal(tnl);
0018   double t1r = 1. / tnl.z();
0019
0020   Surface::RotationType const& rot = surface.rotation();
0021
0022   compute(rot, tn, qh, t1r);
0023 }
0024
0025 JacobianCurvilinearToLocal::JacobianCurvilinearToLocal(const Surface& surface,
0026                                                        const LocalTrajectoryParameters& localParameters,
0027                                                        const GlobalTrajectoryParameters& globalParameters,
0028                                                        const MagneticField& magField)
0029     : theJacobian(ROOT::Math::SMatrixNoInit()) {
0030   // GlobalPoint  x =  globalParameters.position();
0031   // GlobalVector h  = magField.inInverseGeV(x);
0032   GlobalVector h = globalParameters.magneticFieldInInverseGeV();
0033   GlobalVector qh = h * localParameters.signedInverseMomentum();  // changed sign
0034
0035   LocalVector tnl = localParameters.direction();
0036   // GlobalVector tn = surface.toGlobal(tnl); // faster?
0037   GlobalVector tn = globalParameters.momentum() * std::abs(localParameters.signedInverseMomentum());
0038   double t1r = 1. / tnl.z();
0039
0040   Surface::RotationType const& rot = surface.rotation();
0041
0042   compute(rot, tn, qh, t1r);
0043 }
0044
0045 void JacobianCurvilinearToLocal::compute(Surface::RotationType const& rot,
0046                                          GlobalVector const& tn,
0047                                          GlobalVector const& qh,
0048                                          double t1r) {
0049   // Origin: TRSCSD
0050
0051   double cosl = tn.perp();
0052   if (cosl < 1.e-30)
0053     cosl = 1.e-30;
0054   double cosl1 = 1. / cosl;
0055   const GlobalVector un(-tn.y() * cosl1, tn.x() * cosl1, 0.);
0056   const GlobalVector vn(-tn.z() * un.y(), tn.z() * un.x(), cosl);
0057
0058   auto const u = rot.rotate(un.basicVector());
0059   auto const v = rot.rotate(vn.basicVector());
0060
0061   int j = 0, k = 1, i = 2;
0062
0063   double t2r = t1r * t1r;
0064   double t3r = t1r * t2r;
0065
0066   theJacobian(0, 0) = 1.;
0067   for (auto i = 1; i < 5; ++i)
0068     theJacobian(0, i) = 0.;
0069   theJacobian(1, 0) = 0.;
0070   theJacobian(2, 0) = 0.;
0071
0072   theJacobian(1, 1) = -u[k] * t2r;
0073   theJacobian(1, 2) = v[k] * (cosl * t2r);
0074   theJacobian(2, 1) = u[j] * t2r;
0075   theJacobian(2, 2) = -v[j] * (cosl * t2r);
0076
0077   for (auto i = 0; i < 3; ++i) {
0078     theJacobian(3, i) = 0.;
0079     theJacobian(4, i) = 0.;
0080   }
0081
0082   theJacobian(3, 3) = v[k] * t1r;
0083   theJacobian(3, 4) = -u[k] * t1r;
0084   theJacobian(4, 3) = -v[j] * t1r;
0085   theJacobian(4, 4) = u[j] * t1r;
0086
0087   double sinz = un.dot(qh);
0088   double cosz = -vn.dot(qh);
0089   double ui = u[i] * (t3r);
0090   double vi = v[i] * (t3r);
0091   theJacobian(1, 3) = -ui * (v[k] * cosz - u[k] * sinz);
0092   theJacobian(1, 4) = -vi * (v[k] * cosz - u[k] * sinz);
0093   theJacobian(2, 3) = ui * (v[j] * cosz - u[j] * sinz);
0094   theJacobian(2, 4) = vi * (v[j] * cosz - u[j] * sinz);
0095   // end of TRSCSD
0096   //dbg::dbg_trace(1,"Cu2L", localParameters.vector(),di,dj,dk,theJacobian);
0097 }