Pythia Tips

Cross check of the different Pythia available in the market.


-- JuanCastillo - 19 May 2008

Version, manuals and notes

PYTHIA web

Manuals/ Versions tested:

distribution language notes
PYTHIA version 5.720
Last date of change: 29 Nov 1995
JETSET version 7.408
Last date of change: 23 Aug 1995
FORTRAN Pythia manual v. 5.6
PYTHIA version 6.319
Last date of change: 1 Mar 2005
ROOT wrap over the FORTRAN code  
PYTHIA version 6.410
Last date of change: 30 Jan 2007
FORTRAN Release Notes about v 6.1 and later

Installation

FORTRAN code (version 5.720, version 6.410)

  • Create your destination folder and download the corresponding code, following the recommended instructions. Example: tar xzf pythia-6.4.10.tar.gz -C pylib
  • Test that you have all the necessary libraries (see "makefile" section part)

ROOT code

  • Create your PYTHIA folder (henceforth called pythia2)
  • Download a ROOT version (we use 514-00) and place it inside pythia2
  • Untar and unzip this ROOT. Example tar xzf root_v5.18.00.source.tar.gz -C root_v518
  • Install Pythia.
  • Check that you have the necessary environment variables (see next step)
  • Configure ROOT to be installed with Pythia (NOTE: all after configure must be in the same line)

./configure linux --enable-cern --enable-pythia6 
--enable-rfio --enable-mathmore  --enable-roofit 
--enable-asimage --enable-minuit2 --enable-shared 
--enable-xml --enable-soversion  
--with-cern-libdir=$CERNLIB 
--with-pythia6-libdir=$PYTHIA6 

  • Make install ROOT
  • Make ROOT

Generation code structure

First phase: initialization

The general character of the run is determined.

  • Incoming hadrons
  • Energies involved.

It is also possible to select specific final states, and decisions about details in the subsequent generation. This step is finished by a PYINIT call, at which time several variables are initialized in accordance with the values set.

Second phase: loop over events

Each new event generated by a PYEVNT call.
This event may then be analyzed, using information stored in some commonblocks, and the statistics accumulated. Example of common blocks:

COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200)

Purpose: to provide information on latest event generated or, in a few cases, on the accumulated statistics during the run.

COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5)

When an event is generated by a PYEVNT call, it is stored in the common block LUJETS.
Each parton and particle is represented by one line of information (so for 100 events = ). Scheme fixed: status code, parton/particle code, line of mother, lines of first and last daughter (or colour flow information in hard interaction and parton showers), momentum, energy, mass, production vertex and time, and invariant lifetime (if unstable).
Particle Data Group numbering conventions [PDG88].

Final phase: results

This may often be done without the invocation of any PYTHIA routines. From PYSTAT, however, it is possible to obtain a useful list of cross-sections for the different subprocesses.

FORTRAN makefile

We need the next libraries:

  • For Fortran:
FLIBFIL  =  ../flib/libacja.a ../flib/libe877.a
FLIB     = -L../flib -lacja -le877

  • For PAW:
PDFLIB = -L/usr/local/cern/2005-sarge-32/lib/ -lpdflib804 
JETLIB = -L/usr/local/cern/2005-sarge-32/lib/ -ljetset74
PACLIB = -L/usr/local/cern/2005-sarge-32/lib/ -lpacklib
MATLIB = -L/usr/local/cern/2005-sarge-32/lib/ -lmathlib
KERLIB = -L/usr/local/cern/2005-sarge-32/lib/ -lkernlib

FORTRAN variables

  • neven=50 //Number of events (collisions)
  • nmes=max(neven/10,1)
  • mstp(52)=2 //PDFs
  • mstp(51)=3031 //PDFs
  • mstp(51)=4053 //PDFs
  • npart = 9 //number of particles
  • pbeam = 14000.0 //beam energy

FORTRAN Subroutines

This may differ from version to version.
Changes from 5.6 to 6.1:
  • PYTHIA and JETSET merged
  • All JETSET subroutines renamed to begin with PY or contain it. Example:
    - LU = PY
    - UL = PY
    - RLU = PYR
    - KLU = PYK
    - PLU = PYP

Main FORTRAN subroutines

