11institutetext: Pasquale Musella 22institutetext: ETH Institute for Particle
Physics and Astrophysics,
Otto-Stern-Weg 5, Zurich, Switzerland
22email<EMAIL_ADDRESS>33institutetext: Francesco Pandolfi
44institutetext: INFN - Sezione di Roma
Piazzale Aldo Moro 2, Rome, Italy.
44email<EMAIL_ADDRESS>
# Fast and accurate simulation of particle detectors using generative
adversarial networks

Pasquale Musella  Francesco Pandolfi

###### Abstract

Deep generative models parametrised by neural networks have recently started
to provide accurate results in modeling natural images. In particular,
generative adversarial networks provide an unsupervised solution to this
problem. In this work we apply this kind of technique to the simulation of
particle-detector response to hadronic jets. We show that deep neural networks
can achieve high-fidelity in this task, while attaining a speed increase of
several orders of magnitude with respect to traditional algorithms.

###### Keywords:

Generative Adversarial Networks Deep learning High Energy Physics Simulation
Fast simulation Jet images CERN open data.

##  1 Introduction

The extraction of results from high energy physics data crucially relies on
accurate models of particle detectors, and on complex algorithms that infer
the properties of incoming particles from signals recorded in electronic
sensors. Numerical models, based on Monte Carlo methods, are used to simulate
the interaction between elementary particles and matter.
In particular, the Geant4 toolkit geant4_2003  features state-of-the art
models and is employed to simulate particle detectors at the CERN LHC.
Reconstruction algorithms routinely used at collider experiments (see e.g.
ATLAS_exp  and CMS_exp ) are based on estimators of particle trajectories and
energy deposits. This information is subsequently aggregated in order to
reconstruct energy, type and direction of final state particles produced by
the collision of the primary beams.

The CERN LHC complex will undergo a series of upgrades HL_LHC  over the next
ten years that will allow collecting a dataset roughly 30 times larger than
the one currently available. The number of simultaneous interactions per bunch
crossing in such a future dataset will increase by a factor of about 4,
compared to the present levels. It is estimated HSF_whitepaper  that, because
of the larger volume and complexity of the data the shortfall between needs
and bare technology gains is about 4-fold in computing power, if one assumes
constant funding. This gap should therefore be bridged with faster, more
efficient algorithms for particle detector simulation and data reconstruction.

In this work, we develop a generative model parametrised by a deep neural
network, that is capable of predicting the combined effect of particle
detector simulation models and reconstuction algorithms to hadronic jets. The
results are based on samples of simulated hadronic jets produced in proton-
proton collisions at $\sqrt{s}=7$ TeV that was published by the CMS
collaboration on the CERN open data portal. The dataset qcd30 ; qcd50 ; qcd80
; qcd120 ; qcd170 ; qcd300 ; qcd470 ; qcd600 ; qcd800 ; qcd1000 ; qcd1400 ;
qcd1800  is part of the “level 3” category in the the High Energy Physics
(HEP) data preservation classification DPHEP  and it contains the result of
the Geant4 simulation of the CMS detector and of the subsequent data
reconstruction algorithms used by the CMS collaboration.

