GeomPropagators/src/AnalyticalImpactPointExtrapolator.cc

0001 #include "TrackingTools/GeomPropagators/interface/AnalyticalImpactPointExtrapolator.h"
0002 #include "TrackingTools/GeomPropagators/interface/IterativeHelixExtrapolatorToLine.h"
0003 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0004 #include "TrackingTools/TrajectoryState/interface/FreeTrajectoryState.h"
0005 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0006 #include "TrackingTools/AnalyticalJacobians/interface/AnalyticalCurvilinearJacobian.h"
0007 #include "DataFormats/GeometrySurface/interface/Plane.h"
0008 #include "DataFormats/GeometrySurface/interface/PlaneBuilder.h"
0009
0010 // #include "CommonDet/DetUtilities/interface/DetailedDetTimer.h"
0011
0012 AnalyticalImpactPointExtrapolator::AnalyticalImpactPointExtrapolator(const MagneticField* field)
0013     : thePropagator(new AnalyticalPropagator(field, anyDirection)), theField(field) {}
0014
0015 AnalyticalImpactPointExtrapolator::AnalyticalImpactPointExtrapolator(const Propagator& propagator,
0016                                                                      const MagneticField* field)
0017     : thePropagator(propagator.clone()), theField(field) {
0018   thePropagator->setPropagationDirection(anyDirection);
0019 }
0020
0021 TrajectoryStateOnSurface AnalyticalImpactPointExtrapolator::extrapolate(const FreeTrajectoryState& fts,
0022                                                                         const GlobalPoint& vtx) const {
0023   //   static TimingReport::Item& timer = detailedDetTimer("AnalyticalImpactPointExtrapolator");
0024   //   TimeMe t(timer,false);
0025
0026   return extrapolateSingleState(fts, vtx);
0027 }
0028
0029 TrajectoryStateOnSurface AnalyticalImpactPointExtrapolator::extrapolate(const TrajectoryStateOnSurface tsos,
0030                                                                         const GlobalPoint& vtx) const {
0031   if (tsos.isValid())
0032     return extrapolateFullState(tsos, vtx);
0033   else
0034     return tsos;
0035 }
0036
0037 TrajectoryStateOnSurface AnalyticalImpactPointExtrapolator::extrapolateFullState(const TrajectoryStateOnSurface tsos,
0038                                                                                  const GlobalPoint& vertex) const {
0039   //
0040   // first determine IP plane using propagation with (single) FTS
0041   // could be optimised (will propagate errors even if duplicated below)
0042   //
0043   TrajectoryStateOnSurface singleState = extrapolateSingleState(*tsos.freeTrajectoryState(), vertex);
0044   if (!singleState.isValid() || tsos.singleState())
0045     return singleState;
0046   //
0047   // propagate multiTsos to plane found above
0048   //
0049   return thePropagator->propagate(tsos, singleState.surface());
0050 }
0051
0052 TrajectoryStateOnSurface AnalyticalImpactPointExtrapolator::extrapolateSingleState(const FreeTrajectoryState& fts,
0053                                                                                    const GlobalPoint& vertex) const {
0054   //
0055   // initialisation of position, momentum and transverse curvature
0056   //
0057   GlobalPoint x(fts.position());
0058   GlobalVector p(fts.momentum());
0059   double rho = fts.transverseCurvature();
0060   //
0061   // Straight line approximation? |rho|<1.e-10 equivalent to ~ 1um
0062   // difference in transversal position at 10m.
0063   //
0064   double s(0);
0065   if (fabs(rho) < 1.e-10) {
0066     GlobalVector dx(p.dot(vertex - x) / p.mag2() * p);
0067     x += dx;
0068     float sign = p.dot(dx);
0069     s = sign > 0 ? dx.mag() : -dx.mag();
0070   }
0071   //
0072   // Helix case
0073   //
0074   else {
0075     HelixLineExtrapolation::PositionType helixPos(x);
0076     HelixLineExtrapolation::DirectionType helixDir(p);
0077     IterativeHelixExtrapolatorToLine extrapolator(helixPos, helixDir, rho, anyDirection);
0078     if (!propagateWithHelix(extrapolator, vertex, x, p, s))
0079       return TrajectoryStateOnSurface();
0080   }
0081   //
0082   // Define target surface: origin on line, x_local from line
0083   //   to helix at closest approach, z_local along the helix
0084   //   and y_local to complete right-handed system
0085   //
0086   GlobalVector zLocal(p.unit());
0087   GlobalVector yLocal(zLocal.cross(x - vertex).unit());
0088   GlobalVector xLocal(yLocal.cross(zLocal));
0089   Surface::RotationType rot(xLocal, yLocal, zLocal);
0090   PlaneBuilder::ReturnType surface = PlaneBuilder().plane(vertex, rot);
0091   //
0092   // Compute propagated state
0093   //
0094   GlobalTrajectoryParameters gtp(x, p, fts.charge(), theField);
0095   if (fts.hasError()) {
0096     //
0097     // compute jacobian
0098     //
0099     AnalyticalCurvilinearJacobian analyticalJacobian(fts.parameters(), gtp.position(), gtp.momentum(), s);
0100     CurvilinearTrajectoryError cte(
0101         ROOT::Math::Similarity(analyticalJacobian.jacobian(), fts.curvilinearError().matrix()));
0102     return TrajectoryStateOnSurface(gtp, cte, *surface);
0103   } else {
0104     //
0105     // return state without errors
0106     //
0107     return TrajectoryStateOnSurface(gtp, *surface);
0108   }
0109 }
0110
0111 bool AnalyticalImpactPointExtrapolator::propagateWithHelix(const IterativeHelixExtrapolatorToLine& extrapolator,
0112                                                            const GlobalPoint& vertex,
0113                                                            GlobalPoint& x,
0114                                                            GlobalVector& p,
0115                                                            double& s) const {
0116   //
0117   // save absolute value of momentum
0118   //
0119   double pmag(p.mag());
0120   //
0121   // get path length to solution
0122   //
0123   std::pair<bool, double> propResult = extrapolator.pathLength(vertex);
0124   if (!propResult.first)
0125     return false;
0126   s = propResult.second;
0127   //
0128   // get point and (normalised) direction from path length
0129   //
0130   HelixLineExtrapolation::PositionType xGen = extrapolator.position(s);
0131   HelixLineExtrapolation::DirectionType pGen = extrapolator.direction(s);
0132   //
0133   // Fix normalisation and convert back to GlobalPoint / GlobalVector
0134   //
0135   x = GlobalPoint(xGen);
0136   pGen *= pmag / pGen.mag();
0137   p = GlobalVector(pGen);
0138   //
0139   return true;
0140 }