Testing CCQE and 2p2h models in the Neut neutrino interaction generator with published datasets from the MiniBooNE and MINERA experiments
Abstract
The MiniBooNE large axial mass anomaly has prompted a great deal of theoretical work on sophisticated Charged Current QuasiElastic (CCQE) neutrino interaction models in recent years. As the dominant interaction mode at T2K energies, and the signal process in oscillation analyses, it is important for the T2K experiment to include realistic CCQE cross section uncertainties in T2K analyses. To this end, T2K’s Neutrino Interaction Working Group has implemented a number of recent models in NEUT, T2K’s primary neutrino interaction event generator. In this paper, we give an overview of the models implemented, and present fits to published and CCQE cross section measurements from the MiniBooNE and MINERA experiments. The results of the fits are used to select a default cross section model for future T2K analyses, and to constrain the cross section uncertainties of the model. We find a model consisting of a modified relativistic Fermi gas model and multinucleon interactions most consistently describes the available data.
Present address: ]Columbia University, Physics Department, New York, New York 10027, USA Present address: ]University College London, Department of Physics and Astronomy, London, United Kingdom Present address: ]University of Manchester, School of Physics and Astronomy, Manchester, United Kingdom
I Introduction
Charged Current QuasiElastic (CCQE) scattering () is the dominant neutrino interaction process for muon (anti)neutrinos impinging on a nuclear target at neutrino energies on the order of 1 GeV. Because CCQE is a twobody process and the incoming neutrino direction is known for an accelerator experiment, a reasonable approximation of the neutrino energy can be calculated using only the outgoing lepton kinematics. Because of this, CCQE is the preferred signal process for neutrino oscillation experiments which generally require some handle on the incoming neutrino energy to extract neutrino oscillation parameters due to disappearance or appearance in this energy region. However, nuclear effects and interactions which are not distinguishable from CCQE in the final state bias or smear the reconstructed neutrino energy, so a good understanding of these effects is important.
Neutrino interaction generators typically use the relativistic Fermi gas (RFG) model of the nucleus for all neutrinonucleus interactions because of its simplicity. In the RFG model, all possible nucleon momentum states are filled up to the Fermi momentum, there is a constant binding energy required to separate the nucleon from the nucleus, and the neutrino interacts with a single bound nucleon. Neutrinonucleon CCQE scattering for free nucleons is described by the LlewellynSmith formalism Llewellyn Smith (1972), which has been extended to cover neutrinonucleus CCQE scattering in the SmithMoniz RFG model Moniz and Smith (1972), where nucleons bound within the nucleus are described by the RFG nuclear model. The only parameter of the weak current or in the RFG model which is not well constrained by electron scattering data Gallagher et al. (2011); Moniz et al. (1971) is the axial mass, . Results from a global analysis of neutrinodeuterium scattering experiments and pion electroproduction data find GeV/ Budd et al. (2003), which is consistent with other analyses Bernard et al. (2002); Kuzmin et al. (2006, 2008). These results are also consistent with high energy neutrino beam experiments on heavy nuclear targets Lyubushkin et al. (2009).
Recent differential CCQE cross section results from the MiniBooNE collaboration AguilarArevalo et al. (2010, 2013) are significantly higher than expectation, which can only be accounted for in the framework of the SmithMoniz RFG model by inflating the axial mass, giving rise to the term “MiniBooNE large axial mass anomaly”. This came after an earlier large axial mass measurement by K2K Gran et al. (2006), which reported a value of GeV/. Both experiments exhibited not only a largerthanexpected axial mass, but also a supression of low events relative to the expection from the SmithMoniz RFG model. Other experiments using heavy nuclear targets with beam energies in the fewGeV region have also measured cross sections which are consistent with an inflation of the axial mass Abe et al. (2015a, b); Dorman (2009); Wascko (2008), although these results do not paint a coherent picture. More recently, the MINERA experiment Fiorentini et al. (2013); Fields et al. (2013), which is at a somewhat higher neutrino energy than MiniBooNE, has shown good agreement with the SmithMoniz RFG model with = 1.00 GeV/, but requires an enhancement to the transverse component of the cross section, an effect also seen in electronnucleus scattering Bodek et al. (2011). These inconsistent results present a considerable challenge to neutrino oscillation experiments which need to be able to model their signal processes well.
Recent theoretical efforts which have attempted to resolve the “MiniBooNE large axial mass anomaly” have focussed on two main areas: more sophisticated descriptions of the initial ground state of the nucleus; and additional nuclear effects, such as multinucleon interaction models, which go beyond interactions with a single nucleon within the nucleus. The combination of these models would allow for a consistent picture of an axial mass close to 1.00 GeV/, with a suppressed cross section at low and larger cross section at higher relative to a simple RFG model. Comprehensive reviews of available CCQE cross section models can be found in References Gallagher et al. (2011); Morfin et al. (2012); Garvey et al. (2015); AlvarezRuso et al. (2014).
More sophisticated descriptions of the initial state of the nucleus than the RFG model provides are available from a number of authors Benhar and Fabrocini (2000); Ankowski and Sobczyk (2006); Butkevich (2009); Bodek et al. (2014). These models, generically referred to as Spectral Functions (SF), have a more realistic nucleon momentum distribution taking into account the shell structure of the nucleus and correlated pairs of nucleons within the nucleus, and have nonconstant binding energies. Note that although these models include correlations between nucleons in the initial state, they still use the impulse approximation and only consider interactions with a single nucleon. More complex models which go beyond the simple picture of noninteracting fermions are available Leitner et al. (2009); Maieron et al. (2003); Meucci et al. (2004); Lovato et al. (2013); Pandey et al. (2015). However, with the exception of the GiBUU interaction model Buss et al. (2012); Leitner et al. (2009), these are not currently implemented in neutrino interaction generators. In these models, a mean field potential due to the presence of other nucleons within the nucleus is calculated, which will generally depend on the position and momentum of the struck nucleon. These models are not discussed further as they cannot be easily implemented in the NEUT neutrino interaction generator Hayato (2009).
Although alternative nuclear models modify the cross section as a function of the outgoing lepton kinematics significantly, they do not change the total CCQE cross section significantly as a function of neutrino energy Garvey et al. (2015). Additional nuclear effects are also likely to be required to explain the current global dataset. Multinucleon interaction (2p2h) models such as those from Nieves et al. Nieves et al. (2011); Gran et al. (2013) and Martini et al. Martini et al. (2009) go beyond the impulse approximation and include diagrams where two nucleons are involved in the interaction. This adds significant strength to the CCQElike cross section and explains the difference in normalization observed in the MiniBooNE data, which was previously modelled with a large axial mass Nieves et al. (2012); Martini et al. (2011); Nieves et al. (2013); Martini and Ericson (2013). Because these 2p2h models are not twobody processes, they are expected to lead to significant biases in the neutrino energy reconstruction from the outgoing lepton which assumes CCQE kinematics Coloma and Huber (2013); Ankowski et al. (2015). Additionally, the Random Phase Approximation (RPA) is a nuclear screening effect that modifies the propagator for interactions in nuclear matter Jachowicz et al. (1999, 2002); Nieves et al. (2011); Martini et al. (2009) and has a significant effect on the differential cross section as a function of , suppressing the cross section in the low region and enhancing the cross section for GeV. RPA needs to be included, both in interactions with a single nucleon (1p1h) and those from 2p2h calculations, to find good agreement with data. Note that both Nieves and Martini calculations are performed in the context of a local Fermi gas (LFG) model, where the Fermi momentum depends on the local nuclear density, so improvements to the initial state models of the nucleus and improvements to the CCQE interaction models cannot necessarily be combined easily.
Whilst there have been rapid experimental and theoretical developments relating to CCQE cross sections, new nuclear models and nuclear effects have only recently been implemented into neutrino interaction generators or confronted with neutrinonucleus scattering data, and no consistent picture has yet emerged. It is not clear which models fit the global data best, and where the deficiencies now lie, which should be a serious concern for neutrino oscillation experiments. This paper shows the effect of fitting current CCQE and multinucleon models to the MiniBooNE AguilarArevalo et al. (2010, 2013) and MINERA Fiorentini et al. (2013); Fields et al. (2013) datasets to a variety of models implemented in NEUT by members of T2K’s Neutrino Interactions Working Group (NIWG). Previous constraints on the CCQE model produced by the NIWG and used in T2K oscillation analyses only considered an RFG model, and recommended the NEUT default central value for the axial mass = 1.21 GeV based on the value found by the K2K experiment Gran et al. (2006), with an error large enough to cover fits to the MiniBooNE neutrino mode CCQE dataset AguilarArevalo et al. (2010), as is fully described in Reference Abe et al. (2015c). This work improves on the previous situation by including more sophisticated effects proposed to explain the large axial mass anomaly, and by using all of the newly available CCQE data to constrain all model parameters without reference to the default NEUT model.
The models which have been implemented in the NEUT generator are discussed in Section II and Section III discusses crossgenerator validation. Section IV gives a brief overview of the MiniBooNE and MINERA data used in the fit. The nominal NEUT predictions for these datasets are shown in Section V for a variety of models. Section VI discusses the fit procedure. The results of fake data studies and the fit to external data are given in Section VII. In Section VIII we interpret the results and discuss possible implications in cross section and neutrino oscillation analyses and Section IX summarizes the results.
Ii Interaction models in Neut
The motivation for, and an overview of, new CCQE models has already been discussed. This section will briefly outline the important technical details of the models as implemented in NEUT, and highlight any caveats that should be borne in mind when fitting with them. The models used in the fits include the SF model, multinucleon–neutrino interactions, and RPA.
The NEUT implementation of the SF model from Omar Benhar and collaborators Benhar and Fabrocini (2000) is described fully in Reference Furmanski (2015). Although SF is a generic term, in this work it will specifically refer to the Benhar SF. The model information is all encoded in the initial state nucleon distribution shown in Figure 1. Pauli blocking is implemented as a hard cutoff: final state nucleons with momenta less than the Fermi momentum are forbidden. There are two terms in the SF model: a short range correlation term, which extends to higher initial state nucleon momenta, and a mean field term, which contributes the main peak at lower momenta. These terms can be seen in Figure 2, where the twodimensional SF in terms of the removal energy and initial state nucleon momentum has been projected on to the momentum axis. There are three parameters in NEUT which modify the SF as illustrated in Figure 2. The default values for these parameters are given in Table 1. The mean field width and normalization of the correlation term are wellconstrained by electron–scattering data Furmanski (2015) and have little effect on the shape or normalization of the cross section. Thus, they are not considered further in this work. Pauli blocking is modified by changing the Fermi momentum in the fits. It should be noted that in the RFG model, the Fermi momentum defines the Pauli blocking, but also modifies the width of the initial state nucleon distribution. As a result, changing affects a wide range of , whereas changing only affects very low events.
The multinucleon–neutrino (2p2h) model from Nieves et al. Nieves et al. (2011); Gran et al. (2013) has been implemented in NEUT as described in Reference Sinclair (2014). The cross section as a function of neutrino energy and the outgoing lepton kinematics was made available by the authors of the model and is implemented as a series of lookup tables for various nuclear targets and neutrino species. The tables provided had hadronic variables integrated out, so a generic model based on Reference Sobczyk (2012) for simulating the initial and final hadronic states was used for generating NEUT events^{1}^{1}1This model simply enforces energy and momentum conservation, treats initial nucleons as uncorrelated and drawn from a local Fermi gas model, and shares momentum equally between final state nucleons Sobczyk (2012).. The discrepancy between the leptonic and hadronic simulation makes the current NEUT implementation of the Nieves model inadequate for comparisons with experimental measurements of the final state hadrons from CCQE events (such as can be found in Reference Walton et al. (2015)). For this reason, only leptonic measurements are used in this analysis. As the Nieves model is very complex, the current NEUT implementation does not allow fundamental model parameters to be changed. For simplicity, only a simple scaling parameter which changes the normalization of 2p2h events has been considered in this analysis. Note that the Nieves 2p2h model included less decay contributions, where a nucleon excited into a resonance decays without producing a pion Oset and Salcedo (1987); Singh et al. (1998). Contributions from less decay were previously implemented in NEUT and other generators, and have been treated as an intrinsic background in CCQE selections and corrected for. This leads to complications when comparing the full Nieves model to CCQE cross section measurements.
RPA Nieves et al. (2011) is implemented into NEUT as a modification to the 1p1h cross section as a function of and . Figure 3 shows the ratio of the Nieves 1p1h cross section with RPA included over the bare 1p1h cross section; these twodimensional tables of the ratio were supplied by the authors of Reference Nieves et al. (2011) and are used to apply the RPA correction in NEUT. The Nieves RPA calculation uses the local Fermi gas nuclear model, and NEUT only has a global Fermi gas model implemented for 1p1h interactions, but the authors of the RPA calculation have noted Gran et al. (2013) that the same ratio can be applied, with reasonable precision, to a global Fermi gas. Because of the model dependence, the same ratios cannot be applied to modify the 1p1h interactions calculated with a SF model, and no RPA calculation performed in the context of SF nuclear model is available. Two different RPA calculations are available from the same authors, termed relativistic and nonrelativistic, which affect the quenching of the RPA at high ( GeV). The ratio of nonrelativistic to relativistic RPA is shown in Figure 4. Both RPA models are investigated in this analysis as there is no guidance on which model is more physical. The ‘stray’ points in Figures 3 and 4 are artifacts from the authors of the RPA model, who provided the data used to produce these figures. The cause of these artifacts is unknown, but as these points lie outside the kinematically allowed region of (, ) space, they do not affect the RPA implementation in NEUT as no event outside this region can be generated.
With these different ingredients, three distinct candidate CCQE models are available in NEUT, which are all considered in this work:

