File indexing completed on 2023-03-17 11:26:23
0001 #include "TrackingTools/GeomPropagators/interface/StraightLineCylinderCrossing.h"
0002
0003 #include "DataFormats/GeometrySurface/interface/Cylinder.h"
0004 #include "TrackingTools/GeomPropagators/src/RealQuadEquation.h"
0005
0006 #include <cmath>
0007 #include <iostream>
0008 using namespace std;
0009
0010 StraightLineCylinderCrossing::StraightLineCylinderCrossing(const LocalPoint& startingPos,
0011 const LocalVector& startingDir,
0012 const PropagationDirection propDir,
0013 double tolerance)
0014 : theX0(startingPos), theP0(startingDir.unit()), thePropDir(propDir), theTolerance(tolerance) {}
0015
0016 std::pair<bool, double> StraightLineCylinderCrossing::pathLength(const Cylinder& cylinder) const {
0017
0018
0019
0020 double R(cylinder.radius());
0021 PositionType2D xt2d(theX0.x(), theX0.y());
0022
0023
0024
0025 DirectionType2D pt2d(theP0.x(), theP0.y());
0026
0027
0028
0029 RealQuadEquation eq(pt2d.mag2(), 2. * xt2d.dot(pt2d), xt2d.mag2() - R * R);
0030 if (!eq.hasSolution) {
0031
0032
0033
0034
0035
0036
0037
0038
0039
0040
0041 return std::pair<bool, double>(false, 0.);
0042 }
0043
0044
0045
0046 return chooseSolution(eq.first, eq.second);
0047 }
0048
0049 std::pair<bool, double> StraightLineCylinderCrossing::chooseSolution(const double s1, const double s2) const {
0050
0051
0052
0053 if (thePropDir == anyDirection) {
0054 return std::pair<bool, double>(true, (fabs(s1) < fabs(s2) ? s1 : s2));
0055 } else {
0056 int propSign = thePropDir == alongMomentum ? 1 : -1;
0057 if (s1 * s2 < 0) {
0058
0059 return std::pair<bool, double>(true, ((s1 * propSign > 0) ? s1 : s2));
0060 } else if (s1 * propSign > 0) {
0061
0062 return std::pair<bool, double>(true, (fabs(s1) < fabs(s2) ? s1 : s2));
0063 } else {
0064
0065 double shorter = std::min(fabs(s1), fabs(s2));
0066 if (shorter < theTolerance)
0067 return std::pair<bool, double>(true, 0);
0068 else
0069 return std::pair<bool, double>(false, 0.);
0070 }
0071 }
0072 }