\section{Monte Carlo Simulations}

\subsection{Model for Pion Electroproduction in Deuterium}

Because of the role of Fermi momentum in $^2$H, different choices can be made
in the construction of a quasi-free model for pion electroproduction from
deuterium.

\subsubsection{Initial Quasi-Free Model}

The first two analyses of our $^2$H data made use of the model presented in
detail in Ref.~\cite{Gaskellthesis}.  In this model, pion electroproduction is
assumed to follow a quasifree process, where the pion is produced from an
interacting proton (or neutron in the case of $\pi^{-}$ production) with some
initial momentum coming from the Fermi motion in the nucleus.  The cross
section is six fold, so six quantities must be generated, the additional
variable (compared to $^1$H) being the missing mass.  The variables of choice
are the spectrometer in-plane and out-of-plane angles for both the electron and
pion arms (SOS and HMS spectrometer angles), the final electron energy, and the
final pion momentum.  The deuterium wave function is based on a Bonn potential
calculation ~\cite{Koltenukthesis}, while the nucleon momentum is generated
according to the Fermi momentum distribution with the direction being generated
uniformly in cos$\theta_{CM}$ and $\phi_{CM}$.

Due to the fact that the struck nucleon is bound in the nucleus and its energy
is not specified by its momentum, an off-shell prescription was chosen to fully
determine the kinematics.  The ``spectator'' nucleon was taken to be
``on-shell'' in the initial state, while the struck nucleon was taken to be
``off-shell'' with the requirement that the total momentum of the nucleus is
zero, and the total energy is the mass of deuterium, $M_{D}$.
\begin{equation}
E_{struck} = M_{D} - \sqrt{m_{spectator}^2 +p_{Fermi}^2},
\end{equation}
where $m_{spectator}$ is the on-shell spectator nucleon mass and $p_{Fermi}$ is
the momentum of the struck and spectator nucleons.

Another ingredient used in this model is the Jacobian that accounts for the
fact that the events are generated in spectrometer coordinates instead of
physics coordinates.  Thus, one needs to calculate the transformation that
takes the reference system of coordinates for this reaction from the center of
mass to the lab system.  The Jacobian of the cross section transformation is
given in Appendix B of Ref.~\cite{Gaskellthesis}.  The resulting six-fold cross
section can be written as the virtual photon cross section times the virtual
photon flux $\Gamma$,
\begin{equation}
\frac{d^6}{d \sigma_{e} dE_{e} d \Omega_{\pi} dP_{\pi}} = \Gamma
\frac{d\sigma}{d\Omega_{\pi} dP_{\pi}}.
\label{eqn:d6sigma}
\end{equation}

Because the struck nucleon has Fermi momentum, $\phi_{CM}\neq\phi_{\pi q}$ in
this model.  $\phi_{CM}$ cannot be reconstructed for the experimental data,
since it is not possible to uniquely determine the Fermi momentum on an
event-by-event basis.  But this angle is used in the SIMC physics generator,
where the virtual photon cross section is calculated via
\begin{equation}
\frac{d\sigma}{d\Omega_{\pi} dP_{\pi}} = 
\frac{d\sigma_{T}}{d\Omega_{\pi}dP_{\pi}} 
+ \epsilon \frac{d\sigma_{L}}{d\Omega_{\pi} dP_{\pi}} 
+ \sqrt{2\epsilon(1+\epsilon)}\frac{d\sigma_{LT}}
  {d\Omega_{\pi}dP_{\pi}}cos\phi_{CM}
+ \epsilon \frac{d\sigma_{TT}}{d\Omega_{\pi} dP_{\pi}}cos2\phi_{CM},
\label{eqn:dsig_full}
\end{equation} 
where $\sigma_{L}$, $\sigma_{T}$, $\sigma_{LT}$, and $\sigma_{TT}$ are the
usual separated response functions.  Similarly, $\theta_{CM}$ is reconstructed
in the SIMC physics generator using the Fermi momentum.  However, for the
experimental data a different $\theta_{CM}'$ is reconstructed, by boosting to
the photon plus fixed nucleon target system.

