ROOT  6.06/09
Reference Guide
testDistSampler.cxx
Go to the documentation of this file.
1 // program to test distribution sampling
2 
3 #include "Math/Factory.h"
4 #include "Math/DistSampler.h"
5 #include "TF1.h"
6 #include "TH1.h"
7 
8 #include "TStopwatch.h"
9 #include "TRandom.h"
10 #include "TError.h"
11 #include "TCanvas.h"
12 // double Pdf(double x) {
13 // }
14 
15 using namespace ROOT::Math;
16 
17 
18 int testCont1D(int n) {
19 
20  // test gaussian
21 
22  DistSampler * sampler = Factory::CreateDistSampler("Unuran");
23  if (!sampler) return -1;
24  TF1 * f = new TF1("pdf","gaus");
25  f->SetParameters(1,0,1);
26 
27  TH1D * h1 = new TH1D("h1","h1",100,-3,3);
28  TH1D * hr = new TH1D("hr","h2",100,-3,3);
29 
30  sampler->SetFunction(*f,1);
31  bool ret = sampler->Init("AUTO");
32  if (!ret) return -1;
33 
34 
35  TStopwatch w;
36  w.Start();
37  for (int i = 0; i < n; ++i) {
38  h1->Fill(sampler->Sample1D() );
39  }
40  w.Stop();
41  double c = 1.E9/double(n);
42  std::cout << "Unuran sampling - (ns)/call = " << c*w.RealTime() << " " << c*w.CpuTime() << std::endl;
43  new TCanvas("Continous test");
44  h1->SetLineColor(kBlue);
45  h1->Draw();
46 
47  // generate ref histogram
48  w.Start();
49  for (int i = 0; i < n; ++i) {
50  hr->Fill(gRandom->Gaus(0,1) );
51  }
52  w.Stop();
53  std::cout << "TRandom::Gauss sampling - (ns)/call = " << c*w.RealTime() << " " << c*w.CpuTime() << std::endl;
54 
55  hr->Draw("SAME");
56 
57  // do a Chi2 test
58  // switch off printing of info messages from chi2 test
59  gErrorIgnoreLevel = 1001;
60  double prob = h1->Chi2Test(hr,"UU");
61  if (prob < 1.E-6) {
62  std::cerr << "Chi2 test of generated histogram failed" << std::endl;
63  return -2;
64  }
66 
67  return 0;
68 }
69 
70 
71 int testDisc1D(int n) {
72 
73  // test a discrete distribution
74 
75 
76  DistSampler * sampler = Factory::CreateDistSampler("Unuran");
77  if (!sampler) return -1;
78  TF1 * f = new TF1("pdf","TMath::Poisson(x,[0])");
79  double mu = 10;
80  f->SetParameter(0,mu);
81 
82  TH1D * h1 = new TH1D("h2","h1",50,0,50);
83  TH1D * hr = new TH1D("hd","h2",50,0,50);
84 
85  sampler->SetFunction(*f,1);
86  sampler->SetMode(mu);
87  sampler->SetArea(1);
88  bool ret = sampler->Init("DARI");
89  if (!ret) return -1;
90 
91 
92  TStopwatch w;
93  w.Start();
94  for (int i = 0; i < n; ++i) {
95  h1->Fill(sampler->Sample1D() );
96  }
97  w.Stop();
98  double c = 1.E9/double(n);
99  std::cout << "Unuran sampling - (ns)/call = " << c*w.RealTime() << " " << c*w.CpuTime() << std::endl;
100  new TCanvas("Discrete test");
101  h1->SetLineColor(kBlue);
102  h1->Draw();
103 
104  // generate ref histogram
105  w.Start();
106  for (int i = 0; i < n; ++i) {
107  hr->Fill(gRandom->Poisson(mu) );
108  }
109  w.Stop();
110  std::cout << "TRandom::Poisson sampling - (ns)/call = " << c*w.RealTime() << " " << c*w.CpuTime() << std::endl;
111 
112  hr->Draw("SAME");
113 
114  // do a Chi2 test
115  // switch off printing of info messages from chi2 test
116  gErrorIgnoreLevel = 1001;
117  double prob = h1->Chi2Test(hr,"UU");
118  if (prob < 1.E-6) {
119  std::cerr << "Chi2 test of generated histogram failed" << std::endl;
120  return -2;
121  }
122  gErrorIgnoreLevel = 0;
123 
124  return 0;
125 }
126 
127 
128 int testDistSampler(int n = 10000) {
129  int iret = 0;
130  iret |= testCont1D(n);
131  iret |= testDisc1D(n);
132  return iret;
133 }
134 int main() {
135  int iret = testDistSampler();
136  if (iret) std::cerr << "\ntestDistSampler: .... FAILED!" << std::endl;
137  else std::cerr << "\ntestDistSampler: .... OK" << std::endl;
138  return iret;
139 }
140 
virtual void SetMode(double)
set the mode of the distribution (could be useful to some methods) implemented by derived classes if ...
Definition: DistSampler.h:149
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3159
virtual void SetParameters(const Double_t *params)
Definition: TF1.h:439
Double_t RealTime()
Stop the stopwatch (if it is running) and return the realtime (in seconds) passed between the start a...
Definition: TStopwatch.cxx:108
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:56
R__EXTERN Int_t gErrorIgnoreLevel
Definition: TError.h:107
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
Double_t CpuTime()
Stop the stopwatch (if it is running) and return the cputime (in seconds) passed between the start an...
Definition: TStopwatch.cxx:123
int main()
THist< 1, double > TH1D
Definition: THist.h:314
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:75
void SetFunction(Function &func, unsigned int dim)
set the parent function distribution to use for sampling (generic case)
Definition: DistSampler.h:76
TH1F * h1
Definition: legend1.C:5
virtual Double_t Chi2Test(const TH1 *h2, Option_t *option="UU", Double_t *res=0) const
chi^{2} test for comparing weighted and unweighted histograms
Definition: TH1.cxx:1850
virtual void SetLineColor(Color_t lcolor)
Definition: TAttLine.h:54
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2878
int testDisc1D(int n)
Double_t E()
Definition: TMath.h:54
static ROOT::Math::DistSampler * CreateDistSampler(const std::string &samplerType="")
static method to create the distribution sampler class given a string specifying the type Supported s...
Definition: Factory.cxx:167
R__EXTERN TRandom * gRandom
Definition: TRandom.h:62
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:613
virtual double Sample1D()
sample one event in one dimension better implementation could be provided by the derived classes ...
Definition: DistSampler.h:165
virtual void SetArea(double)
set the normalization area of distribution implemented by derived classes if needed ...
Definition: DistSampler.h:153
The Canvas class.
Definition: TCanvas.h:48
Interface class for generic sampling of a distribution, i.e.
Definition: DistSampler.h:61
double f(double x)
int testDistSampler(int n=10000)
virtual bool Init(const char *="")
initialize the generators with the given algorithm Implemented by derived classes who needs it (like ...
Definition: DistSampler.h:105
1-Dim function class
Definition: TF1.h:149
int testCont1D(int n)
Definition: Rtypes.h:61
virtual void SetParameter(Int_t param, Double_t value)
Definition: TF1.h:431
virtual Int_t Poisson(Double_t mean)
Generates a random integer N according to a Poisson law.
Definition: TRandom.cxx:362
const Int_t n
Definition: legend1.C:16
Stopwatch class.
Definition: TStopwatch.h:30