CMS 3D CMS Logo

 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
CircleFromThreePoints.cc
Go to the documentation of this file.
2 
4  const GlobalPoint& mid,
5  const GlobalPoint& outer,
6  double precision) {
7  // move to frame where inner.x() == inner.y() == 0;
8  Vector2D b(mid.x() - inner.x(), mid.y() - inner.y());
9  Vector2D c(outer.x() - inner.x(), outer.y() - inner.y());
10  init(b, c, Vector2D(inner.x(), inner.y()), precision);
11 }
12 
13 void CircleFromThreePoints::init(const Vector2D& b, const Vector2D& c, const Vector2D& offset, double precision) {
14  double b2 = b.mag2();
15  double c2 = c.mag2();
16 
17  double oX(0), oY(0);
18  bool solved = false;
19  if (fabs(b.x()) > fabs(b.y())) { // solve for y first
20  double k = c.x() / b.x();
21  double div = 2 * (k * b.y() - c.y());
22  if (fabs(div) < precision)
23  theCurvature = 0; // if the 3 points lie on a line
24  else {
25  oY = (k * b2 - c2) / div;
26  oX = b2 / (2 * b.x()) - b.y() / b.x() * oY;
27  solved = true;
28  }
29  } else { // solve for x first
30  double k = c.y() / b.y();
31  double div = 2 * (k * b.x() - c.x());
32  if (fabs(div) < precision)
33  theCurvature = 0; // if the 3 points lie on a line
34  else {
35  oX = (k * b2 - c2) / div;
36  oY = b2 / (2 * b.y()) - b.x() / b.y() * oX;
37  solved = true;
38  }
39  }
40  if (solved) {
41  theCurvature = 1. / sqrt(oX * oX + oY * oY);
42  double xC = oX + offset.x();
43  double yC = oY + offset.y();
44  theCenter = Vector2D(xC, yC);
45  // thePhi = acos(xC/sqrt(xC*xC + yC*yC));
46 
47  // if (xC<0.) thePhi = thePhi - PI;
48  // cout << setiosflags(ios::showpoint | ios::fixed);
49  //
50  // cout << "CircleFromThreePoints::init curv = " << theCurvature << endl;
51  // cout << "CircleFromThreePoints::init center prime = " << oX << " " << oY << endl;
52  // cout << "CircleFromThreePoints::init offset = " << offset.x() << " " << offset.y() << endl;
53  // cout << "CircleFromThreePoints::init center = " << theCenter.x()<< " " << theCenter.y() << endl;
54  //
55  // float d = sqrt(theCenter.x()*theCenter.x()+theCenter.y()*theCenter.y());
56  // cout << "CircleFromThreePoints::init dfloat = " << setw(10) << setprecision(5) << d << endl;
57  // cout << "CircleFromThreePoints::init radius = " << 1/theCurvature << endl;
58  }
59 }
const edm::EventSetup & c
T y() const
Definition: PV3DBase.h:60
T sqrt(T t)
Definition: SSEVec.h:19
T y() const
Cartesian y coordinate.
void init(const Vector2D &b, const Vector2D &c, const Vector2D &offset, double precision)
double b
Definition: hdecay.h:118
T mag2() const
The vector magnitude squared. Equivalent to vec.dot(vec)
Basic2DVector< float > Vector2D
static constexpr float b2
T x() const
Definition: PV3DBase.h:59
T x() const
Cartesian x coordinate.