Tuning the GENIE Pion Production Model with MINERvA Data

TL;DR

This paper improves the GENIE neutrino interaction model by tuning it with MINERvA data and other experiments, addressing discrepancies in pion production predictions at few-GeV energies.

Contribution

It introduces specific parameter tunings and corrections to enhance GENIE's agreement with experimental data, including a new correction for low-Q^2 effects.

Findings

Tuned GENIE model better matches MINERvA data

Identified discrepancies between nuclear and nucleon target data

Developed a framework for future model extensions

Abstract

Faced with unresolved tensions between neutrino interaction measurements at few-GeV neutrino energies, current experiments are forced to accept large systematic uncertainties to cover discrepancies between their data and model predictions. In this paper, the widely used pion production model in GENIE is compared to four MINERvA charged current pion production measurements using NUISANCE. Tunings, ie, adjustments of model parameters, to help match GENIE to MINERvA and older bubble chamber data are presented here. We find that scattering off nuclear targets as measured in MINERvA is not in good agreement with scattering off nucleon (hydrogen or deuterium) targets in the bubble chamber data. An additional ad hoc correction for the low-$Q^2$ region, where collective effects are expected to be large, is also presented. While these tunings and corrections improve the agreement of GENIE with the data, the modeling is imperfect. The development of these tunings within the NUISANCE frameworkallows for straightforward extensions to other neutrino event generators and models, and allows omitting and including new data sets as they become available.

Figures

Tables

Table 1. Table 1: Summary of the measurements used in this analysis. W rec subscript 𝑊 rec W_{\mathrm{rec}}  is the true reconstructed hadronic mass assuming the struck nucleon is at rest. None of the measurements veto on activity other than the μ 𝜇 \mu and π 𝜋 \pi in their signal definition, and all selections require 1.5 < E ν < 20 ​ GeV 1.5 subscript 𝐸 𝜈 20 GeV 1.5<E_{\nu}<20~{}\mbox{GeV} .

Channel	${\nu }_{\mu }$CC1${\pi }^{\pm }$ Eberly et al. (2015)	${\nu }_{\mu }$CC$N{\pi }^{\pm }$ McGivern et al. (2016)	${\nu }_{\mu }$CC$1{\pi }^0$ Altinok et al. (2017)	${\overline{\nu }}_{\mu }$CC$1{\pi }^0$ McGivern et al. (2016)
$N_{\mathrm{bins}}$ $p_{\mu }$	8	9	8	9
$N_{\mathrm{bins}}$ ${\theta }_{\mu }$	9	9	9	9
$N_{\mathrm{bins}}$ $T_{\pi }$	7	7	7	7
$N_{\mathrm{bins}}$ ${\theta }_{\pi }$	14	14	11	11
$N_{\mathrm{bins}}$ total	38	39	35	36
Signal definition	1${\pi }^{\pm }$, $\ge$ 0${\pi }^0$	$>0$${\pi }^{\pm }$, $\ge$ 0${\pi }^0$	1${\pi }^0$, 0${\pi }^{\pm }$	1${\pi }^0$, 0${\pi }^{\pm }$
	1${\mu }^{-}$	1${\mu }^{-}$	1${\mu }^{-}$	1${\mu }^{+}$
	$W_{\mathrm{rec}}$	$W_{\mathrm{rec}}$	$W_{\mathrm{rec}}$	$W_{\mathrm{rec}}$
	—	—		—

Table 2. Table 2: Channel by channel contributions to the χ 2 superscript 𝜒 2 \chi^{2} at different stages of the tuning process.

Distribution	Channel	$N_{\mathrm{bins}}$	Default	ANL/BNL	FrAbs Tune	FrInel Tune
$p_{\mu }$ (Rate)	${\nu }_{\mu }$CC1${\pi }^{\pm }$	8	19.1	13.8	12.0	12.3
	${\nu }_{\mu }$CC$N{\pi }^{\pm }$	9	35.4	19.5	26.1	26.8
	${\nu }_{\mu }$CC$1{\pi }^0$	8	11.1	19.6	19.0	19.3
	${\overline{\nu }}_{\mu }$CC$1{\pi }^0$	9	7.4	6.4	6.2	6.3
${\theta }_{\mu }$ (Shape)	${\nu }_{\mu }$CC1${\pi }^{\pm }$	9	7.1	12.4	7.5	7.4
	${\nu }_{\mu }$CC$N{\pi }^{\pm }$	9	4.5	10.4	4.0	4.1
	${\nu }_{\mu }$CC$1{\pi }^0$	9	35.1	71.5	44.5	45.6
	${\overline{\nu }}_{\mu }$CC$1{\pi }^0$	9	9.3	14.0	10.2	10.3
$T_{\pi }$ (Shape)	${\nu }_{\mu }$CC1${\pi }^{\pm }$	7	2.9	2.6	2.5	2.3
	${\nu }_{\mu }$CC$N{\pi }^{\pm }$	7	39.8	34.7	31.2	29.4
	${\nu }_{\mu }$CC$1{\pi }^0$	7	28.3	31.4	30.9	29.9
	${\overline{\nu }}_{\mu }$CC$1{\pi }^0$	7	19.3	17.9	16.6	16.0
${\theta }_{\pi }$ (Shape)	${\nu }_{\mu }$CC1${\pi }^{\pm }$	14	25.4	26.5	13.0	12.6
	${\nu }_{\mu }$CC$N{\pi }^{\pm }$	14	11.7	11.1	6.9	6.2
	${\nu }_{\mu }$CC$1{\pi }^0$	11	13.5	15.0	8.3	8.9
	${\overline{\nu }}_{\mu }$CC$1{\pi }^0$	11	5.7	5.9	3.4	3.5
Total ${\chi }^2$		148	275.6	312.7	242.3	240.7

Table 3. Table 3: Summary of the GENIE dials optimized in this note, their default values and the uncertainties recommended by the GENIE collaboration. We do not use the defaults for M A res superscript subscript 𝑀 A res M_{\mathrm{A}}^{\mathrm{res}} , NormRes NormRes \mathrm{NormRes} and NonRes 1 ​ π 1 𝜋 1\pi and instead impose central values and uncertainties from tunings to ANL and BNL data as described in the text.

Parameter	Default Value	GENIE parameter name
CC Resonant Axial Mass ($M_A^{\mathrm{res}}$)	$1.12\pm 0.22\text{GeV}$	MaCCRES
CC Resonant Normalization ($\mathrm{NormRes}$)	$100\pm 20\text{\%}$	NormCCRES
CC1$\pi$ Nonresonant Normalization (NonRes$1\pi$)	$100\pm 50\text{\%}$	NonRESBGvnCC1pi
		NonRESBGvpCC1pi
		NonRESBGvbarnCC1pi
		NonRESBGvbarpCC1pi
CC2$\pi$ Nonresonant Normalization (NonRes$2\pi$)	$100\pm 50\text{\%}$	NonRESBGvnCC2pi
		NonRESBGvpCC1pi
		NonRESBGvbarnCC1pi
		NonRESBGvbarpCC1pi
Pion Angular Emission ($\pi$-iso)	0 (RS)	Theta_Delta2Npi
Pion Absorption FSI Fraction (FrAbs)	$100\pm 30\text{\%}$	FrAbs_pi
Pion Inelastic FSI Fraction (FrInel)	$100\pm 40\text{\%}$	FrInel_pi

