Measurement of Quasielastic-Like Neutrino Scattering at $\left\sim 3.5$~ GeV on a Hydrocarbon Target

MINERvA presents a new analysis of neutrino induced quasielastic-like interactions in a hydrocarbon tracking target. We report a double-differential cross section using the muon transverse and longitudinal momentum. In addition, differential cross sections as a function of the square of the four-momentum transferred and the neutrino energy are calculated using a quasielastic hypothesis. Finally, an analysis of energy deposited near the interaction vertex is presented. These results are compared to modified GENIE predictions as well as a NuWro prediction. All results use a dataset produced by $3.34\times10^{20}$ protons on target creating a neutrino beam with a peak energy of approximately 3.5 GeV


I. INTRODUCTION
Charged current quasielastic (CCQE) scattering is an important signal process for neutrino oscillation experiments [1] [2][3] [4]. In this process an incoming neutrino exchanges a W boson with a neutron to produce an outgoing charged lepton and a proton. The initial state neutron is not independent but bound inside the nucleus with an initial Fermi momentum. Thus the kinematics of this apparently simple process are significantly modified compared to the predictions from a straightforward two body scattering event. In addition, the presence of the nucleus can cause other processes, for example resonance production and decay followed by pion absorption, to appear quasielastic in a detector. All of these modifications taken together affect not only the total cross section and outgoing hadronic system, but also the outgoing lepton kinematics.
We present in this paper the two-dimensional cross sections for quasielastic-like scattering as a function of muon transverse and longitudinal momentum. Momentum of the final state lepton can be measured precisely and without reference to an interaction model in such neutrino interactions. These muon observables are suitable for comparing to models of inclusive neutrino scattering that can separate events by final state. Using the reconstructed muon angle and momentum the quantities E ν,QE and Q 2 QE are reconstructed. These are the reconstructed neutrino energy (Eq. 1) and square of the four-momentum transferred (Eq. 2) for the limiting case where the initial state nucleon was at rest and there were no nuclear effects. The binding energy, E b , used in this analysis is 34 MeV.
A differential cross section in Q 2 QE is reported which extends to higher and lower Q 2 QE than previously reported in Ref. [5], which was based on a subset of the data presented in this work. Additional studies associated with the vertex energy are presented in this paper, which build on the results of Ref. [5].

II. EXPERIMENT
The MINERvA experiment employs a fine-grained tracking detector for recording neutrino interactions produced by the NuMI beamline at Fermilab [6]. Neutrinos are produced by directing 120 GeV protons from the Main Injector onto a graphite target. The produced charged pions and kaons are focused by two magnetic horns. The charged particles travel along a 675 m long, helium-filled pipe and decay primarily into muon-flavored neutrinos and muons. The polarity of the horns can be switched to produce either a neutrino-dominated or antineutrino-dominated beam at the peak neutrino energy. Approximately 97% of the muon neutrinos that reach the MINERvA detector are produced by pion decay; the remainder come from kaon decay [6,7]. The neutrino beam is inclined at an angle of 58 mrad down rel-ative to the horizontal (z) direction at the MINERvA detector.
The MINERvA detector [8] consists of 120 hexagonal modules. These modules comprise an active tracking volume and a suite of nuclear targets positioned upstream, which are not used for this analysis. This analysis studies the interactions in the scintillator tracking volume which has a fiducial mass of 5.48 tons, and is surrounded by electromagnetic and hadronic calorimeters.
Each tracking module consists of two planes of triangular nested polystyrene scintillator strips with a 1.7 cm strip-to-strip pitch. There are three different plane configurations rotated 0 • and ± 60 • relative to the longitudinal axis of the detector. The modules are arrayed in a stereoscopic combination along the horizontal axis. This enables a three dimensional reconstruction of the neutrino interactions via their charged particles in the final state. The downstream and side electromagnetic calorimeter consists alternating layers of scintillator and 2 mm thick lead planes while alternating scintillaor and 2.54 cm thick steel planes constitute the downstream hadronic calorimeter.
Each scintillator strip is read out by a wavelengthshifting fiber connected to a multi-anode photomultiplier tube. The 3.0 ns timing resolution of the readout electronics is sufficient for separating multiple interactions within a single NuMI beam spill.
The MINOS near detector [9], a scintillator-based magnetized iron spectrometer, is located 2 m downstream of the MINERvA detector. It is used to reconstruct the momentum and charge of muons that originate inside MINERvA and enter MINOS.
This analysis uses data that comprise an integrated 3.34 × 10 20 protons on target received between March 2010 and May 2012, when the focusing horns were set to produce a broad band neutrino beam peaked at 3.5 GeV with a purity of 95% muon neutrinos at the peak.