The SIMC data are then reconstructed in a manner similar to that of the
experimental data.  The simulated yields are calculated and binned versus
$\phi_{\pi q}$, $-t$ and compared to the experimental yields for the same bins.
In order to alleviate problems such as the averaging over the phase space of
each $-t$ bin and the need for L/T separation to have the same values for $W$,
Q$^2$, and $-t$ at high and low $\epsilon$, a model was created such that the
ratio of experimental $Y_{exp}$ to simulated $Y_{sim}$ yields does not depend
(or hardly depends) on $W$, $Q^2$, $-t$, $\theta$ and $\phi$.
\begin{equation}
\left (\frac{d\sigma}{dt}\right )_{exp} = \frac{Y_{exp}}{Y_{sim}} 
\left (\frac{d\sigma}{dt}\right )_{MC}.
\label{eqn:yexp_ysim}
\end{equation} 
The iterative fitting procedure used to create the model for the cross section
is described in Sec.~\ref{sec:iterate}

\subsubsection{Final Quasi-Free Model}

The fact that in the first quasi-free model the $\theta_{CM}$, $\phi$ angles
used in the SIMC physics generator differ from those used in the experimental
data reconstruction leads to a philosophical quandry.  $\phi_{CM}$ and
$\theta_{CM}$ as used in the SIMC physics generator rely on knowledge of the
struck nucleon momentum (if it is moving).  Since this is not accessible
experimentally, how does one iterate the cross section since it is calculated
in terms of these variables?

Given that experimentally we do not know the movement of the struck nucleon,
and in SIMC we (have to) assume something, we ultimately decided to use in the
SIMC physics generator the same $\theta_{CM}$ calculated by boosting to the
photon plus nucleon at rest system, as is used in the experimental data
reconstruction.  In this model, $\phi_{CM}$ equals the experimentally
reconstructed $\phi_{\pi q}$, as was also the case for the $^1$H analysis.
This choice is clear, also for people who later want to look into or use our
cross sections.  This makes the extracted $\sigma_L$, $\sigma_T$,
$\sigma_{LT}$, and $\sigma_{TT}$ effective ones, not directly to be compared to
the ones on $^1$H, but it is better that the influence of off-shell (and
possible other mechanisms in $^2$H) are studied separately, using cross
sections that were determined in a well defined way, than that we incorporate
off shell effects already in some way in the extration of our cross sections
(although the difference in practice may not be large).

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[h]
\begin{center}
\vskip -0.05cm
\includegraphics[angle=0,width=0.6\textwidth]
{eps/dpl_160_63_+0000_var_comparo.eps}
%\vskip -.25cm
\caption{\label{fig:var_comparo} \small Comparison of kinematic variables used
  in the initial SIMC model ($y$-axis) versus the final SIMC model ($x$-axis)
  for the $\pi^+$ setting $Q^2=1.6$ GeV$^2$, $\epsilon=0.63$, 
  $\Theta_{\pi q}=0$.}
\end{center}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In addition to the revised calculation of $\theta_{CM}$ and the use of
$\phi_{\pi q}$ instead of $\phi_{CM}$ in Eqn. \ref{eqn:dsig_full}, the use of a
photon-stationary nucleon C.M. system instead of a photon-moving nucleon system
affects the calculation of $s=W^2$ and $-t$, as well as the Jacobian in
Eqn. \ref{eqn:d6sigma}, in the SIMC physics model.  A comparison of
$\theta_{CM}$, $\phi$ and $W$ used in the two models is shown in
Fig. \ref{fig:var_comparo}.

In parallel with these changes to the physics model, the SIMC event
reconstruction was carefully scrutinized for inconsistencies from the
experimental data reconstruction.  It was discovered that the experimental
reconstruction used $t=(P_n-P_p)^2$, while the default formula used in SIMC was
$t=(P_{\gamma}-P_{\pi})^2$.  In the limit of perfect resolution, no radiative
effects, or smearing, both formulas for $t$ should yield the same result.  But
they don't, because of these effects, and so $t=(P_n-P_p)^2$ was replaced in
all parts of SIMC, including the physics generator.  The differences between
the two definitions of $t$ can be large, as indicated in
Fig. \ref{fig:var_comparo}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\clearpage
\subsection{Simulation of Collimator Pion Punch-Through Events and Missing Mass
Cut Determination}

