LU TP 95–20

August 1995

(revised version of CERN–TH.7112/93)

(first issued December 1993)


[0.5] W5035/W5044

PYTHIA  5.7  and  JETSET  7.4

[5mm] Physics and Manual

[10mm] Torbjörn Sjöstrand

[3mm] Department of Theoretical Physics,

[1mm] University of Lund, Sölvegatan 14A,

[1mm] S-223 62 Lund, Sweden

            *......*        *:::!!:::::::::::*     *::::::!!::::::::::::::*   *::::::::!!::::::::::::::::*  *:::::::::!!:::::::::::::::::*  *:::::::::!!:::::::::::::::::*   *::::::::!!::::::::::::::::*!     *::::::!!::::::::::::::* !!     !! *:::!!:::::::::::*    !!     !!     !* -><- *         !!     !!     !!                !!     !!     !!                !!     !!                       !!     !!        ep             !!     !!                       !!     !!                 pp    !!     !!   e+e-                !!     !!                       !!     !!    

Important note: this is the long writeup of [1mm] T. Sjöstrand, Computer Physics Commun. 82 (1994) 74. [1mm] All references should be to the published version.

Copyright Notice CERNLIB – CERN Program Library Long writeups © Copyright CERN, Geneva 1993 Copyright and any other appropriate legal protection of these computer programs and associated documentation reserved in all countries of the world. These programs or documentation may not be reproduced by any method without prior written consent of the Director-General of CERN or his delegate. Permission for the usage of any programs described herein is granted apriori to those scientific institutes associated with the CERN experimental program or with whom CERN has concluded a scientific collaboration agreement. Requests for information should be addressed to: CERN Program Library Office CERN-CN Division CH-1211 Geneva 23 Switzerland Tel. +41 22 767 4951 Fax. +41 22 767 7155 Bitnet: [email protected] DECnet: VXCERN::CERNLIB (node 22.190) Internet: [email protected]

Trademark notice: All trademarks appearing in this guide are acknowledged as such.


The Pythia and Jetset programs are frequently used for event generation in high-energy physics. The emphasis is on multiparticle production in collisions between elementary particles. This in particular means hard interactions in , and colliders, although also other applications are envisaged. The programs are intended to generate complete events, in as much detail as experimentally observable ones, within the bounds of our current understanding of the underlying physics. Many of the components of the programs represent original research, in the sense that models have been developed and implemented for a number of aspects not covered by standard theory. Although originally conceived separately, the Pythia and Jetset programs today are so often used together that it makes sense to present them here without too much distinction.

Both programs have a long history, and several manuals have come out. The former round of Pythia/Jetset program descriptions appeared in 1987. Meanwhile a large number of additions and changes have been made. Recently a new description therefore appeared in
This is the one and only correct reference to the current versions of Pythia and Jetset. The long writeup that you now have before you is an (unpublished) appendix to the publication above, and need not be separately cited. Instead remember to cite the original literature on the physics topics of particular relevance for your studies. (There is no reason to omit references to good physics papers simply because some of their contents have also been made available as program code.)

Event generators often have a reputation for being ‘black boxes’; if nothing else, this report should provide you with a glimpse of what goes on inside the programs. Some such understanding may be of special interest for new users, who have no background in the field. An attempt has been made to structure the report sufficiently well that many of the sections can be read independently of each other, so you can pick the sections that interest you. I have tried to keep together the physics and the manual sections on specific topics, where practicable, which represents a change of policy compared with previous manual versions. Any feedback on this and other aspects is welcome.

A large number of persons should be thanked for their contributions. Hans-Uno Bengtsson is the originator of the Pythia program, and for many years we worked in parallel on its further development. Mats Bengtsson is the main author of the final-state parton-shower algorithm. Bo Andersson and Gösta Gustafson are the originators of the Lund model, and strongly influenced the early development of the programs. Further comments on the programs have been obtained from users too numerous to be mentioned here, but who are all gratefully acknowledged. To write programs of this size and complexity would be impossible without a strong user feedback.

The moral responsibility for any remaining errors clearly rests with me. However, kindly note that this is a ‘University World’ product, distributed ‘as is’, free of charge, without any binding guarantees. And always remember that the programs do not represent a dead collection of established truths, but rather one of many possible approaches to the problem of multiparticle production in high-energy physics, at the frontline of current research. Be critical!


1 Introduction

Multiparticle production is the most characteristic feature of current high-energy physics. Today, observed particle multiplicities are typically between ten and a hundred, and with future machines this range will be extended upwards. The bulk of the multiplicity is found in jets, i.e. in bunches of hadrons (or decay products of hadrons) produced by the hadronization of quarks and gluons.

The Complexity of High-Energy Processes

To first approximation, all processes have a simple structure at the level of interactions between the fundamental objects of nature, i.e. quarks, leptons and gauge bosons. For instance, a lot can be understood about the structure of hadronic events at LEP just from the ‘skeleton’ process . Corrections to this picture can be subdivided, arbitrarily but conveniently, into three main classes.

Firstly, there are bremsstrahlung-type modifications, i.e. the emission of additional final-state particles by branchings such as or . Because of the largeness of the strong coupling constant , and because of the presence of the triple gluon vertex, QCD emission off quarks and gluons is especially prolific. We therefore speak about ‘parton showers’, wherein a single initial parton may give rise to a whole bunch of partons in the final state. Also photon emission may give sizeable effects in and processes. The bulk of the bremsstrahlung corrections are universal, i.e. do not depend on the details of the process studied, but only on one or a few key numbers, such as the momentum transfer scale of the process. Such universal corrections may be included to arbitrarily high orders, using a probabilistic language. Alternatively, exact calculations of bremsstrahlung corrections may be carried out order by order in perturbation theory, but rapidly the calculations then become prohibitively complicated and the answers correspondingly lengthy.

Secondly, we have ‘true’ higher-order corrections, which involve a combination of loop graphs and the soft parts of the bremsstrahlung graphs above, a combination needed to cancel some divergences. In a complete description it is therefore not possible to consider bremsstrahlung separately, as assumed here. The necessary perturbative calculations are usually very difficult; only rarely have results been presented that include more than one non-‘trivial’ order, i.e. more than one loop. As above, answers are usually very lengthy, but some results are sufficiently simple to be generally known and used, such as the running of , or the correction factor in the partial widths of decay channels. For high-precision studies it is imperative to take into account the results of loop calculations, but usually effects are minor for the qualitative aspects of high-energy processes.

Thirdly, quarks and gluons are confined. In the two points above, we have used a perturbative language to describe the short-distance interactions of quarks, leptons and gauge bosons. For leptons and colourless bosons this language is sufficient. However, for quarks and gluons it must be complemented with a picture for the hadronization process (which can be subdivided into fragmentation and decays), wherein the coloured partons are transformed into jets of colourless hadrons, photons and leptons. This process is still not yet understood from first principles, but has to be based on models. In one sense, hadronization effects are overwhelmingly large, since this is where the bulk of the multiplicity comes from. In another sense, the overall energy flow of a high-energy event is mainly determined by the perturbative processes, with only a minor additional smearing caused by the hadronization step. One may therefore pick different levels of ambition, but in general detailed studies require a detailed modelling of the hadronization process.

The simple structure that we started out with has now become considerably more complex — instead of maybe two final-state partons we have a hundred final particles. The original physics is not gone, but the skeleton process has been dressed up and is no longer directly visible. A direct comparison between theory and experiment is therefore complicated at best, and impossible at worst.

Event Generators

It is here that event generators come to the rescue. In an event generator, the objective strived for is to use computers to generate events as detailed as could be observed by a perfect detector. This is not done in one step, but rather by ‘factorizing’ the full problem into a number of components, each of which can be handled reasonably accurately. Basically, this means that the hard process is used as input to generate bremsstrahlung corrections, and that the result of this exercise is thereafter left to hadronize. This sounds a bit easier than it really is — else this report would be a lot thinner. However, the basic idea is there: if the full problem is too complicated to be solved in one go, try to subdivide it into smaller tasks of manageable proportions. In the actual generation procedure, most steps therefore involve the branching of one object into two, or at least into a very small number, each of which being free to branch in its turn. A lot of bookkeeping is involved, but much is of a repetitive nature, and can therefore be left for the computer to handle.

As the name indicates, the output of an event generator should be in the form of ‘events’, with the same average behaviour and the same fluctuations as real data. In the data, fluctuations arise from the quantum mechanics of the underlying theory. In generators, Monte Carlo techniques are used to select all relevant variables according to the desired probability distributions, and thereby ensure randomness in the final events. Clearly some loss of information is entailed: quantum mechanics is based on amplitudes, not probabilities. However, only very rarely do (known) interference phenomena appear that cannot be cast in a probabilistic language. This is therefore not a more restraining approximation than many others.

Once there, an event generator can be used in many different ways. The five main applications are probably the following:

  • To give physicists a feeling for the kind of events one may expect/hope to find, and at what rates.

  • As a help in the planning of a new detector, so that detector performance is optimized, within other constraints, for the study of interesting physics scenarios.

  • As a tool for devising the analysis strategies that should be used on real data, so that signal-to-background conditions are optimized.

  • As a method for estimating detector acceptance corrections that have to be applied to raw data, in order to extract the ‘true’ physics signal.

  • As a convenient framework within which to interpret the observed phenomena in terms of a more fundamental underlying theory (usually the Standard Model).

Where does a generator fit into the overall analysis chain of an experiment? In ‘real life’, the machine produces interactions. These events are observed by detectors, and the interesting ones are written to tape by the data acquisition system. Afterwards the events may be reconstructed, i.e. the electronics signals (from wire chambers, calorimeters, and all the rest) may be translated into a deduced setup of charged tracks or neutral energy depositions, in the best of worlds with full knowledge of momenta and particle species. Based on this cleaned-up information, one may proceed with the physics analysis. In the Monte Carlo world, the rôle of the machine, namely to produce events, is taken by the event generators described in this report. The behaviour of the detectors — how particles produced by the event generator traverse the detector, spiral in magnetic fields, shower in calorimeters, or sneak out through cracks, etc. — is simulated in programs such as Geant [Bru89]. Traditionally, this latter activity is called event simulation, which is somewhat unfortunate since the same words could equally well be applied to what, here, we call event generation. A more appropriate term is detector simulation. Ideally, the output of this simulation has exactly the same format as the real data recorded by the detector, and can therefore be put through the same event reconstruction and physics analysis chain, except that here we know what the ‘right answer’ should be, and so can see how well we are doing.

Since the full chain of detector simulation and event reconstruction is very time-consuming, one often does ‘quick and dirty’ studies in which these steps are skipped entirely, or at least replaced by very simplified procedures which only take into account the geometric acceptance of the detector and other trivial effects. One may then use the output of the event generator directly in the physics studies.

There are still many holes in our understanding of the full event structure, despite an impressive amount of work and detailed calculations. To put together a generator therefore involved making a choice on what to include, and how to include it. At best, the spread between generators can be used to give some impression of the uncertainties involved. A multitude of approximations will be discussed in the main part of this report, but already here is should be noted that many major approximations are related to the almost complete neglect of the second point above, i.e. of the non-‘trivial’ higher-order effects. It can therefore only be hoped that the ‘trivial’ higher order parts give the bulk of the experimental behaviour. By and large, this seems to be the case; for annihilation it even turns out to be a very good approximation.

The necessity to make compromises has one major implication: to write a good event generator is an art, not an exact science. It is therefore essential not to blindly trust the results of any single event generator, but always to make several cross-checks. In addition, with computer programs of tens of thousands of lines, the question is not whether bugs exist, but how many there are, and how critical their positions. Further, an event generator cannot be thought of as all-powerful, or able to give intelligent answers to ill-posed questions; sound judgement and some understanding of a generator are necessary prerequisites for successful use. In spite of these limitations, the event generator approach is the most powerful tool at our disposal if we wish to gain a detailed and realistic understanding of physics at current or future high-energy colliders.

The Origins of the JETSET and PYTHIA Programs

Over the years, many event generators have appeared. Surveys of generators for physics in general and LEP in particular may be found in [Kle89, Sjö89], for high-energy hadron–hadron () physics in [Ans90, Sjö92, Kno93], and for physics in [HER92]. We refer the reader to those for additional details and references. In this particular report, the two closely connected programs Jetset and Pythia will be described.

Jetset has its roots in the efforts of the Lund group to understand the hadronization process, starting in the late seventies [And83]. The so-called string fragmentation model was developed as an explicit and detailed framework, within which the long-range confinement forces are allowed to distribute the energies and flavours of a parton configuration among a collection of primary hadrons, which subsequently may decay further. This model, known as the Lund string model, or ‘Lund’ for short, contained a number of specific predictions, which were confirmed by data from PETRA and PEP, whence the model gained a widespread acceptance. The Lund string model is still today the most elaborate and widely used fragmentation model at our disposal. It remains at the heart of the Jetset/Pythia programs.

In order to predict the shape of events at PETRA/PEP, and to study the fragmentation process in detail, it was necessary to start out from the partonic configurations that were to fragment. The generation of complete hadronic events was therefore added, originally based on simple exchange and first-order QCD matrix elements, later extended to full exchange with first-order initial-state QED radiation and second-order QCD matrix elements. A number of utility routines were also provided early on, for everything from event listing to jet finding.

By the mid-eighties it was clear that the matrix-element approach had reached the limit of its usefulness, in the sense that it could not fully describe the multijet topologies of the data. (Later on, the use of optimized perturbation theory was to lead to a resurgence of the matrix-element approach, but only for specific applications.) Therefore a parton-shower description was developed [Ben87a] as an alternative to the matrix-element one. The combination of parton showers and string fragmentation has been very successful, and forms the main approach to the description of hadronic events.

In recent years, Jetset has been a fairly stable product, covering the four main areas of fragmentation, final-state parton showers, event generation and general utilities.

The successes of string fragmentation in made it interesting to try to extend this framework to other processes, and explore possible physics consequences. Therefore a number of other programs were written, which combined a process-specific description of the hard interactions with the general fragmentation framework of Jetset. The Pythia program evolved out of early studies on fixed-target proton–proton processes, addressed mainly at issues related to string drawing.

With time, the interest shifted towards higher energies, first to the SPS collider, and later to SSC and LHC, in the context of a number of workshops in the USA and Europe. Parton showers were added, for final-state radiation by making use of the Jetset routine, for initial-state one by the development of the concept of ‘backwards evolution’, specifically for Pythia [Sjö85]. Also a framework was developed for minimum-bias and underlying events [Sjö87a].

Another main change was the introduction of an increasing number of hard processes, within the Standard Model and beyond. A special emphasis was put on the search for the Standard Model Higgs, in different mass ranges and in different channels, with due respect to possible background processes.

The bulk of the machinery developed for hard processes actually depended little on the choice of initial state, as long as the appropriate parton distributions were there for the incoming partons and particles. It therefore made sense to extend the program from being only a generator to working also for and . This process was only completed in 1991, again spurred on by physics workshop activities. Currently Pythia should therefore work equally well for a selection of different possible incoming beam particles.

The tasks of including new processes, and of improving the simulation of already present ones, are never-ending. Work therefore continues apace.

