

\section{CP Violation and Mixing}
\label{sect:cpviolation}
\index{CP violation}

This section discusses how CP violation and mixing work
in the generator. These two topics are closely related, since
a large fraction of CP violating decays occurs via mixing. Exactly what is 
the best way of introducing these two effects into EvtGen is still not clear; 
there are different ways of generating CP violating decays depending on what 
one is interested in.  For example, one might want to generate a sample of events with 
$B\rightarrow a_1^+ \pi^-$, and the $B$'s which tag the flavor of this $B$ 
(e.g., via semileptonic decays) are all $B^0$'s.  On the
other hand, one might want to generate a sample that has the 
correct mixture of $B^0$ and $\bar B^0$ tags, as determined by the decay 
dynamics.  From the point of view of implementation these two situations are 
fairly different, and at the moment we are still working on choosing 
a consistent way of dealing with them.

This section starts out describing how mixing works, and
then proceeds to describe how CP violation is introduced.  Examples
will be provided. CP violating decays are
among the most complicated things this generator does, and 
various ways in which one might get the generation of CP
violating decays wrong will be pointed out.

\subsection{Mixing in the $B^0\bar B^0$ System.}

As is well know, when $\Upsilon(4S)$ decays into the $B^0\bar B^0$ system, 
the $B$'s are produced in a coherent state.  This implies that  
when one of the $B$'s decays, and its flavor is determined, the other
$B$ is known to be in a pure state of the opposite flavor.
This state then evolves (oscillates) independently, with the oscillation frequency
given by $\Delta m /2$, ($\Delta m$ being $m_{2} - m_{1}$, the mass difference between the two $B$ mass 
eigenstates).  Assuming $|B_{0}^{\rm phys}(t = 0)\rangle = |B^{0}\rangle$, its time evolution is given by: 
\begin{equation}
|B_{\rm phys}^{0}(t)\rangle = e^{-\Gamma t/2}e^{i\overline{m} t}\,\left[ |B^0\rangle\,\cos\left(\frac{\Delta m t}{2}\right) + i\frac{q}{p}|\bar B^0\rangle\,\sin\left(\frac{\Delta m t}{2}\right)\right],
\end{equation}
where $\overline{m} \equiv (m_{1} + m_{2})/2$, $q = e^{i\phi_{M}}/\sqrt{2}$, and $p = e^{-i\phi_{M}}/\sqrt{2}$ ($\phi_{M}$ is the mixing angle).  Time $t$ in this formula is the time between the two $B$ decays.

By projecting $|B_{0}^{\rm phys}(t)\rangle$ into a pure $B^0$
or $\bar B^0$ state, one can obtain the amplitudes for the $B$ to decay with
a given flavor.  In turn, squaring them gives the probabilities:
\begin{eqnarray}
P(B_{\rm phys}^{0}\rightarrow B^0) & = & \frac{1}{2} e^{-\Gamma t} (1+\cos(\Delta m t))\\
P(B_{\rm phys}^{0}\rightarrow \bar B^0) & = & \frac{1}{2} e^{-\Gamma t} (1-\cos(\Delta m t))\\
\end{eqnarray}

Integrating these over all times, one obtains the total rates:
\begin{eqnarray}
P(B_{\rm phys}^{0}\rightarrow B^0)  & = & \frac{x_{d}^2}{2\,(1 + x_{d}^{2})}\\
P(B_{\rm phys}^{0}\rightarrow \bar B^0)  & = & \frac{2 + x_{d}^2}{2\,(1 + x_{d}^{2})}\\
\end{eqnarray}
where $x_d \equiv \frac{\Delta m}{\Gamma}$.

Considering in a similar way the time evolution of the coherently produced state ($\Psi(t = 0) = \frac{1}{\sqrt{2}}\,\left [ |B^{0} \overline{B}^{0} \rangle - |\overline{B}^{0} B^{0} \rangle \right ]$ ), one can obtain the ratio:
\begin{equation}
{ N(\Upsilon(4S)\rightarrow B^0B^0)+N(\Upsilon(4S)\rightarrow\bar B^0\bar B^0)
  \over N(\Upsilon(4S)\rightarrow B^0\bar B^0)}= \frac{x_{d}^{2}}{2 + \,x_{d}^{2}}.
