File indexing completed on 2023-10-25 10:05:32
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() {}
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
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
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
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
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
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
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
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
0227 AlgebraicMatrix55 jacobianL2L = jcl.jacobian() * deltaJacobian * jlc.jacobian();
0228
0229
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
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
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
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
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 }