File indexing completed on 2021-02-14 14:31:31
0001 #include "TrackingTools/AnalyticalJacobians/interface/JacobianLocalToCartesian.h"
0002 #include "TrackingTools/AnalyticalJacobians/interface/JacobianCartesianToLocal.h"
0003 #include "TrackingTools/AnalyticalJacobians/interface/JacobianLocalToCurvilinear.h"
0004 #include "TrackingTools/AnalyticalJacobians/interface/JacobianCurvilinearToLocal.h"
0005
0006 #include "DataFormats/GeometrySurface/interface/Surface.h"
0007 #include "DataFormats/TrajectoryState/interface/LocalTrajectoryParameters.h"
0008 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0009
0010 #include "DataFormats/GeometrySurface/interface/Plane.h"
0011
0012 #include "MagneticField/Engine/interface/MagneticField.h"
0013
0014 namespace {
0015
0016 struct M5T : public MagneticField {
0017 M5T() : m(0., 0., 5.) {}
0018 virtual GlobalVector inTesla(const GlobalPoint&) const { return m; }
0019
0020 GlobalVector m;
0021 };
0022
0023 }
0024
0025 #include "FWCore/Utilities/interface/HRRealTime.h"
0026 void st() {}
0027 void en() {}
0028
0029 int main() {
0030
0031
0032
0033 Basic3DVector<float> axis(0.5, 1., 1);
0034
0035 Surface::RotationType rot(axis, 0.5 * M_PI);
0036 std::cout << rot << std::endl;
0037
0038 Surface::PositionType pos(0., 0., 0.);
0039
0040 Plane plane(pos, rot);
0041
0042 GlobalVector g1 = plane.toGlobal(LocalVector(1., 0., 0.));
0043 GlobalVector g2 = plane.toGlobal(LocalVector(0., 1., 0.));
0044 GlobalVector g3 = plane.toGlobal(LocalVector(0., 0., 1.));
0045 AlgebraicMatrix33 Rsub;
0046 Rsub(0, 0) = g1.x();
0047 Rsub(0, 1) = g2.x();
0048 Rsub(0, 2) = g3.x();
0049 Rsub(1, 0) = g1.y();
0050 Rsub(1, 1) = g2.y();
0051 Rsub(1, 2) = g3.y();
0052 Rsub(2, 0) = g1.z();
0053 Rsub(2, 1) = g2.z();
0054 Rsub(2, 2) = g3.z();
0055
0056 std::cout << Rsub << std::endl;
0057
0058 if (rot.xx() != Rsub(0, 0) || rot.xy() != Rsub(1, 0) || rot.xz() != Rsub(2, 0) || rot.yx() != Rsub(0, 1) ||
0059 rot.yy() != Rsub(1, 1) || rot.yz() != Rsub(2, 1) || rot.zx() != Rsub(0, 2) || rot.zy() != Rsub(1, 2) ||
0060 rot.zz() != Rsub(2, 2))
0061 std::cout << " wrong assumption!" << std::endl;
0062
0063 GlobalVector dj(rot.x());
0064 GlobalVector dk(rot.y());
0065 GlobalVector di(rot.z());
0066
0067 GlobalVector un(-1., 1.5, 0.5);
0068 double ui = un.dot(di);
0069 double uj = un.dot(dj);
0070 double uk = un.dot(dk);
0071 std::cout << '\n' << un << std::endl;
0072 std::cout << '\n' << uj << "," << uk << "," << ui << std::endl;
0073 std::cout << rot.rotate(un.basicVector()) << std::endl;
0074 std::cout << rot.rotateBack(un.basicVector()) << '\n' << std::endl;
0075
0076 M5T const m;
0077 LocalTrajectoryParameters tp(1., 1., 1., 0., 0., 1.);
0078 std::cout << tp.vector() << std::endl;
0079 std::cout << tp.charge() << " " << tp.signedInverseMomentum() << std::endl;
0080 GlobalVector tn = plane.toGlobal(tp.momentum()).unit();
0081 GlobalTrajectoryParameters gp(plane.toGlobal(tp.position()), plane.toGlobal(tp.momentum()), tp.charge(), &m);
0082 LocalVector tnl = plane.toLocal(tn);
0083 std::cout << tp.position() << std::endl;
0084 std::cout << tp.momentum() << std::endl;
0085 std::cout << tp.momentum().unit() << std::endl;
0086 std::cout << tp.direction() << std::endl;
0087 std::cout << tnl << std::endl;
0088 std::cout << tn << std::endl;
0089 std::cout << plane.toGlobal(tnl) << std::endl;
0090 std::cout << gp.direction() << std::endl;
0091
0092
0093
0094
0095
0096 {
0097 std::cout << "L ->Cart" << std::endl;
0098 edm::HRTimeType s = edm::hrRealTime();
0099 st();
0100 JacobianLocalToCartesian __attribute__((aligned(16))) jl2c(plane, tp);
0101 en();
0102 edm::HRTimeType e = edm::hrRealTime();
0103 std::cout << e - s << std::endl;
0104 std::cout << jl2c.jacobian() << std::endl;
0105 }
0106
0107
0108 {
0109 std::cout << "L ->Curv from loc" << std::endl;
0110 edm::HRTimeType s = edm::hrRealTime();
0111 st();
0112 JacobianLocalToCurvilinear __attribute__((aligned(16))) jl2c(plane, tp, m);
0113 en();
0114 edm::HRTimeType e = edm::hrRealTime();
0115 std::cout << e - s << std::endl;
0116 std::cout << jl2c.jacobian() << std::endl;
0117 }
0118
0119 {
0120 std::cout << "L ->Curv from loc+glob" << std::endl;
0121 edm::HRTimeType s = edm::hrRealTime();
0122 st();
0123 JacobianLocalToCurvilinear __attribute__((aligned(16))) jl2c(plane, tp, gp, m);
0124 en();
0125 edm::HRTimeType e = edm::hrRealTime();
0126 std::cout << e - s << std::endl;
0127 std::cout << jl2c.jacobian() << std::endl;
0128 }
0129
0130
0131 {
0132 std::cout << "Cart -> Loc" << std::endl;
0133 edm::HRTimeType s = edm::hrRealTime();
0134 st();
0135 JacobianCartesianToLocal __attribute__((aligned(16))) jl2c(plane, tp);
0136 en();
0137 edm::HRTimeType e = edm::hrRealTime();
0138 std::cout << e - s << std::endl;
0139 std::cout << jl2c.jacobian() << std::endl;
0140 }
0141
0142
0143 {
0144 std::cout << "Curv -> Loc from loc" << std::endl;
0145 edm::HRTimeType s = edm::hrRealTime();
0146 st();
0147 JacobianCurvilinearToLocal __attribute__((aligned(16))) jl2c(plane, tp, m);
0148 en();
0149 edm::HRTimeType e = edm::hrRealTime();
0150 std::cout << e - s << std::endl;
0151 std::cout << jl2c.jacobian() << std::endl;
0152 }
0153
0154 {
0155 std::cout << "Curv -> Loc from loc + glob" << std::endl;
0156 edm::HRTimeType s = edm::hrRealTime();
0157 st();
0158 JacobianCurvilinearToLocal __attribute__((aligned(16))) jl2c(plane, tp, gp, m);
0159 en();
0160 edm::HRTimeType e = edm::hrRealTime();
0161 std::cout << e - s << std::endl;
0162 std::cout << jl2c.jacobian() << std::endl;
0163 }
0164
0165 return 0;
0166 }