While Jetset still is formally independent of Pythia, their ties have grown much stronger over the years, and the border-line between the two programs has become more and more artificial. It is no coincidence that the two are presented together here; this way a lot of repetition of common material can be avoided. The price to be paid is that some differences in philosophy will have to be discussed.

About this Report

As we see, Jetset and Pythia started out as very ideologically motivated programs, developed to study specific physics questions in enough detail that explicit predictions could be made for experimental quantities. As it was recognized that experimental imperfections could distort the basic predictions, the programs were made available for general use by experimentalists. It thus became feasible to explore the models in more detail than would otherwise have been possible. As time went by, the emphasis came to shift somewhat, away from the original strong coupling to a specific fragmentation model, towards a description of high-energy multiparticle production processes in general. Correspondingly, the use expanded from being one of just comparing data with specific model predictions, to one of extensive use for the understanding of detector performance, for the derivation of acceptance correction factors, for the prediction of physics at future high-energy accelerators, and for the design of related detectors.

While the ideology may be less apparent, it is still there, however. This is not something unique to the programs discussed here, but inherent in any event generator, or at least any generator that attempts to go beyond the simple parton level skeleton description of a hard process. Do not accept the myth that everything available in Monte Carlo form represents ages-old common knowledge, tested and true. Ideology is present by commissions or omissions in any number of details. Programs like Pythia and Jetset represent a major amount of original physics research, often on complicated topics where no simple answers are available. As a (potential) program user you must be aware of this, so that you can form your own opinion, not just about what to trust and what not to trust, but also how much to trust a given prediction, i.e. how uncertain it is likely to be. Jetset and Pythia are particularly well endowed in this respect, since a number of publications exist where most of the relevant physics is explained in considerable detail. In fact, the problem may rather be the opposite, to find the relevant information among all the possible places. One main objective of the current report is therefore to collect much of this information in one single place. Not all the material found in specialized papers is reproduced, by a wide margin, but at least enough should be found here to understand the general picture and to know where to go for details.

The current report is therefore intended to replace the previous round of published physics descriptions and program manuals [Sjö86, Sjö87, Ben87]. The formal new standard reference is [Sjö94], which is a fairly brief summary of this report — for obvious reasons the full description is too long to be published in its entirety. Further specification could include a statement of the type ‘We use Pythia version X.x and Jetset version Y.y’. (If you are a LaTeX fan, you may want to know that the program names in this report have been generated by the commands \textsc{Jetset} and \textsc{Pythia}.) Kindly do not refer to Jetset/Pythia as ‘unpublished’, ‘private communication’ or ‘in preparation’: such phrases are only creating unnecessary confusion.

In addition, remember that many of the individual physics components are documented in separate publications. If some of these contain ideas that are useful to you, there is every reason to cite them. A reasonable selection would vary as a function of the physics you are studying. The criterion for which to pick should be simple: imagine that a Monte Carlo implementation had not been available. Would you then have cited a given paper on the grounds of its physics contents alone? If so, do not punish the extra effort of turning these ideas into publicly available software. (Monte Carlo manuals are good for nothing in the eyes of many theorists, so often only the acceptance of ‘mainstream’ publications counts.) Here follows a list of some main areas where the Pythia/Jetset programs contain original research:

  • The string fragmentation model [And83].

  • The string effect [And80].

  • Baryon production (diquark/popcorn) [And82, And85].

  • Fragmentation of multiparton systems [Sjö84].

  • Fragmentation effects on determinations [Sjö84a].

  • Initial state parton showers [Sjö85].

  • Final state parton showers [Ben87a].

  • Photon radiation from quarks [Sjö92c]

  • Deep inelastic scattering [And81a, Ben88].

  • Photoproduction [Sch93a] and physics [Sch94a].

  • Parton distributions of the photon [Sch95].

  • Colour flow in hard scatterings [Ben84].

  • Elastic and diffractive cross sections [Sch94].

  • Minijets (multiple parton–parton interactions) [Sjö87a].

  • Rapidity gaps [Dok92].

  • Jet clustering in [Sjö83].

In addition to a physics survey, the current report also contains a complete manual for the two programs. Such manuals have always been updated and distributed jointly with the programs. To a first approximation, we therefore do not have much new to offer here. However, an attempt has been made to group the material more logically according to physics topics than in previous distributions, to tie it closer to the physics description, and to improve the layout and therefore the readability. Any feedback is welcome.

A word of warning may be in place. The program description is fairly lengthy, and certainly could not be absorbed in one sitting. This is not even necessary, since all switches and parameters are provided with sensible default values, based on our best understanding (of the physics, and of what you expect to happen if you do not specify any options). As a new user, you can therefore disregard all the fancy options, and just run the program with a minimum ado. Later on, as you gain experience, the options that seem useful can be tried out. No single user is ever likely to find need for more than a fraction of the total number of possibilities available, yet many of them have been added to meet specific user requests.

In some instances, not even this report will provide you with all the information you desire. You may wish to find out about recent versions of the program, know about related software, pick up a few sample main programs to get going, or get hold of related physics papers. Some such material can be found if you link to my World Wide Web homepage:
and study the contents there.


At all times it should be remembered that this is not a commercial product, developed and supported by professionals. Instead it is a ‘University World’ product, developed by a very few physicists (mainly the current author) originally for their own needs, and supplied to other physicists on an ‘as-is’ basis, free of charge. No guarantees are therefore given for the proper functioning of the programs, nor for the validity of physics results. In the end, it is always up to you to decide for yourself whether to trust a given result or not. Usually this requires comparison either with analytical results or with results of other programs, or with both. Even this is not necessarily foolproof: for instance, if an error is made in the calculation of a matrix element for a given process, this error will be propagated both into the analytical results based on the original calculation and into all the event generators which subsequently make use of the published formulae. In the end, there is no substitute for a sound physics judgement.

This does not mean that you are all on your own, with a program nobody feels responsible for. Attempts are made to check processes as carefully as possible, to write programs that do not invite unnecessary errors, and to provide a detailed and accurate documentation. All of this while maintaining the full power and flexibility, of course, since the physics must always take precedence in any conflict of interests. If nevertheless any errors or unclarities are found, please do communicate them to me, e.g. on phone +46 – 46 – 222 48 16 or e-mail . Every attempt will be made to solve problems as soon as is reasonably possible, given that this support is by one person alone, who also has other responsibilities.

Appendix: The Historical Pythia

While the origin and connotations of the ‘Jetset’ program name should be commonly known, the ‘Pythia’ label may need some explanation.

The myth tells how Apollon, the God of Wisdom, killed the powerful dragon-like monster Python, close to the village of Delphi in Greece. To commemorate this victory, Apollon founded the Pythic Oracle in Delphi, on the slopes of Mount Parnassos. Here men could come to learn the will of the Gods and the course of the future. The oracle plays an important rôle in many of the other Greek myths, such as those of Heracles and of King Oedipus.

Questions were to be put to the Pythia, the ‘Priestess’ or ‘Prophetess’ of the Oracle. In fact, she was a local woman, usually a young maiden, of no particular religious schooling. Seated on a tripod, she inhaled the obnoxious vapours that seeped up through a crevice in the ground. This brought her to a trance-like state, in which she would scream seemingly random words and sounds. It was the task of the professional priests in Delphi to record those utterings and edit them into the official Oracle prophecies, which often took the form of poems in perfect hexameter. In fact, even these edited replies were often less than easy to interpret. The Pythic oracle acquired a reputation for ambiguous answers.

The Oracle existed already at the beginning of the historical era in Greece, and was universally recognized as the foremost religious seat. Individuals and city states came to consult, on everything from cures for childlessness to matters of war. Lavish gifts allowed the temple area to be built and decorated. Many states supplied their own treasury halls, where especially beautiful gifts were on display. Sideshows included the Omphalos, a stone reputedly marking the centre of the Earth, and the Pythic games, second only to the Olympic ones in importance.

Strife inside Greece eventually led to a decline in the power of the Oracle. A serious blow was dealt when the Oracle of Zeus Ammon (see below) declared Alexander the Great to be the son of Zeus. The Pythic Oracle lived on, however, and was only closed by a Roman Imperial decree in 390 ad, at a time when Christianity was ruthlessly destroying any religious opposition. Pythia then had been at the service of man and Gods for a millenium and a half.

The rôle of the Pythic Oracle replies on the course of history is nowhere better described than in ‘The Histories’ by Herodotus [Herbc], the classical and captivating description of the Ancient World at the time of the Great War between Greeks and Persians. Especially famous is the episode with King Croisus of Lydia. Contemplating a war against the upstart Persian Empire, he resolves to ask an oracle what the outcome of a potential battle would be. However, to have some guarantee for the veracity of any prophecy, he decides to send embassies to all the renowned oracles of the known World. The messengers are instructed to inquire the various divinities, on the hundredth day after their departure, what King Croisus is doing at that very moment. From the Pythia the messengers bring back the reply

I know the number of grains of sand as well as the expanse of the sea,
And I comprehend the dumb and hear him who does not speak,
There came to my mind the smell of the hard-shelled turtle,
Boiled in copper together with the lamb,
With copper below and copper above.

The veracity of the Pythia is thus established by the crafty ruler, who had waited until the appointed day, slaughtered a turtle and a lamb, and boiled them together in a copper cauldron with a copper lid. Also the Oracle of Zeus Ammon in the Libyan desert is able to give a correct reply (lost to posterity), while all others fail. King Croisus now sends a second embassy to Delphi, inquiring after the outcome of a battle against the Persians. The Pythia answers

If Croisus passes over the Halys he will dissolve a great Empire.

Taking this to mean he would win, the King collects his army and crosses the border river, only to suffer a crushing defeat and see his Kingdom conquered. When the victorious King Cyrus allows Croisus to send an embassy to upbraid the Oracle, the God Apollon answers through his Prophetess that he has correctly predicted the destruction of a great empire — Croisus’ own — and that he cannot be held responsible if people choose to interpret the Oracle answers to their own liking.

The history of the Pythia program is neither as long nor as dignified as that of its eponym. However, some points of contact exist. You must be very careful when you formulate the questions: any ambiguities will corrupt the reply you get. And you must be even more careful not to misinterpret the answers; in particular not to pick the interpretation that suits you before considering the alternatives. Finally, even a perfect God has servants that are only human: a priest might mishear the screams of the Pythia and therefore produce an erroneous oracle reply; the current author might unwittingly let a bug free in the program Pythia.

2 Physics Overview

In this section we will try to give an overview of the main physics features of Jetset and Pythia, and also to introduce some terminology. The details will be discussed in subsequent sections.

For the description of a typical high-energy event, an event generator should contain a simulation of several physics aspects. If we try to follow the evolution of an event in some semblance of a time order, one may arrange these aspects as follows:

  • Initially two beam particles are coming in towards each other. Normally each particle is characterized by a set of parton distribution functions, which defines the partonic substructure in terms of flavour composition and energy sharing.

  • One shower initiator parton from each beam starts off a sequence of branchings, such as , which build up an initial-state shower.

  • One incoming parton from each of the two showers enters the hard process, where then a number of outgoing partons are produced, usually two. It is the nature of this process that determines the main characteristics of the event.

  • Also the outgoing partons may branch, to build up final-state showers.

  • When a shower initiator is taken out of a beam particle, a beam remnant is left behind. This remnant may have an internal structure, and a net colour charge that relates it to the rest of the final state.

  • The QCD confinement mechanism ensures that the outgoing quarks and gluons are not observable, but instead fragment to colour neutral hadrons.

  • Many of the produced hadrons are unstable and decay further.

Conventionally, only quarks and gluons are counted as partons, while leptons and photons are not. If pushed ad absurdum this may lead to some unwieldy terminology. We will therefore, where it does not matter, speak of an electron or a photon in the ‘partonic’ substructure of an electron, lump branchings together with other ‘parton shower’ branchings such as , and so on. With this notation, the division into the above seven points applies equally well to an interaction between two leptons, between a lepton and a hadron, and between two hadrons.

In the following subsections, we will survey the above seven aspects, not in the same order as given here, but rather in the order in which they appear in the program execution, i.e. starting with the hard process.

2.1 Hard Processes and Parton Distributions

In Jetset, only two hard processes are available. The first and main one is . Here the ‘’ of is used to denote that the photon must be off the mass shell. The distinction is of some importance, since a photon on the mass shell cannot decay. Of course also the can be off the mass shell, but here the distinction is less relevant (strictly speaking, a is always off the mass shell). In the following we may not always use ‘’ consistently, but the rule of thumb is to use a ‘’ only when a process is not kinematically possible for a particle of nominal mass. The quark in the final state of may be , , , , or ; the flavour in each event is picked at random, according to the relative couplings, evaluated at the hadronic c.m. energy. Also the angular distribution of the final pair is included. No parton-distribution functions are needed.

The other Jetset process is a routine to generate and final states, as expected in onium 1 decays such as . Given the current limits on the top mass, toponium will decay weakly much too fast for these processes to be of any interest, so therefore no new applications are expected.

Pythia contains a much richer selection, with close to a hundred different hard processes. These may be classified in many different ways.

One is according to the number of final-state objects: we speak of ‘’ processes, ‘’ ones, ‘’ ones, etc. This aspect is very relevant from a programming point of view: the more particles in the final state, the more complicated the phase space and therefore the whole generation procedure. In fact, Pythia is optimized for and processes. There is currently no generic treatment of processes with three or more particles in the final state, but rather a few different machineries, each tailored to the pole structure of a specific class of graphs. This may be seen as a major limitation, and indeed is so at times. However, often one can come quite far with only one or two particles in the final state, since showers will add the required extra activity. The classification may also be misleading at times, since an -channel resonance is considered as a single particle, even if it is assumed always to decay into two final-state particles. Thus the process is classified as , although the decay treatment of the pair includes the full matrix elements.

Another classification is according to the physics scenario. This will be the main theme of section 8. The following major groups may be distinguished:

  • Hard QCD processes, e.g. .

  • Soft QCD processes, such as diffractive and elastic scattering, and minimum-bias events.

  • Heavy-flavour production, e.g. .

  • Prompt-photon production, e.g. .

  • Photon-induced processes, e.g. .

  • Deep inelastic scattering, e.g. .

  • production, such as the already found in Jetset, or .

  • Standard model Higgs production, where the Higgs is reasonably light and narrow, and can therefore still be considered as a resonance.

  • Gauge boson scattering processes, such as , when the Standard Model Higgs is so heavy and broad that resonant and non-resonant contributions have to be considered together.

  • Non-standard Higgs particle production, within the framework of a two-Higgs-doublet scenario with three neutral and two charged Higgs states.

  • Production of new gauge bosons, such as a .

  • Production of fourth-generation fermions.

  • Leptoquark production.

  • Deviations from Standard Model processes, e.g. due to contact interactions or a strongly interacting gauge boson sector. These scenarios do not always appear as separate processes, but may just be options to some of the processes above.

This is by no means a survey of all interesting physics. Most notable is the absence of supersymmetric particle production and decay, but many other examples could be found. Also, within the scenarios studied, not all contributing graphs have always been included, but only the more important and/or more interesting ones. In many cases, various approximations are involved in the matrix elements coded.

The cross section for a given process is given by


