Semi-Internal Processes
Normally users are expected to implement new processes via the
Les Houches Accord. 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.
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.
Should you actually go ahead, 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 is rather good at handling the phase space of
2 -> 1 and 2 -> 2 processes, is more primitive for
2 -> 3 ones and does not at all address higher multiplicities.
This limits the set of processes that you can implement in this
framework. The produced particles may be resonances, however, so it is
possible to end up with bigger "final" multiplicities through sequential
decays, and to include further matrix-element weighting in those decays.
There are two steps involved in implementing a process:
1) writing a new class, where the matrix elements are implemented,
including information on incoming and outgoing flavours and colours, and
2) making the process available.
We consider these two aspects in turn. An example where it all comes
together is found in main25.cc
.
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
of .h
and .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 SigmaProcess
base class.
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.
method
void initProc()
is called once during initalization, 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 sigmaHat
below.
method
void sigmaKin()
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 sigmaHat
,
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
everything in sigmaHat
.
method
double sigmaHat()
is the key method for cross section calculations and returns a cross section
value, as further described below. It is called when also a preliminary set
of incoming flavours has been picked, in addition to the kinematical ones
already available for sigmaKin
. Typically sigmaHat
is called inside a loop over all allowed incoming flavour combinations,
stored in id1
and 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 incomimg state.
For a 2 -> 1 process, the returned value should be
sigmaHat(sHat), where mH
(= mHat),
sH
(= sHat) and sH2
(= sHat^2)
are available to be used.
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.).
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 variables exist. Instead the obvious ones,
mH, sH, m3, s3, m4, s4, m5, s5
, are complemented by the
four-vectors p3cm, p4cm, p5cm
, from which further invariants
may be calculated. The four-vectors are defined in the cm frame of the
subcollision, with incoming partons along the +-z axis.
In either case, alpha_s and alpha_em have already
been calculated, and are stored in alpS
and alpEM
.
Also other standard variables may be used, like
CoupEW::sin2thetaW()
, and related flavour-dependent
vector and axial couplings in CoupEW
and CKM combinations
in VCKM
.
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
in runBW3
, runBW4
and runBW5
,
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.
method
void setIdColAcol()
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 id1
and id2
,
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, id4
and
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
existing methods.
When the id3Mass()
and id4Mass()
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
the variable swapTU = true
in setIdColAcol()
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
event.
method
double 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.
method
double 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 process[i]
,
iResBeg
<= i
<= iResEnd
,
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.
method
string name()
returns the name of the process, as you want it to be shown in listings.
method
int code()
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 occured in a given event. Numbers below
10000 are reserved for internal PYTHIA use.
method
string inFlux()
this string specifies the combinations of incoming partons that are
allowed for the process under consideration, and thereby which incoming
flavours id1
and id2
the sigmaHat()
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.
* qqbarSame
: a quark and its antiquark;
this is a subset of the above qq
option.
* ff
: any combination of two fermions, two antifermions
or a fermion and an antifermion; is the same as qq
for
hadron beams but also allows processes to work with lepton beams.
* ffbarSame
: a fermion and its antifermion; is the
same as 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.
method
bool convert2mb()
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 false
.
method
int id3Mass(), int id4Mass(), int id5Mass()
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.
method
int resonanceA(), int resonanceB()
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-Wigners. Massless resonances (the gluon and photon) need not be
specified.
method
bool isSChannel()
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
default false
to return true
. See further the
page on couplings and scales.
method
int idTchan1(), int idTchan2()
the 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.
method
double tChanFracPow1(), double tChanFracPow2()
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 tChanFracPow1()
and
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 idTchan1()
or
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 id3Mass()
,
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.
method
useMirrorWeight()
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
weights, true
.
method
int gmZmode()
allows a possibility to override the global mode
WeakZ0:gmZmode
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.
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 MySigma
,
which inherits from Sigma1Process
, Sigma2Process
or 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();
pythia.setSigmaPtr( mySigma);
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 pythia.next()
will therefore generate a mix of the different kinds of processes without
distinction. See also the Program Flow
description.
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
required.
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,
e.g. Compositeness
, while process
would
specify instate and outstate, separated by a 2 (= to), e.g.
ug2u*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 ProcessSelection.xml
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.)
In the 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.