\documentclass{article}

\usepackage{graphicx}

\title{A Program to Estimate Resolution for Charged Particles in GlueX \\
~\\
\large GlueX-doc-1015-v1
}
\author{Mark M.\ Ito \\
Thomas Jefferson National Accelerator Facility \\
12000 Jefferson Avenue \\
Newport News, VA 23606
}
\date{April 2, 2008}

% global parameters
\textwidth=6.5in
\oddsidemargin=0in % use built-in offset of 1 inch for left margin
\evensidemargin=0in % ditto for even pages
\textheight=9in
\topmargin=0in
\headheight=0in % no headers in this document
\headsep=0in

\begin{document}

\bibliographystyle{plain}

\maketitle

\begin{abstract}
This note describes a set of FORTRAN routines ({\tt REZEST}) that
calculates estimates of charged track resolution in transverse momentum
and direction for the GlueX detector geometry. Parameters of that
geometry can be varied to quickly obtain estimates for new
configurations. Since no Monte Carlo is used in the calculations,
results are returned immediately. Instructions on how to obtain and
use the software are provided and a description of the approximations
used is given.
\end{abstract}

\tableofcontents

\section{Introduction}

We often use ``back of the envelope'' estimates of detector
performance to inform detailed design decisions. These calculations
have the advantage of simplicity and provide quick feedback for new
design ideas. The FORTRAN routines described here provide a facility
for making such calculations for GlueX.

The fixed detector geometry concept used is the standard one for
GlueX: a central cylindrical drift chamber (CDC) in a solenoidal field
with a set of planar drift chambers in the forward direction
(FDC). The CDC is assumed to have both axial and stereo layers. The
parameters that can be varied include:

\begin{itemize}
\item $z$ position\footnote{The $z$-direction is along the incident photon beam.} of the target
\item inner and outer radius of the CDC
\item length of the CDC
\item front and back $z$-positions of the FDC
\item amount of material in the FDC and CDC tracking volume
\item amount of material in front of the FDC and CDC
\item magnetic field strength
\item number of axial and stereo position measurements in the CDC
\item stereo angle for stereo layers in the CDC
\item number of position measurements in the FDC
\end{itemize}

Many
approximations are made and generic formulae applied in the
calculations. The numbers that are produced are not a substitute for a
full-scale, hit-based Monte Carlo simulation of the detector. That
having been said, the relative variation of resolution when a particular
parameters are varied should give a good feeling for the effect of
parameter change.

The following assumptions are made:
\begin{itemize}
\item The magnetic field is uniform everywhere.
\item Particles travel in straight lines, independent of momentum.
\item Measurements of track parameters in the FDC are statistically
  independent from those made in the CDC.
\item All position measurements within a detector (FDC or CDC) are
  statistical independent of one another.
\item All positions measurements within a detector are made at
  locations uniformly spaced along the trajectory.
\item All positions measurements within a detector have the same
  resolution.
\end{itemize}

A note on notation: in the following, $\sigma_x$ and $\delta x$ will
be used interchangeable to represent the root mean square deviation of
the quantity $x$.

\section{Using the Software}\label{sec:using}

The principal subroutine is documented in the source code file {\tt
  rezest\_fdc\_cdc.F}:

\small
\begin{verbatim}
      SUBROUTINE REZEST_FDC_CDC(P, LAMBDA, M,
     X     DP_OVER_P, DPHI_TOT, DTHETA_TOT)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C This routine estimates the resolution in GlueX for charged particles
C in tranverse momentum, azimuthal angle, and polar angle.
C
C Input arguments, all REAL*4
C
C     P        Magnitude of total momentum (GeV/c)
C     LAMBDA   Dip angle, difference in polar angle in lab between track
C              and pi/2 (i. e., 90 degrees) (radians)
C     M        Mass of the particle (GeV/c^2)
C
C Output arguments, all REAL*4
C
C     DP_OVER_P  Relative resolution in transverse momentum
C                ("sigma_{p_t}/p_t")
C     DPHI_TOT   Resolution in azimuthal angle ("sigma_phi")
C     DTHETA_TOT Resolution in polar angle ("sigma_theta")
C
C The routine combines the measurements in the FDC and CDC where
C appropriate. Parameters describing the geometry and materials are
C defined in the routine REZEST_COMPONENTS which appears below.
C
\end{verbatim}
\normalsize

