Inclusive and exclusive neutrino--nucleus cross sections and the reconstruction of the interaction kinematics

We present a full kinematic analysis of neutrino-nucleus charged current quasielastic interactions based on the Local Fermi Gas model and the Random Phase Approximation. The model was implemented in the NEUT Monte Carlo framework, which allows us to investigate potentially measurable observables, including hadron distributions. We compare the predictions simultaneously to the most recent T2K and MINERvA charged current (CC) inclusive, CC0$\pi$ and transverse variable results. We pursuit a microscopic interpretation of the relevant reaction mechanisms, with the aim to achieving in neutrino oscillation experiments a correct reconstruction of the incoming neutrino kinematics, free of conceptual biasses. Such study is of the utmost importance for the ambitious experimental program which is underway to precisely determine neutrino properties, test the three-generation paradigm, establish the order of mass eigenstates and investigate leptonic CP violation.


Introduction
The studies of neutrino-nucleus interactions are entering a new stage, motivated by longbaseline experimental programs, in which the statistical uncertainties will diminish and thus the nuclear effects -contributing to the systematical error -have to be kept well under control [1]. The incomplete theoretical knowledge of the neutrino-nucleus interactions influences various stages of experimental analysis. For instance, for the future Hyper-Kamiokande water Cherenkov detector [2], the method for reconstructing the neutrino energy will be mainly based on the kinematics of the outgoing muon, which is the only particle observed, assuming that the reaction-mechanism is two-body charged-current quasielastic (CCQE) dispersion on a bound nucleon. However, the energy range of the neutrino flux produced in the J-PARC facility [3] is such that other physical mechanisms give non negligible contributions to the cross section. In particular, multi-nucleon knockout processes (mainly driven by the excitation of two particle-two hole, 2p2h, components in nuclei) should be taken into account. Since in the latter processes the interaction takes place on a pair of nucleons, the energy balance is different than in the QE case, driven by the excitation of only one nucleon (1p1h). Mismatching the signal coming from these two reaction mechanisms would lead to a bias in the energy reconstruction [4,5]. It is therefore crucial to properly include the 2p2h channel into the Monte Carlo (MC) event generators.

CCQE 1p1h model
We developed a full exclusive CCQE 1p1h MC event generator based on the theoretical scheme developed in [6]. The model is capable to simulate both neutrinos (ν µ n → µ − p) and anti-neutrino 1 (ν µ p → µ + n) reactions for a variety of nuclei: C,O, Al, Ti, Fe and Ca. The original code [6,7] provided the total cross-section for a fixed neutrino energy value, and the CCQE differential cross sections depending on the energy and solid angle of the outgoing charged lepton. The code was included in the NEUT MC generator [8]. The implemented modifications in NEUT keep all the physics of the model: Local Fermi Gas (LFG) nucleondynamics, short and long range Random Phase Approximation (RPA) correlations, Pauli blocking, lepton Coulomb corrections,... and provide an almost fully exclusive cross-section by predicting the hadron kinematics in the first step (weak absorption of the gauge boson) of the reaction.

Implementation of the exclusive CCQE model in the MC
We implement the QE model of Ref. [6] in NEUT MC, so that we could extract both the position of the first step interaction and full hadron kinematics. The model is almost fully exclusive since we compute the cross-section as function of the: • radial position of the interaction in the nucleus.
• outgoing lepton momentum and angle.
• angle between the outgoing proton (neutron) and the transfer momentum direction.
With this information we can obtain the whole event kinematics applying conservation of momentum and energy: • lepton four-momentum p µ .
• final state nucleon four-momentum p µ N .
The model does not include the azimuthal angle of the final hadron with respect to the lepton reaction plane. This angle is generated with a flat probability in the [−π, π] interval.

