Measurement of the muon anti-neutrino double-differential cross section for quasi-elastic scattering on hydrocarbon at~$E_\nu \sim 3.5$ GeV

We present double-differential measurements of anti-neutrino quasi-elastic scattering in the MINERvA detector. This study improves on a previous single differential measurement by using updated reconstruction algorithms and interaction models, and provides a complete description of observed muon kinematics in the form of a double-differential cross section with respect to muon transverse and longitudinal momentum. We include in our signal definition zero-meson final states arising from multi-nucleon interactions and from resonant pion production followed by pion absorption in the primary nucleus. We find that model agreement is considerably improved by a model tuned to MINERvA inclusive neutrino scattering data that incorporates nuclear effects such as weak nuclear screening and two-particle, two-hole enhancements.


I. INTRODUCTION
Although quasielastic neutrino interactions are a key signal process for accelerator-based oscillation experiments, models of these interactions on nuclei have large (∼30%) uncertainties [1,2]. These arise from several sources, including nucleon form factors and final-state interactions wherein the produced particles interact further before exiting the primary nucleus. Final-state interactions can also cause other processes such as resonant pion production to have a zero-meson final state that will appear with a quasielastic topology in a detector. Interactions with multinucleon states can similarly produce zero-meson final states.
These and similar sources of uncertainty on other processes dominate the systematic uncertainty budgets of the current oscillation measurements such as T2K [3] and Nova [4] and will limit the reach of oscillation experiments such as DUNE [5] if not further reduced. Because any one measurement of quasielastic scattering necessarily measures a superposition of these effects, lowering model uncertainties on individual parameters will require many different measurements to untangle the many unknowns.
In this article, we present a critical ingredient in this process: a double-differential measurement of the antineutrino quasielastic (QE) cross section as a function of the transverse and longitudinal momentum of the final-state muon. We include in our measurement events consistent with zero-meson final states arising from resonant pion production followed by pion absorption in the nucleus and from interactions on multinucleon states (frequently referred to as two-particle, two-hole or 2p2h). This ensemble of signal processes, which is defined precisely in Sec. VI, is referred to hereafter as "QE-like." In addition to this primary result, we also present a number of auxiliary measurements including double-differential cross sections as a function of alternate variables, single-differential projections, and comparisons of reconstructed energy near the event vertex to various models. The neutrino energy range of 1.5-15 GeV covered by this measurement is well matched to that of present and future neutrino oscillation experiments with baselines on the 1000 km scale, including MINOS [6], NOvA [7], and DUNE [8].
The measurement described here extends a previous measurement of antineutrino quasielastic scattering by MINERvA [9] and is a companion to similar studies of neutrino scattering [10]. The measurement complements other MINERvA QE-like studies that look in detail at the hadronic component of the final state [11] and that study neutrino interaction cross sections as a function of nuclear mass [12,13].
This article is organized as follows: Section II reviews the current status of neutrino-nucleus QE-like scattering models. Sections III and IV review the MINERvA experiment and simulation. Event reconstruction and selection are discussed in Secs. V and VI. The cross-section extraction procedure and systematic uncertainties are detailed in Secs. VII and VIII. The results and comparisons with models are presented in Sec. IX, and the article is summarized in Sec. X. The Appendixes present additional results, including cross sections with alternate signal definitions; those planning to use the data are encouraged to use the Supplemental Material [14], which provide higher numerical precision.

II. THEORY OF QE-LIKE INTERACTIONS
In charged-current quasielastic (CCQE) scattering on free nucleons, an incoming muon antineutrino interacts with a target proton, exchanging a W AE boson to knock out a neutron and leave a positively charged muon in the final state:ν μ þ p → μ þ þ n: In this case, it is possible to reconstruct certain characteristics of the interaction using only the kinematics of the outgoing charged lepton assuming the initial-state nucleon is at rest. For a nucleon bound within a nucleus, the incoming neutrino energy and the four-momentum transfer Q 2 can be estimated, respectively, as where E b is the initial-state nucleon's binding energy, taken to be 30 MeV, as described in [9], E ν and E μ are the neutrino and muon energy, respectively, p μ and θ μ are the muon's momentum and angle with respect to the neutrino, respectively, and m μ , m n , and m p are the mass of the muon, neutron and proton, respectively. The QE subscript and superscript here and throughout the remainder of this article denote quantities computed under an assumption of a quasielastic hypothesis with the initialstate nucleon at rest.
In the case of a bound nucleon, Fermi motion and nucleon correlations mean that the initial-state nucleon is not at rest, making the QE kinematic variables only estimates of the true values. The final-state interpretation can also be affected, as an ejected nucleon may interact with other nucleons while escaping the nucleus. Other interactions such as resonant pion production can be modified by final-state nuclear effects to have no pions in the final state, thus appearing QE-like. Similarly, interactions with correlated pairs of nucleons can also produce final states that appear quasielastic. All of these nuclear effects can cause quasielastic neutrino interactions on heavy nuclei to differ substantially from those on free nucleons. In this section, we discuss the quasielastic scattering from a free nucleon and several contemporary theories that attempt to model the impact of the nuclear environment. More detail can be found in [15].

A. Quasielastic antineutrino scattering on free nucleons
Because the internal structure of the initial-and finalstate nucleons is governed by the nonperturbative regime of QCD, it is not possible to make a precise ab initio calculation of the neutrino-nucleon quasielastic cross section; it may instead be described by nucleon form factors. In the 1972 review article of Llewellyn-Smith [16], the differential quasielastic cross section is expressed as a function of two vector, one axial-vector, one pseudoscalar and two second-order form factors. All but the axial form factor are known from electron-nucleon scattering measurements. The axial form factor must be taken from neutrino scattering or pion electroproduction measurements and is typically parametrized as a dipole: The value of the axial form factor at Q 2 ¼ 0 has been measured through beta-decay experiments [17,18], leaving one free parameter, the axial mass M A . Deuterium and hydrogen bubble chambers [19,20] have measured the value of M A on free or quasifree nucleons. An average value of M A ¼ 1.014 AE 0.014 GeV/c 2 was extracted by Bodek et al. [21] in 2008. Modern experiments on heavy nuclei have favored higher values of M A [22][23][24], and the discrepancy between these and the deuterium experiments has been attributed to insufficiencies in the nuclear models used to extract the axial mass on heavy nuclei. Alternate parameterizations of the dipole form factor are also available. In particular, a more general "z-expansion" parameterization [25] has been widely adopted in flavor physics and was recently implemented in neutrino event generators.

