Evidence for the decay $B^0\to p\bar{p}\pi^0$

TL;DR

This paper provides the first evidence for the rare decay $B^0 o par{p}\pi^0$ with a measured branching fraction, using data from the Belle experiment, and sets limits on related intermediate decay modes.

Contribution

First evidence for the charmless baryonic decay $B^0 o par{p}\pi^0$ with a measured branching fraction and upper limits on related two-body decays.

Findings

Branching fraction $ imes 10^{-7}$: $5.0\pm1.8 ext{(stat)}\pm0.6 ext{(syst)}$

Significance of 3.1 standard deviations

Upper limit on $B^0 o \Delta^+ar{p}$ and $B^0 o ar{\Delta}^-p$ decays: $1.6 imes10^{-6}$

Abstract

We report a search for the charmless baryonic decay $B^0\to p\bar{p}\pi^0$ with a data sample corresponding to an integrated luminosity of 711~$\rm fb^{-1}$ containing $(772\pm 10)\times 10^6$ $B^0\bar{B}^0$ pairs. The data was collected by the Belle experiment running on the $\Upsilon(4S)$ resonance at the KEKB $e^+e^-$ collider. We measure a branching fraction $\mathcal{B}(B^0\to p\bar{p}\pi^0)= (5.0\pm1.8\pm0.6 )\times 10^{-7}$, where the first uncertainty is statistical and the second is systematic. The signal has a significance of 3.1 standard deviations and constitutes the first evidence for this decay mode. We also search for the intermediate two-body decays $B^{0}\to\Delta^+\bar{p}$ and $B^0\to\bar{\Delta}^-p$, and set an upper limit on the branching fraction: $\mathcal{B}(B^0\to \Delta^+\bar{p})+\mathcal{B}(B^0\to\bar{\Delta}^-p)<1.6\times10^{-6}$ at 90% confidence level.

Tables

Table 1. Table 1: Systematic uncertainties in ℬ ​ ( B 0 → p ​ p ¯ ​ π 0 ) ℬ → superscript 𝐵 0 𝑝 ¯ 𝑝 superscript 𝜋 0 \mathcal{B}(B^{0}\rightarrow p\bar{p}\pi^{0}) . Those listed in the upper section are associated with fitting for the signal yields and are included in the signal significance.

Source	Uncertainty (%)
PDF parametrization	${}_{-3.2}{}^{}{}_{}{}^{+2.9}$
Calibration factor	$11.9$
Fit bias	${}_{-0.0}{}^{}{}_{}{}^{+2.1}$
${\pi }^0$ reconstruction	1.5
Tracking	0.7
Particle identification	$0.6$
Choice of $C_{\mathrm{NN}}$	${}_{-1.1}{}^{}{}_{}{}^{+2.0}$
Incorrectly reconstructed signal events	${}_{-0.8}{}^{}{}_{}{}^{+1.0}$
Number of $B\overline{B}$ pairs	1.4
MC statistics	0.4
Total	${}_{-12.6}{}^{}{}_{}{}^{+12.8}$

Equations

$$C_{\mathrm{NN}}^{\prime }=\ln \left(\frac{C_{\mathrm{NN}}-C_{\mathrm{NN}}^{\min }}{C_{\mathrm{NN}}^{\max }-C_{\mathrm{NN}}}\right),$$

$$\mathcal{L}=e^{-{\sum }_j{Y}_j}\prod _i^N\left(\sum _j{Y}_j{\mathcal{P}}_j(M_{\mathrm{bc}}^i,\Delta E^i,C_{\mathrm{NN}}^{\prime i})\right),$$

$$\mathcal{B}(B^0\to p\bar{p}{\pi }^0)=\frac{Y_{\mathrm{sig}}}{N_{B\overline{B}}\times \epsilon },$$

$$\mathcal{B}(B^0\to p\bar{p}{\pi }^0)=(5.0\pm 1.8\pm 0.6)\times 10^{-7},$$

