source: trunk/documents/UserDoc/DocBookUsersGuides/PhysicsReferenceManual/latex/electromagnetic/standard/MuPgen.tex @ 1222

\section[Gamma Conversion into a Muon - Anti-mu Pair]
 {Gamma Conversion into a Muon - Anti-muon Pair}

The class {\tt G4GammaConversionToMuons} simulates the process of gamma 
conversion into muon pairs.  Given the photon energy and $Z$ and $A$ of the 
material in which the photon converts, the probability for the conversions 
to take place is calculated according to a parameterized total cross section.
Next, the sharing of the photon energy between the $\mu^+$ and $\mu^-$ is
determined.  Finally, the directions of the muons are generated.  Details of
the implementation are given below and can be also found in\,\cite{MuPgen}.

\subsection{Cross Section and Energy Sharing}
In the field of the nucleus, muon pair production on atomic electrons,
$\gamma+e\to e+\mu^+ +\mu^-$, has a threshold of
$2m_\mu(m_\mu+m_e)/m_e\approx 43.9\;{\rm GeV}$ .
%hbu check in MuonBkg.C  in fact 43.905 GeV
Up to several hundred GeV this process has a much lower cross section than 
the corresponding process on the nucleus.  At higher energies, the cross 
section on atomic electrons represents a correction of $\sim 1/Z$ to the 
total cross section.

For the approximately elastic scattering considered here, momentum, but no 
energy, is transferred to the nucleon.  The photon energy is fully shared 
by the two muons according to
\begin{equation}
E_\gamma = E_\mu^+ + E_\mu^-
\end{equation}
or in terms of energy fractions
\[ x_+ = \frac{E_\mu^+}{E_\gamma}, \qquad x_- = \frac{E_\mu^-}{E_\gamma}, \qquad x_+ + x_- = 1\;.
\]
The differential cross section for electromagnetic pair creation of muons in
terms of the energy fractions of the muons is
\begin{equation}
\frac{d\sigma}{d x_+} = 4 \, \alpha \, Z^2 \, r_c^2
\left(1-\frac43\,x_+x_-\right)\log(W)\;,
\label{eq:dSigxPlus}
\end{equation}
where $Z$ is the charge of the nucleus, $r_c$ is the classical radius of the
particles which are pair produced (here muons) and
\begin{equation}
W = W_\infty \;
\frac{1+(D_n\sqrt{e}-2)\,\delta\,/m_\mu}{1+B \, Z^{-1/3}\, \sqrt e \,\delta\,/m_e}
\label{eq:W}
\end{equation}
where
\[
W_\infty = \frac{B \, Z^{-1/3}}{D_n} \, \frac{m_\mu}{m_e} \qquad
\delta = \frac{m_\mu^2}{2\,E_\gamma\,x_+x_-} \qquad 
\sqrt{e}=1.6487\dots .
\]
\begin{eqnarray}
\mbox{For hydrogen} \qquad B=202.4 \qquad D_n=&1.49 & \nonumber\\
\mbox{and for all other nuclei} \qquad B=183 \qquad D_n=&1.54 \, A^{0.27}. & 
\end{eqnarray}
These formulae are obtained from the differential cross section for muon 
bremsstrahlung \cite{Kelner:1995hu} by means of crossing relations.  The
formulae take into account the screening of the field of the nucleus by the
atomic electrons in the Thomas-Fermi model, as well as the finite size of 
the nucleus, which is essential for the problem under consideration.
%
The above parameterization gives good results for $E_\gamma \gg m_\mu$.  The
fact that it is approximate close to threshold is of little practical 
importance.  Close to threshold, the cross section is small and the few low 
energy muons produced will not travel very far.  The cross section 
calculated from Eq.\,(\ref{eq:dSigxPlus}) is positive for 
$E_\gamma > 4 m_\mu$ and
\begin{equation}
x_{\rm min} \leq x \leq x_{\rm max} \quad {\rm with} \quad
x_{\rm min} = \frac12 - \sqrt{\frac{1}{4}-\frac{m_\mu}{E_\gamma}} \qquad
x_{\rm max} = \frac12 + \sqrt{\frac{1}{4}-\frac{m_\mu}{E_\gamma}}\;,
\end{equation}
except for very asymmetric pair-production, close to threshold, which can 
easily be taken care of by explicitly setting $\sigma = 0$ whenever 
$\sigma < 0$.

