---
abstract: 'The [MiniBooNE]{}large axial mass anomaly has prompted a great deal of theoretical work on sophisticated Charged Current Quasi-Elastic (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 [<span style="font-variant:small-caps;">NEUT</span>]{}, T2K’s primary neutrino interaction event generator. In this paper, we give an overview of the models implemented, and present fits to published $\nu_\mu$ and $\overline{\nu}_\mu$ CCQE cross section measurements from the [MiniBooNE]{}and [MINER$\nu$A]{}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.'
author:
- 'C. Wilkinson'
- 'R. Terri'
- 'C. Andreopoulos'
- 'A. Bercellie'
- 'C. Bronner'
- 'S. Cartwright'
- 'P. de Perio'
- 'J. Dobson'
- 'K. Duffy'
- 'A.P. Furmanski'
- 'L. Haegel'
- 'Y. Hayato'
- 'A. Kaboth'
- 'K. Mahn'
- 'K.S. McFarland'
- 'J. Nowak'
- 'A. Redij'
- 'P. Rodrigues'
- 'F.Sánchez'
- 'J.D. Schwehr'
- 'P. Sinclair'
- 'J.T. Sobczyk'
- 'P. Stamoulis'
- 'P. Stowell'
- 'R. Tacik'
- 'L. Thompson'
- 'S. Tobayama'
- 'M.O. Wascko'
- 'J. Żmuda'
bibliography:
- 'niwgpaperdraft.bib'
title: 'Testing CCQE and 2p2h models in the [<span style="font-variant:small-caps;">NEUT</span>]{}neutrino interaction generator with published datasets from the [MiniBooNE]{}and [MINER$\nu$A]{}experiments'
---

Introduction {#sec:intro}
============

Charged Current Quasi-Elastic (CCQE) scattering ($\nu_{\mu} + n \rightarrow p + \mu^{-}$) is the dominant neutrino interaction process for muon (anti)neutrinos impinging on a nuclear target at neutrino energies on the order of 1GeV. Because CCQE is a two-body 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 $\nu_\mu^{\bracketbar}$ disappearance or $\nu_e^{\bracketbar}$ 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 neutrino-nucleus 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. Neutrino-nucleon CCQE scattering for free nucleons is described by the Llewellyn-Smith formalism [@llewellyn-smith], which has been extended to cover neutrino-nucleus CCQE scattering in the Smith-Moniz RFG model [@smith-moniz], 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 [@review; @moniz:fg_escat] is the axial mass, [$M_{\textrm{A}}$]{}. Results from a global analysis of neutrino-deuterium scattering experiments and pion electroproduction data find ${\ensuremath{M_{\textrm{A}}}\xspace}= 1.00\pm0.02$GeV/$c^{2}$ [@bba03], which is consistent with other analyses [@axialStructure; @kuzmin_2006; @kuzmin_2008]. These results are also consistent with high energy neutrino beam experiments on heavy nuclear targets [@NOMAD_CCQE].\
Recent differential CCQE cross section results from the [MiniBooNE]{}collaboration [@mb-ccqe-2010; @mb-ccqe-antinu-2013] are significantly higher than expectation, which can only be accounted for in the framework of the Smith-Moniz 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 [@k2kCCQE], which reported a value of ${\ensuremath{M_{\textrm{A}}}\xspace}=1.20\pm0.12$GeV/$c^{2}$. Both experiments exhibited not only a larger-than-expected axial mass, but also a supression of low-[$Q^{2}$]{}events relative to the expection from the Smith-Moniz RFG model. Other experiments using heavy nuclear targets with beam energies in the few-GeV region have also measured cross sections which are consistent with an inflation of the axial mass [@t2kCCQE; @ingridCCQE; @minosCCQE; @Wascko:2008zz], although these results do not paint a coherent picture. More recently, the [MINER$\nu$A]{}experiment [@minerva-nu-ccqe; @minerva-antinu-ccqe], which is at a somewhat higher neutrino energy than [MiniBooNE]{}, has shown good agreement with the Smith-Moniz RFG model with [$M_{\textrm{A}}$]{}= 1.00GeV/$c^{2}$, but requires an enhancement to the transverse component of the cross section, an effect also seen in electron-nucleus scattering [@tem]. 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.00GeV/$c^{2}$, with a suppressed cross section at low-[$Q^{2}$]{}and larger cross section at higher-[$Q^{2}$]{}relative to a simple RFG model. Comprehensive reviews of available CCQE cross section models can be found in References [@review; @Morfin:2012kn; @garvey_review_2014; @hayato_review_2014].

More sophisticated descriptions of the initial state of the nucleus than the RFG model provides are available from a number of authors [@sf; @ankowski_SF; @butkevich_2009; @eff-sf]. 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 non-constant 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 non-interacting fermions are available [@leitner_2009; @madrid_2003; @meucci_2003; @lovato_2013; @pandey_2014]. However, with the exception of the GiBUU interaction model [@gibuu; @leitner_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 [<span style="font-variant:small-caps;">NEUT</span>]{}neutrino interaction generator [@hayato:neut].

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_review_2014]. 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; @nievesExtension] and Martini [*et al.*]{} [@martini] go beyond the impulse approximation and include diagrams where two nucleons are involved in the interaction. This adds significant strength to the CCQE-like cross section and explains the difference in normalization observed in the [MiniBooNE]{}data, which was previously modelled with a large axial mass [@nieves_MB_nu_2011; @martini_MB_nu_2011; @nieves_MB_anu_2013; @martini_MB_anu_2013]. Because these 2p2h models are not two-body processes, they are expected to lead to significant biases in the neutrino energy reconstruction from the outgoing lepton which assumes CCQE kinematics [@coloma_erec; @ankowski_erec]. Additionally, the Random Phase Approximation (RPA) is a nuclear screening effect that modifies the propagator for interactions in nuclear matter [@Jachowicz:1998fn; @Jachowicz:2002rr; @nieves; @martini] and has a significant effect on the differential cross section as a function of [$Q^{2}$]{}, suppressing the cross section in the low-[$Q^{2}$]{}region and enhancing the cross section for ${\ensuremath{Q^{2}}\xspace}\gtrsim 0.5$GeV$^{2}$. 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 neutrino-nucleus 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]{} [@mb-ccqe-2010; @mb-ccqe-antinu-2013] and [MINER$\nu$A]{} [@minerva-nu-ccqe; @minerva-antinu-ccqe] datasets to a variety of models implemented in [<span style="font-variant:small-caps;">NEUT</span>]{}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 [<span style="font-variant:small-caps;">NEUT</span>]{}default central value for the axial mass [$M_{\textrm{A}}$]{}= 1.21 GeV$/c^{2}$ based on the value found by the K2K experiment [@k2kCCQE], with an error large enough to cover fits to the [MiniBooNE]{}neutrino mode CCQE dataset [@mb-ccqe-2010], as is fully described in Reference [@Abe:2015awa]. 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 [<span style="font-variant:small-caps;">NEUT</span>]{}model.

