GeomPropagators/src/PropagationDirectionChooser.cc

0001 #include "TrackingTools/GeomPropagators/interface/PropagationDirectionChooser.h"
0002 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0003 #include "TrackingTools/GeomPropagators/interface/PropagationExceptions.h"
0004
0005 #include "DataFormats/GeometrySurface/interface/Surface.h"
0006 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0007 #include "DataFormats/GeometrySurface/interface/Plane.h"
0008
0009 #include "MagneticField/Engine/interface/MagneticField.h"
0010
0011 #include <cmath>
0012
0013 PropagationDirection PropagationDirectionChooser::operator()(const FreeTrajectoryState& fts,
0014                                                              const Surface& surface) const {
0015   // understand the surface type (expect cylinder/plane)
0016   // unfortunately dynamic_cast on Sun is broken -- cannot deal with const
0017   // so here we get rid of const
0018   const Surface* sur = (const Surface*)&surface;
0019   const Cylinder* bc = dynamic_cast<const Cylinder*>(sur);
0020   const Plane* bp = dynamic_cast<const Plane*>(sur);
0021   if (bc != nullptr) {
0022     //
0023     // cylinder
0024     //
0025     return (*this)(fts, *bc);
0026   } else if (bp != nullptr) {
0027     //
0028     // plane
0029     //
0030     return (*this)(fts, *bp);
0031   } else {
0032     throw PropagationException("The surface is neither Cylinder nor Plane");
0033   }
0034 }
0035
0036 PropagationDirection PropagationDirectionChooser::operator()(const FreeTrajectoryState& fts, const Plane& plane) const {
0037   // propagation towards arbitrary plane
0038   // need computation of intersection point between plane
0039   // and track (straight line approximation) to set direction
0040   // Copied from BidirectionalPropagator.
0041   PropagationDirection dir = anyDirection;
0042
0043   GlobalPoint x = fts.position();
0044   GlobalVector p = fts.momentum().unit();
0045   GlobalPoint sp = plane.toGlobal(LocalPoint(0., 0.));
0046   GlobalVector v = plane.toGlobal(LocalVector(1., 0., 0.));
0047   GlobalVector w = plane.toGlobal(LocalVector(0., 1., 0.));
0048   AlgebraicMatrix33 a;
0049   a(0, 0) = v.x();
0050   a(0, 1) = w.x();
0051   a(0, 2) = -p.x();
0052   a(1, 0) = v.y();
0053   a(1, 1) = w.y();
0054   a(1, 2) = -p.y();
0055   a(2, 0) = v.z();
0056   a(2, 1) = w.z();
0057   a(2, 2) = -p.z();
0058   AlgebraicVector3 b;
0059   b[0] = x.x() - sp.x();
0060   b[1] = x.y() - sp.y();
0061   b[2] = x.z() - sp.z();
0062
0063   int ifail = !a.Invert();
0064   if (ifail == 0) {
0065     // plane and momentum are not parallel
0066     b = a * b;
0067     // b[2] nb of steps along unit vector parallel to momentum to reach plane
0068
0069     const double small = 1.e-4;  // 1 micron distance
0070     if (fabs(b[2]) < small) {
0071       // already on plane, choose forward as default
0072       dir = alongMomentum;
0073     } else if (b[2] < 0.) {
0074       dir = oppositeToMomentum;
0075     } else {
0076       dir = alongMomentum;
0077     }
0078   }
0079   return dir;
0080 }
0081
0082 PropagationDirection PropagationDirectionChooser::operator()(const FreeTrajectoryState& fts,
0083                                                              const Cylinder& cylinder) const {
0084   // propagation to cylinder with axis along z
0085   // Copied from BidirectionalPropagator.
0086   PropagationDirection dir = anyDirection;
0087
0088   GlobalPoint sp = cylinder.toGlobal(LocalPoint(0., 0.));
0089   if (sp.x() != 0. || sp.y() != 0.) {
0090     throw PropagationException("Cannot propagate to an arbitrary cylinder");
0091   }
0092
0093   GlobalPoint x = fts.position();
0094   GlobalVector p = fts.momentum();
0095
0096   const double small = 1.e-4;  // 1 micron distance
0097   double rdiff = x.perp() - cylinder.radius();
0098
0099   if (fabs(rdiff) < small) {
0100     // already on cylinder, choose forward as default
0101     dir = alongMomentum;
0102   } else {
0103     int rSign = (rdiff >= 0.) ? 1 : -1;
0104     if (rSign == -1) {
0105       // starting point inside cylinder
0106       // propagate always in direction of momentum
0107       dir = alongMomentum;
0108     } else {
0109       // starting point outside cylinder
0110       // choose direction so as to approach cylinder surface
0111       double proj = (x.x() * p.x() + x.y() * p.y()) * rSign;
0112       if (proj < 0.) {
0113         dir = alongMomentum;
0114       } else {
0115         dir = oppositeToMomentum;
0116       }
0117     }
0118   }
0119   return dir;
0120 }