Note that the differential cross section is symmetric in $x_+$ and $x_-$ and
that
\[ x_+ x_- = x - x^2 \] 
where $x$ stands for either $x_+$ or $x_-$.  By defining a constant
\begin{equation}
\sigma_0 = 4 \, \alpha \, Z^2 \, r_c^2 \log(W_\infty)\;
\label{eq:sigma0}
\end{equation}
the differential cross section Eq.\,(\ref{eq:dSigxPlus}) can be rewritten 
as a normalized and symmetric as function of $x$: 
\begin{equation}
\frac{1}{\sigma_0} \, \frac{d\sigma}{dx} = \left[ 1-\frac43 \,(x - x^2) \right]
\, \frac{\log W}{\log W_\infty}\;.
\label{eq:dSigdx}
\end{equation}
This is shown in Fig.\,\ref{plot:dsigdx} for several elements and a wide 
range of photon energies.  The asymptotic differential cross section for 
$E_\gamma \rightarrow \infty$
\[  \frac{1}{\sigma_0} \, \frac{d\sigma_\infty}{dx} = 1-\frac43 \,(x - x^2) \]
is also shown.
\begin{figure}[htpb]
\center\includegraphics[scale=.65]{electromagnetic/standard/MuPgen/dsigdx.eps}
\caption{Normalized differential cross section for pair production as a 
function of $x$, the energy fraction of the photon energy carried by one of 
the leptons in the pair.  The function is shown for three different 
elements, hydrogen, beryllium and lead, and for a wide range of photon 
energies.}
\label{plot:dsigdx}
\end{figure}

\subsection{Parameterization of the Total Cross Section}
The total cross section is obtained by integration of the differential 
cross section Eq.\,(\ref{eq:dSigxPlus}), that is
\begin{equation}
\sigma_{\rm tot}(E_\gamma) = \int_{x_{\rm min}}^{x_{\rm max}} \frac{d\sigma}{d x_+} \, d x_+
= 4 \, \alpha \, Z^2 \, r_c^2 \, \int_{x_{\rm min}}^{x_{\rm max}}\left(1-\frac43\,x_+x_-\right)\log(W) \, d x_+ \;.
\label{eq:sigmatot}
\end{equation}
$W$ is a function of ($x_+, E_\gamma$) and ($Z, A$) of the element 
(see Eq.\,(\ref{eq:W})).  Numerical values of $W$ are given in 
Table\,\ref{tab:W}.

\begin{table}[htbp]\center
\caption{Numerical values of $W$ for $x_+=0.5$ for different elements.}
\label{tab:W}\vskip 1mm
\begin{tabular}{|c|c|c|c|c|} \hline
$E_\gamma$ & W for H & W for Be & W for Cu & W for Pb \\
 GeV      & & & & \\ \hline
     1    & 2.11  & 1.594  & 1.3505 & 5.212 \\
    10    & 19.4  & 10.85  & 6.803 & 43.53  \\
   100    & 191.5 & 102.3  & 60.10 & 332.7  \\
  1000    & 1803  & 919.3  & 493.3 & 1476.1 \\
 10000    & 11427 & 4671   & 1824  & 1028.1 \\
 $\infty$ & 28087 & 8549   & 2607  & 1339.8 \\ \hline
\end{tabular}
\end{table}
%
Values of the total cross section obtained by numerical integration are 
listed in Table\,\ref{tab:sigmatot} for four different elements.  Units are 
in $\mu{\rm barn}\,$, where $1\,\mu{\rm barn} = 10^{-34}\,{\rm m}^2\,$.
%
\begin{table}[htbp]\center
\caption{Numerical values for the total cross section}
\label{tab:sigmatot}\vskip 1mm
\begin{tabular}{|c|c|c|c|c|} \hline
$E_\gamma$ & $\sigma_{\rm tot}$, H & $\sigma_{\rm tot}$, Be & $\sigma_{\rm tot}$, Cu & $\sigma_{\rm tot}$, Pb \\
 GeV      & $\mu{\rm barn}\,$       & $\mu{\rm barn}\,$       & $\mu{\rm barn}\,$       & $\mu{\rm barn}\,$      \\ \hline
     1    & 0.01559 & 0.1515 & 5.047 & 30.22 \\
    10    & 0.09720 & 1.209  & 49.56 & 334.6 \\
   100    & 0.1921  & 2.660  & 121.7 & 886.4 \\
  1000    & 0.2873  & 4.155  & 197.6 & 1476  \\
 10000    & 0.3715  & 5.392  & 253.7 & 1880  \\
 $\infty$ & 0.4319  & 6.108  & 279.0 & 2042  \\ \hline
