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

TL;DR

This paper refines an intranuclear cascade model for pion-nucleus scattering, improving its accuracy and uncertainty estimates, which benefits neutrino experiments like T2K and Super-Kamiokande.

Contribution

It develops a method to tune Final State Interaction models using scattering data, addressing data deficiencies and enabling better neutrino interaction simulations.

Findings

The NEUT model is fitted to extensive pion-nucleus scattering data.

Uncertainty estimates for the model are robust and well-characterized.

The results are directly applicable to neutrino experiments using NEUT.

Abstract

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.

Tables

Table 1. Table 1: Description of the NEUT pion FSI parameters used in the fits described in this work. The overlap in the momentum regions is due to blending of the high and low energy models in NEUT as described in the text.

Parameter	Description	Momentum Region (MeV/$c$)
$f_{\mathrm{ABS}}$	Absorption
$f_{\mathrm{QE}}$	Quasi-elastic scatter
$f_{\mathrm{CX}}$	Single charge exchange
$f_{\mathrm{QEH}}$	Quasi-elastic scatter	$>400$
$f_{\mathrm{CXH}}$	Single charge exchange	$>400$
$f_{\mathrm{INEL}}$	Hadron (N+n$\pi$) production	$>400$

Table 2. Table 2: Summary of π ± superscript 𝜋 plus-or-minus \pi^{\pm} –A scattering data used for this tuning, including beam polarity, nuclear target type(s), momentum range and interaction channel(s). Some of these experiments may have measured data on other target nuclei that are not considered here.

Reference	Polarity	Targets	$p_{\pi }$ [MeV$/c$]	Channel(s)
B. W. Allardyce et al. Allardyce et al. (1973)	${\pi }^{\pm }$	12C, 27Al, 207Pb	710-2000	REAC
A. Saunders et al. Saunders et al. (1996)	${\pi }^{\pm }$	12C, 27Al	116-149	REAC
C. J. Gelderloos et al. Gelderloos et al. (2000)	${\pi }^{-}$	12C, 27Al, 63Cu, 207Pb	531-615	REAC
F. Binon et al. Binon et al. (1969)	${\pi }^{-}$	12C	219-395	REAC
O. Meirav et al. Meirav et al. (1989)	${\pi }^{\pm }$	12C, 16O	128-169	REAC
C. H. Q. Ingram Ingram et al. (1983)	${\pi }^{+}$	16O	211-353	QE
S. M. Levenson et al. Levenson et al. (1983)	${\pi }^{+}$	12C	194-416	QE
M. K. Jones et al. Jones et al. (1993)	${\pi }^{+}$	12C, 208Pb	363-624	QE, CX
D. Ashery et al. Ashery et al. (1981)	${\pi }^{\pm }$	12C, 27Al, 56Fe	175-432	QE, ABS+CX
H. Hilscher et al. Hilscher et al. (1970)	${\pi }^{-}$	12C	156	CX
T. J. Bowles Bowles et al. (1981)	${\pi }^{+}$	16O	128-194	CX
D. Ashery et al. Ashery et al. (1984)	${\pi }^{\pm }$	12C, 16O, 207Pb	265	CX
K. Nakai et al. Nakai et al. (1979)	${\pi }^{\pm }$	27Al, 63Cu	83-395	ABS
E. Bellotti et al. Bellotti et al. (1973a)	${\pi }^{+}$	12C	230	ABS
E. Bellotti et al. Bellotti et al. (1973b)	${\pi }^{+}$	12C	230	CX
I. Navon et al. Navon et al. (1983)	${\pi }^{+}$	12C, 56Fe	128	ABS+CX
R. H. Miller et al. Miller (1957)	${\pi }^{-}$	12C, 207Pb	254	ABS+CX
E. S. Pinzon Guerra et al. Pinzon Guerra et al. (2017)	${\pi }^{+}$	12C	206-295	ABS, CX

Table 3. Table 3: Minimum and maximum values for the NEUT FSI parameters allowed in the finite grid of pre-computed cross section predictions.

