TkSeedGenerator/src/FastCircle.cc

0001 #include "RecoTracker/TkSeedGenerator/interface/FastCircle.h"
0002 #include "DataFormats/Math/interface/AlgebraicROOTObjects.h"
0003
0004 FastCircle::FastCircle(const GlobalPoint& outerHit, const GlobalPoint& middleHit, const GlobalPoint& aVertex)
0005     : theOuterPoint(outerHit),
0006       theInnerPoint(middleHit),
0007       theVertexPoint(aVertex),
0008       theNorm(128.),
0009       theX0(0.),
0010       theY0(0.),
0011       theRho(0.),
0012       theN1(0.),
0013       theN2(0.),
0014       theC(0.),
0015       theValid(true),
0016       theIsLine(false) {
0017   createCircleParameters();
0018 }
0019
0020 FastCircle::FastCircle(const GlobalPoint& outerHit,
0021                        const GlobalPoint& middleHit,
0022                        const GlobalPoint& aVertex,
0023                        double norm)
0024     : theOuterPoint(outerHit),
0025       theInnerPoint(middleHit),
0026       theVertexPoint(aVertex),
0027       theNorm(norm),
0028       theX0(0.),
0029       theY0(0.),
0030       theRho(0.),
0031       theN1(0.),
0032       theN2(0.),
0033       theC(0.),
0034       theValid(true),
0035       theIsLine(false) {
0036   createCircleParameters();
0037 }
0038
0039 namespace {
0040   inline AlgebraicVector3 transform(const GlobalPoint& aPoint, float norm) {
0041     AlgebraicVector3 riemannPoint;
0042
0043     auto p = aPoint.basicVector() / norm;
0044     float R2 = p.perp2();
0045     float fact = 1.f / (1.f + R2);  // let's factorize the common factor out
0046     riemannPoint[0] = fact * p.x();
0047     riemannPoint[1] = fact * p.y();
0048     riemannPoint[2] = fact * R2;
0049
0050     return riemannPoint;
0051   }
0052
0053 }  // namespace
0054
0055 void FastCircle::createCircleParameters() {
0056   AlgebraicVector3 x = transform(theOuterPoint, theNorm);
0057   AlgebraicVector3 y = transform(theInnerPoint, theNorm);
0058   AlgebraicVector3 z = transform(theVertexPoint, theNorm);
0059
0060   AlgebraicVector3 n;
0061
0062   n[0] = x[1] * (y[2] - z[2]) + y[1] * (z[2] - x[2]) + z[1] * (x[2] - y[2]);
0063   n[1] = -(x[0] * (y[2] - z[2]) + y[0] * (z[2] - x[2]) + z[0] * (x[2] - y[2]));
0064   n[2] = x[0] * (y[1] - z[1]) + y[0] * (z[1] - x[1]) + z[0] * (x[1] - y[1]);
0065
0066   double mag2 = n[0] * n[0] + n[1] * n[1] + n[2] * n[2];
0067   if (mag2 < 1.e-20) {
0068     theValid = false;
0069     return;
0070   }
0071   n.Unit();  // reduce n to a unit vector
0072   double c = -(n[0] * x[0] + n[1] * x[1] + n[2] * x[2]);
0073   //  c = -(n[0]*y[0] + n[1]*y[1] + n[2]*y[2]);
0074   //  c = -(n[0]*z[0] + n[1]*z[1] + n[2]*z[2]);
0075
0076   theN1 = n[0];
0077   theN2 = n[1];
0078   theC = c;
0079
0080   if (fabs(c + n[2]) < 1.e-5) {
0081     // numeric limit
0082     // circle is more a straight line...
0083     theValid = false;
0084     theIsLine = true;
0085     return;
0086   }
0087
0088   double x0 = -n[0] / (2. * (c + n[2]));
0089   double y0 = -n[1] / (2. * (c + n[2]));
0090   double rho = sqrt((n[0] * n[0] + n[1] * n[1] - 4. * c * (c + n[2]))) / fabs(2. * (c + n[2]));
0091
0092   theX0 = theNorm * x0;
0093   theY0 = theNorm * y0;
0094   theRho = theNorm * rho;
0095 }