ROOT   Reference Guide
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// http://root.cern.ch/drupal/content/numerical-minimization#multidim_minim
110// http://root.cern.ch/root/html534/tutorials/fit/NumericalMinimization.C.html
111//
113{
114 if (!g) return nullptr; // just a precaution
115 if (g->GetN() < 6) return nullptr; // just a precaution
116
117 // set the TGraph to be fitted (used by the ellipse_TGraph_chi2 function)
119
120 // create minimizer giving a name and (optionally) a name of the algorithm
121#if 0 /* 0 or 1 */
124#elif 0 /* 0 or 1 */
126 ROOT::Math::Factory::CreateMinimizer("Minuit2", "Simplex");
127#elif 1 /* 0 or 1 */
129 ROOT::Math::Factory::CreateMinimizer("Minuit2", "Combined");
130#else /* 0 or 1 */
132 ROOT::Math::Factory::CreateMinimizer("Minuit2", "Scan");
133#endif /* 0 or 1 */
134
135 // set tolerance, etc. ...
136 m->SetMaxFunctionCalls(1000000); // for Minuit and Minuit2
137 m->SetMaxIterations(100000); // for GSL
138 m->SetTolerance(0.001); // edm
139#if 1 /* 0 or 1 */
140 m->SetPrintLevel(1);
141#endif /* 0 or 1 */
142
143 // create function wrapper for minimizer (a IMultiGenFunction type)
145
146 m->SetFunction(f);
147
148 m->Clear(); // just a precaution
149
150 // estimate all initial values (note: good initial values
151 // are CRUCIAL for the minimizing procedure to succeed)
152 double xmin, xmax, ymin, ymax;
153 double x0, y0, a, b, theta;
155 x0 = (xmax + xmin) / 2.0;
156 y0 = (ymax + ymin) / 2.0;
157 a = (ellipse_TGraph->GetX())[0] - x0;
158 b = (ellipse_TGraph->GetY())[0] - y0;
159 theta = ((std::abs(b) > 9999.9 * std::abs(a)) ? 9999.9 : (b / a));
160 a = a * a + b * b;
161 b = a;
162 for (int i = 1; i < ellipse_TGraph->GetN(); i++) {
163 double dx = (ellipse_TGraph->GetX())[i] - x0;
164 double dy = (ellipse_TGraph->GetY())[i] - y0;
165 double d = dx * dx + dy * dy;
166 // try to keep "a" > "b"
167 if (a < d) {
168 a = d;
169 theta = ((std::abs(dy) > 9999.9 * std::abs(dx)) ? 9999.9 : (dy / dx));
170 }
171 if (b > d) b = d;
172 }
173 a = std::sqrt(a); if (!(a > 0)) a = 0.001;
174 b = std::sqrt(b); if (!(b > 0)) b = 0.001;
175 theta = std::atan(theta) * 180.0 / TMath::Pi();
176 if (theta < -45.0) theta += 180.0; // "theta" = -45 ... 135 degrees
177
178 // set the variables to be minimized
179 m->SetVariable(0, "x0", x0, (xmax - xmin) / 100.0);
180 m->SetVariable(1, "y0", y0, (ymax - ymin) / 100.0);
181 m->SetVariable(2, "a", a, a / 100.0);
182 m->SetVariable(3, "b", b, b / 100.0);
183 m->SetVariable(4, "theta", theta, 1);
184
185#if 1 /* 0 or 1 */
186 // set the variables' limits
187 m->SetVariableLimits(0, xmin, xmax);
188 m->SetVariableLimits(1, ymin, ymax);
189 if (theta < 45.0) {
190 if (a < ((xmax - xmin) / 2.0)) a = (xmax - xmin) / 2.0;
191 if (b < ((ymax - ymin) / 2.0)) b = (ymax - ymin) / 2.0;
192 } else {
193 if (a < ((ymax - ymin) / 2.0)) a = (ymax - ymin) / 2.0;
194 if (b < ((xmax - xmin) / 2.0)) b = (xmax - xmin) / 2.0;
195 }
196 m->SetVariableLimits(2, 0, a * 3.0);
197 m->SetVariableLimits(3, 0, b * 3.0);
198 // m->SetVariableLimits(4, theta - 30.0, theta + 30.0); // theta -+ 30
199 m->SetVariableLimits(4, theta - 45.0, theta + 45.0); // theta -+ 45
200#endif /* 0 or 1 */
201
202 // do the minimization
203 m->Minimize();
204
205#if 0 /* 0 or 1 */
206 const double *xm = m->X();
207 std::cout << "Minimum ( "
208 << xm[0] << " , " << xm[1] << " , " // "x0", "y0"
209 << xm[2] << " , " << xm[3] << " , " // "a", "b"
210 << xm[4] << " ) = " // "theta" (in degrees)
211 << m->MinValue() // it's equal to ellipse_TGraph_chi2(xm)
212 << std::endl;
213#endif /* 0 or 1 */
214
215 return m;
216}
217
218//
219// creates a test TGraph with an ellipse
220//
221TGraph *TestGraphRMM(bool randomize = false) {
222 int i;
223
224 // define the test ellipse
225 double x0 = 4; // ellipse's "x" center
226 double y0 = 3; // ellipse's "y" center
227 double a = 2; // ellipse's "semimajor" axis along "x" (> 0)
228 double b = 1; // ellipse's "semiminor" axis along "y" (> 0)
229 double theta = 100; // ellipse's axes rotation angle (-45 ... 135 degrees)
230
231 // gRandom->SetSeed(0);
232 if (randomize) {
233 x0 = 10.0 - 20.0 * gRandom->Rndm();
234 y0 = 10.0 - 20.0 * gRandom->Rndm();
235 a = 0.5 + 4.5 * gRandom->Rndm();
236 b = 0.5 + 4.5 * gRandom->Rndm();
237 theta = 180.0 - 360.0 * gRandom->Rndm();
238 }
239
240 const int n = 100; // number of points
241 double x[n], y[n];
242 double dt = TMath::TwoPi() / double(n);
243 double tmp;
244 theta *= TMath::PiOver2() / 90.0; // degrees -> radians
245 for (i = 0; i < n; i++) {
246 x[i] = a * (std::cos(dt * double(i)) + 0.1 * gRandom->Rndm() - 0.05);
247 y[i] = b * (std::sin(dt * double(i)) + 0.1 * gRandom->Rndm() - 0.05);
248 // rotate the axes
249 tmp = x[i];
250 x[i] = x[i] * std::cos(theta) - y[i] * std::sin(theta);
251 y[i] = y[i] * std::cos(theta) + tmp * std::sin(theta);
252 // shift the center
253 x[i] += x0;
254 y[i] += y0;
255 }
256
257 // create the test TGraph
258 TGraph *g = ((TGraph *)(gROOT->FindObject("g")));
259 if (g) delete g;
260 g = new TGraph(n, x, y);
261 g->SetNameTitle("g", "test ellipse");
262
263 return g;
264}
265
266//
267// "ROOT Script" entry point (the same name as the "filename's base")
268//
269void fitEllipseTGraphRMM(TGraph *g = ((TGraph *)nullptr))
270{
271 if (!g) g = TestGraphRMM(true); // create a "random" ellipse
272
273#if 0 /* 0 or 1 */
274 // create the "ellipse" TF2 (just for fun)
275 TF2 *ellipse = ((TF2 *)(gROOT->GetListOfFunctions()->FindObject("ellipse")));
276 if (ellipse) delete ellipse;
277 ellipse = new TF2("ellipse", ellipse_fcn, -1, 1, -1, 1, 5);
278 ellipse->SetMaximum(2.0); // just for nice graphics
279 ellipse->SetParNames("x0", "y0", "a", "b", "theta");
280 ellipse->SetParameters(0.4, 0.3, 0.2, 0.1, 10);
281#endif /* 0 or 1 */
282
283 // fit the TGraph
285
286#if 1 /* 0 or 1 */
287 // draw everything
288 TCanvas *c = ((TCanvas *)(gROOT->GetListOfCanvases()->FindObject("c")));
289 if (c) { c->Clear(); } else { c = new TCanvas("c", "c"); }
290 c->SetGrid(1, 1);
291 g->Draw("A*");
292 if ( m && (!(m->Status())) ) {
293 const double *xm = m->X();
294 TEllipse *e = new TEllipse(xm[0], xm[1], // "x0", "y0"
295 xm[2], xm[3], // "a", "b"
296 0, 360,
297 xm[4]); // "theta" (in degrees)
298 e->SetFillStyle(0); // hollow
299 e->Draw();
300 }
301 c->Modified(); c->Update(); // make sure it's really drawn
302#endif /* 0 or 1 */
303
304 delete m; // "cleanup"
305
306 return;
307}
308
309// 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 c(i)
Definition RSha256.hxx:101
#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:407
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:117
The Canvas class.
Definition TCanvas.h:23
void Clear(Option_t *option="") override
Remove all primitives from the canvas.
Definition TCanvas.cxx:734
Draw Ellipses.
Definition TEllipse.h:23
virtual void SetMaximum(Double_t maximum=-1111)
Set the maximum value along Y for this function In case the function is already drawn,...
Definition TF1.cxx:3394
virtual void SetParameters(const Double_t *params)
Definition TF1.h:650
virtual void SetParNames(const char *name0="p0", const char *name1="p1", const char *name2="p2", const char *name3="p3", const char *name4="p4", const char *name5="p5", const char *name6="p6", const char *name7="p7", const char *name8="p8", const char *name9="p9", const char *name10="p10")
Set up to 10 parameter names.
Definition TF1.cxx:3462
A 2-Dim function with parameters.
Definition TF2.h:29
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:138
Int_t GetN() const
Definition TGraph.h:130
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:709
Double_t * GetX() const
Definition TGraph.h:137
virtual void Clear(Option_t *="")
Definition TObject.h:119
Double_t Rndm() override
Machine independent random number generator.
Definition TRandom.cxx:552
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
void fitEllipseTGraphRMM(TGraph *g=((TGraph *) nullptr))
ROOT::Math::Minimizer * ellipse_TGraph_minimize(TGraph *g)
double ellipse_TGraph_chi2(const double *x)
Double_t y[n]
Definition legend1.C:17
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