\documentclass[10pt]{article}
\usepackage[margin=1.2in]{geometry}
\usepackage{titling}

\title{The EVENT2 Program}
\date{}
\author{}

\posttitle{\par\end{center}\vspace{-0.5in}}

\begin{document}

\maketitle

\begin{abstract}

This is a collection of notes about the \texttt{EVENT2} program, which implements the Catani-Seymour algorithm for next-to-leading order corrections to two- and three-jet event observables in $e^+e^-$ annihilation. It is not intended to describe the algorithm itself, which will one day be written up. It will hopefully evolve into proper program documentation at some stage though.

\end{abstract}

\section*{Overview}

The program is written somewhat like an event generator | a main routine administers the generation of parton momenta, calls the matrix element routines to calculate their weight, and finally calls a user routine to analyse them, with the momenta stored in something vaguely resembling an event record. At present, the program is designed just to calculate the perturbative coefficients themselves, and does not even contain an $\alpha_S$ calculation, although of course it would be trivial to add this. 

The most important point to note about the algorithm is that the singularities of the $m$-parton real matrix element are cancelled by subtracting several $(\sim m(m-1)(m-2)/2)$ $m-1$-parton counterterms. These are weighted such that when the $m$-parton term diverges, so do one or more $m-1$-parton terms, so that their diverence stays finite. Therefore an `event' consists of many separate sets of momenta with coupled weights, and it is important that the user routine takes this into account. If one is not interested in the errors on the calculated quantities this is trivial | events can be directly histogrammed. However, if error bars are needed it is important that all the sets of momenta are entered into the histogram in one go at the end of the event, so that large positive and negative weights have been cancelled already.

\section*{Phase-space generation}

Although we said that the remaining cross-section is finite, it is not actually. It is guaranteed to have finite integral, but can (and does in general) contain integrable singularities, most importantly square-root singularities. These are notorious for spoiling the convergence of Monte Carlo integrals (their variance is infinite). They must therefore be removed by suitable Jacobian factors. In \texttt{EVENT2} this is done by using a fairly sophisticated multi-channel system.

The phase-space generation looks at first sight something like a parton shower | first a two-parton event is generated, then, using this as a starting point a three-parton event is generated by choosing one of the partons to split into two while the other absorbs its recoil so that all three are on mass-shell. Finally a four-parton event is generated in the same way, choosing one of the two remaining partons to absorb the recoil. However, unlike most parton showers, this can be decribed as a proper transformation of the full four-parton phase-space such that the Jacobian factor can be easily calculated. Different choices of which parton emits and which absorbs the recoil correspond to the different channels of the multi-channel integration. In principle, the convergence might be improved by assigning different weights to different channels. Algorithms even exist to optimize this automatically, but in practice we find that the convergence is perfectly satisfactory simply assigning the same weight to all channels. As usual in multi-channel Monte Carlo, the Jacobian factor is the sum of the factors for all channels, not just the one that was used to generate the event.

The square-root singularities are removed by ensuring that for each variable $x$ in which the full matrix element is singular as $1/x$ (i.e. two-parton invariant mass-squareds and parton energies), $x$ is generated according to \[\frac{dx}{x^{1-1/n}},\] where $n$ is an adjustable parameter. Clearly $n = 1$ gives flat phase-space, $n = 2$ removes square-root singularities, while even larger $n$ pays even more attention to the almost-singular regions. Obviously the appropriate Jacobian factors are applied so that the final result does not depend on the value of $n$, only the size of the errors does. There are actually two adjustable parameters, \texttt{NPOW1} and \texttt{NPOW2}, both defaulting to 2, which control the $2\rightarrow3$ and $3\rightarrow4$ parton steps respectively. For most studies the default values are sufficient, but if for example one was studying the next-to-leading order corrections to a three-jet quantity like thrust in the extreme two-jet-like region, increasing \texttt{NPOW1} would improve the convergence. For technical reasons (certain constants are tabulated rather than being calculated from scratch) \texttt{NPOW1} and \texttt{NPOW2} are integer parameters.

The two-parton event is chosen according to the distribution \[d\cos\theta\,\texttt{COSP}(1+\cos\theta)^2 + \texttt{COSM}(1-\cos\vartheta)^2+\texttt{COSZ},\] where $\cos\theta$ is the angle between the quark and the electron\footnote{The truth is that I haven't checked the parity-violating parts of the code much yet, so it is possible that $\pm\cos\theta$ could be accidently interchanged.}, and \texttt{COSP}, \texttt{COSM} and \texttt{COSZ} are adjustable parameters. Since the matrix elements are not very forward-backward peaked, the exact values are not too critical, but ideally one should set $\texttt{COSZ} = 1$, $\texttt{COSP/M} = 0$ for unorientated matrix elements (see below), $\texttt{COSZ} = 0$, $\texttt{COSP/M} = 1$ (the default) for photon exchange and whatever for the full $Z/\gamma$ exchange.

\section*{Matrix elements}

There are three sets of matrix elements included. None of them is perfect at the moment (but we're working on it!). The first is the original \texttt{ERT} set ($\texttt{METYPE} = 0$), which average over event orientation. Full colour information and most flavour information is available (see below). The second is the set provided for us by the Leiden group ($\texttt{METYPE} = 1$). These are fully oriented matrix elements, but at present only include the photon exchange term. They include the same colour and avour information as the \texttt{ERT} set. The final set is the Giele-Glover matrix elements ($\texttt{METYPE} = 2$ believe it or not), basically ripped out of their \texttt{EERAD} program but with some modifications to use the \texttt{ERT} definitions of the pole terms and a few others associated with summing over parton permutations. No colour or flavour information is available, for technical rather than fundamental reasons. In the long term this is really intended as a cross-check that our code is working rather than as a serious set.

\section*{Flavour information}

For the four-parton matrix element the weight given to the user is always the sum of the contributions from q\={q}gg, qq\={q}\={q} and qq'\={q}\={q}'. The momenta are given in the `event record' in this order. If the cross-sections are needed for the individual processes, it is available (when using the \texttt{ERT} or Leiden matrix elements) from the variables \texttt{GGSUB}, \texttt{QQSUB} and \texttt{QPSUB} respectively, which hold the fraction of the current weight contributed by each subprocess. However, in the four-fermion case this information should be interpreted with a little care, as both sets of matrix elements neglect interference terms from which the contribution would be zero if quarks and antiquarks were averaged over, but which can be non-zero if both quark-antiquark pairs are distinguished. It is hard to imagine many physical observables for which this could be a problem (baryon-baryon (as opposed to baryon-antibaryon) correlations?). Nevertheless I think we should try to include these terms in the Leiden code, just to give a complete description of the process.