The models which have been implemented in the [<span style="font-variant:small-caps;">NEUT</span>]{}generator are discussed in Section \[sec:neut\] and Section \[sec:nuwro\] discusses cross-generator validation. Section \[sec:extdata\] gives a brief overview of the [MiniBooNE]{}and [MINER$\nu$A]{}data used in the fit. The nominal [<span style="font-variant:small-caps;">NEUT</span>]{}predictions for these datasets are shown in Section \[sec:mcpred\] for a variety of models. Section \[sec:fitprod\] discusses the fit procedure. The results of fake data studies and the fit to external data are given in Section \[sec:results\]. In Section \[sec:discussion\] we interpret the results and discuss possible implications in cross section and neutrino oscillation analyses and Section \[sec:summary\] summarizes the results.

Interaction models in [<span style="font-variant:small-caps;">NEUT</span>]{} {#sec: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 [<span style="font-variant:small-caps;">NEUT</span>]{}, 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 probability distribution for initial state protons within an oxygen nucleus for Benhar’s SF model [@sf] as a function of the removal energy ($E_{\mathrm{R}}$) and the magnitude of the nucleon momentum ($|\mathbf{p}|$). The SF is normalized such that the integral of this distribution is 1.[]{data-label="fig:SF_2D_distribution"}](figures/SF_oxygen_2D_new.pdf){width="0.9\columnwidth"}

The [<span style="font-variant:small-caps;">NEUT</span>]{}implementation of the SF model from Omar Benhar and collaborators [@sf] is described fully in Reference [@andy_thesis]. 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 \[fig:SF\_2D\_distribution\]. Pauli blocking is implemented as a hard cut-off: final state nucleons with momenta less than the Fermi momentum ${\ensuremath{p_{\mathrm{F}}}\xspace}^{\mbox{\scriptsize{SF}}}$ 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 \[fig:SF\_dials\], where the two-dimensional 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 [<span style="font-variant:small-caps;">NEUT</span>]{}which modify the SF as illustrated in Figure \[fig:SF\_dials\]. The default values for these parameters are given in Table \[tab:SF\_params\]. The mean field width and normalization of the correlation term are well-constrained by electron–scattering data [@andy_thesis] 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 ${\ensuremath{p_{\mathrm{F}}}\xspace}^{\mbox{\scriptsize{RFG}}}$ affects a wide range of [$Q^{2}$]{}, whereas changing ${\ensuremath{p_{\mathrm{F}}}\xspace}^{\mbox{\scriptsize{SF}}}$ only affects very low [$Q^{2}$]{}events.

The [multinucleon–neutrino]{}(2p2h) model from Nieves [*et al.*]{} [@nieves; @nievesExtension] has been implemented in [<span style="font-variant:small-caps;">NEUT</span>]{}as described in Reference [@peterthesis]. 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 [@multinucleon] for simulating the initial and final hadronic states was used for generating [<span style="font-variant:small-caps;">NEUT</span>]{}events[^1]. The discrepancy between the leptonic and hadronic simulation makes the current [<span style="font-variant:small-caps;">NEUT</span>]{}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 [@minerva_ccqe_p_2014]). For this reason, only leptonic measurements are used in this analysis. As the Nieves model is very complex, the current [<span style="font-variant:small-caps;">NEUT</span>]{}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 $\pi$-less $\Delta$ decay contributions, where a nucleon excited into a $\Delta(1232)$ resonance decays without producing a pion [@oset-fsi; @singh_1998]. Contributions from $\pi$-less $\Delta$ decay were previously implemented in [<span style="font-variant:small-caps;">NEUT</span>]{}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.

[0.9]{} ![The ratio of the [CCQE]{}cross section including the non-relativistic RPA model to the [CCQE]{}cross section without RPA, shown for both muon neutrino and muon antineutrino interactions on carbon. An enhancement of the ratio can be seen at high [$Q^{2}$]{}, and a suppression can be seen at low [$Q^{2}$]{}(and close to the kinematic boundary for antineutrinos). These [$E_{\nu}$]{}and [$Q^{2}$]{}dependent tables are used in [<span style="font-variant:small-caps;">NEUT</span>]{}to apply the RPA model. For these plots, an axial mass value of [$M_{\textrm{A}}$]{}= 1.01 GeV$/c^{2}$ was used.[]{data-label="fig:RPA_ratio_NEUT"}](figures/NU_2D-RPA_ratio.pdf "fig:"){width="\textwidth"}

[0.9]{} ![The ratio of the [CCQE]{}cross section including the non-relativistic RPA model to the [CCQE]{}cross section without RPA, shown for both muon neutrino and muon antineutrino interactions on carbon. An enhancement of the ratio can be seen at high [$Q^{2}$]{}, and a suppression can be seen at low [$Q^{2}$]{}(and close to the kinematic boundary for antineutrinos). These [$E_{\nu}$]{}and [$Q^{2}$]{}dependent tables are used in [<span style="font-variant:small-caps;">NEUT</span>]{}to apply the RPA model. For these plots, an axial mass value of [$M_{\textrm{A}}$]{}= 1.01 GeV$/c^{2}$ was used.[]{data-label="fig:RPA_ratio_NEUT"}](figures/ANU_2D-RPA_ratio.pdf "fig:"){width="\textwidth"}