Table 4. Table 4: Prefit and best fit central values and uncertainties from tuning GENIE to the ANL/BNL pion production measurements. The prefit uncertainties are those recommended by the GENIE collaboration. The tuned values are used as penalty terms with supplied covariance matrix of Fig  5 .

Parameter	GENIE default	ANL/BNL tune
$M_A^{\mathrm{res}}$ [GeV]	$1.12\pm 0.22$	$0.94\pm 0.05$
$\mathrm{NormRes}$ [%]	$100\pm 20$	$115\pm 7$
NonRes$1\pi$ [%]	$100\pm 50$	$43\pm 4$

Table 5. Table 5: Fit results from tuning GENIE parameters in NUISANCE. The “ANL/BNL Value” column shows the contributions when parameters are fixed at values of Table  4 .

Parameter	Default Value	ANL/BNL Value	FrAbs Fit Result	FrInel Result
$M_A^{\mathrm{res}}$ (GeV)	$1.12\pm 0.22$	$0.94\pm 0.05$	$1.07\pm 0.04$	$1.08\pm 0.04$
$\mathrm{NormRes}$ (%)	$100\pm 20$	$115\pm 7$	$94\pm 6$	$92\pm 6$
NonRes$1\pi$ (%)	$100\pm 50$	$43\pm 4$	$44\pm 4$	$44\pm 4$
NonRes$2\pi$ (%)	$100\pm 50$	-	$166\pm 32$	$161\pm 33$
$\pi$-iso	0 = RS	-	1 = Iso (limit)	1 = Iso (limit)
FrAbs (%)	100 $\pm$ 30	-	$109\pm 16$	-
FrInel (%)	100 $\pm$ 40	-	-	$109\pm 24$
MINERvA ${\chi }^2$	275.6	312.7	242.3	240.7
${\chi }_{\text{pen}}^2$	299.3	0.0	9.3	11.1
Total ${\chi }^2$	574.8	312.7	251.6	251.8
$N_{\mathrm{DoF}}$	148	148	145	145

Table 6. Table 6: Individual channel tuning results when the FrAbs dial is treated as the free FSI parameter.

Parameter	${\nu }_{\mu }$CC1${\pi }^{\pm }$	${\nu }_{\mu }$CC$N{\pi }^{\pm }$	${\nu }_{\mu }$CC$1{\pi }^0$	${\overline{\nu }}_{\mu }$CC$1{\pi }^0$
$M_A^{\mathrm{res}}$ (GeV)	0.97 $\pm$ 0.05	$0.97\pm 0.05$	$1.02\pm 0.05$	$0.96\pm 0.05$
$\mathrm{NormRes}$ (%)	110 $\pm$ 7	$110\pm 7$	$104\pm 7$	$111\pm 7$
NonRes$1\pi$ (%)	43 $\pm$ 4	$42\pm 4$	$44\pm 4$	$43\pm 4$
NonRes$2\pi$ (%)	$300$ (limit)	99$\pm$30	$300$ (limit)	$300$ (limit)
$\pi$-iso	1 = Iso (limit)	1 = Iso (limit)	1 = Iso (limit)	1 = Iso (limit)
FrAbs (%)	156 $\pm$ 53	$128\pm 34$	$126\pm 17$	$82\pm 31$
MINERvA ${\chi }^2$	36.6	64.1	92.3	34.6
${\chi }_{\text{pen}}^2$	0.5	0.7	3.2	0.3
Total ${\chi }^2$	37.1	64.8	95.5	34.9
$N_{\mathrm{DoF}}$	35	36	32	33

Table 7. Table 7: Individual channel tuning results when the FrInel dial is treated as the free FSI parameter.

Parameter	${\nu }_{\mu }$CC1${\pi }^{\pm }$	${\nu }_{\mu }$CC$N{\pi }^{\pm }$	${\nu }_{\mu }$CC$1{\pi }^0$	${\overline{\nu }}_{\mu }$CC$1{\pi }^0$
$M_A^{\mathrm{res}}$ (GeV)	0.97 $\pm$ 0.05	$0.97\pm 0.05$	$1.03\pm 0.05$	$0.96\pm 0.05$
$\mathrm{NormRes}$ (%)	109 $\pm$ 7	$108\pm 7$	$103\pm 7$	$112\pm 7$
NonRes$1\pi$ (%)	42 $\pm$ 4	$42\pm 4$	$43\pm 4$	$43\pm 4$
NonRes$2\pi$ (%)	$300$ (limit)	110$\pm$30	$300$ (limit)	$300$ (limit)
$\pi$-iso	1 = Iso (limit)	1 = Iso (limit)	1 = Iso (limit)	1 = Iso (limit)
FrInel (%)	117 $\pm$ 54	$127\pm 33$	0 (limit)	$80\pm 59$
MINERvA ${\chi }^2$	37.1	63.4	86.9	34.9
${\chi }_{\text{pen}}^2$	0.7	1.3	3.4	0.2
Total ${\chi }^2$	37.8	64.7	90.3	35.1
$N_{\mathrm{DoF}}$	35	36	32	33

Table 8. Table 8: Ad hoc low- Q 2 superscript 𝑄 2 Q^{2} suppression model tuning results compared to the tuning results without the low- Q 2 superscript 𝑄 2 Q^{2} suppression.

Parameter	FrAbs Tune	FrAbs + low-$Q^2$ Tune	FrInel Tune	FrInel + low-$Q^2$ Tune
$M_A^{\mathrm{res}}$ $\left(GeV\right)$	$1.07\pm 0.04$	$0.92\pm 0.02$	$1.08\pm 0.04$	$0.93\pm 0.05$
$\mathrm{NormRes}$ (%)	$94\pm 6$	$116\pm 3$	$92\pm 6$	$116\pm 7$
NonRes$1\pi$ (%)	$43\pm 4$	$46\pm 4$	$44\pm 4$	$46\pm 4$
NonRes$2\pi$ (%)	$166\pm 32$	99$\pm$31	$161\pm 33$	$120\pm 32$
$\pi$-iso	$1.0$ (limit)	1.0 (limit)	1.0 (limit)	1.0 (limit)
FrAbs (%)	$109\pm 16$	$48\pm 21$	-	-
FrInel (%)	-	-	$109\pm 24$	$132\pm 27$
Lag. $R_1$	-	$0.32\pm 0.06$	-	$0.37\pm 0.09$
Lag. $R_2$	-	$0.5$ (limit)	-	$0.60\pm 0.16$
MINERvA ${\chi }^2$	242.3	212.2	240.7	215.7
${\chi }_{\text{pen}}^2$	9.3	0.7	11.1	0.5
Total${\chi }^2$	251.6	212.9	251.8	216.2
$N_{\mathrm{DoF}}$	145	143	145	143

Table 9. Table 9: Individual channel FrAbs + low- Q 2 superscript 𝑄 2 Q^{2} tuning results.