\end{tabular}
\end{table}
%
\begin{figure}[htpb]
\center\includegraphics[scale=.7]{electromagnetic/standard/MuPgen/SigTot.eps}
\caption{Total cross section for the Bethe-Heitler process 
$\gamma \rightarrow \mu^+\mu^-$ as a function of the photon energy 
$E_\gamma$ in hydrogen and lead, normalized to the asymptotic cross section 
$\sigma_\infty$.}
\label{plot:SigTot}
\end{figure}

\noindent
Well above threshold, the total cross section rises about linearly in 
$\log(E_\gamma)$ with the slope
\begin{equation}
W_M = \frac{1}{4\, D_n \, \sqrt e \,m_\mu}
\end{equation}
until it saturates due to screening at $\sigma_\infty$.  
Fig.\,\ref{plot:SigTot} shows the normalized cross section where
\begin{equation}
\sigma_\infty = \frac79 \, \sigma_0 \qquad {\rm and} \qquad \sigma_0 = 4 \, \alpha \, Z^2 \, r_c^2 \, \log(W_\infty)\;.
\end{equation}
Numerical values of $W_M$ are listed in Table\,\ref{tab:WM}.

\begin{table}[htbp]\center
\caption{Numerical values of $W_M$.}
\label{tab:WM}\vskip 1mm
\begin{tabular}{|c|c|} \hline
Element & $W_M$ \\
        & 1/GeV \\ \hline
H       & 0.963169 \\
Be      & 0.514712 \\
Cu      & 0.303763 \\
Pb      & 0.220771 \\
\hline
\end{tabular}
\end{table}


The total cross section can be parameterized as
\begin{equation}
\sigma_{\rm par} = \frac{28 \, \alpha \, Z^2 \, r_c^2}{9}  \; \log(1 + W_M C_f E_g)\;,
\label{eq:sigpar}
\end{equation}
with
\begin{equation}
E_g = \left(1-\frac{4 m_\mu}{E_\gamma}\right)^{t}
\left(W_{\rm sat}^{s} + E_\gamma^{s} \right)^{1/s}\;.
\end{equation}
and
\[
W_{\rm sat} = \frac{W_\infty}{W_M} = B \, Z^{-1/3} \, \frac{4\,\sqrt e \,m_\mu^2 }{m_e}\;.
\]
The threshold behavior in the cross section was found to be well 
approximated by $t = 1.479 + 0.00799 D_n $ and the saturation by
$ s = -0.88 $.  The agreement at lower energies is improved using an 
empirical correction factor, applied to the slope $W_M$, of the form
\[
C_f = \left[ 1 + 0.04 \log \left(1+\frac{E_c}{E_\gamma}\right)\right]\;,
\]
where
\[ E_c = \left[ -18.+\frac{4347.}{B \, Z^{-1/3}}\right] \;{\rm GeV}\;.
\]
A comparison of the parameterized cross section with the numerical 
integration of the exact cross section shows that the accuracy of the 
parametrization is better than 2\%, as seen in Fig.\,\ref{plot:SigApRat}.
\begin{figure}[htpb]
\center\includegraphics[scale=.8333]{electromagnetic/standard/MuPgen/SigApRat.eps}
\caption{Ratio of numerically integrated and parametrized total cross 
sections as a function of $E_\gamma$ for hydrogen, beryllium, copper and 
lead.}
\label{plot:SigApRat}
\end{figure}

\subsection{Multi-differential Cross Section and Angular Variables}
The angular distributions are based on the multi-differential cross section
for lepton pair production in the field of the Coulomb center
\[
\frac{d\sigma}{dx_+ \, du_+ \, du_-\,d\varphi} = \frac{4\,Z^2\alpha^3}{\pi}\,\frac{m_\mu^2}{q^4}\,u_+\,u_-
\]
\[
\left\{ \frac{u_+^2+u_-^2}{(1+u_+^2)\,(1+u_-^2)} -2x_+x_- \right.
\]
\begin{equation}
\left.
\left[\frac{u_+^2}{(1+u_+^2)^2}+\frac{u_-^2}{(1+u_-^2)^2}\right]
-\frac{2u_+u_-(1-2x_+x_-)\,\cos\varphi}{(1+u_+^2)\,(1+u_-^2)}
\right\} \,. \\
\label{eq:MultiDiff}
\end{equation}
Here
\begin{equation}
u_\pm = \gamma_\pm \theta_\pm \quad , \qquad \gamma_\pm = \frac{E_\mu^\pm}{m_\mu}
\quad,\qquad q^2=q_{\parallel}^2+q_{\perp}^2\quad, %new
\label{eq:uthetagamma}
\end{equation}
where
\begin{eqnarray}
  q_{\parallel}^2=q_{\min}^2\,(1+x_-u_+^2+x_+u_-^2)^2\,, \nonumber \\
  q_{\perp}^2=m_\mu^2\left[(u_+-u_-)^2+2\,u_+u_-(1-\cos\varphi)
  \right]\,. \,
