GeomPropagators/src/AnalyticalTrajectoryExtrapolatorToLine.cc

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