Measurement of the atmospheric muon charge ratio with the OPERA detector

The OPERA detector at the Gran Sasso underground laboratory (LNGS) was used to measure the atmospheric muon charge ratio in the TeV energy region. We analyzed 403069 atmospheric muons corresponding to 113.4 days of livetime during the 2008 CNGS run. We computed separately the muon charge ratio for single and for multiple muon events in order to select different energy regions of the primary cosmic ray spectrum and to test the charge ratio dependence on the primary composition. The measured charge ratio values were corrected taking into account the charge-misidentification errors. Data have also been grouped in five bins of the"vertical surface energy". A fit to a simplified model of muon production in the atmosphere allowed the determination of the pion and kaon charge ratios weighted by the cosmic ray energy spectrum.


Introduction
Primary cosmic rays (typically protons) impinging on the Earth's atmosphere produce showers of secondary particles which propagate down to the ground level. Most of the interaction products are π and K mesons which in turn decay or interact, depending on their energy and on the air density profile they pass through. The decay of π 0 mesons gives rise to the electromagnetic component of the showers, the decay of π ± and K ± mesons yields mostly muons which are the most penetrating charged particles and therefore the most abundant charged component at sea level. In particular only the most energetic muons can penetrate deep underground. The Gran Sasso laboratory (LNGS) is located at an average depth of 3800 m.w.e. and the minimum muon energy required to reach the underground depth is around 1.5 TeV while the residual underground energy is about 270 GeV averaged over all directions [1].
The muon charge ratio R µ = N µ + /N µ − defined as the number of positive over negative charged muons, results from several contributions: the primary cosmic ray composition (in particular the ratio of protons over heavier primaries), hadronic-interaction features, atmospheric conditions (negligible above a few GeV) and, at very high energy, the contribution of muons from charmed particle decays (prompt muons) [2]. The muon charge ratio at sea level was extensively studied in the past since it is an indicator of important aspects of cosmic rays and particle physics.
An exhaustive compilation of measurements in a wide energy range is reported in Ref. [3]. In the interval from a few hundred MeV to 300 GeV the muon charge ratio R µ stays around 1.27. At higher energies several competing processes can affect its value. Since strong interaction production channels lead to a K + /K − ratio higher than the π + /π − one and the fraction of muons from kaon decays increases with the energy, the muon charge ratio is expected to rise as the energy increases. On the other hand, as the zenith angle increases and hence longer lived mesons have a higher probability to decay in the less dense layers of the high atmosphere the fraction of muons from pion decay increases and the muon charge ratio decreases. We also expect a dependence of the muon charge ratio on the underground muon multiplicity n µ , which is related to the energy of the primary cosmic rays and to their chemical composition. For primaries different from protons the positive charge excess is reduced and so is the muon charge ratio [4].
A simplified model of the atmospheric muon charge ratio is obtained from the muon spectrum [5] (1) where is the primary spectrum of nucleons (evaluated at the muon energy in the atmosphere E µ ) with a spectral index γ + 1 ≃ 2.7. For each of the N par muon parents (π, K, charmed particles etc.) the constants a i and b i contain the kinematical factors for the decay into muons, ǫ i (θ) are the critical energies defined as the energies above which interaction processes dominate over decay. They depend on the ratio m i /τ i (mass over rest lifetime of the muon parent) and on the atmospheric profile density and therefore on the zenith angle θ. A good approximation for ǫ i (θ) which takes into account the Earth's curvature is with cos θ * = 1 − sin 2 θ R e R e + h 2 (3) where R e is the Earth's radius and h is the muon production height. By using eq. 3 the zenith angle is evaluated at the muon production point and not at the detector site. By choosing h = 30 km an agreement within 5% with the precise ǫ i (θ) computation is obtained [6].
The spectrum weighted moments are defined as: where σ ij is the inclusive cross-section for the production of a particle j from the collision of a particle i with a nucleus in the atmosphere and x lab = E j /E i is the energy fraction carried by the secondary particle. Use of the Z factors shows explicitly that particle production in cosmic ray cascades is concentrated in the forward fragmentation region at large x lab . From Eq. 1 the i-th contribution to the muon flux is proportional to Z N i and is suppressed for E µ ≫ ǫ i (θ). Therefore each contribution to the muon charge ratio produced by different muon parents can be disentangled by studying the muon charge ratio as a function of the muon energy E µ .
Considering only π and K parent mesons Eq. 1 can be separated for µ + and µ − where ǫ π ≃ 115 GeV and ǫ K ≃ 850 GeV are the pion and kaon critical energies along the vertical direction, respectively. These energies have to be compared with the corresponding value for charmed particles, ǫ X > 10 7 GeV. The prompt muon component (from charmed particles) is therefore isotropically distributed since the corresponding cos θ * factor is suppressed, at least in the TeV region. Eq. 5 contains most of the aspects already discussed. First we note that the correct variable to describe the evolution of the charge ratio is the product E µ cos θ * , the "vertical surface energy" [7,8]. Moreover the two energy scales which determine the pion and kaon contributions to R µ are the critical energies ǫ π and ǫ K . The evaluation of the muon surface energy E µ depends on the rock depth crossed by the muon to reach the detector and therefore the distribution of E µ cos θ * is related to the shape of the overburden. Measurements of the muon charge ratio at high energies and large zenith angles, corresponding to E µ cos θ * ∼ 0.5 TeV, are given in Ref. [9]. More recent data with large statistics at E µ cos θ * ∼ 1 TeV are presented in Ref. [10]. These results suggest a smooth transition toward the energy region where kaon contribution becomes significant.
The LNGS laboratory is located at E µ cos θ * ≃ 2 TeV, well above the kaon critical energy ǫ K . This allows the measurement of the ratio Z N K + /Z N K − whose value is poorly known in the fragmentation region. This has also a strong impact, for instance on the evaluation of the flux of TeV atmospheric neutrinos, which are dominated by kaon production. OPERA is the first large magnetized detector that can measure the muon charge ratio at the LNGS depth, with an acceptance for cosmic ray muons coming from above A = 599 m 2 · sr (A = 197 m 2 sr for muons crossing the spectrometer sections).
The paper is organized as follows. In Sec. 2 the detector is briefly described, while in Sec. 3 we describe data selection and reconstruction, Monte Carlo simulation and data reduction. In Sec. 4 the muon charge ratio at the LNGS underground depth and its systematic error are evaluated. Finally in Sec. 5 we give R µ as a function of the underground momentum and of the variable E µ cos θ * fitted to Eq. 5.

