PixelSeeding/interface/CircleEq.h

0001 #ifndef RecoTracker_PixelSeeding_CircleEq_h
0002 #define RecoTracker_PixelSeeding_CircleEq_h
0003 /**
0004 | 1) circle is parameterized as:                                              |
0005 |    C*[(X-Xp)**2+(Y-Yp)**2] - 2*alpha*(X-Xp) - 2*beta*(Y-Yp) = 0             |
0006 |    Xp,Yp is a point on the track;                                           |
0007 |    C = 1/r0 is the curvature  ( sign of C is charge of particle );          |
0008 |    alpha & beta are the direction cosines of the radial vector at Xp,Yp     |
0009 |    i.e.  alpha = C*(X0-Xp),                                                 |
0010 |          beta  = C*(Y0-Yp),                                                 |
0011 |    where center of circle is at X0,Y0.                                      |
0012 |                                                                             |
0013 |    Slope dy/dx of tangent at Xp,Yp is -alpha/beta.                          |
0014 | 2) the z dimension of the helix is parameterized by gamma = dZ/dSperp       |
0015 |    this is also the tangent of the pitch angle of the helix.                |
0016 |    with this parameterization, (alpha,beta,gamma) rotate like a vector.     |
0017 | 3) For tracks going inward at (Xp,Yp), C, alpha, beta, and gamma change sign|
0018 |
0019 */
0020
0021 #include <cmath>
0022
0023 template <typename T>
0024 class CircleEq {
0025 public:
0026   CircleEq() {}
0027
0028   constexpr CircleEq(T x1, T y1, T x2, T y2, T x3, T y3) { compute(x1, y1, x2, y2, x3, y3); }
0029
0030   constexpr void compute(T x1, T y1, T x2, T y2, T x3, T y3);
0031
0032   // dca to origin divided by curvature
0033   constexpr T dca0() const {
0034     auto x = m_c * m_xp + m_alpha;
0035     auto y = m_c * m_yp + m_beta;
0036     return std::sqrt(x * x + y * y) - T(1);
0037   }
0038
0039   // dca to given point (divided by curvature)
0040   constexpr T dca(T x, T y) const {
0041     x = m_c * (m_xp - x) + m_alpha;
0042     y = m_c * (m_yp - y) + m_beta;
0043     return std::sqrt(x * x + y * y) - T(1);
0044   }
0045
0046   // curvature
0047   constexpr auto curvature() const { return m_c; }
0048
0049   // alpha and beta
0050   constexpr std::pair<T, T> cosdir() const { return std::make_pair(m_alpha, m_beta); }
0051
0052   // alpha and beta af given point
0053   constexpr std::pair<T, T> cosdir(T x, T y) const {
0054     return std::make_pair(m_alpha - m_c * (x - m_xp), m_beta - m_c * (y - m_yp));
0055   }
0056
0057   // center
0058   constexpr std::pair<T, T> center() const { return std::make_pair(m_xp + m_alpha / m_c, m_yp + m_beta / m_c); }
0059
0060   constexpr auto radius() const { return T(1) / m_c; }
0061
0062   T m_xp = 0;
0063   T m_yp = 0;
0064   T m_c = 0;
0065   T m_alpha = 0;
0066   T m_beta = 0;
0067 };
0068
0069 template <typename T>
0070 constexpr void CircleEq<T>::compute(T x1, T y1, T x2, T y2, T x3, T y3) {
0071   bool noflip = std::abs(x3 - x1) < std::abs(y3 - y1);
0072
0073   auto x1p = noflip ? x1 - x2 : y1 - y2;
0074   auto y1p = noflip ? y1 - y2 : x1 - x2;
0075   auto d12 = x1p * x1p + y1p * y1p;
0076   auto x3p = noflip ? x3 - x2 : y3 - y2;
0077   auto y3p = noflip ? y3 - y2 : x3 - x2;
0078   auto d32 = x3p * x3p + y3p * y3p;
0079
0080   auto num = x1p * y3p - y1p * x3p;  // num also gives correct sign for CT
0081   auto det = d12 * y3p - d32 * y1p;
0082
0083   auto st2 = (d12 * x3p - d32 * x1p);
0084   auto seq = det * det + st2 * st2;
0085   auto al2 = T(1.) / std::sqrt(seq);
0086   auto be2 = -st2 * al2;
0087   auto ct = T(2.) * num * al2;
0088   al2 *= det;
0089
0090   m_xp = x2;
0091   m_yp = y2;
0092   m_c = noflip ? ct : -ct;
0093   m_alpha = noflip ? al2 : -be2;
0094   m_beta = noflip ? be2 : -al2;
0095 }
0096
0097 #endif