Tuning the GENIE Pion Production Model with MINERvA Data

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 still far from perfect. The development of these tunings within the NUISANCE framework means that they can easily be extended to other neutrino event generator models in the future, and compared with, or adapted to include, new datasets very easily.

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. The widely used pion production model in GENIE is compared to four MINERνA charged current pion production measurements using NUISANCE. Tunings, i.e., adjustments of model parameters, to help match GENIE to MINERνA and older bubble chamber data are presented. We find that scattering off nuclear targets as measured in MINERνA is not in good agreement with expectations based upon scattering off nucleon (hydrogen or deuterium) targets in existing bubble chamber data. An additional ad hoc correction for the low-Q 2 region, where collective nuclear effects are expected to be large, is 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 framework allows for straightforward extensions to other neutrino event generators and models, and allows omitting and including new datasets as they become available. DOI: 10.1103/PhysRevD.100.072005

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 (FSIs) 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 FSIs. 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 FSIs 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 Collaboration, and MINERνA experiments [1][2][3][4][5]. 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 dataset. 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 to 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 MINERνA pion production data. The crosssection measurements utilized in this effort have not been reanalyzed or modified in any way. NUISANCE [6] 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 [7,8] neutrino interaction generator is tuned to MINERνA data. Although more sophisticated pion production models exist (e.g., [9][10][11][12]), 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 Sec. 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 Sec. 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 [13,14]. In Sec. V, we tune additional systematic parameters in GENIE to improve agreement with the MINERνA data in combination with the bubble chamber data. In Sec. 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 [15] and MiniBooNE [16]. Finally, in Sec. VII we present our conclusions.