\end{equation}

Currently, mixing is implemented in two decay models: 
$\tt VSS\_MIX$ and $\tt VSS\_BMIX$.  These models both take the mass
difference to be a parameter and use the width assigned in the particle
properties table, i.e., from the evt.pdl file in the BaBar environment.
These models differ in the number of parameters that
the user must supply. The $\tt VSS\_MIX$ model requires that the branching
fractions into mixed and unmixed final states be set by hand, and does not
require that they be consistent with the mass difference and lifetime.
The $\tt VSS_BMIX$ model calculates the branching fractions automatically
and is therefore safer and easier to use. One can use this model as follows:
\begin{verbatim}
Decay Upsilon(4S)
0.5  B+       B-        VSS;
0.5  B0       anti-B0   VSS_BMIX;
Enddecay
\end{verbatim}
Both models give identical results when the $\tt VSS\_MIX$ branching
fractions are set correctly.


\subsection{Mixing and Non-zero Lifetime Differences.}

{\it This section should probably be merged with the previous but 
is currently written separately, mostly due to historical reasons, as
this is only needed for $B_s$ mixing, which was added later.}

In this section $B^0$ and $\bar B^0$ will generically denote any
neutral meson flavor eigenstate. These states are not eigenstates
of the full hamiltonian and hence do not have a definite mass or
lifetime. The states that are eigenstates of the hamiltonian are
\begin{eqnarray}
|B_L\rangle & = & p|B^0\rangle +q|\bar B^0\rangle, \\
|B_H\rangle & = & p|B^0\rangle -q|\bar B^0\rangle.
\end{eqnarray}
The states $B_L$ and $B_H$ have masses, $m_L$ and $m_H$, and widths, $\Gamma_L$
and $\Gamma_H$, respectively such that 
\begin{equation}
|B_{L,H}(t)\rangle=e^{-(\Gamma_{L,T}/2+im_{L,T})t}|B_{L,H}\rangle
\end{equation}
where $B_{L,H}(t)$ is the state after evolving the state $B_{L,H}$ for
a time $t$.

This allows us to evaluate several matrix elements that are useful for
calculating the mixing rates
\begin{eqnarray}
\langle B^0|B^0(t)\rangle & = & {1\over 2p}(\langle B^0|B_L(t)\rangle
                                           +\langle B^0|B_H(t)\rangle)\nonumber\\
                          & = & {1\over 2}(e^{-(\Gamma_L/2+im_L)t}
                                          +e^{-(\Gamma_H/2+im_H)t}),\\
\langle \bar B^0|B^0(t)\rangle & = & {1\over 2p}(\langle \bar B^0|B_L(t)\rangle
                                           +\langle \bar B^0|B_H(t)\rangle)\nonumber\\
                          & = & {q\over 2p}(e^{-(\Gamma_L/2+im_L)t}
                                          -e^{-(\Gamma_H/2+im_H)t}),\\
\langle B^0|\bar B^0(t)\rangle & = & {1\over 2q}(\langle B^0|B_L(t)\rangle
                                           -\langle B^0|B_H(t)\rangle)\nonumber\\
                          & = & {p\over 2q}(e^{-(\Gamma_L/2+im_L)t}
                                          -e^{-(\Gamma_H/2+im_H)t}),\\
\langle \bar B^0|\bar B^0(t)\rangle & = & {1\over 2q}(\langle \bar B^0|B_L(t)\rangle
                                           +\langle \bar B^0|B_H(t)\rangle)\nonumber\\
                          & = & {1\over 2}(e^{-(\Gamma_L/2+im_L)t}
                                          +e^{-(\Gamma_H/2+im_H)t}).
\end{eqnarray}

Now it is straight forward to calculate the mixing rate. Let $\chi$ 
define the probability that a meson produced, at $t=0$, as a $B^0$ will
decay as a $\bar B^0$. It is now easy to see that 
\begin{eqnarray}
\chi & = &{\int_0^\infty |\langle \bar B^0|B^0(t)\rangle|^2dt \over
\int_0^\infty |\langle \bar B^0|B^0(t)\rangle|^2dt +
\int_0^\infty |\langle B^0| B^0(t)\rangle|^2dt }\nonumber\\
  & = & { x^2+y^2 \over x^2+y^2+\vert {p\over q} \vert^2(2+x^2-y^2)} 