$$\mathcal{B}(B^0\to {\Delta }^{+}\bar{p})+\mathcal{B}(B^0\to {\bar{\Delta }}^{-}p)<1.6\times 10^{-6}.$$

Full text

The first observed charmless baryonic $B$ decay was $B^{+}\rightarrow p\bar{p}K^{+}$ Abe:2002ds . Following this first observation, many other charmless baryonic $B$ decays have been found PDG . Except for $B^{+}\rightarrow p\mathchar 28931\relax\pi^{0}~{}{\rm and}~{}p\mathchar 28931\relax\gamma$ decays, all the channels reported to date are entirely reconstructed from charged particles in the final state. A noticeable hierarchy is also observed in the branching fractions of these decays: three-body decays are usually more frequent than their two-body counterparts but less frequent than four-body decays Bevan:2014iga ; Cheng:2014qxa . This phenomenon can be understood in terms of the so-called “threshold effect,” which refers to the fact that the $B$ meson prefers to decay into a di-baryon pair with low invariant mass accompanied by a fast recoil meson Bevan:2014iga ; Hou:2000bz ; Chen:2008sw . This peaking behavior was unexpected, and has led to various speculations about possible mechanisms Bevan:2014iga . Studying additional three-body baryonic decays might provide a better understanding of the dynamics of $B$ decays and the aforementioned threshold-effect. These decays are also useful for $C\!P$ violation studies.

This paper reports a search for a three-body charmless baryonic $B^{0}$ decays to the $p\bar{p}\pi^{0}$ final state charge-conjugate using a data set corresponding to an integrated luminosity of 711 $\rm fb^{-1}$ collected with the Belle detector Belle at the $\mathchar 28935\relax{(4S)}$ resonance at the KEKB asymmetric-energy $e^{+}e^{-}$ (3.5 on 8.0 GeV) collider KEKB . So far, the decay $B^{0}\rightarrow p\bar{p}\pi^{0}$ has not been studied by any experiment. No theoretical prediction for the branching fraction of this process is yet available. A glance at the known branching fractions for $B$ decays PDG shows the three-body charmless baryonic decays to occur in the several times $10^{-6}$ range, indicating that the discovery of the mode $B^{0}\rightarrow p\bar{p}\pi^{0}$ might be possible with the currently available data set.

The Belle detector is a large-solid-angle magnetic spectrometer consisting of a silicon vertex detector (SVD), a 50-layer central drift chamber (CDC), an array of aerogel threshold Cherenkov counters (ACC), a barrel-like arrangement of time-of-flight scintillation counters (TOF), and an electromagnetic calorimeter (ECL) comprising CsI(Tl) crystals. These detector components are located inside a superconducting solenoid coil that provides a 1.5 T magnetic field. An iron flux-return located outside the coil is instrumented (KLM) to detect $K_{L}^{0}$ mesons and to identify muons. Two inner detector configurations were used: a 2.0 cm radius beampipe and a three-layer SVD were used for the first $152\times 10^{6}$ $B\kern 1.79993pt\overline{\kern-1.79993ptB}{}$ pairs of data, while a 1.5 cm radius beampipe, a four-layer SVD, and a small-cell inner drift chamber were used for the remaining $620\times 10^{6}$ $B\kern 1.79993pt\overline{\kern-1.79993ptB}{}$ pairs of data. The detector is described in detail in Ref. Belle . Event selection requirements are optimized using Monte Carlo (MC) simulations. MC events are generated using EvtGen Lange:2001uf , and the detector response is modeled using Geant3 geant3 . Final-state radiation is taken into account using the Photos package Golonka:2005pn .