\label{eq:q2}
\end{eqnarray}
$q^2$ is the square of the momentum ${\bf q}$ transferred to the target 
and $q_{\parallel}^2$ and $q_{\perp}^2$ are the squares of the components 
of the vector ${\bf q}$, which are parallel and perpendicular to the 
initial photon momentum, respectively.
The minimum momentum transfer is 
$q_{\min}=m_\mu^2/(2E_\gamma \, x_+x_-)$.\\

The muon vectors have the components
\begin{equation}
  \begin{array}{rcl}\displaystyle
  {\bf p}_+&=&p_+\,(\;\;\;\sin\theta_+\cos(\varphi_0+\varphi/2)\,,\,\;\;\;\sin\theta_+\sin(\varphi_0+\varphi/2)\,,\,\cos\theta_+)\,,\\ 
  {\bf p}_-&=&p_-\,(     -\sin\theta_-\cos(\varphi_0-\varphi/2)\,,\,     -\sin\theta_-\sin(\varphi_0-\varphi/2)\,,\,\cos\theta_-)\,,
  \end{array}
\label{eq:pvec}
\end{equation}
where $p_{\pm}=\sqrt{E_{\pm}^2-m_\mu^2}$.
The initial photon direction is taken as the $z$-axis.
The cross section of Eq.\,(\ref{eq:MultiDiff}) does not depend on 
$\varphi_0$.  Because of azimuthal symmetry, $\varphi_0$ can simply be 
sampled at random in the interval $(0,\,2\,\pi)$.

Eq.\,(\ref{eq:MultiDiff}) is too complicated for efficient Monte Carlo 
generation.  To simplify, the cross section is rewritten to be symmetric
in $u_+$, $u_-$ using a new variable $u$ and small parameters $\xi,\beta$, 
where $u_\pm=u \pm \xi/2$ and $\beta = u \,\varphi$.  When higher powers 
in small parameters are dropped, the differential cross section in terms 
of $u,\xi,\beta$ becomes
\begin{eqnarray}
\label{mupgen.a}
 \frac{d\sigma}{dx_+ \, d\xi \, d\beta\, u du} & = & \frac{4\,Z^2\alpha^3}{\pi}
\frac{m_\mu^2}{\left(q_{\parallel}^2+m_\mu^2(\xi^2+\beta^2)\right)^2} \\
& &   \left\{\xi^2\left[\frac1{(1+u^2)^2}-2\,x_+x_-\,\frac{(1-u^2)^2}{(1+u^2)^4}\right]+
  \frac{\beta^2(1-2x_+x_-)}{(1+u^2)^2}\right\}\,, \nonumber