RFG+relativistic RPA+2p2h

RFG+nonrelativistic RPA+2p2h

SF+2p2h.
The default values for all variable model parameters are listed in Table 2 and Table 1 for both RFG+RPA+2p2h models and SF+2p2h, respectively.
It should be noted that there are deficiencies for both models as currently implemented in NEUT. The RFG+RPA+2p2h model is very like the full Nieves model as both the RPA and 2p2h calculations are taken from it. However, the Nieves model consistently uses a local Fermi gas, whereas NEUT uses a global Fermi gas model for the 1p1h calculation. Currently there is no ability to vary the value of used in the Nieves model prediction as implemented in NEUT, making the fits slightly inconsistent in this regard^{2}^{2}2The value of the axial mass used for the 2p2h contribution to the cross section was fixed to GeV.. The SF+2p2h model has no RPA correction applied, which is physically inconsistent as the 2p2h enhancement is used (both corrections are due to complications in heavy nuclear targets). As previously noted, no RPA calculation appropriate for a SF model is currently available, so this inconsistency is unavoidable. The nuclear models used for the 1p1h calculation (SF) and the 2p2h calculation (LFG) are also inconsistent, and it has been remarked that the short range correlations included the SF nuclear model may be the same as some contributions to the Nieves 2p2h interaction model, so some contributions may be included twice.
Additionally, the Effective Spectral Function (ESF) Bodek et al. (2014, 2015) has been implemented in NEUT as described in Reference Wilkinson (2015), and is included for comparison with the other nominal models in Section V. The ESF enforces agreement with the longitudinal response function extracted from electron scattering data by modifying the initial state nucleon momentum distribution (using a simple parametrization of the Benhar SF model), and should be used with the Transverse Enhancement Model (TEM), which parametrizes the observed discrepancy between the longitudinal and transverse response functions extracted from electron scattering data as an enhancement to the magnetic form factor Bodek et al. (2011). By construction, the ESF+TEM agrees with elastic electron scattering data, and is extended to neutrino scattering data by modifying the LlewellynSmith interaction formalism for nucleons bound in a nucleus described by the ESF (and with the modified magnetic form factor from the TEM). This model was implemented too late to be a candidate model for the T2K oscillation analysis, and is not considered further in the fitting work described in this paper.
Model parameter  NEUT default value 