Here is the cross section for the hard partonic process, as codified in the matrix elements for each specific process. For processes with many particles in the final state it would be replaced by an integral over the allowed final-state phase space. The are the parton-distribution functions, which describe the probability to find a parton inside beam particle , with parton carrying a fraction of the total momentum. Actually, parton distributions also depend on some momentum scale that characterizes the hard process.

Parton distributions are most familiar for hadrons, such as the proton. Hadrons are inherently composite objects, made up of quarks and gluons. Since we do not understand QCD, a derivation from first principles of hadron parton distributions does not yet exist, although some progress is being made in lattice QCD studies. It is therefore necessary to rely on parametrizations, where experimental data are used in conjunction with the evolution equations for the dependence, to pin down the parton distributions. Several different groups have therefore produced their own fits, based on slightly different sets of data, and with some variation in the theoretical assumptions.

Also for fundamental particles, such as the electron, is it convenient to introduce parton distributions. The function thus parametrizes the probability that the electron that takes part in the hard process retains a fraction of the original energy, the rest being radiated (into photons) in the initial state. Of course, such radiation could equally well be made part of the hard interaction, but the parton-distribution approach usually is much more convenient. If need be, a description with fundamental electrons is recovered for the choice . Note that, contrary to the proton case, electron parton distributions are calculable from first principles, and reduce to the function above for .

The electron may also contain photons, and the photon may in its turn contain quarks and gluons. The internal structure of the photon is a bit of a problem, since the photon contains a point-like part, which is perturbatively calculable, and a vector-meson dominance part, which is not. Normally, the photon parton distributions are therefore parametrized, just as the hadron ones. Since the electron ultimately contains quarks and gluons, hard QCD processes like therefore not only appear in collisions, but also in ones (‘resolved photoproduction’) and in ones (‘doubly resolved 2 events’). The parton distribution function approach here makes it much easier to reuse one and the same hard process in different contexts.

There is also another kind of possible generalization. The two processes , studied in hadron colliders, and , studied in colliders, are really special cases of a common process, , where denotes a fundamental fermion, i.e. a quark, lepton or neutrino. The whole structure is therefore only coded once, and then slightly different couplings and colour prefactors are used, depending on the initial state considered. Usually the interesting cross section is a sum over several different initial states, e.g. and in a hadron collider. This kind of summation is always implicitly done, even when not explicitly mentioned in the text.

2.2 Initial- and Final-State Radiation

In every process that contains coloured and/or charged objects in the initial or final state, gluon and/or photon radiation may give large corrections to the overall topology of events. Starting from a basic process, this kind of corrections will generate , , and so on, final-state topologies. As the available energies are increased, hard emission of this kind is increasingly important, relative to fragmentation, in determining the event structure.

Two traditional approaches exist to the modelling of perturbative corrections. One is the matrix-element method, in which Feynman diagrams are calculated, order by order. In principle, this is the correct approach, which takes into account exact kinematics, and the full interference and helicity structure. The only problem is that calculations become increasingly difficult in higher orders, in particular for the loop graphs. Only in exceptional cases have therefore more than one loop been calculated in full, and often we do not have any loop corrections at all at our disposal. On the other hand, we have indirect but strong evidence that, in fact, the emission of multiple soft gluons plays a significant rôle in building up the event structure, e.g. at LEP, and this sets a limit to the applicability of matrix elements. Since the phase space available for gluon emission increases with the available energy, the matrix-element approach becomes less relevant for the full structure of events at higher energies. However, the perturbative expansion by itself is better behaved at higher energies, owing to the running of . As a consequence, inclusive measurements, e.g. of the rate of well-separated jets, should yield more reliable results.

The second possible approach is the parton-shower one. Here an arbitrary number of branchings of one parton into two (or more) may be put together, to yield a description of multijet events, with no explicit upper limit on the number of partons involved. This is possible since the full matrix-element expressions are not used, but only approximations derived by simplifying the kinematics, and the interference and helicity structure. Parton showers are therefore expected to give a good description of the substructure of jets, but in principle the shower approach has limited predictive power for the rate of well-separated jets (i.e. the 2/3/4/5-jet composition). In practice, shower programs may be patched up to describe the hard-gluon emission region reasonably well, in particular for the annihilation process. Nevertheless, the shower description is not optimal for absolute determinations.

Thus the two approaches are complementary in many respects, and both have found use. However, because of its simplicity and flexibility, the parton-shower option is generally the first choice, while the matrix elements one is mainly used for determinations, angular distribution of jets, triple-gluon vertex studies, and other specialized studies. Obviously, the ultimate goal would be to have an approach where the best aspects of the two worlds are harmoniously married.

2.2.1 Matrix elements

Matrix elements are especially made use of in the Jetset implementation of the process .

For initial-state QED radiation, a first order (unexponentiated) description has been adopted. This means that events are subdivided into two classes, those where a photon is radiated above some minimum energy, and those without such a photon. In the latter class, the soft and virtual corrections have been lumped together to give a total event rate that is correct up to one loop. This approach worked fine at PETRA/PEP energies, but does not do so well for the line shape, i.e. in regions where the cross section is rapidly varying and high precision is strived for.

For final-state QCD radiation, several options are available. The default is the parton-shower one (see below), but the matrix-elements options are also frequently used. In the definition of 3- or 4-jet events, a cut is introduced whereby it is required that any two partons have an invariant mass bigger than some fraction of the c.m. energy. 3-jet events which do not fulfill this requirement are lumped with the 2-jet ones. The first-order matrix-element option, which only contains 3- and 2-jet events therefore involves no ambiguities. In second order, where also 4-jets have to be considered, a main issue is what to do with 4-jet events that fail the cuts. Depending on the choice of recombination scheme, whereby the two nearby partons are joined into one, different 3-jet events are produced. Therefore the second-order differential 3-jet rate has been the subject of some controversy, and Jetset actually contains two different implementations.

By contrast, Pythia does not contain any full higher-order matrix elements, with loop contributions included. There are a few cases where higher-order matrix elements are included at the Born level. Consider e.g. the case of production at a hadron collider, which is contained in the lowest-order process . In an inclusive description, additional jets recoiling against the may be generated by parton showers. Pythia also contains the two first-order processes and . The cross sections for these processes are divergent when the . In this region a correct treatment would therefore have to take into account loop corrections, which are not available in Pythia. Depending on the physics application, one could then use Pythia in one of two ways. In the region of small , the preferred option is lowest-order matrix elements combined with parton showers. For the production of a at large , on the other hand, the shower approach is too imprecise to give the right cross section; additionally the event selection machinery is very inefficient. Here it is advantageous to generate first-order events, and then add showers only to describe additional softer radiation.

2.2.2 Parton showers

The separation of radiation into initial- and final-state showers is arbitrary, but very convenient. There are also situations where it is appropriate: for instance, the process only contains final-state QCD radiation (QED radiation, however, is possible both in the initial and final state), while only contains initial-state QCD one. Similarly, the distinction of emission as coming either from the or from the is arbitrary. In general, the assignment of radiation to a given mother parton is a good approximation for an emission close to the direction of motion of that parton, but not for the wide-angle emission in between two jets, where interference terms are expected to be important.

In both initial- and final-state showers, the structure is given in terms of branchings , specifically , , , , and . Each of these processes is characterized by a splitting kernel . The branching rate is proportional to the integral . The value picked for a branching describes the energy sharing, with daughter taking a fraction and daughter the remaining of the energy. Once formed, the daughters and may in turn branch, and so on.

Each parton is characterized by some virtuality scale , which gives an approximate sense of time ordering to the cascade. In the initial-state shower, values are gradually increasing as the hard scattering is approached, while is decreasing in the final-state showers. Shower evolution is cut off at some lower scale , typically around 1 GeV for QCD branchings. The same cut-off scale is also used to regularize the soft gluon emission divergences in the splitting kernels. From above, a maximum scale is introduced, where the showers are matched to the hard interaction itself. The relation between and the kinematics of the hard scattering is uncertain, and the choice made can strongly affect the amount of well-separated jets.

Despite a number of common traits, the initial- and final-state radiation machineries are in fact quite different, and are described separately below. For historical reasons, the final-state shower is found in Jetset and the initial-state one in Pythia.

Final-state showers are time-like, i.e. partons have . The evolution variable of the cascade is therefore in Jetset associated with the of the branching parton, but this choice is not unique. Starting from , an original parton is evolved downwards in until a branching occurs. The selected value defines the mass of the branching parton, and the of the splitting kernel the parton energy division between its daughters. These daughters may now, in turn, evolve downwards, in this case with maximum virtuality already defined by kinematics, and so on down to the cut-off.

In QCD showers, corrections to the leading-log picture, so-called coherence effects, lead to an ordering of subsequent emissions in terms of decreasing angles. This does not follow automatically from the mass-ordering constraint, but is implemented as an additional requirement on allowed emissions. Photon emission is not affected by angular ordering. It is also possible to obtain non-trivial correlations between azimuthal angles in the various branchings, some of which are implemented as options. Finally, the theoretical analysis strongly suggests the scale choice , and this is the default in the program.

The final-state radiation machinery is applied in the c.m. frame of the hard scattering. The total energy and momentum of the hard-scattering subsystem is preserved, as is the direction of the outgoing partons (in that frame).

In contrast to final-state showers, initial-state ones are space-like. This means that, in the sequence of branchings that lead up from the shower initiator to the hard interaction, particles and have . The ‘side branch’ particle , which does not participate in the hard scattering, may be on the mass shell, or have a time-like virtuality. In the latter case a time-like shower will evolve off it, rather like the final-state radiation described above. To first approximation, the evolution of the space-like main branch is characterized by the evolution variable , which is required to be strictly increasing along the shower, i.e. . Corrections to this picture have been calculated, but are basically absent in Pythia.

Initial-state radiation is handled within the backwards evolution scheme. In this approach, the choice of the hard scattering is based on the use of evolved parton distributions, which means that the inclusive effects of initial-state radiation are already included. What remains is therefore to construct the exclusive showers. This is done starting from the two incoming partons at the hard interaction, tracing the showers ‘backwards in time’, back to the two shower initiators. In other words, given a parton , one tries to find the parton that branched into . The evolution in the Monte Carlo is therefore in terms of a sequence of decreasing space-like virtualities and increasing momentum fractions . Branchings on the two sides are interleaved in a common sequence of decreasing values.

In the above formalism, there is no real distinction between gluon and photon emission. Some of the details actually do differ, as will be explained in the full description.

The initial- and final-state radiation shifts around the kinematics of the original hard interaction. In deep inelastic scattering, this means that the and values that can be derived from the momentum of the scattered lepton do not agree with the values originally picked. In high- processes, it means that one no longer has two jets with opposite and compensating , but more complicated topologies. Effects of any original kinematics selection cuts are therefore smeared out, an unfortunate side-effect of the parton-shower approach.

2.3 Beam Remnants

In a hadron–hadron collision, the initial-state radiation algorithm reconstructs one shower initiator in each beam. This initiator only takes some fraction of the total beam energy, leaving behind a beam remnant which takes the rest. For a proton beam, a quark initiator would leave behind a diquark beam remnant, with an antitriplet colour charge. The remnant is therefore colour-connected to the hard interaction, and forms part of the same fragmenting system. It is further customary to assign a primordial transverse momentum to the shower initiator, to take into account the motion of quarks inside the original hadron, basically as required by the uncertainty principle. This primordial is selected according to some suitable distribution, and the recoil is assumed to be taken up by the beam remnant.

Often the remnant is more complicated, e.g. a initiator would leave behind a proton remnant system in a colour octet state, which can conveniently be subdivided into a colour triplet quark and a colour antitriplet diquark, each of which are colour-connected to the hard interaction. The energy sharing between these two remnant objects, and their relative transverse momentum, introduces additional degrees of freedom, which are not understood from first principles.

Naïvely, one would expect an event to have only one beam remnant, and an event none. This is not always correct, e.g. a interaction in an event would leave behind the and as beam remnants, and a interaction in resolved photoproduction in an event would leave behind one and one in each remnant. Corresponding complications occur for photoproduction in events.

There is another source of beam remnants. If parton distributions are used to resolve an electron inside an electron, some of the original energy is not used in the hard interaction, but is rather associated with initial-state photon radiation. The initial-state shower is in principle intended to trace this evolution and reconstruct the original electron before any radiation at all took place. However, because of cut-off procedures, some small amount may be left unaccounted. Alternatively the user may have chosen to switch off initial-state radiation altogether, but still preserved the resolved electron parton distributions. In either case the remaining energy is given to a single photon of vanishing transverse momentum, which is then considered in the same spirit as ‘true’ beam remnants.

So far we have assumed that each event only contains one hard interaction, i.e. that each incoming particle has only one parton which takes part in hard processes, and that all other constituents sail through unaffected. This is appropriate in or events, but not necessarily so in hadron–hadron collisions. Here each of the beam particles contains a multitude of partons, and so the probability for several interactions in one and the same event need not be negligible. In principle these additional interactions could arise because one single parton from one beam scatters against several different partons from the other beam, or because several partons from each beam take place in separate scatterings. Both are expected, but combinatorics should favour the latter, which is the mechanism considered in Pythia.

The dominant QCD cross sections are divergent for , and drop rapidly for larger . Probably the lowest-order perturbative cross sections will be regularized at small by colour coherence effects: an exchanged gluon of small has a large transverse wave function and can therefore not resolve the individual colour charges of the two incoming hadrons; it will only couple to an average colour charge that vanishes in the limit . In the program, some effective scale is therefore introduced, below which the perturbative cross section is either assumed completely vanishing or at least strongly damped. Phenomenologically, comes out to be a number of the order of 1.5–2.0 GeV.

In a typical ‘minimum-bias’ event one therefore expects to find one or a few scatterings at scales around or a bit above , while a high- event also may have additional scatterings at the scale. The probability to have several high- scatterings in the same event is small, since the cross section drops so rapidly with .

The understanding of multiple interaction is still very primitive, and even the experimental evidence that it exists at all is rather weak. Pythia therefore contains several different options, with a fairly simple one as default. The options differ in particular on the issue of the ‘pedestal’ effect: is there an increased probability or not for additional interactions in an event which is known to contain a hard scattering, compared with one that contains no hard interactions?

2.4 Fragmentation

QCD perturbation theory, formulated in terms of quarks and gluons, is valid at short distances. At long distances, QCD becomes strongly interacting and perturbation theory breaks down. In this confinement regime, the coloured partons are transformed into colourless hadrons, a process called either hadronization or fragmentation. In this paper we reserve the former term for the combination of fragmentation and the subsequent decay of unstable particles.

The fragmentation process has yet to be understood from first principles, starting from the QCD Lagrangian. This has left the way clear for the development of a number of different phenomenological models. Three main schools are usually distinguished, string fragmentation (SF), independent fragmentation (IF) and cluster fragmentation (CF), but many variants and hybrids exist. Being models, none of them can lay claims to being ‘correct’, although some may be better founded than others. The best that can be aimed for is internal consistency, a good representation of existing data, and a predictive power for properties not yet studied or results at higher energies.

Jetset is intimately connected with string fragmentation, in the form of the time-honoured ‘Lund model’. This is the default for all Jetset/Pythia applications, but independent fragmentation options also exist, for applications where one wishes to study the importance of string effects.