\section*{Colour information}

In the \texttt{ERT} and Leiden matrix elements, the variables \texttt{CFSUB}, \texttt{CASUB} and \texttt{TFSUB} give the relative contributions to the current weight from the $C_F$, $C_A$ and $T_f$ terms respectively. The colour factors themselves are removed, i.e. the normalization is \[C_F (C_F \texttt{CFSUB} + C_A\texttt{CASUB} + N_f T_R\texttt{TFSUB}) = 1.\]

\section*{The main routine}

The main routine is called \texttt{EVENT2}. 

\begin{verbatim}
    SUBROUTINE EVENT2(NEV,EM,NFL,USER)
    INTEGER NEV,NFL
    DOUBLE PRECISION EM
    EXTERNAL USER
\end{verbatim}

It takes four arguments,

\begin{description}
    \item[\texttt{NEV}] is the number of events to generate,
    \item[\texttt{EM}] is the total centre-of-mass energy,
    \item[\texttt{NFL}] is the number of avours (used both for the total e+ e! q\={q} normalization and in the second-order cross-section),
    \item[\texttt{USER}] is the name of a user-supplied analysis routine.
\end{description}

The program starts by setting various parameters. Most of these are extremely standard and would never need changing, but at some stage perhaps an improved interface in which the default parameters can be modifi ed should probably be provided. The internal parameters are described below.

\section*{The user routine}

The user must supply a routine to analyze events

\begin{verbatim}
    SUBROUTINE USER(N,NA,ITYPE,P,WEIGHT)
    INTEGER N,NA,ITYPE
    DOUBLE PRECISION P(4,7),WEIGHT
\end{verbatim}

