Angular distributions in Monte Carlo event generation of weak single-pion production

One of the substantial sources of systematic errors in neutrino oscillation experiments that utilize neutrinos from accelerator sources stems from a lack of precision in modeling single-pion production (SPP). Oscillation analyses rely on Monte Carlo event generators (MC), providing theoretical predictions of neutrino interactions on nuclear targets. Pions produced in these processes provide a significant fraction of oscillation signal and background on both elementary scattering and detector simulation levels. Thus, it is of critical importance to develop techniques that will allow us to accommodate state-of-the-art theoretical models describing SPP into MCs. In this work, we investigate various algorithms to implement single-pion production models in Monte Carlo event generators. Based on comparison studies, we propose a novel implementation strategy that combines satisfactory efficiency with high precision in reproducing details of theoretical models predictions, including pion angular distributions. The proposed implementation is model-independent, thereby providing a framework that can include any model for SPP. We have tested the new algorithm with the Ghent Low Energy Model for single-pion production implemented in the NuWro Monte Carlo event generator.


I. INTRODUCTION
Single-pion production (SPP) is one of the main reaction channels relevant for accelerator-based neutrino experiments, where neutrino energies range from a couple of hundred MeVs up to several GeVs [1]. Indeed, in experiments with detectors using Cherenkov radiation, such as T2K [2] and MiniBooNE [3], it is challenging to distinguish neutral pions from electrons. This makes their production the main background for the detection of low-energy electrons. A good understanding of this background is essential for future CP violation measurements in the Hyper-Kamiokande experiment [4] and in attempts to understand the excess of ν e -like events reported by the MiniBooNE collaboration [5]. Moreover, produced pions issue a significant background for other neutrino experiments such as MicroBooNE [6], as it is challenging to distinguish charged pions from muons in Liquid Argon Time Projection Chambers. Regarding oscillation analyses, SPP also contributes to the commonly used CC0π experimental topology [7], provided that the pions get reabsorbed in the nuclear medium or remain otherwise undetected. Furthermore, this interaction channel is itself a part of the signal for oscillation experiments especially with higher-energy neutrino beams such as NOvA [8] and DUNE [9], but also for T2K [10].
Over the past couple of years, the MINERvA, T2K, ArgoNeuT, and MiniBooNE experiments [11][12][13][14][15][16] have collected an increasingly large dataset for (anti-)neutrino-induced single-pion production on nu- * kajetan.niewczas@uwr.edu.pl † alexis.nikolakopoulos@ugent.be clear targets. Subsequently, it has been compared to predictions from several models, revealing significant differences in their description of the data. Moreover, there are apparent tensions between the MiniBooNE, T2K, and MINERvA SPP measurements [17][18][19][20] themselves. Ref. [21] showed that a simultaneous agreement between the results of the ANL and BNL bubble chamber data and the MINERvA experiment could not be reached. Furthermore, it was not possible to provide a consistent description using a single parameter tune for the different SPP channels measured by the latter.
The use of nuclear targets in neutrino oscillation experiments considerably complicates the description of single-pion production because the presence of such a medium affects all of the hadrons in the process. On top of that, final-state interactions (FSI), such as pion absorption or charge exchange pion-nucleon scattering, alter the experimental signal entirely. It is seemingly an intractable problem to provide a detailed microscopic description of FSI over the sizeable phase space of these experiments. For this reason, the FSI are usually treated in an approximate way using intranuclear cascade models [22][23][24] implemented in various Monte Carlo neutrino event generators (MC). An exception is GiBUU that solves the coupled Boltzmann-Uehling-Uhlenbeck (BUU) transport equations instead [25].
A prerequisite for a good description of neutrino-induced single-pion production on nuclei in the factorized approach used in MCs is an accurate model for such scattering off the nucleon. Several models, with varying regions of applicability, have been developed for neutrino-induced SPP off the nucleon [26][27][28][29][30][31][32][33][34][35][36][37][38][39]. However, these models have not readily found their way into Monte Carlo event generators, and if so, without accounting for their full kinematic complexity.
In this work, we perform a detailed study of possible strategies to implement single-pion production models in neutrino event generators. Based on the results of this study, we propose a novel algorithm for the case of SPP on the nucleon target to allow for further progress in the accommodation of information from recent experimental measurements. The algorithm is model-independent, as it only relies on the kinematics of the process and ensures no relevant information is lost on a neutrino-nucleon interaction level. Such a solution allows for any theoretical model to be implemented in MCs, facilitating a comparison of different approaches. Additionally, owing to the separation of the leptonic and hadronic currents, it provides flexibility to modify the former, e.g., with Beyond the Standard Model physics. Furthermore, with appropriately implemented hadronic currents, one can calculate cross sections for charged current, neutral current, and electron-induced SPP with a consistent treatment of both the vector and axial components. We claim that the proposed algorithm will be of great importance for future implementations of neutrino-nucleus single-pion production, dealing with a considerable number of degrees of freedom and hence a critical demand to maintain both numerical efficiency and precision. This paper is structured as follows. In Sec. II, we review the kinematics of lepton-induced single-pion production. Then, in Sec. III, details of the new implementation of SPP in Monte Carlo event generators, as well as the particular numerical tools used, are described. In Sec. IV, we report the results of our study in the context of the implementation performance and physics outcomes. In the last section, Sec. V, we present our conclusions.

