Measurement of TeV atmospheric muon charge ratio with the full OPERA data

The OPERA detector, designed to search for $\nu_{\mu} \to \nu_{\tau}$ oscillations in the CNGS beam, is located in the underground Gran Sasso laboratory, a privileged location to study TeV-scale cosmic rays. For the analysis here presented, the detector was used to measure the atmospheric muon charge ratio in the TeV region. OPERA collected charge-separated cosmic ray data between 2008 and 2012. More than 3 million atmospheric muon events were detected and reconstructed, among which about 110000 multiple muon bundles. The charge ratio $R_{\mu} \equiv N_{\mu^+}/N_{\mu^-}$ was measured separately for single and for multiple muon events. The analysis exploited the inversion of the magnet polarity which was performed on purpose during the 2012 Run. The combination of the two data sets with opposite magnet polarities allowed minimizing systematic uncertainties and reaching an accurate determination of the muon charge ratio. Data were fitted to obtain relevant parameters on the composition of primary cosmic rays and the associated kaon production in the forward fragmentation region. In the surface energy range 1-20 TeV investigated by OPERA, $R_{\mu}$ is well described by a parametric model including only pion and kaon contributions to the muon flux, showing no significant contribution of the prompt component. The energy independence supports the validity of Feynman scaling in the fragmentation region up to $200$ TeV/nucleon primary energy.

To the memory of Prof. G. Giacomelli Abstract The OPERA detector, designed to search for ν µ → ν τ oscillations in the CNGS beam, is located in the underground Gran Sasso laboratory, a privileged location to study TeV-scale cosmic rays. For the analysis here presented, the detector was used to measure the atmospheric muon charge ratio in the TeV region. OPERA collected charge-separated cosmic ray data between 2008 and 2012. More than 3 million atmospheric muon events were detected and reconstructed, among which about 110000 multiple muon bundles. The charge ratio R µ ≡ N µ + /N µ − was measured separately for single and for multiple muon events. The analysis exploited the inversion of the magnet polarity which was performed on purpose during the 2012 Run. The combination of the two data sets with opposite magnet polarities allowed minimizing systematic uncertainties and reaching an accurate determination of the muon charge ratio. Data were fitted to obtain relevant parameters on the composition of primary cosmic rays and the associated kaon production in the forward fragmentation region. In the surface energy range 1-20 TeV investigated by OPERA, R µ is well described by a parametric model including only pion and kaon contributions to the muon flux, showing no significant contribution of the prompt component. The energy independence supports the validity of Feynman scaling in the fragmentation region up to 200 TeV/nucleon primary energy.

