PatternTools/src/TransverseImpactPointExtrapolator.cc

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   // set propagator for bidirectional propagation
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   // set propagator for bidirectional propagation
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   // Compute tip surface
0062   //
0063   if (fabs(tsos.freeState()->transverseCurvature()) < 1.E-6) {
0064     LogDebug("TransverseImpactPointExtrapolator") << "negligeable curvature: using a trick to extrapolate:\n" << tsos;
0065
0066     //0T field probably
0067     //x is perpendicular to the momentum
0068     GlobalVector xLocal = GlobalVector(-tsos.globalMomentum().y(), tsos.globalMomentum().x(), 0).unit();
0069     //y along global Z
0070     GlobalVector yLocal(0., 0., 1.);
0071     //z accordingly
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     // propagate
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   // Compute tip surface
0094   //
0095   if (fabs(fts.transverseCurvature()) < 1.E-6) {
0096     LogDebug("TransverseImpactPointExtrapolator") << "negligeable curvature: using a trick to extrapolate:\n" << fts;
0097
0098     //0T field probably
0099     //x is perpendicular to the momentum
0100     GlobalVector xLocal = GlobalVector(-fts.momentum().y(), fts.momentum().x(), 0).unit();
0101     //y along global Z
0102     GlobalVector yLocal(0., 0., 1.);
0103     //z accordingly
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     // propagate
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 }