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

B.Bourguille b J.Nieves c F.Sánchez111Corresponding author.
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>Institut de Fisica d’Altes Energies (IFAE), The Barcelona Institute of Science
and Technology. Edifici Cn, Universitat Autonoma de Barcelona, Bellaterra
(Barcelona), SpainInstituto de Fisica Corpuscular (IFIC), Centro Mixto CSIC-
Universidad de Valencia, Institutos de Investigacion de Paterna, Apartado
22085, E-46071 Valencia, SpainUniversité de Genève - Faculté des Sciences,
Département de Physique Nucléaire et Corpusculaire (DPNC) 24, Quai Ernest-
Ansermet, CH-1211 Genève 4, Switzerland

###### Abstract

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 kinematic-imbalance 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.

##  1 Introduction

The studies of neutrino-nucleus interactions are entering a new stage,
motivated by long-baseline 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 Alvarez-Ruso:2017oui .
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 Abe:2014oxa , 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 Abe:2012av , extending beyond 10 GeV, 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 Nieves:2012yz ; Alvarez-Ruso:2014bla . It is therefore crucial
to properly include the 2p2h channel into the Monte Carlo (MC) event
generators.

Figure 1:  Left: LFG neutron momentum distribution in carbon. Right: Two-
dimensional neutron momentum distribution correlated to the radial interaction
point inside the nucleus. Distributions are folded with the T2K $\nu-$energy
flux Abe:2012av .

##  2  CCQE $\rm 1p1h$ model

We developed a full exclusive CCQE 1p1h MC event generator based on the
theoretical scheme developed in Nieves:2004wx . The model is capable to
simulate both neutrinos ($\nu_{\mu}n\to\mu^{-}p$) and anti-neutrino222In
general, we will refer to neutrino-induced reactions, unless it is explicitly
mentioned that the discussion is about processes with anti-neutrinos.
($\bar{\nu}_{\mu}p\to\mu^{+}n$) reactions for a variety of nuclei: C,O, Al,
Ti, Fe and Ca. The original code Nieves:2004wx ; Nieves:2011pp  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
Hayato:2009zz . The implemented modifications in NEUT keep all the physics of
the model: Local Fermi Gas (LFG) nucleon-dynamics, 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.

###  2.1  Implementation of the exclusive CCQE model in the MC

We implement the QE model of Ref. Nieves:2004wx  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 and modulus of the target
(hit) nucleon three-momentum from the LFG distribution.

  * •

(anti-)neutrino energy.

  * •

outgoing lepton momentum and angle.

  * •

angle between the outgoing proton (neutron) and the transfer momentum
direction (we test for Pauli blocking at the given radial position).

  * •

The azimuthal angle of the final hadron with respect to the lepton reaction
plane. This angle is generated with a flat probability in the $[-\pi,\pi]$
interval

With this information we can obtain the whole event kinematics applying
conservation of momentum and energy:

  * •

lepton four-momentum $p^{\mu}_{\ell}$.

  * •

target (hit) nucleon four-momentum $p^{\mu}_{N}$.

  * •

final state nucleon four-momentum $p^{\mu}_{N^{\prime}}$.

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 ($\epsilon$). 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.

###  2.2 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 re-interactions. 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_{Fn}=(3\pi^{2}\rho_{n}(r))^{1/3}$ (1)

with $\rho_{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.

Figure 3:  Left: Number of events as function of the $E_{m}^{IA}$ energy for
neutrino scattering off 12C within the LFG model. Both the 1p1h (Eq. (8)) and
2p2h (Eq. (11)) $E_{m}^{IA}$ 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. Center:
Probability to find a neutron in carbon with a momentum ($p_{n}$) for a given
reaction missing energy ($E_{m}^{IA}$) (see Eq. (8)) as predicted by the SF
model Benhar:1989aw ; Benhar:1994hw ; Benhar:2005dj  (contour plot) and for
this implementation of the LFG (box plot). Right: LFG predictions
corresponding to the box plot displayed in the middle panel. In all cases, the
T2K flux Abe:2012av  is used.

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

$\vec{p}_{\nu}+\vec{p}_{N}+\vec{p}_{A-1}=\vec{p}_{\mu}+\vec{p}_{N^{\prime}}+\vec{p}^{\,\prime}_{A-1}$
(2)

where $\vec{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 ($\vec{p}_{N}$) so the total momentum of the initial nucleus
vanishes. Within this approximation, we also consider that the momentum of the
final state nucleus ($\vec{p}^{\,\prime}_{A-1}$) is equal to the residual
momentum of the initial nucleus and both cancel out
($\vec{p}^{\,\prime}_{A-1}=\vec{p}_{A-1}$). Thus the balance of Eq. (2)
reduces to:

$\vec{p}_{\nu}+\vec{p}_{N}=\vec{p}_{\mu}+\vec{p}_{N^{\prime}}$ (3)

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

$E_{\nu}+M_{A}=E_{\mu}+E^{\infty}_{N^{\prime}}+M^{\prime}_{A-1}+\epsilon^{\prime}_{A-1}+T^{\prime}_{A-1}$
(4)

where $M_{A}$ and $M^{\prime}_{A-1}$ are the the ground state masses of the
initial and final nuclei, and $E^{\infty}_{N^{\prime}}$ is the energy of the
asymptotically observed nucleon ($N^{\prime}$). In addition,
$\epsilon^{\prime}_{A-1}>0$ is the excitation energy of the final nucleus,
which average is expected to be between 10 and 20 MeV (see Fig. 2). Finally,
$T^{\prime}_{A-1}$ is the final nucleus kinetic energy, which is very small
(typically of the order of $p_{F_{n}}^{2}/(2M_{A-1})\sim 2$ MeV for carbon
target) and it is approximated to zero in what follows.

The energy balance in Eq. (4) does not apply to cases where any secondary re-
scattering collision changes the energy of the nucleon that emerges after the
weak absorption of the gauge boson333In fact, in these latter situations, the
rupture of the daughter nucleus might occur and the analysis is more
complicated..

The excitation energy can be estimated, in a first approximation, to be the
energy of the hole, within the FG model444Any LFG model implicitly assumes the
existence of a mean-field potential $U=-T_{F}$, which cancels in the
difference of energies, and binds the nucleons.

Figure 4:  Left: Primary proton, created after the absorption of the gauge
boson $W^{+}$ in neutrino processes, momentum distribution in carbon predicted
by the model. Right: Proton momentum distribution predicted by the model as
function of the radial interaction point inside the nucleus. As in Fig. 1,
results have been folded with the T2K neutrino energy flux.

$\epsilon^{\prime}_{A-1}\sim(T_{F}-T_{N})$ (5)

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:

$S_{N}=M^{\prime}_{A-1}-M_{A}+m_{N}$ (6)

where $M_{A}$ and $M^{\prime}_{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 nucleon.
Re-writing Eq. (4) using $S_{N}$, we obtain:

$E_{\nu}+\left(m_{N}+T_{N}-T_{F}\right)=E_{\mu}+S_{N}+E^{\infty}_{N^{\prime}}$
(7)

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 $(\sim 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.

Figure 5:  T2K CC (left) and CC0$\pi$ (right) inclusive double differential
cross-section $d\sigma/dp_{\mu}d\cos{\theta_{\mu}}$. Data, taken from Refs.
Abe:2018uhf  and Abe:2020uub , respectively, are compared to the results
obtained from the present implementation of NEUT with the numerical values for
the $\chi^{2}-$likelihood test compiled in Table 3. Last bin accumulates all
the statistics until the kinematic limit of 30 GeV.

This discussion of the IA energy balance fixes a problem with the nuclear
missing energy ($E_{m}^{IA}$), which appears within the traditional
implementation in MC event generators of the relativistic LFG and RGFG models

$E_{m}^{IA}=E_{\nu}-E_{\mu}-T^{\infty}_{N^{\prime}}=S_{N}+(m_{N^{\prime}}-m_{N})+T_{F}-T_{N}$
(8)

with $E^{\infty}_{N^{\prime}}=T^{\infty}_{N^{\prime}}+m_{N^{\prime}}$. The
value of $E_{m}^{IA}$ becomes negative (non-physical) for some values of
$T_{N}$ when, as it is common, the $T_{F}$ correction is not added555This is
taken into account in some models such as NuWro by adding a constant that
restores the validity of the model. . Equivalently, the problem is caused by
the wrong assumption of taking $E_{N}=\left(m_{N}+T_{N}\right)$, instead of
the correct expression $E_{N}=\left(m_{N}+T_{N}-T_{F}\right)$, which includes
the mean field potential responsible for binding the nuclear system. The
distribution $E_{m}^{IA}$ of energies for a relativistic LFG is depicted in
the left plot of Fig. 3 for neutrino scattering off carbon, where
$S_{n}+(m_{p}-m_{n})=\Delta(^{11}{\rm C})-\Delta(^{12}{\rm C})+\Delta(^{1}{\rm
H})=17.4$ MeV [$\Delta(A_{Z})$ is the mass excess of the $A_{Z}$ nucleus]. The
average $E^{IA}_{m}$ is 28 MeV, very similar to the binding energies used in
RGFG models (25 MeV)Abe:2019arf  or in MINERvA (27.13 MeV)Lu:2018stk . In
addition, we can use this average of $E^{IA}_{m}$ for the LFG model in carbon,
to estimate the average of the excitation energy
$\langle\epsilon^{\prime}_{A-1}\rangle\sim\langle T_{F}-T_{N}\rangle\sim 11$
MeV for this target, using $S_{n}+(m_{p}-m_{n})=17.4$ MeV.

Figure 6:  MINERvA CC (top) and CC0$\pi$ (bottom) inclusive double
differential cross-section $d^{2}\sigma/dp_{\perp}dp_{\parallel}$. Data, taken
from Refs. Filkins:2020xol  and Ruterbories:2018gub , respectively, are
compared to the results obtained from the present implementation of NEUT, with
the $\chi^{2}-$likelihood tests compiled in Table 4.

We pay now attention to the two-dimensional $(p_{N},E_{m}^{IA})$ distribution
shown in the middle and right plots of Fig. 3 following the same
representation as in the Spectral Function (SF) scheme Benhar:1989aw ;
Benhar:1994hw ; Benhar:2005dj . The carbon SF obtained in Benhar:1994hw  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 mass666The
1p1h contribution to the nuclear response function depends on the energy
difference between particle- and hole- nucleons, where the mean-field
potential $-T_{F}(r)$ cancels out., 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. Benhar:1989aw ;
Benhar:1994hw ; Benhar:2005dj . Nevertheless, some differences between both
sets of predictions are visible in Fig. 3, in particular at the edges of the
$(p_{N},E_{m}^{IA})$ contours. This different dependencies might introduce
distinctive differences in 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.

###  2.3  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 Hayato:2009zz .

##  3 Consistent implementation of the CC 2p2h model and the secondary nuclear
collisions

In NEUT Hayato:2009zz , 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.
Nieves:2011pp . In addition to the 1p1h term, the hadron tensor $W^{\mu\nu}$
accounts also for 2p2h contributions evaluated following the LFG scheme of
Ref. Nieves:2011pp , which is fully consistent with the 1p1h implementation
outlined above and based on Nieves:2004wx . 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 Nieves:2011pp . At this
respect, note that in the recent re-computation of Ref. Sobczyk:2020dkn  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 Sobczyk:2020dkn . The produced pair of nucleons is fed into the NEUT
cascade model accounting for the transport of nucleons in the high density
nuclear medium.

Figure 7:  Angular ($\delta\alpha_{T}$,$\delta\phi_{t}$) and transverse
momentum variables ($\vec{\delta}p_{T}$, $\delta p_{T_{X}}$, $\delta
p_{T_{Y}}$) in a reference system, where the neutrino is entering
perpendicular to the plane.

Figure 8:  CC$0\pi 1p$ MINERvA Lu:2018stk ; Cai:2019hpx ; Harewood:2019rzy
(left) and T2K Abe:2018pwo  (right) $\delta\alpha_{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 $\chi^{2}-$values are compiled in Tables
3 and 4.

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:

$\displaystyle E_{\nu}$ $\displaystyle=$ $\displaystyle
E_{\mu}+E^{\infty}_{N^{\prime}_{1}}+E^{\infty}_{N^{\prime}_{2}}+S_{N_{1}N_{2}}$
(9) $\displaystyle-$
$\displaystyle\left(m_{N_{1}}+T_{N_{1}}-T_{F_{1}}\right)-\left(m_{N_{2}}+T_{N_{2}}-T_{F_{2}}\right)$

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

$S_{N_{1}N_{2}}=M^{\prime}_{A-2}-M_{A}+m_{N_{1}}+m_{N_{2}}$ (10)

Finally, $N^{\prime}_{1}$ and $N^{\prime}_{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

$E_{m}^{IA}\Big{|}_{\rm
2p2h}=E_{\nu}-E_{\mu}-T^{\infty}_{N^{\prime}_{1}}-T^{\infty}_{N^{\prime}_{2}}\,,$
(11)

and compared to the 1p1h $E_{m}^{IA}$ values discussed in the previous
section. The excitation of two nucleons leads to a bigger offset of the
$E_{m}^{IA}$ 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.

The remaining neutrino-nucleus interaction channels, mostly with associated
pion production, are simulated based on the existing NEUT Hayato:2009zz  Monte
Carlo event generator. The resonant pion production is based on the Rein-
Sehgal model Rein:1980wg , taking into account eighteen resonances with masses
below 2 GeV and their interference terms, with the axial mass fixed to
$M_{\text{A}}^{\text{RES}}=1.21$ GeV. This model has been compared to
experimental data by the T2K collaboration showing remarkable
agreementAbe:2019arf . Neutral current and charged current coherent pion
production is simulated using the Rein-Sehgal model in Ref. Rein:1982pf . The
CC coherent pion production includes PCAC (partially conserved axial-vector
current) and lepton mass corrections, as discussed in Rein:2006di . DIS
processes are simulated using the GRV98 Gluck:1998xa  parton distribution,
with low-Q2 corrections from the Bodek and Yang model Bodek:2003wc . Secondary
interactions of pions inside the nucleus are simulated using an intra-nuclear
cascade model based on the method developed by Salcedo et al. Salcedo:1987md ,
tuned to external $\pi-$12C data PinzonGuerra:2018rju .

##  4  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 TKI variables to explore the limits of the hadron kinematic
predictions of the model.

###  4.1 Event Simulation and data selection

Events are simulated using the NEUT package with the CCQE and CC2p2h reaction-
mechanisms 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 re-scattering 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.

####  4.1.1 T2K data samples

The neutrino T2K data-sample has different selections obtained from the off-
axis muon neutrino beam, which peaks around 0.6 GeV but it contains a large
energy-tail ranging to the region of tens of GeV. The CC inclusive measurement
considers only the muon production kinematics ignoring all hadronic activity
Abe:2018uhf . The T2K selection criterion for CC0$\pi$ Abe:2020uub  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 TKI variables analysis of the
CC0$\pi 1p$ sample requires the detection of a muon and proton with the
following conditions Abe:2018pwo :

  * •

$0.45~{}{\rm GeV}<|\vec{p}_{p}|<1.0{\rm~{}GeV}$.

  * •

$\cos{\theta_{p}}>0.4$

  * •

$0.25{\rm~{}GeV}<|\vec{p}_{\mu}|<10{\rm~{}GeV}$.

  * •

$\cos{\theta_{\mu}}>-0.6$

with $\theta_{p}$ the polar angle of the outgoing proton, which has a modulus
of the three-momentum $|\vec{p}_{p}|$. T2K presents also a slightly less
restrictive selection criterion for CC0$\pi 1p$ requiring the detection of a
muon and proton with the following conditions Abe:2018pwo :

  * •

$0.45~{}{\rm GeV}<|\vec{p}_{p}|<1.0{\rm~{}GeV}$.

  * •

$\cos{\theta_{p}}>0.4$

Figure 9:  CC$0\pi 1p$ MINERvA Lu:2018stk ; Cai:2019hpx ; Harewood:2019rzy
(left) and T2K Abe:2018pwo  (right) distributions for the TKI angular variable
$\delta\phi_{T}$. Details of the comparison with NEUT results as in Fig. 8.

####  4.1.2 MINERvA data samples

MINERvA, with an average neutrino energy of 3.5 GeV, published different event
selections: the CC inclusive, the CC0$\pi$ that removes events with detected
charged $\pi$’s and electromagnetic activity to eliminate events with
$\pi^{0}$, and the CC$0\pi 1p$ sample which is a sub-sample of CC0$\pi$
requesting the presence of an identified proton in the final state. Another
data sample from MINERvA, the so-called available energy Rodrigues:2015hik ,
requires full simulation of the detector simulation that falls beyond the
capabilities of this work.

The CC inclusive selection in MINERvA requires Filkins:2020xol  a muon with
polar angle ($\theta_{\mu}$) below $20^{o}$. The selection criteria for
CC$0\pi$ demands in addition Ruterbories:2018gub :

  * •

a muon with momentum $|\vec{p}_{\mu}|$ between 1.5 GeV and 10 GeV,

  * •

neither charged, nor neutral pions escaping the nucleus.

  * •

no $\gamma$ with energies about 10 MeV. We check this cut actually do not
affect the MC prediction for CCQE and CC2p2h.

The selection criteria for CC$0\pi 1p$ requires in addition Lu:2018stk ;
Cai:2019hpx ; Harewood:2019rzy :

  * •

a proton with polar angle ($\theta_{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.

###  4.2  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$\pi$ cross-
section that is reported by T2K for a pure carbon target Abe:2020uub  and not
in hydrocarbon as for the CC inclusive and MINERvA.

Component | CH | O | Al | Si | Cl | Ti
---|---|---|---|---|---|---
Weight (%) | 95.02 | 3.18 | 0.26 | 0.27 | 0.55 | 0.69
Nuclei (%) | 96.71 | 2.63 | 0.13 | 0.13 | 0.21 | 0.19
Nucleons (%) | 95.04 | 3.18 | 0.26 | 0.27 | 0.55 | 0.69
Table 1:  MINERvA target nuclear composition in weight fraction and the translation to the fractional composition in nuclei. Last row shows the fraction of nucleons of each of the spices with respect to the total. Data taken from Ref. Aliaga:2013uqz  Component | CH | O | Si | Ti
---|---|---|---|---
Weight (%) | 95.02 | 3.18 | 0.27 | 0.69
Nuclei (%) | 97.04 | 2.64 | 0.12 | 0.19
Nucleons (%) | 95.82 | 3.21 | 0.27 | 0.70
Table 2:  The same as Table 1, but for the T2K experiment. Data taken from
Ref. Amaudruz:2012esa

Figure 10:  CC$0\pi 1p$ MINERvA Lu:2018stk ; Cai:2019hpx ; Harewood:2019rzy
(left) and T2K Abe:2018pwo  (right) distributions for the missing transverse
momentum $|\delta\vec{p}_{T}|$. Details of the comparison with NEUT results as
in Fig. 8.

Figure 11:  MINERvA $\delta p_{T_{x}}$ (top) and $\delta p_{T_{y}}$ (bottom)
distributions Lu:2018stk ; Cai:2019hpx ; Harewood:2019rzy . Left and right
panels show the comparison (details as in Fig. 8) with NEUT in linear and
logarithmic scales, respectively.

###  4.3  Inclusive cross-sections

In this section we compare our model predictions with the CC inclusive and
CC0$\pi$ 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 data-to-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$\pi$ 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 both, the pion re-interactions inside the nucleus and
the primary pion production mechanism.

T2K published both the CC inclusive Abe:2018uhf  and the CC0$\pi$ Abe:2020uub
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$\pi$ 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 Filkins:2020xol
(see top panels of Fig. 6). The cross-section is dominated by DIS, but there
are regions ($1.5~{}\rm{GeV}~{}<p_{\parallel}<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$\pi$ sample from MINERvA is also shown in Fig. 6 (bottom 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_{\parallel}>5.0$ GeV), but also for low ($p_{\parallel}<3.5$ GeV)
longitudinal momenta.

To estimate the agreement data-model, we computed the $\chi^{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 $\chi^{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 $\chi^{2}-$value found in this work is almost a factor of
two smaller than those presented in the experimental paper (NEUT 5.3.2 and
GENIE Andreopoulos:2009rq  2.8.0). The theoretical description of the T2K
CC$0\pi$ data sample is better than that achieved for the inclusive one, with
$\chi^{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 CC$0\pi$ event
selection. The $\chi^{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 Abe:2020uub , since both largely share the physics implementation.

On the other hand, the best-fit $\chi^{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. The description of the MINERvA CC$0\pi$ events,
where DIS has been practically removed, significantly improves, and it is
quantitatively slightly poorer than that seen above for T2K. Nevertheless, the
$\chi^{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 CC$0\pi$ show apparently a larger, though still soft,
dependence on the proton re-scattering probability.

variable | dof | nominal | NrSP: 50% | NrSP: 150%
---|---|---|---|---
CC inclusive | 71 | 110 | 104 | 104
CC$0\pi$ inclusive | 29 | 17 | 17 | 17
$\delta\alpha_{T}$ | 8 | 35 | 29 | 43
$\delta\phi_{T}$ | 8 | 19 | 27 | 36
$\delta p_{T}$ | 8 | 7 | 32 | 13
$\big{|}\Delta\vec{p}\,\big{|}$ | 49 | 261 | 193 | 398
$\Delta\theta$ | 35 | 1021 | 798 | 1322
$\Delta|\vec{p}\,|$ | 49 | 130 | 111 | 177
Table 3:  $\chi^{2}-$likelihood test for T2K variables. The number of degrees of freedom (dof) is given in the second column, while in the another three ones, results obtained for different distributions and three NrSP –proton re-scattering probability– configurations (nominal, 50% and 150% strengths) are collected. variable | dof | nominal | NrSP: 50% | NrSP: 150%
---|---|---|---|---
CC inclusive | 144 | 420 | 435 | 420
CC0$\pi$ | 144 | 208 | 247 | 213
$\delta\alpha_{T}$ | 12 | 21 | 25 | 18
$\delta\phi_{T}$ | 23 | 47 | 105 | 46
$\delta p_{T}$ | 24 | 89 | 155 | 56
$\delta p_{T_{x}}$ | 32 | 63 | 106 | 43
$\delta p_{T_{y}}$ | 33 | 56 | 108 | 40
$|\vec{p}_{\mu}|$ | 32 | 29 | 41 | 26
$\theta_{\mu}$ | 19 | 23 | 24 | 22
$p_{p}$ | 25 | 30 | 33 | 31
$\theta_{p}$ | 26 | 49 | 62 | 43
$p_{n}$ | 24 | 107 | 202 | 86
Table 4:  Same as in Table 3, but for MINERvA experiment.

###  4.4  Transverse kinematic-imbalance (TKI) variables

Figure 12:  Inclusive $|\vec{p}_{\mu}|$ (left) and $\theta_{\mu}$ (right)
distributions from the CC$0\pi 1p$ MINERvA sample Lu:2018stk ; Cai:2019hpx ;
Harewood:2019rzy . Details of the comparison with NEUT results as in Fig. 8.

Figure 13:  MINERVA CC0$\pi$1p differential cross sections Lu:2018stk ;
Cai:2019hpx ; Harewood:2019rzy  in proton kinematics, momentum
($|\vec{p}_{p}|$) and angle ($\theta_{p}$), together with results from the
current implementation of NEUT. Details of the comparison with NEUT as in Fig.
8. Figure 14:  MINERvA CC0$\pi$1p differential cross sections in reconstructed
$|\vec{p}_{n}|$. Data from Ref. Lu:2018stk ; Cai:2019hpx ; Harewood:2019rzy .
We also show results from the current implementation of NEUT, with details of
the comparison as in Fig. 8.

TKI variables Lu:2015hea ; Lu:2015tcr , 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
to the neutrino direction, momentum is computed as

$|\delta\vec{p}_{T}|\equiv|\vec{p}_{T}|=\frac{\left|(\vec{p}_{\mu}+\vec{p}_{p})\times\vec{p}_{\nu}\right|}{|\vec{p}_{\nu}|}$
(12)

In addition, the missing transverse momentum is separated in two components in
relation to the $\nu-\mu$ reaction plane, defined by the neutrino and the
emitted charged lepton. One of the transverse components is contained in the
reaction plane $\delta p_{T_{y}}$, while the another one, $\delta p_{T_{x}}$,
is perpendicular to this plane Cai:2019hpx .

On the other hand, the TKI angular variables read:

$\displaystyle\delta\phi_{T}$ $\displaystyle=$
$\displaystyle\arccos\left[\frac{-\vec{p}_{T_{p}}\cdot\vec{p}_{T_{\mu}}}{|\vec{p}_{T_{p}}||\vec{p}_{T_{\mu}}|}\right]$
(13) $\displaystyle\delta\alpha_{T}$ $\displaystyle=$
$\displaystyle\arccos\left[\frac{-\vec{p}_{T}\cdot\vec{p}_{T_{\mu}}}{|\vec{p}_{T}|\cdot|\vec{p}_{T_{\mu}}|}\right]$
(14)

where $\vec{p}_{T}=\vec{p}_{T_{\mu}}+\vec{p}_{T_{p}}$, with
$\vec{p}_{T_{\mu}}$ (or $\vec{p}_{\perp}$, as used in Fig. 6) and
$\vec{p}_{T_{p}}$ the transverse projections of the muon and proton momenta to
the neutrino direction. This discussion focuses on the QE-like process
$\nu_{\mu}+A\to\mu^{-}+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,
$\delta\vec{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
$|\delta\vec{p}\,|$ gives the magnitude of its recoil momentum, which can be
obtained Furmanski:2016wqo  independently of the unknown incident neutrino
energy. Moreover, assuming perfect balance of momentum in the reaction (see
discussion of Eqs. (2) and (3)), $|\delta\vec{p}\,|$ can be identified to the
neutron target momentum, which is then given in terms of measurable quantities
Lu:2018stk ; Furmanski:2016wqo

$|\vec{p}_{n}|=\sqrt{|\vec{p}_{T}|^{2}+|\vec{p}_{L}|^{2}}$ (15)

with

$|\vec{p}_{L}|=\frac{(M_{A}-E_{\mu}-E_{p}^{\infty}+|\vec{p}_{L_{p}}|)^{2}-(M^{\prime}_{A-1})^{2}-|\vec{p}_{T}|^{2}}{2(M_{A}-E_{\mu}-E_{p}^{\infty}+|\vec{p}_{L_{p}}|+|\vec{p}_{L_{\mu}}|)}$
(16)

where $\vec{p}_{L_{p}}$ and $\vec{p}_{L_{\mu}}$ (or $\vec{p}_{\parallel}$, 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.
Lu:2018stk ; Cai:2019hpx ; Harewood:2019rzy , $M_{A}$ is the mass of the
carbon target, while for $M^{\prime}_{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 neutron777Note, this is the prescription used in the
experimental work, which should be used to compare to the event-distributions
provided in that paper, but it does not correspond to the energy-balance of
Eq. (7) proposed here for genuine QE processes..

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. Nieves:2012yz ):

$E_{\nu}^{\rm
rec}=\frac{m_{p}^{2}-m_{\mu}^{2}+2E_{\mu}(m_{n}-E_{b})-(m_{n}-E_{b})^{2}}{2((m_{n}-E_{b})-E_{\mu}+|\vec{p}_{\mu}|\cos{\theta_{\mu}})}$
(17)

where $\theta_{\mu}$ is the angle of the muon with respect to the average
neutrino direction, $m_{\mu}$ is the mass of the muon, $m_{p}$ is the mass of
the outgoing proton and the constant $E_{b}$ is fixed to 25 MeV in the
experimental T2K results Abe:2018pwo . Once the neutrino energy is known, one
can compute the so-called inferred proton momentum:

$\vec{p}_{p}^{\,\,\rm inf}=\vec{p}^{\,\,\rm rec}_{\nu}-\vec{p}_{\mu}$ (18)

New observables are build by comparing the inferred momentum with the
experimentally measured proton momentum ($\vec{p}_{p}^{\,\infty}$):

$\displaystyle|\Delta\vec{p}\,|$ $\displaystyle=$
$\displaystyle|\vec{p}_{p}^{\,\,\rm inf}-\vec{p}_{p}^{\,\infty}|$ (19)
$\displaystyle\Delta|\vec{p}\,|$ $\displaystyle=$
$\displaystyle|\vec{p}_{p}^{\,\,\rm inf}|-|\vec{p}_{p}^{\,\infty}|$ (20)
$\displaystyle\Delta\theta$ $\displaystyle=$
$\displaystyle\arccos\left[\frac{\vec{p}_{p}^{\,\infty}\cdot\vec{p}_{\nu}}{|\vec{p}_{p}^{\,\infty}||\vec{p}_{\nu}|}\right]-\arccos\left[\frac{\vec{p}_{p}^{\,\,\rm
inf}\cdot\vec{p}_{\nu}}{|\vec{p}_{p}^{\,\,\rm inf}||\vec{p}_{\nu}|}\right]$
(21)

The MINERvA Collaboration reported different cross-sections from its CC$0\pi
1p$ data sample Lu:2018stk ; Cai:2019hpx ; Harewood:2019rzy  :

  * •

Inclusive lepton momentum and angle.

  * •

Reconstructed kinematics such as the target nucleon momentum $|\vec{p}_{n}|$.

  * •

Visible proton momentum.

  * •

TKI angles: $\delta\phi_{T}$ and $\delta\alpha_{T}$.

  * •

Transverse momentum balance: $\delta p_{T}$, $\delta p_{T_{x}}$, $\delta
p_{T_{y}}$.

while T2K reported on slightly different set of variables from its CC$0\pi 1p$
sample Abe:2018pwo :

  * •

Transverse angles (TKI variables with first selection criterion):
$\delta\phi_{T}$ and $\delta\alpha_{T}$, and the missing transverse momentum
$|\delta\vec{p}_{T}|$

  * •

Proton momentum balance (Inferred variables with the second selection
criterion): $|\Delta\vec{p}\,|$, $\Delta\theta$ and $\Delta|\vec{p}\,|$ in
bins of muon momentum and angle.

Experimental and NEUT predictions for the TKI angular variables
$\delta\alpha_{T}$ and $\delta\phi_{T}$, and the missing transverse momentum
$|\delta\vec{p}_{T}|$, both for MINERvA (left) and T2K (right) are shown in
Figs. 8, 9 and 10, respectively. Absolute $\chi^{2}-$values obtained from the
comparison with data are compiled in Tables 3 and 4, respectively. The NEUT
results for the TKI 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 $\chi^{2}-$likelihood
tests. We should point out that NrSP effects change not only the overall
normalization, but also the shape of the distributions. This is clearly
visible, for instance, in the CC$0\pi 1p$ T2K Abe:2018pwo  $\delta\alpha_{T}$
distributions depicted in Fig. 8. There, we see that re-scattering effects
become more important in the region of smallest $\delta\alpha_{T}$ angles.
This observation is in order for some other distributions discussed below in
this subsection.

The transverse momentum components contain different information. The
component $\delta p_{T_{x}}$ is expected to be symmetric around zero, with a
width which depends on the target neutron momentum and further re-scattering
effects, while $\delta p_{T_{y}}$ 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 $\delta p_{T_{x}}$ ($\delta p_{T_{y}}$) distribution.
The $\chi^{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 TKI observables.

MINERvA collaboration also reported $|\vec{p}_{\mu}|$ and $\theta_{\mu}$
distributions from its CC$0\pi 1p$ MINERvA sample Lu:2018stk ; Cai:2019hpx ;
Harewood:2019rzy . 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 $\chi^{2}/dof$ around one (see Table 4)
for nominal NrSP, and some dependence on this latter input as expected when
analysing CC$0\pi 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 additional invisible
protons do not suffer from large momentum transfer biases. The difference
between the two are dominated by low $|\vec{q}\,|$ contributions with feed-
down background cause by NrSP. These tendencies can be observed in Figs. 25
and 26 of Appendix A, where the model estimation for the energy
($q^{0}=E_{\nu}-E_{\mu}$) and momentum
($|\vec{q}\,|=|\vec{p}_{\nu}-\vec{p}_{\mu}|$) transfer distributions for the
MINERvA and T2K CC inclusive, CC$0\pi$ and CC$0\pi$1p event selections are
shown.

The overall agreement with NEUT is good, showing the importance of 2p2h
mechanisms. The results for the TKI variables are good for the MINERvA data,
with statistically acceptable values of $\chi^{2}/dof$ for most of the cases.
The worst comparison is obtained for the reconstructed $|\vec{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
($|\vec{p}_{n}|\geq 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.

Figure 15:  T2K CC0$\pi$1p $|\Delta\vec{p}\,|$ distribution. The panels
correspond to different muon kinematic bins. From left to right and up to
down: $-1<\cos\theta_{\mu}<-0.6$, $-0.6<\cos\theta_{\mu}<0$ with
$|\vec{p}_{\mu}|<250$ MeV, $-0.6<\cos\theta_{\mu}<0$ with
$|\vec{p}_{\mu}|>250$ MeV, $0<\cos\theta_{\mu}<1$ with $|\vec{p}_{\mu}|<250$
MeV, $0<\cos\theta_{\mu}<0.8$ with $|\vec{p}_{\mu}|>250$ MeV,
$0.8<\cos\theta_{\mu}<1$ with $250~{}{\rm MeV}<|\vec{p}_{\mu}|<750$ MeV and
$0.8<\cos\theta_{\mu}<1$ with $|\vec{p}_{\mu}|>750$ MeV. Data taken from Ref.
Abe:2018pwo .  Figure 16:  T2K CC0$\pi$1p $|\Delta\theta|$ distribution. The
panels correspond to the muon kinematic bins specified in Fig. 15. Data taken
from Abe:2018pwo .  Figure 17:  T2K CC0$\pi$1p $\Delta|\vec{p}\,|$
distribution. The panels correspond to the muon kinematic bins specified in
Fig. 15. Data taken from Abe:2018pwo .

The agreement with T2K CC0$\pi$1p data is less impressive and the obtained
$\chi^{2}/dof$, see Table 3, are large for all the three observables
$|\Delta\vec{p}\,|$, $|\Delta\theta|$ and $\Delta|p|$ reported in Abe:2018pwo
, and shown here in Figs. 15–17. The worst situation is found for the
$\Delta\theta$ distribution, with the largest contributions to $\chi^{2}$
produced by the negative $\Delta\theta$ 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 $\chi^{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 $\Delta\theta$). The $\chi^{2}$ figures
for $|\Delta\vec{p}\,|$, and $\Delta|p|$ are also reduced with the re-scaling
factor NrSP=0.5.

In summary from the results of Tables 3 and 4, we conclude that MINERvA TKI
distributions are better described with the 150% enhanced NrSP, while the T2K
ones favor either nominal ($\delta\phi_{T}$ and $\delta p_{T}$) or the 50%
reduced NrSP ($\delta\alpha_{T}$, $\big{|}\Delta\vec{p}\,\big{|}$,
$\Delta\theta$ and $\Delta|\vec{p}\,|$) configurations.

The difference between the agreement found for MINERvA and T2K data samples
might point to an energy dependent deviation. The LFG model could provide a
better approximation to the MINERvA energies than to the T2K ones, which would
be more sensitive to finer details of the low-lying nuclear levels. However,
an overall $\chi^{2}-$analysis is not sufficient to extract robust conclusions
on the energy dependence of the LFG model, and it should also be noted that,
as discussed above, MINERvA and T2K data-sets favor different proton re-
scattering configurations. This could be an indication that the different
agreement exhibited for the MINERvA and T2K data-samples might be also
partially produced by higher energy channels such as 2p2h or resonance
mechanisms. At this respect, it would be very useful to have access to the T2K
$d\sigma/dp_{n}$ cross-section, since this distribution is specially sensitive
to non CC1p1h contributions, as seen in Fig. 14 for MINERvA data.

###  4.5  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$\pi$ samples, while its predictions are
marginally larger for the CC0$\pi$1p data-set. The deviations with 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$\pi$ and the CC0$\pi$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$\pi$1p cross-section is
shown in Table 7, and it is also included in the proportions collected in
Table 6. 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$\pi$ 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$\pi$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 ($\leq$ 450
MeV). A larger re-scattering probability would approach the two results.
Another possible cause of the discrepancy is the larger $|\vec{q}\,|$ values
intrinsic to CC0$\pi$1p with respect to CC0$\pi$ 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.

| Exp [$10^{-39}$ cm2] | Model [$10^{-39}$ cm2] |  $\frac{Exp-Model}{Model}$ [$\%$]
---|---|---|---
MINERvA
CC inclusive | $18.3\pm 1.3$ | $16.7$ | $9.8\pm 7.8$
CC0$\pi$ | $4.64\pm 0.38$ | $3.96$ | $17.1\pm 9.7$
CC0$\pi$1p from $\delta\alpha_{T}$ | $1.85\pm 0.17$ | $1.94$ | $-4.6\pm 9.0$
T2K
CC inclusive | $7.33\pm 0.81$ | $6.01$ | $22\pm 13$
CC0$\pi$ | $4.74\pm 0.59$ | $3.78$ | $25\pm 16$
CC0$\pi$1p from $\delta\alpha_{T}$ | $1.39\pm 0.13$ | $1.46$ | $-4.6\pm 9.0$
CC0$\pi$1p from $|\Delta\vec{p}\,|$ | $3.01\pm 0.34$ | $2.88$ | $5\pm 12$
Table 5:  Integrated cross-sections for the different event selections provided by T2K and MINERvA. The CC0$\pi$1p integrated cross-section has been computed using the $d\sigma/d(\delta\alpha_{T})$ differential distribution for MINERvA and T2K and also $d\sigma/d|\Delta\vec{p}\,|$ for T2K. In all cases, the nominal NrSP configuration is used. | 1p1h (%) | 2p2h (%) | Res (%) | DIS-others (%)
---|---|---|---|---
MINERvA
CC inclusive | 19.0 | 4.7 | 35.4 | 40.9
CC0$\pi$ | 64.8 | 15.2 | 17.8 | 2.1
CC0$\pi$1p from $\delta\alpha_{T}$ | $56.8^{-1.4}_{+1.3}$ | $16.9^{+0.2}_{-0.3}$ | $23.5^{+1.1}_{-1.0}$ | $2.7^{+0.2}_{-0.1}$
T2K
CC inclusive | 46.9 | 6.2 | 31.8 | 15.0
CC0$\pi$ | 78.8 | 10.4 | 10.0 | 0.8
CC0$\pi$1p from $\delta\alpha_{T}$ | $80.2^{-1.0}_{+0.8}$ | $8.4^{+0.2}_{-0.1}$ | $10.8^{+0.7}_{-0.6}$ | $0.6^{+0.0}_{-0.0}$
CC0$\pi$1p from $|\Delta\vec{p}\,|$ | $77.0^{-1.2}_{+0.9}$ | $9.8^{+0.2}_{-0.1}$ | $12.7^{+0.9}_{-0.8}$ | $0.7^{+0.1}_{-0.0}$
Table 6:  Different NEUT contributions to the T2K and MINERvA cross sections given in Table 5. As in this latter table, the two T2K CC0$\pi$1p selections have been considered and all results are obtained using the nominal NrSP. For those samples with explicit proton in the final state (CC0$\pi$1p), we compute the variation of the proportions induced by the use of stronger or weaker proton re-scattering. The differences with respect to the results obtained from the nominal NrSP configuration are shown as an upper (NRsP 150%) and a lower (NRsP 50%) error. As expected, the increase (decrease) in the NRsP strength reduces (increases) the population of CC1p1h in the sample, affecting indirectly the proportions of the other modes.  Nominal Model | NRsP 50% | NRsP 150%
---|---|---
[cm2] | [cm2] | [cm2]
MINERvA
$1.94\times 10^{-39}$ | $2.09\times 10^{-39}$ | $1.81\times 10^{-39}$
| $+7.7\%$ | $-6.7\%$
T2K
$1.46\times 10^{-39}$ | $1.60\times 10^{-39}$ | $1.35\times 10^{-39}$
| $+9.6\%$ | $-7.5\%$
Table 7:  T2K and MINERvA CC0$\pi$1p total cross-section obtained using NRsP
(proton re-scattering probability) configurations of 50% and 150%. We use the
$d\sigma/d(\delta\alpha_{T})$ differential distribution to perform the
numerical integrations. In percentage, we also show the relative variation
with respect to the nominal Model.

##  5  Data vs theoretical predictions in terms of the scaling variable

Figure 18:  Left: CC inclusive MINERvA [
$|\vec{q}\,|=|\vec{p}_{\nu}-\vec{p}_{\mu}|$, $q^{0}=(E_{\nu}-E_{\mu})$] 2D
distribution predicted by the NEUT CC inclusive event generator. The black
solid lines mark fix $\psi^{\prime}$ values across the
$(q^{0},|\vec{q}\,|)-$plane . Right: $\psi^{\prime}$ distribution of events
obtained from the 2D one shown in the left-panel, and separated by the primary
neutrino-nucleon interaction modes.  Figure 19:  Two dimensional
($p_{\parallel},p_{\perp}$) lepton distributions accumulated for different
$\psi^{\prime}-$intervals from the NEUT CC inclusive event generator for
MINERvA flux. The contributions (in percentage) of the different interaction
modes (CC1p1h, C2p2h, CCRes and CCOthers) are also shown in each of the
panels.

Figure 20:  Top: Distribution of the scaling variable $\psi^{\prime}$ versus
the neutrino energy for the CC inclusive MINERvA (left) and T2K (right)
experiments, as predicted by NEUT. Bottom: The distributions shown in the
upper panels, but normalised to unity by column (energy-bin).

Figure 21:  Left: Data/MC ratio $R(\psi^{\prime})$ calculated for the MINERvA
CC inclusive cross-section, and split into the CC1p1h (red), CC2p2h (yellow),
CCRes (gray) and CC-DIS (pink) contributions. Right: $R(\psi^{\prime})$
correlation matrix computed for the MINERvA CC inclusive cross-sections.

Figure 22:  The same as Fig. 21, but for the MINERvA CC0$\pi$ cross-section.

Figure 23:  The same as Fig. 21, but for the T2K CC inclusive cross-section.

Figure 24:  The same as Fig. 21, but for the T2K CC0$\pi$ cross-section.

The scaling variable $\psi^{\prime}$ has been shown to be a powerful tool to
learn details on the neutrino-nucleus interaction Gonzalez-Jimenez:2014eqa .
This variable, proposed originally for electron scattering Donnelly:1999sw ;
Barbaro:1998gu , has been adapted to neutrino-nucleus interactions recently.
The scaling variable ($\psi^{\prime}$) is defined as:

$\psi^{\prime}=\frac{1}{\sqrt{\sqrt{1+\eta_{F}^{2}}-1}}\frac{\lambda-\tau}{\sqrt{(1+\lambda)\tau+\kappa\sqrt{\tau(1+\tau)}}}$
(22)

with

$\displaystyle\eta_{F}=\frac{p_{F}}{m_{n}},\,\kappa=\frac{|\vec{q}\,|}{2m_{n}},\,\lambda=\frac{q^{0}-E_{\rm
shift}}{2m_{n}}\,\text{and},\,\tau=\kappa^{2}-\lambda^{2}$ (23)

where, in addition to variables already introduced, $p_{F}$ is the Fermi
momentum and $E_{\rm shift}$ is a energy shift, which are fixed in this
discussion to 228 MeV and 20 MeV, respectively.

The scaling variable $\psi^{\prime}$ allows to project the energy-momentum
transfer $(q^{0},|\vec{q}\,|)$ 2D-sample of events into a 1D distribution, as
illustrated in Fig. 18 with the NEUT predictions for the CC inclusive MINERvA
sample. The full 2D distribution is shown in the left plot, where we also
display the $\psi^{\prime}-$constant curves. We see that a
$(q^{0},|\vec{q}\,|$) pair determines an unique $\psi^{\prime}$ value, while
the same scaling variable can be constructed out of different energy-momentum
transfer combinations. The $\psi^{\prime}-$distribution of events is depicted
in the right panel. Though $\psi^{\prime}$ is based on muon kinematics, the
underlying contributions from different mechanisms (1p1h, 2p2h Res and DIS)
are better separated since they lead to different $\psi^{\prime}-$behaviors,
as can be appreciated in the figure. Hence, this variable might provide a
method to effectively disentangle the components of nuclear effects in data
and therefore extract valuable constraints on the theoretical model
ingredients. This is the objective of the new analysis proposed in this work,
and the details of which are discussed below.

The scaling variable $\psi^{\prime}$ 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 bins888Experiments normally report $(|\vec{p}_{\mu}|,\theta_{\mu})$
or $(p_{\parallel},p_{\perp})$, $d\sigma_{\rm Exp}/d\vec{p}_{\mu}$, which can
be compared with theoretical MC differential distributions for the same
binning,

$R(\vec{p}_{\mu})=\frac{d\sigma_{\rm Exp}/d\vec{p}_{\mu}}{d\sigma_{\rm
MC}/d\vec{p}_{\mu}}$ (24)

On the other hand within the theoretical MC model, one can associate each
event with a value of the scaling variable $\psi^{\prime}$ (in general, values
of the scaling variable comprised in a certain bin, since both neutrino energy
and muon-kinematics are binned). We propose to express the ratio of data to
theoretical MC predictions (data/MC, in what follows) as function of
$\psi^{\prime}$. Thus, we define the average data/MC ratio as:

$\displaystyle R(\psi^{\prime})$ $\displaystyle=$ $\displaystyle\sum_{{\rm
events\,MC}\,\in\,\psi^{\prime}}\,\sum_{p_{\mu}\theta_{\mu}}f(\vec{p}_{\mu}|\psi^{\prime})R(\vec{p}_{\mu})$
(25) $\displaystyle f(\vec{p}_{\mu}|\psi^{\prime})$ $\displaystyle=$
$\displaystyle\frac{N_{\rm
MC}(\vec{p}_{\mu}|\psi^{\prime})}{N(\psi^{\prime})}$ (26)

with $N(\psi^{\prime})$ the number total of MC generated events, and $N_{\rm
MC}(\vec{p}_{\mu}|\psi^{\prime})$ the number of events with muon kinematics
comprised in the bin around $\vec{p}_{\mu}$, which gives rise to the value
$\psi^{\prime}$ for the scaling variable999Note that Eq. (25) admits a matrix
interpretation of the type $B=MA$, where $A$ and $B$ are ratio vectors and $M$
is a matrix relating $\psi^{\prime}$ and the kinematics of the observed muon.
. Thus, $f(\vec{p}_{\mu}|\psi^{\prime})$ is the fraction of events predicted
by the theoretical model for a given value of $\psi^{\prime}$ and muon
kinematics $\vec{p}_{\mu}$, taking into account the neutrino energy spectrum
of the experiment. The distribution $R(\psi^{\prime})$ is built in such a way
that the ratio between data and the theoretical MC results are weighted
according to the population $f(\vec{p}_{\mu}|\psi^{\prime})$. In the limit in
which $f=1$ for one kinematics-bin and there is no more than one value of
$\psi^{\prime}$ contributing to this bin, $R$ is no more than the ratio of
data to MC for this given value of $\psi^{\prime}$. 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_{\rm MC}(\vec{p}_{\mu}|\psi^{\prime})$
can be split into the different mechanisms (1p1h, 2p2h, Res, DIS,…) considered
in the theoretical approach implemented in the MC. In Fig. 19, we illustrate
the physics content of the 3D $f(\vec{p}_{\mu}|\psi^{\prime})$ transfer matrix
for the CC inclusive MINERvA sample, as predicted by NEUT. In the figure, we
show the number of events $N_{\rm MC}(\vec{p}_{\mu}|\psi^{\prime})$,
accumulated for different $\psi^{\prime}$ intervals. As expected from the top
panel of Fig. 6, the largest contributions to $N_{\rm
MC}(\vec{p}_{\mu}|\psi^{\prime})$, for all $\psi^{\prime}-$regions, are
concentrated in the two dimensional region $[1.5\,{\rm GeV}<p_{\parallel}<4$
GeV] $\times$ $[0.25\,{\rm GeV}<p_{\perp}<1$ GeV]. However, the relative
contributions of the different interaction modes (CC1p1h, C2p2h, CCRes and
CCOthers) change significantly with the $\psi^{\prime}-$bin, following a
pattern consistent with the distribution displayed in the right plot of Fig.
18.

We stress the trivial observation that a $(p_{\parallel},p_{\perp})$ pair does
not unequivocally determine a value of $\psi^{\prime}$, as clearly illustrated
in Fig. 19. This is because the neutrino beam is not monochromatic, and for
each neutrino energy one has a different relation between the
$(p_{\parallel},p_{\perp})$ and $(q^{0},|\vec{q}\,|)$ pairs101010One has
$q^{0}=E_{\nu}-\sqrt{m_{\mu}^{2}+p_{\parallel}^{2}+p_{\perp}^{2}}$ and
$\vec{q}^{\,2}=(E_{\nu}-p_{\parallel})^{2}+p_{\perp}^{2}$. As a consequence in
a $(p_{\parallel},p_{\perp})$ distribution, a fix value of $\psi^{\prime}$
will not be represented by a curve, but instead by a 2D region, with large
overlaps between different $\psi^{\prime}-$regions. We should note, related to
this discussion, that the total number of events for a given $\psi^{\prime}$
[$N(\psi^{\prime})$ in the definition of $f(\vec{p}_{\mu}|\psi^{\prime})$ in
Eq. (26)] receives sizable contributions from a whole interval of neutrino
energies of the incoming broad beam. This is shown in the upper plots of Fig.
20 for both MINERvA and T2K CC inclusive samples. We see that for this
selection, the latter experiment is much more dominated by the 1p1h reaction
mechanism (region of $\psi^{\prime}$ around one ) than MINERvA, for which
resonant and DIS modes, located at higher values of the scaling variable,
become more relevant. Moreover, the dispersion of neutrino energies for
MINERvA is quite significant, while the high-energy tail for T2K is less
important. In the lower plots of Fig. 20, we show the corresponding column
(energy-bin) normalized distributions, where the effect of the neutrino flux
times the neutrino cross-section should largely cancel. We observe an almost
universal pattern, which would be the conditional probability
$P(\psi^{\prime}|E_{\nu})$, corrected only by the detector acceptance effects
that are clearly seen at neutrino energies between 0.2 GeV and 3 GeV.

The errors on $R(\psi^{\prime})$ and the covariance matrix for different
values of $\psi^{\prime}$ 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$\pi$ cross-sections.

###  5.1  Results for $R(\psi^{\prime})$

Minerva CC inclusive $R(\psi^{\prime})$ is shown in the left panel of Fig. 21,
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 $\psi^{\prime}$, 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$\pi$ and DIS channels smoothly balance above
$\psi^{\prime}=1$. We show the $R-$correlation matrix in the right plot of
Fig. 21. The correlation is larger than 90% for any pair of of $\psi^{\prime}$
values. This is a consequence of both the initial correlation between
experimental results and the fact that several $\psi^{\prime}-$values
contributes to the same muon-kinematics bin.

We repeat in Fig. 22 the exercise with the CC0$\pi$ MINERvA data. The tendency
observed is similar to that seen for the CC inclusive. MC predicts smaller
cross-section by a similar amount (10%), though in the CC0$\pi$ case some
structure is observed with a small deficit of MC below $\psi^{\prime}=0$. On
the contrary, the region above $\psi^{\prime}=2$ shows an excess on the MC
predictions with respect to data. This region is dominated by 2p2h and $1\pi$
(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 Rodrigues:2015hik . 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$\pi$ ratio plot.

Figure 25:  MINERvA (up) and T2K (bottom) energy $q^{0}$ (left) and momentum
$|\vec{q}\,|$ (right) transfer distributions for the CC inclusive selections.

Figure 26:  The same as Fig. 25 for CC0$\pi$ selections.

Figure 27:  The same as Fig. 25 for CC0$\pi$1p selections.

The $R(\psi^{\prime})$ results for T2K CC inclusive and CC$0\pi$ samples are
shown in Figs. 23 and 24, 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(\psi^{\prime})\sim 1.2$), while the distribution of the ratio as function
of $\psi^{\prime}$ is again rather flat. There is some reflection of the CCRes
contribution in Fig. 23 around $\psi^{\prime}\sim 1.5$, which might be an
indication of an even smaller single pion cross-section. Nevertheless and
despite the large correlations, the evidence is weak, given the flat behavior
of the ratio and the large errors. Similar low cross-section predictions were
also observed in the T2K CC1$\pi$ results Abe:2019arf . The impact of the
possible CC1$\pi$ miss-modelling is minimised in the case of CC0$\pi$ where
the contribution of CC1$\pi$ 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$\pi$ results which exhibit a visible
decrease below $\psi^{\prime}<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. 23 and 24. The correlations for both CC
inclusive and CC$0\pi$ are smaller than in the case of MINERvA, 60%
correlation between values of $\psi^{\prime}$ below 0 and higher than 1.5.
This reduced correlation gives more credibility to the tendencies of
$R(\psi^{\prime})$.

The large errors ($\approx 10\%$ for T2K and $\approx 5\%$ for MINERvA) are
the consequence of large positive correlations enhancing the experimental
errors.

To understand the large correlations across all values of $\psi^{\prime}$ both
for T2K and MINERvA, we estimated the correlation removing all the off-
diagonal terms in the covariance matrices provided by the experiment. In this
case, the ratios $R(\psi^{\prime})$ 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 (see the discussion above of Figs. 19 and 20). The smallest
observed correlation is reduced from 90% to 40% in the case of MINERvA and
from 60% to 5% in the case of T2K. The better figure in the CC0$\pi$ T2K
sample might be a consequence of the narrow neutrino beam at T2K, which allows
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.

##  6 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, CC$0\pi$ and TKI 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 TKI variables found
here and, in general of its CC0$\pi$ data-sample, does not support a large
modification of the 2p2h contribution as requested by the calorimetric
measurements of that collaboration Rodrigues:2015hik . The latter conclusion
agrees with the findings in Ref. Lu:2018stk  with NuWro Golan:2012wx  and
GiBUU Mosel:2019vhx  event generators (GiBUU results can be found in the
Supplemental Material for Ref. Coplowe:2020yea ). However, we should point out
that the re-weight of the 2p2h strength proposed in Rodrigues:2015hik  is
based on the inclusive sample, which is a superset of the CC0$\pi$ and
CC0$\pi$1p data-sets considered in this work.

We have also proposed a novel comparison between flux-folded data and MC
theoretical predictions accumulated in bins of the scaling variable
$\Psi^{\prime}$, 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.

##  Appendix 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},|\vec{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. 25. The
cut off in $|\vec{q}\,|$ implemented in the CC2p2h model is clearly visible in
the MINERvA distributions. The CC2p2h cross-section at the cutoff
$|\vec{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
$|\vec{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.

Figure 28:  $q^{0}$ vs $p_{\perp}$ for MINERvA CC0$\pi$ sample for the same
$p_{\parallel}$ bins as in Fig. 6. The color code represents the number of
events in arbitrary units.  Figure 29:  $|\vec{q}\,|$ vs $p_{\perp}$ for
MINERvA CC0$\pi$ sample for the same $p_{\parallel}$ bins as in Fig. 6. The
color code represents the number of events in arbitrary units.

The predicted energy and momentum transfer distributions for the MINERvA and
T2K CC0$\pi$ event-selections are shown next in Fig. 26. As in the previous
case, the cut off $|\vec{q}\,|\leq$ 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$\pi$ 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
$|\vec{q}\,|$ similar as that observed for the CC inclusive ones.

The $q^{0}$ and $|\vec{q}\,|$ event distributions used for the TKI variable
analysis (CC0$\pi$1p) are shown in Fig. 27. 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$\pi$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.
27, one can also observe a shift in the value of $|\vec{q}\,|$. Actually, the
mean value of $|\vec{q}\,|$ in the MINERvA (T2K) 1p1h component moves from
0.69 GeV (0.61 GeV) in the CC0$\pi$ data-sample to 0.77 GeV (0.73 GeV) in the
CC0$\pi$1p one. This is a consequence of requesting a proton above 0.45 GeV in
the detector.

Finally in Figs. 28 and 29, we show for the MINERvA CC$0\pi$ sample, the
$q^{0}$ and $|\vec{q}\,|$ distributions as a function of $p_{\perp}$ for the
same $p_{\parallel}$ 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. 28
(note the logarithm scale in the $z-$coordinate). This second enhanced region
is less visible in $|\vec{q}\,|$, see Fig. 29, except for the low longitudinal
momentum bins. The distributions also show that $|\vec{q}\,|$ values are
similar for a given $p_{\perp}$ independently of the longitudinal momentum
except for the case with $p_{\parallel}<3$ GeV. The main differences between
data and MC observed in Fig. 6 as function of $p_{\parallel}$ are most
probably caused by the CCRes and CC2p2h large $q^{0}$ contributions.

###### Acknowledgements.

This work was supported by the Swiss National Foundation Grant No.
200021_85012, the Spanish Ministerio de Economía y Competitividad and the
European Regional Development Fund under contract FIS2017-84038-C2-1-P, the EU
STRONG-2020 project under the program H2020-INFRAIA-2018-1, grant agreement
no. 824093 and by Generalitat Valenciana under contract PROMETEO/2020/023. We
thank the MINERVA collaboration for letting us notice the publication of their
new results.

## References

  * [1] L. Alvarez-Ruso et al.  NuSTEC White Paper: Status and challenges of neutrino–nucleus scattering.  Prog. Part. Nucl. Phys., 100:1–68, 2018.
  * [2] K. Abe et al.  A Long Baseline Neutrino Oscillation Experiment Using J-PARC Neutrino Beam and Hyper-Kamiokande.  12 2014.
  * [3] K. Abe et al.  T2K neutrino flux prediction.  Phys. Rev. D, 87(1):012001, 2013.  [Addendum: Phys.Rev.D 87, 019902 (2013)].
  * [4] J. Nieves, F. Sanchez, I. Ruiz Simo, and M. J. Vicente Vacas.  Neutrino Energy Reconstruction and the Shape of the CCQE-like Total Cross Section.  Phys. Rev., D85:113008, 2012.
  * [5] L. Alvarez-Ruso, Y. Hayato, and J. Nieves.  Progress and open questions in the physics of neutrino cross sections at intermediate energies.  New J. Phys., 16:075015, 2014.
  * [6] J. Nieves, Jose Enrique Amaro, and M. Valverde.  Inclusive quasi-elastic neutrino reactions.  Phys. Rev. C, 70:055503, 2004.  [Erratum: Phys.Rev.C 72, 019902 (2005)].
  * [7] J. Nieves, I. Ruiz Simo, and M. J. Vicente Vacas.  Inclusive Charged–Current Neutrino–Nucleus Reactions.  Phys. Rev., C83:045501, 2011.
  * [8] Yoshinari Hayato.  A neutrino interaction simulation program library NEUT.  Acta Phys. Polon., B40:2477–2489, 2009.
  * [9] O. Benhar, A. Fabrocini, and S. Fantoni.  The Nucleon Spectral Function in Nuclear Matter.  Nucl. Phys. A, 505:267–299, 1989.
  * [10] O. Benhar, A. Fabrocini, S. Fantoni, and I. Sick.  Spectral function of finite nuclei and scattering of GeV electrons.  Nucl. Phys. A, 579:493–517, 1994.
  * [11] Omar Benhar, Nicola Farina, Hiroki Nakamura, Makoto Sakuda, and Ryoichi Seki.  Electron- and neutrino-nucleus scattering in the impulse approximation regime.  Phys. Rev. D, 72:053005, 2005.
  * [12] K. Abe et al.  Measurement of inclusive double-differential $\nu_{\mu}$ charged-current cross section with improved acceptance in the T2K off-axis near detector.  Phys. Rev. D, 98:012004, 2018.
  * [13] K. Abe et al.  Simultaneous measurement of the muon neutrino charged-current cross section on oxygen and carbon without pions in the final state at T2K.  Phys. Rev. D, 101(11):112004, 2020.
  * [14] K. Abe et al.  Measurement of the muon neutrino charged-current single $\pi^{+}$ production on hydrocarbon using the T2K off-axis near detector ND280.  Phys. Rev. D, 101(1):012007, 2020.
  * [15] X.G. Lu et al.  Measurement of final-state correlations in neutrino muon-proton mesonless production on hydrocarbon at $\langle E_{\nu}\rangle=3$ GeV.  Phys. Rev. Lett., 121(2):022504, 2018.
  * [16] A. Filkins et al.  Double-differential inclusive charged-current $\nu_{\mu}$ cross sections on hydrocarbon in MINERvA at $\langle E_{\nu}\rangle\sim$ 3.5 GeV.  Phys. Rev. D, 101(11):112007, 2020.
  * [17] D. Ruterbories et al.  Measurement of Quasielastic-Like Neutrino Scattering at $\langle{E_{\nu}}\rangle\sim 3.5$ GeV on a Hydrocarbon Target.  Phys. Rev. D, 99(1):012004, 2019.
  * [18] J.E. Sobczyk, J. Nieves, and F. Sánchez.  Exclusive-final-state hadron observables from neutrino-nucleus multinucleon knockout.  Phys. Rev. C, 102(2):024601, 2020.
  * [19] T. Cai et al.  Nucleon binding energy and transverse momentum imbalance in neutrino-nucleus reactions (see also supplemental material).  Phys. Rev. D, 101(9):092001, 2020.
  * [20] L. A. Harewood and R. Gran.  Elastic hadron-nucleus scattering in neutrino-nucleus reactions and transverse kinematics measurements (MINERvA technical note 91).  arXiv:1906.10576 2019.
  * [21] K. Abe et al.  Characterization of nuclear effects in muon-neutrino scattering on hydrocarbon with a measurement of final-state kinematics and correlations in charged-current pionless interactions at T2K.  Phys. Rev. D, 98(3):032003, 2018.
  * [22] Dieter Rein and Lalit M. Sehgal.  Neutrino Excitation of Baryon Resonances and Single Pion Production.  Annals Phys., 133:79–153, 1981.
  * [23] Dieter Rein and Lalit M. Sehgal.  Coherent pi0 Production in Neutrino Reactions.  Nucl. Phys. B, 223:29–44, 1983.
  * [24] D. Rein and L.M. Sehgal.  PCAC and the Deficit of Forward Muons in pi+ Production by Neutrinos.  Phys. Lett. B, 657:207–209, 2007.
  * [25] M. Glück, E. Reya, and A. Vogt.  Dynamical parton distributions revisited.  Eur. Phys. J. C, 5:461–470, 1998.
  * [26] A. Bodek and U.K. Yang.  Modeling neutrino and electron scattering cross-sections in the few GeV region with effective LO PDFs.  AIP Conf. Proc., 670(1):110–117, 2003.
  * [27] L. L. Salcedo, E. Oset, M. J. Vicente-Vacas, and C. Garcia-Recio.  Computer Simulation of Inclusive Pion Nuclear Reactions.  Nucl. Phys., A484:557–592, 1988.
  * [28] E. S. Pinzon Guerra et al.  Using world charged $\pi^{\pm}-$nucleus scattering data to constrain an intranuclear cascade model.  Phys. Rev., D99(5):052007, 2019.
  * [29] P.A. Rodrigues et al.  Identification of nuclear effects in neutrino-carbon interactions at low three-momentum transfer.  Phys. Rev. Lett., 116:071802, 2016.  [Addendum: Phys.Rev.Lett. 121, 209902 (2018)].
  * [30] L. Aliaga et al.  Design, Calibration, and Performance of the MINERvA Detector.  Nucl. Instrum. Meth. A, 743:130–159, 2014.
  * [31] P. A. Amaudruz et al.  The T2K Fine-Grained Detectors.  Nucl. Instrum. Meth. A, 696:1–31, 2012.
  * [32] C. Andreopoulos et al.  The GENIE Neutrino Monte Carlo Generator.  Nucl. Instrum. Meth., A614:87–104, 2010.
  * [33] X. G. Lu, D. Coplowe, R. Shah, G. Barr, D. Wark, and A. Weber.  Reconstruction of Energy Spectra of Neutrino Beams Independent of Nuclear Effects.  Phys. Rev., D92(5):051302, 2015.
  * [34] X. G. Lu, L. Pickering, S. Dolan, G. Barr, D. Coplowe, Y. Uchida, D. Wark, M.O. Wascko, A. Weber, and T. Yuan.  Measurement of nuclear effects in neutrino interactions with minimal dependence on neutrino energy.  Phys. Rev. C, 94(1):015503, 2016.
  * [35] Andrew P. Furmanski and Jan T. Sobczyk.  Neutrino energy reconstruction from one muon and one proton events.  Phys. Rev. C, 95(6):065501, 2017.
  * [36] R. Gonzaléz-Jiménez, G.D. Megias, M.B. Barbaro, J.A. Caballero, and T.W. Donnelly.  Extensions of Superscaling from Relativistic Mean Field Theory: the SuSAv2 Model.  Phys. Rev. C, 90(3):035501, 2014.
  * [37] T.W. Donnelly and Ingo Sick.  Superscaling of inclusive electron scattering from nuclei.  Phys. Rev. C, 60:065502, 1999.
  * [38] M.B. Barbaro, R. Cenni, A. De Pace, T.W. Donnelly, and A. Molinari.  Relativistic y - scaling and the Coulomb sum rule in nuclei.  Nucl. Phys. A, 643:137–160, 1998.
  * [39] Tomasz Golan, Cezary Juszczak, and Jan T. Sobczyk.  Final State Interactions Effects in Neutrino-Nucleus Interactions.  Phys. Rev. C, 86:015505, 2012.
  * [40] Ulrich Mosel.  Neutrino event generators: foundation, status and future.  J. Phys. G, 46(11):113001, 2019.
  * [41] D. Coplowe et al.  Supplemental Material For: Probing nuclear effects with neutrino-induced charged-current neutral pion production.  Phys. Rev. D, 102(7):072007, 2020.