Measurement of inclusive charged-current $\nu_\mu$ cross sections as a function of muon kinematics at $\sim6~GeV$ on hydrocarbon

MINERvA presents a new analysis of inclusive charged-current neutrino interactions on a hydrocarbon target. We report single and double-differential cross sections in muon transverse and longitudinal momentum. These measurements are compared to neutrino interaction generator predictions from GENIE, NuWro, GiBUU, and NEUT. In addition, comparisons against models with different treatments of multi-nucleon correlations, nuclear effects, resonant pion production, and deep inelastic scattering are presented. The data recorded corresponds to $10.61\times10^{20}$ protons on target with a peak neutrino energy of approximately 6 GeV. The higher energy and larger statistics of these data extend the kinematic range for model testing beyond previous MINERvA inclusive charged-current measurements. The results are not well modeled by several generator predictions using a variety of input models.


I. INTRODUCTION
Neutrino oscillation experiments [1][2][3][4] depend on neutrino interaction models to correct for detector and nuclear effects. Oscillation experiments at a few GeV of mean neutrino energy plan to use an inclusive chargedcurrent (CC) signal to maximize far detector statistical * Now at Brookhaven National Laboratory † Now at University of Minnesota precision. An important component of these measurements is the identification of the resulting lepton. Accurate prediction of the momentum and angular distributions of the lepton are required to correct the measured rate for efficiency and acceptance in both near and far detectors. Models of neutrino interactions are also used as input to neutrino energy reconstruction; mismodeling of lepton energy is prima facie evidence that neutrino energy reconstruction will be similarly flawed when using that neutrino interaction model as an input. Inclusive cross section measurements have been made on a variety of nuclear targets in the past. Micro-BooNE [5] and T2K [6] have measured double-differential cross sections as a function of muon momentum and angle on argon and hydrocarbon, but at a lower mean neutrino energy than MINERvA [7]. NOMAD [8] as well as MINOS [9] and CCFR [10] made measurements as a function of neutrino energy on carbon and iron respectively. MINERvA has made measurements as a function of neutrino energy using the low-ν method on carbon for both neutrino and anti-neutrino beams [11,12] and as a function of muon transverse and longitudinal momentum in the Low Energy (LE) NuMI beam with a neutrino flux peaked at 3 GeV [13]. The result presented here increases the phase space accepted into the multi-GeV regime and as a result expands the range of transverse momentum from 2.5 to 4.5 GeV and longitudinal momentum from 20 to 60 GeV with a ∼12 times larger sample size and flux normalization uncertainty of approximately 1/2 the size of the previous result. We present here the two-dimensional cross section for the inclusive neutrino scattering as a function of the muon transverse (p t ) and longitudinal momentum (p || ) in the Medium Energy (ME) NuMI beam, which has a neutrino flux peaked near 6 GeV. Figure 1 compares the Low and Medium Energy fluxes used by MINERvA. The muon momentum and angle can be precisely measured. These muon variables are suitable for comparison to exclusive channel measurements and provide a foundation to understand how model predictions combine to form an inclusive cross section prediction. In addition to the twodimensional cross sections, one-dimensional projections, limited to the phase space of the double-differential cross section, are also provided.
Section II describes the experimental setup. Section III describes the simulation of the neutrino interactions, the modifications made to the interaction model, and the simulation of particle propagation through the detector.
The event selection and measurement methods used to extract the differential cross sections are described in Section IV. A description of the sources and determination of systematic uncertainties are presented in Section V. Section VI describes the cross section results while Section VII provides a set of comparisons to multiple neutrino generator predictions as well as modifications to these predictions. Finally, Section VIII provides conclusions that can be drawn from these comparisons.