Local Fermi Gas and nucleon kinematics: implementation of the removal energy correction
The present model utilizes a LFG to describe the nucleus, which provides on one hand a more accurate description of the Fermi momentum and Pauli blocking than those obtained in global FG approaches. On the other hand, it allows to locate the position of the first interaction inside the nucleus, which might affect the strength/relevance of the nuclear reinteractions. We will discuss interactions of neutrinos off carbon, which is the main target material for the most recent neutrino scattering experiments: NOvA, T2K, MINERvA and MiniBooNE. Figure 1 shows that in this nucleus, the interactions mostly occur between 1.5 and 4 fm. The LFG model introduces a relation between the Fermi momentum and the radial position given by the equation: p F n = (3π 2 ρ n (r)) 1/3 (2.1) with ρ n (r), the density of neutrons (protons for anti-neutrino reactions) for a given radial position, r, inside the nucleus. In our MC, we had chosen the neutron (hit nucleon) momentum to be taken between 0 and this local Fermi momentum. The neutron momentum as function of the radial position of the interaction is shown in Fig. 1. We can see that the highest local Fermi momentum is achieved at radius slightly above 1 fm for carbon. At first, the model relies on the Impulse Approximation (IA), where we consider the hit-nucleon is a plane wave state, and the momentum balance reads where p A−1 is the momentum of the remaining (A − 1)−nucleons at the moment of the collision. This momentum should cancel with the hit-nucleon momentum ( p N ) so the total Figure 2. Pictorial representation of the excitation energy scheme. When a nucleon is removed from a deep shell, the hole energy remains until the nucleus is de-excited. The FG removal energy refers only to Fermi level nucleons, which would correspond to leave the daughter (A − 1) nucleus in the ground state within this statistical nuclear model. It is an approximation to the experimental proton and neutron separation energies (S p,n ) of the target nucleus. Knocking out a nucleon from deeper levels will give rise to excited final-nucleus states, equivalent to the binding energy of the nucleon occupying this level ( ). We also show the accumulative occupation numbers when additional shells are considered. The energy scale depicted at the left of the plot has an arbitrary origin, and it is only intended to illustrate the energy differences between shells.
momentum of the initial nucleus vanishes. Within this approximation, we also consider that the momentum of the final state nucleus ( p A−1 ) is equal to the residual momentum of the initial nucleus and both cancel out ( p A−1 = p A−1 ). Thus the balance of Eq. (2.2) reduces to: which is the traditional equation of momentum conservation within the IA model. On the contrary, the IA is broken for the energy balance, due to the need of an energy contribution to transit from the ground state of the target nucleus to a new final nuclear configuration, with the daughter nucleus left in its ground or an excited state or even broken. Actually, the energy conservation equation reads where M A and M A−1 are the the ground state masses of the initial and final nuclei, and E ∞ N is the energy of the asymptotically observed nucleon (N ). In addition, A−1 > 0 is the excitation energy of the final nucleus (see Fig. 2). Finally, T A−1 is the final nucleus kinetic energy, which is very small and it is approximated to zero in what follows. Left: Number of events as function of the E IA m energy for neutrino scattering off 12 C within the LFG model. Both the 1p1h (Eq. (2.8)) and 2p2h (Eq. (3.3)) E IA m distributions are shown by the blue-solid and red-dashed lines, respectively. The gap between the two distributions is caused by the excitation energy of the two holes in the final state. As in Fig. 1, results have been folded with the T2K neutrino energy flux. Right: Probability to find a neutron in carbon with a momentum (p n ) for a given reaction missing energy (E IA m ) (see Eq. (2.8)) as predicted by the SF model [9][10][11] (contour plot) and for this implementation of the LFG (box plot). In both cases, the T2K flux [3] is used.
The energy balance in Eq. (2.4) does not apply to cases where any secondary rescattering collision changes the energy of the nucleon that emerges after the weak absorption of the gauge boson 2 .
The excitation energy can be estimated, in a first approximation, to be the energy of the hole, within the FG model 3 with T F the kinetic energy of the nucleon at the Fermi level for the given radial position, and T N the actual kinetic energy of the knocked out nucleon in the target nucleus. On the other hand, the experimental nucleon separation energy S N can be obtained from the masses of the initial and final nuclei: where M A and M A−1 are the ground-state masses of the initial (A Z ) and final [(A − 1) Z or (A − 1) Z−1 for neutrino or anti-neutrino reactions] nuclei, and m N the mass of the target 2 In fact, in these latter situations, the rupture of the daughter nucleus might occur and the analysis is more complicated. 3 Any LFG model implicitly assumes the existence of a mean-field potential U = −TF , which cancels in the difference of energies, and binds the nucleons. nucleon. Re-writing Eq. (2.4) using S N , we obtain: which reduces to the usual IA energy conservation formula, but with an additional correction: the term S N + T F , which is an approximation of the experimental removal energy and takes into account to some level the excitation of the final state nucleus. In case of the Relativistic Global Fermi Gas (RGFG), T F takes a constant value (∼ 27 MeV), but in the relativistic LFG model, it has a distribution depending on the radial position of the interaction. This dependency introduces several removal energies simulating a continuous distribution of excitation of nuclear states. This discussion of the IA energy balance fixes a problem with the nuclear missing energy (E IA m ), which appears within the traditional implementation in MC event generators of the relativistic LFG and RGFG models The value of E IA m becomes negative (non-physical) for some values of T N when, as it is common, the T F correction is not added 4 . Equivalently, the problem is caused by the wrong assumption of taking E N = (m N + T N ), instead of the correct expression E N = (m N + T N − T F ), which includes the mean field potential responsible for binding the nuclear system. The distribution E IA m of energies for a relativistic LFG is depicted in the left plot of The average E IA m is 28 MeV, very similar to the binding energies used in RGFG models (25 MeV) [14] or in MINERvA (27.13 MeV) [15]. T2K CC (left) and CC0π (right) inclusive double differential cross-section dσ/dp µ d cos θ µ . Data, taken from Refs. [12] and [13], respectively, are compared to the results obtained from the present implementation of NEUT with the numerical values for the χ 2 −likelihood test compiled in Table 3. Last bin accumulates all the statistics until the kinematic limit of 30 GeV.
We pay now attention to the two-dimensional (p N , E IA m ) distribution shown in the right plot of Fig. 3 following the same representation as in the Spectral Function (SF) scheme [9][10][11]. The carbon SF obtained in [10] is comprised of two contributions. The first one is determined by a mean-field description of the nucleus, while the second one takes into account two-nucleon short range correlations, and it is computed within a correlated basis function scheme in isospin-symmetric nuclear matter. The mean-field contribution of the SF modifies the dispersion relation by forcing a set of effective bound masses. This way the value of E IA m is constant for each of the nuclear levels with a broad momentum distribution, which is additionally distorted by the contribution of the correlated part of the SF. The model presented here is based on the LFG approach to the nucleus, where the dispersion relation is fixed to the on-shell target nucleon mass 5 , but with a Fermi level that depends on the spatial position through the local density. Despite its simplicity, the LFG distribution, as shown in Fig. 3, follows a pattern similar to that exhibited by the more realistic one inferred from the SF scheme of Refs. [9][10][11]. Nevertheless, some differences between both sets of predictions are visible in Fig. 3, in particular at the edges of the (p N , E IA m ) contours. This different dependencies might introduce distinctive differences in 5 The 1p1h contribution to the nuclear response function depends on the energy difference between particle-and hole-nucleons, where the mean-field potential −TF (r) cancels out. MINERvA CC (left) and CC0π (right) inclusive double differential cross-section d 2 σ/dp ⊥ dp . Data, taken from Refs. [16] and [17], respectively, are compared to the results obtained from the present implementation of NEUT with the numerical values for the χ 2 −likelihood test compiled in Table 4. the nuclear response when the nucleon target momentum is relevant such as in the case of low energy neutrino interactions and it might explain some of the disagreements discussed later in this work.