III. SIMULATION
The neutrino beam is simulated using a GEANT4based description of the NuMI beamline [6].
Because GEANT4 [10] does not exactly reproduce NA49's hadron production measurements of 158 GeV protons on carbon [11], or other relevant hadron production measurements, the simulation is reweighted as a function of the pion transverse and longitudinal momentum to agree with those data, as described in Ref. [7]. In addition, the flux is constrained by the measurement, in the same beamline, of neutrino scattering off atomic electrons, as described in Ref. [12].
Neutrino interactions are simulated using the GENIE 2.8.4 neutrino event generator [13]. The quasielastic interaction (1p1h) is simulated by the Llewellyn-Smith formalism [14], and the vector form factors are modeled using the BBBA05 model [15]. The dipole form is used for the axial vector form factor, with an axial mass of M A = 0.99 GeV/c 2 . The Rein-Sehgal model [16] is used to simulate resonance production, and the resonant axial mass used is M RES A = 1.12 GeV/c 2 . The DIS cross sections used by GENIE are a leading order model that use the Bodek-Yang prescription [17] for the low Q 2 modification.
The nuclear environment is simulated in GENIE using the relativistic Fermi gas model [18] in combination with the Bodek-Ritchie high momentum tail [19] on the nucleon momentum distribution. The maximum momentum for Fermi motion is assumed to be k F = 0.221 GeV/c. GENIE models intranuclear rescattering, or final state interactions (FSI), of the produced hadrons using the INTRANUKE-hA package [20].
MINERvA makes three sets of corrections to this version of GENIE that are based on improved models or measurements of the underlying processes. First, we reduce the cross section for quasielastic events as a function of energy (q 0 ) and three momentum transfer (q 3 ) ( or four-momentum transfer (Q 2 ) above 3 GeV 2 ) based on the RPA part of the Valencia model [21][22] appropriate for a Fermi gas [23] [24]. Second, we add a prediction for multinucleon (2p2h 1 ) scattering from this same Valencia model [25][26] [27]. However, we know from studies of an inclusive sample of MINERvA data at low momentum transfer that this model significantly underpredicts the measured cross section in a specific region of q 0 and q 3 where 2p2h processes contribute up to one-third of the rate [28]. We therefore enhance the 2p2h prediction from the Nieves model in this region. Integrated over all phase space the rate of 2p2h is increased by 53%. Third, we decrease the non-resonant pion production by 43% to better agree with a fit to measurements on deuterium of that process [29], and take the uncertainties on that process from the same reference rather than use the GENIE-recommended uncertainties. This modified version of the simulation is referred to later in this paper as MINERvA GENIE tune v1 and MnvGENIE-v1 in figures.
The response of the MINERvA detector is simulated using a GEATNT4-based model which uses version 4.9.4p6 GEANT4 [10]. The optical and electronics performance is also simulated, and the absolute energy scale of minimum ionizing energy depositions is set by requiring the average and RMS of the energy deposits by through-going muons to match between the data and simulation as a function of time, as described in Ref. [8]. The absolute energy response of the detector to charged hadrons was determined by placing a scaled-down version of the MINERvA detector in a charged particle test beam, as described in Ref. [30]. In order to simulate the 1 The notation "2p2h" is meant to distinguish hard scattering process with momentum and energy transferred to a pair of nucleons from the standard quasielastic process, "1p1h", where the momentum and energy is transferred to a single nucleon in the hard scattering process.
effects of accidental activity, such as deadtime and reconstruction failures, each simulated neutrino interaction is embedded in an event from the data in both the MIN-ERvA and MINOS detector. These spills from the data are taken from the same periods (typically determined by a configuration change) the Monte Carlo simulates.

IV. SIGNAL DEFINTION
Quasielastic interactions are those where there are one or more nucleons in the final state, but no mesons or excited or heavy baryons. Charged current neutrino interactions where a pion is produced in the primary process but is absorbed in final state interactions appear as quasielastic in the MINERvA detector and are therefore also considered signal events in this analysis.
This analysis uses a quasielastic-like signal, defined as those events with the following final state particles: • One negatively charged muon of angle < 20 • with respect to the neutrino beam Events containing low energy (≤ 10 MeV) photons are accepted because they can arise from nuclear de-excitation in atoms and might not be part of the primary neutrino interaction products. Because of the geometric acceptance (MINOS-MINERvA overlap) and the need for muon charge and momentum identification, a restriction on the muon angle of < 20 • is required.