Parameter	${\nu }_{\mu }$CC1${\pi }^{\pm }$	${\nu }_{\mu }$CC$N{\pi }^{\pm }$	${\nu }_{\mu }$CC$1{\pi }^0$	${\overline{\nu }}_{\mu }$CC$1{\pi }^0$
$M_A^{\mathrm{res}}$ (GeV)	$0.93\pm 0.02$	$0.92\pm 0.02$	$0.96\pm 0.05$	$0.94\pm 0.05$
$\mathrm{NormRes}$ (%)	$115\pm 3$	$117\pm 3$	$114\pm 7$	$115\pm 7$
NonRes$1\pi$ (%)	$43\pm 4$	$43\pm 4$	$45\pm 4$	$43\pm 4$
NonRes$2\pi$ (%)	$300$ (limit)	$70\pm 28$	$300$ (limit)	$300$ (limit)
$\pi$-iso	1 = Iso (limit)	1 = Iso (limit)	1 = Iso (limit)	1 = Iso (limit)
FrAbs (%)	92 $\pm$ 65	$79\pm 40$	$74\pm 22$	$34\pm 35$
Lag. $R_1$	$0.53\pm 0.16$	$0.43\pm 0.13$	$0.21\pm 0.14$	$0.14\pm 0.22$
Lag. $R_2$	$0.50$ (limit)	$0.50$ (limit)	$0.63\pm 0.31$	$1.00$ (limit)
MINERvA ${\chi }^2$	32.2	55.7	71.2	27.7
${\chi }_{\text{pen}}^2$	0.1	0.4	0.5	0.0
Total ${\chi }^2$	32.3	56.1	71.7	27.7
$N_{\mathrm{DoF}}$	33	34	30	31

Table 10. Table 10: Individual channel FrInel + low- Q 2 superscript 𝑄 2 Q^{2} tuning results.

Parameter	${\nu }_{\mu }$CC1${\pi }^{\pm }$	${\nu }_{\mu }$CC$N{\pi }^{\pm }$	${\nu }_{\mu }$CC$1{\pi }^0$	${\overline{\nu }}_{\mu }$CC$1{\pi }^0$
$M_A^{\mathrm{res}}$ (GeV)	$0.93\pm 0.02$	$0.91\pm 0.02$	$0.95\pm 0.05$	$0.94\pm 0.05$
$\mathrm{NormRes}$ (%)	$116\pm 3$	$117\pm 3$	$114\pm 7$	$115\pm 7$
NonRes$1\pi$ (%)	$43\pm 4$	$43\pm 4$	$44\pm 4$	$43\pm 4$
NonRes$2\pi$ (%)	$300$ (limit)	$78\pm 28$	$300$ (limit)	300 (limit)
$\pi$-iso	1 = Iso (limit)	1 = Iso (limit)	1 = Iso (limit)	1 = Iso (limit)
FrInel (%)	$179\pm 63$	$173\pm 37$	$8\pm 125$	$103\pm 57$
Lag. $R_1$	$0.49\pm 0.14$	$0.38\pm 0.13$	$0.25\pm 0.17$	$0.31\pm 0.26$
Lag. $R_2$	$0.50$ (limit)	$0.50$ (limit)	$0.76\pm 0.37$	$1.00$ (limit)
MINERvA ${\chi }^2$	30.8	52.1	69.5	30.9
${\chi }_{\text{pen}}^2$	0.1	0.6	0.2	0.0
Total ${\chi }^2$	30.9	52.7	69.7	30.9
$N_{\mathrm{DoF}}$	33	34	30	31

Table 11. Table 11: Channel by channel contributions to the χ 2 superscript 𝜒 2 \chi^{2} at for the GENIE tunings with and without the low- Q 2 superscript 𝑄 2 Q^{2} correction included.

Distribution	Channel	$N_{\mathrm{bins}}$	FrAbs Tune	FrAbs + low-$Q^2$ Tune	FrInel Tune	FrInel + low-$Q^2$ Tune
$p_{\mu }$ (Rate)	${\nu }_{\mu }$CC1${\pi }^{\pm }$	8	12.0	10.8	12.3	10.9
	${\nu }_{\mu }$CC$N{\pi }^{\pm }$	9	26.1	16.2	26.8	17.9
	${\nu }_{\mu }$CC$1{\pi }^0$	8	19.0	26.2	19.3	26.9
	${\overline{\nu }}_{\mu }$CC$1{\pi }^0$	9	6.2	7.1	6.3	7.2
${\theta }_{\mu }$ (Shape)	${\nu }_{\mu }$CC1${\pi }^{\pm }$	9	7.5	7.4	7.4	7.1
	${\nu }_{\mu }$CC$N{\pi }^{\pm }$	9	4.0	6.3	4.1	5.6
	${\nu }_{\mu }$CC$1{\pi }^0$	9	44.5	20.0	45.6	20.5
	${\overline{\nu }}_{\mu }$CC$1{\pi }^0$	9	10.2	7.0	10.3	6.9
$T_{\pi }$ (Shape)	${\nu }_{\mu }$CC1${\pi }^{\pm }$	7	2.5	2.5	2.3	2.4
	${\nu }_{\mu }$CC$N{\pi }^{\pm }$	7	31.2	28.9	29.4	27.7
	${\nu }_{\mu }$CC$1{\pi }^0$	7	30.9	27.1	29.9	32.0
	${\overline{\nu }}_{\mu }$CC$1{\pi }^0$	7	16.6	15.7	16.0	18.7
${\theta }_{\pi }$ (Shape)	${\nu }_{\mu }$CC1${\pi }^{\pm }$	14	13.0	13.4	12.6	12.6
	${\nu }_{\mu }$CC$N{\pi }^{\pm }$	14	6.9	7.0	6.2	6.3
	${\nu }_{\mu }$CC$1{\pi }^0$	11	8.3	12.2	8.9	9.4
	${\overline{\nu }}_{\mu }$CC$1{\pi }^0$	11	3.4	4.4	3.5	3.7
Total ${\chi }^2$		148	242.3	212.2	240.7	215.7

Equations

$${\chi }^2$$

$$+\sum _{kij}^{N_k}{\Delta }_{k,i}^S(S^{-1}{)}_{ij}{\Delta }_{k,j}^S$$

$${\Delta }_i$$

$${\Delta }_{k,i}^S$$

$${\chi }_{\text{pen}}^2=\sum _{i,j}^{N=3}(x_i-f_i){\left(M^{-1}\right)}_{ij}(x_j-f_j),$$

$$R=\frac{A}{1+\exp \{1-\sqrt{Q^2}/Q_0\}},$$

$$R(Q^2<x_3)$$

$$+\frac{(Q^2-x_1)(Q^2-x_2)}{(x_3-x_1)(x_3-x_2)}.$$

$$w(Q^2)=1-(1-R_1)(1-R(Q^2){)}^2.$$

Full text

I Introduction

In recent years, experimental groups have started to publish neutrino interaction cross-section measurements on nuclear targets in terms of measurable final state particle content, instead of inferred initial interaction channels. This avoids the problem of correcting for complex nuclear effects to make a measurement in terms of the initial interaction channels. For example, events with only a single pion can be produced by the decay of hadronic resonances formed at the primary neutrino interaction, followed by loss of a nucleon from the resonance’s decay as a result of final state interactions (FSI) within the nuclear medium. Such events can also be produced by other sequences of interactions, such as a deep inelastic collision where only a single pion is produced after FSI. A measurement of charged current events with one identified pion in the final state is a benchmark for models, independent of the details of how each model assesses any particular interaction channel’s contributions to that final state. The limitation of giving results in terms of final state particle content is that FSI are important, and result in the contribution of many different interaction channels into a specific final state.

There are tensions between published data from the T2K, MiniBooNE, and MINERvA experiments Mosel and Gallmeister (2017); Alvarez-Ruso et al. (2018); R. González-Jiménez and Jachowicz (2018); Wilkinson et al. (2016); Mahn et al. (2018). These tensions exist in the charged current production of both zero and one pion final states, and a model has yet to emerge that can reliably simulate all experiments at once. This is troubling, as current and future neutrino oscillation experiments require a cross section model which is predictive across the range of energies covered by these experiments and for a variety of targets.