Introduction
Underground experiments detect the penetrating remnants of primary cosmic ray interactions in the atmosphere, namely muons and neutrinos. These are the decay products of charged mesons contained in the particle cascade, mainly pions and kaons. At very high energies also charmed particles are expected to contribute.
The muon charge ratio R µ ≡ N µ + /N µ − , defined as the number of positive over negative charged muons, is studied since many decades. It provides an understanding of the mechanism of multiparticle production in the atmosphere in kinematic regions not accessible to accelerators, as well as information on the primary cosmic ray composition. A charge ratio larger than unity reflects the abundance of protons over heavier nuclei in the primary cosmic radiation. The charge asymmetry is preserved in the secondary hadron production, and consequently in the muon fluxes, due to the steepness of the primary spectrum which enhances the forward fragmentation region [1]. The kaon contribution to the muon flux increases with the muon energy. Since the pro-Corresponding authors. E-mail: mauri@bo.infn.it, sioli@bo.infn.it a Now at Kyungpook National University, Daegu, Korea. b Now at Samsung Changwon Hospital, SKKU, Changwon, Korea. duction of positive kaons is favoured by the associated production Λ K + , the muon charge ratio is expected to rise with energy. Assuming the hypothesis of complete scaling we expect an energy independent charge ratio above the TeV energy region at sea level [1] once the kaon contribution to the muon flux reached its asymptotic value [2]. At higher energies, around O(100) TeV, the heavy flavor contribution, as well as changes in the primary composition, may become significant.
The OPERA experiment, described in detail in Ref. [3], is a hybrid electronic detector/emulsion apparatus, located in the underground Gran Sasso laboratory, at an average depth of 3800 meters of water equivalent. The main physics goal of the experiment is the first observation of neutrino oscillations in direct appearance mode in the ν µ → ν τ channel [4][5][6]. OPERA already reported a first measurement of the atmospheric muon charge ratio at TeV surface energies using the 2008 Run data [7]. Here we present the final results obtained with the complete statistics. OPERA continuously accumulated cosmic ray data with the electronic detectors of the target over the whole year from 2008 up to 2012. However the magnetic spectrometers were active only during the CNGS Physics Runs, being switched off during the CNGS winter shutdowns.
As it was done in Ref. [7], we used the momentum and charge reconstruction obtained via the Precision Trackers (PT) of the OPERA spectrometers [8]. Layers of vertical drift tubes are arranged in PT stations instrumenting the two identical dipole magnets. The momentum and charge information is given by the angle ∆ φ in the bending plane, i.e. the difference between the track directions reconstructed by the two PT stations before and after each magnet arm. For nearly horizontal muons up to four bending angles can be measured in the two dipole magnets.

Data Analysis
The cosmic ray data used for this analysis were collected during the five CNGS Physics Runs between 2008 and 2012. In the first four years (2008-2011) the magnetic field was directed upward in the first arm of both dipoles and in the opposite direction in the second arm (standard polarity, SP). In 2012 the coil currents were reversed and the spectrometer operated in inverted polarity (IP) mode.
A pre-selection was applied in order to select only stable conditions of detector operation. Short periods with increased electronic noise or with any subdetector under test were removed, as well as periods in which the magnets were not in nominal conditions. Details on atmospheric muon event selection, reconstruction and analysis can be found in Refs. [7,9]. In the total SP + IP live time, 3044281 cosmic muon events were recorded. Among them, 113662 are muon bundles, i.e. events with a muon multiplicity n µ greater than 1. To reconstruct the muon charge, the track has to cross at least one magnet arm yielding a measurement of the bending angle ∆ φ by the PT system. This resulted in the reconstruction of momentum and charge for 650492 muons in SP (28.7 % of the total muon events) and 244626 muons in IP (28.9 % of the total muon events).
In order to improve the charge identification purity, the two selection criteria used in Ref. [7] were applied to the data. The first selection is a track quality cut. The ∆ φ bending angle measurement is provided by the PT track reconstruction which is spoiled in events containing a large number of fired tubes, typically due to radiative processes. When the number of PT hits is much larger than the number expected from geometrical considerations [7,9] the event is rejected.
The second selection acts on the charge discrimination power. Events with a bending angle smaller than 3 times the angular resolution were rejected. This corresponds to a maximum detectable momentum up to 1 TeV/c [9]. A further cut was applied to remove a few events with very large deflections (|∆ φ | > 100 mrad), either due to the scattering of low momentum muons (p µ ≤ 5 GeV/c) or mimicked by secondary particles produced in high energy events.
Muons induced by atmospheric neutrinos coming from below were removed from the data set on the basis of timeof-flight measurements. Contributions from muon backscattering or up-going charged particles induced by muons were computed according to Ref. [10] and found to be negligible.
The numbers of single and multiple muons surviving all the selection cuts and used in the computation of the muon charge ratio are reported in Table 2.

