ROOT   Reference Guide
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_t 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//
41 Double_t x0, Double_t y0,
43 Double_t theta) // (in degrees)
44{
45 Double_t 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_t ellipse_fcn(const Double_t *x, const Double_t *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//
93{
94 Double_t v = 0.0;
95 if (!ellipse_TGraph) return v; // just a precaution
96 for (Int_t i = 0; i < ellipse_TGraph->GetN(); i++) {
97 Double_t 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 0; // just a precaution
115 if (g->GetN() < 6) return 0; // 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)
153 Double_t 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_t i = 1; i < ellipse_TGraph->GetN(); i++) {
163 Double_t dx = (ellipse_TGraph->GetX())[i] - x0;
164 Double_t dy = (ellipse_TGraph->GetY())[i] - y0;
165 Double_t 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_t *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//
222 Int_t i;
223
224 // define the test ellipse
225 Double_t x0 = 4; // ellipse's "x" center
226 Double_t y0 = 3; // ellipse's "y" center
227 Double_t a = 2; // ellipse's "semimajor" axis along "x" (> 0)
228 Double_t b = 1; // ellipse's "semiminor" axis along "y" (> 0)
229 Double_t 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_t n = 100; // number of points
241 Double_t x[n], y[n];
243 Double_t tmp;
244 theta *= TMath::PiOver2() / 90.0; // degrees -> radians
245 for (i = 0; i < n; i++) {
246 x[i] = a * (std::cos(dt * Double_t(i)) + 0.1 * gRandom->Rndm() - 0.05);
247 y[i] = b * (std::sin(dt * Double_t(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//
270{
271 if (!g) g = TestGraphRMM(kTRUE); // 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_t *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
ROOT::R::TRInterface & r
Definition: Object.C:4
#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 e(i)
Definition: RSha256.hxx:103
int Int_t
Definition: RtypesCore.h:45
const Bool_t kFALSE
Definition: RtypesCore.h:92
bool Bool_t
Definition: RtypesCore.h:63
double Double_t
Definition: RtypesCore.h:59
const Bool_t kTRUE
Definition: RtypesCore.h:91
float xmin
Definition: THbookFile.cxx:95
float ymin
Definition: THbookFile.cxx:95
float xmax
Definition: THbookFile.cxx:95
float ymax
Definition: THbookFile.cxx:95
double cos(double)
double atan(double)
double sqrt(double)
double sin(double)
#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 corrisponding Minimizer given the string Supported Minimizers types are: ...
Definition: Factory.cxx:63
Documentation for class Functor class.
Definition: Functor.h:400
Abstract Minimizer class, defining the interface for the various minimizer (like Minuit2,...
Definition: Minimizer.h:75
The Canvas class.
Definition: TCanvas.h:23
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:3407
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:644
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:3475
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:132
Int_t GetN() const
Definition: TGraph.h:124
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:669
Double_t * GetX() const
Definition: TGraph.h:131
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom.cxx:552
TGraph * TestGraphRMM(Bool_t randomize=kFALSE)
Double_t ellipse_TGraph_chi2(const Double_t *x)
void fitEllipseTGraphRMM(TGraph *g=((TGraph *) 0))
TGraph * ellipse_TGraph
Double_t ellipse_fcn(Double_t x, Double_t y, Double_t x0, Double_t y0, Double_t a, Double_t b, Double_t theta)
ROOT::Math::Minimizer * ellipse_TGraph_minimize(TGraph *g)
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
auto * m
Definition: textangle.C:8
auto * a
Definition: textangle.C:12