ROOT   Reference Guide
Searching...
No Matches
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"
30
31#include <iostream>
32
33//_____________________________________________________________________________
34Double_t sqr(Double_t x){return x*x;};
35//_____________________________________________________________________________
36//_____________________________________________________________________________
37
38
39Double_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
50class 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
60TH2 * hFoam;
61TH2 * hUnr;
62
63
64Int_t run_foam(int nev){
65 cout<<"--- kanwa started ---"<<endl;
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
123double UCamel2(double * x, double *) {
124 return Camel2(2,x);
125}
126
127int run_unuran(int nev, std::string method = "hitro") {
128 // use unuran
129
130 std::cout << "run unuran " << std::endl;
131
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
173Int_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}
ROOT::R::TRInterface & r
Definition Object.C:4
#define f(i)
Definition RSha256.hxx:104
int Int_t
Definition RtypesCore.h:45
double Double_t
Definition RtypesCore.h:59
#define ClassDef(name, id)
Definition Rtypes.h:325
double exp(double)
R__EXTERN TSystem * gSystem
Definition TSystem.h:559
The Canvas class.
Definition TCanvas.h:23
Definition TCanvas.cxx:708
void Update() override
Definition TCanvas.cxx:2504
A 2-Dim function with parameters.
Definition TF2.h:29
Abstract class representing n-dimensional real positive integrand function.
virtual Double_t Density(Int_t ndim, Double_t *)=0
TFoam is the main class of the multi-dimensional general purpose Monte Carlo event generator (integra...
Definition TFoam.h:21
virtual void GetMCvect(Double_t *)
User may get generated MC point/vector with help of this method.
Definition TFoam.cxx:1176
virtual void MakeEvent()
User method.
Definition TFoam.cxx:1126
virtual void Initialize()
Basic initialization of FOAM invoked by the user.
Definition TFoam.cxx:323
virtual void GetIntegMC(Double_t &, Double_t &)
User method.
Definition TFoam.cxx:1211
virtual void SetnCells(Long_t nCells)
Definition TFoam.h:116
virtual void SetRho(TFoamIntegrand *Rho)
User may use this method to set the distribution object.
Definition TFoam.cxx:1022
virtual void SetPseRan(TRandom *PseRan)
Definition TFoam.h:112
virtual void SetkDim(Int_t kDim)
Definition TFoam.h:115
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:2005
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition TH1.cxx:3350
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition TH1.cxx:3073
2-D histogram with a double per channel (see TH1 documentation)}
Definition TH2.h:292
Service class for 2-Dim histogram classes.
Definition TH2.h:30
Int_t Fill(Double_t)
Invalid Fill method.
Definition TH2.cxx:294
void Divide(Int_t nx=1, Int_t ny=1, Float_t xmargin=0.01, Float_t ymargin=0.01, Int_t color=0) override
Random number generator class based on M.
Definition TRandom3.h:27
This is the base class for the ROOT Random number generators.
Definition TRandom.h:27
virtual void SetSeed(ULong_t seed=0)
Set the random generator seed.
Definition TRandom.cxx:608
Stopwatch class.
Definition TStopwatch.h:28
void Start(Bool_t reset=kTRUE)
Start the stopwatch.
Double_t CpuTime()
Stop the stopwatch (if it is running) and return the cputime (in seconds) passed between the start an...
void Stop()
Stop the stopwatch.
virtual int Load(const char *module, const char *entry="", Bool_t system=kFALSE)
Definition TSystem.cxx:1853
TUnuranMultiContDist class describing multi dimensional continuous distributions.
TUnuran class.
Definition TUnuran.h:79
Double_t y[n]
Definition legend1.C:17
Double_t x[n]
Definition legend1.C:17
const Int_t n
Definition legend1.C:16
TH1F * h1
Definition legend1.C:5
double Dist(void *xp, void *yp)
double dist(Rotation3D const &r1, Rotation3D const &r2)
VecExpr< UnaryOp< Sqr< T >, VecExpr< A, T, D >, T >, T, D > sqr(const VecExpr< A, T, D > &rhs)
constexpr Double_t Pi()
Definition TMath.h:37