V. EVENT RECONSTRUCTION
In order for the quasielastic-like cross sections to be reported as a function of transverse and longitudinal muon momentum, event reconstruction must, at a minimum, measure the momentum, charge, and direction of the muon. If another track consistent with a proton is found, then the proton's kinetic energy and direction can also be measured. Recoil energy is defined as the event activity that is not coming from either the muon or any tracked protons. The absence of pions can be ensured by a cut on the recoil energy, and by requiring no evidence for a Michel electron from the π → µ → e decay chain. Because low energy protons can deposit energy close to the vertex and not be reconstructed as tracks, the recoil energy does not include energy that is less than 150 mm from the vertex.
Every event is required to have a reconstructed primary interaction vertex inside the fiducial volume of the scintillator region of the MINERvA detector. The fiducial volume is defined as the region between scintillator modules 27 and 79 of the detector with an apothem of 850 mm giving a 5.48 metric ton fiducial mass. The tracker composition is described in Sec. VII D. The remaining aspects of event reconstruction are described in the following sections.

A. Muon Track Reconstruction
This analysis requires events to contain one long track that originates in the fiducial volume defined above, traverses the remainder of the MINERvA detector, and extrapolates to a track in MINOS whose momentum and charge are measured. The momentum of the muon as it enters the upstream face of the MINOS detector is measured on the basis of its range or track curvature information, as described in Ref. [9]. The total muon momentum must also include the momentum required to traverse the distance from the event vertex to the downstream end of the MINERvA detector. This momentum is determined by computing the minimum ionization loss for the material traversed.
Long track signatures in MINERvA are categorized as those that possess a contiguous set of clusters for at least 8 planes of the detector. Most of these tracks are made from muons, although sometimes they are made from energetic protons or pions.
The MINERvA detector uses planes of scintillator bars in three different orientations (vertical and ±60 • from vertical) to reconstruct tracks in three dimensions. A long track is found by looking for tracks in at least two of the three views, and for events with a long track found in all three views, they must all be geometrically consistent. Given the geometry of the detector, a muon must traverse at least 9 (11) planes in order to create tracks in two (three) views.
A track matching algorithm ensures that the reconstructed part of the muon track in MINERvA matches with its corresponding momentum-analyzed counterpart in MINOS, by comparing their track timings and extrapolated positions at the back/front of the MIN-ERvA/MINOS detector.
Once a muon track has been identified in a beam spill, an event is formed with a specific (length-corrected) average time for the muon, a vertex, and a set of clusters which will no longer be considered in subsequent reconstruction steps that characterize other aspects of the event. For the remainder of Section V, only clusters that are between 20 ns before and 35 ns after the average muon time are considered. In addition, the energy in each cluster is corrected for the passive material that is upstream of the cluster.
The corrected energy of clusters are further reconstructed into different objects: short tracks, showers, vertex energy, residual energy, and recoil energy. The remainder of this section describes how these objects are formed.

B. Short Track Reconstruction
Low energy (less than 100 MeV) protons or other heavy baryons typically produce short tracks. These traverse just a few planes in the MINERvA detector. There are two tracking algorithms employed for reconstructing short tracks at the neutrino interaction vertex.
Both the algorithms input clusters that are not yet associated with any reconstructed object in the event. The tracking threshold is 4 clusters, which means that particles must have deposited energy in at least 4 different planes to be considered for short tracks.
The first algorithm starts by making track seeds out of sets of at least 4 contiguous planes with clusters near the interaction vertex. An iterative procedure merges these seeds into larger ones based on their relative proximity and angular distributions. A newly made short track is then anchored to the interaction vertex by adding the vertex cluster to the tracks clusters.
The second algorithm performs a scan around the interaction vertex looking for clusters that have the same angular orientations. Most of the short tracks made here consist of contiguous clusters, but some (stub-like, present in just a few planes) can also contain clusters that are not contiguous. The anchoring of the newly made short track to the interaction vertex is based on the proximity of its nearest node to the vertex.