The reconstruction of $B^{0}\rightarrow p\bar{p}\pi^{0}$ proceeds by first reconstructing $\pi^{0}\rightarrow\gamma\gamma$ candidates. An ECL cluster not matched to any track in the CDC is identified as a photon candidate. Such candidates are required to have an energy greater than 50 MeV in the barrel region and greater than 100 MeV in the end-cap regions, where the barrel region covers the polar angle $32^{\circ}<\theta<130^{\circ}$ and the end-cap regions cover the ranges $12^{\circ}<\theta<32^{\circ}$ and $130^{\circ}<\theta<157^{\circ}$. To reject showers produced by neutral hadrons, the energy deposited in the $3\times 3$ array of ECL crystals centered on the crystal with the highest energy must exceed 80% of the energy deposited in the corresponding $5\times 5$ array of crystals. We require that the $\gamma\gamma$ invariant mass be within $20{\mathrm{\,Me\kern-1.00006ptV\!/}c^{2}}$ (about $3.5\sigma$ in resolution) of the $\pi^{0}$ mass PDG . To improve the $\pi^{0}$ momentum resolution, we perform a mass-constrained fit and require that the resulting $\chi^{2}$ be less than 30. This requirement is relatively loose, retaining more than 99% of candidates.

We subsequently combine $\pi^{0}$ candidates with two oppositely charged tracks, identified as a proton-antiproton pair. Such tracks are identified using requirements on the distance of closest approach with respect to the interaction point along the $z$ axis (antiparallel to the $e^{+}$ beam) of $|dz|<3.0$ cm, and in the transverse plane of $dr<0.3$ cm. In addition, charged tracks are required to have a minimum number of SVD hits ($>2$ along the $z$ axis and $>1$ in the transverse direction). Particle identification is achieved using information from the CDC, the TOF, and the ACC subdetectors. This information is combined to form a hadron likelihood $\mathcal{L}_{h}$; a charged track with likelihood ratios of $\mathcal{L}_{p}/(\mathcal{L}_{p}+\mathcal{L}_{K})>0.9$ and $\mathcal{L}_{p}/(\mathcal{L}_{p}+\mathcal{L}_{\pi})>0.9$ is regarded as a proton or antiproton. Furthermore, we reject tracks consistent with either the electron or muon hypothesis. The proton identification efficiency is 75% and the probability for a kaon (pion) to be misidentified as a proton is 6% (2%).

Candidate $B^{0}$ mesons are identified using the beam-energy-constrained mass, $M_{\rm bc}=\sqrt{E_{\rm beam}^{2}-|\vec{p}_{B}c|^{2}}/c^{2}$, and the energy difference $\Delta E=E_{B}-E_{\rm beam}$, where $E_{\rm beam}$ is the beam energy, and $E_{B}$ and $\vec{p}_{B}$ are the reconstructed energy and momentum, respectively, of the $B^{0}$ candidate. All quantities are evaluated in the center-of-mass (CM) frame. To improve the $M_{\rm bc}$ resolution, the momentum $\vec{p}_{B}$ is calculated as $\vec{p}_{B}=\vec{p}_{p}+\vec{p}_{\bar{p}}+\frac{\vec{p}_{\pi^{0}}}{|\vec{p}_{\pi^{0}}|}\sqrt{(E_{\rm beam}-E_{p}-E_{\bar{p}})^{2}/c^{2}-m^{2}_{\pi^{0}}c^{2}}$, where $m_{\pi^{0}}$ is the nominal $\pi^{0}$ mass PDG , $E_{h}$ and $\vec{p}_{h}$ are the energy and momentum of the hadron $h$ ($h=p,~{}\bar{p},~{}\pi^{0}$). In addition, a vertex fit is performed to the charged tracks to form a $B^{0}$ vertex. We require that the $\chi^{2}$ from the fit be less than 200. Events with $M_{\rm bc}>5.25~{}{\mathrm{\,Ge\kern-1.00006ptV\!/}c^{2}}$ and $-0.20~{}\mathrm{\,Ge\kern-1.00006ptV}<\Delta E<0.15~{}\mathrm{\,Ge\kern-1.00006ptV}$ are retained for further analysis. The signal yield is calculated in a smaller region $M_{\rm bc}\in(5.272,5.286)$ ${\mathrm{\,Ge\kern-1.00006ptV\!/}c^{2}}$ and $\Delta E\in(-0.12,+0.06)$ $\mathrm{\,Ge\kern-1.00006ptV}$. In order to reject contributions from charmonium states ($e.g.$, $\eta_{c}$, ${J\mskip-3.0mu/\mskip-2.0mu\psi\mskip 2.0mu}$, $\psi(2S)$, $\chi_{c0}$, $\chi_{c1}$ and $\chi_{c2}$), we apply a “charmonium veto” and exclude the regions of $2.850~{}{\mathrm{\,Ge\kern-1.00006ptV\!/}c^{2}}<m(p\bar{p})<3.128~{}{\mathrm{\,Ge\kern-1.00006ptV\!/}c^{2}}$ and $3.315~{}{\mathrm{\,Ge\kern-1.00006ptV\!/}c^{2}}<m(p\bar{p})<3.735~{}{\mathrm{\,Ge\kern-1.00006ptV\!/}c^{2}}$ from the event sample.

