ROOT   Reference Guide
Searching...
No Matches
foam_demo.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_FOAM
3/// \notebook -nodraw
4/// Demonstrate the TFoam class.
5///
6/// To run this macro type from CINT command line
7///
8/// ~~~{.cpp}
9/// root [0] gSystem->Load("libFoam.so")
10/// root [1] .x foam_demo.C+
11/// ~~~
12///
13/// \macro_code
14///
15/// \author Stascek Jadach
16
17
18#include "Riostream.h"
19#include "TFile.h"
20#include "TFoam.h"
21#include "TH1.h"
22#include "TMath.h"
23#include "TFoamIntegrand.h"
24#include "TRandom3.h"
25
26class TFDISTR: public TFoamIntegrand {
27public:
28 TFDISTR(){};
29 Double_t Density(int nDim, Double_t *Xarg){
30 // Integrand for mFOAM
31 Double_t Fun1,Fun2,R1,R2;
32 Double_t pos1=1e0/3e0;
33 Double_t pos2=2e0/3e0;
34 Double_t Gam1= 0.100e0; // as in JPC
35 Double_t Gam2= 0.100e0; // as in JPC
36 Double_t sPi = sqrt(TMath::Pi());
37 Double_t xn1=1e0;
38 Double_t xn2=1e0;
39 int i;
40 R1=0;
41 R2=0;
42 for(i = 0 ; i<nDim ; i++){
43 R1=R1+(Xarg[i] -pos1)*(Xarg[i] -pos1);
44 R2=R2+(Xarg[i] -pos2)*(Xarg[i] -pos2);
45 xn1=xn1*Gam1*sPi;
46 xn2=xn2*Gam2*sPi;
47 }
48 R1 = sqrt(R1);
49 R2 = sqrt(R2);
50 Fun1 = exp(-(R1*R1)/(Gam1*Gam1))/xn1; // Gaussian delta-like profile
51 Fun2 = exp(-(R2*R2)/(Gam2*Gam2))/xn2; // Gaussian delta-like profile
52 return 0.5e0*(Fun1+ Fun2);
53}
54 ClassDef(TFDISTR,1) //Class of testing functions for FOAM
55};
56ClassImp(TFDISTR)
57
58Int_t foam_demo()
59{
60 TFile RootFile("foam_demo.root","RECREATE","histograms");
61 long loop;
62 Double_t MCresult,MCerror,MCwt;
63 //-----------------------------------------
64 long NevTot = 50000; // Total MC statistics
65 Int_t kDim = 2; // total dimension
66 Int_t nCells = 500; // Number of Cells
67 Int_t nSampl = 200; // Number of MC events per cell in build-up
68 Int_t nBin = 8; // Number of bins in build-up
69 Int_t OptRej = 1; // Wted events for OptRej=0; wt=1 for OptRej=1 (default)
70 Int_t OptDrive = 2; // (D=2) Option, type of Drive =0,1,2 for TrueVol,Sigma,WtMax
71 Int_t EvPerBin = 25; // Maximum events (equiv.) per bin in buid-up
72 Int_t Chat = 1; // Chat level
73 //-----------------------------------------
74 TRandom *PseRan = new TRandom3(); // Create random number generator
75 TFoam *FoamX = new TFoam("FoamX"); // Create Simulator
76 TFoamIntegrand *rho= new TFDISTR();
77 PseRan->SetSeed(4357);
78 //-----------------------------------------
79 cout<<"***** Demonstration Program for Foam version "<<FoamX->GetVersion()<<" *****"<<endl;
80 FoamX->SetkDim( kDim); // Mandatory!!!
81 FoamX->SetnCells( nCells); // optional
82 FoamX->SetnSampl( nSampl); // optional
83 FoamX->SetnBin( nBin); // optional
84 FoamX->SetOptRej( OptRej); // optional
85 FoamX->SetOptDrive( OptDrive); // optional
86 FoamX->SetEvPerBin( EvPerBin); // optional
87 FoamX->SetChat( Chat); // optional
88 //-----------------------------------------
89 FoamX->SetRho(rho);
90 FoamX->SetPseRan(PseRan);
91 FoamX->Initialize(); // Initialize simulator
92 FoamX->Write("FoamX"); // Writing Foam on the disk, TESTING PERSISTENCY!!!
93 //-----------------------------------------
94 long nCalls=FoamX->GetnCalls();
95 cout << "====== Initialization done, entering MC loop" << endl;
96 //-----------------------------------------
97 /*cout<<" About to start MC loop: "; cin.getline(question,20);*/
98 Double_t *MCvect =new Double_t[kDim]; // vector generated in the MC run
99 //-----------------------------------------
100 TH1D *hst_Wt = new TH1D("hst_Wt" , "Main weight of Foam",25,0,1.25);
101 hst_Wt->Sumw2();
102 //-----------------------------------------
103 for(loop=0; loop<NevTot; loop++){
104 /*===============================*/
105 FoamX->MakeEvent(); // generate MC event
106 /*===============================*/
107 FoamX->GetMCvect( MCvect);
108 MCwt=FoamX->GetMCwt();
109 hst_Wt->Fill(MCwt,1.0);
110 if(loop<15){
111 cout<<"MCwt= "<<MCwt<<", ";
112 cout<<"MCvect= ";
113 for ( Int_t k=0 ; k<kDim ; k++) cout<<MCvect[k]<<" "; cout<< endl;
114 }
115 if( ((loop)%100000)==0 ){
116 cout<<" loop= "<<loop<<endl;
117 }
118 }
119
120 //-----------------------------------------
121
122 cout << "====== Events generated, entering Finalize" << endl;
123
124 hst_Wt->Print("all");
125 Double_t eps = 0.0005;
126 Double_t Effic, WtMax, AveWt, Sigma;
127 Double_t IntNorm, Errel;
128 FoamX->Finalize( IntNorm, Errel); // final printout
129 FoamX->GetIntegMC( MCresult, MCerror); // get MC intnegral
130 FoamX->GetWtParams(eps, AveWt, WtMax, Sigma); // get MC wt parameters
131 Effic=0; if(WtMax>0) Effic=AveWt/WtMax;
132 cout << "================================================================" << endl;
133 cout << " MCresult= " << MCresult << " +- " << MCerror << " RelErr= "<< MCerror/MCresult << endl;
134 cout << " Dispersion/<wt>= " << Sigma/AveWt << endl;
135 cout << " <wt>/WtMax= " << Effic <<", for epsilon = "<<eps << endl;
136 cout << " nCalls (initialization only) = " << nCalls << endl;
137 cout << "================================================================" << endl;
138
139 delete [] MCvect;
140 //
141 RootFile.ls();
142 RootFile.Write();
143 RootFile.Close();
144 cout << "***** End of Demonstration Program *****" << endl;
145
146 return 0;
147}
148
149
int Int_t
Definition RtypesCore.h:45
double Double_t
Definition RtypesCore.h:59
#define ClassDef(name, id)
Definition Rtypes.h:325
#define ClassImp(name)
Definition Rtypes.h:364
double sqrt(double)
double exp(double)
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition TFile.h:54
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 SetEvPerBin(Int_t EvPerBin)
Definition TFoam.h:122
virtual void Initialize()
Basic initialization of FOAM invoked by the user.
Definition TFoam.cxx:323
virtual void SetnSampl(Long_t nSampl)
Definition TFoam.h:117
virtual void GetIntegMC(Double_t &, Double_t &)
User method.
Definition TFoam.cxx:1211
virtual const char * GetVersion() const
Definition TFoam.h:127
virtual void GetMCwt(Double_t &)
User may get weight MC weight using this method.
Definition TFoam.cxx:1191
virtual void SetChat(Int_t Chat)
Definition TFoam.h:119
virtual void SetOptDrive(Int_t OptDrive)
Definition TFoam.h:121
virtual void SetnCells(Long_t nCells)
Definition TFoam.h:116
virtual Long_t GetnCalls() const
Definition TFoam.h:131
virtual void SetRho(TFoamIntegrand *Rho)
User may use this method to set the distribution object.
Definition TFoam.cxx:1022
virtual void Finalize(Double_t &, Double_t &)
May be called optionally by the user after the MC run.
Definition TFoam.cxx:1260
virtual void SetOptRej(Int_t OptRej)
Definition TFoam.h:120
virtual void SetPseRan(TRandom *PseRan)
Definition TFoam.h:112
virtual void SetnBin(Int_t nBin)
Definition TFoam.h:118
virtual void SetkDim(Int_t kDim)
Definition TFoam.h:115
virtual void GetWtParams(Double_t, Double_t &, Double_t &, Double_t &)
May be called optionally after the MC run.
Definition TFoam.cxx:1247
1-D histogram with a double per channel (see TH1 documentation)}
Definition TH1.h:618
virtual void Print(Option_t *option="") const
Print some global quantities for this histogram.
Definition TH1.cxx:6975
virtual Int_t Fill(Double_t x)
Increment bin with abscissa X by 1.
Definition TH1.cxx:3350
virtual void Sumw2(Bool_t flag=kTRUE)
Create structure to store sum of squares of weights.
Definition TH1.cxx:8860
virtual Int_t Write(const char *name=0, Int_t option=0, Int_t bufsize=0)
Write this object to the current directory.
Definition TObject.cxx:798
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
constexpr Double_t Sigma()
Stefan-Boltzmann constant in .
Definition TMath.h:277
constexpr Double_t Pi()
Definition TMath.h:37
#define R1(v, w, x, y, z, i)
Definition sha1.inl:134
#define R2(v, w, x, y, z, i)
Definition sha1.inl:137