II. KINEMATICS AND CROSS SECTION
We commence by describing the kinematics of lepton-induced single-pion production, where an incoming lepton with four-momentum k = (E, k) scatters off a nucleon p i by exchange of a single gauge boson with four-momentum q = (ω, q), thereby producing a pion. We denote the four-momenta of the final-state lepton, pion, and recoiling nucleon by k , k π , and p N , and their rest masses by m, M π , and M N , respectively. It is convenient to describe such a process in the hadronic center-of-momentum system (CMS), with the lepton plane defining the x-z plane and the direction of the momentum transfer q defining the z-axis, as depicted in Fig. 1. In the hadronic CMS, for which we denote quantities with a superscript * , the final hadronic system is at rest, meaning k * π = − p * N . We characterize the kinematics by the Lorentz invariants: the invariant hadronic mass W 2 = (q + p i ) 2 = (k π + p N ) 2 and the exchanged four-momentum squared Q 2 = −q 2 = −(k − k ) 2 , along with the produced pion solid angle Ω * π . Within the Born approximation, we can describe the cross section as a contraction of the leptonic and hadronic tensors. The standard calculation of the leptonic tensor for massless incoming leptons yields where η is the metric tensor with signature (+, −, −, −), µναβ is the antisymmetric Levi-Civita tensor ( 0123 = +1), and h is the helicity of the incoming lepton. We define the hadronic tensor as with J µ the hadronic current, and averaging and summation over the spin of the initial and final nucleon are assumed. Respresenting the hadronic current in terms of initial and final state nucleon spinors and the transition operator O µ as one obtains for the hadronic tensor where With these definitions, the cross section is where the coupling constant for the charged current case that we consider in this paper is Using the invariance of the leptonic tensor under rotations of the hadronic plane around q, one can factorize the dependence of the cross section on the azimuthal angle in terms of trigonometric functions, as shown explicitly in Refs. [27,[40][41][42]. Specifically, in the given CMS, we express the cross section as where the functions A, ..., E do not depend on the azimuthal pion angle φ * π . Below, we write them explicitly, in terms of the elements of the leptonic and hadronic tensors, computed for the kinematics of Fig. 1 (with φ * π = 0), and making use of the symmetry properties of L µν , as where H s and H a correspond to the symmetric (real) and antisymmetric (imaginary) parts of the hadronic tensor: For antineutrino interactions, the terms including the imaginary part of the hadronic tensor change sign as all of the off-diagonal terms of the leptonic tensor involving a Lorentz index 2 are purely antisymmetric and proportional to the helicity, while the others are symmetric.
In the context of this work, it is essential to notice that the double-differential cross section d 2 σ/dW dQ 2 and the triple-differential d 3 σ/dW dQ 2 d cos θ * π are entirely determined by the function A as the other contributions disappear after integration over the azimuthal pion angle φ * π . We remark that the presented expressions apply to all electroweak SPP processes, thereby facilitating a consistent treatment of the vector current across electron-and neutrino-induced cases implemented in event generators. Furthermore, a similar separation of the angular dependence is valid for one-nucleon knock-out on a nuclear target or for any semi-leptonic process in which a single on-shell particle defines the hadronic plane for that matter. Thus, similar methods as the ones outlined in the next section should apply to the implementation of microscopic models for exclusive one-nucleon knock-out.

III. MONTE CARLO EVENT GENERATION FOR SINGLE-PION PRODUCTION
The kinematics for weak single-pion production off the nucleon given an incoming neutrino energy, the target nucleon momentum, and an arbitrarily chosen lepton scattering plane, is fully described by four independent variables. In what follows, these quantities are considered to be random variables with a probability distribution defined by Eq. 7. While constructing Monte Carlo event generators, one of the major tasks is to generate these variables efficiently. Here, we discuss several of our approaches, each of them presenting a different trade-off between efficiency, precision, and reliance on precomputed assets.

A. 4D algorithm
The most straightforward approach is to use directly the full cross section formula (7). The available phase space of the independent variables W , where s = (k + p i ) 2 , and the underline marks the quantities calculated in the lepton+hadron center-of-momentum frame: In this approach, for each event, we sample four independent variables, adding a randomly selected lepton scattering plane. We perform the sampling following the order presented in Eq. 14, starting from W , because the range in Q 2 is W -dependent and because one needs both W and Q 2 to specify the hadronic CMS needed to select cos θ * π and φ * π . Such information is enough to generate the full kinematics of an event trial. Each of them has an assigned event weight, given by Eq. 7 multiplied by the Monte Carlo phase space factor The average value of this weight is equal to the total cross section. We obtain the final set of events by applying the accept-reject algorithm on the collection of trials. We will refer to this strategy of generating events as the "4D algorithm".
Although asymptotically correct, we expect this approach to be inefficient, especially with increasing neutrino energies. The efficiency of an accept-reject algorithm depends on the interplay between the distribution's shape and the sampling envelope. Since the former is initially unknown, we choose the latter to be the maximal value of the cross section (Eq. 7), calculated in real-time. As we increase the phase space, we access new regions of low, relative to the envelope, cross section values, leading to event trials with a minimal chance of acceptance. The interplay between the acceptance efficiency and the computation time needed to calculate an event weight are the main features contrasting the proposed algorithms. In the 4D algorithm, for every event trial, we evaluate the value of the cross section given by Eq. 7 once.

B. 3D algorithm
In the second approach, we isolate the dependence of the cross section on the azimuthal pion angle. After performing an integration over φ * π , the differential cross section depends only on the function A, and its explicit dependence on H µν reads where the hadronic tensor elements are functions of three variables: W , Q 2 , cos θ * π . In this case, only three variables are sampled and we attribute event trials with weights obtained from multiplying the results of Eq. 17 by the new Monte Carlo phase space volume As before, the average of the event weights yields the total cross section. We obtain the final set of, yet incomplete, events using the same accept-reject method, with the sampling envelope given by the maximum of Eq. 17. Due to the reduced phase space dimensionality, the accept-reject algorithm for incomplete events, without an assigned value of the pion azimuthal angle, is more efficient.
For already selected events, we sample the variable φ * π using the known probability distribution given, for fixed values of W , Q 2 , cos θ * π , by and its cumulative distribution function: As the derivative of the F (φ * π ) function is known algebraically, its inversion with the Newton method is efficient and converges rapidly. In what follows, we will call this procedure the "3D algorithm".

C. 2D algorithm
The starting point for this approach is the formula obtained from Eq. 17 by integrating out the cos θ * π variable and adopting the notation: As a result, we express the double-differential cross section in terms of 5 combinations of hadronic tensor elements, which depend solely on W and Q 2 . We store their values in the form of lightweight tables. The first step of the '2D algorithm' is to sample a pair of variables (W, Q 2 ) with the probability density defined by Eq. 21. We perform it efficiently, using the precalculated tables with a suitable bilinear interpolation. At this point, we build an incomplete event trial and compute its weight, analogously to the previous approaches, by multiplying the values obtained from Eq. 21 by the Monte Carlo phase space factor We accept the set of incomplete events according to their weights, relative to the maximum of Eq. 21, and only then we assign the values of cos θ * π and φ * π . Such an approach saves a considerable amount of time, avoiding the computation of a full event before applying the accept-reject algorithm.
We select the value of cos θ * π using a probability distribution governed by the function A. To optimize this task, we exploit the smooth character of this function in the region of interest. Having W and Q 2 fixed, we calculate the values of A(cos θ * π ) at k points and approximate as a polynomial of degree k − 1. Then, we obtain the cumulative distribution function as a polynomial of degree k and sample the cos θ * π variable using the inverse sampling method. For (k ≤ 3), we perform the inversion algebraically, while for larger k, numerically, with the bisection method. We have checked that, for most kinematics, the degree of k = 3 provides sufficient precision, while the distributions are almost exact on the whole phase space for degrees k ≥ 7. We will discuss the choice of the optimal value of k in the next section.
Depending on the implementation effort and allowed memory, it is also possible to store in tables hadronic tensor elements H µν (W, Q 2 , cos θ * π ) that allow obtaining the full A(W, Q 2 , cos θ * π ) function. Then, in each event, the maximum of the cos θ * π distribution is given explicitly, and one can sample its value using the accept-reject method. Such an approach enables us to reduce the time-consumption of each trial event further. In what follows, we will denote this approach as the "2D algorithm (table)".
Finally, to finish building the kinematics for the accepted events, we need to sample the variable φ * π . We proceed by repeating the method described in Sec. III B.

D. Numerical tools
To reliably test the performance of the abovementioned sampling algorithms, we performed simulations using the Ghent Low Energy Model (LEM) of single-pion production implemented in the NuWro Monte Carlo event generator. The particular implementation works on a restricted phase space defined by the condition W < 1.5 GeV.

Ghent Low Energy Model of SPP
This single-pion production model is based on the work of Hernández, Nieves, and Valverde (HNV), first presented in Ref. [30] with later improvements of Refs. [43][44][45]. It contains a microscopic description of the SPP at the amplitude level and includes, in addition to the contributions from the ∆(1232) and D 13 (1520) resonances (both direct and crossed channel Feynman diagrams), the lowest-order background diagrams derived from chiral perturbation theory (ChPT). Additionally, it includes a relative phase between the ChPT terms and the dominant partial wave of the ∆-pole, which partially restores unitarity [44].
The Ghent LEM [39] is a custom variant of the model with an independently written code. On top of the standard version, it includes additional s-and u-channel contributions from the spin-1 / 2 resonances P 11 (1440) and S 11 (1535) [46]. Additionally, the same model, working in the relativistic plane wave impulse approximation, was extended to describe neutrino scattering on nuclei [18,19].

NuWro Monte Carlo event generator
NuWro is a versatile Monte Carlo neutrino event generator, which has been developed by the theoretical group of the University of Wroc law since 2005. It is applicable for simulations in the range of neutrino energies covered by the accelerator-based neutrino oscillation experiments, with an upper bound of ∼ 100 GeV. NuWro supports quasielastic, single-pion production, and more inelastic channels (DIS) of neutrino scattering off free nucleons. The neutrino-nucleus interactions are modeled with various nuclear models (e.g., global or local Fermi gas, spectral functions [47,48], or a momentum-dependent nuclear potential [49]) in the impulse approximation, succeeded by final-state interactions of outgoing hadrons simulated using an intranuclear cascade model [22,23]. Moreover, the inclusion of complex nuclear targets enables additional interaction channels such as two-body current processes, coherent pion production, and neutrino scattering off atomic electrons [50]. The code used in this work bases on NuWro version 19.02.2 [51].
The NuWro single-pion production model combines the contribution from the ∆(1232) resonance excitation [52] with a non-resonant background obtained by extrapolating the DIS contribution to lower values of W , blended incoherently in the region W ∈ (1.3, 1.6) GeV [53]. The generated events follow the double-differential cross sections d 2 σ/dW dQ 2 for both the resonant and non-resonant parts.
On top of that, the model obtains the Ω * π distributions using the parametrized ones measured by the BNL bubble chamber experiment [54] for the former, while for the latter, obtains the kinematics using the PYTHIA6 hadronization routines [55]. Alternatively, one can use the parametrization obtained by the ANL experiment [56]. In this work, we refer to this model as "isobar NuWro".
In NuWro, for all of the described single-pion production model implementations, we apply additional optimizations of sampling in the (W, Q 2 ) plane. As it is common for all models and methods presented in this work, this has no impact on our findings nor conclusions.

A. Performance
We have implemented the Ghent LEM in NuWro, applying the five versions of the strategies presented in Sec. III, labeled: "4D alg.", "3D alg.", "2D alg. (k = 3)", "2D alg. (k = 7)", "2D alg. (table)". We summarize their performance in Table I, with four numerical computations: for two neutrino energies E = 1.0, 2.5 GeV, and off both proton and neutron nucleon targets. In the respective columns of these tables, one can find: an average weight σ and its standard deviation s 1M calculated from 1 million trial events, computer time τ needed to   calculate a trial event (before the accept-reject algorithm is applied) in arbitrary units, efficiency of the acceptreject algorithm, and the relative increase of computer time α needed to generate an event with complete kinematics. In the last columns, we present values of S 1M that is a measure of the performance of a given algorithm: an estimate of the time needed to produce a sample of N = 1 × 10 6 events. In a given simulation of efficiency , one has to generate N/ trial events, out of which N events are accepted and require the complete kinematics, while N/ − N are the rejected trial events that require only the weight calculation. Since the computation of a trial event takes time τ and of a complete event τ (1 + α), the overall computer time needed to generate a set of N events becomes Thus, the value of S N depends on three variables: τ , , α, that fully characterize each algorithm. The first, most significant difference between various approaches appears in the values of τ and show that, in the models used, the time needed for generating a trial event is ∼ 20 times smaller for the 2D algorithm off protons, relative to the 3D and 4D algorithms, while off neutrons the difference rises about twice as much. The former stems solely from the computational cost needed to evaluate the hadronic tensor, which is the bottleneck of the Ghent LEM, while the latter comes from the fact that neutrino-induced SPP off the neutron involves two possible final states and both cross sections need to be evaluated to obtain the weight of any of those events.
The differences in for simulations with the same conditions come from the differential cross section shapes as well as the size of the sampled phase space that grows with increasing energy. Due to the dominance of the ∆ ++ resonance, the SPP cross sections for neutrino scattering off the proton target are much more peaked, leading to lower event acceptance efficiency. On the other hand, the differences in efficiencies between particular algorithms within the same simulations come from different dimensionalities of the sampled phase spaces and the fact that the cross sections are not uniform in the additional variables (cos θ * π , φ * π ). Values of the third characteristic variable α represent all of the secondary effort, relative to the event trial computation time, needed to generate the full kinematics of an accepted event. One can see that for the 3D algorithm, in which sampling of the φ * π variable requires to compute the hadronic tensor one additional time, relative to the 4D method, α equals 1.0 and 0.5 for the proton and neutron targets, respectively. The 2D algorithm methods, on top of the φ * π sampling, require additional effort to assign the cos θ * π variable. The increase in α while going from the 2D (table) method to the ones that use polynomial interpolation is almost proportional to the number of times (k = 3, 7, ...) we calculate the hadronic tensor. We expect that one can avoid such behavior using a model implementation that separates the angular dependence algebraically, e.g., in a partial wave expansion, where one can compute the hadronic tensor for different values of cos θ * π at fixed values of Q 2 and W in a much shorter time. However, in general, the cos θ * π dependence is not a priori known. Thus, in this study, we opted to present the most model-independent case.
The resultant performance of all the optimization methods in reducing the total simulation time S N is notable. Considering its execution time and susceptibility to the investigated factors, we conclude that the "2D alg. (table)" method performs best, and in what follows, we use it to generate all of the Monte Carlo simulation results. To strengthen this reasoning, we emphasize that in actual simulations, there is an additional, global computational effort needed to specify the weight of particular interaction channels and initialize the event sampling envelope. In NuWro, we know it as generating test events that require solely an event weight calculation, which is less demanding using the 2D algorithms.

B. Inclusive cross section
To illustrate the accuracy of the implementation of the single-pion production model in NuWro within the "2D alg. (table)" framework, we show several comparisons with the exact results obtained with the original Ghent LEM code. For every presented plot, we compute the Monte Carlo results by averaging over six simulations with 10M events across the whole phase space. The additional band represents a 1σ error on the average. In Fig. 2, we compare the results for the inclusive cross sections as a function of W at fixed Q 2 for electron neutrinos and antineutrinos with an energy of E = 1 GeV, including all possible single-pion production channels. For each value of Q 2 = 0.1, 0.5 GeV 2 , we gathered Monte Carlo events in bins with a width of ∆Q 2 = 0.01 GeV 2 and ∆W = 5 MeV. One can see that the "2D alg. (table)" method provides excellent accuracy. The statistical uncertainty on its results is the smallest for (anti)neutrino reactions on the (neutron)proton, as these are cases with only a single SPP channel accessible. For the other target/helicity combinations, the simulations split the events over two final states, with the one of the higher cross section receiving a larger share, which is reflected in the uncertainty.

C. Angular distributions of the pion
The main strength of the presented approach is the exact implementation of the outgoing pion angular distributions. To illustrate this, in Fig. 3, we plot the cross sections as a function of cos θ * π for values of W = 1230, 1270, 1310 MeV and fixed Q 2 = 0.1 GeV 2 with incoming (anti)neutrino energy E = 1 GeV. We obtained NuWro results in the same way as described in Sec. IV B. Here, we gathered events in bins of ∆Q 2 = 0.01 GeV 2 , ∆W = 5 MeV, and ∆ cos θ * π = 0.04. The obtained Monte Carlo results precisely reproduce the exact model calculations. The shape of the cos θ * π distribution varies with both the interaction channel and kinematics. This behavior is in contrast to the commonly used approach in which the angular dependence of the outgoing pion-nucleon pair is described isotropically or by a distribution independent of the kinematics.
The next comparison, in Fig. 4, concerns the single-differential cross sections as a function of φ * π for electron neutrinos and antineutrinos with an energy of E = 1 GeV, averaged over ∆φ * π = π/25 rad bins. Since the procedure for sampling φ * π is practically exact, shapes of these distributions exemplify the total numerical error propagating from the bilinear and trilinear interpolation of the tabularized information used to sample the values of (W , Q 2 ) and cos θ * π , respectively. Hence, one can interpret this comparison as a good measure of the full accuracy of the proposed algorithm. Regarding the physical results themselves, one immediately notices the asymmetry of dσ/dφ * π around φ * π = π, corresponding to pions produced above or below the lepton scattering plane. As seen in Eq. 7, the D and E functions, which give contributions proportional to sin (φ * π ) and sin (2φ * π ), respectively, are responsible for such behavior. As explained thoroughly in Ref. [42], these asymmetries emerge from relative phase differences between the distinct contributions to the amplitude. Hence, they are not present in models that are described by incoherent sums of resonances, or resonance and background contributions. In the Ghent Low Energy Model, the asymmetry can only arise from the interference between the imaginary part of the resonance propagator (plus the Olsson phases in the case of ∆) and the non-resonant background. Such an asymmetry is also not present in unpolarized electron scattering because both the D function, with the vector-vector contribution proportional to the polarization, and the E function, being a purely vector-axial interference term, disappear in that case.
In Figs. 5 and 6, we show the full two-dimensional Ω * π dependence in the different electron (anti)neutrino-induced SPP channels with E = 1 GeV for fixed Q 2 = 0.1 GeV 2 and W = 1230 MeV, i.e., at the ∆(1232) peak. Here, we average the Monte Carlo results over ∆Q 2 = 0.01 GeV 2 , ∆W = 5 MeV, ∆ cos θ * π = 0.04, and ∆φ * π = π/25 rad bins. Although we performed these NuWro simulations again in the same way as described in Sec. IV B, it is challenging to produce a sufficiently large sample of events to reduce statistical fluctuations. Still, the agreement we find is remarkably good. Analyzing the presented distributions, one can see that the ν(p, pπ + ) and ν(n, nπ − ) cross sections are roughly symmetric with respect to φ * π . These interaction channels only allow isospin 3/2 contributions in the s-channel and are thus dominated ν ‾ e + p -> e + + p + π -FIG. 4: Single-differential cross sections for the ν e -and ν e -induced single-pion production processes as a function of φ * π with incoming energy E = 1 GeV. Solid lines show the Ghent LEM results, while the dashed ones are results of the "2D algorithm (tables)" method implemented in NuWro.
by the ∆(1232) resonance, with minimal impact from the background and thereby minimal interference to generate the asymmetry. The other channels, however, do show a more asymmetric shape as the background contribution grows in relative importance. In Fig. 7, we also present the Ω * π dependence of the "isobar NuWro" model for the same kinematical setup. This model uses angular distributions from the BNL parametrization of Ref. [54], as implemented in NuWro 19.02.2. In the case of neutrino-induced charged pion production off the proton these results are similar to the Ghent LEM, while for the other reaction channels they are quite different. Such behavior originates from the fact that the BNL (ANL) parametrization is obtained from data for the former reaction in the ∆(1232) region. This comparison illustrates that a straightforward application of the same angular distribution to other reaction channels and other phase space regions should be avoided.
Finally, in Fig. 8, we show a shape-only comparison with the single-differential cross sections dσ/d cos θ * π and dσ/dφ * π measured by the ANL [56] and BNL [54] experiments. To obtain this, we performed simulations with the E ANL ∈ (0.2, 6.1) GeV and E BNL ∈ (0.1, 7.5) GeV muon neutrino fluxes off a proton target, and applied a cut on the invariant hadronic mass W < 1.4 GeV.
The theoretical results provide a good agreement for the dσ/d cos θ * π differential cross section, especially in the BNL case, while due to the lack of statistics, the dσ/dφ * π results are not conclusive.

V. CONCLUSIONS
The upcoming precision era of neutrino oscillation experiments requires significant improvements in detecting and modeling more exclusive observables of neutrino scattering. In the case of pion production, a pressing issue is that of their angular distributions. Due to the lack of available data to constrain these quantities, it is of great importance to equip Monte Carlo neutrino event generators with predictions of the most sophisticated theoretical models available. We have made an important step towards this goal, focusing on the (anti)neutrino-induced single-pion production on the nucleon and implementing the Ghent Low Energy Model of Ref. [39] into the NuWro Monte Carlo event generator.
To this end, we investigated various general, model-independent and efficient implementations based on the separation of the pion angular dependence in the hadronic center-of-momentum reference frame. They originate from the idea to sample particular independent kinematic variables needed to build a Monte Carlo event in a specific order, using increasingly differential cross section formulas. In the consecutive steps, such an approach allows performing the time-consuming microscopic model computations only for accepted events and exploits differences in Monte Carlo event generation efficiencies. All of our algorithms start from sampling the (W , Q 2 ) phase space to be able to specify the hadronic CMS. Then, we use approximations that allow us to efficiently choose the value of cos θ * π as well as to exploit the algebraic dependence of the cross section on φ * π . Such approaches provide the flexibility of choosing a different trade-off between efficiency, precision, and reliance on precomputed assets.
To quantify the performance of our implementations, we performed various simulations and measured characteristic quantities for each solution. We conclude that the method labeled "2D alg. (table)", that exploits all the optimizations and mostly relies on precomputed assets, is the most effective solution. We exhaustively checked its accuracy in investigating different multiple-differential cross sections, each time obtaining excellent agreement. The performance of the "2D alg. (k ≥ 3)" methods, which employ a polynomial fit of degree k − 1 for the cos θ * π probability distributions, was intermediate and strongly relied on the choice of the degree k. Although we obtained promising results in the ∆-region for k-values as low as 3, it was necessary to perform simulations with degree k ≥ 7 to reproduce the Ghent LEM predictions over the entire investigated phase space to the percent level. Still, such an approach proved useful, as such a high level of accuracy is not necessary for flux-averaged distributions, and it does not require to precompute nor store the additional tables.
Finally, we compared the new implementation of the Ghent LEM model with the SPP angular distribution data from ANL and BNL bubble chamber experiments and with the results of the nominal single-pion production model of NuWro. We concluded that it is not feasible to use experimental parametrizations for neutrino-induced SPP off the proton in the ∆(1232) region for all channels across the whole phase space. It is of great importance to design Monte Carlo event generators able to provide reliable predictions for such observables.
This work facilitates further studies of nuclear effects in SPP as we can implement more sophisticated models in NuWro. The next important step of this research will be the extension of this implementation framework to single-pion production on the nucleus and investiga-dσ / dW dQ 2 dcos π * dφ π * [10 -47 cm 2 / GeV / GeV 2 / rad] 2D alg. (table) -1 -0.  Yield 2D alg. (table) ANL ν µ + p -> µ -+ p + π + FIG. 8: Pion angular distributions for the neutrino-induced single-pion production on the proton as a function of cos θ * π or φ * π , and data from the ANL [56] and BNL [54] bubble chamber experiments. Solid lines show the Ghent LEM results, while the bins are the results of the "2D algorithm (tables)" method implemented in NuWro. We normalize our cross section predictions to the total experimental yield.