The differences in neutrino fluxes, scattering targets, available phase space and signal definitions between experiments make it difficult to diagnose the exact causes of disagreement within the global data set. In particular, as results must be averaged over the neutrino flux distribution of each experiment, it is difficult to disentangle the energy dependence of an observed deficiency in a particular model, and decide how uncertainties should be propagated in neutrino energy. Tensions between measurements from a single experiment can uncover fundamental problems with a model which should be addressed, before considering the more difficult issue of developing, or empirically tuning a model which fits data from multiple experiments. In this work, we employ published MINERvA pion production data. The cross-section measurements utilized in this effort have not been reanalyzed or modified in any way.

NUISANCE Stowell et al. (2017) was developed to provide the neutrino scattering community with a flexible framework in which various neutrino interaction generators can be validated and empirically tuned to data. Its structure allows for generator tunings to be easily adapted to account for changes in the underlying model or data. In this work, the default pion production model in the GENIE Andreopoulos et al. (2009, 2015) neutrino interaction generator is tuned to MINERvA data. Although more sophisticated pion production models exist (e.g.  Kabirnezhad (2018); Leitner et al. (2009); Nieves et al. (2011); Hernadez and Nieves (2017)), GENIE is widely used by the neutrino scattering community, and its model uncertainties have a central importance to the field. Although the work is only directly applicable to one generator, the methods developed in this paper are easily adaptable to different generators. All the data and methods are publicly available and integrated into the open source NUISANCE framework, facilitating similar studies using other generators and models.

In Section II, the data are reviewed and the goodness-of-fit test statistic is defined for the tuning process. Section III describes the default GENIE pion production model, and reviews comparisons of this model to data. In Section IV, the parameter reweighting package in GENIE is discussed along with the specific parameters tuned therewith. We also discuss other corrections to the GENIE model made to improve agreement with bubble chamber data Rodrigues et al. (2016); Wilkinson et al. (2014). In Section V, we tune additional systematic parameters in GENIE to improve agreement with the MINERvA data in combination with the bubble chamber data. In Section VI, additional low-$Q^{2}$ ad hoc corrections are added to the model to resolve observed tensions, motivated by the need for similar corrections observed at both MINOS Adamson et al. (2015) and MiniBooNE Aguilar-Arevalo et al. (2011). Finally, in Section VII we present our conclusions.

II Data included in the fits

We tune to four of MINERvA’s published charged current pion production measurements taken on a polystyrene scintillator target: $\nu_{\mu}$CC1$\pi^{\pm}$ Eberly et al. (2015), $\nu_{\mu}$CC$N\pi^{\pm}$ McGivern et al. (2016), $\nu_{\mu}$CC$1\pi^{0}$ Altinok et al. (2017) and $\bar{\nu}_{\mu}$CC$1\pi^{0}$ McGivern et al. (2016), summarized in Table 1111In “$\nu_{\mu}$CC$N\pi^{\pm}$ ”, the $N$ indicates one or more identified pions and does not refer to a nucleon.. The MINERvA detector Aliaga et al. (2014) does not determine the polarity of charged pions. The fraction of $\pi^{-}$ in $\nu_{\mu}$CC1$\pi^{\pm}$ sample is small ($\sim$ 1%). Furthermore, the $\nu_{\mu}$CC1$\pi^{\pm}$ and $\nu_{\mu}$CC$N\pi^{\pm}$ signal definition allows for any number of neutral pions. Approximately 3% of the MINERvA $\nu_{\mu}$CC1$\pi^{\pm}$ signal events have at least one neutral pion in the final state. All four analyses include signal definition cuts on the true “reconstructed” mass of the hadronic system assuming the struck nucleon is at rest, $W_{\mathrm{rec}}$, and the true neutrino energy $E_{\nu}$.

The kinematic variable distributions used in this work are the momentum and angle of the outgoing muon with respect to the incoming neutrino beam, $p_{\mu}$ and $\theta_{\mu}$, and the kinetic energy and angle of the outgoing pion with respect to the incoming neutrino beam, $T_{\pi}$ and $\theta_{\pi}$. In the $\nu_{\mu}$CC$N\pi^{\pm}$ channel, where there is at least one $\pi^{\pm}$ in the final state, there is one entry in the distributions of $\theta_{\pi}$ and $T_{\pi}$ for each $\pi^{\pm}$ in an event. The data are reported as efficiency corrected results unfolded to true kinematic variables, which may introduce model dependence. This is notably problematic in regions of low efficiency—present in the charged pion channels at $\theta_{\pi}\sim 90^{\circ}$, $T_{\pi}<50~{}\mbox{MeV}$ and $T_{\pi}>350~{}\mbox{MeV}$, where the signal efficiency is zero Eberly et al. (2015). The pion selection cuts, not present in the signal definition, remove about 50% of signal events, with little dependence upon the muon variables, but a clear impact on the shape of the pion kinematic variables.

The published cross-sections are one dimensional with correlations provided between the bins within each distribution. No correlations are provided between measurements of different final states, or between different one-dimensional projections of the same measurement. These correlations are expected to be large, coming predominantly from flux and detector uncertainties. Additionally, the $\nu_{\mu}$CC1$\pi^{\pm}$ event sample is a subset ($\sim$64%) of the $\nu_{\mu}$CC$N\pi^{\pm}$ event sample, and including both channels introduces a statistical correlation. Not assessing correlations between the distributions, while common practice in this field, is a limitation when tuning models to multiple data sets. It introduces a bias in the $\chi^{2}$ statistic that is difficult to quantify, and requires imposing ad hoc uncertainties Wilkinson et al. (2016) as the test-statistic is not expected to follow a $\chi^{2}$ distribution for the given degrees of freedom.

The covariance matrices contain a flux-dominated normalization component which we expect to be fully correlated across all distributions. To account for the correlated uncertainty, we use the full covariance matrix, $M_{ij}$, for the $p_{\mu}$ distribution and shape-only covariance matrices, $S_{ij}$, for the other three distributions in each of the topologies. Whilst any distribution could set the normalization constraint, the shape of the $p_{\mu}$ distribution for each channel was chosen since it was found to be relatively insensitive to model variations and had good shape agreement with the data. The joint $\chi^{2}$ is therefore defined as the sum of the full $p_{\mu}$ $\chi^{2}$ and shape-only $\theta_{\mu}$, $T_{\pi}$ and $\theta_{\pi}$ $\chi^{2}$’s:

[TABLE]

where $i$ and $j$ are bin indices,

[TABLE]

and $d_{k,i}$ and $m_{k,i}$ are the data and MC values respectively for the $i^{th}$ bin in the $k^{th}$ distribution. The shape-only covariance matrices are provided in the public data release for the $\nu_{\mu}$CC1$\pi^{\pm}$ and $\nu_{\mu}$CC$1\pi^{0}$ measurements, and the method of Ref. Katori (2008) (section 10.6.3) was used to extract them for the $\nu_{\mu}$CC$N\pi^{\pm}$ and $\bar{\nu}_{\mu}$CC$1\pi^{0}$ channels.

III Pion Production in GENIE