\label{eq:MultiDiff2}
\end{eqnarray}
where, in this approximation,
$$
q_{\parallel}^2=q_{\min}^2\,(1+u^2)^2\,.\,
$$
For Monte Carlo generation, it is convenient to replace ($\xi,\beta$) by 
the polar coordinates ($\rho,\psi$) with $\xi=\rho\,\cos\psi$ and 
$\beta=\rho\,\sin\psi$.  Integrating Eq.~\ref{mupgen.a} over $\psi$ and 
using symbolically $du^2$ where $du^2 = 2 u \, du$ yields
\begin{equation}
\label{mupgen.b}
\frac{d\sigma}{dx_+\,d\rho\,du^2} =\frac{4Z^2\alpha^3}{m_\mu^2}\,\frac{\rho^3}{(q_{\parallel}^2/m_\mu^2+\rho^2)^2}
\,\,\left\{\frac{1-x_+x_-}{(1+u^2)^2}-\frac{x_+x_-(1-u^2)^2}{(1+u^2)^4}\right\}.\,
\end{equation}
Integration with logarithmic accuracy over $\rho$ gives
\begin{equation}\label{q4}
\int\!\frac{\rho^3\,d\rho}{(q_{\parallel}^2/m_\mu^2+\rho^2)^2}
\approx \int\limits_{q_{\parallel}/m_\mu}^1\!\frac{d\rho}{\rho}
=\log\left(\frac{m_\mu}{q_{\parallel}}\right)\,.
\end{equation}
Within the logarithmic accuracy, $\log(m_\mu/q_{\parallel})$ can be 
replaced by $\log(m_\mu/q_{\min})$, so that
\begin{equation}
\frac{d\sigma}{dx_+\,du^2}=\frac{4\,Z^2\alpha^3}{m_\mu^2}\,
\left\{\frac{1-x_+x_-}{(1+u^2)^2}-\frac{x_+x_-(1-u^2)^2}{(1+u^2)^4}\right\}\,
\log\left(\frac{m_\mu}{q_{\min}}\right)\,.
\end{equation}
Making the substitution $u^2 = 1/t -1$, $du^2 = -dt \, / t^2$ gives
\begin{equation}
\frac{d\sigma}{dx_+\,dt}=\frac{4\,Z^2\alpha^3}{m_\mu^2}\,
\left[1-2\,x_+x_-+4\,x_+x_-t\,(1-t)\right]\,
\log\left(\frac{m_\mu}{q_{\min}}\right) . \,
\label{eq:sigmadxdt}
\end{equation}
Atomic screening and the finite nuclear radius may be taken into account by 
multiplying the differential cross section determined by 
Eq.\,(\ref{eq:MultiDiff2}) with the factor
\begin{equation}\label{q5}
  \left(F_a(q)-F_n(q)\,\right)^2\,,
\end{equation}
where $F_a$ and $F_n$ are atomic and nuclear form factors.
Please note that after integrating Eq.~\ref{mupgen.b} over $\rho$, the 
$q$-dependence is lost.

\subsection{Procedure for the Generation of Muon - Anti-muon Pairs}

Given the photon energy $E_\gamma$ and $Z$ and $A$ of the material in which 
the $\gamma$ converts, the probability for the conversions to take place is 
calculated according to the parametrized total cross section 
Eq.\,(\ref{eq:sigpar}).  The next step, determining how the photon energy
is shared between the $\mu^+$ and $\mu^-$, is done by generating $x_+$ 
according to Eq.\,(\ref{eq:dSigxPlus}).  The directions of the muons are 
then generated via the auxilliary variables $t,\,\rho,\,\psi$.  In more
detail, the final state is generated by the following five steps, in which
$R_{1,2,3,4,...}$ are random numbers with a flat distribution in the 
interval [0,1]. The generation proceeds as follows.
\\ \\
{\bf 1)} Sampling of the positive muon energy $E_\mu^+ = x_+ \, E_\gamma$. \\
This is done using the rejection technique.
$x_+$ is first sampled from a flat distribution within kinematic limits 
using
\[ x_+ = x_{\rm min} + R_1 (x_{\rm max} - x_{\rm min}) \]
and then brought to the shape of Eq.\,(\ref{eq:dSigxPlus}) by keeping all 
$x_+$ which satisfy 
\[ \left(1-\frac43\,x_+x_-\right)\frac{\log(W)}{\log(W_{\rm max})} < R_2 \,. \]
Here $W_{\rm max}= W(x_+=1/2)$ is the maximum value of $W$, obtained for 
symmetric pair production at $x_+=1/2$.  About 60\% of the events are kept 
in this step.  Results of a Monte Carlo generation of $x_+$ are illustrated 
in Fig.\,\ref{plot:xPlusGen}.  The shape of the histograms agrees with the 
differential cross section illustrated in Fig.\,\ref{plot:dsigdx}.
\begin{figure}[htpb]
\center\includegraphics[scale=.7]{electromagnetic/standard/MuPgen/xPlusGen.eps}
\caption{Histogram of generated $x_+$ distributions for beryllium at three 
different photon energies.  The total number of entries at each energy is 
$10^6$.}
\label{plot:xPlusGen}
\end{figure}
\\ \\
{\bf 2)} Generate $t ( = \frac{1}{\gamma^2 \theta^2 + 1} )$ . \\
The distribution in $t$ is obtained from Eq.(\ref{eq:sigmadxdt}) as
\begin{equation}\label{t}
f_1(t)\,dt=\frac{1-2\,x_+x_-+4\,x_+x_-t\,(1-t)}
{1+C_1/t^2}\,dt\,,\quad 0<t\le 1\,.
\end{equation}
with form factors taken into account by
\begin{equation}\label{C_1}
C_1=
\frac{(0.35\,A^{0.27})^2}{x_+x_-\,E_\gamma/m_\mu }\,.
\end{equation}
In the interval considered, the function $f_1(t)$ will always be bounded 
from above by
\[
\max [f_1(t)]=\frac{1-x_+x_-}{1+C_1}\;.\,
\]
For small $x_+$ and large $E_\gamma$, $f_1(t)$ approaches unity, as shown
in Fig.\,\ref{f1t.eps}.
\begin{figure}[htpb]
\includegraphics[scale=.65]{electromagnetic/standard/MuPgen/f1t_10.eps}\hfill\includegraphics[scale=.65]{electromagnetic/standard/MuPgen/f1t_1000.eps}
  \caption{The function $f_1(t)$ at $E_\gamma = 10\,{\rm GeV}$ (left) and
  $E_\gamma = 1\,{\rm TeV}$ (right) in beryllium for different values of 
  $x_+$.}
