Logo ROOT   master
Reference Guide
unuranFoamTest.C
Go to the documentation of this file.
1 /// \file
2 /// \ingroup tutorial_unuran
3 /// This program must be compiled and executed with Aclic as follows
4 ///
5 /// ~~~{.cpp}
6 /// .x unuranFoamTest.C+
7 /// ~~~
8 ///
9 /// it is an extension of tutorials foam_kanwa.C to compare
10 /// generation of a 2D distribution with unuran and Foam
11 ///
12 /// \macro_code
13 ///
14 /// \author Lorenzo Moneta
15 
16 #include "TH2.h"
17 #include "TF2.h"
18 #include "TSystem.h"
19 #include "TCanvas.h"
20 #include "TMath.h"
21 #include "TRandom3.h"
22 #include "TFoam.h"
23 #include "TFoamIntegrand.h"
24 #include "TStopwatch.h"
25 #include "TROOT.h"
26 
27 
28 #include "TUnuran.h"
29 #include "TUnuranMultiContDist.h"
30 
31 #include <iostream>
32 
33 //_____________________________________________________________________________
34 Double_t sqr(Double_t x){return x*x;};
35 //_____________________________________________________________________________
36 //_____________________________________________________________________________
37 
38 
39 Double_t Camel2(Int_t nDim, Double_t *Xarg){
40 // 2-dimensional distribution for Foam, normalized to one (within 1e-5)
41  Double_t x=Xarg[0];
42  Double_t y=Xarg[1];
43  Double_t GamSq= sqr(0.100e0);
44  Double_t Dist= 0;
45  Dist +=exp(-(sqr(x-1./3) +sqr(y-1./3))/GamSq)/GamSq/TMath::Pi();
46  Dist +=exp(-(sqr(x-2./3) +sqr(y-2./3))/GamSq)/GamSq/TMath::Pi();
47  return 0.5*Dist;
48 }// Camel2
49 
50 class FoamFunction : public TFoamIntegrand {
51  public:
52  virtual ~FoamFunction() {}
53  double Density(int nDim, double * x) {
54  return Camel2(nDim,x);
55  }
56  ClassDef(FoamFunction,1);
57 
58 };
59 
60 TH2 * hFoam;
61 TH2 * hUnr;
62 
63 
64 Int_t run_foam(int nev){
65  cout<<"--- kanwa started ---"<<endl;
66  gSystem->Load("libFoam.so");
67  TH2D *hst_xy = new TH2D("foam_hst_xy" , "FOAM x-y plot", 50,0,1.0, 50,0,1.0);
68  hFoam = hst_xy;
69 
70  Double_t MCvect[2]; // 2-dim vector generated in the MC run
71  //
72  TRandom *PseRan = new TRandom3(); // Create random number generator
73  PseRan->SetSeed(4357);
74  TFoam *FoamX = new TFoam("FoamX"); // Create Simulator
75  FoamX->SetkDim(2); // No. of dimensions, obligatory!
76  FoamX->SetnCells(500); // Optionally No. of cells, default=2000
77  FoamX->SetRho(new FoamFunction() ); // Set 2-dim distribution, included below
78  FoamX->SetPseRan(PseRan); // Set random number generator
79  //
80  // From now on FoamX is ready to generate events
81 
82  // test first the time
83  TStopwatch w;
84 
85  w.Start();
86  FoamX->Initialize(); // Initialize simulator, may take time...
87 
88  //int nshow=5000;
89  int nshow=nev;
90 
91  for(long loop=0; loop<nev; loop++){
92  FoamX->MakeEvent(); // generate MC event
93  FoamX->GetMCvect( MCvect); // get generated vector (x,y)
94  Double_t x=MCvect[0];
95  Double_t y=MCvect[1];
96  //if(loop<10) cout<<"(x,y) = ( "<< x <<", "<< y <<" )"<<endl;
97  hst_xy->Fill(x,y);
98  // live plot
99  if(loop == nshow){
100  nshow += 5000;
101  hst_xy->Draw("lego2");
102  //cKanwa->Update();
103  }
104  }// loop
105  w.Stop();
106 
107  double time = w.CpuTime()*1.E9/nev;
108  cout << "Time using FOAM \t\t " << " \t=\t " << time << "\tns/call" << endl;
109 
110  //
111  hst_xy->Draw("lego2"); // final plot
112  //
113  Double_t MCresult, MCerror;
114  FoamX->GetIntegMC( MCresult, MCerror); // get MC integral, should be one
115  cout << " MCresult= " << MCresult << " +- " << MCerror <<endl;
116  cout<<"--- kanwa ended ---"<<endl;
117 
118  return 0;
119 }//kanwa
120 
121 
122 
123 double UCamel2(double * x, double *) {
124  return Camel2(2,x);
125 }
126 
127 int run_unuran(int nev, std::string method = "hitro") {
128  // use unuran
129 
130  std::cout << "run unuran " << std::endl;
131 
132  gSystem->Load("libUnuran.so");
133 
134  TH2D *h1 = new TH2D("unr_hst_xy" , "UNURAN x-y plot", 50,0,1.0, 50,0,1.0);
135  hUnr= h1;
136 
137  TF2 * f = new TF2("f",UCamel2,0,1,0,1,0);
138 
140 
141  TRandom3 r;
142 
143  TUnuran unr(&r,2); // 2 is debug level
144 
145 
146  // test first the time
147  TStopwatch w;
148 
149  w.Start();
150 
151  // init unuran
152  bool ret = unr.Init(dist,method);
153  if (!ret) {
154  std::cerr << "Error initializing unuran with method " << unr.MethodName() << endl;
155  return -1;
156  }
157 
158  double x[2];
159  for (int i = 0; i < nev; ++i) {
160  unr.SampleMulti(x);
161  h1->Fill(x[0],x[1]);
162 // if (method == "gibbs" && i < 100)
163 // std::cout << x[0] << " , " << x[1] << std::endl;
164  }
165 
166  w.Stop();
167  double time = w.CpuTime()*1.E9/nev;
168  cout << "Time using Unuran " << unr.MethodName() << " \t=\t " << time << "\tns/call" << endl;
169  h1->Draw("lego2");
170  return 0;
171 }
172 
173 Int_t unuranFoamTest(){
174 
175  // visualising generated distribution
176  TCanvas *cKanwa = new TCanvas("cKanwa","Canvas for plotting",600,1000);
177  cKanwa->Divide(1,2);
178  cKanwa->cd(1);
179  int n = 100000;
180 
181 
182  run_unuran(n,"hitro");
183  cKanwa->Update();
184 
185  cKanwa->cd(2);
186 
187  run_foam(n);
188  cKanwa->Update();
189 
190 
191  std::cout <<"\nChi2 Test Results (UNURAN-FOAM):\t";
192  // test chi2
193  hFoam->Chi2Test(hUnr,"UUP");
194 
195  return 0;
196 }
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition: TH1.cxx:3275
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
Random number generator class based on M.
Definition: TRandom3.h:27
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Definition: TStopwatch.cxx:58
TVirtualPad * cd(Int_t subpadnumber=0)
Set current canvas & pad.
Definition: TCanvas.cxx:701
virtual int Load(const char *module, const char *entry="", Bool_t system=kFALSE)
Load a shared library.
Definition: TSystem.cxx:1849
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
#define f(i)
Definition: RSha256.hxx:104
int Int_t
Definition: RtypesCore.h:43
virtual void MakeEvent()
User method.
Definition: TFoam.cxx:1123
Abstract class representing n-dimensional real positive integrand function.
Definition: TFoamIntegrand.h:9
void Stop()
Stop the stopwatch.
Definition: TStopwatch.cxx:77
Double_t x[n]
Definition: legend1.C:17
#define ClassDef(name, id)
Definition: Rtypes.h:322
This is the base class for the ROOT Random number generators.
Definition: TRandom.h:27
virtual Double_t Chi2Test(const TH1 *h2, Option_t *option="UU", Double_t *res=0) const
test for comparing weighted and unweighted histograms
Definition: TH1.cxx:1950
virtual void SetSeed(ULong_t seed=0)
Set the random generator seed.
Definition: TRandom.cxx:597
TH1F * h1
Definition: legend1.C:5
constexpr Double_t Pi()
Definition: TMath.h:38
double Dist(void *xp, void *yp)
ROOT::R::TRInterface & r
Definition: Object.C:4
Service class for 2-Dim histogram classes.
Definition: TH2.h:30
virtual void GetMCvect(Double_t *)
User may get generated MC point/vector with help of this method.
Definition: TFoam.cxx:1173
R__EXTERN TSystem * gSystem
Definition: TSystem.h:556
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:2998
virtual Double_t Density(Int_t ndim, Double_t *)=0
virtual void Initialize()
Basic initialization of FOAM invoked by the user.
Definition: TFoam.cxx:320
TUnuranMultiContDist class describing multi dimensional continuous distributions. ...
A 2-Dim function with parameters.
Definition: TF2.h:29
The Canvas class.
Definition: TCanvas.h:27
virtual void SetnCells(Long_t nCells)
Definition: TFoam.h:116
TUnuran class.
Definition: TUnuran.h:79
double Double_t
Definition: RtypesCore.h:57
Double_t y[n]
Definition: legend1.C:17
virtual void SetPseRan(TRandom *PseRan)
Definition: TFoam.h:112
virtual void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0)
Automatic pad generation by division.
Definition: TPad.cxx:1165
virtual void Update()
Update canvas pad buffers.
Definition: TCanvas.cxx:2433
double exp(double)
Int_t Fill(Double_t)
Invalid Fill method.
Definition: TH2.cxx:294
VecExpr< UnaryOp< Sqr< T >, VecExpr< A, T, D >, T >, T, D > sqr(const VecExpr< A, T, D > &rhs)
const Int_t n
Definition: legend1.C:16
virtual void SetRho(TFoamIntegrand *Rho)
User may use this method to set the distribution object.
Definition: TFoam.cxx:1019
virtual void GetIntegMC(Double_t &, Double_t &)
User method.
Definition: TFoam.cxx:1208
TFoam is the main class of the multi-dimensional general purpose Monte Carlo event generator (integra...
Definition: TFoam.h:21
virtual void SetkDim(Int_t kDim)
Definition: TFoam.h:115
2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:292
Stopwatch class.
Definition: TStopwatch.h:28