\end{eqnarray}
where 
\begin{equation}
x\equiv {\Delta m \over \Gamma},\qquad y\equiv {\Delta\Gamma\over \Gamma}
\end{equation}
and
\begin{equation}
\Gamma\equiv {\Gamma_L+\Gamma_H \over 2},\qquad
\Delta\Gamma \equiv \Gamma_H-\Gamma_L,\qquad
\Delta m \equiv m_H-m_L.
\end{equation}
Similarly the probability that that a meson produced, at $t=0$ 
as a $\bar B^0$ will decay as a $B^0$ is given by
\begin{eqnarray}
\bar \chi & = &{\int_0^\infty |\langle B^0|\bar B^0(t)\rangle|^2dt \over
\int_0^\infty |\langle B^0|\bar B^0(t)\rangle|^2dt +
\int_0^\infty |\langle \bar B^0| \bar B^0(t)\rangle|^2dt }\nonumber\\
  & = & { x^2+y^2 \over x^2+y^2+\vert {q\over p} \vert^2(2+x^2-y^2)}.
\end{eqnarray}
In the limit of $|p/q|=1$ and $y=0$ we have 
\begin{equation}
P_{\pm}(t)=\chi =\bar\chi={x^2\over 2(1+x^2)}
\end{equation}
as expected.
The time distribution, arbitrary normalization, is given by
\begin{equation}
e^{-\Gamma_L t}(1+e^{\Delta\Gamma t}\pm 2e^{\Delta\Gamma t/2}\cos\Delta m t)
\end{equation}
The minus sign corresponds to the case of mixing, and the plus sign
to the case of no mixing.


\subsection{CP-violation.}

The ``first generation'' of models which were implemented in EvtGen
to deal with CP violation did not solve all the problems that are 
involved in simulating the physics of the $B\bar B$ system. To illustrate how 
these models work, let us consider the following example (the decay 
$B\rightarrow \pi^+\pi^-$):
\begin{verbatim}
Alias MYB B0
#
Decay Upsilon(4S)
0.70   MYB B0      VSS;
0.30   MYB B0B     VSS;
Enddecay
#
Decay MYB
1.0   pi+ pi-     SSS_CP alpha dm CP |A| arg(A) |Abar| arg(Abar);
Enddecay
End
\end{verbatim}

The model used in this example is {\tt SSS\_CP}.  Its first argument is {\tt alpha}, is the
relevant CKM angle (in this case it is $\alpha$); the second argument {\tt dm} is
the mass difference of the two $B$ mass eigenstates; the third
argument, {\tt CP}, specifies the CP of the final state, and is either $+1$ or $-1$.
The latter argument could in principle be determined from the particle properties, but at 
the moment EvtGen does not include the C and P quantum numbers in its particle properties
list (this should be fixed). The last 4 arguments are the (complex) amplitudes for 
$B^0$ and $\bar B^0$ to produce the final state.

In this example the $\Upsilon(4S)$ decays with unequal
probabilities to ``MYB'' and $B^0$  and to ``MYB'' and  $\bar B^0$. The 
$B^0$ and $\bar B^0$ simply decay as a $B^0$ and a $\bar B^0$,
respectively, according to the generic decay table.  The interesting
thing is how ``MYB'' decays.  In the {\tt SSS\_CP} model, ``MYB'' will produce a
$\pi^+ \pi^-$ final state, while the other $B$ will provide
a tag of the flavor that is listed in the decay of the $\Upsilon(4S)$ ($B^0$ in 70\% 
of the cases, and $\bar B^0$ in 30\% of the cases).  Therefore, in this example there will
be some number of $B^0$ tags and some number of $\bar B^0$ tags. 
  
If one would like to get a sample of events with only $B^0$ tags, 
one could do it in the following way:
\begin{verbatim}
Alias MYB B0
#
Decay Upsilon(4S)
1.00  MYB B0      VSS;
Enddecay
#
Decay MYB
1.0   pi+ pi-     SSS_CP dm alpha CP |A| arg(A) |Abar| arg(Abar);
Enddecay
End
\end{verbatim}

