Logo ROOT   6.08/07
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 
93  cout << "NO ENCONTRADO!";
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");
268  gPad->Update();
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:439
virtual void Draw(Option_t *option="")
Draw this function with its current attributes.
Definition: TF2.cxx:216
Long64_t Fill(const Double_t *x, Double_t w=1.)
Definition: THnBase.h:153
Class to Wrap a ROOT Function class (like TF1) in a IParamMultiFunction interface of multi-dimensions...
Random number generator class based on the maximally quidistributed combined Tausworthe generator by ...
Definition: TRandom2.h:29
void DoFit(THnSparse *s, TF1 *f, ROOT::Fit::BinData &bd)
Definition: SparseFit4.cxx:162
unsigned int NPoints() const
return number of fit points
Definition: BinData.h:442
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:107
THnSparseT< TArrayD > THnSparseD
Definition: THnSparse.h:232
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:304
int findBin(ROOT::Fit::BinData &bd, const double *x)
Definition: SparseFit4.cxx:76
const FitResult & Result() const
get fit result
Definition: Fitter.h:354
static const double x2[5]
Double_t x[n]
Definition: legend1.C:17
Efficient multidimensional histogram.
Definition: THnSparse.h:52
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.
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
access to the fit configuration (const method)
Definition: Fitter.h:382
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:135
DataOptions : simple structure holding the options on how the data are filled.
Definition: DataOptions.h:28
TRandom2 r(17)
IParamFunction interface (abstract class) describing multi-dimensional parameteric functions It is a ...
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2851
A 3-Dim function with parameters.
Definition: TF3.h:30
Fitter class, entry point for performing all type of fits.
Definition: Fitter.h:94
unsigned int UInt_t
Definition: RtypesCore.h:42
ostream & operator<<(ostream &out, ROOT::Fit::BinData &bd)
Definition: SparseFit4.cxx:57
Double_t E()
Definition: TMath.h:54
Class describing the binned data sets : vectors of x coordinates, y values and optionally error on y ...
Definition: BinData.h:61
A 2-Dim function with parameters.
Definition: TF2.h:33
R__EXTERN TRandom * gRandom
Definition: TRandom.h:66
bool LikelihoodFit(const BinData &data, bool extended=true)
Binned Likelihood fit.
Definition: Fitter.h:181
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:1196
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:307
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:145
you should not use this method at all Int_t Int_t z
Definition: TRolke.cxx:630
unsigned int NDim() const
return coordinate data dimension
Definition: BinData.h:452
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:149
double gaus1D(double *x, double *p)
Definition: SparseFit4.cxx:35
bool showGraphics
Definition: SparseFit4.cxx:28
#define gPad
Definition: TVirtualPad.h:289
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:45
double gaus3D(double *x, double *p)
Definition: SparseFit4.cxx:45
Templated implementation of the abstract base THnSparse.
Definition: THnSparse.h:219
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:1225
void fitSparse3D()
Definition: SparseFit4.cxx:272
tomato 2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:296