The OPERA detector
OPERA is a hybrid experiment with electronic detectors and nuclear emulsions located in Hall C of the underground Gran Sasso Laboratory in central Italy [11]. The main physics goal of the experiment is to observe neutrino flavor oscillations through the appearance of ν τ neutrinos in the ν µ CNGS beam [12]. The detector design was optimized to identify the τ lepton via the topological observation of its decay: this requires a target mass of more than a kton to maximize the neutrino interaction probability and a micrometric resolution to detect the τ decay. To accomplish these requirements the detector concept is based on the Emulsion Cloud Chamber (ECC) technique combined with real-time detection techniques ("electronic detectors").
The ECC basic unit in OPERA is a "brick" made of 56 lead plates, 1 mm thick, providing the necessary mass to cope with the small neutrino cross-section, interleaved with 57 nuclear emulsion films (industrially produced), providing the necessary spatial and angular resolution to identify tau decay topologies. In total, 150000 bricks have been assembled reaching the overall mass of 1.25 kton.
The electronic detectors are used to trigger the neutrino interactions, to locate the brick in which the interaction occurred, to identify muons and measure their charge and momentum.
The detector is composed of two identical parts called supermodules (SM), each one consisting of a target section and a magnetic spectrometer. In the target the bricks are arranged in 29 vertical "walls" transverse to the beam direction interleaved with electronic Target Tracker (TT) Fig. 1: Schematic view of a charged particle crossing one spectrometer. The six PT stations are shown in dark grey; the 24 iron slabs (12 per arm) interleaved with 22 RPC planes are shown in light grey. Each spectrometer arm provides an independent measurement of charge/momentum, provided the track is reconstructed in at least one station (or station doublet) in each side of the arm.
walls. Each TT wall consists of a double layer of 256 scintillator strips, for vertical and horizontal coordinate measurements. The TTs can trigger the data acquisition and locate the brick in which the interaction occurred.
The target section is followed by a magnetic spectrometer (see Fig. 1). A large dipolar iron magnet is instrumented with Resistive Plate Chambers (RPCs). The magnetic field is 1.53 T, directed vertically transverse to the neutrino beam axis. The RPC planes are inserted between the iron slabs: they provide the tracking inside the magnet and the range measurement for stopping muons [13].
The deflection of charged particles in the magnet is measured by six stations of vertical drift tubes, the Precision Trackers (PT), grouped in 3 pairs placed upstream of the first, in between and downstream of the second magnet arms (Fig. 1). Each PT station is formed by four staggered layers of aluminum tubes, 8 m long, with 38 mm outer diameter. The spatial resolution is better than 300 µm in the bending (horizontal) plane: this allows the determination of the muon charge sign with high accuracy, and the momentum measurement with a resolution of better than 20% for momenta < 50 GeV/c for charged particles coming from the CNGS direction [14]. The PT system is triggered by the RPC timing boards with a configuration optimised to collect both beam and cosmic ray muons with high efficiency.
In order to remove ambiguities in the reconstruction of particle trajectories, in particular in multi-track events, each spectrometer is instrumented with additional RPC planes (XPC), with two crossed strip planes rotated by ±42.6 • with respect to the horizontal.
Two RPC planes (VETO) are placed in front of the detector, acting as a veto for charged particles originating from the upstream material (mainly muons from neutrino interactions in the rock).
Finally we emphasize that for this analysis the OPERA detector was used differently from what it was conceived for. This is particularly true for the PT system which was configured and optimized to reconstruct and measure particles traveling along the CNGS direction. The results presented are based on data recorded during the CNGS physics run, from June 18th until November 10th 2008. The detector ran in the standard configuration, with the magnetic field directed along the vertical axis in the first arm of both spectrometers, and in the opposite direction in the second arm. Moreover a sample of cosmic ray muons was collected with the magnetic field switched off in order to improve the alignment between PT stations and to evaluate systematic uncertainties. A limited data sample was obtained inverting the magnet polarity to cross check the charge reconstruction.
The data acquisition was segmented in "extraction periods" of about 12 hours each. Only periods of data taking where all the main detector subsystems ran in stable conditions were considered. They amount to about 78% of the total duration of the run. The total number of selected events is 403069 corresponding to 113.4 days of livetime.