RPA [@nieves] is implemented into [<span style="font-variant:small-caps;">NEUT</span>]{}as a modification to the 1p1h cross section as a function of [$E_{\nu}$]{}and [$Q^{2}$]{}. Figure \[fig:RPA\_ratio\_NEUT\] shows the ratio of the Nieves 1p1h cross section with RPA included over the bare 1p1h cross section; these two-dimensional tables of the ratio were supplied by the authors of Reference [@nieves] and are used to apply the RPA correction in [<span style="font-variant:small-caps;">NEUT</span>]{}. The Nieves RPA calculation uses the [*local*]{} Fermi gas nuclear model, and [<span style="font-variant:small-caps;">NEUT</span>]{}only has a [*global*]{} Fermi gas model implemented for 1p1h interactions, but the authors of the RPA calculation have noted [@nievesExtension] 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 non-relativistic, which affect the quenching of the RPA at high [$Q^{2}$]{}($\gtrsim 0.5$ GeV$^{2}$). The ratio of non-relativistic to relativistic RPA is shown in Figure \[fig:RPA\_types\_ratio\_NEUT\]. Both RPA models are investigated in this analysis as there is no guidance on which model is more physical. The ‘stray’ points in Figures \[fig:RPA\_ratio\_NEUT\] and \[fig:RPA\_types\_ratio\_NEUT\] 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 ([$E_{\nu}$]{}, [$Q^{2}$]{}) space, they do not affect the RPA implementation in [<span style="font-variant:small-caps;">NEUT</span>]{}as no event outside this region can be generated.

[0.9]{} ![The ratio of the non-relativistic RPA correction to the relativistic RPA correction, shown for both muon neutrino and muon antineutrino interactions on carbon. These [$E_{\nu}$]{}and [$Q^{2}$]{}dependent tables are used to reweight [<span style="font-variant:small-caps;">NEUT</span>]{}events from one RPA model to the other. By default, [<span style="font-variant:small-caps;">NEUT</span>]{}events are generated with the non-relativistic RPA model.[]{data-label="fig:RPA_types_ratio_NEUT"}](figures/NU_2D-RPA_types_ratio.pdf "fig:"){width="\textwidth"}

[0.9]{} ![The ratio of the non-relativistic RPA correction to the relativistic RPA correction, shown for both muon neutrino and muon antineutrino interactions on carbon. These [$E_{\nu}$]{}and [$Q^{2}$]{}dependent tables are used to reweight [<span style="font-variant:small-caps;">NEUT</span>]{}events from one RPA model to the other. By default, [<span style="font-variant:small-caps;">NEUT</span>]{}events are generated with the non-relativistic RPA model.[]{data-label="fig:RPA_types_ratio_NEUT"}](figures/ANU_2D-RPA_types_ratio.pdf "fig:"){width="\textwidth"}

With these different ingredients, three distinct candidate CCQE models are available in [<span style="font-variant:small-caps;">NEUT</span>]{}, which are all considered in this work:

1.  RFG+relativistic RPA+2p2h

2.  RFG+non-relativistic RPA+2p2h

3.  SF+2p2h.

The default values for all variable model parameters are listed in Table \[tab:RFG\_params\] and Table \[tab:SF\_params\] 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 [<span style="font-variant:small-caps;">NEUT</span>]{}. 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 [<span style="font-variant:small-caps;">NEUT</span>]{}uses a [*global*]{} Fermi gas model for the 1p1h calculation. Currently there is no ability to vary the value of [$M_{\textrm{A}}$]{}used in the Nieves model prediction as implemented in [<span style="font-variant:small-caps;">NEUT</span>]{}, making the fits slightly inconsistent in this regard[^2]. 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) [@eff-sf; @eff-sf_proc] has been implemented in [<span style="font-variant:small-caps;">NEUT</span>]{}as described in Reference [@my_thesis], and is included for comparison with the other nominal models in Section \[sec:mcpred\]. 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 [@tem]. By construction, the ESF+TEM agrees with elastic electron scattering data, and is extended to neutrino scattering data by modifying the Llewellyn-Smith 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.

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 [@mb-ccqe-wroclaw]. 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 phase-space, so some pragmatism is required.