C. Shower Reconstruction
Showers in MINERvA are defined as activity that is correlated in space but is not track-like. The shower reconstruction routine starts by using in time clusters that are not part of tracks. It combines cluster information from all three of the stereoscopic plane orientations and delivers one or many three dimensional shower objects. The routine first utilizes a peak finding algorithm to select the highest energy clusters in each of the three orientations separately, inside the tracker, the downstream electromagnetic calorimeter, and the hadron calorimeter. These are used for initiating the shower building and are referred to as shower seeds. A particular cluster is considered to be adjacent to the shower seed, if it is present in the next detector plane of the same orientation as the shower seed. The cluster is also required to overlap with the shower seed with respect to its transverse positions. If both the criteria are satisfied, the candidate cluster is added on to the growing shower seed. The shower reconstruction routine follows an iterative procedure, whereby newly formed shower seeds are also merged with each other, depending on their conditions of adjacency and overlapping positions.
The seeds made inside each of the orientations are now matched across the different plane orientations, to trans-form them into three dimensional objects. This growing occurs inside each sub detector of the MINERvA detector. Two seeds are said to be adjacent if they are present in the same or neighboring detector modules and overlap in their transverse dimensions. In this way a three dimensional shower seed is formed. An iterative combining procedure merges these shower seeds to form large, three dimensional shower-like objects.
The energy of the shower is estimated by correcting each cluster in the shower for the passive material traversed just upstream of that cluster.

D. Vertex Energy Reconstruction
The vertex energy in this analysis is defined as the untracked energy within a radius of 150 mm from the muon track vertex. The energy in this region is most likely coming from low energy protons which were not tracked, but could also come from other particles as well. Because the current models of neutrino interactions do not model low energy protons accurately, we collect the energy of clusters in all three views in this region to study it, but we do not make any selection in this analysis based on the energy in this region.

E. Residual Energy Reconstruction
The residual energy is simply the scalar sum of the energy in all the unused clusters in all three views: those that are not already part of tracks, showers or within the primary vertex region. The sum can only be made after the track, shower and vertex energy reconstruction has occurred.

F. Recoil Energy Reconstruction
The "recoil" for this analysis is defined as the activity not associated with the muon and proton candidate tracks, which is at least 150 mm away from the primary vertex. The total recoil energy is determined by taking the calorimetric sum of all the showers and the residual energy.

VI. EVENT SELECTION
In order to be sensitive to both low momentum transfer events where the outgoing nucleons are below tracking threshold, and to take advantage of the additional kinematic handles when a final state proton is identified, this analysis selects events with one negatively charged muon and any number of additional tracked particles originating from the muon vertex. Quasielastic-like interactions with Q 2 < ∼ 0.5 GeV/c 2 , often have soft nucleons in the final state which are below the tracking threshold. Final state protons from interactions with higher Q 2 are energetic enough to be tracked. Quasielastic-like candidate events for this analysis therefore fall into the categories of muon-only, and muon plus N>0 protons. The final cross sections reported in this paper combine these samples, but some reconstructed quantities like vertex energy and recoil energy show distinct differences which are of interest.
In order to be considered in this analysis, an event must have a muon track in MINERvA that is matched into the MINOS detector, and whose charge has been identified in MINOS to be negative. The origin of the muon track must be within the fiducial volume of the tracker volume of the MINERvA detector. Once a muon track is found in a slice of time within the beam spill, the selection process removes activity coming from additional tracks which are not associated with the primary interaction. These tracks occur because of the high intensity of the NuMI beam, and are uncorrelated with the muon track that defines the primary interaction. The activity from tracks whose time difference with respect to the muon time is outside the window of -5 to 10 ns are ignored. The recoil and vertex energy algorithms described above ignore this energy. If the track is within the timing window then particle identification cuts (Sec. VI A) are applied to the track.

A. Charged Pion Rejection
All non-muon tracks passing the above timing cut have a proton particle identification algorithm applied to them. This algorithm compares the charge deposition spectrum to Bethe-Bloch expectations to estimate the momentum as well as provide a hypothesis of particle type. This analysis has a cut varying with Q 2 QE which is relaxed as Q 2 QE increases due to the higher probability of secondary interactions between protons and the detector material.
Another way to reject interactions producing charged pions is to search for candidate Michel electrons that result from the pion to muon to electron decay chain. These electrons are produced later in time than the primary neutrino interaction, and are defined as being either within 200 mm of the muon vertex because they originate from very low energy pions, or within 200 mm of the end points of any tracks passing the timing cut applied above. Michel candidates have at most 70 MeV as constrained by the muon mass, and occur at a later time. Figures 1 and 2 show the multiplicity and energy of the Michel candidates respectively, and Fig. 3 shows the time difference between the initial neutrino interaction and the Michel candidates from the data and simulation. Events with one or more Michel candidates are rejected from the signal sample. . Events with exactly one pion in the final state are labeled by the charge of the pion. If more than one pion is present it is labeled "Nπ". The "other" category is predominately from events with kaons in the final state.