Event reconstruction
The OPERA standard software for beam event reconstruction was complemented with a set of dedicated software tools developed for cosmic ray events. Once the event is tagged as "off-beam", that is outside the CNGS spill window it is classified as cosmic and processed in a dedicated way. This choice was required by the different topologies of beam and cosmic ray events. Beam events come from a well defined direction (the CNGS one) and the reconstruction code is optimized to follow a single long track (the muon escaping from the neutrino-interaction region) on the z-axis. Cosmic ray events come from all directions, they are not generated within the target and a fraction of them (∼ 5% in OPERA) are muon bundles. A brief description of the code follows.
The reference frame is defined to have the z axis along the Hall C longitudinal direction (from north to south), y perpendicular to the floor pointing toward the zenith and x describing a right-handed frame. In this coordinate system, the zenith direction θ is defined by the angle with the y axis, the azimuth direction φ by the angle with the z axis. Event reconstruction is performed separately in the two projected views T xz and T yz . First the event direction is determined by using the Hough transform. Using a Monte Carlo simulation we estimated an angular resolution of better than 0.5 • , both in the θ and φ directions, for single as well as for multiple track events. Then, the direction information is used to subdivide the T xz and T yz views in slices 25 cm wide having the same slant as the reconstructed direction. The hits within the same slice are then processed separately to search for a track "seed" of at least three aligned points. If a seed is found, all the other hits in the corresponding projected view are linked to the selected track according to pre-defined tolerances. Tracks independently reconstructed in each view are then merged together to build the three-dimensional event.
For this analysis, particular attention was devoted to the reconstruction of tracks with the Precision Tracker. A description of the PT system is available in [14], while the reconstruction procedures are detailed in [15]. Here we briefly mention the main steps used to extract charge and momentum from PT hits.
A muon crossing the spectrometer is deflected in the horizontal plane ( Fig. 1). Let φ be the angle between the particle direction and the z axis, then the deflection ∆φ is the difference between the two angles measured at both sides of each magnet arm. Each φ value is obtained by fitting the spatial information provided by the PT system which is made of 8 layers of drift tubes arranged in two "stations". Since cosmic ray tracks are almost uniformly distributed in φ, many of them traverse only one single station 1 of a pair and we refer to them as singlets, otherwise we call them doublets.
Due to the different lever arm the angular resolution for singlets is worse than for doublets. For tracks parallel to the z-axis (φ = 0) it is σ φ ≃ 1 mrad for singlets and ≃ 0.15 mrad for doublets. Since for tilted tracks the number of fired tubes and their mutual distances are larger than for tracks with φ = 0, the errors on the slopes decrease. In order to increase the statistics we decided to consider both singlets and doublets: the percentage of cases in which both angles are reconstructed from doublets is ∼55% of the total, ∼9% are from singlets and the remaining 36% are from mixed configurations, namely cases where an angle is reconstructed from a doublet and the other angle from a singlet.
Considering tracks with φ = 0 and the Multiple Coulomb Scattering (MCS) within one magnet arm the total uncertainty on ∆φ is where d = 0.6 m is the magnetized iron thickness, X 0 = 0.0176 m is the iron radiation length and p is expressed in GeV/c. Since the deflection due to the magnetic field is where B = 1.53 T, the requirement ∆φ/σ ∆φ > 1 provides an estimate for the maximum detectable momentum, p max ≃ 1.25 TeV for doublets, p max ≃ 190 GeV for singlets and p max ≃ 260 GeV for mixed configurations.
For muon momenta p ≪ p max the measurement error can be neglected and the only contributions to the ∆φ uncertainty come from the MCS. In this ideal case the ratio ∆φ B /∆φ MCS ∼ 3.5 corresponds to a chargemisidentification η (defined as the fraction of tracks reconstructed with wrong charge sign) below 10 −3 . In reality there are other effects which spoil the resolution and therefore the charge identification capability, as detailed in Sec. 4.1.
To measure the charge and the momentum of a particle at least one ∆φ angle is needed (for tracks parallel to the z-axis there can be up to 4 independent angles). For each reconstructed ∆φ i the track momentum is computed by using the formula [15] where l = 0.82 m is the total arm length (including RPC gaps),B = Bd/l is the effective magnetic field and dE/dz is the ionization energy loss in the magnetized iron (which depends logarithmically on the muon momentum). The term in the square root takes into account the track slope s yz in the T yz view. The muon charge is determined from the sign of the ∆φ i angle, accounting for the particle arrival direction and the field orientation in the arm. The final muon momentum and charge are computed as the weighted average of the independent measurements.

