AnalyticalJacobians/test/ErrorProp_t.cpp

0001 #include "TrackingTools/GeomPropagators/interface/HelixForwardPlaneCrossing.h"
0002 #include "TrackingTools/GeomPropagators/interface/HelixArbitraryPlaneCrossing.h"
0003 #include "TrackingTools/AnalyticalJacobians/interface/AnalyticalCurvilinearJacobian.h"
0004 #include "TrackingTools/AnalyticalJacobians/interface/JacobianLocalToCurvilinear.h"
0005 #include "TrackingTools/AnalyticalJacobians/interface/JacobianCurvilinearToLocal.h"
0006 #include "TrackingTools/TrajectoryParametrization/interface/GlobalTrajectoryParameters.h"
0007
0008 #include "DataFormats/GeometrySurface/interface/Surface.h"
0009 #include "TrackingTools/TrajectoryParametrization/interface/LocalTrajectoryParameters.h"
0010
0011 #include "DataFormats/GeometrySurface/interface/Plane.h"
0012
0013 #include "MagneticField/Engine/interface/MagneticField.h"
0014
0015 namespace {
0016
0017   struct M5T : public MagneticField {
0018     explicit M5T(double br) : m(br, br, 5.) {}
0019     virtual GlobalVector inTesla(const GlobalPoint&) const { return m; }
0020
0021     GlobalVector m;
0022   };
0023
0024 }  // namespace
0025
0026 #include "FWCore/Utilities/interface/HRRealTime.h"
0027 void st(void* ptr) {
0028   asm("" /* no instructions */
0029       :  /* no inputs */
0030       : /* output */ "rax"(ptr)
0031       : /* pretend to clobber */ "memory");
0032 }
0033 void en(void const* ptr) {
0034   asm("" /* no instructions */
0035       :  /* no inputs */
0036       : /* output */ "rax"(ptr)
0037       : /* pretend to clobber */ "memory");
0038 }
0039
0040 #include <iostream>
0041
0042 int main(int argc, char** argv) {
0043   double br = 0.;
0044   if (argc > 1)
0045     br = 0.1;
0046   M5T const m(br);
0047
0048   Basic3DVector<float> axis(0.5, 1., 1);
0049
0050   Surface::RotationType rot(axis, 0.5 * M_PI);
0051   std::cout << rot << std::endl;
0052
0053   Surface::PositionType pos(0., 0., 0.);
0054
0055   Plane plane(pos, rot);
0056   LocalTrajectoryParameters tpl(-1. / 3.5, 1., 1., 0., 0., 1.);
0057   GlobalVector mg = plane.toGlobal(tpl.momentum());
0058   GlobalTrajectoryParameters tpg(pos, mg, -1., &m);
0059   double curv = tpg.transverseCurvature();
0060   std::cout << curv << " " << mg.mag() << std::endl;
0061   std::cout << tpg.position() << " " << tpg.momentum() << std::endl;
0062
0063   AlgebraicMatrix55 fullJacobian = AlgebraicMatrixID();
0064   AlgebraicMatrix55 deltaJacobian = AlgebraicMatrixID();
0065   GlobalTrajectoryParameters tpg0(tpg);
0066
0067   //HelixForwardPlaneCrossing::PositionType zero(0.,0.,0.);
0068   GlobalPoint zero(0., 0., 0.);
0069   std::cout << std::endl;
0070
0071   {
0072     double totalStep(0.);
0073     double singleStep(1.);
0074     edm::HRTimeType timeFull = 0;
0075     edm::HRTimeType timeInf = 0;
0076     for (int i = 0; i < 10; ++i) {
0077       double s = singleStep;
0078       totalStep += singleStep;
0079
0080       GlobalVector h = tpg.magneticFieldInInverseGeV(tpg.position());
0081       Surface::RotationType rot(Basic3DVector<float>(h), 0);
0082       Plane lplane(zero, rot);
0083       HelixForwardPlaneCrossing::PositionType a(lplane.toLocal(tpg.position()));
0084       HelixForwardPlaneCrossing::DirectionType p(lplane.toLocal(tpg.momentum()));
0085       double lcurv = -h.mag() / p.perp() * tpg.charge();
0086       std::cout << lcurv << " " << p.mag() << std::endl;
0087       HelixForwardPlaneCrossing prop(a, p, lcurv);
0088       LocalPoint x(prop.position(s));
0089       LocalVector dir(prop.direction(s));
0090       std::cout << dir.mag() << std::endl;
0091
0092       GlobalTrajectoryParameters tpg2(
0093           lplane.toGlobal(x), (p.mag() / dir.mag()) * lplane.toGlobal(dir), tpg.charge(), &m);
0094
0095       std::cout << tpg2.position() << " " << tpg2.momentum() << std::endl;
0096       AnalyticalCurvilinearJacobian full;
0097       AnalyticalCurvilinearJacobian delta;
0098       edm::HRTimeType sf = edm::hrRealTime();
0099       st(&tpg2);
0100       full.computeFullJacobian(tpg, tpg2.position(), tpg2.momentum(), h, s);
0101       en(&full.jacobian());
0102       timeFull += (edm::hrRealTime() - sf);
0103       edm::HRTimeType si = edm::hrRealTime();
0104       st(&tpg2);
0105       delta.computeInfinitesimalJacobian(tpg, tpg2.position(), tpg2.momentum(), h, s);
0106       en(&delta.jacobian());
0107       timeInf += (edm::hrRealTime() - si);
0108       std::cout << full.jacobian() << std::endl;
0109       std::cout << std::endl;
0110       std::cout << delta.jacobian() << std::endl;
0111       std::cout << std::endl;
0112       tpg = tpg2;
0113       fullJacobian *= full.jacobian();
0114       deltaJacobian *= delta.jacobian();
0115     }
0116
0117     std::cout << "---------------------------" << std::endl;
0118     std::cout << std::endl;
0119     std::cout << std::endl;
0120     std::cout << "fullJacobian " << timeFull << std::endl;
0121     std::cout << fullJacobian << std::endl;
0122     std::cout << std::endl;
0123     std::cout << "deltaJacobian " << timeInf << std::endl;
0124     std::cout << deltaJacobian << std::endl;
0125     std::cout << std::endl;
0126
0127     AnalyticalCurvilinearJacobian full;
0128     GlobalVector h = tpg0.magneticFieldInInverseGeV(tpg0.position());
0129     edm::HRTimeType s = edm::hrRealTime();
0130     for (int k = 0; k < 50000; ++k) {
0131       st(&tpg0);
0132       full.computeFullJacobian(tpg0, tpg.position(), tpg.momentum(), h, totalStep);
0133       en(&full.jacobian());
0134     }
0135     edm::HRTimeType e = edm::hrRealTime();
0136     std::cout << "one step fullJacobian " << (e - s) / 50000. << " precise timing" << std::endl;
0137     std::cout << full.jacobian() << std::endl;
0138     std::cout << std::endl;
0139
0140     std::cout << "---------------------------" << std::endl;
0141
0142     AlgebraicMatrix55 div;
0143     for (unsigned int i = 0; i < 5; ++i) {
0144       for (unsigned int j = 0; j < 5; ++j) {
0145         //       div(i,j) = fabs(full.jacobian()(i,j))>1.e-20 ? fullJacobian(i,j)/full.jacobian()(i,j)-1 : 0;
0146         div(i, j) = fullJacobian(i, j) - full.jacobian()(i, j);
0147       }
0148     }
0149     std::cout << "Full relative" << std::endl << div << std::endl;
0150     for (unsigned int i = 0; i < 5; ++i) {
0151       for (unsigned int j = 0; j < 5; ++j) {
0152         //       div(i,j) = fabs(full.jacobian()(i,j))>1.e-20 ? deltaJacobian(i,j)/full.jacobian()(i,j)-1 : 0;
0153         div(i, j) = deltaJacobian(i, j) - full.jacobian()(i, j);
0154       }
0155     }
0156     std::cout << "Delta relative" << std::endl << div << std::endl;
0157   }
0158
0159   // for timing no printout
0160
0161   fullJacobian = AlgebraicMatrixID();
0162   deltaJacobian = AlgebraicMatrixID();
0163   tpg = tpg0;
0164   {
0165     double totalStep(0.);
0166     double singleStep(1.);
0167     edm::HRTimeType timeFull = 0;
0168     edm::HRTimeType timeInf = 0;
0169     for (int i = 0; i < 10; ++i) {
0170       double s = singleStep;
0171       totalStep += singleStep;
0172
0173       GlobalVector h = tpg.magneticFieldInInverseGeV(tpg.position());
0174       Surface::RotationType rot(Basic3DVector<float>(h), 0);
0175       Plane lplane(zero, rot);
0176       HelixForwardPlaneCrossing::PositionType a(lplane.toLocal(tpg.position()));
0177       HelixForwardPlaneCrossing::DirectionType p(lplane.toLocal(tpg.momentum()));
0178       double lcurv = -h.mag() / p.perp() * tpg.charge();
0179       HelixForwardPlaneCrossing prop(a, p, lcurv);
0180       LocalPoint x(prop.position(s));
0181       LocalVector dir(prop.direction(s));
0182
0183       GlobalTrajectoryParameters tpg2(
0184           lplane.toGlobal(x), (p.mag() / dir.mag()) * lplane.toGlobal(dir), tpg.charge(), &m);
0185
0186       AnalyticalCurvilinearJacobian full;
0187       AnalyticalCurvilinearJacobian delta;
0188       edm::HRTimeType sf = edm::hrRealTime();
0189       full.computeFullJacobian(tpg, tpg2.position(), tpg2.momentum(), h, s);
0190       timeFull += (edm::hrRealTime() - sf);
0191       edm::HRTimeType si = edm::hrRealTime();
0192       delta.computeInfinitesimalJacobian(tpg, tpg2.position(), tpg2.momentum(), h, s);
0193       timeInf += (edm::hrRealTime() - si);
0194       tpg = tpg2;
0195       fullJacobian *= full.jacobian();
0196       deltaJacobian *= delta.jacobian();
0197     }
0198
0199     AnalyticalCurvilinearJacobian full;
0200     GlobalVector h = tpg0.magneticFieldInInverseGeV(tpg0.position());
0201     edm::HRTimeType s = edm::hrRealTime();
0202     full.computeFullJacobian(tpg0, tpg.position(), tpg.momentum(), h, totalStep);
0203     edm::HRTimeType e = edm::hrRealTime();
0204
0205     std::cout << "---------------------------" << std::endl;
0206     std::cout << std::endl;
0207     std::cout << std::endl;
0208     std::cout << "fullJacobian " << timeFull << std::endl;
0209     std::cout << fullJacobian << std::endl;
0210     std::cout << std::endl;
0211     std::cout << "deltaJacobian " << timeInf << std::endl;
0212     std::cout << deltaJacobian << std::endl;
0213     std::cout << std::endl;
0214
0215     std::cout << "one step fullJacobian " << e - s << std::endl;
0216     std::cout << full.jacobian() << std::endl;
0217     std::cout << std::endl;
0218
0219     std::cout << "---------------------------" << std::endl;
0220
0221     // a long loop to record perf
0222     AlgebraicMatrix55 jjj = full.jacobian();
0223     if (argc > 2) {
0224       for (int kk = 0; kk < 100000; ++kk) {
0225         full.computeFullJacobian(tpg0, tpg.position(), tpg.momentum(), h, totalStep);
0226         jjj = full.jacobian();
0227       }
0228     }
0229   }
0230
0231   //
0232   // plane at end of propagation (equivalent to Curvilinear definition up to rotation in the plane)
0233   //
0234   Basic3DVector<float> newaxis(tpg.momentum().unit());
0235   Surface::RotationType newrot(axis, 0.);
0236   Surface::PositionType newpos(tpg.position());
0237   Plane newplane(newpos, newrot);
0238   //
0239   // total jacobian (local to local)
0240   //
0241   JacobianLocalToCurvilinear jlc(plane, tpl, m);
0242   LocalTrajectoryParameters newlpg(newplane.toLocal(tpg.position()), newplane.toLocal(tpg.momentum()), 1.);
0243   JacobianCurvilinearToLocal jcl(newplane, newlpg, m);
0244   //   AlgebraicMatrix55 jacobianL2L = jcl.jacobian()*full.jacobian()*jlc.jacobian();
0245   AlgebraicMatrix55 jacobianL2L = jcl.jacobian() * deltaJacobian * jlc.jacobian();
0246   //
0247   // redo propagation to target plane (should give identical position ...)
0248   //
0249   HelixArbitraryPlaneCrossing aprop(HelixForwardPlaneCrossing::PositionType(tpg0.position()),
0250                                     HelixForwardPlaneCrossing::DirectionType(tpg0.momentum()),
0251                                     curv);
0252   std::pair<bool, double> as = aprop.pathLength(newplane);
0253   std::cout << "as = " << as.second << std::endl;
0254   GlobalPoint ax(aprop.position(as.second));
0255   GlobalVector ap(aprop.direction(as.second) * tpg0.momentum().perp());
0256   //
0257   // change initial q/p and redo propagation
0258   //
0259   double qpmod = 1. / tpg0.momentum().mag() + 0.01;
0260   GlobalTrajectoryParameters tpg0mod(tpg0.position(), tpg0.momentum().unit() / fabs(qpmod), qpmod > 0 ? 1. : -1., &m);
0261   HelixArbitraryPlaneCrossing aprop2(HelixForwardPlaneCrossing::PositionType(tpg0mod.position()),
0262                                      HelixForwardPlaneCrossing::DirectionType(tpg0mod.momentum()),
0263                                      tpg0mod.transverseCurvature());
0264   std::pair<bool, double> as2 = aprop2.pathLength(newplane);
0265   std::cout << "as2 = " << as2.second << std::endl;
0266   GlobalPoint ax2(aprop2.position(as2.second));
0267   GlobalVector ap2(aprop2.direction(as2.second) * tpg0mod.momentum().perp());
0268   //
0269   // local coordinates after 2nd propagation
0270   LocalVector newDx = newplane.toLocal(ax2 - ax);
0271   LocalVector newDp = newplane.toLocal(ap2 / ap2.z() - ap / ap.z());
0272   std::cout << newDx << std::endl;
0273   std::cout << newDp << std::endl;
0274   //
0275   // prediction of variations with jacobian
0276   //
0277   AlgebraicVector5 dp;
0278   dp(0) = qpmod - 1. / tpg0.momentum().mag();
0279   dp(1) = dp(2) = dp(3) = dp(4) = 0;
0280   AlgebraicVector5 dpPred = jacobianL2L * dp;
0281   std::cout << dpPred << std::endl;
0282   //
0283   // comparison with difference in parameters
0284   //
0285   std::cout << newDx.x() - dpPred[3] << " " << newDx.y() - dpPred[4] << std::endl;
0286   std::cout << newDp.x() - dpPred[1] << " " << newDp.y() - dpPred[2] << std::endl;
0287   std::cout << (1. / ap2.mag() - 1. / ap.mag()) - dpPred[0] << std::endl;
0288   return 0;
0289 }