Systematic uncertainties and unbiased charge ratio
The comparison between the two data sets with opposite magnet polarity (SP and IP) allows checking systematic un- 143628  105278  5252  4533  IP  53575  40086  1785  1740   Table 2 Final statistics for the muon charge ratio measurement; the number of muons surviving the cuts is quoted for both magnet polarity configurations. For muon bundles we provide the total number of muons and not the number of events.
certainties affecting the muon charge ratio. These can be cancelled out using a proper combination of the two data samples (see Appendix A).
The two main sources of systematic uncertainties are due to alignment and charge misidentification.
In principle a different acceptance for µ + and µ − could also contribute to the overall systematic uncertainty. However the symmetry of the detector geometry allows to safely neglect this contribution. An indirect confirmation is given by the compatibility of the charge ratio values computed separately in the two arms of the same magnet, where the magnetic field has opposite directions [9].
Using the SP and IP data sets, we checked the symmetry in the acceptance for each magnet arm. According to the reference frame defined in Ref. [7], where the z-axis points toward the CNGS direction, muons travelling toward the positive z-axis are defined as south-oriented (SO), while muons travelling toward the negative z-axis are defined as northoriented (NO). A muon crosses a magnet arm in one of these two possible "orientations". South-oriented µ + and northoriented µ − are deflected toward west in the first arm in SP mode. The reversals of either the muon incoming orientation or the polarity mode are equivalent ways to exchange the muon bending sign. We computed the ratio A i of the number N i of charge-reconstructed muons in SP mode to the number in IP mode (normalized by their polarity live time), A i = (N i ) SP /(N i ) IP for the two orientations in each magnet arm i. The results are reported in Table 3. The values of A i obtained in one orientation are all compatible with the values obtained in the other orientation, as expected from a charge-symmetric spectrometer. The individual comparison between A i (SO) and A i (NO) for each arm disposes of possible small live time differences among PT stations. The results are consistent with unity within statistical errors.
We have investigated the systematic uncertainty related to the alignment of the PT system. The SP and IP bending angle distributions were compared separately for south-and north-oriented muons in each magnet arm. In case of perfect alignment, the two distributions (normalized by their respective live times) would coincide. In the data, a systematic bending angle shift |δ φ s | ∼ (0.10 ± 0.03) mrad was observed on average (in each magnet arm, for i = 1, . . . , 4: |δ φ s,i | = {< 0.03, 0.07, 0.10, 0.15} mrad). Inverting the muon  Table 3 Ratio between SP and IP numbers of charge-reconstructed muons in each magnet arm. The normalization by the relative polarity live time is globally applied.
orientation, δ φ s,i preserves the absolute value and flips the sign, as expected in case of misalignment. Note that the absolute value is compatible with the alignment systematic uncertainty δ φ syst = 0.08 mrad given in [7]. The observed global shift δ φ s is however an average value. It is a cumulative result of local distortions, tilts and bendings which depend on the muon position, zenith and azimuth. We therefore did not globally correct for δ φ s since the combination of IP and SP data allows to completely remove this systematics at a local level. As detailed in the Appendix A, the unbiased charge ratioR µ is obtained by the normalized sum of µ + over the normalized sum of µ − : where l SP,IP is the respective polarity live time and η is the charge misidentification probability. This combination provides a charge ratio in which the effects induced by misalignments cancel out. Indeed, the last equation is exactly the relation between the reconstructedR µ and the true R µ charge ratio in case of perfect alignment [7]. Inverting this relation, the charge ratio R µ is obtained from the measured R µ corrected by the misidentification probability.
In principle, all the systematic contributions due to misalignment cancel with this combination of SP and IP data. The residual systematic errors which do not cancel are estimated by the difference between the charge ratio values computed separately for SO and NO orientations. Since the alignment bias has opposite sign in the two orientations, we take |R µ (NO) − R µ (SO)| as the systematic uncertainty related to our combination procedure. It was found δ R µ = 0.001 for single muon events and δ R µ = 0.013 for multiple muon events. In the latter the statistical contribution is dominant.
The second source of systematic uncertainty considered is related to the determination of η. The charge misidentification computed with Monte Carlo is η MC = 0.030, nearly independent on the muon momentum in the range 5 GeV/c < ∼ p µ < ∼ 1 TeV/c [9]. We estimated the systematic uncertainty of η using a subsample of experimental data, i.e. the muon tracks reconstructed in both arms of each spectrometer. The probability of wrong charge assignment was evaluated counting the fraction of tracks with different charges, and the experimental η data was derived. The difference between η data and η MC is δ η = 0.007 ± 0.002 [7]. This corre-sponds to a one-sided systematic uncertainty on the charge ratio δ R µ = 0.007.
The final systematic uncertainty is the quadratic sum of the misalignment and the misidentification contributions.

