ROOT  6.07/01
Reference Guide
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Groups Pages
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,
42  Double_t a, Double_t b,
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 //
73 Double_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 //
84 TGraph *ellipse_TGraph = ((TGraph *)0);
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)
118  ellipse_TGraph = g;
119 
120  // create minimizer giving a name and (optionally) a name of the algorithm
121 #if 0 /* 0 or 1 */
123  ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
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;
154  ellipse_TGraph->ComputeRange(xmin, ymin, xmax, ymax);
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];
242  Double_t dt = TMath::TwoPi() / Double_t(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
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:432
void SetMaxIterations(unsigned int maxiter)
set maximum iterations (one iteration can have many function calls)
Definition: Minimizer.h:459
float xmin
Definition: THbookFile.cxx:93
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:3147
virtual Double_t Rndm(Int_t i=0)
Machine independent random number generator.
Definition: TRandom.cxx:512
Float_t theta
Definition: shapesAnim.C:5
return c
float ymin
Definition: THbookFile.cxx:93
Documentation for class Functor class.
Definition: Functor.h:394
#define gROOT
Definition: TROOT.h:344
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
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
TArc * a
Definition: textangle.C:12
const Bool_t kFALSE
Definition: Rtypes.h:92
virtual void SetFillStyle(Style_t fstyle)
Definition: TAttFill.h:52
Int_t GetN() const
Definition: TGraph.h:132
virtual void SetNameTitle(const char *name, const char *title)
Change (i.e. set) all the TNamed parameters (name and title).
Definition: TNamed.cxx:142
double cos(double)
TFile * f
int Status() const
status code of minimizer
Definition: Minimizer.h:438
virtual void Clear()
reset for consecutive minimizations - implement if needed
Definition: Minimizer.h:125
Double_t * GetY() const
Definition: TGraph.h:140
double sqrt(double)
Double_t x[n]
Definition: legend1.C:17
virtual void SetGrid(Int_t valuex=1, Int_t valuey=1)
Definition: TPad.h:326
int d
Definition: tornado.py:11
Abstract Minimizer class, defining the interface for the various minimizer (like Minuit2, Minuit, GSL, etc..) Plug-in's exist in ROOT to be able to instantiate the derived classes like ROOT::Math::GSLMinimizer or ROOT::Math::Minuit2Minimizer via the plug-in manager.
Definition: Minimizer.h:86
virtual double MinValue() const =0
return minimum function value
virtual bool Minimize()=0
method to perform the minimization
static Vc_ALWAYS_INLINE Vector< T > abs(const Vector< T > &x)
Definition: vector.h:450
void Clear(Option_t *option="")
Remove all primitives from the canvas.
Definition: TCanvas.cxx:678
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:640
virtual const double * X() const =0
return pointer to X values at the minimum
Double_t ellipse_TGraph_chi2(const Double_t *x)
Double_t TwoPi()
Definition: TMath.h:45
double sin(double)
TGraph * TestGraphRMM(Bool_t randomize=kFALSE)
virtual void Draw(Option_t *option="")
Draw this ellipse with its current attributes.
Definition: TEllipse.cxx:167
float ymax
Definition: THbookFile.cxx:93
virtual void SetFunction(const ROOT::Math::IMultiGenFunction &func)=0
set the function to minimize
Double_t * GetX() const
Definition: TGraph.h:139
ROOT::R::TRInterface & r
Definition: Object.C:4
SVector< double, 2 > v
Definition: Dict.h:5
TMarker * m
Definition: textangle.C:8
float xmax
Definition: THbookFile.cxx:93
A 2-Dim function with parameters.
Definition: TF2.h:33
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
Double_t Pi()
Definition: TMath.h:44
virtual void SetMaximum(Double_t maximum=-1111)
Set the maximum value along Y for this function In case the function is already drawn, set also the maximum in the helper histogram.
Definition: TF1.cxx:3077
The Canvas class.
Definition: TCanvas.h:48
void SetMaxFunctionCalls(unsigned int maxfcn)
set maximum of function calls
Definition: Minimizer.h:456
double Double_t
Definition: RtypesCore.h:55
void SetTolerance(double tol)
set the tolerance
Definition: Minimizer.h:462
Double_t y[n]
Definition: legend1.C:17
double atan(double)
Draw Ellipses.
Definition: TEllipse.h:44
void fitEllipseTGraphRMM(TGraph *g=((TGraph *) 0))
Double_t PiOver2()
Definition: TMath.h:46
virtual bool SetVariableLimits(unsigned int ivar, double lower, double upper)
set the limits of an already existing variable
Definition: Minimizer.h:217
virtual bool SetVariable(unsigned int ivar, const std::string &name, double val, double step)=0
set a new free variable
A Graph is a graphics object made of two arrays X and Y with npoints each.
Definition: TGraph.h:53
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)
void SetPrintLevel(int level)
set print level
Definition: Minimizer.h:453
virtual void Update()
Update canvas pad buffers.
Definition: TCanvas.cxx:2179
const Bool_t kTRUE
Definition: Rtypes.h:91
const Int_t n
Definition: legend1.C:16
void Modified(Bool_t flag=1)
Definition: TPad.h:407
ROOT::Math::Minimizer * ellipse_TGraph_minimize(TGraph *g)