The thick collimators of the HMS and SOS are very effective in stopping
electrons, but a non-negligible fraction of pions undergo multiple
scattering and ionization energy loss and consequentially end up contributing
to the experimental yield~\cite{Gaskellthesis}.

These pion (hadron) punch-through events have been observed in earlier
experiments, and corrections are needed for a precise yield extraction.  Since
the pions in F$_{\pi}$-1 and F$_{\pi}$-2 are detected in the HMS, the
implementation of the collimator punch-through events was done for only this
arm.  To account for these events, the probability that a pion interacts
hadronically with the collimator (allowing the pion to undergo multiple
scattering and ionization energy loss) was tested ~\cite{Gaskellthesis} and
implemented in the HMS single arm Monte Carlo package as a part of the SIMC
Hall C simulation software.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[h]
\begin{center}
\vskip -0.05cm
\includegraphics[angle=0,width=0.6\textwidth]{eps/MM_pp.eps}
\vskip -.25cm
\caption[Collimator Pion Punch-Through Implementation in SIMC Monte Carlo
Model.]{\label{fig:CPPT} \small Collimator Pion Punch-Through Implementation in
the SIMC Monte Carlo Model for the F$_{\pi}$-1 $\pi^{+}$ data set. The upper
panel shows the normalized experimental yields (full circles), simulation
yields with (red line) and without (black line) collimator punch-through
implementation.  The middle panel shows the ratio of data to simulation with
(full red circles) and without (blue triangles) collimator punch-through
implementation.  The lower panel shows the ratio of the simulated data without
collimator punch-through to the simulated data with this implementation.}
\end{center}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

As can be seen in the lower panel of Fig.~\ref{fig:CPPT}, collimator
punch-through events show up in the long tail of the missing mass distribution.
In extracting the deuterium cross sections, it is desirable to keep as much of
the missing mass tail as possible (up to the two-pion threshold of 1.1 GeV), to
maximize the acceptance of the ``quasifree'' distribution, and to minimize the
systematic uncertainty associated with the missing mass (MM) cut.

After implementing the pion punch-through events in SIMC, we determined the
MM cut value needed to reliably select ``quasielastic'' events
with an electron detected in SOS and a pion detected in HMS. The MM is expected
to approach the ``spectator'' nucleon mass (proton mass $m_{p}$ for $\pi^{-}$
and neutron mass $m_{n}$ for $\pi^{+}$ electroproduction processes). The final
MM cut used in the analysis was $0.875 < MM < 1.03$ GeV.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[ht]
\begin{center}
\vskip -0.05cm
\includegraphics[angle=90,width=8.cm,totalheight=6.cm]{eps/MMcut_Fpi1_unc.ps}
\includegraphics[angle=90,width=8.cm,totalheight=6.cm]{eps/mmcut_unc.ps}
\vskip -0.05cm
\caption[Systematic uncertainties associated with the MM cut.]
{\label{fig:MM_systematics} \small RMS variations in the cross section
associated with the MM cut for F$_{\pi}$-1 (left) and F$_{\pi}$-2 (right).  The
green solid lines denote the average of both high and low $\epsilon$ for
F$_{\pi}$-1 and the average of $\sigma_L$ and $\sigma_T$ for F$_{\pi}$-2, while
the blue solid lines denote the average of each separately.  The dashed lines
indicate the standard deviations of the respective averages.}
\end{center}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The systematic uncertainties associated with this cut were computed in a
slightly different fashion for F$_{\pi}$-1 and F$_{\pi}$-2, but the results
obtained were similar.  For F$_{\pi}$-1, the missing mass cut was varied in
increments of 0.025 GeV from 1.00 to 1.10 GeV for each $-t$ bin at each $Q^2$
and the unseparated cross sections extracted separately at high and low
$\epsilon$.  There was little variation with $-t$ bin, so the rms variations
for each $Q^2$ and charge combination were averaged together.  These variations
are shown in the left panel of Fig. \ref{fig:MM_systematics}, along with the
averages and standard deviations for each $\epsilon$ separately in blue and
together in green.  It is clear that there is no systematic trend with either
$Q^2$ or $\epsilon$.  One might argue that on average the $\pi^-$ cross
sections experience slightly larger variation than the $\pi^+$, but the trend is
weak.  The overall average variation in the cross section due to the MM cut is
$1.34\pm 0.76\%$.  Based on these results, we assigned a scale uncertainty of
1.3\% to $\sigma_{L, T}$, with an additional point-to-point uncertainty of
0.8\%.

