Logo ROOT   6.10/09
Reference Guide
testUnbinGausFit.cxx
Go to the documentation of this file.
1 #include "TMath.h"
2 #include "TSystem.h"
3 #include "TRandom3.h"
4 #include "TTree.h"
5 #include "TROOT.h"
6 
7 #include "Fit/UnBinData.h"
8 #include "Fit/Fitter.h"
9 
10 
11 #include "Math/IParamFunction.h"
12 #include "Math/WrappedTF1.h"
13 #include "Math/WrappedMultiTF1.h"
16 
17 #include "TGraphErrors.h"
18 
19 #include "TStyle.h"
20 
21 
22 #include "Math/DistFunc.h"
23 
24 #include <string>
25 #include <iostream>
26 
27 #include "TStopwatch.h"
28 
29 #define DEBUG
30 
31 // test fit with many dimension
32 
33 const int N = 1; // 1d fit
34 const int NGaus = 3;
35 const int NPar = 8; // sum of 3 gaussians
36 const std::string branchType = "x[1]/D";
37 
38 // for 8 core testing use 1M points
39 const int NPoints = 100000;
40 double truePar[NPar];
41 double iniPar[NPar];
42 //const int nfit = 1;
43 const int strategy = 0;
44 
45 double gausnorm(const double *x, const double *p) {
46 
47  double invsig = 1./std::abs(p[1]);
48  double tmp = (x[0]-p[0]) * invsig;
49  const double sqrt_2pi = 1./std::sqrt(2.* 3.14159 );
50  return std::exp(-0.5 * tmp*tmp ) * sqrt_2pi * invsig;
51 }
52 
53 double gausSum(const double *x, const double *p) {
54 
55  double f = gausnorm(x,p+2) +
56  p[0] * gausnorm(x,p+4) +
57  p[1] * gausnorm(x,p+6);
58 
59  double norm = 1. + p[0] + p[1];
60  return f/norm;
61 }
62 
63 struct MINUIT2 {
64  static std::string name() { return "Minuit2"; }
65  static std::string name2() { return ""; }
66 };
67 
68 // fill fit data structure
70 
71  // fill the unbin data set from a TTree
72  ROOT::Fit::UnBinData * d = 0;
73 
74  // for the large tree access directly the branches
75  d = new ROOT::Fit::UnBinData();
76 
77  unsigned int n = tree->GetEntries();
78 #ifdef DEBUG
79  std::cout << "number of unbin data is " << n << " of dim " << N << std::endl;
80 #endif
81  d->Initialize(n,N);
82  TBranch * bx = tree->GetBranch("x");
83  double vx[N];
84  bx->SetAddress(vx);
85  std::vector<double> m(N);
86  for (int unsigned i = 0; i < n; ++i) {
87  bx->GetEntry(i);
88  d->Add(vx);
89  for (int j = 0; j < N; ++j)
90  m[j] += vx[j];
91  }
92 
93 #ifdef DEBUG
94  std::cout << "average values of means :\n";
95  for (int j = 0; j < N; ++j)
96  std::cout << m[j]/n << " ";
97  std::cout << "\n";
98 #endif
99 
100  return d;
101 }
102 
103 
104 // unbin fit
105 
107 template <class MinType, class T>
108 int DoUnBinFit(T * tree, Func & func, bool debug = false ) {
109 
110  ROOT::Fit::UnBinData * d = FillUnBinData(tree );
111  // need to have done Tree->Draw() before fit
112  //FillUnBinData(d,tree);
113 
114  std::cout << "Filled the fit data " << std::endl;
115  //printData(d);
116 
117 #ifdef DEBUG
118  std::cout << "data size type and size is " << typeid(*d).name() << " " << d->Size() << std::endl;
119 #endif
120 
121 
122 
123  // create the fitter
124  //std::cout << "Fit parameter 2 " << f.Parameters()[2] << std::endl;
125 
126  ROOT::Fit::Fitter fitter;
127  fitter.Config().SetMinimizer(MinType::name().c_str(),MinType::name2().c_str());
128 
129  if (debug)
130  fitter.Config().MinimizerOptions().SetPrintLevel(3);
131  else
132  fitter.Config().MinimizerOptions().SetPrintLevel(1);
133 
134 
135  // set tolerance 1 for tree to be same as in TTTreePlayer::UnBinFIt
136  fitter.Config().MinimizerOptions().SetTolerance(0.01);
137 
138  // set strategy (0 to avoid MnHesse
140 
141 
142  // create the function
143 
144  fitter.SetFunction(func);
145  // need to set limits to constant term
146  fitter.Config().ParSettings(0).SetLowerLimit(0.);
147  fitter.Config().ParSettings(1).SetLowerLimit(0.);
148 
149  if (debug)
150  std::cout << "do fitting... " << std::endl;
151 
152  bool ret = fitter.Fit(*d);
153  if (!ret) {
154  std::cout << " Fit Failed " << std::endl;
155  return -1;
156  }
157  if (debug)
158  fitter.Result().Print(std::cout);
159 
160  // check fit result
161  double chi2 = 0;
162  //if (fitter.Result().Value(0) < 0.5 ) {
163  for (int i = 0; i < NPar; ++i) {
164  double d = (truePar[i] - fitter.Result().Value(i) )/ (fitter.Result().Error(i) );
165  chi2 += d*d;
166  }
167 //}
168 // else {
169 // double truePar2[NPar];
170 // truePar2[0] = 1.-truePar[0];
171 // truePar2[1] = truePar[3];
172 // truePar2[2] = truePar[4];
173 // truePar2[3] = truePar[1];
174 // truePar2[4] = truePar[2];
175 // for (int i = 0; i < N; ++i) {
176 // double d = ( truePar2[i] - fitter.Result().Value(i) )/ (fitter.Result().Error(i) );
177 // chi2 += d*d;
178 // }
179 // }
180  double prob = ROOT::Math::chisquared_cdf_c(chi2,NPar);
181  int iret = (prob < 1.0E-6) ? -1 : 0;
182  if (iret != 0) {
183  std::cout <<"Found difference in fitted values - chi2 = " << chi2
184  << " prob = " << prob << std::endl;
185  fitter.Result().Print(std::cout);
186  }
187 
188  delete d;
189 
190  return iret;
191 
192 }
193 
194 
195 template <class MinType>
196 int DoFit(TTree * tree, Func & func, bool debug = false ) {
197  return DoUnBinFit<MinType, TTree>(tree, func, debug);
198 }
199 // template <class MinType>
200 // int DoFit(TH1 * h1, Func & func, bool debug = false, bool copyData = false ) {
201 // return DoBinFit<MinType, TH1>(h1, func, debug, copyData);
202 // }
203 // template <class MinType>
204 // int DoFit(TGraph * gr, Func & func, bool debug = false, bool copyData = false ) {
205 // return DoBinFit<MinType, TGraph>(gr, func, debug, copyData);
206 // }
207 
208 template <class MinType, class FitObj>
209 int FitUsingNewFitter(FitObj * fitobj, Func & func ) {
210 
211  std::cout << "\n************************************************************\n";
212  std::cout << "\tFit using new Fit::Fitter " << typeid(*fitobj).name() << std::endl;
213  std::cout << "\tMinimizer is " << MinType::name() << " " << MinType::name2() << " func dim = " << func.NDim() << std::endl;
214 
215  int iret = 0;
216  TStopwatch w; w.Start();
217 
218 #ifdef DEBUG
219  std::cout << "initial Parameters " << iniPar << " " << *iniPar << " " << *(iniPar+1) << std::endl;
220  func.SetParameters(iniPar);
221  iret |= DoFit<MinType>(fitobj,func,true );
222  if (iret != 0) {
223  std::cout << "Test failed " << std::endl;
224  }
225 
226 #else
227  for (int i = 0; i < nfit; ++i) {
228  func.SetParameters(iniPar);
229  iret = DoFit<MinType>(fitobj,func, false);
230  if (iret != 0) {
231  std::cout << "Test failed " << std::endl;
232  break;
233  }
234  }
235 #endif
236  w.Stop();
237  std::cout << "\nTime: \t" << w.RealTime() << " , " << w.CpuTime() << std::endl;
238  std::cout << "\n************************************************************\n";
239 
240  return iret;
241 }
242 
243 
244 int testNdimFit() {
245 
246 
247  std::cout << "\n\n************************************************************\n";
248  std::cout << "\t UNBINNED TREE (GAUSSIAN MULTI-DIM) FIT\n";
249  std::cout << "************************************************************\n";
250 
251  TTree t1("t2","a large Tree with gaussian variables");
252  double x[N];
253  Int_t ev;
254  t1.Branch("x",x,branchType.c_str());
255  t1.Branch("ev",&ev,"ev/I");
256 
257  // generate the true parameters
258 // for (int j = 0; j < NGaus; ++j) {
259 // double a = j+1;
260 // double mu = double(j)/NGaus;
261 // double s = 1.0 + double(j)/NGaus;
262 // truePar[3*j] = a;
263 // truePar[3*j+1] = mu;
264 // truePar[3*j+2] = s;
265 // tot += a;
266 // }
267  truePar[0] = 0.2; // % second gaussian
268  truePar[1] = 0.05; // % third gaussian ampl
269  truePar[2] = 0.; // mean first gaussian
270  truePar[3] = 0.5; // s1
271  truePar[4] = 0.; // mean secon gauss
272  truePar[5] = 1;
273  truePar[6] = -3; // mean third gaus
274  truePar[7] = 10;
275 
276 
277 
278  //fill the tree
279  TRandom3 r;
280  double norm = (1+truePar[0] + truePar[1] );
281  double a = 1./norm;
282  double b = truePar[0]/ norm;
283  double c = truePar[1]/ norm;
284  assert(a+b+c == 1.);
285  std::cout << " True amplitude gaussians " << a << " " << b << " " << c << std::endl;
286  for (Int_t i=0;i<NPoints;i++) {
287  for (int j = 0; j < N; ++j) {
288  if (r.Rndm() < a ) {
289  double mu = truePar[2];
290  double s = truePar[3];
291  x[j] = r.Gaus(mu,s);
292  }
293  else if (r.Rndm() < b ) {
294  double mu = truePar[4];
295  double s = truePar[5];
296  x[j] = r.Gaus(mu,s);
297  }
298  else {
299  double mu = truePar[6];
300  double s = truePar[7];
301  x[j] = r.Gaus(mu,s);
302  }
303  }
304 
305  ev = i;
306  t1.Fill();
307 
308  }
309  //t1.Draw("x"); // to select fit variable
310 
311  iniPar[0] = 0.5;
312  iniPar[1] = 0.05;
313  for (int i = 0; i <NGaus; ++i) {
314  iniPar[2*i+2] = 0 ;
315  iniPar[2*i+3] = 1. + 4*i;
316  std::cout << "inipar " << i << " = " << iniPar[2*i+2] << " " << iniPar[2*i+3] << std::endl;
317  }
318 
319  // use simply TF1 wrapper
320  //ROOT::Math::WrappedMultiTF1 f2(*f1);
322 
323  int iret = 0;
324  iret |= FitUsingNewFitter<MINUIT2>(&t1,f2);
325 
326  return iret;
327 }
328 
329 int main() {
330  return testNdimFit();
331 }
int FitUsingNewFitter(FitObj *fitobj, Func &func)
void SetPrintLevel(int level)
set print level
Random number generator class based on M.
Definition: TRandom3.h:27
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
virtual void SetAddress(void *add)
Set address of this branch.
Definition: TBranch.cxx:2148
double iniPar[NPar]
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
void SetTolerance(double tol)
set the tolerance
double truePar[NPar]
virtual Double_t Rndm()
Machine independent random number generator.
Definition: TRandom3.cxx:94
double Error(unsigned int i) const
parameter error by index
Definition: FitResult.h:187
virtual Double_t Gaus(Double_t mean=0, Double_t sigma=1)
Samples a random number from the standard Normal (Gaussian) Distribution with the given mean and sigm...
Definition: TRandom.cxx:235
static std::string name()
double T(double x)
Definition: ChebyshevPol.h:34
virtual Int_t Fill()
Fill all branches.
Definition: TTree.cxx:4383
double Value(unsigned int i) const
parameter value by index
Definition: FitResult.h:180
double gausnorm(const double *x, const double *p)
static std::string name2()
Double_t CpuTime()
Stop the stopwatch (if it is running) and return the cputime (in seconds) passed between the start an...
Definition: TStopwatch.cxx:125
Class describing the unbinned data sets (just x coordinates values) of any dimensions.
Definition: UnBinData.h:42
int Int_t
Definition: RtypesCore.h:41
TArc * a
Definition: textangle.C:12
const int NGaus
ROOT::Math::MinimizerOptions & MinimizerOptions()
access to the minimizer control parameter (non const method)
Definition: FitConfig.h:138
TLatex * t1
Definition: textangle.C:20
const ParameterSettings & ParSettings(unsigned int i) const
get the parameter settings for the i-th parameter (const method)
Definition: FitConfig.h:74
const FitResult & Result() const
get fit result
Definition: Fitter.h:337
double sqrt(double)
int main()
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
Double_t x[n]
Definition: legend1.C:17
double gausSum(const double *x, const double *p)
const int NPoints
int testNdimFit()
const FitConfig & Config() const
access to the fit configuration (const method)
Definition: Fitter.h:365
unsigned int Size() const
return number of fit points
Definition: FitData.h:302
void SetLowerLimit(double low)
set a single lower limit
virtual TBranch * GetBranch(const char *name)
Return pointer to the branch with the given name in this tree or its friends.
Definition: TTree.cxx:4979
bool Fit(const Data &data, const Function &func)
fit a data set using any generic model function If data set is binned a least square fit is performed...
Definition: Fitter.h:129
void SetMinimizer(const char *type, const char *algo=0)
set minimizer type
Definition: FitConfig.h:152
TRandom2 r(17)
ROOT::Math::IParamMultiFunction Func
Fitter class, entry point for performing all type of fits.
Definition: Fitter.h:77
TMarker * m
Definition: textangle.C:8
int nfit
Definition: testFitPerf.cxx:59
virtual Int_t GetEntry(Long64_t entry=0, Int_t getall=0)
Read all leaves of entry and return total number of bytes read.
Definition: TBranch.cxx:1291
void SetFunction(const IModelFunction &func, bool useGradient=false)
Set the fitted function (model function) from a parametric function interface.
Definition: Fitter.cxx:104
virtual void SetParameters(const double *p)=0
Set the parameter values.
double f(double x)
const int N
void Add(double x)
preallocate a data set given size and dimension of the coordinates if a vector already exists with co...
Definition: UnBinData.h:200
double func(double *x, double *p)
Definition: stressTF1.cxx:213
virtual Long64_t GetEntries() const
Definition: TTree.h:381
virtual Int_t Branch(TCollection *list, Int_t bufsize=32000, Int_t splitlevel=99, const char *name="")
Create one branch for each element in the collection.
Definition: TTree.cxx:1660
int DoUnBinFit(T *tree, Func &func, bool debug=false)
void Print(std::ostream &os, bool covmat=false) const
print the result and optionaly covariance matrix and correlations
Definition: FitResult.cxx:446
double chisquared_cdf_c(double x, double r, double x0=0)
Complement of the cumulative distribution function of the distribution with degrees of freedom (upp...
void SetStrategy(int stra)
set the strategy
double f2(const double *x)
ROOT::Fit::UnBinData * FillUnBinData(TTree *tree)
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
bool debug
Definition: tree.py:1
A TTree object has a header with a name and a title.
Definition: TTree.h:78
const std::string branchType
const int NPar
A TTree is a list of TBranches.
Definition: TBranch.h:57
double exp(double)
double norm(double *x, double *p)
Definition: unuranDistr.cxx:40
const Int_t n
Definition: legend1.C:16
WrappedParamFunction class to wrap any multi-dimensional function pbject implementing the operator()(...
virtual unsigned int NDim() const =0
Retrieve the dimension of the function.
int DoFit(TTree *tree, Func &func, bool debug=false)
Stopwatch class.
Definition: TStopwatch.h:28
const int strategy