BMAD refers to a subroutine library\cite{Sagan:2016} %D. Sagan, Nucl. Instrum. & Methods A (2006)
for simulating the dynamics of relativistic beams of charged particles, and an associated 
format for defining beam line components. As long as the guide field components are defined according to BMAD conventions,
the full complement of the analysis tools of the library can be used to investigate the particle dynamics,
including various tracking methods, (Taylor maps, symplectic integration, Runge Kutta), construction of
maps from tracking, evaluation of Twiss parameters, tracking through maps defined by maps or multipoles, etc.
(For details see \href{https://www.classe.cornell.edu/bmad/manual.html}{Bmad-Manual}).

The BMAD formatted representation of the g-2 experiment
is comprised of three distinct branches, 1)M5 beam line, 2)the injection channel
through the backleg iron and inflector, and the 3)storage ring including
DC magnetic field, time dependent quadrupole electric fields and kicker magnetic
fields, and collimator apertures.
In the language
of BMAD, 'branches' are
joined at 'fork' elements and offsets with 'patches'.

The beam lines are assembled as a sequence of elements with fixed length.
The electro-magnetic fields in each element are defined
by field maps,  multipole expansions, or analytic expressions.
Time dependence for pulsed kickers and ramped quadrupoles requires
custom code (as there is no appropriate standard BMAD module).
Specialized custom routines can apply  arbitrary kicker pulse shapes
or quadrupole voltage time dependence.

Particle phase space and spin are tracked through the guide field
using runge-kutta integration.

The program is equipped with a module to generate an initial distribution of muons, or
to read a distribution from a file or both. For example, it is possible to generate a distribution gaussian in transverse phase space,
with longitudinal distribution read from a file. (This is useful for imprinting the 'W' shaped longitudinal distribution on an otherwise gaussian transverse distribution.)

Lattice and phase space parameters are computed with standard BMAD routines,
including twiss parameters, beam moments, closed orbits, betatron tunes, chromaticity,
synchrotron radiation, spin tunes, etc. The Jacobian matrices required to determine
tunes and twiss parameters are based on tracking.

Scattering is somehwat ad hoc. Muons entering the storage volume on trajectories that
pass through quad plates, multple scatter in the plates. Muons that are already inside
the storage volume are lost if their trajectory subsequently intercepts a quad plate or collimator.
If fiber monitors are active, muons multiply scatter in the fibers with energy loss consistent with
measurements.

If the muon decay is turned on, then muons decay with lifetime assigned at birth.
The phase space, spin coordinates and time of the decay at the decay point are recorded.
A positron is created with momentum four-vector appropriately correlated with muon spin.
The positron is tracked until  it strikes a calorimeter plane, where phase space coordinates
and time are recorded, or otherwise leaves the storage
volume.

The walls of the inflector are defined as continuous apertures. Particles that hit the
walls are lost. The ends of the inflector are treated as scattering planes.

\subsection{Quadrupoles}
The quadrupole plates are curved to preserve a fixed radial offset with respect to the
central on momentum design orbit. The curvilinear coordinate system can be represented by
beginning with a full 3-dimensional map, and then extracting an azimuthal
slice of the map (in $r,z$ at fixed $\phi$) or with a fitted multipole expansion of McMillan functions.
The curvature necessarily introduces nonlinearities that are not faithfully included in a 2 dimensional cartesian
expansion. The BMAD code allows several possible ways to specify the quadrupole electric field.
\begin{enumerate}
\item Field map. The field map that is currently incorporated is based on a $\rho, z$ azimuthal slice
  of of a 3 dimensional OPERA map. The quad map can be optionally defined as the sum of individual
  maps for each of the four plates, thus providing the flexibility to set the voltage independently on each plate,
  to model voltage errors, or the scraping configuration.
\item Multipole expansion. The multipole expansion in McMillan functions is derived as a fit to the azimuthal slice
  of the 3 dimensional map. An expansion in cartesian coordinates (that is inconsistent with Maxwell in the curved
  geometry) is another option. The multipole representations are convenient for exploring dependence on
  nonlinearity.
\end{enumerate}

\subsection{Main magnet}
The main magnet is represented as a map or analytic function with uniform field. A unform radial component
can be specified. Measured field errors are incorporated analytically. The azimuthal dependence of the
error field is expanded as a solution to Laplace's equation in cylindrical coordinates in order to ensure
consistency with Maxwell equations.
\subsection{Injection channel}
The magnetic field through the hole in the back leg iron and cryostat and main magnet fringe field is
based on a 3 dimensional OPERA map. A distinct map is computed for the field in the inflector.
The fringe and inflector maps are superimposed as appropriate.

The BMAD model has been used extensively to explore dependencies (capture efficiency, CBO amplitude,
phase space distribution of captured beam, etc.)
on injection trajectory, inflector
orientation, field and aperture, kicker field amplitude, field profile, pulse shape, and timing, and quad voltage.
Effects of quad field errors and misalignment on closed orbit and twiss parameters have been studied.
Simulations informed analysis of measurements with fiber harps. The flexibility to employ various
quadrupole field profiles has enabled study of effects of nonlinearities on decoherence and debunching,
electric field and pitch corrections. Indeed simulation will play a critical role in evaluation of the electric field and pitch systematics.