B. Multi-pion and DIS Rejection
Because of the broad energy range of the NuMI beamline, the MINERvA detector has access to events with large hadronic recoil mass (W). As a result, a large fraction of events can pass the event selection while having large amounts of recoil energy in the detector. An earlier analysis published in Ref. [5] placed a tighter requirement on the recoil energy. This event selection choice introduced a model dependence on the recoil shape of the 2p2h events, and an efficiency dependence is introduced with this tighter selection which is not desirable. Instead, the selection criteria is modified to remove events with recoil energy greater than 500 MeV and the extra background accepted is controlled via control samples described in Sec. VII. Figure 4 is the untracked recoil for  events passing all selection cuts except the recoil cut.
The clustering algorithm to identify showers described in Sec. V C identifies isolated energy deposits that may come from photons as well as neutrons. No more than one isolated shower is allowed in the event. Events with two or more showers are used in the two sideband samples described later in Section VII A.

C. Final Event Sample
After all event selections are applied the analysis has two populations of events: the muon only sample and the muon plus one or more proton track sample. The resulting event rates for each population as well as the combined sample as shown in Fig. 8 for longitudinal momentum (p || ) and transverse momentum(p t ) versions of the two-dimensional results. In addition, a breakdown of the predicted event rate by type is presented in Fig. 9. The event samples after all cuts but before background subtraction are 62776 and 46499 events for the single and multi-track sample, respectively. The predicted purity, defined as the number of signal events per selected events, of these samples are 79.9% and 56.0%, respectively.
The reconstructed proton efficiency as a function of the highest energy proton and angle with respect to the z-axis of the detector is shown in Fig. 10. This plot demonstrates how often quasielastic-like candidates end up in the multi-track events versus the single track event  category.

VII. CROSS SECTION EXTRACTION
To extract cross sections, the background must first be determined and subtracted. In order to reduce model dependence the backgrounds are constrained using independent data samples. After background subtraction the remaining event samples are unfolded to remove detector resolution effects. After that the samples are corrected for detector efficiency. Finally, the event samples are divided by the flux and number of targets to arrive at cross sections. The following subsections describe each of those steps.

A. Background Constraint
This analysis uses three independent control samples for the pion backgrounds: a Michel only sample, a shower only (> 1 shower) sample, and a shower plus Michel sample. These are designed to characterize the single π + backgrounds (Michel), the neutral pion backgrounds (showers), and the multipion backgrounds (Michel plus showers). Figures 11, 12 and 13 show the input distributions for the sideband constraint.
All three samples are fit in bins of transverse muon momentum with binning compatible with the signal sample, with reduced numbers of bins in cases of small sample size. They are fit simultaneously with a standard χ 2 minimization using ROOTs TMinuit class [31]. Three scaling factors are extracted as a function of transverse momentum for the single track sample as well as the multi-track sample. These are then applied to the signal sample and the resulting overall scaling factor is derived and is shown in Fig. 14.
The sideband corrections reduce the pion content at moderate to low transverse momentum. This is consis- tent with previous MINERvA results which indicate a need for an overall reduction in the pion rate as well as preference for even lower rates at low Q 2 QE .

B. Unfolding
The D'Agostini unfolding method [32][33], via the implementation in RooUnfold [34] is used to remove detector resolution effects from the sample. The number of iterations is chosen based on pseudodata studies which take warped input pseudodata and unfold with a different sample's migration matrix. A wide variety of pseudodata samples were used to investigate the number of iterations necessary to unfold to the true pseudodata distribution. The study used model variations ranging from a GENIE 2.8.4 prediction without 2p2h or other modifications to using MINERvA GENIE tune v1 with an additional modification of resonant pion RPA introduced via a MINOS empirical model [35]. This was done by a χ 2 comparison between the unfolded pseudodata and the sample used to generate the pseudodata. In all cases the preferred number of iterations was no more than four, and if the pseudodata and migration model were similar even fewer iterations were required. This result uses four iterations.
During the study it was discovered that certain pseudodata samples would result in large χ 2 at higher numbers of iterations. This was due to fluctuations of particular p t -p || bins with low statistics. At four iterations the probability of the sample producing a large χ 2 was 1 in 1000. We also analyzed the real data sample, using all model variations listed in Tab. I, to more than a hundred iterations and did not see a χ 2 increase. The combination of the very low probability and the performance with the real data gives us confidence the procedure is robust.
The migration matrices for the p t -p || result and their projections are in Figs. 15, 16, 17, 18, and 19. In general the populations are mostly diagonal.

