ROOT   6.10/09 Reference Guide
testNdimFit.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 = 10;
34 const std::string branchType = "x[10]/D";
35 const int NPoints = 100000;
36
37 // const int N = 50;
38 // const std::string branchType = "x[50]/D";
39 // const int NPoints = 10000;
40
41
42
43 double truePar[2*N];
44 double iniPar[2*N];
45 //const int nfit = 1;
46 const int strategy = 0;
47
48 double gausnorm(const double *x, const double *p) {
49
50  double invsig = 1./p[1];
51  double tmp = (x[0]-p[0]) * invsig;
52  const double sqrt_2pi = 1./std::sqrt(2.* 3.14159 );
53  return std::exp(-0.5 * tmp*tmp ) * sqrt_2pi * invsig;
54 }
55
56 double gausnormN(const double *x, const double *p) {
57  double f = 1;
58  for (int i = 0; i < N; ++i)
59  f *= gausnorm(x+i,p+2*i);
60
61  return f;
62 }
63
64 struct MINUIT2 {
65  static std::string name() { return "Minuit2"; }
66  static std::string name2() { return ""; }
67 };
68
69 // fill fit data structure
71
72  // fill the unbin data set from a TTree
73  ROOT::Fit::UnBinData * d = 0;
74
75  // for the large tree access directly the branches
76  d = new ROOT::Fit::UnBinData();
77
78  unsigned int n = tree->GetEntries();
79 #ifdef DEBUG
80  std::cout << "number of unbin data is " << n << " of dim " << N << std::endl;
81 #endif
82  d->Initialize(n,N);
83  TBranch * bx = tree->GetBranch("x");
84  double vx[N];
86  std::vector<double> m(N);
87  for (int unsigned i = 0; i < n; ++i) {
88  bx->GetEntry(i);
90  for (int j = 0; j < N; ++j)
91  m[j] += vx[j];
92  }
93
94 #ifdef DEBUG
95  std::cout << "average values of means :\n";
96  for (int j = 0; j < N; ++j)
97  std::cout << m[j]/n << " ";
98  std::cout << "\n";
99 #endif
100
101  delete tree;
102  tree = 0;
103  return d;
104
105 }
106
107
108 // unbin fit
109
111 template <class MinType, class T>
112 int DoUnBinFit(T * tree, Func & func, bool debug = false ) {
113
114  ROOT::Fit::UnBinData * d = FillUnBinData(tree );
115  // need to have done Tree->Draw() before fit
116  //FillUnBinData(d,tree);
117
118  //std::cout << "data size type and size is " << typeid(*d).name() << " " << d->Size() << std::endl;
119  std::cout << "Fit data size = " << d->Size() << " dimension = " << d->NDim() << std::endl;
120
121
122
123  //printData(d);
124
125  // create the fitter
126  //std::cout << "Fit parameter 2 " << f.Parameters()[2] << std::endl;
127
128  ROOT::Fit::Fitter fitter;
129  fitter.Config().SetMinimizer(MinType::name().c_str(),MinType::name2().c_str());
130
131  if (debug)
132  fitter.Config().MinimizerOptions().SetPrintLevel(3);
133  else
134  fitter.Config().MinimizerOptions().SetPrintLevel(0);
135
136  // set tolerance 1 for tree to be same as in TTTreePlayer::UnBinFIt
137  fitter.Config().MinimizerOptions().SetTolerance(1);
138
139  // set strategy (0 to avoid MnHesse
141
142
143  // create the function
144
145  fitter.SetFunction(func);
146  // need to fix param 0 , normalization in the unbinned fits
147  //fitter.Config().ParSettings(0).Fix();
148
149  bool ret = fitter.Fit(*d);
150  if (!ret) {
151  std::cout << " Fit Failed " << std::endl;
152  return -1;
153  }
154  if (debug)
155  fitter.Result().Print(std::cout);
156
157  // check fit result
158  double chi2 = 0;
159  for (int i = 0; i < N; ++i) {
160  double d = (truePar[i] - fitter.Result().Value(i) )/ (fitter.Result().Error(i) );
161  chi2 += d*d;
162  }
163  double prob = ROOT::Math::chisquared_cdf_c(chi2,N);
164  int iret = (prob < 1.0E-6) ? -1 : 0;
165  if (iret != 0) {
166  std::cout <<"Found difference in fitted values - prob = " << prob << std::endl;
167  if (!debug) fitter.Result().Print(std::cout);
168  for (int i = 0; i < N; ++i) {
169  double d = (truePar[i] - fitter.Result().Value(i) )/ (fitter.Result().Error(i) );
170  std::cout << "par_" << i << " = " << fitter.Result().Value(i) << " true = " << truePar[i] << " pull = " << d << std::endl;
171  }
172
173  }
174
175  delete d;
176
177  return iret;
178
179 }
180
181
182 template <class MinType>
183 int DoFit(TTree * tree, Func & func, bool debug = false ) {
184  return DoUnBinFit<MinType, TTree>(tree, func, debug);
185 }
186 // template <class MinType>
187 // int DoFit(TH1 * h1, Func & func, bool debug = false, bool copyData = false ) {
188 // return DoBinFit<MinType, TH1>(h1, func, debug, copyData);
189 // }
190 // template <class MinType>
191 // int DoFit(TGraph * gr, Func & func, bool debug = false, bool copyData = false ) {
192 // return DoBinFit<MinType, TGraph>(gr, func, debug, copyData);
193 // }
194
195 template <class MinType, class FitObj>
196 int FitUsingNewFitter(FitObj * fitobj, Func & func ) {
197
198  std::cout << "\n************************************************************\n";
199  std::cout << "\tFit using new Fit::Fitter " << typeid(*fitobj).name() << std::endl;
200  std::cout << "\tMinimizer is " << MinType::name() << " " << MinType::name2() << " func dim = " << func.NDim() << std::endl;
201
202  int iret = 0;
203  TStopwatch w; w.Start();
204
205 #ifdef DEBUG
206  std::cout << "initial Parameters " << iniPar << " " << *iniPar << " " << *(iniPar+1) << std::endl;
207  func.SetParameters(iniPar);
208  iret |= DoFit<MinType>(fitobj,func,true );
209
210 #else
211  for (int i = 0; i < nfit; ++i) {
212  func.SetParameters(iniPar);
213  iret = DoFit<MinType>(fitobj,func, false);
214  if (iret != 0) {
215  std::cout << "Fit failed " << std::endl;
216  break;
217  }
218  }
219 #endif
220  w.Stop();
221  std::cout << "\nTime: \t" << w.RealTime() << " , " << w.CpuTime() << std::endl;
222  std::cout << "\n************************************************************\n";
223
224  return iret;
225 }
226
227
228 int testNdimFit() {
229
230
231  std::cout << "\n\n************************************************************\n";
232  std::cout << "\t UNBINNED TREE (GAUSSIAN MULTI-DIM) FIT\n";
233  std::cout << "************************************************************\n";
234
235  TTree * t1 = new TTree("t2","a large Tree with gaussian variables");
236  double x[N];
237  Int_t ev;
238  t1->Branch("x",x,branchType.c_str());
239  t1->Branch("ev",&ev,"ev/I");
240
241  // generate the true parameters
242  for (int j = 0; j < N; ++j) {
243  double mu = double(j)/10.;
244  double s = 1.0 + double(j)/10.;
245  truePar[2*j] = mu;
246  truePar[2*j+1] = s;
247  }
248
249
250  //fill the tree
251  TRandom3 r;
252  for (Int_t i=0;i<NPoints;i++) {
253  for (int j = 0; j < N; ++j) {
254  double mu = truePar[2*j];
255  double s = truePar[2*j+1];
256  x[j] = r.Gaus(mu,s);
257  }
258
259  ev = i;
260  t1->Fill();
261
262  }
263  //t1.Draw("x"); // to select fit variable
264
265
266  for (int i = 0; i <N; ++i) {
267  iniPar[2*i] = 0;
268  iniPar[2*i+1] = 1;
269  }
270
271  // use simply TF1 wrapper
272  //ROOT::Math::WrappedMultiTF1 f2(*f1);
274
275  int iret = 0;
276  iret |= FitUsingNewFitter<MINUIT2>(t1,f2);
277
278  return iret;
279 }
280
281 int main() {
282  return testNdimFit();
283 }
double truePar[2 *N]
Definition: testNdimFit.cxx:43
int DoUnBinFit(T *tree, Func &func, bool debug=false)
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
Definition: TBranch.cxx:2148
unsigned int NDim() const
return coordinate data dimension
Definition: UnBinData.h:280
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
void SetTolerance(double tol)
set the tolerance
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()
Definition: testNdimFit.cxx:65
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
static std::string name2()
Definition: testNdimFit.cxx:66
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
int main()
ROOT::Math::MinimizerOptions & MinimizerOptions()
Definition: FitConfig.h:138
TLatex * t1
Definition: textangle.C:20
ROOT::Fit::UnBinData * FillUnBinData(TTree *tree)
Definition: testNdimFit.cxx:70
const FitResult & Result() const
get fit result
Definition: Fitter.h:337
double sqrt(double)
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
Double_t x[n]
Definition: legend1.C:17
const std::string branchType
Definition: testNdimFit.cxx:34
const int N
Definition: testNdimFit.cxx:33
const FitConfig & Config() const
Definition: Fitter.h:365
unsigned int Size() const
return number of fit points
Definition: FitData.h:302
double iniPar[2 *N]
Definition: testNdimFit.cxx:44
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
const int NPoints
Definition: testNdimFit.cxx:35
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
int testNdimFit()
void SetMinimizer(const char *type, const char *algo=0)
set minimizer type
Definition: FitConfig.h:152
TRandom2 r(17)
Fitter class, entry point for performing all type of fits.
Definition: Fitter.h:77
double gausnorm(const double *x, const double *p)
Definition: testNdimFit.cxx:48
TMarker * m
Definition: textangle.C:8
int nfit
Definition: testFitPerf.cxx:59
double gausnormN(const double *x, const double *p)
Definition: testNdimFit.cxx:56
int FitUsingNewFitter(FitObj *fitobj, Func &func)
virtual Int_t GetEntry(Long64_t entry=0, Int_t getall=0)
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
int DoFit(TTree *tree, Func &func, bool debug=false)
virtual void SetParameters(const double *p)=0
Set the parameter values.
ROOT::Math::IParamMultiFunction Func
double f(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
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)
bool debug
Definition: tree.py:1
A TTree object has a header with a name and a title.
Definition: TTree.h:78
const int strategy
Definition: testNdimFit.cxx:46
A TTree is a list of TBranches.
Definition: TBranch.h:57
double exp(double)
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.
Stopwatch class.
Definition: TStopwatch.h:28