B. Scattering from nuclei
When simulating quasielastic scattering in heavy nuclei, the most commonly used nuclear model is the relativistic Fermi gas (RFG) model proposed by Smith and Moniz [26], in which scattering from a nucleon in a nucleus is treated as if the incoming lepton scatters from an independent nucleon (the "impulse approximation"). However, in the case of the RFG, the target nucleon is not stationary but has a momentum consistent with the Fermi distribution. Thus the cross section for scattering off the nucleus is replaced by an incoherent sum of cross sections for scattering off of individual nucleons, with the remaining nucleus (depleted by one nucleon) as a spectator.
The local Fermi gas (LFG) model is a extension to the RFG model in which a local density approximation [27,28] is used, so that instead of using a constant average field for the whole nucleus, the momentum distribution is dependent on a nucleon's position within the nucleus. This gives a Fermi motion distribution that is not sharply peaked at the momentum limit and is both more natural and reproduces the measured peak of the distribution.
Spectral functions can be also used to improve the relativistic Fermi gas model [29]. The Hamiltonian for a large nucleus is so complicated that it is impractical to try to solve the many-body Schrödinger equation for the entire nucleus. However, if a mean field is used to replace the sum of the individual interactions, a spectral function can be constructed that represents the probability of finding a nucleon with momentum and removal energy within the nucleus. There are a number of non-Fermi-gas approaches [30][31][32][33][34][35].

C. Multinucleon correlations
The models described above, which treat individual nucleons independently, do not fully take into account the nature of the nuclear force, which has a short range with a repulsive core [36]. Interactions between two (or more) spatially close nucleons can give the individual nucleons very high momenta, far above the Fermi momentum. Electron scattering [37] data indicate that approximately 20% of the nucleons in carbon atoms are part of correlated pairs. These experiments have observed ejected nucleons consistent with knocked-out partner nucleons [38,39], with 90% of those pairs found to be proton-neutron pairs, with the remainder being nn and pp pairs. In the case of np pairs, it is expected that charged current QE-like antineutrino scattering on protons within a correlated pair would tend to produce two neutrons-the expected neutron produced by the QE-like interaction, plus the knockedout neutron partner.
The impact of multinucleon states in the initial state has been modeled in a number of ways. Bodek and Ritchie's modification to the relativistic Fermi gas model [40] adds a high-momentum tail to the RFG's momentum distribution, based on the nucleon-nucleon correlation function and fits to data as explained in [41]. While this method attempts to provide a realistic initial momentum distribution, it does not include any model for the ejection of paired nucleons.
Going a step further, Bodek, Budd and Christy [42] have developed a "transverse enhancement model" (TEM). They fit inclusive electron scattering data [43,44], modifying the nucleon magnetic form factors to accommodate the enhancement of the transverse response observed in those data. The resulting form factor modifications plus an unmodified axial vector form factor were then used to predict neutrino-nucleus scattering cross sections, producing results that were consistent with both low-energy data from MiniBooNE [22] and high-energy results from NOMAD [45]. The TEM fit does not attempt to model additional knocked-out nucleons in either the electron or neutrino case. Those empirical fits, when combined with different expressions [46] for how the structure functions should be related, and attaching a two-nucleon knockout are used in the GiBUU generator.
While the approaches described above are empirical models that ascribe effects observed in electron scattering to multinucleon processes and then attempt to predict the effect they might have on neutrino scattering, 2p2h models attempt to predict multinucleon effects in neutrino scattering from first principles. These models consider pairs of nucleons connected by the exchange of virtual pions and rho mesons [47]. There has been a recent dramatic expansion of work on 2p2h models, such as those of Marteau and Martini [48], IFIC Valencia group [49], the SuSA group [34,35], and the Gent group [50]. As of this writing, the IFIC Valencia model is implemented in GENIE, NuWro, and NEUT for the CC process. Empirical versions related to the TEM fits to electron scattering data are also available in GENIE (without twonucleon knockout) and GiBUU (with two-nucleon knockout). Finally, there is a version that implements a simple empirical shape in W and Q2 developed for use in electron scattering codes [51] to generate events as an option in GENIE [52].
Long-range correlations between nucleons are typically modeled using an approach known as the random-phase approximation (RPA) [53]. It is based on the phenomenon, observed in β-decay and muon capture experiments, that the electroweak coupling can be modified by the presence of strongly interacting nucleons in the nucleus, when compared to its free-nucleon coupling strength, similar to the screening of an electric charge in a dielectric. The RPA approach affects cross-section predictions at low energy transfers (and low Q 2 ), where a quenching of the axial current reduces the cross section compared to the RFG prediction. It also introduces a small cross-section enhancement at intermediate Q 2 . Multiple RPA models are available within generators, including those of Nieves [54], Martini [48], Graczyk and Sobczyk [55], and Singh [56]. There is also discussion of the interplay between RPA with the Fermi gas and beyond-the-Fermi-gas models [57,58].

Final-state interactions
Nonquasielastic processes that undergo final-state interactions within the nucleus can have QE-like final states. For example, there are three possible antineutrino chargedcurrent resonant pion production processes: In such events, there is an ∼20% possibility that the pion will be absorbed before it exits the nucleus, leaving a QE-like final state of a single muon and one or more nucleons. Nearly all available models of final-state interactions are intranuclear cascade (INC) models in which final-state hadrons are individually propagated through the nucleus with some probability of undergoing interactions such as absorption or inelastic scattering with the nuclear medium. The details of the interactions vary significantly across different models (typically implemented as part of an event generator such as GENIE [1] or NEUT [59]), but all are tuned to hadron scattering data. Some generators (including GENIE) also provide effective cascade models wherein the cascade of interactions that particles may undergo as they traverse a nucleus is modeled as a single interaction. At least one alternative to INC models exists in the form of a semiclassical nuclear transport model implemented as part of GiBUU [60].

III. MINERvA EXPERIMENT
The MINERvA experiment is situated in the NuMI (Neutrinos at the Main Injector) neutrino beam at Fermilab. The detector and beam are described in detail in [61,62]; this section summarizes their main features, focusing on the components relevant to this study.

A. The NuMI neutrino beam
Fermilab's NuMI beam uses 120 GeV/c protons from the main injector, which impinge on a 1-m-long graphite target. The resulting pions and kaons are focused by a pair of movable parabolic horns. The horn current polarity can be set to focus positively or negatively charged mesons, which decay in a 675-m-long decay pipe, producing muons and neutrinos. An absorber removes any remaining hadrons from the beam and 200 m of rock filter out muons, leaving a beam of primarily neutrinos or antineutrinos, depending on the horn current polarity.
For the low-energy beam configuration used in this work, the peak beam energy was approximately 3 GeV and the horns were configured to focus negative particles. We use data recorded between November 5, 2010 and February 24, 2011, corresponding to 1.020 × 10 20 protons on target (POT). The Monte Carlo simulation described in the next section corresponds to 9.247 × 10 20 POT.

