60 if (!
g)
return ellipse;
61 if (
g->GetN() < 6)
return ellipse;
79 for (i = 0; i <
N; i++) {
80 double x = (
g->GetX())[i] - X0;
81 double y = (
g->GetY())[i] - Y0;
96 S3.
Invert(&tmp); S3 *= -1.0;
98 std::cout <<
"fit_ellipse : linear part of the scatter matrix is singular!" << std::endl;
106 for (i = 0; i < 3; i++) {
108 M[0][i] = M[2][i] / 2.0;
117 std::cout <<
"fit_ellipse : eigenvalues have nonzero imaginary parts!" << std::endl;
121 for (i = 0; i < 3; i++) {
122 tmp = 4.0 * evec[0][i] * evec[2][i] - evec[1][i] * evec[1][i];
123 if (tmp > 0.0)
break;
126 std::cout <<
"fit_ellipse : no min. pos. eigenvalue found!" << std::endl;
134 a1 *= 1.0 / std::sqrt(tmp);
136 std::cout <<
"fit_ellipse : eigenvector for min. pos. eigenvalue is NULL!" << std::endl;
189 std::cout <<
"ConicToParametric : improper input vector length!" << std::endl;
194 double x0 = conic[0];
195 double y0 = conic[1];
199 double B = conic[3] / 2.0;
201 double D = conic[5] / 2.0;
202 double F = conic[6] / 2.0;
205 double J = B * B - A * C;
206 double Delta = A *
F *
F + C * D * D + J *
G - 2.0 * B * D *
F;
207 double I = - (A + C);
210 if (!( (Delta != 0.0) && (J < 0.0) && (
I != 0.0) && (Delta /
I < 0.0) )) {
211 std::cout <<
"ConicToParametric : ellipse (real) specific constraints not met!" << std::endl;
215 x0 += (C * D - B *
F) / J;
216 y0 += (A *
F - B * D) / J;
218 double tmp = std::sqrt((A - C) * (A - C) + 4.0 * B * B);
219 a = std::sqrt(2.0 * Delta / J / (
I + tmp));
220 b = std::sqrt(2.0 * Delta / J / (
I - tmp));
224 tmp = (A - C) / 2.0 / B;
226 if (tmp < 0.0) { theta -= 45.0; }
else { theta += 45.0; }
227 if (A > C) theta += 90.0;
228 }
else if (A > C) theta = 90.0;
231 if (
a <
b) { tmp =
a;
a =
b;
b = tmp; theta -= 90.0; }
233 if (theta < -45.0) theta += 180.0;
234 if (theta > 135.0) theta -= 180.0;
273 for (
int i = 0; i <
n; i++) {
274 x[i] =
a * (std::cos(dt *
double(i)) + 0.1 *
gRandom->
Rndm() - 0.05);
275 y[i] =
b * (std::sin(dt *
double(i)) + 0.1 *
gRandom->
Rndm() - 0.05);
278 x[i] =
x[i] * std::cos(theta) -
y[i] * std::sin(theta);
279 y[i] =
y[i] * std::cos(theta) + tmp * std::sin(theta);
289 g->SetNameTitle(
"g",
"test ellipse");
306 std::cout << std::endl;
307 std::cout <<
"x0 = " << ellipse[0] << std::endl;
308 std::cout <<
"y0 = " << ellipse[1] << std::endl;
309 std::cout <<
"a = " << ellipse[2] << std::endl;
310 std::cout <<
"b = " << ellipse[3] << std::endl;
311 std::cout <<
"theta = " << ellipse[4] << std::endl;
312 std::cout << std::endl;
316 auto c1 =
new TCanvas(
"c1",
"c1", 1000, 800);
321 ellipse[2], ellipse[3],
TMatrixTColumn_const< Double_t > TMatrixDColumn_const
R__EXTERN TRandom * gRandom
A TGraph is an object made of two arrays X and Y with npoints each.
const TMatrixD & GetEigenVectors() const
const TVectorD & GetEigenValuesIm() const
TMatrixT< Element > & Invert(Double_t *det=nullptr)
Invert the matrix and calculate its determinant.
Double_t Rndm() override
Machine independent random number generator.
TVectorT< Element > & ResizeTo(Int_t lwb, Int_t upb)
Resize the vector to [lwb:upb] .
Element Norm2Sqr() const
Compute the square of the 2-norm SUM{ v[i]^2 }.
TVectorD fit_ellipse(TGraph *g)
TGraph * TestGraphDLSF(bool randomize=false)
TVectorD ConicToParametric(const TVectorD &conic)
void fitEllipseTGraphDLSF(TGraph *g=nullptr)
constexpr Double_t PiOver2()
constexpr Double_t TwoPi()