Here, the only available tag (as specified by the $\Upsilon(4S)$ branching fraction) is $B^{0}$.

When the $\Upsilon(4S)$ is decayed and the two $B$'s are produced,
they are assigned lifetimes which are sampled from a pure exponential
distribution, with a mean of the average lifetime specified
in the particle properties list. (If the decay model is one that
includes the effects of mixing, the time distribution is no longer
a pure exponential; see the previous section.  However, as will be
evident below, it doesn't matter what time distribution is used in
the decay of the $\Upsilon(4S)$.)  Since the lifetime distributions
of both the $B$ that decays to the CP channel and the tag $B$ are
no longer exponential, the {\tt SSS\_CP} model has to regenerate 
lifetimes of {\it both} $B$'s. However, the model does not change what
tag one of the $B$ provides. This is important, because it means that this model
cannot enforce the right fraction of $B^0$ and $\bar B^0$ tags.  {\bf The only way to control this is through the decay table}.

The {\tt SSS\_CP} model handles decays to a pair of scalar particles.
Of these, the most important one is probably $B \rightarrow \pi\pi$, but the model also 
works for any other decay with two scalars in the final state, e.g., $K_s\eta'$.
Other models which are implemented in a similar way are {\tt SVS\_CP} and {\tt STS\_CP},
which handle the decays to a scalar and a vector and to a scalar and a 
tensor, respectively. There is also a model for two vectors
in the final state: it is called {\tt SVV\_CP}. One complication with this
model is that there are three partial waves that contribute to the 
rate with different CP. 

All of these models have one feature in common: they generate the $B$ tags with the relative
fractions specified in the decay table. It seems logical that the next step in expanding 
the functionality of the generator should be implementing some way to generate the right 
fraction of $B^0$ and $\bar B^0$ tags according to the amplitudes that are specified for a 
particular decay channel.  Indeed, this would be a very natural thing to do, because given 
the amplitudes, the relative fraction of $B^0$ and $\bar B^0$ tags can be determined.  A new class of models which have been recently put into EvtGen does have this feature.  These new models are particularly important for decays into final states that are not
CP eigenstates, for example $B\rightarrow a_1^+\pi^-$, or decays with direct CP violation (e.g., due to the presence of penguins).

Let us consider in more detail what the problem is.  The time-dependent CP-asymmetry for
the decay $B \rightarrow f$, where $f$ is a CP eigenstate, is defined as:
\begin{equation}
A_f(t) = \frac{\Gamma(B^{0}(t) \rightarrow f) - \Gamma(\overline{B}^{0}(t) \rightarrow f)}{\Gamma(B^{0}(t) \rightarrow f) + \Gamma(\overline{B}^{0}(t) \rightarrow f)}
\end{equation}
If $|A(B^{0} \rightarrow f)|$ = $|A(\overline{B}^{0} \rightarrow f)|$, then the integrated 
asymmetry $\int_{-\infty}^{+\infty}A_f(t)dt = 0$, and the relative fraction of $B^{0}$ ($\bar B^0$) tags is known (1/2).  However, if $|A(B^{0} \rightarrow f)| \neq |A(\overline{B}^{0} \rightarrow f)|$ (which can happen, e.g., because of penguin pollution), $\int_{-\infty}^{+\infty}A_f(t)dt \neq 0$, so the fraction of tags must be somehow determined.


One way to get the correct fraction of tags is to use integrated rates.  For example, for decays into a CP-eigenstate it can be defined as:
\begin{equation}
fr \equiv \frac{\int_{-\infty}^{\infty}\Gamma(\overline{B}(t)^{0} \rightarrow f)dt}{\int_{-\infty}^{\infty}\left[\Gamma(B^{0}(t) \rightarrow f) + \Gamma(\overline{B}^{0}(t) \rightarrow f)\right]dt} 
\end{equation}
so that
\begin{equation}
fr = \frac{|\overline{A}_{f}|^2 \left(1 + |\overline{r}_{f}|^2 + \frac{(1 - |\overline{r}_{f}|^2)}{1+x_{d}^2} \right)}{|\overline{A}_{f}|^2 \left(1 + |\overline{r}_{f}|^2 + \frac{(1 - |\overline{r}_{f}|^2)}{1+x_{d}^2} \right) + |{A}_{f}|^2 \left(1 + |r_{f}|^2 + \frac{(1 - |r_{f}|^2)}{1+x_{d}^2} \right)},
\end{equation}
where $x_{d} \equiv \frac{\Delta m}{\Gamma}$, $r_{f} = e^{2i\,\phi_{M}}\,\frac{\overline{A}}{A}$, and $\overline{r}_{f} = \frac{1}{r_{f}}$.