\label{f1t.eps}
\end{figure}

\begin{figure}[htpb]
% \begin{minipage}[t]{7.5cm}
% \includegraphics[scale=.62]{electromagnetic/standard/MuPgen/f1tgen.eps}
  \center\includegraphics[scale=.8]{electromagnetic/standard/MuPgen/f1tgen.eps}
 \caption{Histograms of generated $t$ distributions for 
  $E_\gamma = 10\,{\rm GeV}$ (solid line) and 
  $E_\gamma = 100\,{\rm GeV}$ (dashed line) with $10^6$ events each.}
 \label{f1t_gen.eps}
% \end{minipage}
\end{figure}
% \hfill
\begin{figure}[htpb]
%  \begin{minipage}[t]{7.5cm}
% { \center\includegraphics[scale=.6]{electromagnetic/standard/MuPgen/PsiGen.eps}
  \center\includegraphics[scale=.8]{electromagnetic/standard/MuPgen/PsiGen.eps}
  \caption{Histograms of generated $\psi$ distributions for beryllium at 
  four different photon energies.}
  \label{plot:PsiGen.eps}
% }
% \end{minipage}
\end{figure}

\noindent
The Monte Carlo generation is done using the rejection technique.  About 
70\% of the generated numbers are kept in this step.  Generated 
$t$-distributions are shown in Fig.\,\ref{f1t_gen.eps}.
\\ \\
{\bf 3)} Generate $\psi$ by the rejection technique using $t$ generated in 
the previous step for the frequency distribution
\begin{equation}\label{q3}
f_2(\psi) =\Big[1-2\,x_+x_-+4\,x_+x_-t\,(1-t)\,(1+\cos(2\psi))\Big]\;, \qquad 0\le\psi\le2\pi\,.
\end{equation}
The maximum of $f_2(\psi)$ is
\begin{equation}
\max [f_2(\psi)]=1-2\,x_+x_-\left[1-4\,t\,(1-t)\right]\,.\,
\end{equation}
Generated distributions in $\psi$ are shown in Fig.\,\ref{plot:PsiGen.eps}.
\\ \\
{\bf 4)} Generate $\rho$. \\
%old The frequency distribution to generate is
%old \begin{equation}
%old f_2(\rho^2) = \frac{1}{\rho^2 + \kappa^2} \qquad  0 \leq \rho^2 \leq 1
%old \end{equation}
%old where $\kappa^2=1/W$ with the value of $W$ calculated in step 1. % according to Eq.\,\ref{eq:W}.
%old The distribution can be generated by direct transformation.
%old By integration
%old \[
%old F_2(x) = \int_0^x f_2(\rho^2) \, d \rho^2 = \log\left(1 + \frac{x}{\kappa^2} \right)
%old \]
%old with the inverse
%old \[
%old F_2^{-1}(x) =  \kappa^2 \, \left(e^x -1\right)\;.
%old \]
%old This is mapped by linear transform $a + b x$ to the range
%old $(0,1)$ of the standard random generator.
%old From
%old \[ x=0 \quad : \qquad F_2^{-1}(a + b x) = \kappa^2 \, \left(e^a -1\right) = 0 \qquad  \Rightarrow \qquad a=0 \]
%old and 
%old \[ x=1 \quad : \qquad F_2^{-1}(a + b x) = \kappa^2 \, \left(e^b -1\right) = 1 \qquad  \Rightarrow \qquad 
%old b=\log\left( 1+1/\kappa^2\right) = \log\left( 1+W^2\right) \; .\]
%old The direct transformation $F_2^{-1}(a + b \, R)$, applied to the flat random distribution $R$ to generate $\rho^2$
%old is therefore
%old \begin{equation}
%old \rho^2 = \kappa^2 \left(e^{R \log\left( 1+1/\kappa^2\right)} -1\right)
%old = \kappa^2 \, \left[ \left( 1+1/\kappa^2\right)^R -1 \right]
%old = \frac{\left( 1+W^2\right)^R -1 }{W^2} \;.
%old \end{equation}
%old Generated distributions in $\rho^2$ are shown in Fig.\,\ref{plot:rho2.eps}.
%old The distribution is, with growing photon energy, increasingly peaked at 0.
%old Note that $\rho$ does not depend on $t$. One could do step 3 before 2 or after step 4.
The distribution in $\rho$ has the form
\begin{equation}\label{rho1}
f_3(\rho)\,d\rho=\frac{\rho^3\,d\rho}{\rho^4+C_2}\,,\quad
0\le\rho\le \rho_{\rm max}\,,\,
\end{equation}
where
\begin{equation}\label{rhomax}
\rho_{\rm max}^2=\frac{1.9}{A^{0.27}}\,\left(\frac{1}{t}-1\right), \,
\end{equation}
and
\begin{equation}\label{C2}
C_2=\frac4{\sqrt{x_+x_-}}\left[\left(\frac{m_\mu}{2E_\gamma x_+x_-\,t}
\right)^2+\left(\frac{m_e}{183 \, Z^{-1/3} \, m_\mu}\right)^2
\right]^2\,.
\end{equation}
The $\rho$ distribution is obtained by a direct transformation applied to 
uniform random numbers $R_i$ according to
\begin{equation}\label{rho}
\rho=\left[C_2(\exp(\beta\,R_i)-1)\right]^{1/4}\,,
\end{equation}
where
\begin{equation}\label{beta}
\beta=\log\left(\frac{C_2+\rho_{\rm max}^4}{C_2}\right)\,.
\end{equation}
Generated distributions of $\rho$ are shown in Fig.\,\ref{plot:rho.eps}

