Notebook

foam_demo¶

Demonstrate the TFoam class.

To run this macro type from ROOT command line

root [0] gSystem->Load("libFoam.so")
 root [1] .x foam_demo.C+

Author: Stascek Jadach
This notebook tutorial was automatically generated with ROOTBOOK-izer from the macro found in the ROOT repository on Friday, March 24, 2023 at 10:47 AM.

In [1]:

%%cpp -d


#include "Riostream.h"
#include "TFile.h"
#include "TFoam.h"
#include "TH1.h"
#include "TMath.h"
#include "TFoamIntegrand.h"
#include "TRandom3.h"

class TFDISTR: public TFoamIntegrand {
public:
   TFDISTR(){};
   Double_t Density(int nDim, Double_t *Xarg){
   // Integrand for mFOAM
   Double_t Fun1,Fun2,R1,R2;
   Double_t pos1=1e0/3e0;
   Double_t pos2=2e0/3e0;
   Double_t Gam1= 0.100e0;  // as in JPC
   Double_t Gam2= 0.100e0;  // as in JPC
   Double_t sPi = sqrt(TMath::Pi());
   Double_t xn1=1e0;
   Double_t xn2=1e0;
   int i;
   R1=0;
   R2=0;
   for(i = 0 ; i<nDim ; i++){
      R1=R1+(Xarg[i] -pos1)*(Xarg[i] -pos1);
      R2=R2+(Xarg[i] -pos2)*(Xarg[i] -pos2);
      xn1=xn1*Gam1*sPi;
      xn2=xn2*Gam2*sPi;
   }
   R1   = sqrt(R1);
   R2   = sqrt(R2);
   Fun1 = exp(-(R1*R1)/(Gam1*Gam1))/xn1;  // Gaussian delta-like profile
   Fun2 = exp(-(R2*R2)/(Gam2*Gam2))/xn2;  // Gaussian delta-like profile
   return 0.5e0*(Fun1+ Fun2);
}
  ClassDef(TFDISTR,1) //Class of testing functions for FOAM
};

Arguments are defined.

In [2]:

TFile RootFile("foam_demo.root","RECREATE","histograms");
long   loop;
Double_t MCresult,MCerror,MCwt;

In [3]:

long NevTot   =     50000;   // Total MC statistics
Int_t  kDim   =         2;   // total dimension
Int_t  nCells   =     500;   // Number of Cells
Int_t  nSampl   =     200;   // Number of MC events per cell in build-up
Int_t  nBin     =       8;   // Number of bins in build-up
Int_t  OptRej   =       1;   // Wted events for OptRej=0; wt=1 for OptRej=1 (default)
Int_t  OptDrive =       2;   // (D=2) Option, type of Drive =0,1,2 for TrueVol,Sigma,WtMax
Int_t  EvPerBin =      25;   // Maximum events (equiv.) per bin in buid-up
Int_t  Chat     =       1;   // Chat level

In [4]:

TRandom *PseRan   = new TRandom3();  // Create random number generator
TFoam   *FoamX    = new TFoam("FoamX");   // Create Simulator
TFoamIntegrand    *rho= new TFDISTR();
PseRan->SetSeed(4357);

In [5]:

cout<<"*****   Demonstration Program for Foam version "<<FoamX->GetVersion()<<"    *****"<<endl;
FoamX->SetkDim(        kDim);      // Mandatory!!!
FoamX->SetnCells(      nCells);    // optional
FoamX->SetnSampl(      nSampl);    // optional
FoamX->SetnBin(        nBin);      // optional
FoamX->SetOptRej(      OptRej);    // optional
FoamX->SetOptDrive(    OptDrive);  // optional
FoamX->SetEvPerBin(    EvPerBin);  // optional
FoamX->SetChat(        Chat);      // optional

*****   Demonstration Program for Foam version 1.02M    *****

In [6]:

