NA49 logo

Input
GEANT
MTSIM
Embeding
Reconstruction
Maching
ROOT MINI-DST
Home page

Last modification:
2003-07-22
e-mail:
cblume@cern.ch korus@cern.ch

An example of a simple simulation

This is an example of a simple simulation which produces for each event one K* resonance. The output consists of two files: the first is a root file with some simple histograms and the second is a data file in ASCII format in which we have the number of the event, the number of produced particles in one event, and for each particle: px, py and pz. Here you can see the source code of this macro and here an example for the output.

#if !defined(__CINT__) || defined(__MAKECINT__)

#include "fstream.h"
#include "TH1.h"
#include "TH2.h"
#include "TF1.h"
#include "TFile.h"
#include "TRandom.h"
#include "TDatime.h"
#include "TSystem.h"
#include "TMath.h"

#endif

void generate_kstar(Int_t nOut = 1, Int_t nevent = 10000)
{

// Definitions
const Int_t kPid = 33; // PID
const Double_t kMass = 0.8961; // Mass
const Double_t kSlope = 0.28; // Slope of the mt-distribution (GeV)
const Double_t kWidth = 1.0; // Width of the y-distribution (Gaussian)
const Double_t kCenter = 2.9; // Center of the y-distribution (Gaussian)
const Double_t kYlow = 0.0; // Range of the y-distribution (Flat)
const Double_t kYhigh = 2.0 * kCenter;
const Int_t kNtrack = 1; // Number of tracks per event
const Double_t kPi = 3.141592654;

TH1F *hY = new TH1F("hY" ,"Rapidity", 12,0.0, 6.0);
TH1F *hPt = new TH1F("hPt" ,"p_{t}" , 12,0.0, 3.0);
TH1F *hPhi = new TH1F("hPhi","#phi" ,100,0.0,360.0);
TH1F *hPx = new TH1F("hPx" ,"p_{x}" ,100,0.0, 5.0);
TH1F *hPy = new TH1F("hPy" ,"p_{y}" ,100,0.0, 5.0);
TH1F *hPz = new TH1F("hPz" ,"p_{z}" ,100,0.0,100.0);
TH2F *hYPt = new TH2F("hYPt","p_{t} vs rapidity",12,0.0,6.0,12,0.0,3.0);

// Gaussian or flat rapidity distribution
Bool_t gaussian = kTRUE;

// The output filenames
Char_t cFileDat[150];
Char_t cFileRoot[150];
sprintf(cFileDat ,"gen_kstar_%02d.d" ,nOut);
sprintf(cFileRoot,"gen_kstar_%02d.root",nOut);

// EXECUTION PART

// Set a new random generator seed
TDatime *datime = new TDatime();
Int_t time = datime->GetTime();
Int_t date = datime->GetDate();
Int_t prid = gSystem->GetPid();
delete datime;
Int_t seed = TMath::Abs(10000 * prid + time - date);
gRandom->SetSeed(seed);

printf(" seed = %d\n",seed);
printf(" nevent = %d\n",nevent);

// The mt-distribution
TF1 *funMt = new TF1("funMt","x*exp(-sqrt([0]*[0]+x*x)/[1])",0.0,5.0);
funMt->SetParameter(0,kMass);
funMt->SetParameter(1,kSlope);

// Open the output file
FILE *fileData = fopen(cFileDat,"w+");

for (Int_t ievent = 0; ievent < nevent; ievent++) {

fprintf(fileData," %4d %6d\n",kNtrack,ievent);

for (Int_t itrack = 0; itrack < kNtrack; itrack++) {

// Rapidity
Double_t y;
if (gaussian) {
y = gRandom->Gaus(kCenter,kWidth);
}
else {
y = gRandom->Rndm() * (kYhigh - kYlow) + kYlow;
}
//Pt
Double_t pt = funMt->GetRandom();
// Phi
Double_t phi = gRandom->Rndm() * 2.0 * kPi;

Double_t px = pt * TMath::Cos(phi);
Double_t py = pt * TMath::Sin(phi);
Double_t pz = TMath::Sqrt(kMass*kMass + pt*pt) * TMath::SinH(y);

hY->Fill(y);
hPt->Fill(pt);
hYPt->Fill(y,pt);
hPhi->Fill(phi*180.0/kPi);
hPx->Fill(px);
hPy->Fill(py);
hPz->Fill(pz);

fprintf(fileData," %2d %8.3f %8.3f %8.3f\n" ,kPid,px,py,pz);

}
}

fclose(fileData);
TFile *fileHist = new TFile(cFileRoot,"RECREATE");
hY->Write();
hPt->Write();
hYPt->Write();
hPhi->Write();
hPx->Write();
hPy->Write();
hPz->Write();
fileHist->Close();

printf("generate done\n");
}