All current models are of a probabilistic and iterative nature. This means that the fragmentation process as a whole is described in terms of one or a few simple underlying branchings, of the type jet hadron + remainder-jet, string hadron + remainder-string, and so on. At each branching, probabilistic rules are given for the production of new flavours, and for the sharing of energy and momentum between the products.

To understand fragmentation models, it is useful to start with the simplest possible system, a colour-singlet 2-jet event, as produced in annihilation. Here lattice QCD studies lend support to a linear confinement picture (in the absence of dynamical quarks), i.e. the energy stored in the colour dipole field between a charge and an anticharge increases linearly with the separation between the charges, if the short-distance Coulomb term is neglected. This is quite different from the behaviour in QED, and is related to the presence of a triple-gluon vertex in QCD. The details are not yet well understood, however.

The assumption of linear confinement provides the starting point for the string model. As the and partons move apart from their common production vertex, the physical picture is that of a colour flux tube (or maybe colour vortex line) being stretched between the and the . The transverse dimensions of the tube are of typical hadronic sizes, roughly 1 fm. If the tube is assumed to be uniform along its length, this automatically leads to a confinement picture with a linearly rising potential. In order to obtain a Lorentz covariant and causal description of the energy flow due to this linear confinement, the most straightforward way is to use the dynamics of the massless relativistic string with no transverse degrees of freedom. The mathematical, one-dimensional string can be thought of as parametrizing the position of the axis of a cylindrically symmetric flux tube. From hadron spectroscopy, the string constant, i.e. the amount of energy per unit length, is deduced to be GeV/fm. The expression ‘massless’ relativistic string is somewhat of a misnomer: effectively corresponds to a ‘mass density’ along the string.

Let us now turn to the fragmentation process. As the and move apart, the potential energy stored in the string increases, and the string may break by the production of a new pair, so that the system splits into two colour-singlet systems and . If the invariant mass of either of these string pieces is large enough, further breaks may occur. In the Lund string model, the string break-up process is assumed to proceed until only on-mass-shell hadrons remain, each hadron corresponding to a small piece of string with a quark in one end and an antiquark in the other.

In order to generate the quark–antiquark pairs which lead to string break-ups, the Lund model invokes the idea of quantum mechanical tunnelling. This leads to a flavour-independent Gaussian spectrum for the of pairs. Since the string is assumed to have no transverse excitations, this is locally compensated between the quark and the antiquark of the pair. The total of a hadron is made up out of the contributions from the quark and antiquark that together form the hadron. Some contribution of very soft perturbative gluon emission may also effectively be included in this description.

The tunnelling picture also implies a suppression of heavy-quark production, . Charm and heavier quarks hence are not expected to be produced in the soft fragmentation, but only in perturbative parton-shower branchings .

When the quark and antiquark from two adjacent string breakings are combined to form a meson, it is necessary to invoke an algorithm to choose between the different allowed possibilities, notably between pseudoscalar and vector mesons. Here the string model is not particularly predictive. Qualitatively one expects a ratio, from counting the number of spin states, multiplied by some wave-function normalization factor, which should disfavour heavier states.

A tunnelling mechanism can also be used to explain the production of baryons. This is still a poorly understood area. In the simplest possible approach, a diquark in a colour antitriplet state is just treated like an ordinary antiquark, such that a string can break either by quark–antiquark or antidiquark–diquark pair production. A more complex scenario is the ‘popcorn’ one, where diquarks as such do not exist, but rather quark–antiquark pairs are produced one after the other. This latter picture gives a less strong correlation in flavour and momentum space between the baryon and the antibaryon of a pair.

In general, the different string breaks are causally disconnected. This means that it is possible to describe the breaks in any convenient order, e.g. from the quark end inwards. One therefore is led to write down an iterative scheme for the fragmentation, as follows. Assume an initial quark moving out along the axis, with the antiquark going out in the opposite direction. By the production of a pair, a meson is produced, leaving behind an unpaired quark . A second pair may now be produced, to give a new meson , etc. At each step the produced hadron takes some fraction of the available energy and momentum. This process may be iterated until all energy is used up, with some modifications close to the end of the string in order to make total energy and momentum come out right.

The choice of starting the fragmentation from the quark end is arbitrary, however. A fragmentation process described in terms of starting at the end of the system and fragmenting towards the end should be equivalent. This ‘left–right’ symmetry constrains the allowed shape of the fragmentation function , where is the fraction of the remaining light-cone momentum (+ for the jet, for the one) taken by each new particle. The resulting ‘Lund symmetric fragmentation function’ has two free parameters, which are determined from data.

If several partons are moving apart from a common origin, the details of the string drawing become more complicated. For a event, a string is stretched from the end via the to the end, i.e. the gluon is a kink on the string, carrying energy and momentum. As a consequence, the gluon has two string pieces attached, and the ratio of gluon to quark string force is 2, a number which can be compared with the ratio of colour charge Casimir operators, . In this, as in other respects, the string model can be viewed as a variant of QCD where the number of colours is not 3 but infinite. Note that the factor 2 above does not depend on the kinematical configuration: a smaller opening angle between two partons corresponds to a smaller string length drawn out per unit time, but also to an increased transverse velocity of the string piece, which gives an exactly compensating boost factor in the energy density per unit string length.

The string will fragment along its length. To first approximation this means that there is one fragmenting string piece between and and a second one between and . One hadron is straddling both string pieces, i.e. sitting around the gluon corner. The rest of the particles are produced as in two simple strings, but strings boosted with respect to the overall c.m. frame. When considered in detail, the string motion and fragmentation is more complicated, with the appearance of additional string regions during the time evolution of the system. These corrections are especially important for soft and collinear gluons, since they provide a smooth transition between events where such radiation took place and events where it did not. Therefore the string fragmentation scheme is ‘infrared safe’ with respect to soft or collinear gluon emission.

For events that involve many partons, there may be several possible topologies for their ordering along the string. An example would be a (the gluon indices are here used to label two different gluon-momentum vectors), where the string can connect the partons in either of the sequences and . The matrix elements that are calculable in perturbation theory contain interference terms between these two possibilities, which means that the colour flow is not always well-defined. Fortunately, the interference terms are down in magnitude by a factor , where is the number of colours, so approximate recipes can be found. In the leading log shower description, on the other hand, the rules for the colour flow are well-defined. A final comment: in the argumentation for the importance of colour flows there is a tacit assumption that soft-gluon exchanges between partons will not normally mess up the original colour assignment; this is likely the case but has not been proven.

2.5 Decays

A large fraction of the particles produced by fragmentation are unstable and subsequently decay into the observable stable (or almost stable) ones. It is therefore important to include all particles with their proper mass distributions and decay properties. Although involving little deep physics, this is less trivial than it may sound: while a lot of experimental information is available, there is also very much that is missing. For charm mesons, it is necessary to put together measured exclusive branching ratios with some inclusive multiplicity distributions to obtain a consistent and reasonably complete set of decay channels, a rather delicate task. For bottom, so far only a rather simple phase-space type of generator has been used for hadronic decays.

Normally it is assumed that decay products are distributed according to phase space, i.e. that there is no dynamics involved in their relative distribution. However, in many cases additional assumptions are necessary, e.g. for semileptonic decays of charm and bottom hadrons one needs to include the proper weak matrix elements. Particles may also be produced polarized and impart a non-isotropic distribution to their decay products. Many of these effects are not at all treated in the program. In fact, spin information is not at all carried along, but has to be reconstructed explicitly when needed.

The normal decay treatment is handled by Jetset, making use of a set of tables where branching ratios and decay modes are stored. In Pythia a separate decay treatment exists, used exclusively for a specific list of particles: , , , , , , , , , , , , and the leptoquark . Together we call these resonances, and contrast the ‘particle decay’ treatment of Jetset with the ‘resonance decay’ one of Pythia. Of course, this is just a matter of terminology: a particle like the could also be called a resonance. What characterizes a (Pythia) resonance is that partial widths and branching ratios are calculated dynamically, as a function of the actual mass of a particle. Therefore not only do branching ratios change between an of nominal mass 100 GeV and one of 200 GeV, but also for a Higgs of nominal mass 200 GeV, the branching ratios would change between an actual mass of 190 GeV and 210 GeV, say. This is particularly relevant for reasonably broad resonances, and in threshold regions. For an approach like this to work, it is clearly necessary to have perturbative expressions available for all partial widths, which is one reason why a corresponding treatment would not be the same for an ordinary hadronic resonance, like the .

The decay products of Pythia resonances are typically quarks, leptons, or other resonances, e.g. or . In decays to quarks, parton showers are automatically added to give a more realistic multijet structure, and one may also allow photon emission off leptons. If the decay products in turn are resonances, further decays are necessary. Often spin information is available in resonance decay matrix elements, contrary to the normal state of affairs in ordinary particle decays. This means that the angular orientations in the two decays of a pair are properly correlated. Occasionally, the information is not available, and then resonances decay isotropically.

The top quark is a special problem. The original machinery is based on the assumption that the is long-lived, so that top hadrons have time to form in the fragmentation process, and afterwards these mesons decay weakly. With current ‘best bet’ mass values, this is not correct, but one should rather consider top decay before fragmentation. Top should then be handled like one of the above resonances. Therefore the program now contains an alternative along these lines, which is the preferred option.

3 Program Overview

This section contains a diverse collection of information. The first part is an overview of previous Jetset and Pythia versions. The second gives instructions for installation of the programs and describes their philosophy: how they are constructed and how they are supposed to be used. It also contains some information on how to read this manual. The third and final part contains several examples of pieces of code or short programs, to illustrate the general style of program usage. The last part is mainly intended as an introduction for completely new users, and can be skipped by more experienced ones.

Since the Jetset and Pythia programs today are so closely connected, and are gradually coalescing, they are presented together in this report. However, they still appear as separate entities, with slightly different style and emphasis.

Jetset is the older of the two, and is at the origin of the whole ‘Lund’ family of event generators. It can be subdivided in two parts. The larger is a generic package for jet fragmentation, particle decays, final-state parton showers, event-analysis routines, and other utitilies. This package can be used in the context of any hard process, provided one is willing to buy the underlying assumption of jet universality, i.e. that the fragmentation process is fundamentally the same whether one is considering an or a event, and that the only differences are to be found in the parton-level processes involved. This package is not only used by all other ‘Lund’ programs, but also by numerous other programs written to study specific processes. The smaller part of Jetset is a generator for annihilation events, according to either a parton-shower or a matrix-element approach. The Jetset program is completely selfcontained.

Pythia is a program made to generate hard or soft processes in collisions between leptons, hadrons and photons, especially at , and colliders. Where Jetset is a loose collection of routines that you can combine as desired, Pythia is a more structured program, where you initially set up what processes you want to study, and thereafter all events will be generated according to this specification. Included is an extensive library of hard subprocess differential cross sections, a library of parton distributions, a process generation machinery, treatment of initial-state showers and beam remnants, and a few odds and ends. Jetset is used for final-state showers, fragmentation and decay, but no other external libraries are needed. An interface to external parton-distribution function libraries is provided, however.

Many programs written by other persons make use of Jetset, and a few also of Pythia. It is not my intention to give a complete list here. A majority of these programs are specific to given collaborations, and therefore not publicly distributed. Below we give a list of a few public programs from the ‘Lund group’, which may have a somewhat wider application. None of them are supported by the current author, so any requests should be directed to the persons mentioned.

  • Ariadne is a generator for dipole emission, written mainly by L. Lönnblad [Pet88]. The dipole provides an alternative formulation of initial- and final-state showers. Jetset or Pythia can be used to generate the hard process and Jetset to do the fragmentation.

  • Aroma is a generator for heavy-flavour processes in leptoproduction, written by G. Ingelman and G. Schuler [Ing88]. It uses Jetset for fragmentation.

  • Fritiof is a generator for hadron–hadron, hadron–nucleus and nucleus–nucleus collisions [Nil87], which makes use of Pythia to generate hard QCD scatterings and of Jetset for fragmentation. Currently H. Pi is responsible for program development.

  • Lepto is a leptoproduction event generator, written mainly by G. Ingelman [Ing80]. It can generate parton configurations in deep inelastic scattering according to a number of possibilities. It makes use of Jetset for fragmentation and additionally has a parton-shower option based on Pythia.

  • Lucifer is a photoproduction generator written by G. Ingelman and A. Weigend [Ing87a]. It is a modification of an earlier version of Pythia and makes use of Jetset.

  • Pompyt is a generator for pomeron interactions written by P. Bruni and G. Ingelman [Bru93]. This program defines parton distributions, flux factors and other aspects specific to the pomeron, which is combined with the standard Pythia machinery for process generation.

  • Twister is a generator for higher-twist processes, written by G. Ingelman [Ing87]. It is a modification of an earlier version of Pythia and makes use of Jetset.

One should also note that a version of Pythia has been modified to include the effects of longitudinally polarized incoming protons. This is the work of St. Güllenstern et al. [Gül93].

3.1 Update History

Both Jetset and Pythia are by now fairly old and well-established programs, but they are still steadily being improved on. While evolution was especially rapid for Jetset in the early days, that program has by now reached a certain level of maturity, and the pace of change has dropped significantly. Pythia, on the other hand, has been continually extended in recent years, and may still see further growth, although most of the basic structure should be in place by now.

In earlier days, before the advent of electronic mail, programs were only infrequently distributed, and version numbers corresponded to distinct new upgrades. Today, the evolutionary process is more continuous and so is the distribution of new versions. In particular, the introduction of a new process or feature is often done on short notice, if no problems of backwards compatibility are involved. With this distribution, the subversion numbers have therefore been expanded to three digits, where the last two give sub-subversions. For every change made in the public file, the sub-subversion number is updated, together with the ‘last date of change’. In most referencing the shorter ‘Jetset version 7.4’ could still be preferable to e.g. ‘Jetset version 7.412’.

For the record, in Tables 1 and 2 we list the official main versions of Jetset and Pythia, respectively, with some brief comments.


Table 1: The main versions of Jetset, with their date of appearance, published manuals, and main changes from previous versions.

No. Date Publ. Main new or improved features
1 Nov 78 [Sjö78] single-quark jets
2 May 79 [Sjö79] heavy-flavour jets
3.1 Aug 79 2-jets in , preliminary 3-jets
3.2 Apr 80 [Sjö80] 3-jets in with full matrix elements,
toponium decays
3.3 Aug 80 softer fragmentation spectrum
4.1 Apr 81 baryon production and diquark fragmentation,
fourth-generation quarks, larger jet systems
4.2 Nov 81 low- physics
4.3 Mar 82 [Sjö82] 4-jets and QFD structure in ,
Jul 82 [Sjö83] event-analysis routines
5.1 Apr 83 improved string fragmentation scheme, symmetric
fragmentation, full 2 order QCD for
5.2 Nov 83 momentum-conservation schemes for IF,
initial-state photon radiation in
5.3 May 84 ‘popcorn’ model for baryon production
6.1 Jan 85 common blocks restructured, parton showers
6.2 Oct 85 [Sjö86] error detection
6.3 Oct 86 [Sjö87] new parton-shower scheme
7.1 Feb 89 new particle codes and common block structure,
more mesons, improved decays, vertex information,
Abelian gluon model, Bose–Einstein effects
7.2 Nov 89 interface to new standard common block,
photon emission in showers
7.3 May 90 [Sjö92d] expanded support for non-standard particles
7.4 Dec 93 [Sjö94] updated particle data and defaults