\begin{description}
    \item[\texttt{N}] is the number of partons in the `event record',
    \item[\texttt{NA}] is the order of $\alpha_S$ (i.e. 0, 1 or 2),
    \item[\texttt{ITYPE}] is the type of contribution (shouldn't be needed for physical observables but provided anyway): 0=tree level, 1=subtraction counterterm, 2=`finite virtual' term,
    \item[\texttt{P}] is the primitive `event record': \texttt{P}($1-4, i$) is the momentum of the $i$th entry. Entries 1-\texttt{N} are the outgoing partons in the order given above for four-parton terms, and q\={q} and q\={q}g for the two- and three-parton terms. Entry 5 is their total four-momentum. Entries 6 and 7 are the incoming electron and positron.
    \item[\texttt{WEIGHT}] is the weight of the event.
\end{description}

The weights are normalized such that their total value is the final cross-section. Thus they do not need to be renormalized at the end of the run. At each order in $\alpha_S$ a factor of $\left(\frac{\alpha_S}{2\pi}\right)^\texttt{NA}$ is extracted. Thus about the simplest imagineable user routine would be something like this...

\begin{verbatim}
     SAVE A,B,C,D
     IF (NA.EQ.0) THEN
       A=A+WEIGHT
     RETURN
     ELSEIF (NA.EQ.1) THEN
       B=B+WEIGHT
     ENDIF

C---CALCULATE THE THRUST OF THE N-PARTON EVENT
    :
    :
     IF (NA.EQ.1) THEN
       C=C+(1-T)*WEIGHT
     ELSEIF (NA.EQ.2) THEN
       D=D+(1-T)*WEIGHT
     ENDIF
\end{verbatim}

At the end of the run, \texttt{A} would be equal to the Born cross-section for $e^+e^-\rightarrow q\bar{q}$ at the requested energy, \texttt{B} would be equal to 2\texttt{A}, indicating that the cross-section through next-to-leading order is \[\sigma=A\left(1+\frac{\alpha_S}{2\pi}2\right).\] \texttt{C} would be equal to 2.10\texttt{A} and \texttt{D} to 45\texttt{A}, indicating that the average value of thrust is given by \[\langle1-T\rangle=\frac{2.10\frac{\alpha_S}{2\pi}+45\left(\frac{\alpha_S}{2\pi}\right)^2}{1+\frac{\alpha_S}{\pi}}.\]

There is a more complex demonstration routine built in to the program, which uses Sjostrand's \texttt{GBOOK} package (similar to \texttt{HBOOK}) to calculate the distribution of various event shapes, as well as their average values. It is recommended that you create a new user routine of your own, rather than modifying the existing one.

The whole of the \texttt{EVENT2} program is Lorentz invariant. The momenta can be specified in any frame in routine \texttt{GENTWO}, which generates the original q\={q} pair, and the choice will carry through the whole procedure. However, the default frame is the $e^+e^-$ centre-of-mass, so it should be safe to assume that in the user routine. For die-hard dogmatists it is possible to write any event shape in a Lorentz-invariant way though, so it should be possible to avoid any assumption about the frame used.

\section*{The common blocks}

There are three common blocks. Two are intended to be internal, while the third gives extra information to the user routine. \texttt{CONCOM} holds all the constants used by the program,

\begin{verbatim}
    INTEGER METYPE,NF
    DOUBLE PRECISION CF,CA,TR,PI,PISQ,HF,CUTOFF,CQ,ONF
    COMMON /CONCOM/ CF,CA,TR,PI,PISQ,HF,CUTOFF,CQ,ONF,NF,METYPE
\end{verbatim}

All values are set at the beginning of \texttt{EVENT2}. At present there is no possibility to change them without recompiling. They are not copied to the \texttt{EERAD} common blocks, so if using the Giele-Glover matrix elements, some care is needed. The first few should be reasonably obvious, (with $\texttt{TR} = \texttt{HF} = \frac{1}{2}$).

\texttt{CUTOFF} is an infrared cut on the phase-space. If any pair of partons has $m^2_{ij} < \texttt{CUTOFF} Q^2$ the whole event is thrown away. Strictly speaking, the algorithm does not need an infrared cutoff | all integrals are finite, so it can be set to zero. However, for very small $m_{ij}$ one ends up taking the difference between two very large numbers and the numerical inaccuracy of the machine becomes important. The default value \texttt{CUTOFF} = $10^{-8}$ works for the Leiden matrix elements. For \texttt{ERT}, it could even be decreased to about $10^{-10}$ if necessary, but the Giele-Glover matrix elements seem to be extremely unreliable and something more like $10^{-6} $ is necessary. We have explicitly checked by working in quadruple precision that this is a numerical problem not a physics problem or a bug. A study of the \texttt{CUTOFF}-dependence is needed, but has not yet been done.

\texttt{CQ} holds the coupling constant, $\texttt{CQ} = \sum_q^{N_f}e^2_q$. It will eventually become an array once Z exchange is included. \texttt{NF} is just a copy of the user input \texttt{NFL} and \texttt{ONF} holds \texttt{1/NF} for convenience. Setting \texttt{NFL = 0} is safely dealt with\footnote{Actually it isn't at the moment. This is currently the only known bug.}. \texttt{METYPE} was already described above.

\texttt{SAMCOM} holds all the constants associated with importance sampling.

\begin{verbatim}
    INTEGER NPOW1,NPOW2
    DOUBLE PRECISION POWINT(10),XPOW1,XPOW2,COSP,COSM,COSZ
    COMMON /SAMCOM/ POWINT,XPOW1,XPOW2,COSP,COSM,COSZ,NPOW1,NPOW2
\end{verbatim}

The last five were already described above. \texttt{XPOW}$i$ is simply shorthand for 1-1/\texttt{NPOW}$i$, while \texttt{POWINT} pretabulates the definite integral \[\texttt{POWINT}(n)=\int_0^1\,dx\sqrt{1- x^n}=\frac{\sqrt{(\pi)}\Gamma(1/n)}{2^{2/n}n\Gamma(1/2+17n)}.\] Since it only extends to 10, \texttt{NPOW}$i$ can only be up to 10, but this is easily big enough for most imagineable applications.

\texttt{SUBCOM} holds the sub-process contributions to the current weight (for the $\alpha_S^2$ terms only) as described earlier,

\begin{verbatim}
    DOUBLE PRECISION CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB
    COMMON /SUBCOM/ CFSUB,CASUB,TFSUB,GGSUB,QQSUB,QPSUB
\end{verbatim}

\section*{EERAD common blocks}

The \texttt{EERAD} common blocks are set in routine \texttt{setup}. They are reasonably well explained in the comments. A couple of points worth noting are: \texttt{ymin} has to be set to a non-zero value even though it is never used; \texttt{as} is the strong coupling constant, which has to be set to $2\pi$ to match \texttt{EVENT2}'s normalization; \texttt{iorder} has to be left at 0 otherwise the subtraction terms will not match the real phase-space; there is no check that \texttt{nf} is the same as that used by \texttt{EVENT2}, but the integrals will not converge if it is not. \texttt{setup} prints lots of spurious messages, which can be ignored.

\section*{Timing}

To get errors similar to the LEP1 yellow book, we need about 20 million events\footnote{This refers to the middle of phase-space, we find significantly better convergence than them going into the two-jet region of phase-space.}. The main program, using ERT/Leiden/Giele-Glover matrix elements takes about 19/32/145 cpu hours to do this on CERNSP. In addition, the user routine takes about 2 cpu hours to calculate the thrust distribution, and 14 cpu hours for the demonstration routine, which calculates the distributions of C- and D-parameters, thrust, $y_{3JADE}$ and energy-energy correlation.

\section*{The future}

The program is not yet complete. The main problem is the reliance on Giele-Glover code. We are gradually moving away from this, but still need them at present because we have not yet implemented the Leiden matrix elements for Z exchange. The eventual aim is that they will only be a cross-check of the Leiden code. The eventual aim is to also work for arbitrary currents including W exchange and polarized leptons, but that's still a way o yet. These are the main changes that can be anticipated. Smaller cosmetic changes can also be expected once work gets more seriously under way on the DIS program, to make the different programs as similar as possible.

\section*{Modification Log}

\begin{description}
    \item[0.0] 18/01/96 First prerelease version put on web page.
    \item[0.1] 26/02/96 Bug fix in event orientation.
    \item[0.2] 15/03/96 Removed need for Lorentz boosts.
\end{description}

\end{document}