Results
The charge ratio of single muons impinging on the apparatus was computed combining the two polarity data sets according to Eq. 1. After the correction for charge misidentification and detector misalignment, the final measurement with the complete 5-year statistics yields the result: The charge ratio of multiple muon events was computed using all the muon charges reconstructed in events with n µ > 1. It is not computed within the bundle itself, but summing up all the positive and the negative charges belonging to the bundle subsample. The result after polarity combination and correction for misidentification is significantly lower than the single muon value: R µ (n µ > 1) = 1.098 ± 0.023(stat.) +0.015 −0.013 (syst.) The smaller value of the charge ratio for multiple muon events originates from two effects. First, as pointed out in [7], the multiple muon sample naturally selects heavier primaries, thus a neutron enriched primary beam ( A 3.4 for single muons, A 8.5 for bundles). Second, the selection of muon bundles biases the Feynman-x distribution towards the central region (x F E secondary /E primary → 0), in which the sea quark contribution to secondary particle production becomes relevant [9]. Both processes cause a decrease in the charge ratio.
The single muon charge ratio was projected at the Earth surface using a Monte Carlo based unfolding technique for the muon energy E µ [9]. As a first attempt, only pion and kaon contributions to the total muon flux are considered. We used the analytic approximation described in [7] to infer the fractions of charged mesons decaying into a positive muon, f π + and f K + . This approach does not yet consider any energy dependence of the proton excess in the primary composition. In this case the muon flux and charge ratio depend on the vertical surface energy E µ cos θ * , where θ * is the zenith angle at the muon production point [11].
R µ is computed as a function of the vertical surface muon energy, binned according to the energy resolution, which  Fig. 1 The muon charge ratio measured by OPERA as a function of the vertical surface energy E µ cos θ * (black points). Our data are fitted together with the L3+C [15] data (open triangles). The fit result is shown by the continuous line. The dashed, dotted and dash-dot lines are, respectively, the fit results with the inclusion of the RQPM [21], QGSM [21] and VFGS [22] models for prompt muon production in the atmosphere. The vertical inner bars denote the statistical uncertainty, the full bars show the total uncertainty. Results from other experiments, MINOS Near and Far Detectors [16,17], CMS [18] and Utah [19], are shown for comparison.
is of the order of d(log 10 E µ /GeV) 0.15 in a logarithmic scale [9]. In each bin the two polarity data sets are combined and the obtained value is corrected for the charge misidentification. The two contributions to the systematic uncertainty are computed and added in quadrature. The results are shown in Fig. 1, together with data from other experiments (L3+C [15], MINOS Near and Far Detectors [16,17], CMS [18] and Utah [19]). The information for each of the four E µ cos θ * bins are presented in Table 4: the energy range, the most probable value of the energy distribution in the bin, the average zenith angle, the charge ratio R µ , the statistical and systematic uncertainties.
Following the procedure described in [7], we fitted our data and those from [15] (for the high and low energy regions) in order to infer the fractions f π + and f K + . In this approach, the atmospheric charged kaon/pion production ratio R K/π had to be fixed. For this, we took the weighted average of experimental values reviewed in [20], R K/π = 0.127. The fit yields f π + = 0.5512 ± 0.0014 and f K + = 0.705 ± 0.014, corresponding to a muon charge ratio from pion decay R π = 1.2281 ± 0.0007 and a muon charge ratio from kaon decay R K = 2.39 ± 0.07.
Taking into account various models for charm production, namely RQPM [21], QGSM [21] and VFGS [22], the positive pion and kaon fractions obtained from the fit are unchanged within statistical errors. The results are shown in Fig. 1. The prompt muon component does not significantly contribute to R µ up to E µ cos θ * < ∼ 10 TeV.
Recently, an enlightening analytic description of the muon charge ratio considering an explicit dependence on the relative proton excess in the primary cosmic rays, δ 0 = (p − n)/(p + n), was presented in [2]: Here p and n fluxes are defined as where the index i runs over the primary ions (H, He, CNO, Mg-Si, Fe) and E N is the primary nucleon energy. The contributions from decays of pions and kaons are included in the kinematic factors A i , B i , ε i (i = π, K) described in [2,11]. An analogous contribution from charm decay is foreseen at high energies but still not observed. The spectrum weighted moments Z i j [2] are contained in β and α K : Isospin symmetry allows expressing the pion contribution in terms of f π + , where Here α π is obtained replacing the subscript K with the subscript π in α K . We extracted from the data the composition parameter δ 0 and the factor Z pK + related to the associated production Λ K + in the forward region. The Z pK + moment is still poorly known and its predicted value considerably differs for different Monte Carlo codes [12,13].
In Eq. 2 the charge ratio does not exclusively depend on the vertical surface energy. Since the spectra of primary nuclei have different spectral indices, the parameter δ 0 depends on the primary nucleon energy E N . In the energy range of interest the approximation E N 10 × E µ can be used [2].
The correct way of taking into account the different dependencies is to simultaneously fit Eq. 2 as a function of the two variables (E µ , cos θ * ). In each (E µ , cos θ * ) bin the data sets with opposite polarities are combined andR µ is corrected for the charge misidentification.
The pion moments Z pπ + and Z pπ − were set to the values reported in [2], since the fraction of positive pions in the atmosphere f π + = 0.5512 ± 0.0014 derived in this work is robust and consistent with previous measurements [16,17] and with the Z Nπ values based on fixed target data [14]. The moment Z pK − was also set to the value given in [2], since for K − there is no counterpart of the associated production Λ K + . On the other hand K − are equally produced in K + K − pairs by protons and neutrons (Z pK − Z nK − ).  Table 4 The charge ratio in bins of E µ cos θ * . Here reported are the energy bin range, the most probable value of the energy distribution in the bin (MPV, evaluated using the full Monte Carlo simulation described in [9]), the average zenith angle, the charge ratio and the statistical and systematic uncertainties.  Fig. 2 Our measurement of the muon charge ratio as a function of the surface energy E µ (black points). The two-dimensional fit in (E µ , cos θ * ) yields a measurement of the composition parameter δ 0 and of the factor Z pK + . The fit result is projected on the average OPERA zenith cos θ * 0.7 and shown by the continuous line. Results from other experiments, L3+C (only for 0.675 < cos θ < 0.75) [15], MI-NOS Near and Far Detectors [16,17], CMS [18] and Utah [19], are also shown for comparison.
A linear energy dependence in logarithmic scale of the parameter δ 0 was assumed, δ 0 = a + b log 10 (E N /GeV/nucleon), as suggested by direct measurements of the primary composition and by the Polygonato model [23]. We fixed the slope at b = −0.035 which was obtained fitting the values reported in [2].
We made a two-dimensional fit of OPERA and L3+C data as a function of (E µ , cos θ * ) to Eq. 2 with δ 0 and Z pK + as free parameters. The fit yields the composition parameter at the average energy measured by OPERA E µ = 2 TeV (corresponding to E N ≈ 20 TeV/nucleon) δ 0 ( E µ ) = 0.61± 0.02 and the factor Z pK + = 0.0086 ± 0.0004.
The result of the fit in two variables (E µ , cos θ * ) is projected on the average OPERA zenith cos θ * 0.7 and is shown in Fig. 2 together with the measured charge ratio as a function of the surface muon energy. The energy independence of the charge ratio above the TeV supports the validity of the Feynman scaling in the fragmentation region.

