AnalyticalJacobians/src/JacobianLocalToCurvilinear.cc

0001 #include "TrackingTools/AnalyticalJacobians/interface/JacobianLocalToCurvilinear.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 JacobianLocalToCurvilinear::JacobianLocalToCurvilinear(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 hq = h * localParameters.signedInverseMomentum();  // changed sign
0015
0016   LocalVector tnl = localParameters.direction();
0017   GlobalVector tn = surface.toGlobal(tnl);
0018
0019   Surface::RotationType const& rot = surface.rotation();
0020
0021   compute(rot, tnl, tn, hq);
0022 }
0023
0024 JacobianLocalToCurvilinear::JacobianLocalToCurvilinear(const Surface& surface,
0025                                                        const LocalTrajectoryParameters& localParameters,
0026                                                        const GlobalTrajectoryParameters& globalParameters,
0027                                                        const MagneticField& magField)
0028     : theJacobian(ROOT::Math::SMatrixNoInit()) {
0029   GlobalVector h = globalParameters.magneticFieldInInverseGeV();
0030   GlobalVector hq = h * localParameters.signedInverseMomentum();  // changed sign
0031
0032   LocalVector tnl = localParameters.direction();
0033   GlobalVector tn = surface.toGlobal(tnl);  // globalParameters.momentum().unit();
0034
0035   Surface::RotationType const& rot = surface.rotation();
0036
0037   compute(rot, tnl, tn, hq);
0038 }
0039
0040 void JacobianLocalToCurvilinear::compute(Surface::RotationType const& rot,
0041                                          LocalVector const& tnl,
0042                                          GlobalVector const& tn,
0043                                          GlobalVector const& hq) {
0044   // Origin: TRSDSC
0045
0046   GlobalVector dj(rot.x());
0047   GlobalVector dk(rot.y());
0048
0049   // rotate coordinates because of wrong coordinate system in orca
0050   double tvwX = tnl.z(), tvwY = tnl.x(), tvwZ = tnl.y();
0051   double cosl = tn.perp();
0052   if (cosl < 1.e-30)
0053     cosl = 1.e-30;
0054   double cosl1 = 1. / cosl;
0055
0056   GlobalVector un(-tn.y() * cosl1, tn.x() * cosl1, 0.);
0057   double uj = un.dot(dj);
0058   double uk = un.dot(dk);
0059   double sinz = -un.dot(hq);
0060
0061   GlobalVector vn(-tn.z() * un.y(), tn.z() * un.x(), cosl);
0062   double vj = vn.dot(dj);
0063   double vk = vn.dot(dk);
0064   double cosz = vn.dot(hq);
0065
0066   theJacobian(0, 0) = 1.;
0067   for (auto i = 1; i < 5; ++i)
0068     theJacobian(0, i) = 0.;
0069
0070   theJacobian(1, 0) = 0.;
0071   theJacobian(2, 0) = 0.;
0072
0073   theJacobian(1, 1) = tvwX * vj;
0074   theJacobian(1, 2) = tvwX * vk;
0075   theJacobian(2, 1) = tvwX * uj * cosl1;
0076   theJacobian(2, 2) = tvwX * uk * cosl1;
0077
0078   for (auto i = 0; i < 3; ++i) {
0079     theJacobian(3, i) = 0.;
0080     theJacobian(4, i) = 0.;
0081   }
0082
0083   theJacobian(3, 3) = uj;
0084   theJacobian(3, 4) = uk;
0085   theJacobian(4, 3) = vj;
0086   theJacobian(4, 4) = vk;
0087
0088   theJacobian(1, 3) = tvwY * sinz;
0089   theJacobian(1, 4) = tvwZ * sinz;
0090   theJacobian(2, 3) = tvwY * (cosz * cosl1);
0091   theJacobian(2, 4) = tvwZ * (cosz * cosl1);
0092   // end of TRSDSC
0093
0094   //dbg::dbg_trace(1,"Loc2Cu", localParameters.vector(),x,dj,dk,theJacobian);
0095 }