Predictions of the final state hadron kinematics
In the left panel of Fig. 4 we show the momentum distribution of the primary proton, created after the absorption of the gauge boson in neutrino processes, as predicted by the model presented in this work. In the right plot of the figure, we show the proton momentum correlated to the radial position of the primary interaction. Low energy protons are produced close to the outer surface of the nucleus having a large probability to survive nuclear re-scattering. The maximum momentum of the proton is limited by the energy of the neutrino, but the lowest values are determined by Pauli blocking, which is also function of the radial position of the interaction. The fact that the Pauli effects become less relevant at large radii (4 fm) allows the proton momentum to have values close to zero contrary to less sophisticated models such as the global FG [8].

Consistent implementation of the CC 2p2h model
In NEUT [8], the events are generated according to the distribution of the outgoing lepton, i.e. using the weight given by the value of the double-differential inclusive cross section, expressed as the contraction of lepton and hadron tensors, as given for instance in Eq. (2) of Ref. [7]. In addition to the 1p1h term, the hadron tensor W µν accounts also for 2p2h contributions evaluated following the LFG scheme of Ref. [7], which is fully consistent with the 1p1h implementation outlined above and based on [6]. It is computed, for (anti-)neutrino reactions, separately for proton-neutron and proton-proton (neutron-neutron) final states to provide isospin dependent final states. The location of the interaction vertex in the nucleus is chosen according to the density profile, and the initial state nucleons are picked below the Fermi level corresponding to the radial position following the LFG model recipe. The outgoing nucleons at the weak vertex are distributed according to the available phase-space. This is because all the hadron-variables are integrated out in the calculation of the inclusive lepton cross sections carried out in [7]. At this respect, note that in the recent re-computation of Ref. [18] some of these integrations have been undone, opening the possibility to improve on this phase-space prescription. The final state nucleons are generated uniformly in the center of the mass of the hadronic system and boosted to the laboratory rest frame. Next, their momenta are tested against the local Fermi level to implement Pauli blocking. This procedure neglects the dynamics of underlying nuclear model and produces a symmetric distribution of outgoing nucleons [18]. The produced pair of nucleons is fed into the NEUT cascade model accounting for the transport of nucleons in the high density nuclear medium. We have introduced the same removal energy corrections than in the case of the 1p1h, but we need to take into account the fact that two nucleons are removed from the nucleus. This provides an energy balance equation: where we have neglected the kinetic energy of the daughter nucleus, have contemplated the possibility of different (isospin) Fermi levels for the hit nucleons N 1 and N 2 , and MINERvA [15] (left) and T2K [19] (right) δα T distributions compared to the predictions of the implementation of NEUT in this work. The simulation has been done for three configurations of the NrSP (nominal 50% and 150% strengths) and the obtained absolute χ 2 −values are compiled in Tables 3 and 4.
Finally, N 1 and N 2 are the two outgoing nucleons asymptotically observed. In the left plot of Fig. 3, we showed the distribution of missing energies from the 2p2h mechanism and compared to the 1p1h E IA m values discussed in the previous section. The excitation of two nucleons leads to a bigger offset of the E IA m values and a longer energy tail compared to that of the 1p1h distribution. The immediate consequence of this implementation is that the average E IA m is around 45 MeV which is larger than the one of the CC1p1h and larger than previous implementations of CC2p2h models, where the typical CC1p1h E IA m were implemented. This correction will reduce the overall 2p2h cross-section for low energy transfers.
NEUT MC model and the secondary nuclear collisions The remaining neutrinonucleus interaction channels, mostly with associated pion production, are simulated based on the existing NEUT [8] Monte Carlo event generator. The resonant pion production is based on the Rein-Sehgal model [20], taking into account eighteen resonances with masses below 2 GeV and their interference terms, with the axial mass fixed to M RES A = 1.21 GeV. This model has been compared to experimental data by the T2K collaboration showing remarkable agreement [14]. Neutral current and charged current coherent pion production is simulated using the Rein-Sehgal model in Ref. [21]. The CC coherent pion production includes PCAC (partially conserved axial-vector current) and lepton mass corrections, as discussed in [22]. DIS processes are simulated using the GRV98 [23] parton distribution, with low-Q 2 corrections from the Bodek and Yang model [24]. Secondary interactions of pions inside the nucleus are simulated using an intra-nuclear cascade model based on the method developed by Salcedo et al. [25], tuned to external π− 12 C data [26].

Comparison to experimental data
In this section we discuss the comparison of the predictions from the model with recent data from MINERvA and T2K cross sections with no pions in the final state. The implementation of the model inside NEUT allows us to make a direct comparison with the experimental cross-sections since all the interaction channels are considered including, the transport of the nucleons and pions inside of the nucleus. The data selected include inclusive muon kinematics and transverse variables to explore the limits of the hadron kinematic predictions of the model.

Event Simulation and data selection
Events are simulated using the NEUT package with the CCQE and CC2p2h reactionmechanisms described above. We take the fluxes from the experiment releases according to their best understanding. The simulation is done for three configurations of the nucleon rescattering probability (NrSP): nominal, and 50% and 150% strengths. This is only applied to the proton re-scattering while pions are kept to their nominal NEUT values. We select events according to the particles emitted by the nucleus after the interaction taking into account the event acceptance of the experiments as described in their published designs.

