12 February 1993
Updates to
PYTHIA 5.6 and JETSET 7.3
Physics and Manual
Torbjorn Sjostrand
Theory Division, CERN
CH-1211 Geneva 23
Switzerland
--------------------------------------------------
1) Introduction.
Since the PYTHIA/JETSET programs are being updated frequently,
it is also important to keep the documentation up to date.
The big manual that appeared as CERN-TH.6488/92 described the
status as of 20 May 1992. The LaTeX file with this manual will
be updated, though less frequently than the programs. (One
revised version has appeared, which documents changes made
up until 6 September.) Further, it is not very economical to
have to get a new copy of a long manual just to see if any
interesting new features have been added recently. Therefore
I have here collected the main changes that have taken place
after 20 May 1992. This file will be updated as regularly as
the programs themselves. It is not intended to be as complete
as the ordinary manual, but should be sufficient for the
intended purpose.
--------------------------------------------------
2) Changes in JETSET 7.3
16 June: LUHEPC modified to avoid problems with event history
when pileup events are generated in PYTHIA.
17 June: add option 16 to LUEDIT. With a CALL LUEDIT(16), some
missing daughter pointers in K(I,4) amd K(I,5) will be
reconstructed from the mother pointers in K(I,3).
17 June: when 3-jet events are generated according to the scalar
gluon option of LUEEVT, the full angular distribution (including
weak effects) is included according to the formulae in Ref. Lae80.
26 June: a techni-eta is introduced as particle 38, with default
mass PMAS(38,1) equal to 350 GeV and the three (dominant) decay
modes b + bbar, t + tbar and g + g. So as not to mess up the ordering
of existing decay channels, these modes are put in lines 1201 - 1203.
The techni-eta has zero spin, is neutral, singlet under weak SU(2),
but a colour octet.
10 July: introduce code 110 for a diffractive state of a pi0, a rho0,
or a gamma.
29 July: in LUSHOW do not set history pointers in K(I,4) and K(I,5)
for photon.
31 July: skip calculation of momentum shift in LUBOEI for pairs
already so close that numerical precision problems arise.
2 September: modify expression for maximum weight in tau decay in
LUDECY to make generation more efficient.
17 October: introduce extra parameter to determine action of LUCELL
calls. Cells with transverse energy below a threshold PARU(58) are
not considered, i.e. the energy in those cells is set vanishing.
By default PARU(58)=0., i.e. the old behaviour is retained.
Also correct mistake in LUCLUS, which gave a wrong value for the
minimum distance between two clusters as stored in PARU(63) (does not
affect the algorithm as such, only this complementary piece of
information).
28 October: introduce possibility for b quark to develop parton
shower in a top hadron decay, i.e. T -> W + b + spectator. Is
intended for scenarios where the top quark first fragments and only
later decays; cf. PYTHIA note of 14 October. This possibility is
regulated by switch MSTJ(27), with default value 0 (for backwards
compatibility). The possibilities are:
MSTJ(27) = 0 : (default) no parton shower, i.e. old behaviour.
= 1 : b parton shower; the W is thereafter assumed to decay
isotropically.
= 2 : b parton shower; the W thereafter decays anisotropically,
as in option 0, but with W momentum reduced when the
b parton shower initiator is off the mass shell.
Differences between options 1 and 2 are minor; in principle 2 is the
better but 1 is better matched to what is done when one uses option
MSTP(48)=1 and lets top decay before fragmentation.
29 October: correct a minor bug in LUZDIS (alternative fragmentation
function f(z) = (1-z)**c could not be reached).
3 November: change the default values for heavy flavour fragmentation
functions to correspond to the Peterson parametrizations with values
epsilon_c = 0.06, epsilon_b = 0.006 and epsilon_t = 0.00001. This
corresponds to a change of PARJ(54) = -0.06, PARJ(55) = -0.006, and
PARJ(56) - PARJ(58) = -0.00001. To enable these values, it is
recommended to use MSTJ(11)=3.
--------------------------------------------------
3) Changes in PYTHIA 5.6
2 June: modification of the scheme for preservation of (x, Q2) in
leptoproduction. As before, this is achieved by a post facto
modification of kinematics. The main problem is that modification may
often be impossible, due to kinematical constraints. Earlier such
events were simply thrown, which meant some bias in the generated
distributions. In the new scheme, the modification is done somewhat
differently, so that the failure rate has gone down significantly.
In addition, in case of failure, up to five new tries are made from
the same hard scattering kinematics, i.e. nominal (x, Q2) values,
but with new initial and final state shower evolution. The
preservation option should therefore be more useful.
22 June: introduce a call to LUEDIT(16) in PYEVNT so as to reconstruct
a few lost parts of event history (see JETSET, 17 June).
22 June: introduce parametrizations of the gamma-p total cross-section
in PYXTOT. The variable MSTP(30) can be used to select between the
two alternatives, as follows.
MSTP(30) = 1 : (default) the fit by Landshoff, see LP-HEP proceedings.
= 2 : the CERN-HERA fit, see Review of Particle Properties.
See further 8 July.
23 June: correct trivial bug which gave too high rate for the decay
of supersymmetric Higgs particles to tau+ tau-.
24 June: introduce (simple) description of multiple interactions in
gamma-p events, using the same formalism as for p-p events. Thus you
can use MSEL=1 to generate a combination of high-pT and low-pT
events, where a high-pT event may actually contain several
parton-parton interactions. The switches MSTP(81) and so on work
exactly as in p-p events. The main limitation is that this only
works for a photon beam, i.e. the formalism can not be used for
photoproduction events in e-p scattering. (The main reason is that
the multiple interactions formalism currently needs to be initialized
at a fixed c.m. energy of the hadronic system, while the e-p
event generation gives a spectrum of hadronic c.m. energies.)
29 June: implement process 149, g + g -> eta_techni (see also JETSET,
26 June). Three decay modes are included, b + bbar, t + tbar and g + g.
If kinematically allowed, the t + tbar decay mode always dominates.
The two main free parameters are the techni-eta mass and the decay
constant F_pi. The latter appears inversely quadratically in all
the partial widths. Since the total cross-section is proportional to
the g + g partial width, F_pi also affects the size of the
cross-section. F_pi is stored in PARP(46), and has the default value
123 GeV. Further information about the model implemented can be found
in T. Appelquist and G. Triantaphyllou, Yale University preprint,
in preparation.
3 July: expand process 77 so it also works for like sign longitudinal
W scattering, W+ W+ -> W+ W+ and W- W- -> W- W-.
6 July: introduce further option MSTP(22)=4 for parton shower
evolution scale in deep inelastic scattering chosen to be
Q2 * (1-x) * max(1, ln(1/x)) (see Ingelman, HERA workshop).
6 July: modify simulation of charged Higgs to include the down-type
part of quark couplings in process 161, and to include the running of
down-type mass in processes 143 and 161 (if requested by MSTP(37)).
7 July: correct an error in the MSTP(13)=2 (non-default) option of
PYSTEL, whereby the structure functions of partons inside a photon
inside an electron were overestimated by about 10% - 20%.
7 July: add an option 3 to a PYFRAM call, whereby you can boost a
leptoproduction event to the hadronic c.m. frame. Mainly intended
for deep inelastic scattering, but can also be used in photoproduction.
Note that the lepton, and any photons radiated by it, are not removed
from the event record; if you only want to study the hadrons, you have
to remove it yourself.
8 July: add options MSTP(30)=0 for gamma-p events and MSTP(31)=0 for
hadron-hadron ones, wherein you may put your own values for
cross-sections, according to the following:
PARP(21) - PARP(26) : (D=6*0.) cross-sections and slope parameters
that you can set yourself in the options MSTP(30)=0 and MSTP(31)=0,
respectively. If these options are used, at least PARP(21) and
PARP(25) must be assigned nonvanishing values, or the program is
likely to crash. All values are taken to refer to the c.m. energy of
the current run, i.e. the energy dependence of parameters need not be
specified (unlike the generic parametrizations found in the programs).
PARP(21) : total cross-section sigma_tot in mb.
PARP(22) : elastic cross-section sigma_el in mb.
PARP(23) : single diffractive cross-section sigma_sd in mb.
PARP(24) : double diffractive cross-section sigma_dd in mb.
PARP(25) : nuclear slope parameter B in GeV^(-2).
PARP(26) : curvature parameter C in GeV^(-4).
Note added January 1993: when the VDM description of gamma-p events
was expanded on 23 September, this section was not expanded accordingly.
In particular, with the new description it is also necessary to
specify how cross-sections are distributed among rho0, omega and phi
VDM states. After the PYINIT call, the cross-sections in VINT(101) -
VINT(107) should therefore be split among the new variables VINT(241) -
VINT(267) in the relative rho0/omega/phi proportions. In other words,
one must select VINT(241) + VINT(251) + VINT(261) = VINT(101), etc.
Further, one should put VINT(248) = VINT(258) = VINT(268) = VINT(121),
and MINT(101) = 3 if the beam partile (in the PYINIT call) is a photon,
while MINT(102) = 3 if the target particle is a photon.
8 July: include elastic and diffractive cross-sections for gamma-p
interactions. You can therefore use MSEL=2, as for p-p, to get
a mixture of elastic, single diffractive, double diffractive and
ordinary non-diffractive events. It is not possible, in the same
run, to generate also the direct contributions, where the photon
interacts pointlike (see MSTP(14)).
A summary of the current status:
For gamma-p events, where the photon acts like a hadron in the
vector dominance model sense, the default parametrization of the
total cross-section is taken from Landshoff. The cross-section for
the elastic process gamma p to rho^0 p has been parametrized as
sigma_el(s) = sigma_tot(s) * 0.058 * (s / 1 GeV^2)^0.08,
based on the numbers predicted by Schuler and Terron.
For lack of any other information, the single diffractive cross-section
has been assumed equal to the elastic one, and the double diffractive
to a third of this. The latter follows from the relation suggested by
Goulianos
sigma_dd(s) = sigma_sd^2(s) / (3 * sigma_el(s)) .
In elastic scattering, the t distribution is always parametrized as
d(sigma)/dt = const. * exp(B * t + C * t^2).
where B is the slope parameter and C the curvature. For gamma-p, one
assumes
B = 9.2 GeV^(-2) + 0.5 GeV^(-2) * ln(0.25 GeV^(-2) * s) ,
which is a simplified version of the form given by Schuler and Terron,
while the curvature parameter C has been assumed vanishing.
10 July: replace the existing elastic and diffractive kinematics
generation machinery by a new one, as follows:
For elastic scattering, t is generated according to the formula given
above. If the curvature C is vanishing or negative, then t is allowed
to go up to the kinematical limit. For a positive C, i.e. for a
distribution which eventually diverges, the distribution is cut off at
t = -B/(2C).
For diffractive scattering, the distribution in t is the
square-root of the distribution for elastic events, i.e. B and C
are rescaled by a factor 1/2. For each diffractive system an
additional factor dM^2/M^2 is multiplied on to give the mass M
of the excited state. Initially the two or three variables t,
M_1^2 and M_2^2 are selected independently of each other, within
the full allowed range. However, a multiplet of variables which is
internally kinematically inconsistent is then rejected and a new
one generated in its place. This leads to a slight shift of the
average abs(t) value upwards and a non-negligible suppression of
the high-mass tails of the diffractive systems.
A rho^ formed by gamma to rho^0 in elastic or diffractive scattering
is normally longitudinally polarized, and therefore its decay angular
distribution in rho^0 to pi^+ pi^- is taken to be proportional to
sin^2(theta), where the reference axis is given by the rho^0 direction
of motion.
The mass spectrum of the diffractive system is assumed to start
0.2 GeV above the mass of the incoming particle, or 0.2 GeV above the
rho^ mass for an incoming gamma. This value is stored in PARP(102)
and could be changed, but must be at least a pi^0 mass (plus some
saftey margin). A light diffractive system, with a mass less than 1 GeV
(PARP(103)) above the mass of the incoming particle, is taken to decay
isotropically into a two-body state. Single-resonance diffractive
states, such as a Delta^+, are therefore not explicily generated, but
are assumed described in an average, smeared-out sense.
A more massive diffractive system is treated as before, i.e. the
'Pomeron' can either interact with a quark or with a gluon, and thus
pull out either a single or a double string (see MSTP(101)).
6 August: modify PYREMN so that check on remnant mass is not made when
there is no remnant.
7 August: modify PYREMN so that energy sharing in resolved photon
remnant is as for meson; see MSTP(92) etc.
8 August: modify PYSTFU and PYSTEL so that the PDFLIB photon structure
functions may be used. To use this feature, put MSTP(56)=2 and select
PDFLIB set in MSTP(55). As before, you need to enable lines beginning
with C! to allow the proper calls to be made.
12 August: replace PDFLIB calls to STRUCTF by calls to STRUCTM, to
allow for a difference between up and down sea quarks. Also allow
option MSTP(57)=0 (structure functions evaluated at lower cutoff scale)
for PDFLIB.
31 August: modify PYKMAP to protect from numerical precision problem in
2 -> 3 processes.
23 September: introduce new parametrizations for total and partial
cross-sections, the former according to Donnachie-Landshoff
(CERN-TH.6635/92), the latter according to Schuler-Sjostrand (paper
in preparation). The partial cross-sections include elastic, single
diffractive on either of the two sides, and double diffractive.
The new parametrizations are automatically default for
gamma-p events, MSTP(30)=1, while the old parametrizations are
retained for the alternative, MSTP(30)=2. For backwards compatibility
the default is unchanged for hadron-hadron collisions, and the new
option is added as MSTP(31)=6. In principle you are recommended
to use this option, since e.g. it reproduces the Tevatron total and
elastic cross-sections much better. You should still be warned that
some parameters, e.g. for multiple interactions, may need to be retuned
somewhat for this alternative.
In addition, the treatment of diffractive system is further improved,
in two respects.
(1) The slope parameter B is evaluated in Pomeron language, which gives
predictions that depend on the nature of the incoming particles, on the
c.m. energy, and on the mass(es) of the diffractive system(s),
separately for elastic, single and double diffractive events.
(2) In photoproduction, the photon may be represented either by a
rho0, an omega or a phi, in the appropriate proportions. The couplings
are inversely proportional to the expressions f_V^2/(4pi), which are
stored in PARP(161) - PARP(163) (default values 2.20, 23.6, 18.4) for
rho0, omega and phi, respectively.
Further, the default value of PARP(102) is changed from 0.2 to 0.28
GeV; this represents the minimum mass allowed for a diffractive system
above the mass of the incoming particle.
For internal use, a number of MINT and VINT variables are now used,
including MINT(92), MINT(101) - MINT(107), VINT(107) and VINT(241) -
VINT(271).
25 September: introduce two new options MSTP(14) = 2 and 3 for the
structure of a photon. The options now include:
MSTP(14) (D=1) : nature of photon beam.
= 0 : direct, i.e. without internal structure.
= 1 : resolved, i.e. all interactions where a photon structure
function is needed.
= 2 : resolved, VDM only, i.e. where the photon structure function
can be taken to be an appropriately scaled-down version of
the rho0 (= pi0) one.
= 3 : resolved, anomalous only, i.e. the additional perturbatively
calculable part of the photon structure function initiated by
branchings gamma -> q + qbar.
In principle, options 2 and 3 together correspond to option 1, modulo
the fact that existing structure function parametrizations do not give
a 1-to-1 correspondence. However, the separate treatments in 2 and 3
allows a more precise modelling e.g. of beam jet structure, than can be
obtained with the more rough option 1. This will be described in a
paper by Schuler and Sjostrand, in preparation. However, it should be
remembered that large uncertainties remain. A full description of
photon interactions therefore require separate runs, either two with
MSTP(14) = 0 and 1, or three with MSTP(14) = 0, 2 and 3. Kindly note
that the current description only allows for one incoming photon, i.e.
does not yet cover the gamma-gamma case, where e.g. one photon is
direct and another is resolved (according to either of the options).
The structure function routine PYSTFU is modified to call the pi+
structure function for MSTP(14)=2, and turn it into a pi0 one which
then is scaled down. For MSTP(14)=3 a new routine PYSTAG is called by
PYSTFU. This routine (plus help routines PYSTGS and PYDILN), adapted
from G. Schuler, contains the anomalous part of the photon structure
function. A nominal lower cutoff scale P0 is needed, at which the
anomalous structure function is assumed vanishing. This cutoff P0 is
stored in PARP(15), and has the default value 0.5 GeV. For charm and
bottom flavours the threshold is assumed to appear at
Q^2 = PARP(16) * m_q^2. By default PARP(16) = 1. The actual scale P^2,
at which the photon is resolved, is logarithmically distributed between
the lower cutoff P0^2 and the scale Q^2 of the hard interaction. The
current P^2 value is stored in VINT(272). In the initial state
radiation machinery, an option MSTP(66) = 1 is included by default,
wherein the P^2 scale is used as an additional lower cutoff scale for
QCD radiation. The option MSTP(66) = 0 will always do backwards
evolution down to the standard cutoff PARP(62).
1 October: introduce routine PYSTFL, which allows a modification of the
low-Q^2 and low-x behaviour of the proton structure function. This
option is accessed with the switch value MSTP(57)=2. It can be used
in conjunction with any internal or PDFLIB proton structure function.
Its main effect is to make structure functions vanish in the limit
Q^2 -> 0.
5 October: introduce option to dampen the cross-section of an arbitrary
hard 2 -> 2 process in the limit pT -> 0, in the same way as is done
with QCD processes in options MSTP(82) >= 2 of the multiple interactions
scenario. More precisely, the cross-section is multiplied by
pT^4/(pT0^2 +pT^2)^2, where pT0 = PARP(82). In addition alpha_S is
evaluated at the scale pT0^2 + pT^2 rather than pT^2. This possibility
is switched on with MSTP(85) = 1; default is MSTP(85) = 0, i.e. off.
14 October: introduce option to allow top quarks to decay before they
fragment. To use this option, put MSTP(48)=1. Additionally you have to
change to MDCY(6,1)=1. Default is MSTP(48)=0, MDCY(6,1)=0, i.e. top
quarks fragment into top hadrons. Currently only the decay mode
t -> b + W is implemented. Note that, in the new mode, the top decays
to a W with code +-24, while for top hadron decay the separate code
89 (Wvirt) has been used. This is of importance e.g. if one wishes to
select only semileptonic decays. The treatment of top is still not that
of the standard resonances, in that cross-section for top production
is not reduced when the subsequent W decay is limited to a few modes;
you have to multiply with the respective branching rations yourself.
Also kindly note that spin correlations are not included, i.e. the t
and the W both decay isotropically in their respective rest frames.
28 October: As consequence of JETSET change one can now also include
parton shower evolution in PYTHIA. The following summarizes the recent
top changes:
If a top is lighter than about 120 GeV, then the top quark is
sufficiently long-lived that top hadrons T are formed, which
subsequently decay. For heavier tops, the decay is faster, and
therefore hadrons have no time to be formed. Of course there is
a transition region where the issue is not so clear - where things
will work out differently in each event. So, in principle, one
should try both scenarios to check there are no serious disagreement
between them. So far as I can tell, they indeed seem to match well.
Earlier, only the first alternative was available. Now you can use
MSTP(48) to switch between the two.
MSTP(48) : (D=0) flag to be set.
= 0: as before; top quark fragments to top hadrons which then
decay. Recommended for top mass below 120 GeV.
= 1: new; top decays to W + b, and the b subsequently fragments.
Recommended for top mass above 120 GeV. You must also change
MDCY(6,1)=1 to use this option, see below.
Depending on which alternative you choose above, there are other
switches that can be set. The first deals with QCD radiation off
the b quark. This was negligible when the top was assumed to weigh
below 100 GeV, but is more significant for the heavier tops now
used. The main effect of such gluon emission is to soften the
momentum of the B quarks produced, and add a few extra soft pions.
In the option MSTP(48)=0, you can switch on this emission with
MSTJ(27) = 1 or 2. The default is MSTJ(27)=0, i.e. no emission,
for reasons of backwards compatibility, so you are recommended to
change. When MSTP(48) = 1, QCD radiation is automatically included
as part of standard, and need not be selected. (It can be switched
off, together with other final state radiation, by putting
MSTP(71) = 0.) Some theoretical ambiguity exists in the modelling
of this gluon emission, but what is offered should be a reasonable
first approximation.
The t -> b + W is followed by the W decay. In principle this decay
carries some spin information, i.e. is not quite isotropic. This
information is not available in the MSTP(48) = 1 option i.e. here
the W decays isotropically. For MSTP(48) = 0 it is made use of
when MSTJ(27) = 2, while MSTJ(27) = 1 implies isotropic W decay.
In practice, effects are small, to the best of my knowledge.
For MSTP(48) = 0 top hadrons are produced and subsequently decay.
Therefore the decay channels of the top are set for the top hadron
pseudoparticle 86. For MSTP(48) = 1 no hadrons are formed, and
it is instead the t quark itself, with code 6, that decays. Since
this was not originally foreseen, you must change to MDCY(6,1)=1.
In scheme MSTP(48) = 0 the decay W is in fact the Wvirt with code
89, while the standard W with code 24 is used for MSTP(48) = 1.
The separate code 89 is related to the necessity to be more
sophisticated if the t mass is below or around the W one, and is
therefore now more of a historical curiosity. It implies, however,
that you have to select the W decay channels differently in the
two scenarios.
If you print cross-sections with PYSTAT(1), then remember that these
do not include the top and W branching ratios, for the case that
you only allow some channels. You therefore have to multiply with
these yourself. This holds true in either option.
So far, I have not implemented the potential t -> b + H mode. In
principle, this is not difficult, but did not get top priority.
It is recommended to switch to Peterson fragmentation for the heavy
quarks and keep the standard Lund one for the other hadrons. This
you do by putting MSTJ(11) = 3. I have included as defaults the
recommended epsilon values epsilon_c = 0.06, epsilon_b = 0.006 and
epsilon_t = 0.00001 in the new Jetset version, so you need not set
these yourself. The c and b values are typical LEP fit results,
with QCD emission included. Note that there is some mismatch of
Q2 scales between e+e- and pp production of heavy flavours, so the
epsilon values above are not holy numbers, but should be a good first
guess.
29 October: trace how the position of a t is moved by shower evolution,
so that option MSTP(48)=1 takes effect also e.g. in Z0 -> t + tbar.
14 December: fix numerical precision problem in PYDILN.
13 January: correct mass for gamma + p -> phi + p in PYXTOT.
Constrain acos evaluation in PYI3AU to avoid math library bug on SUN.
29 January: possibility to run PYTHIA with varying energies from one
event to the next. Here comes a summary of the new features:
It is now possible to use PYTHIA in a mode where the energy can be
varied from one event to the next, without the need to reinitialize
with a new PYINIT call. This allows a significant speed-up of
execution, although it is not as fast as running at a fixed
energy. It can not be used for everything - we will come to
the fine print at the end - but should be applicable for most tasks.
The master switch to access this possibility is in MSTP(171).
By default it is off, so you must set MSTP(171)=1 before
initialization. There are two submodes of running,
MSTP(172)= 1 or 2. In the former mode, PYTHIA will generate an
event at the requested energy. This means that you have to know
which energy you want beforehand. In the latter mode, PYTHIA
will often return without having generated an event - with flag
MSTI(61)=1 to signal that - and you are then requested to give
a new energy. The energy spectrum of accepted events will then,
in the end, be your naive input spectrum weighted with the
cross-section of the processes you study. We will come back to
this.
The energy can be varied, whichever frame is given in the PYINIT
call. In addition two new frames are allowed, to increase the
flexibility. The summary is:
CALL PYINIT(FRAME,BEAM,TARGET,WIN)
FRAME : frame of experiment.
= 'CMS' : WIN gives nominal c.m. energy used at initialization.
For each event to be generated thereafter, PARP(171) should
be filled by the fractional energy of that one, i.e.
E_cm = PARP(171)*WIN. Therefore PARP(171) should normally
be smaller than unity.
= 'FIXT' : WIN gives initialization beam energy, with target at
rest. For each event PARP(171) should be filled by the
fractional beam energy of that one, i.e. E_beam = PARP(171)*WIN.
Therefore PARP(171) should normally be smaller than unity.
= 'USER' : give in beam and target three-momenta in
P(1,1) - P(1,3) and P(2,1) - P(2,3), respectively. These
variables have to be given at initialization and for each
new event. The particle masses are assumed to be the nominal
ones and the particle energies are calculated by the program.
The c.m. energy for events to be generated should normally
be smaller than the one used at initialization.
= 'FOUR' : as 'USER', except that you should give P(1,1) - P(1,4)
and P(2,1) - P(2,4), i.e. also the energies of the two
incoming particles. These particles need not be on the
mass shell.
= 'FIVE' : as 'USER', except that you should give P(1,1) - P(1,5)
and P(2,1) - P(2,5), i.e. also the energies and masses of the
two incoming particles. These particles need not be on the
mass shell. For a particle with m**2 = E**2 - p**2 < 0, one
should give P(I,5) = -sqrt(-m**2). Since the mass should be
uniquely given by the energy and momentum, the only advantage
of this option compared with 'FOUR' is that a small mass
may be better known by the user than the program can calculate
itself, due to the limited precision.
To illustrate the use of the MSTP(172)=2 facility, consider the case of
beamstrahlung in e+e- linear colliders. This is just for convenience;
what is said here can be translated easily into other situations.
Assume that the beam spectrum is given by D(z), where z is the
fraction retained by the original e after beamstrahlung. Therefore
0 <= z <=1 and the integral of D(z) is unity. This is not perfectly
general; one could imagine branchings e- -> e- + gamma -> e- + e+ + e-,
which gives a multiplication in the number of beam particles. This
could either be expressed in terms of a D(z) with integral larger than
unity or in terms of an increased luminosity. We will assume the latter,
and use D(z) properly normalized. Given a nominal s = 4*E_beam**2,
the actual s' after beamstrahlung is given by s' = z_1 * z_2 * s.
For a process with a cross-section sigma(s) the total cross-section is
sigma_tot = Int_0^1 Int_0^1 D(z_1) D(z_2) *
* sigma(s' = z_1 * z_2 * s) dz_1 dz_2.
The cross-section sigma may in itself be an integral over a number
of additional phase space variables. If the maximum of the differential
cross-section is known, a correct procedure to generate events is
1) pick z_1 and z_2 according to D(z_1) and D(z_2);
2) pick a set of phase space variables of the process, for the given
s' of the event;
3) evaluate sigma(s') and compare with sigma_max;
4) if event is rejected, then return to 1) to generate new variables;
5) else continue generation to give complete event.
You as a user are assumed to take care of step 1), and present the
resulting kinematics with incoming e+ and e- of varying energy.
Thereafter PYTHIA will do steps 2 - 5, and either return an event or
put MSTI(61)=1 to signal failure.
The maximization procedure of PYTHIA does search in phase space to
find sigma_max, but it does not vary the s' energy in this process.
Therefore the maximum search in the PYINIT call should be performed
where the cross-section is largest. For processes with increasing
cross-section as function of energy this means at the largest energy
that will ever be encountered, i.e. s' = s in the case above. This is
the 'standard' case, but often one encounters other behaviours, where
more complicated procedures are needed. One such case would be the
process e+ e- -> Z* -> Z + H, which is known to have a cross-section
which increases near the threshold but is decreasing asymptotically.
If one already knows that the maximum, for a given Higgs mass,
appears at 300 GeV, say, then the PYINIT call should be made with that
energy, even if subsequently one will be generating events for a 500
GeV collider.
In general, it may be necessary to modify the selection of z_1 and
z_2 and assign a compensating event weight. For instance, consider
a process with a 1/s cross-section. Then the sigma_tot expression
above may be rewritten as
sigma_tot = Int_0^1 Int_0^1 D(z_1)/z_1 D(z_2)/z_2 *
* z_1 z_2 sigma(s' = z_1 * z_2 * s) dz_1 dz_2.
The expression z_1 z_2 sigma(s') is now flat in s', i.e. not only
can sigma_max be found at a convenient energy
such as the maximum one, but additionally the PYTHIA generation
efficiency (the likelihood of surviving step 4)) is greatly
enhanced. The price to be paid is that z has to be selected
according to D(z)/z rather than according to D(z). Note that
D(z)/z is not normalized to unity. One therefore needs to define
I_D = Int_0^1 D(z)/z and a properly normalized
D'(z) = (1/I_D) * D(z)/z.
Then
sigma_tot = Int_0^1 Int_0^1 D'(z_1) D'(z_2) *
* (I_D)**2 z_1 z_2 sigma(s' = z_1 * z_2 * s) dz_1 dz_2.
Therefore the proper event weight is WT = (I_D)**2 z_1 z_2.
This weight should be stored, for each event, in PARP(173).
Additionally you must put MSTP(173)=1 for the program to
make use of weights at all. Often D(z) are not known analytically;
therefore I_D is also not known beforehand, but may have to be
evaluated (by you) during the course of the run. Then you should
just use WT = z_1 z_2 in PARP(173) and do the overall normalization
yourself in the end. Only the cross-sections are affected by this,
the events themselves still are properly distributed in s' and
internal phase space.
Above it has been assumed tacitly that D(z) -> 0 for z -> 0. If not,
D(z)/z is divergent, and it is not possible to define a
properly normalized D'(z) = D(z)/z. If the cross-section is truly
diverging like 1/s, then a D(z) which is nonvanishing for z -> 0
does imply an infinite total cross-section, whichever way things
are considered. In cases like that, it may be necessary to impose
a lower cut on z, based on some physics or detector consideration.
The most difficult cases are those with a very narrow and high
peak, such as the Z0. One could initialize at the energy of maximum
cross-section and use D(z) as is, but efficiency might turn out to
be very low. One might then be tempted to do more complicated
transforms of the kind illustrated above. As a rule it is then
convenient to work in the variables tau_z = z_1 * z_2 and
y_z = (1/2) ln (z_1/z_2), cf. section 7.2 of the manual.
Clearly, the better the behaviour of the cross-section can be
modelled in the choice of z_1 and z_2, the better the overall
event generation efficieny. Even under the best of circumstances,
the efficiency will still be lower than for runs with fix energy.
There is also a non-negligible time overhead for using variable
energies in the first place, from kinematics reconstruction
and (in part) from the phase space selection. One should therefore
not use variable energies when not needed, and not use a large
range of energies sqrt(s') if in the end only a smaller range is of
experimental interest.
The following new common block variables therefore exist:
MSTP(171) (D=0) possibility of variable energies from one event to
the next.
= 0 : no; i.e. the energy is fixed at the initialization call.
= 1 : yes; i.e. a new energy has to be given for each new event.
MSTP(172) (D=2) options for variable energy generation, applicable
when MSTP(171)=1.
= 1 : an event is generated at the requested energy, i.e.
internally a loop is performed over possible event
configurations until one is accepted. If the requested c.m.
energy of an event is below PARP(2) the run is aborted.
Cross-section information can not be trusted with this
option, since it depends on how you decided to pick
the requested energies.
= 2 : only one event configuration is tried. If that is accepted,
the event is generated in full. If not, no event is
generated, and the status code MSTI(61)=1 is returned.
You are then expected to give a new energy, looping until an
acceptable event is found. No event is generated if the
requested c.m. energy is below PARP(2), instead status code
MSTI(61) is set =1 to signal the failure. In principle,
cross-sections should come out correctly with this option.
MSTP(173) (D=0) possibility for user to give in event weight to
compensate for a biased choice of beam spectrum.
= 0 : no, i.e. event weight is unity.
= 1 : yes; weight to be given for each event in PARP(173).
PARP(171) : to be set, event-by-event, when variable energies are
allowed, i.e. MSTP(171)=1. If PYINIT is called with
FRAME='CMS' (='FIXT'), PARP(171) multiplies the c.m. energy
(beam energy) used at initialization. For the options 'USER',
'FOUR' and 'FIVE', PARP(171) is dummy.
PARP(173) : event weight to be given by user when MSTP(173)=1.
MSTI(61) : status flag set when events are generated. It is only of
interest for runs with variable energies, MSTP(171)=1, with the
option MSTP(172)=2.
= 0 : an event has been generated.
= 1 : no event was generated, either because the c.m. energy was
too low or because the Monte Carlo phase space point
selection machinery rejected the trial point. A new energy
is to be picked by the user.
Also a few additional internal variables have been introduced:
MINT(111) : the frame given in PYINIT call, 0 - 5 for
'NONE', 'CMS', 'FIXT', 'USER', 'FOUR' and 'FIVE', respectively.
VINT(289) : squared c.m. energy found in PYINIT call.
VINT(290) : the WIN argument of a PYINIT call.
VINT(291) - VINT(300) : the two five-vectors of the two incoming
particles, as reconstructed in PYINKI or PYEVKI.
This facility may be combined with most other aspects of the program.
For instance, it is possible to simulate beamstrahlung as above and
still include bremsstrahlung with MSTP(11)=1. Further, one may
multiply the overall event weight of PARP(173) with a
kinematics-dependent weight given by PYEVWT, although it is not
recommended (since the chances of making a mistake are also
multiplied). However, a few things do NOT work.
* It is not possible to use pile-up events, i.e. you must have
MSTP(131)=0.
* The possibility of giving in your own cross-section optimization
coefficients, option MSTP(121) = 2, would require more input
than enumerated in the manual, and this option should therefore
not be used. You can still use MSTP(121) = 1, however.
* The multiple interactions scenario with MSTP(82) >= 2 only works
approximately for energies different than the initialization one.
If the c.m. energy spread is smaller than a factor 2, say, the
approximation should be reasonable, but if the spread is larger
one may have to subdivide into subruns of different energy bins.
The initialization should be made at the largest energy to be
encountered - whenever multiple interactions are possible (i.e.
for incoming hadrons and resolved photons) this is where
cross-sections are largest anyway, and so is no further constraint.
There is no simple possibility to change PARP(82) during the
course of the run, i.e. a energy-independent pT0 must be assumed.
The default option MSTP(82) = 1 works fine, i.e. does not suffer
from the constraints above. If so desired pTmin=PARP(81) can be
set differently for each event, as a function of c.m. energy.
11 February: inclusion of J/Psi in the VDM component of
photoproduction. It appears on equal footing with rho0, omega and
phi when a gamma-p event is generated with MSTP(14)=2. The
coupling f_J/Psi**2/4pi is stored in PARP(164). The following
points should be noted.
* The parametrization of the J/Psi + p total cross section is uncertain.
We use the phi-p cross-section as a reference and scale it down
by about a factor 10, in agreement with current estimates, but
there is an uncertainty by maybe a factor of 2 in this number.
What has been measured experimentally (at low energies) is only
the elastic gamma + p -> J/Psi + p cross-section, which is
proportional to the sigma_tot**2/f_J/Psi**2. In other words,
only a combination of sigma_tot and f_J/Psi is measured, which
means that a change of sigma_tot is allowed if f_J/Psi is rescaled.
* At HERA energies the VDM cross section goes up by about 1.4 mb by
the inclusion of this new component. This is so little that it does
not affect any of the discussion of photoproduction properties
in general, but it gives a non-negligible contribution to charm
production in general. If you only want to simulate this component,
i.e. switch off the rho0, omega and phi contributions, you can
do this by setting PARP(161) - PARP(163) to be very large, e.g.
1E10.
* The fraction of elastic J/Psi production is much lower than is the
case for rho0. This is expected, based on the smallness of the
total cross section (noted above) and the optical theorem. A J/Psi
obviously may appear in the final state of elastic and single
diffractive events (where the proton becomes a diffractive system).
A diffractive J/Psi state normally would produce e.g. D + Dbar
mesons, but if the mass of the diffractive state is below
this threshold, the program will allow diffractive J/Psi + pi0
states. Therefore also double diffractive events may contain J/Psi.
* For the nondiffractive events, one would need knowledge about the
J/Psi structure functions to do a good job. For lack of anything
else, I assume the same valence spectrum as for pions, which
certainly means the c quarks are too soft. Therefore one should be
careful not to trust the results too much, in particular not at
high transverse momenta. Most are produced in low-pT events,
where structure functions are not explicitly used, and where the
uncertainty is smaller.
* The J/Psi decays isotropically, i.e. no attempt has been made to
include the possibility of j/Psi polarization.
The manual should be corrected as follows.
PARP(164) : (D=11.5) the coupling f_J/Psi**2/4pi used in the VDM model
when MSTP(14)=2.
VINT(271) - VINT(280) : cross sections for the J/Psi VDM component,
parallelling the current description of VINT(241) - VINT(270).
VINT(281) : what used to be VINT(271), i.e. ratio between total
gamma+X and total pi0+X cross section.
VINT(282) : what used to be VINT(272), i.e. scale P**2 at which an
anomalous photon is resolved.