Table 2: The main versions of Pythia, with their date of appearance, published manuals, and main changes from previous versions.

No. Date Publ. Main new or improved features
1 Dec 82 [Ben84] synthesis of predecessors Compton, Highpt and
3.3 Feb 84 [Ben84a] scale-breaking parton distributions
3.4 Sep 84 [Ben85] more efficient kinematics selection
4.1 Dec 84 initial- and final-state parton showers, and
4.2 Jun 85 multiple interactions
4.3 Aug 85 , , and processes
4.4 Nov 85 , , processes
4.5 Jan 86 production, diffractive and elastic events
4.6 May 86 angular correlation in resonance pair decays
4.7 May 86 and processes
4.8 Jan 87 [Ben87] variable impact parameter in multiple interactions
4.9 May 87 process
5.1 May 87 massive matrix elements for heavy quarks
5.2 Jun 87 intermediate boson scattering
5.3 Oct 89 new particle and subprocess codes, new common block
structure, new kinematics selection, some
lepton–lepton and lepton–hadron interactions,
new subprocesses
5.4 Jun 90 -dependent widths, resonances not on the mass shell,
new processes, new parton distributions
5.5 Jan 91 improved and , several new processes
5.6 Sep 91 [Sjö92d] reorganized parton distributions, new processes,
user-defined external processes
5.7 Dec 93 [Sjö94] new total cross sections, photoproduction, top decay

All versions preceding Jetset 7.3 and Pythia 5.6 should now be considered obsolete, and are no longer supported. For stable applications, the earlier combination Jetset 6.3 and Pythia 4.8 could still be used, however.

Jetset version 7 and Pythia version 5 have been evolved in parallel, so some of the processes added in later versions of Pythia make use of particle data only found in Jetset from that time onwards. Although it would be possible to combine Pythia 5.7 with Jetset 7.3, e.g., it is not recommended. From the current versions onwards, checks have therefore been introduced to detect the use of (potentially) incompatible subversions, with warnings issued at initialization if that should be the case.

Previous versions of the manuals have contained detailed lists of modifications from one version to the next, see e.g. [Sjö92d]. Below we only reproduce the updates that appear with the most recent versions of the programs. Some of them were introduced in later editions of Pythia 5.6 with Jetset 7.3, while others are completely new. If nothing is explicitly said, these changes do not affect backwards compatibility, but only add new features.

3.1.1 Updates in JETSET 7.4

Changes from version 7.3 to 7.4 are not so large, although the impact of the updated particle data and parameter default values may need to be studied.

  • Particle data have been updated in accordance with the 1992 Review of Particle Properites [PDG92]. (As usual, with a free interpretation of inconsistencies, unclarities and other gaps in the knowledge.) Changes are especially drastic for charm and bottom. In the bottom sector the decay properties are now given individually for , , , and , i.e. the generic data for ‘pseudoparticle’ 85 are only used for other weakly decaying baryons.

  • Also a few other Standard Model parameters have been updated, such as the and masses and widths, and the CKM matrix elements.

  • Fragmentation and parton shower parameters have been modified to reflect current LEP knowledge [LEP90], i.e. a minor retuning starting from an average of the ‘best’ parameter values obtained by the four LEP collaborations. Bose-Einstein effects are still left out. Flavour composition is unchanged, except for a suppression of production. Affected by the change are MSTJ(11), PARJ(21), PARJ(23), PARJ(24), PARJ(41), PARJ(42), PARJ(54)PARJ(58) and PARJ(81).

  • Several other default values have been changed for switches and parameters in the , parton shower and fragmentation parts of the programs. These changes are intended to reflect our current best knowledge. See MSTJ(26), MSTJ(27), MSTJ(41), MSTJ(46), MSTJ(50), MSTJ(110), PARJ(26), and PARJ(121)PARJ(125).

  • A common title page for Jetset and Pythia has been introduced with the LULOGO routine. Sub-subversion numbers are also given.

  • Several options have been added for the LUSHOW shower routine. See MSTJ(41), MSTJ(47) and MSTJ(50).

  • A quark produced in the decay of a top hadron is allowed to radiate according to the standard parton shower scheme.

  • The scalar gluon option contains the full electroweak angular distribution of 3-jet events.

  • The LUCOMP routine has been modified. Among other things, the now appears as a separate compressed code, further codes for diffractive states have been added to the current list, and the pomeron (reggeon, ) has been added as particle 29 (28, 38).

  • A minimum threshold for calorimeter cell energy has been introduced for the LUCELL routine.

  • All obsolescent features of the Fortran 90 standard have been removed, i.e. the program should work well either with a Fortran 77 compiler or with a Fortran 90 one.

  • A few minor errors have been corrected.

The following changes have been made since the beginning of 1994, i.e. since the original distribution 7.400:

  • Jetset version 7.401, 11 February 1994:

    • Protect against overflow in LUZDIS (needed on some machines).

  • Jetset version 7.402, 7 April 1994:

    • New option to suppress either hard or soft radiation in LUSHOW, see MSTJ(40).

    • A generic interface to an external decay library has been introduced, see MSTJ(28) and SUBROUTINE LUTAUD.

    • In a few places, a dot has been moved from the end of one line to the beginning of the next continuation line, or the other way around, to keep together tokens such as .EQ. or .AND., since some debuggers may otherwise complain.

    • A source of (harmless) division by zero in LUSHOW has been removed.

  • Jetset version 7.403, 15 July 1994:

    • Leptons and photons which are unrelated to the system feeling the Bose-Einstein effects do not have their energies and momenta changed in the global rescaling step of LUBOEI. (Example: events, where one decays leptonically; before these lepton momenta could be slightly changed, but now not.) Further, the LUBOEI routine has been changed to avoid an unintentional gap in the limits of the very first bin.

    • The option LUEDIT(16) (used e.g. from PYEVNT) has been improved with a more extensive search for missing daughter pointers.

    • The KLU(I,16) procedure for finding rank has been rewritten to work in the current Jetset version, which it did not before. However, note that it will only work for MSTU(16)=2. As a general comment, the options 14–17 of KLU were written at a time when possible event histories were less complex, and can not be guaranteed always to work today.

  • Jetset version 7.404, 25 August 1994:

    • LUSHOW has been corrected, so that if , or quarks (or or quarks masked as or ones) are given with masses that vary from event to event (a Breit-Wigner shape, e.g.), the current mass rather than the nominal mass is used to define the cut-off scales of parton shower evolution.

    • LULOGO has been modified to take into account that a new Pythia/Jetset description has been published in
      T. Sjöstrand, Computer Phys. Commun. 82 (1994) 74
      and is from now on the standard reference to these two programs.

  • Jetset version 7.405, 27 January 1995:

    • LUCELL has been corrected, in that in the option with smearing of energy rather than transverse energy, the conversion factor between the two was applied in the wrong direction.

    • LUSHOW has been corrected in one place where the PMTH array was addressed with the wrong order of the indices. This affected quark mass corrections in the matching to the three-jet matrix elements.

    • An additional check has been included in LUBOEI that there are at least two particles involved in the Bose-Einstein effects. (No problem except in some bizarre situations.)

  • Jetset version 7.406, 20 February 1995:

    • A new option has been added for the behaviour of the running in ULALEM. This is not added as a true physics scenario, but only to produce results with a given, fixed value for the hard events, while still keeping the conventional value in the limit. See MSTU(101), PARU(103), PARU(104). Additionally, the constant has been added to the parameter list, see PARU(105).

    • The LULOGO routine has been updated to reflect my change of affiliation.

  • Jetset version 7.407, 21 June 1995:

    • Header and LULOGO have been updated with respect to phone number and WWW access.

    • The PHEP and VHEP variables in the /HEPEVT/ common block are now assumed to be in DOUBLE PRECISION, in accord with the proposed LEP 2 workshop addendum to the standard.

    • In LUTEST a missing decimal point on the energy check has been reinstated.

    • In LUINDF an expression has been protected against vanishing denominator.

  • Jetset version 7.408, 23 August 1995:

    • Check against division by zero in LUSHOW.

3.1.2 Updates in PYTHIA 5.7

The updates from version 5.6 to 5.7 are all minor, and just about any program that ran with version 5.6 will also work with Pythia 5.7. However, as for Jetset, it should be noted that some important default values have been changed.

  • New parametrizations of the total cross sections of hadronic reactions, based in Donnachie–Landshoff [Don92], which replace the old ones.

  • New parametrizations of elastic and single and double diffractive cross sections of hadronic reactions, based on Schuler–Sjöstrand [Sch94, Sch93a], which replace the old ones. Also the slope parameters, the diffractive mass distributions and other aspects of the event generation have been changed accordingly.

  • A possibility to give own total, elastic and diffractive cross sections.

  • The single diffractive cross section has been split into its two constituents, and . As a consequence, the diffractive subprocess codes 92–94 have received changed meaning.

  • A new common block PYINT7 has been added for the expanded total cross section information, and this information has been partly removed from other common blocks.

  • A much extended description of photoproduction physics, with the possibility to simulate separately VMD, anomalous and direct processes [Sch93, Sch93a].

  • The selection of proton parton distributions that come with the program has been updated with the CTEQ2 ones, while some others have been removed. New default is the leading-order fit CTEQ2L.

  • Since the Pdflib library now has been expanded to contain also parton-distribution functions for the photon, the interfaces to the Pakpdf and Phopdf libraries have been removed. In addition, the interface to Pdflib has been modified, and is now for appropriate for Pdflib version 4.

  • An extension of hadron parton distributions into the low- and low- region [Sch93a].

  • The top quark can be made to decay before it has time to fragment. In view of the current best estimate for the top mass, this is the expected behaviour, and is therefore now default. Further, a parton shower is allowed to evolve in the top decay. Also fourth generation quarks are allowed to decay before they fragment, and so on.

  • It is possible to call PYEVNT with energies that vary from one event to the next, without the need to reinitialize.

  • Improved scheme for post-factor conservation of and in deep inelastic scattering.

  • Processes 15, 19, 30 and 35 have been expanded to cover production in addition to the one, with full interference.

  • New process 80, .

  • New process 110, .

  • New process 149, .

  • New option for initial state radiation to restrict angular range of emission in accordance with coherence considerations.

  • Some options have been added or removed, and default values have been changed. This includes KFIN (top parton distributions off by default), MSTP(7), MSTP(11), MSTP(14), MSTP(23), MSTP(30) (removed), MSTP(31),  MSTP(34), MSTP(45), MSTP(48), MSTP(49), MSTP(62), MSTP(67), MSTP(101), PARP(13), PARP(81), PARP(82), PARP(47) and PARP(101).


  • A few minor errors have been corrected.