B. MINERvA detector
The MINERvA detector is composed of an inner detector (ID) and an outer detector. The most upstream portion of the ID consists of active scintillator planes interspersed with passive nuclear targets. This region is used for studies of the A dependence of neutrino interaction cross sections but is not used in the work described here. Immediately downstream of the nuclear target region is the central tracker, followed by electromagnetic and hadronic calorimeters (ECAL and HCAL, respectively).
The tracker is composed of 124 active scintillator planes, each consisting of 127 strips of doped polystyrene scintillator with a titanium dioxide coating. The strips have a triangular cross section 17 mm (height) by 33 mm (base) and vary in length between 122 and 245 cm depending on the position within the plane. Each scintillator plane is installed in one of three orientations, X, U or V. In the X orientation, the strips are vertical. Strips in the U or V planes are oriented at AE60°with respect to the strips in the X planes. Each module in the active tracker region consists of two planes of scintillator strips, alternating between UX and VX configurations. A 2-mm-thick lead collar covers the outermost 15 cm of each plane, forming the side electromagnetic calorimeter.
The downstream ECAL consists of ten modules that are similar to tracker modules, except that the 2-mm-thick lead collar is replaced by a 2-mm-thick sheet of lead covering the plane. The 20 HCAL modules, downstream of the ECAL, each contain only one plane of scintillator, followed by a 2.54-cm-thick plane of steel.
Light produced in the scintillator is collected by a 1.2 mm diameter wavelength-shifting optical fiber inserted in a hole passing along the length of the strip and transmitted by the optical fibers to Hamamatsu H8804MOD-2 photomultiplier tubes (PMTs), as described in [61]. The full detector has 507 PMTs, each of which consists of 64 pixels. The PMTs are read out via a data acquisition system that is described in detail in [63]. Raw PMT counts are transformed into estimated energy deposited in the strip via the calibration chain also described in [63]. The intrinsic time resolution is 4 ns.

C. The MINOS Near Detector
The MINOS Near Detector [6] is located 2 m downstream of MINERvA and is used to measure the charge and momentum of muons exiting the back of MINERvA. The 1 kT MINOS detector is composed of 2.54-cmthick steel planes, interspersed with 1-cm-thick layers of scintillator. The scintillator planes are formed from 4.1cm-wide scintillator strips, with orientation of the strips alternating between þ45°and −45°to the vertical in successive planes. The first 120 planes are instrumented for fine sampling; in this region, every fifth steel plane is followed by a fully instrumented scintillator plane, while all other steel planes are followed by a partially instrumented scintillator plane. The coarse-sampling region, further downstream, has only the fully instrumented scintillator every five planes; there are no partial scintillator planes in this region. The MINOS detector is magnetized by a coil that runs in a loop passing through the detector, generating a toroidal field with an average strength of 1.3 T.

A. Beam flux simulation
MINERvA's simulation chain begins with G4Numi [64], a GEANT4 [65] based simulation of the NuMI beam line from primary proton beam to the MINERvA detector. The FTFP_BERT inelastic scattering model of GEANT version 4.9.2.p03 is used. This raw simulation is found to disagree with existing hadroproduction data from the NA49 [66] and other experiments [67,68] and is therefore corrected so that both differential and total interaction cross sections in the simulation match these external data sets. Version 1 of the PPFX package is used to implement these corrections [69].
We also use neutrino-electron scattering data collected in the MINERvA detector with the beam line in neutrino mode (focusing positive pions) as an independent constraint on the flux model, as described in [70]. This constraint lowers the predicted neutrino flux by 2%-4% depending on the neutrino energy. While an equivalent measurement is not available for the antineutrino running mode due to low statistics for theν − e process in that configuration, the known correlations between the neutrino and antineutrino fluxes are used to apply this constraint to the antineutrino flux distribution. As shown in Fig. 1, applying the constraint results in a 1%-3% decrease in the antineutrino flux prediction.

B. Neutrino event generation
MINERvA uses a modified version of the GENIE neutrino interaction event generator [71] version 2.8.4 to model physics processes within the primary interaction nucleus. Simulated event distributions using this generator, with data constraints described in Sec. VII, are used to estimate background levels, resolution effects, acceptance and efficiency.
GENIE models the nucleus using the relativistic Fermi gas model [26] incorporating the Bodek-Ritchie highmomentum tail [40] that simulates short-range correlations. For carbon, the maximum momentum for Fermi motion is taken as k F ¼ 0.221 GeV/c, and Pauli blocking is also included. Quasielastic cross sections follow Llewellyn-Smith's prescription. Vector form factors are modeled by default using the BBBA05 model [72]. For the axial vector form factor f A , a dipole form is used, with f A ð0Þ ¼ 1.2670 and axial mass M A ¼ 0.99 GeV/c 2 [73].
GENIE uses the Rein-Sehgal model [74] to simulate baryon resonance production, which provides cross sections for 16 different resonance states. The resonant axial mass M RES A is taken to be 1.12 GeV/c 2 . Deep-inelastic scattering (DIS) cross sections are calculated with an effective leading order model with a low-Q 2 modification from Bodek and Yang [75]. Hadronic showering is modeled with the Andreopoulos-Gallagher-Kehayias-Yang (AGKY) model [76]. The Bodek-Yang model also describes other low-energy nonresonant pion production processes. Rescattering of nucleons and pions in the nucleus is simulated using the INTRANUKE-hA intranucleon hadron cascade package [77]. While the resonant interactions described earlier account for the majority of pion production, other inelastic processes, as described by Bodek and Yang [75] are also possible. In particular, nonresonant pion production followed by finalstate interaction (FSI) can produce a QE-like signature.
In addition to the basic processes simulated in GENIE 2.8.4 we also apply three additional corrections. First, we reweight quasielastic events as a function of the energy and 3-momentum transfers q 0 and q 3 to include the random phase approximation model as predicted by the Valencia model of Nieves, Amaro, and Valverde [54] and implemented for MINERvA [78]. Figure 2 shows the Q 2 dependence of this correction. Second, QE-like interactions on multinucleon pairs are simulated using the Valencia IFIC model. We modify this model to match MINERvA inclusive neutrino scattering data reported in [79], which enhances this contribution by approximately 60%. 1 Finally, the normalization of nonresonant pion production is reduced to 43% of the default GENIE 2.8.4 prediction, based on a fit to pion-production data on deuterium from bubble-chamber experiments at Argonne and Brookhaven National Laboratories [80]. We reduce the uncertainty on the normalization of this process to 5%, based on the same data fit. This modified version of GENIE 2.8.4 is hereafter referred to as MINERvA -tuned GENIE.

C. Detector simulation
The GEANT4 toolkit [81] v4.9.4p02 with the QGSP_BERT physics list is used to simulate propagation through the material of the detector and has been validated using a scaled-down version of the MINERvA detector in a test beam [82]. The optical and electronics systems are also simulated, which allows the energy depositions recorded by GEANT4 to be converted to a simulated readout that can be analyzed as if it were MINERvA data. These simulated data are overlaid with actual data to include the effects of multiple neutrino interactions, noise and dead time, which  The IFIC model has a cutoff of q 3 < 1.2 GeV. While this is marginally within the MINERvA kinematic range, the bulk of the distribution is near 0.5 AE 0.2 GeV and the contribution above 1.2 GeV would be negligible. is a result of the ∼150 ns digitization window following activity above threshold, during which additional deposits will not be recorded.