Parameter	Grid min.	Grid max.	Step size
$f_{\mathrm{QE}}$	0.1	1.7	0.1
$f_{\mathrm{ABS}}$	0.35	1.95	0.1
$f_{\mathrm{CX}}$	0.1	1.6	0.1
$f_{\mathrm{INEL}}$	0.2	2.6	0.2
$f_{\mathrm{QEH}}$	0.8	2.8	0.2
$f_{\mathrm{CXH}}$	1.2	2.4	0.2

Table 4. Table 4: Post-fit FSI parameters and the minimum χ 2 superscript 𝜒 2 \chi^{2} value with fixed normalization parameters using the TMultiDimFit and GNU-Octave interpolation methods.

Parameter	Best fit $\pm$ 1$\sigma$
	TMultiDimFit	GNU-Octave
$f_{\mathrm{QE}}$	1.07 $\pm$ 0.06	1.07 $\pm$ 0.04
$f_{\mathrm{ABS}}$	1.50 $\pm$ 0.08	1.40 $\pm$ 0.06
$f_{\mathrm{CX}}$	0.69 $\pm$ 0.05	0.70 $\pm$ 0.04
$f_{\mathrm{INEL}}$	0.89 $\pm$ 0.20	1.00 $\pm$ 0.15
$f_{\mathrm{QEH}}$	1.90 $\pm$ 0.25	1.82 $\pm$ 0.11
${\chi }^2$/n.d.o.f	150.74/59	149.03/59

Table 5. Table 5: Post-fit NEUT FSI parameters and the minimum χ 2 superscript 𝜒 2 \chi^{2} /n.d.o.f values obtained for fits with floating normalization parameters and the specified external data selections.

Parameter	Best fit $\pm$ 1$\sigma$
	Carbon-only	Light nuclei	All nuclei
$f_{\mathrm{QE}}$	1.07 $\pm$ 0.07	1.08 $\pm$ 0.07	1.08 $\pm$ 0.07
$f_{\mathrm{ABS}}$	1.24 $\pm$ 0.05	1.25 $\pm$ 0.05	1.26 $\pm$ 0.05
$f_{\mathrm{CX}}$	0.79 $\pm$ 0.05	0.80 $\pm$ 0.04	0.80 $\pm$ 0.04
$f_{\mathrm{INEL}}$	0.63 $\pm$ 0.27	0.71 $\pm$ 0.21	0.70 $\pm$ 0.20
$f_{\mathrm{QEH}}$	2.16 $\pm$ 0.34	2.14 $\pm$ 0.24	2.13 $\pm$ 0.22
${\chi }^2$/n.d.o.f	18.36/27	40.14/44	53.48/59

Table 6. Table 6: Post-fit NEUT FSI parameters and the minimum χ 2 superscript 𝜒 2 \chi^{2} /n.d.o.f values obtained for fits with floating normalization parameters and the five data sets with strongest pulls in Figure  8 removed.

Parameter	Best fit $\pm$ 1$\sigma$
	Carbon-only	Light nuclei	All nuclei
$f_{\mathrm{QE}}$	1.07 $\pm$ 0.08	1.09 $\pm$ 0.06	1.09 $\pm$ 0.06
$f_{\mathrm{ABS}}$	1.23 $\pm$ 0.05	1.25 $\pm$ 0.04	1.26 $\pm$ 0.04
$f_{\mathrm{CX}}$	0.79 $\pm$ 0.05	0.81 $\pm$ 0.04	0.80 $\pm$ 0.04
$f_{\mathrm{INEL}}$	0.63 $\pm$ 0.27	0.71 $\pm$ 0.21	0.71 $\pm$ 0.19
$f_{\mathrm{QEH}}$	2.15 $\pm$ 0.34	2.14 $\pm$ 0.24	2.13 $\pm$ 0.21
${\chi }^2$/n.d.o.f	17.95/26	35.55/42	45.13/54

Table 7. Table 7: Ψ Ψ \Psi for various scaling factors.

Scale	Mean	RMS
$1/1$	0.39	2.86
$1/9$	0.00	2.02
$1/25$	-0.11	1.28
$1/49$	-0.14	1.09
$1/81$	-0.09	1.00

Table 8. Table 8: Post-fit FSI parameters after error scaling

Parameter	Best fit $\pm$ 1$\sigma$
$f_{\mathrm{QE}}$	1.07 $\pm$ 0.31
$f_{\mathrm{ABS}}$	1.40 $\pm$ 0.43
$f_{\mathrm{CX}}$	0.70 $\pm$ 0.30
$f_{\mathrm{INEL}}$	1.00 $\pm$ 1.10
$f_{\mathrm{QEH}}$	1.82 $\pm$ 0.86

Equations

$$\frac{\rho (r)}{{\rho }_0}=\frac{1+w\frac{r^2}{c^2}}{1+\exp \left(\frac{r-c}{\alpha }\right)},$$

$$\begin{array}{r}{\chi }^2(f_{\mathrm{FSI}})=\sum _i^{\mathrm{Datasets}}\left(\sum _j^{n_i}\frac{1}{n_i}{\left(\frac{{\sigma }_j^{\mathrm{Data}}-{\lambda }_i^{-1}{\sigma }_j^{\mathrm{NEUT}}(f_{\mathrm{FSI}})}{\Delta {\sigma }_j^{\mathrm{Data}}}\right)}^2+{\left(\frac{{\lambda }_i-1}{{ϵ}_i}\right)}^2\right)\\ +\sum _{i,j}^{10}\left({\sigma }_i^{\mathrm{DUET}}-{\sigma }_i^{\mathrm{NEUT}}(f_{\mathrm{FSI}})\right)(V_{ij}^{-1}{)}^{\mathrm{DUET}}\left({\sigma }_j^{\mathrm{DUET}}-{\sigma }_j^{\mathrm{NEUT}}(f_{\mathrm{FSI}})\right),\end{array}$$

$$\Psi =\frac{{\sigma }_j^{\mathrm{Best}\mathrm{fit}}(f_{\mathrm{FSI}})-{\sigma }_j^{\mathrm{Data}}}{\Delta {\sigma }_j^{\mathrm{MC}}(\xi \Delta f_{\mathrm{FSI}})}.$$

Full text

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. Alvarez-Ruso et al. (2014); Garvey et al. (2015); Mosel (2016); Katori and Martini (2018); Alvarez-Ruso et al. (2018); Mahn et al. (2018); Betancourt et al. (2018)). 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 re-interaction 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$\pi$), such as T2K Abe et al. (2011a) and the planned Hyper-Kamiokande Hyp (2016). CC0$\pi$ events are dominated by charged-current quasi-elastic (CCQE) interactions $\nu_{l}+n\rightarrow l^{-}+p$ and $\bar{\nu}_{l}+p\rightarrow 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 Abe et al. (2015). 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$\pi$ 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 migrate events from one topology to another. Hadronic interactions inside (FSI) and outside (SI) the target nucleus are related, as discussed in Refs Dytman (2009); Golan et al. (2012); Salcedo et al. (1988); Oset and Salcedo (1987), 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 $\pi^{\pm}$–A scattering model used in NEUT Hayato (2009)—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 Böhlen et al. (2014); Asai (2006); Andreopoulos et al. (2015a); Golan et al. (2012). 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 Buss et al. (2012). The procedures used to evaluate the impact of these cascade model uncertainties are described elsewhere Abe et al. (2015); de Perio (2014) and will not be discussed here. The updates to the NEUT cascade model described have been incorporated in recent T2K analyses Abe et al. (2018).

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 $\pi^{\pm}$–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 Hayato (2009); de Perio (2014). 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:

[TABLE]