1.01 GeV/  
Fermi momentum,  209 MeV/ 
Meanfield width  200 MeV/ (Benhar nominal Benhar and Fabrocini (2000)) 
Norm. of the  Benhar nominal Benhar and Fabrocini (2000) 
correlation term  (correlated tail 20% of total) 
2p2h normalization  100% Nieves model Nieves et al. (2011); Gran et al. (2013) 
Axial form factor  Dipole 
Vector form factors  BBBA05 Bradford et al. (2006) 
Model parameter  NEUT default value 

1.01 GeV/  
Fermi momentum,  217 MeV/ 
RPA  Nieves relativistic or 
nonrelativistic model Nieves et al. (2011)  
2p2h normalization  100% Nieves model Nieves et al. (2011); Gran et al. (2013) 
Axial form factor  Dipole 
Vector form factors  BBBA05 Bradford et al. (2006) 
We note that both of our candidate models are expected to break down at low momentum transfer because they do not include nuclear effects such as nuclear excitations and collective resonances. In other analyses which fit models to CCQE data, bins which are dominated by low momentum transfer events are excluded Juszczak et al. (2010). In this analysis we have not followed any such bin masking procedure. Arguably, to obtain a realistic value of the model parameters, one should only fit the model in its stated region of validity. However, the main focus of this analysis is to obtain central values and errors for the T2K oscillation analysis, where the cross section model is used for all regions of phasespace, so some pragmatism is required.
Iii NuWro as a validation tool for new interaction models
The NuWro Monte Carlo generator for neutrino interactions has been developed over the past 10 years at the University of Wrocław Golan et al. (2012). It was the first MC generator to have an implementation of the Benhar SF Benhar and Fabrocini (2000) and the Nieves 2p2h model included Nieves et al. (2011); Gran et al. (2013), and served as the benchmark for the NEUT development of both models. The implementation of the SF model in NuWro was based on the code written for Reference Ankowski and Sobczyk (2008) and subsequently optimized for NuWro. The Nieves model implementation in NuWro used a series of lookup tables for the 2p2h crosssection as a function of leptonic variables for various nuclear targets and neutrino species so is very similar as in NEUT, although it has since been improved to use a more general formalism which depends on a number of nuclear response functions which can be extracted from the Nieves code, and therefore reduces the number of lookup tables required. The same generic model Sobczyk (2012) was used to simulate the initial and final hadronic states in NuWro as was used in NEUT. For both the SF and Nieves 2p2h models, NuWro and NEUT are in good agreement, which provides a useful validation of the NEUT implementations of these models.
Iv External datasets
Four datasets are used in the CCQE fits presented in this work: the MiniBooNE neutrino AguilarArevalo et al. (2010) (2010) and antineutrino AguilarArevalo et al. (2013) (2013) results; and the MINERA neutrino Fiorentini et al. (2013) (2013) and antineutrino Fields et al. (2013) (2013) results. All experimental details and information about these results, which is reproduced here, are taken from the references cited above unless otherwise stated.
The singledifferential cross section results are given in terms of , the fourmomentum transfer calculated from lepton kinematics under the quasielastic hypothesis, which is calculated using the equations:
(1) 
(2) 
where is the muon energy; , and are the neutron, proton and muon masses, respectively; and where is the binding energy of carbon assumed in the analysis^{3}^{3}3Note that the binding energy is just the value assumed when calculating , so we must use the same value as the experiments when producing comparable distributions, but it need not be consistent with the binding energy used in our simulation.. For both MiniBooNE datasets and the MINERA neutrino dataset, MeV; for the MINERA antineutrino dataset, MeV.
In the MiniBooNE analysis, is calculated from the unfolded and distributions. The MINERA analysis unfolds the distribution calculated with the reconstructed and values. The errors on the distributions for both experiments include the uncertainties relating to the muon reconstruction, so should cover the difference in the method used to produce the cross section results. We note that the main results of our analysis use the MiniBooNE doubledifferential results only, so there is no possible tension from differences between the methods used to produce distributions.
iv.1 MiniBooNE neutrino
The MiniBooNE CCQE data has been released as a doubledifferential cross section as a function of , where is the kinetic energy of the outgoing muon and is the angle between the incoming neutrino and outgoing muon. Differential cross sections were also released as a function of or , but the doubledifferential result was preferred as it is contains the most information and has minimal model dependence. The MiniBooNE data release included central values for each bin and the diagonal elements of the shapeonly covariance matrix; correlations between bins were not released. Additionally, the overall normalization uncertainly was given as 10.7% for neutrino running.
The MiniBooNE CCQE cross sections are released as both CCQEcorrected, and CCQElike measurements. The CCQElike sample is obtained by selecting events in which a muon was detected with no pions, but no requirement was made on the proton. The CCQEcorrected measurement is produced by subtracting background events (where the primary interaction was not CCQE) based on the NUANCE Casper (2002) generator prediction. The dominant background is CC1, and a dedicated sample was used to tune the NUANCE prediction which was used in the background subtraction. It should be noted that the NUANCE CC1 simulation included less decay. The published signal purity for the neutrino dataset is 77%.
CCQElike results are less model dependent than CCQEcorrected results (as they do not rely on the experiment’s own MC correction strategy), but make the analysis dependent on the simulation of the background in the MiniBooNE detector, which cannot be tuned to the MiniBooNE data in the same way MiniBooNE’s background model could be. CCQEcorrected results are used in this analysis. A downside of using the CCQEcorrected data is the explicit subtraction of less decay events in the MiniBooNE analysis, which forms part of the Nieves multinucleon–neutrino prediction which we treat as signal in the analysis. Unfortunately, there is no obvious way to account for this effect, so we ignore it for the analysis presented. We note that Nieves et al. also used the CCQEcorrected dataset to compare to their full models Nieves et al. (2012, 2013).
iv.2 MiniBooNE antineutrino
The MiniBooNE antineutrino data has been released in the same format as the neutrino mode data. Again, the doubledifferential CCQEcorrected results are used. The overall normalization uncertainty was given as 13.0% for antineutrino running. This is likely to be strongly correlated with the normalization uncertainty for the neutrino mode data, as the uncertainly comes mostly from the flux normalization uncertainty. However, as this information was not released, no correlation is assumed in this analysis.
The correction strategy for the antineutrino dataset is more complicated than for the neutrino mode sample because of the relatively high contamination in the beam, which is the largest background in the antineutrino CCQE sample (MiniBooNE is an unmagnetized detector). There is also a large CC1 background, the analogue of the CC1 contamination in the neutrino dataset. Two properties are used to measure the background AguilarArevalo et al. (2011): 8% of induced CC interactions produce no decay electron due to muonnucleus capture; and the induced CC1 events can be identified independently of induced CC1 as most mesons are absorbed. Unfortunately, this property makes CC1 a bigger background to the CCQE analysis in antineutrino mode, and means that there is no sample with which to directly tune the CC1 production from the NUANCE resonance model, so the neutrino mode CC1 has to be used (as was done for the neutrino mode sample). Other backgrounds are subtracted using the NUANCE interaction model after some tuning and corrections. As a result of the two large backgrounds in the antineutrino sample, the purity of the CCQElike sample is 61%, making the correction larger than for the neutrino mode sample.
iv.3 MinerA
The CCQE datasets from MINERA are released as CCQEcorrected singledifferential fluxaveraged cross section as a function of , where the flux has been averaged over the region GeV. There is an additional requirement that GeV, with as defined in Equation 1. Covariance matrices and central values have been released to perform fits to both shapeonly and absolutely normalized neutrino and antineutrino datasets. In this work, the absolutely normalized distributions have been used in the fit.
The correction strategy for the MINERA data is to fit the relative normalizations of simulated background distributions to the data in terms of the recoil energy, energy deposited outside a vertex region (the recoil region), and then subtract the predicted background from the CCQElike sample. The published purity for the neutrino dataset ranges from 65% at low to 40% at high (with an overall purity of 49%). The purity for the antineutrino dataset is given as 77%. The purity is lower for the neutrino analysis because events with a proton from the initial interaction are more complicated to reconstruct than those with a neutron^{4}^{4}4The antineutrino analysis has an additional cut requiring no additional (other than the muon) tracks from the vertex, and allows only one isolated energy shower, whereas the neutrino mode analysis allows two Fields et al. (2013); Fiorentini et al. (2013)..
In the MINERA CCQE analyses, the efficiency for selecting events with is very low because the MINOS near detector, downstream of MINERA, is used to tag muons. This introduces a small model dependence on the results because an RFG model was used to correct for the unsampled region of phasespace. The MINERA collaboration subsequently released a distribution where the cross section is measured for CCQE events with . As this dataset is less modeldependent, it has been used in the fits, and will be consistently used in this analysis. MINERA also made crosscorrelations between the neutrino and antineutrino datasets available in a data release after the publication of their CCQE papers. The correlation matrices released include both shape and normalization errors, but it is possible to extract shapeonly correlation matrices using the method given in Reference Katori (2008). The full matrix including both shape and normalization errors included is shown in Figure 5.
V Monte Carlo prediction
For each of the four experimental results included in the fit, one million CCQE and 2p2h events were generated with NEUT for each model using the default parameters given in Tables 2 and 1 and the published flux for each dataset. The flux averaged cross section predictions were produced using the following method:

For each event apply experimentspecific cuts and, if the event passes, calculate the relevant reconstructed quantity and fill the 1D or 2D event rate histogram.

Calculate the event rate by integrating the MC event rate histogram (flux cross section).

Integrate the published flux histogram to get the average flux.

Scale the filled histogram by the event rate divided by the average flux to get the flux averaged cross section per target nucleon.

Divide the content of each bin by the bin width.
The default predictions for a variety of models available in NEUT, as well as the data, are shown in Figures 6, 7 and 8 for the MINERA, MiniBooNE singledifferential and MiniBooNE doubledifferential samples, respectively. The Nieves 2p2h contribution is also shown on these plots for reference.
To produce a meaningful nominal for the MiniBooNE datasets, it is necessary to fit the MiniBooNE normalization parameters. The single and doubledifferential plots shown in Fig. 7 and 8 are scaled according to the MiniBooNE normalization parameter at the best fit point. The best fit values of the pull parameters and are given in Table 3. Additionally, the nominal predictions for the MiniBooNE doubledifferential datasets, without the scaling factor applied, are shown in Figure 9, which are easier to interpret by eye.
Fit type  

Neutrino 1D  RFG  0.7320.007  — 
SF+2p2h  0.7410.007  —  
RPA+2p2h  0.7600.007  —  
ESF+TEM  0.8040.008  —  
Antineutrino 1D  RFG  —  0.8050.011 
SF+2p2h  —  0.8260.011  
RPA+2p2h  —  0.7740.010  
ESF+TEM  —  0.8030.011  
Neutrino 2D  RFG  0.7250.011  — 
SF+2p2h  0.7560.011  —  
RPA+2p2h  0.7600.011  —  
ESF+TEM  0.8270.012  —  
Antineutrino 2D  RFG  —  0.8080.015 
SF+2p2h  —  0.8380.015  
RPA+2p2h  —  0.8020.015  
ESF+TEM  —  0.8330.015 
Note that the doubledifferential cross section plots shown in Figures 9 have been rebinned. In the distributions released by MiniBooNE, and used in the fits, there are 20 bins uniformly distributed between 1 and 1. For ease of presentation, these have been rebinned and the results are shown in eight slices of varying sizes, where merged bins have been averaged and their errors combined in quadrature.
Vi Fit procedure
All minimizations are performed using the MIGRAD algorithm of the MINUIT package James and Roos (1975), using the statistic:
(3) 
where are the model parameters varied in the fit, and are the diagonals of the MiniBooNE shapeonly covariance matrices for the neutrino and antineutrino results, is the crosscovariance matrix provided by MINERA, and and are the normalization parameters for MiniBooNE neutrino and antineutrino datasets, with published normalization uncertainties of (10.7%) and (13.0%)^{5}^{5}5Note that the MINERA normalization uncertainty is included in the covariance matrix, so also contributes a penalty term to the fit..
Fits to individual datasets only include the relevant terms from the definition in Equation 3, and fits to single MINERA datasets neglect crosscorrelations (the summation is only over the relevant eight bins).
vi.1 Parameter GoodnessofFit (PGoF) test
Standard goodnessoffit tests, such as the Pearson test used as an example here, test the agreement between prediction and data; however, some issues can arise with their use in global fits, as discussed in Reference Maltoni and Schwetz (2003). The basic problem is that much of the data will have limited power to constrain any one parameter, but agree well with the prediction regardless of the parameter values. These data will add little to the , but contribute another degree of freedom. Thus the found may be deceptively good despite not agreeing well with those parts of the dataset that actually have power to constrain key parameters. It is also possible that a dataset with a large number of datapoints (such as MiniBooNE) which does agree well with a model may hide disagreements with other datasets included in a global fit for which fewer datapoints are available (such as MINERA); again, the key problem is a dilution of the .
This problem is worsened in the case of datasets for which correlations between datapoints have not been included, where the DOF can be much less than 1, such as is the case for MiniBooNE. Looking at the Pearson test statistic is not very illuminating when fitting to both MiniBooNE and MINERA datasets.
The PGoF is a more rigorous test proposed in Reference Maltoni and Schwetz (2003) for fitting to global datasets and has been used extensively in sterile neutrino literature Conrad et al. (2013); Kopp et al. (2013), where there are often contradictory results coming from different experiments, and the fitters are fitting to many different experiments which are sensitive to different parameters. It is also referred to as the Likelihood Ratio Test in both statistics and other HEP literature. The PGoF test statistic is given by
(4) 
where are the parameters floated in the fits, is the number of datasets, is the minimum in a fit to all subsets of the data, and is the minimum obtained in a fit to the th subset of the data. The PGoF test statistic forms a distribution with the number of degrees of freedom
where and are the number of degrees of freedom for each fit.
The aim of the PGoF is to test the compatibility of the different datasets in the framework of the model. Put simply, it tests whether the best fit parameter values to subsets of the data pulls the fit parameters far from the best fit values found when fitting to all of the data. If different subsets favor very different values, then those subsets are not compatible in the framework of the model (though individually each may be able to find parameter combinations which produce a good fit).
A further advantage of the PGoF test in the situation where some of the data lacks correlations is that the number of degrees of freedom come from the number parameters varied in the fits, not from the number of bins that dataset contributes, which mitigates against the issue.
The PGoF test still assumes that the datasets follow a distribution, but allows for a lower effective number of degrees of freedom. This assumption is not quantitatively correct due to the aforementioned lack of correlations in the MiniBooNE CCQE data. The values returned should be taken with the caution that they highlight tensions between datasets, but are not to be interpreted in the same manner as they would if all correlations were reported.
Vii Fit results
vii.1 Fake data studies
The fitter was validated in two ways. Firstly, Asimov fake datasets Cowan et al. (2011) were produced to estimate the size of the errors that would be produced from the fit and used as a sanity check of the real fit results. The Asimov tests also provide a very basic test of the fitting framework developed for this analysis. Secondly, pull studies were performed to check that the definition given in Equation 3 is an unbiased estimator of the parameter central values and errors. For all parameters, the biases were less than 10% across the entire parameter range allowed in the fit, so we conclude that the fitter behaves well.
vii.2 Combined fit
The results for the combined fits to all four datasets are given for both relativistic and nonrelativistic RFG+RPA+2p2h models and the SF+2p2h model in Table 4. The best fit distributions are compared with data for MINERA in Figure 10, and for MiniBooNE in Figure 11. Relativistic RPA is used in the figures, as this was the best fit of the two RPA models available. In the legends of these figures, each line is given two values, the contribution from that dataset to the in the combined fit, and the total in parentheses. Note that in Figure 10, the contributions from MINERA are calculated for the individual datasets, which necessarily ignores crosscorrelations and makes these numbers slightly misleading. Explicitly, due to crosscorrelations, so the values shown in the figure should be treated with caution.
Fit type  /  (GeV/)  2p2h norm. (%)  (MeV/)  