Charmless hadronic decays suffer from large amount of continuum background, arising from light quark production ($e^{+}e^{-}\rightarrow q\overline{}q,~{}q=u,d,s,c$). To suppress this background, we use a multivariate analyzer based on a neural network (NN) Feindt:2006pm that distinguishes jet-like continuum events from more spherical $B\kern 1.79993pt\overline{\kern-1.79993ptB}{}$ events. The NN uses the following input variables: the cosine of the angle between the thrust axis Brandt:1964sa of the $B^{0}$ candidate and the thrust axis of the rest of the event; the cosine of the angle between the $B^{0}$ thrust axis and the $+z$ axis; the cosine of the angle between the $+z$ axis and the $B^{0}$ candidate flight direction; a set of 18 modified Fox-Wolfram moments SFW ; the ratio of the second to zeroth (unmodified) Fox-Wolfram moments; the separation along the $z$ axis between the two $B$ vertices; and the $B$-flavor tagging information Kakuno:2004cf . All but for the last two quantities are evaluated in the CM frame. The NN is trained using MC simulated signal events and $q\overline{}q$ background events. The NN generates a single output variable ($C_{\rm NN}$) that ranges from $-1$ for background-like events to $+1$ for signal-like events. We require $C_{\rm NN}>-0.5$, which rejects approximately 86% of the $q\overline{}q$ background while retaining 94% of the signal. We then translate $C_{\rm NN}$ to a new variable

[TABLE]

where $C^{\rm min}_{\rm NN}=-0.5$ and $C_{\rm NN}^{\rm max}=1.0$. This translation is advantageous as the $C^{\prime}_{\rm NN}$ distribution for both signal and background is well described by a sum of Gaussian functions.

After applying all selection criteria, approximately 7% of the events have multiple $B^{0}$ candidates. For these events, we retain the candidate having the smallest sum of $\chi^{2}$ values obtained from the $\pi^{0}\rightarrow\gamma\gamma$ mass-constrained fit and the $B^{0}$ vertex-constrained fit. According to MC simulation, this criterion selects the correct $B^{0}$ candidate in 83% of multiple-candidate events.

We measure the signal yield by performing an unbinned extended maximum likelihood fit to the variables $M_{\rm bc}$, $\Delta E$, and $C^{\prime}_{\rm NN}$. The likelihood function is defined as

[TABLE]

where $Y_{j}$ is the yield of component $j$; $\mathcal{P}_{j}(M_{\rm bc}^{i},\Delta E^{i},C^{\prime i}_{\rm NN})$ is the probability density function (PDF) of component $j$ for event $i$; $j$ runs over all signal and background components; and $i$ runs over all events in the sample ($N$). The background components consist of continuum events, $b\rightarrow c$ (generic $B$) processes, and rare charmless processes. The latter two backgrounds are small compared to the continuum events and are studied using MC simulations. The rare charmless background shows a peaking structure in the $M_{\rm bc}$ distribution, most of which arises from $B^{+}\rightarrow p\bar{p}\rho^{+}$ decays. As correlations among the variables $M_{\rm bc}$, $\Delta E$, and $C^{\prime}_{\rm NN}$ are found to be small, the three-dimensional PDFs $\mathcal{P}_{j}(M_{\rm bc}^{i},\Delta E^{i},C^{\prime i}_{\rm NN})$ are factorized into the product of separate one-dimensional PDFs.