This analysis begins with version 2.12.6 of GENIE, which is close to what is used by MINERvA, T2K, NOvA and MicroBooNE. We use the Smith-Moniz relativistic Fermi gas (RFG) model Smith and Moniz (1972) with an added high momentum tail as per Bodek and Ritchie Bodek and Ritchie (1981). The Valencia random phase approximation screening Nieves et al. (2004) is applied as a weight to quasielastic events. The two-particle two-hole process is simulated using the Valencia model Nieves et al. (2011); Gran et al. (2013). MINERvA currently uses a modification of v2.8.4 Altinok et al. (2017); Patrick et al. (2018); Mislivec et al. (2018) with an increased rate for the Valencia two-particle two-hole process; that modification is not used here. An important difference in single pion production between v2.8.x and v2.12.x is the angular distributions of single pion events in the Rein-Sehgal model, discussed below. A sample of 2.5 million events were generated using the MINERvA flux predictions Aliaga et al. (2016), a polystyrene target and the official GENIE 2.12.6 splines Andreopoulos et al. (2017).

To simulate pion production, GENIE uses the Rein-Sehgal (RS) model Rein and Sehgal (1981) with a hadronic invariant mass cut of $W\leq 1.7~{}\mbox{GeV}$. Of the 18 resonances in the RS model, the $\Delta(1600)$ and $N(1990)$ were not included due to their unclear experimental status at the time of implementation. Resonance-resonance and resonance-nonresonance interference terms are not included. Lepton mass terms are only included in calculating phase space limits and are neglected when calculating the cross sections. A discussion of the limitations of this simplification can be found in Ref. Graczyk and Sobczyk (2008). In earlier versions—including v2.8.4—the pion-nucleon distribution was isotropic in the resonance rest frame, but was changed in 2.12.x. Here we use the non-isotropic model as our default and reweight to the isotropic distribution, explained later. The RS nonresonant background is not used by GENIE; rather, a deep inelastic scattering (DIS) model is extended to cover that invariant mass region. The DIS model uses the Bodek-Yang parametrization Bodek and Yang (2003), and the AGKY model to describe hadronization Yang et al. (2009). In the AGKY model, the KNO model Koba et al. (1972) is used for $W\leq 2.3~{}\mbox{GeV}$ and PYTHIA T. Sjostrand and Skands (2008) is used for $W\geq 3.0~{}\mbox{GeV}$, with a smooth transition in between the two, implemented by randomly selecting the results of one model or the other for each event.

In addition to pion production on a single nucleon, it is also possible for a neutrino to produce a pion by scattering coherently off the nucleus. GENIE uses the Rein-Sehgal coherent pion production model Rein and Sehgal (1983, 2007), including the effect of lepton masses in the cross-section calculation. MINERvA has found that the RS coherent pion production model needs to be suppressed by $\sim 50\%$ at $T_{\pi}$ $<500~{}\mbox{MeV}$ to agree with data Mislivec et al. (2018). This correction also moves the shape of the $T_{\pi}$ spectrum closer to the predictions of the Berger-Sehgal coherent model Berger and Sehgal (2009). The $\nu_{\mu}$CC1$\pi^{\pm}$ channel has a small contribution from coherent production in the lowest $Q^{2}$ bins but the inclusion of this suppression has only a small effect on the MC predictions. To maintain a model similar to that currently being used by MINERvA, this suppression is included in the analysis presented in Sections IV and onwards.

The “hA Intranuke” effective cascade model Dytman (2009) is used to model pion and nucleon FSI. In this model, the effect of intranuclear scattering is parameterized as a single cascade step applied to each particle emanating from the primary interaction. This model steps hadrons through a nucleus of radius $r\sim A^{1/3}$ and a nuclear density function derived from electron scattering data. The hadron’s mean free path is determined from tabulated hadron-proton and hadron-neutron cross sections Arndt et al. (2006). The probability to interact with the nucleus is high; it is e.g. $\sim$73% for a pion from an $E_{\nu}$ = 3 GeV quasielastic event in carbon. When a FSI occurs, the possible interactions (absorption, pion production, knockout, charge exchange, elastic scatter) are chosen according to their proportions for iron.

Default GENIE predictions separated by nucleon level interaction channels for the MINERvA data are shown in Fig. 1. The shape of the $p_{\mu}$ distributions agree well with the data for all four measurements. However, the model overestimates the cross section for $\pi^{\pm}$ production and as a result the $\chi^{2}$ for the $\nu_{\mu}$CC1$\pi^{\pm}$ and $\nu_{\mu}$CC$N\pi^{\pm}$, given in the fourth column (“Default”) of Table 2, are large. The model overestimates $\theta_{\mu}$ below $<5^{\circ}$ in the $\pi^{0}$ channels, although it does correctly predict the shape of the $\theta_{\mu}$ distribution in the $\pi^{\pm}$ channels. The model underestimates the production rate at large $\theta_{\mu}$ in $\nu_{\mu}$CC$1\pi^{0}$. The shape of the $T_{\pi}$ distribution is in larger disagreement for $\nu_{\mu}$CC$N\pi^{\pm}$ data than for $\nu_{\mu}$CC1$\pi^{\pm}$. Since the $\nu_{\mu}$CC$N\pi^{\pm}$ distributions summed over all identified $\pi^{\pm}$, redistributing kinetic energy between $\pi^{\pm}$ in events with more than one $\pi^{\pm}$ could resolve some of this tension. The $\pi^{0}$ channels are under-predicted at low $T_{\pi}$. Finally, GENIE predictions are too high in magnitude at $\theta_{\pi}\approx 50^{\circ}$ in both the $\nu_{\mu}$CC1$\pi^{\pm}$ and $\nu_{\mu}$CC$N\pi^{\pm}$ channels, and the prediction has the wrong shape in the $\nu_{\mu}$CC1$\pi^{\pm}$ channel. Comparisons using the transport theory based GiBUU model Buss et al. (2012) show similar shape disagreements despite GiBUU’s use of an advanced semiclassical cascade model to simulate FSI Mosel and Gallmeister (2017).

Each of the measurements are shown as MC/data ratio distributions in Fig. 2. Similar comparisons between the MiniBooNE and MINERvA experiments are found in Ref. Mahn et al. (2018). The shape-only data sets ($\theta_{\mu}$, $\theta_{\pi}$, $T_{\pi}$) were normalized to match the data before the ratio was taken and the error bars in Fig. 2 reflect the extracted shape-only uncertainties on the data, so that the distributions reflect their contributions to the total $\chi^{2}$.

IV Tunable parameters in the GENIE Model

The GENIE event generator allows assessment of systematic uncertainties through the GENIE reweighting package. A large number of event weighting “dials” are included to allow model uncertainties to be evaluated. The dials adjusted in this note are summarized in Table 3 and are chosen because of their connection to the kinematic variables and interaction modes studied herein.

Experiments often use variations in the charged-current resonant axial mass, $M_{\mathrm{A}}^{\mathrm{res}}$, as a systematic uncertainty which varies both the normalization and $Q^{2}$ shape of resonant interactions along with variations in a total resonant cross-section normalization dial, $\mathrm{NormRes}$. Variations in $\mathrm{NormRes}$ approximates the behaviour of varying $F_{\mathrm{A}}(0)$ in the axial form factor in the Rein-Sehgal model. Since low $\theta_{\mu}$ correlates with low $Q^{2}$, variations in $M_{\mathrm{A}}^{\mathrm{res}}$ have the largest effect on the shape of the muon angular distributions as shown in Fig. 3, and have a small effect on the $\theta_{\pi}$ spectrum.