The systematic uncertainty due to the missing mass cut is evaluated in a
similar fashion for F$_{\pi}$-2, except that the cut was varied in
increments/decrements of 0.010 GeV and the four separated cross sections
$\sigma_{L, T, LT, TT}$ were extracted in each case.  The averages
and standard deviations for $\sigma_{L, T}$ are shown in the right panel
Figure~\ref{fig:MM_systematics}, separately in blue and their combined average
in green.  The overall average uncertainty is $\approx 0.7\%$ in all cases. The
standard deviation for $\sigma_{L}$ is 0.6\% while for $\sigma_{T}$ it is s
0.9\%. The combined standard deviation is 0.8\%.  Based on these results, we
assigned a scale uncertainty of 0.7\% to $\sigma_{L, T}$, with an additional
point-to-point uncertainty of 0.6-1.5\%.

\subsection{Selection of $t$ Bins for F$_{\pi}$-1 $\pi^-$ Data}

Because of problems encountered during the running of the F$_{\pi}$-1
experiment, not all of the planned HMS settings were acquired for the low
$\epsilon$ $\pi^-$ data at $Q^2$=0.6 and 1.0 GeV$^2$.  This led to inadequate
$\phi$ coverage at larger $-t$, degrading the reliability of the L/T
separations for those kinematics.  To determine the appropriate range in $-t$
for which reliable separated cross sections could be determined in these cases,
the following study was performed.  Since radiative processes can allow
otherwise excluded kinematic regions to otherwise be populated, the $\phi$
coverage of each $-t$ bin was simulated with electron radiation turned on and
off.  Fig.~\ref{fig:rad_on_off} shows a representative example at $Q^2$=0.6
GeV$^2$.  Based on the criterion that it is necessary to have
$\phi_{norad}>270^o$ for reliable L/T separations, it was decided that the
region $-t>0.07$ GeV$^2$ for the $Q^2$=0.6 GeV$^2$ data (corresponding to the
last $-t$ bin in the left panel of Fig.~\ref{fig:rad_on_off}), and the
region $-t>0.13$ for the $Q^2$=1.0 GeV$^2$ data (corresponding to the last two
$-t$ bins in the right panel of Fig.~\ref{fig:rad_on_off}) had inadequate
$\phi$ coverage.  In order to provide some constraint to the Monte Carlo
model fitting procedure (discussed in the next section), these regions were
retained in the analysis.  But reliable separated cross sections for these
regions cannot be determined, and so will not be quoted.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[ht]
\begin{center}
\vskip -0.05cm
\hskip -0.05cm
\includegraphics[angle=0,width=7.75cm,totalheight=11.5cm]{eps/phi_coverage_mn_060_37.ps}
\hskip 0.05cm
\includegraphics[angle=0,width=7.75cm,totalheight=11.5cm]{eps/phi_coverage_mn_100_33.ps}
\vskip -0.05cm
\caption[Simulated $\phi$ coverage for the $Q^2$=0.6 GeV$^2$, low $\epsilon$
settings.]{\label{fig:rad_on_off} \small Simulated $\phi$ coverage for the
$Q^2$=0.6 GeV$^2$, low $\epsilon$ settings.  Due to poor $\phi$ coverage
($<270^O$), separated cross sections cannot be quoted for the last $-t$ bin for
$Q^2$=0.6 GeV$^2$ (left panel) and the last two $-t$ bins for $Q^2$=1.0 GeV$^2$
(right panel).}
\end{center}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Iteration of the Monte Carlo cross section model
\label{sec:iterate}}

