GeomPropagators/src/StraightLinePropagator.cc

0001 #include "TrackingTools/GeomPropagators/interface/StraightLinePropagator.h"
0002
0003 #include "DataFormats/Math/interface/AlgebraicROOTObjects.h"
0004 #include "DataFormats/GeometrySurface/interface/Plane.h"
0005 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0006 #include "TrackingTools/GeomPropagators/interface/PropagationExceptions.h"
0007
0008 std::pair<TrajectoryStateOnSurface, double> StraightLinePropagator::propagateWithPath(const FreeTrajectoryState& fts,
0009                                                                                       const Plane& plane) const {
0010   // propagate parameters
0011   LocalPoint x;
0012   LocalVector p;
0013   double s = 0;
0014   bool parametersOK = propagateParametersOnPlane(fts, plane, x, p, s);
0015   if (!parametersOK)
0016     return std::make_pair(TrajectoryStateOnSurface(), 0.);
0017
0018   // compute propagated state
0019   if (fts.hasError()) {
0020     return std::make_pair(propagatedState(fts, plane, jacobian(s), x, p), s);
0021   } else {
0022     // return state without errors
0023     return std::make_pair(TSOS(LocalTrajectoryParameters(x, p, fts.charge()), plane, theField), s);
0024   }
0025 }
0026
0027 std::pair<TrajectoryStateOnSurface, double> StraightLinePropagator::propagateWithPath(const FreeTrajectoryState& fts,
0028                                                                                       const Cylinder& cylinder) const {
0029   // propagate parameters
0030   GlobalPoint x;
0031   GlobalVector p;
0032   double s = 0;
0033   bool parametersOK = propagateParametersOnCylinder(fts, cylinder, x, p, s);
0034   if (!parametersOK)
0035     return std::make_pair(TrajectoryStateOnSurface(), 0.);
0036
0037   // compute propagated state
0038   if (fts.hasError()) {
0039     return std::make_pair(propagatedState(fts, cylinder, jacobian(s), x, p), s);
0040   } else {
0041     // return state without errors
0042     return std::make_pair(TSOS(GlobalTrajectoryParameters(x, p, fts.charge(), theField), cylinder), s);
0043   }
0044 }
0045
0046 TrajectoryStateOnSurface StraightLinePropagator::propagatedState(const FTS& fts,
0047                                                                  const Surface& surface,
0048                                                                  const AlgebraicMatrix55& jacobian,
0049                                                                  const LocalPoint& x,
0050                                                                  const LocalVector& p) const {
0051   if (fts.hasError()) {
0052     // propagate error
0053     TSOS tmp(fts, surface);
0054     const AlgebraicSymMatrix55& eLocal = tmp.localError().matrix();
0055     AlgebraicSymMatrix55 lte = ROOT::Math::Similarity(jacobian, eLocal);
0056     LocalTrajectoryError eloc(lte);
0057     LocalTrajectoryParameters ltp(x, p, fts.charge());
0058     return TSOS(ltp, eloc, surface, theField);
0059   } else {
0060     // return state without errors
0061     return TSOS(LocalTrajectoryParameters(x, p, fts.charge()), surface, theField);
0062   }
0063 }
0064
0065 TrajectoryStateOnSurface StraightLinePropagator::propagatedState(const FTS& fts,
0066                                                                  const Surface& surface,
0067                                                                  const AlgebraicMatrix55& jacobian,
0068                                                                  const GlobalPoint& x,
0069                                                                  const GlobalVector& p) const {
0070   if (fts.hasError()) {
0071     // propagate error
0072     TSOS tmp(fts, surface);
0073     const AlgebraicSymMatrix55& eLocal = tmp.localError().matrix();
0074     AlgebraicSymMatrix55 lte = ROOT::Math::Similarity(jacobian, eLocal);
0075     LocalTrajectoryError eloc(lte);
0076
0077     TSOS tmp2(tmp.localParameters(), eloc, surface, theField);
0078     GlobalTrajectoryParameters gtp(x, p, fts.charge(), theField);
0079     return TSOS(gtp, tmp2.cartesianError(), surface);
0080   } else {
0081     // return state without errors
0082     return TSOS(GlobalTrajectoryParameters(x, p, fts.charge(), theField), surface);
0083   }
0084 }
0085
0086 AlgebraicMatrix55 StraightLinePropagator::jacobian(double& s) const {
0087   //Jacobian for 5*5 local error matrix
0088   AlgebraicMatrix55 j = AlgebraicMatrixID();  //Jacobian
0089
0090   double dir = (propagationDirection() == alongMomentum) ? 1. : -1.;
0091   if (s * dir < 0.)
0092     return j;
0093
0094   j(3, 1) = s;
0095   j(4, 2) = s;
0096
0097   return j;
0098 }
0099
0100 bool StraightLinePropagator::propagateParametersOnCylinder(
0101     const FTS& fts, const Cylinder& cylinder, GlobalPoint& x, GlobalVector& p, double& s) const {
0102   GlobalPoint sp = cylinder.toGlobal(LocalPoint(0., 0.));
0103   if (sp.x() != 0. || sp.y() != 0.) {
0104     throw PropagationException("Cannot propagate to an arbitrary cylinder");
0105   }
0106
0107   x = fts.position();
0108   p = fts.momentum();
0109   s = cylinder.radius() - x.perp();
0110
0111   double dir = (propagationDirection() == alongMomentum) ? 1. : -1.;
0112   if (s * dir < 0.)
0113     return false;
0114
0115   double dt = s / p.perp();
0116   x = GlobalPoint(x.x() + p.x() * dt, x.y() + p.y() * dt, x.z() + p.z() * dt);
0117
0118   return true;
0119 }
0120
0121 bool StraightLinePropagator::propagateParametersOnPlane(
0122     const FTS& fts, const Plane& plane, LocalPoint& x, LocalVector& p, double& s) const {
0123   //Do extrapolation in local frame of plane
0124   //  LocalPoint sp = plane.toLocal(plane.position());
0125   x = plane.toLocal(fts.position());
0126   p = plane.toLocal(fts.momentum());
0127   s = -x.z();  // sp.z() - x.z(); local z of plane always 0
0128
0129   //double dir = (propagationDirection() == alongMomentum) ? 1. : -1.;
0130   //if(s*dir < 0.) return false;
0131   if ((p.x() != 0 || p.y() != 0) && p.z() == 0 && s != 0)
0132     return false;
0133
0134   x = LocalPoint(x.x() + (p.x() / p.z()) * s, x.y() + (p.y() / p.z()) * s, x.z() + s);
0135
0136   return true;
0137 }