AnalyticalJacobians/test/Jacobian_t.cpp

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 }  // namespace
0024
0025 #include "FWCore/Utilities/interface/HRRealTime.h"
0026 void st() {}
0027 void en() {}
0028
0029 int main() {
0030   // GlobalVector xx(0.5,1.,1.);
0031   // GlobalVector yy(-1.,0.5,1.);
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   // verify cart to curv and back....
0093   // AlgebraicMatrixID();
0094
0095   // L ->Cart
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   // L -> Curv
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   // Cart -> Loc
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   // Curv -> Loc
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 }