File indexing completed on 2024-05-10 02:21:30
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 }
0025
0026 #include "FWCore/Utilities/interface/HRRealTime.h"
0027 void st(void* ptr) {
0028 asm(""
0029 :
0030 : "rax"(ptr)
0031 : "memory");
0032 }
0033 void en(void const* ptr) {
0034 asm(""
0035 :
0036 : "rax"(ptr)
0037 : "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
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
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
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
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
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
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
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
0245 AlgebraicMatrix55 jacobianL2L = jcl.jacobian() * deltaJacobian * jlc.jacobian();
0246
0247
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
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
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
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
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 }