The PDF of signal events consists of two parts: one for candidates that are correctly reconstructed, and one for those incorrectly reconstructed, $i.e.$, at least one daughter originates from the other (tag-side) $B$. For the former case, the $M_{\rm bc}$ and $\Delta E$ distributions are modeled with Gaussian and Crystal Ball (CB) Skwarnicki:1986xj functions, respectively, while the $C^{\prime}_{\rm NN}$ distribution is modeled with a sum of Gaussian and bifurcated Gaussian functions having a common mean. The peak positions and resolutions of the $M_{\rm bc}$, $\Delta E$, and $C^{\prime}_{\rm NN}$ PDFs are adjusted to account for data-MC differences observed in a high-statistics control sample of $B^{0}\rightarrow\kern 1.99997pt\overline{\kern-1.99997ptD}{}^{0}(\rightarrow K^{+}\pi^{-})\pi^{0}$ decays. For the latter case, the correlated two-dimensional $M_{\rm bc}$-$\Delta E$ distribution is modeled with a non-parametric PDF Cranmer:2000du , and the $C^{\prime}_{\rm NN}$ component is modeled with a Gaussian function. The fraction of incorrectly reconstructed decays ($\sim 4\%$ in the signal region) is taken from MC simulation. For the rare charmless background, the $C^{\prime}_{\rm NN}$ component is modeled with a bifurcated Gaussian function. The $M_{\rm bc}$ and $\Delta E$ components are modeled by a joint two-dimensional non-parametric PDF. We model the $M_{\rm bc}$, $\Delta E$, and $C^{\prime}_{\rm NN}$ distributions of continuum background with an ARGUS Albrecht:1990am function having its endpoint fixed to $5.29~{}{\mathrm{\,Ge\kern-1.00006ptV\!/}c^{2}}$, a first-order polynomial, and a sum of two Gaussians having a common mean, respectively. For the generic $B$ background, we use a bifurcated Gaussian function to model the $C^{\prime}_{\rm NN}$ shape, while the similar shapes as of continuum background are used to model the $M_{\rm bc}$ and $\Delta E$ distributions. In addition to the fitted yields $Y_{j}$, all shape parameters for continuum background are also floated. All other parameters are fixed to the corresponding MC values.

The projections of the fit are shown in Fig. 1.

From the fit, we extract $40.5\pm 14.2$ signal events, $1490.3\pm 34.5$ continuum, $100.6\pm 35.0$ generic $B$, and $6.5\pm 10.1$ rare charmless background events in the $M_{\rm bc}-\Delta E$ signal region. The resulting branching fraction is calculated as

[TABLE]

where $Y_{\rm sig}$ represents the extracted signal yield, $N_{B\kern 1.25995pt\overline{\kern-1.25995ptB}{}}=(772\pm 11)\times 10^{6}$ is the total number of $B\kern 1.79993pt\overline{\kern-1.79993ptB}{}$ events, $\varepsilon=(10.53\pm 0.04)\%$ is the reconstruction efficiency. The efficiency is corrected to account for possible differences in particle identification (PID) and $\pi^{0}$ detection efficiencies between data and simulations. In Eq. (3) we assume equal production of $B^{0}\kern 1.79993pt\overline{\kern-1.79993ptB}{}^{0}$ and $B^{+}B^{-}$ pairs at the $\mathchar 28935\relax{(4S)}$ resonance. The result is

[TABLE]

where the first uncertainty is statistical and the second is systematic. This is the first measurement of this branching fraction.