Dials are available to vary the normalization of the nonresonant $1\pi$ production channels in GENIE (e.g. NonRESBGvnCC1pi, NonRESBGvpCC1pi) but each dial introduces similar modifications to the predictions. To reduce the number of free parameters in the fit described in Section V, these dials were grouped into a single background scaling for nonresonant $1\pi$ production, NonRes$1\pi$, following the approach in Ref. Rodrigues et al. (2016); Wilkinson et al. (2014). A similar treatment was also applied to nonresonant $2\pi$ production, NonRes$2\pi$, with the neutrino and antineutrino related parameters assumed to be 100% correlated in both cases. The effects of varying the nonresonant contributions are shown in Fig. 4. Variations in the NonRes$2\pi$ dial introduce a large change in normalization for the $\nu_{\mu}$CC$N\pi^{\pm}$ channel and has a minor effect in the other single pion channels as the fraction of multi-$\pi$ events is small.

Reanalysis of data from ANL and BNL bubble chambers has provided a tuning of GENIE’s single pion production model on free nucleons. The work showed that a small shift in $M_{\mathrm{A}}^{\mathrm{res}}$ was required to model the low-$Q^{2}$ region and a large suppression of the nonresonant $\pi$ production ($-54\%$) was required to match the observed cross sections of $\pi^{+}$ and $\pi^{0}$ production. The reanalysis used the measured ratios of the rates of single $\pi$ production to CCQE measurements to cancel errors in the flux. We note that by using CCQE data multiple times, they introduce hidden correlations which may have a small effect on the postfit uncertainties. However, as the single pion statistical uncertainties at ANL Radecky and others (ANL) and BNL Kitagaki and others (BNL) were magnitudes higher than the CCQE statistical uncertainty Barish and others (ANL); Kitagaki and others (BNL), the effect was neglected in that work, and is also neglected here. The resulting parameter tunes shown in Table 4 and Fig. 5 have been partially adopted by MINERvA and NOvA which both apply the nonresonant rescaling of 43% but leave the other parameters unchanged.

Fig. 6 shows MINERvA data and the predictions of GENIE when its output has been reweighted to reflect the parameter changes of Table 4. The channel-by-channel contributions to the $\chi^{2}$ are given in the fifth column (“ANL/BNL”) of Table 2. Incorporating the parameter changes improves the total normalization agreement in the $p_{\mu}$ distributions for $\nu_{\mu}$CC1$\pi^{\pm}$ and $\nu_{\mu}$CC$N\pi^{\pm}$. The $\chi^{2}$ for the $p_{\mu}$ distribution is also improved in the $\bar{\nu}_{\mu}$CC$1\pi^{0}$ channel, even though the ANL/BNL data is from neutrino interactions only. The $\chi^{2}$ for the $p_{\mu}$ distribution in the $\nu_{\mu}$CC$1\pi^{0}$ channel is somewhat worse as the parameter tunes reduce the predicted nucleon $\nu_{\mu}$CC$1\pi^{0}$ cross section. The modification of $M_{\mathrm{A}}^{\mathrm{res}}$ shifts the $\theta_{\mu}$ predictions to lower values, increasing the $\chi^{2}$ contributions. The $T_{\pi}$ and $\theta_{\pi}$ distributions change mostly by normalization, having a smaller effect on the $\chi^{2}$. The overall agreement of GENIE with MINERvA data is not improved by incorporating the ANL/BNL information. Indeed, the total $\chi^{2}$ increases, largely because of the $\chi^{2}$ contributions from the $\theta_{\mu}$ distributions.

GENIE provides a dial that influences the resonances’ decay into the pion-nucleon system in the resonance rest frame, $\pi$-iso, and allows events to be reweighted continuously between the default anisotropic distribution ($\pi$-iso = 0) and the isotropic distribution ($\pi$-iso = 1). The Adler angle222The Adler angle is the angle between the pion and the three momentum transfer in the resonance rest frame Adler (1968). is highly sensitive to the $\pi$-iso parameter and has been measured by neutrino induced pion production experiments on single nucleons, such as ANL Radecky and others (ANL), BNL Kitagaki and others (BNL), BEBC Allen et al. (1980); Allasia et al. (1983) and FNAL Bell. et al. (1978). Nucleon data strongly prefers an anisotropic process, as shown in Fig. 7. Nonetheless, $\pi$-iso has some impact, albeit one that does depend on how FSI are modelled, on the shape of MINERvA $\theta_{\pi}$ and $T_{\pi}$ distributions, seen in the bottom of Fig. 7, and was therefore included in this work.

The GENIE hA model for FSI has uncertainties from the $\pi-A$ cross-section data to which the model was tuned. The total $\pi-A$ cross section has a stronger constraint than each of the individual interaction cross sections, so GENIE provides dials to vary the fractional contribution of each component. The available fractional dials are pion absorption (FrAbs), pion inelastic scattering (FrInel), pion elastic scattering, pion charge exchange and pion production.

V Tuning the GENIE Model

Fig. 6 and Table 2 show the unsatisfactory agreement of the GENIE prediction against MINERvA data. The disagreement worsens after incorporating the prior constraint from ANL and BNL data; this correction, based on nucleon data, is inadequate. This section describes fits that improve the agreement with MINERvA data. The parameters $M_{\mathrm{A}}^{\mathrm{res}}$, $\mathrm{NormRes}$ and NonRes$1\pi$ are included in the fits with a penalty term added to the $\chi^{2}$ from the ANL and BNL data. The penalty term uses the covariance, $M$, shown in Fig. 5:

[TABLE]

where $x_{i}$ are the parameter values $i$ at each iteration of the fit, and $f_{i}$ are the parameter values from the fit to ANL and BNL data. The GENIE default model is strongly disfavoured with $\chi^{2}_{\textrm{pen}}$ = 299.3, but changing NonRes$1\pi$ to 43% while leaving all the other parameters at their default values reduces the $\chi^{2}_{\textrm{pen}}$ to 21.8, showing that the largest tension is due to the NonRes$1\pi$ parameter.

The $\pi$-iso dial is allowed to vary in the range 0–1 in the fit, corresponding to a continuous variation between an RS angular distribution and an isotropic distribution for $\Delta(1232)$ decay. To avoid the normalization of the $\nu_{\mu}$CC$N\pi^{\pm}$ measurement pulling parameters in the $\nu_{\mu}$CC1$\pi^{\pm}$ model, the NonRes$2\pi$ dial was allowed to vary between 0-300% of the nominal value.

When varying one of the five hA pion FSI dials, GENIE automatically adjusts the remaining parameters to preserve the total pion cross-section and maintain agreement with pion-nucleus scattering data. This “cushion” technique introduced instabilities in the $\chi^{2}$ surface, so it was not possible to include multiple pion FSI parameters in a simultaneous fit. Instead we performed fits with only one of the FSI parameters floating. No $\chi^{2}$ penalty terms were added for the FSI dial in either tuning: the parameters were driven solely by MINERvA data. The charge exchange and pion production dials had small contributions to the overall $\chi^{2}$ for the selected data, forcing the parameters to be inflated beyond +3$\sigma$ of GENIE’s recommendation, with large post-fit uncertainties. Furthermore, the pion elastic scattering parameter is strongly constrained by external data, so its 1$\sigma$ variation has a small impact on the MINERvA distributions. The non-FSI fit parameters’ (e.g. $M_{\mathrm{A}}^{\mathrm{res}}$) central values and uncertainties all agreed for the five fits. Here we present the results from the FrAbs and FrInel fits.