T2K data samples
The neutrino T2K data-sample has different selections. The CC inclusive measurement considers only the muon production kinematics ignoring all hadronic activity [12]. The T2K selection criterion for CC0π [13] requires no charged or neutral pions in the final state. Based on this selection, T2K provides double differential cross-sections for carbon and oxygen nuclear targets. The comparison is performed only on carbon data to keep a common nuclear target across the different measurements and experiments.
The T2K selection criterion for the single transverse variables analysis of the CC0π1p sample requires the detection of a muon and proton with the following conditions [19]: • cos θ µ > −0.6 with θ p the polar angle of the outgoing proton, which has a modulus of the three-momentum | p p |. T2K presents also a slightly less restrictive selection criterion for CC0π1p requiring the detection of a muon and proton with the following conditions [19]: • 0.45 GeV < | p p | < 1.0 GeV.

MINERvA data samples
MINERvA published different event selections: the CC inclusive, the CC0π that removes events with detected charged π's and electromagnetic activity to eliminate events with π 0 , and the CC0π1p sample which is a sub-sample of CC0π requesting the presence of an identified proton in the final state. Another data sample from MINERvA, the so-called available energy [27], requires full simulation of the detector simulation that falls beyond the capabilities of this work.
The CC inclusive selection in MINERvA requires [16] a muon with polar angle angle (θ µ ) below 20 o . The selection criteria for CC0π requires [17] in addition: • a muon with momentum | p µ | between 1.5 GeV and 10 GeV, • neither charged, nor neutral pions escaping the nucleus.
• no γ with energies about 10 MeV. We check this cut actually do not affect the MC prediction for CCQE and CC2p2h.
The selection criteria for CC0π1p requires in addition [15]: • a proton with polar angle (θ p ) below 70 o .
• a proton with momentum between 0.45 GeV and 1.2 GeV.
These two conditions are applied to protons leaving the nucleus after the nucleon re-scattering (NrS).

Target composition
MINERvA and T2K targets are composed materials made of several components: CH, O, Al,.. Proportions in weight and nuclear content are given in Tables 1 and 2. To simulate the experimental composition we take the approximation of selecting only the main 3 components: C,H and O in proportions given in the tables. This actually has an appreciable effect on the selection mainly because of the large Fermi momentum and nuclear radius of the oxygen affecting the nucleon transport in the nucleus, accounting for secondary collisions. Both experiments have the same target material (plastic scintillators) with very similar composition. The correction to the cross-section prediction introduced by be oxygen contribution is estimated by our models to be at the order of a percent. The exception to this treatment is the CC0π cross-section that is reported by T2K for a pure carbon target [13] and not in hydrocarbon as for the CC inclusive and MINERvA.

Inclusive cross-sections
In this section we compare our model predictions with the CC inclusive and CC0π measurements. The selection of events in the first reaction relies on the muon kinematics and ignore the hadronic component of the interaction. This sample allows to make a datato-model comparison with a reduced bias from selection and detector acceptances, but it relies on the proper description of other interaction channels in the Monte Carlo. The CC0π sample requires that there are no pions observed in the final state. This selection, reduces the contribution of channels beyond CC1p1h and CC2p2h but it is affected by the correct modelling of the pion re-interactions inside the nucleus. T2K published both the CC inclusive [12] and the CC0π [13] cross-sections as a double differential distributions as function of the muon momentum and angle. The comparisons of these data-samples and the NEUT model predictions are shown in the left and right panels of Fig. 5, respectively. The CC inclusive spectrum shows a sizeable contribution from resonant and deep inelastic scattering (DIS) channels, while for the CC0π such components are considerably reduced, see also Table 6. The remaining CCRes contribution comes from pion absorption in the nucleus. The model predicts reasonably well the tendencies in the experimental data.
MINERvA published the muon 2D longitudinal and transverse momentum distributions for the CC inclusive cross-section measurements [16] (see left panels of Fig. 6). The cross-section is dominated by DIS, but there are regions (1.5 GeV < p < 5.0 GeV) where the QE and 2p2h contributions are relevant. The agreement with the present version of NEUT in the regions where CCQE is relevant is good, while the high longitudinal momentum histograms show cross-section predictions lower than data. The CC0π sample from MINERvA is also shown in Fig. 6 (right panels). It can be observed that the predictions are qualitatively similar to those found for the MINERvA CC inclusive case, lower than the data for large (p > 5.0 GeV), but also for low (p < 3.5 GeV) longitudinal momenta.
To estimate the agreement data-model, we computed the χ 2 −merit function between the model predictions and data using the full error covariance matrices provided by the T2K and MINERvA experiments. The results are compiled in Tables 3 and 4, respectively. The comparison is done for three different scale factors of the proton-nucleon cross-sections accounting for secondary collisions. Absolute χ 2 −values for the T2K CC-inclusive reaction are close to the dof and they show very little variations between the three examined scenarios. This is because re-scattering effects should not affect the inclusive cross section, and the observed differences should be produced by MC fluctuations (note that the three simulations are statistically independent). The χ 2 −value found in this work is almost a factor of two smaller than those presented in the experimental paper. The theoretical description of the T2K CC0π data sample is better than that achieved for the inclusive Figure 11. MINERvA δp Tx (top) and δp Ty (bottom) distributions [28]. Left and right panels show the comparison (details as in Fig. 8) with NEUT in linear and logarithmic scales, respectively. one, with χ 2 −value around 17 for 29 degrees of freedom. In this case, the different NrSP assumptions do not practically alter the results, as expected since the kinematics of the proton is not used in the CC0π event selection. The χ 2 value achieved with the present scheme is among the best reported in the T2K paper and it approaches the model "NEUT 5.4.1 LFG" in Table V of [13], since both largely share the physics implementation.
On the other hand, the best-fit χ 2 values for the CC inclusive MINERvA data sample, collected in Table 4, are among the best five reported in the experimental paper. Since this is dominated by the DIS cross-section, as discussed in Fig. 6, this result also confronts the prediction for this reaction channel in NEUT.  Table 4. Same as in Table 3, but for MINERvA experiment.
quantitatively slightly poorer than that seen above for T2K. Nevertheless, the χ 2 value is among the best two reported in the experimental paper and significantly better than any of the models including 2p2h contributions. NEUT predictions for MINERvA CC0π show apparently a larger, though still soft, dependence on the proton re-scattering probability.

