Normally users are expected to implement new processes via the
Les Houches Accord, or by supplying
Les Houches Event Files. Then
you do all flavour, colour and phase-space selection externally,
before your process-level events are input for further processing
by PYTHIA. However, it is also possible to implement a
new process in exactly the same way as the internal PYTHIA
ones, thus making use of the internal phase space selection machinery
to sample an externally provided cross-section expression.
The MadGraph5 program [Alw11] allows you to do exactly that,
i.e. it can be used to generate C++ code that can be linked into
the existing PYTHIA framework, see
Should you decide to go ahead on your own,
this page gives a brief summary how to do that. If you additionally
want to introduce a new resonance species, with its own internal
width calculations, you will find further instructions
here. It is strongly
recommended to shop around for a similar process that has already
been implemented, and to use that existing code as a template.
Look for processes with the same combinations of incoming flavours
and colour flows, rather than the shape of the cross section itself.
With a reasonable such match the task should be of medium difficulty,
without it more demanding.
PYTHIA's internal phase-space generators are rather good at handling
the phase space of 2 → 1 and 2 → 2 processes,
are more primitive for 2 → 3 ones and do not at all
address higher multiplicities. An option is therefore also provided
for external phase-space generators to be used, which must then be
encapsulated to inherit from PYTHIA's
- The Cross Section Calculation
- The Cross Section Class
- Access to a process
- Implementing an external phase-space generator
PhaseSpace base class
(or one of its derivatives).
The set of processes that can be implemented in this framework is
therefore in principle unlimited, though the user must supply external
phase-space generators for non-trivial 2 → 3 processes
and all higher 2 → n multiplicities.
Note, however, that the produced particles may be resonances, so it is
possible to end up with bigger "final" multiplicities through sequential
decays, also with the internal phase-space generators, and to include
further matrix-element weighting in those decays.
For processes using PYTHIA's internal phase-space generators, there
are three steps involved in implementing a process:
We consider these aspects in turn. An example where it all comes
together is found in
- making use of the PYTHIA-provided kinematics information to
calculate the relevant cross section,
- writing a new class, where the matrix elements are implemented,
including information on incoming and outgoing flavours and colours,
- making the process available.
For processes for which an external phase-space generator will be
used, step 1 above changes to writing a new class, where the
phase-space generator is implemented, and making use of that to
calculate the relevant cross section. There are no example programs
illustrating how to do this yet, but the methodology is described
below, under "Implementing an external phase-space generator".
The Cross Section Calculation
The key method for the cross section calculation is
SigmaProcess::sigmaHat(), described below. At the point when
it is called, the kinematics has already been set up, and from these
phase space variables the differential cross section is to be calculated.
For a 2 → 1 process, the returned value should be
mH (= mHat),
sH (= sHat) and
sH2 (= sHat^2)
are available to be used. Incoming partons are massless. Overload the
convertM2() method below if you instead plan to return
For a 2 → 2 process, instead d(sigmaHat)/d(tHat)
should be returned, based on provided
mH, sH, sH2, tH, tH2, uH, uH2, m3, s3, m4, s4 and
pT2 values (
s3 = m3*m3 etc.). Incoming
partons are massless. Overload the
below if you instead plan to return |M|^2.
For a 2 → 3 process, instead |M|^2 should be
returned, with normalization such that |M|^2 / (2 sHat) integrated
over the three-body phase space gives the cross section. Here no standard
set of Mandelstam-style variables exists. Instead the obvious ones,
mH, sH, m3, s3, m4, s4, m5, s5, are complemented by the
p3cm, p4cm, p5cm, from which further invariants
may be calculated. The four-vectors are defined in the CM frame of the
subcollision, with massless incoming partons along the +-z axis.
In either case, alpha_s and alpha_em have already
been calculated, and are stored in
Also other standard variables may be used, like
CoupEW::sin2thetaW(), and related flavour-dependent
vector and axial couplings in
CoupEW and CKM combinations
In case some of the final-state particles are resonances, their
squared masses have already been selected according to a Breit-Wigner
with a linearly running width Gamma(m) = Gamma(m_0) * m / m_0.
More precisely, the mass spectrum is weighted according to
w_BW(m^2) d(m^2), where
w_BW(m^2) = (1/pi) * (m * Gamma(m)) / ( (m^2 - m_0^2)^2 + (m * Gamma(m))^2 ) .
If you would like to have another expression, the above weights are stored
respectively. If you divide out one of these factors, you just remain with
a phase space selection d(m^2) for this particle,
and can multiply on your desired shape factor instead. Unfortunately, the
Monte Carlo efficiency will drop if your new mass distribution differs
dramatically from the input one. Therefore it does make sense to adjust the
database value of the width to be slightly (but not too much) broader
than the distribution you have in mind. Also note that, already by default,
the wings of the Breit-Wigner are oversampled (with a compensating lower
internal weight) by partly sampling like (a + b/m^2 + c/m^4) d(m^2),
where the last term is only used for gamma^*/Z^0.
As alternative to the kinematics variables defined above, also the two
pME, for masses and
four-momenta, respectively, can be used for cross-section calculations.
Here indices 0 and 1 are the two incoming beams, and 2 and onwards the
outgoing particles. Note that this differs by one step from the normal
internal labeling, where slot 0 is left empty. The four-momenta are
defined in the rest frame of the subcollision, with the incoming partons
along the +-z direction. The kinematics need not agree with the
"correct" one stored in the event record, for three reasons.
1) Gauge invariance forces matrix-element calculations to use
the same masses for incoming as outgoing legs of a particle species,
say b quarks. Therefore the kinematics of the two incoming
partons is recalculated, relative to the normal event record, to put
the partons on the mass shell. (Note that initial masses is a technical
issue, not the correct physics picture: the incoming partons are likely
to be spacelike virtual rather than on the mass shell.)
2) In principle each fermion flavour has to be treated separately,
owing to a different mass. However, in many cases fermions can be
assumed massless, which speeds up the calculations, and further gains
occur if then different flavours can use the same cross-section
expression. In MadGraph the default is that fermions up to and including
the c quark and the mu lepton are considered massless,
while the b quark and the tau lepton are considered
massive. This can be modified however, and below we provide four flags
that can be used to consider the "borderline" fermions either as
massless or as massive when matrix elements are evaluated, to match the
assumptions made for the matrix elements themselves.
3) For 2 → 2 and 2 → 3 processes of massive
identical particles (or antiparticles) in the final state, such as
t tbar or W^+ W^-, the kinematics is here adjusted
so that the two or three particles have the same mass, formed as a
suitable average of the actual Breit-Wigner-distributed masses. This
allows the evaluation of matrix-element expressions that only have
meaning if the two/three have the same mass.
Thus the mass array
mME and the four-momentum array
pME present values both for initial- and final-state
particles based on these mass principles suited for matrix-element input.
Note that these variables therefore differ from the kinematics stored in
the event record proper, where incoming fermions are always massless and
outgoing resonances have independent Breit-Wigner mass distributions.
The conversion from the normal to the special kinematics is done
by calling the
setupForME() method. This you have to do
yourself in the
SigmaHat() member of your derived class.
Alternatively it could be done in
SigmaKin(), i.e. before
the loop over incoming flavours, but then these would be considered
massless. The identity of final-state particles is obtained from the
methods. Should the conversion to
pME not work,
setupForME() will return
false, and then the cross section should be put zero.
default = off)
Let the c quark be massive or not in the kinematics set up for
external matrix-element evaluation.
default = on)
Let the b quark be massive or not in the kinematics set up for
external matrix-element evaluation.
default = off)
Let the mu lepton be massive or not in the kinematics set up for
external matrix-element evaluation.
default = on)
Let the tau lepton be massive or not in the kinematics set up for
external matrix-element evaluation.
The Cross Section Class
The matrix-element information has to be encoded in a new class.
The relevant code could either be put before the main program in the
same file, or be stored separately, e.g. in a matched pair
.cc files. The latter may be more
convenient, in particular if the cross sections are lengthy, or if you
intend to build up your own little process library, but of course
requires that these additional files are correctly compiled and linked.
The class has to be derived either from
Sigma1Process, for 2 → 1 processes, from
Sigma2Process, for 2 → 2 ones, or from
Sigma3Process, for 2 → 3 ones. (The
Sigma0Process class is used for elastic, diffractive
and minimum-bias events, and is not recommended for use beyond that.)
These are in their turn derived from the
The class can implement a number of methods. Some of these are
compulsory, others strongly recommended, and the rest are to be
used only when the need arises to override the default behaviour.
The methods are:
A constructor for the derived class obviously must be available.
Here you are quite free to allow a list of arguments, to set
the parameters of your model, or even to create a set of closely
related but distinct processes. For instance, g g → Q Qbar,
Q = c or b, is only coded once, and then the
constructor takes the quark code (4 or 5) as argument,
to allow the proper amount of differentiation.
A destructor is only needed if you plan to delete the process
before the natural end of the run, and require some special behaviour
at that point. If you call such a destructor you will leave a pointer
dangling inside the
Pythia object you gave it in to,
if that still exists.
is called once during initialization, and can then be used to set up
parameters, such as masses and couplings, and perform calculations
that need not be repeated for each new event, thereby saving time.
This method needs not be implemented, since in principle all
calculations can be done in
is called once a kinematical configuration has been determined, but
before the two incoming flavours are known. This routine can therefore
be used to perform calculations that otherwise might have to be repeated
over and over again in
sigmaHat below. For instance
a flavour-independent cross section calculation for a q g
initial state would be repeated 20 times in
five times for the five quark flavours allowed in the incoming beams,
times twice to include antiquarks, times twice since the (anti)quark
could be in either of the two beams. You could therefore calculate the
result once only and store it as a private data member of the class.
It is optional whether you want to use this method, however, or put
is the key method for cross section calculations and returns a cross section
value, as described in the previous section. It is called when also a
preliminary set of incoming flavours has been picked, in addition to the
kinematical ones already available for
sigmaHat is called inside a loop over all allowed
incoming flavour combinations, stored in
id2, with fixed kinematics, as already illustrated above.
The sum over the different flavour combinations provides the total
cross section, while their relative size is used to make a selection of
a specific incoming state.
to be called by the user from inside
sigmaKin()) to setup alternative kinematics
pME arrays, better
suited for matrix-element calculations. See the end of the previous
section for a more detailed description. Should the method return
false then the conversion did not work, and
sigmaKin()) should be set to
is called only once an initial state and a kinematical configuration has
been picked. This routine must set the complete flavour information and
the colour flow of the process. This may involve further random choices,
between different possible final-state flavours or between possible
competing colour flows. Private data members of the class may be used to
retain some information from the previous steps above.
When this routine is called the two incoming flavours have already
been selected and are available in
whereas the one, two or three outgoing ones either are fixed for a given
process or can be determined from the instate (e.g. whether a W^+
or W^- was produced). There is also a standard method in
VCKM to pick a final flavour from an initial one with CKM
mixing. Once you have figured out the value of
id3 and, the case being,
id5, you store these values permanently by a call
setId( id1, id2, id3, id4, id5), where the last two may be
omitted if irrelevant.
Correspondingly, the colours are stored with
setColAcol( col1, acol1, col2, acol2, col3, acol3, col4, acol4,
col5, acol5), where the final ones may be omitted if irrelevant.
Les Houches style colour tags are used, but starting with number 1
(and later shifted by the currently requested offset). The
input is grouped particle by particle, with the colour index before the
anticolour one. You may need to select colour flow dynamically, depending
on the kinematics, when several distinct possibilities exist. Trivial
operations, like swapping colours and anticolours, can be done with
methods have been used, the order of the outgoing particles may be
inconsistent with the way the tHat and uHat
variables have been defined. A typical example would be a process like
q g → q' W with tHat defined between incoming and
outgoing quark, but where
id3Mass() = 24 and so the
process is to be stored as q g → W q'. One should then put
swapTU = true in
for each event where the tHat and uHat variables
should be swapped before the event kinematics is reconstructed. This
variable is automatically restored to
false for each new
double SigmaProcess::weightDecayFlav( Event& process)
is called to allow a reweighting of the simultaneous flavour choices of
resonance decay products. Is currently only used for the
q qbar → gamma*/Z^0 gamma*/Z^0 process, and will likely not
be of interest for you.
double SigmaProcess::weightDecay( Event& process, int iResBeg, int iResEnd)
is called when the basic process has one or several resonances, after each
set of related resonances in
has been allowed to decay. The calculated weight, to be normalized
to the range between 0 and 1, is used to decide whether to accept the
decay(s) or try for a new decay configuration. The base-class version of
this method returns unity, i.e. gives isotropic decays by default.
This method may be called repeatedly for a single event. For instance, in
q qbar → H^0 Z^0 with H^0 → W^+ W^-, a first call
would be made after the H^0 and Z^0 decays, and then
depend only on the Z^0 decay angles since the H^0
decays isotropically. The second call would be after the W^+ W^-
decays and then involve correlations between the four daughter fermions.
returns the name of the process, as you want it to be shown in listings.
returns an integer identifier of the process. This has no internal function,
but is only intended as a service for the user to rapidly (and hopefully
uniquely) identify which process occurred in a given event. Numbers below
10000 are reserved for internal PYTHIA use.
this string specifies the combinations of incoming partons that are
allowed for the process under consideration, and thereby which incoming
calls will be looped over. It is always possible to pick a wider flavour
selection than strictly required and then put to zero cross sections in
the superfluous channels, but of course this may cost some extra execution
time. Currently allowed options are:
gg: two gluons.
qg: one (anti)quark and one gluon.
qq: any combination of two quarks, two antiquarks or
a quark and an antiquark.
qqbar: any combination of a quark and an antiquark;
a subset of the
qqbarSame: a quark and its antiquark;
a subset of the
ff: any combination of two fermions, two antifermions
or a fermion and an antifermion; is the same as
hadron beams but also allows processes to work with lepton beams.
ffbar: any combination of a fermion and an antifermion;
is the same as
qqbar for hadron beams but also allows processes
to work with lepton beams.
ffbarSame: a fermion and its antifermion; is the
qqbarSame for hadron beams but also allows processes
to work with lepton beams.
ffbarChg: a fermion and an antifermion that combine
to give charge +-1.
fgm: a fermion and a photon (gamma).
ggm: a gluon and a photon.
gmgm: two photons.
it is assumed that cross sections normally come in dimensions such that
they, when integrated over the relevant phase space, obtain the dimension
GeV^-2, and therefore need to be converted to mb. If the cross section
is already encoded as mb then
convert2mb() should be
overloaded to instead return
it is assumed that 2 → 1 cross sections are encoded as
sigmaHat(sHat), and 2 → 2 ones as
d(sigmaHat)/d(tHat) in the
convertM2() is overloaded to instead return
true then the return value is instead assumed to be the
squared matrix element |M|^2, and
SigmaProcess::sigmaHatWrap(...) converts to
sigmaHat(sHat) or d(sigmaHat)/d(tHat), respectively.
This switch has no effect on 2 → 3 processes, where
|M|^2 is the only allowed input anyway.
are the one, two or three final-state flavours, where masses are to be
selected before the matrix elements are evaluated. Only the absolute value
should be given. For massless particles, like gluons and photons, one need
not give anything, i.e. one defaults to 0. The same goes for normal light
quarks, where masses presumably are not implemented in the matrix elements.
Later on, these quarks can still (automatically) obtain constituent masses,
once a u, d or s flavour has been selected.
are the codes of up to two s-channel resonances contributing to
the matrix elements. These are used by the program to improve the phase-space
selection efficiency, by partly sampling according to the relevant
Breit-Wigner distributions. Massless resonances (the gluon and photon)
need not be specified.
normally the choice of renormalization and factorization scales in
2 → 2 and 2 → 3 processes is based on the
assumption that t- and u-channel exchanges dominates the
cross section. In cases such as f fbar → gamma* → f' fbar'
a 2 → 2 process actually ought to be given scales as a
2 → 1 one, in the sense that it proceeds entirely through
an s-channel resonance. This can be achieved if you override the
false to return
true. See further the
page on couplings and scales.
normally no intermediate state is shown in the event record for
2 → 2 and 2 → 3 processes. However, in case
idSChannel is overloaded to return a nonzero value,
an intermediate particle with that identity code is inserted into the
event record, to make it a 2 → 1 → 2 or
2 → 1 → 3
process. Thus if both
are overloaded, a process will behave and look like it proceeded through
a resonance. The one difference is that the implementation of the
matrix element is not based on the division into a production and a
decay of an intermediate resonance, but is directly describing the
transition from the initial to the final state.
there are two different 3-body phase-space selection machineries,
of which the non-QCD one is default. If you overload this method
instead the QCD-inspired machinery will be used. The differences
between these two is related to which
phase space cuts
can be set, and also that the QCD machinery assumes (almost) massless
the non-QCD 2 → 3 phase space selection machinery is rather
primitive, as already mentioned. The efficiency can be improved in
processes that proceed though t-channel exchanges, such as
q qbar' → H^0 q qbar' via Z^0 Z^0 fusion, if the
identity of the t-channel-exchanged particles on the two side
of the event are provided. Only the absolute value is of interest.
in the above kind of 2 → 3 phase-space selection, the
sampling of pT^2 is done with one part flat, one part weighted
like 1 / (pT^2 + m_R^2) and one part like
1 / (pT^2 + m_R^2)^2. The above values provide the relative
amount put in the latter two channels, respectively, with the first
obtaining the rest. Thus the sum of
tChanFracPow2() must be below unity. The final results
should be independent of these numbers, but the Monte Carlo efficiency
may be quite low for a bad choice. Here m_R is the mass of the
exchanged resonance specified by
idTchan2(). Note that the order of the final-state
listing is important in the above q qbar' → H^0 q qbar' example,
i.e. the H^0 must be returned by
since it is actually the pT^2 of the latter two that are
selected independently, with the first pT then fixed
by transverse-momentum conservation.
in 2 → 3 processes the phase space selection used here
involves a twofold ambiguity basically corresponding to a flipping of
the positions of last two outgoing particles. These are assumed equally
likely by default,
false, but for processes proceeding entirely
through t-channel exchange the Monte Carlo efficiency can be
improved by making a preselection based on the relative propagator
allows a possibility to override the global mode
for a specific process. The global mode normally is used to switch off
parts of the gamma^*/Z^0 propagator for test purposes. The
above local mode is useful for processes where a Z^0 really is
that and nothing more, such as q qbar → H^0 Z^0. The default
value -1 returned by
gmZmode() ensures that the global
mode is used, while 0 gives full gamma^*/Z^0 interference,
1 gamma^* only and 2 Z^0 only.
Access to a process
Once you have implemented a class, it is straightforward to make use of
it in a run. Assume you have written a new class
which inherits from
Sigma3Process, which in their turn inherit from
SigmaProcess. You then create an instance of this class
and hand it in to a
pythia object with
SigmaProcess* mySigma = new MySigma();
If an external phase-space generator should be used for this process
(see "Implementing an external phase-space generator" below),
this should be specified as a second argument in the call to
setSigmaPtr(), as in:
pythia.setSigmaPtr( new mySigma(), new myPhaseSpaceGenerator() );
If you have several processes you can repeat the procedure any number
of times. When
pythia.init() is called these processes
are initialized along with any internal processes you may have switched on,
and treated in exactly the same manner. The
will therefore generate a mix of the different kinds of processes without
distinction. See also the Program Flow
If the code should be of good quality and general usefulness, it would
be simple to include it as a permanently available process in the
standard program distribution. The final step of that integration ought to
be left for the PYTHIA authors, but here is a description of what is
A flag has to be defined, that allows the process to be switched on;
by default it should always be off. The name of the flag should be
chosen of the type
model:process. Here the
model would be related to the general scenario considered,
specify instate and outstate, separated by a 2 (= to), e.g.
When several processes are implemented and "belong together" it is
also useful to define a
model:all switch that affects
all the separate processes.
The flags should normally be stored in the
file or one of its daughters for a specific kind of processes. This is to
make them easily found by users. You could create and use your own
.xml file, so long as you then add that name to the
list of files in the
Index.xml file. (If not,
the flags would never be created and the program would not work.)
ProcessContainer.c file, the
SetupContainers::init() method needs to be expanded to
create instances of the processes switched on. This code is fairly
repetitive, and should be easy to copy and modify from the code
already there. The basic structure is
(i) check whether a process is requested by the user and, if so,
(ii) create an instance of the matrix-element class,
(iii)create a container for the matrix element and its associated
phase-space handling, and
(iv) add the container to the existing process list.
Two minor variations are possible. One is that a set of related
processes are lumped inside the the same initial check, i.e. are
switched on all together. The second is that the matrix-element
constructor may take arguments, as specified by you (see above).
If so, the same basic matrix element may be recycled for a set of
related processes, e.g. one for a composite u and one for
a composite d. Obviously these variations may be combined.
Implementing an external phase-space generator
An external phase-space generator can be interfaced by encapsulating
it within a class inheriting from PYTHIA's
base class. The following three virtual methods must be defined:
// Determine how phase space should be sampled.
virtual bool setupSampling();
// Select a trial event kinematics.
virtual bool trialKin(bool inEvent = true, bool repeatSame = false);
// Construct final (accepted) event kinematics.
virtual bool finalKin();
Optionally, a further virtual method is available to specify whether
beam particles are resolved in partons or scatter directly,
// Inform whether beam particles are resolved or scatter directly.
virtual bool isResolved();
with default return value
setupSampling() step the main point is to determine
the upper estimate of the cross section integrated over the allowed
phase space regions, and this should be stored in
The ratio between the correct cross section and its upper estimate
is a measure of the phase-space selection efficiency, and the purpose
of this step is to optimize the sampling accordingly. To this end any
convenient set of phase-space variables may be chosen. The
x2H varables should be used to
denote the incoming parton momentum fractions, however, to be used in
trialKin() intermediate step the same set of internal
variables can be used, and fed into the
to evaluate the cross section in the given phase space point, multiplied
by the integrated cross section. This value is to be stored in
sigmaNw, and the ratio
will be used to determine whether the trial event is accepted or not.
finalKin() step the output is more standardized.
The key values are the ones stored in the
pH arrays, the former for masses and the latter for
four-momenta. Here the first two slots represent the two incoming
partons and the subsequent ones up to ten outgoing particles. Other
particle properties, like the number of final-state particles, their
identities and colours, and more, are defined by the
A tailor-made 2 → 3
generator could be defined, e.g., by starting from the code for PYTHIA's
PhaseSpace2to3tauycyl base class, which provides
a specific representation of 3-parton phase space, used for generic
2 → 3 processes in PYTHIA. The virtual functions
described could then be redefined to generate a different sampling of
3-parton phase space. One example of this is provided by the existing
PhaseSpace2to3yyycyl class, which PYTHIA uses for
massless QCD processes. Note the interplay between the phase-space
variables, generated and saved here, and how they are used by the
matrix-element codes. For general processes, the user can define
samplings in terms of their own phase-space parametrizations, as long
as the corresponding matrix elements use the same variables to
evaluate the cross-section expressions.