II. EXPERIMENT
The MINERvA experiment employs a fine-grained tracking detector for recording neutrino interactions produced by the NuMI beamline at Fermilab [14,15]. Neutrinos are created by directing 120 GeV protons from the Main Injector onto a graphite target. The resulting charged pions and kaons are focused by two magnetic horns. A neutrino-dominated or anti-neutrinodominated beam is produced by switching the polarity of the horns. This analysis uses data from neutrinodominated beam.
The MINERvA detector [7] consists of 120 hexagonal modules that create an active tracking volume preceded by a set of passive nuclear targets. This result includes only those interactions in the active tracking volume with a fiducial mass of 5.48 tons. The active target volume is surrounded by electromagnetic and hadronic calorimeters.
Each tracking module is made of two planes. Each plane is comprised of triangular polystyrene scintillator strips with a 1.7 cm strip-to-strip pitch. To allow for better three-dimensional reconstruction in a highmultiplicity environment, planes are oriented in three different directions, 0 • and ± 60 • relative to the vertical axis of the detector. The downstream and side electromagnetic calorimeter consists of alternating layers of scintillator and 2 mm thick lead planes. The downstream and side hadronic calorimeters consists of alternating scintillator and 2.54 cm thick steel planes.
Multi-anode photomultiplier tubes read out the scintillator strips via wavelength-shifting fibers. The timing resolution measured by thoroughgoing muons is 3.0 ns and sufficient to separate multiple interactions within a single NuMI beam spill. Muons that originate in MIN-ERvA are analyzed by the MINOS near detector [16], a magnetized spectrometer composed of scintillator and iron and located 2 m downstream of the MINERvA detector. The requirement that muons are analyzed in MINOS restricts this analysis to muons with p || > 1.5 GeV/c and θ µ < 20 • , which means a restricted acceptance for events 8 . This analysis uses data that correspond to 10.61 × 10 20 protons on target (POT), received between September 2013 and February 2017 while the horn polarity was set to focus positively charged particles, creating a beam that is predominantly muon neutrinos.

III. SIMULATION
A GEANT4-based simulation of the NuMI beamline is used to predict the neutrino flux. To improve the prediction, the simulation is reweighted as a function of pion kinematics to correct for differences between the GEANT4 [17] 1 prediction and hadron production measurements of 158 GeV protons on carbon from the NA49 experiment [18] and other relevant hadron production measurements. A description of this procedure is found in Ref. [15]. In addition, an in situ measurement of neutrino scattering off atomic electrons is used, as described in Ref. [19], to constrain the flux prediction.
Neutrino interactions are simulated using the GENIE neutrino event generator [20] version 2.12.6. Quasielastic (1p1h) interactions are simulated using the Llewellyn-Smith formalism [21] with the vector form factors modeled using the BBBA05 model [22]. The axial vector form factor uses the dipole form with an axial mass of M A = 0.99 GeV/c 2 . Resonance production is simulated using the Rein-Sehgal model [23] with an axial mass of M RES A = 1.12 GeV/c 2 . Higher invariant mass interactions are simulated using a leading order model for deep inelastic scattering (DIS) with the Bodek-Yang prescription [24] for the modification at low momentum transfer squared, Q 2 .
A relativistic Fermi gas model [25] is used with an additional Bodek-Ritchie high momentum tail [26] to account for nucleon-nucleon short range correlations. 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 [27].
To better describe MINERvA data, a variety of modifications to the interaction model are made. To better simulate quasielastic events, the cross section is modified as a function of energy and three momentum transfer based on the random phase approximation (RPA) part of the Valencia model [28,29] appropriate for a Fermi gas [30,31]. Multi-nucleon scattering (two-particle two-hole or "2p2h") is simulated by the same Valencia model [32][33][34], but the cross section is increased in specific regions of energy and three momentum transfer based on fits to MINERvA data [35] in a lower energy beam configuration. Integrated over all phase space, the rate of 2p2h is increased by 50% over the nominal prediction. Based on fits done in Ref. [36], we decrease the non-resonant pion production by 43% and reduce the uncertainty compared to the base GENIE model uncertainties. This modified version of the simulation is referred to as MINERvA Tune v1.
The response of the MINERvA detector is simulated using GEANT4 [17] version 4.9.3.p6 with the QGSP BERT physics list. The optical and electronics performance is also simulated. Through-going muons are used to set the absolute energy scale of minimum ionizing energy depositions by requiring the average and RMS of energy deposits match between data and simulation as a function of time. A full description is found in Ref. [7]. Measurements using a charged particle test beam [37] and a scaled-down version of the MINERvA detector set the absolute energy response to charged hadrons. The effects of accidental activity are simulated by overlaying hits in both MINERvA and MINOS from data corresponding to random beam spills appropriate to the time periods in the simulation.

