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