II. DATA INCLUDED IN THE FITS
We tune to four of MINERνA's published charged current pion production measurements taken on a polystyrene scintillator target: ν μ CC1π AE [17], ν μ CCNπ AE [18], ν μ CC1π 0 [19], andν μ CC1π 0 [18], summarized in Table I. 1 The MINERνA detector [20] does not determine the polarity of charged pions. The fraction of π − in the ν μ CC1π AE sample is small (∼1%). Furthermore, the ν μ CC1π AE and ν μ CCNπ AE signal definition allows for any number of neutral pions. Approximately 3% of the MINERνA ν μ CC1π AE 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 rec , and the true neutrino energy E ν .
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 μ and θ μ , and the kinetic energy and angle of the outgoing pion with respect to the incoming neutrino beam, T π and θ π . In the ν μ CCNπ AE channel, where there is at least one π AE in the final state, there is one entry in the distributions of θ π and T π for each π AE 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 θ π ∼ 90°, T π < 50 MeV, and T π > 350 MeV, where the signal efficiency is zero [17]. The pion selection cuts, not present in the signal definition, remove about 50% of the 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 onedimensional projections of the same measurement. These correlations are expected to be large, coming predominantly from flux and detector uncertainties. Additionally, the ν μ CC1π AE event sample is a subset (∼64%) of the ν μ CCNπ AE event sample, and including both channels introduces a statistical correlation. Not assessing correlations between the distributions, while a common practice in this field, is a limitation when tuning models to multiple datasets. It introduces a bias in the χ 2 statistic that is difficult to quantify, and requires imposing ad hoc uncertainties [4] as the test statistic is not expected to follow a χ 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 μ distribution and shape-only covariance matrices, S ij , for the other three distributions in each of the topologies. While any distribution could set the normalization constraint, the shape of the p μ 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 χ 2 is therefore defined as the sum of the full p μ χ 2 and shape-only θ μ , T π , and θ π χ 2 's: where i and j are bin indices, Signal definition GeV not applicable not applicable θ μ < 25°not applicable and d k;i and m k;i are the data and MC values, respectively, for the ith bin in the kth distribution. The shape-only covariance matrices are provided in the public data release for the ν μ CC1π AE and ν μ CC1π 0 measurements, and the method of Ref. [21] (Sec. 10.6.3) was used to extract them for the ν μ CCNπ AE andν μ CC1π 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 MINERνA, T2K, NOνA, and MicroBooNE. We use the Smith-Moniz relativistic Fermi gas model [22] with an added high momentum tail as per Bodek and Ritchie [23]. The Valencia random phase  12.6 default model predictions compared to MINERνA data. Colors correspond to particle content at the nucleon interaction. "Other" is dominated by coherent pion production. "MC shape" shows the total MC prediction after it has been normalized to match the total data normalization. In the case of the shape-only distributions (θ μ , T π , θ π ) the shape-only χ 2 =N bins values are shown. All cross sections are per nucleon.
approximation screening [24] is applied as a weight to quasielastic events. The two-particle two-hole process is simulated using the Valencia model [11,25]. MINERνA currently uses a modification of v2.8.4 [19,26,27] 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 × 10 6 events were generated using the MINERνA flux predictions [28], a polystyrene target, and the official GENIE 2.12.6 splines [29].
To simulate pion production, GENIE uses the Rein-Sehgal (RS) model [30] with a hadronic invariant mass cut of W ≤ 1.7 GeV. Of the 18 resonances in the RS model, the Δð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. [31]. 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 nonisotropic 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 [32], and the Andreopoulos-Gallagher-Kehayias-Yang (AGKY) model to describe hadronization [33]. In the AGKY model, the Koba-Nielsen-Olesen model [34] is used for W ≤ 2.3 GeV and PYTHIA [35] is used for W ≥ 3.0 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 [36,37], including the effect of lepton masses in the cross-section calculation. MINERνA has found that the RS coherent pion production model needs to be suppressed by ∼50% at T π < 500 MeV to agree with the data [27]. This correction also moves the shape of the T π spectrum closer to the predictions of the Berger-Sehgal coherent model [38]. The ν μ CC1π AE 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 MINERνA, this suppression is included in the analysis presented from Sec. IV onwards.
The "hA Intranuke" effective cascade model [39] is used to model pion and nucleon FSIs. In this model, the effect of intranuclear scattering is parametrized as a single cascade step applied to each particle emanating from the primary interaction. This model steps hadrons through a nucleus of radius r ∼ 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 hadronneutron cross sections [40]. The probability to interact with the nucleus is high; it is, e.g., ∼73% for a pion from an E ν ¼ 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 MINERνA data are shown in Fig. 1. The shape of the p μ distributions agree well with the data for all four measurements. However, the model overestimates the cross section for π AE production, and as a result the χ 2 for the ν μ CC1π AE and ν μ CCNπ AE , given in the fourth column ("Default") of Table II, are large. The model overestimates θ μ below < 5°in the π 0 channels, although it does correctly predict the shape of the θ μ distribution in the π AE channels. The model underestimates the production rate at large θ μ in ν μ CC1π 0 . The shape of the T π distribution is in larger disagreement for ν μ CCNπ AE data than for ν μ CC1π AE . Since the ν μ CCNπ AE distributions are summed over all identified π AE , redistributing kinetic energy between π AE in events with more than one π AE could resolve some of this tension. The π 0 channels are underpredicted at low T π . Finally, GENIE predictions are too high in magnitude at θ π ≈ 50°in both the ν μ CC1π AE and ν μ CCNπ AE channels, and the prediction has the wrong shape in the ν μ CC1π AE channel. Comparisons using the transport theory based Giessen Boltzmann-Uehling-Uhlenbeck (GiBUU) model [41] show similar shape disagreements despite GiBUU's use of an advanced semiclassical cascade model to simulate FSIs [1]. Each of the measurements are shown as MC/data ratio distributions in Fig. 2. Similar comparisons between the MiniBooNE and MINERνA experiments are found in Ref. [5]. The shape-only datasets (θ μ , θ π , T π ) 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 χ 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 paper are summarized in Table III 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 res A , 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, NormRes. Variations in NormRes approximate the behavior of varying F A ð0Þ in the axial form factor in the Rein-Sehgal model. Since low θ μ correlates with low Q 2 , variations in M res A have the largest effect on the shape of the muon angular distributions as shown in Fig. 3, and have a small effect on the θ π spectrum.
Dials are available to vary the normalization of the nonresonant 1π 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 Sec. V, these dials were grouped into a single background scaling for nonresonant 1π production, NonRes1π, following the approach in Refs. [13,14]. A similar treatment was also applied to nonresonant 2π production, NonRes2π, 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 NonRes2π dial introduce a large change in normalization for the ν μ CCNπ AE channel and has a minor effect in the other single pion channels as the fraction of multi-π 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 res A was required to model the low-Q 2 region and a large suppression of the nonresonant π production (−54%) was required to match the observed cross sections of π þ and π 0 production. The reanalysis used the measured ratios of the rates of single π production to charged current quasielastic (CCQE) measurements to cancel errors in the flux. We note that by using CCQE data multiple times, they introduce TABLE III. Summary of the GENIE dials optimized in this paper, their default values, and the uncertainties recommended by the GENIE Collaboration. We do not use the defaults for M res A , NormRes, and NonRes1π 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 res A )  hidden correlations which may have a small effect on the postfit uncertainties. However, as the single pion statistical uncertainties at ANL [42] and BNL [43] were magnitudes higher than the CCQE statistical uncertainty [44,45], the effect was neglected in that work, and is also neglected here. The resulting parameter tunes shown in Table IV and Fig. 5 have been partially adopted by MINERνA and NOνA which both apply the nonresonant rescaling of 43% but leave the other parameters unchanged. Figure 6 shows MINERνA data and the predictions of GENIE when its output has been reweighted to reflect the parameter changes of Table IV. The channel-by-channel contributions to the χ 2 are given in the fifth column ("ANL/BNL") of Table II. Incorporating the parameter changes improves the total normalization agreement in the p μ distributions for ν μ CC1π AE and ν μ CCNπ AE . The χ 2 for the p μ distribution is also improved in theν μ CC1π 0 channel, even though the ANL/BNL data is from neutrino interactions only. The χ 2 for the p μ distribution in the ν μ CC1π 0 channel is somewhat worse as the parameter tunes reduce the predicted nucleon ν μ CC1π 0 cross section. The modification of M res A shifts the θ μ predictions to lower values, increasing the χ 2 contributions. The T π and θ π distributions change mostly by normalization, having a smaller effect on the χ 2 ". The overall agreement of GENIE with MINERνA data is not improved by incorporating the ANL/BNL information. Indeed, the total χ 2 increases, largely because of the χ 2 contributions from the θ μ distributions. GENIE provides a dial that influences the resonances' decay into the pion-nucleon system in the resonance rest frame, π-iso, and allows events to be reweighted continuously between the default anisotropic distribution (π-iso ¼ 0) and the isotropic distribution (π-iso ¼ 1). The Adler angle 2 is highly sensitive to the π-iso parameter and has been measured by neutrino induced pion production experiments on single nucleons, such as ANL [42], BNL [43], BEBC [47,48], and FNAL [49]. Nucleon data strongly prefer an anisotropic process, as shown in Fig. 7. Nonetheless, π-iso has some impact, albeit one that does depend on how FSI are modeled, on the shape of MINERνA θ π and T π distributions, seen at the bottom of Fig. 7, and was therefore included in this work. The GENIE hA model for FSIs has uncertainties from the π − A cross-section data to which the model was tuned. The total π − 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. Figure 6 and Table II show the unsatisfactory agreement of the GENIE prediction against MINERνA 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 MINERνA data. The parameters M res A , NormRes, and NonRes1π are included in the fits with a penalty term added to the χ 2 from the ANL and BNL data. The penalty term uses the covariance, M, shown in Fig. 5:

V. TUNING THE GENIE MODEL
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 disfavored with χ 2 pen ¼ 299.3, but changing NonRes1π to 43% while leaving all the other parameters at their default values reduces the χ 2 pen to 21.8, showing that the largest tension is due to the NonRes1π parameter.
The π-iso dial is allowed to vary in the range 0 to 1 in the fit, corresponding to a continuous variation between a RS angular distribution and an isotropic distribution for Δð1232Þ decay. To avoid the normalization of the ν μ CCNπ AE  measurement pulling parameters in the ν μ CC1π AE model, the NonRes2π dial was allowed to vary between 0% and 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 χ 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 χ 2 penalty terms were added for the FSI dial in either tuning: the parameters  FIG. 6. GENIE ANL/BNL single pion tuning model predictions compared to MINERνA data. The distributions have been weighted to the ANL/BNL tuning parameter set, and have had the coherent pion correction applied. Colors correspond to particle content at the nucleon interaction. "Other" is dominated by coherent pion production. "MC shape" shows the total MC prediction after it has been normalized to match the total data normalization. In the case of the shape-only distributions (θ μ , T μ , θ π ) the shape-only χ 2 =N bins values are shown. All cross sections are per nucleon.
were driven solely by MINERνA data. The charge exchange and pion production dials had small contributions to the overall χ 2 for the selected data, forcing the parameters to be inflated beyond þ3σ of GENIE's recommendation, with large postfit uncertainties. Furthermore, the pion elastic scattering parameter is strongly constrained by external data, so its 1σ variation has a small impact on the MINERνA distributions. The non-FSI fit parameters' (e.g., M res A ) 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 [50] 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 χ 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 V. Figure 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 χ 2 and best-fit parameter values.
Comparing to the results of the ANL and BNL reanalysis, larger values of M res A and smaller values of NormRes were found by the fit, pulling the parameters closer to GENIE nominal. The NonRes1π parameter is strongly bound by the bubble chamber data and the MINERνA data did little to improve on this constraint. The penalty term contributed to the χ 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 postfit correlation matrices are provided in Fig. 10. The ANL/BNL input correlations are largely maintained in our fit.
Tables VI and VII show the results when individual MINERνA 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 χ 2 is steered by the bubble chamber χ 2 penalty. The individual channel fits also found values at the 300% limit for NonRes2π dial, except       ν μ CC1π 0 channels, but not in the ν μ CC1π AE orν μ CC1π 0 channels. Furthermore, the ν μ CC1π 0 shows the strongest χ 2 penalty, indicating tension with the ANL/BNL prior. Given the different kinematic regions covered by the channels (see Table I) 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 θ π distributions. Disagreements in the θ π spectrum are clearly seen in the data/MC ratios of Figs. 8 and 9, and the large χ 2 values observed for the ν μ CCNπ AE and ν μ CC1π 0 channels. The individual χ 2 contributions in the joint tuning best fit, shown in sixth and seventh columns ("FrAbs tune" and "FrInel tune") of Table II, show that not all distributions in all channels benefit from the model variations, as the default GENIE fits have a better χ 2 for some distributions.
In particular, the ν μ CC1π 0 channel distributions have worse agreement after the tuning, with only the θ π distribution improving in χ 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 χ 2 results, Figs. 8 and 9 show that there are still unresolved shape disagreements in both the T π and θ μ kinematics.
The tension between MINERνA's nuclear data and the constraints from the ANL and BNL nucleon data is difficult to confidently pinpoint; the lack of lepton mass effects [51], modification to the resonance propagator in the nucleus [52,53], missing diagrams describing the nonresonant background contributions [9], dynamical coupled channels [54], interactions on correlated initial states, and the pion FSI model [1] 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. Figure 11 (not used for any tuning) shows the Q 2 distributions observed at MINERνA for our tunes. The data are below the predictions of the tunes of Sec. V at low values of Q 2 . There are also differences at low θ μ , as shown in Figs. 8 and 9. Measurements of ν μ CC1π AE and ν μ CC1π 0 interactions on mineral oil at MiniBooNE have shown a data/MC shape discrepancy for the RS implementation in the NUANCE model [55,56] in both Q 2 and cos θ μ distributions [16,57]. In the MINOS quasielastic analysis [15] on iron, which used NeuGen [58], a similar disagreement was observed when studying pion production dominated sidebands. Indeed, concerns about low-Q 2 modeling date back almost a decade [31]. The data from MINOS and MiniBooNE experiments and the MINERνA data on polystyrene 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 θ μ modeling by introducing a Q 2 -dependent correction to the model.
The MINOS Collaboration suppression was expressed as where the free parameters A ¼ 1.010 and Q 0 ¼ 0.156 GeV were empirically extracted from bin-by-bin fits in Q 2 to the data, and a hard cutoff at Q 2 < 0.7 GeV 2 was imposed. We chose an empirical function so that the shape of the suppression preferred by each of the MINERνA 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 three points ðx i ; R i Þ i¼1;2;3 between 0.0 < x < 1.0 and 0.0 < R < 1.0, where x ≡ Q 2 . Motivated the ANL/BNL curves in Fig. 11, the correction is assumed to approach unity as Q 2 approaches 0.7 GeV 2 , providing the constraint ðx 3 ; R 3 Þ ¼ ð0.7 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 GeV 2 ; 1.0): This interpolation function is then used to calculate the correction for each event as 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 GeV 2 so that R 2 describes the curvature at the center point of the correction. Expressing the weights with Eqs. (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 Eq. (7) ensures that wðQ 2 Þ → 1.0 as x → 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 VIII. 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 ν μ CC1π AE and ν μ CCNπ AE 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 MINERνA data when FrAbs (FrInel) is taken as a free parameter. As anticipated, the predictions now have better agreement with the data in regard to the θ μ distribution, and the χ 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 res A and 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. Figure 11 shows the comparison of all our models directly against MINERνA data in Q 2 .    Although the tuning sees improvement in the χ 2 for the ν μ CC1π 0 andν μ CC1π 0 distributions, the ν μ CC1π AE and ν μ CCNπ AE distributions get worse, hinting at tensions in the charged and neutral pion production channels. Tables IX and X show the results of the fits to individual channels, and Table XI shows the breakdown of contributions to the χ 2 from the individual channels. The best-fit χ 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 VIII (or by comparing  Tables VI and IX) 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ν μ CC1π 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 MINERνA data in the ν μ CC1π AE , ν μ CCNπ AE , ν μ CC1π 0 , andν μ CC1π 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 μ , θ μ , T π , and θ π . Parameter fits were performed with either the fraction of pions absorbed or the fraction of pions inelastically scattered in FSIs as a floating parameter, with broadly similar conclusions for the two cases.
The results of the fit (see Table V) 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 MINERνA nuclear target data and the light-target bubble chamber datasets 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 MINERνA 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 VI and VII). Because the four channels cover different kinematic regions (see Table I) 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 MINERνA datasets.
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 θ μ and Q 2 distributions, tensions remain. In particular, fits to individual MINERνA 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 datasets, 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 MINERνA 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.