\subsection{Getting the code}

Two methods:

\begin{enumerate}
\item Get the tar ball from
\begin{center}
{\tt http://www.jlab.org/$\sim$marki/misc/rezest.tar}
\end{center}
\item Check it out from the subversion repository with the command
\begin{center}
{\tt svn checkout
 https://halldsvn.jlab.org/repos/trunk/home/marki/gluex/rezest}
\end{center}
\end{enumerate}

\subsection{Building the files}

There is a simple makefile in the directory:

\begin{verbatim}
> cd <rezest directory>
> make
gfortran -g   -c -o rezest.o rezest.F
gfortran -g   -c -o rezest_fdc_cdc.o rezest_fdc_cdc.F
ar rcv librezest.a rezest.o rezest_fdc_cdc.o
a - rezest.o
a - rezest_fdc_cdc.o
gfortran -g   -c -o rezest_point.o rezest_point.F
gfortran -o rezest_point rezest_point.o -L. -lrezest
gfortran -g   -c -o rezest_point_comp.o rezest_point_comp.F
gfortran -o rezest_point_comp rezest_point_comp.o -L. -lrezest
\end{verbatim}

This creates three files that you care about:

\begin{enumerate}
\item {\tt librezest.a}: the object library
\item {\tt rezest\_point}: a binary
\item {\tt rezest\_point\_comp}: a binary 
\end{enumerate}

\subsection{Using the files}

\subsubsection{\tt librezest.a}

Link this into your own program to retrieve resolution values. It
contains the routine {\tt rezest\_fdc\_cdc} described in
Section~\ref{sec:using}.

\subsubsection{\tt rezest\_point}

This stand-alone binary will accept three arguments on standard input
(space separated on a single line) and will write the resolution
values on standard output.

\begin{description}
\item[Input:] ~ 
\begin{enumerate}
\item total magnitude of momentum in GeV/$c$
\item dip angle in radians (polar angle referenced from $90^\circ$ in
  the lab; $\theta_{\rm dip} = \pi/2 - \theta$ there $\theta$ is the standard
  polar angle)
\item mass in GeV/$c^2$ 
\end{enumerate}
\item[Output:] ~
\begin{enumerate}
\item relative resolution in transverse momentum ($\delta p_t/p_t$)
\item resolution in azimuthal angle ($\delta\phi$)
\item resolution in polar angle ($\delta\theta$)
\end{enumerate}
\end{description}

For example:

\begin{verbatim}
> rezest_point
1.0 1.22 0.139
  2.8597742E-02  1.6185496E-02  8.1181811E-04
\end{verbatim}
says that a 1~GeV/$c$ particle traveling in a direction 1.22 radians
forward of the transverse direction with a mass of 139~MeV has an
estimated $\delta p_t/p_t = 2.9\%$, $\delta\phi = 16$~mR, and
$\delta\theta = 0.8$~mR.

\subsubsection{\tt rezest\_point\_comp}

Same as {\tt rezest\_point} except that the individual components of the
resolution are listed. Note that the first three numbers being output
are the same as for {\tt rezest\_point}.

Output:

\begin{enumerate}
\item relative resolution in transverse momentum
\item resolution in azimuthal angle
\item resolution in polar angle
\item resolution in curvature ($k = 1/R$) due to multiple
  scattering in the CDC in inverse meters
\item same for FDC
\item resolution in curvature due to position resolution in the CDC
\item same for FDC
\item total resolution in curvature in the CDC
\item same for the FDC
\item resolution in azimuthal angle ($\phi$) due to multiple scattering in
  the CDC in radians
\item same for FDC
\item resolution in azimuthal angle due to position resolution in the CDC
\item same for FDC
\item resolution in azimuthal angle due to curvature resolution in the CDC
\item same for FDC
\item total resolution in azimuthal angle in the CDC
\item same for the FDC
\item resolution in polar angle ($\theta$) due to multiple scattering in the
  CDC in radians
\item same for FDC
\item resolution in polar angle due to position resolution in the CDC
  in radians
\item same for FDC
\item total resolution in polar angle in the CDC
\item same for the FDC 
\end{enumerate}

For example:

\begin{verbatim}
> rezest_point_comp
1.0 1.22 0.139
  2.8597742E-02  1.6185496E-02  8.1181811E-04  4.2887099E-02
  2.2642065E-02  1.3919008E-02  0.1605156      4.5089267E-02
  0.1621047      2.5197803E-03  3.1950074E-04  5.0111138E-04
  2.4516270E-03  1.7020071E-02  4.7496673E-02  1.7212879E-02
  4.7560975E-02  2.5197803E-03  3.1950074E-04  7.7741774E-04
  7.9118297E-04  2.6369814E-03  8.5325917E-04
\end{verbatim}
(All values in the output actually appear on one line.)

\section{How the Estimates Are Done}

In this section the various formulae and concepts used to produce the
estimates are described.

\subsection{Transverse Momentum Resolution}

The formulae used to estimate transverse are taken from the Particle
Data Group's Review of Particle Physics\cite{Yao:2006px}. Using (for
the most part) their notation and descriptions, for a particle with
charge $q$ of momentum $p$ in a uniform magnetic field $B$ with a
pitch angle (dip angle) $\lambda$
\begin{equation}\label{eq:p-and-B}
p\cos\lambda = (0.3)qBR
\end{equation}
where $r$ is the radius of curvature in the projection of the
trajectory onto the bend plane, $p$ is in GeV/c, $B$ is in Tesla, and
$R$ is in meters. In the remainder of this note it will be assumed that
$q = 1$. The curvature $k = {1 \over R}$.  The variance of $k$ has two
contributions,
\begin{equation}
(\delta k)^2 = (\delta k_{\rm res})^2 + (\delta k_{\rm ms})^2 .
\end{equation}
\begin{equation}
\delta k_{\rm res} = {\epsilon \over L^{\prime 2}} \sqrt{720 \over N + 4}
\end{equation}
where $\epsilon$ is the position resolution in meters, $L^{\prime}$
is the projected length of the track onto the bending plane in meters
and $N$ is the number of measurements.
\begin{equation}
\delta k_{\rm ms} =
{(0.016~{\rm GeV}/c)z \over Lp\beta\cos^2\lambda} \sqrt{N_{\rm rl}}
\end{equation}
where $N_{\rm rl}$ is the number of radiation lengths in the
detector and $L$ is the total track length in the detector.

\subsection{Curvature vs.\ Momentum}

The estimates in the previous section are given in terms of the error
on curvature. These must be converted to momentum. Rewriting
Eq.~\ref{eq:p-and-B}

\begin{equation}
p_t = (0.3)BR = {(0.3)B \over k}
\end{equation}
So
\begin{equation}
{dp_t \over dk} = -{(0.3)B \over k^2} = -(0.3)BR^2
\end{equation}
and since
\begin{equation}
\delta p_t = \left\vert dp_t \over dk \right\vert \delta k ,
\end{equation}
we have
\begin{equation}
{\delta p_t \over p_t} = R\delta k = {\delta k \over k} .
\end{equation}

\subsection{Error on Slope and $y$-intercept of a Straight-Line Fit, Equally
Spaced Measurements}

To estimate the error due to position resolution on the direction of a
fitted track, we use the error on the slope of a straight line fitted
to the same number of measurements as the fundamental input. This
approach ignores the fact that the trajectory is in fact curved, but
is being used to estimate the error in the angle and not the angle
itself.

From Ref.~\cite{Bevington+92}
\begin{equation}
\chi^2 = \sum\left[{1 \over \sigma_i^2}(y_i - a - bx_i)^2\right]
\end{equation}
For equal errors ($\sigma_i = \sigma$ $\forall i$,),
$\chi^2$ is minimized for $a$ and $b$ when
\begin{equation}
a = {1 \over \Delta^\prime}\left(\sum{x_i^2}\sum{y_i}
- \sum{x_i}\sum{x_iy_i}\right)
\end{equation}
\begin{equation}
b = {1 \over \Delta^\prime}\left(n\sum{x_iy_i} - \sum{x_i}\sum{y_i}\right)
\end{equation}
where
\begin{equation}
\Delta^\prime = n\sum{x_i^2} - \left(\sum{x_i}\right)^2
\end{equation}
The variance of the parameters $a$ and $b$ are
\begin{equation}
\sigma_a^2 \approx {\sigma^2 \over \Delta^\prime}\sum{x_i^2}
\end{equation}
and
\begin{equation}
\sigma_b^2 \approx {n\sigma^2 \over \Delta^\prime} .
\end{equation}
For $n$ equally spaced measurements spanning the interval $[0,L]$,
\begin{equation}
x_i = {L(i-1) \over n-1}
\end{equation}
and using expressions from Appendix~\ref{app:algebra}, we get
\begin{equation}
\sigma_a^2 = {2\sigma^2(2n-1) \over n(n+1)}
\end{equation}
\begin{equation}
\sigma_b^2 = {12\sigma^2(n-1) \over L^2n(n+1)}
\end{equation}

\begin{figure}
\centering\includegraphics[height=2in]{slope.png}
\caption{Generic straight line through origin.}\label{fig:slope}
\end{figure}

We need to translate an error in slope to an error in angle.
In Fig.~\ref{fig:slope}, $\theta = \tan^{-1}b$ so
\begin{equation}
\delta\theta = \left| d\theta \over db \right| \delta b
= {\delta b\over \sec^2 \theta}.
\end{equation}

\subsection{Angular Error Due to Multiple Coulomb Scattering}

\begin{figure}
\centering\includegraphics[width=5.0in]{ms_pdg.png}
\caption{Quantities used to describe multiple Coulomb scattering. The
  particle is incident in the plane of the figure.}\label{fig:ms-pdg}
\end{figure}

Again from Ref.~\cite{Yao:2006px}, we define
\begin{equation} \theta_0 = \theta^{\rm rms}_{\rm plane} = {1
\over
 \sqrt 2}\theta^{\rm rms}_{\rm space} .
\end{equation}
The central angular distribution is approximately Gaussian with a
width given by
\begin{equation}
\theta_0 = {(13.6~{\rm MeV}) \over \beta c p} z \sqrt{x/X_0}
\left[
1 + 0.038 \ln(x/X_0)
\right] .
\end{equation}
As seen in Fig.~\ref{fig:ms-pdg}The angle $\Psi_{\rm plane}$ is the
angle between an unscattered trajectory and a line drawn from the
entrance point of the detector to the exit point. We use this as an
approximation to the contribution of multiple scattering to both the
azimuthal and polar angles.
\begin{equation}
\Psi^{\rm rms}_{\rm plane} = {1 \over \sqrt 3} \theta_0 .
\end{equation}

\subsection{Contribution to Azimuthal Angle Resolution \\
  from Curvature Resolution}

\begin{figure}
\centering\includegraphics[height=4in]{curve_error_in_phi.png}
\caption{The concept for estimating the effect of curvature resolution
  on azimuthal angle resolution.}
\end{figure}

Curvature $k = 1/R$ and direction in plane is measured at some
position rotated from the vertex by an angle $\alpha$ about the center
of curvature, not at the vertex. To infer the azimuthal angle $\phi$
at the vertex, track must be swum backward through angle
$\alpha$. Determination of $\alpha$ depends on $R$ and thus on $k$. An
error in $\alpha$ translates directly into an error in $\phi$.
\begin{equation}
\sin{\alpha \over 2} = {r_{\rm mid} \over 2R} = {r_{\rm mid}k \over 2}
\end{equation}
so
\begin{equation}
\delta\alpha = r_{\rm mid} \delta k \sec{\alpha \over 2}
\end{equation}
where $\alpha = 2\sin^{-1}(r_{\rm mid}/2R)$
We take as an approximation $r_{\rm mid} = (r_{\rm in} + r_{\rm out})/2$.

Note that a side effect of using a straight-line approximation for the
trajectory is that in some cases $r_{\rm out}$ can be greater than
$2R$, {\it i.\ e.}, the curving track never reaches the outer radius
obtained from the straight-line approximation. In those cases, the
outer radius is set to $2R$. This avoids arguments to the inverse sine
greater than one.

\subsection{Geometry of the CDC}

\begin{figure}
\centering\includegraphics[width=5.0in]{CDC_geometry.png}
\caption{Geometry of the CDC}\label{fig:CDC}
\end{figure}

Fig.~\ref{fig:CDC} shows a schematic of the CDC chamber displaying
some of the parameters that serve as input to the estimate:
\begin{description}
\item[$r_{\rm min}$] minimum radius of the CDC tracking volume
\item[$r_{\rm max}$] maximum radius of the CDC tracking volume
\item[$z_{\rm CDC}$] $z$-coordinate distance from center of target to
  downstream end of the CDC
\end{description}
Other parameters are:
\begin{description}
\item[$r_{\rm min,stereo}$] minimum radius of the CDC stereo layers
\item[$r_{\rm max,stereo}$] maximum radius of the CDC stereo layers
\item[$n_{\rm rl,FDC}$] number of radiation lengths measured transverse to
  the tracking layers, {\it i.~e.}, at $90^\circ$ in the lab ($n_{\rm
    rl} = x/X_0$)
\item[$n_{\rm m,FDC}$] number of position measurements for a track which
  passes through all layers of the CDC (sum of number of axial and
  stereo layers)
\item[$n_{\rm m,stereo}$] number of position measurements in stereo layers
\end{description}

For the CDC, $L_{\rm det} = r_{\rm max,CDC} - r_{\rm min,CDC}$ is the
radial thickness of the active volume. Then $L^\prime = L_{\rm det}$
and $L = L_{\rm det}/\cos\lambda$.

For tracks that exit the end of the CDC, the straight-line track
approximation is used to scale the number of measurements, the number
of radiation lengths and the transverse length of the track
($L^\prime$). In fact the straight-line approximation is used to
determine whether this scaling is necessary or not.

\subsection{Geometry of the FDC}

\begin{figure}
\centering\includegraphics[width=5.0in]{FDC_geometry.png}
\caption{Geometry of the FDC.}\label{fig:FDC}
\end{figure}

Fig.~\ref{fig:FDC} shows a schematic of the FDC chamber displaying
some of the parameters that serve as input to the estimate:
\begin{description}
\item[$z_{\rm min}$] minimum $z$ of the FDC tracking volume
\item[$z_{\rm max}$] maximum $z$ of the FDC tracking volume
\item[$r_{\rm FDC}$] outer radius of the FDC
\end{description}
Other parameters are:
\begin{description}
\item[$n_{\rm rl,FDC}$] number of radiation lengths measured transverse to
  the tracking layers, {\it i.~e.}, in the forward direction ($n_{\rm
    rl} = x/X_0$)
\item[$n_{\rm m,FDC}$] number of one-dimensional position measurements
  for a track which passes through all layers of the
  FDC\footnote{Since not all layers of the FDC will yield position
    measurements in a particular chosen dimension, the number of
    effective measurement planes will be only a fraction of the total
    measurement planes present in the chamber. While this is true of
    the CDC as well, the stereo angle is assumed to be so small that
    stereo layers are assumed to contribute fully to position
    determination in the $r$-$\phi$ plane (an approximation).}
\end{description}

For the FDC $L_{\rm det} = z_{\rm max} - z_{\rm min}$ is the length of
the active volume along the beam direction. Then
$L^\prime = L_{\rm det}/\tan\lambda$ and $L = L_{\rm det}/\sin\lambda$.

For tracks that exit the side of the FDC, like the CDC, the
straight-line track approximation is used both to determine if the
track is ``exiting early'' and to scale the number of measurements,
the number of radiation lengths and the transverse length of the track
($L^\prime$).
\subsection{Combining the CDC and the FDC}

For each quantity, the measurements in the CDC and FDC are
statistically independent. They are therefore combined using:
\begin{equation}
\sigma_{\rm total} = 
{1
\over
\sqrt{{1 \over \sigma_{\rm CDC}^2} + {1 \over \sigma_{\rm FDC}^2}}
}
\end{equation}
to get the overall resolution.

\section{Results}

Table~\ref{tab:figs} is a guide to plots which show results of
estimates of resolutions for the parameters in the design at this
writing. In each figure, the upper left graph shows resolution
contributions from the CDC, the upper right shows those from the FDC,
and the lower left shows the combination of the CDC and FDC. The lower
right contains a legend. In each of the upper plots (of the CDC and
FDC alone), the contribution from position resolution is shown in
yellow and that from multiple scattering is shown in blue. In
addition, for the azimuthal resolution figures, there is a
contribution from curvature error shown in magenta. CDC total
resolution is displayed in red. This red curve is the same data for
the CDC alone and in the combined. Likewise, the blue curve in the FDC
alone graph and in the combined is the total resolution from the FDC
and displays the same data in both cases.  In addition, in each of the
combined graphs at lower left there is a comparison with the results
of HDGEANT for the same angle or momentum as appropriate.

\begin{table}[h]
\centering\begin{tabular}{|lccc|}
\hline
& \bf $\delta p_t/p_t$ & \bf $\delta\theta$ & \bf $\delta\phi$ \\
resolution as a function total momentum
  & Fig.~\ref{fig:pscan-pt}
  & Fig.~\ref{fig:pscan-polar}
  & Fig.~\ref{fig:pscan-azimuth} \\ 
resolution as a function of polar angle
  & Fig.~\ref{fig:thetascan-pt}
  & Fig.~\ref{fig:thetascan-polar}
  & Fig.~\ref{fig:thetascan-azimuth} \\ 
\hline
\end{tabular}
\caption{Figures displaying resolution estimates.}\label{tab:figs}
\end{table}

\begin{figure}[ht]
\centering\includegraphics[height=4.1in]{figure/pscan_pt_20.png}
\caption{Estimates for resolution in relative transverse momentum as a
  function of total momentum at $20^\circ$.}\label{fig:pscan-pt}
\end{figure}

\begin{figure}[ht]
\centering\includegraphics[height=4.1in]{figure/pscan_polar_20.png}
\caption{Estimates for resolution in polar angle as a
  function of total momentum at $20^\circ$.}\label{fig:pscan-polar}
\end{figure}

\begin{figure}[ht]
\centering\includegraphics[height=4.1in]{figure/pscan_azimuth_20.png}
\caption{Estimates for resolution in azimuthal angle as a
  function of total momentum at $20^\circ$.}\label{fig:pscan-azimuth}
\end{figure}

\begin{figure}[ht]
\centering\includegraphics[height=4.1in]{figure/thetascan_pt_1_0.png}
\caption{Estimates for resolution in relative transverse momentum as a
  function of polar angle at $20^\circ$.}\label{fig:thetascan-pt}
\end{figure}

\begin{figure}[ht]
\centering\includegraphics[height=4.1in]{figure/thetascan_polar_1_0.png}
\caption{Estimates for resolution in polar angle as a
  function of polar angle at $20^\circ$.}\label{fig:thetascan-polar}
\end{figure}

\begin{figure}[ht]
\centering\includegraphics[height=4.1in]{figure/thetascan_azimuth_1_0.png}
\caption{Estimates for resolution in azimuthal angle as a
  function of polar angle at $20^\circ$.}\label{fig:thetascan-azimuth}
\end{figure}

The plots show reasonable agreement with the HDGEANT results.

\appendix

\section{Some useful equations}\label{app:algebra}

This appendix contains some intermediate results in the calculation of the error on the slope in the case of equally spaced measurements.

The well-known expressions
\begin{equation}
\sum_{i=1}^n i = {n(n+1) \over 2}
\end{equation}
and
\begin{equation}
\sum_{i=1}^n i^2 = {n(n+1)(2n+1) \over 6}
\end{equation}
allow us to obtain the following results for equally space measurements:
\begin{equation}
\sum x_i = {nL \over 2}
\end{equation}
\begin{equation}
\sum x_i^2 = {L^2n(2n-1) \over 6(n-1)}
\end{equation}
and for the parameter in the straight-line-fit error expressions we get
\begin{equation}
\Delta^\prime = {L^2n^2(n+1) \over 12(n-1)} .
\end{equation}


\bibliography{rezest}

\end{document}