Transverse variables
Transverse variables [29,30], depicted in Fig. 7, lead to observable distributions with minimal dependence on neutrino energy, which provide direct constraints on nuclear effects in (anti-)neutrino-nucleus interactions. Obviously, the measurement of the distribution of these variables require the observation of the outgoing proton (neutrino reactions) or the neutron (anti-neutrino reactions) in the final state. The missing transverse, with respect Figure 12. Inclusive | p µ | (left) and θ µ (right) distributions from the CC0π1p MINERvA sample [15]. Details of the comparison with NEUT results as in Fig. 8. Figure 13. MINERVA CC0π1p differential cross sections [15] in proton kinematics, momentum (| p p |) and angle (θ p ), together with results from the current implementation of NEUT. Details of the comparison with NEUT as in Fig. 8.
to the neutrino direction, momentum is computed as In addition, the missing transverse momentum is separated in two components in relation to the ν − µ reaction plane, defined by the neutrino and the emitted charged lepton. One MINERvA CC0π1p differential cross sections in reconstructed | p n |. Data from Ref. [15]. We also show results from the current implementation of NEUT, with details of the comparison as in Fig. 8.
of the transverse components is contained in the reaction plane δp Ty , while the another one, δp Tx , is perpendicular to this plane.
On the other hand, the transverse angular variables read: where p T = p Tµ + p Tp , with p Tµ (or p ⊥ , as used in Fig. 6) and p Tp the transverse projections of the muon and proton momenta to the neutrino direction. This discussion focuses on the QE-like process ν µ +A → µ − +p+X, where X is a final-state hadronic system consisting of the nuclear remnant with possible additional protons but without pions that indicate resonant or other processes. There is an imbalance, δ p, between the initial neutrino momentum and the sum of final-state lepton and hadron momenta as a result of nuclear effects. Under the assumption that X is just the remnant nucleus, with (A − 1) nucleons, then |δ p | gives the magnitude of its recoil momentum, which can be obtained [31] independently of the unknown incident neutrino energy. Moreover, assuming perfect balance of momentum in the reaction (Eq. (2.3)), |δ p | can be identified to the neutron target momentum, which is then given in terms of measurable quantities [15] | where p Lp and p Lµ (or p , as used in Fig. 6) denote the projections of the corresponding three-vectors on the direction of the incoming neutrino. For the MINERvA measurement of Ref. [15], M A is the mass of the carbon target, while for M A−1 is taken (M A − m n + E b ), with E b the nucleon binding energy that is fixed to 27.13 MeV and m n is the mass of the neutron 6 . The T2K collaboration considered additional distributions based on the neutrino energy reconstruction formula used in neutrino oscillation experiments. This prescription assumes a genuine QE event where the target nucleon is at rest, and the neutrino energy is reconstructed assuming energy and momentum conservation (see for instance Ref. [4]): • Proton momentum balance (Inferred variables with the second selection criterion): |∆ p |, ∆θ and ∆| p | in bins of muon momentum and angle.
Experimental and NEUT predictions for the transverse angular variables δα T and δφ T , and the missing transverse momentum |δ p T |, both for MINERvA (left) and T2K (right) are shown in Figs. 8, 9 and 10, respectively. Absolute χ 2 −values obtained from the comparison with data are compiled in Tables 3 and 4, respectively. The NEUT results for transverse variables are in general in an acceptable agreement with both MINERvA and T2K data, though the MINERvA distributions are better described. As expected, we observe some significant dependence on the proton re-scattering probability, which is reflected in the χ 2 −likelihood tests.
The transverse momentum components contain different information. The component δp Tx is expected to be symmetric around zero, with a width which depends on the target neutron momentum and further re-scattering effects, while δp Ty can be asymmetric due to leading effect of the neutrino boost. Results are shown in Fig. 11. The tendency is very well described by the model with the long tails dominated by 2p2h, resonant and, DIS mechanisms in the δp Tx (δp Ty ) distribution. The χ 2 comparison, see Table 4 shows an excellent agreement with a preference for an increase in the NrSP, which reduces the contribution of the CCQE channel by reducing the probability of proton tagging in detectors. This tendency is shared by most of the other transverse observables.
MINERvA collaboration also reported inclusive | p µ | and θ µ distributions from its CC0π1p MINERvA sample [15]. The comparison of these data with the current implementation of NEUT is shown in Fig. 12. We find a quite good description of these two event distributions, with χ 2 /dof around one (see Table 4) for nominal NrSP, and some dependence on this latter input as expected when analysing CC0π1p data-samples.
The experimental results with a visible proton in the final state are biased towards high momentum transfers since the proton should have at least 450 MeV to be detected. The typical proton momentum and angle distributions in MINERvA are shown in Fig. 13. On the contrary, the samples with no visible protons do not suffer from large momentum transfer biases. The difference between the two are dominated by low | q | contributions with feed-down background cause by NrSP. These tendencies can be observed in Figs. 22 and 23 of Appendix A, where the model estimation for the energy (q 0 = E ν − E µ ) and momentum (| q | = | p ν − p µ |) transfer distributions for the MINERvA and T2K CC inclusive, CC0π and CC0π1p event selections are shown.
The overall agreement with NEUT is good, showing the importance of 2p2h mechanisms. The results for transverse variables is good for the MINERvA data, with statistically acceptable values of χ 2 /dof for most of the cases. The worst comparison is obtained for the reconstructed | p n | variable, where a large discrepancy is observed in the region around 0.3 GeV (see Fig. 14). This is at the transition from the CC1p1h dominated cross-section to the one dominated by resonance and CC2p2h mechanisms. This is actually the most distinctive difference in all the comparisons of this work and a nice reference observable to try model variations. As it is is shown in Fig. 14 and Table 4, the variation of the proton re-scattering probability does no alleviate the discrepancy. Since, the high momentum Figure 15. T2K CC0π1p |∆ p | distribution. The panels correspond to different muon kinematic bins. From left to right and up to down: −1 < cos θ µ < −0.6, −0.6 < cos θ µ < 0 with | p µ | < 250 MeV, −0.6 < cos θ µ < 0 with | p µ | > 250 MeV, 0 < cos θ µ < 1 with | p µ | < 250 MeV, 0 < cos θ µ < 0.8 with | p µ | > 250 MeV, 0.8 < cos θ µ < 1 with 250 MeV < | p µ | < 750 MeV and 0.8 < cos θ µ < 1 with | p µ | > 750 MeV. Data taken from Ref. [19].
(| p n | ≥ 0.5 GeV) region is well reproduced by the model, the discrepancy seems to be led by the transition, either from non described tails in the CC1p1h, which might come from high energy neutron target components predicted by realistic SFs, or by a miss-representation of resonant or CC2p2h models. In any case, it seems that a re-weight of the cross-section will not improve the agreement.
The agreement with T2K CC0π1p data is less impressive and the obtained χ 2 /dof , see Table 3, are large for all the three observables |∆ p |, |∆θ| and ∆|p| reported in [19], and shown here in Figs. 15-17. The worst situation is found for the ∆θ distribution, with the largest contributions to χ 2 produced by the negative ∆θ bins. There, the number of events is always very small independent of the muon-kinematics, and the present model fails to properly describe those data, though one should bear in mind that these bins have a negligible weight in the totally integrated cross section. Indeed, if these bins are removed, the merit-function χ 2 is reduced to 66 for 30 degrees of freedom. The agreement is slightly better for the re-scaling factor NrSP=0.5 that reduces the scattering of the outgoing proton in the nucleus (54 for negative values of ∆θ). The χ 2 figures for |∆ p |, and ∆|p| are also reduced with the re-scaling factor NrSP=0.5. The difference between the agreement found for MINERvA and T2K data samples might point to an energy dependent deviation. The LFG model might provide a better approximation to the MINERvA energies than to the T2K ones, which might be more sensitive to finer details of the low-lying nuclear levels. Figure 16. T2K CC0π1p |∆θ| distribution. The panels correspond to the muon kinematic bins specified in Fig. 15. Data taken from [19].