Generative adversarial networks (GANs) goodfellow_generative_2014  are pairs
of neural networks, a generative model and a discriminative one, that are
trained concurrently as players of a minimax game. The task of the generative
network is to produce, starting from a latent space with a fixed distribution,
samples that the discriminative model tries to separate from samples drawn
from a target dataset. It can be shown goodfellow_generative_2014  that with
this kind of setup the generator is able to learn the distribution of the
target dataset, provided that the generative and discriminative models have
enough capacity (it should be noted that this condition has actually been
proven to hold for mixtures of neural networks gan_mix , not for single ones;
the interested reader may find additional information about GANs convergence
in gan_mix ; f_gan ; w_gan1 ; w_gan2  and references therein.

Since they were first proposed, GANs have been applied to an increasingly
large number of problems in machine learning, mostly dealing with natural-
image data, but not only. Applications of adversarial networks were also
proposed in the context of HEP, mainly with two purposes: training of robust
discriminators that are insensitive to systematic effects, or uncorrelated
from observables used for signal extraction Louppe:2016ylz ; Shimmin:2017mfk ;
nips17_systs , and for event generation and detector simulation
deOliveira:2017pjk ; Paganini:2017hrr ; nips17_calo . Similar applications
were proposed in the context of cosmic ray experiments Erdmann:2018kuh .

We use GANs to train a generative model that generates the reconstruction-
level detector response to a hadronic jet, conditionally on its particle level
content. We represent hadronic jets as “gray-scale” images of fixed size
centred on the jet axis, where the pixel intensity reflects the fraction of
jet energy deposited in the corresponding geometrical cell. The architecture
of the networks and the problem formulation, that can be classified as a
domain mapping one, are based on the image-to-image translation work described
in pix2pix . We introduce a few differences to tailor the approach to the
generation of jet images: we explicitly model the set of non-empty pixels in
the generated images, which are much sparser than in natural images; we
enforce a good modelling of the total pixel intensity, though the combined use
of feature matching improved_gans  and of a dedicated adversarial classifier;
and we condition the generator on a number of auxiliary features. The use of
feature matching on physics-inspired features was already introduced in the
HEP context in Paganini:2017hrr , while conditional GANs are common in image
generation problems (see, e.g. cgan  and references in pix2pix ). Previous
work on the handling of sparsity in hadronic jet images was proposed in
Paganini:2017hrr  and deOliveira:2017rwa , and it was based on the engineering
of high level auxiliary features.

Related work has recently been presented in Paganini:2017hrr  and
Erdmann:2018kuh . The authors of Paganini:2017hrr  proposed to use a deep
convolutional neural network to simulate calorimeter showers, thus aiming at
modelling the particle interaction with the detector medium. The solution that
we explore here allows the largest reduction in computation time, by
predicting directly the objects used at analysis level, and thus reproducing
the output of both detector simulation and reconstruction algorithms. This
philosophy is similar to that of the parametrised detectors simulations
delphes  that are often used in HEP for phenomenological studies, and that are
very limited in accuracy. We show that using a deep neural network model
allows attaining accuracies that are comparable to that of the full simulation
and reconstruction chain. The approach of Erdmann:2018kuh , that studied the
application of GANs to the generation of air-showers, is more similar to ours,
as it aims at predicting the patterns reconstructed by the detectors,
conditionally on the energy and type of the primary particles.

##  2 Inputs and problem formulation

For this study we use simulated samples of hadronic jets produced by the CMS
collaboration and published on the CERN open data portal. In particular we
take hadronic jets produced in proton-proton collisions at $\sqrt{s}=7$ TeV.
These events feature state-of-the-art characteristics in terms of simulation
and reconstruction algorithms in HEP. The events were generated with the
PYTHIA6 event generator pythia6 , the CMS detector response was simulated
using Geant4 geant4_2003 . Concurrent proton-proton interactions (“pile-up”)
were simulated, roughly reproducing the LHC running conditions of 2011. The
samples contain the results of the full CMS reconstruction chain CMS_exp .

In the input dataset, hadronic jets were clustered with the anti-kt algorithm
antikt , using the FastJet library fastjet  and a distance parameter of 0.5.
We used two sets of hadronic jets: those clustered from the list of stable
particles produced by PYTHIA6, and those clustered from the list of
reconstructed particle candidates. In the following we term “particle-level
jets” the former, and “reconstructed jets” the latter.

Following standard practices in collider experiments, we employ a cylindrical
system of coordinates. The origin of the coordinate system is set to the
centre of the CMS detector, the z axis is chosen to be parallel to the beam
line, and the x axis is chosen to point towards the centre of the LHC ring. We
indicate as $\phi$ and $\theta$ the azimuthal and polar angles, respectively,
and we define the pseudorapidity $\eta$ as $\log(\cot(\theta/2))$.

In each event, we select particle-level jets with a transverse component of
the momentum above 20 GeV and with an absolute value of the pseudorapidity
below 2.5. A search is then performed in the reconstructed jet collection to
find reconstructed jets satisfying
$\sqrt{\Delta\eta^{2}+\Delta\phi^{2}}<0.3$, where $\Delta\eta$ and
$\Delta\phi$ are, respectively, the difference between pseudo-rapidity and
azimuth of the particle-level and reconstructed jets. Pairs of reconstructed
and particle-level jets satisfying these conditions are considered in this
study. As we are interested in the particle content of jets, no jet energy
calibrations are applied to the reconstructed jets.

Jet images are constructed by opening a window of size
$\Delta\eta\times\Delta\phi=0.3\times 0.3$ around the particle-level jet axis.
The window is split into $32\times 32$ identical square pixels and the
intensity associated to each pixel is proportional to the total transverse
momentum of the jet components contained in it, divided by the transverse
momentum of the particle-level jet. Roughly 80-90% of the jet energy is
contained in the jet window that we considered, and the jet components not
contained in the window are used to fill the closest pixel of the image
border. The choice of studying the central part of the jet was aimed at
limiting the data size, while keeping sufficiently high spacial resolution.
Increasing the window to the $0.5\times 0.5$ region would have in fact almost
doubled the image size and thus considerably increased the training time. This
technical limitation will have to be addressed by future work, but is not
expected to significantly change the conclusion of this work.

Figure 1:  Distribution of the jet transverse momentum at particle level for
the hadronic jets in the input dataset. The blue histogram shows the
distribution as coming from the dataset, while the red one shows the
distribution obtained after applying the weights used for training.

Each pair of jet images is further associated to four auxiliary features: the
transverse momentum ($p^{gen}_{T}$), pseudo-rapidity ($\eta^{gen}$) and
azimuthal angle ($\phi^{gen}$) of the particle-level jet, and the number of
pile-up interactions ($n_{PU}$) that were simulated in the event under
consideration. No further preprocessing or standardisation was performed on
the input data, which is available on the Zenodo information server data .

The distribution of $p^{gen}_{T}$ for the selected hadronic jets is shown by
the blue histogram in figure 1. The shape of the distributions is multi-modal
because the original dataset was split in several bins in the scale of the
parton-parton interaction. A set of weighting factors, as a function of
$p^{gen}_{T}$ was applied to obtain a falling distribution. The distribution
obtained after applying such weighting factors is shown by the red histogram.

###  2.1 Notation and learning setup

We adopt the following notation: we denote by $\mathbf{x}$ and $\mathbf{y}$
the jet image at particle and reconstruction level, respectively, we indicate
with $\mathbf{c}$ the set of auxiliary features, and we use $\mathbf{z}$ to
denote a latent space of uniformly distributed noise. With this notation, our
problem can be formulated as follows.

  * •

Given:

    * $(\mathbf{c},\mathbf{x})\sim p_{cx}(\mathbf{c},\mathbf{x})$

    * $\mathbf{y}\sim p_{y}(\mathbf{y}|\mathbf{x},\mathbf{c})$

    * $\mathbf{z}\sim p_{z}(\mathbf{z})=U(\mathbf{z})$;

where $p_{cx}$ and $p_{y}$ are the input-data distributions, and $U$ indicates
the uniform distribution.

  * •

We want to construct a function $G$, such that

    * $\mathbf{y^{\prime}}=G(\mathbf{z},\mathbf{x},\mathbf{c})\sim p_{y}(\mathbf{y^{\prime}}|\mathbf{x},\mathbf{c})$.

We note that the $\mathbf{c}$ and $\mathbf{x}$ variables sets play, from a
mathematical point of view, an identical role in the problem, as we want to
condition the output of $G$ is conditional on the union of the two. The reason
to separate them is mostly conceptual, as $\mathbf{x}$ represents the image
data associated with a jet, while $\mathbf{c}$ parametrises information about
jet kinematics and the environment. Furthermore, the two sets of variables are
treated differently by the neural network architecture, as detailed in the
following.

The function $G$ is a generative model that approximates the combined response
of particle detector simulation and reconstruction algorithms to hadronic
jets. Following the GAN paradigm, we look for a solution to this problem by
introducing a discriminative model $D$ and setting-up a minimax game between
the two models, with value function $V(G,D)$ defined as:

$\begin{array}[]{ll}V(G,D)~{}=&\min_{G}\max_{D}\\{\mathbb{E}_{(\mathbf{c},\mathbf{x},\mathbf{y})}\left[\log(D(\mathbf{y},\mathbf{x},\mathbf{c}))\right]\\\\[4.30554pt]
&~{}~{}+~{}\mathbb{E}_{(\mathbf{c},\mathbf{x},\mathbf{z})}\left[\log(1-D(G(\mathbf{z},\mathbf{x},\mathbf{c}),\mathbf{x},\mathbf{c}))\right]\large\\}\end{array}$
(1)

While the problem could in principle be solved using the GAN setup alone, we
inject additional information in order to stabilise and speed-up the
convergence. In particular, we take into account two facts:

  1. 1.

$G$ and $\mathbf{y}$ should match on average;

  2. 2.

$\mathbf{y}$ is very sparse (on average, roughly 3% of the pixels have non-
zero values).

The authors of pix2pix  show that the first requirement can be efficiently
satisfied by adding an $L_{2}$\- or $L_{1}$-norm term to the loss function. We
adopt this approach, using in particular an $L_{2}$-norm term.

To explicitly take into account the second requirement, we modify the
structure of the generator, by increasing the depth of its output: one channel
is used to model the pixel intensity, while a second channel, to which we
refer as a “soft-mask”, models the probability of a pixel to be non-zero. We
denote the two channels as $G_{0}$ and $G_{1}$, and we modify the generative
model loss function by adding a term of this form:

$\begin{array}[]{ll}\lambda&{\rm{\cal
L}_{P}}(G)~{}=\lambda~{}\mathbb{E}_{\mathbf{y}=0}\left[-\log(G_{0}(\mathbf{z},\mathbf{x},\mathbf{c}))\right]\\\
&+~{}\lambda~{}\mathbb{E}_{\mathbf{y}>0}\left[-\log(1-G_{0}(\mathbf{z},\mathbf{x},\mathbf{c}))+\left(G_{1}(\mathbf{z},\mathbf{x},\mathbf{c})-\mathbf{y}\right)^{2}/2\right]\\\
\end{array}$ (2)

where $\lambda$ is the associated hyperparameter.

To generate images, we sample the soft-mask probabilities to create a “hard-
mask” binary stochastic layer
$G_{1}^{\prime}(\mathbf{z})=\mathbf{1}_{\mathbf{z}<G_{1}}$, where $\mathbf{1}$
is the so-called indicator function. The GAN value function in eq. 4 becomes
$V(G_{2},D)$, where $G_{2}$ is defined as $G_{0}\cdot G_{1}^{\prime}$, and the
differentiability is preserved by replacing $G_{1}^{\prime}$ with $G_{1}$
during back-propagation hinton_bsn .

Finally, we enforce a good modelling of the total image intensity, which is
proportional to the reconstructed jet energy, with two additions:

  * •

we add an extra term to the generative model loss function, proportional to
the mean squared error of the total image intensities:

$\tau{\rm{\cal
L}_{T}}(G_{2})~{}=\tau~{}\mathbb{E}\left[\left(I(G_{2}(\mathbf{z},\mathbf{x},\mathbf{c}))-I(\mathbf{y})\right)^{2}\right]$
(3)

where the operator $I$ computes the total intensity of the jet images and
$\mu$ is the associated hyperparameter.

  * •

We introduce a second discriminative model $D_{T}$ that receives as input the
total reconstructed jet image intensities, and the auxiliary features
$\mathbf{c}$, and we set-up an additional minimax game, whose importance is
controlled by the $\mu$ hyperparameter, with value function
$V_{T}(G_{2},D_{T})$, defined as below:

$\begin{array}[]{ll}\mu&V_{T}(G_{2},D_{T})=\\\\[4.30554pt]
&\mu\min_{G_{2}}\max_{D_{T}}\\{\mathbb{E}\left[\log(D_{T}(I(\mathbf{y}),I(\mathbf{x}),\mathbf{c}))\right]\\\\[4.30554pt]
&~{}~{}+~{}\mathbb{E}\left[\log(1-D_{T}(I(G_{2}(\mathbf{z},\mathbf{x},\mathbf{c})),I(\mathbf{x}),\mathbf{c}))\right]\large\\}\end{array}$
(4)

###  2.2 Model architecture

Figure 2:  Graphical summary of the generator network. Boxes shows the data
representation, while arrows represent operations. Input and output nodes
annotated with bold-text labels. On the right side of the diagram, white and
blue boxes are concatenated before applying up-convolution operations.
Convolutions applied to the $\mathbf{x}$ node and its daughters have a filter
size of 3x3, while a 1x1 convolution is applied to the $\mathbf{c}$ node. All
up-convolutions use 3x3 filters.

The generative model and the discriminative model $D$ are implemented as
convolutional neural networks lecun-98 . For the generator we adopt the so
called “U-net” architecture Unet , that consists of an encoding section
followed by a decoding one, with additional skip connections linking encoding
and decoding layers with the same spatial dimension. The input images are
first fed into a batch normalisation layer BN , which allows running the
network on non-standardised inputs. Afterwards, we use 5 encoding layers and 5
decoding ones. Each encoding layer consists of a convolutional unit, followed
by a batch-normalisation one and by a leaky ReLU activation LReLU , with a
slope of 0.2 in the negative domain. The encoding filters size is chosen to be
of 3$\times$3, with a stride of 2 in order to reduce the representation width.
The number of filters is set to 16 for the first layer and it is doubled at
each step. The decoding layers comprise a concatenation unit to implement the
skip connections, followed by up-convolutional units, batch normalisation and
leaky ReLU activation ones. Drop-out units are also employed in the first two
layers of the decoding section. At each step in the decoding section, the
depth of the representation is halved, while its width is doubled. This is
achieved by decreasing the number of convolutional filters, while
appropriately choosing the stride and padding parameters. Auxiliary
conditional features are injected in the architecture as follows: they are
first passed through a batch normalisation unit and then into a 1$\times$1
convolution unit whose output matches the depth of the last encoding layer;
the 1$\times$1 convolutional unit output is subsequently concatenated with
that of the last encoding convolution. Noise can be injected into the
architecture at the same level. However, we obtained better results by feeding
noise only in the form of a stochastic sampling of the output soft-mask.
Figure 2 shows a graphical summary of the generator architecture.

The discriminative model $D$ uses 4 layers, comprising 3$\times$3
convolutional filter units, batch normalisation units and leaky ReLU
activation ones. Stride and padding are tuned in such a way that the largest
field of view of the convolution layers is of 13$\times$13 pixels. The model
acts as a “patch-GAN” pix2pix , i.e. it is only sensitive to the local
structure of the jet images. The convolutional layers are followed by a fully-
connected layer with a sigmoid activation function. The auxiliary variables
are treated similarly to what is done in the generative model and are injected
at the input of the fully connected layer.

The model $D_{T}$ takes as input the total intensities for $\mathbf{y}$ (or
$G_{2}$) as well as $\mathbf{c}$, and is parametrised as a feed-forward fully
connected neural network with 4 layers, using dense units, batch normalisation
and leaky ReLU activations. The fully connected layers have widths of
64-64-32-16 and are followed by an output layer with a sigmoid activation
function.

Figure 3:  Average jet images obtained at the particle-level (“gen”), after
detector simulation and reconstruction (“reco”), and those predicted by the
generative model (“pred”).

We did not perform a formal optimisation of the neural networks architecture,
but we picked a particular set of values after exploring the parameter space
in terms of width and depth of the networks, based on two factors: the
performance of the model and the computational times required.

Figure 4:  Joint distribution of the particle-level (“gen”) and reconstruction
level pixel intensities. The left plot shows the distribution of the input
data, while the right plot is obtained using the generative model.

A software package implementing the neural networks’ instantiation and
training is openly available in code .

##  3 Results

The models were trained on two million jet images extracted from the dataset
described in section 2. The TensorFlow tf  framework, and the Keras keras
high level interface were used to implement and train the models. Computing
resources from the Piz Daint Cray supercomputer located at the Swiss Centre
for Supercomputing were used to obtain the results that we present here. Two
independent Adam adam  optimisers were employed for parameter sets of the
generative and discriminative models. NVIDIA Pascal P100 GPUs were used to
accelerate the computations and the models were trained for 10-20 epochs,
which were sufficient to achieve convergence. The training time for these
models was around 1 hour per epoch. Inference ran at roughly 100Hz on Intel
Xeon CPUs and at roughly 10kHz on NVIDIA Pascal P100 GPUs, which are to be
compared to the typical time scale for event simulation and reconstruction,
i.e. $10^{-1}$–$10^{-2}$ Hz.

Figure 5:  Aggregated pixel intensities for different rings in
$\Delta\eta$–$\Delta\phi$. Blue histograms are obtained from the input data,
while red ones are obtained using the generative model.

Figure 6:  Evolution of the aggregated pixel intensities for different rings
in $\Delta\eta$–$\Delta\phi$ as a function of the particle level jet
transverse momentum. Solid lines represent the median of the distribution,
filled regions show the inter-quartile range, while dashed lines mark the 10%
and 90% quantiles. Blue lines are obtained from the input data, while red ones
are obtained using the generative model.

Figure 3 shows the average jet image obtained at particle-level by aggregating
the $\mathbf{y}$ values, as well as those obtained at reconstruction-level by
aggregating either the $\mathbf{x}$ values or the $G_{2}$ ones. One can
observe that the average effect of the detector simulation and reconstruction
algorithms is to generally spread the jet energy away from the core. We see
that this general trend is correctly reproduced by our set-up.

Figure 4 shows the correlation between the particles and detector-level pixel
intensities for the input dataset and for the predicted images, obtained from
a test sample of 100000 images. Three sub-populations can be observed in the
distributions:

  * •

well measured jet components populate the diagonal;

  * •

errors in the position reconstruction of the jet components lead to energy
migration between close-by pixels and thus contribute to the sub-populations
located close to the horizontal and vertical axes;

  * •

non-reconstructed jet components manifest as empty reconstruction-level pixels
that correspond to non-empty ones at particle level and therefore to the sub-
population located along the horizontal axis

  * •

pile-up effects lead to non-empty pixels in the reconstruction level image in
correspondence to empty particle-level pixels and so contribute to the sub-
population close to the vertical axis.

As can be seen from the figure, our set-up (right panel) achieves a good
modelling of the relative weight of the three sub-populations, which result
from a non-trivial set of effects. The sub-populations located along the
diagonal, the horizontal axis and the vertical one account for about 40%, 25%
and 25% of the non-empty pixels, respectively and our algorithm is able to
reproduce such numbers with a relative accuracy of roughly 30%.

In figure 5 we show the distributions of the total pixel intensities obtained
integrating over rings in $\Delta R$ centred on the particle-level jets axis.
Blue histograms are obtained from the input dataset, while red ones show the
results of the generative model. Figure 6 shows the evolution of the
aggregated pixel intensities distribution as a function of the particle-level
jet transverse momentum. These results show that our set-up allows good
modelling of hadronic jet structure over more than two orders of magnitude in
jet transverse momentum.

We further investigate the goodness of the learned model by evaluating its
ability to reproduce high level jet features that are typically used in
physics analyses. We concentrate, in particular on two sets of variables:

  1. 1.

variables used in the context of quark/gluon discrimination;

  2. 2.

jet substructure variables used in the context of merged jets discrimination.

Figure 7:  Distribution of high level variables used for quark/gluon
discrimination (first two rows) and merged jets tagging (last row). Blue
histograms are obtained from the input data, while red ones are obtained using
the generative model.

Figure 8:  Evolution of the quark/gluon (first two rows) discrimination and
merged jet tagging (last row) variables as a function of the particle level
jet transverse momentum. Solid lines represent the median of the distribution,
filled regions show the inter-quartile range, while dashed lines mark the 10%
and 90% quantiles. Blue lines are obtained from the input data, while red ones
are obtained using the generative model.

From the first set, we choose the so-called major and minor axes, i.e. the
square root of the eigenvalues of the $\eta$-$\phi$ covariance matrix of the
jet image, and the $p_{T}D$ variable, i.e. the ratio between the square root
of the second and first non-central moment of the pixel intensities
Pandolfi2013 . From the second set, we choose the ratio between the 2- and
1-subjettines subjett  and that between the 2- and 3-subjettiness. The
subjettiness variables were computed using the FastJet package fastjet1 ;
fastjet2 , approximating each jet as a set of mass-less particles with
energies and directions obtained from the pixel intensities and positions.

Figure 7 shows the distribution of the quark/gluon and merged jets
discrimination variables that we considered aggregated over the test dataset,
while figure 8 shows the evolution of the distributions as a function of the
transverse momentum of the jet at particle level. The level at which these
variables are predicted by our set-up is good (in general, the marginal
densities and the value of the distributions quantiles are predicted with an
accuracy of 5-10%), even though some mismodelling can be observed for the
quark-gluon discriminating variables. In particular, mismodelling at the
10-20% level in the marginal density, and in the $p^{gen}_{T}$ dependence of
the 75% and 90% quantiles can be observed for the $p_{T}D$ and the major and
minor axes. These kind of disagreements point to the fact that the correlation
between the number of non-empty pixels and their energy sharing is not
perfectly modelled.

###  3.1 Discussion

The results that we discussed above represent a step forward in terms of
accuracy of fast simulation systems proposed in the context of collider
detector physics. We believe that three main aspects contributed to this:

  * •

the use of a generative model that is designed to handle spatial correlations
well, and the use of a conditioning space (i.e. that of particle-level images)
that encodes large amounts of spatial information;

  * •

the explicit handling of the sparsity through the soft-mask layer;

  * •

the use of physics-driven constraints on the total intensity of the jet
images.

The method that we outline here has potential application in conjunction with
fast detector simulation models, or parametrised ones delphes ; cms_fast_sim .
In this context, being able to accurately predict the output of simulation and
reconstruction algorithms for objects like hadronic jets, which are ubiquitous
at the LHC, would allow to save large amounts of the computing power, by
reducing the cost of producing simulated events samples.

While a relatively stable set-up was established by tuning the model hyper-
parameters, this aspect of the work is not yet completely satisfactory, as the
region of hyper-parameters space that lead to satisfactory results was found
to be relatively narrow. A review of the loss function structure, possibly
incorporating the use of alternative formulations of the GAN game f_gan ;
w_gan1 ; w_gan2 , and of the model training strategy in general, will be
important to allow streamlining the method, and will the subject of future
work.

Furthermore, the current approach introduces noise only through the stochastic
generation of the set of active pixels. Our attempts at injecting noise at a
more fundamental level in the generative model structure have been
unsuccessful so far. Similar problems have been reported by researchers
working on natural image generation (see e.g. pix2pix ). A more extensive
investigation of the handling of noise will also be subject of subsequent
work.

##  4 Conclusions

We have reported on a method that uses deep neural networks to learn the
response of particle detectors simulation and reconstruction algorithms. The
method is based on generative adversarial networks and it was applied to the
generation of hadronic jet images at the CERN LHC.

We trained a generative model to reproduce the combined response of state-of-
the-art simulation and reconstruction algorithms. This was possible thanks to
the exploitation of the open datasets published by the CMS collaboration under
the HEP data preservation initiative.

Starting from proposals made for natural image processing, we devised a hybrid
set-up based on the combined use of generative adversarial networks and
analytic loss functions that is able to take into account the conditioning on
auxiliary variables and physics-driven constraints on the generation process.

Our method allows reducing the computation time required to obtain
reconstruction-level hadronic jets from particle-level jets by several orders
of magnitude, while achieving a very good accuracy in reproducing the
simulation and reconstruction algorithms response. The model is in particular
capable of reproducing the evolution of the reconstructed jet shapes as a
function of several conditional variables. Physics-driven high level features
commonly used for merged jets tagging and quark/gluon discrimination, and
their evolution, are also generally well modelled.

The results obtained with this work represent a promising step forward towards
the development of fast and accurate simulation systems that will be crucial
for the future of collider experiments in high energy physics.

## Acknowledgements

We thank the CMS collaboration for publishing state of the art simulated data
under the open access policy. We strongly support this initiative and believe
that it will be crucial to spark developments of new algorithms from which the
HEP community as a whole can profit. We thank Dr. M. Donegà, Prof. G.
Dissertori, and Dr. M. Pierini, for their careful review of this manuscript.
This work was supported by the Swiss Centre for Supercomputing under the
project D78. The final publication is available at
https://doi.org/10.1007/s41781-018-0015-y.

## Conflicts of interest

On behalf of all authors, the corresponding author states that there is no
conflict of interest.

## References

  * (1) S. Agostinelli, et al., Geant4-a simulation toolkit, Nuclear Instruments and Methods in Physics Research Section A 506(3), 250 (2003).  DOI 10.1016/S0168-9002(03)01368-8
  * (2) G. Aad, et al., The ATLAS Experiment at the CERN Large Hadron Collider, JINST 3, S08003 (2008).  DOI 10.1088/1748-0221/3/08/S08003
  * (3) S. Chatrchyan, et al., The CMS Experiment at the CERN LHC, JINST 3, S08004 (2008).  DOI 10.1088/1748-0221/3/08/S08004
  * (4) G. Apollinari, et al., High-Luminosity Large Hadron Collider (HL-LHC) (2017).  DOI 10.23731/CYRM-2017-004
  * (5) A.A. Alves, Jr, et al., A Roadmap for HEP Software and Computing R&D for the 2020s.  Tech. rep. (2017).  URL https://arxiv.org/abs/1712.06982
  * (6) /QCD_Pt-30to50_TuneZ2_7TeV_pythia6/Summer11LegDRPU_S13_START53_LV6-v1/AODSIM.  DOI 10.7483\/OPENDATA.CMS.Q3BX.69VQ
  * (7) /QCD_Pt-50to80_TuneZ2_7TeV_pythia6/Summer11LegDRPU_S13_START53_LV6-v1/AODSIM.  DOI 10.7483\/OPENDATA.CMS.84VC.RU8W
  * (8) /QCD_Pt-80to120_TuneZ2_7TeV_pythia6/Summer11LegDRPU_S13_START53_LV6-v1/AODSIM.  DOI 10.7483\/OPENDATA.CMS.PUTE.7H2H
  * (9) /QCD_Pt-120to170_TuneZ2_7TeV_pythia6/Summer11LegDRPU_S13_START53_LV6-v1/AODSIM.  DOI 10.7483\/OPENDATA.CMS.QJND.HA88
  * (10) /QCD_Pt-170to300_TuneZ2_7TeV_pythia6/Summer11LegDRPU_S13_START53_LV6-v1/AODSIM.  DOI 10.7483\/OPENDATA.CMS.WKRR.DCJP
  * (11) /QCD_Pt-300to470_TuneZ2_7TeV_pythia6/Summer11LegDRPU_S13_START53_LV6-v1/AODSIM.  DOI 10.7483\/OPENDATA.CMS.X3XQ.USQR
  * (12) /QCD_Pt-470to600_TuneZ2_7TeV_pythia6/Summer11LegDRPU_S13_START53_LV6-v1/AODSIM.  DOI 10.7483\/OPENDATA.CMS.BKTD.SGJX
  * (13) /QCD_Pt-600to800_TuneZ2_7TeV_pythia6/Summer11LegDRPU_S13_START53_LV6-v1/AODSIM.  DOI 10.7483\/OPENDATA.CMS.EJT7.KSAY
  * (14) /QCD_Pt-800to1000_TuneZ2_7TeV_pythia6/Summer11LegDRPU_S13_START53_LV6-v1/AODSIM.  DOI 10.7483\/OPENDATA.CMS.S3D5.KF2C
  * (15) /QCD_Pt-1000to1400_TuneZ2_7TeV_pythia6/Summer11LegDRPU_S13_START53_LV6-v1/AODSIM.  DOI 10.7483\/OPENDATA.CMS.96U2.3YAH
  * (16) /QCD_Pt-1400to1800_TuneZ2_7TeV_pythia6/Summer11LegDRPU_S13_START53_LV6-v1/AODSIM.  DOI 10.7483\/OPENDATA.CMS.RC9V.B5KX
  * (17) /QCD_Pt-1800_TuneZ2_7TeV_pythia6/Summer11LegDRPU_S13_START53_LV6-v1/AODSIM.  DOI 10.7483\/OPENDATA.CMS.CX2X.J3KW
  * (18) R. Mount, et al., Data Preservation in High Energy Physics, Intermediate report of the ICFA-DPHEP Study Group (2009).  URL https://arxiv.org/abs/0912.0255
  * (19) I.J. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xu, D. Warde-Farley, S. Ozair, A. Courville, Y. Bengio, Generative Adversarial Networks (2014).  URL https://arxiv.org/abs/1406.2661
  * (20) S. Arora, R. Ge, Y. Liang, T. Ma, Y. Zhang, in _Proceedings of the 34th International Conference on Machine Learning_ , _Proceedings of Machine Learning Research_ , vol. 70, ed. by D. Precup, Y.W. Teh (PMLR, International Convention Centre, Sydney, Australia, 2017), _Proceedings of Machine Learning Research_ , vol. 70, pp. 224–232.  URL http://proceedings.mlr.press/v70/arora17a.html
  * (21) S. Nowozin, B. Cseke, R. Tomioka, f-GAN: Training Generative Neural Samplers using Variational Divergence Minimization (2016).  URL https://arxiv.org/abs/1606.00709
  * (22) I. Gulrajani, F. Ahmed, M. Arjovsky, V. Dumoulin, A.C. Courville, Improved training of wasserstein gans, CoRR abs/1704.00028 (2017).  URL http://arxiv.org/abs/1704.00028
  * (23) M. Arjovsky, S. Chintala, L. Bottou, Wasserstein GAN (2017).  URL https://arxiv.org/abs/1701.07875
  * (24) G. Louppe, M. Kagan, K. Cranmer, Learning to pivot with adversarial networks (2016).  URL https://arxiv.org/abs/1611.01046
  * (25) C. Shimmin, P. Sadowski, P. Baldi, E. Weik, D. Whiteson, E. Goul, A. Søgaard, Decorrelated Jet Substructure Tagging using Adversarial Neural Networks, Phys. Rev. D96(7), 074034 (2017).  DOI 10.1103/PhysRevD.96.074034
  * (26) V. Estrade, et al., in _NIPS 2017 - workshop Deep Learning for Physical Sciences_ (Long Beach, United States, 2017), pp. 1–5
  * (27) L. de Oliveira, M. Paganini, B. Nachman, Learning Particle Physics by Example: Location-Aware Generative Adversarial Networks for Physics Synthesis, Comput. Softw. Big Sci. 1(1), 4 (2017).  DOI 10.1007/s41781-017-0004-6
  * (28) M. Paganini, L. de Oliveira, B. Nachman, Accelerating Science with Generative Adversarial Networks: An Application to 3D Particle Showers in Multilayer Calorimeters, Phys. Rev. Lett. 120(4), 042003 (2018).  DOI 10.1103/PhysRevLett.120.042003
  * (29) F. Carminati, et al., in _NIPS 2017 - workshop Deep Learning for Physical Sciences_ (Long Beach, United States, 2017), pp. 1–5
  * (30) M. Erdmann, L. Geiger, J. Glombitza, D. Schmidt, Generating and refining particle detector simulations using the Wasserstein distance in adversarial networks (2018).  URL https://arxiv.org/abs/1802.03325
  * (31) P. Isola, J. Zhu, T. Zhou, A.A. Efros, Image-to-image translation with conditional adversarial networks, CoRR abs/1611.07004 (2016).  URL http://arxiv.org/abs/1611.07004
  * (32) T. Salimans, I.J. Goodfellow, W. Zaremba, V. Cheung, A. Radford, X. Chen, Improved techniques for training gans, CoRR abs/1606.03498 (2016).  URL http://arxiv.org/abs/1606.03498
  * (33) M. Mirza, S. Osindero, Conditional generative adversarial nets, CoRR abs/1411.1784 (2014).  URL http://arxiv.org/abs/1411.1784
  * (34) L. de Oliveira, M. Paganini, B. Nachman, in _18th International Workshop on Advanced Computing and Analysis Techniques in Physics Research (ACAT 2017) Seattle, WA, USA, August 21-25, 2017_ (2017)
  * (35) J. de Favereau, C. Delaere, P. Demin, A. Giammanco, V. Lemaître, A. Mertens, M. Selvaggi, DELPHES 3, A modular framework for fast simulation of a generic collider experiment, JHEP 02, 057 (2014).  DOI 10.1007/JHEP02(2014)057
  * (36) T. Sjostrand, S. Mrenna, P.Z. Skands, PYTHIA 6.4 Physics and Manual, JHEP 05, 026 (2006).  DOI 10.1088/1126-6708/2006/05/026
  * (37) M. Cacciari, G.P. Salam, G. Soyez, The anti- k t jet clustering algorithm, Journal of High Energy Physics 2008(04), 063 (2008).  URL http://stacks.iop.org/1126-6708/2008/i=04/a=063
  * (38) M. Cacciari, G.P. Salam, G. Soyez, FastJet User Manual, Eur. Phys. J. C72, 1896 (2012).  DOI 10.1140/epjc/s10052-012-1896-2
  * (39) The data used in this work is available through the following electronic reference.  DOI 10.5281/zenodo.1467678.  URL http://doi.org/10.5281/zenodo.1467678
  * (40) H. Geoffry, et al.  Neural networks for machine learning.  www.coursera.org
  * (41) Y. LeCun, L. Bottou, Y. Bengio, P. Haffner, Gradient-based learning applied to document recognition, Proceedings of the IEEE 86(11), 2278 (1998)
  * (42) O. Ronneberger, P. Fischer, T. Brox, U-net: Convolutional networks for biomedical image segmentation, CoRR abs/1505.04597 (2015).  URL http://arxiv.org/abs/1505.04597
  * (43) S. Ioffe, C. Szegedy, Batch normalization: Accelerating deep network training by reducing internal covariate shift, CoRR abs/1502.03167 (2015).  URL http://arxiv.org/abs/1502.03167
  * (44) A.L. Maas, A.Y. Hannun, A.Y. Ng, in _Proc. icml_ , vol. 30 (2013), vol. 30, p. 3
  * (45) A software package related to this work is available through the following electronic reference.  DOI 10.5281/zenodo.1467665.  URL http://doi.org/10.5281/zenodo.1467665
  * (46) M. Abadi, P. Barham, J. Chen, Z. Chen, A. Davis, J. Dean, M. Devin, S. Ghemawat, G. Irving, M. Isard, M. Kudlur, J. Levenberg, R. Monga, S. Moore, D.G. Murray, B. Steiner, P.A. Tucker, V. Vasudevan, P. Warden, M. Wicke, Y. Yu, X. Zhang, Tensorflow: A system for large-scale machine learning, CoRR abs/1605.08695 (2016).  URL http://arxiv.org/abs/1605.08695
  * (47) F. Chollet, et al.  Keras.  https://keras.io (2015)
  * (48) D.P. Kingma, J. Ba, Adam: A method for stochastic optimization, CoRR abs/1412.6980 (2014).  URL http://arxiv.org/abs/1412.6980
  * (49) F. Pandolfi, Search for the Standard Model Higgs boson in the H $\to$ ZZ $\to\ell^{+}\ell^{-}q\bar{q}$ decay channel at CMS on 4.6 fb-1 of 7 TeV proton-proton collision data, The European Physical Journal Plus 128(10), 117 (2013).  DOI 10.1140/epjp/i2013-13117-x
  * (50) J. Thaler, K. Van Tilburg, Identifying boosted objects with n-subjettiness, Journal of High Energy Physics 2011(3), 15 (2011).  DOI 10.1007/JHEP03(2011)015
  * (51) M. Cacciari, G.P. Salam, G. Soyez, FastJet User Manual, Eur. Phys. J. C72, 1896 (2012).  DOI 10.1140/epjc/s10052-012-1896-2
  * (52) M. Cacciari, G.P. Salam, Dispelling the $N^{3}$ myth for the $k_{t}$ jet-finder, Phys. Lett. B641, 57 (2006).  DOI 10.1016/j.physletb.2006.08.037
  * (53) R. Rahmat, R. Kroeger, A. Giammanco, The fast simulation of the cms experiment, Journal of Physics: Conference Series 396(6), 062016.  DOI 10.1088/1742-6596/396/6/062016

##  Appendix A Comparison with traditional fast simulation methods

In this appendix we compare the results obtained using our generative neural
networks with methods for fast simulation traditionally used for
phenomenological studies in high energy physics. These packages parametrise
particle detectors responses, on a particle by particle basis, in terms of
gaussian energy and momentum resolutions and binomial reconstruction
probabilities. In particular, we compared the results obtained with our GAN
with the output of the DELPHES delphes  package. In order to do that, we used
version 3.4.1 of the program and the standard CMS datacard included in the
package.

In figures 9 and 10 we show the same quantities as shown in figures 5 and 6,
respectively, but with the addition of the results obtained with the DELPHES
gaussian-smearing model (shown in gray). As can be seen the simple gaussian-
smearing approach fails to describe the fine structure of particle jets, and
shows that our neural network based approach provides results of superior
quality, compared to traditional fast simulation methods.

Figure 9:  Aggregated pixel intensities for different rings in
$\Delta\eta$–$\Delta\phi$. Blue histograms are obtained from the input data;
red ones are obtained using the generative model; gray histograms are obtained
with a gaussian-smearing model.

Figure 10:  Evolution of the aggregated pixel intensities for different rings
in $\Delta\eta$–$\Delta\phi$ as a function of the particle level jet
transverse momentum. Solid lines represent the median of the distribution,
filled regions show the inter-quartile range, while dashed lines mark the 10%
and 90% quantiles. Blue lines are obtained from the input data; red ones are
obtained using the generative model; gray ones using the gaussian-smearing
model.

##  Appendix B Differential characterisation of the image translation

In this appendix we report a differential comparison of the jet images
obtained at particle level, reconstruction level, and through our generative
model. In figures 11 and 12 we show the same quantities as shown in figures 5
and 6, respectively, but with the addition of results obtained using particle
level jets.

Figure 11:  Aggregated pixel intensities for different rings in
$\Delta\eta$–$\Delta\phi$. Blue histograms are obtained from the full
simulation and reconstruction chain; red ones are obtained using the
generative model; gray histograms show the quantities obtained before detector
simulation.

Figure 12:  Evolution of the aggregated pixel intensities for different rings
in $\Delta\eta$–$\Delta\phi$ as a function of the particle level jet
transverse momentum. Solid lines represent the median of the distribution,
filled regions show the inter-quartile range, while dashed lines mark the 10%
and 90% quantiles. Blue histograms are obtained from the full simulation and
reconstruction chain; red ones are obtained using the generative model; gray
histograms show the quantities obtained before detector simulation.


# OBSERVATION OF $C\\!P$ VIOLATION IN CHARM DECAYS AT LHCb

F. BETTI on behalf of the LHCb collaboration  Università di Bologna,
Dipartimento di Fisica e Astronomia,
Istituto Nazionale di Fisica Nucleare - Sezione di Bologna,
viale Berti Pichat 6/2, Bologna (40127), Italy

###### Abstract

A search for charge-parity ($C\\!P$) violation in $D^{0}\to K^{-}K^{+}$ and
$D^{0}\to\pi^{-}\pi^{+}$ decays is reported, using $pp$ collision data
corresponding to an integrated luminosity of 5.9 $\mathrm{fb}^{-1}$ collected
at a center-of-mass energy of 13 TeV with the LHCb detector. The flavor of the
$D^{0}$ meson is determined from the charge of the pion in $D^{*}(2010)^{+}\to
D^{0}\pi^{+}$ decays or from the charge of the muon in $\kern
1.79993pt\overline{\kern-1.79993ptB}{}\to D^{0}\mu^{-}\bar{\nu}_{\mu}X$
decays. The difference between the $C\\!P$ asymmetries in $D^{0}\to
K^{-}K^{+}$ and $D^{0}\to\pi^{-}\pi^{+}$ decays is measured to be $\Delta
A_{C\\!P}=[-18.2\pm 3.2\,(\rm stat.)\pm 0.9\,(\rm syst.)]\times 10^{-4}$ for
$\pi$-tagged and $\Delta A_{C\\!P}=[-9\pm 8\,(\rm stat.)\pm 5\,(\rm
syst.)]\times 10^{-4}$ for $\mu$-tagged $D^{0}$ mesons. The combination with
previous LHCb results leads to

$\Delta A_{C\\!P}=(-15.4\pm 2.9)\times 10^{-4},$

where the uncertainty includes both statistical and systematic contributions.
The measured value differs from zero by more than five standard deviations.
This is the first observation of $C\\!P$ violation in the decay of charm
hadrons.

##  1 Introduction

The noninvariance of fundamental interactions under the combined action of
charge conjugation ($C$) and parity ($P$) transformations, so-called $C\\!P$
violation, is a necessary condition for the dynamical generation of the baryon
asymmetry of the universe [1]. $C\\!P$ violation is included in the Standard
Model (SM) of particle physics through an irreducible complex phase in the
Cabibbo-Kobayashi-Maskawa (CKM) quark-mixing matrix [2, 3]. Several
experiments established the presence of $C\\!P$ violation in weak interactions
in the $K$\- and $B$-meson systems [4, 5, 6, 7, 8, 9, 10, 11, 12], and all
results are well interpreted within the CKM formalism. However, the size of
$C\\!P$ violation in the SM is too small to account for the observed matter-
antimatter asymmetry [13, 14, 15], suggesting the existence of beyond-the-SM
sources of $C\\!P$ violation.

Despite decades of experimental searches, the observation of $C\\!P$ violation
in the charm sector has not yet been achieved. Because of the presence of low-
energy strong-interaction effects, theoretical predictions of the size of
$C\\!P$ violation in charm decays are difficult to compute reliably, and the
asymmetries are expected to be of the order of $10^{-4}$–$10^{-3}$ in
magnitude [16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32,
33, 34].

Searches for $C\\!P$ violation in ${{D}^{0}}\\!\to{{K}^{-}}{{K}^{+}}$ and
${{D}^{0}}\\!\to{{\pi}^{-}}{{\pi}^{+}}$ modes111The inclusion of charge-
conjugate decay modes is implied throughout except in asymmetry definitions.
have been performed by the BaBar [35], Belle [36], CDF [37, 38] and LHCb [39,
40, 41, 42, 43] collaborations, which measured values of $C\\!P$ asymmetries
consistent with zero within a precision of a few per mille. This document
presents a measurement of the difference of the time-integrated $C\\!P$
asymmetries in ${{D}^{0}}\\!\to{{K}^{-}}{{K}^{+}}$ and
${{D}^{0}}\\!\to{{\pi}^{-}}{{\pi}^{+}}$ decays, performed using $p$ $p$
collision data collected with the LHCb detector between 2015 and 2018 at a
center-of-mass energy of 13 TeV, corresponding to an integrated luminosity of
5.9$\text{\,fb}^{-1}$.

The time-dependent $C\\!P$ asymmetry, $A_{{C\\!P}}(f;\,t)$, between states
produced as ${D}^{0}$ or ${\kern
1.99997pt\overline{\kern-1.99997ptD}{}}{}^{0}$ mesons decaying to a $C\\!P$
eigenstate $f$ at time $t$ is defined as

$A_{{C\\!P}}(f;\,t)\equiv\frac{\Gamma({{D}^{0}}(t)\to f)-\Gamma({{\kern
1.99997pt\overline{\kern-1.99997ptD}{}}{}^{0}}(t)\to
f)}{\Gamma({{D}^{0}}(t)\to f)+\Gamma({{\kern
1.99997pt\overline{\kern-1.99997ptD}{}}{}^{0}}(t)\to f)},$ (1)

where $\Gamma$ denotes the time-dependent rate of a given decay. For
$f={{K}^{-}}{{K}^{+}}$ or $f={{\pi}^{-}}{{\pi}^{+}}$, $A_{{C\\!P}}(f;\,t)$ can
be expressed in terms of a direct component associated to $C\\!P$ violation in
the decay amplitude and another component associated to $C\\!P$ violation in
${D}^{0}$–${\kern 1.99997pt\overline{\kern-1.99997ptD}{}}{}^{0}$ mixing or in
the interference between mixing and decay. The corresponding time-integrated
asymmetry, $A_{{C\\!P}}(f)$, can be written to first order in the
${D}^{0}$–${\kern 1.99997pt\overline{\kern-1.99997ptD}{}}{}^{0}$ mixing
parameters as [37, 44]

$A_{{C\\!P}}(f)\approx a_{{C\\!P}}^{\rm dir}(f)-\frac{\langle
t(f)\rangle}{\tau({{D}^{0}})}\,A_{\Gamma}(f),$ (2)

where $\langle t(f)\rangle$ denotes the mean decay time of ${{D}^{0}}\to f$
decays in the reconstructed sample, $a_{{C\\!P}}^{\rm dir}(f)$ is the direct
$C\\!P$ asymmetry, $\tau({{D}^{0}})$ the ${D}^{0}$ lifetime and
$A_{\Gamma}(f)$ the asymmetry between the ${{D}^{0}}\to f$ and ${{\kern
1.99997pt\overline{\kern-1.99997ptD}{}}{}^{0}}\to f$ effective decay widths
[45, 46]. Taking $A_{\Gamma}$ to be independent of the final state [19, 47,
48], the difference between $C\\!P$ asymmetries in
${{D}^{0}}\\!\to{{K}^{-}}{{K}^{+}}$ and
${{D}^{0}}\\!\to{{\pi}^{-}}{{\pi}^{+}}$ decays is

$\displaystyle\Delta A_{{C\\!P}}$ $\displaystyle\equiv$ $\displaystyle
A_{{C\\!P}}({{K}^{-}}{{K}^{+}})-A_{{C\\!P}}({{\pi}^{-}}{{\pi}^{+}})$ (3)
$\displaystyle\approx$ $\displaystyle\Delta a_{{C\\!P}}^{\rm
dir}-\frac{\Delta\langle t\rangle}{\tau({{D}^{0}})}\,A_{\Gamma},$

where $\Delta a_{{C\\!P}}^{\rm dir}\equiv a_{{C\\!P}}^{\rm
dir}(K^{-}K^{+})-a_{{C\\!P}}^{\rm dir}(\pi^{-}\pi^{+})$ and $\Delta\langle
t\rangle$ is the difference of the mean decay times $\langle
t({{K}^{-}}{{K}^{+}})\rangle$ and $\langle t({{\pi}^{-}}{{\pi}^{+}})\rangle$.

The ${D}^{0}$ mesons considered in this analysis are produced in two ways:
promptly at a $p$ $p$ collision point (primary vertex, PV) in the strong
${{D}^{*}}(2010)^{+}\\!\to{{D}^{0}}{{\pi}^{+}}$ decay (hereafter
${{D}^{*}}(2010)^{+}$ is referred to as ${{D}^{*+}}$) or at a vertex displaced
from any PV in semileptonic ${{\kern
1.79993pt\overline{\kern-1.79993ptB}{}}}\\!\to{{D}^{0}}\mu^{-}\bar{\nu}_{\mu}X$
decays, where ${{\kern 1.79993pt\overline{\kern-1.79993ptB}{}}}$ denotes a
hadron containing a $b$ quark and $X$ stands for additional particles. The
flavor at production of ${D}^{0}$ mesons from ${{D}^{*+}}$ decays is obtained
from the charge of the accompanying pion ($\pi$-tagged), whereas that of
${D}^{0}$ mesons from semileptonic $b$-hadron decays is determined from the
charge of the accompanying muon ($\mu$-tagged). The raw asymmetries measured
for $\pi$-tagged and $\mu$-tagged ${D}^{0}$ decays are defined as

$\displaystyle\begin{split}A_{\rm
raw}^{\pi{\mbox{-}}\rm{tagged}}(f)&\equiv\frac{N\left({{D}^{*+}}\to{{D}^{0}}(f)\pi^{+}\right)-N\left({{D}^{*-}}\to{{\kern
1.99997pt\overline{\kern-1.99997ptD}{}}{}^{0}}(f)\pi^{-}\right)}{N\left({{D}^{*+}}\to{{D}^{0}}(f)\pi^{+}\right)+N\left({{D}^{*-}}\to{{\kern
1.99997pt\overline{\kern-1.99997ptD}{}}{}^{0}}(f)\pi^{-}\right)},\\\ A_{\rm
raw}^{\mu{\mbox{-}}\rm{tagged}}(f)&\equiv\frac{N({{\kern
1.79993pt\overline{\kern-1.79993ptB}{}}}\to{{D}^{0}}(f)\,\mu^{-}\bar{\nu}_{\mu}X)-N({B}\to{{\kern
1.99997pt\overline{\kern-1.99997ptD}{}}{}^{0}}(f)\,\mu^{+}\nu_{\mu}X)}{N({{\kern
1.79993pt\overline{\kern-1.79993ptB}{}}}\to{{D}^{0}}(f)\,\mu^{-}\bar{\nu}_{\mu}X)+N({B}\to{{\kern
1.99997pt\overline{\kern-1.99997ptD}{}}{}^{0}}(f)\,\mu^{+}\nu_{\mu}X)},\end{split}$
(4)

where $N$ is the measured signal yield for each given decay. These can be
approximated as

$\displaystyle\begin{split}A_{\rm raw}^{\pi{\mbox{-}}\rm{tagged}}(f)&\approx
A_{{C\\!P}}(f)+A_{\rm D}(\pi)+A_{\rm P}({{D}^{*}}),\\\ A_{\rm
raw}^{\mu{\mbox{-}}\rm{tagged}}(f)&\approx A_{{C\\!P}}(f)+A_{\rm
D}(\mu)+A_{\rm P}({B}),\end{split}$ (5)

where $A_{\rm D}(\pi)$ and $A_{\rm D}(\mu)$ are detection asymmetries due to
different reconstruction efficiencies between positive and negative tagging
particles, whereas $A_{\rm P}({{D}^{*}})$ and $A_{\rm P}({B})$ are the
production asymmetries of ${D}^{*}$ mesons and $b$ hadrons, arising from the
hadronization of charm and beauty quarks in $p$ $p$ collisions. The involved
terms, averaged over phase space for selected events are
${\mathcal{O}}(10^{-2})$ or less [49, 50, 51, 52], hence the approximations in
Eqs. 5 are valid up to corrections of ${\mathcal{O}}(10^{-6})$. The values of
the detection and production asymmetries are independent of the final state
$f$, and thus cancel in the difference, if the kinematic distributions of the
two channels are equal, resulting in

$\Delta A_{{C\\!P}}=A_{\rm raw}({{K}^{-}}{{K}^{+}})-A_{\rm
raw}({{\pi}^{-}}{{\pi}^{+}}).$ (6)

The relation between $\Delta A_{{C\\!P}}$ and the measurable raw asymmetries
in ${{K}^{-}}{{K}^{+}}$ and ${{\pi}^{-}}{{\pi}^{+}}$ makes the determination
of $\Delta A_{{C\\!P}}$ largely insensitive to systematic uncertainties.

##  2 Selection

The LHCb detector is a single-arm forward spectrometer designed for the study
of particles containing $b$ or $c$ quarks [53, 54]. The magnetic-field
polarity of the dipole magnet used by the LHCb tracking system is reversed
periodically during data taking to mitigate the differences of reconstruction
efficiencies of particles with opposite charges, though the analysis presented
in this document is expected to be insensitive to such effects.

The online event selection is performed by a trigger, which consists of a
hardware stage based on information from the calorimeter and muon systems,
followed by two software stages. ${D}^{0}$ candidates are fully reconstructed
in the second software stage using kinematic, topological and particle-
identification (PID) criteria. In the $\mu$-tagged sample, ${D}^{0}$
candidates are combined with muons to form $B$ candidates, under the
requirement that they are consistent with originating from a common vertex. In
addition, requirements on the invariant mass of the ${{D}^{0}}\mu$ system,
$m({{D}^{0}}\mu)$, and on the corrected mass222The corrected mass is defined
as $m_{\rm
corr}\equiv\sqrt{m({{D}^{0}}\mu)^{2}+p_{\perp}({{D}^{0}}\mu)^{2}}+p_{\perp}({{D}^{0}}\mu)$
[55], where $p_{\perp}({{D}^{0}}\mu)$ is the momentum of the ${{D}^{0}}\mu$
system transverse to the flight direction of the $b$ hadron. are applied in
the $\mu$-tagged sample.

In certain kinematic regions very large raw asymmetries, up to 100%, occur
because, for a given magnet polarity, low-momentum particles of one charge at
small or large polar angles in the horizontal plane may be deflected out of
the detector or into the LHC beam pipe, whereas particles with the other
charge are more likely to remain within the acceptance. For this reason, in
the offline selection, fiducial requirements are imposed to exclude kinematic
regions characterized by large detection asymmetries for the tagging particle.
About 35% and 10% of the selected candidates are rejected by these fiducial
requirements for the $\pi$-tagged and $\mu$-tagged samples, respectively. For
$\pi$-tagged ${D}^{0}$ mesons, a requirement on the ${D}^{0}$
$\chi^{2}_{\text{IP}}$ is applied333The $\chi^{2}_{\text{IP}}$ is defined as
the difference between the $\chi^{2}$ of the PV reconstructed with and without
the considered particle. to suppress the background of ${D}^{0}$ mesons
produced in $B$ decays, and PID requirements on the ${D}^{0}$ decay products
are tightened. The ${D}^{0}$ and pion candidates are combined to form
${D}^{*+}$ candidates by requiring a good fit quality of the ${D}^{*+}$
vertex, that is constrained to coincide with the nearest PV [56]. The
invariant mass of ${D}^{0}$ candidates is required to lie within a range of
about $\pm 3$ standard deviations around the known ${{D}^{0}}$ mass. For
$\mu$-tagged mesons, in order to suppress the combinatorial background due to
random combinations of charged kaon or pion pairs not originating from a
${D}^{0}$ decay, the $B$ candidates are further filtered using a dedicated
boosted decision tree (BDT) that uses variables related to the topology and
the kinematics of the reconstructed decay. A veto in the invariant mass of the
$\mu^{\mp}\pi^{\pm}$ ($\mu^{\mp}K^{\pm}$) pair, where the pion (kaon) is given
the muon mass hypothesis, is applied to suppress background from $b$-hadron
decays to ${c}{\overline{c}}\pi^{\pm}X$ (${c}{\overline{c}}K^{\pm}X$), where
the ${c}{\overline{c}}$ resonance decays to a pair of muons.

The data sample includes events with multiple ${D}^{*+}$ and $B$ candidates,
that are mostly due to a common reconstructed ${D}^{0}$ meson combined with
different tagging particles. The fractions of events with multiple candidates
are about 10% and 0.4% in the $\pi$-tagged and $\mu$-tagged samples,
respectively. When multiple candidates are present in the event, only one is
kept randomly.

Since the detection and production asymmetries are expected to depend on the
kinematics of the reconstructed particles, the possible difference between the
kinematic distributions of reconstructed ${D}^{*+}$ or $B$ candidates and of
the tagging pions or muons in the $K^{-}K^{+}$ and $\pi^{-}\pi^{+}$ decay
modes may induce an incomplete cancellation in the difference in Eq. 6. Hence,
a small correction to the $K^{-}K^{+}$ sample is applied by means of a
weighting procedure: for the $\pi$-tagged sample, the ratio between the three-
dimensional background-subtracted distributions of pseudorapidity, transverse
momentum and azimuthal angle of the ${D}^{*+}$ meson in the $K^{-}K^{+}$ and
$\pi^{-}\pi^{+}$ modes is taken and candidate-by-candidate weights are
calculated. An analogous procedure is followed for the $\mu$-tagged sample,
where ${D}^{0}$ distributions are used in place of those of the ${D}^{*+}$
meson. It is then checked _a posteriori_ that the distributions of the same
variables for tagging pions and muons are also equalized by the weighting. The
application of the weights results in a small variation of $\Delta
A_{{C\\!P}}$, below $10^{-4}$ for both the $\pi$-tagged and $\mu$-tagged
samples.

Figure 1: Mass distributions of selected (top) $\pi^{\pm}$-tagged and (bottom)
$\mu^{\pm}$-tagged candidates for (left) ${{K}^{-}}{{K}^{+}}$ and (right)
${{\pi}^{-}}{{\pi}^{+}}$ final states of the ${D}^{0}$-meson decays, with fit
projections overlaid.

##  3 Measurement of the Asymmetries

For each decay mode, simultaneous least-square fits to the binned mass
distributions of ${D}^{*+}$ and ${D}^{*-}$ candidates for the $\pi$-tagged
sample, or ${D}^{0}$ and ${\kern
1.99997pt\overline{\kern-1.99997ptD}{}}{}^{0}$ candidates for the $\mu$-tagged
sample, are performed to obtain the raw asymmetries of signal and background
components, which are free parameters of the fits.

In the analysis of the $\pi$-tagged sample the fits are performed to the
$m({{D}^{0}}{{\pi}^{+}})$ and $m({{\kern
1.99997pt\overline{\kern-1.99997ptD}{}}{}^{0}}{{\pi}^{-}})$ distributions,
that are defined using the known value of the ${D}^{0}$ mass [37]. The signal
mass model consists of the sum of three Gaussian functions and a Johnson
$S_{{U}}$ function [57], whereas the combinatorial background is described by
an empirical function of the form
$[m({{D}^{0}}{{\pi}^{+}})-m({{D}^{0}})-m({{\pi}^{+}})]^{\alpha}e^{\beta\,m({{D}^{0}}{{\pi}^{+}})}$.
All the parameters of the models are free to be adjusted by the fit and are
shared among positive and negative tags, except for the mean values of the
Gaussian functions, which are different to take into account small shifts in
the raw mass measurements between opposite tags.

In the analysis of the $\mu$-tagged sample, the fits are performed to the
$m({{D}^{0}})$ distributions. The signal is described by the sum of two
Gaussian functions convolved with a truncated power-law function accounting
for final-state photon radiation effects, while the combinatorial background
is modeled by an exponential function. A small contribution from
${{D}^{0}}\\!\to{{K}^{-}}{{\pi}^{+}}$ decays with a misidentified kaon or pion
is visible and is modeled as the tail of a Gaussian function. The fit
parameters are shared among positive and negative tags, except for the mean
values of the Gaussian functions.

Fits are performed to subsamples of data split according to magnet polarities
and years of data taking. The final results are obtained by averaging the
partial $\Delta A_{{C\\!P}}$ values corresponding to each subsample, which are
found to be in good agreement. Performing single fits the overall $\pi$-tagged
and $\mu$-tagged samples gives small differences of the order of a few
$10^{-5}$. Figure 1 displays the $m({{D}^{0}}{{\pi}^{+}})$ and $m({{D}^{0}})$
distributions corresponding to the entire samples. The $\pi$-tagged
($\mu$-tagged) signal yields are approximately $44$ ($9$) million
${{D}^{0}}\\!\to{{K}^{-}}{{K}^{+}}$ decays and $14$ ($3$) million
${D}^{0}\rightarrow\pi^{-}\pi^{+}$ decays.

##  4 Systematic Uncertainties

Several sources of systematic uncertainties affecting the measurement are
considered and studied independently for the $\pi$-tagged and $\mu$-tagged
samples. In the case of $\pi$-tagged decays, the dominant systematic
uncertainty is related to the knowledge of the signal and background mass
models. It is evaluated by generating pseudoexperiments according to the
baseline fit model, then fitting both baseline and alternative models to those
data and considering the difference between the resulting values of $\Delta
A_{{C\\!P}}$. A value of $0.6\times 10^{-4}$, corresponding to the largest
observed variation, is assigned as a systematic uncertainty. A similar study
with pseudoexperiments is also performed with the $\mu$-tagged sample and a
value of $2\times 10^{-4}$ is found.

In the case of $\mu$-tagged decays, the main systematic uncertainty is due to
the possibility that the ${D}^{0}$ flavor is not tagged correctly by the muon
charge because of misreconstruction. The probability of wrongly assigning the
${D}^{0}$ flavor (mistag) is measured on a large sample of $\mu$-tagged
${{D}^{0}}\\!\to{{K}^{-}}{{\pi}^{+}}$ decays by comparing the charges of kaon
and muon candidates. Mistag rates are found to be at the percent level and
compatible for positively and negatively tagged decays, and the corresponding
systematic uncertainty is estimated to be $4\times 10^{-4}$.

Systematic uncertainties of $0.2\times 10^{-4}$ and $1\times 10^{-4}$
accounting for the knowledge of the weights used in the kinematic weighting
procedure are assessed for $\pi$-tagged and $\mu$-tagged decays, respectively.
A fraction of ${D}^{0}$ mesons from $B$ decays (secondary decays) is still
present in the final $\pi$-tagged sample even after the requirement that the
${D}^{0}$ trajectory points back to the PV. Possible different levels of
contamination from secondary decays in ${{D}^{0}}\\!\to{{K}^{-}}{{K}^{+}}$ and
${{D}^{0}}\\!\to{{\pi}^{-}}{{\pi}^{+}}$ samples may bias the value of $\Delta
A_{{C\\!P}}$ because of an incomplete cancellation of the production
asymmetries of $b$ hadrons. The fractions of secondary decays are estimated by
performing a fit to the distribution of the ${D}^{0}$-candidate impact
parameter in the plane transverse to the beam direction, and the corresponding
systematic uncertainty is estimated to be $0.3\times 10^{-4}$. A systematic
uncertainty associated to the presence of background components peaking in
$m({{D}^{0}}\pi)$ and not in $m({{D}^{0}})$ is determined by fits to the
$m({{D}^{0}})$ distributions after the removal of the signal window
requirement, where these components are modeled using fast simulation. The
main sources are the ${{D}^{0}}\\!\to{{K}^{-}}{{\pi}^{+}}{{\pi}^{0}}$ decay
for the ${{K}^{+}}{{K}^{-}}$ mode, and the
${{D}^{0}}\to{{\pi}^{-}}{\mu^{+}}{{\nu}_{\mu}}$ and
${{D}^{0}}\to{{\pi}^{-}}{e^{+}}{{\nu}_{e}}$ decays for the
${{\pi}^{+}}{{\pi}^{-}}$ mode. Yields and raw asymmetries of the peaking-
background components measured from the fits are then used as inputs to
pseudoexperiments performed to evaluate the corresponding effects on the
determination of $\Delta A_{{C\\!P}}$, resulting in a systematic uncertainty
of $0.5\times 10^{-4}$.

In the case of $\mu$-tagged decays, the fractions of reconstructed ${\kern
1.79993pt\overline{\kern-1.79993ptB}{}}$ decays can be slightly different
between the $K^{-}K^{+}$ and $\pi^{-}\pi^{+}$ decay modes, which could lead to
a small bias in $\Delta A_{{C\\!P}}$. Using the LHCb measurements of the
$b$-hadron production asymmetries [49], the associated systematic uncertainty
on $\Delta A_{{C\\!P}}$ is estimated to be $1\times 10^{-4}$. The combination
of a difference in the $B$ reconstruction efficiency as a function of the
decay time between the ${{D}^{0}}\\!\to{{K}^{-}}{{K}^{+}}$ and
${{D}^{0}}\\!\to{{\pi}^{-}}{{\pi}^{+}}$ modes and the presence of neutral
$B$-meson oscillations may also cause an imperfect cancellation of $A_{\rm
P}(B)$ in $\Delta A_{{C\\!P}}$, and the related systematic uncertainty is
estimated to be $2\times 10^{-4}$.

The total systematic uncertainties on $\Delta A_{{C\\!P}}$ are given by the
sum in quadrature of all individual contributions, and are equal to $0.9\times
10^{-4}$ and $5\times 10^{-4}$ for the $\pi$-tagged and $\mu$-tagged samples,
respectively. A summary of all systematic uncertainties is reported in Table
1.

Table 1: Systematic uncertainties on $\Delta A_{{C\\!P}}$ for $\pi$\- and $\mu$-tagged decays (in $10^{-4}$). The total uncertainties are obtained as the sums in quadrature of the individual contributions. Source |  $\pi$-tagged |  $\mu$-tagged
---|---|---
Fit model | 0.6 | 2
Mistag | – | 4
Weighting | 0.2 | 1
Secondary decays | 0.3 | –
Peaking background | 0.5 | –
$B$ fractions | – | 1
$B$ reco. efficiency | – | 2
Total | 0.9 | 5

Numerous additional robustness checks are carried out. The measured value of
$\Delta A_{{C\\!P}}$ is studied as a function of several geometrical and
kinematic variables. Furthermore, the total sample is split into subsamples
taken in different run periods within the years of data taking, also
distinguishing different magnet polarities. No evidence for unexpected
dependences of $\Delta A_{{C\\!P}}$ is found in any of these tests. A check
using more stringent PID requirements is performed, and all variations of
$\Delta A_{{C\\!P}}$ are found to be compatible within statistical
uncertainties. An additional check concerns the measurement of $\Delta A_{\rm
bkg}$, which is the difference of the background raw asymmetries in
${{K}^{-}}{{K}^{+}}$ and ${{\pi}^{-}}{{\pi}^{+}}$ final states. The prompt
background is mainly composed of genuine ${D}^{0}$ candidates paired with
unrelated pions originating from the PV, so $\Delta A_{\rm bkg}$ is expected
to be compatible with zero. A value of $\Delta A_{\rm bkg}=(-2\pm 4)\times
10^{-4}$ is obtained.

##  5 Results

The measured differences of time-integrated $C\\!P$ asymmetries of
${{D}^{0}}\\!\to{{K}^{-}}{{K}^{+}}$ and ${{D}^{0}}\to{{\pi}^{-}}{{\pi}^{+}}$
decays are [58]

$\displaystyle\Delta A_{{C\\!P}}^{\pi{\mbox{-}}\rm{tagged}}$
$\displaystyle=\left[-18.2\pm 3.2\,\rm{(stat.)}\pm
0.9\,\rm{(syst.)}\right]\times 10^{-4},$ $\displaystyle\Delta
A_{{C\\!P}}^{\mu{\mbox{-}}\rm{tagged}}$ $\displaystyle=\left[-9\pm
8\,\rm{(stat.)}\pm 5\,\rm{(syst.)}\right]\times 10^{-4},$

both in good agreement with world averages [59] and previous LHCb results [42,
41].

The full combination with previous LHCb measurements [42, 41] gives the
following value of $\Delta A_{{C\\!P}}$

$\Delta A_{{C\\!P}}=\left(-15.4\pm 2.9\right)\times 10^{-4},$

where the uncertainty includes statistical and systematic contributions. The
significance of the deviation from zero corresponds to 5.3 standard
deviations. This is the first observation of $C\\!P$ violation in the decay of
charm hadrons.

As shown in Eq. 3, the interpretation of $\Delta A_{{C\\!P}}$ in terms of
direct $C\\!P$ violation and $A_{\Gamma}$ requires knowledge of the difference
of reconstructed mean decay times for ${{D}^{0}}\\!\to{{K}^{-}}{{K}^{+}}$ and
${{D}^{0}}\\!\to{{\pi}^{-}}{{\pi}^{+}}$ decays normalized to the ${D}^{0}$
lifetime. The values corresponding to the present measurements, using the
world average of the ${D}^{0}$ lifetime [60], are $\Delta\\!\left\langle
t\right\rangle^{\pi{\mbox{-}}\rm{tagged}}/\tau({{D}^{0}})=0.135\pm 0.002$ and
$\Delta\\!\left\langle
t\right\rangle^{\mu{\mbox{-}}\rm{tagged}}/\tau({{D}^{0}})=-0.003\pm 0.001$,
whereas that corresponding to the full combination is $\Delta\left\langle
t\right\rangle/\tau({{D}^{0}})=0.115\pm 0.002$. The uncertainties include
statistical and systematic contributions. By using the LHCb average [45, 46]
$A_{\Gamma}=(-2.8\pm 2.8)\times 10^{-4}$, from Eq. 3 it is possible to derive

$\Delta a_{{C\\!P}}^{\rm dir}=\left(-15.7\pm 2.9\right)\times 10^{-4}.$

As expected, $\Delta A_{{C\\!P}}$ is primarily sensitive to direct $C\\!P$
violation.

In summary, this document reports the first observation of a nonzero $C\\!P$
asymmetry in charm decays, using large samples of
${{D}^{0}}\\!\to{{K}^{-}}{{K}^{+}}$ and
${{D}^{0}}\\!\to{{\pi}^{-}}{{\pi}^{+}}$ decays collected with the LHCb
detector. The result is consistent with, although in magnitude at the upper
end of, SM expectations. In the next future, further measurements with charmed
particles, along with possible theoretical improvements, will help clarify the
present physics picture, to establish whether this result is consistent with
the SM or indicates the presence of new physics processes in the up-quark
sector.

## References

## References

  * [1] A. D. Sakharov, Pisma Zh. Eksp. Teor. Fiz. 5 (1967) 32
  * [2] N. Cabibbo, Phys. Rev. Lett. 10 (1963) 531
  * [3] M. Kobayashi and T. Maskawa, Prog. Theor. Phys. 49 (1973) 652
  * [4] J. H. Christenson, J. W. Cronin, V. L. Fitch, and R. Turlay, Phys. Rev. Lett. 13 (1964) 138
  * [5] KTeV collaboration, A. Alavi-Harati et al., Phys. Rev. Lett. 83 (1999) 22, arXiv:hep-ex/9905060
  * [6] NA48 collaboration, A. Lai et al., Eur. Phys. J. C22 (2001) 231, arXiv:hep-ex/0110019
  * [7] BaBar collaboration, B. Aubert et al., Phys. Rev. Lett. 87 (2001) 091801, arXiv:hep-ex/0107013
  * [8] Belle collaboration, K. Abe et al., Phys. Rev. Lett. 87 (2001) 091802, arXiv:hep-ex/0107061
  * [9] BaBar collaboration, B. Aubert et al., Phys. Rev. Lett. 93 (2004) 131801, arXiv:hep-ex/0407057
  * [10] Belle collaboration, Y. Chao et al., Phys. Rev. Lett. 93 (2004) 191802, arXiv:hep-ex/0408100
  * [11] LHCb collaboration, R. Aaij et al., Phys. Rev. Lett. 110 (2013) 221601, arXiv:1304.6173
  * [12] LHCb collaboration, R. Aaij et al., Phys. Lett. B712 (2012) 203, Erratum ibid. B713 (2012) 351, arXiv:1203.3662
  * [13] A. G. Cohen, D. B. Kaplan, and A. E. Nelson, Ann. Rev. Nucl. Part. Sci. 43 (1993) 27, arXiv:hep-ph/9302210
  * [14] A. Riotto and M. Trodden, Ann. Rev. Nucl. Part. Sci. 49 (1999) 35, arXiv:hep-ph/9901362
  * [15] W.-S. Hou, Chin. J. Phys. 47 (2009) 134, arXiv:0803.1234
  * [16] M. Golden and B. Grinstein, Phys. Lett. B222 (1989) 501
  * [17] F. Buccella et al., Phys. Rev. D51 (1995) 3478, arXiv:hep-ph/9411286
  * [18] S. Bianco, F. L. Fabbri, D. Benson, and I. Bigi, Riv. Nuovo Cim. 26N7 (2003) 1, arXiv:hep-ex/0309021
  * [19] Y. Grossman, A. L. Kagan, and Y. Nir, Phys. Rev. D75 (2007) 036008, arXiv:hep-ph/0609178
  * [20] M. Artuso, B. Meadows, and A. A. Petrov, Ann. Rev. Nucl. Part. Sci. 58 (2008) 249, arXiv:0802.2934
  * [21] J. Brod, A. L. Kagan, and J. Zupan, Phys. Rev. D86 (2012) 014023, arXiv:1111.5000
  * [22] H.-Y. Cheng and C.-W. Chiang, Phys. Rev. D85 (2012) 034036, Erratum ibid. D85 (2012) 079903, arXiv:1201.0785
  * [23] H.-Y. Cheng and C.-W. Chiang, Phys. Rev. D86 (2012) 014014, arXiv:1205.0580
  * [24] H.-n. Li, C.-D. Lu, and F.-S. Yu, Phys. Rev. D86 (2012) 036012, arXiv:1203.3120
  * [25] E. Franco, S. Mishima, and L. Silvestrini, JHEP 05 (2012) 140, arXiv:1203.3131
  * [26] D. Pirtskhalava and P. Uttayarat, Phys. Lett. B712 (2012) 81, arXiv:1112.5451
  * [27] T. Feldmann, S. Nandi, and A. Soni, JHEP 06 (2012) 007, arXiv:1202.3795
  * [28] J. Brod, Y. Grossman, A. L. Kagan, and J. Zupan, JHEP 10 (2012) 161, arXiv:1203.6659
  * [29] G. Hiller, M. Jung, and S. Schacht, Phys. Rev. D87 (2013) 014024, arXiv:1211.3734
  * [30] Y. Grossman and D. J. Robinson, JHEP 04 (2013) 067, arXiv:1211.3361
  * [31] B. Bhattacharya, M. Gronau, and J. L. Rosner, Phys. Rev. D85 (2012) 054014, arXiv:1201.2351
  * [32] S. Müller, U. Nierste, and S. Schacht, Phys. Rev. Lett. 115 (2015) 251802, arXiv:1506.04121
  * [33] A. Khodjamirian and A. A. Petrov, Phys. Lett. B774 (2017) 235, arXiv:1706.07780
  * [34] F. Buccella, A. Paul, and P. Santorelli, arXiv:1902.05564
  * [35] BaBar collaboration, B. Aubert et al., Phys. Rev. Lett. 100 (2008) 061803, arXiv:0709.2715
  * [36] Belle collaboration, M. Starič et al., Phys. Lett. B670 (2008) 190, arXiv:0807.0148
  * [37] CDF collaboration, T. Aaltonen et al., Phys. Rev. D85 (2012) 012009, arXiv:1111.5023
  * [38] CDF collaboration, T. Aaltonen et al., Phys. Rev. Lett. 109 (2012) 111801, arXiv:1207.2158
  * [39] LHCb collaboration, R. Aaij et al., Phys. Rev. Lett. 108 (2012) 111602, arXiv:1112.0938
  * [40] LHCb collaboration, R. Aaij et al., Phys. Lett. B723 (2013) 33, arXiv:1303.2614
  * [41] LHCb collaboration, R. Aaij et al., JHEP 07 (2014) 041, arXiv:1405.2797
  * [42] LHCb collaboration, R. Aaij et al., Phys. Rev. Lett. 116 (2016) 191601, arXiv:1602.03160
  * [43] LHCb collaboration, R. Aaij et al., Phys. Lett. B767 (2017) 177, arXiv:1610.09476
  * [44] M. Gersabeck et al., J. Phys. G39 (2012) 045005, arXiv:1111.6515
  * [45] LHCb collaboration, R. Aaij et al., JHEP 04 (2015) 043, arXiv:1501.06777
  * [46] LHCb collaboration, R. Aaij et al., Phys. Rev. Lett. 118 (2017) 261803, arXiv:1702.06490
  * [47] A. L. Kagan and M. D. Sokoloff, Phys. Rev. D80 (2009) 076008, arXiv:0907.3917
  * [48] D.-S. Du, Eur. Phys. J. C50 (2007) 579, arXiv:hep-ph/0608313
  * [49] LHCb collaboration, R. Aaij et al., Phys. Lett. B774 (2017) 139, arXiv:1703.08464
  * [50] LHCb collaboration, R. Aaij et al., Phys. Lett. B728 (2014) 607, arXiv:1308.1048
  * [51] LHCb collaboration, R. Aaij et al., Phys. Lett. B718 (2013) 902, arXiv:1210.4112
  * [52] LHCb collaboration, R. Aaij et al., Phys. Lett. B713 (2012) 186, arXiv:1205.0897
  * [53] LHCb collaboration, A. A. Alves Jr. et al., JINST 3 (2008) S08005
  * [54] LHCb collaboration, R. Aaij et al., Int. J. Mod. Phys. A30 (2015) 1530022, arXiv:1412.6352
  * [55] E653 collaboration, K. Kodama et al., Phys. Rev. Lett. 66 (1991) 1819
  * [56] W. D. Hulsbergen, Nucl. Instrum. Meth. A552 (2005) 566, arXiv:physics/0503191
  * [57] N. L. Johnson, Biometrika 36 (1949) 149
  * [58] LHCb collaboration, R. Aaij et al., arXiv:1903.08726, accepted by Phys. Rev. Lett.
  * [59] Heavy Flavor Averaging Group, Y. Amhis et al., Eur. Phys. J. C77 (2017) 895, arXiv:1612.07233, updated results and plots available at https://hflav.web.cern.ch
  * [60] Particle Data Group, M. Tanabashi et al., Phys. Rev. D98 (2018) 030001