RFG+rel.RPA+2p2h  97.8/228  1.150.03  2712  2235  0.790.03  0.780.03 
RFG+nonrel.RPA+2p2h  117.9/228  1.070.03  3412  2255  0.800.04  0.750.03 
SF+2p2h  97.5/228  1.330.02  0 (at limit)  2344  0.810.02  0.860.02 
It is clear from Figures 10 and 11 that MiniBooNE is not completely dominating the fits, as might be expected given the large number of bins in each of the MiniBooNE datasets. Indeed, these fits exploit the fact that, without correlations, . It is also clear that neither model fits all of the datasets perfectly at the best fit point, which is not reflected by the reduced values of 97.5/228 and 97.8/228 for the SF+2p2h and RFG+rel.RPA+2p2h models, respectively. As the MiniBooNE public data release lacks bintobin correlations, the contributions are not as large as would be expected for the number of bins contributed. This may explain why so many theoretical models are able to find good agreement with the MiniBooNE CCQE data.
In all fits, it was observed that the MiniBooNE normalization values tended to be suppressed for both neutrino and antineutrino datasets indicating that the MC underestimated the published data by 20–30% (10–20%) for the RFG+RPA+2p2h (SF+2p2h) models. It is not possible to accurately determine the favored MINERA normalization as the normalization uncertainty is included in the published covariance matrix, but the output distributions show that the MC normalization is approximately equal to the data normalization.
Because of the large pulls on the MiniBooNE normalization parameters, shapeonly fits were also performed (see Reference Wilkinson (2015) for further details). It was found that the best fit parameters were not significantly changed, indicating that there is not a significant bias to the other parameters caused by the large pulls on the MiniBooNE normalization parameters. A recent reanalysis of the MINERA flux Golan (2015) results in an increase to the normalization of previous MINERA cross section results including the CCQE samples used in this analysis. Although these updated datasets are not included in this work, we note that as the results were found to be largely unchanged in a shapeonly fit, the main results will not be significantly affected. Additionally, results from fits to individual datasets, and to various combinations of datasets can be found in Reference Wilkinson (2015).
vii.3 PGoF results
Using the PGoF test defined in Section VI.1, it is possible to test the compatibility between different subsets of the data. Tables 6, 7, and 8 show a breakdown of the four datasets used in the the combined fits for each initial CCQE model assumption. The Standard Goodness of Fit (SGoF) for each row is determined using Pearson’s test, where is found by minimizing the function given in Equation 3, including only the terms for the relevant datasets. The PGoF test is found by subtracting for each of the constituent datasets from the minimum of the combined dataset. The formulae for calculating the PGoF test statistic are given explicitly in Table 5. The for each dataset is again determined by minimizing the function given in Equation 3 with only the relevant terms included.
All  

MINERA  
MiniBooNE  
MA vs MB  
vs 
In each fit, the , 2p2h normalization, , and any MiniBooNE normalization terms are allowed to float.
/  SGoF (%)  /  PGoF (%)  

All  117.9/228  100.00  25.3/6  0.03 
MINERA  30.3/13  0.42  0.4/3  93.09 
MiniBooNE  65.7/212  100.00  3.4/3  33.09 
69.1/142  100.00  12.7/3  0.53  
46.1/83  99.97  10.4/3  1.55  
MA vs MB  117.9/228  100.00  21.9/3  0.01 
vs  117.9/228  100.00  2.6/3  45.12 
One subtlety must be kept in mind when analysing the results in Tables 6, 7, and 8: the PGoF test is only appropriate for statistically independent datasets. This makes the interpretation difficult for MINERA, where crosscorrelations are provided and used in the fits. Whenever a subset of data includes both MINERA and datasets, the fits include crosscorrelations, but if only one dataset is included, they are not. This means that two of the rows in each table give slightly unreliable results: “MINERA”, and “ vs ”. In each case, the function for the combined dataset includes crosscorrelations, and the functions for the subdivided dataset does not. The issue is most obvious in Table 7, where the “ vs ” row gives a negative PGoF . These values are still useful as a comparison between models and to give a rough idea of compatibility between datasets, but the exact values must be treated with caution.
/  SGoF (%)  /  PGoF (%)  

All  97.8/228  100.00  17.9/6  0.66 
MINERA  23.4/13  3.74  1.0/3  79.03 
MiniBooNE  58.6/212  100.00  2.0/3  57.69 
62.6/142  100.00  16.1/3  0.11  
38.5/83  100.00  6.1/3  10.75  
MA vs MB  97.8/228  100.00  15.9/3  0.12 
vs  97.8/228  100.00  3.3/3  100.00 
/  SGoF (%)  /  PGoF (%)  

All  97.5/228  100.00  41.1/6  0.00 
MINERA  12.6/13  47.75  1.0/3  79.49 
MiniBooNE  50.2/212  100.00  6.5/3  8.92 
54.8/142  100.00  25.1/3  0.00  
34.1/83  100.00  8.5/3  3.61  
MA vs MB  97.5/228  100.00  34.6/3  0.00 
vs  97.5/228  100.00  8.5/3  3.59 
The PGoF test highlights the incompatibility of the various datasets within the framework of the SF+2p2h and both RFG+RPA+2p2h models, despite the apparent goodness of fit when only considering . The level of agreement given in the final column of Tables 6, 7, and 8 should be interpreted as the level of agreement between the datasets included in that row. For example, it is clear that for all models considered the agreement found between the MINERA and MiniBooNE datasets (which include both neutrino and antineutrino samples) have the lowest level of agreement as shown by the “MA vs MB” row. In contrast, the level of agreement between the neutrino and antineutrino datasets (which include the MINERA and MiniBooNE samples) show relatively good agreement, indicating that fitting to the neutrino and antineutrino datasets separately produces similar best fit parameter values.
It is clear from Table 8 that the SF+2p2h model does not fit the various datasets well, the poor PGoF statistics indicate that the datasets favor very different parameter values when fit separately. This is particularly true for any fits involving the MiniBooNE neutrino dataset, though there is no a priori reason to exclude this dataset and improve the fit results. The PGoF tests for RFG+RPA+2p2h using both relativistic and nonrelativistic RPA, shown in Tables 6 and 7, show much better compatibility between experiments than SF+2p2h. There is still a considerable amount of tension, which is largely due to differences between MINERA and MiniBooNE. Because of the relatively poor consistency between datasets for the SF+2p2h model compared with RFG+rel.RPA+2p2h, the latter model is a better choice as the default model for T2K oscillation analyses.
vii.4 Rescaling parameter errors
Assuming Gaussian statistics, 1 errors on a single fit parameter are defined by the parameter value for which Olive et al. (2014). MINUIT uses this assumption when calculating the errors at the minimum, which were included with the best fit values for the combined fit in Table 4. However, as well as motivating the use of the PGoF test, the lack of bin correlations from MiniBooNE also means that Gaussian statistics no longer work as expected when estimating parameter errors.
There is a large body of literature looking at how this problem affects fits to parton density distributions, where global fits include a large number of datasets, many of which did not provide bin correlations Pumplin et al. (2001); Stump et al. (2001); Collins and Pumplin (2001). A summary of the work of one PDF fitting group is given in Reference Pumplin et al. (2001) and was used as a guide here. Their solution for producing reasonable parameter error estimates is to inflate the value of the used to define the 1 parameter errors, although no generic solution is offered for defining that value. In the case of the PDF fits in Reference Pumplin et al. (2001), the used was very large, 100, although it should be kept in mind that many more datasets are used in that fit than in the current work.
The PGoF gives a value for the incompatibility between the datasets: how much the increases between the best fit points of each experiment and the best fit point for the combined dataset. The PGoF value can therefore be used as a measure of how much the errors have to be inflated to cover the difference between the best fit parameter values from the combined fit and the best fit values of individual datasets, this is shown explicitly in Equation 5, where the value used to define the 1 error is given by , and the rescaling parameter is given by .
(5) 
Note that this PGoF rescaling procedure does not modify the correlations between parameters, it simply rescales the error on each parameter.
Fit type  /  (GeV/)  2p2h (%)  (MeV/) 