\begin{figure}[htpb]
% \begin{minipage}[t]{7.5cm}
% \includegraphics[scale=.6]{electromagnetic/standard/MuPgen/rho.eps}
\center\includegraphics[scale=.8]{electromagnetic/standard/MuPgen/rho.eps}
\caption{Histograms of generated $\rho$ distributions for beryllium at 
 two different photon energies.  The total number of entries at each energy
 is $10^6$.}
 \label{plot:rho.eps}
% \end{minipage}
\end{figure}
% \hfill
\begin{figure}[htpb]
% \begin{minipage}[t]{7.5cm}
% \includegraphics[scale=.6]{electromagnetic/standard/MuPgen/thetaPlus.eps}
\center\includegraphics[scale=.8]{electromagnetic/standard/MuPgen/thetaPlus.eps}
\caption{Histograms of generated $\theta_+$ distributions at different photon energies.}
\label{plot:thetaPlus}
% \end{minipage}
\end{figure}

\noindent
{\bf 5)} Calculate $\theta_+,\theta_-$ and $\varphi$ from $t, \rho,\psi$ with
\begin{equation}
\gamma_\pm = \frac{E_\mu^\pm}{m_\mu} \qquad {\rm and} \qquad
u=\sqrt{\frac1t-1}\,.
\label{eq:gammau}
\end{equation}
according to
\begin{equation}\label{s2}
\theta_+=
\frac{1}{\gamma_+}\,\left(u +\frac{\rho}{2}\,\cos\psi\right)\,,\quad \theta_-=
\frac{1}{\gamma_-}\,\left(u -\frac{\rho}{2}\,\cos\psi\right)\, \quad {\rm and}
\quad \varphi=\frac{\rho}{u} \, \sin\psi\, . \,
\end{equation}
The muon vectors can now be constructed from Eq.\,(\ref{eq:pvec}), where 
$\varphi_0$ is chosen randomly between 0 and $2\pi$.
Fig.\,\ref{plot:thetaPlus} shows distributions of $\theta_+$ at different 
photon energies (in beryllium).  The spectra peak around $1/\gamma$ as 
expected.

