ROOT   6.10/09 Reference Guide
SparseFit4.cxx
Go to the documentation of this file.
1 #include "THnSparse.h"
2 #include "TH2.h"
3 #include "TF2.h"
4 #include "TF3.h"
5 #include "TCanvas.h"
6 #include "TApplication.h"
7 #include "TMath.h"
8
9 #include "Fit/SparseData.h"
10 #include "HFitInterface.h"
11 #include "Fit/Fitter.h"
12 #include "Math/WrappedMultiTF1.h"
13
14 #include "TRandom2.h"
15
16 #include <iostream>
17 #include <iterator>
18 #include <algorithm>
19
20 #include <list>
21 #include <vector>
22
23 #include <cmath>
24 #include <cassert>
25
26 using namespace std;
27
28 bool showGraphics = false;
29
30 TRandom2 r(17);
31 Int_t bsize[] = { 10, 10, 10 };
32 Double_t x_min[] = { 0., 0., 0. };
33 Double_t x_max[] = { 10., 10., 10. };
34
35 double gaus1D(double *x, double *p)
36 {
37  return p[0]*TMath::Gaus(x[0],p[1],p[2]);
38 }
39
40 double gaus2D(double *x, double *p)
41 {
42  return p[0]*TMath::Gaus(x[0],p[1],p[2]) * TMath::Gaus(x[1],p[3],p[4]);
43 }
44
45 double gaus3D(double *x, double *p)
46 {
47  return p[0] * TMath::Gaus(x[0],p[1],p[2])
48  * TMath::Gaus(x[1],p[3],p[4])
49  * TMath::Gaus(x[2],p[5],p[6]);
50 }
51
52 double pol2D(double *x, double *p)
53 {
54  return p[0]*x[0]+ p[1]*x[0]*x[0] + p[2]*x[1] + p[3]*x[1]*x[1] + p[4];
55 }
56
57 ostream& operator << (ostream& out, ROOT::Fit::BinData& bd)
58 {
59  const unsigned int ndim( bd.NDim() );
60  const unsigned int npoints( bd.NPoints() );
61  for ( unsigned int i = 0; i < npoints; ++i )
62  {
63  double value = 0, error = 0;
64  const double *x = bd.GetPoint(i, value, error);
65  for ( unsigned int j = 0; j < ndim; ++j )
66  {
67  out << " x[" << j << "]: " << x[j];
68  }
69  out << " value: " << value;
70  out << " error: " << error;
71  out << endl;
72  }
73  return out;
74 }
75
76 int findBin(ROOT::Fit::BinData& bd, const double *x)
77 {
78  const unsigned int ndim = bd.NDim();
79  const unsigned int npoints = bd.NPoints();
80
81  for ( unsigned int i = 0; i < npoints; ++i )
82  {
83  double value1 = 0, error1 = 0;
84  const double *x1 = bd.GetPoint(i, value1, error1);
85  bool thisIsIt = true;
86  for ( unsigned int j = 0; j < ndim; ++j )
87  {
88  thisIsIt &= fabs(x1[j] - x[j]) < 1E-15;
89  }
90  if ( thisIsIt ) return i;
91  }
92
94  copy(x, x+ndim, ostream_iterator<double>(cout, " " ));
95  cout << endl;
96
97  return -1;
98 }
99
101 {
102  const unsigned int ndim = bd1.NDim();
103  const unsigned int npoints = bd1.NPoints();
104
105  bool equals = true;
106
107  cout << "Equals" << endl;
108
109  for ( unsigned int i = 0; i < npoints && equals; ++i )
110  {
111  double value1 = 0, error1 = 0;
112  const double *x1 = bd1.GetPoint(i, value1, error1);
113
114  int bin = findBin(bd2, x1);
115
116  double value2 = 0, error2 = 0;
117  const double *x2 = bd2.GetPoint(bin, value2, error2);
118
119  equals &= ( value1 == value2 );
120  cout << " v: " << equals;
121  equals &= ( error1 == error2 );
122  cout << " e: " << equals;
123  for ( unsigned int j = 0; j < ndim; ++j )
124  {
125  equals &= fabs(x1[j] - x2[j]) < 1E-15;
126  cout << " x[" << j << "]: " << equals;
127  }
128
129  cout << " bd1: ";
130  std::copy(x1, &x1[ndim], ostream_iterator<double>(cout, " "));
131  cout << " value:" << value1;
132  cout << " error:" << error1;
133
134  cout << " bd2: ";
135  std::copy(x2, &x2[ndim], ostream_iterator<double>(cout, " "));
136  cout << " value:" << value2;
137  cout << " error:" << error2;
138
139  cout << " equals: " << equals;
140
141  cout << endl;
142  }
143
144  return equals;
145 }
146
147 void fillSparse(THnSparse* s, TF1* f, int nEvents = 5)
148 {
149  const unsigned int ndim = s->GetNdimensions();
150
151  for ( Int_t e = 0; e < nEvents; ++e ) {
152  Double_t *points = new Double_t[ndim];
153  for ( UInt_t i = 0; i < ndim; ++ i )
154  points[i] = r.Uniform( x_min[0] * .9 , x_max[0] * 1.1 );
155  double value = gRandom->Poisson( f->EvalPar(points));
156  s->Fill(points, value );
157  cout << value << " " << s->GetNbins() << endl;
158  delete [] points;
159  }
160 }
161
163 {
164  ///////////////// CREATE THE SPARSE DATA
165  cout << "DoFit: dim = " << s->GetNdimensions() << " - Retrieving the Sparse Data Structure" << endl;
166  //ROOT::Fit::SparseData d(s);
168  ROOT::Fit::FillData(d, s, 0);
169  d.GetBinData(bd);
170
171  ///////////////// FITS
172  // Fit preparation
173  bool ret;
174  ROOT::Fit::Fitter fitter;
177  fitter.Config().SetMinimizer("Minuit2");
178
180  opt.fUseEmpty = true;
181  opt.fIntegral = true;
182
183  /////////////////////////////////////////////////////////////////////////
184  cout << "\n ******* Likelihood with BinData and NoZeros *******" << endl;
185  ROOT::Fit::BinData bdNoZeros;
186  d.GetBinDataNoZeros(bdNoZeros);
187  ret = fitter.LikelihoodFit(bdNoZeros, if2, true);
188  fitter.Result().Print(std::cout);
189  if (!ret)
190  std::cout << "Fit Failed " << std::endl;
191
192  /////////////////////////////////////////////////////////////////////////
193  cout << "\n ******* Likelihood with BinData with Zeros *******" << endl;
194  ROOT::Fit::BinData bdWithZeros(opt);
195  d.GetBinDataIntegral(bdWithZeros);
196  ret = fitter.LikelihoodFit(bdWithZeros, if2, true);
197  fitter.Result().Print(std::cout);
198  if (!ret)
199  std::cout << "Fit Failed " << std::endl;
200
201  /////////////////////////////////////////////////////////////////////////
202 }
203
205 {
206
207  std::cout << "1D SPARSE FIT \n" << std::endl;
208
209  const unsigned int ndim = 1;
210
211  THnSparseD* s1 = new THnSparseD("mfND-s1", "s1-Title", ndim, bsize, x_min, x_max);
212
213  TF1* f = new TF1("func1D", gaus1D, x_min[0] - 2, x_max[0] + 2, 3);
214  f->SetParameters(10., 5., 2.);
215
216  fillSparse(s1,f,5);
217
218  cout << "1D Fit : Retrieving the Sparse Data Structure" << endl;
219  //ROOT::Fit::SparseData d(s1);
221  ROOT::Fit::FillData(d, s1, 0);
222 }
223
225 {
226
227  std::cout << "2D SPARSE FIT \n" << std::endl;
228  const unsigned int ndim = 2;
229
230  THnSparseD* s1 = new THnSparseD("mfND-s1", "s1-Title", ndim, bsize, x_min, x_max);
231
232  TF2* f = new TF2("func2D", gaus2D, x_min[0],x_max[0], x_min[1], x_max[1], 5);
233  f->SetParameters(40,5,2,5,1);
234
235  for (int ix=1; ix <= bsize[0]; ++ix) {
236  for (int iy=1; iy <= bsize[1]; ++iy) {
237  double x = s1->GetAxis(0)->GetBinCenter(ix);
238  double y = s1->GetAxis(1)->GetBinCenter(iy);
239
240  double value = gRandom->Poisson( f->Eval(x,y) );
241  if ( value )
242  {
243  double points[] = {x,y};
244  s1->Fill( points, value );
245  }
246  }
247  }
248
250  DoFit(s1, f, bd);
251
252  TH2D* h2 = new TH2D("2D Blanked Hist Fit", "h1-title",
253  bsize[0], x_min[0], x_max[0],
254  bsize[1], x_min[1], x_max[1]);
255  //cout << "Filling second histogram" << endl;
256  for ( unsigned int i = 0; i < bd.NPoints(); ++i )
257  {
258  const double* x;
259  double value, error;
260  x = bd.GetPoint(i, value, error);
261  value = (value)?value:-10;
262  h2->Fill(x[0], x[1], value);
263  }
264
265  if (showGraphics) {
266  h2->Draw("colZ");
267  f->Draw("SAME");
269  }
270 }
271
273 {
274  const unsigned int ndim = 3;
275
276  std::cout << "3D SPARSE FIT \n" << std::endl;
277
278  THnSparseD* s1 = new THnSparseD("mfND-s1", "s1-Title", ndim, bsize, x_min, x_max);
279
280  TF3* f = new TF3("func3D", gaus3D,
281  x_min[0],x_max[0],
282  x_min[1],x_max[1],
283  x_min[2],x_max[2],
284  7);
285  f->SetParameters(100,5,2,5,1,5,2);
286
287  for (int ix=1; ix <= bsize[0]; ++ix) {
288  for (int iy=1; iy <= bsize[1]; ++iy) {
289  for (int iz=1; iz <= bsize[2]; ++iz) {
290  double x = s1->GetAxis(0)->GetBinCenter(ix);
291  double y = s1->GetAxis(1)->GetBinCenter(iy);
292  double z = s1->GetAxis(2)->GetBinCenter(iz);
293  double value = gRandom->Poisson( f->Eval(x,y,z) );
294
295  if ( value )
296  {
297  double points[] = {x,y,z};
298  s1->Fill( points, value );
299  }
300  }
301  }
302  }
303
305  DoFit(s1, f, bd);
306 }
307
308
309 int main(int argc, char** argv)
310 {
311
312  bool do3dfit = false;
313
314  // Parse command line arguments
315  for (Int_t i=1 ; i<argc ; i++) {
316  std::string arg = argv[i] ;
317  if (arg == "-g") {
318  showGraphics = true;
319  }
320  if (arg == "-v") {
321  showGraphics = true;
322  //verbose = true;
323  }
324  if (arg == "-3d") {
325  do3dfit = true;
326  }
327  if (arg == "-h") {
328  cerr << "Usage: " << argv[0] << " [-g] [-v]\n";
329  cerr << " where:\n";
330  cerr << " -g : graphics mode\n";
331  cerr << " -v : verbose mode\n";
332  cerr << " -3d : 3d fit";
333  cerr << endl;
334  return -1;
335  }
336  }
337
338
339  TApplication* theApp = 0;
340
341  if ( showGraphics )
342  theApp = new TApplication("App",&argc,argv);
343
344  // fitSparse1D();
345  fitSparse2D();
346  if (do3dfit) fitSparse3D();
347
348  cout << "C'est fini!" << endl;
349
350  if ( showGraphics ) {
351  theApp->Run();
352  delete theApp;
353  theApp = 0;
354  }
355
356  return 0;
357 }
358
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:588
virtual void Draw(Option_t *option="")
Draw this function with its current attributes.
Definition: TF2.cxx:217
Long64_t Fill(const Double_t *x, Double_t w=1.)
Definition: THnBase.h:143
Random number generator class based on the maximally quidistributed combined Tausworthe generator by ...
Definition: TRandom2.h:27
void DoFit(THnSparse *s, TF1 *f, ROOT::Fit::BinData &bd)
Definition: SparseFit4.cxx:162
int Int_t
Definition: RtypesCore.h:41
int main(int argc, char **argv)
Definition: SparseFit4.cxx:309
int equals(Double_t n1, Double_t n2, double ERRORLIMIT=1.E-10)
Definition: UnitTesting.cxx:24
STL namespace.
Long64_t GetNbins() const
Definition: THnSparse.h:91
THnSparseT< TArrayD > THnSparseD
Definition: THnSparse.h:216
double pol2D(double *x, double *p)
Definition: SparseFit4.cxx:52
const double * GetPoint(unsigned int ipoint, double &value) const
retrieve at the same time a pointer to the coordinate data and the fit value More efficient than call...
Definition: BinData.h:347
int findBin(ROOT::Fit::BinData &bd, const double *x)
Definition: SparseFit4.cxx:76
const FitResult & Result() const
get fit result
Definition: Fitter.h:337
static const double x2[5]
Double_t x[n]
Definition: legend1.C:17
Efficient multidimensional histogram.
Definition: THnSparse.h:36
Int_t bsize[]
Definition: SparseFit4.cxx:31
bool operator==(ROOT::Fit::BinData &bd1, ROOT::Fit::BinData &bd2)
Definition: SparseFit4.cxx:100
virtual void Run(Bool_t retrn=kFALSE)
Main application eventloop. Calls system dependent eventloop via gSystem.
Class to Wrap a ROOT Function class (like TF1) in a IParamMultiFunction interface of multi-dimensions...
void fitSparse1D()
Definition: SparseFit4.cxx:204
Double_t x_min[]
Definition: SparseFit4.cxx:32
virtual Double_t GetBinCenter(Int_t bin) const
Return center of bin.
Definition: TAxis.cxx:464
void fillSparse(THnSparse *s, TF1 *f, int nEvents=5)
Definition: SparseFit4.cxx:147
const FitConfig & Config() const
Definition: Fitter.h:365
point * points
Definition: X3DBuffer.c:20
const int nEvents
Definition: testRooFit.cxx:42
void fitSparse2D()
Definition: SparseFit4.cxx:224
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
void SetMinimizer(const char *type, const char *algo=0)
set minimizer type
Definition: FitConfig.h:152
void FillData(BinData &dv, const TH1 *hist, TF1 *func=0)
fill the data vector from a TH1.
TAxis * GetAxis(Int_t dim) const
Definition: THnBase.h:125
DataOptions : simple structure holding the options on how the data are filled.
Definition: DataOptions.h:28
TRandom2 r(17)
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2851
A 3-Dim function with parameters.
Definition: TF3.h:28
Fitter class, entry point for performing all type of fits.
Definition: Fitter.h:77
unsigned int UInt_t
Definition: RtypesCore.h:42
ostream & operator<<(ostream &out, ROOT::Fit::BinData &bd)
Definition: SparseFit4.cxx:57
Class describing the binned data sets : vectors of x coordinates, y values and optionally error on y ...
Definition: BinData.h:53
A 2-Dim function with parameters.
Definition: TF2.h:29
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
constexpr Double_t E()
Definition: TMath.h:74
bool LikelihoodFit(const BinData &data, bool extended=true)
Binned Likelihood fit.
Definition: Fitter.h:164
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 Double_t Eval(Double_t x, Double_t y=0, Double_t z=0, Double_t t=0) const
Evaluate this function.
Definition: TF1.cxx:1197
static const double x1[5]
double f(double x)
double Double_t
Definition: RtypesCore.h:55
Double_t y[n]
Definition: legend1.C:17
you should not use this method at all Int_t Int_t Double_t Double_t Double_t e
Definition: TRolke.cxx:630
THist< 2, double, THistStatContent, THistStatUncertainty > TH2D
Definition: THist.hxx:316
virtual Double_t Uniform(Double_t x1=1)
Returns a uniform deviate on the interval (0, x1).
Definition: TRandom.cxx:606
Int_t GetNdimensions() const
Definition: THnBase.h:135
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
void Print(std::ostream &os, bool covmat=false) const
print the result and optionaly covariance matrix and correlations
Definition: FitResult.cxx:446
1-Dim function class
Definition: TF1.h:150
double gaus1D(double *x, double *p)
Definition: SparseFit4.cxx:35
bool showGraphics
Definition: SparseFit4.cxx:28
double gaus2D(double *x, double *p)
Definition: SparseFit4.cxx:40
This class creates the ROOT Application Environment that interfaces to the windowing system eventloop...
Definition: TApplication.h:39
unsigned int NPoints() const
return number of fit points
Definition: FitData.h:294
double gaus3D(double *x, double *p)
Definition: SparseFit4.cxx:45
Templated implementation of the abstract base THnSparse.
Definition: THnSparse.h:203
unsigned int NDim() const
return coordinate data dimension
Definition: FitData.h:310
virtual Int_t Poisson(Double_t mean)
Generates a random integer N according to a Poisson law.
Definition: TRandom.cxx:362
Double_t x_max[]
Definition: SparseFit4.cxx:33
Int_t Fill(Double_t)
Invalid Fill method.
Definition: TH2.cxx:292
virtual Double_t EvalPar(const Double_t *x, const Double_t *params=0)
Evaluate function with given coordinates and parameters.
Definition: TF1.cxx:1226
void fitSparse3D()
Definition: SparseFit4.cxx:272
tomato 2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:290