Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
fitEllipseTGraphRMM.cxx
Go to the documentation of this file.
1//
2// The "fitEllipseTGraphRMM" macro uses the "ROOT::Math::Minimizer"
3// interface for fitting an ellipse to a set of data points from a TGraph
4//
5// To try this macro, in a ROOT prompt, do:
6// .L fitEllipseTGraphRMM.cxx // or ".L fitEllipseTGraphRMM.cxx++"
7// fitEllipseTGraphRMM(TestGraphRMM());
8// for (int i=0; i<10; i++) { fitEllipseTGraphRMM(); gSystem->Sleep(333); }
9//
10// Last update: Thu Jul 31 18:00:00 UTC 2014
11//
12// Changes:
13// 2014.07.31 - (initial version)
14//
15
16#include "TROOT.h"
17#include "TMath.h"
18#include "TRandom.h"
19#include "TGraph.h"
20#include "TF2.h"
21#include "TCanvas.h"
22#include "TEllipse.h"
23#include "Math/Minimizer.h"
24#include "Math/Factory.h"
25#include "Math/Functor.h"
26
27#include <cmath>
28#include <iostream>
29
30//
31// ellipse_fcn calculates the "normalized distance" from the ellipse's center
32//
33// x, y = point for which one wants to calculate the "normalized distance"
34// x0 = ellipse's "x" center
35// y0 = ellipse's "y" center
36// a = ellipse's "semimajor" axis along "x" (> 0)
37// b = ellipse's "semiminor" axis along "y" (> 0)
38// theta = ellipse's axes rotation angle (-45 ... 135 degrees)
39//
40double ellipse_fcn(double x, double y,
41 double x0, double y0,
42 double a, double b,
43 double theta) // (in degrees)
44{
45 double v = 9999.9;
46 if ((a == 0.0) || (b == 0.0)) return v; // just a precaution
47 // shift the center
48 x -= x0;
49 y -= y0;
50 // un-rotate the axes
51 theta *= TMath::Pi() / 180.0; // degrees -> radians
52 v = x;
53 x = x * std::cos(theta) + y * std::sin(theta);
54 y = y * std::cos(theta) - v * std::sin(theta);
55 // "scale" axes
56 x /= a;
57 y /= b;
58 // calculate the "normalized distance"
59 v = x * x + y * y;
60 v = std::sqrt(v);
61 return v;
62}
63
64//
65// x[0] = "x"
66// x[1] = "y"
67// params[0] = ellipse's "x" center ("x0")
68// params[1] = ellipse's "y" center ("y0")
69// params[2] = ellipse's "semimajor" axis along "x" ("a" > 0)
70// params[3] = ellipse's "semiminor" axis along "y" ("b" > 0)
71// params[4] = ellipse's axes rotation angle ("theta" = -45 ... 135 degrees)
72//
73double ellipse_fcn(const double *x, const double *params)
74{
75 return ellipse_fcn(x[0], x[1], // "x", "y"
76 params[0], params[1], // "x0", "y0"
77 params[2], params[3], // "a", "b"
78 params[4]); // "theta" (in degrees)
79}
80
81//
82// the TGraph to be fitted (used by the ellipse_TGraph_chi2 function below)
83//
85//
86// x[0] = ellipse's "x" center ("x0")
87// x[1] = ellipse's "y" center ("y0")
88// x[2] = ellipse's "semimajor" axis along "x" ("a" > 0)
89// x[3] = ellipse's "semiminor" axis along "y" ("b" > 0)
90// x[4] = ellipse's axes rotation angle ("theta" = -45 ... 135 degrees)
91//
92double ellipse_TGraph_chi2(const double *x)
93{
94 double v = 0.0;
95 if (!ellipse_TGraph) return v; // just a precaution
96 for (int i = 0; i < ellipse_TGraph->GetN(); i++) {
97 double r = ellipse_fcn((ellipse_TGraph->GetX())[i], // "x"
98 (ellipse_TGraph->GetY())[i], // "y"
99 x[0], x[1], // "x0", "y0"
100 x[2], x[3], // "a", "b"
101 x[4]); // "theta" (in degrees)
102 r -= 1.0; // ellipse's "radius" in "normalized coordinates" is always 1
103 v += r * r;
104 }
105 return v;
106}
107
108///
109/// https://root.cern.ch/doc/master/NumericalMinimization_8C.html
110///
112{
113 if (!g) return nullptr; // just a precaution
114 if (g->GetN() < 6) return nullptr; // just a precaution
115
116 // set the TGraph to be fitted (used by the ellipse_TGraph_chi2 function)
118
119 // create minimizer giving a name and (optionally) a name of the algorithm
120#if 0 /* 0 or 1 */
122 ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
123#elif 0 /* 0 or 1 */
125 ROOT::Math::Factory::CreateMinimizer("Minuit2", "Simplex");
126#elif 1 /* 0 or 1 */
128 ROOT::Math::Factory::CreateMinimizer("Minuit2", "Combined");
129#else /* 0 or 1 */
131 ROOT::Math::Factory::CreateMinimizer("Minuit2", "Scan");
132#endif /* 0 or 1 */
133
134 // set tolerance, etc. ...
135 m->SetMaxFunctionCalls(1000000); // for Minuit and Minuit2
136 m->SetMaxIterations(100000); // for GSL
137 m->SetTolerance(0.001); // edm
138#if 1 /* 0 or 1 */
139 m->SetPrintLevel(1);
140#endif /* 0 or 1 */
141
142 // create function wrapper for minimizer (a IMultiGenFunction type)
144
145 m->SetFunction(f);
146
147 m->Clear(); // just a precaution
148
149 // estimate all initial values (note: good initial values
150 // are CRUCIAL for the minimizing procedure to succeed)
151 double xmin, xmax, ymin, ymax;
152 double x0, y0, a, b, theta;
154 x0 = (xmax + xmin) / 2.0;
155 y0 = (ymax + ymin) / 2.0;
156 a = (ellipse_TGraph->GetX())[0] - x0;
157 b = (ellipse_TGraph->GetY())[0] - y0;
158 theta = ((std::abs(b) > 9999.9 * std::abs(a)) ? 9999.9 : (b / a));
159 a = a * a + b * b;
160 b = a;
161 for (int i = 1; i < ellipse_TGraph->GetN(); i++) {
162 double dx = (ellipse_TGraph->GetX())[i] - x0;
163 double dy = (ellipse_TGraph->GetY())[i] - y0;
164 double d = dx * dx + dy * dy;
165 // try to keep "a" > "b"
166 if (a < d) {
167 a = d;
168 theta = ((std::abs(dy) > 9999.9 * std::abs(dx)) ? 9999.9 : (dy / dx));
169 }
170 if (b > d) b = d;
171 }
172 a = std::sqrt(a); if (!(a > 0)) a = 0.001;
173 b = std::sqrt(b); if (!(b > 0)) b = 0.001;
174 theta = std::atan(theta) * 180.0 / TMath::Pi();
175 if (theta < -45.0) theta += 180.0; // "theta" = -45 ... 135 degrees
176
177 // set the variables to be minimized
178 m->SetVariable(0, "x0", x0, (xmax - xmin) / 100.0);
179 m->SetVariable(1, "y0", y0, (ymax - ymin) / 100.0);
180 m->SetVariable(2, "a", a, a / 100.0);
181 m->SetVariable(3, "b", b, b / 100.0);
182 m->SetVariable(4, "theta", theta, 1);
183
184#if 1 /* 0 or 1 */
185 // set the variables' limits
186 m->SetVariableLimits(0, xmin, xmax);
187 m->SetVariableLimits(1, ymin, ymax);
188 if (theta < 45.0) {
189 if (a < ((xmax - xmin) / 2.0)) a = (xmax - xmin) / 2.0;
190 if (b < ((ymax - ymin) / 2.0)) b = (ymax - ymin) / 2.0;
191 } else {
192 if (a < ((ymax - ymin) / 2.0)) a = (ymax - ymin) / 2.0;
193 if (b < ((xmax - xmin) / 2.0)) b = (xmax - xmin) / 2.0;
194 }
195 m->SetVariableLimits(2, 0, a * 3.0);
196 m->SetVariableLimits(3, 0, b * 3.0);
197 // m->SetVariableLimits(4, theta - 30.0, theta + 30.0); // theta -+ 30
198 m->SetVariableLimits(4, theta - 45.0, theta + 45.0); // theta -+ 45
199#endif /* 0 or 1 */
200
201 // do the minimization
202 m->Minimize();
203
204#if 0 /* 0 or 1 */
205 const double *xm = m->X();
206 std::cout << "Minimum ( "
207 << xm[0] << " , " << xm[1] << " , " // "x0", "y0"
208 << xm[2] << " , " << xm[3] << " , " // "a", "b"
209 << xm[4] << " ) = " // "theta" (in degrees)
210 << m->MinValue() // it's equal to ellipse_TGraph_chi2(xm)
211 << std::endl;
212#endif /* 0 or 1 */
213
214 return m;
215}
216
217//
218// creates a test TGraph with an ellipse
219//
220TGraph *TestGraphRMM(bool randomize = false)
221{
222 // define the test ellipse
223 double x0 = 4; // ellipse's "x" center
224 double y0 = 3; // ellipse's "y" center
225 double a = 2; // ellipse's "semimajor" axis along "x" (> 0)
226 double b = 1; // ellipse's "semiminor" axis along "y" (> 0)
227 double theta = 100; // ellipse's axes rotation angle (-45 ... 135 degrees)
228
229 // gRandom->SetSeed(0);
230 if (randomize) {
231 x0 = 10.0 - 20.0 * gRandom->Rndm();
232 y0 = 10.0 - 20.0 * gRandom->Rndm();
233 a = 0.5 + 4.5 * gRandom->Rndm();
234 b = 0.5 + 4.5 * gRandom->Rndm();
235 theta = 180.0 - 360.0 * gRandom->Rndm();
236 }
237
238 const int n = 100; // number of points
239 double x[n], y[n];
240 double dt = TMath::TwoPi() / double(n);
241 double tmp;
242 theta *= TMath::PiOver2() / 90.0; // degrees -> radians
243 for (int i = 0; i < n; i++) {
244 x[i] = a * (std::cos(dt * double(i)) + 0.1 * gRandom->Rndm() - 0.05);
245 y[i] = b * (std::sin(dt * double(i)) + 0.1 * gRandom->Rndm() - 0.05);
246 // rotate the axes
247 tmp = x[i];
248 x[i] = x[i] * std::cos(theta) - y[i] * std::sin(theta);
249 y[i] = y[i] * std::cos(theta) + tmp * std::sin(theta);
250 // shift the center
251 x[i] += x0;
252 y[i] += y0;
253 }
254
255 // create the test TGraph
256 TGraph *g = static_cast<TGraph *>(gROOT->FindObject("g"));
257 if (g) delete g;
258 g = new TGraph(n, x, y);
259 g->SetNameTitle("g", "test ellipse");
260
261 return g;
262}
263
264//
265// "ROOT Script" entry point (the same name as the "filename's base")
266//
268{
269 if (!g) g = TestGraphRMM(true); // create a "random" ellipse
270
271 // fit the TGraph
273
274 // draw everything
275 auto c1 = new TCanvas("c1","c1", 1000, 800);
276 c1->SetGrid(1, 1);
277 g->Draw("A*");
278 if ( m && (!(m->Status())) ) {
279 const double *xm = m->X();
280 TEllipse *e = new TEllipse(xm[0], xm[1], // "x0", "y0"
281 xm[2], xm[3], // "a", "b"
282 0, 360,
283 xm[4]); // "theta" (in degrees)
284 e->SetFillStyle(0); // hollow
285 e->Draw();
286 }
287 c1->Modified();
288 c1->Update(); // make sure it's really drawn
289
290 delete m; // "cleanup"
291}
292
293// end of file fitEllipseTGraphRMM.cxx by Silesius Anonymus
#define d(i)
Definition RSha256.hxx:102
#define b(i)
Definition RSha256.hxx:100
#define f(i)
Definition RSha256.hxx:104
#define g(i)
Definition RSha256.hxx:105
#define a(i)
Definition RSha256.hxx:99
#define e(i)
Definition RSha256.hxx:103
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
float xmin
float ymin
float xmax
float ymax
#define gROOT
Definition TROOT.h:406
R__EXTERN TRandom * gRandom
Definition TRandom.h:62
static ROOT::Math::Minimizer * CreateMinimizer(const std::string &minimizerType="", const std::string &algoType="")
static method to create the corresponding Minimizer given the string Supported Minimizers types are: ...
Definition Factory.cxx:63
Documentation for class Functor class.
Definition Functor.h:47
Abstract Minimizer class, defining the interface for the various minimizer (like Minuit2,...
Definition Minimizer.h:119
The Canvas class.
Definition TCanvas.h:23
Draw Ellipses.
Definition TEllipse.h:23
A TGraph is an object made of two arrays X and Y with npoints each.
Definition TGraph.h:41
Double_t * GetY() const
Definition TGraph.h:140
Int_t GetN() const
Definition TGraph.h:132
virtual void ComputeRange(Double_t &xmin, Double_t &ymin, Double_t &xmax, Double_t &ymax) const
Compute the x/y range of the points in this graph.
Definition TGraph.cxx:733
Double_t * GetX() const
Definition TGraph.h:139
virtual void Clear(Option_t *="")
Definition TObject.h:119
Double_t Rndm() override
Machine independent random number generator.
Definition TRandom.cxx:559
void fitEllipseTGraphRMM(TGraph *g=nullptr)
TGraph * TestGraphRMM(bool randomize=false)
double ellipse_fcn(double x, double y, double x0, double y0, double a, double b, double theta)
TGraph * ellipse_TGraph
ROOT::Math::Minimizer * ellipse_TGraph_minimize(TGraph *g)
https://root.cern.ch/doc/master/NumericalMinimization_8C.html
double ellipse_TGraph_chi2(const double *x)
Double_t y[n]
Definition legend1.C:17
return c1
Definition legend1.C:41
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
constexpr Double_t PiOver2()
Definition TMath.h:51
constexpr Double_t Pi()
Definition TMath.h:37
constexpr Double_t TwoPi()
Definition TMath.h:44
TMarker m
Definition textangle.C:8