 ROOT   Reference Guide
Searching...
No Matches
unuranFoamTest.C File Reference

## Detailed Description

This program must be compiled and executed with Aclic as follows.

.x unuranFoamTest.C+

it is an extension of tutorials foam_kanwa.C to compare generation of a 2D distribution with unuran and Foam

#include "TH2.h"
#include "TF2.h"
#include "TSystem.h"
#include "TCanvas.h"
#include "TMath.h"
#include "TRandom3.h"
#include "TFoam.h"
#include "TFoamIntegrand.h"
#include "TStopwatch.h"
#include "TROOT.h"
#include "TUnuran.h"
#include <iostream>
//_____________________________________________________________________________
Double_t sqr(Double_t x){return x*x;};
//_____________________________________________________________________________
//_____________________________________________________________________________
Double_t Camel2(Int_t nDim, Double_t *Xarg){
// 2-dimensional distribution for Foam, normalized to one (within 1e-5)
Double_t x=Xarg;
Double_t y=Xarg;
Double_t GamSq= sqr(0.100e0);
Double_t Dist= 0;
Dist +=exp(-(sqr(x-1./3) +sqr(y-1./3))/GamSq)/GamSq/TMath::Pi();
Dist +=exp(-(sqr(x-2./3) +sqr(y-2./3))/GamSq)/GamSq/TMath::Pi();
return 0.5*Dist;
}// Camel2
class FoamFunction : public TFoamIntegrand {
public:
virtual ~FoamFunction() {}
double Density(int nDim, double * x) {
return Camel2(nDim,x);
}
ClassDef(FoamFunction,1);
};
TH2 * hFoam;
TH2 * hUnr;
Int_t run_foam(int nev){
cout<<"--- kanwa started ---"<<endl;
TH2D *hst_xy = new TH2D("foam_hst_xy" , "FOAM x-y plot", 50,0,1.0, 50,0,1.0);
hFoam = hst_xy;
Double_t MCvect; // 2-dim vector generated in the MC run
//
TRandom *PseRan = new TRandom3(); // Create random number generator
PseRan->SetSeed(4357);
TFoam *FoamX = new TFoam("FoamX"); // Create Simulator
FoamX->SetkDim(2); // No. of dimensions, obligatory!
FoamX->SetnCells(500); // Optionally No. of cells, default=2000
FoamX->SetRho(new FoamFunction() ); // Set 2-dim distribution, included below
FoamX->SetPseRan(PseRan); // Set random number generator
//
// From now on FoamX is ready to generate events
// test first the time
w.Start();
FoamX->Initialize(); // Initialize simulator, may take time...
//int nshow=5000;
int nshow=nev;
for(long loop=0; loop<nev; loop++){
FoamX->MakeEvent(); // generate MC event
FoamX->GetMCvect( MCvect); // get generated vector (x,y)
Double_t x=MCvect;
Double_t y=MCvect;
//if(loop<10) cout<<"(x,y) = ( "<< x <<", "<< y <<" )"<<endl;
hst_xy->Fill(x,y);
// live plot
if(loop == nshow){
nshow += 5000;
hst_xy->Draw("lego2");
//cKanwa->Update();
}
}// loop
w.Stop();
double time = w.CpuTime()*1.E9/nev;
cout << "Time using FOAM \t\t " << " \t=\t " << time << "\tns/call" << endl;
//
hst_xy->Draw("lego2"); // final plot
//
Double_t MCresult, MCerror;
FoamX->GetIntegMC( MCresult, MCerror); // get MC integral, should be one
cout << " MCresult= " << MCresult << " +- " << MCerror <<endl;
cout<<"--- kanwa ended ---"<<endl;
return 0;
}//kanwa
double UCamel2(double * x, double *) {
return Camel2(2,x);
}
int run_unuran(int nev, std::string method = "hitro") {
// use unuran
std::cout << "run unuran " << std::endl;
TH2D *h1 = new TH2D("unr_hst_xy" , "UNURAN x-y plot", 50,0,1.0, 50,0,1.0);
hUnr= h1;
TF2 * f = new TF2("f",UCamel2,0,1,0,1,0);
TUnuran unr(&r,2); // 2 is debug level
// test first the time
w.Start();
// init unuran
bool ret = unr.Init(dist,method);
if (!ret) {
std::cerr << "Error initializing unuran with method " << unr.MethodName() << endl;
return -1;
}
double x;
for (int i = 0; i < nev; ++i) {
unr.SampleMulti(x);
h1->Fill(x,x);
// if (method == "gibbs" && i < 100)
// std::cout << x << " , " << x << std::endl;
}
w.Stop();
double time = w.CpuTime()*1.E9/nev;
cout << "Time using Unuran " << unr.MethodName() << " \t=\t " << time << "\tns/call" << endl;
h1->Draw("lego2");
return 0;
}
Int_t unuranFoamTest(){
// visualising generated distribution
TCanvas *cKanwa = new TCanvas("cKanwa","Canvas for plotting",600,1000);
cKanwa->Divide(1,2);
cKanwa->cd(1);
int n = 100000;
run_unuran(n,"hitro");
cKanwa->Update();
cKanwa->cd(2);
run_foam(n);
cKanwa->Update();
std::cout <<"\nChi2 Test Results (UNURAN-FOAM):\t";
// test chi2
hFoam->Chi2Test(hUnr,"UUP");
return 0;
}
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(Rotation3D const &r1, Rotation3D const &r2)
constexpr Double_t Pi()
Definition TMath.h:37

Definition in file unuranFoamTest.C.