The following changes have been made since the beginning of 1994, i.e. since the original distribution 5.700:

  • Pythia version 5.701, 27 January 1994:

    • The machinery to handle interactions is expanded to the level already available for . This in particular means that a number of new options appear for MSTP(14). Affected are also MINT(105), MINT(107), MINT(108), MINT(109), VINT(282) (removed), VINT(283) and VINT(284). Parametrizations are introduced for meson–meson total, elastic and diffractive cross sections, needed for the VMD part of the photon. The treatment of cross sections for hard processes, of initial state radiation, of beam remnants and of other aspects are also expanded to cover the new possibilities. A first study of the relevant physics aspects is found in [Sch94a].

    • An option is introduced to modify the scale of the anomalous part of the photon parton distributions, see MSTP(59) and PARP(59).

    • Correction of an error, where the generation of jet and low- events could give incorrect cross section information with PYSTAT(1) at low energies. The event generation itself was correct. (The error was introduced as a consequence of allowing variable energies.)

    • A rejection is introduced for top events where the top mass (selected according to a Breit-Wigner) is too low to allow the decay into a on the mass shell.

    • The correction of a few other minor bugs, probably harmless.

  • Pythia version 5.702, 13 February 1994:

    • The interface to Pdflib has been modified to reflect that TMAS should no longer be set except in first PDFSET call. (Else a huge amount of irrelevant warning messages are generated by Pdflib.)

    • The STOP statement in a few dummy routines has been modifed to avoid irrelevant compilation warning messages on IBM mainframes.

    • A few labels have been renumbered.

  • Pythia version 5.703, 22 February 1994:

    • Removal of a bug in PYRESD, which could give (under some specific conditions) errors in the colour flow.

  • Pythia version 5.704, 7 April 1994:

    • Process 11 has been corrected, for the part that concerns anomalous couplings (contact interactions) in the process. The error was present in the expression for and obvious permutations, while , and the others were correct.

    • The option with post-facto conservation in deep inelastic scattering can give infinite loops when applied to process 83, in particular if one asks for the production of a top. (Remember that the standard DIS kinematics is defined for massless quarks.) Therefore the switch MSTP(23) has been modifed so that by default only process 10 is affected.

    • PYRESD is modified to ensure isotropic angular distributions in the decays of the top or a fourth generation particle, i.e. in . This may not be the correct distribution but, unless explicit knowledge exists for a given process, this should always be the default.

    • In processes 16, 20, 31 and 36 the propagator has been modified to include -dependent widths in the Breit-Wigner shape. The most notable effect is a suppression of the low-mass tail of the mass spectrum.

    • When Pdflib is used, PDFSET is now only called whenever a different structure function is requested. For events therefore only one call is made, while interactions still involve a call to PDFSET for each STRUCTM one, since and structure functions have to be called alternatingly. MINT(93) is used to keep track of latest structure function called.


    • A number of minor errors have been corrected.

  • Pythia version 5.705, 15 July 1994:

    • A completely new possibility to have Pythia mix different allowed processes (direct, VMD and anomalous) in and interactions. This option can be accessed with MSTP(14)=10. The relevant physics description and programming details may be found in sections 7.7.2 and 8.3.2. This facility is still not definitive, in that it is hoped to gradually enhance it with further features. The cross-section output of the PYSTAT has been expanded to reflect the further subdivision of the total cross section.

    • The new facility above has required a major restructuring of some of the code: the routine PYEVKI has been removed, new routines PYINBM, PYINPR and PYSAVE created, and some material has been moved to or from PYINIT, PYINRE and PYINKI. New variables include MSTI(9), MINT(121), MINT(122), MINT(123) and VINT(285).

    • The GRV leading-order dynamically generated parton distributions for the and have been included as options, see MSTP(51) and MSTP(53).

    • A parametrization of the homogeneous solution to the anomalous photon parton distributions have been added as an option, see MSTP(56)=3.

    • The treatment of the anomalous photon component can be modified with the new switch MSTP(15) and variable PARP(17); at the same time MSTP(59) and PARP(59) have been removed. The new options are mainly intended for comparative studies and should not normally be touched.

    • The option MSTP(92)=5 for beam remnant treatment erroneously missed some statements which now have been inserted. Further, new options have been added for the beam remnant splitting of momentum between a hadron and a quark/diquark jet, where MSTP(94) should now be used rather than MSTP(92).

    • In PYDIFF the recoiling gluon energy is calculated in a numerically more stable fashion.

  • Pythia version 5.706, 25 August 1994:

    • New processes 167 and 168, and , respectively, have been introduced. These contact interaction production processes of excited quarks complement the quark–gluon fusion ones in processes 147 and 148, and obey the same general rules, see section 8.5.5.

    • The option MSTP(57)=3 now also allows a dampening of parton distributions.

    • A few minor errors have been corrected.

  • Pythia version 5.707, 20 October 1994:

    • A major bug discovered in processes 121 and 122 (and thus also affecting 181, 182, 186 and 187), : the kinematics was incorrectly handed on to the Kunszt matrix elements. This affected the default option , but effects were especially dramatic when the alternative was used. The choice of appropriate scale for structure functions introduces a further uncertainty in cross sections for the processes above. So long as only quarks are considered, the mass is a reasonable choice, but for the alternative this is presumably too low. Therefore new options have been introduced in MSTP(39), with the default behaviour changed.

    • Another important bug corrected in the calculation of the reduction of cross section when decay modes are forced. This occured when both and produced a , and and decay modes were set differently.

  • Pythia version 5.708, 25 October 1994:

    • A few further places changed to make processes 181, 182, 186 and 187 work (see version 5.707 above).

  • Pythia version 5.709, 26 October 1994:

    • The matrix element for has been replaced, using the formulae of
      D. Bardin, M. Bilenky, D. Lehner, A. Olchevski and T. Riemann, CERN-TH.7295/94,
      but with the dependence on the variable not integrated out (D. Bardin, private communication). This avoids some problems encountered in the old expressions when one or both ’s were far off the mass shell.

    • Change in calls to Pdflib, so that the input is always at least the of the respective set.

    • Extra protection against infinite loops in PYSSPA.

  • Pythia version 5.710, 27 January 1995:

    • The dimensions of the HGZ array in PYRESD has been expanded to avoid accidental writing outside the bounds.

    • VINT(41)-VINT(66) are saved and restored in PYSCAT, for use in low- events, when beam remnant treatment has failed (with nonzero MINT(57)).

    • The routine PYSTGH has been replaced by the routine PYSTHG. This contains an improved parametrization of the homogeneous evolution of an anomalous photon from some given initial scale. The argument NF of the PYSTGH routine has been removed; now is always automatically converted to the relevant -flavour value from its 4-flavour one, at flavour thresholds.

  • Pythia version 5.711, 20 February 1995:

    • New possibilities have been added to switch between electroweak couplings being expressed in terms of a running or in terms of a fixed Fermi constant . This affects both decay widths and process cross sections, in the routines PYINRE, PYRESD, PYWIDT and PYSIGH. See MSTP(8) for details; default corresponds to old behaviour.

    • The option MSTP(37)=1, with running quark masses in couplings to Higgs bosons, only works when is allowed to run (so one can define a value). Therefore a check has been introduced in PYWIDT and PYSIGH that the option MSTP(37)=1 is only executed if additionally MSTP(2).

    • Some non-physics changes have been made in the RKBBV and STRUCTM codes so as to avoid some (in principle harmless) compiler warnings.

  • Pythia version 5.712, 15 March 1995:

    • A serious error has been corrected in the MSTP(173)=1 option, i.e. when the program is run with user-defined weights that should compensate for a biased choice of variable beam energies. This both affected the relative admixture of low- and high- events and the total cross section obtained by Monte Carlo integration. (PYRAND changed.)

    • In order to improve the flexibility and efficiency of the variable-energy option, the user should now set PARP(174) before the PYINIT call, and thereafter not change it. This allows PARP(173) weights of arbitrary size. (PYRAND and PYMAXI changed.)

    • MSTI(5) (and MINT(5)) are now changed so they count the number of successfully generated events, rather than the number of tries made. This change only affects runs with variable energies, MSTP(171)=1 and MSTP(172)=2, where MSTI(61)=1 signals that a user-provided energy has been rejected in the weighting. This change also affects PARI(2), which becomes the cross section per fully generated event. (PYEVNT changed.)

    • The option MSTP(14)=10 has now been extended so that it also works for deep inelastic scattering of an electron off a (real) photon, i.e. subprocess 10. What is obtained is a mixture of the photon acting as a vector meson and it acting as an anomalous state. This should therefore be the sum of what can be obtained with MSTP(14)=2 and =3. It is distinct from MSTP(14)=1 in that different sets are used for the parton distributions — in MSTP(14)=1 all the contributions to the photon distributions are lumped together, while they are split in VMD and anomalous parts for MSTP(14)=10. Also the beam remnant treatment is different, with a simple Gaussian distribution (at least by default) for MSTP(14)=1 and the VMD part of MSTP(14)=10, but a powerlike distribution between PARP(15) and for the anomalous part of MSTP(14)=10. (PYINIT, PYINPR and PYSTAT changed.)
      To access this option for and as incoming beams, it is only necessary to set MSTP(14)=10 and keep MSEL at its default value. Unlike the corresponding option for and , no cuts are overwritten, i.e. it is still the responsability of the user to set these appropriately. Those especially appropriate for DIS usage are CKIN(21)-CKIN(22) or CKIN(23)-CKIN(24) for the range (former or latter depending on which side is the incoming real photon), and CKIN(35)-CKIN(36) for the range. A further new option has been added (in PYKLIM) to set the range as well, see CKIN(39)-CKIN(40).
      A warning about the usage of Pdflib for photons. So long as MSTP(14)=1, i.e. the photon is not split up, Pdflib is accessed by MSTP(56)=2 and MSTP(55) the parton distribution set, as described in the manual. However, when the VMD and anomalous pieces are split, the VMD part is based on a rescaling of pion distributions by VMD factors (except for the SaS sets, that already come with a separate VMD piece). Therefore, to access Pdflib for MSTP(14)=10, it is not correct to set MSTP(56)=2 and a photon distribution in MSTP(55). Instead, one should put MSTP(56)=2, MSTP(54)=2 and a pion distribution code in MSTP(53), while MSTP(55) has no function. The anomalous part is still based on the SaS parametrization, with PARP(15) as main free parameter.

    • A change has been made in PYREMN to reduce the possibility of infinite loops.

  • Pythia version 5.713, 22 March 1995:

    • The SaS parton distributions of the photons are now available, see [Sch95]. There are four new sets. These differ in that two use a  GeV and two a  GeV, and in that two use the DIS and two the conventions for the dominant non-leading contributions. (However, the fits are formally still leading-order, in that not all next-to-leading contributions have been included.) New default is the SaS 1D set. Furthermore, for the definition of , additional terms appear that do not form part of the parton distributions itself. To partly take this into account, an additional doubling of the possibilities has been included. These eight possibilites can be accesed with MSTP(55). The default value of PARP(15) has been changed from 0.5 to 0.6 GeV, for consistency with SaS 1D.
      The generic routine PYSTFU has been rewritten to handle the interfacing. The old routines PYSTAG, PYSTGS, PYDILN and PYSTHG have been removed. Instead the routines of the SaSgam library have been inserted. In order to avoid any clashes, the routines SAS*** have been renamed PYG***. Thus new routines are PYGGAM, PYGVMD, PYGANO, PYGBEH and PYGDIR. The common block SASCOM is renamed PYINT8. If you want to use the parton distributions for standalone purposes, you are encouraged to use the original SaSgam routines rather than going the way via the Pythia adaptations.

    • PYDOCU has been corrected so that PARI(2) refers to the full cross section for and processes, rather than that of the latest subprocess considered.

    • An additional check has been inserted into PYREMN.

  • Pythia version 5.714, 22 March 1995:

    • Some minor modifications to PYSTFU and PYGGAM in the wake of the changes of the previous version.

  • Pythia version 5.715, 24 April 1995:

    • An unfortunate choice of default values has been corrected: the old MSTP(3)=2 value implied that was entirely based on the value of the proton structure function; also e.g. for annihilation events. Thus the in PARJ(81) was overwritten, i.e. did not keep the value required by standard phenomenology, which typically gave too narrow jets. (While switching to MSTP(3)=1 it worked fine.) In the modified option MSTP(3)=2 this has been corrected, to better agree with user expectations. Since further changes were made in version 5.716, we refer below for additional comments.

    • The form for PTMANO, the for anomalous processes, as used in PYINPR when processes are mixed for or events, has been updated to match (as well as can be expected) the SaS 1D photon distributions.

  • Pythia version 5.716, 30 June 1995:

    • The strategy for the changes to in version 5.715 above have been modified for better transparency. Now PARJ(81) is used for resonance decays (including e.g. decay, from which it is determined), and PARP(72) for other timelike showers. PARJ(81) is not overwritten for MSTP(3)=2, but only for =3. Changes affect PYINIT, PYEVNT and PYRESD.

    • A new multiplicative factor has been introduced for the scale choice of the hard scattering in PYSIGH, affecting parton distributions and , see PARP(34).

    • PYREMN has been corrected for occasional too large boost factors.

    • An error in PYSIGH for process 148 has been corrected.

    • The MSTP(62)=1 option of PYSSPA is modified to avoid division by zero.

    • Header has been updated with WWW-information.

  • Pythia version 5.717, 23 August 1995:

    • MIN1, MIN2, MAX1, MAX2, MINA and MAXA in PYSIGH have had an extra M prefixed to avoid confusion with Fortran functions.

    • Protect against MDCY(0,1) being accessed in PYSIGH.

    • Protect against THB=0 in PYRAND.

    • Protect against YSTMAX-YSTMIN = 0 in PYSIGH.

    • Check for moved leptoquark at beginning of PYRESD just like for other particles with colour.

3.2 Program Installation

Several ‘authorized’ sources of the programs exist. The ‘master copy’ of the programs is the one found on my World Wide Web homepage
There you have: jetset74.f the Jetset code, pythia57.f the Pythia code, pythia57.tex this common Pythia/Jetset manual, and update57.notes plain text update notes to the manual. In addition to these, one may also find older versions of the program and manuals, sample main programs and other pieces of related software, and other physics papers.

The lack of stable versions may make it less convenient to rely on the above files. New versions are introduced in the general distribution of the CERN program library, maybe once a year. These versions are better checked before release, and should be useful for most applications. However, clearly, they may be less up-to-date. Read the CERN Computer Newsletter for announcements. Copies of the programs are also available via anonymous ftp, e.g. from the asisftp server at CERN.

The programs are written entirely in standard Fortran 77, and should run on any machine with such a compiler. To a first approximation, program compilation should therefore be straightforward.

Unfortunately, experience with many different compilers has been uniform: the options available for obtaining optimized code actually produce erroneous code (e.g. operations inside DO loops are moved out before them, where some of the variables have not yet been properly set). Therefore the general advice is to use a low optimization level. Note that this is often not the default setting.

SAVE statements have been included in accordance with the Fortran standard. Since most ordinary machines take SAVE for granted, this part is not particularly well tried out, however.

All default settings and particle and process data are stored in BLOCK DATA LUDATA for Jetset and BLOCK DATA PYDATA for Pythia. These subprograms must be linked for a proper functioning of the other routines. On some machines this is not done automatically but must be forced by you, in particular if Jetset and Pythia are maintained as libraries from which routines are to be loaded only when they are needed. In this connection we note that the library approach does not give any significant space advantages over a loading of the packages as a whole, since a normal run will call on most of the routines anyway, directly or indirectly.

Since most machines in current use are 32-bit ones, this is the precision normally assumed. A few pieces of code have therefore had to be written in double precision. As a rule of thumb, double-precision variables have as first character D, but there are a few exceptions.

For applications at very high energies, such as LHC, the use of single precision for any real variable is a problem. It might then be necessary to rewrite the program completely, i.e. to have a declaration IMPLICIT DOUBLE PRECISION(A-H,O-Z) at the beginning of each subprogram, and to change all real constants to double precision. Needless to say, the latter is a major undertaking. In some cases, shortcuts are available. On the IBM, for instance, the AUTODBL compiler option for automatic precision doubling works fine, provided only that an even number of integers precede real numbers in common blocks. In Jetset you therefore need to introduce an additional integer variable (NPAD, say) directly after N in the LUJETS common block, and in Pythia an additional integer (MSEPAD) after MSEL in the PYSUBS common block. Some pieces of code will then actually run in quadruple precision.

A test program, LUTEST, is included in the Jetset package. It is disguised as a subroutine, so you have to run a main program

      CALL LUTEST(1)

This program will generate six hundred events of different types, under a variety of conditions. If Jetset has not been properly installed, this program is likely to crash, or at least generate a number of erroneous events. This will then clearly be marked in the output, which otherwise will just contain a few sample event listings and a table of the number of different particles produced. To switch off the output of normal events and final table, use LUTEST(0) instead of LUTEST(1). The final tally of errors detected should read 0.

In exactly the same vein, a test program PYTEST comes with the Pythia package. You then have to run a program

      CALL PYTEST(1)

As before the alternative PYTEST(0) will give a less extensive listing. No errors should appear during execution.

3.3 Program Philosophy

The Monte Carlo programs are built as slave systems, i.e. you, the user, have to supply the main program. From this the various subroutines are called on to execute specific tasks, after which control is returned to the main program. Some of these tasks may be very trivial, whereas the ‘high-level’ routines by themselves may make a large number of subroutine calls. Many routines are not intended to be called directly by you, but only from higher-level routines such as LUEXEC, LUEEVT, PYINIT or PYEVNT.

Basically, this means that there are three ways by which you communicate with the programs. First, by setting common block variables, you specify the details of how the programs should perform specific tasks, i.e. which subprocesses should be generated (for Pythia), which particle masses should be assumed, which coupling constants used, which fragmentation scenarios, and so on with hundreds of options and parameters. Second, by calling subroutines you tell the programs to generate events according to the rules established above. Normally there are few subroutine arguments, and those are usually related to details of the physical situation, such as what c.m. energy to assume for events. Third, you can either look at the common block LUJETS to extract information on the generated event, or you can call on various functions and subroutines to analyse the event further for you.

It should be noted that, while the physics content is obviously at the centre of attention, the Jetset/Pythia package also contains a very extensive setup of auxiliary service routines. The hope is that this will provide a comfortable working environment, where not only events are generated, but where you also linger on to perform a lot of the subsequent studies. Of course, for detailed studies, it may be necessary to interface the output directly to a detector simulation program.

The general rule is that all routines have names that are six characters long, beginning with LU for Jetset routines and PY for Pythia ones. Real-valued functions in Jetset begin with UL instead. There are three exceptions to both the length and the initial character rules: KLU, PLU and RLU. The former two functions are strongly coupled to the K and P matrices in the LUJETS common block, the latter uses R to emphasize the rôle as a random-number generator. Also common block names are six characters long and start with LU or PY.

On the issue of initialization, Jetset and Pythia behave quite differently. Most Jetset routines work without any initialization (except for the one implied by the presence of BLOCK DATA LUDATA, see above), i.e. each event and each task stand on their own. Current common block values are used to perform the tasks in specific ways, and those rules can be changed from one event to the next (or even within the generation of one and the same event) without any penalty. The random-number generator is initialized at the first call, but usually this is transparent. Therefore the two Jetset routines LUEEVT (and some of the routines called by it) and LUONIA are basically the only ones to contain some elements of initialization, where there are a few advantages if events are generated in a coherent fashion, but even here the penalty for not doing it is small.