Monte Carlo simulation
Two different Monte Carlo simulations have been used in order to meet two different and conflicting requirements.  . A 3σ cut of the Gaussian fit to Monte Carlo events where secondary particle production was switched off, was applied to the rescaled data (see text).
to correct for detector effects. This task was accomplished with a code which generates muons directly at the detector level fast enough to produce the required statistical amount of events. On the other side, information on cosmic ray primaries and the links between underground and surface variables required the use of a complete generator which simulates the full cosmic ray cascade in the atmosphere.
The first Monte Carlo code (MC1 hereafter) is based on a fast parameterized event generator developed in the framework of the MACRO experiment [16] and adapted for OPERA. The generator considers the MACRO primary cosmic ray composition model [17] and, for each primary mass, it extracts from built-in probability tables the underground muon multiplicity. This choice ensures a self-consistency on the predicted muon flux underground since the composition model and the probability tables were obtained using the same hadronic interaction model. This means that the systematic errors on the primary composition and on the interaction model cancel and the Monte Carlo generator predicts the correct muon flux in the Hall B of LNGS. The event direction (θ, φ) is sampled, together with the muon radial distribution with respect to the shower axis and the muon underground momentum. Finally, the event kinematics is processed within the OPERA standard software chain, from the generation of the relevant physical processes in the electronic detectors up to the reconstruction level. The atmospheric muon charge ratio is introduced by hand (R µ = 1.4). Taking into account the livetime normalization, the ratio between OPERA data and Monte Carlo MC1 predictions is The difference from unity of the ratio is mainly due to the subdetectors efficiency mismatch between data and detector simulation. The second Monte Carlo program (MC2) is based on the more detailed and physics-inspired Monte Carlo program described in Ref. [18], based on the FLUKA code [19]. Here all the main physical processes are implemented, from the primary cosmic ray interactions and shower propagation in the atmosphere up to the muon transport in the overburden. The primary composition model is described in Ref. [20] based on a global fit of several experimental observations. This code is predictive of the muon charge ratio and allows the study of the relations between underground and surface muon energies. Due to its complexity the event production with MC2 is statistically limited.
In the following with the term "Monte Carlo" we refer to MC1 unless otherwise specified.