Conclusions
The atmospheric muon charge ratio R µ was measured with the complete statistics accumulated along the five years of data taking. The combination of the two data sets collected with opposite magnet polarities allows reaching the most accurate measurement in the high energy region to date. The underground charge ratio was evaluated separately for single and for multiple muon events. For single muons, the integrated R µ value is R µ (n µ = 1) = 1.377 ± 0.006(stat.) +0.007 −0.001 (syst.) while for muon bundles R µ (n µ > 1) = 1.098 ± 0.023(stat.) +0.015 −0.013 (syst.) The integral value and the energy dependence of the charge ratio for single muons are compatible with the expectation from a simple model [2,14] which takes into account only pion and kaon contributions to the atmospheric muon flux. We extracted the fractions of charged pions and kaons decaying into positive muons, f π + = 0.5512±0.0014 and f K + = 0.705 ± 0.014.
Considering the composition dependence embedded in Eq. 2, we inferred a proton excess in the primary cosmic rays δ 0 = 0.61 ± 0.02 at the energy E N ≈ 20 TeV/nucleon and a spectrum weighted moment Z pK + = 0.0086 ± 0.0004.
The observed behaviour of R µ as a function of the surface energy from ∼ 1 TeV up to 20 TeV (about 200 TeV/nucleon for the primary particle) shows no deviations from a simple parametric model taking into account only pions and kaons as muon parents, supporting the hypothesis of limiting fragmentation up to primary energies/nucleon around 200 TeV. We are indebted to our technical collaborators for the excellent quality of their work over many years of design, prototyping and construction of the detector and of its facilities.