C. Efficiency Correction
The efficiency correction in the p t -p || space is shown in Fig. 20. The zero efficiency at lower p || for higher p t bins is the effect of the muon angle requirement in the signal definition. Bins where this requirement still accepts a very small component of signal events near bins edges are not reported in the final cross section result due to extremely small efficiencies and event rates. The efficiency in these boundary regions were kept as smooth as resolutions would allow.
In addition to the overall efficiency of the signal, the generator level (QE, 2p2h, resonant) components of the signal efficiency are also presented in Fig. 21, and show that the relative efficiencies for each predicted component are similar. The cutoff at p t of 1.2 GeV is due to the FIG. 14: Combined scaling factor for π ± ,π 0 , and N π backgrounds. Top is the sample with only the muon reconstructed. The bottom is the sample with the muon and proton candidates reconstructed.
restriction of the 2p2h model having a cutoff at q 3 at 1.2 GeV.
The Q 2 QE efficiency is shown in Fig. 22. The reduction at Q 2 QE greater than 0.1 GeV 2 is due to the tracking thresholds and PID performance. The small increase of efficiency at 0.75 GeV 2 is due to the removal of the PID requirements, see Sec. VI A.
The E ν,QE efficiency is shown in Fig. 23. The lower efficiency at lower E ν,QE is due to the MINOS match requirement, but in general the sample is 50% efficient or better. The point of inflection at 8 GeV occurs when the population of tracks move from the forward fine grained region of the MINOS detector to the less instrumented spectrometer downstream.

D. Normalization
The analysis uses the flux described in Ref. [7]. In addition, the each result is scaled by the number of nucleons in the fiducial volume, including hydrogen, with a mix of 88.51% carbon, 8.18% hydrogen, 2.5% oxygen, 0.47% titanium, 0.2% chlorine, 0.07% aluminum, and 0.07% silicon. The number of neutrons in the region is 1.49×10 30 .
The p t -p || and Q 2 QE results are flux averaged results. The normalization factor is the result of the integral of the predicted flux from 0 to 120 GeV. This corresponds to 2.877 × 10 −8 cm −2 per proton on target. For the E ν,QE result the normalization is done in the same manner as Ref. [36]. For each bin of E ν,QE the events are normalized by the integral, in true E ν , of the flux predic-   tion with the same limits. This results in a cross section which is not a true total cross section, but is a welldefined quantity.

VIII. SYSTEMATIC UNCERTAINTIES
The systematic uncertainties for the p t -p || , Q 2 QE , and E ν,QE are shown in Figs. 24, 25, and 26 respectively.
Flux is typically the dominant uncertainty in most of the phase space, at the 8-9% level for flux averaged results and up to 12% in the falling edge of the beam focusing regions and higher neutrino energy regions for non-fluxintegrated results. Another large systematic uncertainty is the muon reconstruction which dominates in particular regions of phase space due to the peaked nature of the muon momentum spectrum.
The flux uncertainty is largest except at very large transverse momentum. The FSI model in particular plays a larger role in this region at the 8% level. This is also a region where backgrounds are larger than signal and the sideband constraint is limited by statistics.
To evaluate the uncertainty of the enhancement to 2p2h described in Sec. III, the application of the 2p2h enhancement is varied in three ways: the enhancement is applied only to the np nucleon pairs in the 2p2h model, the enhancement is applied only to the nn nucleon pairs, and the enhancement is applied instead to the 1p1h model. Each of these variations consistently describes the data in Ref. [28], but attributes it to a different hard scattering process. The measurement of quasielasticlike scattering is largely immune to these variations because the primary effect of these variations is to modify the recoil spectrum in the detector. However, the selection only rejects events with greater than 500 MeV of visible recoil energy, much larger than is typical for 1p1h or 2p2h interaction.
The FSI model uncertainties are evaluated using the  standard GENIE reweighting infrastructure. Due to the correlated nature of the parameters we do not use the uncertainty associated with inelastic pion interactions, which is correlated with the pion absorption uncertainty. The uncertainty associated with the mean free path of pions and inelastic interactions of nucleons starts to contribute at the several per cent level beyond p t >1 GeV.
The interaction model uncertainties are evaluated using the standard GENIE reweighting infrastructure. No source of uncertainty dominates for most of the phase space. The expection is non-resonant backgrounds where an interaction on a proton results in two pions, and two high-twist modifications in the Bodek-Yang model (A,B) which are individually greater than three percent at p t >1.5 GeV.
The systematic uncertainty resulting from the appli-  cation of RPA to the CCQE events is broken into two kinematic regions as described in the Ref. [22] analysis of the model from Ref. [21]. Variations of uncertainties on the input parameters to the relativistic model are used as one source of uncertainty, though the non-relativistic version of the model is used technically to implement this. Above 0.4 GeV 2 , the uncertainty is 60% of the difference between the two on the high side, and no modification on the low side. For Q 2 ≤ 0.4 GeV 2 , the uncertainty is set by the comparison of the model with muon capture data on multiple nuclei [24,37]. This is set at 25% change in the strength of suppression due to the poor description of muon capture on carbon.
The other category includes a large number of systematics including reconstruction efficiencies, calorimetric response to single particles constrained by Ref. [30], inelastic interaction cross sections for pions, protons, neutrons in the detector material, and detector modeling. A subset of these-proton reconstruction efficiency, particle identification by dE dx in the scintillator, and MINERvA -MINOS matching efficiencies-contribute significantly at large p t .