NuWro as a validation tool for new interaction models {#sec:nuwro}
=====================================================

The NuWro Monte Carlo generator for neutrino interactions has been developed over the past $\sim$10 years at the University of Wrocław [@Golan:2012wx]. It was the first MC generator to have an implementation of the Benhar SF [@sf] and the Nieves 2p2h model included [@nieves; @nievesExtension], and served as the benchmark for the [<span style="font-variant:small-caps;">NEUT</span>]{}development of both models. The implementation of the SF model in NuWro was based on the code written for Reference [@Ankowski:2007uy] and subsequently optimized for NuWro. The Nieves model implementation in NuWro used a series of lookup tables for the 2p2h cross-section as a function of leptonic variables for various nuclear targets and neutrino species so is very similar as in [<span style="font-variant:small-caps;">NEUT</span>]{}, 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 [@multinucleon] was used to simulate the initial and final hadronic states in NuWro as was used in [<span style="font-variant:small-caps;">NEUT</span>]{}. For both the SF and Nieves 2p2h models, NuWro and [<span style="font-variant:small-caps;">NEUT</span>]{}are in good agreement, which provides a useful validation of the [<span style="font-variant:small-caps;">NEUT</span>]{}implementations of these models.

External datasets {#sec:extdata}
=================

Four datasets are used in the CCQE fits presented in this work: the [MiniBooNE]{}neutrino [@mb-ccqe-2010] (2010) and antineutrino [@mb-ccqe-antinu-2013] (2013) results; and the [MINER$\nu$A]{}neutrino [@minerva-nu-ccqe] (2013) and antineutrino [@minerva-antinu-ccqe] (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 single-differential cross section results are given in terms of [$Q^{2}_{\textrm{QE}}$]{}, the four-momentum transfer calculated from lepton kinematics under the quasi-elastic hypothesis, which is calculated using the equations: $${\ensuremath{E_{\nu}}\xspace}^{\mathrm{QE}} = \frac{2M'_{n}E_{\mu}-(M'^{2}_{n}+m_{\mu}^{2}-M^{2}_{p})}{2(M'_{n}-E_{\mu}+\sqrt{E_{\mu}^{2}-m_{\mu}^{2}}\cos{\theta_{\mu}})},
  \label{eq:ccqe-enuqe}$$ $${\ensuremath{Q^{2}_{\textrm{QE}}}\xspace}= -m_{\mu}^{2}+2{\ensuremath{E_{\nu}}\xspace}^{\mathrm{QE}}(E_{\mu}-\sqrt{E_{\mu}^{2}-m_{\mu}^{2}}\cos{\theta_{\mu}}),
  \label{eq:ccqe-q2qe}$$ where $E_\mu$ is the muon energy; $M_{n}$, $M_{p}$ and $m_{\mu}$ are the neutron, proton and muon masses, respectively; and $M'_{n} = M_{n} - V$ where $V$ is the binding energy of carbon assumed in the analysis[^3]. For both [MiniBooNE]{}datasets and the [MINER$\nu$A]{}neutrino dataset, $V = 34$MeV; for the [MINER$\nu$A]{}antineutrino dataset, $V = 30$MeV.

In the [MiniBooNE]{}analysis, [$Q^{2}_{\textrm{QE}}$]{}is calculated from the unfolded $T_{\mu}$ and $\cos\theta_{\mu}$ distributions. The [MINER$\nu$A]{}analysis unfolds the [$Q^{2}_{\textrm{QE}}$]{}distribution calculated with the reconstructed $p_{\mu}$ and $\cos\theta_{\mu}$ values. The errors on the [$Q^{2}_{\textrm{QE}}$]{}distributions for both experiments include the uncertainties relating to the muon reconstruction, so should cover the difference in the method used to produce the [$Q^{2}_{\textrm{QE}}$]{}cross section results. We note that the main results of our analysis use the [MiniBooNE]{}double-differential results only, so there is no possible tension from differences between the methods used to produce [$Q^{2}_{\textrm{QE}}$]{}distributions.

[MiniBooNE]{}neutrino
---------------------

The [MiniBooNE]{}CCQE data has been released as a double-differential cross section as a function of $(T_{\mu}, \cos\theta_{\mu})$, where $T_{\mu}$ is the kinetic energy of the outgoing muon and $\theta_{\mu}$ is the angle between the incoming neutrino and outgoing muon. Differential cross sections were also released as a function of [$Q^{2}_{\textrm{QE}}$]{}or ${\ensuremath{E_{\nu}}\xspace}^{\mathrm{QE \; RFG}}$, but the double-differential 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 shape-only 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 CCQE-corrected, and CCQE-like measurements. The CCQE-like sample is obtained by selecting events in which a muon was detected with no pions, but no requirement was made on the proton. The CCQE-corrected measurement is produced by subtracting background events (where the primary interaction was not CCQE) based on the NUANCE [@nuanceMC] generator prediction. The dominant background is CC1$\pi^{+}$, 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$\pi^{+}$ simulation included $\pi$-less $\Delta$ decay. The published signal purity for the neutrino dataset is 77%.

CCQE-like results are less model dependent than CCQE-corrected 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. CCQE-corrected results are used in this analysis. A downside of using the CCQE-corrected data is the explicit subtraction of $\pi$-less $\Delta$ 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 CCQE-corrected dataset to compare to their full models [@nieves_MB_nu_2011; @nieves_MB_anu_2013].

[MiniBooNE]{}antineutrino
-------------------------

The [MiniBooNE]{}antineutrino data has been released in the same format as the neutrino mode data. Again, the double-differential CCQE-corrected 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 $\nu_{\mu}$ contamination in the $\bar{\nu}_{\mu}$ beam, which is the largest background in the antineutrino CCQE sample ([MiniBooNE]{}is an unmagnetized detector). There is also a large CC1$\pi^{-}$ background, the analogue of the CC1$\pi^{+}$ contamination in the neutrino dataset. Two properties are used to measure the $\nu_{\mu}$ background [@AguilarArevalo:2011sz]: 8% of $\nu_{\mu}$-induced CC interactions produce no decay electron due to muon-nucleus capture; and the $\nu_{\mu}$-induced CC1$\pi^{+}$ events can be identified independently of $\bar{\nu}_{\mu}$-induced CC1$\pi^{-}$ as most $\pi^{-}$ mesons are absorbed. Unfortunately, this property makes CC1$\pi^{-}$ a bigger background to the CCQE analysis in antineutrino mode, and means that there is no sample with which to directly tune the CC1$\pi^{-}$ production from the NUANCE resonance model, so the neutrino mode CC1$\pi^{+}$ 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 CCQE-like sample is 61%, making the correction larger than for the neutrino mode sample.

[MINER$\nu$A]{}
---------------

The CCQE datasets from [MINER$\nu$A]{}are released as CCQE-corrected single-differential flux-averaged cross section as a function of [$Q^{2}_{\textrm{QE}}$]{}, where the flux has been averaged over the region $1.5 \leq {\ensuremath{E_{\nu}}\xspace}\leq 10$GeV. There is an additional requirement that $1.5 \leq {\ensuremath{E_{\nu}}\xspace}^{\mathrm{QE}} \leq 10$GeV, with ${\ensuremath{E_{\nu}}\xspace}^{\mathrm{QE}}$ as defined in Equation \[eq:ccqe-enuqe\]. Covariance matrices and central values have been released to perform fits to both shape-only 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 [MINER$\nu$A]{}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 CCQE-like sample. The published purity for the neutrino dataset ranges from 65% at low [$Q^{2}_{\textrm{QE}}$]{}to 40% at high [$Q^{2}_{\textrm{QE}}$]{}(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].

In the [MINER$\nu$A]{}CCQE analyses, the efficiency for selecting events with $\theta_{\mu} > 20^{\circ}$ is very low because the MINOS near detector, downstream of [MINER$\nu$A]{}, 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 phase-space. The [MINER$\nu$A]{}collaboration subsequently released a distribution where the cross section is measured for CCQE events with $\theta_{\mu} \leq 20^{\circ}$. As this dataset is less model-dependent, it has been used in the fits, and will be consistently used in this analysis. [MINER$\nu$A]{}also made cross-correlations 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 shape-only correlation matrices using the method given in Reference [@teppeiThesis]. The full matrix including both shape and normalization errors included is shown in Figure \[fig:cross-correlation\].

![Cross-correlation matrix including both shape and normalization uncertainties for the [MINER$\nu$A]{}neutrino and antineutrino samples. The eight neutrino and eight antineutrino bins shown here correspond to the eight [$Q^{2}_{\textrm{QE}}$]{}bins from the [MINER$\nu$A]{}datasets.[]{data-label="fig:cross-correlation"}](figures/MINERvA_cross-correlations.pdf){width="0.9\columnwidth"}

Monte Carlo prediction {#sec:mcpred}
======================

For each of the four experimental results included in the fit, one million CCQE and 2p2h events were generated with [<span style="font-variant:small-caps;">NEUT</span>]{}for each model using the default parameters given in Tables \[tab:RFG\_params\] and \[tab:SF\_params\] and the published flux for each dataset. The flux averaged cross section predictions were produced using the following method:

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

2.  Calculate the event rate by integrating the MC event rate histogram (flux $\times$ cross section).

3.  Integrate the published flux histogram to get the average flux.

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

5.  Divide the content of each bin by the bin width.

The default predictions for a variety of models available in [<span style="font-variant:small-caps;">NEUT</span>]{}, as well as the data, are shown in Figures \[fig:MIN\_20deg\_nominal\], \[fig:MB\_1D\_nominal\] and \[fig:MB\_2D\_nominal\] for the [MINER$\nu$A]{}, [MiniBooNE]{}single-differential and [MiniBooNE]{}double-differential samples, respectively. The Nieves 2p2h contribution is also shown on these plots for reference.

To produce a meaningful nominal $\chi^{2}$ for the [MiniBooNE]{}datasets, it is necessary to fit the [MiniBooNE]{}normalization parameters. The single and double-differential plots shown in Fig. \[fig:MB\_1D\_nominal\] and \[fig:MB\_2D\_nominal\] are scaled according to the [MiniBooNE]{}normalization parameter at the best fit point. The best fit values of the pull parameters $\lambda^{\mathrm{MB}}_{\nu}$ and $\lambda^{\mathrm{MB}}_{\bar{\nu}}$ are given in Table \[tab:MB\_nominal\_fits\]. Additionally, the nominal predictions for the [MiniBooNE]{}double-differential datasets, without the scaling factor applied, are shown in Figure \[fig:MB\_2D\_nominal\_UNSCALED\], which are easier to interpret by eye.

[0.9]{} ![Nominal model predictions for the [MINER$\nu$A]{}datasets with [$M_{\textrm{A}}$]{}$ = 1.01$GeV/$c^{2}$ and all other model parameters at their default values. The relativistic RPA model is shown.[]{data-label="fig:MIN_20deg_nominal"}](figures/MINERvA_numu_1D_20deg_nominal-models.pdf "fig:"){width="\textwidth"}

[0.9]{} ![Nominal model predictions for the [MINER$\nu$A]{}datasets with [$M_{\textrm{A}}$]{}$ = 1.01$GeV/$c^{2}$ and all other model parameters at their default values. The relativistic RPA model is shown.[]{data-label="fig:MIN_20deg_nominal"}](figures/MINERvA_numubar_1D_20deg_nominal-models.pdf "fig:"){width="\textwidth"}

[0.9]{} ![Nominal model predictions for the [MiniBooNE]{}single-differential datasets with [$M_{\textrm{A}}$]{}$ = 1.01$GeV/$c^{2}$ and all other model parameters at their default values. The relativistic RPA calculation is shown. Normalization parameters are applied as given in Table \[tab:MB\_nominal\_fits\].[]{data-label="fig:MB_1D_nominal"}](figures/MB_numu_1D_nominal-models.pdf "fig:"){width="\textwidth"}

[0.9]{} ![Nominal model predictions for the [MiniBooNE]{}single-differential datasets with [$M_{\textrm{A}}$]{}$ = 1.01$GeV/$c^{2}$ and all other model parameters at their default values. The relativistic RPA calculation is shown. Normalization parameters are applied as given in Table \[tab:MB\_nominal\_fits\].[]{data-label="fig:MB_1D_nominal"}](figures/MB_numubar_1D_nominal-models.pdf "fig:"){width="\textwidth"}

[0.8]{} ![image](figures/MB_nu_2D_nominal-models.pdf){width="\textwidth"}

[0.8]{} ![image](figures/MB_antinu_2D_nominal-models.pdf){width="\textwidth"}

[0.8]{} ![image](figures/MB_nu_2D_nominal-models_UNSCALED.pdf){width="\textwidth"}

[0.8]{} ![image](figures/MB_antinu_2D_nominal-models_UNSCALED.pdf){width="\textwidth"}

Note that the double-differential cross section plots shown in Figures \[fig:MB\_2D\_nominal\_UNSCALED\] have been rebinned. In the distributions released by [MiniBooNE]{}, and used in the fits, there are 20 $\cos \theta_{\mu}$ bins uniformly distributed between -1 and 1. For ease of presentation, these have been rebinned and the results are shown in eight $\cos \theta_{\mu}$ slices of varying sizes, where merged bins have been averaged and their errors combined in quadrature.

Fit procedure {#sec:fitprod}
=============

All minimizations are performed using the MIGRAD algorithm of the MINUIT package [@minuit], using the $\chi^{2}$ statistic: $$\begin{aligned}
  \chi^{2}(\vec{\mathbf{x}}) &= \left\lbrack\sum^{N}_{k=0} \left(\frac{\nu_{k}^{DATA}-\lambda_{\nu}^{-1}\nu_{k}^{MC}(\vec{\mathbf{x}})}{\sigma_{k}} \right)^{2} 
  + \left(\frac{\lambda_{\nu} - 1}{\varepsilon_{\nu}} \right)^{2}\right \rbrack \rightarrow \mathrm{{MiniBooNE\xspace}~\nu} \notag\\
  &+ \left\lbrack\sum^{M}_{l=0} \left(\frac{\nu_{l}^{DATA}-\lambda_{\bar{\nu}}^{-1}\nu_{l}^{MC}(\vec{\mathbf{x}})}{\sigma_{l}} \right)^{2} 
  + \left(\frac{\lambda_{\bar{\nu}} - 1}{\varepsilon_{\bar{\nu}}} \right)^{2}\right\rbrack \rightarrow \mathrm{{MiniBooNE\xspace}~\bar{\nu}}\notag\\
  &+ \left\lbrack\sum_{i=0}^{16} \sum_{j=0}^{16} \left(\nu_{i}^{DATA} - \nu_{i}^{MC}(\vec{\mathbf{x}})\right)V_{ij}^{-1}\left(\nu_{j}^{DATA} - \nu_{j}^{MC}(\vec{\mathbf{x}})\right)\right\rbrack \rightarrow \mathrm{MINER}\nu\mathrm{A}
\label{eq:chi2_def}\end{aligned}$$ where $\vec{\mathbf{x}}$ are the model parameters varied in the fit, $\sigma_k$ and $\sigma_l$ are the diagonals of the [MiniBooNE]{}shape-only covariance matrices for the neutrino and antineutrino results, $V_{ij}$ is the cross-covariance matrix provided by [MINER$\nu$A]{}, and $\lambda_{\alpha}$ and $\lambda_{\beta}$ are the normalization parameters for [MiniBooNE]{}neutrino and antineutrino datasets, with published normalization uncertainties of $\varepsilon_{\alpha}$ (10.7%) and $\varepsilon_{\beta}$ (13.0%)[^5].

Fits to individual datasets only include the relevant terms from the $\chi^{2}$ definition in Equation \[eq:chi2\_def\], and fits to single [MINER$\nu$A]{}datasets neglect cross-correlations (the summation is only over the relevant eight bins).

Parameter Goodness-of-Fit (PGoF) test {#sec:pgof}
-------------------------------------

Standard goodness-of-fit tests, such as the Pearson $\chi^{2}_{\mathrm{min}}$ 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 [@pgof-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 $\chi^{2}$, but contribute another degree of freedom. Thus the $\chi^{2}_{\mathrm{min}}$ 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 [MINER$\nu$A]{}); again, the key problem is a dilution of the $\chi^{2}$.

This problem is worsened in the case of datasets for which correlations between datapoints have not been included, where the $\chi^{2}/$DOF can be much less than 1, such as is the case for [MiniBooNE]{}. Looking at the Pearson $\chi^{2}_{\mathrm{min}}$ test statistic is not very illuminating when fitting to both [MiniBooNE]{}and [MINER$\nu$A]{}datasets.

The PGoF is a more rigorous test proposed in Reference [@pgof-2003] for fitting to global datasets and has been used extensively in sterile neutrino literature [@conrad12GF; @kopp13], 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 $$\chi_{\mbox{\scriptsize{PGoF}}}^{2}(\vec{\mathbf{x}}) = \chi^{2}_{\mathrm{tot}}(\vec{\mathbf{x}}) - \sum^{D}_{r=1} \chi^{2}_{r, \: \mathrm{min}}(\vec{\mathbf{x}}),
  \label{eq:pgof-test-statistic}$$ where $\vec{\mathbf{x}}$ are the parameters floated in the fits, $D$ is the number of datasets, $\chi^{2}_{\mathrm{tot}}$ is the minimum $\chi^{2}$ in a fit to all $D$ subsets of the data, and $\chi^{2}_{r, \: min}$ is the minimum $\chi^2$ obtained in a fit to the $r$th subset of the data. The PGoF test statistic forms a $\chi^{2}$ distribution with the number of degrees of freedom $$P_{\mbox{\scriptsize{PGoF}}} = \sum_{r=1}^{D}P_{r} - P_{\mathrm{tot}},$$ where $P_r$ and $P_{\mathrm{tot}}$ 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 $\chi^{2}_{\mathrm{min}}/{\ensuremath{N_{\mathrm{DOF}}}\xspace}<< 1$ issue.

The PGoF test still assumes that the datasets follow a $\chi^{2}$ 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 $p$-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.

Fit results {#sec:results}
===========

Fake data studies {#sec:ccqe_fake_data}
-----------------

The fitter was validated in two ways. Firstly, Asimov fake datasets [@asimov] 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 $\chi^{2}$ definition given in Equation \[eq:chi2\_def\] 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.

Combined fit {#sec:combined-fit}
------------

The results for the combined fits to all four datasets are given for both relativistic and non-relativistic RFG+RPA+2p2h models and the SF+2p2h model in Table \[tab:combinedSummary\]. The best fit distributions are compared with data for [MINER$\nu$A]{}in Figure \[fig:MIN\_20deg\_jointFit\], and for [MiniBooNE]{}in Figure \[fig:MB\_2D\_jointFit\]. 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 $\chi^{2}$ values, the contribution from that dataset to the $\chi^{2}_{\mathrm{min}}$ in the combined fit, and the total $\chi^{2}_{\mathrm{min}}$ in parentheses. Note that in Figure \[fig:MIN\_20deg\_jointFit\], the contributions from [MINER$\nu$A]{}are calculated for the individual datasets, which necessarily ignores cross-correlations and makes these numbers slightly misleading. Explicitly, $\chi^{2}_{\mathrm{M{\nu}A \: total}} \neq \chi^{2}_{\mathrm{M{\nu}A \: \nu}} + \chi^{2}_{\mathrm{M{\nu}A \: \bar{\nu}}}$ due to cross-correlations, so the values shown in the figure should be treated with caution.

It is clear from Figures \[fig:MIN\_20deg\_jointFit\] and \[fig:MB\_2D\_jointFit\] 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, $\chi^{2}_{\mathrm{MB}} \approx \chi^{2}_{\mathrm{M{\nu}A}}$. 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 $\chi^{2}$ 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 bin-to-bin correlations, the $\chi^{2}_{\mathrm{MB}}$ 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 [MINER$\nu$A]{}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.

[0.9]{} ![Comparison of the best fit from the combined fits detailed in Table \[tab:combinedSummary\] with the [MINER$\nu$A]{}datasets used in the fit. The $\chi^{2}$ values in the legend are the contribution from each dataset at the best fit point (and the total $\chi^{2}_{\mathrm{min}}$ for the combined fit).[]{data-label="fig:MIN_20deg_jointFit"}](figures/Comparison_MIN_1D_numu.pdf "fig:"){width="\textwidth"}

[0.9]{} ![Comparison of the best fit from the combined fits detailed in Table \[tab:combinedSummary\] with the [MINER$\nu$A]{}datasets used in the fit. The $\chi^{2}$ values in the legend are the contribution from each dataset at the best fit point (and the total $\chi^{2}_{\mathrm{min}}$ for the combined fit).[]{data-label="fig:MIN_20deg_jointFit"}](figures/Comparison_MIN_1D_numubar.pdf "fig:"){width="\textwidth"}

[0.8]{} ![image](figures/Comparison_MB_2D_nu.pdf){width="\textwidth"}

[0.8]{} ![image](figures/Comparison_MB_2D_antinu.pdf){width="\textwidth"}

Because of the large pulls on the [MiniBooNE]{}normalization parameters, shape-only fits were also performed (see Reference [@my_thesis] 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 [MINER$\nu$A]{}flux [@minerva_new_flux] results in an increase to the normalization of previous [MINER$\nu$A]{}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 shape-only 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 [@my_thesis].

PGoF results {#sec:pgof_results}
------------

Using the PGoF test defined in Section \[sec:pgof\], it is possible to test the compatibility between different subsets of the data. Tables \[tab:PGoF\_RFG\_non-relativistic\], \[tab:PGoF\_RFG\_relativistic\], and \[tab:PGoF\_SF\] 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 $\chi^{2}_{\mathrm{min}}$ test, where $\chi^{2}_{\mathrm{min}}$ is found by minimizing the function given in Equation \[eq:chi2\_def\], including only the terms for the relevant datasets. The PGoF test is found by subtracting $\chi^{2}_{\mathrm{min}}$ for each of the constituent datasets from the minimum of the combined dataset. The formulae for calculating the PGoF test statistic $\chi^2_{\mathrm{PGoF}}$ are given explicitly in Table \[tab:PGoF\_definitions\]. The $\chi^{2}_{\mathrm{min}}$ for each dataset is again determined by minimizing the function given in Equation \[eq:chi2\_def\] with only the relevant terms included.

In each fit, the [$M_{\textrm{A}}$]{}, 2p2h normalization, [$p_{\mathrm{F}}$]{}, and any [MiniBooNE]{}normalization terms are allowed to float.

One subtlety must be kept in mind when analysing the results in Tables \[tab:PGoF\_RFG\_non-relativistic\], \[tab:PGoF\_RFG\_relativistic\], and \[tab:PGoF\_SF\]: the PGoF test is only appropriate for statistically independent datasets. This makes the interpretation difficult for [MINER$\nu$A]{}, where cross-correlations are provided and used in the fits. Whenever a subset of data includes both [MINER$\nu$A]{}$\nu$ and $\bar{\nu}$ datasets, the fits include cross-correlations, but if only one dataset is included, they are not. This means that two of the rows in each table give slightly unreliable results: “[MINER$\nu$A]{}”, and “$\nu$ vs $\bar{\nu}$”. In each case, the $\chi^{2}$ function for the combined dataset includes cross-correlations, and the $\chi^{2}$ functions for the subdivided dataset does not. The issue is most obvious in Table \[tab:PGoF\_RFG\_relativistic\], where the “$\nu$ vs $\bar{\nu}$” row gives a negative PGoF $\chi^{2}$. 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.

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 $\chi^{2}_{\mathrm{min}}/{\ensuremath{N_{\mathrm{DOF}}}\xspace}$. The level of agreement given in the final column of Tables \[tab:PGoF\_RFG\_non-relativistic\], \[tab:PGoF\_RFG\_relativistic\], and \[tab:PGoF\_SF\] 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 [MINER$\nu$A]{}and [MiniBooNE]{}datasets (which include both neutrino and antineutrino samples) have the lowest level of agreement as shown by the “M$\nu$A vs MB” row. In contrast, the level of agreement between the neutrino and antineutrino datasets (which include the [MINER$\nu$A]{}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 \[tab:PGoF\_SF\] 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 non-relativistic RPA, shown in Tables \[tab:PGoF\_RFG\_non-relativistic\] and \[tab:PGoF\_RFG\_relativistic\], show much better compatibility between experiments than SF+2p2h. There is still a considerable amount of tension, which is largely due to differences between [MINER$\nu$A]{}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.

Rescaling parameter errors {#sec:scale-errors}
--------------------------

Assuming Gaussian statistics, 1$\sigma$ errors on a single fit parameter are defined by the parameter value for which $\chi^{2} = \chi^{2}_{\mathrm{min}} + 1$ [@pdg2014]. 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 \[tab:combinedSummary\]. 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 [@pumplin2000_pdfs; @stump_2001_pdfs; @collins_2001_pdfs]. A summary of the work of one PDF fitting group is given in Reference [@pumplin2000_pdfs] and was used as a guide here. Their solution for producing reasonable parameter error estimates is to inflate the value of the $\Delta \chi^{2}$ used to define the 1$\sigma$ parameter errors, although no generic solution is offered for defining that value. In the case of the PDF fits in Reference [@pumplin2000_pdfs], the $\Delta \chi^{2}$ used was very large, $\sim$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 $\chi^{2}$ 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 \[eq:pgof-rescaling\], where the value used to define the 1$\sigma$ error is given by $\Delta \chi^{2}$, and the rescaling parameter is given by $r$. $$\begin{aligned}
  \Delta \chi^{2} &= \chi^{2}_{\mbox{\scriptsize{PGoF}}}/{\ensuremath{N_{\mathrm{DOF}}}\xspace}\notag \\
  r &= \sqrt{\chi^{2}_{\mbox{\scriptsize{PGoF}}}/{\ensuremath{N_{\mathrm{DOF}}}\xspace}}
  \label{eq:pgof-rescaling}\end{aligned}$$ Note that this PGoF rescaling procedure does not modify the correlations between parameters, it simply rescales the error on each parameter.

There is some ambiguity over which PGoF statistic to use, the ‘All’ or ‘M$\nu$A vs MB’ row of Table \[tab:PGoF\_RFG\_relativistic\], with $\chi^{2}_{\mbox{\scriptsize{PGoF}}}/{\ensuremath{N_{\mathrm{DOF}}}\xspace}$ values of 17.9/6 and 15.9/3 respectively. The more conservative value is from the ‘M$\nu$A vs MB’ (because the greatest differences are between experiments, not between neutrino and antineutrino running), so this is used[^6]. To be explicit, we multiply the parameter errors from MINUIT by $r = \sqrt{15.9/3}\approx2.3$ based on this statistic, as shown in Table \[tab:final-errors\]. It can be seen from Table \[tab:final-errors\] that the 2p2h normalization is strongly suppressed and, even with the rescaled error, is nearly 3$\sigma$ 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 [@mb-ccqe-wroclaw; @mb-ccqe-2010; @mb-ccqe-antinu-2013], it is still significantly higher than the value of ${\ensuremath{M_{\textrm{A}}}\xspace}\simeq 1$ GeV$/c^{2}$ found by fitting to light target data and pion electroproduction data [@bbba05], and the inflated 1$\sigma$ error does not cover this difference.

Discussion of the fit results {#sec:discussion}
=============================

The results from the fits presented in Section \[sec:combined-fit\] show that none of the models which are currently available in [<span style="font-variant:small-caps;">NEUT</span>]{}describe all of the [CCQE]{}data adequately. In particular, there is a significant difference between [MiniBooNE]{}and [MINER$\nu$A]{}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 [$M_{\textrm{A}}$]{}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:2015awa] 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 [@bbba05]. 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 [$M_{\textrm{A}}$]{}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 \[tab:final-errors\]. This suppression is driven by [MINER$\nu$A]{}, which would completely suppress the 2p2h component of the model if [MiniBooNE]{}were not included in the fit [@my_thesis]. 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 [<span style="font-variant:small-caps;">NEUT</span>]{}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 CCQE-corrected cross section results from both [MINER$\nu$A]{}and [MiniBooNE]{}have part of the 2p2h signal region removed as a background ($\pi$-less $\Delta$ 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; @martini]. Future cross section measurements should be encouraged to focus on exclusive final states (CC$0\pi$) 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 [<span style="font-variant:small-caps;">NEUT</span>]{}is not equivalent to the full Nieves model because the 1p1h component in [<span style="font-variant:small-caps;">NEUT</span>]{}uses a global, rather than local, Fermi gas nuclear model. Using a global Fermi gas model results in greater interaction strength in the low [$Q^{2}$]{}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]{}$\chi^{2}$ 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 \[tab:PGoF\_SF\]), so should not be discounted completely.

The lack of reported [MiniBooNE]{}correlations and non-Gaussian 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$\sigma$ parameter errors cover the disagreement between the [MINER$\nu$A]{}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 bin-to-bin 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 non-relativistic RPA. The best fit parameters and uncertainties of the model are given in the second row of Table \[tab:final-errors\], and are correlated according to the matrix shown in Figure \[fig:final-correlation\].

![Correlation matrix for the best fit RFG+rel.RPA+2p2h model parameters.[]{data-label="fig:final-correlation"}](figures/RPA_ALL_2D_covar-noMBNORMs.pdf){width="0.9\columnwidth"}

Summary {#sec:summary}
=======

In this paper, we have shown how T2K’s NIWG uses previously published CCQE datasets from the [MiniBooNE]{}and [MINER$\nu$A]{}experiments to test CCQE+2p2h models in the [<span style="font-variant:small-caps;">NEUT</span>]{}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 ${\ensuremath{M_{\textrm{A}}}\xspace}= 1.15 \pm 0.03$GeV/$c^2$, the normalization on the 2p2h model $27\pm12$%, and ${\ensuremath{p_{\mathrm{F}}}\xspace}=223\pm5$MeV/$c$. 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 [@t2k_numubar_disapp_2015].

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 [@review]. The new models on the market are not perfect, and their implementation into [<span style="font-variant:small-caps;">NEUT</span>]{}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 [$M_{\textrm{A}}$]{}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 [MINER$\nu$A]{}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 $\pi$-less $\Delta$ decay from the CCQE-corrected 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.

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 [MINER$\nu$A]{}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 [<span style="font-variant:small-caps;">NEUT</span>]{}. We are grateful to T. Katori and P. Litchfield for helpful comments on the manuscript. We thank the J-PARC 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 HERA-B 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.

[^1]: This 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 [@multinucleon].

[^2]: The value of the axial mass used for the 2p2h contribution to the cross section was fixed to ${\ensuremath{M_{\textrm{A}}}\xspace}= 1.01$ GeV$/c^{2}$.

[^3]: Note that the binding energy $V$ is just the value assumed when calculating [$Q^{2}_{\textrm{QE}}$]{}, so we must use the same value as the experiments when producing comparable [$Q^{2}_{\textrm{QE}}$]{}distributions, but it need not be consistent with the binding energy used in our simulation.

[^4]: The 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 [@minerva-antinu-ccqe; @minerva-nu-ccqe].

[^5]: Note that the [MINER$\nu$A]{}normalization uncertainty is included in the covariance matrix, so also contributes a penalty term to the fit.

[^6]: It should also be noted that this rescaling procedure more than covers the difference between neutrino and antineutrino datasets.