The starting pion electroproduction cross section model used in the simulation
was taken from previous experimental data.  The individual response functions
$\sigma_{L}$, $\sigma_{T}$, $\sigma_{LT}$, and $\sigma_{TT}$ were then
optimized to follow the $-t$ and $Q^2$ dependences of the $^2$H experimental
data through the procedure described below.  The dependence on $W$ was assumed
to follow the phase space factor $1/(W^2-M_{p}^2)^2$.

First, the ratios of experimental to simulated data were computed for each
($-t, \phi$) bin, separately for each HMS central angle setting.  Then these
ratios were combined in a weighted average over the experimental $\phi$
acceptance.  Using Eqn.~\ref{eqn:yexp_ysim}, the unseparated total experimental
cross sections were computed for each $-t$ bin and the $\phi$ dependences were
fitted, extracting simultaneously all four separated response functions
$\sigma_{L, T, LT, TT}$ versus $t$.  Fitting all four response functions
simultaneously is advantageous, as it allows one to calculate all four terms of
the cross section, including the associated uncertainties.  As discussed in
Ref.~\cite{Tanjathesis}, this method of cross section extraction relies on the
assumption that the variations between simulation and experimental data cross
section are small, and no modification of the model functional form is
necessary.

In order to ensure a deliable description of the data by the SIMC physics
model, the model cross section was optimized using a standalone FORTRAN program
which incorporates a MINUIT routine for fitting corrections to the individual
response functions in a $\chi^2$ minimization fashion ~\cite{MINUIT}.  For each
iteration, the experimental cross section was calculated using
Eqn. \ref{eqn:yexp_ysim}, and the four response functions were fitted
simultaneously from the $\phi$-dependences at both $\epsilon$ settings.  The
separated response functions $\sigma_L$, $\sigma_T$, $\sigma_{LT}$ and
$\sigma_{TT}$ are in fact the fit parameters.  Finally, the separated cross
sections were fitted to the model cross section in a multiple parameter fit
using the MINUIT routine described above, which was then used as an input to
the next iterative calculation of the experimental cross section.  In this way,
the $Q^2$ and $-t$ and dependences of each of the cross section terms were
iteratively converged to the data.

The entire iteration procedure can be summarized in five basic steps:
\begin{itemize}
\item{use parameter values $p_{i=1,...,12}$ to generate simulated yields,
  cut and binned in the same manner as the experimental data [simc,
  sim\_yields.kumac],}