Appendix A: Combination of data sets
A systematic shift of the bending angle distribution biases the integral value of the muon charge ratio. Moreover, since the unfolding of the surface muon energy is based on the underground muon momentum, a curvature bias has an important effect on the bin-to-bin migration matrix, i.e. the probability of measuring a surface energy E i at a true energy E j .
Due to misalignment there are in principle two different migration matrices U + and U − for each magnet polarity. Given the symmetry of the detector, the exchange of the magnet polarity is equivalent to the exchange of the charge sign (see Sect. 2.1), thus U + SP = U − IP and coherently U + IP = U − SP . In general, with a curvature bias that shifts the bending angle distribution, a different charge misidentification η + and η − for positive and negative muons is expected for both standard and inverted magnet polarity. Given the symmetry of the detector, the relations η + SP = η − IP and η − SP = η + IP are valid. However we verified that after the application of the second selection criterion (the bending angle cut) the charge misidentification η is insensitive to the charge sign. We applied a rigid curvature bias δ φ s and observed that the bin construction clearly separates positive and negative bins. Therefore a symmetric misidentification η = η + = η − is assumed. Each energy bin content N ± is the integral of the true charged muon flux Φ µ convolved with the migration matrix U and corrected for the charge misidentification. For the standard polarity SP we have: where E 1 , E 2 are the lower and upper bounds of the reconstructed energy bin. The positive flux contribution can be rewritten in terms of the true charge ratio R µ and the negative flux: Writing the same equations for the inverted polarity IP, the symmetries described above are taken into account: Thanks to the symmetric detector setup, the data combination able to cancel the misalignment systematic errors is the ratio (N + SP + N + IP )/(N − SP + N − IP ), where the numbers are normalized by the respective polarity live times. Indeed, writing the integrands only, we obtain: Thus the unbiased charge ratio is given by the normalized sum of µ + over the normalized sum of µ − : The last equation is exactly the relation between the recon-structedR µ and true R µ charge ratio in case of perfect alignment [7].