In Pythia, on the other hand, a sizeable amount of initialization is performed in the PYINIT call, and thereafter the events generated by PYEVNT all obey the rules established at that point. Therefore common block variables that specify methods to be used have to be set before the PYINIT call and then not be changed afterwards, with few exceptions. Of course, it is possible to perform several PYINIT calls in the same run, but there is a significant time overhead involved, so this is not something one would do for each new event.

Apart from writing a title page, giving a brief initialization information, printing error messages if need be, and responding to explicit requests for listings, all tasks of the programs are performed ‘silently’. All output is directed to unit MSTU(11), by default 6, and it is up to you to set this unit open for write. The only exceptions are RLUGET, RLUSET and LUUPDA where, for obvious reasons, the input/output file number is specified at each call. Here you again have to see to it that proper read/write access is set.

The programs are extremely versatile, but the price to be paid for this is having a large number of adjustable parameters and switches for alternative modes of operation. No single user is ever likely to need more than a fraction of the available options. Since all these parameters and switches are assigned sensible default values, there is no reason to worry about them until the need arises.

Unless explicitly stated (or obvious from the context) all switches and parameters can be changed independently of each other. One should note, however, that if only a few switches/parameters are changed, this may result in an artificially bad agreement with data. Many disagreements can often be cured by a subsequent retuning of some other parameters of the model, in particular those that were once determined by a comparison with data in the context of the default scenario. For example, for annihilation, such a retuning could involve one QCD parameter ( or ), the longitudinal fragmentation function, and the average transverse fragmentation momentum.

The programs contain a number of checks that requested processes have been implemented, that flavours specified for jet systems make sense, that the energy is sufficient to allow hadronization, that the memory space in LUJETS is large enough, etc. If anything goes wrong that the program can catch (obviously this may not always be possible), an error message will be printed and the treatment of the corresponding event will be cut short. In serious cases, the program will abort. As long as no error messages appear on the output, it may not be worthwhile to look into the rules for error checking, but if but one message appears, it should be enough cause for alarm to receive prompt attention. Also warnings are sometimes printed. These are less serious, and the experienced user might deliberately do operations which go against the rules, but still can be made to make sense in their context. Only the first few warnings will be printed, thereafter the program will be quiet. By default, the program is set to stop execution after ten errors, after printing the last erroneous event.

It must be emphasized that not all errors will be caught. In particular, one tricky question is what happens if an integer-valued common block switch or subroutine/function argument is used with a value that is not defined. In some subroutine calls, a prompt return will be expedited, but in most instances the subsequent action is entirely unpredictable, and often completely haywire. The same goes for real-valued variables that are assigned values outside the physically sensible range. One example will suffice here: if PARJ(2) is defined as the suppression factor, a value will not give more profuse production of than of , but actually a spillover into production. Users, beware!

3.4 Manual Conventions

In the manual parts of this report, some conventions are used. All names of subprograms, common blocks and variables are given in upper-case ‘typewriter’ style, e.g. MSTP(111)=0. Also program examples are given in this style.

If a common block variable must have a value set at the beginning of execution, then a default value is stored in one of the block data subprograms LUDATA and PYDATA. Such a default value is usually indicated by a ‘(D=…)’ immediately after the variable name, e.g.

  • (D=1) choice of fragmentation scheme.

All variables in the Jetset common blocks (with very few exceptions, clearly marked) can be freely changed from one event to the next, or even within the treatment of one single event. In the Pythia common blocks the situation is more complicated. The values of many switches and parameters are used already in the PYINIT call, and cannot be changed after that. The problem is mentioned in the preamble to the afflicted common blocks, which in particular means /PYPARS/ and /PYSUBS/. For the variables which may still be changed from one event to the next, a ‘(C)’ is added after the ‘(D=…)’ statement.

Normally, variables internal to the program are kept in separate common blocks and arrays, but in a few cases such internal variables appear among arrays of switches and parameters, mainly for historical reasons. These are denoted by ‘(R)’ for variables you may want to read, because they contain potentially interesting information, and by ‘(I)’ for purely internal variables. In neither case may the variables be changed by you.

In the description of a switch, the alternatives that this switch may take are often enumerated, e.g.

  • (D=1) choice of fragmentation scheme.

    • no jet fragmentation at all.

    • string fragmentation according to the Lund model.

    • independent fragmentation, according to specification in MSTJ(2) and MSTJ(3).

If you then use any value other than 0, 1 or 2, results are unpredictable. The action could even be different in different parts of the program, depending on the order in which the alternatives are identified.

It is also up to you to choose physically sensible values for parameters: there is no check on the allowed ranges of variables. We gave an example of this at the end of the preceding section.

Subroutines you are expected to use are enclosed in a box at the point where they are defined:


This is followed by a description of input or output parameters. The difference between input and output is not explicitly marked, but should be obvious from the context. In fact, the event-analysis routines of section 15.4 return values, while all the rest only have input variables.

Routines that are only used internally are not boxed in. However, we use boxes for all common blocks, so as to enhance the readability.

3.5 Getting Started with JETSET

As a first example, assume that you want to study the production of 2-jet systems at 20 GeV energy. To do this, write a main program

      CALL LU2ENT(0,2,-2,20.)
      CALL LULIST(1)

and run this program, linked together with Jetset. The routine LU2ENT is specifically intended for storing two entries (jets or particles). The first argument (0) is a command to perform fragmentation and decay directly after the entries have been stored, the second and third that the two entries are (2) and (), and the last that the c.m. energy of the pair is 20 GeV. When this is run, the resulting event is stored in the LUJETS common block. This information can then be read out by you. No output is produced by LU2ENT itself, except for a title page which appears once for every Jetset/Pythia run.

Instead the second command, to LULIST, provides a simple visible summary of the information stored in LUJETS. The argument (1) indicates that the short version should be used, which is suitable for viewing the listing directly on an 80-column terminal screen. It might look as shown here.

                        Event listing (summary)

  I  particle/jet KS     KF orig   p_x     p_y     p_z      E       m

  1  (u)       A  12      2    0   0.000   0.000  10.000  10.000   0.006
  2  (u~)      V  11     -2    0   0.000   0.000 -10.000  10.000   0.006
  3  (string)     11     92    1   0.000   0.000   0.000  20.000  20.000
  4  (rho+)       11    213    3   0.098  -0.154   2.710   2.856   0.885
  5  (rho-)       11   -213    3  -0.227   0.145   6.538   6.590   0.781
  6  pi+           1    211    3   0.125  -0.266   0.097   0.339   0.140
  7  (Sigma0)     11   3212    3  -0.254   0.034  -1.397   1.855   1.193
  8  (K*+)        11    323    3  -0.124   0.709  -2.753   2.968   0.846
  9  p~-           1  -2212    3   0.395  -0.614  -3.806   3.988   0.938
 10  pi-           1   -211    3  -0.013   0.146  -1.389   1.403   0.140
 11  pi+           1    211    4   0.109  -0.456   2.164   2.218   0.140
 12  (pi0)        11    111    4  -0.011   0.301   0.546   0.638   0.135
 13  pi-           1   -211    5   0.089   0.343   2.089   2.124   0.140
 14  (pi0)        11    111    5  -0.316  -0.197   4.449   4.467   0.135
 15  (Lambda0)    11   3122    7  -0.208   0.014  -1.403   1.804   1.116
 16  gamma         1     22    7  -0.046   0.020   0.006   0.050   0.000
 17  K+            1    321    8  -0.084   0.299  -2.139   2.217   0.494
 18  (pi0)        11    111    8  -0.040   0.410  -0.614   0.751   0.135
 19  gamma         1     22   12   0.059   0.146   0.224   0.274   0.000
 20  gamma         1     22   12  -0.070   0.155   0.322   0.364   0.000
 21  gamma         1     22   14  -0.322  -0.162   4.027   4.043   0.000
 22  gamma         1     22   14   0.006  -0.035   0.422   0.423   0.000
 23  p+            1   2212   15  -0.178   0.033  -1.343   1.649   0.938
 24  pi-           1   -211   15  -0.030  -0.018  -0.059   0.156   0.140
 25  gamma         1     22   18  -0.006   0.384  -0.585   0.699   0.000
 26  gamma         1     22   18  -0.034   0.026  -0.029   0.052   0.000
                 sum:  0.00        0.000   0.000   0.000  20.000  20.000

(A few blanks have been removed between the columns to make it fit into the format of this text.) Look in the particle/jet column and note that the first two lines are the original and , where ‘bar’ is actually written ‘’ to save space in longer names. The parentheses enclosing the names, ‘(u)’ and ‘(u~)’, are there as a reminder that these jets actually have been allowed to fragment. The jets are still retained so that event histories can be studied. Also note that the KF (flavour code) column contains 2 in the first line and in the second. These are the codes actually stored to denote the presence of a and a , cf. the LU2ENT call, while the names written are just conveniences used when producing visible output. The A and V near the end of the particle/jet column indicate the beginning and end of a string (or cluster, or independent fragmentation) parton system; any intermediate entries belonging to the same system would have had an I in that column. (This gives a poor man’s representation of an up-down arrow, .)

In the orig (origin) column, the zeros indicate that and are two initial entries. The subsequent line, number 3, denotes the fragmenting string system as a whole, and has origin 1, since the first parton of this string system is entry number 1. The particles in lines 4–10 have origin 3 to denote that they come directly from the fragmentation of this string. In string fragmentation it is not meaningful to say that a particle comes from only the quark or only the one. It is the string system as a whole that gives a , a , a , a , a , a , and a . Note that some of the particle names are again enclosed in parentheses, indicating that these particles are not present in the final state either, but have decayed further. Thus the in line 13 and the in line 14 have origin 5, as an indication that they come from the decay of the in line 5. Only the names not enclosed in parentheses remain at the end of the fragmentation/decay chain, and are thus experimentally observable. The actual status code used to distinguish between different classes of entries is given in the KS column; codes in the range 1–10 correspond to remaining entries, and those above 10 to those that have fragmented or decayed.

The columns with p_x, p_y, p_z, E and m are quite self-explanatory. All momenta, energies and masses are given in units of GeV, since the speed of light is taken to be . Note that energy and momentum are conserved at each step of the fragmentation/decay process (although there exist options where this is not true). Also note that the axis plays the rôle of preferred direction, along which the original partons are placed. The final line is intended as a quick check that nothing funny happened. It contains the summed charge, summed momentum, summed energy and invariant mass of the final entries at the end of the fragmentation/decay chain, and the values should agree with the input implied by the LU2ENT arguments. (In fact, warnings would normally appear on the output if anything untoward happened, but that is another story.)

The above example has illustrated roughly what information is to be had in the event record, but not so much about how it is stored. This is better seen by using a 132-column format for listing events. Try e.g. the following program

      CALL LU3ENT(0,1,21,-1,30.,0.9,0.7)
      CALL LULIST(2)
      CALL LUEDIT(3)
      CALL LULIST(2)

where a 3-jet event is generated in the first line and listed in the second. This listing will contain the numbers as directly stored in the common block LUJETS


For particle I, K(I,1) thus gives information on whether or not a jet or particle has fragmented or decayed, K(I,2) gives the particle code, K(I,3) its origin, K(I,4) and K(I,5) the position of fragmentation/decay products, and P(I,1)–P(I,5) momentum, energy and mass. The number of lines in current use is given by N, i.e. 1 I N. The V matrix contains decay vertices; to view those LULIST(3) has to be used. It is important to learn the rules for how information is stored in LUJETS.

The third line in the program illustrates another important point about Jetset: a number of routines are available for manipulating the event record after the event has been generated. Thus LUEDIT(3) will remove everything except stable charged particles, as shown by the result of the second LULIST call. More advanced possibilities include things like sphericity or clustering routines.

Apart from the input arguments of subroutine calls, control on the doings of Jetset may be imposed via the LUDAT1, LUDAT2, LUDAT3 and LUDAT4 common blocks. Here sensible default values are always provided. A user might want to switch off all particle decays by putting MSTJ(21)=0 or increase the ratio in fragmentation by putting PARJ(2)=0.40, to give but two examples. It is by exploring the possibilities offered here that Jetset can be turned into an extremely versatile tool, even if all the nice physics is already present in the default values.

As a final, semirealistic example, assume that the spectrum of particles is to be studied in 91.2 GeV annihilation events, where is to be defined with respect to the sphericity axis. Using the HBOOK package (version 4, watch out for version- or installation-specific differences) for histogramming, a complete program might look like

C...Common blocks.

C...Reserve histogram memory and book histograms.
      CALL HLIMIT(10000)
      CALL HBOOK1(1,’pT spectrum of pi+’,100,0.,5.,0.)

C...Number of events to generate. Loop over events.
      DO 110 IEVT=1,NEVT

C...Generate event. List first one.
      CALL LUEEVT(0,91.2)

C...Find sphericity axis and rotate event so sphericity along z axis.
      CALL LUEDIT(31)

C...Loop over all particles, but skip if not pi+.
      DO 100 I=1,N
      IF(K(I,2).NE.211) GOTO 100

C...Calculate pT and fill in histogram.
      CALL HF1(1,PT,1.)

C...End of particle and event loops.

C...Normalize histogram properly and list it.
      CALL HOPERA(1,’+’,1,1,20./NEVT,0.)


Study this program, try to understand what happens at each step, and run it to check that it works. You should then be ready to look at the relevant sections of this report and start writing your own programs.

3.6 Getting Started with PYTHIA

A Pythia run has to be more strictly organized than a Jetset one, in that it is necessary to initialize the generation before events can be generated, and in that it is not possible to change switches and parameters freely during the course of the run. A fairly precise recipe for how a run should be structured can therefore be given.

Thus, the usage of Pythia can be subdivided into three steps.

  • The initialization step. It is here that all the basic characteristics of the coming generation are specified. The material in this section includes the following.

    • Common blocks, at least the following, and maybe some more:

    • Selection of required processes. Some fixed ‘menus’ of subprocesses can be selected with different MSEL values, but with MSEL=0 it is possible to compose ‘à la carte’, using the subprocess numbers. To generate processes 14, 18 and 29, for instance, one needs

    • Selection of kinematics cuts in the CKIN array. To generate hard scatterings with 5 GeV 10 GeV, for instance, use


      Unfortunately, initial- and final-state radiation will shift around the kinematics of the hard scattering, making the effects of cuts less predictable. One therefore always has to be very careful that no desired event configurations are cut out.

    • Definition of underlying physics scenario, e.g. top mass.

    • Selection of parton-distribution sets, definitions, and all other details of the generation.

    • Switching off of generator parts not needed for toy simulations, e.g. fragmentation for parton level studies.

    • Initialization of the event generation procedure. Here kinematics is set up, maxima of differential cross sections are found for future Monte Carlo generation, and a number of other preparatory tasks carried out. Initialization is performed by PYINIT, which should be called only after the switches and parameters above have been set to their desired values. The frame, the beam particles and the energy have to be specified.

            CALL PYINIT(’CMS’,’p’,’pbar’,1800.)
    • Any other initial material required by the user, e.g. histogram booking.

  • The generation loop. It is here that events are generated and studied. It includes the following tasks:

    • Generation of the next event, with

            CALL PYEVNT
    • Printing of a few events, to check that everything is working as planned, with

            CALL LULIST(1)
    • An analysis of the event for properties of interest, either directly reading out information from the LUJETS common block or making use of a number of utility routines in Jetset.

    • Saving of events on tape, or interfacing to detector simulation.

  • The finishing step. Here the tasks are:

    • Printing a table of deduced cross sections, obtained as a by-product of the Monte Carlo generation activity, with the command

            CALL PYSTAT(1)
    • Printing histograms and other user output.