IV. CROSS SECTION EXTRACTION
A sample of neutrino charged-current interactions is extracted by requiring the track identified as being from a muon to be matched between MINERvA and MINOS, and to be negatively charged. In addition, the reconstructed interaction vertex must be within a specified fiducial volume. To avoid model dependence introduced by correcting for kinematic regions without acceptance, we only report results for charged-current cross section where the muon angle with respect to the neutrino direction is less than 20 • , the muon p t is less than 4.5 GeV, and the muon p || is between 1.5 GeV and 60 GeV.
Using these criteria, a sample of 4,105,696 interactions was selected. The simulation predicts an average selection efficiency of 64% in the p t -p || phase space, where the efficiency loss is due to the MINERvA-MINOS geometric acceptance. After all selection cuts, the sample in muon transverse and longitudinal momentum space is shown in Fig. 2, decomposed into predicted components. Events are labeled by categories within GENIE except for events given a DIS label. To explore how contributions from DIS events rely on the validity of the neutrino-quark scattering model with different "depth" of the inelasticity, DIS is divided into two categories. "True DIS" events are those events where the invariant mass of the hadronic system, W , is greater than 2.0 GeV/c 2 and Q 2 greater than 1.0 GeV 2 /c 4 . "Soft DIS" represents the remainder of the GENIE DIS events. While the cut defining the 'True DIS' regime is defined by Q 2 and hadronic invariant mass W values, and could therefore be used for model comparisons, the 'Soft DIS' definition is not. The modeling of inelastic events below the 'True DIS' region can vary widely across generators, in terms of the kinematic coverage of the resonance model and handling of non-resonant contributions in the resonance region [38]. No two generators handle these aspects in exactly the same way, so the 'Soft DIS' label here is relevant only to GENIE simulations. A final category is "other CC", which contains CC events not belonging to the other categories, such as coherent charged pion production. The background category con-   tains charged-current events from other neutrino flavors and anti-muon type neutrinos as well as neutral-current interactions. A total of 8655 (0.2%) background events are predicted. Backgrounds at p || less than 2.5 GeV and p t less than 0.4 GeV are primarily from neutral-current interactions where a pion was reconstructed as a muon in MINOS. These backgrounds have a maximum contribution of 10 percent of the predicted event rate at these low p || and p t , and are typically much smaller. Backgrounds at high p t are mostly anti-neutrino contamination due to muon charge misidentification which accounts for about one percent of the sample in the highest p t bin.
The predicted background contributions are subtracted from the sample. Detector resolution effects (see Figs. 4 and 3) are then removed using the D'Agostini unfolding method [39,40], via the implementation in RooUnfold [41]. To understand the necessary regularization strength, 10 different model predictions as pseudodata were unfolded using the MINERvA Tune v1. These fake data were derived by reweighting the 2p2h strength, QE RPA, non-resonsant pion production reweight, resonant pion production. Many of the models used appear in Table I. The unfolded models were then compared to their true distributions via a χ 2 test taking full consideration of correlations. The optimal number of iterations was determined when the χ 2 approached one per degree of freedom and was not changing as a function of the number of iterations. Only statistical uncertainties were considered in determining the number of iterations. In all variations the required number of iterations was no more than 10. In addition, a fit to the data in reconstructed p t -p || was performed, and the MINERvA Tune v1 prediction reweighted to the data as an additional fake data sample.A fit is performed as a function of p t -p || to provide a reweight value forcing the MINERvA Tune v1 to better agree with the data. This reweighted prediction is used as another fake data sample. Reweighting is done in true kinematic quantities, although the fit is done with reconstructed quantities, and propagated through the Monte Carlo detector response prediction. Based on these studies, the data was unfolded using 10 iterations.
Finally, the sample is corrected for efficiency and acceptance. The selection efficiency is shown in Fig. 5. The large efficiency in the 6-7 GeV p || and highest p t bin is due to a fractionally large sample of muons with generated angle greater than 20 degrees passing event selection and appearing in this bin.
The efficiency-corrected distribution is then divided by the integral of the flux with neutrino energies between 0 and 100 GeV averaged over the fiducial volume, which is 6.32 × 10 −8 ± 3.9% per cm 2 per proton on target, and the number of nucleons in the fiducial volume, 3.23 × 10 30 ±1.4%, with a mass fraction of 88.51% carbon, 8.18% hydrogen, 2.5% oxygen, 0.47% titanium, 0.2% chlorine, 0.07% aluminum, and 0.07% silicon.