Data reduction
A set of data quality cuts were made in order to isolate a clean sample of reconstructed muon events. First at least one reconstructed ∆φ angle is required for each event (acceptance cut). Then events with a large number of PT hits potentially dangerous for the muon charge determination are removed (clean PT cut). This typically occurs when some drift tubes are fired by secondary particles (δ-rays, showers etc.) and the best χ 2 track could result from a fake tube configuration. In order to evaluate the maximum number of fired tubes/track allowed by geometrical considerations a special version of MC1 switching off delta ray and secondary particle production was run. By naming M and N the number of fired tubes from Monte Carlo simulation and experimental data, respectively, we derived the functional form M = M (φ), a six-order polynomial shown in Fig. 2a. M (φ) was used to rescale the experimental distribution N as N ′ = N − M (φ) (Fig. 2b).
We considered only tracks with N ′ < 3σ (one-sided cut), where σ is the standard deviation of the Gaussian fit to N ′ . We verified by visual inspection that events rejected by the latter cut are characterized by a large number of additional fired tubes in the neighbourhood of the correct ones.
A further cut was applied on the ∆φ angle (deflection cut). Events having a ∆φ smaller or compatible with the experimental resolution were rejected. On the basis of the plot shown in Fig. 3 where as expected for small deflection values R µ → 1, only events with ∆φ/σ ∆φ > 3 were selected. The effect of this cut is visible in Fig. 4 in which the ∆φ distribution is shown before and after its application. In these plots, experimental data (black points) are plotted with the corresponding Monte Carlo distributions split in the two regions corresponding to positive particles (q true > 0) and negative particles (q true < 0). The chargemisidentification η corresponds to the overlapping region of the two distributions. Averaged over all the event samples η is reduced from 0.080±0.002 to 0.030±0.001 by this cut. The source of events with large ∆φ angles and reconstructed with wrong charge-sign was investigated. A visual scan of Monte Carlo events confirms the hypothesis that they are due to secondary particles in the neighbourhood of the true muon track: if the two tracks are very close, it may happen that the track reconstructed with the best χ 2 is the wrong one. A further selection ∆φ < 100 mrad was used to reject these fake tracks, with a small impact on the statistics. This last selection affects the sample with p µ < ∼ 5 GeV/c.  Fig. 3: Dependence of the measured charge ratio R µ on the deflection angle expressed in units of experimental resolutions. A cut at ∆φ/σ ∆φ > 3 was applied in the data analysis. Note that the fitted value R µ = 1.342 ± 0.015 was obtained with the bins indicated in the plot (the first 3 bins have not been used).
The muon charge ratio was computed separately for single muon events (i.e. event multiplicity n µ = 1) and multiple muon events (n µ > 1). Single muon events are selected by requiring single tracks in each projected view merged in the three dimensional space. Multiple muon events are selected by requiring a muon multiplicity ≥ 2 in both views, with tracks identified and unambiguously merged in 3D space.
Tab. 1 lists the number of events remaining at each stage of the selection process. Note that data and MC1 event rates are absolute (given in day −1 ) and not normalized one to the other. Also note that the clean PT cut has a stronger impact on data reduction and that the effect on the experimental data is different from that of Monte Carlo. This was expected since the percentage of events with PT digits not related to the muon track is intrinsically larger in the experimental data. The clean PT cut was tuned in order to be left with a clean data sample at the expense of a considerable loss of statistics.