V. EVENT RECONSTRUCTION
Calibrated energy depositions in the scintillator strips (referred to subsequently as "hits") are reconstructed into antineutrino interaction candidates through a series of steps. First, the ensemble of hits collected over the 10-μs-long NuMI beam spill are grouped into time slices corresponding to individual neutrino interactions. Hits within the same time slice are then collected into clusters that are adjacent in strip space and contained within the same scintillator plane. The position of the cluster is taken to be the energy-weighted average of the hit (strip) positions; the cluster time is set to the time of the highestenergy hit.

A. Track reconstruction
Track reconstruction begins by collecting clusters within a single time slice into "seeds" containing three clusters in consecutive planes of the same (X, U or V) orientation that fit to a straight line. Seeds are merged into track candidates within each view (X, U and V), and candidates are formed into three-dimensional tracks, which are fitted with a Kalman filter routine [83,84], in combination with additional untracked clusters in planes adjacent to the track. This allows tracks to be extrapolated through areas of high activity (such as a hadron shower). This algorithm is then repeated to identify additional charged particles in the event until no further tracks are identified.
A similar reconstruction algorithm is performed in parallel in the MINOS detector, where time slices are selected by looking at hits clustered in space and time. The hits in a given time slice are then formed into clusters, which are grouped into tracks if their positions are correlated. Each track's path is then estimated using a Kalman filter; unlike in MINERvA, MINOS tracks curve due to the detector's magnetic field. For tracks stopping within the detector and not entering the coil, the track's momentum is estimated via range; otherwise, the momentum is estimated via curvature through the Kalman fit. For the data considered here, MINOS's magnet was configured to focus positive muons.
Once tracks have been formed in both MINERvA and MINOS, they are then matched between the two detectors. MINOS tracks are matched to MINERvA muons when activity is measured in the last five planes of MINERvA, and a track starts in the first four planes of MINOS within 200 ns of the MINERvA track time. The MINERvA track is extrapolated forwards to the first MINOS plane, and the MINOS track is extrapolated back to the last plane of MINERvA. The point of closest approach between the two tracks is required to be less than 40 cm. The final step of track reconstruction is known as muon "cleaning." MINOS-MINERvA matched tracks are deemed to be muons. Energy beyond the expected deposition of a minimum-ionizing particle is removed from the muon track and added to the ensemble of unmatched clusters considered for further reconstruction.

B. Recoil energy reconstruction
We refer to final-state energy not associated with the muon track as "recoil energy." In this study we consider only energy deposited in the tracker and ECAL portions of the detector and further require recoil cluster times to be between 20 ns before and 35 ns after the path-lengthcorrected average time of clusters on the muon track. We also exclude all clusters likely to be due to PMT cross talk and clusters within 10 cm of the muon vertex from the recoil energy sum, to minimize dependence on simulations of energy near the vertex, which are sensitive to details of final-state and multinucleon interactions. Energy in all remaining clusters is summed and calorimetrically corrected: where C sd i is a calorimetric constant obtained from the simulation for subdetector i that corrects for the passive material fraction in that subdetector (1.22 for the tracker and 2.013 for the ECAL).

VI. EVENT SELECTION
Before identifying selection criteria for isolating signal events, it is necessary to clearly define what is meant by "signal." For MINERvA's first studies of quasielastic scattering [9,10], we attempted to measure events in which the underlying neutrino-nucleon interaction was quasielastic, regardless of how those events were modified by final-state interactions. Several other experiments have recently published measurements [22,24,85] of QElike events with a final state of an appropriately charged muon, plus nucleons. In this case, resonant pion production events where the pion is absorbed become part of the signal to be measured. However, in MINERvA's scintillator tracker, which is able to resolve proton tracks above a kinetic energy of 120 MeV and to detect the energy of lower-energy particles, this definition is not ideal. For this study, we define our signal to be events that are antineutrino charged-current events occurring in the MINERvA tracker fiducial volume, have post-FSI final states without mesons, prompt photons above nuclear deexcitation energies, heavy baryons, or protons above our proton tracking kinetic energy threshold of 120 MeV, and include a muon emitted at an angle with respect to the beam of less than 20°, 1.5 GeV < p k < 15 GeV and p T < 1.5 GeV (matching the region where tracks can be reconstructed in both MINERvA and MINOS with wellreconstructed momentum). This is similar to the QE-like (often called CC0Pi) definitions used by other experiments [22,86], modified slightly to correspond to MINERvA's acceptance, which is poor for events with high angle muons, very low or very high momentum muons, but able to reject high momentum protons. We also report alternate results where the signal definition consists of interactions that were initially generated in GENIE as quasielastic (that is, no resonant or deep inelastic scatters but including scatters from nucleons in correlated pairs with zero-meson final states), regardless of the final-state particles produced.
We begin the event selection by identifying time slices containing at least one track reconstructed in the MINERvA detector and matched to a track in the MINOS detector as described in Sec. V. This provides a high-purity sample of charged-current events. To isolate antineutrino event candidates, we further require that the charge-momentum ratio (q/p) returned by the MINOS Kalman fit be positive. Because we also require no visible proton in the final state, the remaining neutrino contamination in our samples is quite low-0.6% in the simulation-and is accounted for in the acceptance calculation. Because MINERvA experiences some dead time after an event has been recorded, we further require that no more than one strip immediately upstream of the track vertex (projected along the track direction) or immediately adjacent to these strips be dead at the time of the neutrino event. This reduces candidates arising from through-going muons generated upstream of the detector to less than 0.1%. We require the reconstructed interaction vertex to be within the fiducial volume of our detector; the vertex must be within a hexagon of apothem 850 mm and fall within modules 27-80, inclusive, corresponding to 108 tracking planes. We also require our reconstructed muon longitudinal momentum to be less than 15 GeV. This removes very energetic, forward-going muons that have poor energy reconstruction in MINOS.
To reduce backgrounds from non-QE-like events, we require that no tracks other than the muon track be reconstructed between 20 ns before and 35 ns after the muon track (the same time window used for recoil energy reconstruction). This reduces backgrounds from events with charged pions, particularly at high Q 2 QE where the recoil cut described below is very loose, while the narrow time window minimizes the likelihood that signal events are rejected due to overlapping neutrino interactions.
Charged pions and high-energy protons do not always leave reconstructable tracks; they do, however, deposit clusters of energy in the detector. We therefore consider recoil energy, reconstructed as described in Sec. V and shown in Fig. 3. We find that the purity of the QE-like sample depends on both the recoil energy and on the Q 2 QE of the interaction, with high Q 2 QE interactions having larger recoil (see Fig. 4). We therefore apply a Q 2 QE -dependent cut on the recoil energy:

VII. CROSS-SECTION EXTRACTION
The double-differential cross section versus variables x and y in bin ði; jÞ is constructed using Recoil energy for data (points) and simulation (colors, POT normalized to data). In the simulation, signal events include both quasielastic events (purple), 2p2h scatters (green) and resonant or DIS events (pink and red) with a QE-like signature. The backgrounds consist of quasielastic and 2p2h events with a non-QE-like signature (hatched purple and green) and non-QE events (resonant and DIS) without a QE-like signature (hatched pink and red). All cuts except the recoil cut are applied.
where N αβ is the number of data events reconstructed in bin ðα; βÞ, N bkgd αβ is the estimated number of background events reconstructed in bin ðα; βÞ, U is an unfolding operation transforming reconstructed bin ðα; βÞ to true bin ði; jÞ, ϵ ij is the product of reconstruction efficiency and detector acceptance for events in true bin ði; jÞ, Φ is the flux of incoming antineutrinos (either integrated or for the given bin, as described later), T is the number of scattering targets (here, the number of nucleons), and Δx i (Δy j ) is the width of bin i (j).
We report our primary cross-section measurement in bins of muon transverse (p T ) and longitudinal momentum (p k ) with respect to the neutrino beam direction. We choose these as our primary results as they are quantities that we have directly measured. For comparison with other experiments, we also report auxiliary measurements versus Q 2 QE and E QE ν , both reconstructed in the quasielastic hypothesis from the muon kinematics [see Eqs. (2) and (3)]. The bin boundaries are shown in Table I.
Two bins at highest p T and lowest p k and four bins at highest Q 2 QE and lowest E QE ν are not reported due to poor acceptance in those regions. Note that, as Q 2 QE and E QE ν are reconstructed from the muon kinematics, they are both functions of both p k and p T . Figure 5 shows lines of constant Q 2 QE and E QE ν , projected onto the p k /p T phase space. For most of the region considered by this analysis, E QE ν correlates fairly well with p k and Q 2 QE with p T . This simplification breaks down at high p T and low p k . For both versions of the double-differential cross sections, we also report projections onto each axis, resulting in onedimensional distributions of p T , p k , Q 2 QE , and E QE ν . For the single-differential cross section versus E QE ν , we report a flux-weighted cross section, where each bin has been divided by the flux integrated over the energy range of that bin only, rather than the entire antineutrino flux integrated over all energies. Note that care must be taken in interpreting this quantity, as E QE ν does not correspond exactly to true antineutrino energy.
A total of 17 621 interactions pass our reconstruction cuts for data. Distributions of these events versus muon p T in bins of p k are shown in Fig. 6.

A. Background subtraction
The term N bkgd αβ in Eq. (9) refers to the estimated number of reconstructed data events that correspond to background processes. Recall that our QE-like signal, explained in Sec. VI, is defined as having a final state containing a μ þ , any number of neutrons, any number of protons with less than 120 MeV kinetic energy, and no pions, other hadrons, or prompt photons. Thus, background events in our sample could, for example, correspond to resonant events with pions that did not make a track and that generated recoil distributions that fell within our cuts. Figure 7 shows p T , p k , E QE ν , and Q 2 QE distributions in the data and simulation, with the latter subdivided into signal and background. Backgrounds in this analysis arise primarily from events involving charged pions. MINERvA's charged pion production analysis [87] suggests that GENIE 2.8.4 overpredicts the rate of resonant pion production. We therefore use a data-driven fitting procedure to constrain the backgrounds predicted by GENIE. Since the constraint can in principle be different in each p T /p k bin, the fit would ideally be done separately in each bin. However, the limited statistics of our data sample caused attempts to fit each bin separately to fail. The fits are instead performed separately for five regions of the p T /p k phase spaces, chosen by combining p T /p k bins with similar background shapes.
For each of the five regions, the recoil energy, after all other cuts, is compared for data and for signal and background Monte Carlo. The TFractionFitter tool, part   Table I. of the ROOT framework [88], is used to perform a fractional fit of the simulation to data, where the relative normalization of the signal and background distributions is allowed to vary. The shapes of the distributions are not varied. Figure 8 shows the recoil distributions in data and (area-normalized) simulation for one of the five regions of p T /p k , before and after tuning the signal and background fractions. In each bin, a scale is extracted corresponding to the factor by which the background fraction was rescaled relative to the nominal simulation to give the best fit. The estimated background fraction in each bin of the data distribution corresponds to the background fraction of the Monte Carlo in that bin, multiplied by this scale factor.
The scales for the p T versus p k regions are shown in Table II. In most cases, as suggested by [89], the simulation is found to predict too high a fraction of background events. Figure 9 shows the signal fraction as a function of the muon kinematic variables. After background subtraction, the signal data sample has estimated 14 839 events.

B. Unfolding
Detector smearing is corrected using a migration matrix that describes the relationship between true and reconstructed bins of p T and p k . The migration matrix for our simulated reconstructed QE-like signal distribution is shown in Fig. 10. The x axis indicates bins in the reconstructed variables, where the bins of p k are repeated for each bin of p T . The y axis indicates bins in the true variables, arranged in the same way. Thus any events on the diagonal were reconstructed in the correct bin of both p k and p T . An event reconstructed in the wrong bin of p k (but the right p T bin) will be displayed in another bin in the same subplot; one reconstructed in the wrong p T bin will appear in a different subplot.
We use the iterative method of D'Agostini [90], as implemented in the ROOT package RooUnfold [91], with four iterations. The unfolding procedure was validated using an ensemble test, in which ten data-sized subsamples of the simulation were selected and warped by an adjustment of the quasielastic axial mass by AE25%. These samples were then unfolded using the migration matrix generated from the full unwarped simulation; the warped simulation was recovered within four iterations.

C. Efficiency and acceptance correction
The unfolded distributions are then corrected for detector acceptance and reconstruction efficiency. The most significant effect on acceptance is from the requirement that QE and E QE ν . Here the MINERvA -tuned GENIE simulation is absolutely normalized to the POT of the data sample, but the background corrections described in this section have not yet been applied. The Q 2 QE distribution is shown on a log scale to highlight the high Q 2 QE region. In the simulation, signal events include both quasielastic events (purple), 2p2h scatters (green) and resonant or DIS events (pink and red) with a QE-like signature. The backgrounds consist of quasielastic and 2p2h events with a non-QElike signature (hatched purple and green) and non-QE events (resonant and DIS) without a QE-like signature (hatched pink and red).
final-state muons are matched in MINOS, limiting the muon's angle with respect to the beam line to a maximum of 20°. The MINOS-match requirement also limits our ability to accept muons with low longitudinal momentum ≲1.5 GeV/c which will stop in MINERvA or not produce enough activity to be analyzed in the MINOS spectrometer. The largest source of inefficiency is due to the Q 2 QEdependent E recoil cut.
We estimate the product of acceptance and efficiency using the full MINERvA -tuned GENIE+GEANT4 simulation: where N generated and reconstructed ij is the number of simulated events generated in p T bin i and p k bin j that also pass all reconstruction cuts (except the fiducial cuts on position and muon angle) and N generated ij is the total number of events generated in p T bin i and p k bin j. Figure 11 shows the product of efficiency and acceptance versus p k and p T . The low acceptance at high p T and low p k is due to the MINOS match requirement and angle cut. The efficiency also decreases at higher energies, where interactions are more likely to include large amounts of recoil energy and may be vetoed by our Q 2 QE -dependent E recoil cut. The overall efficiency × acceptance of the sample is 52.5%.