V. SYSTEMATIC UNCERTAINTIES
Systematic uncertainties in this analysis fall under three different categories: flux, detector response, and neutrino interaction model uncertainties. The uncertainties from individual sources are evaluated by reextracting the cross section using modified simulations. The size of each modification is related to the uncertainty in each source. Flux uncertainty, a typical leading uncertainty in neutrino cross section measurements, is below 4% for almost all the phase space because of the flux constraint that comes from a measurement of neutrino-electron scattering in the same beam [19]. The normalization uncertainty of 1.4% corresponds to the uncertainty in the number of target nucleons and is based on material assays and weight measurements of productionquality scintillator planes.
Uncertainty in the detector response to hadrons is evaluated using shifts determined by in situ measurements of a smaller version of the detector in a test beam [37]. Uncertainties in inelastic interaction cross sections for particles in the detector material are independently varied based on data-Monte Carlo differences between GEANT particle cross sections and world data on neutrons [42][43][44][45], pions [46][47][48][49], and protons [50][51][52].
Muon reconstruction uncertainty is dominated by the muon energy scale uncertainty, which is constrained by a fit to the reconstructed neutrino energy distribution for low recoil neutrino charged-current events, whose cross section is known to be flat as a function of neutrino energy. In the low recoil fit procedure, we include the model uncertainties that are used in this and other MIN-ERvA results, including uncertainty on the 2p2h process informed by MINERvA data [35]. Because the low recoil cut is at 800 MeV, this results in a small uncertainty in the fit due to cross section modeling [53].The resulting uncertainty on the muon energy scale is 1%. The fit causes the flux and muon energy scale uncertainties to be correlated and that correlation is propagated to the final result. Uncertainty in the matching efficiency is from imperfect modeling of the efficiency loss from accidental activity in the MINOS near detector when matching muon tracks from MINERvA to MINOS. This last efficiency is also determined by a data-simulation comparison as a function of instantaneous neutrino beam intensity.
Interaction model uncertainties are evaluated using the standard GENIE reweighting infrastructure [20,54]. Because this is an inclusive analysis with very low backgrounds and few selection cuts, model uncertainties are never the dominant uncertainty in any p t -p || bin. These uncertainties are most significant at the highest p t bins where the geometric acceptance changes dramatically and at low p t bins where the backgrounds are the largest. The reconstruction of the muon is largely unaffected by the hadronic shower. The efficiency of the selection as a function of p t -p || -hadronic energy was evaluated and found to be flat as a function of hadronic energy.
The fractional uncertainties in the one-dimensional projections are shown in Figs. 6 and 7. The fractional uncertainties in the two-dimensional result are shown in Fig. 8. The dominant uncertainties are the muon energy scale uncertainty and the flux normalization. The muon energy scale uncertainty has the largest effect on the cross section measurement at the rising and falling edges of the peak of the muon momentum spectrum, where the slope between bins is largest. The muon momentum peaks at approximately p || = 5 GeV and p t = 0.6 GeV. It should be noted there are particular regions where the pion re-interaction probability uncertainty, grouped in the hadronic response systematic, is large. This is due to interactions in which the reconstructed muon was not the primary muon, but was instead due to a high energy pion. This is also how a population of neutral current interactions populate the lowest p || bins at low p t .