Integrated cross-section
We have performed the numerical integrals of the differential distributions for the six event samples examined in this work. The obtained cross sections, both from data and from the model predictions are compiled in Table 5. The results show a common tendency of the NEUT model to predict lower cross-sections for the inclusive and the CC0π samples, while its predictions are marginally larger for the CC0π1p data-set. The deviations with 1.39 ± 0.13 1.46 −4.6 ± 9.0 Table 5. Integrated cross-sections for the different event selections provided by T2K and MIN-ERvA. The CC0π1p integrated cross-section has been computed using the dσ/d(δα T ) differential distribution.
the theoretical approach for both experiments are similar, except for the case of the CC inclusive. This is expected due to the DIS cross-section, and the very different proportion predicted for this channel for the MINERvA and T2K experiments, see Table 6. The change in the Exp to Model ratio observed from the CC0π and the CC0π1p may be a consequence of both the proton momentum prediction below the detector detection thresholds and the proton NrSP. The effect of the NRsP on the integrated CC0π1p cross-section is shown in Table 7. The number of events with visible protons increases when reducing the NrSP and vice-versa. The effects on MINERvA are slightly reduced with respect to those found for T2K due to the larger proton momentum expected at higher neutrino energies. Even if a reduced NrSP choice will bring the numerical values closer to the differences seen for the CC0π samples, we should be cautious. The effect of the CCRes model is not apparent in the results, but however, 23% in MINERvA and 12.5% in T2K of the CC0π1p events come from CCRes according to our model. The modelling of the CCRes should also take into account the absorption of the emerging pions by the nucleus. If, on the contrary, we assume that the CCRes is well simulated, the results point to a deficit in the model prediction of low momentum protons (≤ 450 MeV). A larger re-scattering probability would approach the two results. Another possible cause of the discrepancy is the larger | q | values intrinsic to CC0π1p with respect to CC0π since the proton should be emitted with momentum greater that 450 MeV. Nevertheless, given the experimental uncertainties also included in Table 5, we conclude that the present model implemented in NEUT leads to reasonable integrated cross sections for all six samples considered in this work.