D. Flux and target number correction
To convert an acceptance-corrected distribution to a cross section, we divide by the number of nucleons, the total number of POT producing the neutrino beam, and the estimated antineutrino flux per POT. These are summarized in Table III.
The NuMI beam's flux prediction is explained in detail in [64] and is summarized in Sec. IV. For distributions in the p T /p k phase space, we report flux-integrated cross sections. We do the same for the single-differential cross section dσ/dQ 2 QE . We integrate over the entire available flux range of 0-100 GeV, to get a total integrated flux of 2.295 × 10 −8 cm −2 per proton on target.
For the cross section as a function of E QE ν /Q 2 QE , one can create an approximate flux-weighted cross section, where the number of events in each E QE ν bin is normalized by the flux in the corresponding E ν bin. This is not a true total cross section, since E QE ν is not the true neutrino energy, except in the case of quasielastic scatters off of hydrogen. However E QE ν is closely correlated to E true ν (see Fig. 12), making the flux-weighted cross section a close approximation of the total cross section versus energy.  The target for an antineutrino quasielastic scatter [Eq. (1)] is a proton. For QE-like scattering, it is possible that a scattering process could originate on a neutron (e.g.ν μ n → μ þ Δ − ), where the resonance decays Δ − → nπ − and the pion is absorbed, or on a nucleon pair. We use the total number of nucleons in the fiducial volume as the target number normalization. The fiducial volume is made up of a combination of polystyrene, doping agents, epoxy and light-tight coating. The predominant material is polystyrene, which is composed of equal parts carbon and hydrogen. A full summary of the composition of the MINERvA tracker is available in [61]. We estimate that the fiducial volume used in this analysis contains 3.23 × 10 30 nucleons, of which 1.76 × 10 30 are protons.

VIII. SYSTEMATIC UNCERTAINTIES
Systematic uncertainties on the cross-section measurements arise from many sources. To assess these systematic uncertainties, we vary parameters in the simulation within their uncertainties and recalculate the cross sections using new estimates of efficiency, backgrounds, unsmearing, flux and target number corrections. The difference between this new cross section and the original result is taken to be the systematic uncertainty on the   cross section due to that source. The sources of systematic uncertainty are discussed below. Systematic uncertainties in each category, and total systematic uncertainties, are available in Appendix A.

A. Flux uncertainties
The simulation of the NuMI flux and its uncertainties are described in detail in [69]. Uncertainties in the antineutrino flux arise primarily from uncertainties in hadron production rates and in parameters that control the alignment of the NuMI focusing system, such as the position of the focusing horns. These uncertainties are constrained with both external data and with a MINERvA measurement of elastic neutrino scattering on electrons [70]. The total uncertainty in the focusing peak is approximately 8% and rises to 11% at the falling edge of the focusing peak, where beam focusing uncertainties are large.
Muon tracking efficiencies in MINERvA and MINOS are measured by reconstructing tracks in one detector, extrapolating to the other detector, and observing the fraction of tracks matched in both detectors in data and in the simulation. The simulation is corrected for small discrepancies between tracking efficiencies, 0.5% AE0.25% for MINERvA and 0.5 (2.5)% AE0.25 (1.25)% for MINOS for muons with momentum greater than (less than) 3.0 GeV.
Potential angular reconstruction biases are estimated both by cutting tracks in half and comparing the reconstructed angles of both halves, as well as studies of forward-going events such as neutrino-electron scattering and low hadronic recoil events. These studies limit additional angular smearing or bias in the data relative to the simulation to below 1 mrad.
Smearing of reconstruction vertices causes some events within the fiducial volume to be misreconstructed outside the fiducial volume, and vice versa. We estimate the uncertainty due to this effect by smearing reconstructed vertices in the simulation by 0.9 mm in x, 1.25 mm in y and 1 cm in z; this results in a negligible change in measured cross section. The MINERvA tracker consists of primarily scintillator strips, with smaller portions of epoxy, tape, reflective coating and wavelength shifting fibers. The total uncertainty on the mass of the tracker is 1.4% [61].

C. Model uncertainties
Models used in the simulation include various parameters that carry uncertainties. These include uncertainties in signal, background and final-state interaction models. Most of these are evaluated using the reweighting prescription and parameter uncertainties recommended by the GENIE Collaboration [1]. These parameters are listed in Table IV, along with the amount by which they are varied and the approximate effect on the cross sections. GENIE uncertainties that change particle fates cannot be modeled using the reweighting method. In this case, we generate an alternative simulated sample in which these parameters, including the effective nuclear radius, formation zone and hadronization model, have been adjusted.
The RPA correction described in Sec. IV B is applied when calculating the central values. The correction is varied within the uncertainties shown in Fig. 2.
Similarly, the addition of the 2p2h process is estimated by adding events from correlated pairs as described in Sec. IV and Ref. [93]. The uncertainty on this is determined by using several variations of the tuning procedure, including fits that allow interactions on pp pairs, np pairs, or single-nucleon interactions to be tuned. The differences in cross section obtained using these three variants from the standard simulation are added in quadrature as a systematic error due to the 2p2h model. Table V summarizes the effects of these variations on the extracted cross sections.

D. Recoil reconstruction uncertainties
Several sources of uncertainties can affect the reconstruction of recoil energy, which can in turn change the background estimates and efficiencies used to estimate the cross sections. Quasielastic antineutrino events have a hadronic final state consisting of a neutron. In order for neutrons to deposit recoil energy in the detector, they must undergo an interaction, with resulting charged particles (usually protons) then depositing visible energy. The most significant source of uncertainty associated with neutrons is due to the GEANT4 neutron interaction model. To evaluate this uncertainty, we vary the mean free path of neutrons in the detector, with the variations spanning discrepancies between GEANT4 and thin target neutron scattering data on copper, iron and carbon [94][95][96][97][98][99][100][101].
Energy response of protons has been measured in the MINERvA test beam detector [82]. To propagate uncertainties on this measurement to the cross sections, we shift simulated recoil energy deposited by protons by uncertainties derived from comparisons of the test beam measurements and GEANT4. The variation depends on the proton energy: 4% below 50 MeV and 3.5% above 50 MeV. The proton response affects our event rate measurement by less than 1% across our whole phase space, as the track cut removes many protons, and only a small amount of those that remain pass into our selected sample by making this shift. The pion calorimetric response has also been constrained by test beam studies to an accuracy of 4%, for pions with a kinetic energy between 400 and 1900 MeV. We thus separate our pions into two categories-"constrained" within this energy range and "unconstrained" outside of it. Pions within the constrained range have their energy fraction varied by AE4%, while others have it varied by AE5%. The pion response has only a minor effect (< 1% across our whole phase space) on our cross sections.
For the other particles (electromagnetic and kaons), we vary the recoil by AE3%. This uncertainty was derived by observing the energy response for Michel electrons (electrons from muon decay), which have a well-known energy spectrum. This change mainly affects the Q 2 (p T ) shape and contributes its maximum of around 1% uncertainty at low Q 2 .
These uncertainties are dominated by neutron interaction modeling, which ranges from 2-6%; the other uncertainties are less than 1%.