FoamX->SetRho(rho);
FoamX->SetPseRan(PseRan);
FoamX->Initialize(); // Initialize simulator
FoamX->Write("FoamX");     // Writing Foam on the disk, TESTING PERSISTENCY!!!

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

In [7]:

long nCalls=FoamX->GetnCalls();
cout << "====== Initialization done, entering MC loop" << endl;

====== Initialization done, entering MC loop

In [8]:

/*cout<<" About to start MC loop: ";  cin.getline(question,20);*/
Double_t *MCvect =new Double_t[kDim]; // vector generated in the MC run

In [9]:

TH1D  *hst_Wt = new TH1D("hst_Wt" ,  "Main weight of Foam",25,0,1.25);
hst_Wt->Sumw2();

In [10]:

for(loop=0; loop<NevTot; loop++){
/*===============================*/
   FoamX->MakeEvent();           // generate MC event
/*===============================*/
   FoamX->GetMCvect( MCvect);
   MCwt=FoamX->GetMCwt();
   hst_Wt->Fill(MCwt,1.0);
   if(loop<15){
      cout<<"MCwt= "<<MCwt<<",  ";
      cout<<"MCvect= ";
      for ( Int_t k=0 ; k<kDim ; k++) cout<<MCvect[k]<<" "; cout<< endl;
   }
   if( ((loop)%100000)==0 ){
      cout<<"   loop= "<<loop<<endl;
   }
}

MCwt= 1,  MCvect= 0.68053985 0.69250597 
   loop= 0
MCwt= 1,  MCvect= 0.39743987 0.22347379 
MCwt= 1,  MCvect= 0.41862392 0.37268423 
MCwt= 1,  MCvect= 0.33221191 0.37801703 
MCwt= 1,  MCvect= 0.32221499 0.25437954 
MCwt= 1,  MCvect= 0.61444622 0.60520452 
MCwt= 1,  MCvect= 0.30018061 0.38244034 
MCwt= 1,  MCvect= 0.76521983 0.777539 
MCwt= 1,  MCvect= 0.78407102 0.69301713 
MCwt= 1,  MCvect= 0.72028183 0.66087924 
MCwt= 1,  MCvect= 0.73255425 0.64108329 
MCwt= 1,  MCvect= 0.24614236 0.23288176 
MCwt= 1,  MCvect= 0.34828797 0.51862775 
MCwt= 1,  MCvect= 0.44052705 0.39626735 
MCwt= 1,  MCvect= 0.65947295 0.59166751

In [11]:

cout << "====== Events generated, entering Finalize" << endl;

hst_Wt->Print("all");
Double_t eps = 0.0005;
Double_t Effic, WtMax, AveWt, Sigma;
Double_t IntNorm, Errel;
FoamX->Finalize(   IntNorm, Errel);     // final printout
FoamX->GetIntegMC( MCresult, MCerror);  // get MC intnegral
FoamX->GetWtParams(eps, AveWt, WtMax, Sigma); // get MC wt parameters
Effic=0; if(WtMax>0) Effic=AveWt/WtMax;
cout << "================================================================" << endl;
cout << " MCresult= " << MCresult << " +- " << MCerror << " RelErr= "<< MCerror/MCresult << endl;
cout << " Dispersion/<wt>= " << Sigma/AveWt << endl;
cout << "      <wt>/WtMax= " << Effic <<",    for epsilon = "<<eps << endl;
cout << " nCalls (initialization only) =   " << nCalls << endl;
cout << "================================================================" << endl;

delete [] MCvect;

