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(void* ptr) {
0027   asm("" /* no instructions */
0028       :  /* no inputs */
0029       : /* output */ "rax"(ptr)
0030       : /* pretend to clobber */ "memory");
0031 }
0032 void en(void const* ptr) {
0033   asm("" /* no instructions */
0034       :  /* no inputs */
0035       : /* output */ "rax"(ptr)
0036       : /* pretend to clobber */ "memory");
0037 }
0038
0039 int main() {
0040   // GlobalVector xx(0.5,1.,1.);
0041   // GlobalVector yy(-1.,0.5,1.);
0042
0043   Basic3DVector<float> axis(0.5, 1., 1);
0044
0045   Surface::RotationType rot(axis, 0.5 * M_PI);
0046   std::cout << rot << std::endl;
0047
0048   Surface::PositionType pos(0., 0., 0.);
0049
0050   Plane plane(pos, rot);
0051
0052   GlobalVector g1 = plane.toGlobal(LocalVector(1., 0., 0.));
0053   GlobalVector g2 = plane.toGlobal(LocalVector(0., 1., 0.));
0054   GlobalVector g3 = plane.toGlobal(LocalVector(0., 0., 1.));
0055   AlgebraicMatrix33 Rsub;
0056   Rsub(0, 0) = g1.x();
0057   Rsub(0, 1) = g2.x();
0058   Rsub(0, 2) = g3.x();
0059   Rsub(1, 0) = g1.y();
0060   Rsub(1, 1) = g2.y();
0061   Rsub(1, 2) = g3.y();
0062   Rsub(2, 0) = g1.z();
0063   Rsub(2, 1) = g2.z();
0064   Rsub(2, 2) = g3.z();
0065
0066   std::cout << Rsub << std::endl;
0067
0068   if (rot.xx() != Rsub(0, 0) || rot.xy() != Rsub(1, 0) || rot.xz() != Rsub(2, 0) || rot.yx() != Rsub(0, 1) ||
0069       rot.yy() != Rsub(1, 1) || rot.yz() != Rsub(2, 1) || rot.zx() != Rsub(0, 2) || rot.zy() != Rsub(1, 2) ||
0070       rot.zz() != Rsub(2, 2))
0071     std::cout << " wrong assumption!" << std::endl;
0072
0073   GlobalVector dj(rot.x());
0074   GlobalVector dk(rot.y());
0075   GlobalVector di(rot.z());
0076
0077   GlobalVector un(-1., 1.5, 0.5);
0078   double ui = un.dot(di);
0079   double uj = un.dot(dj);
0080   double uk = un.dot(dk);
0081   std::cout << '\n' << un << std::endl;
0082   std::cout << '\n' << uj << "," << uk << "," << ui << std::endl;
0083   std::cout << rot.rotate(un.basicVector()) << std::endl;
0084   std::cout << rot.rotateBack(un.basicVector()) << '\n' << std::endl;
0085
0086   M5T const m;
0087   LocalTrajectoryParameters tp(1., 1., 1., 0., 0., 1.);
0088   std::cout << tp.vector() << std::endl;
0089   std::cout << tp.charge() << " " << tp.signedInverseMomentum() << std::endl;
0090   GlobalVector tn = plane.toGlobal(tp.momentum()).unit();
0091   GlobalTrajectoryParameters gp(plane.toGlobal(tp.position()), plane.toGlobal(tp.momentum()), tp.charge(), &m);
0092   LocalVector tnl = plane.toLocal(tn);
0093   std::cout << tp.position() << std::endl;
0094   std::cout << tp.momentum() << std::endl;
0095   std::cout << tp.momentum().unit() << std::endl;
0096   std::cout << tp.direction() << std::endl;
0097   std::cout << tnl << std::endl;
0098   std::cout << tn << std::endl;
0099   std::cout << plane.toGlobal(tnl) << std::endl;
0100   std::cout << gp.direction() << std::endl;
0101
0102   // verify cart to curv and back....
0103   // AlgebraicMatrixID();
0104
0105   int N = 1000000;
0106
0107   // L =>Cart
0108   {
0109     std::cout << "L =>Cart ";
0110     edm::HRTimeType s = edm::hrRealTime();
0111     for (int i = 0; i < N; ++i) {
0112       st(&tp);
0113       JacobianLocalToCartesian __attribute__((aligned(16))) jl2c(plane, tp);
0114       en(&jl2c.jacobian());
0115     }
0116     edm::HRTimeType e = edm::hrRealTime();
0117     std::cout << double(e - s) / N << std::endl;
0118     JacobianLocalToCartesian __attribute__((aligned(16))) jl2c(plane, tp);
0119     std::cout << jl2c.jacobian() << std::endl;
0120   }
0121
0122   // L => Curv
0123   {
0124     std::cout << "L =>Curv from loc ";
0125     edm::HRTimeType s = edm::hrRealTime();
0126     for (int i = 0; i < N; ++i) {
0127       st(&tp);
0128       JacobianLocalToCurvilinear __attribute__((aligned(16))) jl2c(plane, tp, m);
0129       en(&jl2c.jacobian());
0130     }
0131     edm::HRTimeType e = edm::hrRealTime();
0132     std::cout << double(e - s) / N << std::endl;
0133     JacobianLocalToCurvilinear __attribute__((aligned(16))) jl2c(plane, tp, m);
0134     std::cout << jl2c.jacobian() << std::endl;
0135   }
0136
0137   {
0138     std::cout << "L =>Curv from loc+glob ";
0139     edm::HRTimeType s = edm::hrRealTime();
0140     for (int i = 0; i < N; ++i) {
0141       st(&tp);
0142       JacobianLocalToCurvilinear __attribute__((aligned(16))) jl2c(plane, tp, gp, m);
0143       en(&jl2c.jacobian());
0144     }
0145     edm::HRTimeType e = edm::hrRealTime();
0146     std::cout << double(e - s) / N << std::endl;
0147     JacobianLocalToCurvilinear __attribute__((aligned(16))) jl2c(plane, tp, gp, m);
0148     std::cout << jl2c.jacobian() << std::endl;
0149   }
0150
0151   // Cart => Loc
0152   {
0153     std::cout << "Cart => Loc ";
0154     edm::HRTimeType s = edm::hrRealTime();
0155     for (int i = 0; i < N; ++i) {
0156       st(&tp);
0157       JacobianCartesianToLocal __attribute__((aligned(16))) jl2c(plane, tp);
0158       en(&jl2c.jacobian());
0159     }
0160     edm::HRTimeType e = edm::hrRealTime();
0161     std::cout << double(e - s) / N << std::endl;
0162     JacobianCartesianToLocal __attribute__((aligned(16))) jl2c(plane, tp);
0163     std::cout << jl2c.jacobian() << std::endl;
0164   }
0165
0166   // Curv => Loc
0167   {
0168     std::cout << "Curv => Loc from loc ";
0169     edm::HRTimeType s = edm::hrRealTime();
0170     for (int i = 0; i < N; ++i) {
0171       st(&tp);
0172       JacobianCurvilinearToLocal __attribute__((aligned(16))) jl2c(plane, tp, m);
0173       en(&jl2c.jacobian());
0174     }
0175     edm::HRTimeType e = edm::hrRealTime();
0176     std::cout << double(e - s) / N << std::endl;
0177     JacobianCurvilinearToLocal __attribute__((aligned(16))) jl2c(plane, tp, m);
0178     std::cout << jl2c.jacobian() << std::endl;
0179   }
0180
0181   {
0182     std::cout << "Curv => Loc from loc + glob ";
0183     edm::HRTimeType s = edm::hrRealTime();
0184     for (int i = 0; i < N; ++i) {
0185       st(&tp);
0186       JacobianCurvilinearToLocal __attribute__((aligned(16))) jl2c(plane, tp, gp, m);
0187       en(&jl2c.jacobian());
0188     }
0189     edm::HRTimeType e = edm::hrRealTime();
0190     std::cout << double(e - s) / N << std::endl;
0191     JacobianCurvilinearToLocal __attribute__((aligned(16))) jl2c(plane, tp, gp, m);
0192     std::cout << jl2c.jacobian() << std::endl;
0193   }
0194
0195   return 0;
0196 }