The NUISANCE interface to MINUIT2 James and Roos (1975) was used to perform the fits. At each iteration, the GENIE-ReWeight package was used on an event-by-event basis to update the MC predictions before the total $\chi^{2}$ was calculated. The uncertainties in the fitted parameters were determined using the HESSE routine in MINUIT2. The best fit results from the joint tuning are shown in Table 5. Fig. 8 shows the ratios of the best fit prediction to the data for all four kinematic variables of interest when the pion absorption FSI parameter (FrAbs) is floated in the fit; Fig. 9 is the same, but when the inelastic scattering FSI parameter (FrInel) is floated. Notably, the two FSI fits are very similar in both minimum $\chi^{2}$ and best-fit parameter values.

Comparing to the results of the ANL and BNL reanalysis, larger values of $M_{\mathrm{A}}^{\mathrm{res}}$ and smaller values of $\mathrm{NormRes}$ were found by the fit, pulling the parameters closer to GENIE nominal. The NonRes$1\pi$ parameter is strongly bound by the bubble chamber data and the MINERvA data did little to improve on this constraint. The penalty term contributed to the $\chi^{2}$ by 9.3 for the FrAbs fit and 11.1 for the FrInel fit. This is a significant improvement over the default, but indicates that there is mild tension between the nucleon and nuclear data. The post-fit correlation matrices are provided in Fig 10. The ANL/BNL input correlations are largely maintained in our fit.

Tables 6 and 7 show the results when individual MINERvA data were tuned in separate fits. Since three of the four channels were removed in these fits the constraint from data is weakened and the total $\chi^{2}$ is steered by the bubble chamber $\chi^{2}$ penalty. The individual channel fits also found values at the 300% limit for NonRes$2\pi$ dial, except in the $\nu_{\mu}$CC$N\pi^{\pm}$ channel, where the result was unchanged by the fit. Only the $\nu_{\mu}$CC$N\pi^{\pm}$ channel has a significant contribution from nonresonant 2$\pi^{\pm}$ production. In the other fits, the parameter is largely unconstrained and has little impact on the fitted distributions. The $\chi^{2}$ per degree of freedom is indicative of a poor fit in the $\nu_{\mu}$CC$N\pi^{\pm}$ and $\nu_{\mu}$CC$1\pi^{0}$ channels, but not in the $\nu_{\mu}$CC1$\pi^{\pm}$ or $\bar{\nu}_{\mu}$CC$1\pi^{0}$ channels. Furthermore, the $\nu_{\mu}$CC$1\pi^{0}$ shows the strongest $\chi^{2}$ penalty, indicating tension with the ANL/BNL prior. Given the different kinematic regions covered by the channels (see Table 1) and the different physics (e.g. fraction of coherent pion production) it is difficult to infer what combination of effects are at work. Isotropic emission was preferred in all fits, driven by the $\theta_{\pi}$ distributions. Disagreements in the $\theta_{\pi}$ spectrum are clearly seen in the data/MC ratios of Fig. 8 and 9, and the large $\chi^{2}$ values observed for the $\nu_{\mu}$CC$N\pi^{\pm}$ and $\nu_{\mu}$CC$1\pi^{0}$ channels.

The individual $\chi^{2}$ contributions in the joint tuning best fit, shown in sixth and seventh columns (“FrAbs Tune” and “FrInel Tune”) of Table 2, show that not all distributions in all channels benefit from the model variations, as the default GENIE fits have a better $\chi^{2}$ for some distributions. In particular, the $\nu_{\mu}$CC$1\pi^{0}$ channel distributions have worse agreement after the tuning, with only the $\theta_{\pi}$ distribution improving in $\chi^{2}$, whereas all channels benefit from the shift to isotropic emission. While there is an overall improvement over the ANL/BNL tune when comparing the combined $\chi^{2}$ results, Figs. 8 and 9 show that there are still unresolved shape disagreements in both the $T_{\pi}$ and $\theta_{\mu}$ kinematics.

The tension between MINERvA’s nuclear data and the constraints from ANL and BNL nucleon data is difficult to confidently pinpoint; the lack of lepton mass effects Berger and Sehgal (2007), modification to the resonance propagator in the nucleus Singh et al. (1998); Oset and Salcedo (1987), missing diagrams describing the non-resonant background contributions Kabirnezhad (2018), dynamical coupled channels Nakamura et al. (2017), interactions on correlated initial states, and the pion FSI model Mosel and Gallmeister (2017) are all part of an incomplete list of possible culprits.

VI Ad hoc $Q^{2}$ Supression

Further modifications beyond the standard GENIE dials are required to resolve the observed tensions. Fig. 11 (not used for any tuning) shows the $Q^{2}$ distributions observed at MINERvA in for our tunes. The data is below the predictions of the tunes of Section V at low values of $Q^{2}$. There are also also differences at low $\theta_{\mu}$, as shown in Fig. 8 and 9. Measurements of $\nu_{\mu}$CC1$\pi^{\pm}$ and $\nu_{\mu}$CC$1\pi^{0}$ interactions on mineral oil at MiniBooNE have shown a data/MC shape discrepancy for the RS implementation in the NUANCE model Casper (2002); Przewlocki (2009) in both $Q^{2}$ and $\cos\theta_{\mu}$ distributions Aguilar-Arevalo et al. (2011); Aquilar-Arevalo et al. (2011). In the MINOS quasielastic analysis Adamson et al. (2015) on iron, which used NeuGen Gallagher (2002), a similar disagreement was observed when studying pion production dominated sidebands. Indeed, concerns about low-$Q^{2}$ modeling date back almost a decade Graczyk and Sobczyk (2008). The data from MINOS and MiniBooNE experiments and the MINERvA data on CH studied herein suggest that the RS implementation common to each of the generators needs to be suppressed at low $Q^{2}$. Collective effects, which are usually modeled in the random phase approximation, are known to affect the $Q^{2}$ distribution of neutrino-nucleus reactions at low $Q^{2}$. Motivated by these considerations, we attempted to improve the $\theta_{\mu}$ modeling by introducing a $Q^{2}$-dependent correction to the model.

The MINOS collaboration suppression was expressed as

[TABLE]

where the free parameters $A=1.010$ and $Q_{0}=0.156~{}\mbox{GeV}$ were empirically extracted from bin-by-bin fits in $Q^{2}$ to the data, and a hard cut-off at $Q^{2}<0.7~{}\mbox{GeV}^{2}$ was imposed.

We chose an empirical function so that the shape of the suppression preferred by each of the MINERvA channels could be extracted. The empirical correction function is applied to events with a resonance decay inside the nucleus giving rise to a pion. Our suppression term is defined by choosing 3 points $(x_{i},R_{i})_{i=1,2,3}$ between $0.0<x<1.0$ and $0.0<R<1.0$, where $x\equiv Q^{2}$. Motivated the ANL/BNL curves in Fig. 11, the correction is assumed to approach unity as $Q^{2}$ approaches $0.7~{}\mbox{GeV}^{2}$, providing the constraint $(x_{3},R_{3})=(0.7~{}\mbox{GeV}^{2},1.0)$. Lagrange interpolation is used to derive a curvature from $R_{2}$ by assuming a simple interpolation between the points $(x_{1},0.0)$, $(x_{2},R_{2})$, and $(0.7~{}\mbox{GeV}^{2},1.0)$:

[TABLE]

This interpolation function is then used to calculate the correction for each event as

[TABLE]