Alignment of the PT system
The measurement of the muon charge is strongly affected by the alignment precision of the PT system. Misalignment effects have "global" or "local" contributions. To correct for global effects, which are the dominant ones, each station is treated as an independent rigid body and relative rotations and translations of one station with respect to the others are searched for. The local misalignment contribution takes into account possible distortions or bendings within each station.
A first alignment campaign was carried out with a theodolite to measure the position of the PT walls in the OPERA coordinate system. Recently a more refined alignment using cosmic ray muons was performed. The alignment procedure was carried out in two steps a) PT stations forming a doublet were aligned with the whole data sample, since the space in between has no magnetic field and tracks do not suffer any deflection; b) each doublet (pair) treated as a unit, separated by the iron magnet arm, was aligned using special runs with the magnetic field switched off. This procedure allowed aligning two PT stations within a doublet with a spatial accuracy of ∼0.1 mm and an angular accuracy of ∼0.1 mrad, and to align two doublets with an angular accuracy of 0.2 mrad.
Local effects, such as bendings or distortions, contribute at the second order level and due to the present limited statistics have not been corrected for. However in Sec. 4.2 we provide an estimate of the systematic uncertainty on R µ introduced by these effects.   Table 1: Progressive reduction of the number of events per day after each selection cut, for data (left) and for MC1 (right). The effect of data reduction is also shown by reporting the fraction of events referred to the original sample (f 1 ) and to the previous cut (f 2 ). The total number of experimental events surviving the cuts is 44370. misidentification η and the unfolded charge ratio are reported. The η value, defined as the fraction of tracks reconstructed with wrong charge sign, was extracted from the Monte Carlo simulation. Once η is known, the unfolded charge ratio is obtained according to the formula (see Appendix A): The single muon sample was subdivided into three classes: tracks reconstructed exclusively as doublets, tracks reconstructed exclusively as singlets and as mixed. We verified that the fraction of these classes for experimental data and for Monte Carlo simulation are compatible: 54.8% (doublets), 9.0% (singlets) and 36.2% (mixed) for real data to be compared to 52.5%, 10.0% and 37.5% for Monte Carlo simulation (the errors are < ∼ 0.5%). The final charge ratio value for single muon events, integrated over all the classes, is: The same procedure was applied to multiple muon events. We selected events with n µ > 1 and reconstructed the charge of muons crossing the spectrometer section. Events were classified in this category provided that more than one muon was reconstructed in the detector even though only one charge was measured. In other words, the muon multiplicity is used to "tag" events generated by heavier and more energetic primaries. The possibility to compute the muon charge ratio within the same event is presently excluded by the lack of statistics of high multiplicity muon bundles. The charge ratio is R meas µ (n µ > 1) = 919/753 = 1.22 ± 0.06 and the corresponding unfolded value, obtained from Eq. 10 R unf µ (n µ > 1) = 1.23 ± 0.06.
This value is 2.4σ away from the value for single muon events, consistent with the hypothesis of dilution of R µ due to the neutron enhancement in the primary nuclei. Tab. 3 gives information obtained with MC2 on some variables of single muon events and muon bundles in the OPERA detector. In particular, are given the average pri-mary mass number A , the average primary energy/nucleon E/A , the fraction of Hydrogen nuclei over the total (H fraction), the ratio of protons over neutrons in the primary radiation N p /N n and finally the measured muon charge ratio R unf µ .

