ROOT   6.08/07 Reference Guide
GAMinTutorial.cxx
Go to the documentation of this file.
1 #include "TH1.h"
2 #include "TF1.h"
3 #include "Fit/BinData.h"
4 #include "Fit/Chi2FCN.h"
5
6 #include "Math/WrappedMultiTF1.h"
7 #include "Math/Minimizer.h"
9 #include "Math/Factory.h"
10 #include "HFitInterface.h"
11
12 #include "TMath.h"
13
14 #include "TApplication.h"
15 #include "TLegend.h"
16 #include "TCanvas.h"
17 #include "TStyle.h"
18 #include "TStopwatch.h"
19
20 //#define DEBUG
21
24  return par[0] + par[1]*x[0];
25 }
26
27 // Gaussian Peak function
29  return par[0]*TMath::Gaus(x[0],par[1], par[2]);
30 }
31
32 // Sum of background and peak function
34  return background(x,par) + gaussianPeak(x,&par[2]) + gaussianPeak(x,&par[5]);
35 }
36
37 // We'll look for the minimum at 2 and 7
38 double par0[8] = { 1, 0.05, 10 , 2, 0.5 , 10 , 7 , 1. };
39 const int ndata = 10000;
40 const double gAbsTolerance = 0.1;
41 int gVerbose = 0;
42 bool showGraphics = false;
43
44 using std::cout;
45 using std::endl;
46
47 int GAMinimize(ROOT::Math::IMultiGenFunction& chi2Func, double& xm1, double& xm2)
48 {
49  // minimize the function
51  if (min == 0) {
52  cout << "Error creating minimizer " << endl;
53  return -1;
54  }
55
56
57  // Set the parameters for the minimization.
58  min->SetMaxFunctionCalls(1000000);
59  min->SetMaxIterations(100000);
61  min->SetPrintLevel(gVerbose);
62  min->SetFunction(chi2Func);
63  ROOT::Math::GeneticMinimizerParameters params; // construct with default values
64  params.fNsteps = 100; // increset number of steps top 100 (default is 40)
65  params.fPopSize = 200; // increset number of steps top 100 (default is 40)
66  params.fSeed = 111; // use a fixed seed
67
68  min->SetParameters(params);
69
70
71
72  // initial values of the function
73  double x0[8];
74  std::copy(par0, par0 + 8, x0);
75  x0[3] = xm1;
76  x0[6] = xm2;
77
78  for (unsigned int i = 0; i < chi2Func.NDim(); ++i) {
79 #ifdef DEBUG
80  cout << "set variable " << i << " to value " << x0[i] << endl;
81 #endif
82  if ( i == 3 || i == 6 )
83  min->SetLimitedVariable(i,"x" + ROOT::Math::Util::ToString(i),x0[i], 0.1,2,8);
84  else
85  min->SetLimitedVariable(i,"x" + ROOT::Math::Util::ToString(i),x0[i], 0.1,x0[i]-2,x0[i]+2);
86  }
87
88
89  // minimize
90  if ( !min->Minimize() ) return 1;
91  double minval = min->MinValue();
92
93  // show the results
94  cout << "--------------------------------------" << endl;
95  cout << "GeneticMinimizer(" << xm1 << "," << xm2 << ") : ";
96  cout << "chi2 min value = " << minval;
97  cout << " minimum values m1 = " << min->X()[3] << " m2 = " << min->X()[6] << endl;
98  if (gVerbose) {
99  std::cout << "All Fit parameters : ";
100  for (unsigned int i = 0; i < chi2Func.NDim(); ++i)
101  cout << "x(" << i << ") = " << min->X()[i] << " ";
102  cout << endl;
103  }
104  xm1 = min->X()[3];
105  xm2 = min->X()[6];
106
107  return 0;
108 }
109
110
111 const double* Min2Minimize(ROOT::Math::IMultiGenFunction& chi2Func, double &xm1, double &xm2) {
112
113  // minimize the function
114  ROOT::Math::Minimizer * min = ROOT::Math::Factory::CreateMinimizer("Minuit2", "Migrad");
115  if (min == 0) {
117  if (min == 0) {
118  cout << "Error creating minimizer " << endl;
119  exit(-1);
120  }
121  }
122
123  // Set the minimizer options
124  min->SetMaxFunctionCalls(1000000);
125  min->SetMaxIterations(100000);
127  min->SetPrintLevel(gVerbose);
128  min->SetFunction(chi2Func);
129
130  // initial values of the function
131  double x0[8];
132  std::copy(par0, par0 + 8, x0);
133  x0[3] = xm1;
134  x0[6] = xm2;
135
136  for (unsigned int i = 0; i < chi2Func.NDim(); ++i) {
137 #ifdef DEBUG
138  cout << "set variable " << i << " to value " << x0[i] << endl;
139 #endif
140  min->SetVariable(i,"x" + ROOT::Math::Util::ToString(i),x0[i], 0.1);
141  }
142
143  // minimize
144  if ( !min->Minimize() ) exit(1);
145  double minval = min->MinValue();
146
147  // show the results
148  cout << "--------------------------------------" << endl;
149  cout << "Minuit2Minimizer(" << xm1 << "," << xm2 << ") : ";
150  cout << "chi2 min value = " << minval;
151  cout << " minimum values - m1 = " << min->X()[3] << " m2 = " << min->X()[6] << endl;
152  if (gVerbose) {
153  std::cout << "All Fit parameters : ";
154  for (unsigned int i = 0; i < chi2Func.NDim(); ++i)
155  cout << "x(" << i << ") = " << min->X()[i] << " ";
156  cout << endl;
157  }
158  xm1 = min->X()[3];
159  xm2 = min->X()[6];
160
161  return min->X();
162 }
163
165 {
166  double x1=0, x2=0;
167
168
169  // create a TF1 from fit function and generate histogram
170  TF1 * fitFunc = new TF1("fitFunc",fitFunction,0,10,8);
171  fitFunc->SetParameters(par0);
172
173  // Create a histogram filled with random data from TF1
174  TH1D * h1 = new TH1D("h1","",100, 0, 10);
175  for (int i = 0; i < ndata; ++i) {
176  h1->Fill(fitFunc->GetRandom() );
177  }
178
179  // perform the fit
181  ROOT::Fit::FillData(d,h1);
182  ROOT::Math::WrappedMultiTF1 f(*fitFunc);
183
184  // Create the function for fitting.
186
187  // Minimizer a first time with Minuit2
188  TStopwatch t;
189  t.Start();
190  const double* xmin1 = Min2Minimize(chi2Func, x1, x2);
191  t.Stop();
192  cout << "Minuit2 minimization Time :\t " << t.RealTime() << endl;
193
194
195  fitFunc->SetParameters(&xmin1[0]);
196  fitFunc->SetLineColor(kBlue+3);
197  TF1 * fitFunc0 = new TF1;
198  fitFunc->Copy(*fitFunc0);
199
200  if (showGraphics) {
201  h1->Draw();
203  h1->SetFillColor(kYellow-5);
204  gStyle->SetOptStat(0);
205  fitFunc0->Draw("same");
206  }
207
208
209
210  // Look for an approximation with a Genetic Algorithm
211  t.Start();
212  GAMinimize(chi2Func, x1,x2);
213  t.Stop();
214  cout << "Genetic minimization Time :\t " << t.RealTime() << endl;
215
216
217  // Minimize a second time with Minuit2
218  t.Start();
219  const double* xmin2 = Min2Minimize(chi2Func, x1, x2);
220  t.Stop();
221  cout << "Second Minuit2 minimization Time :\t " << t.RealTime() << endl;
222
223  fitFunc->SetParameters(&xmin2[0]);
224  if (showGraphics) {
225  fitFunc->SetLineColor(kRed+3);
226  fitFunc->DrawCopy("same");
227
228  TLegend* legend = new TLegend(0.61,0.72,0.86,0.86);
232  legend->Draw();
233
235  }
236
237  return 0;
238  // Min2Minimize will exit with a different value if there is any
239  // error.
240 }
241
242 int main(int argc, char **argv)
243 {
244  // Parse command line arguments
245  for (Int_t i=1 ; i<argc ; i++) {
246  std::string arg = argv[i] ;
247  if (arg == "-g") {
248  showGraphics = true;
249  }
250  if (arg == "-v") {
251  gVerbose = 1;
252  }
253  if (arg == "-vv") {
254  gVerbose = 3;
255  }
256  if (arg == "-h") {
257  std::cout << "Usage: " << argv[0] << " [-g] [-v]\n";
258  std::cout << " where:\n";
259  std::cout << " -g : graphics mode\n";
260  std::cout << " -v : verbose mode\n";
261  std::cout << " -vv : very verbose mode\n";
262  std::cout << std::endl;
263  return -1;
264  }
265  }
266  TApplication* theApp = 0;
267  if (showGraphics)
268  theApp = new TApplication("App",&argc,argv);
269
270  int status = GAMinTutorial();
271
272  if (showGraphics) {
273  theApp->Run();
274  delete theApp;
275  theApp = 0;
276  }
277
278  return status;
279 }
const int ndata
double par[1]
Definition: unuranDistr.cxx:38
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3125
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:439
virtual Double_t GetRandom()
Return a random number following this function shape.
Definition: TF1.cxx:1798
void SetMaxIterations(unsigned int maxiter)
set maximum iterations (one iteration can have many function calls)
Definition: Minimizer.h:459
virtual bool Minimize()
method to perform the minimization
Double_t RealTime()
Stop the stopwatch (if it is running) and return the realtime (in seconds) passed between the start a...
Definition: TStopwatch.cxx:110
This class displays a legend box (TPaveText) containing several legend entries.
Definition: TLegend.h:27
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
virtual void Copy(TObject &f1) const
Copy this F1 to a new F1.
Definition: TF1.cxx:770
int GAMinimize(ROOT::Math::IMultiGenFunction &chi2Func, double &xm1, double &xm2)
virtual TF1 * DrawCopy(Option_t *option="") const
Draw a copy of this function with its current attributes.
Definition: TF1.cxx:1114
Definition: Rtypes.h:61
Class to Wrap a ROOT Function class (like TF1) in a IParamMultiFunction interface of multi-dimensions...
virtual double MinValue() const
return minimum function value
R__EXTERN TStyle * gStyle
Definition: TStyle.h:418
virtual void Draw(Option_t *option="")
Draw this legend with its current attributes.
Definition: TLegend.cxx:373
Double_t background(Double_t *x, Double_t *par)
double fitFunc(double *x, double *p)
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
virtual void Draw(Option_t *option="")
Draw this function with its current attributes.
Definition: TF1.cxx:1086
int gVerbose
Definition: Rtypes.h:61
TLegend * legend
Definition: pirndm.C:35
const double gAbsTolerance
static const double x2[5]
Double_t fitFunction(Double_t *x, Double_t *par)
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
Double_t x[n]
Definition: legend1.C:17
virtual void Run(Bool_t retrn=kFALSE)
Main application eventloop. Calls system dependent eventloop via gSystem.
double par0[8]
Abstract Minimizer class, defining the interface for the various minimizer (like Minuit2, Minuit, GSL, etc..) Plug-in&#39;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
Double_t gaussianPeak(Double_t *x, Double_t *par)
virtual const double * X() const =0
return pointer to X values at the minimum
TH1F * h1
Definition: legend1.C:5
virtual void SetLineColor(Color_t lcolor)
Set the line color.
Definition: TAttLine.h:46
Chi2FCN class for binnned fits using the least square methods.
Definition: Chi2FCN.h:68
virtual void SetFunction(const ROOT::Math::IMultiGenFunction &func)=0
set the function to minimize
void FillData(BinData &dv, const TH1 *hist, TF1 *func=0)
fill the data vector from a TH1.
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2851
virtual void SetFillColor(Color_t fcolor)
Set the fill area color.
Definition: TAttFill.h:42
const double * Min2Minimize(ROOT::Math::IMultiGenFunction &chi2Func, double &xm1, double &xm2)
virtual unsigned int NDim() const =0
Retrieve the dimension of the function.
Class describing the binned data sets : vectors of x coordinates, y values and optionally error on y ...
Definition: BinData.h:61
void SetParameters(const GeneticMinimizerParameters &params)
tomato 1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:618
Double_t Gaus(Double_t x, Double_t mean=0, Double_t sigma=1, Bool_t norm=kFALSE)
Calculate a gaussian function with mean and sigma.
Definition: TMath.cxx:452
virtual const double * X() const
return pointer to X values at the minimum
void SetMaxFunctionCalls(unsigned int maxfcn)
set maximum of function calls
Definition: Minimizer.h:456
static const double x1[5]
double f(double x)
double Double_t
Definition: RtypesCore.h:55
TLegendEntry * AddEntry(const TObject *obj, const char *label="", Option_t *option="lpf")
Add a new entry to this legend.
Definition: TLegend.cxx:280
void SetTolerance(double tol)
set the tolerance
Definition: Minimizer.h:462
int GAMinTutorial()
virtual void SetFunction(const ROOT::Math::IMultiGenFunction &func)
set the function to minimize
int main(int argc, char **argv)
virtual bool SetVariable(unsigned int ivar, const std::string &name, double val, double step)=0
set a new free variable
std::string ToString(const T &val)
Utility function for conversion to strings.
Definition: Util.h:42
1-Dim function class
Definition: TF1.h:149
void SetOptStat(Int_t stat=1)
The type of information printed in the histogram statistics box can be selected via the parameter mod...
Definition: TStyle.cxx:1257
THist< 1, double, THistStatContent, THistStatUncertainty > TH1D
Definition: THist.hxx:301