Unscaled  97.8/228  1.150.03  2712  2235 
PGoF scaling  1.150.06  2727  22311 
There is some ambiguity over which PGoF statistic to use, the ‘All’ or ‘MA vs MB’ row of Table 7, with values of 17.9/6 and 15.9/3 respectively. The more conservative value is from the ‘MA vs MB’ (because the greatest differences are between experiments, not between neutrino and antineutrino running), so this is used^{6}^{6}6It should also be noted that this rescaling procedure more than covers the difference between neutrino and antineutrino datasets.. To be explicit, we multiply the parameter errors from MINUIT by based on this statistic, as shown in Table 9. It can be seen from Table 9 that the 2p2h normalization is strongly suppressed and, even with the rescaled error, is nearly 3 away from the Nieves nominal model prediction. It is also clear that although the axial mass value preferred in the fit is not as strongly inflated as in fits to MiniBooNE data alone Juszczak et al. (2010); AguilarArevalo et al. (2010, 2013), it is still significantly higher than the value of GeV found by fitting to light target data and pion electroproduction data Bradford et al. (2006), and the inflated 1 error does not cover this difference.
Viii Discussion of the fit results
The results from the fits presented in Section VII.2 show that none of the models which are currently available in NEUT describe all of the CCQE data adequately. In particular, there is a significant difference between MiniBooNE and MINERA data which forces the model parameters to compromise between the two, as well as a large change in the normalization for the MiniBooNE datasets. Although the value obtained from the fit to the RFG+rel.RPA+2p2h model is lower than that obtained from past fits of the RFG model to MiniBooNE data alone (see Reference Abe et al. (2015c) as an example), it is still inconsistent with that obtained in global fits to light target bubble chamber data or high energy heavy target data Bradford et al. (2006). Additionally, the data requires a large supression of the nominal 2p2h model. The SF+2p2h model, in which the 2p2h component is completely suppressed, requires an inflated value to fit the data. This is unsurprising as the SF model alone does not significantly change the total cross section. Including an RPA calculation appropriate for the SF model is likely to reduce the tension with the 2p2h model, and is likely to change this conclusion significantly, this work will be revisited when such a calculation is available. Both fits also initially imply that there may need to be additional interactions used that may mimic CCQE interactions or change the shape of the distributions through additional, but currently unmodelled, effects in the nucleus.
There seems to be a shape problem with the 2p2h model, which leads to the suppression of the 2p2h cross section when 2p2h normalization is allowed to vary in the fit. Recall that at the best fit point 2p2h is suppressed to 27% of the Nieves nominal value as shown in Table 9. This suppression is driven by MINERA, which would completely suppress the 2p2h component of the model if MiniBooNE were not included in the fit Wilkinson (2015). The SF+2p2h model shows the same disagreement, and 2p2h is completely removed at the best fit point. The inability to change the shape of the 2p2h prediction is clearly a deficiency in the current NEUT implementation of the Nieves 2p2h model, and much better agreement might be found if more parameters could be varied in the model. A further significant issue is that the CCQEcorrected cross section results from both MINERA and MiniBooNE have part of the 2p2h signal region removed as a background (less decay). It is not clear how this issue should be treated in the fits, and has simply been ignored here, as indeed has been done by the 2p2h model builders Nieves et al. (2011); Martini et al. (2009). Future cross section measurements should be encouraged to focus on exclusive final states (CC) rather than initial state processes (CCQE), which will avoid such an issue in the future. As previously remarked, the RFG+RPA+2p2h model implemented in NEUT is not equivalent to the full Nieves model because the 1p1h component in NEUT uses a global, rather than local, Fermi gas nuclear model. Using a global Fermi gas model results in greater interaction strength in the low region than with a local Fermi gas, which is also where the 2p2h model contributes most interaction strength. It is possible that the 2p2h shape issue is due to a conflict with the 1p1h model, and that a more consistent LFG+RPA+2p2h model would resolve this issue.
Although both the RFG+rel.RPA+2p2h model and SF+2p2h model give reasonable agreement with data at the best fit point, it is difficult to trust standard goodness of fit tests as the lack of MiniBooNE correlations means that Gaussian statistics no longer work correctly. An alternative measure of the goodness of fit, the PGoF, was used to try to improve the situation. Although the PGoF procedure still assumes Gaussian statistics, it highlights disagreements within the combined dataset by dividing the dataset into subsamples. These disagreements are completely hidden by the standard goodness of fit tests because the MiniBooNE contribution is so low relative to the number of degrees of freedom it contributes to the fit. The PGoF showed that there was considerably better agreement between the best fit parameter values obtained in fits to subsamples of the data for the RFG+rel.RPA+2p2h fit, which gives some confidence to the fit result. For the SF+2p2h model, the fits to subsamples of the data pulled to drastically different parameter values at the best fit points, which is highly undesirable behaviour if the fit results are to be used as prior uncertainties in oscillation analyses, and indicates that the model is a bad fit to the global dataset. But the SF+2p2h model can fit individual datasets well (as is clear in Table 8), so should not be discounted completely.
The lack of reported MiniBooNE correlations and nonGaussian behaviour of the test statistic also means that standard parameter error estimation does not work, and returns smaller parameter errors than are reasonable given the level of disagreement between the datasets used in the fit. An unrealistically tight constraint on cross section parameters would lead to biases in the near detector fit for T2K. To circumvent this problem a PGoF error inflation procedure was defined to ensure that the 1 parameter errors cover the disagreement between the MINERA and MiniBooNE datasets. This is a conservative approach, but as no model seems able to describe all of the available data, such an approach was necessary. Such ad hoc procedures are necessary when incomplete information is available from some of the datasets included in the fit. The lack of information about bintobin correlations for the MiniBooNE datasets significantly complicates this analysis, and may significantly change the results. We note that in the literature, many statements about how well various models agree with the MiniBooNE datasets are made which assume Gaussian statistics. It is clear that an appreciation of this issue is important for future model comparisons, and that the availability of complete information for new cross section results will be critical for building consistent CCQE models.
For T2K, the results of this fit are part of a larger set of cross section model systematic uncertainties recommended by the NIWG that can be used as prior inputs for the oscillation analyses and various cross section analyses. In this case, the model used for these analyses is the RFG+rel.RPA+2p2h model, since the SF+2p2h model is disfavored in the fits and relativistic RPA is preferred over nonrelativistic RPA. The best fit parameters and uncertainties of the model are given in the second row of Table 9, and are correlated according to the matrix shown in Figure 12.
Ix Summary
In this paper, we have shown how T2K’s NIWG uses previously published CCQE datasets from the MiniBooNE and MINERA experiments to test CCQE+2p2h models in the NEUT neutrino interaction generator. For each model, the parameters that describe the data are fit, with both the SGoF and PGoF used to select the model that best describes the data. In this case, the RFG+rel.RPA+2p2h model is considered the best candidate, with GeV/, the normalization on the 2p2h model %, and MeV/. Tensions between the two experiments require an error scaling procedure outlined by the PGoF test, with the final result providing prior inputs into various future T2K analyses. This is the first time a comprehensive analysis has been performed with these models using a neutrino interaction generator and published, and the first time that such models will be used in an oscillation analysis with full detector simulations Abe et al. (2015d).
Moving away from the RFG model for CCQE interactions is an ambitious step for a neutrino experiment as it is a departure from the standard which has been used for decades Gallagher et al. (2011). The new models on the market are not perfect, and their implementation into NEUT and other neutrino interaction generators will always have technical foibles. However, further theoretical development of these models requires the engagement of the experimental community, and so using them in our simulations is essential to move the field onwards. It is also clear that the current approach of inflating is inadequate, and something better must be done in order to make precision measurements of neutrino oscillation parameters.
The fitting framework developed by the NIWG for this analysis is extensible and the general method for producing cross section errors developed in this work will be used with new CCQE models and datasets in future, and with new cross section channels entirely, to continue to contrain systematic errors for T2K oscillation and cross section analyses. The results from the CCQE fits presented here will also help inform the future model development required to fit the data. It is clear that alternative 2p2h models and fundamental parameters in the 2p2h model should be investigated, to see whether the disagreement with the 2p2h shape is telling us something meaningful about the Nieves model. It is also probable that the current RPA model is too inflexible, and this is partially responsible for the disagreement between MiniBooNE and MINERA data. Both of these problems may relate to the fact that for several years, the only data available for theorists to use for building models against was from the MiniBooNE neutrino dataset, which is difficult to use due to the lack of correlations and the explicit subtraction of less decay from the CCQEcorrected result. Converging on a new CCQE model which adequately describes all current and future data is likely to require several iterations between experimentalists and model builders. Confronting all the available models with a variety of data, as has been done in this analysis, and including these models in full Monte Carlo simulations, as will be done in T2K with the output from this analysis, is an important step in this cycle from the experimental side.
Acknowledgements.
The authors would like to thank the members of the T2K Collaboration and the authors of the NuWro generator for their help and support in this analysis. We thank the MINERA and MiniBooNE collaborations for assistance in understanding their results and studies beyond the published work. We are grateful to J. Nieves and his collaborators for providing the code required to implement their 2p2h and RPA models in NEUT. We are grateful to T. Katori and P. Litchfield for helpful comments on the manuscript. We thank the JPARC staff for superb accelerator performance and the CERN NA61/SHINE collaboration for providing valuable particle production data. We acknowledge the support of MEXT, Japan; NSERC, NRC and CFI, Canada; National Science Centre (NCN), Poland; MINECO and ERDF funds, Spain; SNSF and SER, Switzerland; STFC, UK; and DOE, USA. We also thank CERN for the UA1/NOMAD magnet, DESY for the HERAB magnet mover system, NII for SINET4, the WestGrid and SciNet consortia in Compute Canada, GridPP, UK, and the Emerald High Performance Computing facility in the Centre for Innovation, UK.References
 Llewellyn Smith (1972) C. H. Llewellyn Smith, Phys. Rep. 3, 261 (1972).
 Moniz and Smith (1972) E. J. Moniz and R. A. Smith, Nucl. Phys. B101, 605 (1972).
 Gallagher et al. (2011) H. Gallagher, G. Garvey, and G. P. Zeller, Annu. Rev. Nucl. Part. Sci. 61, 355 (2011).
 Moniz et al. (1971) E. J. Moniz, I. Sick, R. R. Whitney, J. R. Ficenec, R. D. Kephart, and W. P. Trower, Phys. Rev. Lett. 26, 445 (1971).
 Budd et al. (2003) H. S. Budd, A. Bodek, and J. Arrington, in 2nd International Workshop on NeutrinoNucleus Interactions in the Few GeV Region (NuInt 02) Irvine, California, December 1215, 2002 (2003) arXiv:hepex/0308005 [hepex] .
 Bernard et al. (2002) V. Bernard, L. Elouadrhiri, and U.G. Meißner, Journal of Physics G: Nuclear and Particle Physics 28, R1 (2002).
 Kuzmin et al. (2006) K. S. Kuzmin, V. V. Lyubushkin, and V. A. Naumov, Acta Phys. Polon. B37, 2337 (2006), arXiv:hepph/0606184 [hepph] .
 Kuzmin et al. (2008) K. S. Kuzmin, V. V. Lyubushkin, and V. A. Naumov, Eur. Phys. J. C54, 517 (2008), arXiv:0712.4384 [hepph] .
 Lyubushkin et al. (2009) V. Lyubushkin et al., The European Physical Journal C 63, 355 (2009).
 AguilarArevalo et al. (2010) A. AguilarArevalo et al. (MiniBooNE Collaboration), Phys. Rev. D81, 092005 (2010), arXiv:1002.2680 [hepex] .
 AguilarArevalo et al. (2013) A. AguilarArevalo et al. (MiniBooNE Collaboration), Phys. Rev. D88, 032001 (2013), arXiv:1301.7067 [hepex] .
 Gran et al. (2006) R. Gran et al. (K2K), Phys. Rev. D74, 052002 (2006), arXiv:hepex/0603034 [hepex] .
 Abe et al. (2015a) K. Abe et al. (T2K), Phys. Rev. D92, 112003 (2015a), arXiv:1411.6264 [hepex] .
 Abe et al. (2015b) K. Abe et al. (T2K), Phys. Rev. D91, 112002 (2015b), arXiv:1503.07452 [hepex] .
 Dorman (2009) M. Dorman (MINOS Collaboration), AIP Conf. Proc. 1189, 133 (2009).
 Wascko (2008) M. O. Wascko (SciBooNE), Proceedings, 23rd International Conference on Neutrino Physics and Astrophysics (Neutrino 2008), J. Phys. Conf. Ser. 136, 042028 (2008).
 Fiorentini et al. (2013) G. Fiorentini et al. (MINERvA Collaboration), Phys. Rev. Lett. 111, 022502 (2013), arXiv:1305.2243 [hepex] .
 Fields et al. (2013) L. Fields et al. (MINERvA Collaboration), Phys. Rev. Lett. 111, 022501 (2013), arXiv:1305.2234 [hepex] .
 Bodek et al. (2011) A. Bodek, H. S. Budd, and E. Christy, Eur. Phys. J. C 71, 1726 (2011).
 Morfin et al. (2012) J. G. Morfin, J. Nieves, and J. T. Sobczyk, Adv. High Energy Phys. 2012, 934597 (2012), arXiv:1209.6586 [hepex] .
 Garvey et al. (2015) G. T. Garvey, D. A. Harris, H. A. Tanaka, R. Tayloe, and G. P. Zeller, Phys. Rept. 580, 1 (2015), arXiv:1412.4294 [hepex] .
 AlvarezRuso et al. (2014) L. AlvarezRuso, Y. Hayato, and J. Nieves, New J. Phys. 16, 075015 (2014), arXiv:1403.2673 [hepph] .
 Benhar and Fabrocini (2000) O. Benhar and A. Fabrocini, Phys. Rev. C62, 034304 (2000), arXiv:nuclth/9909014 [nuclth] .
 Ankowski and Sobczyk (2006) A. M. Ankowski and J. T. Sobczyk, Phys. Rev. C 74, 054316 (2006).
 Butkevich (2009) A. Butkevich, Phys. Rev. C80, 014610 (2009), arXiv:0904.1472 [nuclth] .
 Bodek et al. (2014) A. Bodek, M. Christy, and B. Coopersmith, Eur. Phys. J. C74, 3091 (2014), arXiv:1405.0583 [hepph] .
 Leitner et al. (2009) T. Leitner, O. Buss, L. AlvarezRuso, and U. Mosel, Phys. Rev. C 79, 034601 (2009).
 Maieron et al. (2003) C. Maieron, M. Martinez, J. Caballero, and J. Udias, Phys. Rev. C68, 048501 (2003), arXiv:nuclth/0303075 [nuclth] .
 Meucci et al. (2004) A. Meucci, C. Giusti, and F. D. Pacati, Nucl. Phys. A739, 277 (2004), arXiv:nuclth/0311081 [nuclth] .
 Lovato et al. (2013) A. Lovato, S. Gandolfi, R. Butler, J. Carlson, E. Lusk, S. C. Pieper, and R. Schiavilla, Phys. Rev. Lett. 111, 092501 (2013), arXiv:1305.6959 [nuclth] .
 Pandey et al. (2015) V. Pandey, N. Jachowicz, T. Van Cuyck, J. Ryckebusch, and M. Martini, Phys. Rev. C92, 024606 (2015), arXiv:1412.4624 [nuclth] .
 Buss et al. (2012) O. Buss, T. Gaitanos, K. Gallmeister, H. van Hees, M. Kaskulov, et al., Phys. Rept. 512, 1 (2012), arXiv:1106.1344 [hepph] .
 Hayato (2009) Y. Hayato, Acta Physica Polonica B 40, 2477 (2009).
 Nieves et al. (2011) J. Nieves, I. R. Simo, and M. J. V. Vacas, Phys. Rev. C 83, 045501 (2011).
 Gran et al. (2013) R. Gran, J. Nieves, F. Sanchez, and M. J. Vicente Vacas, Phys. Rev. D88, 113007 (2013), arXiv:1307.8105 [hepph] .
 Martini et al. (2009) M. Martini, M. Ericson, G. Chanfray, and J. Marteau, Phys. Rev. C80, 065501 (2009), arXiv:0910.2622 [nuclth] .
 Nieves et al. (2012) J. Nieves, I. Ruiz Simo, and M. Vicente Vacas, Phys. Lett. B707, 72 (2012), arXiv:1106.5374 [hepph] .
 Martini et al. (2011) M. Martini, M. Ericson, and G. Chanfray, Phys. Rev. C84, 055502 (2011), arXiv:1110.0221 [nuclth] .
 Nieves et al. (2013) J. Nieves, I. Ruiz Simo, and M. Vicente Vacas, Phys. Lett. B721, 90 (2013), arXiv:1302.0703 [hepph] .
 Martini and Ericson (2013) M. Martini and M. Ericson, Phys. Rev. C87, 065501 (2013), arXiv:1303.7199 [nuclth] .
 Coloma and Huber (2013) P. Coloma and P. Huber, Phys. Rev. Lett. 111, 221802 (2013), arXiv:1307.1243 [hepph] .
 Ankowski et al. (2015) A. M. Ankowski, O. Benhar, P. Coloma, P. Huber, C.M. Jen, C. Mariani, D. Meloni, and E. Vagnoni, Phys. Rev. D92, 073014 (2015), arXiv:1507.08560 [hepph] .
 Jachowicz et al. (1999) N. Jachowicz, S. Rombouts, K. Heyde, and J. Ryckebusch (ISOLDE), Phys. Rev. C59, 3246 (1999).
 Jachowicz et al. (2002) N. Jachowicz, K. Heyde, J. Ryckebusch, and S. Rombouts, Phys. Rev. C65, 025501 (2002).
 Abe et al. (2015c) K. Abe et al. (T2K Collaboration), Phys. Rev. D91, 072010 (2015c), arXiv:1502.01550 [hepex] .
 Ciofi degli Atti and Simula (1996) C. Ciofi degli Atti and S. Simula, Phys. Rev. C53, 1689 (1996), arXiv:nuclth/9507024 [nuclth] .
 Furmanski (2015) A. Furmanski, Chargedcurrent Quasielasticlike neutrino interactions at the T2K experiment, Ph.D. thesis, University of Warwick (2015).
 Sinclair (2014) P. Sinclair, Implementation of a multinucleon neutrino interaction simulation and comparison with T2K data, Ph.D. thesis, Imperial College London (2014).
 Sobczyk (2012) J. T. Sobczyk, Phys. Rev. C 86, 015504 (2012).
 Walton et al. (2015) T. Walton et al. (MINERvA), Phys. Rev. D91, 071301 (2015), arXiv:1409.4497 [hepex] .
 Oset and Salcedo (1987) E. Oset and L. Salcedo, Nuclear Physics A 468, 631 (1987).
 Singh et al. (1998) S. K. Singh, M. J. VicenteVacas, and E. Oset, Phys. Lett. B416, 23 (1998), [Erratum: Phys. Lett.B423,428(1998)].
 Bodek et al. (2015) A. Bodek, M. E. Christy, and B. Coopersmith, “Effective spectral function for quasielastic scattering on nuclei from H to Pb,” (2015), arXiv:1409.8545 [nuclth] .
 Wilkinson (2015) C. Wilkinson, Constraining neutrino interaction uncertainties for oscillation experiments, Ph.D. thesis, University of Sheffield (2015).
 Bradford et al. (2006) R. Bradford, A. Bodek, H. Budd, and J. Arrington, Nuclear Physics B  Proceedings Supplements 159, 127 (2006).
 Juszczak et al. (2010) C. Juszczak, J. T. Sobczyk, and J. Żmuda, Phys. Rev. C 82, 045502 (2010).
 Golan et al. (2012) T. Golan, C. Juszczak, and J. T. Sobczyk, Phys. Rev. C86, 015505 (2012), arXiv:1202.4197 [nuclth] .
 Ankowski and Sobczyk (2008) A. M. Ankowski and J. T. Sobczyk, Phys. Rev. C77, 044311 (2008), arXiv:0711.2031 [nuclth] .
 Casper (2002) D. Casper, Nucl. Phys. Proc. Suppl. 112, 161 (2002), arXiv:hepph/0208030 [hepph] .
 AguilarArevalo et al. (2011) A. AguilarArevalo et al. (MiniBooNE), Phys. Rev. D84, 072005 (2011), arXiv:1102.1964 [hepex] .
 Katori (2008) T. Katori, A Measurement of the muon neutrino charged current quasielastic interaction and a test of Lorentz violation with the MiniBooNE experiment, Ph.D. thesis, Indiana U. (2008).
 James and Roos (1975) F. James and M. Roos, Comput. Phys. Commun. 10, 343 (1975).
 Maltoni and Schwetz (2003) M. Maltoni and T. Schwetz, Phys. Rev. D68, 033020 (2003), arXiv:hepph/0304176 [hepph] .
 Conrad et al. (2013) J. Conrad, C. Ignarra, G. Karagiorgi, M. Shaevitz, and J. Spitz, Adv. High Energy Phys. 2013, 163897 (2013), arXiv:1207.4765 [hepex] .
 Kopp et al. (2013) J. Kopp, P. A. N. Machado, M. Maltoni, and T. Schwetz, JHEP 05, 050 (2013), arXiv:1303.3011 [hepph] .
 Cowan et al. (2011) G. Cowan, K. Cranmer, E. Gross, and O. Vitells, Eur. Phys. J. C71, 1554 (2011), arXiv:1007.1727 [physics.dataan] .
 Golan (2015) T. Golan (Presented at the 10th international workshop on neutrinonucleus interactions in the fewGeV region (NuInt2015), 2015).
 Olive et al. (2014) K. Olive et al. (Particle Data Group), Chin. Phys. C38, 090001 (2014).
 Pumplin et al. (2001) J. Pumplin, D. Stump, and W. Tung, Phys. Rev. D65, 014011 (2001), arXiv:hepph/0008191 [hepph] .
 Stump et al. (2001) D. Stump, J. Pumplin, R. Brock, D. Casey, J. Huston, et al., Phys. Rev. D65, 014012 (2001), arXiv:hepph/0101051 [hepph] .
 Collins and Pumplin (2001) J. C. Collins and J. Pumplin, (2001), arXiv:hepph/0105207 [hepph] .
 Abe et al. (2015d) K. Abe et al. (T2K), (2015d), arXiv:1512.02495 [hepex] .