Systematic uncertainty on R µ
The main sources of systematic uncertainties in the determination of R unf µ are related to the alignment accuracy of the PT system and to the determination of the η value.
The systematic uncertainty due to misalignment effects was evaluated in different ways. A given offset ∆φ → ∆φ+δφ can be directly propagated in the algorithm which computes the charge ratio to evaluate R µ → R µ + δR µ . The δφ = 0.2 mrad uncertainty on the alignment accuracy obtained with magnets off (Sec. 3.5) corresponds to δR µ = 0.03. However a more powerful procedure was used to better estimate this systematics. We considered all muon tracks crossing both arms of each spectrometer, thus providing two independent deflection values ∆φ per spectrometer for the same muon track. With perfect alignment and neglecting the energy loss the difference δ∆φ = ∆φ arm1 − ∆φ arm2 should be peaked at zero. The two distributions, one for each spectrometer, are shown in Fig. 5 together with a Gaussian fit to the central part of the distributions, where the effects of muon energy loss in the magnet iron are negligible. The two peaks are at 0.08 mrad and -0.07 mrad respectively, ∼2 standard deviations away from zero. A misalignment of 0.08 mrad produces an error on the charge ratio δR µ ≃ 0.015. We quote this number as the limiting alignment accuracy of each doublet with respect to the other. This number is conservative since it assumes that all four arms are affected independently from the same uncertainty. In reality only the the outer two doublets of each magnet contribute to this error, since a given offset in the central doublet cancels the systematic uncertainty for ∆φ arm1 and ∆φ arm2 .
A further test which also incorporates local effects consists in comparing the values R i µ (i=1,...,4) in each magnet arm. The average difference from the mean value i |R i µ − R µ |/4 = 0.017 is within the statistical accuracy of each δR i µ = 0.03. Another consistency check exploited a small data sample (∼9 days of livetime) obtained after inverting the polarity of the magnetic field. Running with inverted magnetic polarity could in principle cancel the systematic error related to misalignment effects. The result is R inverted µ = 1.36 ± 0.04, corresponding to the unfolded value R inverted µ = 1.39 ± 0.04. Even if the statistical error is larger than the systematic error quoted above, the result is in good agreement with the value obtained with normal polarity.
The charge-misidentification η was previously estimated using Monte Carlo simulations. As already discussed the value is larger than what is expected from multiple scattering alone. The difference is ascribed to the inclusion of spurious effects, such as the production of secondary particles near the muon trajectory, timing errors, and other second order effects not reproducible with the Monte Carlo program. Therefore we expect that the systematic uncertainty on η is one-sided, being η real ≥ η MC . To estimate this difference η was evaluated using experimental data for a subsample of events. We considered all muon tracks crossing both arms of each spectrometer, which provide two independent deflections ∆φ of the same muon track. In this case, the probability that the two deflection angles have opposite sign is p = 2η(1 − η) and therefore η = 1 − √ 1 − 2p. This formula neglects the correlation between the two ∆φ angles, since they are built using a common track (the one in between the two arms). The correct η(p) relation was derived using a Monte Carlo simulation applied to the experimental and simulated data. It was found for the case of doublets η data = 0.018 ± 0.002 and η MC = 0.012±0.002. Considering doublets and mixed configuration together, we found η data = 0.026±0.002 and η MC = 0.019 ± 0.002. The difference δη = 0.007 was used as an estimate of the systematic uncertainty on η which corresponds to δR µ = 0.007.
The final systematic error is taken as the quadratic sum of its contributions and it is assumed be the same for single and for multiple muon events: 5 R µ as a function of p µ and E µ cos θ * The underground muon momentum p µ was computed using Eq. 8. The muon charge ratio as a function of p µ is shown in Fig. 6, where the widths of the horizontal error bars correspond approximately to the (average) muon momentum resolution. A linear fit gives a 0 = 1.29 ± 0.06 and a 1 = 0.05 ± 0.03 with χ 2 /dof = 13.7/15. The data are also compatible with the hypothesis of a constant charge ratio, since the fit to a constant yields a 0 = 1.379 ± 0.015 with χ 2 /dof = 16.2/16 and therefore ∆χ 2 /dof = 2.47/1 (corresponding to ∼1.6 sigma). The muon energy at the surface (E µ ) is directly related to the underground residual energy (E µ ≃ p µ ) and to the rock amount crossed by the muon to reach the detector level. In fact, the energy loss of high energy muons in the rock is usually expressed as where h is the rock depth while the two energy-dependent parameters α and β are the contributions of the ionization energy loss and the radiative processes, respectively. Eq. 15 can be integrated to obtain the approximate formula which connects the surface and underground muon energies. However, Eq. 16 is valid only on average. The "resolution" dE µ = E rec µ − E true µ is dominated by the statistical   Table 3: Primary cosmic ray information for single and multiple muon events (see text). Reported numbers were obtained with MC2 and with the composition model fitted in [20]. Only statistical errors are quoted. Systematic uncertainties related to the composition model dominate and can be inferred from the cited reference (δ A ≃ 1). In the last column the measured (and unfolded) charge ratios are given.  Table 4: Main information for the five bins in E µ cos θ * . From left to right: the energy range and average value, the number of muons reconstructed with positive and negative charges, the unfolded charge ratio, the statistical and systematic errors.
fluctuations due to the discrete processes described by the term β [21]. We evaluated E µ with a full Monte Carlo simulation to build the table E µ = f (h, p µ ). For this purpose the code MC2 was used since it contains a detailed description of the muon flux at the surface and the muon transport in the Gran Sasso rock. The (h, p µ ) plane was divided into 10×10 equally-spaced bins and in each bin the average E µ value was computed. The binning was chosen coarse enough to have a large statistical sample in each bin without affecting the resolution dE µ , which is of the order of 0.15 in the logarithmic scale. The surface muon charge ratio was computed as a function of the variable E µ cos θ * binned according to the resolution. Finally, the experimental values were corrected in each bin for the corresponding charge mis-identification and shown in Fig. 7 with black points (only single muon events were considered). The present statistics does not allow to draw conclusions about the highest energy data point shown in the figure.
Tab. 4 gives some information for each of the five bins considered: the energy range and average value, the statistical sample, the unfolded charge ratio, the statistical and systematic errors. The latter were evaluated computing in each bin the two contributions discussed in Sec. 4.2.
In Fig. 7 are shown for comparison the data from other experiments for which we could recover information on the E µ cos θ * variable. For the low energy region we took data from Ref. [22] and Ref. [23] (we choose data points with uncertainties δR µ <0.05) while in the high energy region the data are from Ref. [9] and Ref. [10]. For the latter, since the angular information were not provided in the paper, we plotted the R µ integrated value in correspondence of the E µ cos θ * value given in Ref. [8]. We also report a recent result from Ref. [24] where the vertical muon charge ratio is given in the range 1-3 TeV (average value 1.3 TeV).
Finally, we fit our data to Eq. 5, using a procedure similar to what is described in Ref. [10]. We rewrite Eq. 5 in the form: where R Kπ = Z N K /Z N π and f π + = 1−f π − = Z N π + /Z N π (and similarly for kaons). f π + and f K + were left free to vary while we fixed the kinematical parameters a π = 0.674, a K = 0.246, b π = 1.061, b K = 1.126 and the fraction of kaons over pions in the atmosphere R Kπ = 0.149 [5]. The fit of R µ = φ µ + /φ µ − takes into account data from [22] and [23] for the low energy region and data from this work at higher energies. The fit yields the values f π + = 0.5514±0.0014 and f K + = 0.680±0.015 which correspond to a ratio R π = Z N π + /Z N π − = 1.229±0.001 and R K = Z N K + /Z N K − = 2.12±0.03 for pions and kaons respectively. The result of the fit is shown in Fig. 7 as a continuous line. The contribution of the prompt muon component to R µ was evaluated for three different charm production models: the phenomenological non-perturbative models RQPM and QGSM [2] and the semi-empirical model from Volkova et al. [25]. In [2] the prompt muon flux and charge  ratios are parametrized as a function of the muon energy. The results of the fit extended to include the prompt contribution as predicted by these models, are shown in Fig.  7. The pion and kaon charge ratios obtained from the fit are unchanged within the statistical errors.

