Neutrino Energy Reconstruction from Semi-inclusive Samples

PHYSICAL REVIEW C 105, 025502 (2022) Neutrino energy reconstruction from semi-inclusive samples R. González-Jiménez ,1M. B. Barbaro ,2,3J. A. Caballero ,4,5T. W. Donnelly ,6N. Jachowicz ,7 G. D. Megias ,4,8K. Niewczas ,7,9A. Nikolakopoulos,7J. W. Van Orden,10 and J. M. Udías 1 1Grupo de Física Nuclear, Departamento de Estructura de la Materia, Física Térmica y Electrónica and IPARCOS, Facultad de Ciencias Físicas, Universidad Complutense de Madrid, CEI Moncloa, Madrid 28040, Spain 2Dipartimento di Fisica, Università di Torino and INFN, Sezione di Torino, Via P. Giuria 1, 10125 Torino, Italy 3Université Paris-Saclay, CNRS/IN2P3, IJCLab, 91405 Orsay, France 4Departamento de Física Atómica, Molecular y Nuclear, Universidad de Sevilla, 41080 Sevilla, Spain 5Instituto de Física Teórica y Computacional Carlos I, Granada 18071, Spain 6Center for Theoretical Physics, Laboratory for Nuclear Science and Department of Physics, Massachusetts Institute of Technology, Cambridge, Massachusetts 02139, USA 7Department of Physics and Astronomy, Ghent University, Proeftuinstraat 86, B-9000 Gent, Belgium 8Research Center for Cosmic Neutrinos, Institute for Cosmic Ray Research, University of Tokyo, Kashiwa, Chiba 277-8582, Japan 9Institute of Theoretical Physics, University of Wrocław, Plac Maxa Borna 9, 50-204, Wrocław, Poland 10Department of Physics, Old Dominion University, Norfolk, Virginia 23529, USA and Jefferson Laboratory, 12000 Jefferson Avenue, Newport News, Virginia 23606, USA (Received 4 April 2021; accepted 4 February 2022; published 23 February 2022) We study neutrino-nucleus charged-current reactions on finite nuclei for the situation in which an outgoing muon and a proton are detected in coincidence; i.e., we focus on semi-inclusive cross sections. We limit our attention to one-body current interactions (quasielastic scattering) and assess the impact of different nuclear effects in the determination of the neutrino energy. We identify the regions in phase space where the neutrino energy can be reconstructed relatively well and study whether the cross section in those regions is significant. Our results indicate that it is possible to filter more than 50% of all events according to the muon and proton kinematics, so that for the DUNE and T2K fluxes the neutrino energy can be determined with uncertainties of less than 1% and 3%, respectively. Furthermore, we find that the reconstructed neutrino energy does not depend strongly on how one treats the final-state interactions and is not much affected by the description of the initial state. On the other hand, the estimations of the uncertainty on the neutrino energy show important sensitivity to the modeling of the initial state. DOI: 10.1103/PhysRevC.105.025502 I. INTRODUCTION The T2K and MINERvA experiments and those based on Liquid Argon Time Projection Chamber (LArTPC) detectors have shown their capabilities to measure the final-state lepton (μ±or e±) and to identify one or more charged particles in coincidence [1–5]. Future experiments, such as DUNE [6], will incorporate an enhanced tracking capability for hadrons in the final state. Also, it is worth mentioning the SK-Gd project [7], that improves the detection and identification capabilities of neutrons by adding Gd salts to the Super Kamiokande water tank. Compared to inclusive experiments, where only the final lepton is detected, the additional information about the hadrons in the final state, namely, semi-inclusive scattering, will improve the reconstruction of the incoming neutrino energy. In Refs. [8–11] explorations of the possibilities arising from the extended knowledge of the final state, specifically focusing on events where there is simultaneous detection of the lepton and a nucleon, were presented. In Ref. [10], it was proposed to study the average neutrino energy corresponding to a given semi-inclusive event E=dEEφ(E)d6σ(E) dldkldNdp N dEφ(E)d6σ(E) dldkldNdp N ,(1) where d6σ(E) dldkldNdp Nis the sixfold differential cross section for a fixed neutrino energy Eand fixed muon and final nucleon kinematics. φ(E) is a given flux distribution, normalized as dEφ(E)=1. Furthermore, the standard deviation for the average neutrino energy can be obtained from the first and second statistical moments E=E2−E2,(2) where E2=dE E2φ(E)d6σ(E) dldkldNdp N dEφ(E)d6σ(E) dldkldNdp N .(3) Thus, provided that expressions for the flux and cross section are known, and given that the 4-momenta of the lepton 2469-9985/2022/105(2)/025502(16) 025502-1 ©2022 American Physical Society R. GONZÁLEZ-JIMÉNEZ et al. PHYSICAL REVIEW C 105, 025502 (2022) and proton in the final state are both measured, the average neutrino energy will be defined up to E±E.Fromnow on, we will refer to the average energy [Eq. (1)] as reconstructed energy, since it is an estimator of the most likely neutrino energy associated with that event. Notice that in the kinematical regions where the mere detection of the final lepton and nucleon determines the incoming neutrino energy with good accuracy, any adequate neutrino energy estimator will provide results very similar to the average energy derived from Eq. (1). The impact of different assumptions for the nuclear models involved on the neutrino energies and uncertainties were studied in Ref. [10]. To summarize the conclusions in said reference, we note the following: (1) The reconstructed neutrino energy Edepends only moderately on the nuclear model introduced in Eq. (1). (2) The corresponding uncertainty of the reconstructed energy does depend on the nuclear model, but it may be relatively low for a large fraction of the events. The basis for these observations is the fact that the neutrino energy gets essentially blurred by the missing energy of the nuclear system, that is, the energy required to knock out the observed nucleon, while the nuclear recoil is generally very small. In light nuclei, such as 12Cor 16O, for many events the main contribution will come from the nucleons in the p shell(s), for which the missing energy is a rather well-known quantity. In this paper, we extend the previous analysis presented in Ref. [10] by examining the whole phase space and scrutinize in detail the potential for model-independent neutrino-energy determinations. We consider semi-inclusive reactions involving an incident neutrino followed by detection of a charged lepton and a nucleon in the final state together with no produced pions; that is, we focus on events of the type CC1μ1p0π, having chosen in the present work to emphasize muons and protons in the final state. As discussed later, this selection of events does not mean that one and only one nucleon is assumed to be present in the final state, only that at least one is present. Indeed, depending on the kinematics chosen, there must be other nucleons beyond the one actually detected. The study presented in Ref. [10] is extended here to include various models that treat the issue of hadronic final-state interactions. We discuss a typical situation, that is, we make specific choices for the measured 4-momenta, in order to orient the reader to the basic characteristics of semi-inclusive reactions before going on to analyze a broad region of the full phase space. We first introduce the semi-inclusive kinematics and cross section in general terms (Sec. II), and then particularize for the quasielastic (QE) interaction (Sec. II A). In Sec. III,thesemiinclusive cross section is studied for a fixed set of kinematics. Full phase space results are shown in Sec. IV. In Secs. IV A and IVB, we assess the effect of final-state interactions and the description of the initial state on the neutrino energy determination. In Sec. IV C, we show the regions of phase space where the neutrino energy is reconstructed with the lowest error. Finally, we draw our conclusions in Sec. V. II. KINEMATICS AND CROSS SECTION For the discussion that follows, we will assume that the final-state lepton (here a muon is assumed) with 4-momentum (El,kl) and a nucleon (here a proton is assumed) with 4momentum (EN,pN) are detected in coincidence. No other particles are assumed to be detected, although, depending on the specific kinematics assumed, they must be present (see below). We work in the laboratory frame where the target nucleus is at rest, the incoming neutrino momentum is along ˆz, and the lepton kinematical variables are contained in the ˆx-ˆzplane. The angle between the incident neutrino and the outgoing lepton is θl, while in the chosen coordinate system the polar and azimuthal angles that specify the direction of the outgoing nucleon are θNand φN, respectively. The magnitude of the nucleon’s 3-momentum is given by pN=|pN|. Apart from the detected nucleon, the hadronic final state contains an undetected hadronic system having missing 4-momentum (EB,pB), namely, a total energy of EBand a missing 3momentum pB≡pm. If one denotes by qthe 3-momentum transferred from the leptons to the hadronic system, one has pm=q−pN.(4) The undetected hadronic system has invariant mass MB(M0 B at threshold with MB⩾M0 B) and total energy EB=TB+MB=(MB)2+pm2,(5) which defines the kinetic energy of the unobserved final-state system, TB.FromEq.(5) one has EB=E−El−TN+M0 A−mN,(6) where M0 Ais the target ground-state mass and TN=EN−mN is the kinetic energy of the detected nucleon. This leads to an expression for the so-called missing energy, Em=MB−M0 B+Es=E−El−TN−TB,(7) where Es=M0 B+mN−M0 Ais the separation energy and the (typically very small) recoil kinetic energy difference has been neglected. Clearly, if one knew the missing-energy Emthen the incident neutrino energy Ewould also be known. The magnitude of this missing momentum pmis given by pm=k2+kl2+p2 N−2kklcos θl−2kpNcos θN +2klpN(cos θlcos θN+sin θlsin θNcos φN)1 2.(8) Depending on the specific kinematics, i.e., the value of the missing energy, the residual system may be the daughter nucleus in its ground state (this defines the threshold for the semi-inclusive reaction to become possible), or it may be in a discrete excited state (these are states in the residual nucleus that lie below the threshold where a second nucleon can be ejected), and, while they de-excite by γdecay, that process is slow on the nuclear timescale and thus these states may be treated effectively as stationary states. Then, at a well-defined threshold, a second nucleon must be emitted (this is not optional: there are no nuclear states involving one nucleon and a residual bound nucleus above this point); and so on with more particles in the final state in addition to the one special nucleon that is assumed to be detected. At even larger missing energy 025502-2 NEUTRINO ENERGY RECONSTRUCTION FROM SEMI- … PHYSICAL REVIEW C 105, 025502 (2022) (roughly 140 MeV), pion production becomes possible (still with the lepton and one nucleon assumed to be detected) and beyond where more particles may be present in the undetected part of the final-state system. We consider neutrino energy distributions from DUNE and T2K fluxes. The flux-averaged semi-inclusive cross section for this process is given by d6σ dkldldp NdN=dEφ(E)d6σ(E) dkldldp NdN .(9) With the neutrino energy and missing energy related through energy-momentum conservation, it is convenient to change the integration variable from Eto Em, thereby writing the expression for the semi-inclusive cross section in a more familiar way, found for instance in inclusive and exclusive electron scattering. The semi-inclusive cross section is then given by d6σ dkldldp NdN=dEmφ(E)Fk2 lp2 NM∗ B (2π)5EBfrec μνHμν, (10) with frec = 1−E−ˆ kν·(pN+kl) EB (11) being the nuclear recoil factor, where ˆ kνis a unit vector along the neutrino beam line. Fis given by F=GF √22 cos2θc,(12) where GFis the Fermi constant and θcis the Cabibbo mixing angle. The charged-current lepton tensor is μν =2 EElKi,μKf,ν +Ki,ν Kf,μ −gμν KiKf−ihμναβ Kα iKβ f, (13) with h=+1 for antineutrinos and −1 for neutrinos. The terms Kiand Kfrepresent the four-momenta of the initial and final leptons involved in the process. Up to this point the discussion is general; i.e., it does not depend on the nuclear model or the reaction channel. All the complexity regarding the hadronic part of the interaction is in the hadron tensor Hμν , which is discussed in Sec. II A. Equation (6) tells us that the neutrino energy can be reconstructed very well from a semi-inclusive sample of events when it is dominated by a narrow and well-known missingenergy region. This should be the case for QE scattering, where the neutrino scatters elastically from a bound nucleon such that the missing energy is of the order of the binding energy of the nucleon. In the following, we focus on the QE reaction and we describe the content of the different nuclear models employed in this work. A. Nuclear models for quasielastic scattering We focus on neutrino-induced charge-current QE scattering, with one boson exchanged between lepton and one-body hadron currents, within the impulse approximation. Thus, we do not consider meson-exchange currents (MEC) nor the processes where real pions may be produced in the final state. In this work, we restrict our attention to 16O, although analyses along these lines can easily be performed for other nuclei [12]. We will describe the initial state as a set of relativistic mean-field (RMF) wave functions that correspond to different shells labeled by the relativistic quantum number κ.The hadron tensor for a given shell is Hμν κ=ρκ(Em) mj,sNJμ κ,mj,sN(Q,PN)∗Jν κ,mj,sN(Q,PN), where Q=Ki−Kf,ρκ(Em) is the missing-energy density, Jμ κ,mj,sNis the hadron current in momentum space defined as Jμ κ,mj,sN(PN,Q)=dpsN(pN,q+p)Oμ(Q)mj κ(p),(14) mjis the third component of the total angular momentum jof the bound nucleon, and sNis the spin projection of the final nucleon. functions are relativistic independentparticle wave functions describing the bound and scattered nucleon and Oμis the usual boson-nucleon-nucleon operator, for which we use its CC2 form (see Refs. [13–15] for details). In the framework of a pure shell model, one simply has ρκ(Em)=δ(Em−Eκ m), where Eκ mis the energy eigenvalue for a given shell. This missing-energy distribution of the shell model, however, is only a first approximation to the one in real nuclei, that for the valence shells has been measured in quite a number of electron scattering experiments. Generally speaking, effects beyond mean field, such as shortand long-range correlations, modify the actual missing-energy distribution predicted by the shell model. In this work, as already done in Ref. [10], we also take as a reference the spectral function formalism. The spectral function incorporates the probability of finding the nucleon in the initial state with certain energy and momentum. It includes the depletion of the occupation of the shell-model states and the appearance of nucleons at deeper (namely, higher) missing energy, in both cases due to correlations, both long and short ranged. For the purpose of this work, we will consider the Rome spectral function [16,17] as a fair representation of the missing energy and missing-momentum distribution of the nucleons in the target nuclei, as measured in electron scattering experiments. Other spectral function calculations are available in the literature [18,19]; however, since all of these have been constrained to some extent to reproduce the (e,ep) electron scattering cross sections, their results will not differ much from the ones here. The spectral function is easy to incorporate in a fully factorized, plane-wave calculation (as in Refs. [10,20]), with the exclusive 6-differential cross section given by d6σ dkldldp NdN=KS(Em,pm)σνN.(15) Kis a function containing kinematic factors, S(Em,pm)isthe initial-state spectral function, and σνNis the charged-current elastic neutrino-nucleon cross section for an off-shell nucleon with initial momentum pm(see Refs. [10,17,20] for details). 025502-3 R. GONZÁLEZ-JIMÉNEZ et al. PHYSICAL REVIEW C 105, 025502 (2022) FIG. 1. (a) The Rome spectral function (integrated over pm) as a function of the missing energy. Our parametrization of the background is represented by the dashed line. (b) Momentum distributions from the Rome spectral function and from our representation. In this work, however, we use a representation of the spectral function amenable to our relativistic distorted-wave and unfactorized calculations that is sufficient to achieve the goals of this work. Thus, we divide the missing-energy phase space into several regions. In the lowest-energy region, below the two-nucleon threshold, we will have the p-shell states, with an energy dependence not given by δfunctions but rather by the energy distribution seen in the spectral function [Fig. 1(a)]. We identify each shell with a different region of missingenergy following the analysis presented in Ref. [10]. The high Emand pmpart of the spectral function due to correlations is accounted for by introducing an swave [21,22], broad in momentum space (narrow in rspace, approximately 0.85 fm), that is fitted to reproduce the momentum distribution of the spectral function in this region. The momentum distribution obtained in this representation, compared to the one of the Rome spectral function, is shown in Fig. 1(b). The specific regions and occupation numbers are summarized in Table I. Above the two-nucleon knockout threshold, the independentparticle shells and the background coexist. To account for it, we have parameterized the missing-energy profile of the background in the region 25 <Em<100 MeV [dashed blue line in Fig. 1(a)]. In the region Em>100 MeV, we assume that there is only background, which is well described by an exponential TABLE I. Correspondence between missing-energy regions and shells in oxygen. In the last column are the occupation numbers. Em(MeV) Shells 16O 0–16.5p1/21.51 16.5–25 p3/23.47 25–100 s1/2+backg. 2.22 s1/21.62 backg. 0.60 100–300 backg. 0.80 fall-off. The explicit expressions for these functions are given in the Appendix. The representation of the initial state we use in our modeling essentially contains the same (albeit somewhat simplified) missing-energy and momentum structure of the Rome spectral function. Indeed, despite the fact that our calculations are unfactorized (due to the relativistic effects [23,24] and eventually FSI (final-state interactions) [14]), our cross section results, when FSI are neglected, are within few percent of the ones obtained with the fully factorized calculation based on the spectral function approach of the Rome model (see, e.g., Fig. 2), showing that ingredients preventing factorization and negative-energy components have a small effect on the FIG. 2. Single-differential cross sections for the DUNE (a) and T2K (b) fluxes with the models discussed in the text. 025502-4 NEUTRINO ENERGY RECONSTRUCTION FROM SEMI- … PHYSICAL REVIEW C 105, 025502 (2022) cross sections computed here [25]. This means that our results are representative of what MC (Monte Carlo) event generator based on the spectral function+factorized calculations, even considering FSI, would produce. To summarize, in our rationale (impulse approximation) the spectral function is a reasonably realistic representation of the initial state (energy and momentum) of the nucleon which will contribute, for the reaction at hand, to a final state in which we will see at least one knocked-out nucleon. We want to emphasize that in this way, we incorporate the experimental constraints provided by electron scattering experiments on missing-energy and momentum distribution in the initial nucleon, certainly much better than within any Fermi gas approach or a pure shell model. In the following, we discuss the description of the final state that we incorporate in our calculations. In a way, the final state is to a large extent determined by the experimental signal that is to be described: Is there only a proton and no other hadrons in the final state? Or, at the other extreme, does the experimental signal contain every event for which at least one proton is seen? Generally, the actual experimental situation will be a combination of these two extreme, simplified, cases. We will look for some representation of these situations, in order to study the effects of the different definitions of the final state on the reliability of the determination of the neutrino energy. (1) Real relativistic optical potential (rROP) or energydependent relativistic mean-field (ED-RMF). In this case, the final-state nucleon is a solution of the Dirac equation with a real potential, and the absence of an imaginary part in the potential means that no flux is lost. We will use both the rROP and the ED-RMF, the difference between them being in the relativistic mean-field potential seen by the final nucleon [26]. In the rROP case, we use the real part of the energydependent A-independent oxygen (EDAI-O) optical potential [27], while the ED-RMF is the RMF potential (the same as for the bound state) but multiplied by a phenomenological function that weakens the potential for increasing nucleon momenta [28]. The nucleon wave functions within the ED-RMF model are eigenstates of the same Hamiltonian, and therefore, orthogonality between initial and final state is satisfied; i.e., Pauli blocking is consistently incorporated. This orthogonality is not as good in the rROP model; therefore, one should be cautious when the momentum of the nucleon is smaller than approximately pN<300 MeV. For pNlarger than around 400 MeV, the overlap between the initial and final states is negligible and hence, orthogonality is not an issue. These approaches, in the pure shell-model case, have been shown to successfully describe inclusive scattering data for both neutrinos and electrons [26,28]. In the case of semi-inclusive scattering, as considered here, these models would provide an estimate of the situation in which the hadronic final-state signature consists of at least one proton. There may be other nucleons that arise from correlations in the initial state or other hadrons, such as nucleons or pions, produced during the interaction of the knocked-out nucleon(s) with the residual system. The scenarios with two or more nucleons in the final configurations necessarily arise from kinematics in which the missing energy, Em, is above the two-nucleon knockout threshold. (2) Elastic-only channel in the FSI, represented by a complex relativistic optical potential (ROP). The whole ROP which is fitted to reproduce elastic protonnucleus scattering, and that contains real and imaginary parts, is employed in this case. Hence, this calculation allows us to estimate the probability that the (primary) nucleon knocked out during the interaction with the boson propagates through the residual system with elastic scattering only. This primary nucleon does not knock out other nucleons or create new hadrons in its way out, nor does it lose energy in any way apart from elastic recoils. The angle of the nucleon can change, though. This situation can be considered equivalent to running the cascade models retaining only elastic interactions. However, the calculation presented here is of course a fully quantum mechanical one. The loss of flux implied by the imaginary part of the potential would lead to a strong underestimation of the inclusive cross section, in which the outgoing nucleon remains undetected. Thus, the ROP estimation would be more in line with an experimental signature of having one proton detected and no other hadron. Additional hadrons, however, could appear due to correlations in the initial state and subsequent FSI of the secondary nucleon. The “one proton and one only” signature might be enforced in the calculation by keeping the missing-energy below the two-nucleon emission threshold, for the initial state, and the elastic (full optical potential) condition in the final state. These ingredients, but with the pure shell model, have been widely applied to analyze exclusive (e,ep) data on different target nuclei [13,14,29]forwhich there is certainty that there is only one proton in the final state. The theoretical prediction is scaled to the data by the spectroscopic factor and the agreement with data is outstanding. As in this work we use the spectral function representation, the scale factor is already included in the theory. Were we to use these ingredients to analyze the exclusive data, we would find very good agreement for the p1/2shell. For the p3/2shell, one has to consider that with the binning of the spectral function used here, the small states located at missing energies around the main p3/2state are all summed up. When all these states are taken into account, the agreement is good, at least for low-to-moderate values of Em. (3) Relativistic plane-wave impulse approximation (RPWIA): The final nucleon is described by a relativistic plane wave, so final-state interactions and Pauli blocking effects are neglected. Although RPWIA is an oversimplified description of the process that is not suitable for some experimental situations, it is a very 025502-5 R. GONZÁLEZ-JIMÉNEZ et al. PHYSICAL REVIEW C 105, 025502 (2022) common first-order theoretical estimate that can help in improving our understanding of the dynamical properties involved in semi-inclusive processes. We include it here as a reference. (4) In addition to the previous models, and as a further reference, we present here our calculations in the factorized spectral function approach (SFA) described by Eq. (15). The SFA is completely or partially implemented in some of the Monte Carlo event generators used in neutrino experiments [30,31] and thus it is useful for reference. In Fig. 2, we compare the predictions of the models described above for the single-differential cross sections as a function of the final nucleon momentum for the DUNE and T2K fluxes. As expected, our representation of the spectral function calculation in RPWIA and the factorized SFA results are essentially identical and yield the largest cross sections, due to the absence of Pauli blocking and of the distortion, that shifts the strength to regions kinematically suppressed [28,32]. The ratios of the peak values of the rROP results to those of the RPWIA in the two panels in the figure are 0.84 [Fig. 2(a)] and 0.82 [Fig. 2(b)]. The rROP and ED-RMF results are very similar both in shape and magnitude. The ROP result has a very similar shape to that of the rROP, albeit smaller in magnitude as a consequence of the restriction to elastic-only propagation for the knocked-out proton: The ratios at the peaks are 0.71 [Fig. 2(a)] and 0.75 [Fig. 2(b)]. These results are consistent with those found in previous studies, e.g., Ref. [15]. The main idea in this work is to compare the accuracy in the energy determination performed based on cross sections and kinematics (missing-energy and momentum) given by these models, which either correspond approximately to possible experimental samples or signatures or to common ingredients of MC event generators. We will first do so for a fixed final-state kinematics (Sec. III) and then extend the analysis to the full phase space (Sec. IV). III. SELECTED EVENTS AND TRAJECTORIES To get a deeper understanding of the procedures used to reconstruct the neutrino energy and the definition of its error, we analyze in detail the semi-inclusive cross section for a fixed set of kinematics. After this discussion, which should help in understanding the basic characteristics of semi-inclusive reactions of the type considered in the present study as well as those involving electron scattering, specifically, (e,ep) reactions, in the following section we present an analysis that extends over the full phase space. From Eqs. (7) and (8), it is clear that for fixed values of the observable parameters the values of Emand pmare determined for each value of E. Thus, for fixed values of the observables, the integral over Ein Eq. (9) follows a curve or trajectory in the pm-Emplane. In Fig. 3, trajectories are shown for selected kinematics for the detected particles, namely, El=3800 MeV, θl=7deg,TN=140 MeV, φN=180 deg. This represents a “typical” situation and clearly many others could have been chosen to illustrate the basic behavior to be expected. What FIG. 3. The Em−pmtrajectories are shown for selected “typical” kinematics: El=3800 MeV, θl=7deg,TN=140 MeV, φN= 180 deg. Each line corresponds to a different value of the proton scattering angle θN(in degrees). Here, we plot the Rome spectral function as a background to allow one to easily identify the different regions of the spectral function that are crossed by the trajectories. is varied here is the polar angle for the detected proton, θN. As stated above, as one goes along a given trajectory the neutrino energy Evaries. It starts at the lower boundary which defines the threshold for the semi-inclusive reaction to occur and grows as the missing energy increases. The semi-inclusive cross section thereby produced for a particular neutrino flux is then obtained as a line integral along the specific trajectory. In effect, each event where a muon and proton are detected in coincidence corresponds to a specific trajectory. Only were one to have a monoenergetic neutrino beam would a point in the Em-pmplane be selected; however, with broad-band beams a weighted line integral is required. That said, one still sees a striking pattern to the behavior one should expect when performing such line integrals. The strength in the Rome spectral function, which should provide a good starting point for the characteristics to be expected in semi-inclusive reactions, is extremely localized. One sees the largest concentration of strength where the pshells are located (at around Em=20 MeV) with less where the broad sshell is located (at around Em=50 MeV); at still larger values of Em(and pm) the Rome group spectral function does have some strength, although it is spread over a wide region in the Em-pmplane and is too small to be seen in this representation. Furthermore, we note that pion production cannot occur until one reaches Em∼mπand that it is not appreciable until Em∼m−mN∼300 MeV. Of course, in the present work we are considering only events that have no pions. We do have some knowledge about this generic behavior of the distribution of strength from inclusive electron scattering, (e,e). Inclusive scattering corresponds to performing integrals over specific regions in the Em-pmplane [33–35]. A very similar pattern is expected in that case and what is found for 16Ois that somewhat over 50% of the inclusive cross section stems from the pshells, about 25% comes from the s-shell region and the rest comes from a broad region 025502-6 NEUTRINO ENERGY RECONSTRUCTION FROM SEMI- … PHYSICAL REVIEW C 105, 025502 (2022) FIG. 4. The sixfold differential cross section computed with the DUNE flux and the SFA model is shown as a function of the missing energy Emon linear (a) and semilog (b) scales. Muon and proton variables are fixed to El=3800 MeV, θl=7 deg, TN=140 MeV, φN=180 deg, as in Fig. 3. at higher missing-energy. This is borne out in semi-inclusive (e,ep) electron scattering studies, although only a few data exist in that case. Thus, we expect this general picture to be the case for CCνreactions. Indeed, in the semi-inclusive case that forms the focus of the present work, we expect that the line integrals discussed above are at their largest when the trajectories cross the peaks in the p-shell region together with somewhat reduced strength coming from the s-shell region and a negligible amount arising from the higher Emregion. Given that this is the case, we can then also expect that the optical model based approach discussed above is a reasonable one, whereas were the high-Emregion to be important this would not be obvious. In passing we note that while other choices of two variables to replace Emand pmcan of course be made, the generic behavior seen here strongly suggests that the present choice is a good one and that other choices may not reflect the highly localized nature of the nuclear response. Let us next see what the various models yield for the weighted cross sections. In Fig. 4, we represent the integrand of Eq. (9) along the trajectories shown in Fig. 3, i.e., the sixfold differential cross section for fixed muon and proton kinematics as a function of the missing energy. By varying the angle θN, the cross section changes its magnitude and its shape but, in general, we observe a profile that resembles the ρ(Em) function used in our model [Fig. 1(a)], namely, two prominent peaks corresponding to the pshells, a wide belly for the sshell, and a background that extends up to high missing energies. Clearly, as expected, the p-shell strength is largest, the s-shell strength is smaller, and the high-Em strength is completely negligible, being down by several orders of magnitude. By examining these results in the light of the trajectories shown in Fig. 3we see that the general behavior we expect to occur is borne out. For example, the trajectories for θN=60 and 80 deg both pass through the p-shell region near its peak. However, one trajectory intersects the s-shell region more than the other one does and this results in relatively different amounts from the sshell compared with the pshell. Or, if one has events that correspond to large values of θNwith the chosen kinematics introduced above, then the cross sections are very small, as they should be, since neither the pnor s-shell regions are crossed. Thus, we are interested in events that fulfill two conditions: (1) The neutrino energy needs to be reconstructed rather well. For this to happen, most of the strength should be concentrated in a small Emregion. (2) The cross section is large and, hence, the probability of finding events around such kinematics is high. Given the structure of the spectral function, it seems very likely that for events in which the pshells dominate these requirements will be fulfilled. Another reason that makes the “p-shell dominated events” particularly interesting is that in 16Othe pshells are below the two-nucleon emission threshold. Hence, this is the part of the spectral function that can be determined experimentally from exclusive (e,ep) experiments and, therefore, is well constrained. By looking at Fig. 3, the trajectories of greatest interest are those in the region around 80 <θ N<100 deg; however, for larger θNthe size of the cross section starts to decrease rapidly. To quantify and better assess these results, for each of the curves shown in Fig. 3, we present in Table II the TABLE II. The kinematics chosen here are the following: El= 3800 MeV, θl=7deg,TN=140 MeV, and φN=180 deg and different θN(first column). We show the mean neutrino energy E, its relative 1σerror E/E, the cross section (c.s.), and the weight of the pshells [r, defined in Eq. (16)]. We have used the DUNE flux and the SFA model. θNEE/Ec.s. p-shell (deg) (MeV) (%) (10−42 cm2/MeV2) weight (r) 40 3976 0.80 0.036 0.60 50 3970 0.53 0.30 0.65 60 3974 0.51 0.65 0.48 70 3980 0.51 0.56 0.29 80 3965 0.44 0.68 0.76 90 3964 0.56 0.24 0.86 100 3981 1.2 0.026 0.66 110 4043 1.6 0.0038 0.11 120 4061 1.5 0.0014 0.0025 025502-7 R. GONZÁLEZ-JIMÉNEZ et al. PHYSICAL REVIEW C 105, 025502 (2022) FIG. 5. 2D histogram with the number of events in bins of the relative error (%) and p-shell weight. The DUNE and T2K fluxes were employed in panels (a) and (b), respectively. The calculations correspond to the rROP model, but similar results are found with the other approaches. Brighter areas correspond to bins with more events. values of E,E/E, the cross section (integrated over Em), and the weight of the p-shell region (0 <Em<25 MeV).1 For the particular kinematics studied here, the error in the reconstructed energy is small in all cases. Also, one does not see a simple correlation between p-shell dominance and small error. For example, for θN=40 deg the error is larger than for θN=60 deg, while the weight of the pshells is smaller in the latter case. However, we should not draw general conclusions from the study of just one particular choice of kinematics. Thus, in the next section we address these and other questions in a more systematic way by analyzing the whole phase space. IV. FULL PHASE-SPACE RESULTS In this section, we extend the previous analysis (restricted to very particular kinematics) to the whole phase space, using the DUNE and T2K neutrino fluxes, peaked around E≈2.5 GeV and E≈0.6 GeV, respectively. Also, we study the dependence of the outcomes upon the different models described in Sec. IV. With that purpose in mind, we populate the whole phase space with a few millions of events distributed following the sixfold differential semi-inclusive cross section [Eq. (9)] given by each of the models mentioned above, and for each event we compute the average neutrino energy (E) and its error (E) according to Eqs. (1) and (2). We stress that the analysis presented here does not take into account non-QE interactions that can contribute to the one muon–one proton sample. Furthermore, it does not account for detector efficiency and resolution. For these reasons, 1We compute the weight of the pshells with respect to the full cross section as r=25 0dEmφ(E)d6σ(E) dldkldNdp N 300 0dEmφ(E)d6σ(E) dldkldNdp N ,(16) with Emin MeV. the error in the reconstructed neutrino energy reported here should be understood as a lowest intrinsic bound. Also, the reconstructed energy will likely change when non-QE processes are explicitly considered in the cross section. It is important to realize that the uncertainty in the neutrino energy for these semi-inclusive experiments is basically given by the relevant range of missing energy. We can see that for QE processes (only nuclear excitations and no pions produced), from about Em>200 MeV the cross section has fallen by several orders of magnitude, and hence the neutrino energy is effectively restricted to a limited region of Em. Thus, in absolute terms, DUNE and T2K events have similar errors in the neutrino energy, while in relative terms, the errors for the DUNE events are smaller, owing to the larger neutrino energies in the DUNE flux. This is clearly shown in Fig. 5 that presents a 2D histogram with the number of events in bins of E <E>and p-shell weight (r). We see that the bulk of the events (intense yellow-orange regions) concentrates in a small region corresponding to 0.6<r<0.9 and E <E><1% (1 <E<3%) for the DUNE (T2K) flux. This means that for most of these small-error events, it is likely that the detected proton was ejected from a pshell. As noted earlier, if one neglects the nuclear recoil, the missing energy is trivially defined from the reconstructed neutrino. Accordingly, we define the reconstructed missing energy as Em=E−El−TN.(17) In Fig. 6, we show the single-differential cross section as a function of Emfor the five models and the DUNE and T2K neutrino fluxes. All models show similar shapes and significant differences in the size of the cross section, in line with the results discussed in Fig. 2. It is interesting that for all cases one observes a clear dominance of the Em<40 MeV region, corresponding to the p-shell region. In Fig. 7(a), we show the cumulative distributions as a function of the relative error in the reconstructed neutrino energy. We observe that around 50% of the events have an error lower than 1% for the DUNE flux (thicker lines). For the T2K flux (thinner lines), the relative errors are somewhat 025502-8 NEUTRINO ENERGY RECONSTRUCTION FROM SEMI- … PHYSICAL REVIEW C 105, 025502 (2022) FIG. 6. Single-differential cross section as a function of the reconstructed missing-energy defined in Eq. (17), for the DUNE (a) and T2K (b) fluxes. larger and for around 47% of the events E E<3%. Although for all models we see a similar trend, it seems that the elasticonly ROP estimation has a slightly better capability for good reconstruction of the neutrino energy. In this case, we find that around 60% (51%) of the events have a neutrino energy uncertainty below 1% (3%) for the DUNE (T2K) flux. This is somewhat expected, as this calculation corresponds to a signal very much enriched in just one proton events. The trade-off is, of course, that fewer events will qualify in the first place. It is also interesting to look at these cumulative distributions as a function of the absolute errors. These are presented in Fig. 7(b). For the two fluxes, the majority of the expected events (around 80%) has an error between 15 and 40 MeV, consistent with a majority of events coming from the p-shells. To end this section, we comment that the estimations we made here from our theoretical prescription not only allow one to predict how many events there may be in an experiment leading to a given uncertainty, but they also allow one to identify where in the phase space these events lie, as we show in Sec. IV C. But first, we will study the dependence of these predictions on the ingredients of the models. A. Dependence of the reconstructed energy and its error on the final-state interactions To study how the reconstructed neutrino energy Edepends on the model, we compute for each event (that is, for each set of muon and proton kinematics) a systematic error that quantifies the deviation of the predictions between a given model and a reference one. We choose the rROP model as reference, although the conclusions are independent of this choice. Thus, for each event we compute EFSI,i =1 2|ErROP −Ei|,(18) where the index irefers to RPWIA, ROP, ED-RMF, or factorized SFA. Notice that in this section, all the results share the same description of the initial state, except for the case of SFA, for which it is just slightly different due to factorization, the absence of negative-energy components, and the representation of the spectral function used in our calculation. Hence, by comparing the results of these models we are actually evaluating the impact of FSI (and Pauli blocking effects) on the quality of the neutrino energy reconstruction; this is the reason for the index FSI in the previous equation. The cumulative distributions as a function of EFSI are shown in Fig. 8. These results clearly show that the FSI uncertainty, evaluated via Eq. (18), is very small. EED-RMF and EFSI,ROP are, as expected, very small since the cross section shapes in these models are similar to the reference one, rROP. The largest EFSI is found for the SFA model. This is due to the differences in the initial state, the factorization assumption, and the lack of FSI and Pauli blocking. In any case, for nearly 98% (90%) of the events, the FSI error remains FIG. 7. (a) Cumulative distributions as a function of the relative error (in percentage) in the reconstructed neutrino energy [error given by Eq. (2)]. The thinner (thicker) lines are the results for the DUNE (T2K) flux. (b) As for panel (a), except that the cumulative distributions are represented as a function of the absolute errors (in MeV). 025502-9 R. GONZÁLEZ-JIMÉNEZ et al. PHYSICAL REVIEW C 105, 025502 (2022) [1] R. Acciarri, C. Adams, J. Asaadi, B. Baller, T. Bolton, C. Bromberg, F. Cavanna, E. Church, D. Edmunds, A. Ereditato, S. Farooq, B. Fleming, H. Greenlee, G. Horton-Smith, C. James, E. Klein, K. Lang, P. Laurens, R. Mehdiyev, B. Page, O. Palamara, K. Partyka, G. Rameika, B. Rebel, M. Soderberg, J. Spitz,A.M.Szelc,M.Weber,T.Yang,andG.P.Zeller,Phys. Rev. D 90, 012008 (2014). [2] T. Cai, X.-G. Lu, L. A. Harewood, C. Wret, F. Akbar, D. A. Andrade, M. V. Ascencio, L. Bellantoni, A. Bercellie, M. Betancourt et al. (MINERνA Collaboration), Phys. Rev. D 101, 092001 (2020). [3] X.-G. Lu, M. Betancourt, T. Walton, F. Akbar, L. Aliaga, O. Altinok, D. A. Andrade, M. Ascencio, L. Bellantoni et al. (MINERvA Collaboration), Phys. Rev. Lett. 121, 022504 (2018). [4] T. Walton, M. Betancourt, L. Aliaga, O. Altinok, A. Bodek, A. Bravar, H. Budd, M. J. Bustamante, A. Butkevich et al. (MINERvA Collaboration), Phys. Rev. D 91, 071301 (2015). [5] K. Abe, J. Amey, C. Andreopoulos, L. Anthony, M. Antonova, S. Aoki, A. Ariga, Y. Ashida, Y. Azuma, S. Ban et al. (T2K Collaboration), Phys.Rev.D98, 032003 (2018). [6] R. Acciarri, M. A. Acero, M. Adamowski, C. Adams, P. Adamson, S. Adhikari, Z. Ahmad, C. H. Albright, T. Alion, E. Amador et al. (DUNE Collaboration), arXiv:1601.05471 [physics.ins-det]. [7] C. Simpson, K. Abe, C. Bronner, Y. Hayato, M. Ikeda, H. Ito, K. Iyogi, J. Kameda, Y. Kataoka, Y. Kato et al.,Astrophys. J. 885, 133 (2019). [8] U. Mosel, O. Lalakulich, and K. Gallmeister, Phys. Rev. Lett. 112, 151802 (2014). [9] A. P. Furmanski and J. T. Sobczyk, Phys. Rev. C 95, 065501 (2017). [10] J. W. Van Orden and T. W. Donnelly, Phys. Rev. C 100, 044620 (2019). [11] L. Munteanu, S. Suvorov, S. Dolan, D. Sgalaberna, S. Bolognesi, S. Manly, G. Yang, C. Giganti, K. Iwamoto, and C. Jesús-Valls, Phys.Rev.D101, 092003 (2020). [12]G.D.Megias,M.B.Barbaro,J.A.Caballero,J.E.Amaro, T. W. Donnelly, I. R. Simo, and J. W. V. Orden, J. Phys. G: Nucl. Part. Phys. 46, 015104 (2018). [13] J. M. Udías, P. Sarriguren, E. Moya de Guerra, E. Garrido, and J. A. Caballero, Phys. Rev. C 48, 2731 (1993). [14] J. M. Udías, P. Sarriguren, E. Moya de Guerra, E. Garrido, and J. A. Caballero, Phys. Rev. C 51, 3246 (1995). [15] M. C. Martínez, P. Lava, N. Jachowicz, J. Ryckebusch, K. Vantournhout, and J. M. Udías, Phys. Rev. C 73, 024607 (2006). [16] O. Benhar, A. Fabrocini, S. Fantoni, and I. Sick, Nucl. Phys. A 579, 493 (1994). [17] O. Benhar, N. Farina, H. Nakamura, M. Sakuda, and R. Seki, Phys.Rev.D72, 053005 (2005). [18] K. Amir-Azimi-Nili, J. Udias, H. Müther, L. Skouras, and A. Polls, Nucl. Phys. A 625, 633 (1997). [19] M. V. Ivanov, A. N. Antonov, G. D. Megias, J. A. Caballero, M. B. Barbaro, J. E. Amaro, I. Ruiz Simo, T. W. Donnelly, and J. M. Udías, Phys. Rev. C 99, 014610 (2019). [20] J. M. Franco-Patino, J. Gonzalez-Rosa, J. A. Caballero, and M. B. Barbaro, Phys.Rev.C102, 064626 (2020). [21] C. Colle, W. Cosyn, J. Ryckebusch, and M. Vanhalst, Phys. Rev. C89, 024603 (2014). [22] J. Ryckebusch, M. Vanhalst, and W. Cosyn, J. Phys. G: Nucl. Part. Phys. 42, 055104 (2015). [23] J. Caballero, T. Donnelly, E. Moya de Guerra, and J. Udías, Nucl. Phys. A 632, 323 (1998). [24] S. Gardner and J. Piekarewicz, Phys. Rev. C 50, 2822 (1994). [25] M. Ivanov, R. González-Jiménez, J. Caballero, M. Barbaro, T. Donnelly, and J. Udías, Phys. Lett. B 727, 265 (2013). [26] R. González-Jiménez, M. B. Barbaro, J. A. Caballero, T. W. Donnelly, N. Jachowicz, G. D. Megias, K. Niewczas, A. Nikolakopoulos, and J. M. Udías, Phys.Rev.C101, 015503 (2020). [27] E. D. Cooper, S. Hama, B. C. Clark, and R. L. Mercer, Phys. Rev. C 47, 297 (1993). [28] R. González-Jiménez, A. Nikolakopoulos, N. Jachowicz, and J. M. Udías, Phys. Rev. C 100, 045501 (2019). [29] J. M. Udías, J. A. Caballero, E. Moya de Guerra, J. R. Vignote, and A. Escuderos, Phys. Rev. C 64, 024614 (2001). [30] K. Abe, R. Akutsu, A. Ali, C. Alt, C. Andreopoulos, L. Anthony, M. Antonova, S. Aoki, A. Ariga, T. Arihara et al. (T2K Collaboration), Nature (London) 580, 339 (2020). [31] NuWro official repository, https://github.com/NuWro/nuwro. [32] A. Nikolakopoulos, N. Jachowicz, N. Van Dessel, K. Niewczas, R. González-Jiménez, J. M. Udías, and V. Pandey, Phys. Rev. Lett. 123, 052501 (2019). [33] D. B. Day, J. S. McCarthy, T. W. Donnelly, and I. Sick, Annu. Rev. Nucl. Part. Sci. 40, 357 (1990). [34] O. Moreno, T. W. Donnelly, J. W. Van Orden, and W. P. Ford, Phys. Rev. D 90, 013014 (2014). [35] J. E. Amaro, M. B. Barbaro, J. A. Caballero, R. GonzálezJiménez, G. D. Megias, and I. R. Simo, J. Phys. G: Nucl. Part. Phys. 47, 124001 (2020). [36] V. Pandey, N. Jachowicz, M. Martini, R. González-Jiménez, J. Ryckebusch, T. Van Cuyck, and N. Van Dessel, Phys.Rev.C 94, 054609 (2016). [37] J. Van Orden and T. Donnelly, Ann. Phys. 131, 451 (1981). [38] A. De Pace, M. Nardi, W. Alberico, T. Donnelly, and A. Molinari, Nucl. Phys. A 741, 249 (2004). [39] J. Amaro, M. Barbaro, J. Caballero, T. Donnelly, and A. Molinari, Phys. Rep. 368, 317 (2002). [40] M. Martini, M. Ericson, G. Chanfray, and J. Marteau, Phys. Rev. C 80, 065501 (2009). [41] J. Nieves, I. R. Simo, and M. J. Vicente Vacas, Phys. Rev. C 83, 045501 (2011). [42] T. Van Cuyck, N. Jachowicz, R. González-Jiménez, J. Ryckebusch, and N. Van Dessel, Phys.Rev.C95, 054611 (2017). [43] https://www.jlab.org/physics/hall-b/clas. 025502-16