Logo ROOT  
Reference Guide
pythiaExample.C
Go to the documentation of this file.
1/// \file
2/// \ingroup tutorial_pythia
3/// Using Pythia6 with ROOT
4///
5/// To make an event sample (of size 100) do
6///
7/// ~~~{.cpp}
8/// shell> root
9/// root [0] .L pythiaExample.C
10/// root [1] makeEventSample(1000)
11/// ~~~
12///
13/// To start the tree view on the generated tree, do
14///
15/// ~~~{.cpp}
16/// shell> root
17/// root [0] .L pythiaExample.C
18/// root [1] showEventSample()
19/// ~~~
20///
21/// The following session:
22///
23/// ~~~{.cpp}
24/// shell> root
25/// root [0] .x pythiaExample.C(500)
26/// ~~~
27///
28/// will execute makeEventSample(500) and showEventSample()
29///
30/// Alternatively, you can compile this to a program
31/// and then generate 1000 events with
32///
33/// ~~~{.cpp}
34/// ./pythiaExample 1000
35/// ~~~
36///
37/// To use the program to start the viewer, do
38///
39/// ~~~{.cpp}
40/// ./pythiaExample -1
41/// ~~~
42///
43/// NOTE 1: To run this example, you must have a version of ROOT
44/// compiled with the Pythia6 version enabled and have Pythia6 installed.
45/// The statement `gSystem->Load("$HOME/pythia6/libPythia6");` (see below)
46/// assumes that the directory containing the Pythia6 library
47/// is in the pythia6 subdirectory of your $HOME. Locations
48/// that can specify this, are:
49///
50/// ~~~{.cpp}
51/// Root.DynamicPath resource in your ROOT configuration file
52/// (/etc/root/system.rootrc or ~/.rootrc).
53/// Runtime load paths set on the executable (Using GNU ld,
54/// specified with flag `-rpath').
55/// Dynamic loader search path as specified in the loaders
56/// configuration file (On GNU/Linux this file is
57/// etc/ld.so.conf).
58/// For Un*x: Any directory mentioned in LD_LIBRARY_PATH
59/// For Windows: Any directory mentioned in PATH
60/// ~~~
61///
62/// NOTE 2: The example can also be run with ACLIC:
63///
64/// ~~~{.cpp}
65/// root > gSystem->Load("libEG");
66/// root > gSystem->Load("$ROOTSYS/../pythia6/libPythia6"); //change to your setup
67/// root > gSystem->Load("libEGPythia6");
68/// root > .x pythiaExample.C+
69/// ~~~
70///
71/// \macro_code
72///
73/// \author Christian Holm Christensen
74
75#include "TApplication.h"
76#include "TPythia6.h"
77#include "TFile.h"
78#include "TError.h"
79#include "TTree.h"
80#include "TClonesArray.h"
81#include "TH1.h"
82#include "TF1.h"
83#include "TStyle.h"
84#include "TLatex.h"
85#include "TCanvas.h"
86#include "Riostream.h"
87#include <cstdlib>
88using namespace std;
89
90
91#define FILENAME "pythia.root"
92#define TREENAME "tree"
93#define BRANCHNAME "particles"
94#define HISTNAME "ptSpectra"
95#define PDGNUMBER 211
96
97// This function just load the needed libraries if we're executing from
98// an interactive session.
99void loadLibraries()
100{
101#ifdef R__MACOSX
102 gSystem->Load("libpythia6_dummy");
103#endif
104}
105
106// nEvents is how many events we want.
107int makeEventSample(Int_t nEvents)
108{
109 // Load needed libraries
110 loadLibraries();
111
112 // Create an instance of the Pythia event generator ...
113 TPythia6* pythia = new TPythia6;
114
115 // ... and initialise it to run p+p at sqrt(200) GeV in CMS
116 pythia->Initialize("cms", "p", "p", 200);
117
118 // Open an output file
119 TFile* file = TFile::Open(FILENAME, "RECREATE");
120 if (!file || !file->IsOpen()) {
121 Error("makeEventSample", "Couldn;t open file %s", FILENAME);
122 return 1;
123 }
124
125 // Make a tree in that file ...
126 TTree* tree = new TTree(TREENAME, "Pythia 6 tree");
127
128 // ... and register a the cache of pythia on a branch (It's a
129 // TClonesArray of TMCParticle objects. )
130 TClonesArray* particles = (TClonesArray*)pythia->GetListOfParticles();
131 tree->Branch(BRANCHNAME, &particles);
132
133 // Now we make some events
134 for (Int_t i = 0; i < nEvents; i++) {
135 // Show how far we got every 100'th event.
136 if (i % 100 == 0)
137 cout << "Event # " << i << endl;
138
139 // Make one event.
140 pythia->GenerateEvent();
141
142 // Maybe you want to have another branch with global event
143 // information. In that case, you should process that here.
144 // You can also filter out particles here if you want.
145
146 // Now we're ready to fill the tree, and the event is over.
147 tree->Fill();
148 }
149
150 // Show tree structure
151 tree->Print();
152
153 // After the run is over, we may want to do some summary plots:
154 TH1D* hist = new TH1D(HISTNAME, "p_{#perp} spectrum for #pi^{+}",
155 100, 0, 3);
156 hist->SetXTitle("p_{#perp}");
157 hist->SetYTitle("dN/dp_{#perp}");
158 char expression[64];
159 sprintf(expression,"sqrt(pow(%s.fPx,2)+pow(%s.fPy,2))>>%s",
160 BRANCHNAME, BRANCHNAME, HISTNAME);
161 char selection[64];
162 sprintf(selection,"%s.fKF==%d", BRANCHNAME, PDGNUMBER);
163 tree->Draw(expression,selection);
164
165 // Normalise to the number of events, and the bin sizes.
166 hist->Sumw2();
167 hist->Scale(3 / 100. / hist->Integral());
168 hist->Fit("expo", "QO+", "", .25, 1.75);
169 TF1* func = hist->GetFunction("expo");
170 func->SetParNames("A", "-1/T");
171 // and now we flush and close the file
172 file->Write();
173 file->Close();
174
175 return 0;
176}
177
178// Show the Pt spectra, and start the tree viewer.
179int showEventSample()
180{
181 // Load needed libraries
182 loadLibraries();
183
184 // Open the file
185 TFile* file = TFile::Open(FILENAME, "READ");
186 if (!file || !file->IsOpen()) {
187 Error("showEventSample", "Couldn;t open file %s", FILENAME);
188 return 1;
189 }
190
191 // Get the tree
192 TTree* tree = (TTree*)file->Get(TREENAME);
193 if (!tree) {
194 Error("showEventSample", "couldn't get TTree %s", TREENAME);
195 return 2;
196 }
197
198 // Start the viewer.
199 tree->StartViewer();
200
201 // Get the histogram
202 TH1D* hist = (TH1D*)file->Get(HISTNAME);
203 if (!hist) {
204 Error("showEventSample", "couldn't get TH1D %s", HISTNAME);
205 return 4;
206 }
207
208 // Draw the histogram in a canvas
209 gStyle->SetOptStat(1);
210 TCanvas* canvas = new TCanvas("canvas", "canvas");
211 canvas->SetLogy();
212 hist->Draw("e1");
213 TF1* func = hist->GetFunction("expo");
214
215 char expression[64];
216 sprintf(expression,"T #approx %5.1f", -1000 / func->GetParameter(1));
217 TLatex* latex = new TLatex(1.5, 1e-4, expression);
218 latex->SetTextSize(.1);
219 latex->SetTextColor(4);
220 latex->Draw();
221
222 return 0;
223}
224
225void pythiaExample(Int_t n=1000) {
226 makeEventSample(n);
227 showEventSample();
228}
229
230#ifndef __CINT__
231int main(int argc, char** argv)
232{
233 TApplication app("app", &argc, argv);
234
235 Int_t n = 100;
236 if (argc > 1)
237 n = strtol(argv[1], NULL, 0);
238
239 int retVal = 0;
240 if (n > 0)
241 retVal = makeEventSample(n);
242 else {
243 retVal = showEventSample();
244 app.Run();
245 }
246
247 return retVal;
248}
249#endif
250
#define e(i)
Definition: RSha256.hxx:103
int Int_t
Definition: RtypesCore.h:45
void Error(const char *location, const char *msgfmt,...)
Use this function in case an error occurred.
Definition: TError.cxx:187
R__EXTERN TStyle * gStyle
Definition: TStyle.h:412
R__EXTERN TSystem * gSystem
Definition: TSystem.h:559
int main(int argc, char *argv[])
Definition: cef_main.cxx:54
This class creates the ROOT Application Environment that interfaces to the windowing system eventloop...
Definition: TApplication.h:39
virtual void SetTextColor(Color_t tcolor=1)
Set the text color.
Definition: TAttText.h:43
virtual void SetTextSize(Float_t tsize=1)
Set the text size.
Definition: TAttText.h:46
The Canvas class.
Definition: TCanvas.h:23
An array of clone (identical) objects.
Definition: TClonesArray.h:29
1-Dim function class
Definition: TF1.h:213
virtual void SetParNames(const char *name0="p0", const char *name1="p1", const char *name2="p2", const char *name3="p3", const char *name4="p4", const char *name5="p5", const char *name6="p6", const char *name7="p7", const char *name8="p8", const char *name9="p9", const char *name10="p10")
Set up to 10 parameter names.
Definition: TF1.cxx:3482
virtual Double_t GetParameter(Int_t ipar) const
Definition: TF1.h:512
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition: TFile.h:54
static TFile * Open(const char *name, Option_t *option="", const char *ftitle="", Int_t compress=ROOT::RCompressionSetting::EDefaults::kUseCompiledDefault, Int_t netopt=0)
Create / open a file.
Definition: TFile.cxx:4011
virtual TObjArray * GetListOfParticles() const
Definition: TGenerator.h:176
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:618
virtual void SetXTitle(const char *title)
Definition: TH1.h:413
virtual TFitResultPtr Fit(const char *formula, Option_t *option="", Option_t *goption="", Double_t xmin=0, Double_t xmax=0)
Fit histogram with function fname.
Definition: TH1.cxx:3893
virtual Double_t Integral(Option_t *option="") const
Return integral of bin contents.
Definition: TH1.cxx:7829
virtual TF1 * GetFunction(const char *name) const
Return pointer to function with name.
Definition: TH1.cxx:8894
virtual void Draw(Option_t *option="")
Draw this histogram with options.
Definition: TH1.cxx:3074
virtual void SetYTitle(const char *title)
Definition: TH1.h:414
virtual void Scale(Double_t c1=1, Option_t *option="")
Multiply this histogram by a constant c1.
Definition: TH1.cxx:6566
virtual void Sumw2(Bool_t flag=kTRUE)
Create structure to store sum of squares of weights.
Definition: TH1.cxx:8863
To draw Mathematical Formula.
Definition: TLatex.h:18
virtual void Draw(Option_t *option="")
Default Draw method for all objects.
Definition: TObject.cxx:197
void SetLogy(Int_t value=1) override
Set Lin/Log scale for Y.
Definition: TPad.cxx:5947
TPythia is an interface class to F77 version of Pythia 6.2
Definition: TPythia6.h:84
void Initialize(const char *frame, const char *beam, const char *target, float win)
Calls PyInit with the same parameters after performing some checking, sets correct title.
Definition: TPythia6.cxx:437
void GenerateEvent()
generate event and copy the information from /HEPEVT/ to fPrimaries
Definition: TPythia6.cxx:298
void SetOptStat(Int_t stat=1)
The type of information printed in the histogram statistics box can be selected via the parameter mod...
Definition: TStyle.cxx:1589
virtual int Load(const char *module, const char *entry="", Bool_t system=kFALSE)
Load a shared library.
Definition: TSystem.cxx:1855
A TTree represents a columnar dataset.
Definition: TTree.h:79
const Int_t n
Definition: legend1.C:16
Definition: file.py:1
Definition: tree.py:1