Two functions in the new version of EvtOtherB return this fraction (fractB0CP for decays into a CP eigenstate and fractB0nonCP for decays into non-CP eigenstates).  However, analytical calculations of the integrals can be done only in the simplest cases.  For example, for the $B \rightarrow \rho \pi$ decay one would have to integrate over the Dalitz plot as well as over time to get the total rates, which is very hard to do!  A better solution is to use the so-called acceptance-rejection method, which is effectively equivalent to numerical integration.  This method has been used in the ``second generation'' models which introduce direct CP-violation into EvtGen.

  The new models that are available are: {\tt SSS\_CP\_PNG} for the $B \rightarrow \pi^{+} \pi^{-}$ decay with penguin contributions; {\tt BTO2PI\_CP\_ISO} for the three isospin-related $B\rightarrow \pi \pi$ modes; {\tt BTOKPI\_CP\_ISO} for the four isospin-related $B \rightarrow K \pi$ modes; and {\tt SVS\_CP\_ISO} for the five isospin-related generic $scalar \rightarrow (pseudo)scalar + vector$ modes.  There is another decay model for the latter case, which also returns the correct fraction of $B^0$ ($\bar B^0$) tags.  It is called {\tt SVS\_NONCPEIGEN}.  Here is an example of how this model can be used for the $B \rightarrow a_1 \pi$ decay:
\begin{verbatim}
Alias MYB B0
Decay Upsilon(4S)
1.000 MYB anti-B0                         VSS;
Enddecay
Decay MYB
1.000 a_1-  pi+          SVS_NONCPEIGEN 1.22 0.51e12 0.0 
                                                       1.0 0.0 3.0 0.1 
                                                       3.0 0.1 1.0 0.0;
Enddecay
End
\end{verbatim}

The first parameter of this model is the corresponding CKM angle (in this case, it is $\alpha$); the second parameter is $\Delta m$; the third parameter is the ``flip'' which determines the fraction of ``MYB'' $\rightarrow f$ to ``MYB'' $\rightarrow \overline{f}$ decays, where the state specified in the decay table is considered the ``$f$'' state.  If the flip is set to 0, ``MYB'' always decays into the $f$ state.  The next four parameters are the absolute values and phases for the amplitudes $B^{0} \rightarrow f$ and $\bar B^{0} \rightarrow f$, respectively, and the last four parameters are those for the amplitudes $B^{0} \rightarrow \overline{f}$ and $\overline{B}^{0} \rightarrow \overline{f}$.

It is important not to confuse the ``flip'' parameter with the number of tags.  For example, one would like to generate 10,000 $B\overline{B}$ events with the flip parameter set to 0.5.  This means that one would get 5,000 ``MYB'' $\rightarrow f$ decays and 5,000 ``MYB'' $\rightarrow \overline{f}$ decays, but within each 5,000 event sample the number of $B^{0}$ and $\bar B^{0}$ tags is determined by the amplitudes specified in the decay file! 

As before, ``MYB'' in this example decays into a CP-mode, while the other B provides a tag.  However, the model itself generates the correct fraction of $B^{0}$ ($\bar B^{0}$) tags, which is determined from the last eight parameters (the amplitudes), and is independent of the $\Upsilon(4S)$ branching fractions.



\subsection{Special CP models}

 The models described above are meant to be generic models that can
be applied to different decays with appropriate arguments. All of the
above models are two body decays. There is, however, many other decays
that are not two body decays and these can not be treated quite as
simply since there will in general be interference between various
resonances. As an example consider the decay $B\rightarrow \pi^+\pi^-\pi^0$
which is dominated by the $\rho(770)$ resonances.

 To handle the $B\rightarrow 3\pi$ and the $B\rightarrow 4\pi$