Data vs theoretical predictions in terms of the scaling variable
The scaling variable ψ has been shown to be a powerful tool to learn details on the neutrinonucleus interaction [32]. This variable, proposed originally for electron scattering [33,34], has been adapted to neutrino-nucleus interactions recently. The scaling variable (ψ ) is  Table 7. T2K and MINERvA CC0π1p total cross-section obtained using NRsP (proton rescattering probability) configurations of 50% and 150%. As in Table 5, we use the dσ/d(δα T ) differential distribution to perform the numerical integrations. In percentage, we also show the relative variation with respect to the nominal Model. defined as: where, in addition to variables already introduced, p F is the Fermi momentum and E shift is a energy shift, which are fixed in this discussion to 228 MeV and 20 MeV, respectively. The scaling variable ψ is not accessible to traditional neutrino experiments, since the neutrino energy is not measurable in an event by event basis to compute momentum and energy transfers. Neutrino experiments usually report their flux-averaged cross-section results as function of muon-kinematics bins 7 , dσ Exp /d p µ , which can be compared with  Figure 19. The same as Fig. 18, but for the MINERvA CC0π cross-section. theoretical MC differential distributions for the same binning, On the other hand within the theoretical MC model, one can associate each event with a value of the scaling variable ψ (in general, values of the scaling variable comprised in a certain bin, since both neutrino energy and muon-kinematics are binned). We propose  Figure 21. The same as Fig. 18, but for the T2K CC0π cross-section.
to express the ratio of data to theoretical MC predictions (data/MC, in what follows) as function of ψ . Thus, we define the average data/MC ratio as: with N (ψ ) the number total of MC generated events, and N MC ( p µ |ψ ) the number of events with muon kinematics comprised in the bin around p µ , which gives rise to the value ψ for the scaling variable. Thus, f ( p µ |ψ ) is the fraction of events predicted by the theoretical model for a given value of ψ and muon kinematics p µ , taking into account the neutrino energy spectrum of the experiment. The distribution R(ψ ) is built in such a way that the ratio between data and the theoretical MC results are weighted according to the population f ( p µ |ψ ). In the limit in which f = 1 for one kinematics-bin and there is no more than one value of ψ contributing to this bin, R is no more than the ratio of data to MC for this given value of ψ . This is not the case in most of the experimental bins, but we still expect that some of the deviations from measurements are accumulated in the corresponding value of the scaling variable.
In addition, the number of events N MC ( p µ |ψ ) can be split into the different mechanisms (1p1h1, 2p2h, Res, DIS,...) considered in the theoretical approach implemented in the MC.
The errors on R(ψ ) and the covariance matrix for different values of ψ can be propagated using the above definition and the covariance matrix from the experiments. We have applied this algorithm to T2K and MINERvA inclusive and CC0π cross-sections.

5.1
Results for R(ψ ) Minerva CC inclusive R(ψ ) is shown in the left panel of Fig. 18, where the contributions of the different reaction channels are also shown. The total distribution shows a remarkable pattern, with the ratio almost constant, independent of ψ , around 1.1 (i.e. the MC prediction is some 10% lower than the data). This means that the proportion of the different reaction channels is well balanced. The 2p2h contribution is very small and it is difficult that this mechanism could significantly influence the overall picture, but both 1π and DIS channels smoothly balance above ψ = 1. We show the R−correlation matrix in the right plot of Fig. 18. The correlation is larger than 90% for any pair of of ψ values. This is a consequence of both the initial correlation between experimental results and the fact that several ψ −values contributes to the same muon-kinematics bin. We repeat in Fig. 19 the exercise with the inclusive CC0π MINERvA data. The tendency observed is similar to that seen for the CC inclusive. MC predicts smaller crosssection by a similar amount (10%), though in the CC0π case some structure is observed with a small deficit of MC below ψ = 0. On the contrary, the region above ψ = 2 shows an excess on the MC predictions with respect to data. This region is dominated by 2p2h and 1π (Res) channels. It is important to notice that these results do not call for a large modification of the 2p2h contribution as requested by the calorimetric measurements in MINERvA [27]. The bin to bin correlation is shown in the right panel of the figure. The observed correlations are smaller than in the previous case, with some regions below 70%. This might explain the appearance of some structure in the CC0π ratio plot.
The R(ψ ) results for T2K CC inclusive and CC0π samples are shown in Figs. 20 and 21, respectively. The observed patterns in both cases are similar to those discussed above for the MINERvA experiment. Our model predicts smaller cross-section than the measurements by about 20% (R(ψ ) ∼ 1.2), while the distribution of the ratio as function of ψ is again rather flat. There is some reflection of the CCRes contribution in Fig. 20 around ψ ∼ 1.5, which might be an indication of an even smaller single pion cross-section. Similar low cross-section predictions were also observed in the T2K CC1π results [14]. The impact of the possible CC1π miss-modelling is minimised in the case of CC0π where the contribution of CC1π becomes as the same level than CC 2p2h, see Table 6. The region around the CCQE peak shows a flat dependency with the scaling variable contrary to the MINERvA CC0π results which exhibit a visible decrease below ψ < 0. It is also observed that the prediction is balanced between the different reaction channels. The correlation matrices are shown in the right panels of Figs. 20 and 21. The correlations for both CC inclusive and CC0π are smaller than in the case of MINERvA, 60% correlation between values of ψ below 0 and higher than 1.5. This reduced correlation gives more credibility to the tendencies of R(ψ ).
The large errors (≈ 10% for T2K and ≈ 5% for MINERvA) are the consequence of To understand the large correlations across all values of ψ both for T2K and MIN-ERvA, we estimated the correlation removing all the off-diagonal terms in the covariance matrices provided by the experiment. In this case, the ratios R(ψ ) for different values of the scaling variable can be correlated only because receive contributions from the same muon-kinematics bin for different neutrino energies, contained in the non-monochromatic beam. The smallest observed correlation is reduced from 90% to 40% in the case of MIN-ERvA and from 60% to 5% in the case of T2K. The better figure in CC0π T2K sample might be a consequence of the the lower non-CC0π backgrounds, the narrow neutrino beam at T2K, but also from the fact that the low energy will allow to separate better the regions dominated by 1p1h and the resonant and DIS components. The experimental correlations errors come mostly from flux uncertainties, statistical and systematic experimental errors but also from bin to bin migrations in the extraction of the cross-section. Improvement in flux determination, larger statistics to reduce the bin size and select the proper data representation might improve the conclusions of this study.