The signal significance is calculated as $\sqrt{-2\ln(\mathcal{L}_{0}/\mathcal{L}_{\rm max})}$, where $\mathcal{L}_{0}$ is the likelihood value when the signal yield is fixed to zero, and $\mathcal{L}_{\rm max}$ is the likelihood value of the nominal fit. To include systematic uncertainties in the significance, we convolve the likelihood distribution with a Gaussian function whose width is set to the total systematic uncertainty that affects the signal yield. The resulting significance is 3.1 standard deviations. Thus, our measurement constitutes the first evidence for this decay mode.

The systematic uncertainty in ${\cal B}(B^{0}\rightarrow p\bar{p}\pi^{0})$ arises from several sources, as listed in Table 1. The uncertainty due to the fixed parameters in the PDF is estimated by varying them individually according to their statistical uncertainties. For each variation, the branching fraction is recalculated, and the difference with the nominal value is taken as the systematic uncertainty associated with that parameter. The smoothing parameters of the non-parametric functions are also varied. The differences in the fit results are included as systematic uncertainties. We add all uncertainties in quadrature to obtain the overall uncertainty due to PDF parametrization. The uncertainties due to errors in the calibration factors used to account for data-MC differences in the signal PDF are evaluated separately but in a similar manner. To test the stability of our fitting procedure, we generate and fit a large ensemble of pseudoexperiments. We find a potential fit bias of $+2.1$%. We attribute this bias to neglecting small correlations among the fitted observables. We assign a 1.5% systematic uncertainty due to $\pi^{0}$ reconstruction; this is determined from a study of $\tau^{-}\rightarrow\pi^{-}\pi^{0}\nu_{\tau}$ decays Ryu:2014vpc . The systematic uncertainty due to the track reconstruction efficiency is 0.35% per track, as determined from a study of partially reconstructed $D^{*+}\rightarrow D^{0}\pi^{+},D^{0}\rightarrow K^{0}_{S}\pi^{+}\pi^{-}$ decays. A 0.6% systematic uncertainty is assigned due to the particle identification efficiency of the proton-antiproton pair; this is determined from a study of $\Lambda\rightarrow p\pi^{-}$ decays. We determine the systematic uncertainty due to the $C_{\rm NN}$ selection by applying different $C_{\rm NN}$ criteria and comparing the results with that of the $C_{\rm NN}$ nominal selection. The uncertainty due to the estimated fraction of incorrectly reconstructed signal events is obtained by varying this fraction by $\pm 50\%$. The systematic uncertainty due to the total number of $B\kern 1.79993pt\overline{\kern-1.79993ptB}{}$ pairs is 1.4%, and the uncertainty due to MC used to evaluate the reconstructed efficiency is 0.4%. The total systematic uncertainty is obtained by adding each source in quadrature, as they are assumed to be uncorrelated.

Figure 2 shows the background-subtracted and efficiency-corrected distribution of $m(p\bar{p})$, where the charmonium veto is removed. For the background subtraction, we use the sPlot technique Pivk:2004ty , with $M_{\rm bc}$, $\Delta E$ and $C^{\prime}_{\rm NN}$ as the discriminating variables.

As expected, an enhancement near threshold is visible. The background-subtracted distributions of $m(p\pi^{0})$ and $m(\bar{p}\pi^{0})$ are shown in Fig. 3.

No obvious structure is seen in these distributions.

We also search for the intermediate two-body decay $B^{0}\rightarrow\Delta^{+}(\rightarrow p\pi^{0})\bar{p}$. Events with $m(p\pi^{0})<1.4$ ${\mathrm{\,Ge\kern-1.00006ptV\!/}c^{2}}$ are selected for this search. No significant signal is observed in this mass range. We set an upper limit on the branching fraction of $B^{0}\rightarrow\Delta^{+}\bar{p}$ at 90% confidence level (C.L.) using a Bayesian approach. The limit is obtained by integrating the likelihood function from zero to infinity; the value that corresponds to 90% of this total area is taken as the 90% C.L. upper limit. We include the systematic uncertainty in the calculation by convolving the likelihood distribution with a Gaussian function whose width is set equal to the total systematic uncertainty of $\mathcal{B}(B^{0}\rightarrow p\bar{p}\pi^{0})$. As we do not know the flavor of the $B$ meson at decay, we express our result as a sum of final states containing either a $\Delta^{+}$ or a $\bar{\Delta}^{-}$. The result is