VI. RESULTS
Three results are presented: two single-differential cross sections and a double-differential cross section using the p t and p || of the muon. The single-differential cross sections are shown in Figs. 9 and 10. The ratio of data to MINERvA Tune v1 for the single-differential results is shown in Figs. 11 and 12. The double-differential cross section as a function of p t and p || is shown in Fig. 13. The single-differential cross sections are derived from the two-dimensional result which means the additional phase space restrictions of the two-dimensional spaces are incorporated. The cross section versus p t includes a restriction of 1.5 ≤ p || ≤ 60 GeV/c, while the cross section versus p || includes a restriction of p t ≤ 4.5 GeV/c. The difference between data and MIN-ERvA Tune v1, shown in Fig.  11, is due to a mismodeling of the cross section as function of muon kinematics combined with the angular acceptance and muon energy scale. A study was performed by correcting the prediction to match the data in the p t projection. The p || prediction after this correction was compared to the data and was found to be consistent in both normalization and shape within the muon energy scale uncertainty. Figure 14 shows the ratio of data to simulation. For p || between 3 and 15 GeV and low values of p t , the cross section is overpredicted. In this region the dominant process is resonant pion production which has been previously measured by MINERvA [55][56][57][58][59][60], MiniBooNE [61], and T2K [62]. Many of these measurements, including a MINOS measurement [63], indicate the need to reduce the predicted cross section at low Q 2 which corresponds here to regions of low p t . A set of comparisons against a variety of resonant pion production model modifications is shown in Sec. VII. At p t > 0.85 GeV/c the Monte Carlo prediction consistently underpredicts the data by 10-25%. This high p t region is dominated by the "True DIS" process for p || > 6 GeV/c. Poorly understood neutrino DIS nuclear effects could contribute to the underprediction in this region of kinematics.

VII. COMPARISONS
In this section the extracted cross sections are compared to a variety of predictions. Four different groups of models and the data are presented as a ratio with respect to MINERvA Tune v1. Each group investigates a different aspect of the MINERvA Tune v1. No model combination completely describes the data. We also compare to a group of predictions from different neutrino event generators.
The first of these groups of comparisons, see Fig. 15, considers alterations to MINERvA Tune v1 by adding or subtracting selected changes that were made to the original default GENIE prediction. The cases plotted are: (a) GENIE v2.12.6 with Valencia 2p2h [33,34,64,65], (GENIE 2.12.6); (b) Case a) but with a reduction in the non-resonant pion production of 43% [36], removed and a reweighting of resonant pion production to the MK model [69], (MK Model); (c) Case a) but with a reweighting to reduce the resonant pion cross section at low Q 2 according to a MI-NOS parameterization [63], (Pion LowQ2-MINOS); (d) Case a) but with a reweighting to reduce the resonant pion cross section at low Q 2 according to a MINERvA parameterization [70], (MINERvA Tune v2). Figure 18 shows a comparison between various neutrino generator predictions including GiBUU [71,72] for two different versions, two different nuclear models from NuWro [73,74], and NEUT [75]. 2 Table I gives the χ 2 statistics comparing data to the various model predictions listed above. The χ 2 is calculated using Eq. 1 and summing over i, j which takes into account the covariance between bins. The standard χ 2 calculation assumes the underlying uncertainty is normally distributed. This assumption is not correct as some normally distributed uncertainties manifest as lognormal when propagated to the final result. An example is the flux normalization which introduces uncertainty in the measurement by division. To understand the span of these differences we report both the standard χ 2 and a log-normal version. The model with the lowest standard and log-normal χ 2 is the NuWro prediction with a local Fermi gas (LFG) nuclear model, which like the NEUT local Fermi gas prediction, does well at predicting the data between 4 <p || < 8 GeV and p t > 1.5 GeV . The effect of Peele's Pertinent Puzzle [77][78][79] is clearly shown by the effect of the GiBUU v2019 prediction (or the difference in the peak region between MINERvA Tune v1 and GE-NIE v2.12.6 with Valencia 2p2h in Figs. 10 and 9) which, by eye, has a different normalization than the data. This effect occurs when the dominant uncertainty of a result is a highly correlated normalization uncertainty. The lognormal χ 2 increases significantly for models where the dominant difference with data is the normalization. Because of the flux is a dominant uncertainty in this analysis the log-normal χ 2 is a better estimation, but the ordering of models within each group of model modifications is the same using either estimator.
The MINERvA Tune v1 agrees better with the data than GENIE v2.12.6 with Valencia 2p2h. The inclusion or removal of components can either improve or degrade the data Monte Carlo agreement. The inclusion of the non-resonant pion reduction results in poorer agreement than the base model. GENIE v2.12.6 with Valencia 2p2h and QE RPA and a non-resonant pion production reduction is the best combination to predict the data while the model with the low recoil enhancement has the largest χ 2 . The MINERvA Tune v1 is supported by a variety of exclusive measurements in a 3 GeV neutrino focused beam [36,80,81] indicate a need for the all modifications but when compared against data using the 6 GeV beam [82] GENIE v2.12.6 with Valencia 2p2h and QE RPA and a non-resonant pion production reduction is less discrepant. The modifications to DIS have mild modifications to the overall prediction. Models which further modify the resonant pion cross section improve the prediction. The data prefer a low Q 2 type suppression for resonant pions. To understand the effect of the modifications in detail the χ 2 is broken down into contributions for each kinematic bin. Figure 19 shows the bin-by-bin contributions to the overall χ 2 . The best prediction, based on overall χ 2 , from each model group from Sec. VII is compared against the MINERvA Tune v1. The metric used is the difference in χ 2 on a bin-by-bin basis. The value for the i th bin is the result of the calculation shown in Eqs. 1 -2 where x is the cross section and V is the measurement covariance matrix. Due to the anti-correlations introduced by the unfolding procedure some large bin-to-bin anticorrelations will appear in this metric. Negative values indicate an improvement due to the model with respect to MINERvA Tune v1. GENIE 2.12.6 with Valencia 2p2h and QE RPA and a non-resonant pion reduction, MINERvA Tune v1 with low Q 2 pion suppression models, and NuWro with a local Fermi gas roughly improve in the same regions of phase space. The low Q 2 pion suppression modification appears to be an overcorrection for the lowest p t -p || bins, which is evident in Fig. 17. The modification to the DIS model is different than the other three modifications. The nCTEQ15 model improves agreement in regions of increasing p t as a function of p || . This is a kinematic boundary between the "Soft DIS" and "True DIS".