PYINIT(FRAME,BEAM,TARGET,WIN)

Purpose: to initialize the generation procedure.
FRAME : 'CMS','FIXT','USER','NONE'
BEAM, TARGET : 'p+' 'p~-' 'n0' 'n~0' 'pi+' 'pi-' 'e-' 'e+' 'nu_e' 'nu_e~' 'mu-' 'mu+' 'nu_mu' 'nu_mu~' 'gamma' 'Lambda0' 'Sigma-' 'Sigma0' 'Sigma+' 'Xi-' 'Xi0' 'Omega-'
WIN : related to energy of system, exact meaning depends on FRAME
Examples:
     CALL PYINIT ('CMS','p','n',pbeam)%BR%
     CALL PYINIT ('CMS','n','n',pbeam)%BR%
     CALL PYINIT ('CMS','p~-','p',pbeam)%BR%
     CALL PYINIT ('FIXT','p','p',pbeam)%BR%

PYEVNT

Purpose: to generate one event of the type specified by the PYINIT call. Called inside loops.

PYSTAT(MSTAT)

It prints out cross-sections statistics, decay widths, branching ratios, status codes and parameter values. PYSTAT may be called at any time, after the PYINIT call, e.g. at the end of the run, or not at all. MSTAT : specification of information desired.

  • 1 : event kind of and corresponding cross-sections. Numbers already including PYKCUT.
  • 2 : resonances, particle codes (KF), allowed decay channels.
  • 3 : hard interaction flavours KFIN(I,J) for beam and target particles.
  • 4 : kinematical cuts CKIN(I) in the current run.
  • 5 : status codes MSTP(I) and parameters PARP(I) used in the current run.

ISUB and MSEL

This codes control (using switch or if analogues) some physics processes. See documentation for more details.

  • QCD Jets. MSEL = 1,2 ISUB = 11,12,13,28,53,68
  • Heavy Flavours. MSEL=4,5,6,7,8,35,36,37,38. ISUB = 81,82,83,84,85
  • Minimum Bias. MSEL = 1,2. ISUB = 91,92,93,95
  • Prompt Photons. MSEL = 10. ISUB = 14,18,29,114,115
  • Single Z/W Production. MSEL = 11,12,13,14,(21). ISUB = 1,2,15,16,30,31,35,36,131,(141)
  • Neutral Higgs Production. MSEL = 16. ISUB = 3,5,8,24,26,(71,72,73,76,77),102,103,111,112,113,121,122,123,124
  • gamma/Z/W Pairs. MSEL = 15. ISUB = 19,20,22,23,25,69,70,71,72,73,76,77
  • New Gauge Bosons. MSEL = 21,22,24. ISUB = 141,142,144
  • Extended Higgs scenario. MSEL = 23. ISUB = 3,24,26,102,103,121,122,123,124,(141),143,151,152,153,156,157,158,161,171,172,173,174,176,177,178,179
  • Leptoquarks. MSEL = 25. ISUB = 145,162,163,164
  • Compositeness and Anomalous Couplings. ISUB = 11,12,20,165,166
  • Excited fermions. ISUB = 147,148

Analyzing PYTHIA ROOT output files

Converting output PAW files onto ROOT files

Source: ROOT manual

You can use h2root to convert your HBOOK/PAW histograms or ntuple files (output of the FORTRAN code) into ROOT files. To use this program, after login in ROOT, you type the shell script command:

h2root <hbookfile> <rootfile>

If you do not specify the second parameter, a file name is automatically generated for you. If hbookfile is of the form file.hbook, then the ROOT file will be called file.root.

In our case, ROOT trees are filled with leaves after the n-tupla

data chtags/'id','idpare','nrpare','idgran','nrgran','vx','vy','vz','px','py','pz','en','mass','rap'/

ROOT files from a PAW file in this way generated will be called from here "ROOT-translated".
ROOT-translated file unveiled:
ROOT-translated file unveiled

Using a TSelector

Source: ROOT manual

Open a ROOT session.
You know the name of the tree stored on in (if not, just browse inside the file using TBrowser). We call it h777.

  1. ) Load the file: TFile f("myroot.root")
  2. ) Create a selector h777->MakeSelector("h777analysis")
  3. ) Modify your selector
  4. ) Create a chain TChain chain2("h777")
  5. ) Load the file(s) onto the chain chain2.Add("myroot.root")
  6. ) Process the selector chain2.Process("h777analysis.C")

