File indexing completed on 2024-04-06 12:31:35
0001 #include "TrackingTools/PatternTools/interface/TransverseImpactPointExtrapolator.h"
0002 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0003 #include "DataFormats/GeometrySurface/interface/Surface.h"
0004
0005 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0006 #include "DataFormats/TrajectorySeed/interface/PropagationDirection.h"
0007 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0008 #include "DataFormats/GeometryVector/interface/Point2DBase.h"
0009 #include "DataFormats/GeometryVector/interface/Vector2DBase.h"
0010 #include "DataFormats/GeometrySurface/interface/PlaneBuilder.h"
0011
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013
0014 TransverseImpactPointExtrapolator::TransverseImpactPointExtrapolator() : thePropagator(nullptr) {}
0015
0016 TransverseImpactPointExtrapolator::TransverseImpactPointExtrapolator(const MagneticField* field)
0017 : thePropagator(new AnalyticalPropagator(field, anyDirection)) {}
0018
0019 TransverseImpactPointExtrapolator::TransverseImpactPointExtrapolator(const Propagator& u) : thePropagator(u.clone()) {
0020 thePropagator->setPropagationDirection(anyDirection);
0021 }
0022
0023 TrajectoryStateOnSurface TransverseImpactPointExtrapolator::extrapolate(const FreeTrajectoryState& fts,
0024 const GlobalPoint& vtx) const {
0025 return doExtrapolation(fts, vtx, *thePropagator);
0026 }
0027
0028 TrajectoryStateOnSurface TransverseImpactPointExtrapolator::extrapolate(const TrajectoryStateOnSurface tsos,
0029 const GlobalPoint& vtx) const {
0030 if (!tsos.isValid())
0031 return tsos;
0032 return doExtrapolation(tsos, vtx, *thePropagator);
0033 }
0034
0035 TrajectoryStateOnSurface TransverseImpactPointExtrapolator::extrapolate(const FreeTrajectoryState& fts,
0036 const GlobalPoint& vtx,
0037 const Propagator& p) const {
0038
0039
0040
0041 std::unique_ptr<Propagator> p_cloned = SetPropagationDirection(p, anyDirection);
0042 return doExtrapolation(fts, vtx, *(p_cloned.get()));
0043 }
0044
0045 TrajectoryStateOnSurface TransverseImpactPointExtrapolator::extrapolate(const TrajectoryStateOnSurface tsos,
0046 const GlobalPoint& vtx,
0047 const Propagator& p) const {
0048 if (!tsos.isValid())
0049 return tsos;
0050
0051
0052
0053 std::unique_ptr<Propagator> p_cloned = SetPropagationDirection(p, anyDirection);
0054 return doExtrapolation(tsos, vtx, *(p_cloned.get()));
0055 }
0056
0057 TrajectoryStateOnSurface TransverseImpactPointExtrapolator::doExtrapolation(const TrajectoryStateOnSurface tsos,
0058 const GlobalPoint& vtx,
0059 const Propagator& p) const {
0060
0061
0062
0063 if (fabs(tsos.freeState()->transverseCurvature()) < 1.E-6) {
0064 LogDebug("TransverseImpactPointExtrapolator") << "negligeable curvature: using a trick to extrapolate:\n" << tsos;
0065
0066
0067
0068 GlobalVector xLocal = GlobalVector(-tsos.globalMomentum().y(), tsos.globalMomentum().x(), 0).unit();
0069
0070 GlobalVector yLocal(0., 0., 1.);
0071
0072 GlobalVector zLocal(xLocal.cross(yLocal));
0073
0074 const Surface::PositionType& origin(vtx);
0075 Surface::RotationType rotation(xLocal, yLocal, zLocal);
0076 ReferenceCountingPointer<Plane> surface = PlaneBuilder().plane(origin, rotation);
0077
0078 return p.propagate(*tsos.freeState(), *surface);
0079 } else {
0080 ReferenceCountingPointer<Plane> surface =
0081 tipSurface(tsos.globalPosition(), tsos.globalMomentum(), 1. / tsos.transverseCurvature(), vtx);
0082
0083
0084
0085 return p.propagate(tsos, *surface);
0086 }
0087 }
0088
0089 TrajectoryStateOnSurface TransverseImpactPointExtrapolator::doExtrapolation(const FreeTrajectoryState& fts,
0090 const GlobalPoint& vtx,
0091 const Propagator& p) const {
0092
0093
0094
0095 if (fabs(fts.transverseCurvature()) < 1.E-6) {
0096 LogDebug("TransverseImpactPointExtrapolator") << "negligeable curvature: using a trick to extrapolate:\n" << fts;
0097
0098
0099
0100 GlobalVector xLocal = GlobalVector(-fts.momentum().y(), fts.momentum().x(), 0).unit();
0101
0102 GlobalVector yLocal(0., 0., 1.);
0103
0104 GlobalVector zLocal(xLocal.cross(yLocal));
0105
0106 const Surface::PositionType& origin(vtx);
0107 Surface::RotationType rotation(xLocal, yLocal, zLocal);
0108 ReferenceCountingPointer<Plane> surface = PlaneBuilder().plane(origin, rotation);
0109
0110 return p.propagate(fts, *surface);
0111 } else {
0112 ReferenceCountingPointer<Plane> surface =
0113 tipSurface(fts.position(), fts.momentum(), 1. / fts.transverseCurvature(), vtx);
0114
0115
0116
0117 return p.propagate(fts, *surface);
0118 }
0119 }
0120
0121 ReferenceCountingPointer<Plane> TransverseImpactPointExtrapolator::tipSurface(const GlobalPoint& position,
0122 const GlobalVector& momentum,
0123 const double& signedTransverseRadius,
0124 const GlobalPoint& vertex) const {
0125 LogDebug("TransverseImpactPointExtrapolator") << position << "\n"
0126 << momentum << "\n"
0127 << "signedTransverseRadius : " << signedTransverseRadius << "\n"
0128 << vertex;
0129
0130 typedef Point2DBase<double, GlobalTag> PositionType2D;
0131 typedef Vector2DBase<double, GlobalTag> DirectionType2D;
0132
0133 PositionType2D x0(position.x(), position.y());
0134 DirectionType2D t0(-momentum.y(), momentum.x());
0135 t0 = t0 / t0.mag();
0136
0137 PositionType2D xc(x0 + signedTransverseRadius * t0);
0138
0139 DirectionType2D vtxDirection(xc.x() - vertex.x(), xc.y() - vertex.y());
0140 double vtxDistance = vtxDirection.mag();
0141
0142 const Surface::PositionType& origin(vertex);
0143 GlobalVector xLocal(vtxDirection.x() / vtxDistance, vtxDirection.y() / vtxDistance, 0.);
0144 if (vtxDistance < fabs(signedTransverseRadius)) {
0145 LogDebug("TransverseImpactPointExtrapolator") << "Inverting the x axis.";
0146 xLocal = -xLocal;
0147 }
0148 GlobalVector yLocal(0., 0., 1.);
0149 GlobalVector zLocal(xLocal.cross(yLocal));
0150 if (zLocal.dot(momentum) < 0.) {
0151 LogDebug("TransverseImpactPointExtrapolator") << "Inverting the y,z frame.";
0152 yLocal = -yLocal;
0153 zLocal = -zLocal;
0154 }
0155 Surface::RotationType rotation(xLocal, yLocal, zLocal);
0156
0157 LogDebug("TransverseImpactPointExtrapolator") << "plane center: " << origin << "\n"
0158 << "plane rotation axis:\n"
0159 << xLocal << "\n"
0160 << yLocal << "\n"
0161 << zLocal << "\n"
0162 << "x0: " << x0 << "\n"
0163 << "t0: " << t0 << "\n"
0164 << "xc: " << xc << "\n"
0165 << "vtxDirection: " << vtxDirection;
0166
0167 return PlaneBuilder().plane(origin, rotation);
0168 }