Foam_Kanwa

This program can be execute from the command line as folows:

root -l foam_kanwa.C

Author: Stascek Jadach
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Saturday, June 19, 2021 at 10:42 AM.

In [1]:
%%cpp -d
#include "Riostream.h"
#include "TFoam.h"
#include "TCanvas.h"
#include "TH2.h"
#include "TMath.h"
#include "TFoamIntegrand.h"
#include "TRandom3.h"

In [2]:
%%cpp -d
Double_t sqr(Double_t x){
   return x*x;
}

In [3]:
%%cpp -d
Double_t Camel2(Int_t nDim, Double_t *Xarg){
   // 2-dimensional distribution for Foam, normalized to one (within 1e-5)
   Double_t x=Xarg[0];
   Double_t y=Xarg[1];
   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;
}
In [4]:
cout<<"--- kanwa started ---"<<endl;
TH2D  *hst_xy = new TH2D("hst_xy" ,  "x-y plot", 50,0,1.0, 50,0,1.0);
Double_t MCvect[2]; // 2-dim vector generated in the MC run
--- kanwa started ---
In [5]:
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->SetRhoInt(Camel2);  // Set 2-dim distribution, included below
FoamX->SetPseRan(PseRan);  // Set random number generator
FoamX->Initialize();       // Initialize simulator, may take time...
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
F                                                                              F
F                   ****************************************                   F
F                    ******      TFoam::Initialize    ******                   F
F                   ****************************************                   F
F                                                      FoamX                   F
F    Version =           1.02M   =                   Release date:  2005.04.10 F
F       kDim =          2 =  Dimension of the hyper-cubical space              F
F     nCells =        500 =  Requested number of Cells (half of them active)   F
F     nSampl =        200 =  No of MC events in exploration of a cell          F
F       nBin =          8 =  No of bins in histograms, MC exploration of cell  F
F   EvPerBin =         25 =  Maximum No effective_events/bin, MC exploration   F
F   OptDrive =          2 =  Type of Driver   =1,2 for Sigma,WtMax             F
F     OptRej =          1 =  MC rejection on/off for OptRej=0,1                F
F   MaxWtRej =             1.1   =       Maximum wt in rejection for wt=1 evts F
F                                                                              F
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
2222222222222222222222222222222222222222222222222
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
F                                                                              F
F                    ***  TFoam::Initialize FINISHED!!!  ***                   F
F     nCalls =      99800 =            Total number of function calls          F
F     XPrime =       1.3929609   =     Primary total integral                  F
F     XDiver =      0.39362177   =     Driver  total integral                  F
F   mcResult =      0.99933914   =     Estimate of the true MC Integral        F
F                                                                              F
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF

visualising generated distribution

In [6]:
TCanvas *cKanwa = new TCanvas("cKanwa","Canvas for plotting",600,600);
cKanwa->cd();

From now on foamx is ready to generate events

In [7]:
int nshow=5000;
for(long loop=0; loop<100000; loop++){
   FoamX->MakeEvent();            // generate MC event
   FoamX->GetMCvect( MCvect);     // get generated vector (x,y)
   Double_t x=MCvect[0];
   Double_t y=MCvect[1];
   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
(x,y) =  ( 0.68053985, 0.69250597 )
(x,y) =  ( 0.39743987, 0.22347379 )
(x,y) =  ( 0.41862392, 0.37268423 )
(x,y) =  ( 0.33221191, 0.37801703 )
(x,y) =  ( 0.32221499, 0.25437954 )
(x,y) =  ( 0.61444622, 0.60520452 )
(x,y) =  ( 0.30018061, 0.38244034 )
(x,y) =  ( 0.76521983, 0.777539 )
(x,y) =  ( 0.78407102, 0.69301713 )
(x,y) =  ( 0.72028183, 0.66087924 )
In [8]:
hst_xy->Draw("lego2");  // final plot
cKanwa->Update();
In [9]:
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;
 MCresult= 1.0002172 +- 0.00081511086
--- kanwa ended ---

Draw all canvases

In [10]:
%jsroot on
gROOT->GetListOfCanvases()->Draw()