GeomPropagators/src/StraightLineCylinderCrossing.cc

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   // radius of cylinder and transversal position relative to axis
0019   //
0020   double R(cylinder.radius());
0021   PositionType2D xt2d(theX0.x(), theX0.y());
0022   //
0023   // transverse direction
0024   //
0025   DirectionType2D pt2d(theP0.x(), theP0.y());
0026   //
0027   // solution of quadratic equation for s - assume |theP0|=1
0028   //
0029   RealQuadEquation eq(pt2d.mag2(), 2. * xt2d.dot(pt2d), xt2d.mag2() - R * R);
0030   if (!eq.hasSolution) {
0031     /*
0032     double A=   pt2d.mag2();
0033     double B=   2.*xt2d.dot(pt2d);
0034     double C=   xt2d.mag2()-R*R;
0035     cout << "A= " << pt2d.mag2()
0036      << " B= " << 2.*xt2d.dot(pt2d)
0037      << " C= " << xt2d.mag2()-R*R
0038      << " D= " << B*B - 4*A*C
0039      << endl;
0040     */
0041     return std::pair<bool, double>(false, 0.);
0042   }
0043   //
0044   // choice of solution and verification of direction
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   // follows the logic implemented in HelixBarrelCylinderCrossing
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       // if different signs return the positive one
0059       return std::pair<bool, double>(true, ((s1 * propSign > 0) ? s1 : s2));
0060     } else if (s1 * propSign > 0) {
0061       // if both positive, return the shortest
0062       return std::pair<bool, double>(true, (fabs(s1) < fabs(s2) ? s1 : s2));
0063     } else {
0064       // if both negative, check if the smaller (abs value) is smaller then tolerance
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 }