IX. RESULTS
The analysis provides five different cross section results: a double differential cross section in p t -p || with comparisons to MINERvA GENIE tune v1 and additional modifications to the model as well as a comparison to NuWro [38], and four single differential results. In addition, the double differential cross section is projected into Q 2 QE which is highly correlated with p t and E ν,QE which is correlated with p || . An analysis measuring the reconstructed energy near the vertex is presented. Furthermore, an extraction of the double differential cross section with a modified signal definition can be found in Sec. XI.

A. Double Differential Cross Section
The double differential cross section of p t versus p || is shown in Fig. 27. The effect of the 2p2h enhancement and RPA reweighting are clearly visible and improve agreement with the data. The MINERvA GENIE tune v1 agrees with the data except in regions of low and high p t where the data cross section is smaller than predicted.
A comprehensive set of variations of the MIN-ERvA GENIE tune v1 is provided both to demonstrate the effect of each component and to provide guidance on some choices which appear to better simulate the data. Variations applied to the base GENIE model include the application of RPA to CCQE (Fig. 28), the addition of 2p2h (Fig. 29), the combination of applying RPA and adding 2p2h (Fig. 30). In addition, a comparison to a different neutrino generator prediction, NuWro [38] revision 17.09, is shown in Fig. 31. The application of an empirical resonant pion low Q 2 suppression whose form is based on MINOS data [35] is applied as an addition to the MINERvA GENIE tune v1. Table I provides a breakdown of χ 2 for the various model variations. Comparisons of models to data are presented under two calculations: standard and lognormal χ 2 . A pathology referred to as Peelle's-Pertinent-Puzzle(PPP) [39][40] [41] occurs when there are highly correlated systematic uncertainies and a large source of normalization uncertainty. This is observed when the overall scale of the model is well below or above the data.
A standard χ 2 calculation assumes the underlying uncertainty is normally distributed, but this isn't the case for all sources of uncertainty in the cross section extraction. Sources of uncertainty which come from multiplicative operations are log-normal. Background subtraction is a source of additive uncertainty while efficiency corrections, normalization factors, such as flux, are multiplicative. Statistical uncertainties are typically additive, but the unfolding procedure makes this less clear. We present both the linear and log versions of the χ 2 as a demonstration of just how different the best fit is under these two fitting prescriptions.
In the case of linear χ 2 a model with the best fit is clearly not modeling the data as seen in Fig. 28. While the log-normal case is more consistent with our intuition, it still shows disagreement. In the double differential cross section there are two regions of clear data-Monte Carlo difference, which drive the overall poor χ 2 . These are more clearly shown in the single differential cross sections.