HOWTO: Modify your selector. A good idea is to create a TEventList ,a list of selected events (entries, that are maybe not corresponding to collisions but to particles) in a TTree.

Using macros

A chain must be created and filled.
Depending on if cuts are applied or not, the chain is read in a way or in another.
For cutting,usually TCut objects used. TCut in ROOT documentation

  • DisplayQuarks.C: Arguments: none. 1 ROOT file analyzed. 2 cuts applied. Plots drawed via TChain.Draw().Quark production is displayed. Quarks are identified via its PID label.
  • DisplayAll.C: Arguments: string filename. 1 ROOT file analyzed. No cuts applied. Plots drawed via TChain.Draw(). Draw("id"), Draw("eventnum"), Draw("dwinlow"), Draw("dwinup"), Draw("mass"), Draw("rap") in the same canva. Used for testing a proper filling.

  • MultDistAllFiles.C: Arguments: 1 file with filenames, output ROOT filename. For each file, it plots a multiplicity distribution using another function MultDist, that gives an average charged mutiplicity and its error. With the averages of all files, it plots a $N_{ch}$ vs ${\sqrt{s}}$ dependence.
  • AverMult.C: Arguments: number of files. It ask for each file, it plots a multiplicity distribution using another function MultDist, that prints onto screen $N_{ch}$. Version previous to MultDistAllFiles.C
  • MultDistOneFile.C: Arguments: 1 filename. A TTree is declared an branches are set from it.It plots a multiplicity distribution. Basically it is the function MultDist used above that returns an average charged mutiplicity and its error.
  • ParticleComp.C: Arguments: string fileone, string filetwo. Two chains are created from the ROOT files which names (or paths) are given as strings. After them, 2 TTrees and branches for reading PID of the stored particles. 3 plots are drawn: two for PID (one per file) and 1 for the ratio. Each plot is having a TPaveText with the relevant information. One output txt file with name of the particle and the compared production is created.
  • ParticleProd.C: Arguments: string filename. It's doing what the previous macro is doing but only for one file.
  • ProductionPIDs.C: Arguments: string filename.It plots in a canva particle production after PIDs for different sectors (9).

  • MomentumDistribution.C: Member content directly using chain.Draw("leave1:leave2") or similar. All events are read. No cuts are applied. The histograms are default ones. Used for tests on content of the file.
  • NumberOfEventsCut.C: Create a chain of ROOT-translated files and draw tree members. Cut on number of events. Histogram characteristics can be defined.
  • Multi1.C: Plots the number of particles produced in one event vs number of events with that number.

Hot Topics & Error messages

Documentation:

From ROOT forum.

  • Handling diffractive and elastic scattering events Pydiff
  • Search results for Pythia6 Pythia6
  • Illustrates the use of Pythia6 in ROOT macro

PDG code 9902210 unknown

Running for MB event generator, all the energies.

Warning in <TParticle::SetPdg>: PDG code 9902210 unknown from TDatabasePDG

TODO: Identificate the origin of this mistake (check if it's appearing in another tunings).

Solution: Cut on PID, not allowing it bigger than 500000.

Zero Energy-pz

This creates a crash while calculating the rapidity.
All the energies, but more probable @ high energies.

Ex: Event 4, (74 particles) at 12000.000 GeV center-of-mass energy
proton 0
uu_1 0
eta' 0
pi+ 0

TODO: Investigate the frequency of this events.

Solution: give an arbitrary big rapidity (10) to all the particles with Energy-pz = 0.
Topic attachments
I Attachment Action Size Date Who Comment
pythia_root_leaves.pngpng pythia_root_leaves.png manage 16.7 K 2008-03-04 - 16:59 JuanCastillo ROOT-translated file unveiled
Topic revision: r4 - 2009-10-21, JuanCastillo
 
This site is powered by FoswikiCopyright © by the contributing authors. All material on this collaboration platform is the property of the contributing authors.
Ideas, requests, problems regarding GSI Wiki? Send feedback
Imprint (in German)
Privacy Policy (in German)