The most probable values are $\theta_+\sim m_\mu/E_\mu^+ = 1 / \gamma_+$. In the small angle
approximation used here, the values of $\theta_+$ and $\theta_-$
can in principle be any positive value from 0 to $\infty$.
In the simulation, this may lead (with a very small probability, of the 
order of $m_\mu/E_\gamma$) to unphysical events in which $\theta_+$ or
$\theta_-$ is greater than $\pi$.  To avoid this, a limiting angle 
$\theta_{\rm cut}=\pi$ is introduced, and the angular sampling repeated, 
whenever $\max(\theta_+,\,\theta_-)>\theta_{\rm cut}$ \, .

\begin{figure}[htpb]
\center\includegraphics[scale=.65]{electromagnetic/standard/MuPgen/Fig1.eps}
\caption{Angular distribution of positive (or negative) muons.
%old $E_\gamma=10\:{\rm GeV}$, $x_+=0.3$; iron.
The solid curve represents
the results of the exact calculations. The histogram is the simulated
distribution. The angular distribution for pairs created in the field
of the Coulomb centre (point-like target) is shown by the dashed curve
for comparison.}
\label{plot:Fig1}
\end{figure}

\begin{figure}[htpb]
% \begin{minipage}[t]{7.5cm}
% { \center\includegraphics[scale=.65]{electromagnetic/standard/MuPgen/Fig2.eps}
\center\includegraphics[scale=0.8]{electromagnetic/standard/MuPgen/Fig2.eps}
 \caption{Angular distribution in logarithmic scale.  The curve corresponds 
  to the exact calculations and the histogram is the simulated 
  distribution.}
 \label{plot:Fig2}
% }
% \end{minipage}
\end{figure}
% \hfill
\begin{figure}[htpb]
% \begin{minipage}[t]{7.5cm}  
% \includegraphics[scale=.65]{electromagnetic/standard/MuPgen/Fig3.eps}
\center\includegraphics[scale=0.8]{electromagnetic/standard/MuPgen/Fig3.eps}
\caption{Distribution of the difference of transverse momenta of positive 
 and negative muons (with logarithmic x-scale).}
  \label{plot:Fig3}
% \end{minipage}
\end{figure}

Figs.\,\ref{plot:Fig1},\ref{plot:Fig2} and \ref{plot:Fig3} show
distributions of the simulated angular characteristics of muon pairs in 
comparison with results of exact calculations.  The latter were obtained 
by means of numerical integration of the squared matrix elements with 
respective nuclear and atomic form factors.  All these calculations were 
made for iron, with $E_\gamma=10\,{\rm GeV}$ and $x_+=0.3$.  As seen from 
Fig.\,\ref{plot:Fig1}, wide angle pairs (at low values of the argument in 
the figure) are suppressed in comparison with the Coulomb center 
approximation.  This is due to the influence of the finite nuclear size
which is comparable to the inverse mass of the muon.  Typical angles of 
particle emission are of the order of
$1/\gamma_\pm=m_\mu/E_\mu^\pm$ (Fig.\,\ref{plot:Fig2}). Fig.\,\ref{plot:Fig3}
illustrates the influence of the momentum transferred to the target on the 
angular characteristics of the produced pair.  In the frame of the often 
used model which neglects target recoil, the pair particles would be 
symmetric in transverse momenta, and coplanar with the initial photon.


\subsection{Status of this document}
28.05.02 created by H.Burkhardt. \\
01.12.02 re-worded by D.H. Wright \\

\begin{latexonly}

\begin{thebibliography}{99}

\bibitem{MuPgen}
H.~Burkhardt, S.~Kelner, and R.~Kokoulin, ``Monte Carlo Generator for Muon 
Pair Production''. CERN-SL-2002-016 (AP) and CLIC Note 511, May 2002.

\bibitem{Kelner:1995hu}
S.~R. Kelner, R.~P. Kokoulin, and A.~A. Petrukhin, ``About cross section 
for high energy muon bremsstrahlung,''. Moscow Phys. Eng. Inst. 024-95, 1995.
  31pp.

\end{thebibliography}

\end{latexonly}

\begin{htmlonly}

\subsection{Bibliography}

\begin{enumerate}
\item H.~Burkhardt, S.~Kelner, and R.~Kokoulin, 
``Monte Carlo Generator for Muon 
Pair Production''. CERN-SL-2002-016 (AP) and CLIC Note 511, May 2002.

\item S.~R. Kelner, R.~P. Kokoulin, and A.~A. Petrukhin, ``About cross section 
for high energy muon bremsstrahlung,''. Moscow Phys. Eng. Inst. 024-95, 1995.
  31pp.
\end{enumerate}

\end{htmlonly}