B. Single Differential Cross Sections
The one dimensional projections of the p t -p || result are shown in Fig. 32.
A single differential cross section in E ν,QE , which is correlated with p || , is provided in Fig. 33. The quasielastic-like MiniBooNE result [42] is shown for comparison at lower energies. The reader should again keep in mind this result is not strictly an absolute cross section, σ(E ν ), but dσ dEν,QE . Even at a given E ν,QE value, two experiments will measure different cross sections due to different fluxes, detector acceptances, and signal definitions. In particular, the θ µ requirement of the signal definition in this analysis will have a large effect at low E ν,QE , demonstrated by the curve labeled "MnvGENIE v1 -No Angle Cut". A single differential cross section of Q 2 QE is presented in Fig. 34 with the same signal model curves as the double differential result. The addition of the default GENIE 2.8.4 prediction is provided for reference. The effect of the 2p2h enhancement is strongest at moderate Q 2 QE . The prediction disagrees with the data at very low Q 2 QE and high Q 2 QE . Except at the very edges of the Q 2 QE distribution, the rate and shape of the Q 2 QE distribution is well described confirming the anomalous M A evaluated by [43,44] for the same range of Q 2 QE was due to unmodeled multi-nucleon effects not readily available at the time of those analyses.
For Q 2 QE < 0.01 GeV 2 , one-third of the predicted rate comes from resonant pion production which enter the sample via FSI. The initial state model is missing the effect of a suppression at low Q 2 , which the background constraint from Sec. VII indicates is necessary to describe the sidebands. Previous MINERvA results have also measured cross sections less than the prediction in  comparable phase space [28][45] [46] [47]. In addition, MI-NOS measured a similar effect [35].
The region where Q 2 QE > 2 GeV 2 is predicted to be ∼70% true CCQE events. The remaining rate comes from resonant pion events with a small contribution from DIS interactions. The prediction of a larger cross section at high Q 2 QE could be due to any combination of the following: the dipole form of the axial form factor or final state interactions or initial state nucleon kinematics.
A variety of data is available to constrain the axial form factor such as data on deuterium [48] [57] have data for Q 2 > 2 GeV 2 , and they are statistics limited. The z-expansion model [60] [61] provides a different parameterization of the axial form factor, extracted from the deuterium data listed above, suggests that the current data permits significant variation from the dipole model at large Q 2 .

C. Vertex Energy
This analysis ignores energy within 150 mm of the interaction vertex to ensure it is not strongly sensitive to modeling the particle content in this region. However, the MINERvA GENIE tune v1 attempts to provide a more descriptive model of the data and it is important to see how this tune predicts the vertex energy. The prior MINERvA investigation [5] preferred a spectrum with 25% more protons. At the time a 2p2h model was not part of the reference model. The vertex energy distribution within a 150 mm region around the vertex with the MINERvA GENIE tune v1 is shown for both the 1-track and multi-track samples in Fig. 35, respectively. Backgrounds have been constrained in these samples. The ratio plot with respect to default GENIE 2.8.4 with full systematic errors are shown in Fig. 36. It is clear without the 2p2h and RPA modeling the data would be poorly described, but the data also favors the empirical enhancement used in the MINERvA GENIE tune v1.
As described in Sec. VIII the effect of the 2p2h enhancement has a systematic uncertainty derived by three different applications of the fit to various potential contributors, np-pair 2p2h, nn-pair 2p2h, and QE-only. Figure 37 shows the effect of these variations on the vertex energy distribution.
The sample has enough events to further break these vertex energy distributions into bins of p t , shown in Figs. 38 and 39. Regions with noticeable differences between the simulation and data include low p t with large vertex energy. Overall, the single track sample has a χ 2 of 355 per 247 degrees of freedom, while the multi-track sample has a χ 2 of 195 per 104 degrees of freedom, so both samples have significant disagreements with the MIN-ERvA GENIE tune v1.
Events with no second track reconstructed and p t <0.4 GeV 2 show a prediction of more events with large vertex energy than seen in data. This is also seen in the multi-track events at low vertex energy for p t <0.4 GeV 2 . The predicted fraction of the event rate by different signal and background processes is shown in Fig. 40. The regions of Monte Carlo excess correspond to regions of the vertex energy where resonant pion production contributes more to the signal.

X. CONCLUSIONS
This analysis measures a double differential cross section with respect to the longitudinal and transverse momentum of the muon for quasielastic-like events. A suite of various additional processes and models were added to GENIE and are compared to the data.
The MINERvA GENIE tune v1 models the data well except low and high p t . A low p t , and in turn low Q 2 QE , the addition of a low Q 2 suppression to resonant events would better replicate the data. At high p t a plausible explanation of the Monte Carlo data difference is the current model of the axial form factor does not work in this region.
Finally, a detailed look at the energy deposited near the interaction vertex shows very good agreement with the MINERvA GENIE tune v1 for overall vertex energy but deviates when separated into bins of p t . These results are consistent with the previous MINERvA result [5], and demonstrate that the enhancement of 2p2h processes provides a model for such additional low energy protons.

XI. APPENDIX: QUASIELASTIC RESULT
A similar analysis was done with a different signal definition to provide a measurement for predictions which cannot produce a post-FSI signal. The model dependence of this result appears in the cross section modeling and FSI systematic uncertainties which are much larger than the quasielastic-like result. This result applies the same techniques as described in Sec. VII. The signal def-inition for the result described in this section is: true quasielastic events, including 2p2h events, with a muon angle less than 20 degrees.
The p t p || differential cross section is shown in Fig. 44   [GeV] ,QE ν