%\documentclass[11pt]{cernrep}
%\usepackage{graphicx}
%\usepackage{floatflt}  
%\usepackage{fancybox}  %allows for shadow boxes etc.
%\usepackage{color}     %allows use of colors (RBG only, I think)
%\usepackage{multirow}
%\usepackage{myext}
%\usepackage{here}
%
% ========================= RDF DEFINITIONS =========================
\newcommand*{\eg}{{\it e.g.},\ }
\newcommand*{\ie}{{\it i.e.},\ }
\newcommand*{\etal}{{\it et al.\ }}
\newcommand*{\gev}{{\rm\,GeV}}
\newcommand*{\gevc}{{\rm\,GeV/c}}
\newcommand*{\mevc}{{\rm\,MeV/c}}
\newcommand*{\tev}{{\rm\,TeV}}
\newcommand*{\pt}{$p_T$}
\newcommand*{\bigpt}{$P_T$}
\newcommand*{\UE}{``underlying event"}
\newcommand*{\BBR}{``beam-beam remnants"}
\newcommand*{\MB}{Min-Bias}
\newcommand*{\pthard}{$p_T({\rm hard})$}
\newcommand*{\ptcut}{$p_T\!>\!0.5\,{\rm GeV/c}$}
\newcommand*{\etacut}{$|\eta|\!<\!1$}
\newcommand*{\hardcut}{$p_T\!({\rm hard})>3{\rm\,GeV/c}$}
\newcommand*{\hardzero}{$p_T\!({\rm hard})>0{\rm\,GeV/c}$}
\newcommand*{\ptchj}{$P_T\!({\rm chgjet}\#1)$}
\newcommand*{\aveN}{$\langle\! N_{\rm chg}\!\rangle$}
\newcommand*{\avePT}{$\langle P_T{\rm sum}\!\rangle$}
\newcommand*{\etaphi}{$\eta$-$\phi$}
%======================END RDF DEFINITIONS ==========================
%
\def\bea{\begin{eqnarray}}
\def\eea{\end{eqnarray}}
\def\bes{\begin{eqnarray*}}
\def\ees{\end{eqnarray*}}
\newcommand{\MEG}{MEG}
\def\q2{Q^2}              
\def\qt{Q_T}
\newcommand{\SHG}{SHG}

%\begin{document}
%\bibliographystyle{unsrt}
\newenvironment{myitemize}[2]
{ \begin{small}
  {\large\tt\bf #1} { \sc #2 } \vspace{.4cm} % CsB \vspace{-.2cm}
  \begin{itemize}
  \setlength{\itemsep}{2pt} % CsB \setlength{\itemsep}{-4pt}
  \setlength{\parsep}{0pt}
  \setlength{\topsep}{0pt}
  \setlength{\partopsep}{0pt}
}
{\end{itemize}\end{small}}
\newcommand{\myitem}[3]{ \item {\footnotesize\tt #1}~{\bf\tt #2}: {\it
    #3} }

\section{Monte Carlo\protect\footnote{Section coordinators:
I. Hinchliffe, J. Huston}$^{,}$~\protect\footnote{Contributing authors: 
C. Bal\'azs, E.\ Boos, M.\ Dobbs, W.\ Giele, 
I.\ Hinchliffe, Rick Field, J.\ Huston, V.\ Ilyin, 
J.\ Kanzaki, B. Kersevan, K.\ Kato, Y.\ Kurihara, 
L.\ L\"onnblad, K.\ Mazumudar, M.\ Mangano, S.\ Mrenna, 
F.\ Paige, I. Puljak, E.\ Richter-Was,
M.\ Seymour, T.\ Sj\"ostrand,  M. T\"onnesmann. B.\ Webber, D.\ Zeppenfeld}}
\label{sec:mc,qcdsm}

\subsection{Introduction}

The Monte-Carlo intergroup focused on four main issues: the problem of
interfacing partonic event generators to showering Monte-Carlos, an
implementation using this interface to calculate backgrounds which are poorly
simulated by the showering Monte-Carlos alone, a comparison of the {\tt 
HERWIG} and {\tt PYTHIA} parton  shower models with the predictions of a soft
gluon resummation program (ResBos),  and studies of the underlying events at
hadron colliders and how well they are modeled by the Monte-Carlo generators.





Section~\ref{interface} discusses a strategy whereby generic Fortran
common blocks are presented for use by High Energy
Physics event generators for the transfer of event configurations from
parton level generators to showering and hadronization event
generators.

Section~\ref{acer} discusses the {\bf AcerMC} Monte Carlo Event
Generator which is dedicated to the generation
of the Standard Model background processes at $pp$ LHC collisions. 
The program itself provides a library of the massive matrix elements 
and phase space modules for the generation of 
selected processes: $q \bar q \to W(\to \ell \nu) b \bar b$, $q \bar q \to W(\to \ell \nu) t \bar t$,
$gg, q \bar q \to Z/\gamma^*(\to \ell \ell) b \bar b$,
$gg, q \bar q \to Z/\gamma^*(\to \ell \ell, \nu \nu, b \bar b) t \bar t$, 
the QCD $gg, q \bar q \to t \bar t b \bar b$ and EW  $gg \to (Z/W/\gamma \to)
t \bar t b \bar b$. The hard process event, 
generated with these modules, can be completed by the initial and final 
state radiation, hadronization and decays, simulated with either the
{\tt PYTHIA} or {\tt HERWIG} Monte Carlo Event Generators. Interfaces to both
these generators are provided in the distribution version. The {\bf AcerMC} also uses several other external libraries: {\tt CERNLIB}, {\tt HELAS}, {\tt VEGAS}.


In Section~\ref{Huston}, predictions for the $p_T$ distribution for a Higgs particle were generated using two approaches: (1) a soft-gluon resummation technique, using the program {\tt ResBos} and (2) a parton shower technique using the Monte Carlo programs  {\tt HERWIG} and {\tt PYTHIA}. An understanding of the 
kinematics of the Higgs boson, or of any other Standard Model or non-
Standard Model particle, and the characteristics of any jets associated with its production, is of great interest for physics at the Tevatron, the LHC or any future hadron colliders. The transverse momentum distribution of the Higgs boson depends primarily on the details  of the soft gluon emission from the initial state partons. The effects of these soft gluon emissions can be described either by a
resummation calculation or by a Monte Carlo parton shower formalism. Comparisons between the two techniques for several different Higgs masses and for several center-of-mass energies show relatively good agreement between the {\tt ResBos} predictions and those of {\tt HERWIG} and recent versions of {\tt PYTHIA}. 




In Section~\ref{field}, the behavior of the \UE\ in hard scattering proton-antiproton collisions at $1.8\tev$ is studied and 
compared with the QCD Monte-Carlo models. The ``hard scattering" component consists of the outgoing two ``jets" plus initial and final-state radiation. The \UE\ is everything 
except the two outgoing hard scattered ``jets" and consists of the \BBR\ plus possible contributions 
from the ``hard scattering" arising from initial and final-state radiation.  In addition multiple parton 
scattering might contribute to the \UE.  The data indicate that
neither {\tt ISAJET} or {\tt HERWIG} produce 
enough charged particles (with \ptcut) from the ``beam-beam remnant" component and that 
{\tt ISAJET} produces too many charged particles from initial-state radiation.  The ``tuning" of {\tt ISAJET} and 
{\tt PYTHIA} to fit the \UE\ is explored.






In addition to these activities, discussions were held with the
``standard model'' and ``beyond the standard model'' groups in order
to assess their needs. Many of these needs are expressed in terms of
the desire to have new processes added to the showering Monte-Carlo
event generators. The interface discussed in Section \ref{interface}
should make this easier for users to implement the processes themselves, obviating
the need to have the Monte-Carlo authors do it. Among the important
items mentioned are:



\begin{itemize}
\item The production of $WZ$ final states should include the
  contributions from virtual $W$ and $Z$ as these are vital for SUSY
  searches at the Tevatron
\item  Resummation of the $log(m_H/m_b)$  and $log(m_H/p_t)$  terms in
      the $gg\to H b \overline{b}$ process.
\item  Radion phenomenology
\item A better understanding of the $gg\to H$ process and in
  particular the effect of a jet veto
\item Correct $\tau$ polarization in Higgs decay.
\end{itemize}


\subsection{Generic User Process Interface for Event Generators}
\label{interface}


Modularization of High Energy Particle Physics event generation is
becoming increasingly useful as the complexity of Monte Carlo programs
grows. To accommodate this trend, several authors of popular Monte
Carlo and matrix element 
programs attending the {\it Physics at TeV Colliders Workshop}
in Les Houches, 2001 have agreed on a generic format for the transfer
of parton level event configurations from matrix element event generators 
(\MEG) to showering and hadronization event generators (\SHG).

\begin{center} \color{blue}
\begin{tabular}{rcl}
  CompHEP\cite{Pukhov:1999gg}     & &        \\
  Grace\cite{Ishikawa:1993qr}     & & {\tt HERWIG}\cite{Corcella:2001bw} \\
  MadGraph\cite{Stelzer:1994ta}   & \hspace{1cm} $\Rightarrow$
                                  \hspace{1cm} & \tt ISAJET\cite{Baer:1999sp} \\
%%%       SusyGen\cite{Katsanevas:1998fb} 
  VecBos\cite{Berends:1991ax}     & & {\tt PYTHIA}\cite{Sjostrand:2000wi} \\
  WbbGen\cite{Caravaglios:1999yr} & & \ldots \\
  \ldots    & &        \\
\end{tabular}
\color{black} \end{center}

Events generated this way are customarily called user (or
user-defined) processes, to distinguish them from the internal
processes that come with the \SHG. Specific solutions are already in
use, including an interface of WbbGen with
{\tt HERWIG}~\cite{Mangano:2001xp} and an interface of CompHEP with
{\tt PYTHIA}~\cite{Belyaev:2000wn}---that experience is exploited here.


%% I know of atleast two others:   VecBos --> {\tt HERWIG} (see joey's mail)
%%                             and Grace --> {\tt PYTHIA}  (see spires)

Since the specification of the user process interface in May 2001 at
Les Houches, the interface has been (or is being) implemented in a
number of \MEG\ and \SHG\ programs.  An implementation has been
included in {\tt PYTHIA}~6.2, described (with an example) in
Ref.~\cite{Sjostrand:2001yu}.  A {\tt HERWIG} implementation is in progress
and will appear in Version~6.5.  \MEG\ implementations exist
for the Madison Collection of User Processes (MADCUP)~\cite{Cranmer:2002kc},
ALPGEN~\cite{Caravaglios:1999yr,Mangano:2002ap},
and CompHEP~\cite{Ilyin:2002vi} (to be available publicly soon). An
implementation in MadGraph is in preparation.  Other \MEG\
implementations include Refs.~\cite{Kersevan:2002dd,Dobbs:2001dq}.

The user process interface discussed here is not intended as a
replacement for \mbox{ {\tt HEPEVT}}~\cite{Sjostrand:1989tj}, 
which is the standard Fortran
common block for interface between generators and analysis/detector
simulation. The user process common blocks address the communication
between two event generators only, a \MEG\ one and a \SHG\ one, and not the 
communication of event generators with the outside world. 



In the course of a normal event generation run, this communication occurs 
at two stages: (1) at initialization, to establish the basic parameters 
of the run as a whole and (2) for each new event that is to be 
transferred from the \MEG\ to the \SHG. Each of these two stages here
corresponds to its own Fortran common block.\footnote{ An interface in 
C++ has been developed in Ref.~\cite{Dobbs:2001ck} and contains similar 
information content as that discussed here.} These common blocks are 
described in detail in the next two sections, followed by some examples.

One can also foresee that each stage will be associated with its own
subroutine, called from the \SHG, where information is put in the
respective common block, based on output from the \MEG. The details of
these subroutines are likely to be specific to a given \MEG\ and may
also be specific to a given \SHG.  The subroutine names {\tt UPINIT} and
{\tt UPEVNT} (each with no arguments) were chosen for the
{\tt PYTHIA}~6.2 implementation.  They are intended to
be generic (the usual {\tt PY} prefixes are omitted), and only dummy
versions are packaged with the program.  It is recommended that other
\SHG\ authors use the same dummy routine names (with zero arguments)
such that for simple cases which do not require intervention 
`by hand', \MEG\ authors will be able to interface several \SHG s with a
single set of routines. Example routines are presented in the {\tt PYTHIA}
documentation~\cite{Sjostrand:2001yu}.

In general, a user process run may consist of a number of subprocesses,
each denoted by a unique integer identifier. If the user
wishes to have the \SHG\ unweight events using acceptance-rejection and/or
mix together events from different processes, then the user process
author will need to supply a subroutine that is able to return an event of
the requested subprocess type to the \SHG\ on demand. 
The author may choose to organize the subroutine to
generate the event `on the fly', or to read the event from a file
stream (with a separate file stream for each subprocess).
The \SHG\ will also need information about the subprocess cross section
and/or maximum event weight to select which process is generated next
and for acceptance-rejection. This information will need to be
known from the onset (and could, for example, be determined in advance
from an initialization run). Alternatively, the user may already have 
a proper mixture of subprocesses from the \MEG\ and only wish the \SHG\
to process events in the order they are handed in. We therefore allow for 
several different event weight models.


If extra information is needed for a specific user implementation,
then a implementation-specific common block should be
created. The meaning of the user process common block
elements should not be overloaded, as this would defeat the generic
purpose.



The descriptions in this paper are intended for event generator
authors and may appear complex---most of the details will be
transparent to the casual user.

\subsubsection{`User Process' Run Information}

        The run common block contains information which pertains to a
        collection of events. \\
        In general this information is process dependent and it is
        impossible to include everything in a generic common
        block. Instead only the most general information is included
        here, and it is expected that users will have to intervene
        `by hand' for many cases (i.e.\ a user may need to specify
        which cutoffs are used to avoid singularities, which jet 
        clustering algorithm has been used to recombine partons in a
        next-to-leading-order calculation, the effective parton masses, 
        \ldots).



\begin{verbatim}
      integer     MAXPUP
      parameter ( MAXPUP=100 )
      integer     IDBMUP, PDFGUP, PDFSUP, IDWTUP, NPRUP, LPRUP
      double precision EBMUP, XSECUP, XERRUP, XMAXUP
      common /HEPRUP/ IDBMUP(2), EBMUP(2), PDFGUP(2), PDFSUP(2),
     +                IDWTUP, NPRUP, XSECUP(MAXPUP), XERRUP(MAXPUP), 
     +                XMAXUP(MAXPUP), LPRUP(MAXPUP)
\end{verbatim}
%      CHARACTER*28 CHPRUP
%      COMMON/HEPCUP/CHPRUP(MAXPUP)

\vspace{.5cm}

\begin{myitemize}{HEPRUP}{ `User Process' Run Common Block }

   \myitem{parameter} { MAXPUP=100 }
        {maximum number of different processes to be interfaced at one
        time} 

\hspace{-.2cm} {\bf \underline{ Beam Information } } \\
        Beam particle 1 (2) is defined as traveling along the +Z (--Z)
        direction.
   \myitem{integer} { IDBMUP(2) }
        {ID of beam particle 1 and 2 according to the Particle Data
        Group convention~\cite{Garren:2000st}}
   \myitem{double} { EBMUP(2) }
        {energy in GeV of beam particles 1 and 2}
   \myitem{integer} { PDFGUP(2) }
        { the author group for beam 1 and 2, according to the Cernlib
        PDFlib~\cite{Plothow-Besch:1993qj} specification }
   \myitem{integer} { PDFSUP(2) }
        { the PDF set ID for beam 1 and 2, according to the Cernlib
        PDFlib specification } \\
    For $e^+e^-$ or when the \SHG\ defaults are
        to be used, PDFGUP=--1, PDFSUP=--1 should be specified. \\
    The PDFlib enumeration of PDFs is sometimes out of date, but it is the
        only unique integer labels for PDFs available. In the case
        where a PDF not included in PDFlib is being used, this
        information will have to be passed `by hand'.

\hspace{-.2cm} {\bf \underline{ Process Information } }

   \myitem{integer} { IDWTUP }
        {master switch dictating how the event weights
        (XWGTUP) are interpreted} \\
        The user is expected to pick the most appropriate event weight 
        model for a run,
        given the \MEG\ input at hand and the desired output. Normally 
        the \SHG\ should be able to handle all of the models.\\   
        A summary of the IDWTUP switches is presented in
        Table~\ref{t_IDWTUP}.
        \begin{itemize}
        \setlength{\itemsep}{-4pt} \setlength{\parsep}{0pt}
        \setlength{\topsep}{0pt} \setlength{\partopsep}{0pt} 

        \item[+1] Events are weighted on input and the \SHG\ is asked
        to produce events with weight +1 as output. XSECUP and XERRUP
        need not be provided, but are calculated by the \SHG. XWGTUP
        is a dimensional quantity, in pb, with a mean value converging
        to the cross section of the process.  The \SHG\ selects the
        next subprocess type to be generated, based on the relative
        size of the XMAXUP(i) values. The user-supplied interface routine
        must return an event of the requested type on demand from the
        \SHG, and the maximum weight XMAXUP (or a reasonable
        approximation to it) must be known from the onset.
        A given event is accepted with a probability
        XWGTUP/XMAXUP(i). In case of rejection, a new event type and a
        new event are selected.  If XMAXUP(i) is chosen too low, such
        that XWGTUP violates the XMAXUP(i), the \SHG\ will issue a
        warning and update XMAXUP(i) with the new maximum weight
        value.
        If events of some types are already available unweighted, then
        a correct mixing of these processes is ensured by putting
        XWGTUP = XMAXUP(i). In this option also the internal \SHG\
        processes are available, and can be mixed with the external
        ones.  All weights are positive definite. $k$-factors may be
        included on an event by event basis by the user process 
        interface by re-scaling the XWGTUP for each event.

        \item[--1] Same as above (IDWTUP=+1), but the event weights
        may be either positive or negative on input, and the \SHG\
        will produce events with weight +1 or --1 as 
        output.\footnote{Negative-weight events may occur e.g.\ in 
        next-to-leading-order calculations. They should cancel against
        positive-weight events in physical distributions. The details
        of this cancellation are rather subtle when considered in
        the context of showers and hadronization, however, and a
        proper treatment would require more information than discussed 
        here. The negative-weight options should therefore be used with 
        some caution, and the negative-weight events should be a 
        reasonably small fraction of the total event sample.}
        A given event would be accepted with a probability
        $|$XWGTUP$|$/$|$XMAXUP(i)$|$ and assigned weight {\tt
        sign}(1,XWGTUP), where the {\tt sign} function transfers the
        sign of XWGTUP onto 1. A
        physics process with alternating cross section sign must
        be split in two IDPRUP types,\footnote{The motivation for this
        requirement is best understood with a simple example: imagine
        two subprocesses with the same cross section. The first
        process includes events with both positive and negative event
        weights such that two events out of three have weight +1 and the
        third --1. All events from the second process have positive
        weight +1. In this scenario these two processes should be
        `mixed' with proportions 3:1 to account for the cancellations
        that occur for the first process. The proportions for the
        mixing are communicated to the \SHG\ by supplying the positive
        and negative contributions to the cross section separately.
        }
        based on the sign of
        the cross section, such that all events of a particular IDPRUP
        have the same event weight sign.  Also the XMAXUP(i) values must
        be available for these two IDPRUP types separately, so that
        $|$XMAXUP(i)$|$ gives the relative mixing of event types, with
        event acceptance based on $|$XWGTUP$|$/$|$XMAXUP(i)$|$.

        \item[+2] Events are weighted on input and the \SHG\ is asked
        to produce events with weight +1 as output. The \SHG\ selects
        the next subprocess type to be generated, based on the
        relative size of the XSECUP(i) values.  The user-supplied
        interface routine must return an event of the requested type
        on demand from the \SHG. The cross sections XSECUP(i) must be
        known from the onset.
%though it may be poorly known at the beginning and gradually 
%updated by the user process
%       This would imply the presence of a separate file of events,
%        alternatively a separate generation algorithm, for each
%        process type.  
        A given event is accepted with a probability
        XWGTUP/XMAXUP(i).  In case of rejection, a new event of
        the same type would be requested.  In this scenario only the
        ratio XWGTUP/XMAXUP(i) is of significance. If events of
        some types are already available unweighted, then a correct
        mixing of these processes is ensured by putting XWGTUP =
        XMAXUP(i). A $k$-factor can be applied
        to each process by re-scaling the respective XSECUP(i)
        value at the beginning of the run, but cannot be given
        individually for each event. In this option also the internal
        \SHG\ processes are available, and can be mixed with the user
        processes.

        \item[--2] Same as above (IDWTUP=+2), but the event weights
        may be either positive or negative on input, and the \SHG\
        will produce events with weight +1 or --1 as output.  A
        physics process with alternating cross section sign must
        therefore be split in two IDPRUP types, based on the sign of
        the cross section, such that all events of a particular IDPRUP
        have the same event weight sign.  Also the XSECUP(i) and
        XMAXUP(i) values must be available for these two IDPRUP types
        separately, so that $|$XSECUP(i)$|$ gives the relative mixing
        of event types, with event acceptance based on
        $|$XWGTUP$|$/$|$XMAXUP(i)$|$ and the total cross section of
        the two IDPRUP types combined given by XSECUP(i)+XSECUP(j).

        \item[+3] Events are unweighted on input such that all events
        come with unit weight XWGTUP=+1. The \SHG\ will only ask for
        the next event.  If any mixing or unweighting is desired, it
        will have to be performed by the user process interface.  The 
        \SHG\ will not reject any events (unless it encounters other 
        kinds of problems).  If a $k$-factor is desired, it is the
        responsibility of the user process interface. When events are  
        read sequentially from an already existing file, this would    
        imply one common $k$-factor for all processes. In this option it 
        is not possible to mix with internal \SHG\ processes.


        \item[--3] Same as above (IDWTUP=+3), but the event weights
        may be either +1 or --1 on input. A single process identifier
        (IDPRUP) may include events with both positive and negative
        event weights.

        \item[+4] Same as (IDWTUP=+3), but events are weighted
        on input and the average of the event weights (XWGTUP) is
        the cross section in pb. When histogramming results on analyzed
        events, these weights would have to be used.
        The \SHG\ will only ask for the next event
        and will not perform any mixing or unweighting. Neither XSECUP
        nor XMAXUP needs to be known or supplied. In this
        option it is not possible to mix with internal \SHG\
        processes.

        \item[--4] Same as (IDWTUP=+4), but event weights may be either
        positive or negative on input and the average of the event
        weights (XWGTUP) is the cross section. A single process identifier
        (IDPRUP) may include events with both positive and negative
        event weights.

        \end{itemize} 

\begin{table}
\begin{center}
\begin{tabular}{c c c c c} \hline
        & event selection& control of        & XWGTUP &        \\
 IDWTUP & according to   & mixing/unweighting& input & output \\ \hline
%        
+1      & XMAXUP         & \SHG\        & + weighted     & +1 \\
--1     & XMAXUP         & \SHG\        & $\pm$ weighted & $\pm$1 \\
+2      & XSECUP         & \SHG\        & + weighted     & +1 \\
--2     & XSECUP         & \SHG\        & $\pm$ weighted & $\pm$1 \\
+3      & ---            & user interface & +1             & +1 \\
--3     & ---            & user interface & $\pm$1         & $\pm$1 \\
+4      & ---            & user interface & + weighted     & + weighted \\
--4     & ---            & user interface & $\pm$ weighted & $\pm$ weighted \\
\hline
\end{tabular}
\caption{\label{t_IDWTUP} Summary of the options available for the
master weight switch IDWTUP.}
\end{center}
\end{table}

   \myitem{integer} { NPRUP }
        { the number of different user subprocesses } \\
        i.e.\ LPRUP and other arrays will have NPRUP entries, LPRUP(1:NPRUP)

   \myitem{double} { XSECUP(J) }
        { the cross section for process J in pb } \\
        This entry is mandatory for IDWTUP=$\pm$2.

   \myitem{double} { XERRUP(J) }
        { the statistical error associated with the cross section 
        of process J in pb }\\
        It is not expected that this information will be used by the
        \SHG, except perhaps for printouts.

   \myitem{double} { XMAXUP(J) }{ the maximum XWGTUP for process J } \\
        For the case of weighted events (IDWTUP=$\pm$1,$\pm$2), this entry is
        mandatory---though it need not be specified to a high degree
        of accuracy. If too small a number is specified, the \SHG\
        will issue a warning and increase XMAXUP(J).
        This entry has no meaning for IDWTUP=$\pm$3,$\pm$4. 

   \myitem{integer} { LPRUP(J) } {a listing of all user process IDs
        that can appear in IDPRUP of HEPEUP for this run} \\
        When
        communicating between the user process and \SHG, the LPRUP
        code will be used. Example: if LPRUP(1)=1022, then the \SHG\
        will ask for an event of type 1022, not 1.

%   \myitem{character*28} { CHPRUP(J) }
%       { character string of process name (optional) }
%       This entry may be used by the \SHG s for print statements.
%       This entry is contained in its own common block, HEPCUP.

\end{myitemize}

\begin{center}{\rule[1mm]{5in}{1mm}}\end{center} %%%%%line

\subsubsection{`User Process' Event Information}

\begin{verbatim}
      integer MAXNUP
      parameter ( MAXNUP=500 )
      integer NUP, IDPRUP, IDUP, ISTUP, MOTHUP, ICOLUP
      double precision XWGTUP, SCALUP, AQEDUP, AQCDUP, 
     +                 PUP, VTIMUP, SPINUP
      common /HEPEUP/ NUP, IDPRUP, XWGTUP, SCALUP, AQEDUP, AQCDUP,
     +                IDUP(MAXNUP), ISTUP(MAXNUP), MOTHUP(2,MAXNUP),
     +                ICOLUP(2,MAXNUP), PUP(5,MAXNUP), VTIMUP(MAXNUP),
     +                SPINUP(MAXNUP)
\end{verbatim}

\vspace{.5cm}

\begin{myitemize}{HEPEUP}{ `User Process' Event Common Block }

  \myitem{parameter} { MAXNUP=500 }
        { maximum number of particle entries }
%         Parton level processes will typically have only a few
%         entries. A relatively large number is specified to also allow for
%         the possibility of interfacing at the stage between the parton
%         shower and hadronization steps of event generation.

  \myitem{integer} { NUP } 
        {number of particle entries in this event} \\
        An event with NUP=0 denotes the case where the user process is unable
        to provide an event of the type requested by the \SHG\
        (i.e.\ if the user process
        is providing events to the \SHG\ by reading them sequentially
        from a file and the end of the file is reached).

  \myitem{integer} { IDPRUP } 
        {ID of the process for this event} \\
        The process ID's are not intended to be generic. The entry is
        a hook which the event generators can use to
        translate into their own scheme, or use in print statements
        (e.g.\ so that cross section information can be shown per
        process). \\
        When IDWTUP$=\pm1,\pm2$ the next process to be generated 
        is selected by the \SHG, and so IDPRUP is set by the \SHG. For 
        IDWTUP$=\pm3,\pm4$ the process is selected by the \MEG, and
        IDPRUP is set by the \MEG.

  \myitem{double} { XWGTUP }
        { event weight } \\
        \underline{weighted events}:
        if the user process supplies
        weighted events and the \SHG\ is asked to produce 
        unweighted events, this number
        will be compared against XMAXUP in the run common block HEPRUP
        for acceptance-rejection. \\
        \underline{unweighted events}:
        if the user process supplies events which have already been
        unweighted, this number should be set to +1 (-1 for negative
        weight events in e.g. a NLO calculation).

        The precise definition of XWGTUP depends on the master weight
        switch IDWTUP in the run common block. More information is
        given there.

  \myitem{double} { SCALUP } 
        {scale of the event in GeV, as used for calculation of  PDFs} \\
        If the scale has not been defined, this should be denoted
        by setting the scale to --1.

  \myitem{double} { AQEDUP } 
        {the QED coupling $\alpha_\mathrm{QED}$ used for this event 
        (e.g.\ $\frac{1}{128}$) }

  \myitem{double} { AQCDUP } 
        {the QCD coupling $\alpha_\mathrm{QCD}$ used for this event }

        When $\alpha_\mathrm{QED}$ and/or $\alpha_\mathrm{QCD}$ is not
        relevant for the process, or in the case where the user
        process prefers to let the \SHG\ use its defaults,
        AQEDUP=--1 and/or AQCDUP=--1 should be specified.
      
\hspace{-.2cm} {\bf \underline{ ID, Status, and Parent-Child History } }
  \myitem{integer} { IDUP(I) } 
        {particle ID according to Particle Data Group 
        convention~\cite{Garren:2000st}} \\
        undefined (and possibly non-physical) ``particles'' should be
        assigned IDUP=0 (i.e.\ the $WZ$ particle in the example given
        in the MOTHUP description below)

  \myitem{integer} { ISTUP(I) } {status code} 
        \begin{itemize}
        \setlength{\itemsep}{-4pt} \setlength{\parsep}{0pt}
        \setlength{\topsep}{0pt} \setlength{\partopsep}{0pt} 
        \item[--1] Incoming particle
        \item[+1] Outgoing final state particle
        \item[--2] Intermediate space-like propagator defining an 
        $x$ and $Q^2$ which should be preserved
        \item[+2] Intermediate resonance, Mass should be preserved
        \item[+3] Intermediate resonance, for documentation only\footnote{
                Treatment of ISTUP(I)=+3 entries may be generator 
                dependent (in particular see
                Ref.~\cite{Sjostrand:2001yu} 
                for the special treatment in {\tt PYTHIA}).}
        \item[--9] Incoming beam particles at time $t=-\infty$
        \end{itemize} 
        The recoil from a parton shower (including
        photon emission) needs to be absorbed by particles in the
        event. Without special instructions, this can alter the mass
        of intermediate particles. The ISTUP flag +2 allows the user
        process to specify which intermediate states should have their
        masses preserved, i.e.\ for $e^+e^- \rightarrow Z^0 h^0
        \rightarrow q \bar{q} b \bar{b}$, the $Z^0$ and $h^0$ 
        would be flagged with ISTUP=+2. \\
        The primary application of the ISTUP=--2 status code is deep
        inelastic scattering (a negative number is chosen for this
        status code because the propagator in some sense can be
        thought of as incoming). See the example below. \\
        The status code ISTUP=--9 specifying incoming beams is not
        needed in most cases because the beam particle energy and identity
        is contained in the HEPRUP run information common block. The
        primary application of this status code will be non-collinear
        beams and/or runs for which the beam energy varies event by
        event (note that using the --9 status code to vary the machine
        energy may produce problems if the \SHG\
        is asked to combine separate processes).  The
        use of ISTUP=--9 entries is optional, and is only necessary
        when the information in HEPRUP is insufficient. If entries
        with ISTUP=--9 are specified, this information will over-ride
        any information in HEPRUP.

  \myitem{integer} { MOTHUP(2,I) } 
        {index of first and last mother} \\
        For decays, particles will normally have only one mother.\\ In
        this case either MOTHUP(2,I)=0 or MOTHUP(2,I)=MOTHUP(1,I).
        Daughters of a $2\rightarrow n$ process have 2 mothers.  This
        scheme does not limit the number of mothers, but in practice
        there will likely never be more than 2 mothers per particle.

        The history (intermediate particles) will be used by the
        \SHG s to decipher which
        combinations of particles should have their masses fixed and
        which particle decays should be ``dressed'' by the parton
        shower. Example: for $q\bar{q}'\rightarrow W^-Zg\rightarrow l^-
        \nu l^+ l^- g$, intermediate ``particles'' $WZ$, $W$, and $Z$
        could be specified with ISTUP=+2. Here the $WZ$ ``particle''
        would have its own entry in the common block with IDUP=0.  The
        showering generator would preserve the invariant masses of
        these ``particles'' when absorbing the recoil of the parton
        shower.

        In a case like $e^+e^- \rightarrow \mu^+\mu^-\gamma$
        proceeding via a
        $\gamma^*/Z^0$, where the matrix element contains an
        interference term between initial and final-state emission,
        this ambiguity in the parent-child history of the $\gamma$ has
        to be resolved explicitly by the user process.

  \hspace{-.2cm} {\bf \underline{ Color Flow } } \\
        A specific choice of color flow for a particular event is often
        unphysical, due to interference effects. However,
        \SHG s require a specific color state from which to begin
        the shower---it is the responsibility of the user process to
        provide a sensible choice for the color flow of a particular event.

       
  \myitem{integer} { ICOLUP(1,I) }
        {integer tag for the color flow line passing through the color of
        the particle}
  \myitem{integer} { ICOLUP(2,I) }
        {integer tag for the color flow line passing through the
        anti-color of the particle} \\
        The tags can be viewed as numbering the different color lines in
        the $N_C\rightarrow \infty$ limit.
        The color/anti-color of a particle are defined with respect to
        the physical time order of the process so as to allow a unique
        definition of color flow also through intermediate particles.\\
        This scheme is chosen because it has the fewest ambiguities, and
        when used with the history information, it supports Baryon
        number violation (an example is given below). \\
        To avoid confusion it is recommended that integer tags larger than 
        MAXNUP (i.e.\ 500) are used. The actual value of the tag has no meaning
        beyond distinguishing the lines in a given process.
      
  \hspace{-.2cm} {\bf \underline{ Momentum and Position } }
  \myitem{double} { PUP(5,I) }
        { lab frame momentum $(P_x, P_y, P_z, E, M)$ of particle in GeV } \\
        The mass is the `generated mass' for this particle, 
        $M^2=E^2-|\vec{p}|^2$ (i.e.\ not
        necessarily equal to the on-shell mass). The mass may 
        be negative, which denotes negative $M^2$ (i.e.\
        $M=2$ implies $M^2=4$ whereas $M=-2$ implies $M^2=-4$). \\
        Both $E$ and $M$ are needed for numerical reasons, the user
        should explicitly calculate and provide each one.

  \myitem{double} { VTIMUP(I) }
        { invariant lifetime $c\tau$ (distance from production to decay) in
        mm} \\
        Combined with the directional information from the momentum,
        this is enough to determine vertex locations.  Note that this
        gives the distance of travel for the particle from birth to
        death, in this particular event, and not its distance from the
        origin.

\hspace{-.2cm} {\bf \underline{ Spin / Helicity} }
  \myitem{double} { SPINUP(I) }
        { cosine of the angle between the spin-vector of particle I
          and the 3-momentum of the decaying particle, 
        specified in the lab frame } \\
        This scheme is neither general nor complete, but is chosen as
        the best compromise.
% since external legs to the partonic
%        subprocess will often be massless with a given chirality.
        The main foreseen application is $\tau$'s with a specific
        helicity. Typically a   
        relativistic $\tau^-$ ($\tau^+$) from a $W^-$ ($W^+$) has 
        helicity --1 (+1) (though this might be changed by the boost
        to the lab frame), so SPINUP(I)= --1 (+1). The use of a floating
        point number allows for the extension to the non-relativistic
        case. 
        Unknown or unpolarized particles should be given SPINUP(I)=9. The
        lab frame is the frame in which the four-vectors are
        specified.
        % 9 is used instead of 999 to allow for packing.        

\end{myitemize}

%
% {\tt HERWIG} uses separations at the femto scale (i.e. for b-decays the
%       separation in space WRT to the b-vertex is the thing of
%       importance, but since the b-vertex can be displaced at the 1mm 
%       level, machine precision (even for double) is not enough)
%
%
% NFUP, KFUP, Q2UP (which served this purpose in {\tt PYTHIA}::PYUPPR) are not
% needed.

%\item Communication with Tauola/Photos is not addressed. (!)
%\item A generic file read/write routine could be supplied.
%\item Special treatment is normally needed for Deep Inelastic
%Scattering, wherein the recoil from the parton shower should not
%affect the incoming electron. The showering generators will be able to
%recognize these cases by the identity of the incoming particles.

%By allowing the incident quark to be both incoming and a daughter to
%the $\gamma^*$, the history will allow the showering event generator
%to know that the recoil from the parton shower should not alter 
%$\hat{s}$ and $\hat{t}$ of the interaction.
%\begin{center}{\rule[1mm]{5in}{1mm}}\end{center} %%%%%line

\subsection*{Example: hadronic $t \bar{t} $ production }
%
\parbox{1.4in}{\includegraphics[height=1in]{example_tt.eps}}
%\begin{center}
\parbox{4.5in}{{\tiny \begin{tabular}{ccccccc}
I & ISTUP(I) & IDUP(I)   & MOTHUP(1,I) & MOTHUP(2,I) 
                                        & ICOLUP(1,I) & ICOLUP(2,I) \\ \hline
1 &--1 &  21~($g$)               & 0   & 0         & 501 & 502     \\
2 &--1 &  21~($g$)               & 0   & 0         & 503 & 501     \\
%3 & +2 &  0\footnotesize{ (undefined)}
%                                 &  1 & 2        & 503 & 502     \\
3 & +2 & --6~($\bar{t}$)         &  1 & 2        &  0  & 502     \\
4 & +2 &  6~($t$)                &  1 & 2        & 503 &  0      \\
5 & +1 & --5~($\bar{b}$)         &  3 & 3        &  0  & 502     \\
6 & +1 & --24~($W^-$)            &  3 & 3        &  0  &  0      \\
7 & +1 &  5~($b$)                &  4 & 4        & 503 &  0      \\
8 & +1 &  24~($W^+$)             &  4 & 4        &  0  &  0      \\
\hline
\end{tabular}}}
%\end{center}
%

The $t$ and $\bar{t}$ are given ISTUP=+2, which informs the \SHG\
to preserve their invariant masses when
showering and hadronizing the event. An intermediate s-channel gluon
has been drawn in the diagram, but since this graph cannot be usefully
distinguished from the one with a t-channel top exchange,
an entry has not been included for it in the event record.

The definition of a line as `color' or `anti-color' depends on the
orientation of the graph. This ambiguity is resolved  by defining 
color and anti-color according to the physical time order. 
A quark will always have its color tag ICOLUP(1,I) filled,
but never its anti-color tag ICOLUP(2,I). The reverse is true for an
anti-quark, and a gluon will always have information in both
ICOLUP(1,I) and ICOLUP(2,I) tags.

Note the difference in the treatment by the parton shower of the
above example, and an identical final state, where the intermediate
particles are not specified:
\begin{center}{\tiny 
\begin{tabular}{ccccccc}
I & ISTUP(I) & IDUP(I)   & MOTHUP(1,I) & MOTHUP(2,I) 
                                        & ICOLUP(1,I) & ICOLUP(2,I) \\ \hline
1 &--1        &  21~($g$)       &  0 & 0        & 501 & 502     \\
2 &--1        &  21~($g$)       &  0 & 0        & 503 & 501     \\
3 & +1        & --5~($\bar{b}$) &  1 & 2        &  0  & 502     \\
4 & +1        & --24~($W^-$)    &  1 & 2        &  0  &  0      \\
5 & +1        &  5~($b$)        &  1 & 2        & 503 &  0      \\
6 & +1        &  24~($W^+$)     &  1 & 2        &  0  &  0      \\
\hline
\end{tabular}}
\end{center}
%
In this case the parton shower will evolve the $b,\bar{b}$ without
concern for the invariant mass of any pair of particles. 
Thus the parton shower may alter the invariant mass of the $Wb$
system (which may be undesirable if the $Wb$ was generated from a top
decay).

%The above example is ofcourse unphysical, in reality we only know that
%we started with an initial state $w^+e^-

\subsection*{Example: $gg \rightarrow gg$ }
%
\parbox{1.4in}{\includegraphics[height=1.2in]{example_gg.eps}}
\parbox{4.5in}{\tiny \begin{tabular}{ccccccc}
I & ISTUP(I) & IDUP(I) & MOTHUP(1,I) & MOTHUP(2,I) 
                                & ICOLUP(1,I) & ICOLUP(2,I) \\ \hline
1 &--1        &  21~($g$)    &  0 & 0        &   501 & 502   \\
2 &--1        &  21~($g$)    &  0 & 0        &   502 & 503   \\
3 & +1        &  21~($g$)    &  1 & 2        &   501 & 504   \\
4 & +1        &  21~($g$)    &  1 & 2        &   504 & 503   \\
\hline
\end{tabular}}

\subsection*{Example: Baryon number violation in decays}
%
\parbox{1.4in}{\includegraphics[height=1.2in]{example_baryon.eps}}
\parbox{4.5in}{\tiny \begin{tabular}{ccccccc}
I & ISTUP(I) & IDUP(I)   & MOTHUP(1,I) & MOTHUP(2,I) 
                                        & ICOLUP(1,I) & ICOLUP(2,I) \\ \hline
1 &--1  &  11~($e^-$)                            &  0 & 0    &  0   &   0   \\
2 &--1  & --11~($e^+$)                           &  0 & 0    &  0   &   0   \\
3 & +2  &  23~($Z^0$)                            &  1 & 2    &  0   &   0   \\
4 & +2  &--1000002~($\stackrel{\sim}{\bar{u}}$)  &  3 & 3    &  0   &  501  \\
5 & +2  &  1000002~($\stackrel{\sim}{u}$)        &  3 & 3    & 501  &   0   \\
6 & +1  &  1~($d$)                               &  4 & 4    & 502  &   0   \\
7 & +1  &  1~($d$)                               &  4 & 4    & 503  &   0   \\
8 & +1  & --1~($\bar{d}$)                        &  5 & 5    &  0   &  504  \\
9 & +1  & --1~($\bar{d}$)                       &  5 & 5    &  0   &  505  \\
\hline
\end{tabular}
}%

Three ``dangling'' color lines intersect at the vertex joining the
$\stackrel{\sim}{\bar{q}},q,q'$ 
(and $\stackrel{\sim}{q},\bar{q},\bar{q}'$), which
corresponds to a Baryon number source (sink) of +1 (-1), and will be
recognizable to the \SHG s.

\subsection*{Example: Baryon number violation in production}
%
\parbox{1.5in}{\includegraphics[height=1.2in]{example_baryon2.eps}}
\parbox{4.5in}{\tiny \begin{tabular}{ccccccc}
I & ISTUP(I) & IDUP(I)               & MOTHUP(1,I) & MOTHUP(2,I) 
                                        & ICOLUP(1,I) & ICOLUP(2,I) \\ \hline
1 &--1 &  1~($d$)                               & 0 & 0   & 501  & 0        \\
2 &--1 &  1~($d$)                               & 0 & 0   & 502  & 0        \\
3 & +2 & --1000002~($\stackrel{\sim}{\bar{u}}$) & 1 & 2   &  0   & 503      \\
4 & +1 & --2~($\bar{u}$)                        & 3 & 3   &  0   & 503      \\
5 & +1 &  1000022~($\stackrel{\sim}{\chi}^0$)   & 3 & 3   &  0   & 0        \\
\hline
\end{tabular}}

\subsection*{Example: deep inelastic scattering }
%
\label{example_dis.eps}
\parbox{1.5in}{\includegraphics[height=1.2in]{example_dis.eps}}
\parbox{4.5in}{\tiny \begin{tabular}{ccccccc}
I & ISTUP(I) & IDUP(I) & MOTHUP(1,I) & MOTHUP(2,I) 
                                & ICOLUP(1,I) & ICOLUP(2,I) \\ \hline
1 &--1        &  11~($e^-$)    &  0 &  0       &   0   &  0    \\
2 &--1        &  21~($g$)      &  0 &  0       &  503  & 501   \\
3 &--2        &  22~($\gamma$) &  1 &  0       &   0   &  0    \\
4 & +1        &  11~($e^-$)    &  1 &  0       &   0   &  0    \\
5 & +1        &  -4~($\bar{c}$)&  2 &  3       &   0   & 501   \\
6 & +1        &  4~($c$)       &  2 &  3       &  502  &  0    \\
7 & +1        &  21~($g$)      &  2 &  3       &  503  & 502   \\
\hline
\end{tabular}}

For DIS, the $x$ and $q^2$ of the $\gamma$ should not be altered by
the parton shower, so the $\gamma$ is given ISTUP=--2. We have not
specified the internal quark and gluon lines which will be dressed by
the parton shower, such that the partonic event configuration may be
drawn as follows, \\
\parbox{1.6in}{\includegraphics[height=1.5in]{example_dis2.eps}}
\parbox{4.5in}{If information about the quark and gluon propagators is desired (i.e\
for human readability), then those entries may be included with status
code ISTUP=+3.}

\subsection{The Monte Carlo Event Generator AcerMC }
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\label{acer}

\unboldmath

Despite the existence of a large repertoire of  processes implemented in
universal generators such as {\tt PYTHIA} or {\tt HERWIG}, a number of  Standard Model
background processes crucial for studying expected physics potential of the LHC
experiments are still missing.  For some of these processes, the matrix element
expressions are rather lengthy and/or complex, and to achieve a reasonable generation
efficiency, it is necessary to tailor the phase-space selection procedure to the dynamics
of the process.  The practical solution could  therefore be to produce a choice of dedicated
matrix-element-based generators with standardized interfaces (such as the one proposed in
section \ref{interface}) to the more general Monte Carlo programs such as {\tt HERWIG} or {\tt PYTHIA},  which are then used to complete the event generation.

The {\bf AcerMC} Monte Carlo Event Generator \cite{Kersevan:2002dd} follows up on this idea.
It is dedicated to the simulation of  Standard Model background processes in 
LHC collisions. The program itself provides a library of the massive matrix elements 
and phase space modules for the generation of a few selected $2 \to 4$ processes.
The hard process event, generated with these modules, can be completed by the addition of 
initial and final state radiation, hadronization and decays, simulated with either 
{\tt PYTHIA~6.2} \cite{Sjostrand:2001yu} or {\tt HERWIG 6.3} \cite{Corcella:2001bw} 

Interfaces of {\tt AcerMC 1.0} to both  {\tt PYTHIA~6.2} and {\tt HERWIG 6.3} are prepared following the standard proposed in Section
\ref{interface},
and are provided in the distribution version. The {\tt AcerMC 1.0} also uses
several other external libraries: {\tt CERNLIB}, {\tt HELAS}, {\tt VEGAS}.
The matrix element codes have been derived with the help of the {\tt MADGRAPH/HELAS} package.
The  typical  efficiency achieved for the generation of unweighted events 
is up to {\bf 30\%}, rather impressively high given the complicated topology of the 
implemented processes.

The very first version of this library was  interfaced to {\tt PYTHIA 6.1 } within the
standard of the so called  {\it external processes} in {\tt PYTHIA}
\cite{Kersevan:2001ab, Kersevan:2001cd, Kersevan:2001ef}.
Since then, after upgrading to the {\tt AcerMC 1.0} standard, the efficiencies have
significantly improved due to an additional optimization step in the phase space
generation.  Also in the new version, the interface standard has been changed from {\tt PYTHIA 6.1 } to {\tt
PYTHIA~6.2 } conventions, an interface to the {\tt HERWIG~6.3} generator was introduced, and the
native {\tt AcerMC 1.0} calculations of the $\alpha_{QED}$ and $\alpha_{QCD}$ couplings
were coded to allow for consistent benchmarking between results obtained with {\tt
PYTHIA~6.2 } and {\tt HERWIG~6.3}.  We also added the $q \bar q \to W(\to \ell \nu) t \bar t$ and
$gg, q \bar q \to Z/\gamma^*(\to \ell \ell, \nu \nu, b \bar b) t \bar t$ and electro-weak 
$gg \to (Z/W/\gamma \to) b \bar b t \bar t$  processes, which have been
implemented for the first time in the {\tt AcerMC 1.0} library presented here.

It is not necessarily the case that the the lowest order
 matrix element calculations  
for a given topology represents the total expected background of a given type. This is
particularly true concerning the heavy flavour content of an event.  The heavy flavour in a
given event might  occur in the hard process of a much simpler topology, as the effect
of including higher order QCD corrections via the shower mechanism. This is  the case, for example, 
for the presence of b-quarks  in  inclusive Z-boson or W-boson production. 
$Wb \bar b$ or $Zb \bar b$ final states can be calculated through higher order radiative corrections to inclusive W and Z production (through parton showering), or with the use of explicit $Wb \bar b$ and $Zb \bar b$ matrix elements. 
%Inclusive  production  has
%a total cross-section orders of magnitude higher than the  (matrix-element-based)
%$Wb \bar b$ or $Zb \bar b$ production.  
The matrix-element-based calculation
itself is a very good reference point to compare with parton shower approaches using
different fragmentation/hadronization models.  It also helps to  study matching
procedures between calculations at a fixed $\alpha_{QCD}$ order and parton shower
approaches. For exclusively hard topologies matrix-element-based calculations
usually represent a much better approximation than the parton shower ones.

The physics processes implemented in the {\tt AcerMC 1.0} library represent a set of important
Standard Model background processes.  These processes are all  key backgrounds for
discoveries in  channels characterized by the presence of  heavy flavour jets and/or
multiple isolated leptons \cite{ATLAS:1999ab}.  The  Higgs boson searches,  $t \bar t
H$, $ZH, WH$ (with $H \to b \bar b$),  $gg \to H$ with ($H\to ZZ^* \to 4 \ell$), and  $b\bar
b h/H/A$ (with $h/H/A \to \mu \mu $) are the most obvious examples of such channels.  We will
briefly discuss the physics motivations for the  processes of interest and the 
implementations that are available:

\boldmath {\bf $Z/\gamma^*(\to \ell \ell)  b \bar b$ production} \unboldmath 
has, over the last several years,  has been recognized as one of the most substantial backgrounds for  several Standard Model (SM) and Minimal Supersymmetric
Standard Model (MSSM) Higgs boson decay modes, as well as for the observability of  SUSY particles.  There is a rather wide spectrum of {\it regions of interest} for this
background.  In all cases, the leptonic $Z/\gamma^*$ decay is relevant; events with
di-lepton invariant mass below, at, or  above  the Z-boson mass could be of interest.  This process enters an analysis
either by the accompanying b-quarks being tagged as b-jets, or by the presence of leptons
from the b-quark semi-leptonic decays in these events. 

A good understanding of those backgrounds and the availability of a credible Monte Carlo
generator which allows the study of the  expected acceptances for different final states
topologies, is crucial. Despite a large effort expended at the time of the
Aachen Workshop \cite{Aachen1990}, such well established Monte Carlo generators were missing
for several years\footnote{
The matrix element for the $gg \to Z b \bar b \to b \bar b \ell \ell$ production has been
published in \cite{Kleiss1990} and, at the  time of the Aachen Workshop, implemented into the
{\tt EUROJET} Monte Carlo \cite{vanEijk:1994ab}.  Since that generator did not allow for the
possibility for having a fully generated hadronic event, with initial and final state
radiation and hadronization, it was
interfaced to {\tt PYTHIA 5.6} \cite{Sjostrand:1994yb} Monte Carlo for the analyses presented in \cite{DellaNegra1990}.  This program, however, is no longer
supported.
The same matrix element has been directly implemented into {\tt PYTHIA 5.7}
\cite{Sjostrand:1994yb}. However, with this implementation the algorithm for the phase space
generation has never worked credibly and thus it was removed from  
{\tt PYTHIA 6.1} \cite{Sjostrand:2001yu}.}.
Recently, the massless matrix elements for the $gg, q \bar q \to Z b \bar b$ processes have been
implemented in the general purpose Monte Carlo program {\tt MCFM} \cite{Campbell:2000bg}.  In
that implementation radiative corrections to this process have also been addressed. 
The massive matrix elements, with an interface to {\tt PYTHIA 6.1}, became 
available in \cite{Kersevan:2001cd}. The {\bf AcerMC} library discussed here
includes an even more efficient implementation of the algorithm presented
in \cite{Kersevan:2001cd}.
The same process is also implemented in the very recent  
version of {\tt HERWIG (6.3)} \cite{Corcella:2001bw}.

\boldmath {\bf $Z/\gamma^*(\to \ell \ell, \nu \nu, b \bar b)  t \bar t$ production}
 \unboldmath at the LHC is an irreducible background to a Higgs search in  association with a 
top-quark pair \cite{Gunion:1994ab}. With the  $Z/\gamma^*(\to b \bar b)$ decay, it is also
an irreducible resonant background to a Higgs search in the $t \bar t H$  production channel
in which the Higgs boson decays to a b-quark pair \cite{Richterwas:1998ab}. 

\boldmath {\bf $W (\to \ell \nu) b \bar b$ production} \unboldmath at the LHC is
recognized as a substantial irreducible 
background for both   Standard Model (SM) and Minimal Supersymmetric
Standard Model (MSSM) Higgs boson searches,  in the associated production mode
$WH$, followed by the decay $H \to b \bar b$.
The massive matrix element for the $q \bar q \to W g^*(\to b \bar b)$ process
 has been calculated \cite{Mangano:1993kq} and interfaced with the {\tt HERWIG 5.6} Monte Carlo \cite{Marchesini:1988cf,Knowles:1988vs,Catani:1991rr} several years ago.
A more recent implementation of the $Wbb$ + multi-jet final states 
is available from \cite{Mangano:2002ab}. Recently, the massless matrix element
has been implemented in the general purpose Monte Carlo program
{\tt MCFM} \cite{Campbell:2000bg}, where the radiative corrections
to this process are also addressed. Yet another implementation of the
$q \bar q \to W ( \to \ell \nu)g^* (\to b \bar b)$ massive matrix elements,
with an interface to {\tt PYTHIA 6.1} became available in \cite{Kersevan:2001ab}.
The {\bf AcerMC} library discussed here  includes a more efficient implementation of
the algorithm presented in \cite{Kersevan:2001ab}. 

\boldmath {\bf $W (\to \ell \nu) t \bar t$ production } \unboldmath 
at the LHC, has,  to our knowledge,  not been implemented
in any publicly available code so far. It is of interest
\footnote{ We thank M. A. Mangano for bringing
this process to our attention and for providing benchmark numbers for verifying 
the total cross-section} because it contributes 
an overwhelming background \cite{Mazzucato:2002ab} to the measurement of the 
Standard Model Higgs  self-couplings at LHC in the most promising 
channel $pp \to HH \to WWWW$.

\boldmath {\bf $t \bar t b \bar b$ production  } \unboldmath
at  the LHC is a dominant irreducible background for both Standard Model (SM) and Minimal
Super-symmetric Standard Model (MSSM) Higgs boson searches in  associated production, $ t
\bar t H$, followed by the decay $H \to b \bar b$.  The potential for the observability of
this channel has been carefully studied and documented in \cite{ATLAS:1999ab} and
\cite{Richterwas:1998ab}.  The proposed analysis requires identifying four b-jets,
reconstruction of both top-quarks in the hadronic and leptonic modes and the visibility of a
peak in the invariant mass distribution of the remaining b-jets.  The irreducible $t \bar
t b \bar b$ background contributes about 60-70\% of the total background from  $t \bar
t$ events ($t \bar t b \bar b$, $t \bar t b j$, $t \bar t j j$).
In the {\bf AcerMC} library, we have implemented both QCD and EW processes leading to the
$t \bar t b \bar b$ final state, namely $gg, q \bar q \to t \bar t b \bar b$ and
$gg \to (Z/W/\gamma \to) t \bar t b \bar b$. The contribution from  EW processes, 
never studied thus far, is surprisingly important in the mass range of the b-quark system 
around 120~GeV, see \cite{Kersevan:2002dd}.
It would seem that the analyses documented in \cite{ATLAS:1999ab}
and \cite{Richterwas:1998ab} might need revisiting.

These complete the list of the native {\bf AcerMC 1.0} processes implemented so far (see
Table~\ref{AcerMC:Tab1}).  Having all these different production processes implemented in
a consistent framework, which can also be directly used  for generating standard {\tt
PYTHIA} or {\tt HERWIG} processes, represents a very convenient environment for several
phenomenological studies related to LHC physics.

\begin{Tabhere}
\vspace{-0.4cm}
\newcommand{\lstrut}{{$\strut\atop\strut$}}
  \caption {\em Matrix-element-based processes implemented  in the{\bf AcerMC} library.
}\label{AcerMC:Tab1}
\vspace{2mm} 
\begin{center}
\begin{tabular}{|c|c|c|} \hline \hline
Process id & Process specification & Efficiency for generation \\ 
           &                       &  of  unweighted events \\ 
\hline\hline
 1      &  $g + g  \to   t \bar t b \bar b $   &  20.2 \%  \\
\hline
 2      &  $q + \bar q  \to   t \bar t b \bar b $   &  26.3 \%  \\
\hline
 3      &  $q + \bar q  \to W(\to \ell \nu)  b \bar b $   &  33.0 \%  \\
\hline
 4      &  $q + \bar q  \to W(\to \ell \nu)  t \bar t $   &  21.0 \%  \\
\hline
 5      &  $g + \bar g  \to Z(\to \ell \ell)  b \bar b $   &  39.0 \%  \\
\hline
 6      &  $q + \bar q  \to Z(\to \ell \ell)  b \bar b $   &   31.7 \% \\
\hline
 7      &  $g +  g  \to Z(\to \ell \ell, \nu \nu, b \bar b)  t \bar t $   &   28.2 \%  \\
\hline
 8      &  $q + \bar q  \to Z(\to \ell \ell, \nu \nu, b \bar b)  t \bar t $   &   34.6 \%  \\
\hline
 9      &  $g +  g  \to (Z/W/\gamma \to) b \bar b  t \bar t $   &   11.2 \%  \\
\hline \hline
\end{tabular}
\end{center}
\end{Tabhere}

\boldmath
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Monte Carlo algorithm} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\unboldmath

{\bf AcerMC 1.0} produces unweighted events with colour flow information using the
{\tt MADGRAPH/HELAS} \cite{Stelzer:1994ta} package, {\tt PDFLIB} \cite{Plothow-Besch:1993qj} and either the native
or {\tt PYTHIA 6.2}/{\tt HERWIG 6.3} coded running couplings $\alpha_s$ and $\alpha_{\rm
QED}$ (user's choice), for matrix element calculation and native (multi-channel based)
phase space generation procedures.  The generated events are then passed to either the {\tt
PYTHIA 6.2} or {\tt HERWIG 6.3} event generators, where the fragmentation and
hadronization procedures, as well as the initial and final state radiation, are added and
final unweighted events are produced.
\vspace{0.2cm}

\boldmath
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\noindent{\bf The Matrix Element Calculation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\unboldmath
\vspace{0.2cm}

The {\sf FORTRAN}-coded squared matrix elements of the processes were obtained by using the {\tt
MADGRAPH/HELAS} \cite{Stelzer:1994ta} package, taking properly into account the masses and
helicity contributions of the participating particles. The particle masses, charges and
coupling values that were passed to {\tt MADGRAPH} were taken from the interfaced libraries
({\tt PYTHIA/HERWIG}) to preserve the internal consistency of the event generation
procedure. In addition, the (constant) coupling values of $\alpha_s$ and $\alpha_{\rm
QED}$ were replaced with the appropriate running functions that were either taken from the
interfaced generators or provided by the {\bf AcerMC} code according to user settings.
In addition, a slightly modified {\tt MADGRAPH/HELAS} code was used for obtaining the colour
flow information of the implemented processes.
\vspace{0.2cm}

\boldmath
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\noindent{\bf The Four Fermion Phase Space Generation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\unboldmath
\vspace{0.2cm}

The four-fermion phase space corresponding to the processes discussed was modeled
using the importance sampling technique based on the procedures implemented in
the $e^+e^-$ event generators {\tt FERMISV} \cite{Hilgart:1993xu}, {\tt EXCALIBUR}
\cite{Berends:1995xn} and {\tt NEXTCALIBUR} \cite{Berends:2000fj}. For each implemented
process, a sequence of different kinematic diagrams ({\it channels}) modeling
the expected event topologies was constructed and the relative weights between the
contributions of each sampling channel was subsequently obtained by using a
multi-channel self-optimising approach \cite{Kleiss:1994qy}. Eventually, additional
smoothing of the phase space was obtained by using a modified {\tt VEGAS}
routine to improve the generation efficiency. 

The procedure of multi-channel importance sampling used in the event generation can
briefly be outlined as follows. An analytically integrable function
$g(\vec{\Phi})$, which aims to approximate the peaking behaviour of the differential
cross-section dependence on various kinematic quantities, is introduced into the
differential cross-section equation as:
\begin{equation}
d\sigma = s(\vec{\Phi})\, d\vec{\Phi} = \frac{s(\vec{\Phi})}{g(\vec{\Phi})} \cdot
g(\vec{\Phi})\, d\vec{\Phi} = w(\vec{\Phi})\,g(\vec{\Phi})\, d\vec{\Phi},
\end{equation}
where  $d\vec{\Phi}$ denotes the (four-)particle phase space and 
$s(\vec{\Phi})$ summarizes the matrix element, flux and structure functions,
all of which depend on the chosen phase space point.  The function $g(\vec{\Phi})$
is required to be unitary, i.e. a normalized probability density:
\begin{equation}
\int g(\vec{\Phi}) d\vec{\Phi} = 1.
\end{equation}
Since the peaking behaviour of $s(\vec{\Phi})$ can be very complex, due to the several possible
topologies introduced by a large number of contributing Feynman diagrams, 
the function $g(\vec{\Phi})$
is composed of a weighted sum of several channels $g_i(\vec{\Phi})$, each of which is adapted to a
certain event topology:
\begin{equation}
\int g(\vec{\Phi})  = \sum_i \alpha_i \cdot g_i(\vec{\Phi}). 
\end{equation}
The values of the relative weights $\alpha_i$ are determined from a multi-channel
self-optimization procedure  in order to minimize the variance of the 
weights, $w(\vec{\Phi})$ \cite{Kleiss:1994qy}. The phase space points are than sampled from
the function $g(\vec{\Phi})$, first by randomly choosing a channel $i$ according to the
relative frequencies $\alpha_i$, and then deriving the required four momenta from the
chosen $g_i(\vec{\Phi})$ using unitary algorithms. 

The modeling of the kinematic channels has heavily relied on the procedures developed in the {\tt
NEXTCALIBUR} program \cite{Berends:2000fj}; nevertheless, many additions and improvements
were made. The detailed description of the implementations of the four-momenta sampling in all
existing kinematic channels is omitted for the sake of brevity; an example of the
extended/added procedures used in {\bf AcerMC}, as given below, should serve as a representative
illustration,  For further details on the applied method and
unitary algorithms the reader is referred to the original papers
(e.g. \cite{Berends:2000fj,Hilgart:1993xu}).  

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\noindent {\bf Example: Breit-Wigner Function with s-dependent Width}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\vspace{0.2cm}

In some topologies of  processes involving $W^\pm$ or $Z^0$ bosons, a bias of the matrix
element towards large values in the high $s^*_{W/Z}$ region is evident. This, in turn, means
that a more accurate description of the tail of the $s^*_W$ distribution is needed. Consequently, the
Breit-Wigner sampling function was replaced by:
\begin{equation}
BW_s(s^*_W)=\frac{s^*_W}{(s^*_W - M_W^2)^2 + M_W^2 \Gamma_W^2},
\label{e:bwnew}
\end{equation}
which is proportional to the (more accurate) Breit-Wigner function with an $s^*_W$ dependent
width (W in the above formula denotes either a $W^\pm$ or a $Z^0$ boson). 

After some calculation the whole unitary procedure can thus be listed as follows:
\begin{itemize}
\item Introduce a new variable  $\eta = (s^*_W-M_W^2)/(M_W\;\Gamma_W)$. The integral of
the above function thus gives:
\begin{equation}
F(\eta)=  \{ \frac{M_W^2}{M_W \Gamma_W} \cdot {\rm atan}(\eta) \}  + 
\{ \frac{1}{2} \{ \log(\eta^2+1) \}, =
F_1(\eta) + F_2(\eta)
\label{e:bwint}
\end{equation}
\item Calculate the kinematic limits $\eta_{\rm min}$ and $\eta_{\rm max}$.
\item Calculate the {\it normalisation} factors $\Delta_1=F_1(\eta_{\rm
max})-F_1(\eta_{\rm min})$, $\Delta_2=F_2(\eta_{\rm max})-F_2(\eta_{\rm min})$
and $\Delta_s = \Delta_1 + \Delta_2$; the term $\Delta_2$ can actually be
negative and thus does not represent a proper normalisation.

\item Obtain a (pseudo-)random number $\rho_1$. 
\item If $\rho_1 \leq \Delta_2/\Delta_s$ then:
\begin{itemize}
\item Obtain a (pseudo-)random number $\rho_2$;
\item Construct $\eta$ as:
\begin{eqnarray*}
X & = & \Delta_2 \cdot \rho_2 + F_2(\eta_{\rm min}),\\
\eta & = & \sqrt(e^{2X}-1),
\end{eqnarray*}
which is the inverse of the (normalized) cumulant $(F_2(\eta)-F_2(\eta_{\rm
min}))/\Delta_2$.
\item Note that the condition $\rho_1 \leq \Delta_2/\Delta_s$ can be fulfilled
only if $\Delta_2 \geq 0$, which means that $\eta_{\rm max}$ is positive and
greater than $\eta_{\rm min}$.
\end{itemize}
\item Conversely, if $\rho_1 > \Delta_2/\Delta_s$ then:
\begin{itemize}
\item Obtain a (pseudo-)random number $\rho_2$;
\item Construct $\eta$ as:
\begin{eqnarray*}
X & = & \Delta_1 \cdot \rho_2 + F_1(\eta_{\rm min}),\\
\eta & = & \tan(\frac{M_W \Gamma_W}{M_W^2}\cdot X)
\end{eqnarray*}
which is the inverse of the (normalized) cumulant $(F_1(\eta)-F_1(\eta_{\rm
min}))/\Delta_1$.
\item If the obtained $\eta$ is less than zero then calculate the normalized
probability densities:
\begin{eqnarray*}
P_1 & = & \frac{1}{\Delta_1} \cdot \{\frac{M_W^2}{M_W \Gamma_W} \cdot
\frac{1}{1+\eta^2} \} \\
P_s & = &  \frac{1}{\Delta_s} \cdot \{\frac{M_W^2}{M_W \Gamma_W} \cdot
\frac{1}{1+\eta^2} + \frac{\eta}{1+\eta^2} \}
\end{eqnarray*}
\item Obtain a (pseudo-)random number $\rho_3$;
\item If $\rho_3 > P_s/P_1$ map $\eta \to -\eta$.
\item If the new $\eta$ falls outside the kinematic limits $[\eta_{\rm
min},\eta_{\rm max}]$ the event is rejected.
\item Note also that the last mapping can only occur if the original $\eta$ was
negative, since $P_s < P_1$ only in the region $\eta < 0$. 
\item Calculate the value of $s^*_W$ using the inverse of $\eta$ definition:
\begin{equation}
s^*_W = (M_W\;\Gamma_W) \cdot \eta + M_W^2
\end{equation}
\end{itemize}
The weight corresponding to the sampled value $\eta$ is exactly:
\begin{equation}
\Delta_s \cdot \frac{(s^*_W - M_W^2)^2 + M_W^2 \Gamma_W^2}{s^*_W},
\end{equation}
which is the (normalized) inverse of Equation \ref{e:bwnew} as requested.
\end{itemize}
\begin{figure}
\vspace{-0.8cm}
\begin{center}
\mbox{
     \includegraphics[width=7cm]{acph_2.eps} 
     \includegraphics[width=7cm]{acph_1.eps} 
}
\end{center}
\vspace{-0.8cm}
\caption{\em {\bf Left} Comparisons of the two functional terms contributing to 
${\rm BW}_s(\eta)$ given by Equation \ref{e:bwnew}. Note that the scaling factor $A$ is
chosen in view of making the contributions more transparent; it is much too small compared
to the real case of $W^\pm/Z^0$ bosons. {\bf Right} Comparison of the (normalized)
distributions of differential cross-section for the process $q \bar q \to W b \bar b$
(dashed) and sampling functions (solid line) with respect to the variables obtained by
importance sampling, as described in the text.\label{f:bwcomp}
}
\vspace{0.2cm}
\end{figure}

Using the above re-sampling procedure, the whole approach remains completely unitary,
i.e. no events are rejected when either there are no limits set on the value of $\eta$ or they are
symmetric, $|\eta_{\rm min}|=\eta_{\rm max}$. In the contrary case, a small fraction of
sampling values is rejected.

As it turns out in subsequent generator level studies, this generation procedure
provides a much better agreement with the  differential distributions than the
{\it usual} (width independent) Breit-Wigner; an example obtained for the  $q \bar q \to W b
\bar b$ process is shown in Figure~\ref{f:bwcomp}. The evident consequence is
that the unweighting efficiency is substantially improved due to the reduction
of the event weights in the high $s^*_W$ region.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\noindent {\bf Modified {\tt VEGAS} Algorithm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\vspace{0.2cm}

Using the  multi-channel approach previously described, the total generation (unweighting)
efficiency amounts to about $3-10\%$ depending on the complexity of the chosen process.
 In order to further improve the efficiency, a set of modified {\tt VEGAS}
\cite{Lepage:1978sw} routines was used as a (pseudo-)random number generator for
sampling the peaking quantities in each kinematic channel. After training
all the sampling grids (of dimensions 4-7, depending on the kinematic channel), the
generation efficiency increased to values of above $20\%$. The motivation for this approach
is, that for  unitary algorithms, only a very finite set of simple sampling functions is
available, since the functions have to have simple analytic integrals for which an inverse
function also exists. Consequently, the non-trivial kinematic distributions can not be
adequately described by simple functions  in the whole sampling domain (e.g. the
$\tau$ distribution, c.f. Figure~\ref{f:tauveg}) and some additional smoothing might be
necessary.
\begin{figure}
\vspace{-1.1cm}
\begin{center}
\mbox{
     \includegraphics[width=7cm]{ttbb_tau.eps}
     \includegraphics[width=7cm]{ttbb_tauveg.eps}
}
\end{center}
\vspace{-0.4cm}
\caption{\em
Comparison between the sampling distribution for the $\tau = \hat{s}/s \in [\tau_{\rm
min},1]$ variable before and after the application of modified {\tt VEGAS} \cite{Lepage:1978sw}
smoothing procedure (see text).\label{f:tauveg}}
\end{figure}
\vspace{0.5cm}


 In addition, the random number distributions should, due to the applied
importance sampling, have a reasonably flat behaviour to be approached by an adaptive
algorithm such as {\tt VEGAS}\footnote{At this point  a disadvantage of using the
adaptive algorithms of the {\tt VEGAS} type should be stressed, namely that these are
burdened with the need of training them on usually very large samples of events before
committing them to event generation.}. The principal modification of {\tt VEGAS},
besides adapting it to function as a (pseudo-)random number generator instead of
the usual {\it integrator}, was based on the discussions \cite{Ohl:1999qm,Jadach:1999vf} that in the case
of event generation, i.e.the unweighting of events to  weight one, reducing the maximal
value of event weights is in principle of higher importance than achieving the
minimal weight variance. Thus, the learning algorithm was modified accordingly.
By observing the distributions of the event weights before and after the inclusion
of the modified {\tt ac-VEGAS} algorithm (Fig.~\ref{f:weights}), it is evident
that {\tt ac-VEGAS} quite efficiently clusters the weights at lower values.
\begin{figure}
%\vspace{-.8cm}
\begin{center}
\mbox{
     \includegraphics[width=8cm]{zbb_wei.eps}
}
\end{center}
\vspace{-0.5cm}
\caption{\em
The distribution of event weights using  the Multi-Channel approach only (dotted histogram)
 and after application of the {\tt VEGAS} (dashed histogram) and {\tt ac-VEGAS} (full histogram) 
algorithms in the $ g g \to (Z^0 \to) l \bar{l} b \bar{b}$ process.\label{f:weights} 
}
\end{figure}

\boldmath
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\noindent{\bf Colour Flow Information}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\unboldmath
\vspace{0.2cm}

Before the generated events are passed to {\tt PYTHIA/HERWIG} for further treatment,
additional information on the colour flow/connection of the event has to be obtained.  To
provide an illustration, the method of the colour flow determination is described for the 
processes $ g g \to t \bar{t} b \bar{b}$.

For the process $ g g \to t \bar{t} b \bar{b}$ six colour flow configurations are
possible, as shown in Figure~\ref{f:cols}. With 36 Feynman diagrams contributing to the
process and at least half of them participating in two or more colour flow configurations,
calculations by hand would prove to be very tedious. Consequently, a slightly modified
colour matrix summation procedure from {\tt MADGRAPH} \cite{Stelzer:1994ta} was used to
determine the colour flow combinations of the diagrams and the corresponding colour
factors. The derived squared matrix elements for the separate colour flow combinations
$|{\mathcal{M}_{\rm flow}}|^2$ were used as sampling weights on an event-by-event basis to
decide on a colour flow configuration of the event before it is passsed on to {\tt PYTHIA/HERWIG}
for showering and fragmentation. The procedure was verified to give identical results
for the colour flow combinations and corresponding colour factors when applied to
the processes published in \cite{Bengtsson:1985yx}.
\begin{figure}
\vspace{-0.6cm}
\begin{center}
\mbox{
      \includegraphics[width=8cm]{ttbb_cols.eps}
]}
\end{center}
\vspace{-0.6cm}
\caption{\em 
A diagrammatic representation of the six colour flow configurations in the process
$g g \to t \bar{t} b \bar{b}$.\label{f:cols}
}
\end{figure}

\boldmath
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{How to use the package} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\unboldmath

The {\bf AcerMC 1.0} package consists of the library of the matrix-element-based generators
for the selected processes, interfaces to the {\tt PYTHIA 6.2} and {\tt HERWIG 6.3}
generators, and two main programs: {\tt demo\_hw.f} and {\tt demo\_py.f}. The  makefiles provided
allow the user to build executables with either of these generators:
{\tt demo\_hw.exe} or {\tt demo\_py.exe}.

There are two steering input files: {\tt run.card} and {\tt acermc.card}, which have the
same form for both executables. 
The {\tt run.card} file provides switches for modifying the generated process, the number of events, 
the structure functions, and predefined option for hadronisation/fragmentation, random number, etc..
The  {\tt acermc.card} file provides switches for modifying more specialized settings for
the  {\bf AcerMC 1.0} library itself. Once the user decides on the setup for the generated process,
only the {\tt run.card} very likely needs to be modified for the job submission.



The same executables can be also used for generating standard  {\tt PYTHIA 6.2}/{\tt HERWIG 6.3}
processes. Examples of how to run such jobs are provided as well,
in respectively {\tt demo\_hw.f} and {\tt demo\_py.f}.
If the user requires that the  {\bf AcerMC} library is not used,  $ttH$ production will 
be generated with  {\tt demo\_py.exe}, and the  {\tt HERWIG 6.3} implementation of the $Zbb$ production
will be generated with {\tt demo\_hw.exe}. Only in this case will the {\tt run.card} file be
read, so the user should  implemente her/his steering there or create another xxx.card file,
and add the respective code to the {\tt demo\_xx.f}.

The {\bf AcerMC 1.0} matrix-element-based generators are very highly optimized, using
multi-channel optimization and additional improvement with the {\tt VEGAS} grid.  The
generation modules require three kinds of the input data to perform the generation of
unweighted events: {\bf (1)} A file containing the list of the values of relative channel
weights obtained by the multi-channel optimization.  {\bf (2)} A file containing the
pre-trained {\tt VEGAS} grid.  {\bf (3)} A file containing the maximum weight $wt_{\rm
max}$, $\epsilon$-cutoff maximum weight $wt_{\rm max}^\epsilon$ and the 100 events with
the highest weights. This means that in the case of changing the default running conditions,
such as the structure functions or centre-of-mass energy, in order to recover the initial efficiency for
event generation, the user should repeat the process of preparation of the internal data files
with the inputs for the phase-space generator modules.  Pre-trained data sets obtained
using $\sqrt{s}=14\;$ TeV, {\tt PYTHIA} default $\alpha_s(Q^2)$ and $\alpha_{\rm
QED}(Q^2)$ and the CTEQ5L (parametrized) structure function set are already provided for each
implemented process.

The number of required input files might at first glance seem large, considering
that many event generators do not require any input files for operation; the
difference is not so much in the complexity of the phase space generation as
in the fact that many event generators require a {\it warming run} instead.
That is,  before the generation of unweighted events is performed a certain number of
weighted events (typically of the order of $10^4$) is generated in order to
obtain the relative multi-channel weights (in case multi-channel phase space
generation is used) and/or the optimized {\tt VEGAS} grid and/or an estimate of
the maximal weight. Such an approach can have an advantage when event generation
is very fast and the phase space regions with the highest weights are well known
(as done for the $2 \to 2$ processes in {\tt PYTHIA}); on the other hand, when
the phase space topology of the process is more complex and the event generation
is comparatively slow, generating a relatively small number of e.g. $10^4$
weighted events {\it every time} a generator is started can become CPU wasteful
and/or inaccurate in terms of maximum weight estimation.

 Reasonably accurate
estimation of the latter is mostly crucial for correct event unweighting; event
generators using the {\it  warming-up} method for maximal weight search often find still 
higher weights during the normal run and reset the maximal weight
accordingly. In this case however, the statistically correct approach would be to
reject all events generated beforehand and start the event generation anew.
This is almost never implemented due to the CPU consumption and the possibility
of hitting a weak singularity. With a small pre-sampled set the
generator can, however, badly under-estimate the maximum weight and a large number
of events can be accepted with a too-high probability. The only hope
of obtaining correct results  in such cases is that the weight {\it plateau} will be hit
sufficiently early in the event generation process. Consequently, such an approach
can be very dangerous when generating  small numbers of events\footnote{{\it Small}
is a somewhat relative quantifier, since the size of an representative sample 
should depend on the phase space dimension, i.e. the number of particles in the
final state. For example, with 4 particles in the final state, $10^5$ events can still
be considered  relatively small statistics.}.



In contrast to the {\it warming-up} approach, we have decided that
using separate training runs with large numbers of weighted events to obtain the
optimized grids and maximum weight estimates is preferable. In case the user wants
to produce data sets for non-default settings, this can easily be done
by configuring the switches in the {\tt acermc.card}


\boldmath
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Outlook and conclusions} 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\unboldmath

We have presented here the {\bf AcerMC 1.0} Monte Carlo Event Generator,which is based on a library
of  matrix-element-based generators and an interface to the universal event generator
{\tt PYTHIA~6.2} and {\tt HERWIG 6.3}. The interface is based on the standard proposed
in section \ref{interface}. 

The presented library fulfils following goals:
\begin{itemize}
\item
It provides a possibility to generate a few Standard Model background processes which
were recognised as being very dangerous for searches for  {\it New Physics} at LHC,
and for which generation was either unavailable or not straightforward thus far.
\item
Although the hard process event is generated with a matrix-element-based generator, the provided
interface allows a complete event to be generated  with either {\tt PYTHIA~6.2} or {\tt HERWIG 6.3}
initial and final state radiation, multiple interaction, hadronization, fragmentation 
and decays.
\item
The interface can be also used for studying systematic differences between 
 {\tt PYTHIA~6.2} or {\tt HERWIG~6.3} predictions for the underlying QCD processes.
\end{itemize}


These complete the list of the native {\bf AcerMC} processes implemented so far is:
$q \bar q \to W(\to \ell \nu) b \bar b$,
 $gg, q \bar q \to Z/\gamma^*(\to \ell \ell) b \bar b$, 
QCD $gg, q \bar q \to t \bar t b \bar b$ and EW $gg \to (Z/W/\gamma^* \to) t \bar t b \bar b$,
 $q \bar q \to W(\to \ell \nu) t \bar t$ and
 $gg, q \bar q \to Z/\gamma^*(\to \ell \ell, \nu \nu, b \bar b) t \bar t$.
We plan to extend this crucial list of processes,
gradually in the near future.

\subsection{Comparisons of Higgs Boson Properties with Soft Gluon Emission:
  Analytic and Parton Showering Methods}\label{Huston}

In the near future, experiments at the Tevatron and the LHC will
search for evidence of both the Higgs boson and new phenomena 
that supersede the Standard Model.
Important among the tools that will be used in these searches are
event generators based on parton showering (PS-EG's).  The most
versatile and popular of these are the Monte Carlos {\tt HERWIG}
\cite{Corcella:2001bw,Marchesini:1988cf,Knowles:1988vs,Catani:1991rr}, {\tt ISAJET}  \cite{Baer:1999sp}, and {\tt PYTHIA} 
\cite{Sjostrand:2000wi,Sjostrand:2001yu}.
PS-EG's are useful because they accurately describe the emission of
multiple soft gluons (which is, in effect, an all orders problem in
QCD) and also allow a direct connection with non--perturbative models
of hadronization.  In the parton shower, energy--momentum is conserved
at every step, and realistic predictions can be made for arbitrary
physical quantities.
However, the prediction of the  total cross section is only accurate to leading
order accuracy, and, thus, can demonstrate a sizable dependence on
the choice of scale used for the parton distribution functions (PDF's)
or coupling constants (particularly $\alpha_s$).  Also, in general,
they do not accurately describe kinematic configurations where a hard
parton is emitted at a large angle with respect to other partons.  In
distinction to PS-EG's are certain analytic calculations which account
for multiple soft gluon emission and higher order corrections to the
hard scattering simultaneously using resummation.  
The resummation technique systematically includes towers of logarithms
which are formally of the same order using the renormalization group.
These calculations,
however, integrate over the kinematics of the soft gluons, and, thus,
are limited in their predictive power.  They can, for example,
describe the kinematics of a heavy gauge boson produced in hadron
collision, but cannot predict the number or distribution of jets that
accompany it.  

Much recent work has focused on 
correcting the parton shower predictions to reproduce the
hard emission limit, where the exact leading order
matrix element gives an accurate description,
with work proceeding on 
extending this correction to
next-to-leading order~\cite{Miu:1998ju,Corcella:1999gs,Chen:2001ci,Collins:2000gd}. 
In order to match this precision, it is also
important to verify that the PS-EG programs correctly reproduce
the expected logarithmic structure in their simulation of multiple,
soft gluon emission.
The best approximation we have of the
expected logarithmic structure is represented by resummation calculations.

Since 
Standard Model Higgs boson production is a primary
focus of the physics program
at the Tevatron and  the LHC, and since several PS-EG and analytic resummation
predictions exist, we have chosen $gg\to H$ production as
a benchmark for evaluating the consistency and accuracy of 
these two approaches.
In particular, the transverse momentum of the Higgs boson 
$Q_T^H$ depends primarily
on the details of the soft gluon emission from the initial state
partons. 
The $gg$ initial state is particularly interesting, since the large color
charge may emphasize any differences that might exist in parton
shower/resummation implementations.  
Furthermore, for the $gg$ induced process, the details of non--perturbative
physics (e.g. intrinsic $k_T$) are less important \cite{Balazs:2000sz,Catani:2000zg}.

To this end, we have compared Higgs production using {\tt HERWIG},
{\tt PYTHIA} (several recent versions) and ResBos \cite{Balazs:1997xd} for 
Higgs masses of 125 and 500 GeV at center-of-mass energies of 1.96 TeV 
(${{\rm p}\bar{\rm p}}$), 14 and 40 TeV (${\rm pp}$).  
The two different masses and three
different center of mass energies provide a wide variation of kinematics 
that test the showering/resummation processes.
This work extends our results published in \cite{Djouadi:2000gu,Catani:2000zg}, 
\cite{Balazs:2000wv} and \cite{Balazs:2000sz}.

\subsubsection{Parton Showers}

PS-EG's are based on the factorization theorem \cite{Bodwin:1985ft}, which,
roughly, states that physical observables in any sensible gauge theory
are the product of short--distance functions and long--distance
functions.  The short--distance functions are calculable in
perturbation theory.  The long--distance functions are fit at a scale,
but their evolution to any other scale is also calculable in
perturbation theory.

A standard application of the factorization theorem is to describe
heavy boson production at a hadron collider to a fixed order in
$\alpha_s$.  The production cross section is obtained by convoluting
the partonic subprocesses evaluated at the scale $Q$ with the PDF's
evaluated at $Q$.  The partons involved in the hard collision must be
sufficiently virtual to be resolved inside the proton, and a natural
choice for the scale $Q$ is the mass of the heavy boson \cite{Bengtsson:1986gz}.
Prior to the hard collision, however, the partons are not resolvable
in the proton (i.e. the proton is intact) and have virtualities at a
much lower scale $Q_0$ of the order of 1 GeV.  The connection between
the partons at the low scale $Q_0$ and those at the high scale $Q$ is
described by the DGLAP evolution equations \cite{Gribov:1972rt,Altarelli:1977zs,Dokshitzer:1977sg}.  
The DGLAP
equations include the most important kinematic configurations of the
splittings $a \to b c$, where $a,b$ and $c$ represent different types
of partons in the hadron ($q,g$, etc.).  Starting from a measurement
of the PDF's at a low scale $Q_0$, a solution of the DGLAP equations
yields the PDF's at the hard scale $Q$.  Equivalently, starting with a
parton $c$ involved in a hard collision, it is also possible to
determine probabilistically which splittings generated $c$.  In the
process of de--evolving parton $c$ back to the valence quarks in the
proton, a number of spectator partons (e.g. parton $b$ in the
branching $a\to bc$) are resolved.  These partons constitute a shower
of soft jets that accompany the heavy boson, and influence its
kinematics.

The shower described above occurs with unit probability and does not
change the total cross section for heavy boson production calculated
at the scale $Q$ \cite{Odorico:1980gg}.  The showering can be attached to the
hard--scattering process based on a probability distribution {\it
  after} the hard scattering has been selected.  Once kinematic cuts
are applied, the transverse momentum and rapidity of the heavy boson
populate regions never accessed by the differential partonic cross
section calculated at a fixed order.  This is consistent, since the
fixed--order calculation was inclusive and
was never intended to describe the detailed kinematics of the
heavy boson.  The parton shower, in effect, resolves the structure of
the inclusive state of partons denoted as $X$.  In practice, the fixed
order partonic cross section (without showering) can still be used to
describe properties of the decay leptons as long as the measurable is
not highly correlated with the heavy boson kinematics.

Here, we review parton showering schematically.  More details can be
found, for example, in Ref. \cite{Ellis:1996yz}.  First, for simplicity, consider the case
of final state or forward showering, where the parton virtuality $Q$
evolves forward to the low scale $Q_0$.  The basis for developing a
probabilistic picture of final state showering is the DGLAP equation
for the fragmentation functions: \bea Q {\partial \over \partial Q}
D_a(x,Q) = \int_{x}^{1-\epsilon}{dz \over z}{\alpha_{abc}(z,Q) \over
  \pi}
\hat{P}_{a\to bc}(z) D_b(x/z,Q) \nonumber \\
-D_a(x,Q) \int_{x}^{1-\epsilon}{dz}{\alpha_{abc}(z,Q) \over \pi}
\hat{P}_{a\to bc}(z), \eea where $\hat{P}_{a\to bc}$ is an
unregularized splitting function, $\alpha_{abc}$ is the coupling times
color factor, and $\epsilon$ is a cutoff.  The equation can be
rewritten as \bea {\partial \over \partial \ln Q^2} \left(
  D_a(x,Q)/\Delta(Q) \right) = \int_{x}^{1}{dz \over
  z}{\alpha_{abc}(z,Q) \over 2\pi} \hat{P}_{a\to bc}(z)
(D_b(x/z,Q)/\Delta(Q)) \nonumber \eea or, after integrating both sides
of the expression, \bea D_a(x,t') = D_a(x,t) \Delta(t') +
\int_{t'}^{t}\int_{x}^{1} dt'' {dz\over z} {\Delta(t') \over
  \Delta(t'')} {\alpha_{abc}(z,t'') \over 2\pi} \hat{P}_{a\to bc}(z)
D_b(x/z,t''), \eea where $t=\ln Q^2,$ with similar definitions for
$t'$ and $t''$.  
The function 
\bea 
   \Delta(t^\prime) = \exp\left(
  -\int_{t_0}^{t'}\int_{\epsilon}^{1-\epsilon} dt'' dz
  {\alpha_{abc}(z,t'') \over 2\pi} \hat{P}_{a\to bc}(z) \right) 
\eea
is called the Sudakov form factor, 
and gives the probability of evolving from the scale $t^\prime$ to $t_0$ 
with no resolvable branchings, 
where  $t_0$ is a cutoff scale for the showering.
The Sudakov $\Delta(t^\prime)$ contains all the information necessary to reconstruct a
shower, since it encodes the change in virtuality of a parton until a
resolvable showering occurs.
Showering is reduced to iterative solutions of the equation
$r=\Delta(t^\prime)/\Delta(t'')$, where $r$ is a random number uniformly
distributed in the interval $[0,1]$, until a solution for $t^\prime$ is
found which is below the cutoff $t_0$.  

For the case of initial state radiation, several modifications are
necessary.  The fragmentation function is replaced by a parton
distribution function, and the evolution proceeds backwards from a
large, negative scale $-|Q^2|$ to a small, negative cutoff scale
$-|Q_0^2|$.  There are two equivalent formulations of backwards
showering based on the probabilities 
\bea 
  \exp\left(-\int_{t'}^{t}\int_{\epsilon}^{1-\epsilon} dt'' dz
  {\alpha_{abc}(z,t'') \over 2\pi} \hat{P}_{a\to bc}(z) {x'
    f_a(x',t')\over x f_b(x,t') } \right), x'=x/z, 
\label{algo1}
\eea \cite{Sjostrand:1985xi},
and 
\bea 
{\Delta(t') \over f_b(x,t')} {f_a(x,t'') \over \Delta(t'') }
\label{algo2}
\eea 
\cite{Marchesini:1984bm}.
%
After choosing the change in virtuality, a particular backwards
branching is selected from the probability function based on their
relative weights (a summation over all possible branchings $a\to bc$
is implied these expressions), and the splitting variable is a solution
to the equation 
\bea 
  \int_\epsilon^{x/x'} {dz\over z}
  \hat{P}_{a\to bc}(z)f(x/z,t') = r \int_\epsilon^{1-\epsilon} {dz\over
  z} \hat{P}_{a\to bc}(z)f(x/z,t'), 
\eea 
where $r$ is a random number.
The details of how a full shower is reconstructed in {\tt PYTHIA},
 for example, can be found in
Ref.~\cite{Sjostrand:1994yb}.  The structure of the shower can be complex: the
transverse momentum of the heavy boson is built up from the whole
series of splittings and boosts, and is known only at the end of the
shower, after the final boost.


The PS-EG formulation outlined above is fairly independent of the hard
scattering process considered.  Only the initial choice of partons and
possibly the high scale differs.  Therefore, this formalism can be
applied universally to many different scattering problems.  In effect,
soft gluons are not sensitive to the specifics of the hard scattering,
only the color charge of the incoming partons.  This statement is true
to leading logarithm.

The parton showering of {\tt PYTHIA} obeys a strict ordering in
virtuality: the parton that initiates a hard scattering has a larger
(negative for initial state showers) virtuality than any other parton in the shower.  Parton
showers in {\tt HERWIG} proceed via a coherent branching process in which
a strict angular ordering is imposed on sequential gluon emissions:
the evolution variable is not virtuality, but a generalized virtuality
$\xi$.  
For an initial state shower, with parton splitting $c\to ba$, and
where $a$ has the largest virtuality, the variable $\xi=(p_b \cdot
p_c)/(E_bE_c)$.
At all values of $x$, the coherent branching algorithm
correctly sums the leading logarithmic contributions. At large $x\sim 1$, it
also sums the next-to-leading order contributions \cite{Catani:1991rr}, with
an appropriate definition of the splitting kernel.
(The exact definition of LL and NLL will be given later.)
Because of
the demonstrated importance of coherence effects, {\tt PYTHIA} includes
an additional veto on showers which are not also angular--ordered.
Note, however, that this does not make the two schemes equivalent --
some late emissions in a {\tt HERWIG} shower can have virtuality larger
than previous emissions.

\subsubsection{Analytic results}

At hadron colliders, the partonic cross sections for heavy boson
production can receive substantial corrections at higher orders in
$\alpha_s$.  This affects not only the total production rate, but also 
the kinematics of the heavy boson.  At leading order, 
the heavy boson has a $\delta(Q_T^2)$ distribution in $Q_T^2$.  At
next--to--leading order, the real emission of a single gluon generates
a contribution to $d\sigma/dQ_T^2$ that behaves as
$Q_T^{-2}\alpha_s(Q_T^2)$ and $Q_T^{-2}\alpha_s(Q_T^2)\ln(Q^2/Q_T^2)$
while the soft, and virtual corrections are
proportional to $-\delta(Q_T^2)$.  At higher orders, the most singular
terms follow the pattern of
$\alpha_s^n(Q_T^2)\sum_{m=0}^{2n-1}\ln^m(Q^2/Q_T^2)$ $=\alpha_s^n
L\equiv V^n$.  The logarithms arise from the incomplete cancellation
of the virtual and real QCD corrections. This cancellation becomes
complete for the integrated spectrum, where the real gluon can become
arbitrarily soft and/or collinear to other partons.  The pattern of
singular terms suggest that perturbation theory should be performed in
powers of $V^n$ instead of $\alpha_s^n$.  This reorganization of the
perturbative series is called resummation.

The first studies of soft gluon emission resummed the leading
logarithms \cite{Dokshitzer:1978yd,Parisi:1979se}, leading to a suppression of the cross
section at small $Q_T$.  The suppression underlines the importance of
including sub--leading logarithms \cite{Ellis:1981sj}.  The most rigorous
approach to the problem of multiple gluon emission is the
Collins--Soper--Sterman (CSS) formalism for transverse momentum
resummation \cite{Collins:1982zc}, which resums all of the important
logarithms.  This is achieved after a Fourier transformation with
respect to $Q_T$ in the transverse coordinate $b$, so that the series 
involving the delta function and terms $V^n$ simplifies to the form of an
exponential.  Hence, the soft gluon emission is resummed or
exponentiated in this $b$--space formalism.  The
Fourier transformation is the result of expressing the
transverse--momentum conserving delta functions
$\delta^{(2)}(\vec{Q}_T - \sum \vec{k}_{T_i})$ in their Fourier
representation.  Also, the exponentiation is accomplished through the
application of the renormalization group equation.  
Despite the successes of the $b$--space
formalism, there are several drawbacks. Most notable for the present
study is that it integrates out the soft gluon dynamics and does not
have a simple physical interpretation.

The CSS formalism was used by its authors to predict both the total
cross section to NLO and the kinematic distributions of the
heavy boson to resummed NLL order \cite{Collins:1985kg} at hadron colliders.  A
similar treatment was presented using the AEGM formalism \cite{Altarelli:1984pt},
that does not involve a Fourier transform, but is evaluated directly
in transverse momentum $Q_T$ space.  When expanded in $\alpha_s$, the 
two formalisms are equivalent to the NNNL order, and agree with
the NLO fixed order calculation of the total cross section \cite{Altarelli:1985kp}.
A more detailed numerical comparison of the two predictions can be
found in Ref. \cite{Arnold:1991yk}.  The AEGM formalism has been
reinvestigated, and an approximation to the $b$--space formalism has
been developed in $Q_T$--space which retains much of its predictive 
features \cite{Ellis:1998ii}.

In the $b$--space formalism, the differential cross section of the
heavy boson produced in association with soft gluons is: 
\bea
 {d\sigma(h_1 h_2 \to B^{(*)}X) \over dQ^2\,dQ^2_T\,dy} = {1\over (2
 \pi)^2} \int_{}^{} d^2 b \, e^{i {\vec b}\cdot{\vec Q_T}}
 \widetilde{W}(b,Q,x_1,x_2) + Y(Q_T,Q,x_1,x_2).  
\label{Eq:CSSFormula}
\eea 
where $Q$, $Q_T$ and $y$ describe the kinematics of the off--shell
heavy boson $B^{(*)}$. The function $Y$ is free of $\ln(Q^2/Q_T^2)$ and 
corrects for the soft gluon approximation in the high $Q_T$ region. 
The function $\widetilde W$ has the form: 
\bea 
 {\widetilde W}(b,Q,x_1,x_2) = e^{-S(b,Q)} 
 \left( C_{il} \otimes f_{l/h_1} \right)(x_1,b)
 \left( C_{jl} \otimes f_{l/h_2} \right)(x_2,b)
% \Bigg(C \otimes f\Bigg)(x_1,b)
% \Bigg(C \otimes f\Bigg)(x_2,b) 
 H_{ij}(Q,y), 
\eea 
where 
\bea & &
 S(b,Q,C_1,C_2) = \int_{C_1^2/b^2}^{C_2^2Q^2} {d {\bar \mu}^2\over
 {\bar \mu}^2} \left[ \ln\left({C_2^2Q^2\over {\bar \mu}^2}\right)
 A\big(\alpha_s({\bar \mu})\big) + B\big(\alpha_s({\bar \mu})\big)
 \right], 
\eea 
and 
\bea 
 \left( C_{jl} \otimes f_{l/h_1} \right)
 (x_1,\mu) = \int_{x_1}^{1} {d \xi_1 \over \xi_1} \, C_{jl}\left( {x_1 \over
 \xi_1}, C_1, C_2, \mu=C_3/b\right) f_{l/h_1}\left(\xi_1, \mu=C_3/b\right).  
\eea 
In these expressions, $C_1$, $C_2$ and $C_3$ are renormalization scales, 
$H$ is a function that describes the hard scattering, and $A$, $B$ and 
$C$ are calculated perturbatively in powers of $\alpha_s$: 
\bes
 \{A,B\}=\sum_{n=1}^{\infty} \left({\alpha_s(\mu)\over\pi}\right)^n
 \{A^{(n)},B^{(n)}\} ~~~ {\rm and} ~~~
  C_{ij} =\sum_{n=0}^{\infty} \left({\alpha_s(\mu)\over\pi}\right)^n C^{(n)}.
\ees 
The functions $C^{(n)}$ are mostly responsible for the change in the total
production cross section at higher orders.  In fact, $\left( C\otimes f\right)$ 
is simply a redefinition of the parton distribution
function obtained by convoluting the standard ones with an
ultraviolet--safe function. These generalized parton distributions encode
both the longitudinal momentum fraction and the transverse recoil of the 
initial state parton.

We remove $C_1, C_2$ and $C_3$ from these expressions by
choosing their canonical values \cite{Collins:1982zc}, which also removes
large logarithms from the expansion of the resummed expression. 
At leading order in $C_{ij}$, the resummed part of the expression for the 
production of an on--shell heavy boson simplifies considerably to: 
\bea
  {d\sigma(h_1 h_2 \to B^{(*)}X) \over dQ^2_T} = \sigma_0 \int_{}^{} {d^2 b
  \over (2 \pi)^2} \, e^{i {\vec b}\cdot{\vec Q_T}} e^{-S(b,Q)}
  {f(x_1,b) f(x_2,b) \over f(x_1,Q) f(x_2,Q)}, 
\eea 
where
\bes 
  \sigma_0 = \kappa \int_{}^{} {dx_1 \over x_1} f(x_1,Q) f(x_2,Q), 
\ees 
and $\kappa$ contains physical constants.  
The expression contains two factors, the total cross section
at leading order $\sigma_0$, and a cumulative probability function in
$Q_T^2$ that describes the transverse momentum of the heavy boson
(the total integral over $Q_T^2$ transforms
$e^{i\vec{b}\cdot\vec{Q}_T}$ to $\delta^{(2)}(\vec{b})$).  Except for
the complication of the Fourier transform, the term $e^{-S/2}
f(x,b)/f(x,Q)$ is analogous to $\Delta(Q) f(x,Q')/\Delta(Q') f(x,Q)$
of the PS-EG of Eq.~(\ref{algo2}).

Equation (\ref{Eq:CSSFormula}), which is formulated in $b$--space, has a 
similar structure in $Q_T$--space.
This is surprising, since the $b$--space result depends critically on
the conservation of total transverse momentum.  To NNNL accuracy,
however, the $Q_T$ space expression agrees exactly with the $b$--space
expression, and has the form \cite{Ellis:1998ii}: 
\bea 
  {d\sigma(h_1 h_2 \to
  B^{(*)}X) \over dQ^2\,dQ^2_T\,dy} = {d\over
  dQ_T^2}\widetilde{W}(Q_T,Q,x_1,x_2) + Y(Q_T,Q,x_1,x_2).  
\eea 
Again ignoring $Y$, we can rewrite this expression as: 
\bea {d\sigma(h_1 h_2
  \to B^{(*)}X) \over dQ_T^2 } = \sigma_1 \left({d\over dQ_T^2}
  \left[e^{-S(Q_T,Q)} {\left(C\otimes f\right)(x_1,Q_T) \left(C\otimes
        f\right)(x_2,Q_T) \over \left(C\otimes f\right)(x_1,Q)
      \left(C\otimes f\right)(x_2,Q)} \right]\right), \eea \bes
\sigma_1 = \kappa \int_{}^{} {dx_1 \over x_1} \left(C\otimes
  f\right)(x_1,Q) \left(C\otimes f\right)(x_2,Q).  \ees The factor
$\sigma_1$ is the total cross section to a fixed order, while the rest
of the function yields the probability that the heavy boson has a
transverse momentum $Q_T$.

\subsubsection{Methodology for comparison}\label{sec:resumps,mc,qcdsm}

To make a comparison between the distributions from analytic and parton
showering calculations, we must quantify the differences in theoretical input
and identify what approximations have been made in each one.   >From the
discussion of the analytic resummation calculations, we see that the emission
of multiple soft gluons is described best by perturbation theory not in terms
of $\alpha_s$, but in powers of $\alpha_s$ times logarithms of large  numbers.
Typically, the logarithms are classified  according to the orders of $\alpha_s$
and the power of logarithms in  $(i)$ the Sudakov exponent, or in $(ii)$ the
perturbative expansion of  the resummed $\qt$ distribution.

In classification $(i)$\footnote{See also sect. \ref{sec:low-qt,resum,qcdsm}}, 
it is argued that terms with $A^{(n)}$ are leading
compared to  terms with $B^{(n)}$ because the former are enhanced by a large
logarithm.  Also, $X^{(n)}$ is leading compared to $X^{(n+1)}$ because the
latter is  suppressed by $\alpha_s$ for $X = A, B$ or $C$. The comparison of
$C$ to $A$ and  $B$ is somewhat {\it ad hoc}, because $C$ does not appear in
the Sudakov in all  resummation schemes. In this approach, the lowest order
resummed result (LL) has  only $A^{(1)}$ and $C^{(0)}$. Additionally, the next
order (NLL) has $A^{(2)}$, $B^{(1)}$ and $C^{(0)}$, because the first two
are suppressed by an $\alpha_s$  and the last by the inverse of a large
logarithm when compared to the leading  terms. This is an elegant and simple
classification scheme, because it  relies only on the structure of the Sudakov
exponent. On the other hand, it has an {\it ad hoc} mixing $\alpha_s$ and
$\log^{-1}$--suppressed terms, and  depends on the resummation scheme in
dealing with $C^{(n)}$.

In classification $(ii)$, the resummed cross section is expanded and reorganized in 
towers of logarithms as 
  \bea
  \frac{d\sigma}{dQ_T^2} \sim 
\begin{array}{llll}
\alpha_s^1 ( L_1^1 ) &   +      &          &   \\
\alpha_s^2 ( L_2^1   &+ L_2^3 ) & +        &   \\
\alpha_s^3 ( L_3^1   &+ L_3^3   & + L_3^5 )& + \cdots,
\end{array}
  \eea
where $L_i$ are various linear combinations of $\ln(Q/Q_T)$.
Remarkably, each column of logs receives contribution only from $A^{(n)}$, 
$B^{(n)}$ and $C^{(n-1)}$ ($n$ numbers the columns). Since $X^{(n)}$ 
and $X^{(n+1)}$ are clearly ordered by $\alpha_s$ ($X = A,B,C$),
the first column is called the leading tower of logs, i.e. LL, the next NLL, etc. 
This is a more--involved classification, but closer to the spirit of the 
resummation, and more precise. In the rest of this work, we use the naming 
convention $(ii)$ when comparing to the analytic resummation.
Thus, when we say that the coherent branching algorithm in {\tt HERWIG} correctly sums the 
leading logarithmic contributions, we are referring to the $A^{(1)}$ and $B^{(1)}$ terms.

\begin{figure}[!ht]
%%  \resizebox{8.1cm}{8cm}{
%%    \includegraphics*{fig_qt_logs.eps} }
\centering
  \resizebox{10cm}{10cm}{ 
\includegraphics*{fig_qt_abc.eps}
    }
    \caption{\em (top) Absolute $Q_T$ distribution of the Higgs boson from an analytic resummation
calculation to different orders of accuracy in perturbation theory;
      (left, bottom) Same as (top), but normalized to the same total
cross section; (right, bottom) Same as (left,bottom), but on a logarithmic scale.  
All curves are for the production of a 125 GeV Higgs boson
at the LHC based on the CTEQ4M parton distribution functions.}
    \label{fig:abc}
\end{figure}

Figure~\ref{fig:abc} demonstrates the different predictions
for different choices of organizing the perturbative expansion.
The top frame shows the absolute distributions of $Q_T^H$ when the
perturbative
coefficients include: $A^{(1,2)},~B^{(1,2)},~C^{(0,1)}$ (solid red),
$A^{(1,2)},~B^{(1)}$ and $C^{(0,1)}$ (long-dashed blue), $A^{(1)},~B^{(1)}$
and
$C^{(0)}$ (solid light blue), and only $A^{(1)}$ and $C^{(0)}$ (short-dashed
black).
The change from LL to NLL results in an increase of the normalization of the
cross section as well as a slight shift of the peak towards higher
transverse
momentum, but the shapes are very similar. The lower left frame shows the
same
distributions, but all normalized to the same production rate. This
eliminates
the normalization dependence and focuses on the variation of the shape. The
lower
right frame, shown in logarithmic scale, amplifies the differences at high
$Q_T$'s.

\subsubsection{Numerical comparisons}

In the following plots, we show several numerical comparisons of ResBos,
{\tt HERWIG} and three versions of {\tt PYTHIA}.  Why three versions of {\tt PYTHIA}?
This reflects the evolution of the program.  
The version {\tt PYTHIA}-5.7 was used for almost all LHC Monte Carlo analyses for Higgs production.
The showering in {\tt PYTHIA}-5.7, as described earlier, is virtuality ordered, 
%
with an additional requirement of angular ordering.  No other kinematic restrictions
are imposed.  In later versions, beginning with {\tt PYTHIA}-6, each parton emission
was required to match the kinematic constraints of the NLO process, i.e. for
the first emission in the parton shower accompanying $gg\to H$,
the kinematic constraints $\hat{s}+\hat{t}+\hat{u}=M_H^2$ must be satisfied.
Since $\hat{s}$ and $\hat{t}$ are the two variables of the virtuality--ordered
shower, a cut on $\hat{u}$ restricts their range.  The plots labeled 
by {\tt PYTHIA}6.1 demonstrate the effect of the  $\hat{u}$ cut.
% 
Finally, in {\tt PYTHIA}-6.2, an additional hard matrix element correction was
applied to the parton shower.  In this approach, the maximum virtuality of 
the shower is increased from its nominal value at $Q=m_H$ up to the largest
kinematically allowed value.  Furthermore, each parton emission is corrected
by the ratio of the exact matrix element squared at NLO to the approximate
matrix element squared given by the parton shower approximation.
Thus, showers generated using {\tt PYTHIA}-6.2 should have a closer agreement
to the ResBos predictions at high $Q_T^H$ (where ResBos relies on
the NLO prediction).  In contrast, showers induced in {\tt HERWIG} will still
have a cutoff set by $m_H$.  Note, however, that the {\tt PYTHIA} prediction
still has the LO normalization, and a full rescaling is necessary.  
This may be appropriate if the $K$--factor correction to the LO prediction
is the same as the one for the NLO emission.  

Figures~\ref{fig:tev125} and \ref{fig:lhc125} show the predicted $Q_T$ 
distributions for
production of a Standard Model Higgs boson with mass $m_H=125$ GeV
at collider energies of 1.96, and 14 TeV. 
For all plots, the distributions have been normalized to the same cross
section (ResBos NLL); without this normalization, the
PS-EG predictions would be about a factor of 2 lower than the ResBos
curves.  
Two
ResBos curves are shown: LL (including $A^{(1)}$, $B^{(1)}$ and $C^{(0)}$) and
NLL (also including $A^{(2)}$, $B^{(2)}$ and $C^{(1)}$ ). The inclusion of the NLL terms leads to a
slightly harder $Q_T$ distribution, as discussed previously. The
ResBos curves appear close to the {\tt HERWIG} predictions, and somewhat less
close to the predictions of {\tt PYTHIA} (versions 6.1 and after).

In general, the PS-EG predictions are in fair agreement with
the analytic resummation ones for low $Q_T$, where multiple, soft
gluon emission is the most important.
The agreement of {\tt HERWIG} with the ResBos curves becomes even
better if the shape comparison are made by normalizing the
cross sections in the low $Q_T$ region alone, away from the effects of
the exact matrix element for Higgs plus jet. 
This is illustrated in Figures~\ref{fig:lhc125}.
The {\tt PYTHIA}-6+ predictions peak at a noticeably lower value
of $Q_T$ than either ResBos or {\tt HERWIG}.
A striking feature of the plots is the change
induced in {\tt PYTHIA}-6.1 relative to {\tt PYTHIA}-5.7,
indicating the importance of kinematic constraints.
Note also that the average Higgs transverse momentum increases with
increasing center-of-mass energy, due to the increasing phase space
available for gluon emission. 
At high $Q_T$, the exact matrix element for Higgs plus
jet, present both in the most recent version of {\tt PYTHIA} as well as
ResBos, correctly describes Higgs production at transverse momenta
on the order of the Higgs mass or larger, while a pure parton
showering description of the high $Q_T$ end is inadequate. 
The change observed in 6.2,
only visible at large values of $Q_T$, is the
result of the 
matrix element (Higgs + jet) corrections to the parton shower.

The final comparison plot, 
Figure~\ref{fig:lhc500},  shows the $Q_T$ distributions for a 500 GeV
Higgs boson
generated through $gg$ fusion at 14 TeV (pp) and 40 TeV (pp).
Comparisons are made of the ResBos(LL and NLL), {\tt PYTHIA} and
{\tt HERWIG}  predictions. As in
Figures~\ref{fig:tev125} and \ref{fig:lhc125} the distributions have been scaled to have the same
total cross section. The average transverse momentum for a 500 GeV
Higgs is noticeably larger than that for a 125 GeV Higgs in all of the
predictions, as expected, since the hard scale for the
process is $m_H$.  
The agreement between the {\tt HERWIG} (version 6.1 and
later) predictions and the ResBos  curves is better than for the
predictions with a 125 GeV mass Higgs. The agreement
between {\tt HERWIG} and the ResBos curves remains very good.
For high $Q_T$, the {\tt PYTHIA} 6.2 prediction, with matrix element
corrections, agrees with the ResBos prediction.

\begin{figure}
\centering
  \resizebox{15cm}{10cm}{
    \includegraphics*[0,50][650,600]{higgs_125_tevatron_all.eps} }
    \caption{\em Transverse momentum $Q_T$ distributions of a 125 GeV 
Higgs boson produced at the Tevatron for {\tt HERWIG}, different versions
of {\tt PYTHIA} and different perturbative orders of ResBos.  The
bottom plot is for an expanded $Q_T$ range.  All curves are
normalized to the same total cross section.}
    \label{fig:tev125}
\end{figure}

%\begin{figure}
%  \resizebox{15cm}{15cm}{
%    \includegraphics[0,0][650,650]{higgs_150_tevatron_all.eps} }
%    \caption{\em Same as Fig.~\ref{fig:tev125}, but for $m_H=150$ GeV.}
%    \label{fig:tev150}
%\end{figure}

\begin{figure}
\resizebox{8cm}{10cm}{
\includegraphics*[0,50][600,650]{higgs_125_14tev_all.eps}}
\resizebox{8cm}{10cm}{
\includegraphics[0,50][600,650]{higgs_125_14tev_all_09.eps}}
\caption{\em Comparison of light Higgs boson production at
the LHC. The right hand plots have  {\tt PYTHIA} 5.7, 6.1 and 
{\tt HERWIG} distributions normalized
to 90\% of the ResBos NLL cross section. }
    \label{fig:lhc125}
\end{figure}

\begin{figure}
\resizebox{8cm}{10cm}{
\includegraphics[0,0][600,600]{higgs_500_14tev_all.eps}}
\resizebox{8cm}{10cm}{\includegraphics[0,0][600,600]{higgs_500_40tev_all.eps}}
\caption{\em Comparison of Higgs boson production at
the LHC and 40 TeV $m_H=500$ GeV.}
    \label{fig:lhc500}
\end{figure}

\subsubsection{Properties of showers}

In the previous section, there were notable differences
in the predictions of {\tt PYTHIA} and {\tt HERWIG}, with 
{\tt HERWIG} giving a superior description of the Higgs
boson properties for low $Q_T^H$.  Therefore, we have
investigated several basic properties of the parton emissions
to determine the cause of this discrepancy.  It is already
known that the {\tt HERWIG} shower is of the coherent kind, whereas
{\tt PYTHIA} is virtuality ordered with approximate angular ordering
superimposed.  On the other hand, the same was true in
{\tt PYTHIA}5.7, but closer agreement with {\tt HERWIG} is
obtained after making the $\hat{u}$ cut.  Perhaps a more
careful analysis of the {\tt HERWIG} kinematics will result
in even better agreement.

Figure~\ref{fig:showering} shows
a comparison of showering properties for light Higgs
boson production at the Tevatron.  
For the purpose of these plots, the matrix element corrections
in {\tt PYTHIA} have been turned off.
The leftmost plot shows
the transverse momentum of the Higgs boson and the first
parton emission for {\tt PYTHIA} (red) and {\tt HERWIG} (green).
Note that the leading emission describes the Higgs boson $Q_T$ from
the full shower very well down to fairly low values $\sim 10$ GeV.
The main discrepancy between {\tt PYTHIA} and {\tt HERWIG} occurs in
the same kinematic region, exactly where the properties of
several emissions become important.  

The rightmost plot shows the largest (negative) virtuality
of a parton in the shower.  
The blue curve shows the effect of having no $\hat{u}$ cut, as
in the older version of {\tt PYTHIA}.  Clearly, {\tt PYTHIA} without
the $\hat{u}$ cut does
not have the $Q^{-2}$ behavior expected from the derivative
of the Sudakov.
While the agreement between {\tt PYTHIA}-6.2 and {\tt HERWIG} is markedly
improved, {\tt PYTHIA} still has a residual enhancement near the upper
scale $Q=m_H$ in this example.  This requires further investigation.

\begin{figure}[!ht]
\resizebox{8cm}{8cm}{
\includegraphics*{pt_first.eps}}
\resizebox{8cm}{8cm}{
\includegraphics*{virt.eps}}
\caption{\em Comparison of showering properties for light Higgs
boson production at the Tevatron.  The leftmost plot shows
the transverse momentum of the Higgs boson and the first
parton emission for {\tt PYTHIA} (red) and {\tt HERWIG} (green).
The rightmost plot shows the largest (negative) virtuality
of parton in the shower.  The blue curve shows the effect of having no $\hat{u}$ cut.}
\label{fig:showering}
\end{figure}

To summarize:-

\begin{itemize}
\item   The differences between ResBos LL and NLL are in the direction expected and are relatively
subtle.

\item   The two newer versions of {\tt PYTHIA}  and {\tt HERWIG} both approximately agree with  the predictions
of ResBos LL/NLL, with the {\tt HERWIG} shape agreement  being somewhat  better in the low
$Q_T$ region.

\item   The agreement of {\tt HERWIG} with ResBos becomes better as:
        $(a)$    the center of mass energy increases and
        $(b)$    the Higgs mass  increases
        the agreement seems to be better with LL.
\end{itemize}

\subsection{Studies of underlying events using CDF data}
\label{field}

Fig.~\ref{rdf_fig1} illustrates the way QCD Monte-Carlo models simulate a proton-antiproton collision in which 
a "hard" $2$-to-$2$ parton scattering with transverse momentum, \pthard, has occurred.  The resulting event contains 
particles that originate from the two outgoing partons ({\it plus initial and final-state radiation}) and particles 
that come from the breakup of the proton and antiproton (\ie \BBR).  The ``hard scattering" component consists 
of the outgoing two ``jets" plus initial and final-state radiation. The \UE\ is everything 
except the two outgoing hard scattered ``jets" and consists of the \BBR\ plus possible contributions 
from the ``hard scattering" arising from initial and final-state radiation. 

\begin{figure}
\begin{center}
\includegraphics[scale=0.8]{Field_fig1.eps}\end{center}
\caption{\em Illustration of the way the QCD Monte-Carlo models simulate a proton-antiproton collision in 
which a hard $2$-to-$2$ parton scattering with transverse momentum, \pthard, has occurred.  The 
resulting event contains particles that originate from the two outgoing partons 
({\it plus initial and final-state radiation}) and particles that come from the breakup of the 
proton and antiproton ({\it\BBR}).  The ``hard scattering" component consists of the outgoing two ``jets" 
plus initial and final-state radiation. The \UE\ is everything except the two outgoing hard scattered ``jets" 
and consists of the \BBR\ plus possible contributions from the "hard scattering" arising from 
initial and final-state radiation.
\label{rdf_fig1}}
\end{figure}

The \BBR\ are what is left over after a parton is knocked out of each of the initial two beam 
hadrons.  It is the reason hadron-hadron collisions are more ``messy" than electron-positron annihilations and no one 
really knows how it should be modeled.  For the QCD Monte-Carlo models the \BBR\ are an 
important component of the \UE.  Also, it is possible that multiple parton scattering contributes to the 
\UE.  {\tt PYTHIA} \cite{Sjostrand:2000wi} models the \UE\ in proton-antiproton collision 
by including multiple parton interactions. In addition to the hard $2$-to-$2$ parton-parton scattering and the \BBR, 
sometimes there is a second ``semi-hard" $2$-to-$2$ parton-parton scattering that contributes particles to the \UE.  


Of course, from a certain point of view there is no such thing as an \UE\ in a proton-antiproton 
collision.  There is only an ``event" and one cannot say where a given particle in the event originated.  On the other 
hand, hard scattering collider ``jet" events have a distinct topology.  On the average, the outgoing hadrons 
``remember" the underlying the $2$-to-$2$ hard scattering subprocess.  An average hard scattering event consists of a 
collection (or burst) of hadrons traveling roughly in the direction of the initial beam particles and two collections of 
hadrons (\ie ``jets") with large transverse momentum.  The two large transverse momentum ``jets" are roughly back 
to back in azimuthal angle.  One can use the topological structure of hadron-hadron collisions to study the 
\UE\ \cite{Field:2001ab,Field:2000ab,Tano:2001ab}. The ultimate goal is to understand the physics of the \UE, but since it is very 
complicated and involves both non-perturbative as well as perturbative QCD it seems unlikely that this will happen 
soon.  In the mean time, we would like to tune the QCD Monte-Carlo models to do a better job fitting the \UE.  
The \UE\ is an unavoidable background to most collider observables.  To find ``new" physics 
at a collider it is crucial to have Monte-Carlo models that simulate accurately ``ordinary" hard-scattering collider 
events.  This report  will compare collider observables that are sensitive to the \UE\ with the QCD 
Monte-Carlo model predictions of {\tt PYTHIA} 6.115 \cite{Sjostrand:2000wi}, {/tt HERWIG} 5.9 \cite{Marchesini:1988cf,Knowles:1988vs,Catani:1991rr}, and {\tt ISAJET} 7.32 \cite{Paige:1986vk}
and discuss the tuning of {\tt ISAJET} and {\tt PYTHIA}.

\begin{figure}
\begin{center}
\includegraphics[scale=0.8]{Field_fig3.eps}\end{center}
\caption{\em Illustrations of correlations in azimuthal angle $\Delta\phi$ relative to the direction of the 
leading charged jet in the event, chgjet\#1.  The angle $\Delta\phi=\phi-\phi_{\rm chgjet\#1}$  
is the relative azimuthal angle between charged particles and the direction of chgjet\#1.  The``toward" 
region is defined by $|\Delta\phi|<60^\circ$ and \etacut, while the ``away" 
region is $|\Delta\phi|>120^\circ$ and \etacut.   The ``transverse" region is defined by 
$60^\circ<|\Delta\phi|<120^\circ$  and \etacut.  Each region has an 
area in \etaphi\ space of $4\pi/3$.
On an event by event basis, we define ``transMAX (``transMIN") to 
be the maximum (minimum) of the two ``transverse" pieces, $60^\circ<\Delta\phi<120^\circ$ and \etacut, 
and $60^\circ<-\Delta\phi<120^\circ$ and \etacut. ``TransMAX" and ``transMIN" each have an area in \etaphi\ space 
of $2\pi/3$.  The sum of  ``TransMAX" and ``transMIN" is the total ``transverse" region with area $4\pi/3$.
\label{rdf_fig3}}
\end{figure}

\subsubsection{The ``Transverse" Region}

In a proton-antiproton collision large transverse momentum outgoing partons manifest themselves, in the laboratory, 
as a clusters of particles ({\it both charged and neutral}) traveling in roughly the same direction.  These clusters are 
referred to as ``jets".  In this analysis we examine only the charged particle component of ``jets".  Our philosophy in 
comparing the QCD Monte-Carlo models with data is to select a region where the data is very ``clean" so that ``what 
you see is what you get" ({\it almost}).  Hence, we consider only charged particles measured by the CDF central tracking 
chamber (CTC) in the region \ptcut\ and \etacut, where the track finding efficiency is high and uniform 
(estimated to be $92\%$ efficient) and we restrict ourselves to charged particle jets with transverse momentum less than 
$50\gevc$.  The data presented here are uncorrected.  Instead the theoretical Monte-Carlo models are corrected for 
the track finding efficiency by removing, on the average, $8\%$ of the charged particles.  The theory curves have an 
error ({\em statistical plus systematic}) of about $5\%$.   Thus, to within $10\%$ ``what you see is what you get".

\begin{figure}
\begin{center}
\includegraphics[scale=0.5]{Field_fig5.eps}\end{center}
\caption{\em Data on the average number of charged particles (\ptcut, \etacut) in the ``transverse'' region defined in 
Fig.~\ref{rdf_fig3} as a function of transverse momentum of the leading charged jet compared with the QCD 
Monte-Carlo predictions of {\tt HERWIG 5.9}, {\tt ISAJET} 7.32, and {\tt PYTHIA} 6.115 with their default parameters and 
with \hardcut. Each point corresponds to the \aveN\  in a $1\gevc$ bin. The solid (open) points are the 
\MB\ (JET20) data. The theory curves are corrected for the track finding efficiency and have an 
error ({\it statistical plus systematic}) of around $5\%$.
\label{rdf_fig5}}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[scale=0.5]{Field_fig6.eps}\end{center}
\caption{\em As for Fig~\ref{rdf_fig5} except that the average scalar
  $P_T$ is shown. \label{rdf_fig6}
}
\end{figure}

Charged particle ``jets" are defined as clusters of charged particles (\ptcut, \etacut) in ``circular regions" of 
\etaphi\ space with radius $R = 0.7$.   Every charged particle in the event is assigned to a ``jet", 
with the possibility that some jets might consist of just one charged particle.  The transverse momentum of a 
charged jet, $P_{T}\!({\rm chgjet})$, is the {\it scalar} \pt\ sum of the particles in the jet.  We use the 
direction of the leading charged particle jet to define correlations in azimuthal angle, $\Delta\phi$.  
The angle $\Delta\phi=\phi-\phi_{\rm chgjet\#1}$  is the relative azimuthal angle between a charged 
particle and the direction of chgjet\#1.  The``toward" region is defined by $|\Delta\phi|<60^\circ$ and \etacut, 
while the ``away" region is $|\Delta\phi|>120^\circ$ and \etacut.   The ``transverse" region is defined by 
$60^\circ<|\Delta\phi|<120^\circ$  and \etacut.  The three regions ``toward", ``transverse", and ``away" are 
shown in Fig.~\ref{rdf_fig3}.  Each region has an area in \etaphi\ space of $4\pi/3$.  As illustrated in 
Fig.~\ref{rdf_fig3}, the ``toward" region contains the leading charged particle jet, while the ``away" region, 
on the average, contains the ``away-side" jet.  The ``transverse" region is perpendicular to the plane of the 
hard $2$-to-$2$ scattering and is therefore very sensitive to the \UE. 

\begin{figure}
\begin{center}
\includegraphics[scale=0.5]{Field_fig7.eps}\end{center}
\caption{\em Data from Fig.~\ref{rdf_fig5} compared with the 
QCD Monte-Carlo predictions of {\tt ISAJET} 7.32 (default parameters and \hardcut). The predictions of {\tt ISAJET} are divided 
into two categories: charged particles that arise from the break-up of the beam and target ({\it beam-beam remnants}), 
and charged particles that result from the outgoing jets plus initial and final-state radiation 
({\it hard scattering component}). The theory curves are corrected for the track finding efficiency and have an 
error ({\it statistical plus systematic}) of around $5\%$.\label{rdf_fig7}
}
\end{figure}

Fig.~\ref{rdf_fig5} and Fig.~\ref{rdf_fig6} compare the ``transverse" \aveN\ and the ``transverse" \avePT, 
respectively, with the QCD 
Monte-Carlo predictions of {\tt HERWIG}, {\tt ISAJET}, and {\tt PYTHIA} 6.115 with their default parameters and \hardcut. The solid 
points are \MB\ data and the open points are the JET20 data. The JET20 data connect 
smoothly to the \MB\ data and allow us to study observables over the range $0.5 < P_T\!({\rm chgjet}\#1) < 50\gevc$.  
The average number of charged particles in the ``transverse" region doubles in going from \ptchj $=1.5\gevc$ to $2.5\gevc$ 
and then forms an approximately constant ``plateau" for \ptchj $>5\gevc$. If we 
suppose that the \UE\ is uniform in azimuthal angle $\phi$ and pseudo-rapidity $\eta$, the observed $2.3$ charged 
particles at \ptchj $=20\gevc$ translates to $3.8$ charged particles per unit pseudo-rapidity with \ptcut\
(multiply by $3$ to get $360^\circ$, divide by $2$ for the two units of pseudo-rapidity, 
multiply by $1.09$ to correct for the track 
finding efficiency).  We know that if we include all $p_T > 50$ MeV/c there are, on the average, about four charged 
particles per unit rapidity in a ``soft" proton-antiproton collision at 1.8 TeV \cite{Abe:1990td}.  The data in 
Fig.~\ref{rdf_fig5} imply that in the 
\UE\ of a hard scattering there are, on the average, about $3.8$ charged particles per unit rapidity with \ptcut!  
Assuming a charged particle \pt\ distribution of $e^{-2p_T}$ (see Fig.~\ref{rdf_fig25}) implies that there are 
roughly $10$ 
charged particles per unit pseudo-rapidity with $p_T >  0$ in the \UE\ (factor of e).  Since we examine 
only those charge particles with \ptcut, we cannot accurately extrapolate to low \pt, however, it is clear that 
the \UE\ has a charge particle density that is at least a factor of two larger than the four charged 
particles per unit rapidity seen in ``soft" proton-antiproton collisions at this energy.  

The \MB\ data were collected with a very ``loose" trigger.  The CDF \MB\ trigger requirement removes 
elastic scattering and most of the single and double diffraction events, but keeps essentially all
the ``hard-scattering" events.  
In comparing with the QCD Monte-Carlo models we do require that the models satisfy the CDF \MB\ trigger, 
however, for \ptchj $> 5\gevc$ essentially all the generated events satisfy the \MB\ trigger (\ie the 
\MB\ trigger is unbiased for large \pt\ ``jets").  If we had enough \MB\ events we would not need the JET20 data, 
but because of the fast fall-off of the cross section we run out of statistics at around \ptchj $=20\gevc$ (that 
is why the \MB\ data errors get large at around $20\gevc$).  The JET20 data were collected by requiring at least 
$20\gev$ of energy ({\it charged plus neutral}) in a cluster of calorimeter cells.  We do not use the calorimeter 
information, but instead look only at the charged particles measured in the CTC in the same way we do for the \MB\
data.  The JET20 data is, of course, biased for low \pt\ jets and we do not show the JET20 data below 
\ptchj\ around $20\gevc$.  At large \ptchj\ values the JET20 data becomes unbiased and, in fact, we 
know this occurs at around $20\gevc$ because it is here that it agrees with the ({\it unbiased}) \MB\ data.  

\begin{figure}
\begin{center}
\includegraphics[scale=0.5]{Field_fig10.eps}\end{center}
\caption{\em The ``hard scattering" component ({\it outgoing jets plus initial and final-state radiation}) of the number 
of charged particles (\ptcut, \etacut) in the ``transverse" region defined in Fig.~\ref{rdf_fig3} as a function of 
the transverse momentum of the leading charged jet from the QCD Monte-Carlo predictions of 
{\tt HERWIG 5.9}, {\tt ISAJET} 7.32, and {\tt PYTHIA} 6.115 with their default parameters and with \hardcut. The curves are corrected 
for the track finding efficiency and have an error ({\it statistical plus systematic}) of around $5\%$.\label{rdf_fig10}
}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[scale=0.6]{Field_fig12.eps}\end{center}
\caption{\em Data on the average number of ``transMAX" and ``transMIN" charged particles (\ptcut, \etacut) as a function 
of the transverse momentum of the leading charged jet compared with the QCD Monte-Carlo predictions of 
{\tt HERWIG 5.9}, {\tt ISAJET} 7.32, and {\tt PYTHIA} 6.115 with their default parameters and with \hardcut. The solid (open) points are 
the \MB\ (JET20) data. The theory curves are corrected for the track finding efficiency and have an error 
({\it statistical plus systematic}) of around $5\%$.\label{rdf_fig12}
}
\end{figure}

We expect the ``transverse" region to be composed predominately of particles that arise from the break-up of the 
beam and target and from initial and final-state radiation. This is clearly the case for the QCD Monte-Carlo models 
as can be seen in Figs.~7-9, where the predictions for the ``transverse" region are divided into two categories: charged 
particles that arise from the break-up of the beam and target ({\it beam-beam remnants}), and charged particles that 
result from the outgoing jets plus initial and final-state radiation ({\it hard scattering component}). For {\tt PYTHIA} the 
``beam-beam remnant" contribution includes contributions from multiple parton scattering.  It is interesting to see 
that in the QCD Monte-Carlo models it is the \BBR\ that are producing the approximately constant ``plateau".  
The contributions from initial-state and final-state radiation increase as \ptchj\ 
increases.  In fact, for {\tt ISAJET} it is the sharp rise in the initial-state radiation component that is causing the 
disagreement with the data for \ptchj $>20\gevc$.  The hard scattering component of {\tt HERWIG} and {\tt PYTHIA} 
does not rise nearly as fast as the hard scattering component of {\tt ISAJET}.  

There are two reasons why the hard scattering component of {\tt ISAJET} is different from {\tt HERWIG} and {\tt PYTHIA}.  The 
first is due to different fragmentation schemes.  {\tt ISAJET} uses independent fragmentation, which produces too many 
soft hadrons when partons begin to overlap.  The second difference arises from the way the QCD Monte-Carlos 
produce ``parton showers".  {\tt ISAJET} uses a leading-log picture in which the partons within the shower are ordered 
according to their invariant mass.  Kinematics requires that the invariant mass of daughter partons be less than the 
invariant mass of the parent.  {\tt HERWIG} and {\tt PYTHIA} modify the leading-log picture to include ``color coherence 
effects" which leads to ``angle ordering" within the parton shower.  Angle ordering produces less high \pt\ radiation 
within a parton shower which is what is seen in Fig.~\ref{rdf_fig10}.

Of course, the origin of an outgoing particle (``beam-beam remnant" or ``hard-scattering") is not an experimental 
observable.  Experimentally one cannot say where a given particle comes from.  However, we do know the origins of 
particles generated by the QCD Monte-Carlo models and Figs.~7-9 show the 
composition of the ``transverse" region as predicted by {\tt ISAJET}, {\tt HERWIG}, and {\tt PYTHIA}. 

\begin{figure}
\begin{center}
\includegraphics[scale=0.6]{Field_fig14.eps}\end{center}
\caption{\em Data on the average number of ``transMAX" charged particles (\ptcut, \etacut) as a function of the 
transverse momentum of the leading charged jet compared with the QCD Monte-Carlo predictions of 
{\tt HERWIG 5.9} (default parameters and \hardcut). The predictions of {\tt HERWIG} are divided into two categories: 
charged particles that arise from the break-up of the beam and target ({\it beam-beam remnants}), and charged particles 
that result from the outgoing jets plus initial and final-state radiation ({\it hard scattering component}).  The theory 
curves are corrected for the track finding efficiency and have an error ({\it statistical plus systematic}) of 
around $5\%$.\label{rdf_fig14}
}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[scale=0.5]{Field_fig15.eps}\end{center}
\caption{\em As for Fig~\ref{rdf_fig14} except that ``transMIN" is shown. \label{rdf_fig15}
}
\end{figure}

\subsubsection{Maximum and Minimum ``Transverse" Regions}

We now break up the ``transverse" region into two pieces.  As illustrated in Fig.~\ref{rdf_fig3}, on an event 
by event basis, we 
define ``transMAX" (``transMIN") to be the maximum (minimum) of the two ``transverse" pieces, 
$60^\circ<\Delta\phi<120^\circ$, \etacut, and $60^\circ<-\Delta\phi<120^\circ$, \etacut.  Each has an area 
in \etaphi\ space of $2\pi/3$ and what we previously referred 
to as the ``transverse" region is the sum of  ``transMAX" and ``transMIN".  One expects that ``transMAX"  will pick 
up more of the initial and final state radiation while ``transMIN" should be more sensitive to the ``beam-beam 
remnant" component of the \UE.  Furthermore, one expects that the ``beam-beam remnant" component 
will nearly cancel in the difference, ``transMAX" minus ``transMIN".  If this is true then the difference, ``transMAX" 
minus ``transMIN", would be more sensitive to the ``hard scattering" component (\ie initial and final-state radiation).  
I believe that this idea was first proposed by Bryan Webber and then implemented in a paper by Jon Pumplin \cite{Pumplin:1998ix} and then investigated in CDF by V. Tano~\cite{Tano:2001ab}. .
Fig.~\ref{rdf_fig12}  show the data on the \aveN\ for 
the``transMAX" and ``transMIN" region as a function of the \ptchj\ compared with QCD Monte-Carlo predictions of 
{\tt HERWIG}, {\tt ISAJET}, and {\tt PYTHIA} with their default parameters and \hardcut.
The data on  \avePT, show  similar behaviour. 
 Fig.~\ref{rdf_fig14} and  Fig.~\ref{rdf_fig15}, 
 show the data 
on \aveN\  for ``transMAX",  and ``transMIN", 
compared with QCD Monte-Carlo predictions of {\tt HERWIG}.   The predictions of {\tt HERWIG} are divided into two 
categories: charged particles that arise from the break-up of the beam and target ({\it beam-beam remnants}), 
and charged 
particles that result from the outgoing jets plus initial and final-state radiation ({\it hard scattering component}).  
It is 
clear from these plots that in the QCD Monte-Carlo models the ``transMAX" is more sensitive to the ``hard scattering 
component" of the \UE\ while ``transMIN" is more sensitive to the \BBR, especially 
at large values of \ptchj.   For example, for {\tt HERWIG} at \ptchj $=40\gevc$ the hard scattering 
component makes up $62\%$ of the ``transMAX" \aveN\ with $38\%$ coming from the \BBR.  On the 
other hand, the hard scattering component makes up only $42\%$ of the ``transMIN" \aveN\ with $58\%$ coming from 
the \BBR\ at \ptchj $=40\gevc$.  Taking difference between ``tansMAX" and ``transMIN" does 
not completely remove the ``beam-beam remnant" component, but reduces it to only about $32\%$ at \ptchj $=40\gevc$.

\begin{figure}
\begin{center}\includegraphics[scale=0.35]{Field_fig17.eps}\includegraphics[scale=0.35]{Field_fig18.eps}\end{center}
\caption{\em Data on the transverse momentum distribution of charged particles (\ptcut, \etacut) in the ``transverse'' region 
for \ptchj $>2\gevc$, $5\gevc$, and $30\gevc$, where chgjet\#1 is the leading charged particle jet. Each point 
corresponds to $dN_{chg}/dp_T$ and the integral of the distribution gives the average number of charged particles 
in the transverse region, $\langle N_{chg}({\rm transverse})\rangle$. The data are compared with the QCD Monte-Carlo 
model predictions of {\tt HERWIG 5.9} (default parameters and \hardcut).  The theory curves are corrected for the track 
finding efficiency and have an error ({\it statistical plus
  systematic}) of around $5\%$. Right plot is  compared with the QCD Monte-Carlo 
model predictions of {\tt ISAJET} 7.32
}
\label{rdf_fig17}
\end{figure}

\begin{figure}
\begin{center}\includegraphics[scale=0.35]{Field_fig19.eps}\includegraphics[scale=0.35]{Field_fig20.eps}\end{center}
\caption{\em The data from Fig.~\ref{rdf_fig17}
for \ptchj $> 5\gevc$ compared with the QCD Monte-Carlo model
predictions of {\tt HERWIG 5.9} (left) and {\tt ISAJET} (right)
(default parameters and \hardcut). 
\label{rdf_fig19}
}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[scale=0.35]{Field_fig21.eps}\end{center}
\caption{\em Data from Fig.~\ref{rdf_fig17}
for \ptchj $> 30\gevc$ compared with the QCD Monte-Carlo model predictions of {\tt HERWIG 5.9} 
(default parameters and \hardcut).  The theory curves are corrected for the track finding efficiency and have an 
error ({\it statistical plus systematic}) of around $5\%$. The solid curve is the total (``hard scattering" plus 
\BBR) and the dashed curve shows the contribution arising from the break-up of the beam 
particles (\BBR).\label{rdf_fig21}
}
\end{figure}

\subsubsection{The Transverse Momentum Distribution in the ``Transverse" Region}

Fig.~\ref{rdf_fig17} shows the data on the transverse momentum distribution of charged particles in the ``transverse" 
region for \ptchj $> 2\gevc$, $5\gevc$, and $30\gevc$. Each point corresponds to $dN_{\rm chg}/dp_T$ and the integral 
of the distribution gives the average number of charged particles in the ``transverse" region, \aveN, shown in 
Fig.~\ref{rdf_fig5}.  Fig.~\ref{rdf_fig5} shows only mean values, while Fig.~17 shows the 
distribution from which the mean is computed.  In Fig.~\ref{rdf_fig17} the data are compared with the QCD hard 
scattering Monte-Carlo models predictions {\tt HERWIG} and {\tt ISAJET} Since 
these distributions fall off sharply as \pt\ increases, it is essentially only the first few points at low \pt\ that 
determine the mean.  The approximately constant plateau seen in Fig.~\ref{rdf_fig5} is a result of the low \pt\ points in 
Fig.~\ref{rdf_fig17}  not changing much as \ptchj\ changes.  However, the high \pt\ points do 
increase considerably as \ptchj\ increases. This effect cannot be seen by simply examining the average number 
of  ``transverse" particles.  Fig.~\ref{rdf_fig17} shows the growth of 
the ``hard scattering component" at large \pt\ in the "transverse region" (\ie three or more hard scattering jets).

For the QCD Monte-Carlo models the \pt\ distribution in the ``transverse" region, at low values of \ptchj, is 
dominated by the ``beam-beam remnant" contribution with very little ``hard scattering" component.  This can be seen 
in Fig.~\ref{rdf_fig19}  which shows both the ``beam-beam remnant" component together 
with the total overall predictions of {\tt HERWIG} and {\tt ISAJET}, respectively, for \ptchj $>5\gevc$.  For the QCD 
Monte-Carlo models the \pt\ distribution in the ``transverse" region, at low values of \ptchj, measures directly 
the \pt\ distribution of the \BBR.  Both {\tt ISAJET} and {\tt HERWIG} have the wrong \pt\ dependence in the ``transverse" region 
due to a ``beam-beam remnant" component that falls off too rapidly as \pt\ increases.  It is, of course, understandable 
that the Monte-Carlo models might be slightly off on the parameterization of the \BBR.  This component cannot be 
calculated from perturbation theory and must be determined from data.

\begin{figure}
\begin{center}
\includegraphics[scale=0.375]{pttrack_max_min_07.eps}
\includegraphics[scale=0.375]{pttrack_max-min_07.eps}
\end{center}
\caption{\em On the left is shown data on the average {\it scalar} \pt\ sum of charged particles ($p_T > 0.4\gevc$, \etacut) within the 
maximum (MAX) and minimum (MIN) ``transverse cones" versus the transverse energy of 
the leading ({\it highest $E_T$}) ``calorimeter jet" compared with the QCD Monte-Carlo model predictions of {\tt HERWIG} and {\tt PYTHIA}.On the right is shown the difference (MAX-MIN) for each event.   
\label{rdf_fig24}
}
\end{figure}


Fig.~\ref{rdf_fig21} shows both the ``beam-beam remnant" component together with the overall prediction of {\tt HERWIG} for 
\ptchj $>30\gevc$.  Here the QCD Monte-Carlo models predict a large ``hard scattering" component 
corresponding to the production of more than two large \pt\ jets. {\tt HERWIG}, {\tt ISAJET}, and {\tt PYTHIA} all do well at 
describing the high \pt\ tail of this distribution.  However, Fig.~\ref{rdf_fig17} shows that {\tt ISAJET} produces too 
many charged particles at low \pt\ which comes from an overabundance of soft particles produced in the 
hard scattering. Fig.~\ref{rdf_fig17} shows that the large rise in the transverse charged multiplicity 
from the hard scattering component of {\tt ISAJET} seen in Fig.~\ref{rdf_fig7} comes from soft particles.  This is 
to be expected from a model that employs independent fragmentation such as {\tt ISAJET}.  Independent fragmentation 
does not differ much from color string or cluster fragmentation for the hard particles, but independent fragmentation 
produces too many soft particles.

\begin{figure}
\begin{center}
\includegraphics[scale=0.4]{Field_fig25.eps}\includegraphics[scale=0.4]{Field_fig26.eps}\end{center}
\caption{\em Data from Fig.~\ref{rdf_fig17}
for \ptchj $> 5\gevc$ compared with the QCD Monte-Carlo model predictions of a ``tuned" version of {\tt ISAJET} 7.32
(CUTJET$=12\gevc$, \hardcut).  For the ``tuned" version of {\tt ISAJET} the \pt\ distribution of the \BBR\ is generated 
according to $e^{-bp_T}$, where $b=2/(\gevc)$.  The theory curves are corrected for the track finding efficiency 
and have an error ({\it statistical plus systematic}) of around $5\%$. The solid curve is the total 
(``hard scattering" plus \BBR) and the dashed curve shows the contribution arising from the break-up of the 
beam particles (\BBR).\label{rdf_fig25}
}
\end{figure}

Note that the transverse momentum distribution of the ``beam-beam remnant" component for both {\tt HERWIG} and {\tt ISAJET}  
does not change in going from \ptchj $>5\gevc$ to \ptchj $>30\gevc$.  

\subsubsection{``Transverse Regions" Versus ``Transverse Cones"}

In a complementary CDF analysis, Valeria Tano \cite{Tano:2001ab} has studied the \UE\ in hard scattering processes 
by defining ``transverse cones" instead of ``transverse regions".   
The ``transverse cones" (with 
radius in \etaphi\ space of $R = 0.7$) are located at the same pseudo-rapidity as the leading jet but 
with azimuthal angle $\Delta\phi=+90^\circ$ and $\Delta\phi=-90^\circ$ relative to the leading ``jet".    
In the cone analysis the ``jet" is a ``calorimeter jet" ({\it charged plus neutrals}) defined using 
the standard CDF cluster algorithm.  Maximum (MAX) and minimum (MIN) 
``transverse" cones are determined, on an event-by-event basis, similar to the ``transMAX" and ``transMIN" regions 
described in Section III.  Each ``transverse cone" has an area in \etaphi\ space of $\pi R^2 = 0.49\pi$ 
(compared with $0.67\pi$).
Fig.~\ref{rdf_fig24} shows data at $1.8\tev$ on the average {\it scalar} \pt\ sum of charged particles 
($p_T > 0.4\gevc$, \etacut) within the 
MAX and MIN ``transverse cones" and the value of MAX-MIN for each event versus the transverse energy of the leading ({\it highest $E_T$}) ``calorimeter jet". (A similar analysis has been carried out at 630 GeV and a publication reporting the results from both energies is in preparation.) The data is 
compared with  QCD hard scattering Monte-Carlo models predictions from  {\tt HERWIG} and {\tt PYTHIA}.  The ``transverse cone" analysis 
covers the range $50 < E_T({\rm calorimeter\ jet}\#1) < 300\gev$, while the ``transverse region" analysis examines only 
charged particles and covers the range $0 <$\ptchj $<50\gevc$.  One cannot directly compare the two 
analysis, but if one scales the low $E_T({\rm jet}\#1)$ points in Fig.~\ref{rdf_fig24} by the ratio of 
areas $0.67\pi/0.49\pi=1.36$, one gets 
approximate agreement with the high \ptchj\ points.  Fig.~\ref{rdf_fig24} indicates that both {\tt HERWIG} and {\tt PYTHIA} correctly describe the MIN cone distributions but that {\tt PYTHIA} generates too much energy for the MAX cone. Both analyses together provide a  
good handle on the \UE\ in hard scattering processes.

\subsubsection{Tuning the Models to Fit the ``Underlying Event"}

\subsubsection{Tuning {\tt ISAJET}}
{\tt ISAJET} generates the \pt\ distribution of the \BBR\ according the power-law distribution 
$1/(1+P_T^2/b)^4$, where $b$ is chosen to give a mean \pt\ of primary particles (\ie before decay) of 
$450\mevc$. Fig.~\ref{rdf_fig19} indicates that this yields a \pt\ distribution that falls off too rapidly.  
Since one does not know a priori how to parameterize the \pt\ distribution of the \BBR, it is interesting to 
see if we can modify {\tt ISAJET} to do a better job at fitting the data.  Fig.~\ref{rdf_fig25}  
shows the \pt\ distribution of the ``beam-beam remnant" component together with the total overall predictions of 
a ``tuned" version {\tt ISAJET}.   For the ``tuned" version, {\tt ISAJET} is modified to generate the \pt\ distribution of 
the \BBR\ according to an exponential distribution of the form $e^{-bp_T}$, where $b=2/(\gevc)$. Also, for the 
tuned version of {\tt ISAJET} the parameter CUTJET is increased from its default value of $6\gevc$ to $12\gevc$ in 
order to reduce the amount of initial-state radiation.

Fig.~\ref{rdf_fig27} compares the data on the average number of charged particles in the ``transverse" region 
with the ``tuned" version of {\tt ISAJET}, where the predictions for the ``transverse� region are divided into 
two categories: charged particles that arise from the break-up of the beam and target ({\it beam-beam remnants}), 
and charged particles that result from the outgoing jets plus initial and final-state radiation 
({\it hard scattering component}).  If one compares Fig.~\ref{rdf_fig27} with Fig.~\ref{rdf_fig7} one sees that 
the tuned version of {\tt ISAJET} has a larger ``beam-beam remnant" plateau and less particles from 
the ``hard-scattering" component (\ie initial-state radiation).  The ``tuned" version of {\tt ISAJET} does a much better
job fitting the \UE\ as can be seen by comparing Fig.~\ref{rdf_fig25} with Fig.~\ref{rdf_fig17}.  However,  
Fig.~\ref{rdf_fig27} shows that there are still too many charged particles in the ``transverse" region at large
\ptchj\ values.

\begin{figure}
\begin{center}
\includegraphics[scale=0.5]{Field_fig27.eps}\end{center}
\caption{\em Data on the average number of charged particles (\ptcut, \etacut) in the ``transverse" region defined 
in Fig.~\ref{rdf_fig3} as a function of the transverse momentum of the leading charged jet compared with the 
QCD Monte-Carlo predictions of a ``tuned" version of {\tt ISAJET} 7.32 (CUTJET$=12\gevc$, \hardcut).  For the ``tuned'' version 
of {\tt ISAJET} the \pt\ distribution of the \BBR\ is generated according to $e^{-bp_T}$, where $b=2/(\gevc)$.
The predictions are divided into two categories: charged particles that arise from the break-up of the beam and 
target ({\it beam-beam remnants}), and charged particles that result from the outgoing jets plus initial and 
final-state radiation ({\it hard scattering component}).
\label{rdf_fig27}
}
\end{figure}

\subsubsection{Tuning {\tt PYTHIA}}
Now that we have constructed collider observables that are sensitive to the \UE\ we would like to tune 
the multiple parton interaction parameters of {\tt PYTHIA} to fit the data.  There are many tunable parameters.  Here we 
consider only the parameters given in Table~1.  The default values of the parameters are given in 
Table~2.  Note that the {\tt PYTHIA} default values sometimes change as the version changes \footnote{ The latest versions of PYTHIA (6.120 and higher) include additional parameters that allow one to adjust the 
energy dependence of multiple parton interactions.}

\begin{table}
\begin{center}
\caption{\em {\tt PYTHIA} multiple parton scattering parameters.}
\label{rdf_table1}
\begin{tabular}{||c|c|l||}
\hline\hline
{\bf Parameter} & {\bf Value} & {\bf Description} \\
\hline\hline
MSTP(81) &      $0$  &  Multiple-Parton Scattering off \\
\hline
           &  $1$  &  Multiple-Parton Scattering on  \\
\hline
MSTP(82) &      $1$  &  Multiple interactions assuming the same probability, \\
& & with an abrupt cut-off $P_T{\rm min}$=PARP(81) \\
\hline
           &  $3$  &  Multiple interactions assuming a varying impact parameter \\
& & and a hadronic matter overlap consistent with a \\
& & single Gaussian matter distribution, with a smooth turn-off $P_{T0}$=PARP(82) \\
\hline
           &  $4$  &  Multiple interactions assuming a varying impact parameter \\
& & and a hadronic matter overlap consistent with a \\
& & double Gaussian matter distribution (governed by PARP(83) and PARP(84)) \\
& & with a smooth turn-off $P_{T0}$=PARP(82) \\
\hline\hline
\end{tabular}
\end{center}
\end{table}

\begin{table}
\begin{center}
\caption{\em Default values for some of the multiple parton scattering parameters of {\tt PYTHIA}.}
\label{rdf_table2}
\begin{tabular}{||c|c|c||}
\hline\hline
{\bf Parameter} & {\bf {\tt PYTHIA} 6.115} & {\bf {\tt PYTHIA} 6.125} \\
\hline\hline
MSTP(81) &      $1$  &  $1$ \\
\hline
MSTP(82) &  $1$  &  $1$  \\
\hline
PARP(81) &      $1.4\gevc$  & $1.9\gevc$ \\
\hline
PARP(82) &      $1.55\gevc$  & $2.1\gevc$ \\
\hline
PARP(83) &      $0.5$  & $0.5$ \\
\hline
PARP(84) &      $0.2$  & $0.2$ \\
\hline\hline
\end{tabular}
\end{center}
\end{table}

\begin{figure}
\begin{center}
\includegraphics[scale=0.5]{Field_fig28.eps}\\
\includegraphics[scale=0.5]{Field_fig29.eps}\end{center}
\caption{\em Data on the average number of charged particles (\ptcut, \etacut) in the ``transverse'' region 
as a function of the transverse momentum of the leading charged jet compared with the QCD Monte-Carlo predictions of 
{\tt PYTHIA} 6.115 with different structure functions and different multiple parton interaction parameters and 
with \hardzero. The theory curves are corrected for the track finding efficiency and have an 
error ({\it statistical plus systematic}) of around $5\%$.\label{rdf_fig28}
}
\end{figure}

Fig.~\ref{rdf_fig28}  shows data on the average number of charged 
particles in the ``transverse" region compared with 
the QCD Monte-Carlo predictions of {\tt PYTHIA} 6.115 with different structure functions and different multiple parton 
interaction parameters and with \hardzero.  For {\tt PYTHIA} the amount of multiple parton scattering depends 
on the parton distribution functions (\ie the structure functions) and hence the number of particles produced in the 
``transverse" region (\ie the \UE) changes if one changes the structure functions.  {\tt HERWIG} and 
{\tt ISAJET} do not include multiple parton scattering and for them the number of particles in the ``transverse" is 
essentially independent of the choice of structure functions.

\begin{figure}
\begin{center}
\includegraphics[scale=0.5]{Field_fig30.eps}\end{center}
\caption{\em Data on the average number of ``transMAX" and ``transMIN" charged particles (\ptcut, \etacut) as a function 
of the transverse momentum of the leading charged jet defined compared with the QCD Monte-Carlo predictions of 
{\tt PYTHIA} 6.115 ({\it tuned version}, CTEQ4L, MSTP(82) = 4, PARP(82) = $1.4\gevc$, \hardzero). The theory curves are 
corrected for the track finding efficiency and have an error ({\it statistical plus systematic}) of around $5\%$. \label{rdf_fig30}
}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[scale=0.5]{Field_fig32.eps}\end{center}
\caption{\em Data on the average {\it scalar} \pt\ sum of ``transMAX" and ``transMIN" charged particles (\ptcut, \etacut) 
as a function of the transverse momentum of the leading charged jet defined compared with the QCD Monte-Carlo 
predictions of {\tt PYTHIA} 6.115 ({\it tuned version}, CTEQ4L, MSTP(82) = 4, PARP(82) = $1.4\gevc$, \hardzero). The theory 
curves are corrected for the track finding efficiency and have an error ({\it statistical plus systematic}) of 
around $5\%$. \label{rdf_fig32}
}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[scale=0.4]{Field_fig34.eps}\end{center}
\caption{\em Data from Fig.~\ref{rdf_fig17} compared with the QCD 
Monte-Carlo model predictions of {\tt PYTHIA} 6.115 ({\it tuned version}, CTEQ4L, MSTP(82) = 4, PARP(82) = $1.4\gevc$, \hardzero).  
The theory curves are corrected for the track finding efficiency and have an error ({\it statistical plus systematic}) 
of around $5\%$.\label{rdf_fig34}
}
\end{figure}

Figs.~\ref{rdf_fig30}-\ref{rdf_fig34} show the results of a ``tuned" version of 
{\tt PYTHIA} 6.115 with MSTP(82) = $4$ and PARP(82) = $2.4\gevc$ 
using the CTEQ4L structure functions.  One must first choose a structure function and then tune the multiple parton 
scattering parameters for that structure function.  In generating the {\tt PYTHIA} curves in 
Figs.~30-35  we have taken 
\hardzero.  In general the perturbative $2$-to-$2$ parton scattering subprocesses diverge as \pthard\ goes to 
zero.  {\tt PYTHIA} regulates these divergences using the same cut-off parameters that are used to regulate the multiple 
parton scattering cross section (see Table~1).  This allows for the possibility of using {\tt PYTHIA} to simultaneously 
describe both ``soft" and ``hard" collisions.  Most of the CDF \MB\ events are ``soft", with less than $3\%$ of the 
events having \ptchj $>5\gevc$.  There is no clear separation between ``soft" and ``hard" collisions, but 
roughly speaking \ptchj $<2\gevc$ corresponds to ``soft" \MB\ collisions and demanding \ptchj $>5\gevc$ 
assures a ``hard" collision.  Figs.~\ref{rdf_fig30}-\ref{rdf_fig32} show that the ``tuned" version of 
{\tt PYTHIA} with \hardzero\ describes fairly well the transition between ``soft" and ``hard" collisions.  
The QCD Monte-Carlo models with \hardcut\ cannot describe the data for \ptchj $<3\gevc$ 
(see Fig.~\ref{rdf_fig5} and Fig.~\ref{rdf_fig6}), whereas {\tt PYTHIA} 
with \hardzero\ seems to do a good job on the ``transverse" observables as  \ptchj\ goes to zero.
Fig.~\ref{rdf_fig34} shows the data on the transverse momentum distribution of charged particles in the ``transverse" 
region compared with the ``tuned" version of {\tt PYTHIA} 6.115 (CTEQ4L, MSTP(82) = $4$, PARP(82) = $2.4\gevc$).  The fit is 
not perfect, but it is much better than the {\tt HERWIG} prediction shown in Fig.~\ref{rdf_fig17}.  Multiple parton 
scattering produces 
more large \pt\ particles in the ``transverse" region, which is what is needed to fit the data.  
The \pt\ distribution in the ``transverse" region, at low values of \ptchj, 
for the ``tuned" version of {\tt PYTHIA} is also 
dominated by the ``beam-beam remnant" contribution as is the case for {\tt HERWIG} (see Fig.~\ref{rdf_fig19}).  However, for 
{\tt PYTHIA} the ``beam-beam remnant" component includes contributions from multiple parton scattering, which results 
in a less steep \pt\ distribution. Also, unlike {\tt ISAJET} and {\tt HERWIG}  for {\tt PYTHIA} 
the ``beam-beam remnant" component increases as \ptchj\ increases due to multiple parton scattering.

\subsubsection{Tuning {\tt HERWIG}}
The latest version of {\tt HERWIG} includes a multiple-parton scattering option \cite{Butterworth:1996zw}, which Jon Butterworth 
has tuned to fit some of the data presented here \cite{Butterworth:2001ab}.

\subsubsection{Summary}

The \UE\ in a hard scattering process is a complicated and interesting object which involves aspects of 
both non-perturbative and perturbative QCD.  Studying the ``transMAX" and ``transMIN" pieces of the ``transverse" 
region provides additional information not contained in the sum.  In the QCD Monte-Carlo models the various 
components that make up the \UE\ are weighted differently in ``transMAX" and ``transMIN" terms.  
The ``transMAX" term preferentially selects the ``hard component" of the \UE\ ({\it outgoing jets plus 
initial and final-state radiation}) while the ``transMIN" term preferentially selects the ``beam-beam remnant" 
component.  Unfortunately one cannot cleanly isolate a single component of the \UE\ since all 
components contribute to both  ``transMAX", ``transMIN", and to the difference.  However, requiring the 
Monte-Carlo models to fit both ``transMAX" and ``transMIN" (or the sum and difference) puts additional constraints 
on the way the generators model the \UE.  

{\tt ISAJET} ({\it with independent fragmentation}) produces too many ({\it soft}) particles in the \UE\ with the 
wrong dependence on \ptchj. {\tt HERWIG} and {\tt PYTHIA} modify the leading-log picture to include ``color coherence 
effects" which leads to ``angle ordering" within the parton shower and they do a better job describing the ``underlying 
event".  Both {\tt ISAJET} and {\tt HERWIG} have the too steep of a \pt\ dependence of the ``beam-beam remnant" component 
of the \UE\ and hence do not have enough \BBR\ with \ptcut.  A modified version of {\tt ISAJET} in which the \BBR\ are
generated with an exponential distribution of the form $e^{-bp_T}$, where $b=2/(\gevc)$ improves the fit to the data.
{\tt PYTHIA} with multiple parton scattering does the best job at fitting the data.

The increased activity in the \UE\ of a hard scattering over that observed in  ``soft" collisions cannot be 
explained solely by initial-state radiation. Multiple parton interactions provide a natural way of explaining the 
increased activity in the \UE\ in a hard scattering.  A hard scattering is more likely to occur when 
the ``hard cores" of the beam hadrons overlap and this is also when the probability of a multiple parton interaction is 
greatest.  For a soft grazing collision the probability of a multiple parton interaction is small. However, multiple 
parton interactions are very sensitive to the parton structure functions (PDF).  You must first decide on a particular 
PDF and then tune the multiple parton interactions to fit the data.  

One should not take the ``tuned" version of {\tt PYTHIA} 6.115 (CTEQ4L, MSTP(82) = $4$, PARP(82) = $2.4\gevc$) 
presented here too seriously.  It is encouraging that it describes fairly well the "transverse" region over the 
range $0<$\ptchj$<50\gevc$ including the transition from ``soft" to ``hard" collisions.  However, it is still not 
quite right.  For example, it does not reproduce very well the multiplicity distribution of ``soft" collisions.  
More work needs to be done in tuning the Monte-Carlo models.  In addition, more work needs to be 
done before one can say for sure that the multiple parton interaction approach is correct.   {\tt HERWIG} without multiple 
parton scattering is not that far off the data.  Maybe we simply need to change and improve the way the 
Monte-Carlo models handle the ``beam-beam remnant" component.

\section*{Acknowledgments}

Work (ERW) supported in part by the European Community's Human Potential
Programme under contract HPRN-CT-2000-00149 Physics at Colliders. 
The work was supported in part by the Director, Office of Energy

Research, Office of High Energy and Nuclear Physics of the U.S. Department of Energy under Contracts
DE--AC03--76SF00098 (IH) and DE-AC02-98CH10886 (FP).  Accordingly, the U.S.
Government retains a nonexclusive, royalty-free license to publish or
reproduce the published form of this contribution, or allow others to
do so, for U.S. Government purposes.

%\bibliography{biblifile}
%\begin{thebibliography}{99}
%%\cite{Pukhov:1999gg}
%\bibitem{Pukhov:1999gg}
%A.~Pukhov {\it et al.},
%
%%``CompHEP: A package for evaluation of Feynman diagrams and
%%integration over multi-particle phase space. User's manual for version 33,''
%hep-ph/9908288.
%%%CITATION = HEP-PH 9908288;%%
%
%%\cite{Ishikawa:1993qr}
%\bibitem{Ishikawa:1993qr}
%T.~Ishikawa, T.~Kaneko, K.~Kato, S.~Kawabata, Y.~Shimizu and H.~Tanaka
%		   [MINAMI-TATEYA group Collaboration],
%%``GRACE manual: Automatic generation of tree amplitudes in
%% Standard Models: Version 1.0,''
%KEK-92-19.
%
%%\cite{Stelzer:1994ta}
%\bibitem{Stelzer:1994ta}
%T.~Stelzer and W.~F.~Long,
%%``Automatic generation of tree level helicity amplitudes,''
%Comput.\ Phys.\ Commun.\  {\bf 81}, 357 (1994)
%[hep-ph/9401258].
%%%CITATION = HEP-PH 9401258;%%
%
%%\cite{Berends:1991ax}
%\bibitem{Berends:1991ax}
%F.~A.~Berends, H.~Kuijf, B.~Tausk and W.~T.~Giele,
%%``On the production of a W and jets at hadron colliders,''
%Nucl.\ Phys.\ B {\bf 357}, 32 (1991).
%%%CITATION = NUPHA,B357,32;%%
%
%%  %\cite{Katsanevas:1998fb}
%%  \bibitem{Katsanevas:1998fb}
%%  S.~Katsanevas and P.~Morawitz,
%%  %``SUSYGEN-2.2: A Monte Carlo event generator for MSSM sparticle    
%%  % production at e+ e- colliders,''
%%  Comput.\ Phys.\ Commun.\  {\bf 112}, 227 (1998)
%%  [hep-ph/9711417].
%%  %%CITATION = HEP-PH 9711417;%%
%
%%\cite{Caravaglios:1999yr}
%\bibitem{Caravaglios:1999yr}
%F.~Caravaglios, M.~L.~Mangano, M.~Moretti and R.~Pittau,
%%``A new approach to multi-jet calculations in hadron collisions,''
%Nucl.\ Phys.\ B {\bf 539}, 215 (1999)
%[hep-ph/9807570].
%%%CITATION = HEP-PH 9807570;%%
%
%\bibitem{alpgen} M.L. Mangano, M. Moretti, F. Piccinini, R. Pittau and
%A. Polosa, (to appear) refer to {\tt
%http://cern.ch/mlm/alpgen/alpgen.html}.
%
%%\cite{Corcella:2001bw}
%\bibitem{Corcella:2001bw}
%G.~Corcella {\it et al.},
%%``{\tt HERWIG} 6: An event generator for hadron emission reactions with  interfering gluons (including supersymmetric processes),''
%JHEP {\bf 0101}, 010 (2001)
%[hep-ph/0011363].
%%%CITATION = HEP-PH 0011363;%%
%
%%                                     Isajet
%%\cite{Baer:1999sp}
%\bibitem{Baer:1999sp}
%H.~Baer, F.~E.~Paige, S.~D.~Protopopescu and X.~Tata,
%%``ISAJET 7.48: A Monte Carlo event generator for p p, anti-p p, and  
%% e+ e- reactions,''
%hep-ph/0001086.
%%%CITATION = HEP-PH 0001086;%%
%
%%\cite{Sjostrand:2001wi}
%\bibitem{Sjostrand:2001wi}
%
%T.~Sjostrand, P.~Eden, C.~Friberg, L.~Lonnblad, G.~Miu, S.~Mrenna and E.~Norrbin,
%%``High-energy-physics event generation with {\tt PYTHIA} 6.1,''
%Comput.\ Phys.\ Commun.\  {\bf 135}, 238 (2001)
%[hep-ph/0010017].
%%%CITATION = HEP-PH 0010017;%%
%
%%\cite{Mangano:2001xp}
%\bibitem{Mangano:2001xp}
%M.~L.~Mangano, M.~Moretti and R.~Pittau,
%hep-ph/0108069.
%%%CITATION = HEP-PH 0108069;%%
%
%%\cite{Belyaev:2000wn}
%\bibitem{Belyaev:2000wn}
%A.~S.~Belyaev {\it et al.},
%%``CompHEP-{\tt PYTHIA} interface: Integrated package for the collision events  generation based on exact matrix elements,''
%hep-ph/0101232.
%%%CITATION = HEP-PH 0101232;%%
%
%\bibitem{Sjostrand:2001yu}
%T. Sj\"ostrand, L. L\"onnblad and S. Mrenna, LU~TP~01-21, 
%[hep-ph/0108264].
%\bibitem{MadCUP} K.~Cranmer, T.~Figy, W.~Quayle, D.~Rainwater, and
%D.~Zeppenfeld, refer to  {\tt 
%http://pheno.physics.wisc.edu/Software/MadCUP/}.
%
%\bibitem{comhepIMP} V.A.~Ilyin, ``Interfacing parton and hadron level
%generators'', presented at the {\it Physics and Detectors for a 90 to
%800 GeV Linear Collider Second Workshop of the Extended ECFA/DESY
%Study} in Saint Malo, France, April 2002.
%
%\bibitem{Dobbs:2001dq}
%M.~Dobbs,
%%``Phase space veto method for next-to-leading order event generators
%%  in hadronic collisions,''
%to appear in Phys.\ Rev.\ D, [hep-ph/0111234].
%%%CITATION = HEP-PH 0111234;%%
%
%\bibitem{Kersevan:2002dd}
%B.~P.~Kersevan and E.~Richter-Was,
%%``The Monte Carlo event generator AcerMC version 1.0 with interfaces
%%  to {\tt PYTHIA} 6.2 and {\tt HERWIG} 6.3,''
%hep-ph/0201302 submitted to Comp. Phys. Commuc.
%%%CITATION = HEP-PH 0201302;%%
%
%\bibitem{HEPEVT}  T.~Sj\"{o}strand et al., in ``Z physics at LEP 1'',
%  eds. G.~Altarelli, R.~Kleiss and C.~Verzegnassi, 
%	 CERN 89-08 (Geneva, 1989), Vol. 3, p. 327.
%
%%\cite{Dobbs:2001ck}
%\bibitem{Dobbs:2001ck}
%M.~Dobbs and J.~B.~Hansen,
%%``The HepMC C++ Monte Carlo event record for High Energy Physics,''
%Comput.\ Phys.\ Commun.\  {\bf 134}, 41 (2001).
%%%CITATION = CPHCB,134,41;%%
%
%%\cite{Garren:2000st}
%\bibitem{Garren:2000st}
%L.~Garren, I.~G.~Knowles, T.~Sjostrand and T.~Trippe,
%%``Monte Carlo particle numbering scheme: in Review of Particle Physics (RPP 2000),''
%Eur.\ Phys.\ J.\ C {\bf 15}, 205 (2000).
%%%CITATION = EPHJA,C15,205;%%
%
%%\cite{Plothow-Besch:1993qj}
%\bibitem{Plothow-Besch:1993qj}
%H.~Plothow-Besch,
%%``PDFLIB: A Library of all available parton density functions of the
%%nucleon, the pion and the photon and the corresponding alpha-s calculations,''
%Comput.\ Phys.\ Commun.\  {\bf 75}, 396 (1993);
%CERN Program Library Long Writeup W5051 (2000);
%refer to  {\tt http://consult.cern.ch/writeup/pdflib/}.
%%%CITATION = CPHCB,75,396;%%
%
%\bibitem{ATL-COM-PHYS-2001-013}
%B. Kersevan and E. Richter-Was, 
%{\it The $q \bar q \to W(\to \ell \nu)g^*(\to  b \bar b)$  process: matrix element
%implementation to {\tt PYTHIA} 6.1}, ATLAS Physics Communication, ATL-COM-PHYS-2001-013 (2001).
%
%\bibitem{ATL-COM-PHYS-2001-014}
%B. Kersevan and E. Richter-Was, {\it The $gg, q \bar q \to
% Z/\gamma^* (\to \ell \ell) b \bar b$ process: matrix element implementation
% to {\tt PYTHIA} 6.1}, ATLAS Physics Communication, ATL-COM-PHYS-2001-014 (2001).
%
%\bibitem{ATL-COM-PHYS-2001-025}
%B. Kersevan and E. Richter-Was, {\it The $gg, q \bar q \to
% t \bar t b \bar b$ process: matrix element implementation
% to {\tt PYTHIA} 6.1}, ATLAS Physics Communication, ATL-COM-PHYS-2001-025 (2001).
%
%\bibitem{ATL-PHYS-TDR}
%ATLAS Collaboration, {\it ATLAS Detector and Physics Performance TDR},
%CERN-LHCC/99-15 (1999). 
%
%\bibitem{Aachen1990}
%Proceedings of the Large Hadron Collider Workshop, Aachen, 1990, 
%edited by G. Jarlskog and D. Rein, CERN 90-10/ECFA 90-133.
%
%\bibitem{Kleiss1990}
%B. van Eijk and R. Kleiss, in \cite{Aachen1990}, page 184.
%
%\bibitem{EUROJET}
%B. van Eijk et al., EUROJET Monte Carlo, Nucl. Phys. {\bf B292} (1987) 1.
%
%\bibitem{DellaNegra1990}
%M. Della Negra et al., in \cite{Aachen1990}, page 509.
%
%\bibitem{Pythia57}
%T. Sjostrand, Comp. Phys Comm. {\bf 82} (1994) 74.
%
%\bibitem{Sjostrand:2001yu[was:Pythia61]}
%T. Sjostrand et al.,{\it High energy physics generation with {\tt PYTHIA} 6.1},
%eprint hep-ph/0010017, LU-TP 00-30, October 2000.
%
%\bibitem{RKEllis2000}
%J.Campbell and R. K. Ellis, Phys. Rev. {\bf D62} (2000) 114012.
%
%\bibitem{Gunion94}
%J. Gunion, Phys. Rev. Lett. {\bf 72} (1994) 199.
%
%\bibitem{ATL-PHYS-98-132}
%E. Richter-Was and M. Sapinski, ATLAS Internal Note, ATL-PHYS-98-132 (1998).
%
%\bibitem{Wbb93}
%M. L. Mangano, Nucl. Phys. {\bf B405} (1993) 536.
%
%\bibitem{Mangano:2002ab[was:Wbbgen]}
%M. L. Mangano, M. Moretti and R. Pittau,{\it  WBBGEN:  program for the generation
%of $W b \bar b$ + multijet final states in hadronic collisions},
%{\it URL:mlm.home.cern.ch/mlm/wbb/wbb.html}.
%
%\bibitem{Note-Federica}
%F. Mazzucato, {\it Studies on the Standard Model self-couplings},
%ATLAS Internal Note in preparation.
%
%\bibitem{Madgraph}
%T. Stelzer and W. F. Long, Comput. Phys. Commun. {\bf 81} (1994) 357. 
%
%\bibitem{PDFLIB}
%H. Plothow-Besch, Comput.Phys.Commun. {\bf 75} (1993) 396. 
%
%\bibitem{Kleiss1994}
%R. Kleiss and R. Pittau, Comput. Phys. Commun. {\bf 83} (1994) 141.
%
%\bibitem{nextcalibur}
%F. A. Berends, C. G. Papadopoulos and R. Pittau, 
%{\it NEXTCALIBUR: a four fermion generator for electron-positron collision},
%hep-ph/0011031. 
%
%\bibitem{fermisv}
%J. Hilgart, R. Kleiss, F. Le Dibider, Comp. Phys. Comm. {\bf 75} (1993) 191.
%
%\bibitem{excal}
% F.A. Berends, R. Kleiss, Comput. Phys. Commun. {\bf 85} (1996) 11.
%
%\bibitem{vegas}
% G.P. Lepage, J. Comput. Phys. {\bf 27} (1978) 192.
%
%\bibitem{cflows}
% H.-U. Bengtsson, Comput. Phys. Commun. {\bf 31} (1984) 323.
%
%\bibitem{ohl}
%T. Ohl, {\it Electroweak Gauge Bosons at Future Electron-Positron Colliders}, preprint
%hep-ph/9911437
%
%\bibitem{foam}
%S. Jadach, Comput. Phys. Commun.  {\bf 130} (2000)  244.
%
%\bibitem{ATL-COM-PHYS-2001-032}
%B. Kersevan and E. Richter-Was, {\it What is the $W b \bar b$ background at LHC?...},
%ATLAS Internal Note, ATL-COM-PHYS-2001-032 (2001).
%
%% M.~Bengts%son, T.~Sj\"ostrand, 
%%and M.~van~Zijl, Z. Phys. {\bf C32}, 67 (1986); T.~Sj\"ostrand and M.~van~Zijl%, Phys. Rev. {\bf 
%%D36}, 2019 (1987).
%
%\bibitem{Balazs:1997xd}
%C.~Balazs and C.~P.~Yuan,
%%``Soft gluon effects on lepton pairs at hadron colliders,''
%Phys.\ Rev.\ D {\bf 56}, 5558 (1997)
%[arXiv:hep-ph/9704258].
%%%CITATION = HEP-PH 9704258;%%
%
%\bibitem{Djouadi:2000gu}
%A.~Djouadi {\it et al.},
%%``The Higgs working group: Summary report,''
%arXiv:hep-ph/0002258.
%%%CITATION = HEP-PH 0002258;%%
%
%\bibitem{Balazs:2000wv}
%C.~Balazs and C.~P.~Yuan,
%%``Higgs boson production at the LHC with soft gluon effects,''
%Phys.\ Lett.\ B {\bf 478}, 192 (2000)
%[arXiv:hep-ph/0001103].
%%%CITATION = HEP-PH 0001103;%%
%
%\bibitem{Balazs:2000sz}
%C.~Balazs, J.~Huston and I.~Puljak,
%%``Higgs production: A comparison of parton showers and resummation,''
%Phys.\ Rev.\ D {\bf 63}, 014021 (2001)
%[arXiv:hep-ph/0002032].
%%%CITATION = HEP-PH 0002032;%%
%
%\bibitem{factor}
%G.~T.~Bodwin, S.~J.~Brodsky and G.~P.~Lepage,
%%``Factorization Of The Drell-Yan Cross-Section,''
%ANL-HEP-CP-85-32
%
%\bibitem{ISR}
%M.~Bengtsson, T.~Sjostrand and M.~van Zijl,
%%``Initial State Radiation Effects On W And Jet Production,''
%Z.\ Phys.\ C {\bf 32}, 67 (1986).
%
%\bibitem{DGLAP1}
%Gribov, V. N. and Lipatov, L. N. Yad. Fiz {\bf 15}, (1972).
%
%\bibitem{DGLAP2}
%Altarelli, Guido and Parisi, G. Nucl. Phys.,
% {\bf B126},298 (1977).
%
%\bibitem{DGLAP3} Dokshitzer, Yuri L. Sov. Phys. JETP",
%   {\bf 46},641 (1977).
%
%\bibitem{odorico} Odorico, R. Nucl. Phys.,
%     { \bf B172},157 (1980).
%
%\bibitem{kbook}Ellis, R. K. and Stirling, W. James and Webber, B. R.,
%  "QCD and collider physics. Cambridge
%  Monogr. Part. Phys. Nucl. Phys. Cosmol. (1996).
%
%\bibitem{Sjostrand:1985xi,[was:back1]} T.~Sj\"ostrand, Phys. Lett. {\bf 157B}, 321 (1985).
%
%\bibitem{back2}
%G.~Marchesini and B.~R.~Webber,
%%``Simulation Of QCD Jets Including Soft Gluon Interference,''
%Nucl.\ Phys.\ B {\bf 238}, 1 (1984).
%
%\bibitem{Catani:1991rr}
%S.~Catani, B.~R.~Webber and G.~Marchesini,
%%``QCD Coherent Branching And Semiinclusive Processes At Large X,''
%Nucl.\ Phys.\ B {\bf 349}, 635 (1991).
%
%\bibitem{DDT} Dokshitzer, Yuri L. and Diakonov, Dmitri and Troian, S. I.
%  Phys. Lett.,
%{\bf B79}, 269 (1978).
%
%\bibitem{Parisi}
%G.~Parisi and R.~Petronzio,
%%``Small Transverse Momentum Distributions In Hard Processes,''
%Nucl.\ Phys.\ B {\bf 154} (1979) 427.
%
%\bibitem{logs}
%S.~D.~Ellis, N.~Fleishon and W.~J.~Stirling,
%%``Logarithmic Approximations, Quark Form-Factors And Quantum Chromodynamics,''
%Phys.\ Rev.\ D {\bf 24}, 1386 (1981).
%
%\bibitem{backtoback}
%J.~C.~Collins and D.~E.~Soper,
%%``Back To Back Jets In QCD: Comparison With Experiment,''
%Phys.\ Rev.\ Lett.\  {\bf 48} (1982) 655.
%
%\bibitem{cssWZ}
%J.~C.~Collins, D.~E.~Soper and G.~Sterman,
%%``Transverse Momentum Distribution In Drell-Yan Pair And W And Z Boson Production,''
%Nucl.\ Phys.\ B {\bf 250}, 199 (1985).
%
%\bibitem{AEGM}
%G.~Altarelli, R.~K.~Ellis, M.~Greco and G.~Martinelli,
%%``Vector Boson Production At Colliders: A Theoretical Reappraisal,''
%Nucl.\ Phys.\ B {\bf 246}, 12 (1984).
%
%\bibitem{AEM85}
%G.~Altarelli, R.~K.~Ellis and G.~Martinelli,
%%``Vector Boson Production At Present And Future Colliders,''
%Z.\ Phys.\ C {\bf 27}, 617 (1985).
%
%\bibitem{Arnold-Kauffman} P~B.~Arnold and R.~P.~Kauffman,
%%``W And Z Production At Next-To-Leading Order: From Large Q(T) To Small,''
%Nucl.\ Phys.\ B {\bf 349}, 381 (1991).
%
%\bibitem{resum2}
%R.~K.~Ellis and S.~Veseli,
%Nucl.\ Phys.\ B {\bf 511}, 649 (1998)
%
%\bibitem{field_prd} {\it Charged Jet Evolution and the Underlying Event in Proton-Antiproton Collisions at 1.8 TeV}, The CDF 
%Collaboration (T. Affolder et al.), FERMILAB-PUB-01/211-E ({\it to appear in Phys. Rev. D.}).
%
%\bibitem{field_dpf} {\it The Underlying Event in Large Transverse Momentum Charged Jet and Z-boson Production at 1.8 TeV}, talk 
%presented by Rick Field at DPF2000, Columbus, OH, August 11, 2000.
%
%\bibitem{tano} {\it A Comparison of the Underlying Event in Jet and \MB\ Events}, talk presented by Joey Huston at DPF2000, 
%Columbus, OH, August 11, 2000.  {\it The Underlying Event in Jet and Minimum Bias Events at the Tevatron}, talk 
%presented by Valeria Tano at ISMD2001, Datong, China, September 1-7, 2001.
%
%\bibitem{Marchesini:1988cf,Knowles:1988vs,Catani:1991rr[was:herwig]} G.~Marchesini and B.~R.~Webber, Nucl. Phys {\bf B310}, 461 (1988); I.~G.~ Knowles, 
%Nucl. Phys. {\bf B310}, 571 (1988); S.~Catani, G.~Marchesini, and B.~R.~Webber, Nucl. Phys. {\bf B349}, 635 (1991).
%
%\bibitem{Paige:1986vk,[was:isajet]} F.~Paige and S.~Protopopescu, BNL Report, BNL38034, 1986 ({\it unpublished}),
%version 7.32.
%
%\bibitem{CDF3} F.~Abe \etal, Phys. Rev., {\bf D41}, 2230 (1990).
%
%\bibitem{pumplin} {\it Hard Underlying Event Corrections to Inclusive Jet Cross-Sections}, Jon Pumplin, Phys. Rev. D57, 5787 (1998).
%
%\bibitem{butterworth1} {\it Multiparton Interactions in Photoproduction at HERA}, J.~M.~Butterworth, J.~R.~Forshaw, 
%and M.~H.~Seymour, Zeit. f\"ur Phys. {\bf C72}, 637 (1996).
%
%\bibitem{butterworth2} Talk presented by Jon Butterworth at Snowmass
%  2001 ({\it unpublished}).
%
%\end{thebibliography}

%\end{document}

%
