A Simple C++ Simulation For Beginners

Phenomenology in physics, mostly deals with the simulation of events and obtaining data from simulations to compare it with real time event datas. Obviously during event processing, it’s not necessary to give extra effort to visuality. So one should not confuse it with visual simulations. We’r just making event based calculations. Therefore you can ask what exactly do we simulate? or can any calculation be a simulation? Notice that in scientific experiments, you always need a satisfactory amount of statistics. So basically you should have a scenario for gathering statistics in simulations. Here i’d like to present a calculation of “Pi” number as a simulation sample. Here we are collecting statistics via producing random numbers which is included a circle with r=1.

Calculation: Pi number
The method: Monte Carlo Simulation
Fundamental Formulas: (pi)r 2 and x2+y2 = 1 (Note that radius of circle is unit 1.)

1-) Write below code and compile it writing “g++ pi.cpp -o pi.x”

#include<iostream>
#include<math.h>
#include<stdlib.h>
#include<time.h>

using namespace std;
int main(){
	int jmax=1000; // maximum value of HIT number. (Length of output file)
	int imax=1000; // maximum value of random numbers for producing HITs.
	double x,y;    // Coordinates
	int hit;       // storage variable of number of HITs
	srand(time(0));
	for (int j=0;j<jmax;j++){
		hit=0;
		x=0; y=0;
		for(int i=0;i<imax;i++){
			x=double(rand())/double(RAND_MAX);
			y=double(rand())/double(RAND_MAX);
		if(y<=sqrt(1-pow(x,2))) hit+=1; }          //Choosing HITs according to analytic formula of circle
	cout<<""<<4*double(hit)/double(imax)<<endl; }  // Print out Pi number
}

2-) To understand the code: We have just 2 loops here. The inner loop produce random number (<1) and uses these numbers for coordinates x,y. “If” condition increases hit number if this (x,y) point locates in the area of quarter circle (Look at the figure below.)

Hit Production Area

The outer loop resets our variables and print out Pi number according to formula: Area of Quarter Circle/ Area of Square = (1/4)πr2/r2 = (1/4)π= accepted hits / total hits = hits / imax.

3-) Run it as “./pi.x > pi.dat”

4-) Draw output file.

Pi Graph

I used below “root macro” to read and convert it to a .root file.

{
	gROOT->Reset();
	ifstream in;
	in.open("pi.dat");
	Float_t x; Int_t nlines = 0;
	TFile *f = new TFile("pi1.root","RECREATE");
	TH1F *h1 = new TH1F("h1","pi_grafik",100,2.5,4.0);
	TNtuple *ntuple = new TNtuple("ntuple","pi","x");
	for (nlines=0; nlines<10000; nlines++) {
		in >> x;
		if (!in.good()) {break;}
		if (nlines < 5) {printf("x=%5f\n",x);}
		h1->Fill(x);
		ntuple->Fill(x);
		nlines++;
	}
	in.close();
	f->Write();
	printf("%d deger bulundu\n",nlines);
	h1->SetXTitle("pi");
	h1->SetYTitle("Olay");
	h1->Draw();
}
  • Paste above root macro in a C file and name it as “pintuple.c”
  • Open your root analysis program in the same directory you saved pintuple.c : “root”
  • Execute the file: “root> .x pintuple.c”
  • You’ll get an ntuple file called “pi1.root”
  • Write “TBrowser g” in root.
  • Open pi1.root file and you’ll get the above histogram. Congratulations 🙂
Advertisements

Tags: ,

About Kuday

Pure Mind

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s

%d bloggers like this: