Using world charged pion--nucleus scattering data to constrain an intranuclear cascade model

The NEUT intranuclear cascade model is described and fit to a large body of \pipm--nucleus scattering data. Methods are developed to deal with deficiencies in the available historical data, and robust uncertainty estimates are produced. The results are compared to a variety of simulation packages, and the data itself. This work provides a method for tuning Final State Interaction models, which are of particular interest to neutrino experiments that operate in the few-GeV energy region, and provides results which can be used directly by the T2K and Super-Kamiokande collaborations, for whom NEUT is the primary simulation package.

The NEUT intranuclear cascade model is described and fit to a large body of π ± -nucleus scattering data. Methods are developed to deal with deficiencies in the available historical data, and robust uncertainty estimates are produced. The results are compared to a variety of simulation packages, and the data itself. This work provides a method for tuning Final State Interaction models, which are of particular interest to neutrino experiments that operate in the few-GeV energy region, and provides results which can be used directly by the T2K and Super-Kamiokande collaborations, for whom NEUT is the primary simulation package.

I. INTRODUCTION
There is ongoing work dedicated to understanding or measuring neutrino interactions in the 0.1-10 GeV energy range aimed to reduce systematic uncertainties for current and future neutrino oscillation experiments (see, for example, recent reviews in Refs. [1][2][3][4][5][6][7]). The key problem is providing the best estimate of the incoming neutrino energy (as their oscillations are energy dependent) using measurable quantities in a detector: namely the outgoing particle content and kinematics. After a neutrino interacts with nucleon(s) within the nucleus, or quark(s) within a nucleon, the interaction products are propagated out of the nuclear medium before being measurable in a detector. Nucleons and pions have a considerable reinteraction probability when propagated through the nucleus, which can change the outgoing particle type, number and kinematics. These effects are known in the neutrino scattering community as Final State Interactions (FSI).
The effect of FSI is particularly important for oscillation analyses that use kinematic neutrino energy reconstruction, which selects neutrino events with a charged lepton and no pions in the final state (CC0π), such as T2K [8] and the planned Hyper-Kamiokande [9]. CC0π events are dominated by charged-current quasi-elastic (CCQE) interactions ν l + n → l − + p andν l + p → l + + n. For such events, one can reconstruct the incoming neutrino energy perfectly, provided that the incoming neutrino beam direction is known, and the initial state struck nucleon is at rest [10]. In a nuclear environment, the neutrino energy reconstruction is smeared by the initial state nucleon momentum, as they are not at rest inside the nucleus. Additionally, biases are introduced to the reconstructed neutrino energy by introducing events in the CC0π sample that are not CCQE, and so do not have the simple kinematic relationship to the neutrino energy. Events where a pion is produced at the primary neutrino interaction vertex but is later absorbed by the nucleus through FSI are a large contribution to this source of bias, thus FSI in carbon and oxygen are of great importance to the T2K oscillation analysis. Pions which interact outside the nucleus, known as secondary interactions (SI), can result in mis-reconstructing the event, and can 2 migrate events from one topology to another. Hadronic interactions inside (FSI) and outside (SI) the target nucleus are related, as discussed in Refs [11][12][13][14], although the relationship is not trivial. FSI and SI are among the nuclear effects that dominate systematic uncertainties for long-baseline neutrino oscillation analyses.
In this work, we describe the tuning of the π ± -A scattering model used in NEUT [15]-the primary interaction simulation package used by the T2K and Super-Kamiokande experiments-to the world dataset and describe its relation to the intranuclear semi-classical cascade models for FSI and SI. Although this work focuses on a single model, as it was motivated by the needs of the T2K experiment, the methods and the general conclusions developed herein are applicable to other cascade models [12,[16][17][18]. Comparisons are made between the results from this tuning and various other models, including the GiBUU simulation package-which goes beyond the semi-classical approximations made by cascade models [19]. The procedures used to evaluate the impact of these cascade model uncertainties are described elsewhere [10,20] and will not be discussed here. The updates to the NEUT cascade model described have been incorporated in recent T2K analyses [21].
Section II provides an overview of the NEUT cascade model and the parameters used in this tuning work. Section III provides a summary of the external scattering data sets. Section IV describes the fit strategy and methods. Section V presents the best fit values for the FSI parameters and their correlation matrix. Section VI provides an overview and set of comparisons to other models for π ± -A interactions. Finally, Section VII presents our conclusions and outlook.

II. NEUT CASCADE
Pion FSI are simulated in NEUT using a semi-classical intranuclear cascade model [15,20]. After a neutrino interaction occurs, the starting positions of the interaction products are chosen randomly from a nuclear density profile in the form of a three-parameter Fermi model (Woods-Saxon potential) described in Equation 1 and shown in Figure 1: where r is the distance from the center of the nucleus, c is the nuclear radius, α is the surface thickness, and w modifies the shape of the potential. The parameters c, α and w were determined from an analysis of electron-scattering measurements [22]. Note that the mass number (A) dependence of the model is encoded in these parameters.
In the case of oxygen, the two-parameter Fermi model is used (w = 0). The initial kinematic information of outgoing particles is taken directly from the neutrino interaction simulation and is fed into the cascade model. The initial interaction products are propagated "classically", in finite steps within the nuclear medium. The steps are in space and were chosen as dx = R N /100, where R N is the size of the nucleus and is defined such that ρ(R N )/ρ 0 ≈ 10 −4 , which for carbon is ∼2.5 times the nuclear radius from Ref. [22]. The step spacing was chosen such that the probability of two or more interactions at each step was negligible. R N acts to stop the cascade once the nuclear density, and thus the reinteraction probability, becomes negligible. The probabilities for each interaction type is calculated at each step and a random number generator is used to determine which, if any, interaction takes place. The cascade continues until the particle is absorbed or its position exceeds R N . Particles produced through inelastic FSI are also included into the cascade model from the relevant production point(s). The particles are propagated independently and are only subject to the potential from the nucleus, not from other nucleons and pions. Figure  2 shows an example cartoon of the intranuclear cascade 3 mechanism.
For low momentum pions (defined in NEUT as p π < 500 MeV/c), tables computed from the Oset et al. model [13] are used to determine quasi-elastic, single charge exchange, and absorption interaction probabilities inside the nucleus, and relate them to extra-nuclear π ± -A measurements. This model involves a many-body calculation in infinite nuclear matter with a local density approximation included. The π ± -A scattering is represented as a wave in a complex optical potential. Contributions from individual reaction channels are obtained by separating the real and complex parts of the potential and calculating the corresponding Feynman diagrams. For high momentum pions (p π > 500 MeV/c), the interaction probabilities are calculated from π ± -A scattering data off free proton and deuterium compiled by the PDG [23]. The two models are blended linearly in the 400 < p π < 500 MeV/c region to avoid discontinuities.
The most commonly used generators and simulation toolkits, including GENIE [18], NuWro [12], FLUKA [16], and Geant4 [17], have similar cascade models implemented, although the details vary. The notable exception is GiBUU [19], where the Boltzmann-Uehling-Uhlenbeck transport equations are solved for a more complete description of the nuclear medium. In GiBUU, the particles experience the potential from both the excited nuclear remnant and the other particles produced in the initial interaction. Comparisons between the results of this work, the data, and a variety of these alternate simulations are presented in Section VI.
The model is parameterized by the scaling factors summarized in Table I, henceforth referred to simply as "FSI parameters" (f FSI ). Each parameter scales the corresponding microscopic probability of a π ± interaction at each step of the cascade, except for f CX , which scales the charge exchange fraction of low momentum quasi-elastic (QE) scattering.
A reweighting scheme is used to allow the propagation of variations of these parameters to T2K analyses [10,20]. The scheme uses the information in the cascade for each event and re-calculates the interaction probabilities for each step with varied parameters to obtain a new value for the event probability. The FSI weight is defined as the ratio between the varied and nominal total event probabilities. However, the reweighting procedure is not exact due to the factorization into individual parameter variations, so was not used in this work.

III. SUMMARY OF EXTERNAL SCATTERING DATA
Although T2K analyses are predominantly interested in light nuclear targets (carbon and oxygen), there are many heavier targets in the T2K detectors [24,25]. As such, this work includes external data taken with carbon, oxygen, aluminium, iron, copper, and lead targets.
We used data from π + and π − beams over a momen- tum range from 60-2000 MeV/c. The interaction channels are defined exclusively from the number of pions in the final state, with any number of nucleons. This allows for direct comparisons between the external measurements of cross sections and the NEUT predicted cross sections. The following interaction channels were used: • Absorption (ABS): the incident pion is absorbed by the nucleus, and thus no pions are observed in the final state.
• Quasi-elastic Scattering (QE): only one pion is observed in the final state, and it has the same charge as the incident beam. The interaction is with a nucleon within the nucleus which differentiates QE from elastic scattering, where the interaction is with the nucleus as a whole. There is no requirement for the struck nucleon to be observed in the final state, as it may undergo FSI itself and not be observable.
• Single Charge Exchange (CX): the charged pion interacts with the nucleus and a single π 0 is observed in the final state.
• Absorption + Single Charge Exchange (ABS+CX): sum of ABS and CX (e.g., a charged pion is present in the initial state, but none are observable in the final state).
• Reactive (REAC): sum of ABS + CX + QE, double charge exchange, and hadron production. Double charge exchange events have a single pion in the final state which has the opposite charge to the incident beam. Hadron production events have final states with more than one pion.
These channels do not have a one-to-one correspondence to the parameters defined in Table I, but describe experimentally distinguishable interactions. Elastic scattering is experimentally defined as events with a very small scattering angle, where there is a very small momentum transfer to the target nucleus. Measurements of the elastic and total (reactive + elastic) cross sections are not used for this tuning as NEUT does not simulate elastic π ± -A scattering.
There are some changes to the datasets used, compared to earlier iterations of this work, described in Ref [10,20]. In addition to the new results from the DUET collaboration, Refs. [30,35] have been added. Measurements of exclusive pion production [44,45], double charge exchange [46] and diffractive production [47] are not used, as more inclusive channels are preferred here. Ref. [48] was taken on nuclear targets other than those identified in Table II, and so is not included here. Some measurements of ABS [33,49,50] were removed as those were performed with the goal of understanding multi-nucleon correlations and concentrated on final states with multiple protons instead of the more inclusive cross sections considered here. Measurements of the total cross section were not used [51,52], and measurements of reactive cross section where an optical model was used to separate the reactive and the elastic components from the total cross section were neglected [53][54][55] as it introduces model dependence. Finally, Ref. [56] has been neglected as the normalization of the results is inconsistent with other measurements, discussed in Ref. [26].

IV. FIT STRATEGY
The goal of the fit is two-fold: 1. Find the set of FSI parameters that best fit the external scattering data listed in Section III.
2. Set uncertainties for the FSI parameters that span the errors from the external data, and extract correlations between the parameters Given the computational cost of the NEUT cas-cade, predictions for the relevant π ± -A cross sections, σ NEUT j (f FSI ), as a function of the FSI parameters, f FSI , were pre-computed for a finite grid of FSI parameters. The FSI parameters are rescaling parameters of the nominal NEUT FSI prediction, where a parameter value of 1 is the nominal NEUT value, and are not fundamental physics parameters in their own right. The minimum and maximum values allowed for the FSI parameters and the step sizes are summarized in Table III. The predictions were only calculated for values of p π for which data was available to reduce the computational load. The predictions were calculated by running the NEUT cascade with a mono-energetic pion beam using the piscat simulation of NEUT [15,20]. The best fit is found by minimizing the χ 2 -statistic defined by where n i represents the number of data points in each data set except DUET, σ Data j and ∆σ Data j are the external data set cross sections and their respective uncertainties, σ DUET i are the DUET data, and (V ij ) DUET is the DUET covariance matrix shown in Figure 4. Note that every channel measured by the uncorrelated datasets (n i ) is treated as a single, uncorrelated point, in the fit. The normalization parameters λ i were included as penalty terms for each data set (except DUET). The uncertainty for the normalization parameters ( i ) was assigned to be 40% following the representative correlations in the DUET data sets as seen in Figure 4, and was an ad-hoc choice made for this analysis. These λ i normalization parameters include a fully correlated component between the datapoints of each dataset, to allow the effect of no correlations to be investigated. Note that they are not used as fit parameters for all fits described in this work.
The dependence of the χ 2 in Equation 2 to the f CXH parameter was found to be very weak, so it was decided to fix this parameter to its nominal value for the fit. An additional parameter designed to uniformly scale all the microscopic interaction probabilities (f ALL ) was also investigated, but similarly, it was found to have little impact on the fit and was consequently not included.
The minimization was performed using the MIGRAD algorithm of the MINUIT package [57]. The advantage of using this algorithm is that it provides both best fit parameters and their correlations. The difficulty is that the algorithm requires a χ 2 surface with smooth first and second order derivatives. Two different interpolation methods were used to smooth the finite pre-computed grid, to determine if biases were introduced by either.
1. TMultiDimFit The TMultiDimFit [58]. class of ROOT [59] was used to obtain a polynomial expression for the χ 2 grid in terms of the FSI parameters. The best-fit polynomial function contained up to fourth-degree polynomials with 53 terms in total, including cross-terms. A comparison of the bestfit polynomial function to the finite grid reported a χ 2 /n.d.o.f of 0.29.
2. GNU-Octave n-dim splines The interpn function of GNU-Octave [60] was used to obtain a multidimensional spline interpolation of the χ 2 grid. Cubic splines are evaluated around the requested point.
In general, it is difficult to compare multi-dimensional distributions and the two interpolation methods. Their one-dimensional projections of each is shown for illustrative purposes in Figure 5. The f CXH parameter caused problems with the convergence of the TMultiDimFit parameterization due to its low constraining power (there is very little available relevant data), and thus was not included for this comparison. The interpolation methods are found to be consistent and no significant biases are expected.

V. FIT RESULTS
The fit was carried out in two ways. Firstly, the normalization parameters for each dataset, λ i from Equation 2, were fixed to their nominal value of 1.0, with results shown in Section V A. Secondly, the λ i 's were treated as constrained parameters in the fit, as described in Section V B. The former strategy forms the main result of this work, whereas the latter was used as a cross check.

A. Fixed Normalization Parameters
The best fit FSI parameters, while keeping the normalization parameters fixed, are presented in Table IV for both interpolation methods. The spread in the results from the methods was covered by the uncertainties on the fitted parameters. The χ 2 between the two methods, using the covariance matrix calculated for the GNU-Octave fit, was 4.14 for 5 degrees of freedom. The minimum χ 2 values from each method are in agreement and are shown in the last row of Table IV, along with the number of degrees of freedom in the fit. This confirms that the interpolation methods are not introducing biases, and indicates that the χ 2 did not have a local minimum that would affect the minimization process. Thus, no additional uncertainties due to the interpolation method choice were deemed necessary. Figures 6 and 7 show the covariance matrices obtained from MINUIT using each interpolation method. Stronger correlations across the FSI parameters are observed when using the TMultiDimFit interpolation. This can be understood as the effect of the polynomial parameterization, which inherently carries strong correlations from the large number of cross-terms. For this reason, it was decided to use the GNU-Octave interpolation for the final results.

B. Varying Normalization Parameters
To investigate the dependence of the χ 2 and fit results on each data set, the normalization parameters (λ i in Equation 2) were allowed to vary. This keeps the n.d.o.f the same, but the χ 2 is expected to reduce as correlated changes in the datasets from a single measurement introduce a smaller χ 2 penalty (e.g., if there is a flux modeling issue with a dataset). This approach has been used elsewhere when fitting to datasets with missing covariance information [61,62]. Including these parameters as pull terms, rather than adding correlations between datapoints in a covariance matrix, gives greater insight into the behaviour of the fitter. For instance, a dataset with datapoints that all pull strongly on the fitted FSI parameters in the same direction in Section V A would have a normalization parameter largely deviating from 1.0.
The fit was performed for three selections of external data sets: data on carbon nuclei only; data on light nuclei (carbon, oxygen, aluminum); and data on light and heavy nuclei (carbon, oxygen, aluminum, iron, copper, and lead). This was done to investigate how the cascade model scales with increasing nuclear size, and the effect on the FSI parameters. Table V shows the best fit FSI parameters for each case using the GNU-Octave method. The agreement across the three cases indicates that the model is able to consistently describe all the data, and that the fit is well behaved. Figure 8 shows the fitted normalization parameters for the three selections of external data. The normalization parameters roughly follow a Gaussian distribution, and the results from each of the three cases essentially overlap each other. The varying sizes of the post-fit errors is understood to be a consequence of assigning the same 40% pre-fit error to all the parameters, as different channels have been experimentally measured with different levels  f of precision. For instance, inelastic and reactive processes tend to have much smaller uncertainties and the model does much better there. As a further cross-check, the fit was repeated after removing the five data sets whose post-fit normalization parameters were more than 1σ away from the nominal  value of 1.0 in Figure 8. Table VI shows the post-fit FSI parameters for that study. The effect on all parameters and the remaining normalizations was found to be negligible.

C. Drawing Error Envelopes
The constraints of the FSI parameters from the fit were used to form the 1σ variations of the macroscopic scattering cross sections to allow for comparisons with external data. This was done as: 1. Generate a random correlated throw from the obtained covariance matrix (e.g., Figure 7) through Cholesky decomposition. Throws outside the finite grid defined in Table III are discarded since the interpolation does not apply.
2. Draw the probability distribution function (PDF) for each momentum value for which σ NEUT (f FSI ) has been calculated using the GNU-Octave spline interpolation. Figures 9 and 10 show examples of these distributions for two combinations of interaction channel and momentum.
3. Fit a Gaussian function to the PDFs.
4. Use the mean and variance from the Gaussian for each value of momentum to draw the 1σ error envelopes. Figure 11 shows the resulting error bands for the π + -12 C cross sections, obtained using the constraints from the correlation matrix in Figure 7 and the procedure described above. The coverage of the error is insufficient when using the constraints from Table IV, and was also true for the data on more nuclei. This is dealt with in Section V D.

D. Error Inflation
The problem of lack of coverage arises from fitting to datasets for which no covariances are provided. Here we follow the same procedure as described in Ref. [61], which builds on an error inflation procedure defined in Ref. [62]. We inflate the fit uncertainties such that 68% of the data is indeed covered at 1σ, as one would expect from the fit, whilst keeping the post-fit central values and correlations of the FSI parameters. Explicitly, we scale the cross section uncertainties obtained from the fit, ∆σ MC j (f FSI ), by some scaling factor ξ, and calculate the level of agreement between the pion-nucleus cross section corresponding to the best fit FSI parameters, σ MC j (f FSI ), and the data, σ Data j , given the inflated uncertainty on  11. Comparison of the available π + -12 C cross section external data with the NEUT best fit (solid black line) and 1σ band (red), and the best fit (dashed black line) and 1σ band (blue) obtained from a previous iteration of this fit [10,20].
the pion-nucleus cross section, ∆σ MC j (ξ∆f FSI ), using the figure of merit Ψ, By increasing the scaling factor ξ until the distribution of Ψ plotted for all data points has an RMS of 1, we enforce the expected coverage in a naive way-e.g. 68% of the data are covered by the post-fit uncertainties at 1σ, neglecting correlations between data points, and between different points in the MC. This approximation almost certainly results in over coverage, and very conservative inflated uncertainty bands, but given the lack of information available in the fit, rigorous definitions of coverage are not possible and we prefer over-to undercoverage. Table VII shows the Ψ distribution mean and RMS for various values of ξ. A linear fit to these RMS values was performed, and the scaling value for which the RMS was equal to 1, ξ = 57.0, was found. Such a large scaling factor is not very surprising given the large number of datasets without correlations and is consistent with the results found elsewhere [62].
The final post-fit values of the FSI parameters after applying scaling to the χ 2 surface are shown in Table VIII. As the correlation among parameters is not affected by the scaling procedure the correlation matrix presented in Figure 7 remains the same. The scaled uncertainty bands for π + -12 C scattering is shown in Figure 12 and  can be compared with the unscaled uncertainty bands shown in Figure 11. The results are compared to earlier work [10,20] dedicated to constraining the FSI parameters to data on T2K. Comparisons of the old and new results with all considered π ± -A combinations is found in Ref. [63]. We note that after scaling of the uncertainties, there is good agreement found between the fit results with fixed and varied normalization parameters, shown in Tables V and VIII The modeling of π ± -A interactions is fundamental for a complete description of neutrino-nucleus interactions. All complete neutrino event generators include a model for these Final State Interactions. In this section, a variety of available models are compared to the tuned NEUT model with scaled uncertainties, presented in this work, and the data. The models can be divided into the following three categories:

Effective Models
• genie hA is a simple, data-driven, effective model [11,18]. Instead of a full cascade model, the total cross section for each scattering process within the nucleus is calculated-effectively reducing the cascade to a single step. The model is normalized to pion/nucleon-iron scattering data and cross sections for targets other than iron are obtained by scaling by A 2/3 [11]. The hA model is the default in genie and has been used by most MINERvA and NOvA analyses published to date [64,65].
• The genie hA 2014 model is a development version of the hA model and is included here for completeness. It includes a wider range of π ± -A data to reduce the need for A 2/3 extrapolation [66].

2.
Cascade models all follow the same general principles as the NEUT cascade model described in Section II. The implementation, input data, and theoretical approximations often differ to reflect the motivations and priorities of the simulation packages.
• genie hN is an alternative full cascade model for FSI available in genie [11,18]. Only data on free nucleons is used as an input. The development version of this model, genie hN2015, is used for these comparisons. Work is ongoing [67] to incorporate the Oset et al. model [13] as is used in NEUT and NuWro. • The PEANUT (Pre-Equilibrium Approach to NUclear Thermalization) model is an intra-nuclear cascade model implemented in fluka [16,68]. Similar to NEUT, it uses the Oset et al. model [13] to describe the absorptive width of the optical potential for pion momenta below 300 MeV/c. π ± -free nucleon cross sections are used to describe elastic, quasi-elastic, and charge exchange interactions. • The Bertini cascade model of Geant 4.9.4 [69] is part of the QGSP_BERT physics list. It is valid for pion momenta below 9.9 GeV/c.
It also handles all other long-live hadrons. A detailed treatment of pre-equilibrium and evaporation physics is included, relevant at momenta below 200 MeV where the cascade model approach begins to fail, as the de Broglie wavelength of the probe is roughly the same as the distance between nucleons in the target nucleus. The CERN-HERA compilation of hadron-nucleus elementary cross section data is used as input [70].
• The nuwro event generator uses cascade model [12] based on the Oset et al. model [13], as in NEUT. It introduces a phenomenological treatment of the formation zone (or time) effect, which can interplay with the pion absorption probability in a non-trivial momentum dependent way.

Transport models
• The Giessen Boltzmann-Uehling-Uhlenbeck (GiBUU) model is an implementation of transport model for nuclear reactions [19]. It describes the dynamical evolution of the interacting nuclear system through a coupled set of semi-classical kinetic equations, while taking into account the hadronic potentials and the equation of state of nuclear matter within the Boltzmann-Uehling-Uhlenbeck (BUU) theory.
Thin target particle gun simulations were used to produce microscopic cross sections for each of the models described above, and used to compare to data here, with the exception of GiBUU, where the the predictions for REAC and ABS of π ± − 12 C and 63 Cu were taken directly from Ref. [19]. To allow for a consistent comparison, the interactions channels were defined using only the final state particles, as described in Section II.
Comparisons between the tuned NEUT results obtained in Section V, the data used in the tuning (summarized in Table II), and the alternative models described above, are shown for carbon in Figures 13 and 14. Similar plots for oxygen, aluminium, iron, copper and lead are found in Appendix A, including ratios of the models to data.
The NEUT prediction with scaled error bands covers ∼2/3 of the available data points as intended. In the energy regions where there is data, most of the models are in reasonable agreement. However, there are noticeable differences between the models outside the data range.
Most notably the NEUT model predicts a much larger cross section for the single charge exchange channel at higher momenta. Given that the agreement is recovered in the reactive channel, this is likely to be caused by differences in the tagging of hadron production events (π ± -A→ π + π − X where the π + is absorbed). The effective GENIE models also predict a much larger π ± absorption cross section at higher momenta. Additionally, most of  FIG. 13. Comparison of the available π + -12 C cross section external data with the NEUT best fit and its 1σ error band obtained in this work, and other models.  FIG. 14. Comparison of the available π − -12 C cross section external data with the NEUT best fit and its 1σ error band obtained in this work, and other models.
the models fail to properly reproduce the data for the absorption cross section at low momenta for heavy nuclei (see Appendix A). One possible explanation is the lack of a model for the Coulomb attraction felt by the negatively charge pion, which is included for the GiBUU model, the only model which obtains good agreement with the data in this region.

VII. SUMMARY AND OUTLOOK
The cascade model used to simulate pion interactions in NEUT has been tuned to a variety of available data on several nuclear targets, and uncertainties have been produced. These have recently been used in analyses by the T2K and Super-Kamiokande collaborations, for whom NEUT is the primary neutrino simulation software. The uncertainties produced in this work are reduced with respect to previous attempts to constrain the pion interaction parameters [10,20] by ∼50%, which should correspond to a similar reduction in the FSI and SI uncertainties, currently a large uncertainty in the T2K oscillation analysis [10,71]. Additionally, the methods developed can be applied for tuning other similar cascade simulation packages.
One of the additions to this fit over previous attempts is the new DUET data [43,72]. DUET made measurements specifically targeted at T2K's needs and pro-vided more information than available from most historical measurements from the 1950's-1990's. To reduce the uncertainties further for future experiments more modern data is needed. Fortunately, there are dedicated experiments [73,74] in charged particle test beams aiming to make such measurements on argon to fulfill this need for the future DUNE oscillation program [75].
One limitation of this work is neglecting the measured pion kinematic distributions, for which only very limited data is available [32,76]. NEUT comparisons to these distributions can be found in Ref. [20]. Although this would require significant development of the NEUT cascade model, having such data from new π ± -A scattering experiments would provide a stringent constraint on models in the future.
Elizabeth II Graduate Scholarship in Science & Technology, the Japan Society for the Promotion of Science (JSPS) Postdoctoral Fellowship for Research in Japan, the H2020 Grant No. RISE-GA644294-JENNIFER, the Alfred P. Sloan Foundation and the DOE Early Career program, USA. Computations were performed on the GPC supercomputer at the SciNet HPC Consortium [77]. SciNet is funded by: the Canada Foundation for Innovation under the auspices of Compute Canada; the Government of Ontario; Ontario Research Fund -Research Excellence; and the University of Toronto. Ratios of the π + -12 C, 16 O, 27 Al data and the GEANT4 (black), GENIE (red), NuWro (blue), and GiBUU (magenta) predictions to the NEUT best fit for the five interaction channels used in the fit. The green band represents the ±1σ band.     63 Cu, 207 Pb data and the GEANT4 (black), GENIE (red), NuWro (blue), and GiBUU (magenta) predictions to the NEUT best fit. The green band represents the ±1σ band.