VIII. CONCLUSIONS
This paper presents the double-differential and singledifferential cross sections as a function of the muon transverse and longitudinal momenta using data recorded by the MINERvA detector in the NuMI beamline. The high statistics and high neutrino energy of these new data, combined with a low flux uncertainty compared to previous measurements [13] mean that new regions of kinematic phase space can be examined to unprecedented precision. The data are compared to model predictions that reweight different components in the GENIE prediction or to external generator predictions. Some modifications to the GENIE prediction are inspired by measurements of previous exclusive or restricted phase space measurements on earlier data taken by MINERvA. Other modifications to the GENIE prediction (MK, AMU, nCTEQ15, and nCTEQν) represent replacements of a particular set of interaction channels. Finally, generators other than GENIE provide a different set of nuclear models, particle transport, and interaction channel models. None of these predictions describe the data well based on χ 2 tests. Similar models were least discrepant when compared with the inclusive double-differential cross section measured [13] in the 3 GeV neutrino focused beam.
The single-and double-differential measurements provide indication of the need for a low Q 2 suppression (low p t ) for the resonant pion production channel. In addition, most of the components in MINERvA Tune v1 are favored, while the low recoil enhancement from increased 2p2h production is disfavored. Comparing the bin-by-bin χ 2 , in Fig. 19, the application of the least discrepant DIS model (nCTEQ15) indicates large changes in the region between "Soft DIS" and "True DIS". Overall, NuWro with a local Fermi gas nuclear model best describes the data. The result cannot differentiate the specific source of mismodeling in regions where all the underlying processes contribute to the prediction. Other methods, either via exclusive, semi-inclusive, or inclusive measurements using other kinematic variables are needed to investigate these complex regions. This work represents an important benchmark that can be used to validate future ensembles of models tuned to agree with exclusive results.