[TABLE]

This is the first such limit and is in agreement with the theoretical predictions Chernyak:1990ag ; Cheng:2001tr .

In summary, using the full set of Belle data, we report a measurement of the branching fraction for $B^{0}\rightarrow p\bar{p}\pi^{0}$ decays. We obtain $\mathcal{B}(B^{0}\rightarrow p\bar{p}\pi^{0})=(5.0\pm 1.8\pm 0.6)\times 10^{-7}$, where the first uncertainty is statistical and the second is systematic. The significance of this result is 3.1 standard deviations, and thus this measurement constitutes the first evidence for this decay. We also search for the intermediate two-body decays $B^{0}\rightarrow\Delta^{+}\bar{p}$ and $B^{0}\rightarrow\bar{\Delta}^{-}p$, and set an upper limit on the branching fraction, $\mathcal{B}(B^{0}\rightarrow\Delta^{+}\bar{p})+\mathcal{B}(B^{0}\rightarrow\bar{\Delta}^{-}p)<1.6\times 10^{-6}$ at 90% C.L.

We thank the KEKB group for the excellent operation of the accelerator; the KEK cryogenics group for the efficient operation of the solenoid; and the KEK computer group, and the Pacific Northwest National Laboratory (PNNL) Environmental Molecular Sciences Laboratory (EMSL) computing group for strong computing support; and the National Institute of Informatics, and Science Information NETwork 5 (SINET5) for valuable network support. We acknowledge support from the Ministry of Education, Culture, Sports, Science, and Technology (MEXT) of Japan, the Japan Society for the Promotion of Science (JSPS), and the Tau-Lepton Physics Research Center of Nagoya University; the Australian Research Council including grants DP180102629, DP170102389, DP170102204, DP150103061, FT130100303; Austrian Science Fund under Grant No. P 26794-N20; the National Natural Science Foundation of China under Contracts No. 11435013, No. 11475187, No. 11521505, No. 11575017, No. 11675166, No. 11705209; Key Research Program of Frontier Sciences, Chinese Academy of Sciences (CAS), Grant No. QYZDJ-SSW-SLH011; the CAS Center for Excellence in Particle Physics (CCEPP); the Shanghai Pujiang Program under Grant No. 18PJ1401000; the Ministry of Education, Youth and Sports of the Czech Republic under Contract No. LTT17020; the Carl Zeiss Foundation, the Deutsche Forschungsgemeinschaft, the Excellence Cluster Universe, and the VolkswagenStiftung; the Department of Science and Technology of India; the Istituto Nazionale di Fisica Nucleare of Italy; National Research Foundation (NRF) of Korea Grants No. 2015H1A2A1033649, No. 2016R1D1A1B01010135, No. 2016K1A3A7A09005 603, No. 2016R1D1A1B02012900, No. 2018R1A2B3003 643, No. 2018R1A6A1A06024970, No. 2018R1D1 A1B07047294; Radiation Science Research Institute, Foreign Large-size Research Facility Application Supporting project, the Global Science Experimental Data Hub Center of the Korea Institute of Science and Technology Information and KREONET/GLORIAD; the Polish Ministry of Science and Higher Education and the National Science Center; the Grant of the Russian Federation Government, Agreement No. 14.W03.31.0026; the Slovenian Research Agency; Ikerbasque, Basque Foundation for Science, Spain; the Swiss National Science Foundation; the Ministry of Education and the Ministry of Science and Technology of Taiwan; and the United States Department of Energy and the National Science Foundation.