====== Events generated, entering Finalize
TH1.Print Name  = hst_Wt, Entries= 50000, Total sum= 49993
 fSumw[0]=0, x=-0.025, error=0
 fSumw[1]=0, x=0.025, error=0
 fSumw[2]=0, x=0.075, error=0
 fSumw[3]=0, x=0.125, error=0
 fSumw[4]=0, x=0.175, error=0
 fSumw[5]=0, x=0.225, error=0
 fSumw[6]=0, x=0.275, error=0
 fSumw[7]=0, x=0.325, error=0
 fSumw[8]=0, x=0.375, error=0
 fSumw[9]=0, x=0.425, error=0
 fSumw[10]=0, x=0.475, error=0
 fSumw[11]=0, x=0.525, error=0
 fSumw[12]=0, x=0.575, error=0
 fSumw[13]=0, x=0.625, error=0
 fSumw[14]=0, x=0.675, error=0
 fSumw[15]=0, x=0.725, error=0
 fSumw[16]=0, x=0.775, error=0
 fSumw[17]=0, x=0.825, error=0
 fSumw[18]=0, x=0.875, error=0
 fSumw[19]=0, x=0.925, error=0
 fSumw[20]=0, x=0.975, error=0
 fSumw[21]=49991, x=1.025, error=223.587
 fSumw[22]=1, x=1.075, error=1
 fSumw[23]=1, x=1.125, error=1
 fSumw[24]=0, x=1.175, error=0
 fSumw[25]=0, x=1.225, error=0
 fSumw[26]=7, x=1.275, error=2.64575
TH1.Print Name  = TFoamMaxwt_hst_Wt1, Entries= 76619, Total sum= 76619
TH1.Print Name  = TFoamMaxwt_hst_Wt2, Entries= 76619, Total sum= 54981.3
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
F                                                                              F
F                   ****************************************                   F
F                    ******     TFoam::Finalize       ******                   F
F                   ****************************************                   F
F     NevGen =      76619 = Number of generated events in the MC generation    F
F     nCalls =     176419 = Total number of function calls                     F
F                   ----------------------------------------                   F
F      AveWt =      0.71759331   =     Average MC weight                       F
F      WtMin =    1.472245e-19   =     Minimum MC weight (absolute)            F
F      WtMax =       4.1676265   =     Maximum MC weight (absolute)            F
F                   ----------------------------------------                   F
F     XPrime =       1.3929609   =     Primary total integral, R_prime         F
F     XDiver =      0.39362177   =     Driver  total integral, R_loss          F
F                   ----------------------------------------                   F
F      IntMC =      0.99957943 +-      0.00115528  = Result of the MC Integral F
F   mCerelat =    0.0011557661   =     Relative error of the MC integral       F
F  <w>/WtMax =       0.7211993   =              MC efficiency, acceptance rate F
F  Sigma/<w> =      0.31991763   =              MC efficiency, variance/ave_wt F
F      WtMax =           0.995   =              WtMax(esp= 0.0005)             F
F      Sigma =      0.22957075   =              variance of MC weight          F
F <OveW>/<W> =   0.00018763091   =              Contrib. of events wt>MaxWtRej F
F                                                                              F
FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFF
TH1.Print Name  = TFoamMaxwt_hst_Wt1, Entries= 76619, Total sum= 76619
TH1.Print Name  = TFoamMaxwt_hst_Wt2, Entries= 76619, Total sum= 54981.3
================================================================
 MCresult= 0.99957943 +- 0.00115528 RelErr= 0.0011557661
 Dispersion/<wt>= 0.31991763
      <wt>/WtMax= 0.7211993,    for epsilon = 0.0005
 nCalls (initialization only) =   99800
================================================================

In [12]:

RootFile.ls();
RootFile.Write();
RootFile.Close();
cout << "***** End of Demonstration Program  *****" << endl;

return 0;

TFile**		foam_demo.root	histograms
 TFile*		foam_demo.root	histograms
  OBJ: TH1D	hst_Wt	Main weight of Foam : 0 at: 0x7f5f82f92f50
  KEY: TFoam	FoamX;1	General purpose self-adapting Monte Carlo event generator
  KEY: TProcessID	ProcessID0;1	45513a12-ca31-11ed-9471-942c8a89beef
***** End of Demonstration Program  *****