where $r$ is the distance from the center of the nucleus, $c$ is the nuclear radius, $\alpha$ is the surface thickness, and $w$ modifies the shape of the potential. The parameters $c,\alpha$ and $w$ were determined from an analysis of electron-scattering measurements Jager et al. (1974). 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_{\mathrm{N}}/100}$, where $R_{\mathrm{N}}$ is the size of the nucleus and is defined such that $\rho(R_{\mathrm{N}})/\rho_{0}\approx 10^{-4}$, which for carbon is $\sim$2.5 times the nuclear radius from Ref. Jager et al. (1974). The step spacing was chosen such that the probability of two or more interactions at each step was negligible. $R_{\mathrm{N}}$ acts to stop the cascade once the nuclear density, and thus the re-interaction 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_{\mathrm{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 mechanism.

For low momentum pions (defined in NEUT as $p_{\pi}<500$ MeV$/c$), tables computed from the Oset et al. model Salcedo et al. (1988) are used to determine quasi-elastic, single charge exchange, and absorption interaction probabilities inside the nucleus, and relate them to extra-nuclear $\pi^{\pm}$–A measurements. This model involves a many-body calculation in infinite nuclear matter with a local density approximation included. The $\pi^{\pm}$–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_{\pi}>500$ MeV$/c$), the interaction probabilities are calculated from $\pi^{\pm}$–A scattering data off free proton and deuterium compiled by the PDG Group (2010). The two models are blended linearly in the $400<p_{\pi}<500$ MeV$/c$ region to avoid discontinuities.

The most commonly used generators and simulation toolkits, including GENIE Andreopoulos et al. (2015a), NuWro Golan et al. (2012), FLUKA Böhlen et al. (2014), and Geant4 Asai (2006), have similar cascade models implemented, although the details vary. The notable exception is GiBUU Buss et al. (2012), 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 1, henceforth referred to simply as “FSI parameters” ($f_{\mathrm{FSI}}$). Each parameter scales the corresponding microscopic probability of a $\pi^{\pm}$ interaction at each step of the cascade, except for $f_{\mathrm{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 Abe et al. (2015); de Perio (2014). 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 Abe et al. (2011b, 2012). As such, this work includes external data taken with carbon, oxygen, aluminium, iron, copper, and lead targets.

We used data from $\pi^{+}$ and $\pi^{-}$ beams over a momentum 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 $\pi^{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 1, 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 $\pi^{\pm}$–A scattering.

Table 2 lists the external data, specifying the channels, targets, and $p_{\pi}$ range measured by each experiment. Most of the datasets used date from the 1950’s–1990’s, and contain limited information: in particular, covariances between datapoints in each channel, and between channels, are not available. A notable exception is the recent measurement of the ABS and CX cross sections by the DUET collaboration Pinzon Guerra et al. (2017). These measurements were partially motivated by the need to cover the lack of ABS measurements in the $\Delta$(1232) resonance region, which is of particular importance for T2K. Figure 3 shows the cross sections measured by DUET, and Figure 4 shows the covariance matrix provided between all DUET datapoints.

There are some changes to the datasets used, compared to earlier iterations of this work, described in Ref Abe et al. (2015); de Perio (2014). In addition to the new results from the DUET collaboration, Refs. Meirav et al. (1989); Hilscher et al. (1970) have been added. Measurements of exclusive pion production Grion et al. (1989); Rahav et al. (1991), double charge exchange Wood et al. (1992) and diffractive production Cronin et al. (1957) are not used, as more inclusive channels are preferred here. Ref. Rowntree et al. (1999) was taken on nuclear targets other than those identified in Table 2, and so is not included here. Some measurements of ABS Jones et al. (1993); Giannelli et al. (2000); Ransome et al. (1991) 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 Wilkin et al. (1973); Clough et al. (1974), 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 Takahashi et al. (1995); Carroll et al. (1976); Aoki et al. (2007) as it introduces model dependence. Finally, Ref. Crozon et al. (1964) has been neglected as the normalization of the results is inconsistent with other measurements, discussed in Ref. Allardyce et al. (1973).

IV Fit Strategy

The goal of the fit is two-fold:

Find the set of FSI parameters that best fit the external scattering data listed in Section III. 2. 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 cascade, predictions for the relevant $\pi^{\pm}$–A cross sections, $\sigma_{j}^{\mathrm{NEUT}}(f_{\mathrm{FSI}})$, as a function of the FSI parameters, $f_{\mathrm{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 3. 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 Hayato (2009); de Perio (2014).

The best fit is found by minimizing the $\chi^{2}$-statistic defined by

[TABLE]

where $n_{i}$ represents the number of data points in each data set except DUET, $\sigma_{j}^{\mathrm{Data}}$ and $\Delta\sigma_{j}^{\mathrm{Data}}$ are the external data set cross sections and their respective uncertainties, $\sigma_{i}^{\mathrm{DUET}}$ are the DUET data, and $(V_{ij})^{\mathrm{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 $\lambda_{i}$ were included as penalty terms for each data set (except DUET). The uncertainty for the normalization parameters ($\epsilon_{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 $\lambda_{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 $\chi^{2}$ in Equation 2 to the $f_{\mathrm{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_{\mathrm{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 James and Roos (1975). 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 $\chi^{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.

TMultiDimFit The TMultiDimFit tmu (2017). class of ROOT Bru (1997) was used to obtain a polynomial expression for the $\chi^{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 best-fit polynomial function to the finite grid reported a $\chi^{2}$/n.d.o.f of 0.29. 2. 2.

GNU-Octave n-dim splines The interpn function of GNU-Octave oct was used to obtain a multi-dimensional spline interpolation of the $\chi^{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_{\mathrm{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, $\lambda_{i}$ from Equation 2, were fixed to their nominal value of 1.0, with results shown in Section V.1. Secondly, the $\lambda_{i}$’s were treated as constrained parameters in the fit, as described in Section V.2. The former strategy forms the main result of this work, whereas the latter was used as a cross check.

V.1 Fixed Normalization Parameters

The best fit FSI parameters, while keeping the normalization parameters fixed, are presented in Table 4 for both interpolation methods. The spread in the results from the methods was covered by the uncertainties on the fitted parameters. The $\chi^{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 $\chi^{2}$ values from each method are in agreement and are shown in the last row of Table 4, 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 $\chi^{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.

V.2 Varying Normalization Parameters

To investigate the dependence of the $\chi^{2}$ and fit results on each data set, the normalization parameters ($\lambda_{i}$ in Equation 2) were allowed to vary. This keeps the n.d.o.f the same, but the $\chi^{2}$ is expected to reduce as correlated changes in the datasets from a single measurement introduce a smaller $\chi^{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 Wilkinson et al. (2016); Pumplin et al. (2001). 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.1 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 5 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 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$\sigma$ away from the nominal value of 1.0 in Figure 8. Table 6 shows the post-fit FSI parameters for that study. The effect on all parameters and the remaining normalizations was found to be negligible.

V.3 Drawing Error Envelopes

The constraints of the FSI parameters from the fit were used to form the 1$\sigma$ variations of the macroscopic scattering cross sections to allow for comparisons with external data. This was done as:

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 3 are discarded since the interpolation does not apply. 2. 2.

Draw the probability distribution function (PDF) for each momentum value for which $\sigma^{\mathrm{NEUT}}(f_{\mathrm{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. 3.

Fit a Gaussian function to the PDFs. 4. 4.

Use the mean and variance from the Gaussian for each value of momentum to draw the 1$\sigma$ error envelopes.

Figure 11 shows the resulting error bands for the $\pi^{+}$–12C 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 4, and was also true for the data on more nuclei. This is dealt with in Section V.4.

V.4 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. Wilkinson et al. (2016), which builds on an error inflation procedure defined in Ref. Pumplin et al. (2001). We inflate the fit uncertainties such that 68% of the data is indeed covered at 1$\sigma$, 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, $\Delta\sigma_{j}^{\mathrm{MC}}(f_{\mathrm{FSI}})$, by some scaling factor $\xi$, and calculate the level of agreement between the pion-nucleus cross section corresponding to the best fit FSI parameters, $\sigma_{j}^{\mathrm{MC}}(f_{\mathrm{FSI}})$, and the data, $\sigma_{j}^{\mathrm{Data}}$, given the inflated uncertainty on the pion-nucleus cross section, $\Delta\sigma_{j}^{\mathrm{MC}}(\xi\Delta f_{\mathrm{FSI}})$, using the figure of merit $\Psi$,

[TABLE]

By increasing the scaling factor $\xi$ until the distribution of $\Psi$ 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$\sigma$, 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 under-coverage. Table 7 shows the $\Psi$ distribution mean and RMS for various values of $\xi$. A linear fit to these RMS values was performed, and the scaling value for which the RMS was equal to 1, $\xi$ = 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 Pumplin et al. (2001).

The final post-fit values of the FSI parameters after applying scaling to the $\chi^{2}$ surface are shown in Table 8. 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 $\pi^{+}$–12C 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 Abe et al. (2015); de Perio (2014) dedicated to constraining the FSI parameters to data on T2K. Comparisons of the old and new results with all considered $\pi^{\pm}$–A combinations is found in Ref. Pinzon Guerra (2017). 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 5 and 8.

VI Model Comparisons

The modeling of $\pi^{\pm}$–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 Andreopoulos et al. (2015a); Dytman (2009). 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 A2/3 Dytman (2009). The hA model is the default in genie and has been used by most MINERvA and NOvA analyses published to date Aliaga et al. (2014); Aurisano et al. (2015).

• •

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 $\pi^{\pm}$–A data to reduce the need for A2/3 extrapolation Andreopoulos et al. (2015b). 2. 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 Andreopoulos et al. (2015a); Dytman (2009). 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 Perdue and Dytman (2018) to incorporate the Oset et al. model Salcedo et al. (1988) 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 Böhlen et al. (2014); Ferrari et al. (2005). Similar to NEUT, it uses the Oset et al. model Salcedo et al. (1988) to describe the absorptive width of the optical potential for pion momenta below 300 MeV$/c$. $\pi^{\pm}$–free nucleon cross sections are used to describe elastic, quasi-elastic, and charge exchange interactions.

• •

The Bertini cascade model of Geant 4.9.4 Wright and Kelsey (2015) 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 Flaminio et al. (1983).

• •

The nuwro event generator uses cascade model Golan et al. (2012) based on the Oset et al. model Salcedo et al. (1988), 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. 3. 3.

Transport models

• •

The Giessen Boltzmann-Uehling-Uhlenbeck (GiBUU) model is an implementation of transport model for nuclear reactions Buss et al. (2012). 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 $\pi^{\pm}-^{12}$C and 63Cu were taken directly from Ref. Buss et al. (2012). 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 2), 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 $\sim$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 ($\pi^{\pm}$-A$\rightarrow\pi^{+}\pi^{-}X$ where the $\pi^{+}$ is absorbed). The effective GENIE models also predict a much larger $\pi^{\pm}$ absorption cross section at higher momenta. Additionally, most of 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 Abe et al. (2015); de Perio (2014) by $\sim$50%, which should correspond to a similar reduction in the FSI and SI uncertainties, currently a large uncertainty in the T2K oscillation analysis Abe et al. (2015, 2017). 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 Ieki et al. (2015); Pinzon Guerra et al. (2017). DUET made measurements specifically targeted at T2K’s needs and provided 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 Cavanna et al. (2014); Abi et al. (2017) in charged particle test beams aiming to make such measurements on argon to fulfill this need for the future DUNE oscillation program Acciarri et al. (2015).

One limitation of this work is neglecting the measured pion kinematic distributions, for which only very limited data is available Levenson et al. (1983); Fujii et al. (2001). NEUT comparisons to these distributions can be found in Ref. de Perio (2014). Although this would require significant development of the NEUT cascade model, having such data from new $\pi^{\pm}$–A scattering experiments would provide a stringent constraint on models in the future.

Acknowledgements.

The authors would like to thank the members of the T2K collaboration and the authors of the NEUT generator for their help and support. We acknowledge the support of MEXT, Japan; NSERC (Grant No. SAPPJ-2014-00031), NRC and CFI, Canada; CEA and CNRS/IN2P3, France; SNSF and SERI, Switzerland; STFC, UK; and DOE, USA. In addition, participation of individual researchers and institutions has been further supported by funds from the Ontario Graduate Scholarship, the Queen 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 Loken et al. (2010). 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.

Appendix A Additional model comparisons

In Section VI, a comparison of the tuned NEUT results obtained in Section V, the data used in the fits described in this work (summarized in Table 2), and a variety of alternative intra-nuclear rescattering models were compared for $\pi^{\pm}$–carbon scattering (see Figures 13 and 14). Here, the similar plots are presented for: oxygen — Figures 15 and 16; aluminium — Figures 17 and 18; iron — Figures 19 and 20; copper — Figures 21 and 22; lead — Figures 23 and 24.

In addition, ratios of those plots with respect to the best fit NEUT predictions obtained in this work, are presented for reference for $\pi^{+}$–12C,16O,27Al ($\pi^{-}$-12C,16O,27Al) scattering in Figure 25 (Figure 26), and for $\pi^{+}$–56Fe,63Cu,207Pb ($\pi^{-}$-56Fe,63Cu,207Pb) scattering in Figure 25 (Figure 26).