IX. RESULTS
Double-differential cross sections versus p T and p k are shown in Fig. 13. Double-differential cross sections versus E QE ν (E true ν ) and Q 2 QE are shown in Fig. 14 (Fig. 15). In each case, simulated cross sections are also plotted, where the simulation uses the MINERvA -tuned GENIE model described in Sec. IV. Results corrected to a quasielastic, rather than QE-like, signal definition are available in Appendix B.
The MINERvA -tuned GENIE model agrees well with MINERvA data, in spite of the fact that the tune was made to an independent (neutrino rather than antineutrino) data set. In the following sections, we compare these results with many alternate models and discuss the impact of the individual components of the MINERvA tune. Figure 16 shows the data and the different components included in the MnvGENIE model. Events generated as pure 1p1h CCQE are shown in green; the resonance (and a A. Comparisons to alternate GENIE models Figure 17 shows the measured differential cross sections and several variations on the GENIE model; single-differential projections versus p T , p k , E QE ν , and Q 2 QE are shown in Fig. 18 for the same variations. In particular, the MINERvA-tuned GENIE (MnvGENIE in figures) model is GENIE with RPA and MINERvAtuned 2p2h effects added. 2 Other permutations of RPA, 2p2h and MINERvA -tuned 2p2h are also shown.  Figure 19 shows the ratio of the measured differential cross section to the MINERvA-tuned GENIE model. Table VI shows the χ 2 for 58 degrees of freedom for the models shown in that figure and for additional theoretical models described in Sec. IX.
A standard χ 2 comparison for these models relative to the data using the statistical uncertainties derived from the data gives a best agreement (the green curve with RPA but no 2p2h contribution) with the curve that lies furthest away from the data. This is due to the dominance of multiplicative normalization uncertainties in the covariance matrix, which leads to the well-known pathology of Peelle's pertinent puzzle [102,103]. This effect is well documented in the nuclear cross-section literature [104] and, in the limit of pure multiplicative uncertainties, a χ 2 derived from the log of the cross section is preferred to one derived from the cross section itself, since in the former case the multiplicative factors are normally distributed.
In addition, the χ 2 statistic is known to have biases when uncertainties estimated from the counting statistics of individual data points are used. For statistical uncertainties estimated to be proportional to ffiffiffiffiffi N i p , points that fluctuate downwards are given smaller estimated uncertainties and hence greater weight in the χ 2 calculation. For normalization uncertainties, the effect is even greater with the uncertainty being directly proportional to N i . The normalization uncertainties in these data are highly correlated from bin to bin and substantial relative to the other uncertainties. For this reason, we report the χ 2 using both the cross section itself (linear) and the log of the cross section (log-normal) in Table VI in the next section. The lowest log-normal χ 2 is for the MINERvA -tuned GENIE model with default 2p2h and the RPA correction which appears as the orange curve in Figs. 17-19.

B. Comparisons to NuWro models
We have also compared the data to several models available in the NuWro event generator. Table VI summarizes the agreement between the data and these models while Figs. 20 and 22 show the comparisons. NuWro also includes models of 2p2h and RPA and additionally includes an implementation of the transverse enhancement model describe in Sec. II. The relativistic Fermi gas nuclear model is labeled GFG (global Fermi gas) to distinguish it from an alternate LFG nuclear model. A spectral function model is also available in NuWro and included in the comparisons.
All of the NuWro models have higher χ 2 values than the MINERvA -tuned GENIE model. Even when comparing very similar primary interaction models (e.g. default GENIE 3 and NuWro GFG without RPA or 2p2h) between the two generators, the agreement with data is quite different. We believe this is due to the different FSI models used by NuWro and GENIE, which impact the predicted contribution to the cross section from events that include a pion that is absorbed before exiting the primary nucleus. Of the NuWro models, the preferred model includes RPA and 2p2h contributions, as is also the case with GENIE variants described above.

Figures 23 and 24 show the cross sections versus E QE
ν , corrected to cross section or proton, compared to results from MiniBooNE [85] and NOMAD [45]. The MiniBooNE cross sections quoted are the average of their reported cross sections on mineral oil and their estimated rates on pure carbon, as our scintillator target lies approximately halfway between those two compositions. NOMAD is only shown for the true-QE assumption as they only quote results for that process. We note that the caveats discussed in Sec. VII should be taken into account when comparing these results to other experimentsnamely that this is an approximation of the energydependent cross section and that the approximations will have different impact on these results than those measured in beams with different neutrino energy spectra. Although the MINERvA data points show a small dip in the TABLE VI. This table summarizes the χ 2 values for comparisons of the differential cross section dσ 2 /dp T dp k to a wide variety of models. The top five are GENIE variations illustrated in Fig. 17 while the bottom seven are variations of the NuWro event generator illustrated in Figs. 20 and 21. The MINERvA -tuned GENIE model includes RPA and MINERvA -tuned 2p2h and was used in the extraction of the cross section. The χ 2 is for 58 degrees of freedom.

Model
Conventional Our "default" GENIE includes a correction to the single nonresonant pion production rate based on bubble chamber inputs discussed in Sec. IV. This has little effect on the CCQE-like cross-section prediction, since very few CCQE-like events arise from nonresonant pion production. ∼4-6 GeV region, they are consistent within uncertainties with models that predict a smoothly rising cross section in this region (see Fig. 18). The MINERvA neutrino flux [69] changes rapidly in the 4-6 GeV region with uncertainties that are dominated by the neutrino beam focusing system (see the distribution of uncertainties versus E QE in Fig. 26 of Appendix A).

D. Vertex energy distributions
Because interactions on multinucleon pairs are expected to include additional low-energy nucleons compared to standard QE interactions, reconstructed energy near the interaction vertex is useful for judging the efficacy of 2p2h models. Figure 25 shows energy reconstructed in scintillator strips that are within 100 mm of the interaction vertex, in the sample used to produce the cross sections discussed earlier, but before background subtraction and efficiency, flux and target number corrections. Also shown are the expected distributions for default and MINERvA -tuned GENIE and ratios to MINERvA -tuned GENIE for the data and several GENIE variants. Models that omit a 2p2h component have very poor agreement with the data, but the case for RPA suppression is not as strong. The model that lies closest to the data is the MINERvA -tuned GENIE with the RPA correction omitted. This is in conflict with similar conclusions drawn from cross sections versus kinematic distributions, indicating that while MINERvA -tuned GENIE is a definite improvement over default GENIE 2.8.4, more improvements are needed to properly simulate the hadronic component of anti-neutrino QElike interactions.

X. CONCLUSION
We have presented a measurement of a QE-like cross section for antineutrino scattering on scintillator. The signal definition requires no charged pions in the final state and no protons with kinetic energies above 120 MeV. This variant of the QE-like definition allows us to include quasielastic scatters off of NN pairs in the nucleus in our signal definition and closely matches the actual sensitivity of our detector to low-energy protons. The main result is presented as a function of muon kinematics p T and p k . We also present an energydependent flux normalized cross section in terms of the neutrino energy and 4-momentum transfer squared as calculated from the muon kinematics using a quasielastic assumption. These data are compared to a large number of models for antineutrino QE-like scattering. In particular, we have applied corrections to the default GENIE 2.8.4 scattering model for the random phase approximation and have added 2p2h processes that have been tuned to the observed recoil distributions in an independent MINERvA neutrino scattering sample. This MINERvA -tuned model agrees better with our data than default GENIE both visually (Fig. 19) and when a log-normal χ 2 is calculated, as is more appropriate when multiplicative uncertainties are significant. Moreover, comparisons of reconstructed energy near the interaction vertex between these data and various models indicate poor agreement with models that do not include 2p2h. In conclusion, addition of RPA and 2p2h effects to the simulation substantially improves agreement with the MINERvA QE-like data over default GENIE. Addition of either RPA or 2p2h alone is not sufficient. However, substantial discrepancies between the improved model and data remain, indicating that more model development is needed. This is the first double-differential measurement of quasielastic or QE-like scattering cross sections for antineutrinos in this energy range, which is very similar to the expected spectrum of the DUNE experiment, and will be an essential component in the development and tuning of models used in future neutrino oscillation measurements.

ACKNOWLEDGMENTS
This document was prepared by members of the MINERvA Collaboration using the resources of the Fermi National Accelerator Laboratory (Fermilab), a U.S. Department of Energy, Office of Science, HEP User Facility. Fermilab is managed by Fermi Research Alliance, LLC (FRA), acting under Contract No. DE-AC02-07CH11359. These resources included support for the MINERvA construction project, and support for construction also was granted by the United States National Science Foundation under Grant No. PHY-0619727 and by the University of Rochester. Support for participating scientists was provided by NSF and DOE (USA), by CAPES and CNPq (Brazil), by CoNaCyT (Mexico), by Proyecto Basal FB 0821, CONICYT PIA ACT1413, Fondecyt 3170845 and 11130133 (Chile), by DGI-PUCP and UDI/VRI-IGI-UNI (Peru), the UK Science and Technology Facilities Council (UK), and by the Latin American Center for Physics (CLAF). We thank the MINOS Collaboration for use of its near detector data. Finally, we thank the staff of Fermilab for support of the beam line, the detector and computing infrastructure. Figure 26 shows a summary of the fractional systematic uncertainties from each category of uncertainty described in Sec. VIII for the one-dimensional cross sections. The model uncertainties have been further subdivided into those primarily affecting the signal models (quasielastic and 2p2h), background models, and final-state interactions. Figures 27 and 28 show fractional uncertainties on the double-differential cross sections. Figure 29 shows the correlation matrix for all systematic uncertainties.  The main focus of the analysis was the calculation of CCQE-like double-differential cross sections shown above, which correspond to our measurement for the signal definition described in Sec. VI. As an extension to the analysis, however, we also calculated a true CCQE cross section. Recall that, for the CCQE-like cross section, our signal corresponded to interactions with a CCQE-like final state, even if that final state was generated by a resonant or DIS interaction followed by FSI. For the true CCQE definition, our signal corresponds only to events where the initial interaction was quasielastic, even if FSI created final-state particles such as pions that mimicked a nonquasielastic interaction. The signal also includes 2p2h events where a CCQE interaction takes place on a correlated pair.

APPENDIX A: SUMMARY OF SYSTEMATIC UNCERTAINTIES
The true CCQE double-differential cross sections are shown in Fig. 30 (d 2 σ/dp T dp k ) and 31 [dσðE QE ν Þ/dQ 2 QE ], while one-dimensional projections are shown in Fig. 32. Also shown in these figures are the predictions in our default simulation, which includes both CCQE and 2p2h contributions. a. p T and p k The differential cross section with respect to muon parallel and transverse momentum, d 2 σ/dp T dp k , is shown in Table VII. The statistical uncertainty on these measurements are shown in Table VIII, and the systematic uncertainty in Table IX.
Table X shows the differential cross section dσ/dp T , generated by projecting the two-dimensional measurement onto the p T axis. Table XI shows the differential cross section dσ/dp k , generated by projecting the twodimensional measurement onto the p k axis. Both tables also include statistical and systematic uncertainties.  Tables XXV, XXVI,  XXX and XXXI.   TABLE VII. Measured double-differential CCQE-like cross section d 2 σ/dp T dp k . Units are 10 −41 cm 2 /GeV 2 /nucleon. Columns represent bins of p T (GeV); rows are bins of p k (GeV).

CCQE cross sections
The tables in this section list our cross-section measurements for the true quasielastic signal definition explained in Sec. B. In each case, we show a table of values corresponding to the cross section, followed by a table of   TABLE VIII. Statistical uncertainty on the measured doubledifferential CCQE-like cross section d 2 σ/dp T dp k . Units are 10 −41 cm 2 /GeV 2 /nucleon. Columns represent bins of p T (GeV); rows are bins of p k (GeV).    a. p T and p k The differential cross section with respect to muon parallel and transverse momentum, d 2 σ/dp T dp k , is shown in Table XXII. The statistical uncertainty on these measurements is in Table XXIII, and the systematic uncertainty in Table XXIV.
Table XXV shows the differential cross section dσ/dp T , generated by projecting the two-dimensional measurement onto the p T axis. Table XXVI shows the differential cross Statistical uncertainty on the measured doubledifferential true CCQE cross section d 2 σ/dp T dp k . Units are 10 −41 cm 2 /GeV 2 /nucleon. Columns represent bins of p T (GeV); rows are bins of p k (GeV). section dσ/dp k , generated by projecting the two-dimensional measurement onto the p k axis. Both tables also include statistical and systematic uncertainties.
b. E QE ν and Q 2

QE
The reconstructed energy-dependent true CCQE cross section versus Q 2 QE , dσðE QE ν Þ/dQ 2 QE , is shown in Table XXVII. The statistical uncertainty on these measurements is in Table XXVIII, and the systematic uncertainty in Table XXIX. By projecting, we can get the total reconstructed energy-dependent cross section σðE QE ν Þ shown in Table XXX and the single-differential cross section dσ/dQ 2 QE shown in Table XXXI.  XXIV. Systematic uncertainty on the measured double-differential true CCQE cross section d 2 σ/dp T dp k . Units are 10 −41 cm 2 /GeV 2 /nucleon. Columns represent bins of p T (GeV); rows are bins of p k (GeV).

Bin
Cross section Statistical uncertainty Systematic uncertainty