\item{determine the ratios of experimental to simulated yields,
${Y_{exp}}\{Y_{sim}$ [average\_ratios.f, plot\_r.kumac]}
\item{calculate the unseparated experimental cross sections for each 
$(\phi, -t)$ bin using Eqn.~\ref{eqn:yexp_ysim} [calc\_xsect.f,
    plot\_unsep.kumac]}, 
\item{fit the $\phi$-dependences at both epsilon to determine the separated
cross sections [xfit\_in\_phi.f, plot\_sep.kumac]},
\item{fit the $t$-dependences of the separated cross sections to determine
  the set of parameter values to be used in the next iteration [xfit\_in\_t.f]}
\end{itemize}
Through this process, some parameters were held fixed while others were allowed
to vary.  After several iterations, the parameter values are reasonably
converged and the kinematic dependences of the various parameters and the
${Y_{exp}}/{Y_{sim}}$ ratios were investigated and the parameterizations
and starting parameter values modified where necessary to yield a better fit.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[h]
\begin{center}
\vskip -0.2cm
\hskip -0.5cm
\includegraphics[angle=0,width=8.7cm,totalheight=15.cm]
{eps/exp_kinematics_plw_corrected_p29.ps}
\hskip -0.7cm
\includegraphics[angle=0,width=8.7cm,totalheight=15.cm]
{eps/exp_kinematics_plw_corrected_p30.ps}
\vskip -1.5cm
\caption{\label{fig:Kinematics_MM}
\small Experimental (black points) to simulated (red) data comparison for
various physics and HMS/SOS optical transport quantities for the $\pi^+$
setting $Q^2=1.0$ GeV$^2$, $\epsilon=0.65$, $\Theta_{\pi q}=-4^o$.}
\end{center}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

In addition to the ${Y_{exp}}/{Y_{sim}}$ ratios versus $-t, \phi, \theta_{CM}$,
the normalized experimental yields and simulated yields for a variety of
physics and HMS and SOS optical transport matrix variables were compared during
the iteration process, in order to asses the quality of match between the data
and the simulation.  Some examples are shown in Fig.~\ref{fig:Kinematics_MM}.
The iteration procedure was repeated until satisfactory agreement between the
experimental and simulated distributions was obtained, the values of
$\sigma_{L,T,LT,TT}$ (and the associated fit parameters) were stable in
subsequent iterations, and the most important parameters not change too
much versus $Q^2$.  

The final parameterizations used for the four separated response functions in
the $\pi^+$ analysis are:
\begin{eqnarray}
\frac{d\sigma_{L}}{dt} &=& g(W) \biggl(p_{1}+p_{2} ln(Q^2)\biggr)
e^{(p_{3}+p_{4}ln(Q^2))(-t)},\nonumber \\
\frac{d\sigma_{T}}{dt} &=& g(W) \biggl( p_{5}+p_{6}ln(Q^2)+
\biggl(p_{7}+p_{8}ln(Q^2)\biggr)\biggl(
\frac{|t|-|t_{ave}|}{|t_{ave}|}\biggr)\biggr),\nonumber \\
\frac{d\sigma_{LT}}{dt} &=& g(W) p_{9} e^{p_{10}(-t)} sin\theta_{CM},\nonumber \\
\frac{d\sigma_{TT}}{dt} &=& g(W) f(t) 
\frac{p_{11}}{Q^2 }e^{-Q^2} sin^2\theta_{CM}, \nonumber \\
\label{eqn:mc_model_pl}
\end{eqnarray}
where $g(W)=1/(W^2-M_{p}^2)^2$ is the phase space factor (used in previous
analyses of electroproduction data described in Ref.~\cite{Brauel77,Bebek78}),
$f(t)=-t/(-t-m_{\pi}^2)^2$ is the pion pole factor, $|t_{ave}|$ is the average
$-t$ value for a given kinematic setting, given by $|t_{ave}|=(0.105+0.04
ln(Q^2)Q^2$, and $p_{i=1,...,12}$ are the fit parameters.

For the $\pi^-$ analysis, a slightly different parameterization yielded a
better fit:
\begin{eqnarray}
\frac{d\sigma_{L}}{dt} &=& g(W) \biggl(p_{1}+p_{2} ln(Q^2)\biggr)
e^{(p_{3}+p_{4}ln(Q^2))(-t)},\nonumber \\
\frac{d\sigma_{T}}{dt} &=& g(W) \biggl( p_{5}+\frac{p_{6}}{Q^4+0.1}+
\biggl(p_{7}+p_{8}ln(Q^2)\biggr)\biggl(
\frac{|t|-|t_{ave}|}{|t_{ave}|}\biggr)\biggr),\nonumber \\
\frac{d\sigma_{LT}}{dt} &=& g(W) p_{9} e^{p_{10}(-t)} sin\theta_{CM},\nonumber \\
\frac{d\sigma_{TT}}{dt} &=& g(W) f(t) \biggl(
\frac{p_{11}}{Q^2 }+\frac{p_{12}}{Q^4+0.2}\biggr) sin^2\theta_{CM}, \nonumber \\
\label{eqn:mc_model_mn}
\end{eqnarray}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[h]
\begin{center}
\vskip -0.2cm
\includegraphics[angle=90,width=8.7cm]{eps/param_q2_it27.ps}
\caption{\label{fig:pl_it27}
\small Final parameter values for the $F_{\pi}-1$ $\pi^+$ data, plotted versus
$Q^2$.  The parameters not shown were fixed to the following values:
$p_2=-785$, $p_4=21.2$, $p_6=-34$, $p_8=-60$, $p_9=p_{12}=0$, $p_{10}=-20$.
Note that to better investigate the low $Q^2$-dependence, the $Q^2=0.60$ setting
was subdivided into $Q^2=0.55$, 0.65 GeV$^2$ settings and iterated separately
from the others.
}
\end{center}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[h]
\begin{center}
\vskip -0.2cm
\includegraphics[angle=90,width=8.7cm]{eps/param_q2_it29.ps}
\caption{\label{fig:mn_it29}
\small Final parameter values for the $F_{\pi}-1$ $\pi^-$ data, plotted versus
$Q^2$.  The parameters not shown were fixed to the following values:
$p_2=-205$, $p_4=10$, $p_5=0$, $p_8=45$, $p_9=-200$, $p_{10}=-20$, $p_{11}=0$.
$p_{12}$ varies strongly with $Q^2$, but the contribution of TT is small.
}
\end{center}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{figure}[h]
\begin{center}
\vskip -1.0cm
\includegraphics[angle=0,width=10cm,totalheight=15.cm]{eps/plot_unsep_mn_p1.ps}
\vskip -3.0cm
\caption{\label{fig:Phifit} \small Unseparated experimental $\pi^-$ cross
sections $\frac{d^2\sigma}{dt d\phi}$ as a function of azimuthal angle $\phi$
at $Q^2=0.60$ GeV$^2$ high (blue) and low (black) $\epsilon$.  The curves shown
represent the fit of the measured values of the cross section to
Eqn.~\ref{eqn:dsig_full}.  In this fit, all four parameters $\sigma_{L, T, LT,
TT}$ are extracted simultaneously, separately for each $-t$ bin.
}
\end{center}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

The final sets of parameter values for the $F_{\pi}-1$ data are shown in
Figs. \ref{fig:pl_it27}, \ref{fig:mn_it29}.  A representative example of the
measured experimental cross sections as a function of $\phi$ for a particular
$-t$ bin and for low and high $\epsilon$ is shown in Fig.~\ref{fig:Phifit}.

\subsection{SIMC Model Dependence}

The SIMC model dependence was estimated by modifying the fit parameters and
investigating the resulting differences (uniterated) on the separated cross
sections.  For the $\pi^-$ settings, LT and TT were separately set to zero
(keeping other parameters at their final values).  The variation of the
separated cross sections in four additional late-stage iterations, including
different parameterizations with non-pathologic behavior, were also
investigated.  For each $Q^2, t$ bin, the error-weighted mean, the standard
deviations of each iteration and of the mean with respect to the official
value, and their error-weighted fractional differences were computed.  The
model-dependences of the L,T cross sections were similar to each other, and did
not exhibit any clear $Q^2$ or $t$ dependence.  From these variations, a
point-to-point uncertainty of 0.4\%, an $\epsilon$-uncorrelated, $t$-correlated
(which could affect the $t$-dependence) uncertainty of 1.4-3.9\%, and a scale
uncertainty of 0.3-2.0\% was computed.  These variations are comparable or
smaller than the statistical uncertainties in these cross sections (per
$t$-bin) of 5-10\%.

The procedure for the $\pi^+$ settings was similar, except that LT=0 in the
default model.  In this case, separated cross sections with $p_9=\pm 100$ were
computed.  In addition, the variation of the separated cross sections in seven
additional late-stage iterations was investigated.  From these variations, a
point-to-point uncertainty of 0.4\%, an $\epsilon$-uncorrelated, $t$-correlated
(which could affect the $t$-dependence) uncertainty of 0.7-3.3\%, and a scale
uncertainty of 0.3-1.7\% was computed.  These uncertainties are expected to be
smaller than the $\pi^-$ model uncertainties because of the substantially more
complete $\phi$-coverage, especially at $Q^2=0.6$, 1.0 GeV$^2$.

In contrast to the L,T model dependences, the variation in the LT, TT
interference response functions is strongly $t$-dependent.  The model
dependences are least for the smallest $-t$ bin at each $Q^2$, and increases
for the larger $-t$ bins.  These higher $-t$ bins have relatively poorer
statistics as well as incomplete $\phi$ coverage at low $\epsilon$ [as well as
for high $\epsilon$ at $\pi^-$ $Q^2=0.6$, 1.0 GeV$^2$].  As a result, they are
computed separately for each $t$-bin.
