PixelTriplets/test/fastDPHI_t.cpp

0001 // this test documents the derivation of the fast deltaphi used in gpu doublet code..
0002 //
0003 //
0004 //
0005 #include <cmath>
0006 #include <algorithm>
0007 #include <numeric>
0008 #include <cassert>
0009
0010 /**
0011 | 1) circle is parameterized as:                                              |
0012 |    C*[(X-Xp)**2+(Y-Yp)**2] - 2*alpha*(X-Xp) - 2*beta*(Y-Yp) = 0             |
0013 |    Xp,Yp is a point on the track (Yp is at the center of the chamber);      |
0014 |    C = 1/r0 is the curvature  ( sign of C is charge of particle );          |
0015 |    alpha & beta are the direction cosines of the radial vector at Xp,Yp     |
0016 |    i.e.  alpha = C*(X0-Xp),                                                 |
0017 |          beta  = C*(Y0-Yp),                                                 |
0018 |    where center of circle is at X0,Y0.                                      |
0019 |    Alpha > 0                                                                |
0020 |    Slope dy/dx of tangent at Xp,Yp is -alpha/beta.                          |
0021 | 2) the z dimension of the helix is parameterized by gamma = dZ/dSperp       |
0022 |    this is also the tangent of the pitch angle of the helix.                |
0023 |    with this parameterization, (alpha,beta,gamma) rotate like a vector.     |
0024 | 3) For tracks going inward at (Xp,Yp), C, alpha, beta, and gamma change sign|
0025 |
0026 */
0027
0028 template <typename T>
0029 class FastCircle {
0030 public:
0031   FastCircle() {}
0032   FastCircle(T x1, T y1, T x2, T y2, T x3, T y3) { compute(x1, y1, x2, y2, x3, y3); }
0033
0034   void compute(T x1, T y1, T x2, T y2, T x3, T y3);
0035
0036   T m_xp;
0037   T m_yp;
0038   T m_c;
0039   T m_alpha;
0040   T m_beta;
0041 };
0042
0043 template <typename T>
0044 void FastCircle<T>::compute(T x1, T y1, T x2, T y2, T x3, T y3) {
0045   bool flip = std::abs(x3 - x1) > std::abs(y3 - y1);
0046
0047   auto x1p = x1 - x2;
0048   auto y1p = y1 - y2;
0049   auto d12 = x1p * x1p + y1p * y1p;
0050   auto x3p = x3 - x2;
0051   auto y3p = y3 - y2;
0052   auto d32 = x3p * x3p + y3p * y3p;
0053
0054   if (flip) {
0055     std::swap(x1p, y1p);
0056     std::swap(x3p, y3p);
0057   }
0058
0059   auto num = x1p * y3p - y1p * x3p;  // num also gives correct sign for CT
0060   auto det = d12 * y3p - d32 * y1p;
0061   if (std::abs(det) == 0) {
0062     // and why we flip????
0063   }
0064   auto ct = num / det;
0065   auto sn = det > 0 ? T(1.) : T(-1.);
0066   auto st2 = (d12 * x3p - d32 * x1p) / det;
0067   auto seq = T(1.) + st2 * st2;
0068   auto al2 = sn / std::sqrt(seq);
0069   auto be2 = -st2 * al2;
0070   ct *= T(2.) * al2;
0071
0072   if (flip) {
0073     std::swap(x1p, y1p);
0074     std::swap(al2, be2);
0075     al2 = -al2;
0076     be2 = -be2;
0077     ct = -ct;
0078   }
0079
0080   m_xp = x1;
0081   m_yp = y1;
0082   m_c = ct;
0083   m_alpha = al2 - ct * x1p;
0084   m_beta = be2 - ct * y1p;
0085 }
0086
0087 // compute curvature given two points (and origin)
0088 float fastDPHI(float ri, float ro, float dphi) {
0089   /*
0090   x3=0 y1=0 x1=0;
0091   y3=ro
0092   */
0093
0094   // auto x2 = ri*dphi;
0095   // auto y2 = ri*(1.f-0.5f*dphi*dphi);
0096
0097   /*
0098   auto x1p = x1-x2;
0099   auto y1p = y1-y2;
0100   auto d12 = x1p*x1p + y1p*y1p;
0101   auto x3p = x3-x2;
0102   auto y3p = y3-y2;
0103   auto d32 = x3p*x3p + y3p*y3p;
0104   */
0105
0106   /*
0107   auto x1p = -x2;
0108   auto y1p = -y2;
0109   auto d12 = ri*ri;
0110   auto x3p = -x2;
0111   auto y3p = ro-y2;
0112   auto d32 = ri*ri + ro*ro - 2.f*ro*y2;
0113   */
0114
0115   // auto rat = (ro -2.f*y2);
0116   // auto det =  ro - ri - (ro - 2.f*ri -0.5f*ro)*dphi*dphi;
0117
0118   //auto det2 = (ro-ri)*(ro-ri) -2.*(ro-ri)*(ro - 2.f*ri -0.5f*ro)*dphi*dphi;
0119   // auto seq = det2 +  dphi*dphi*(ro-2.f*ri)*(ro-2.f*ri);    // *rat2;
0120   // auto seq = (ro-ri)*(ro-ri) +  dphi*dphi*ri*ro;
0121
0122   // and little by little simplifing and removing higher over terms
0123   // we get
0124   auto r2 = (ro - ri) * (ro - ri) / (dphi * dphi) + ri * ro;
0125
0126   // d2 = (ro-ri)*(ro-ri)/(4.f*r2 -ri*ro);
0127   // return -2.f*dphi/std::sqrt(seq);
0128
0129   return -1.f / std::sqrt(r2 / 4.f);
0130 }
0131
0132 #include <iostream>
0133
0134 template <typename T>
0135 bool equal(T a, T b) {
0136   //  return float(a-b)==0;
0137   return std::abs(float(a - b)) < std::abs(0.01f * a);
0138 }
0139
0140 int n = 0;
0141 void go(float ri, float ro, float dphi, bool print = false) {
0142   ++n;
0143   float x3 = 0.f, y3 = ro;
0144   float x2 = ri * sin(dphi);
0145   float y2 = ri * cos(dphi);
0146
0147   FastCircle<float> c(0, 0, x2, y2, x3, y3);
0148
0149   auto cc = fastDPHI(ri, ro, dphi);
0150   if (print)
0151     std::cout << c.m_c << ' ' << cc << std::endl;
0152   assert(equal(c.m_c, cc));
0153 }
0154
0155 int main() {
0156   go(4., 7., 0.1, true);
0157
0158   for (float r1 = 2; r1 < 15; r1 += 1)
0159     for (float dr = 0.5; dr < 10; dr += 0.5)
0160       for (float dphi = 0.02; dphi < 0.2; dphi += 0.2)
0161         go(r1, r1 + dr, dphi);
0162
0163   std::cout << "done " << n << std::endl;
0164   return 0;
0165 };