where $R_{1}$ defines the magnitude of the correction function at the intercept, $x_{1}=0.0$. $x_{2}$ is chosen to be $Q^{2}=0.35~{}\mbox{GeV}^{2}$ so that $R_{2}$ describes the curvature at the centre point of the correction. Expressing the weights with Equations 6 and 7 ensures that the magnitude at $x_{2}$ always lies between $R_{1}$ and 1.0, avoiding parameter sets with large unphysical peaks in the correction function. Additionally, the squared term in Equation 7 ensures that $w(Q^{2})\rightarrow 1.0$ as $x\rightarrow x_{3}$, avoiding discontinuous steps in the weighting function at $x_{3}$. The fitted parameters $R_{1}$ and $R_{2}$ were limited to $0.0<R_{1}<1.0$ and $0.5<R_{2}<1.0$ to avoid extraneous solutions, e.g. double peaks.

The fit results are shown in Table 8. The correction from the fit with FrAbs taken as a free parameter are compared to the MINOS low-$Q^{2}$ correction in Fig. 12. Our fits obtain a suppression factor that is similar to the MINOS one, with almost identical suppression at $Q^{2}$ = 0, albeit with less curvature, particularly in the $\nu_{\mu}$CC1$\pi^{\pm}$ and $\nu_{\mu}$CC$N\pi^{\pm}$ channels. The correction factors from the fit with FrInel or FrAbs as free parameters give similar results.

The correlation matrices for the fits including a $Q^{2}$ dependent suppression are provided in Fig 13. Again, the ANL/BNL input prior covariance is maintained. The parameters largely correlate in the same way for the FrAbs and FrInel fit, and for the FrInel fit the $R_{1}$ and $R_{2}$ parameters are negatively correlated.

Figure 14 (Fig. 15) shows the ratio of the resulting fits to the MINERvA data when FrAbs (FrInel) is taken as a free parameter. As anticipated, the predictions now have better agreement with the data in regards to the $\theta_{\mu}$ distribution, and the $\chi^{2}$ values are improved by the introduction of our ad hoc low-$Q^{2}$ correction. Other fit parameters are for the most part unchanged by the introduction of the low-$Q^{2}$ correction. Furthermore, $M_{\mathrm{A}}^{\mathrm{res}}$ and $\mathrm{NormRes}$ are closer to their values when fitting ANL and BNL data, indicating the $Q^{2}$ correction alleviates the tension between nucleon and nuclear modeling. Fig. 11 shows the comparison of all our models directly against MINERvA data in $Q^{2}$. Although the tuning sees improvement in the $\chi^{2}$ for the $\nu_{\mu}$CC$1\pi^{0}$ and $\bar{\nu}_{\mu}$CC$1\pi^{0}$ distributions, the $\nu_{\mu}$CC1$\pi^{\pm}$ and $\nu_{\mu}$CC$N\pi^{\pm}$ distributions get worse, hinting at tensions in the charged and neutral pion production channels.

Tables 9 and 10 show the results of the fits to individual channels, and Table 11 shows the breakdown of contributions to the $\chi^{2}$ from the individual channels. The best fit $\chi^{2}$ value was significantly improved for each channel tuning when using a low-$Q^{2}$ suppression and the extracted parameters were consistent with the ANL/BNL tunings. Pion kinematic distributions are not improved, and in some cases are slightly worse, as a result of including the low-$Q^{2}$ suppression. It is clear from Table 8 (or by comparing Tables 6 and 9) that the low-$Q^{2}$ suppression has a similar effect in the fit to the FrAbs parameter. When the low-$Q^{2}$ suppression is introduced, FrAbs tends to consistently lower values. It is also clear that the $\bar{\nu}_{\mu}$CC$1\pi^{0}$ channel favors stronger low-$Q^{2}$ suppression than the other channels.

VII Concluding Remarks

We have adjusted the parameters of the GENIE model that are important for pion production to match MINERvA data in the $\nu_{\mu}$CC1$\pi^{\pm}$, $\nu_{\mu}$CC$N\pi^{\pm}$, $\nu_{\mu}$CC$1\pi^{0}$ and $\bar{\nu}_{\mu}$CC$1\pi^{0}$ channels, using the NUISANCE framework. We incorporate existing results which informs the GENIE model using ANL and BNL bubble chamber data from scattering off protons and deuterons. Fits of selected GENIE model parameters were done using the kinematic distributions $p_{\mu}$, $\theta_{\mu}$, $T_{\pi}$ and $\theta_{\pi}$. Parameter fits were performed with either the fraction of pions absorbed or the fraction of pions inelastically scattered in FSI as a floating parameter, with broadly similar conclusions for the two cases.

The results of the fit (see Table 5) show that the tuning improves the GENIE pion production model significantly, but tensions remain. The pull on the ANL/BNL prior demonstrates a tension between MINERvA nuclear target data and the light-target bubble-chamber data sets used to make the prior, indicating a deficiency in the GENIE nuclear model which cannot be fixed by modifying the available reweighting dials. Additionally, fitting to individual MINERvA pion production channels produces different best-fit parameters, demonstrating that GENIE cannot describe the different exclusive channels in a consistent manner with the available dials (shown in Tables 6 and 7). Because the four channels cover different kinematic regions (see Table 1) and contain different physics (e.g. different coherent pion production contributions or nonresonant processes), it is difficult to pinpoint the origin of the discrepancy between the model and the different MINERvA data sets.

Following experimental hints of discrepancies at low-$Q^{2}$ for a variety of cross-section measurements on nuclear targets, an additional empirical low-$Q^{2}$ suppression was introduced and the fits were repeated. Although the data showed a preference for a strong suppression at low-$Q^{2}$ and the agreement improved for $\theta_{\mu}$ and $Q^{2}$ distributions, tensions remain. In particular, fits to individual MINERvA channels still produced different results, and favor different parameter values for the low-$Q^{2}$ suppression.

The main conclusion of this work is that current neutrino experiments operating in the few–GeV region should think critically about single pion production models and uncertainties, as the Monte Carlo models which are currently widely used in the field are unable to explain multiple data sets, even when they are from a single experiment.

A key strength of this analysis is its development within the NUISANCE framework, allowing it to be easily repeated with alternate model assumptions, neutrino interaction generators, and different data. The developments presented here will be used in future iterations of this work, as the MINERvA collaboration works towards a GENIE model that provides a good description of all their available data, and can be easily applied to other measurements and experiments.

Acknowledgements

P.S., L.P., and C.V.C.W. would like to thank the UK Science and Technology Facilities Council (STFC) for Ph.D. funding support. P.S. acknowledges the Fermilab Neutrino Physics Center for the scholarship that funded this work and thanks Fermilab and the MINERvA collaboration for their hospitality during this work. C.W. acknowledges the support of the Swiss National Science Foundation and SERI.

This document was prepared by members of the MINERvA Collaboration using the resources of the Fermi National Accelerator Laboratory (Fermilab), a U.S. Department of Energy, Office of Science, HEP User Facility. Fermilab is managed by Fermi Research Alliance, LLC (FRA), acting under Contract No. DE-AC02-07CH11359. These resources included support for the MINERvA construction project, and support for construction also was granted by the United States National Science Foundation under Award No. PHY-0619727 and by the University of Rochester. Support for participating scientists was provided by NSF and DOE (USA); by CAPES and CNPq (Brazil); by CoNaCyT (Mexico); by Proyecto Basal FB 0821, CONICYT PIA ACT1413, Fondecyt 3170845 and 11130133 (Chile); by CONCYTEC, DGI-PUCP, and IDI/IGI-UNI (Peru); and by the Latin American Center for Physics (CLAF); NCN Opus Grant No. 2016/21/B/ST2/01092 (Poland). We thank the MINOS Collaboration for use of its near detector data. Finally, we thank the staff of Fermilab for support of the beam line, the detector, and computing infrastructure.