Conclusions
We have presented an exclusive final state model to describe CC1p1h interactions. The approach is based on a LFG picture of the nucleus and uses a consistent implementation of the removal energy, that provides an estimation of the excitation of the final nuclear system. The model has been included in NEUT to profit from the existing simulation of 2p2h, pion production and DIS mechanisms and on the transport simulation of the hadrons inside the nucleus after the interaction. Predictions are simultaneously compared to the most recent T2K and MINERvA inclusive, CC0π and transverse variable results, showing an acceptable agreement with the data from both experiments. Results from T2K suffer from low statistics, but they also show worse agreement with the model predictions. This might be an indication of some energy dependency that is not properly accounted by this implementation. The correct modelling of the energy removal reduces the amount of interactions at low q 2 and the total CC1p1h cross-section facilitating the agreement with the experimental results. On the other hand, the overall good description of the MINERvA transverse variables found here and, in general of its CC0π data-sample, does not support a large modification of the 2p2h contribution as requested by the calorimetric measurements of that collaboration.
We have also proposed a novel comparison between flux-folded data and MC theoretical predictions accumulated in bins of the scaling variable Ψ , which can be used to signal possible deficiencies of theoretical schemes.
A microscopic interpretation of the relevant reaction mechanisms becomes essential in neutrino oscillation experiments in order to achieve a correct reconstruction of the incoming neutrino kinematics, free of conceptual biasses. Studies, as the one presented in this work, are of the utmost importance for the ambitious experimental program which is underway to precisely determine neutrino properties, test the three-generation paradigm, establish the order of mass eigenstates and investigate leptonic CP violation.

A Appendix: Comparison of momentum and energy transfer distributions
The different experimental selection criteria might bias the momentum and energy transfers. The comparison of the different accessible (q 0 , | q |)−phase space provide some indications of possible deficits in the models. The model estimation for the energy and momentum transfer distributions for the MINERvA and T2K CC inclusive event selections are shown in Fig. 22. The cut off in | q | implemented in the CC2p2h model is clearly visible in the MINERvA distributions. The CC2p2h cross-section at the cutoff | q | = 1.3 GeV represents around 10% of the total one for this momentum transfer. Although, we find similar drops for MINERvA and T2K, there are more contributions to the cross-section in the first experiment above | q | = 1.3 GeV, and hence we expect that the implementation of this cutoff to make a bigger impact in the total cross-section determination for MINERvA. As expected, the resonant, 2p2h and DIS contributions are larger for the MINERvA energies, see Table 6.
The predicted energy and momentum transfer distributions for the MINERvA and T2K CC0π event-selections are shown next in Fig. 23. As in the previous case, the cut off | q | ≤ 1.3 GeV for the CC2p2h contribution is clearly visible, both for MINERvA and T2K distributions, with the drop at the cutoff amounting around 10% of the total cross-section As expected, the CC0π event-sample has smaller contamination from DIS and resonant processes, which leads to a small bias in the event selection. Moreover, the CC1p1h contributions in this data selection show a dependence on q 0 and | q | similar as that observed for the CC inclusive ones.
The q 0 and | q | event distributions used for the transverse variable analysis (CC0π1p) are shown in Fig. 24. The significant reduction of the CCRes and CCOther contributions is evident in the plots and it can be also seen in Table 6. The detector acceptance cut employed for the MINERvA CC0π1p selection of events is observed as a change of slope around q 0 = 0.6 GeV, clearly visible in the 1p1h distribution. In the case of T2K, the hard cutoff is not visible and the distributions of momentum and energy transfers are narrower and shifted towards lower values, with smaller contamination from Res and DIS-Others components (see also Table 6). In Fig. 24, one can also observe a shift in the value of | q |. Actually, the mean value of | q | in the MINERvA (T2K) 1p1h component moves from 0.69 GeV (0.61 GeV) in the CC0π data-sample to 0.77 GeV (0.73 GeV) in the CC0π1p one. This is a consequence of requesting a proton above 0.45 GeV in the detector.
Finally in Figs. 25 and 26, we show for the MINERvA CC0π sample, the q 0 and | q | distributions as a function of p ⊥ for the same p binning as in Fig. 6. A high occupancy region corresponding to the CC1p1h contribution is clearly observed. A second one, mainly due to CC2p2h and CCRes events, is also visible in q 0 , see Fig. 25 (note the logarithm scale in the z−coordinate). This second enhanced region is less visible in | q |, see Fig. 26, except for the low longitudinal momentum bins. The distributions also show that | q | values are similar for a given p ⊥ independently of the longitudinal momentum except for the case with p < 3 GeV. The main differences between data and MC observed in Fig. 6 as function of p are most probably caused by the CCRes and CC2p2h large q 0 contributions.