\section{Algorithm} \label{sect:algorithm} \index{algorithm} To illustrate how the event selection algorithm works consider the decay $B\rightarrow D^{*} \tau \bar\nu$, $D^* \rightarrow D\pi$, and $\tau \to \pi \nu$. The general case is a straight forward generalization of this example. The decay amplitude can be written as \def\AmpB{A^{B\rightarrow D^*\tau\nu}_{\lambda_{D^*}\lambda_{\tau}}} \def\AmpDstar{A^{D^*\rightarrow D\pi}_{\lambda_{D^*}}} \def\Amptau{A^{\tau\rightarrow \pi\nu}_{\lambda_{\tau}}} \begin{equation} A=\sum_{\lambda_{D^*}\lambda_{\tau}}\AmpB \times \AmpDstar \times \Amptau, \label{Eq:ampeq} \end{equation} where $\lambda^{\phantom '}_{D^*}$ and $\lambda^{\phantom '}_{\tau}$ label the states of spin degrees of freedom of the $D^{*}$ and the $\tau$, respectively. Thus, $\AmpB$ represents the decay amplitude for $B \to D^* \tau \nu$ for the six different combinations of $D^*$ and $\tau$ states. A possible implementation of Eq.~\ref{Eq:ampeq} is to generate kinematics according to phase space for the entire decay chain and to calculate the probability, the amplitude squared, which is used in an accept-reject algorithm. This approach has two serious limitations. First, the maximum probability of the decay chain must be known. This is logicistally difficult given the large number of potential decay chains in $B$ decays. Second, for long decay chains the accept-reject algorithm can be very inefficient as the entire chain must be regenerated if the event is rejected. We have implemented an algorithm that generates a decay chain as a sequence of sub-decays, thus avoiding both of these limitations. First the decay of the $B$ is considered. Kinematics are generated according to phase space and the probability is calculated \begin{equation} P_B=\sum_{\lambda^{\phantom '}_{D^*}\lambda^{\phantom '}_{\tau}}|\AmpB|^2. \end{equation} The kinematics are regenerated until the event passes an accept-reject algorithm based on $P_B$. After decaying the $B$ we form the spin density matrix \begin{equation} \rho^{D^*}_{\lambda^{\phantom '}_{D^*}\lambda^{'}_{D^*}}=\sum_{\lambda^{\phantom '}_{\tau}}\AmpB [{A^{B\rightarrow D^*\tau\nu}_{\lambda^{'}_{D^*}\lambda^{\phantom '}_{\tau}}}]^*, \end{equation} which describes a $D^*$ from the $B \to D^* \tau \nu$ decay after summing over the degrees of freedom for the $\tau$. To generate the $D^* \to D \pi$ decay, proceed as with the $B$, including also $\rho^{D^*}$ \begin{equation} P_{D^*}= {{1} \over {{\rm Tr}\,\rho^{D^*}}} \sum_{\lambda^{\phantom '}_{D^*} \lambda^{'}_{D^*}} \rho^{D^*}_{\lambda^{\phantom '}_{D^*}\lambda^{'}_{D^*}} \AmpDstar [{A^{D^*\rightarrow D\pi}_{\lambda^{'}_{D^*}}}]^*, \label{Eq:probdstar} \end{equation} where the scale factor, $1/{\rm Tr}\,\rho^{D^*}$, is proportional to the decay rate, and does not affect the angular distributions. This scale factor makes the maximum decay probability of each sub-decay independent of the full decay chain. Finally, we decay the $\tau$. We form the density matrix \begin{equation} \tilde\rho^{D^*}_{\lambda^{\phantom '}_{D^*}\lambda^{'}_{D^*}}=\AmpDstar [{A^{D^*\rightarrow D\pi}_{\lambda^{'}_{D^*}}}]^*, \end{equation} which encapsulates the information about the $D^*$ decay needed to properly decay the $\tau$ with the full correlations between all kinematic variables in the decay. Using the $\tilde\rho^{D^*}$ matrix we calculate the spin density matrix of the $\tau$ \begin{equation} \rho^{\tau}_{\lambda^{\phantom '}_{\tau}\lambda^{'}_{\tau}}=\sum_{\lambda^{\phantom '}_{D^*}\lambda^{'}_{D^*}} \tilde\rho^{D^*}_{\lambda^{\phantom '}_{D^*}\lambda^{'}_{D^*}} \AmpB[{A^{B\rightarrow D^*\tau\nu}_{\lambda^{'}_{D^*}\lambda^{'}_{\tau}}}]^*. \end{equation} As in the other decays, kinematics are generated according to phase space and the accept-reject is based on the probability calculated as in Eq.~\ref{Eq:probdstar}, replacing $D^*$ with $\tau$. The algorithm was illustrated above using an example which should convey the idea. In general consider the decay \begin{equation} A\rightarrow B_1B_2...B_N \end{equation} where the amplitudes are denoted by \begin{equation} A^{A\rightarrow B_1B_2...B_N}_{\lambda_A\lambda_{B_1}\lambda_{B_2}...\lambda_{B_N}}. \end{equation} The probability for, $P_A$, for this decay is given by \begin{equation} P_A=\sum_{\lambda_A\lambda'_{A}\lambda_{B_1}...\lambda_{B_N}} \rho^{A}_{\lambda_A\lambda'_{A}} A^{A\rightarrow B_1B_2...B_N}_{\lambda_A\lambda_{B_1}\lambda_{B_2}...\lambda_{B_N}} [A^{A\rightarrow B_1B_2...B_N}_{\lambda'_{A}\lambda_{B_1}\lambda_{B_2}...\lambda_{B_N}}]^* . \end{equation} The forward spin-density matrix $\rho^{B_i}$, given that $B_j$, $j