To illustrate this structure, imagine a toy example, where one wants to simulate the production of a 300 GeV Higgs particle. In Pythia, a program for this might look something like the following.

C...Common blocks.

C...Number of events to generate. Switch on proper processes.

C...Select t and H masses and kinematics cuts in mass.

C...For simulation of hard process only: cut out unnecessary tasks.

C...Initialize and list partial widths.
      CALL PYINIT(’CMS’,’p’,’p’,16000.)
      CALL PYSTAT(2)

C...Book histograms.
      CALL HLIMIT(10000)
      CALL HBOOK1(1,’Higgs mass’,50,275.,325.,0.)

C...Generate events. Look at first few.
      DO 200 IEV=1,NEV

C...Loop over particles to find Higgs and histogram its mass.
      DO 100 I=1,N
  100 IF(K(I,2).EQ.25) HMASS=P(I,5)
      CALL HF1(1,HMASS,1.)

C...Print cross sections and histograms.
      CALL PYSTAT(1)


Here 102, 123 and 124 are the three main Higgs production graphs , , and , and MSUB(ISUB)=1 is the command to switch on process ISUB. Full freedom to combine subprocesses ‘à la carte’ is ensured by MSEL=0; ready-made ‘menus’ can be ordered with other MSEL numbers. The PMAS commands set the masses of the top quark and the Higgs itself, and the CKIN variables the desired mass range of the Higgs — a Higgs with a 300 GeV nominal mass actually has a fairly broad Breit–Wigner type mass distribution. The MSTP switches that come next are there to modify the generation procedure, in this case to switch off initial- and final-state radiation, multiple interactions among beam jets, and fragmentation, to give only the ‘parton skeleton’ of the hard process. The PYINIT call initializes Pythia, by finding maxima of cross sections, recalculating the Higgs decay properties (which depend on the Higgs mass), etc. The decay properties can be listed with PYSTAT(2).

Inside the event loop, PYEVNT is called to generate an event, and LULIST(1) to list the event. The information used by LULIST(1) is the event record, stored in the common block LUJETS. Here one finds all produced particles, both final and intermediate ones, with information on particle species and event history (K array), particle momenta (P array) and production vertices (V array). In the loop over all particles produced, 1 through N, the Higgs particle is found by its code, K(I,2)=25, and its mass is stored in P(I,5).

After all events have been generated, PYSTAT(1) gives a summary of the number of events generated in the various allowed channels, and the inferred cross sections.

In the run above, a typical event listing might look like the following.

                         Event listing (summary)

    I  particle/jet     KF    p_x      p_y      p_z       E        m

    1  !p+!           2212    0.000    0.000 8000.000 8000.000    0.938
    2  !p+!           2212    0.000    0.000-8000.000 8000.000    0.938
    3  !g!              21   -0.505   -0.229   28.553   28.558    0.000
    4  !g!              21    0.224    0.041 -788.073  788.073    0.000
    5  !g!              21   -0.505   -0.229   28.553   28.558    0.000
    6  !g!              21    0.224    0.041 -788.073  788.073    0.000
    7  !H0!             25   -0.281   -0.188 -759.520  816.631  300.027
    8  !W+!             24  120.648   35.239 -397.843  424.829   80.023
    9  !W-!            -24 -120.929  -35.426 -361.677  391.801   82.579
   10  !e+!            -11   12.922   -4.760 -160.940  161.528    0.001
   11  !nu_e!           12  107.726   39.999 -236.903  263.302    0.000
   12  !s!               3  -62.423    7.195 -256.713  264.292    0.199
   13  !c~!             -4  -58.506  -42.621 -104.963  127.509    1.350
   14  (H0)             25   -0.281   -0.188 -759.520  816.631  300.027
   15  (W+)             24  120.648   35.239 -397.843  424.829   80.023
   16  (W-)            -24 -120.929  -35.426 -361.677  391.801   82.579
   17  e+              -11   12.922   -4.760 -160.940  161.528    0.001
   18  nu_e             12  107.726   39.999 -236.903  263.302    0.000
   19  s         A       3  -62.423    7.195 -256.713  264.292    0.199
   20  c~        V      -4  -58.506  -42.621 -104.963  127.509    1.350
   21  ud_1      A    2103   -0.101    0.176 7971.328 7971.328    0.771
   22  d         V       1   -0.316    0.001  -87.390   87.390    0.010
   23  u         A       2    0.606    0.052   -0.751    0.967    0.006
   24  uu_1      V    2203    0.092   -0.042-7123.668 7123.668    0.771
                sum:  2.00     0.00     0.00     0.00 15999.98 15999.98

The above event listing is abnormally short, in part because some columns of information were removed to make it fit into this text, in part because all initial- and final-state QCD radiation, all non-trivial beam jet structure, and all fragmentation was inhibited in the generation. Therefore only the skeleton of the process is visible. In lines 1 and 2 one recognizes the two incoming protons. In lines 3 and 4 are incoming partons before initial-state radiation and in 5 and 6 after — since there is no such radiation they coincide here. Line 7 shows the Higgs produced by fusion, 8 and 9 its decay products and 10–13 the second-step decay products. Up to this point lines give a summary of the event history, indicated by the exclamation marks that surround particle names (and also reflected in the K(I,1) code, not shown). From line 14 onwards come the particles actually produced in the final states, first in lines 14–16 particles that subsequently decayed, which have their names surrounded by brackets, and finally the particles and jets left in the end, including beam remnants. Here this also includes a number of unfragmented jets, since fragmentation was inhibited. Ordinarily, the listing would have gone on for a few hundred more lines, with the particles produced in the fragmentation and their decay products. The final line gives total charge and momentum, as a convenient check that nothing unexpected happened. The first column of the listing is just a counter, the second gives the particle name and information on status and string drawing (the A and V), the third the particle-flavour code (which is used to give the name), and the subsequent columns give the momentum components.

One of the main problems is to select kinematics efficiently. Imagine for instance that one is interested in the production of a single with a transverse momentum in excess of 50 GeV. If one tries to generate the inclusive sample of events, by the basic production graphs , then most events will have low transverse momenta and will have to be discarded. That any of the desired events are produced at all is due to the initial-state generation machinery, which can build up transverse momenta for the incoming and . However, the amount of initial-state radiation cannot be constrained beforehand. To increase the efficiency, one may therefore turn to the higher-order processes and , where already the hard subprocess gives a transverse momentum to the . This transverse momentum can be constrained as one wishes, but again initial- and final-state radiation will smear the picture. If one were to set a cut at 50 GeV for the hard-process generation, those events where the was given only 40 GeV in the hard process but got the rest from initial-state radiation would be missed. Not only therefore would cross sections come out wrong, but so might the typical event shapes. In the end, it is therefore necessary to find some reasonable compromise, by starting the generation at 30 GeV, say, if one knows that only rarely do events below this value fluctuate up to 50 GeV. Of course, most events will therefore not contain a above 50 GeV, and one will have to live with some inefficiency. It is not uncommon that only one event out of ten can be used, and occasionally it can be even worse.

If it is difficult to set kinematics, it is often easier to set the flavour content of a process. In a Higgs study, one might wish, for example, to consider the decay , with each or . It is therefore necessary to inhibit all other and decay channels, and also to adjust cross sections to take into account this change, all of which is fairly straightforward. However, if one wanted to consider instead the decay , with a meson producing a lepton, not only would there then be the problem of different leptonic branching ratios for different :s (which means that fragmentation and decay treatments would no longer decouple), but also that of additional pair production in parton-shower evolution, at a rate that is unknown beforehand. In practice, it is therefore impossible to force decay modes in a consistent manner.

4 Monte Carlo Techniques

Quantum mechanics introduces a concept of randomness in the behaviour of physical processes. The virtue of event generators is that this randomness can be simulated by the use of Monte Carlo techniques. In the process, the program authors have to use some ingenuity to find the most efficient way to simulate an assumed probability distribution. A detailed description of possible techniques would carry us too far, but in this section some of the most frequently used approaches are presented, since they will appear in discussions in subsequent sections. Further examples may be found e.g. in [Jam80].

First of all one assumes the existence of a random number generator. This is a (Fortran) function which, each time it is called, returns a number in the range between 0 and 1, such that the inclusive distribution of numbers is flat in the range, and such that different numbers are uncorrelated. The random number generator that comes with Jetset is described at the end of this section, and we defer the discussion until then.

4.1 Selection From a Distribution

The situation that is probably most common is that we know a function which is non-negative in the allowed range . We want to select an ‘at random’ so that the probability for a given is proportional to . Here might be a fragmentation function, a differential cross section, or any of a number of distributions.

One does not have to assume that the integral of is explicitly normalized to unity: by the Monte Carlo procedure of picking exactly one accepted value, normalization is implicit in the final result. Sometimes the integral of does carry a physics content of its own, as part of an overall weight factor we want to keep track of. Consider, for instance, the case when represents one or several phase-space variables and a differential cross section; here the integral has a meaning of total cross section for the process studied. The task of a Monte Carlo is then, on the one hand, to generate events one at a time, and, on the other hand, to estimate the total cross section. The discussion of this important example is deferred to section 7.4.

If it is possible to find a primitive function which has a known inverse , an can be found as follows (method 1):


The statement of the first line is that a fraction of the total area under should be to the left of . However, seldom are functions of interest so nice that the method above works. It is therefore necessary to use more complicated schemes.

Special tricks can sometimes be found. Consider e.g. the generation of a Gaussian . This function is not integrable, but if we combine it with the same Gaussian distribution of a second variable , it is possible to transform to polar coordinates


and now the and distributions may be easily generated and recombined to yield . At the same time we get a second number , which can also be used. For the generation of transverse momenta in fragmentation, this is very convenient, since in fact we want to assign two transverse degrees of freedom.

If the maximum of is known, in the range considered, a hit-or-miss method will always yield the correct answer (method 2):

  • select an with even probability in the allowed range, i.e. ;

  • compare a (new) with the ratio ; if , then reject the value and return to point 1 for a new try;

  • otherwise the most recent value is retained as final answer.

The probability that is proportional to ; hence the correct distribution of retained values. The efficiency of this method, i.e. the average probability that an will be retained, is . The method is acceptable if this number is not too low, i.e. if does not fluctuate too wildly.

Very often does have narrow spikes, and it may not even be possible to define an . An example of the former phenomenon is a function with a singularity just outside the allowed region, an example of the latter an integrable singularity just at the and/or borders. Variable transformations may then be used to make a function smoother. Thus a function which blows up as for , with an close to 0, would instead be roughly constant if transformed to the variable .

The variable transformation strategy may be seen as a combination of methods 1 and 2, as follows. Assume the existence of a function , with over the range of interest. Here is picked to be a ‘simple’ function, such that the primitive function and its inverse are known. Then (method 3):

  • select an according to the distribution , using method 1;


  • otherwise the most recent value is retained as final answer.

This works, since the first step will select with a probability and the second retain this choice with probability . The total probability to pick a value is then just the product of the two, i.e. .

If has several spikes, method 3 may work for each spike separately, but it may not be possible to find a that covers all of them at the same time, and which still has an invertible primitive function. However, assume that we can find a function , such that over the range considered, and such that the functions each are non-negative and simple, in the sense that we can find primitive functions and their inverses. In that case (method 4):

  • select an at random, with relative probability given by the integrals

  • for the selected, use method 1 to find an , i.e.


  • otherwise the most recent value is retained as final answer.

This is just a trivial extension of method 3, where steps 1 and 2 ensure that, on the average, each value picked there is distributed according to : the first step picks with relative probability , the second with absolute probability (this is one place where one must remember to do normalization correctly); the product of the two is therefore and the sum over all gives back .

We have now arrived at an approach that is sufficiently powerful for a large selection of problems. In general, for a function which is known to have sharp peaks in a few different places, the generic behaviour at each peak separately may be covered by one or a few simple functions , to which one adds a few more to cover the basic behaviour away from the peaks. By a suitable selection of the relative strengths of the different ’s, it is possible to find a function that matches well the general behaviour of , and thus achieve a reasonable Monte Carlo efficiency.

The major additional complication is when is a multidimensional variable. Usually the problem is not so much itself, but rather that the phase-space boundaries may be very complicated. If the boundaries factorize it is possible to pick phase-space points restricted to the desired region. Otherwise the region may have to be inscribed in a hyper-rectangle, with points picked within the whole hyper-rectangle but only retained if they are inside the allowed region. This may lead to a significant loss in efficiency. Variable transformations may often make the allowed region easier to handle.

There are two main methods to handle several dimensions, each with its set of variations. The first method is based on a factorized ansatz, i.e. one attempts to find a function which is everywhere larger than , and which can be factorized into , where . Here each may in its turn be a sum of functions , as in method 4 above. First, each is selected independently, and afterwards the ratio is used to determine whether to retain the point.

The second method is useful if the boundaries of the allowed region can be written in a form where the maximum range of is known, the allowed range of only depends on , that of only on and , and so on until , whose range may depend on all the preceding variables. In that case it may be possible to find a function that can be integrated over through to yield a simple function of , according to which is selected. Having done that, is selected according to a distribution which now depends on , but with through integrated over. In particular, the allowed range for is known. The procedure is continued until is reached, where now the function depends on all the preceding values. In the end, the ratio is again used to determine whether to retain the point.

4.2 The Veto Algorithm

The ‘radioactive decay’ type of problems is very common, in particular in parton showers, but it is also used, e.g. in the multiple interactions description in Pythia. In this kind of problems there is one variable , which may be thought of as giving a kind of time axis along which different events are ordered. The probability that ‘something will happen’ (a nucleus decay, a parton branch) at time is described by a function , which is non-negative in the range of values to be studied. However, this naïve probability is modified by the additional requirement that something can only happen at time if it did not happen at earlier times . (The original nucleus cannot decay once again if it already did decay; possibly the decay products may decay in their turn, but that is another question.)

The probability that nothing has happened by time is expressed by the function and the differential probability that something happens at time by . The basic equation then is


For simplicity, we shall assume that the process starts at time , with .

The above equation can be solved easily if one notes that :


and thus


With this is nothing but the textbook formulae for radioactive decay. In particular, at small times the correct decay probability, , agrees well with the input one, , since the exponential factor is close to unity there. At larger , the exponential gives a dampening which ensures that the integral of never can exceed unity, even if the integral of does. The exponential can be seen as the probability that nothing happens between the original time 0 and the final time . In the parton-shower language, this is (almost) the so-called Sudakov form factor.

If has a primitive function with a known inverse, it is easy to select values correctly:


which has the solution


If is not sufficiently nice, one may again try to find a better function , with for all . However to use method 3 with this would not work, since the method would not correctly take into account the effects of the exponential term in . Instead one may use the so-called veto algorithm:

  • start with and ;

  • add 1 to and select