Conclusions
The atmospheric muon charge ratio R µ = N µ + /N µ − was measured using the spectrometers of the OPERA underground detector. We analyzed four months of data taken to be compared to R unf µ (n µ > 1) = 1.23 ± 0.06 for muon bundles. This difference of about ∼2.4σ supports the hypothesis of the decrease of the muon charge ratio with increasing primary mass. This is the first indication of such an effect which provides a further handle for the correct understanding and modelling of the secondary production in the atmosphere.
The underground muon charge ratio is consistent with past measurements in a similar energy region. Data suggest a slight increase of R µ with the underground muon momentum, although a fit to a constant charge ratio cannot be excluded.
The dependence of R µ on the vertical surface energy E µ cos θ * shows an increase in the region 1-3 TeV and it is compatible with a model which considers only the π and K contributions to the muon charge ratio. A fit of the low energy data and our data with a simplified description of the atmospheric muon flux provides a value of the pion and kaon charge ratios R π = Z N π + /Z N π − = 1.229±0.001 and R K = Z N K + /Z N K − = 2.12±0.03, respectively. The inclusion of the prompt muon component does not modify the fit results. It is however intriguing to observe that our measurement lies in the region where the charmed particle production may start to give an observable contribution to the muon charge ratio. A larger statistical sample or an experimental measurement with a new detector at very large depths could shed light on the region E µ cos θ * > ∼ 10 TeV. The data collected by OPERA at the end of its scientific program will allow to improve the measurement in this energy region.  Fig. 7: R µ values measured by OPERA in bins of E µ cos θ * (black points). Also plotted are the data in the low energy region from MINOS-ND [22] and L3+C [23] and in the high energy region from Utah [9], MINOS [10] and LVD [24] experiments. The result of the fit of OPERA and L3+C data to Eq. 17 is shown by the continuous line. The dashed, dotted and dash-dot lines are, respectively, the fit results with the inclusion of the RQPM, QGSM [2] and VFGS [25] models for prompt muon production in the atmosphere.