decays two special models have been written. These models are 
called {\tt BTO3PI\_CP} and {\tt BTO4PI\_CP}. The second
of model {\tt BTO4PI\_CP} needs more work interms of understanding the
conventions of the phases between the resonances that contributes.

\newpage

\subsection{Unified implementation of CP-violating and mixing}

This section describes a {\it proposal} for how to unify and 
simplify the simulation of CP-violating decays and mixing in EvtGen, this 
should allow writing models that works independently of
how the neutral $B$ meson is produced, i.e. in a coherent pair in an
$\Upsilon(4S)$ decay or incoherently e.g. in a $p\bar p$ collision.

\subsubsection{CP violation}

The decays that we are considering here are 
common final states to a neutral $B$ meson and its antiparticle,
e.g. $B\to J/\psi K_S$
or $B\to J/\psi K^*$. At the $\Upsilon(4S)$ in the absence of
direct CP violation these decays can be simplified as we can make the 
assumption that you have an equal number of $B^0$ and $\bar B^0$
tags, or equivalent that there is no time integrated asymmetry.
However, we want to write the code such that no assumptions are
made about the numbers of $B^0$ and $\bar B^0$ tags.
For $B$ mesons that are produced incoherently we always have to
consider the asymmetries in the rates for an initial $B^0$ or
$\bar B^0$. 

To allow generating these asymmetries correctly these models
will be allowed to modify the flavor of the $B$ meson that is decayed.
For the $\Upsilon(4S)$ this means changing the the flavor of 'the
other' $B$ meson in the $\Upsilon(4S)$ decay, the tag $B$. 
For incoherent produced $B$
mesons this means that if a specific flavor is given to 
EvtGen it can return a decay of the opposite flavor. If
e.g. EvtGen is given a $B^{**}+$ to decay and it decays it to
$B^0\pi^+$ the decay of the $B^0$ might decide to flip the 
flavor of the decay chain to allow generating the right 
time integrated asymmetry. This means that after EvtGen is
done the code that called EvtGen will get back a $B^{**}-$
and needs to be prepared to handle this, e.g. by applying CP
to the decay chain that produced the $B^{**}+$. Also note the
meaning of branching fractions in the decay table, for modes
that are common to a $B$ and its anti-particle, the $\bar B$
the above implies that the branching fraction listed for the
$B$ and the $\bar B$ are the average branching fractions.

Practically, this is implemented in the {\tt EvtCPUtil}
class, which will handle the cases of coherent and incoherently
produced $B$ mesons.

How is this different from the current implementation, V00-09-38?
Not very much, there are already some models, e.g. the 
{\tt SVS\_NONCPEIGEN}
that has the full functionality that allows for direct CP violation 
and will generate the right mixture of $B$ and $\bar B$ tags.
All CP models should convert to this more general model. For 
applying this to $B$ mesons that are produced incoherent, i.e.
with a definite flavor, we need to modify the {\tt EvtCPUtil::otherB()}
method to allow the change of the flavor.

\subsubsection{Mixing}

In the $\Upsilon(4S)$ system mixing is handled by the decay of the
$\Upsilon(4S)$ meson via models such as {\tt VSS\_BMIX}. This 
implementation seems sufficient at this point. Note that we
are not considering the case of the $\Upsilon(5S)$ at this point.

For incoherently produced neutral $B$ mesons, $B^0$ or $B_s$, we assume that
the flavor is tagged and mixing will be generated by adding a 
decay like $B^0\to \bar B^0$ where the mixed particle have 
zero lifetime. Models that simulate CP violation, or in general
models that handle decays common to the $B$ and the $\bar B$ 
decay will remove any mixing, this is done in {\tt EvtCPUtil}. 
(Models that do not make any assumption about
the flavor of the $B$ should probably be using a 'neutral $B$' particle
instead of a specific flavor, but at BABAR we have not implemented 
this as it would be very confusing at this point to a lot of
analysis code that uses Monte Carlo truth.)

This means that models that simulates CP-violating decays should not
need to turn of mixing, it will practically be ignored. 

The mixing of incoherently produced  $B$ mesons is not done in a 
model, this causes some problems as it is hard to control the
mixing, e.g., turning it of or modifying the parameters. We need to
invent a way to control this.


